Merge new AOI methods into full mode workflow (#94)

Fixes #36 

* Add workflow for removing aoi losses and calculating q_absorbed

* The engine needs to fit vf calculator as well

* Speed up calculation by ignoring "null" surfaces

* ~2x speed up

* Add "is_empty" property to surfaces

* Speed up vf methods by ignoring null surfaces

* Use "is_empty" flag in aoi methods as well

* Test engine calculated absorbed values with diffuse AOI losses

* Check that VF calculator uses reflectivity matrix for diffuse faoi

* Test for pv engine using faoi functions

* Update both run functions to accept vf calculator params

* run_timeseries_engine works but problem with run_timeseries_parallel

* Fix so that can run faoi_functions with parallel mode

* Try to fix error in Python 2

pvfactors/engine.py

@@ -164,6 +164,9 @@ def fit(self, timestamps, DNI, DHI, solar_zenith, solar_azimuth,
         # Add timeseries irradiance results to pvarray
         self.irradiance.transform(self.pvarray)
 
+        # Fit VF calculator
+        self.vf_calculator.fit(self.n_points)
+
         # Skip timesteps when:
         #    - solar zenith > 90, ie the sun is down
         #    - DNI or DHI is negative, which does not make sense
@@ -193,17 +196,18 @@ def run_full_mode(self, fn_build_report=None):
 
         # Get the irradiance modeling matrices
         # shape = n_surfaces, n_timesteps
-        irradiance_mat, _, invrho_mat, _ = \
+        irradiance_mat, rho_mat, invrho_mat, _ = \
             self.irradiance.get_full_ts_modeling_vectors(pvarray)
 
-        # Calculate view factors
+        # --- Calculate view factors
         # shape = n_surfaces, n_surfaces, n_timesteps
         ts_vf_matrix = self.vf_calculator.build_ts_vf_matrix(pvarray)
         pvarray.ts_vf_matrix = ts_vf_matrix
         # Reshape for broadcasting and inverting
         # shape = n_timesteps, n_surfaces, n_surfaces
         ts_vf_matrix_reshaped = np.moveaxis(ts_vf_matrix, -1, 0)
 
+        # --- Solve mathematical problem
         # Build matrix of inverse reflectivities
         # shape = n_surfaces, n_surfaces
         invrho_mat = np.diag(invrho_mat[:, 0])
@@ -220,18 +224,32 @@ def run_full_mode(self, fn_build_report=None):
         # shape = n_surfaces, n_timesteps
         qinc = np.dot(invrho_mat, q0)
 
-        # Derive other irradiance terms
+        # --- Derive other irradiance terms
         # shape = n_surfaces, n_timesteps
         isotropic_mat = ts_vf_matrix[:-1, -1, :] * irradiance_mat[-1, :]
         reflection_mat = qinc[:-1, :] - irradiance_mat[:-1, :] - isotropic_mat
 
-        # Update surfaces with values: the list is ordered by index
+        # --- Calculate AOI losses and absorbed irradiance
+        rho_mat = np.tile(rho_mat[:, 0], (rho_mat.shape[0], 1)).T
+        # shape [n_surfaces + 1, n_surfaces + 1, n_timestamps]
+        vf_aoi_matrix = (self.vf_calculator
+                         .build_ts_vf_aoi_matrix(pvarray, rho_mat))
+        pvarray.ts_vf_aoi_matrix = vf_aoi_matrix
+        # shape [n_surfaces, n_surfaces]
+        irradiance_abs_mat = (
+            self.irradiance.get_summed_components(pvarray, absorbed=True))
+        # Calculate absorbed irradiance
+        qabs = (np.einsum('ijk,jk->ik', vf_aoi_matrix, q0)[:-1, :]
+                + irradiance_abs_mat)
+
+        # --- Update surfaces with values: the lists are ordered by index
         for idx_surf, ts_surface in enumerate(pvarray.all_ts_surfaces):
             ts_surface.update_params(
                 {'q0': q0[idx_surf, :],
                  'qinc': qinc[idx_surf, :],
                  'isotropic': isotropic_mat[idx_surf, :],
-                 'reflection': reflection_mat[idx_surf, :]})
+                 'reflection': reflection_mat[idx_surf, :],
+                 'qabs': qabs[idx_surf, :]})
 
         # Return report if function was passed
         report = None if fn_build_report is None else fn_build_report(pvarray)

pvfactors/geometry/pvarray.py

@@ -69,6 +69,7 @@ def __init__(self, axis_azimuth=None, gcr=None, pvrow_height=None,
 
         # These attributes will be updated by the engine
         self.ts_vf_matrix = None
+        self.ts_vf_aoi_matrix = None
 
     @classmethod
     def init_from_dict(cls, pvarray_params, param_names=None):

pvfactors/geometry/timeseries.py

@@ -249,6 +249,12 @@ def u_vector(self):
                     np.array([-self.n_vector[1, :], self.n_vector[0, :]]))
         return u_vector
 
+    @property
+    def is_empty(self):
+        """Check if surface is "empty" by checking if its length is always
+        zero"""
+        return np.nansum(self.length) < DISTANCE_TOLERANCE
+
 
 class TsLineCoords(object):
     """Timeseries line coordinates class: will provide a helpful shapely-like

pvfactors/irradiance/base.py

@@ -118,7 +118,7 @@ def get_summed_components(self, pvarray, absorbed=True):
             for component in list_components:
                 value += ts_surface.get_param(component)
             irradiance_mat.append(value)
-        return irradiance_mat
+        return np.array(irradiance_mat)
 
     def update_ts_surface_sky_term(self, ts_surface, name_sky_term='sky_term'):
         """Update the 'sky_term' parameter of a timeseries surface.

pvfactors/run.py

@@ -25,6 +25,7 @@ def run_timeseries_engine(fn_build_report, pvarray_parameters,
                           fast_mode_pvrow_index=None,
                           fast_mode_segment_index=None,
                           irradiance_model_params=None,
+                          vf_calculator_params=None,
                           ghi=None):
     """Run timeseries simulation without multiprocessing. This is the
     functional approach to the :py:class:`~pvfactors.engine.PVEngine` class.
@@ -77,6 +78,9 @@ def run_timeseries_engine(fn_build_report, pvarray_parameters,
     irradiance_model_params : dict, optional
         Dictionary of parameters that will be passed to the irradiance model
         class as kwargs at instantiation (Default = None)
+    vf_calculator_params : dict, optional
+        Dictionary of parameters that will be passed to the VF calculator
+        class as kwargs at instantiation (Default = None)
     ghi : array-like, optional
         Global horizontal irradiance values [W/m2] (Default = None)
 
@@ -87,18 +91,17 @@ class as kwargs at instantiation (Default = None)
         function
     """
 
-    # Prepare irradiance model inputs
-    irradiance_model_params = ({} if irradiance_model_params is None
-                               else irradiance_model_params)
+    # Prepare input parameters
+    irradiance_model_params = irradiance_model_params or {}
+    vf_calculator_params = vf_calculator_params or {}
     # Instantiate classes and engine
     irradiance_model = cls_irradiance(**irradiance_model_params)
+    vf_calculator = cls_vf(**vf_calculator_params)
     pvarray = cls_pvarray.init_from_dict(pvarray_parameters)
-    vf_calculator = cls_vf()
     eng = cls_engine(pvarray, irradiance_model=irradiance_model,
                      vf_calculator=vf_calculator,
                      fast_mode_pvrow_index=fast_mode_pvrow_index,
                      fast_mode_segment_index=fast_mode_segment_index)
-
     # Fit engine
     eng.fit(timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt,
             surface_azimuth, albedo, ghi=ghi)
@@ -118,7 +121,8 @@ def run_parallel_engine(report_builder, pvarray_parameters,
                         cls_irradiance=HybridPerezOrdered,
                         cls_vf=VFCalculator, fast_mode_pvrow_index=None,
                         fast_mode_segment_index=None,
-                        irradiance_model_params=None, n_processes=2,
+                        irradiance_model_params=None,
+                        vf_calculator_params=None, n_processes=2,
                         ghi=None):
     """Run timeseries simulation using multiprocessing. Here, instead of a
     function that will build the report, the users will need to pass a class
@@ -173,6 +177,9 @@ def run_parallel_engine(report_builder, pvarray_parameters,
     irradiance_model_params : dict, optional
         Dictionary of parameters that will be passed to the irradiance model
         class as kwargs at instantiation (Default = None)
+    vf_calculator_params : dict, optional
+        Dictionary of parameters that will be passed to the VF calculator
+        class as kwargs at instantiation (Default = None)
     n_processes : int, optional
         Number of parallel processes to run for the calculation (Default = 2)
     ghi : array-like, optional
@@ -192,9 +199,9 @@ class (or object)
     if np.isscalar(albedo):
         albedo = albedo * np.ones(len(dni))
 
-    # Prepare irradiance model inputs
-    irradiance_model_params = ({} if irradiance_model_params is None
-                               else irradiance_model_params)
+    # Prepare class input parameters
+    irradiance_model_params = irradiance_model_params or {}
+    vf_calculator_params = vf_calculator_params or {}
 
     # Fix: np.array_split doesn't work well on pd.DatetimeIndex objects
     if isinstance(timestamps, pd.DatetimeIndex):
@@ -216,6 +223,7 @@ class (or object)
     folds_fast_mode_pvrow_index = [fast_mode_pvrow_index] * n_processes
     folds_fast_mode_segment_index = [fast_mode_segment_index] * n_processes
     folds_irradiance_model_params = [irradiance_model_params] * n_processes
+    folds_vf_calculator_params = [vf_calculator_params] * n_processes
     folds_ghi = ([ghi] * n_processes if ghi is None else
                  np.array_split(ghi, n_processes))
     report_indices = list(range(n_processes))
@@ -228,7 +236,8 @@ class (or object)
                        folds_cls_engine, folds_cls_irradiance, folds_cls_vf,
                        folds_fast_mode_pvrow_index,
                        folds_fast_mode_segment_index,
-                       folds_irradiance_model_params, folds_ghi,
+                       folds_irradiance_model_params,
+                       folds_vf_calculator_params, folds_ghi,
                        report_indices))
 
     # Start multiprocessing
@@ -272,7 +281,7 @@ class (or object)
         solar_zenith, solar_azimuth, surface_tilt, surface_azimuth,\
         albedo, cls_pvarray, cls_engine, cls_irradiance, cls_vf, \
         fast_mode_pvrow_index, fast_mode_segment_index, \
-        irradiance_model_params, ghi, idx = args
+        irradiance_model_params, vf_calculator_params, ghi, idx = args
 
     report = run_timeseries_engine(
         report_builder.build, pvarray_parameters,
@@ -283,6 +292,7 @@ class (or object)
         fast_mode_pvrow_index=fast_mode_pvrow_index,
         fast_mode_segment_index=fast_mode_segment_index,
         irradiance_model_params=irradiance_model_params,
+        vf_calculator_params=vf_calculator_params,
         ghi=ghi)
 
     return report, idx

pvfactors/tests/conftest.py

@@ -138,6 +138,7 @@ def fn_report_example():
     def fn_report(pvarray): return {
         'qinc_front': pvarray.ts_pvrows[1].front.get_param_weighted('qinc'),
         'qinc_back': pvarray.ts_pvrows[1].back.get_param_weighted('qinc'),
+        'qabs_back': pvarray.ts_pvrows[1].back.get_param_weighted('qabs'),
         'iso_front': pvarray.ts_pvrows[1]
         .front.get_param_weighted('isotropic'),
         'iso_back': pvarray.ts_pvrows[1].back.get_param_weighted('isotropic')}

pvfactors/tests/test_engine.py

@@ -41,6 +41,10 @@ def test_pvengine_float_inputs_iso(params):
         pvarray.ts_pvrows[1].front.get_param_weighted('qinc'), 1099.6948573)
     np.testing.assert_almost_equal(
         pvarray.ts_pvrows[2].front.get_param_weighted('qinc'), 1102.76149246)
+    # Check absorbed
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[1].front.get_param_weighted('qabs'),
+        1099.6948573 * 0.99)
 
 
 def test_pvengine_float_inputs_perez(params):
@@ -80,6 +84,13 @@ def test_pvengine_float_inputs_perez(params):
     np.testing.assert_almost_equal(
         pvarray.ts_pvrows[1].back.get_param_weighted('qinc'),
         116.49050349491208)
+    # Check absorbed irradiance
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[2].front.get_param_weighted('qabs'),
+        1112.37717553 * 0.99)
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[1].back.get_param_weighted('qabs'),
+        116.49050349491208 * 0.97)
 
 
 def test_pvengine_ts_inputs_perez(params_serial,
@@ -115,6 +126,8 @@ def test_pvengine_ts_inputs_perez(params_serial,
         report['iso_front'], [42.816637, 42.780206])
     np.testing.assert_array_almost_equal(
         report['iso_back'], [1.727308, 1.726535])
+    np.testing.assert_array_almost_equal(
+        report['qabs_back'], report['qinc_back'] * 0.97)
 
 
 def test_run_fast_mode_isotropic(params):
@@ -545,7 +558,7 @@ def test_check_direct_shading_continuity():
     np.testing.assert_allclose(out, expected_out)
 
 
-def test_create_with_rho_init(params, pvmodule_canadian):
+def test_create_engine_with_rho_init(params, pvmodule_canadian):
     """Check that can create PV engine with rho initialization
     from faoi functions"""
     # Create inputs
@@ -559,3 +572,51 @@ def test_create_with_rho_init(params, pvmodule_canadian):
     # Check that rho values are the ones calculated
     np.testing.assert_allclose(engine.irradiance.rho_front, 0.02900688)
     np.testing.assert_allclose(engine.irradiance.rho_back, 0.02900688)
+
+
+def test_create_engine_with_rho_init(params, pvmodule_canadian):
+    """Run PV engine calcs with faoi functions for AOI losses"""
+
+    # Irradiance inputs
+    timestamps = dt.datetime(2019, 6, 11, 11)
+    DNI = 1000.
+    DHI = 100.
+
+    irradiance_model = HybridPerezOrdered()
+    pvarray = OrderedPVArray.init_from_dict(params)
+    faoi_fn = faoi_fn_from_pvlib_sandia(pvmodule_canadian)
+    vfcalculator = VFCalculator(faoi_fn_front=faoi_fn, faoi_fn_back=faoi_fn)
+    eng = PVEngine(pvarray, irradiance_model=irradiance_model,
+                   vf_calculator=vfcalculator)
+
+    # Make sure aoi methods are available
+    assert eng.vf_calculator.vf_aoi_methods is not None
+
+    # Fit engine
+    eng.fit(timestamps, DNI, DHI,
+            params['solar_zenith'],
+            params['solar_azimuth'],
+            params['surface_tilt'],
+            params['surface_azimuth'],
+            params['rho_ground'])
+
+    # Run timestep
+    pvarray = eng.run_full_mode(fn_build_report=lambda pvarray: pvarray)
+    # Checks
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[0].front.get_param_weighted('qinc'),
+        1110.1164773159298)
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[1].front.get_param_weighted('qinc'), 1110.595903991)
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[2].front.get_param_weighted('qinc'), 1112.37717553)
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[1].back.get_param_weighted('qinc'),
+        116.49050349491208)
+    # Check absorbed irradiance: calculated using faoi functions
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[2].front.get_param_weighted('qabs'),
+        [1099.1251094])
+    np.testing.assert_almost_equal(
+        pvarray.ts_pvrows[1].back.get_param_weighted('qabs'),
+        [114.4690984])