From 9c5bc1c2524b48206312f7bd1cc400f589211967 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 18:38:35 -0500
Subject: [PATCH 01/17] vectorize _vf_ground_sky_2d across surface_tilt

also a few other micro-optimizations
---
 pvlib/bifacial/infinite_sheds.py | 22 +++++------------
 pvlib/bifacial/utils.py          | 42 ++++++++++++++++++--------------
 2 files changed, 30 insertions(+), 34 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index 91b33551bc..711251130b 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -38,17 +38,10 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
 
     Returns
     -------
-    fgnd_sky : float
+    fgnd_sky : numeric
         Integration of view factor over the length between adjacent, interior
-        rows. [unitless]
-    fz : ndarray
-        Fraction of distance from the previous row to the next row. [unitless]
-    fz_sky : ndarray
-        View factors at discrete points between adjacent, interior rows.
-        [unitless]
-
+        rows.  Shape matches that of ``surface_tilt``. [unitless]
     """
-    # TODO: vectorize over surface_tilt
     # Abuse utils._vf_ground_sky_2d by supplying surface_tilt in place
     # of a signed rotation. This is OK because
     # 1) z span the full distance between 2 rows, and
@@ -57,12 +50,9 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    fz_sky = np.zeros((len(rotation), npoints))
-    for k, r in enumerate(rotation):
-        vf, _ = utils._vf_ground_sky_2d(z, r, gcr, pitch, height, max_rows)
-        fz_sky[k, :] = vf
+    fz_sky, _ = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
     # calculate the integrated view factor for all of the ground between rows
-    return np.trapz(fz_sky, z, axis=1)
+    return np.trapz(fz_sky, z, axis=0)
 
 
 def _poa_ground_shadows(poa_ground, f_gnd_beam, df, vf_gnd_sky):
@@ -469,7 +459,7 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
         on the surface that is not reflected away. [unitless]
 
     npoints : int, default 100
-        Number of points used to discretize distance along the ground.
+        Number of discretization points for calculating integrated viewfactors.
 
     Returns
     -------
@@ -701,7 +691,7 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         etc. A negative value is a reduction in back irradiance. [unitless]
 
     npoints : int, default 100
-        Number of points used to discretize distance along the ground.
+        Number of discretization points for calculating integrated viewfactors.
 
     Returns
     -------
diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index bac8923536..8e83e66535 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -4,7 +4,7 @@
 """
 import numpy as np
 from pvlib.tools import sind, cosd, tand
-
+import time
 
 def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth):
     """
@@ -104,7 +104,7 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
         Position on the ground between two rows, as a fraction of the pitch.
         x = 0 corresponds to the point on the ground directly below the
         center point of a row. Positive x is towards the right. [unitless]
-    rotation : float
+    rotation : numeric
         Rotation angle of the row's right edge relative to row center.
         [degree]
     gcr : float
@@ -124,26 +124,32 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
         Fraction of sky dome visible from each point on the ground. [unitless]
     wedge_angles : array
         Angles defining each wedge of sky that is blocked by a row. Shape is
-        (2, len(x), 2*max_rows+1). ``wedge_angles[0,:,:]`` is the
-        starting angle of each wedge, ``wedge_angles[1,:,:]`` is the end angle.
-        [degree]
+        (2, len(x), len(rotation), 2*max_rows+1). ``wedge_angles[0,:,:,:]``
+        is the starting angle of each wedge, ``wedge_angles[1,:,:,:]`` is the
+        end angle. [degree]
     """
-    x = np.atleast_1d(x)  # handle float
+    # handle floats:
+    st = time.perf_counter()
+    x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
+    rotation = np.atleast_1d(rotation)[np.newaxis, :, np.newaxis]
     all_k = np.arange(-max_rows, max_rows + 1)
     width = gcr * pitch / 2.
+    distance_to_row_centers = (all_k - x) * pitch
+    dy = width * sind(rotation)
+    dx = width * cosd(rotation)
     # angles from x to right edge of each row
-    a1 = height + width * sind(rotation)
-    b1 = (all_k - x[:, np.newaxis]) * pitch + width * cosd(rotation)
-    phi_1 = np.degrees(np.arctan2(a1, b1))
+    a1 = height + dy
+    b1 = distance_to_row_centers + dx
+    phi_1 = np.arctan2(a1, b1)  # dimensions: (x, rotation, row)
     # angles from x to left edge of each row
-    a2 = height - width * sind(rotation)
-    b2 = (all_k - x[:, np.newaxis]) * pitch - width * cosd(rotation)
-    phi_2 = np.degrees(np.arctan2(a2, b2))
-    phi = np.stack([phi_1, phi_2])
-    swap = phi[0, :, :] > phi[1, :, :]
-    # swap where phi_1 > phi_2 so that phi_1[0,:,:] is the lesser angle
+    a2 = height - dy
+    b2 = distance_to_row_centers - dx
+    phi_2 = np.arctan2(a2, b2)  # dimensions: (x, rotation, row)
+    phi = np.stack([phi_1, phi_2])  # dimensions: (row edge, x, rotation, row)
+    swap = phi_1 > phi_2
+    # swap where phi_1 > phi_2 so that phi[0,:,:,:] is the lesser angle
     phi = np.where(swap, phi[::-1], phi)
     # right edge of next row - left edge of previous row
-    wedge_vfs = 0.5 * (cosd(phi[1, :, 1:]) - cosd(phi[0, :, :-1]))
-    vf = np.sum(np.where(wedge_vfs > 0, wedge_vfs, 0.), axis=1)
-    return vf, phi
+    wedge_vfs = 0.5 * (np.cos(phi[1, :, :, 1:]) - np.cos(phi[0, :, :, :-1]))
+    vf = np.sum(np.clip(wedge_vfs, a_min=0., a_max=None), axis=-1)
+    return vf, np.degrees(phi)

From af68f8478ed5045bd091c4d33c405ced87b30ae5 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 19:12:28 -0500
Subject: [PATCH 02/17] whatsnew

---
 docs/sphinx/source/whatsnew/v0.9.5.rst | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst
index bca0f7b8f5..d9c8ab7ea2 100644
--- a/docs/sphinx/source/whatsnew/v0.9.5.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.5.rst
@@ -20,12 +20,14 @@ Enhancements
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
 * Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
   support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
-
 * Added an optional ``model`` parameter to
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa`
   to enable use of the hay-davies sky diffuse irradiance model
   instead of the default isotropic model. (:pull:`1668`)
+* Improved execution speed of the infinite sheds model by 10-25% for
+  simulations with time series surface orientation. (:issue:`1680`, :pull:`1682`)
+
 
 Bug fixes
 ~~~~~~~~~
@@ -74,3 +76,4 @@ Contributors
 * Michael Deceglie (:ghuser:`mdeceglie`)
 * Saurabh Aneja (:ghuser:`spaneja`)
 * John Moseley (:ghuser:`johnMoseleyArray`)
+* Areeba Turabi (:ghuser:`aturabi`)

From e5bdf5a58724c7426cdc36149438603647a98b6a Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 19:24:19 -0500
Subject: [PATCH 03/17] remove straggler instrumentation

---
 pvlib/bifacial/utils.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index 8e83e66535..ca54083f1c 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -4,7 +4,6 @@
 """
 import numpy as np
 from pvlib.tools import sind, cosd, tand
-import time
 
 def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth):
     """
@@ -129,7 +128,6 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
         end angle. [degree]
     """
     # handle floats:
-    st = time.perf_counter()
     x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
     rotation = np.atleast_1d(rotation)[np.newaxis, :, np.newaxis]
     all_k = np.arange(-max_rows, max_rows + 1)

From fb1d943f3b019ec3db9554f9e6b81aace0c54bd7 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 19:27:52 -0500
Subject: [PATCH 04/17] docstring improvement

---
 pvlib/bifacial/utils.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index ca54083f1c..a1bb9c9e2c 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -119,8 +119,9 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
 
     Returns
     -------
-    vf : numeric
-        Fraction of sky dome visible from each point on the ground. [unitless]
+    vf : array
+        Fraction of sky dome visible from each point on the ground.
+        Shape is (len(x), len(rotation)). [unitless]
     wedge_angles : array
         Angles defining each wedge of sky that is blocked by a row. Shape is
         (2, len(x), len(rotation), 2*max_rows+1). ``wedge_angles[0,:,:,:]``

From bb82373c958a69d409023cc5af7eff62c66ece14 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 19:33:14 -0500
Subject: [PATCH 05/17] stickler

---
 pvlib/bifacial/infinite_sheds.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index 711251130b..634154bccd 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -50,7 +50,8 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    fz_sky, _ = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
+    fz_sky, _ = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height,
+                                        max_rows)
     # calculate the integrated view factor for all of the ground between rows
     return np.trapz(fz_sky, z, axis=0)
 

From ee421a49f760fbbb84cb7242f492bdf38ed7b9ca Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 2 Mar 2023 21:42:34 -0500
Subject: [PATCH 06/17] update failing test

---
 pvlib/tests/bifacial/test_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/tests/bifacial/test_utils.py b/pvlib/tests/bifacial/test_utils.py
index 388a222394..7b6240769a 100644
--- a/pvlib/tests/bifacial/test_utils.py
+++ b/pvlib/tests/bifacial/test_utils.py
@@ -35,7 +35,7 @@ def test_system_fixed_tilt():
     c22 = (-2 - sqr3) / np.sqrt(1.25**2 + (2 + sqr3)**2)  # right edge row 0
     c23 = (0 - sqr3) / np.sqrt(1.25**2 + (0 - sqr3)**2)  # right edge row 1
     vf_2 = 0.5 * (c23 - c22 + c21 - c20)  # vf at point 1
-    vfs_ground_sky = np.array([vf_0, vf_1, vf_2])
+    vfs_ground_sky = np.array([[vf_0], [vf_1], [vf_2]])
     return syst, pts, vfs_ground_sky
 
 

From 0d80d7d12d2f9c4178577ab8c90b70ec5c008adf Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Fri, 3 Mar 2023 12:58:51 -0500
Subject: [PATCH 07/17] remove "per-point" language

---
 pvlib/bifacial/infinite_sheds.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index 634154bccd..f09d2afb91 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -12,8 +12,8 @@
 def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
                          pitch, max_rows=10, npoints=100):
     """
-    Integrated and per-point view factors from the ground to the sky at points
-    between interior rows of the array.
+    Integrated view factor from the ground to the sky underneath
+    interior rows of the array.
 
     Parameters
     ----------

From 9bf93601d7e1ff954313e2a0239d7f23a91d2f8b Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Wed, 8 Mar 2023 12:59:38 -0500
Subject: [PATCH 08/17] sprinkle in some `del`s

---
 pvlib/bifacial/utils.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index a1bb9c9e2c..391e35194f 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -128,6 +128,12 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
         is the starting angle of each wedge, ``wedge_angles[1,:,:,:]`` is the
         end angle. [degree]
     """
+    # This function creates large float64 arrays of size
+    # (2*len(x)*len(rotation)*len(max_rows)) or ~100 MB for
+    # typical time series inputs.  This function uses `del` to
+    # allow python to recapture intermediate variables and
+    # reduce peak memory usage.
+
     # handle floats:
     x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
     rotation = np.atleast_1d(rotation)[np.newaxis, :, np.newaxis]
@@ -136,19 +142,30 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
     distance_to_row_centers = (all_k - x) * pitch
     dy = width * sind(rotation)
     dx = width * cosd(rotation)
+
     # angles from x to right edge of each row
     a1 = height + dy
     b1 = distance_to_row_centers + dx
     phi_1 = np.arctan2(a1, b1)  # dimensions: (x, rotation, row)
+    del a1, b1  # reduce max memory usage
+
     # angles from x to left edge of each row
     a2 = height - dy
     b2 = distance_to_row_centers - dx
     phi_2 = np.arctan2(a2, b2)  # dimensions: (x, rotation, row)
+    del a2, b2  # reduce max memory usage
+
     phi = np.stack([phi_1, phi_2])  # dimensions: (row edge, x, rotation, row)
     swap = phi_1 > phi_2
+    del phi_1, phi_2  # reduce max memory usage
+
     # swap where phi_1 > phi_2 so that phi[0,:,:,:] is the lesser angle
     phi = np.where(swap, phi[::-1], phi)
+    del swap  # reduce max memory usage
+
     # right edge of next row - left edge of previous row
     wedge_vfs = 0.5 * (np.cos(phi[1, :, :, 1:]) - np.cos(phi[0, :, :, :-1]))
     vf = np.sum(np.clip(wedge_vfs, a_min=0., a_max=None), axis=-1)
+    del wedge_vfs  # reduce max memory usage
+
     return vf, np.degrees(phi)

From 0bc108719b99e95e722ae326eacdaf92b7362690 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Wed, 8 Mar 2023 17:19:31 -0500
Subject: [PATCH 09/17] contorted version using numpy's `out` parameter

---
 pvlib/bifacial/infinite_sheds.py   |  3 +-
 pvlib/bifacial/utils.py            | 45 +++++++++++++-----------------
 pvlib/tests/bifacial/test_utils.py |  8 +++---
 3 files changed, 24 insertions(+), 32 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index f09d2afb91..f6634a3602 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -50,8 +50,7 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    fz_sky, _ = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height,
-                                        max_rows)
+    fz_sky = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
     # calculate the integrated view factor for all of the ground between rows
     return np.trapz(fz_sky, z, axis=0)
 
diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index 391e35194f..3831a76fef 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -122,17 +122,12 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
     vf : array
         Fraction of sky dome visible from each point on the ground.
         Shape is (len(x), len(rotation)). [unitless]
-    wedge_angles : array
-        Angles defining each wedge of sky that is blocked by a row. Shape is
-        (2, len(x), len(rotation), 2*max_rows+1). ``wedge_angles[0,:,:,:]``
-        is the starting angle of each wedge, ``wedge_angles[1,:,:,:]`` is the
-        end angle. [degree]
     """
     # This function creates large float64 arrays of size
     # (2*len(x)*len(rotation)*len(max_rows)) or ~100 MB for
-    # typical time series inputs.  This function uses `del` to
-    # allow python to recapture intermediate variables and
-    # reduce peak memory usage.
+    # typical time series inputs.  This function makes heavy
+    # use of numpy's out parameter to avoid allocating new
+    # memory.  Unfortunately that comes at the cost of readability.
 
     # handle floats:
     x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
@@ -143,29 +138,27 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
     dy = width * sind(rotation)
     dx = width * cosd(rotation)
 
+    phi = np.empty((2, x.shape[0], rotation.shape[1], len(all_k)))
+
     # angles from x to right edge of each row
     a1 = height + dy
-    b1 = distance_to_row_centers + dx
-    phi_1 = np.arctan2(a1, b1)  # dimensions: (x, rotation, row)
-    del a1, b1  # reduce max memory usage
+    phi[0] = distance_to_row_centers + dx
+    np.arctan2(a1, phi[0], out=phi[0])
 
     # angles from x to left edge of each row
     a2 = height - dy
-    b2 = distance_to_row_centers - dx
-    phi_2 = np.arctan2(a2, b2)  # dimensions: (x, rotation, row)
-    del a2, b2  # reduce max memory usage
-
-    phi = np.stack([phi_1, phi_2])  # dimensions: (row edge, x, rotation, row)
-    swap = phi_1 > phi_2
-    del phi_1, phi_2  # reduce max memory usage
+    phi[1] = distance_to_row_centers - dx
+    np.arctan2(a2, phi[1], out=phi[1])
 
-    # swap where phi_1 > phi_2 so that phi[0,:,:,:] is the lesser angle
-    phi = np.where(swap, phi[::-1], phi)
-    del swap  # reduce max memory usage
+    # swap angles so that phi[0,:,:,:] is the lesser angle
+    phi.sort(axis=0)
 
     # right edge of next row - left edge of previous row
-    wedge_vfs = 0.5 * (np.cos(phi[1, :, :, 1:]) - np.cos(phi[0, :, :, :-1]))
-    vf = np.sum(np.clip(wedge_vfs, a_min=0., a_max=None), axis=-1)
-    del wedge_vfs  # reduce max memory usage
-
-    return vf, np.degrees(phi)
+    next_edge = phi[1, :, :, 1:]
+    np.cos(next_edge, out=next_edge)
+    prev_edge = phi[0, :, :, :-1]
+    np.cos(prev_edge, out=prev_edge)
+    np.subtract(next_edge, prev_edge, out=next_edge)
+    np.clip(next_edge, a_min=0., a_max=None, out=next_edge)
+    vf = np.sum(next_edge, axis=-1) / 2
+    return vf
diff --git a/pvlib/tests/bifacial/test_utils.py b/pvlib/tests/bifacial/test_utils.py
index 7b6240769a..f9851d5083 100644
--- a/pvlib/tests/bifacial/test_utils.py
+++ b/pvlib/tests/bifacial/test_utils.py
@@ -79,10 +79,10 @@ def test__unshaded_ground_fraction(
 def test__vf_ground_sky_2d(test_system_fixed_tilt):
     # vector input
     ts, pts, vfs_gnd_sky = test_system_fixed_tilt
-    vfs, _ = utils._vf_ground_sky_2d(pts, ts['rotation'], ts['gcr'],
-                                     ts['pitch'], ts['height'], max_rows=1)
+    vfs = utils._vf_ground_sky_2d(pts, ts['rotation'], ts['gcr'],
+                                  ts['pitch'], ts['height'], max_rows=1)
     assert np.allclose(vfs, vfs_gnd_sky, rtol=0.1)  # middle point vf is off
     # test with singleton x
-    vf, _ = utils._vf_ground_sky_2d(pts[0], ts['rotation'], ts['gcr'],
-                                    ts['pitch'], ts['height'], max_rows=1)
+    vf = utils._vf_ground_sky_2d(pts[0], ts['rotation'], ts['gcr'],
+                                 ts['pitch'], ts['height'], max_rows=1)
     assert np.isclose(vf, vfs_gnd_sky[0])

From 6f490ea363dad00f6ea72c228b7bb1ff4dcade92 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Wed, 8 Mar 2023 17:32:12 -0500
Subject: [PATCH 10/17] even more out

---
 pvlib/bifacial/utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index 3831a76fef..becaa89697 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -142,12 +142,12 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
 
     # angles from x to right edge of each row
     a1 = height + dy
-    phi[0] = distance_to_row_centers + dx
+    np.add(distance_to_row_centers, dx, out=phi[0])
     np.arctan2(a1, phi[0], out=phi[0])
 
     # angles from x to left edge of each row
     a2 = height - dy
-    phi[1] = distance_to_row_centers - dx
+    np.subtract(distance_to_row_centers, dx, out=phi[1])
     np.arctan2(a2, phi[1], out=phi[1])
 
     # swap angles so that phi[0,:,:,:] is the lesser angle

From 79b48491de27c0b683c6ccdbabca64646d77130f Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 9 Mar 2023 09:04:05 -0500
Subject: [PATCH 11/17] add vectorize kwarg to infinite_sheds

---
 docs/sphinx/source/whatsnew/v0.9.5.rst |  8 ++++---
 pvlib/bifacial/infinite_sheds.py       | 33 ++++++++++++++++++++------
 2 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst
index d9c8ab7ea2..0e2c6d8b03 100644
--- a/docs/sphinx/source/whatsnew/v0.9.5.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.5.rst
@@ -25,9 +25,11 @@ Enhancements
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa`
   to enable use of the hay-davies sky diffuse irradiance model
   instead of the default isotropic model. (:pull:`1668`)
-* Improved execution speed of the infinite sheds model by 10-25% for
-  simulations with time series surface orientation. (:issue:`1680`, :pull:`1682`)
-
+* Added an optional ``vectorize`` parameter to 
+  :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
+  :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` which,
+  when set to ``True``, enables faster calculation but with increased
+  memory usage. (:issue:`1680`, :pull:`1682`)
 
 Bug fixes
 ~~~~~~~~~
diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index f6634a3602..d85aedc73e 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -10,7 +10,7 @@
 from pvlib.irradiance import beam_component, aoi, haydavies
 
 def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
-                         pitch, max_rows=10, npoints=100):
+                         pitch, max_rows=10, npoints=100, vectorize=False):
     """
     Integrated view factor from the ground to the sky underneath
     interior rows of the array.
@@ -35,6 +35,9 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
         Maximum number of rows to consider in front and behind the current row.
     npoints : int, default 100
         Number of points used to discretize distance along the ground.
+    vectorize : bool, default False
+        If True, vectorize the view factor calculation across ``surface_tilt``.
+        This increases speed with the cost of increased memory usage.
 
     Returns
     -------
@@ -50,7 +53,14 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    fz_sky = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
+    # calculate the integrated view factor for all of the ground between rows
+    if vectorize:
+        fz_sky = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
+    else:
+        fz_sky = np.zeros((npoints, len(rotation)))
+        for k, r in enumerate(rotation):
+            vf = utils._vf_ground_sky_2d(z, r, gcr, pitch, height, max_rows)
+            fz_sky[:, k] = vf[:, 0]  # remove spurious rotation dimension
     # calculate the integrated view factor for all of the ground between rows
     return np.trapz(fz_sky, z, axis=0)
 
@@ -391,7 +401,7 @@ def _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt,
 def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
                        solar_azimuth, gcr, height, pitch, ghi, dhi, dni,
                        albedo, model='isotropic', dni_extra=None, iam=1.0,
-                       npoints=100):
+                       npoints=100, vectorize=False):
     r"""
     Calculate plane-of-array (POA) irradiance on one side of a row of modules.
 
@@ -461,6 +471,10 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
     npoints : int, default 100
         Number of discretization points for calculating integrated viewfactors.
 
+    vectorize : bool, default False
+        If True, vectorize the view factor calculation across ``surface_tilt``.
+        This increases speed with the cost of increased memory usage.
+
     Returns
     -------
     output : dict or DataFrame
@@ -527,7 +541,8 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
     # method differs from [1], Eq. 7 and Eq. 8; height is defined at row
     # center rather than at row lower edge as in [1].
     vf_gnd_sky = _vf_ground_sky_integ(
-        surface_tilt, surface_azimuth, gcr, height, pitch, max_rows, npoints)
+        surface_tilt, surface_azimuth, gcr, height, pitch, max_rows, npoints,
+        vectorize)
     # fraction of row slant height that is shaded from direct irradiance
     f_x = _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt,
                            surface_azimuth, gcr)
@@ -600,7 +615,7 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                    gcr, height, pitch, ghi, dhi, dni,
                    albedo, model='isotropic', dni_extra=None, iam_front=1.0,
                    iam_back=1.0, bifaciality=0.8, shade_factor=-0.02,
-                   transmission_factor=0, npoints=100):
+                   transmission_factor=0, npoints=100, vectorize=False):
     """
     Get front and rear irradiance using the infinite sheds model.
 
@@ -693,6 +708,10 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
     npoints : int, default 100
         Number of discretization points for calculating integrated viewfactors.
 
+    vectorize : bool, default False
+        If True, vectorize the view factor calculation across ``surface_tilt``.
+        This increases speed with the cost of increased memory usage.
+
     Returns
     -------
     output : dict or DataFrame
@@ -746,14 +765,14 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
         gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
         albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_front,
-        npoints=npoints)
+        npoints=npoints, vectorize=vectorize)
     # back side POA irradiance
     irrad_back = get_irradiance_poa(
         surface_tilt=backside_tilt, surface_azimuth=backside_sysaz,
         solar_zenith=solar_zenith, solar_azimuth=solar_azimuth,
         gcr=gcr, height=height, pitch=pitch, ghi=ghi, dhi=dhi, dni=dni,
         albedo=albedo, model=model, dni_extra=dni_extra, iam=iam_back,
-        npoints=npoints)
+        npoints=npoints, vectorize=vectorize)
 
     colmap_front = {
         'poa_global': 'poa_front',

From 3d3eef0aabf1c98f2cce9656aa42a10a52dddef4 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 9 Mar 2023 09:04:50 -0500
Subject: [PATCH 12/17] benchmark both vectorize=True and vectorize=False

---
 benchmarks/benchmarks/infinite_sheds.py | 26 ++++++++++++++++---------
 1 file changed, 17 insertions(+), 9 deletions(-)

diff --git a/benchmarks/benchmarks/infinite_sheds.py b/benchmarks/benchmarks/infinite_sheds.py
index 573cf19d7a..755c2ec3ef 100644
--- a/benchmarks/benchmarks/infinite_sheds.py
+++ b/benchmarks/benchmarks/infinite_sheds.py
@@ -10,7 +10,11 @@
 
 class InfiniteSheds:
 
-    def setup(self):
+    # benchmark variant parameters (run both vectorize=True and False)
+    params = [True, False]
+    param_names = ['vectorize']
+
+    def setup(self, vectorize):
         self.times = pd.date_range(start='20180601', freq='1min',
                                    periods=1440)
         self.location = location.Location(40, -80)
@@ -38,7 +42,7 @@ def setup(self):
                 gcr=self.gcr
             )
 
-    def time_get_irradiance_poa_fixed(self):
+    def time_get_irradiance_poa_fixed(self, vectorize):
         infinite_sheds.get_irradiance_poa(
             surface_tilt=self.surface_tilt,
             surface_azimuth=self.surface_azimuth,
@@ -51,10 +55,11 @@ def time_get_irradiance_poa_fixed(self):
             dhi=self.clearsky_irradiance['dhi'],
             dni=self.clearsky_irradiance['dni'],
             albedo=self.albedo,
-            npoints=self.npoints
+            npoints=self.npoints,
+            vectorize=vectorize,
         )
 
-    def time_get_irradiance_poa_tracking(self):
+    def time_get_irradiance_poa_tracking(self, vectorize):
         infinite_sheds.get_irradiance_poa(
             surface_tilt=self.tracking['surface_tilt'],
             surface_azimuth=self.tracking['surface_azimuth'],
@@ -67,10 +72,11 @@ def time_get_irradiance_poa_tracking(self):
             dhi=self.clearsky_irradiance['dhi'],
             dni=self.clearsky_irradiance['dni'],
             albedo=self.albedo,
-            npoints=self.npoints
+            npoints=self.npoints,
+            vectorize=vectorize,
         )
 
-    def time_get_irradiance_fixed(self):
+    def time_get_irradiance_fixed(self, vectorize):
         infinite_sheds.get_irradiance(
             surface_tilt=self.surface_tilt,
             surface_azimuth=self.surface_azimuth,
@@ -83,10 +89,11 @@ def time_get_irradiance_fixed(self):
             dhi=self.clearsky_irradiance['dhi'],
             dni=self.clearsky_irradiance['dni'],
             albedo=self.albedo,
-            npoints=self.npoints
+            npoints=self.npoints,
+            vectorize=vectorize,
         )
 
-    def time_get_irradiance_tracking(self):
+    def time_get_irradiance_tracking(self, vectorize):
         infinite_sheds.get_irradiance(
             surface_tilt=self.tracking['surface_tilt'],
             surface_azimuth=self.tracking['surface_azimuth'],
@@ -99,5 +106,6 @@ def time_get_irradiance_tracking(self):
             dhi=self.clearsky_irradiance['dhi'],
             dni=self.clearsky_irradiance['dni'],
             albedo=self.albedo,
-            npoints=self.npoints
+            npoints=self.npoints,
+            vectorize=vectorize,
         )

From 8445aa58c3a1e97fd987444025cbc5f63190bfc7 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 9 Mar 2023 09:10:27 -0500
Subject: [PATCH 13/17] improve comments

---
 pvlib/bifacial/utils.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index becaa89697..782ceb09f7 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -127,7 +127,9 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
     # (2*len(x)*len(rotation)*len(max_rows)) or ~100 MB for
     # typical time series inputs.  This function makes heavy
     # use of numpy's out parameter to avoid allocating new
-    # memory.  Unfortunately that comes at the cost of readability.
+    # memory.  Unfortunately that comes at the cost of some
+    # readability: because arrays get reused to avoid new allocations,
+    # variable names don't always match what they hold.
 
     # handle floats:
     x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
@@ -142,6 +144,7 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
 
     # angles from x to right edge of each row
     a1 = height + dy
+    # temporarily store one leg of the triangle in phi:
     np.add(distance_to_row_centers, dx, out=phi[0])
     np.arctan2(a1, phi[0], out=phi[0])
 
@@ -153,11 +156,15 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
     # swap angles so that phi[0,:,:,:] is the lesser angle
     phi.sort(axis=0)
 
-    # right edge of next row - left edge of previous row
+    # now re-use phi's memory again, this time storing cos(phi).
     next_edge = phi[1, :, :, 1:]
     np.cos(next_edge, out=next_edge)
     prev_edge = phi[0, :, :, :-1]
     np.cos(prev_edge, out=prev_edge)
+    # right edge of next row - left edge of previous row, again
+    # reusing memory so that the difference is stored in next_edge.
+    # Note that the 0.5 view factor coefficient is applied after summing
+    # as a minor speed optimization.
     np.subtract(next_edge, prev_edge, out=next_edge)
     np.clip(next_edge, a_min=0., a_max=None, out=next_edge)
     vf = np.sum(next_edge, axis=-1) / 2

From 3f46b6a198871442b8da2a4bb949ba420f4173e1 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 9 Mar 2023 09:21:06 -0500
Subject: [PATCH 14/17] test coverage

---
 pvlib/tests/bifacial/test_infinite_sheds.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/pvlib/tests/bifacial/test_infinite_sheds.py b/pvlib/tests/bifacial/test_infinite_sheds.py
index 404b0914c1..66eced25ad 100644
--- a/pvlib/tests/bifacial/test_infinite_sheds.py
+++ b/pvlib/tests/bifacial/test_infinite_sheds.py
@@ -42,7 +42,8 @@ def test_system():
     return syst, pts, vfs_ground_sky
 
 
-def test__vf_ground_sky_integ(test_system):
+@pytest.mark.parametrize("vectorize", [True, False])
+def test__vf_ground_sky_integ(test_system, vectorize):
     ts, pts, vfs_gnd_sky = test_system
     # pass rotation here since max_rows=1 for the hand-solved case in
     # the fixture test_system, which means the ground-to-sky view factor
@@ -50,7 +51,7 @@ def test__vf_ground_sky_integ(test_system):
     vf_integ = infinite_sheds._vf_ground_sky_integ(
         ts['rotation'], ts['surface_azimuth'],
         ts['gcr'], ts['height'], ts['pitch'],
-        max_rows=1, npoints=3)
+        max_rows=1, npoints=3, vectorize=vectorize)
     expected_vf_integ = np.trapz(vfs_gnd_sky, pts)
     assert np.isclose(vf_integ, expected_vf_integ, rtol=0.1)
 
@@ -262,7 +263,8 @@ def test__backside_tilt():
     assert np.allclose(back_az, np.array([0., 330., 90., 180.]))
 
 
-def test_get_irradiance():
+@pytest.mark.parametrize("vectorize", [True, False])
+def test_get_irradiance(vectorize):
     # singleton inputs
     solar_zenith = 0.
     solar_azimuth = 180.
@@ -282,7 +284,7 @@ def test_get_irradiance():
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back,
         bifaciality=0.8, shade_factor=-0.02, transmission_factor=0,
-        npoints=npoints)
+        npoints=npoints, vectorize=vectorize)
     expected_front_diffuse = np.array([300.])
     expected_front_direct = np.array([700.])
     expected_front_global = expected_front_diffuse + expected_front_direct
@@ -300,11 +302,11 @@ def test_get_irradiance():
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back,
         bifaciality=0.8, shade_factor=-0.02, transmission_factor=0,
-        npoints=npoints)
+        npoints=npoints, vectorize=vectorize)
     result_front = infinite_sheds.get_irradiance_poa(
         surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         gcr, height, pitch, ghi, dhi, dni,
-        albedo, iam=iam_front)
+        albedo, iam=iam_front, vectorize=vectorize)
     assert isinstance(result, pd.DataFrame)
     expected_poa_global = pd.Series(
         [1000., 500., result_front['poa_global'][2] * (1 + 0.8 * 0.98),

From d53197b2f5bc586fac936e193a5385ee11636647 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Thu, 9 Mar 2023 09:21:54 -0500
Subject: [PATCH 15/17] stickler

---
 pvlib/bifacial/infinite_sheds.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index d85aedc73e..bc64ab5c5e 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -53,9 +53,9 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    # calculate the integrated view factor for all of the ground between rows
     if vectorize:
-        fz_sky = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height, max_rows)
+        fz_sky = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height,
+                                         max_rows)
     else:
         fz_sky = np.zeros((npoints, len(rotation)))
         for k, r in enumerate(rotation):

From 63f6675d3d4cfcf9210224c2612a5aa8c885ba79 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Wed, 15 Mar 2023 19:31:04 -0400
Subject: [PATCH 16/17] Apply suggestions from code review

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 pvlib/bifacial/infinite_sheds.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index bc64ab5c5e..5b8d7fd974 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -12,7 +12,7 @@
 def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
                          pitch, max_rows=10, npoints=100, vectorize=False):
     """
-    Integrated view factor from the ground to the sky underneath
+    Integrated view factor to the sky from the ground underneath
     interior rows of the array.
 
     Parameters
@@ -469,7 +469,7 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
         on the surface that is not reflected away. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated viewfactors.
+        Number of discretization points for calculating integrated view factors.
 
     vectorize : bool, default False
         If True, vectorize the view factor calculation across ``surface_tilt``.
@@ -706,7 +706,7 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         etc. A negative value is a reduction in back irradiance. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated viewfactors.
+        Number of discretization points for calculating integrated view factors.
 
     vectorize : bool, default False
         If True, vectorize the view factor calculation across ``surface_tilt``.

From 514a1e98ccdb3ae51581de1f01a94cb7dce5a6b5 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Wed, 15 Mar 2023 19:35:33 -0400
Subject: [PATCH 17/17] stickler

---
 pvlib/bifacial/infinite_sheds.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py
index 5b8d7fd974..78e1a0f2a0 100644
--- a/pvlib/bifacial/infinite_sheds.py
+++ b/pvlib/bifacial/infinite_sheds.py
@@ -469,7 +469,8 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
         on the surface that is not reflected away. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated view factors.
+        Number of discretization points for calculating integrated view
+        factors.
 
     vectorize : bool, default False
         If True, vectorize the view factor calculation across ``surface_tilt``.
@@ -706,7 +707,8 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         etc. A negative value is a reduction in back irradiance. [unitless]
 
     npoints : int, default 100
-        Number of discretization points for calculating integrated view factors.
+        Number of discretization points for calculating integrated view
+        factors.
 
     vectorize : bool, default False
         If True, vectorize the view factor calculation across ``surface_tilt``.