Changes to Difference constraints

src/mssm/src/python/constraints.py

@@ -38,7 +38,8 @@ class Constraint:
   Line 3 shows how the constraint can be absorbed for model fitting by first computing new bases (f_i(x) - f_{i-1}(x)) and then estimating
   the coefficients based on those (this is implemented in mgcv's ``smoothCon`` function when setting ``sparse.cons=2``). Note that the constrained
   model, just like when using the QR-based or dropping approach, requires dropping one of the k original coefficients for estimation since only k-1
-  coefficients are allowed to vary.
+  coefficients are allowed to vary. The same sum-to-zero constraint can be achieved by fixing one of the central bases in the original model to it's neighbor (e.g., setting b_2 = -b_3) or
+  by setting b_1= -b_2 - b_3 - b_4 - b_5. mssm fixes one of the central bases to it's neighbor.
 
   With a B-splines basis, it would be necessary that b_1 = b_2 = b_3 = b_4 = b_5 for the model to fit a constant f(x) over all values of x (Eilers & Marx, 1996).
   In the constrained model this is no longer possible due to b_1 = -b_2, effectively removing the constant from the space of functions that f(x) can approximate.

src/mssm/src/python/formula.py

@@ -773,7 +773,7 @@ def __absorb_constraints(self):
                elif term_constraint == ConstType.DROP:
                   sterm.Z.append(Constraint(int(matrix_term_v.shape[1]/2),ConstType.DROP))
                elif term_constraint == ConstType.DIFF:
-                  sterm.Z.append(Constraint(None,ConstType.DIFF))
+                  sterm.Z.append(Constraint(int(matrix_term_v.shape[1]/2),ConstType.DIFF))
             else:
                if vi == 0:
                   matrix_term = matrix_term_v
@@ -789,7 +789,7 @@ def __absorb_constraints(self):
             elif term_constraint == ConstType.DROP:
                sterm.Z.append(Constraint(int(matrix_term.shape[1]/2),ConstType.DROP))
             elif term_constraint == ConstType.DIFF:
-                  sterm.Z.append(Constraint(None,ConstType.DIFF))
+               sterm.Z.append(Constraint(int(matrix_term.shape[1]/2),ConstType.DIFF))
 
     def build_penalties(self):
       """Builds the penalties required by ``self.__terms``. Called automatically whenever needed. Call manually only for testing."""
@@ -1401,8 +1401,9 @@ def build_smooth_term_matrix(ci,n_j,has_scale_split,sti,sterm,var_map,var_mins,v
             # Applies difference re-coding for sum-to-zero coefficients.
             # Based on smoothCon in mgcv(2017). See constraints.py
             # for more details.
-            matrix_term_v = np.diff(matrix_term_v)
-
+            matrix_term_v = np.diff(np.concatenate((matrix_term_v[:,sterm.Z[vi].Z:matrix_term_v.shape[1]],matrix_term_v[:,:sterm.Z[vi].Z]),axis=1))
+            matrix_term_v = np.concatenate((matrix_term_v[:,matrix_term_v.shape[1]-sterm.Z[vi].Z:],matrix_term_v[:,:matrix_term_v.shape[1]-sterm.Z[vi].Z]),axis=1)
+
       if vi == 0:
          matrix_term = matrix_term_v
       else:
@@ -1414,7 +1415,8 @@ def build_smooth_term_matrix(ci,n_j,has_scale_split,sti,sterm,var_map,var_mins,v
       elif sterm.Z[0].type == ConstType.DROP:
          matrix_term = np.delete(matrix_term,sterm.Z[0].Z,axis=1)
       elif sterm.Z[0].type == ConstType.DIFF:
-         matrix_term = np.diff(matrix_term)
+         matrix_term = np.diff(np.concatenate((matrix_term[:,sterm.Z[0].Z:matrix_term.shape[1]],matrix_term[:,:sterm.Z[0].Z]),axis=1))
+         matrix_term = np.concatenate((matrix_term[:,matrix_term.shape[1]-sterm.Z[0].Z:],matrix_term[:,:matrix_term.shape[1]-sterm.Z[0].Z]),axis=1)
 
    m_rows, m_cols = matrix_term.shape
    #print(m_cols)

src/mssm/src/python/penalties.py

@@ -57,7 +57,8 @@ def diff_pen(n,constraint,m=2):
         S = np.delete(np.delete(S,Z,axis=1),Z,axis=0)
         D = np.delete(D,Z,axis=0)
      elif constraint.type == ConstType.DIFF:
-        D = np.diff(D,axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017)
+        D = np.diff(np.concatenate((D[Z:D.shape[0],:],D[:Z,:]),axis=0),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017)
+        D = np.concatenate((D[D.shape[0]-Z:,:],D[:D.shape[0]-Z,:]),axis=0) 
         S = D @ D.T
 
   S = scp.sparse.csc_array(S)
@@ -113,7 +114,9 @@ def TP_pen(S_j,D_j,j,ks,constraint):
         S_TP = scp.sparse.csc_array(np.delete(np.delete(S_TP.toarray(),Z,axis=1),Z,axis=0))
         D_TP = scp.sparse.csc_array(np.delete(D_TP.toarray(),Z,axis=0))
      elif constraint.type == ConstType.DIFF:
-        D_TP = np.diff(D_TP.toarray(),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017)
+        D_TP = D_TP.toarray()
+        D_TP = np.diff(np.concatenate((D_TP[Z:D_TP.shape[0],:],D_TP[:Z,:]),axis=0),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017)
+        D_TP = np.concatenate((D_TP[D_TP.shape[0]-Z:,:],D_TP[:D_TP.shape[0]-Z,:]),axis=0) 
         S_TP = D_TP @ D_TP.T
         S_TP = scp.sparse.csc_array(S_TP)
         D_TP = scp.sparse.csc_array(D_TP)