Merge pull request #108 from palatej/develop

STL, linear filters

jdplus-incubator-base/jdplus-filters-base-parent/jdplus-filters-base-r/src/main/java/jdplus/filters/base/r/LinearFilters.java

@@ -0,0 +1,60 @@
+/*
+ * Copyright 2022 National Bank of Belgium
+ * 
+ * Licensed under the EUPL, Version 1.2 or – as soon they will be approved 
+ * by the European Commission - subsequent versions of the EUPL (the "Licence");
+ * You may not use this work except in compliance with the Licence.
+ * You may obtain a copy of the Licence at:
+ * 
+ * https://joinup.ec.europa.eu/software/page/eupl
+ * 
+ * Unless required by applicable law or agreed to in writing, software 
+ * distributed under the Licence is distributed on an "AS IS" basis,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the Licence for the specific language governing permissions and 
+ * limitations under the Licence.
+ */
+package jdplus.filters.base.r;
+
+import java.util.function.DoubleFunction;
+import jdplus.toolkit.base.api.math.Complex;
+import jdplus.toolkit.base.api.math.matrices.Matrix;
+import jdplus.toolkit.base.core.data.DataBlock;
+import jdplus.toolkit.base.core.math.linearfilters.FiniteFilter;
+import jdplus.toolkit.base.core.math.matrices.FastMatrix;
+
+/**
+ *
+ * @author palatej
+ */
+@lombok.experimental.UtilityClass
+public class LinearFilters {
+
+    /**
+     * 
+     * @param weights Weights of the filter, from its lower bound to the upper bound
+     * @param lbound Lower bound
+     * @param start First frequency 
+     * @param step Step between two successive frequencies
+     * @param n Number of steps
+     * @return A matrix with the used frequencies in the first column, the gain
+     * of the filter in the second column and the phase of the filter in the
+     * third column
+     * 
+     */
+    public Matrix frequencyResponse(double[] weights, int lbound, double start, double step, int n){
+        FiniteFilter ff=FiniteFilter.of(weights, lbound);
+        DoubleFunction<Complex> fn = ff.frequencyResponseFunction();
+        FastMatrix fr=FastMatrix.make(n, 3);
+        double f=start;
+        for (int i=0; i<n; ++i){
+            Complex cur=fn.apply(f);
+            DataBlock row = fr.row(i);
+            row.set(0, f);
+            row.set(1, cur.abs());
+            row.set(2, cur.arg());
+            f+=step;
+        }
+        return fr;
+    }
+}

jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/main/java/jdplus/stl/base/r/StlDecomposition.java

@@ -31,8 +31,8 @@
  */
 @lombok.experimental.UtilityClass
 public class StlDecomposition {
-    
-    public Matrix stl(double[] data, int period, boolean mul, int swindow, int twindow, int lwindow, int sjump, int tjump, int ljump, int nin, int nout, double weightThreshold, String weightsFunction, boolean legacy) {
+
+    public Matrix stl(double[] data, int period, boolean mul, int swindow, int twindow, int lwindow, int sdegree, int tdegree, int ldegree, int sjump, int tjump, int ljump, int nin, int nout, double weightThreshold, String weightsFunction, boolean legacy) {
         if (nin < 1) {
             nin = 1;
         }
@@ -45,32 +45,44 @@ public Matrix stl(double[] data, int period, boolean mul, int swindow, int twind
         if (twindow == 0) {
             twindow = LoessSpec.defaultTrendWindow(period, swindow);
         }
-        if (lwindow == 0){
-            lwindow=period%2 == 0 ? period+1 : period;
+        if (lwindow == 0) {
+            lwindow = period % 2 == 0 ? period + 1 : period;
+        }
+        if (legacy) {
+            return stl_legacy(data, period, mul, swindow, twindow, lwindow, sdegree, tdegree, ldegree, sjump, tjump, ljump, nin, nout, weightThreshold, weightsFunction);
+        } else {
+            return stl_new(data, period, mul, swindow, twindow, lwindow, sdegree, tdegree, ldegree, sjump, tjump, ljump, nin, nout, weightThreshold, weightsFunction);
         }
+    }
+
+    private Matrix stl_legacy(double[] data, int period, boolean mul, int swindow, int twindow, int lwindow, int sdegree, int tdegree, int ldegree, int sjump, int tjump, int ljump, int nin, int nout, double weightThreshold, String weightsFunction) {
+
         StlLegacySpec spec = StlLegacySpec.defaultSpec(period, swindow, true);
         spec.setMultiplicative(mul);
-        spec.setLegacy(legacy);
+        spec.setLegacy(true);
         spec.setNi(nin);
         spec.setNo(nout);
         spec.setNs(swindow);
+        spec.setSdeg(sdegree);
         spec.setNsjump(sjump);
         spec.setNt(twindow);
+        spec.setTdeg(tdegree);
         spec.setNtjump(tjump);
         spec.setNl(lwindow);
+        spec.setLdeg(ldegree);
         spec.setNljump(ljump);
         spec.setWthreshold(weightThreshold);
         spec.setWfn(WeightFunction.valueOf(weightsFunction).asFunction());
         StlLegacy stl = new StlLegacy(spec);
         DoubleSeq y = DoubleSeq.of(data).cleanExtremities();
-        
+
         int n = y.length();
         if (!stl.process(y)) {
             return null;
         }
-        
+
         FastMatrix M = FastMatrix.make(n, 7);
-        
+
         M.column(0).copy(y);
         M.column(2).copyFrom(stl.getTrend(), 0);
         M.column(3).copyFrom(stl.getSeason(), 0);
@@ -90,7 +102,56 @@ public Matrix stl(double[] data, int period, boolean mul, int swindow, int twind
         }
         return M;
     }
-
+
+    private Matrix stl_new(double[] data, int period, boolean mul, int swindow, int twindow, int lwindow, int sdegree, int tdegree, int ldegree, int sjump, int tjump, int ljump, int nin, int nout, double weightThreshold, String weightsFunction) {
+        LoessSpec tspec = LoessSpec.builder()
+                .window(twindow)
+                .jump(tjump)
+                .degree(tdegree)
+                .build();
+        LoessSpec sspec = LoessSpec.builder()
+                .window(swindow)
+                .jump(sjump)
+                .degree(sdegree)
+                .build();
+        LoessSpec lspec = LoessSpec.builder()
+                .window(lwindow)
+                .jump(ljump)
+                .degree(ldegree)
+                .build();
+
+        StlSpec spec = StlSpec.builder()
+                .multiplicative(mul)
+                .trendSpec(tspec)
+                .seasonalSpec(SeasonalSpec.builder()
+                        .period(period)
+                        .seasonalSpec(sspec)
+                        .lowPassSpec(lspec)
+                        .build())
+                .innerLoopsCount(nin)
+                .outerLoopsCount(nout)
+                .robustWeightThreshold(weightThreshold)
+                .robustWeightFunction(WeightFunction.valueOf(weightsFunction))
+                .build();
+
+        RawStlKernel stl = new RawStlKernel(spec);
+        DoubleSeq y = DoubleSeq.of(data).cleanExtremities();
+
+        int n = y.length();
+        RawStlResults rslt = stl.process(DoubleSeq.of(data));
+
+        FastMatrix M = FastMatrix.make(n, 7);
+
+        M.column(0).copy(y);
+        M.column(2).copy(rslt.getTrend());
+        M.column(3).copy(rslt.getSeasonal());
+        M.column(4).copy(rslt.getIrregular());
+        M.column(5).copy(rslt.getFit());
+        M.column(6).copy(rslt.getWeights());
+        M.column(1).copy(rslt.getSa());
+        return M;
+    }
+
     private int max(int[] v) {
         int m = v[0];
         for (int i = 1; i < v.length; ++i) {
@@ -100,23 +161,23 @@ private int max(int[] v) {
         }
         return m;
     }
-    
+
     public Matrix mstl(double[] data, int[] periods, boolean mul, int[] swindow, int twindow, int nin, int nout, boolean nojump, double weightThreshold, String weightsFunction) {
         if (periods == null || (swindow != null && periods.length != swindow.length)) {
             return null;
         }
         if (twindow == 0) {
             twindow = LoessSpec.defaultTrendWindow(max(periods));
         }
-        
+
         MStlSpec.Builder builder = MStlSpec.builder()
                 .innerLoopsCount(nin)
                 .outerLoopsCount(nout)
                 .multiplicative(mul)
                 .trendSpec(LoessSpec.of(twindow, 1, nojump))
                 .robustWeightThreshold(weightThreshold)
                 .robustWeightFunction(WeightFunction.valueOf(weightsFunction));
-        
+
         if (swindow == null) {
             for (int i = 0; i < periods.length; ++i) {
                 builder.seasonalSpec(SeasonalSpec.createDefault(periods[i], nojump));
@@ -129,17 +190,17 @@ public Matrix mstl(double[] data, int[] periods, boolean mul, int[] swindow, int
             for (int i = 0; i < periods.length; ++i) {
                 builder.seasonalSpec(new SeasonalSpec(periods[i], swindow[i], nojump));
             }
-            
+
         }
         MStlSpec spec = builder.build();
-        
+
         MStlKernel stl = MStlKernel.of(spec);
         DoubleSeq y = DoubleSeq.of(data).cleanExtremities();
         stl.process(y);
-        
+
         int n = y.length();
         FastMatrix M = FastMatrix.make(n, 6 + periods.length);
-        
+
         M.column(0).copyFrom(stl.getY(), 0);
         M.column(1).copy(M.column(0));
         M.column(2).copyFrom(stl.getTrend(), 0);
@@ -157,22 +218,22 @@ public Matrix mstl(double[] data, int[] periods, boolean mul, int[] swindow, int
         M.column(j).copyFrom(stl.getWeights(), 0);
         return M;
     }
-    
+
     public Matrix istl(double[] data, int[] periods, boolean mul, int[] swindow, int[] twindow, int nin, int nout, boolean nojump, double weightThreshold, String weightsFunction) {
         if (periods == null || (swindow != null && periods.length != swindow.length)) {
             return null;
         }
         if (twindow != null && twindow.length != periods.length) {
             return null;
         }
-        
+
         IStlSpec.Builder builder = IStlSpec.builder()
                 .innerLoopsCount(nin)
                 .outerLoopsCount(nout)
                 .multiplicative(mul)
                 .robustWeightThreshold(weightThreshold)
                 .robustWeightFunction(WeightFunction.valueOf(weightsFunction));
-        
+
         for (int i = 0; i < periods.length; ++i) {
             SeasonalSpec sspec;
             if (swindow == null) {
@@ -189,16 +250,16 @@ public Matrix istl(double[] data, int[] periods, boolean mul, int[] swindow, int
                 tspec = LoessSpec.of(twindow[i], 1, nojump);
             }
             builder.periodSpec(new IStlSpec.PeriodSpec(tspec, sspec));
-            
+
         }
         IStlSpec spec = builder.build();
-        
+
         DoubleSeq y = DoubleSeq.of(data).cleanExtremities();
         MStlResults rslt = IStlKernel.process(y, spec);
-        
+
         int n = y.length();
         FastMatrix M = FastMatrix.make(n, 6 + periods.length);
-        
+
         M.column(0).copy(y);
         M.column(1).copy(M.column(0));
         M.column(2).copy(rslt.getTrend());
@@ -216,13 +277,13 @@ public Matrix istl(double[] data, int[] periods, boolean mul, int[] swindow, int
         M.column(++j).copy(rslt.getWeights());
         return M;
     }
-    
+
     public double[] loess(double[] y, int window, int degree, int jump) {
         LoessSpec spec = LoessSpec.of(window, degree, jump, null);
         LoessFilter filter = new LoessFilter(spec);
         double[] z = new double[y.length];
         filter.filter(IDataGetter.of(y), null, IDataSelector.of(z));
         return z;
     }
-    
+
 }

jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/test/java/jdplus/stl/base/r/StlDecompositionTest.java

@@ -28,7 +28,14 @@ public void testFilter() {
 
     @Test
     public void testStl() {
-        Matrix decomp = StlDecomposition.stl(Data.ABS_RETAIL, 12, false, 0, 0, 0, 0, 0, 0, 2, 5, 0.1, WeightFunction.TRICUBE.name(), false);
+        Matrix decomp = StlDecomposition.stl(Data.ABS_RETAIL, 12, false, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 5, 0.1, WeightFunction.TRICUBE.name(), false);
+//        System.out.println(decomp);
+        assertTrue(null != decomp);
+    }
+
+    @Test
+    public void testStl_legacy() {
+        Matrix decomp = StlDecomposition.stl(Data.ABS_RETAIL, 12, false, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 5, 0.1, WeightFunction.TRICUBE.name(), true);
 //        System.out.println(decomp);
         assertTrue(null != decomp);
     }

pom.xml

@@ -59,7 +59,7 @@
 
         <!-- dependencies -->
         <jd2.version>2.2.4</jd2.version>
-        <jdplus-main.version>3.2.1</jdplus-main.version>
+        <jdplus-main.version>3.2.2-SNAPSHOT</jdplus-main.version>
         <protobuf.version>3.21.12</protobuf.version>
     </properties>