Correct condition 5 in detect_clearksy (#596)

* Fixed inconsistencies between the detect_clearsky function and the original reference.

* Minor cleanup.

* Revert samples_per_interval change, rename to clarify

* Style fixes

* Style fixes

* whatsnew

* Variable name edits, delete extinct commented code lines

* Document change to API of detect_clearsky

Author: cwhanse

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

@@ -18,6 +18,7 @@ API Changes
   and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`)
 * Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem`
   and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``
+* Changed key names for `components` returned from :py:func:`pvlib.clearsky.detect_clearsky`
 
 
 Enhancements
@@ -35,6 +36,7 @@ Bug fixes
 ~~~~~~~~~
 * Fix when building documentation using Matplotlib 3.0 or greater.
 * Fix and improve :func:`~pvlib.solarposition.hour_angle` (:issue:`598`)
+* Fix error in :func:`pvlib.clearsky.detect_clearsky` (:issue:`506`)
 
 
 Testing
@@ -45,6 +47,7 @@ Testing
 Contributors
 ~~~~~~~~~~~~
 * Will Holmgren (:ghuser:`wholmgren`)
-* Cliff Hansen (:ghuser:`cwhanse`)
 * Leland Boeman (:ghuser:`lboeman`)
+* Ben Ellis (:ghuser:`bhellis725`)
+* Cliff Hansen (:ghuser:`cwhanse`)
 * Mark Mikofski (:ghuser:`mikofski`)

pvlib/clearsky.py

@@ -680,36 +680,37 @@ def detect_clearsky(measured, clearsky, times, window_length,
         raise NotImplementedError('algorithm does not yet support unequal '
                                   'times. consider resampling your data.')
 
-    samples_per_window = int(window_length / sample_interval)
+    intervals_per_window = int(window_length / sample_interval)
 
     # generate matrix of integers for creating windows with indexing
     from scipy.linalg import hankel
-    H = hankel(np.arange(samples_per_window),                   # noqa: N806
-               np.arange(samples_per_window-1, len(times)))
+    H = hankel(np.arange(intervals_per_window),                   # noqa: N806
+               np.arange(intervals_per_window - 1, len(times)))
 
     # calculate measurement statistics
     meas_mean = np.mean(measured[H], axis=0)
     meas_max = np.max(measured[H], axis=0)
-    meas_slope = np.diff(measured[H], n=1, axis=0)
+    meas_diff = np.diff(measured[H], n=1, axis=0)
+    meas_slope = np.diff(measured[H], n=1, axis=0) / sample_interval
     # matlab std function normalizes by N-1, so set ddof=1 here
     meas_slope_nstd = np.std(meas_slope, axis=0, ddof=1) / meas_mean
-    meas_slope_max = np.max(np.abs(meas_slope), axis=0)
     meas_line_length = np.sum(np.sqrt(
-        meas_slope*meas_slope + sample_interval*sample_interval), axis=0)
+        meas_diff * meas_diff +
+        sample_interval * sample_interval), axis=0)
 
     # calculate clear sky statistics
     clear_mean = np.mean(clearsky[H], axis=0)
     clear_max = np.max(clearsky[H], axis=0)
-    clear_slope = np.diff(clearsky[H], n=1, axis=0)
-    clear_slope_max = np.max(np.abs(clear_slope), axis=0)
+    clear_diff = np.diff(clearsky[H], n=1, axis=0)
+    clear_slope = np.diff(clearsky[H], n=1, axis=0) / sample_interval
 
     from scipy.optimize import minimize_scalar
 
     alpha = 1
     for iteration in range(max_iterations):
         clear_line_length = np.sum(np.sqrt(
-            alpha*alpha*clear_slope*clear_slope +
-            sample_interval*sample_interval), axis=0)
+            alpha * alpha * clear_diff * clear_diff +
+            sample_interval * sample_interval), axis=0)
 
         line_diff = meas_line_length - clear_line_length
 
@@ -718,7 +719,8 @@ def detect_clearsky(measured, clearsky, times, window_length,
         c2 = np.abs(meas_max - alpha*clear_max) < max_diff
         c3 = (line_diff > lower_line_length) & (line_diff < upper_line_length)
         c4 = meas_slope_nstd < var_diff
-        c5 = (meas_slope_max - alpha*clear_slope_max) < slope_dev
+        c5 = np.max(np.abs(meas_slope -
+                           alpha * clear_slope), axis=0) < slope_dev
         c6 = (clear_mean != 0) & ~np.isnan(clear_mean)
         clear_windows = c1 & c2 & c3 & c4 & c5 & c6
 
@@ -749,13 +751,21 @@ def rmse(alpha):
 
     if return_components:
         components = OrderedDict()
-        components['mean_diff'] = c1
-        components['max_diff'] = c2
-        components['line_length'] = c3
-        components['slope_nstd'] = c4
-        components['slope_max'] = c5
-        components['mean_nan'] = c6
+        components['mean_diff_flag'] = c1
+        components['max_diff_flag'] = c2
+        components['line_length_flag'] = c3
+        components['slope_nstd_flag'] = c4
+        components['slope_max_flag'] = c5
+        components['mean_nan_flag'] = c6
         components['windows'] = clear_windows
+
+        components['mean_diff'] = np.abs(meas_mean - alpha * clear_mean)
+        components['max_diff'] = np.abs(meas_max - alpha * clear_max)
+        components['line_length'] = meas_line_length - clear_line_length
+        components['slope_nstd'] = meas_slope_nstd
+        components['slope_max'] = (np.max(
+            meas_slope - alpha * clear_slope, axis=0))
+
         return clear_samples, components, alpha
     else:
         return clear_samples