added autoapproximate functionality for cos basis

Author: passscoed

src/pyANOVAapprox/__init__.py

@@ -8,6 +8,8 @@
 from scipy.optimize import bisect
 from scipy.sparse.linalg import lsqr
 from scipy.special import erf
+import csv
+import os
 
 # from sklearn.metrics import roc_auc_score
 

src/pyANOVAapprox/approx.py

@@ -187,6 +187,9 @@ def __init__(
     ):
 
         self.X = X
+        
+        if X.shape[0] != y.shape[0]:
+            raise ValueError("X and y have different lengths.")
 
         setting = approx_setting(
             parent=self,
@@ -384,7 +387,7 @@ def _autoapproximate(
         B,
         maxiter,
         solver,
-        verbose,
+        verbosity,
         solver_max_iter,
         solver_weights,
         solver_verbose,
@@ -398,8 +401,16 @@ def _autoapproximate(
 
         D = dict([(u, tuple([1.0] * len(u))) for u in setting.U])
         t = dict([(u, tuple([1.0] * len(u))) for u in setting.U])
+        
+        if verbosity>3:
+            if not os.path.exists("log"):
+                os.mkdir("log")
+            with open('log/log.csv', 'w', newline='') as csvfile:
+                csv.writer(csvfile, delimiter=',').writerow(["it"] + setting.U)
 
         for idx in range(maxiter):
+            if verbosity > 0:
+                print("===== Iteration ", str(idx + 1), " =====")
             bw = compute_bandwidth(B, D, t)
             if setting.N is not None:
                 self.addSetting(setting)
@@ -409,9 +420,14 @@ def _autoapproximate(
             setting.N = [bw[i] for i in setting.U]
             self.addTrafo()
 
-            if verbose:
-                print("bw in iteration", str(idx + 1), "are", str(bw))
-
+            if verbosity > 0:
+                for i in setting.U:
+                    print("bw in", str(i), ":", bw[i])
+                #print("bw in iteration", str(idx + 1), "are", str(bw))
+                print()
+            if verbosity > 3:
+                with open('log/log.csv', 'a', newline='') as csvfile:
+                    csv.writer(csvfile, delimiter=',').writerow(["bw in it"+str(idx+1)] + [str(bw[i]) for i in setting.U])
             self.approximate(
                 lam=lam,
                 solver=solver,
@@ -421,16 +437,25 @@ def _autoapproximate(
                 tol=solver_tol,
             )
 
-            D, t = self.estimate_rates(lam=lam, verbose=verbose)
-            if verbose:
-                print(
-                    "estimated rates in iteration",
-                    str(idx + 1),
-                    "are D =",
-                    str(D),
-                    "and t =",
-                    str(t),
-                )
+            D, t = self.estimate_rates(lam=lam, verbosity=verbosity)
+            if verbosity > 1:
+                for i in setting.U:
+                    print("estimated rates for", str(i), ": D = ", str(D[i]), "and t = ", str(t[i]) )
+                print()
+                
+                #print(
+                #    "estimated rates in iteration",
+                #    str(idx + 1),
+                #    "are D =",
+                #    str(D),
+                #    "and t =",
+                #    str(t),
+                #)
+            if verbosity > 3:
+                with open('log/log.csv', 'a', newline='') as csvfile:
+                    wr = csv.writer(csvfile, delimiter=',')
+                    wr.writerow(["D in it"+str(idx+1)] + [str(D[i]) for i in setting.U])
+                    wr.writerow(["t in it"+str(idx+1)] + [str(t[i]) for i in setting.U])
         return D, t
 
     def autoapproximate(
@@ -440,7 +465,7 @@ def autoapproximate(
         B=None,
         maxiter=2,
         solver="lsqr",
-        verbose=False,
+        verbosity=0,
         solver_max_iter=50,
         solver_weights=None,
         solver_verbose=False,
@@ -462,7 +487,7 @@ def autoapproximate(
                 B=B,
                 maxiter=maxiter,
                 solver=solver,
-                verbose=verbose,
+                verbosity=verbosity,
                 solver_max_iter=solver_max_iter,
                 solver_weights=solver_weights,
                 solver_verbose=solver_verbose,
@@ -476,7 +501,7 @@ def autoapproximate(
                     B=B,
                     maxiter=maxiter,
                     solver=solver,
-                    verbose=verbose,
+                    verbosity=verbosity,
                     solver_max_iter=solver_max_iter,
                     solver_weights=solver_weights,
                     solver_verbose=solver_verbose,

src/pyANOVAapprox/bandwidth.py

@@ -13,17 +13,22 @@ def getfcu(ghat, u):
     return fcu
 
 
-def getaxissum(ghat, u, j):
+def getaxissum(ghat, u, j, system):
     fcu = getfcu(ghat, u)
     idx = [s.u for s in ghat.settings].index(u)
     bws = ghat.settings[idx].bandwidths
 
     fcuj = np.sum(abs(fcu) ** 2, axis=tuple([i for i in range(len(u)) if i != j]))
 
-    fcuj = (
-        fcuj[np.mod(range(int(bws[j] / 2), bws[j]), bws[j] - 1)]
-        + fcuj[range(int(bws[j] / 2) - 1, -1, -1)]
-    )
+    if system == "exp":
+        fcuj = (
+            fcuj[np.mod(range(int(bws[j] / 2), bws[j]), bws[j] - 1)]
+            + fcuj[range(int(bws[j] / 2) - 1, -1, -1)]
+        )
+    elif system == "cos":
+        pass
+    else:
+        raise ValueError("For this basis is estimate rates not implemented")
     # fcuj = np.concatenate(fcuj[math.ceil(bws[j]/2):] + [fcuj[0]]) + fcuj[math.ceil(bws[j]/2)-1::-1]
     return fcuj
 
@@ -106,7 +111,7 @@ def most_common_value(data):
     return t[max_index]
 
 
-def estimate_rates(self, lam, settingnr=None, verbose=False):
+def estimate_rates(self, lam, settingnr=None, verbosity=0):
     us = [s.u for s in self.getTrafo(settingnr).settings]
     nhat = GroupedCoefficients(
         self.getTrafo(settingnr).settings,
@@ -118,25 +123,28 @@ def estimate_rates(self, lam, settingnr=None, verbose=False):
 
     mcl = np.exp(most_common_value(np.log(abs(self.getFc(settingnr)[lam].data))))
     threshold = 100 * mcl**2
-
-    if verbose:
-        ps = []
+    
+    if verbosity>5:
+        if not os.path.exists(os.path.join("log","figures")):
+            os.mkdir(os.path.join("log","figures"))
+        num = 0
+        
+    system = self.getTrafo(settingnr).system
 
     for u in us:
         if len(u) == 0:
             continue
 
-        if verbose:
+        if verbosity>5:
             fig, ax = plt.subplots()
             ax.set_xscale("log")
             ax.set_yscale("log")
             ax.set_title(str(u))
-            ps.append(ax)
 
         for j in range(len(u)):
-            axissum = getaxissum(self.getFc(settingnr)[lam], u, j)
+            axissum = getaxissum(self.getFc(settingnr)[lam], u, j, system)
 
-            axissumnum = getaxissum(nhat, u, j)
+            axissumnum = getaxissum(nhat, u, j, system)
             idx = next(
                 (
                     i
@@ -146,7 +154,7 @@ def estimate_rates(self, lam, settingnr=None, verbose=False):
                 None,
             )
 
-            if verbose:
+            if verbosity>5:
 
                 y = np.cumsum(axissum[::-1])[::-1]
                 ax.plot(
@@ -166,18 +174,20 @@ def estimate_rates(self, lam, settingnr=None, verbose=False):
                 D[u][j] = Duj
                 t[u][j] = -tuj / 2
 
-                if verbose:
+                if verbosity>5:
                     x = np.arange(1, idx + 1)
                     ax.plot(
                         x,
                         D[u][j] * x ** (-2 * t[u][j]),
                         linewidth=2,
                         #                        color=j
                     )
-    if verbose:
-        plt.figure(figsize=(12, 9))
-        for ax in ps:
-            plt.sca(ax)
-        plt.show()
+        if verbosity>5:
+            #plt.figure(figsize=(12, 9))
+            #for ax in ps:
+            #    plt.sca(ax)
+            fig.savefig(os.path.join("log","figures", str(num).strip()+ "_rates_" + str(u).strip() + ".png"))
+            num = num + 1
+            plt.close(fig) 
 
     return D, t