Dev

gfosd/components/__init__.py

@@ -2,9 +2,10 @@
                                    SumAbs,
                                    SumHuber,
                                    SumQuantile,
-                                   SumCard)
+                                   SumCard,
+                                   SumSquareReference)
 from gfosd.components.inequality_constraints import Inequality
-from gfosd.components.basis_constraints import Basis, Periodic
+from gfosd.components.basis_constraints import Basis, Periodic, Fourier
 from gfosd.components.finite_set import FiniteSet, Boolean
 from gfosd.components.aggregate import Aggregate
 from gfosd.components.equality_constraints import (FirstValEqual,

gfosd/components/aggregate.py

@@ -37,7 +37,10 @@ def prepare_attributes(self, T, p=1):
         self._Pz = sp.block_diag([
             c._Pz for c in self._gf_list
         ])
-        self._make_q()
+        # same logic as above but for q
+        qx = self._gf_list[g_ix]._q[:self._gf_list[g_ix].x_size]
+        qz = np.concatenate([c._q[self._gf_list[g_ix].x_size:] for c in self._gf_list])
+        self._q = np.concatenate([qx, qz])
         self._make_r()
         self._gx = self._make_gx()
         self._gz = self._make_gz()
@@ -53,6 +56,9 @@ def _make_P(self):
             c._Pz for c in self._gf_list
         ])
 
+    def _make_r(self):
+        self._r = np.sum([c._r for c in self._gf_list])
+
     def _make_gx(self):
         gx = []
         for ix, component in enumerate(self._gf_list):

gfosd/components/base_graph_class.py

@@ -94,10 +94,10 @@ def _make_P(self, size):
         return sp.dok_matrix(2 * (size,))
 
     def _make_q(self):
-        self._q = None
+        self._q = np.zeros(self.x_size + self.z_size)
 
     def _make_r(self):
-        self._r = None
+        self._r = 0
 
     def _make_g(self, size):
         return []

gfosd/components/basis_constraints.py

@@ -11,10 +11,11 @@
 
 """
 
-
 import numpy as np
 import scipy.sparse as sp
 from gfosd.components.base_graph_class import GraphComponent
+from gfosd.components.utilities import make_basis_matrix, make_regularization_matrix
+
 
 class Basis(GraphComponent):
     def __init__(self, basis, penalty=None, *args, **kwargs):
@@ -51,8 +52,8 @@ def _make_g(self, size):
         else:
             # typically 'abs', 'huber', or 'quantile'
             g = [{'g': self._penalty,
-                         'args': {'weight': self.weight},
-                         'range': (0, size)}]
+                  'args': {'weight': self.weight},
+                  'range': (0, size)}]
         return g
 
     def _make_P(self, size):
@@ -65,8 +66,32 @@ def _make_P(self, size):
             P = sp.dok_matrix(2 * (size,))
         return P
 
+
 class Periodic(Basis):
     def __init__(self, period, *args, **kwargs):
         self._period = period
         M = sp.eye(period)
         super().__init__(M, *args, **kwargs)
+
+
+class Fourier(Basis):
+    def __init__(self, num_harmonics, length, periods, standing_wave=False, trend=False, max_cross_k=None,
+                 custom_basis=None, weight=1, **kwargs):
+        _B = make_basis_matrix(
+            num_harmonics,
+            length,
+            periods,
+            standing_wave,
+            trend,
+            max_cross_k,
+            custom_basis)
+        _D = make_regularization_matrix(
+            num_harmonics,
+            weight,
+            periods,
+            standing_wave,
+            trend,
+            max_cross_k,
+            custom_basis
+        )
+        super().__init__(basis=_B, penalty=_D, weight=weight, **kwargs)

gfosd/components/equality_constraints.py

@@ -53,6 +53,34 @@ def _make_B(self):
     def _make_c(self):
         self._c = np.concatenate([np.atleast_1d(self._c[-1]),
                                   [self._last_val]])
+
+
+class ValsEqual(GraphComponent):
+    def __init__(self, indices=0, value=0, *args, **kwargs):
+        self.indices = np.atleast_1d(indices)
+        self._val = value
+        super().__init__(*args, **kwargs)
+        # always retain helper variable
+        self._has_helpers = True
+
+    def _make_A(self):
+        super()._make_A()
+        super()._make_B()
+        super()._make_c()
+        self._A = sp.bmat([
+           [self._A.tocsr()[-1]],
+            [sp.dok_matrix((len(self.indices), self._A.shape[1]))]
+        ])
+
+    def _make_B(self):
+        self._B = sp.bmat([
+            [self._B.tocsr()[-1]], #XXXXX got to here! need to update this matrix generating function to identify all indices that are to be set equal
+            [sp.coo_matrix(([1], ([0], [self._B.shape[1]-1])), shape=(len(self.indices), self._B.shape[1]))]
+        ])
+
+    def _make_c(self):
+        self._c = np.concatenate([np.atleast_1d(self._c[-1]),
+                                  [self.val] * len(self.indices)])
 
 class AverageEqual(GraphComponent):
     def __init__(self, value=0, period=None, *args, **kwargs):

gfosd/components/sums.py

@@ -1,6 +1,8 @@
+import numpy as np
 import scipy.sparse as sp
 from gfosd.components.base_graph_class import GraphComponent
 
+
 class SumSquare(GraphComponent):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -9,17 +11,65 @@ def __init__(self, *args, **kwargs):
     def _make_P(self, size):
         return self.weight * 2 * sp.eye(size)  # note the (1/2) in canonical form!
 
+
+class SumSquareReference(GraphComponent):
+    """
+    phi(x;v) = || x - v ||_2^2 = ||x||_2^2 -2v^Tx + ||v||_2^2
+    """
+
+    def __init__(self, vec, mask=None, *args, **kwargs):
+        self._vec = vec
+        self._mask = mask
+        super().__init__(*args, **kwargs)
+        if self._mask is not None:
+            self._has_helpers = True
+
+    def _set_z_size(self):
+        if self._mask is None:
+            self._z_size = (self._T - self._diff) * self._p
+        else:
+            self._z_size = self._mask.shape[0]
+
+    def _make_P(self, size):
+        return self.weight * 2 * sp.eye(size)  # note the (1/2) in canonical form!
+
+    def _make_q(self):
+        if self._mask is None:
+            qx = self.weight * (-2) * self._vec
+            self._q = np.r_[qx, np.zeros(self.z_size)]
+        else:
+            qz = (-2) * self.weight * self._mask @ self._vec
+            self._q = np.r_[np.zeros(self.x_size), qz]
+
+    def _make_r(self):
+        if self._mask is None:
+            self._r = self.weight * np.sum(np.square(self._vec))
+        else:
+            self._r = self.weight * np.sum(np.square(self._mask @ self._vec))
+
+    def _make_A(self):
+        self._A = self._mask
+
+
 class SumAbs(GraphComponent):
-    def __init__(self, *args, **kwargs):
+    def __init__(self, weighting_mat=None, *args, **kwargs):
+        self.weighting_mat = weighting_mat
         super().__init__(*args, **kwargs)
         return
 
+    def _make_B(self):
+        if self.weighting_mat is None:
+            self._B = sp.eye(self._A.shape[0], self.z_size) * -1
+        else:
+            self._B = self.weighting_mat
+
     def _make_g(self, size):
         g = [{'g': 'abs',
               'args': {'weight': self.weight},
               'range': (0, size)}]
         return g
 
+
 class SumHuber(GraphComponent):
     def __init__(self, M=1, *args, **kwargs):
         self._M = M
@@ -32,6 +82,7 @@ def _make_g(self, size):
               'range': (0, size)}]
         return g
 
+
 class SumQuantile(GraphComponent):
     def __init__(self, tau, *args, **kwargs):
         self.tau = tau
@@ -44,6 +95,7 @@ def _make_g(self, size):
               'range': (0, size)}]
         return g
 
+
 class SumCard(GraphComponent):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)