From 38922c8cd3b0a1bd97a5698ee65e99b35a35495d Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 25 Feb 2022 08:33:36 -0500 Subject: [PATCH 1/8] proof of concept sparse linear solve --- pvfactors/engine.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index 9d8b252..b3e059f 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -2,11 +2,24 @@ timeseries simulations.""" import numpy as np +from scipy.sparse import csc_matrix +from scipy.sparse.linalg import spsolve from pvfactors.viewfactors import VFCalculator from pvfactors.irradiance import HybridPerezOrdered from pvfactors.config import DEFAULT_RHO_FRONT, DEFAULT_RHO_BACK +def _sparse_solve(A, b): + xs = [] + for A_slice, b_slice in zip(A, b.T): + A_sparse = csc_matrix(A_slice) + b_sparse = csc_matrix(b_slice).T + x = spsolve(A_sparse, b_sparse) + xs.append(x) + x = np.stack(xs).T + return x + + class PVEngine(object): """Class putting all of the calculations together into simple workflows. @@ -218,12 +231,8 @@ def run_full_mode(self, fn_build_report=None): # Subtract matrices: will rely on broadcasting # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1 a_mat = invrho_ts_diag - ts_vf_matrix_reshaped - # Calculate inverse, requires specific shape - # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1 - inv_a_mat = np.linalg.inv(a_mat) - # Use einstein sum to get final timeseries radiosities - # shape = n_surfaces + 1, n_timesteps - q0 = np.einsum('ijk,ki->ji', inv_a_mat, irradiance_mat) + # solve the linear system a_mat * q0 = irradiance_mat for q0 + q0 = _sparse_solve(a_mat, irradiance_mat) # Calculate incident irradiance: will rely on broadcasting # shape = n_surfaces + 1, n_timesteps qinc = np.einsum('ijk,ki->ji', invrho_ts_diag, q0) From 4450fcc7f97d9af424c7b812456c395a00941bf6 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Fri, 25 Feb 2022 10:38:06 -0500 Subject: [PATCH 2/8] a nice sprinkling of dels --- pvfactors/engine.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index b3e059f..c5d6ed1 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -228,29 +228,36 @@ def run_full_mode(self, fn_build_report=None): invrho_ts_diag = np.zeros(shape_system) for idx_surf in range(shape_system[1]): invrho_ts_diag[:, idx_surf, idx_surf] = invrho_mat[idx_surf, :] + del invrho_mat # Subtract matrices: will rely on broadcasting # shape = n_timesteps, n_surfaces + 1, n_surfaces + 1 a_mat = invrho_ts_diag - ts_vf_matrix_reshaped + del ts_vf_matrix_reshaped # solve the linear system a_mat * q0 = irradiance_mat for q0 q0 = _sparse_solve(a_mat, irradiance_mat) + del a_mat # Calculate incident irradiance: will rely on broadcasting # shape = n_surfaces + 1, n_timesteps qinc = np.einsum('ijk,ki->ji', invrho_ts_diag, q0) + del invrho_ts_diag # --- Derive other irradiance terms # shape = n_surfaces, n_timesteps isotropic_mat = ts_vf_matrix[:-1, -1, :] * irradiance_mat[-1, :] + del ts_vf_matrix reflection_mat = qinc[:-1, :] - irradiance_mat[:-1, :] - isotropic_mat - + del irradiance_mat # --- Calculate AOI losses and absorbed irradiance # Create tiled reflection matrix of # shape = n_surfaces + 1, n_surfaces + 1, n_timestamps rho_ts_tiled = np.moveaxis(np.tile(rho_mat.T, (shape_system[1], 1, 1)), -1, 0) + del rho_mat # Get vf AOI matrix # shape [n_surfaces + 1, n_surfaces + 1, n_timestamps] vf_aoi_matrix = (self.vf_calculator .build_ts_vf_aoi_matrix(pvarray, rho_ts_tiled)) + del rho_ts_tiled pvarray.ts_vf_aoi_matrix = vf_aoi_matrix # Get absorbed irradiance matrix # shape [n_surfaces, n_timestamps] @@ -259,6 +266,8 @@ def run_full_mode(self, fn_build_report=None): # Calculate absorbed irradiance qabs = (np.einsum('ijk,jk->ik', vf_aoi_matrix, q0)[:-1, :] + irradiance_abs_mat) + del irradiance_abs_mat + del vf_aoi_matrix # --- Update surfaces with values: the lists are ordered by index for idx_surf, ts_surface in enumerate(pvarray.all_ts_surfaces): From 077a160440094809673395257c99f3ae931db079 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Sat, 26 Feb 2022 08:49:31 -0500 Subject: [PATCH 3/8] document sparse solve helper function --- pvfactors/engine.py | 35 +++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index c5d6ed1..27b1db3 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -9,14 +9,41 @@ from pvfactors.config import DEFAULT_RHO_FRONT, DEFAULT_RHO_BACK -def _sparse_solve(A, b): +def _sparse_solve_3D(A, b): + """ + Solve the linear system A*x=b for x, where A is sparse (contains mostly + zeros). + + For large matrices with many zeros, this is faster and uses less memory + than the dense analog `np.linalg.solve`. + + Parameters + ---------- + A : np.ndarray of shape (M, N, N) + First dimension is timestamps, second and third are surface dimensions + b : np.ndarray of shape (M, N) + First dimension is timestamps, second is surfaces + + Returns + ------- + x : np.ndarray of shape (M, N) + First dimension is timestamps, second is surfaces + """ + # implementation notes: + # - csc_matrix seemed to be the fastest option of the various + # sparse matrix formats in scipy.sparse + # - unfortunately the sparse matrix formats are 2-D only, so + # iteration across the time dimension is required + # - scipy 1.8.0 added new sparse arrays (as opposed to sparse matrices); + # they are 2-D only at time of writing, but in the future if they + # become N-D then it may make sense to use them here. xs = [] - for A_slice, b_slice in zip(A, b.T): + for A_slice, b_slice in zip(A, b): A_sparse = csc_matrix(A_slice) b_sparse = csc_matrix(b_slice).T x = spsolve(A_sparse, b_sparse) xs.append(x) - x = np.stack(xs).T + x = np.stack(xs) return x @@ -234,7 +261,7 @@ def run_full_mode(self, fn_build_report=None): a_mat = invrho_ts_diag - ts_vf_matrix_reshaped del ts_vf_matrix_reshaped # solve the linear system a_mat * q0 = irradiance_mat for q0 - q0 = _sparse_solve(a_mat, irradiance_mat) + q0 = _sparse_solve_3D(a_mat, irradiance_mat.T).T del a_mat # Calculate incident irradiance: will rely on broadcasting # shape = n_surfaces + 1, n_timesteps From e0e1e3e37d7bedba78ea201d7f790fe87adf9089 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Sat, 26 Feb 2022 08:49:41 -0500 Subject: [PATCH 4/8] add scipy dependency --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index 38a4094..1af9f4f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ pvlib>=0.9.0,<0.10.0 shapely>=1.6.4.post2,<2 +scipy>=1.2.0 matplotlib future six From 5268d5fdb7659791cba34ad497e76c95f358533e Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Sat, 26 Feb 2022 09:09:54 -0500 Subject: [PATCH 5/8] add basic sparse solve test --- pvfactors/tests/test_engine.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/pvfactors/tests/test_engine.py b/pvfactors/tests/test_engine.py index c5f1991..2ff8e8c 100644 --- a/pvfactors/tests/test_engine.py +++ b/pvfactors/tests/test_engine.py @@ -1,4 +1,4 @@ -from pvfactors.engine import PVEngine +from pvfactors.engine import PVEngine, _sparse_solve_3D from pvfactors.geometry.pvarray import OrderedPVArray from pvfactors.irradiance import IsotropicOrdered, HybridPerezOrdered from pvfactors.irradiance.utils import breakup_df_inputs @@ -709,3 +709,13 @@ def test_engine_variable_albedo(params, df_inputs_clearsky_8760): expected_bfg_after_aoi = 14.670709 np.testing.assert_allclose(bfg, expected_bfg) np.testing.assert_allclose(bfg_after_aoi, expected_bfg_after_aoi) + + +def test__sparse_solve_3d(): + """Verify the sparse solver returns same results as np.linalg.solve""" + A = np.random.rand(5, 3, 3) + b = np.random.rand(5, 3) + x_dense = np.linalg.solve(A, b) + x_sparse = _sparse_solve_3D(A, b) + assert x_sparse.shape == b.shape + np.testing.assert_allclose(x_dense, x_sparse) From 12f93e5895a475b358394a2d8da3cd74f4804eaa Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Sat, 26 Feb 2022 09:12:31 -0500 Subject: [PATCH 6/8] what's new --- docs/sphinx/whatsnew/v1.5.3.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/sphinx/whatsnew/v1.5.3.rst b/docs/sphinx/whatsnew/v1.5.3.rst index e8bbad4..4fb0b83 100644 --- a/docs/sphinx/whatsnew/v1.5.3.rst +++ b/docs/sphinx/whatsnew/v1.5.3.rst @@ -6,12 +6,17 @@ v1.5.3 (TBD) Enhancements ------------ -* Add python3.10 to test configuration (ghpull:`129`) +* Add python 3.10 to test configuration (:ghpull:`129`) +* :py:meth:`pvfactors.engine.PVEngine.run_full_mode` is faster and uses less + memory for simulations with large ``n_pvrows`` and many timestamps (:ghpull:`140`) Requirements ------------ +* ``scipy>=1.2.0`` is now a direct requirement (before now it was only + indirectly required through pvlib) (:ghpull:`140`) + Bug Fixes ---------- From 0d2e6f6d78bc59e0082b3e988f40ffdcb0d4841b Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Sun, 27 Feb 2022 10:54:34 -0500 Subject: [PATCH 7/8] edits from review --- pvfactors/engine.py | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index 27b1db3..b823198 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -9,7 +9,7 @@ from pvfactors.config import DEFAULT_RHO_FRONT, DEFAULT_RHO_BACK -def _sparse_solve_3D(A, b): +def _sparse_solve_3D(A: np.ndarray, b: np.ndarray) -> np.ndarray: """ Solve the linear system A*x=b for x, where A is sparse (contains mostly zeros). @@ -28,15 +28,17 @@ def _sparse_solve_3D(A, b): ------- x : np.ndarray of shape (M, N) First dimension is timestamps, second is surfaces + + Notes + ----- + - csc_matrix seemed to be the fastest option of the various + sparse matrix formats in scipy.sparse + - unfortunately the sparse matrix formats are 2-D only, so + iteration across the time dimension is required + - scipy 1.8.0 added new sparse arrays (as opposed to sparse matrices); + they are 2-D only at time of writing, but in the future if they + become N-D then it may make sense to use them here. """ - # implementation notes: - # - csc_matrix seemed to be the fastest option of the various - # sparse matrix formats in scipy.sparse - # - unfortunately the sparse matrix formats are 2-D only, so - # iteration across the time dimension is required - # - scipy 1.8.0 added new sparse arrays (as opposed to sparse matrices); - # they are 2-D only at time of writing, but in the future if they - # become N-D then it may make sense to use them here. xs = [] for A_slice, b_slice in zip(A, b): A_sparse = csc_matrix(A_slice) @@ -230,6 +232,14 @@ def run_full_mode(self, fn_build_report=None): report Saved results from the simulation, as specified by user's report function. If no function is passed, nothing will be returned. + + Notes + ----- + This function allocates several large arrays, some of which are only + needed for part of the function. To reduce peak memory usage, this + function uses the `del` statement to allow intermediate arrays to be + garbage collected early and the underlying memory to be reused sooner + than if the arrays were allowed to go out of scope naturally. """ # Get pvarray pvarray = self.pvarray From 480b350f5f5171297f79c7860f3a6cfb40f90585 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Thu, 3 Mar 2022 15:00:06 -0500 Subject: [PATCH 8/8] bugfix for empty inputs --- pvfactors/engine.py | 4 ++++ pvfactors/tests/test_engine.py | 23 +++++++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index b823198..e398a19 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -39,6 +39,10 @@ def _sparse_solve_3D(A: np.ndarray, b: np.ndarray) -> np.ndarray: they are 2-D only at time of writing, but in the future if they become N-D then it may make sense to use them here. """ + if b.size == 0: + # prevent ValueError from np.stack call below. + # it's empty, so values don't matter -- just shape (and dtype?) + return np.zeros_like(b) xs = [] for A_slice, b_slice in zip(A, b): A_sparse = csc_matrix(A_slice) diff --git a/pvfactors/tests/test_engine.py b/pvfactors/tests/test_engine.py index 2ff8e8c..624320e 100644 --- a/pvfactors/tests/test_engine.py +++ b/pvfactors/tests/test_engine.py @@ -719,3 +719,26 @@ def test__sparse_solve_3d(): x_sparse = _sparse_solve_3D(A, b) assert x_sparse.shape == b.shape np.testing.assert_allclose(x_dense, x_sparse) + + +def test_engine_empty_inputs(params, df_inputs_clearsky_8760): + """Empty inputs yields empty outputs""" + df_inputs = df_inputs_clearsky_8760.iloc[:0] + timestamps = df_inputs.index + dni = df_inputs.dni.values + dhi = df_inputs.dhi.values + solar_zenith = df_inputs.solar_zenith.values + solar_azimuth = df_inputs.solar_azimuth.values + surface_tilt = df_inputs.surface_tilt.values + surface_azimuth = df_inputs.surface_azimuth.values + + # Run engine + pvarray = OrderedPVArray.init_from_dict(params) + eng = PVEngine(pvarray) + eng.fit(timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, + surface_azimuth, params['rho_ground']) + qinc = eng.run_full_mode( + fn_build_report=lambda pvarray: (pvarray.ts_pvrows[1] + .back.get_param_weighted('qinc'))) + + assert len(qinc) == 0