Merge pull request #158 from slacgismo/loss-dev

Work on signal decompositions

solardatatools/_osd_signal_decompositions.py

@@ -27,21 +27,34 @@
 import numpy as np
 
 from gfosd import Problem
-from gfosd.components import SumAbs, SumSquare, SumCard, SumQuantile, Aggregate, AverageEqual,\
-    Periodic, Inequality, FirstValEqual, LastValEqual, NoCurvature, NoSlope
+from gfosd.components import (
+    SumAbs,
+    SumSquare,
+    SumCard,
+    SumQuantile,
+    Aggregate,
+    AverageEqual,
+    Periodic,
+    Inequality,
+    FirstValEqual,
+    LastValEqual,
+    NoCurvature,
+    NoSlope,
+    Fourier,
+)
 
 
 def _osd_l2_l1d1_l2d2p365(
-        signal,
-        w0=10,
-        w1=50,
-        w2=1e5,
-        return_all=False,
-        yearly_periodic=False,
-        solver="QSS",
-        use_ixs=None,
-        sum_card=False,
-        verbose=False
+    signal,
+    w0=10,
+    w1=50,
+    w2=1e5,
+    return_all=False,
+    yearly_periodic=False,
+    solver="QSS",
+    use_ixs=None,
+    sum_card=False,
+    verbose=False,
 ):
     """
     Used in: solardatatools/algorithms/time_shifts.py
@@ -70,7 +83,7 @@ def _osd_l2_l1d1_l2d2p365(
     :return: A tuple with two 1d numpy arrays containing the two signal component estimates
     """
     if solver != "QSS":
-        sum_card=False
+        sum_card = False
 
     if sum_card:
         # Scale objective
@@ -84,17 +97,18 @@ def _osd_l2_l1d1_l2d2p365(
         w2 /= 1e4
 
     c1 = SumSquare(weight=w0)
-
-    c2 = SumSquare(weight=w2, diff=2)
+    T = len(signal)
+    c2 = Fourier(3, T, 365.2425, weight=1e-3)
 
     if sum_card:
         c3 = SumCard(weight=w1, diff=1)
     else:
         c3 = SumAbs(weight=w1, diff=1)
 
     if len(signal) >= 365:
-        c2 = Aggregate([SumSquare(weight=w2, diff=2), AverageEqual(0, period=365), Periodic(365)])
-        if yearly_periodic and not sum_card: # SumCard does not work well with Aggregate class
+        if (
+            yearly_periodic and not sum_card
+        ):  # SumCard does not work well with Aggregate class
             c3 = Aggregate([c3, Periodic(365)])
         elif yearly_periodic and sum_card:
             print("Cannot use Periodic Class with SumCard.")
@@ -104,7 +118,7 @@ def _osd_l2_l1d1_l2d2p365(
     problem = Problem(signal, classes, use_set=use_ixs)
     problem.decompose(solver=solver, verbose=verbose, eps_rel=1e-6, eps_abs=1e-6)
 
-    s_error =  problem.decomposition[0]
+    s_error = problem.decomposition[0]
     s_seas = problem.decomposition[1]
     s_hat = problem.decomposition[2]
 
@@ -115,15 +129,14 @@ def _osd_l2_l1d1_l2d2p365(
 
 
 def _osd_tl1_l2d2p365(
-        signal,
-        use_ixs=None,
-        tau=0.75,
-        w0=1,
-        w1=500,
-        yearly_periodic=True,
-        return_all=False,
-        solver="OSQP",
-        verbose=False
+    signal,
+    use_ixs=None,
+    tau=0.75,
+    w0=1,
+    yearly_periodic=True,
+    return_all=False,
+    solver="OSQP",
+    verbose=False,
 ):
     """
     Used in:
@@ -151,10 +164,8 @@ def _osd_tl1_l2d2p365(
     :return: A tuple with three 1d numpy arrays containing the three signal component estimates
     """
     c1 = SumQuantile(tau=tau, weight=w0)
-    c2 = SumSquare(weight=w1, diff=2)
-
-    if len(signal) > 365 and yearly_periodic:
-        c2 = Aggregate([c2, Periodic(365)])
+    T = len(signal)
+    c2 = Fourier(3, T, 365.2425, weight=1e-3)
 
     classes = [c1, c2]
 
@@ -173,12 +184,12 @@ def _osd_l1_l1d1_l2d2p365(
     signal,
     use_ixs=None,
     w0=2e-6,  # l1 term, scaled
-    w1=40e-6, # l1d1 term, scaled
-    w2=6e-3, # seasonal term, scaled
+    w1=40e-6,  # l1d1 term, scaled
+    w2=6e-3,  # seasonal term, scaled
     return_all=False,
     solver=None,
     sum_card=False,
-    verbose=False
+    verbose=False,
 ):
     """
     Used in solardatatools/algorithms/capacity_change.py
@@ -202,29 +213,25 @@ def _osd_l1_l1d1_l2d2p365(
     :param verbose: Sets verbosity
     :return: A tuple with three 1d numpy arrays containing the three signal component estimates
     """
-    if solver!="QSS":
-        sum_card=False
+    if solver != "QSS":
+        sum_card = False
 
     c1 = SumAbs(weight=w0)
-
-    c2 = SumSquare(weight=w2, diff=2)
+    T = len(signal)
+    c2 = Fourier(3, T, 365.2425, weight=1e-3)
+    # TODO: evaluate the weight used here
     if len(signal) >= 365:
-        c2 = Aggregate([c2,
-                        AverageEqual(0, period=365),
-                        Periodic(365)
-                        ])
+        pass
     else:
-        w1 /= 5 # PWC weight needs adjusting when dataset is short
+        w1 /= 5  # PWC weight needs adjusting when dataset is short
 
     if sum_card:
         c3 = SumCard(weight=w1, diff=1)
     else:
         c3 = SumAbs(weight=w1, diff=1)
 
     # Linear term to describe yearly degradation of seasonal component
-    c4 =  Aggregate([NoCurvature(),
-                     FirstValEqual(0)
-                     ])
+    c4 = Aggregate([NoCurvature(), FirstValEqual(0)])
 
     classes = [c1, c2, c3, c4]
 
@@ -242,12 +249,7 @@ def _osd_l1_l1d1_l2d2p365(
 
 
 def _osd_l2_l1d2_constrained(
-        signal,
-        w0=1,
-        w1=5,
-        return_all=False,
-        solver="OSQP",
-        verbose=False
+    signal, w0=1, w1=5, return_all=False, solver="OSQP", verbose=False
 ):
     """
     Used in solardatatools/algorithms/clipping.py
@@ -265,11 +267,7 @@ def _osd_l2_l1d2_constrained(
      and the weight
     """
     c1 = SumSquare(weight=w0)
-    c2 = Aggregate([
-        SumAbs(weight=w1, diff=2),
-        FirstValEqual(0),
-        LastValEqual(1)
-    ])
+    c2 = Aggregate([SumAbs(weight=w1, diff=2), FirstValEqual(0), LastValEqual(1)])
 
     classes = [c1, c2]
 

solardatatools/algorithms/loss_factor_analysis.py

@@ -16,6 +16,8 @@
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
+import scipy.sparse as sp
+from scipy.signal import sawtooth, find_peaks
 from gfosd import Problem
 import gfosd.components as comp
 from gfosd.components.base_graph_class import GraphComponent
@@ -132,7 +134,7 @@ def estimate_degradation_rate(
             if verbose:
                 progress.update()
             tau = np.random.uniform(0.85, 0.95)
-            weight = np.random.uniform(0.5, 5)
+            weight = np.random.uniform(0.1, 1)
             good_ixs = np.arange(len(self.use_ixs))[self.use_ixs]
             num_values = int(len(good_ixs) * fraction_hold)
             # select random index locations for removal
@@ -142,7 +144,18 @@ def estimate_degradation_rate(
             new_ixs = ~np.logical_or(~self.use_ixs, msk)
             self.use_ixs = new_ixs
             # remake modified problem
-            self.problem = self.make_problem(tau=tau, weight_soiling_stiffness=weight)
+            self.problem = self.make_problem(
+                tau=tau,
+                weight_soiling_stiffness=weight,
+                num_harmonics=self.user_settings["num_harmonics"],
+                deg_type=self.user_settings["deg_type"],
+                include_soiling=self.user_settings["include_soiling"],
+                weight_seasonal=self.user_settings["weight_seasonal"],
+                weight_soiling_sparsity=self.user_settings["weight_soiling_sparsity"],
+                use_capacity_change_labels=self.user_settings[
+                    "use_capacity_change_labels"
+                ],
+            )
             # run signal decomposition and shapley attribution
             self.estimate_losses()
             # record running results
@@ -297,16 +310,21 @@ def plot_waterfall(self):
         fig = waterfall_plot(data, index)
         return fig
 
-    def plot_decomposition(self, figsize=(16, 8.5)):
+    def plot_decomposition(self, plot_capacity_component=True, figsize=(16, 8.5)):
         """
         Creates a figure with subplots illustrating the estimated signal components found through decomposition
 
         :param figsize: size of figure (tuple)
         :return: matplotlib figure
         """
-        _fig_decomp = self.problem.plot_decomposition(
-            exponentiate=True, figsize=figsize
-        )
+        if plot_capacity_component:
+            _fig_decomp = self.problem.plot_decomposition(
+                exponentiate=True, figsize=figsize
+            )
+        else:
+            _fig_decomp = self.problem.plot_decomposition(
+                exponentiate=True, figsize=figsize, skip=1
+            )
         _ax = _fig_decomp.axes
         _ax[0].plot(
             np.arange(len(self.energy_data))[~self.use_ixs],
@@ -315,12 +333,19 @@ def plot_decomposition(self, figsize=(16, 8.5)):
             marker=".",
             ls="none",
         )
-        _ax[0].set_title("weather and system outages")
-        _ax[1].set_title("capacity changes")
-        _ax[2].set_title("soiling")
-        _ax[3].set_title("degradation")
-        _ax[4].set_title("baseline")
-        _ax[5].set_title("measured energy (green) and model minus weather")
+        if plot_capacity_component:
+            _ax[0].set_title("weather and system outages")
+            _ax[1].set_title("capacity changes")
+            _ax[2].set_title("soiling")
+            _ax[3].set_title("degradation")
+            _ax[4].set_title("baseline")
+            _ax[5].set_title("measured energy (green) and model minus weather")
+        else:
+            _ax[0].set_title("weather and system outages")
+            _ax[1].set_title("soiling")
+            _ax[2].set_title("degradation")
+            _ax[3].set_title("baseline")
+            _ax[4].set_title("measured energy (green) and model minus weather")
         plt.tight_layout()
         return _fig_decomp
 
@@ -408,10 +433,11 @@ def make_problem(
         deg_type="linear",
         include_soiling=True,
         weight_seasonal=10e-2,
-        weight_soiling_stiffness=1e0,
+        weight_soiling_stiffness=5e-1,
         weight_soiling_sparsity=1e-2,
         weight_deg_nonlinear=10e4,
         deg_rate=None,
+        use_capacity_change_labels=True,
     ):
         """
         Constuct the signal decomposition problem for estimation of loss factors in PV energy data.
@@ -446,18 +472,18 @@ def make_problem(
         c1 = comp.SumQuantile(tau=tau)
         # Smooth periodic term
         length = len(self.log_energy)
-        periods = [365.2425]  # average length of a year in days
-        _B = make_basis_matrix(num_harmonics, length, periods)
-        _D = make_regularization_matrix(num_harmonics, weight_seasonal, periods)
-        c2 = comp.Basis(basis=_B, penalty=_D)
+        period = 365.2425  # average length of a year in days
+        c2 = comp.Fourier(num_harmonics, length, period, weight=weight_seasonal)
         # Soiling term
         if include_soiling:
+            sawtooth_dict = make_sawtooth_dictionary(length)
             c3 = comp.Aggregate(
                 [
                     comp.Inequality(vmax=0),
-                    comp.SumAbs(weight=weight_soiling_stiffness, diff=2),
-                    comp.SumQuantile(
-                        tau=0.98, weight=10 * weight_soiling_sparsity, diff=1
+                    comp.Basis(
+                        basis=sawtooth_dict,
+                        penalty="abs",
+                        weight=weight_soiling_stiffness,
                     ),
                     comp.SumAbs(weight=weight_soiling_sparsity),
                 ]
@@ -488,8 +514,8 @@ def make_problem(
             )
         elif deg_select.value == "none":
             c4 = SetEqual(val=np.zeros(length))
-        # capacity change term — leverage previous analysis from SDT pipeline
-        if self.capacity_change_labels is not None:
+        # capacity change term: leverage previous analysis from SDT pipeline
+        if use_capacity_change_labels and self.capacity_change_labels is not None:
             basis_M = np.zeros((length, len(set(self.capacity_change_labels))))
             for lb in set(self.capacity_change_labels):
                 slct = np.array(self.capacity_change_labels) == lb
@@ -692,3 +718,27 @@ def _make_A(self):
 
     def _make_c(self):
         self._c = self._val
+
+
+def make_sawtooth_dictionary(T):
+    ks = np.arange(2, 32)
+    phases = [0, np.pi]
+    dictionary = [-1 * np.ones(T).reshape((-1, 1))]
+    for _k in ks:
+        for ph in phases:
+            dictionary.append(make_st(_k, ph, T))
+    dictionary = np.concatenate(dictionary, axis=1)
+    return sp.lil_array(dictionary)
+
+
+def make_st(k, phase, t):
+    start = -phase
+    end = k * 2 * np.pi - phase
+    wf = sawtooth(np.linspace(start, end, t), 0) / 2 - 0.5
+    peaks = find_peaks(wf)
+    indices = np.r_[np.s_[None], peaks[0], np.s_[None]]
+    num_segments = len(peaks[0]) + 1
+    out = np.zeros((t, num_segments))
+    for s in np.arange(num_segments):
+        out[indices[s] : indices[s + 1], s] = wf[indices[s] : indices[s + 1]]
+    return out