def detect_clearsky(measured, clearsky, times, window_length,<EOL>mean_diff=<NUM_LIT>, max_diff=<NUM_LIT>,<EOL>lower_line_length=-<NUM_LIT:5>, upper_line_length=<NUM_LIT:10>,<EOL>var_diff=<NUM_LIT>, slope_dev=<NUM_LIT:8>, max_iterations=<NUM_LIT:20>,<EOL>return_components=False):
<EOL>deltas = np.diff(times.values) / np.timedelta64(<NUM_LIT:1>, '<STR_LIT>')<EOL>unique_deltas = np.unique(deltas)<EOL>if len(unique_deltas) == <NUM_LIT:1>:<EOL><INDENT>sample_interval = unique_deltas[<NUM_LIT:0>]<EOL><DEDENT>else:<EOL><INDENT>raise NotImplementedError('<STR_LIT>'<EOL>'<STR_LIT>')<EOL><DEDENT>intervals_per_window = int(window_length / sample_interval)<EOL>from scipy.linalg import hankel<EOL>H = hankel(np.arange(intervals_per_window),                   <EOL>np.arange(intervals_per_window - <NUM_LIT:1>, len(times)))<EOL>meas_mean = np.mean(measured[H], axis=<NUM_LIT:0>)<EOL>meas_max = np.max(measured[H], axis=<NUM_LIT:0>)<EOL>meas_diff = np.diff(measured[H], n=<NUM_LIT:1>, axis=<NUM_LIT:0>)<EOL>meas_slope = np.diff(measured[H], n=<NUM_LIT:1>, axis=<NUM_LIT:0>) / sample_interval<EOL>meas_slope_nstd = np.std(meas_slope, axis=<NUM_LIT:0>, ddof=<NUM_LIT:1>) / meas_mean<EOL>meas_line_length = np.sum(np.sqrt(<EOL>meas_diff * meas_diff +<EOL>sample_interval * sample_interval), axis=<NUM_LIT:0>)<EOL>clear_mean = np.mean(clearsky[H], axis=<NUM_LIT:0>)<EOL>clear_max = np.max(clearsky[H], axis=<NUM_LIT:0>)<EOL>clear_diff = np.diff(clearsky[H], n=<NUM_LIT:1>, axis=<NUM_LIT:0>)<EOL>clear_slope = np.diff(clearsky[H], n=<NUM_LIT:1>, axis=<NUM_LIT:0>) / sample_interval<EOL>from scipy.optimize import minimize_scalar<EOL>alpha = <NUM_LIT:1><EOL>for iteration in range(max_iterations):<EOL><INDENT>clear_line_length = np.sum(np.sqrt(<EOL>alpha * alpha * clear_diff * clear_diff +<EOL>sample_interval * sample_interval), axis=<NUM_LIT:0>)<EOL>line_diff = meas_line_length - clear_line_length<EOL>c1 = np.abs(meas_mean - alpha*clear_mean) < mean_diff<EOL>c2 = np.abs(meas_max - alpha*clear_max) < max_diff<EOL>c3 = (line_diff > lower_line_length) & (line_diff < upper_line_length)<EOL>c4 = meas_slope_nstd < var_diff<EOL>c5 = np.max(np.abs(meas_slope -<EOL>alpha * clear_slope), axis=<NUM_LIT:0>) < slope_dev<EOL>c6 = (clear_mean != <NUM_LIT:0>) & ~np.isnan(clear_mean)<EOL>clear_windows = c1 & c2 & c3 & c4 & c5 & c6<EOL>clear_samples = np.full_like(measured, False, dtype='<STR_LIT:bool>')<EOL>clear_samples[np.unique(H[:, clear_windows])] = True<EOL>previous_alpha = alpha<EOL>clear_meas = measured[clear_samples]<EOL>clear_clear = clearsky[clear_samples]<EOL>def rmse(alpha):<EOL><INDENT>return np.sqrt(np.mean((clear_meas - alpha*clear_clear)**<NUM_LIT:2>))<EOL><DEDENT>alpha = minimize_scalar(rmse).x<EOL>if round(alpha*<NUM_LIT>) == round(previous_alpha*<NUM_LIT>):<EOL><INDENT>break<EOL><DEDENT><DEDENT>else:<EOL><INDENT>import warnings<EOL>warnings.warn('<STR_LIT>'<EOL>% max_iterations, RuntimeWarning)<EOL><DEDENT>if isinstance(measured, pd.Series):<EOL><INDENT>clear_samples = pd.Series(clear_samples, index=times)<EOL><DEDENT>if return_components:<EOL><INDENT>components = OrderedDict()<EOL>components['<STR_LIT>'] = c1<EOL>components['<STR_LIT>'] = c2<EOL>components['<STR_LIT>'] = c3<EOL>components['<STR_LIT>'] = c4<EOL>components['<STR_LIT>'] = c5<EOL>components['<STR_LIT>'] = c6<EOL>components['<STR_LIT>'] = clear_windows<EOL>components['<STR_LIT>'] = np.abs(meas_mean - alpha * clear_mean)<EOL>components['<STR_LIT>'] = np.abs(meas_max - alpha * clear_max)<EOL>components['<STR_LIT>'] = meas_line_length - clear_line_length<EOL>components['<STR_LIT>'] = meas_slope_nstd<EOL>components['<STR_LIT>'] = (np.max(<EOL>meas_slope - alpha * clear_slope, axis=<NUM_LIT:0>))<EOL>return clear_samples, components, alpha<EOL><DEDENT>else:<EOL><INDENT>return clear_samples<EOL><DEDENT>
Detects clear sky times according to the algorithm developed by Reno
and Hansen for GHI measurements [1]. The algorithm was designed and
validated for analyzing GHI time series only. Users may attempt to
apply it to other types of time series data using different filter
settings, but should be skeptical of the results.

The algorithm detects clear sky times by comparing statistics for a
measured time series and an expected clearsky time series.
Statistics are calculated using a sliding time window (e.g., 10
minutes). An iterative algorithm identifies clear periods, uses the
identified periods to estimate bias in the clearsky data, scales the
clearsky data and repeats.

Clear times are identified by meeting 5 criteria. Default values for
these thresholds are appropriate for 10 minute windows of 1 minute
GHI data.

Parameters
----------
measured : array or Series
    Time series of measured values.
clearsky : array or Series
    Time series of the expected clearsky values.
times : DatetimeIndex
    Times of measured and clearsky values.
window_length : int
    Length of sliding time window in minutes. Must be greater than 2
    periods.
mean_diff : float, default 75
    Threshold value for agreement between mean values of measured
    and clearsky in each interval, see Eq. 6 in [1].
max_diff : float, default 75
    Threshold value for agreement between maxima of measured and
    clearsky values in each interval, see Eq. 7 in [1].
lower_line_length : float, default -5
    Lower limit of line length criterion from Eq. 8 in [1].
    Criterion satisfied when
    lower_line_length < line length difference < upper_line_length
upper_line_length : float, default 10
    Upper limit of line length criterion from Eq. 8 in [1].
var_diff : float, default 0.005
    Threshold value in Hz for the agreement between normalized
    standard deviations of rate of change in irradiance, see Eqs. 9
    through 11 in [1].
slope_dev : float, default 8
    Threshold value for agreement between the largest magnitude of
    change in successive values, see Eqs. 12 through 14 in [1].
max_iterations : int, default 20
    Maximum number of times to apply a different scaling factor to
    the clearsky and redetermine clear_samples. Must be 1 or larger.
return_components : bool, default False
    Controls if additional output should be returned. See below.

Returns
-------
clear_samples : array or Series
    Boolean array or Series of whether or not the given time is
    clear. Return type is the same as the input type.

components : OrderedDict, optional
    Dict of arrays of whether or not the given time window is clear
    for each condition. Only provided if return_components is True.

alpha : scalar, optional
    Scaling factor applied to the clearsky_ghi to obtain the
    detected clear_samples. Only provided if return_components is
    True.

References
----------
[1] Reno, M.J. and C.W. Hansen, "Identification of periods of clear
sky irradiance in time series of GHI measurements" Renewable Energy,
v90, p. 520-531, 2016.

Notes
-----
Initial implementation in MATLAB by Matthew Reno. Modifications for
computational efficiency by Joshua Patrick and Curtis Martin. Ported
to Python by Will Holmgren, Tony Lorenzo, and Cliff Hansen.

Differences from MATLAB version:

    * no support for unequal times
    * automatically determines sample_interval
    * requires a reference clear sky series instead calculating one
      from a user supplied location and UTCoffset
    * parameters are controllable via keyword arguments
    * option to return individual test components and clearsky scaling
      parameter