Fix to WVM methodology to track changes in MATLAB version.

docs/sphinx/source/whatsnew/v0.9.0.rst

@@ -117,6 +117,8 @@ Bug fixes
 * Reindl model fixed to generate sky_diffuse=0 when GHI=0.
   (:issue:`1153`, :pull:`1154`)
 * Update GFS product names for GFS v16. (:issue:`1202`, :pull:`1203`)
+* Corrected methodology error in scaling.py WVM. Tracks with fix in MATLAB ver.
+  (:issue:`1206`, :pull:`1213`)
 
 Testing
 ~~~~~~~
@@ -148,3 +150,4 @@ Contributors
 * Joshua Stein (:ghuser:`jsstein`)
 * Tony Lorenzo (:ghuser:`alorenzo175`)
 * Damjan Postolovski (:ghuser:`dpostolovski`)
+* Joe Ranalli (:ghuser:`jranalli`)

pvlib/scaling.py

@@ -159,7 +159,12 @@ def latlon_to_xy(coordinates):
 
 def _compute_wavelet(clearsky_index, dt=None):
     """
-    Compute the wavelet transform on the input clear_sky time series.
+    Compute the wavelet transform on the input clear_sky time series. Uses a
+    top hat wavelet [-1,1,1,-1] shape, based on the difference of successive
+    centered moving averages. Smallest scale (filter size of 2) is a degenerate
+    case that resembles a Haar wavelet. Returns one level of approximation
+    coefficient (CAn) and n levels of detail coefficients (CD1, CD2, ...,
+    CDn-1, CDn).
 
     Parameters
     ----------
@@ -174,7 +179,8 @@ def _compute_wavelet(clearsky_index, dt=None):
     Returns
     -------
     wavelet: numeric
-        The individual wavelets for the time series
+        The individual wavelets for the time series. Format follows increasing
+        scale (decreasing frequency): [CD1, CD2, ..., CDn, CAn]
 
     tmscales: numeric
         The timescales associated with the wavelets in seconds [s]
@@ -185,7 +191,7 @@ def _compute_wavelet(clearsky_index, dt=None):
     Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
     Energy, vol. 4, no. 2, pp. 501-509, 2013.
 
-    [3] Wavelet Variability Model - Matlab Code:
+    [2] Wavelet Variability Model - Matlab Code:
     https://pvpmc.sandia.gov/applications/wavelet-variability-model/
     """
 
@@ -209,31 +215,37 @@ def _compute_wavelet(clearsky_index, dt=None):
 
     # Compute wavelet time scales
     min_tmscale = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet timescale
-    max_tmscale = int(12 - min_tmscale)  # maximum wavelet timescale
+    max_tmscale = int(13 - min_tmscale)  # maximum wavelet timescale
 
     tmscales = np.zeros(max_tmscale)
     csi_mean = np.zeros([max_tmscale, len(cs_long)])
+    # Skip averaging for the 0th scale
+    csi_mean[0, :] = cs_long.values.flatten()
+    tmscales[0] = 1
     # Loop for all time scales we will consider
-    for i in np.arange(0, max_tmscale):
-        j = i+1
-        tmscales[i] = 2**j * dt  # Wavelet integration time scale
-        intvlen = 2**j  # Wavelet integration time series interval
+    for i in np.arange(1, max_tmscale):
+        tmscales[i] = 2**i * dt  # Wavelet integration time scale
+        intvlen = 2**i  # Wavelet integration time series interval
         # Rolling average, retains only lower frequencies than interval
+        # Produces slightly different end effects than the MATLAB version
         df = cs_long.rolling(window=intvlen, center=True, min_periods=1).mean()
         # Fill nan's in both directions
         df = df.fillna(method='bfill').fillna(method='ffill')
         # Pop values back out of the dataframe and store
         csi_mean[i, :] = df.values.flatten()
+        # Shift to account for different indexing in MATLAB moving average
+        csi_mean[i, :] = np.roll(csi_mean[i, :], -1)
+        csi_mean[i, -1] = csi_mean[i, -2]
 
-    # Calculate the wavelets by isolating the rolling mean frequency ranges
+    # Calculate detail coefficients by difference between successive averages
     wavelet_long = np.zeros(csi_mean.shape)
     for i in np.arange(0, max_tmscale-1):
         wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :]
-    wavelet_long[max_tmscale-1, :] = csi_mean[max_tmscale-1, :]  # Lowest freq
+    wavelet_long[-1, :] = csi_mean[-1, :]  # Lowest freq (CAn)
 
     # Clip off the padding and just return the original time window
     wavelet = np.zeros([max_tmscale, len(vals)])
     for i in np.arange(0, max_tmscale):
-        wavelet[i, :] = wavelet_long[i, len(vals)+1: 2*len(vals)+1]
+        wavelet[i, :] = wavelet_long[i, len(vals): 2*len(vals)]
 
     return wavelet, tmscales

pvlib/tests/test_scaling.py

@@ -48,21 +48,24 @@ def positions():
 @pytest.fixture
 def expect_tmscale():
     # Expected timescales for dt = 1
-    return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+    return [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
 
 
 @pytest.fixture
 def expect_wavelet():
     # Expected wavelet for indices 5000:5004 for clear_sky_index above (Matlab)
-    return np.array([[-0.025, 0.05, 0., -0.05, 0.025],
-                     [0.025, 0., 0., 0., -0.025],
-                     [0., 0., 0., 0., 0.]])
+    e = np.zeros([13, 5])
+    e[0, :] = np.array([0, -0.05, 0.1, -0.05, 0])
+    e[1, :] = np.array([-0.025, 0.05, 0., -0.05, 0.025])
+    e[2, :] = np.array([0.025, 0., 0., 0., -0.025])
+    e[-1, :] = np.array([1, 1, 1, 1, 1])
+    return e
 
 
 @pytest.fixture
 def expect_cs_smooth():
     # Expected smoothed clear sky index for indices 5000:5004 (Matlab)
-    return np.array([1., 1.0289, 1., 0.9711, 1.])
+    return np.array([1., 1., 1.05774, 0.94226, 1.])
 
 
 def test_latlon_to_xy_zero():
@@ -94,7 +97,7 @@ def test_compute_wavelet_series(clear_sky_index, time,
     csi_series = pd.Series(clear_sky_index, index=time)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_series_numindex(clear_sky_index, time,
@@ -103,21 +106,29 @@ def test_compute_wavelet_series_numindex(clear_sky_index, time,
     csi_series = pd.Series(clear_sky_index, index=dtindex)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_array(clear_sky_index,
                                expect_tmscale, expect_wavelet):
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_array_invalid(clear_sky_index):
     with pytest.raises(ValueError):
         scaling._compute_wavelet(clear_sky_index)
 
 
+def test_compute_wavelet_dwttheory(clear_sky_index, time,
+                                   expect_tmscale, expect_wavelet):
+    # Confirm detail coeffs sum to original signal
+    csi_series = pd.Series(clear_sky_index, index=time)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(np.sum(wavelet, 0), csi_series)
+
+
 def test_wvm_series(clear_sky_index, time, positions, expect_cs_smooth):
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, positions, cloud_speed)