diff --git a/jdplus-incubator-base/jdplus-filters-base-parent/jdplus-filters-base-core/src/test/java/jdplus/filters/base/core/FSTFilterTest.java b/jdplus-incubator-base/jdplus-filters-base-parent/jdplus-filters-base-core/src/test/java/jdplus/filters/base/core/FSTFilterTest.java index 94257b5..ebce51d 100644 --- a/jdplus-incubator-base/jdplus-filters-base-parent/jdplus-filters-base-core/src/test/java/jdplus/filters/base/core/FSTFilterTest.java +++ b/jdplus-incubator-base/jdplus-filters-base-parent/jdplus-filters-base-core/src/test/java/jdplus/filters/base/core/FSTFilterTest.java @@ -33,7 +33,7 @@ public void testSymmetric() { .nleads(4) .build(); FSTFilter.Results rslt = ff.make(.1 * i, 0, false); - System.out.println(DoubleSeq.of(rslt.getFilter().weightsToArray())); + // System.out.println(DoubleSeq.of(rslt.getFilter().weightsToArray())); } } @@ -45,11 +45,11 @@ public void testASymmetric() { .nleads(i) .build(); FSTFilter.Results rslt = ff.make(1, 0, true); - System.out.print(rslt.getS()); - System.out.print('\t'); - System.out.print(rslt.getF()); - System.out.print('\t'); - System.out.println(DoubleSeq.of(rslt.getFilter().weightsToArray())); +// System.out.print(rslt.getS()); +// System.out.print('\t'); +// System.out.print(rslt.getF()); +// System.out.print('\t'); +// System.out.println(DoubleSeq.of(rslt.getFilter().weightsToArray())); } } diff --git a/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-core/src/test/java/jdplus/highfreq/base/core/extendedairline/ExtendedAirlineMappingTest.java b/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-core/src/test/java/jdplus/highfreq/base/core/extendedairline/ExtendedAirlineMappingTest.java index 05eaa6e..6fbd1ce 100644 --- a/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-core/src/test/java/jdplus/highfreq/base/core/extendedairline/ExtendedAirlineMappingTest.java +++ b/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-core/src/test/java/jdplus/highfreq/base/core/extendedairline/ExtendedAirlineMappingTest.java @@ -44,7 +44,7 @@ public void testAR_int() { ExtendedAirlineMapping mapping = ExtendedAirlineMapping.of(spec); DoubleSeq p = mapping.getDefaultParameters(); ArimaModel m = mapping.map(p); - System.out.println(m); +// System.out.println(m); DoubleSeq np = mapping.parametersOf(m); assertEquals(p, np); } @@ -61,7 +61,7 @@ public void testAR_noint() { ExtendedAirlineMapping mapping = ExtendedAirlineMapping.of(spec); DoubleSeq p = mapping.getDefaultParameters(); ArimaModel m = mapping.map(p); - System.out.println(m); +// System.out.println(m); DoubleSeq np = mapping.parametersOf(m); assertTrue(p.distance(np) < 1.e-9); } @@ -78,7 +78,7 @@ public void test_int() { ExtendedAirlineMapping mapping = ExtendedAirlineMapping.of(spec); DoubleSeq p = mapping.getDefaultParameters(); ArimaModel m = mapping.map(p); - System.out.println(m); +// System.out.println(m); DoubleSeq np = mapping.parametersOf(m); assertEquals(p, np); } @@ -95,7 +95,7 @@ public void test_noint() { ExtendedAirlineMapping mapping = ExtendedAirlineMapping.of(spec); DoubleSeq p = mapping.getDefaultParameters(); ArimaModel m = mapping.map(p); - System.out.println(m); +// System.out.println(m); DoubleSeq np = mapping.parametersOf(m); assertTrue(p.distance(np) < 1.e-9); } diff --git a/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-r/src/test/java/jdplus/highfreq/base/r/SplinesTest.java b/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-r/src/test/java/jdplus/highfreq/base/r/SplinesTest.java index 639ef72..8bf1790 100644 --- a/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-r/src/test/java/jdplus/highfreq/base/r/SplinesTest.java +++ b/jdplus-incubator-base/jdplus-highfreq-base-parent/jdplus-highfreq-base-r/src/test/java/jdplus/highfreq/base/r/SplinesTest.java @@ -43,6 +43,7 @@ import jdplus.toolkit.base.core.ssf.ISsfLoading; import jdplus.toolkit.base.core.ssf.StateStorage; import jdplus.toolkit.base.core.timeseries.calendars.HolidaysUtility; +import tck.demetra.data.Data; /** * @@ -55,7 +56,7 @@ public class SplinesTest { static { DoubleSeq y; try { - InputStream stream = ExtendedAirlineMapping.class.getResourceAsStream("/births.txt"); + InputStream stream = Data.class.getResourceAsStream("/births.txt"); Matrix edf = MatrixSerializer.read(stream); y = edf.column(0); } catch (IOException ex) { @@ -101,7 +102,7 @@ public static void main(String[] args) { //StateItem l = AtomicModels.localLevel("l", .01, false, Double.NaN); StateItem l = AtomicModels.localLinearTrend("l", .01, .01, false, false); StateItem sw = AtomicModels.seasonalComponent("sw", "HarrisonStevens", 7, .01, false); - StateItem sy = AtomicModels.dailySplineComponent("sy", 1968, 2024, pos, 0, .01, false); + StateItem sy = AtomicModels.dailySplines("sy", 1968, pos, 0, .01, false); StateItem reg = AtomicModels.timeVaryingRegression("reg", X, 0.01, false); StateItem n = AtomicModels.noise("n", .01, false); ModelEquation eq = new ModelEquation("eq1", 0, true); @@ -149,7 +150,7 @@ public static void main(String[] args) { public static void main2(String[] args) { DoubleSeq y; try { - InputStream stream = ExtendedAirlineMapping.class.getResourceAsStream("/usclaims.txt"); + InputStream stream = Data.class.getResourceAsStream("/usclaims.txt"); Matrix us = MatrixSerializer.read(stream); y = us.column(0); } catch (IOException ex) { @@ -168,7 +169,7 @@ public static void main2(String[] args) { // StateItem l = AtomicModels.localLevel("l", .01, false, Double.NaN); StateItem l = AtomicModels.localLinearTrend("l", .01, .01, false, false); // StateItem sw = AtomicModels.seasonalComponent("sw", "HarrisonStevens", 7, .01, false); - StateItem sy = AtomicModels.weeklySplineComponent("sy", 1967, 2025, pos, 6, 0, .01, false); + StateItem sy = AtomicModels.regularSplines("sy", 1967, new double[]{1,5,10,15,20,30,40,50}, 0, .01, false); // StateItem reg=AtomicModels.timeVaryingRegression("reg", X, 0.01, false); StateItem n = AtomicModels.noise("n", .01, false); ModelEquation eq = new ModelEquation("eq1", 0, true); diff --git a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-core/src/test/java/jdplus/stl/base/core/StlLegacyTest.java b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-core/src/test/java/jdplus/stl/base/core/StlLegacyTest.java index 40e4ca0..7abf51c 100644 --- a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-core/src/test/java/jdplus/stl/base/core/StlLegacyTest.java +++ b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-core/src/test/java/jdplus/stl/base/core/StlLegacyTest.java @@ -66,7 +66,7 @@ public void testDefault() { spec.setNs(13); spec.setNt(23); spec.setNi(2); - spec.setNo(1); + spec.setNo(0); spec.setNsjump(0); spec.setNtjump(0); spec.setNljump(0); diff --git a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/main/java/jdplus/stl/base/r/StlDecomposition.java b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/main/java/jdplus/stl/base/r/StlDecomposition.java index d538698..95f03e2 100644 --- a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/main/java/jdplus/stl/base/r/StlDecomposition.java +++ b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/main/java/jdplus/stl/base/r/StlDecomposition.java @@ -31,7 +31,7 @@ */ @lombok.experimental.UtilityClass public class StlDecomposition { - + 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; @@ -54,9 +54,9 @@ public Matrix stl(double[] data, int period, boolean mul, int swindow, int twind 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(true); @@ -75,14 +75,14 @@ private Matrix stl_legacy(double[] data, int period, boolean mul, int swindow, i 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); @@ -102,7 +102,7 @@ private Matrix stl_legacy(double[] data, int period, boolean mul, int swindow, i } 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) @@ -119,7 +119,7 @@ private Matrix stl_new(double[] data, int period, boolean mul, int swindow, int .jump(ljump) .degree(ldegree) .build(); - + StlSpec spec = StlSpec.builder() .multiplicative(mul) .trendSpec(tspec) @@ -133,25 +133,29 @@ private Matrix stl_new(double[] data, int period, boolean mul, int swindow, int .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()); + if (!rslt.getWeights().isEmpty()) { + M.column(6).copy(rslt.getWeights()); + } else { + M.column(6).set(1); + } 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) { @@ -161,7 +165,7 @@ 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; @@ -169,7 +173,7 @@ public Matrix mstl(double[] data, int[] periods, boolean mul, int[] swindow, int if (twindow == 0) { twindow = LoessSpec.defaultTrendWindow(max(periods)); } - + MStlSpec.Builder builder = MStlSpec.builder() .innerLoopsCount(nin) .outerLoopsCount(nout) @@ -177,7 +181,7 @@ public Matrix mstl(double[] data, int[] periods, boolean mul, int[] swindow, int .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)); @@ -190,17 +194,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); @@ -218,7 +222,7 @@ 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; @@ -226,14 +230,14 @@ public Matrix istl(double[] data, int[] periods, boolean mul, int[] swindow, int 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) { @@ -250,16 +254,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()); @@ -277,7 +281,7 @@ 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); @@ -285,5 +289,5 @@ public double[] loess(double[] y, int window, int degree, int jump) { filter.filter(IDataGetter.of(y), null, IDataSelector.of(z)); return z; } - + } diff --git a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/test/java/jdplus/stl/base/r/StlDecompositionTest.java b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/test/java/jdplus/stl/base/r/StlDecompositionTest.java index 7c18cf2..936794d 100644 --- a/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/test/java/jdplus/stl/base/r/StlDecompositionTest.java +++ b/jdplus-incubator-base/jdplus-stl-base-parent/jdplus-stl-base-r/src/test/java/jdplus/stl/base/r/StlDecompositionTest.java @@ -28,14 +28,14 @@ public void testFilter() { @Test public void testStl() { - 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); + Matrix decomp = StlDecomposition.stl(Data.ABS_RETAIL, 12, false, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 0, 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); + Matrix decomp = StlDecomposition.stl(Data.ABS_RETAIL, 12, false, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 0, 0.1, WeightFunction.TRICUBE.name(), true); // System.out.println(decomp); assertTrue(null != decomp); } diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/AtomicModels.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/AtomicModels.java index 7d01760..c0b7675 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/AtomicModels.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/AtomicModels.java @@ -16,6 +16,9 @@ */ package jdplus.sts.base.core.msts; +import jdplus.sts.base.core.splines.DailySpline; +import jdplus.sts.base.core.splines.RegularSpline; +import jdplus.sts.base.core.splines.SplineData; import jdplus.sts.base.core.msts.internal.ArItem; import jdplus.sts.base.core.msts.internal.ArItem2; import jdplus.sts.base.core.msts.internal.ArimaItem; @@ -39,12 +42,9 @@ import jdplus.sts.base.core.msts.internal.VarSeasonalComponentItem; import jdplus.sts.base.core.msts.internal.VarNoiseItem; import jdplus.toolkit.base.api.math.matrices.Matrix; -import jdplus.sts.base.core.msts.internal.RegularSplineItem; import jdplus.sts.base.core.msts.internal.SplineItem; import jdplus.sts.base.core.msts.internal.VarRegressionItem; -import jdplus.toolkit.base.core.ssf.sts.splines.DailySpline; -import jdplus.toolkit.base.core.ssf.sts.splines.SplineData; -import jdplus.toolkit.base.core.ssf.sts.splines.WeeklySpline; +import jdplus.toolkit.base.api.data.DoubleSeq; /** * @@ -149,15 +149,27 @@ public StateItem periodicComponent(String name, double period, int[] k, double c return new PeriodicItem(name, period, k, cvar, fixedvar); } - public StateItem regularSplineComponent(String name, int[] pos, int startpos, double cvar, boolean fixedvar) { - return new RegularSplineItem(name, pos, startpos, cvar, fixedvar); + public StateItem regularSplines(String name, double period, int startpos, double cvar, boolean fixedvar) { + RegularSpline rs = RegularSpline.of(period); + SplineData sd = new SplineData(rs); + return new SplineItem(name, sd, startpos, cvar, fixedvar); } - - public StateItem dailySplineComponent(String name, int startYear, int endYear, int[] pos, int startpos, double cvar, boolean fixedvar) { - return new SplineItem(name, SplineData.of(new DailySpline(startYear, pos), endYear-startYear+1), startpos, cvar, fixedvar); + + public StateItem regularSplines(String name, double period, int nnodes, int startpos, double cvar, boolean fixedvar) { + RegularSpline rs = RegularSpline.of(period, nnodes); + SplineData sd = new SplineData(rs); + return new SplineItem(name, sd, startpos, cvar, fixedvar); + } + + public StateItem regularSplines(String name, double period, double[] nodes, int startpos, double cvar, boolean fixedvar) { + RegularSpline rs = RegularSpline.of(period, DoubleSeq.of(nodes)); + SplineData sd = new SplineData(rs); + return new SplineItem(name, sd, startpos, cvar, fixedvar); } - public StateItem weeklySplineComponent(String name, int startYear, int endYear, int[] pos, int startPos, int modelStart, double cvar, boolean fixedvar) { - return new SplineItem(name, SplineData.of(new WeeklySpline(startYear, startPos, pos), endYear-startYear+1), modelStart, cvar, fixedvar); + public StateItem dailySplines(String name, int startYear, int[] pos, int startpos, double cvar, boolean fixedvar) { + DailySpline rs = new DailySpline(startYear, pos); + SplineData sd = new SplineData(rs); + return new SplineItem(name, sd, startpos, cvar, fixedvar); } } diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/CompositeModelEstimation.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/CompositeModelEstimation.java index ef31f68..e8c2149 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/CompositeModelEstimation.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/CompositeModelEstimation.java @@ -60,6 +60,7 @@ public static CompositeModelEstimation estimationOf(CompositeModel model, FastMa .build(); MstsMapping mapping = model.mapping(); monitor.process(data, mapping, parameters == null ? null : DoubleSeq.of(parameters)); + rslt.model = model; rslt.likelihood = monitor.getLikelihood(); rslt.ssf = monitor.getSsf(); rslt.cmpPos = rslt.getSsf().componentsPosition(); @@ -72,6 +73,7 @@ public static CompositeModelEstimation estimationOf(CompositeModel model, FastMa public static CompositeModelEstimation computationOf(CompositeModel model, FastMatrix data, DoubleSeq fullParameters, boolean marginal, boolean concentrated) { CompositeModelEstimation rslt = new CompositeModelEstimation(); + rslt.model = model; rslt.data = data; rslt.fullParameters = fullParameters.toArray(); MstsMapping mapping = model.mapping(); @@ -97,6 +99,7 @@ public static CompositeModelEstimation computationOf(CompositeModel model, FastM private double[] fullParameters, parameters; private String[] parametersName, cmpName; private StateStorage smoothedStates, filteredStates, filteringStates; + private CompositeModel model; public StateStorage getSmoothedStates() { if (smoothedStates == null) { @@ -230,6 +233,76 @@ public StateStorage getFilteringStates() { return filteringStates; } + private static int find(String[] ss, String s) { + for (int i = 0; i < ss.length; ++i) { + if (ss[i].equalsIgnoreCase(s)) { + return i; + } + } + return -1; + } + + public FastMatrix getSmoothedComponents(int p) { + ModelEquation equation = model.getEquation(p); + StateStorage ss = getSmoothedStates(); + int nr = ss.size(); + int nc = equation.getItemsCount(); + FastMatrix C = FastMatrix.make(nr, nc); + int[] cmpDim = ssf.componentsDimension(); + + for (int i = 0; i < nc; ++i) { + ModelEquation.Item item = equation.getItem(i); + String cmp = item.getCmp(); + int pos = find(cmpName, cmp); + if (pos >= 0) { + ISsfLoading loading = item.getLoading(); + int start = cmpPos[pos], end = start + cmpDim[pos]; + DataBlock col = C.column(i); + for (int j = 0; j < nr; ++j) { + col.set(j, loading.ZX(j, ss.a(j).range(start, end))); + } + LoadingInterpreter li = item.getC(); + if (li != null) { + col.mul(li.value()); + } + } + } + return C; + } + + public FastMatrix getSmoothedComponentVariance(int eq) { + ModelEquation equation = model.getEquation(eq); + StateStorage ss = getSmoothedStates(); + int nr = ss.size(); + int nc = equation.getItemsCount(); + FastMatrix C = FastMatrix.make(nr, nc); + int[] cmpDim = ssf.componentsDimension(); + + for (int i = 0; i < nc; ++i) { + ModelEquation.Item item = equation.getItem(i); + String cmp = item.getCmp(); + int pos = find(cmpName, cmp); + if (pos >= 0) { + ISsfLoading loading = item.getLoading(); + int start = cmpPos[pos]; + DataBlock col = C.column(i); + for (int j = 0; j < nr; ++j) { + double v = loading.ZVZ(j, ss.P(j).extract(start, cmpDim[pos], start, cmpDim[pos])); + col.set(j, v < 0 ? 0 : v); + } + LoadingInterpreter li = item.getC(); + if (li != null) { + col.mul(li.value()*li.value()); + } + } + } + return C; + } + + public CompositeModel getModel() { + return model; + } + public DoubleSeq signal(int obs, int[] cmps) { if (obs >= data.getColumnsCount()) { return null; diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/RegularSplineItem.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/RegularSplineItem.java deleted file mode 100644 index 88a696b..0000000 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/RegularSplineItem.java +++ /dev/null @@ -1,108 +0,0 @@ -/* - * Copyright 2023 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.sts.base.core.msts.internal; - -import jdplus.sts.base.core.msts.StateItem; -import jdplus.toolkit.base.api.data.DoubleSeq; -import jdplus.sts.base.core.msts.MstsMapping; -import jdplus.sts.base.core.msts.VarianceInterpreter; -import java.util.Collections; -import java.util.List; -import jdplus.sts.base.core.msts.ParameterInterpreter; -import jdplus.toolkit.base.core.ssf.ISsfLoading; -import jdplus.toolkit.base.core.ssf.StateComponent; -import jdplus.toolkit.base.core.ssf.sts.RegularSplineComponent; - -/** - * - * @author palatej - */ -public class RegularSplineItem extends StateItem { - - private final VarianceInterpreter v; - private final int startpos; - private final RegularSplineComponent.Data data; - - public RegularSplineItem(String name, int[] pos, int startpos, double cvar, boolean fixedvar) { - super(name); - v = new VarianceInterpreter(name + ".var", cvar, fixedvar, true); - this.data = RegularSplineComponent.Data.of(pos); - this.startpos = startpos; - } - - private RegularSplineItem(RegularSplineItem item) { - super(item.name); - this.data = item.data; - this.startpos = item.startpos; - v = item.v.duplicate(); - } - - @Override - public RegularSplineItem duplicate() { - return new RegularSplineItem(this); - } - - @Override - public void addTo(MstsMapping mapping) { - mapping.add(v); - mapping.add((p, builder) -> { - double var = p.get(0); - builder.add(name, RegularSplineComponent.stateComponent(data, var), RegularSplineComponent.loading(data, startpos)); - return 1; - }); - } - - @Override - public List parameters() { - return Collections.singletonList(v); - } - - @Override - public StateComponent build(DoubleSeq p) { - double var = p.get(0); - return RegularSplineComponent.stateComponent(data, var); - } - - @Override - public int parametersCount() { - return 1; - } - - @Override - public ISsfLoading defaultLoading(int m) { - if (m > 0) { - return null; - } - return RegularSplineComponent.loading(data, startpos); - } - - @Override - public int defaultLoadingCount() { - return 1; - } - - @Override - public int stateDim() { - return data.getDim(); - } - - @Override - public boolean isScalable() { - return !v.isFixed(); - } - -} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/SplineItem.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/SplineItem.java index a5853e5..8be8987 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/SplineItem.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/msts/internal/SplineItem.java @@ -16,6 +16,8 @@ */ package jdplus.sts.base.core.msts.internal; +import jdplus.sts.base.core.splines.SplineComponent; +import jdplus.sts.base.core.splines.SplineData; import jdplus.sts.base.core.msts.StateItem; import jdplus.toolkit.base.api.data.DoubleSeq; import jdplus.sts.base.core.msts.MstsMapping; @@ -25,8 +27,6 @@ import jdplus.sts.base.core.msts.ParameterInterpreter; import jdplus.toolkit.base.core.ssf.ISsfLoading; import jdplus.toolkit.base.core.ssf.StateComponent; -import jdplus.toolkit.base.core.ssf.sts.splines.SplineComponent; -import jdplus.toolkit.base.core.ssf.sts.splines.SplineData; /** * @@ -38,7 +38,7 @@ public class SplineItem extends StateItem { private final SplineData data; private final int startpos; - public SplineItem(String name, SplineData data, int startpos,double cvar, boolean fixedvar) { + public SplineItem(String name, SplineData data, int startpos, double cvar, boolean fixedvar) { super(name); v = new VarianceInterpreter(name + ".var", cvar, fixedvar, true); this.data = data; diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/CubicSplines.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/CubicSplines.java new file mode 100644 index 0000000..d0f3558 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/CubicSplines.java @@ -0,0 +1,385 @@ +/* + * To change this license header, choose License Headers in Project Properties. + * To change this template file, choose Tools | Templates + * and open the template in the editor. + */ +package jdplus.sts.base.core.splines; + +import jdplus.toolkit.base.api.data.DoubleSeq; +import jdplus.toolkit.base.api.math.MathException; +import java.util.Arrays; +import java.util.function.DoubleUnaryOperator; +import jdplus.toolkit.base.core.data.DataBlock; + +/** + * + * @author Jean Palate + */ +@lombok.experimental.UtilityClass +public class CubicSplines { + + public static abstract class Spline implements DoubleUnaryOperator { + + final double[] b; + final double[] c; + final double[] d; + final double[] xa, ya; + + Spline(DoubleSeq x, DoubleSeq y) { + this.xa = x.toArray(); + this.ya = y.toArray(); + int size = xa.length - 1; + b = new double[size]; + c = new double[size + 1]; + d = new double[size]; + } + + protected double interpolate(double x, int index) { + double delx = x - xa[index]; + return ya[index] + delx * (b[index] + delx * (c[index] + delx * d[index])); + } + + int size() { + return xa.length; + } + } + + private static class NaturalSpline extends Spline { + + @Override + public double applyAsDouble(double x) { + if (x < xa[0]) { + return extrapolate0(x); + } else if (x > xa[xa.length - 1]) { + return extrapolate1(x); + } else { + int pos = Arrays.binarySearch(xa, x); + if (pos >= 0) { + return ya[pos]; + } else { + return interpolate(x, -pos - 2); + } + } + } + + private double extrapolate0(double x) { + double df = b[0]; + return ya[0] + (x - xa[0]) * df; + } + + private double extrapolate1(double x) { + int n = xa.length - 1; + double dx = xa[n] - xa[n - 1], dx2 = dx * dx; + double df = b[n - 1] + 2 * c[n - 1] * dx + 3 * d[n - 1] * dx2; + return ya[n] + (x - xa[n]) * df; + } + + private NaturalSpline(DoubleSeq X, DoubleSeq Y) { + super(X, Y); + } + } + + public static Spline natural(DoubleSeq X, DoubleSeq Y) { + if (X.length() != Y.length()) { + throw new IllegalArgumentException(); + } + NaturalSpline spline = new NaturalSpline(X, Y); + + int psize = spline.size() - 1; + int sys_size = psize - 1; + /* linear system is sys_size x sys_size */ + double[] xa = spline.xa, ya = spline.ya; + double[] g = new double[psize]; + double[] diag = new double[psize]; + double[] offdiag = new double[psize]; + + for (int i = 0; i < sys_size; i++) { + double h_i = xa[i + 1] - xa[i]; + double h_ip1 = xa[i + 2] - xa[i + 1]; + double ydiff_i = ya[i + 1] - ya[i]; + double ydiff_ip1 = ya[i + 2] - ya[i + 1]; + double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0; + double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0; + offdiag[i] = h_ip1; + diag[i] = 2.0 * (h_ip1 + h_i); + g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i); + } + + if (sys_size == 1) { + spline.c[1] = g[0] / diag[0]; + } else { + DoubleSeq g_vec = DoubleSeq.of(g, 0, sys_size); + DoubleSeq diag_vec = DoubleSeq.of(diag, 0, sys_size); + DoubleSeq offdiag_vec = DoubleSeq.of(offdiag, 0, sys_size - 1); + DataBlock solution_vec = DataBlock.of(spline.c, 1, sys_size + 1); + solveTriDiag(diag_vec, offdiag_vec, g_vec, solution_vec); + + } + + for (int i = 0; i < psize; i++) { + double h_i = xa[i + 1] - xa[i]; + spline.b[i] = (ya[i + 1] - ya[i]) / h_i - h_i / 3.0 * (2.0 * spline.c[i] + spline.c[i + 1]); + spline.d[i] = (spline.c[i + 1] - spline.c[i]) / (3.0 * h_i); + } + return spline; + } + + private static class PeriodicSpline extends Spline { + + private double translate(double x) { + int n = size() - 1; + if (x < xa[0] || x > xa[n]) { + double xc = x - xa[0]; + double dx = xa[n] - xa[0]; + double m = Math.floor(xc / dx); + xc -= m * dx; + return xc + xa[0]; + } else { + return x; + } + } + + @Override + public double applyAsDouble(double x) { + double xc = translate(x); + int pos = Arrays.binarySearch(xa, xc); + if (pos >= 0) { + return ya[pos]; + } else { + return interpolate(xc, -pos - 2); + } + } + + private PeriodicSpline(DoubleSeq X, DoubleSeq Y) { + super(X, Y); + } + + } + + public Spline periodic(DoubleSeq X, DoubleSeq Y) { + int n = Y.length(); + if (n != X.length() || Y.get(0) != Y.get(n - 1)) { + throw new IllegalArgumentException(); + } + Spline spline = new PeriodicSpline(X, Y); + + int psize = spline.size() - 1; + /* Engeln-Mullges + Uhlig "n" */ + int sys_size = psize; + /* linear system is sys_size x sys_size */ + + double[] xa = spline.xa, ya = spline.ya; + + if (sys_size == 2) { + /* solve 2x2 system */ + + double h0 = xa[1] - xa[0]; + double h1 = xa[2] - xa[1]; + + double A = 2.0 * (h0 + h1); + double B = h0 + h1; + + double g0 = 3.0 * ((ya[2] - ya[1]) / h1 - (ya[1] - ya[0]) / h0); + double g1 = 3.0 * ((ya[1] - ya[2]) / h0 - (ya[2] - ya[1]) / h1); + + double det = 3.0 * (h0 + h1) * (h0 + h1); + spline.c[1] = (A * g0 - B * g1) / det; + spline.c[2] = (-B * g0 + A * g1) / det; + spline.c[0] = spline.c[2]; + + } else { + double[] g = new double[psize]; + double[] diag = new double[psize]; + double[] offdiag = new double[psize]; + + for (int i = 0; i < sys_size - 1; i++) { + double h_i = xa[i + 1] - xa[i]; + double h_ip1 = xa[i + 2] - xa[i + 1]; + double ydiff_i = ya[i + 1] - ya[i]; + double ydiff_ip1 = ya[i + 2] - ya[i + 1]; + double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0; + double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0; + offdiag[i] = h_ip1; + diag[i] = 2.0 * (h_ip1 + h_i); + g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i); + } + int ilast = sys_size - 1; + double h_i = xa[ilast + 1] - xa[ilast]; + double h_ip1 = xa[1] - xa[0]; + double ydiff_i = ya[ilast + 1] - ya[ilast]; + double ydiff_ip1 = ya[1] - ya[0]; + double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0; + double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0; + offdiag[ilast] = h_ip1; + diag[ilast] = 2.0 * (h_ip1 + h_i); + g[ilast] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i); + DoubleSeq g_vec = DoubleSeq.of(g, 0, sys_size); + DoubleSeq diag_vec = DoubleSeq.of(diag, 0, sys_size); + DoubleSeq offdiag_vec = DoubleSeq.of(offdiag, 0, sys_size); + DataBlock solution_vec = DataBlock.of(spline.c, 1, sys_size + 1); + + solveCyclicTriDiag(diag_vec, offdiag_vec, g_vec, solution_vec); + spline.c[0] = spline.c[psize]; + } + for (int i = 0; i < psize; i++) { + double h_i = xa[i + 1] - xa[i]; + + spline.b[i] = (ya[i + 1] - ya[i]) / h_i - h_i / 3.0 * (2.0 * spline.c[i] + spline.c[i + 1]); + spline.d[i] = (spline.c[i + 1] - spline.c[i]) / (3.0 * h_i); + } + return spline; + } + + /* for description natural method see [Engeln-Mullges + Uhlig, p. 92] + * + * diag[0] offdiag[0] 0 ..... + * offdiag[0] diag[1] offdiag[1] ..... + * 0 offdiag[1] diag[2] + * 0 0 offdiag[2] ..... + */ + public void solveTriDiag(DoubleSeq diag, DoubleSeq offdiag, DoubleSeq b, DataBlock x) { + int N = diag.length(); + double[] gamma = new double[N]; + double[] alpha = new double[N]; + double[] c = new double[N]; + double[] z = new double[N]; + + /* Cholesky decomposition + A = L.D.L^t + lower_diag(L) = gamma + diag(D) = alpha + */ + alpha[0] = diag.get(0); + gamma[0] = offdiag.get(0) / alpha[0]; + + if (alpha[0] == 0) { + throw new MathException(MathException.DIVBYZERO); + } + + for (int i = 1; i < N - 1; i++) { + alpha[i] = diag.get(i) - offdiag.get(i - 1) * gamma[i - 1]; + gamma[i] = offdiag.get(i) / alpha[i]; + if (alpha[i] == 0) { + throw new MathException(MathException.DIVBYZERO); + } + } + if (N > 1) { + alpha[N - 1] = diag.get(N - 1) - offdiag.get(N - 2) * gamma[N - 2]; + } + /* update RHS */ + z[0] = b.get(0); + for (int i = 1; i < N; i++) { + z[i] = b.get(i) - gamma[i - 1] * z[i - 1]; + } + for (int i = 0; i < N; i++) { + c[i] = z[i] / alpha[i]; + } + + /* backsubstitution */ + x.set(N - 1, c[N - 1]); + if (N >= 2) { + for (int i = N - 2, j = 0; j <= N - 2; j++, i--) { + x.set(i, c[i] - gamma[i] * x.get(i + 1)); + } + } + } + + public void solveCyclicTriDiag(DoubleSeq diag, DoubleSeq offdiag, DoubleSeq b, DataBlock x) { + int N = diag.length(); + double[] delta = new double[N]; + double[] gamma = new double[N]; + double[] alpha = new double[N]; + double[] c = new double[N]; + double[] z = new double[N]; + double sum = 0.0; + + /* factor */ + if (N == 1) { + x.set(0, b.get(0) / diag.get(0)); + return; + } + + alpha[0] = diag.get(0); + if (alpha[0] == 0) { + throw new MathException(MathException.DIVBYZERO); + } + gamma[0] = offdiag.get(0) / alpha[0]; + delta[0] = offdiag.get(N - 1) / alpha[0]; + + for (int i = 1; i < N - 2; i++) { + alpha[i] = diag.get(i) - offdiag.get(i - 1) * gamma[i - 1]; + gamma[i] = offdiag.get(i) / alpha[i]; + delta[i] = -delta[i - 1] * offdiag.get(i - 1) / alpha[i]; + if (alpha[i] == 0) { + throw new MathException(MathException.DIVBYZERO); + } + } + + for (int i = 0; i < N - 2; i++) { + sum += alpha[i] * delta[i] * delta[i]; + } + + alpha[N - 2] = diag.get(N - 2) - offdiag.get(N - 3) * gamma[N - 3]; + gamma[N - 2] = (offdiag.get(N - 2) - offdiag.get(N - 3) * delta[N - 3]) / alpha[N - 2]; + alpha[N - 1] = diag.get(N - 1) - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2]; + + /* update */ + z[0] = b.get(0); + for (int i = 1; i < N - 1; i++) { + z[i] = b.get(i) - z[i - 1] * gamma[i - 1]; + } + sum = 0.0; + for (int i = 0; i < N - 2; i++) { + sum += delta[i] * z[i]; + } + z[N - 1] = b.get(N - 1) - sum - gamma[N - 2] * z[N - 2]; + for (int i = 0; i < N; i++) { + c[i] = z[i] / alpha[i]; + } + + /* backsubstitution */ + x.set(N - 1, c[N - 1]); + x.set(N - 2, c[N - 2] - gamma[N - 2] * x.get(N - 1)); + if (N >= 3) { + for (int i = N - 3, j = 0; j <= N - 3; j++, i--) { + x.set(i, c[i] - gamma[i] * x.get(i + 1) - delta[i] * x.get(N - 1)); + } + } + } + + /* Perform a binary search natural an array natural values. + * + * The parameters index_lo and index_hi provide an initial bracket, + * and it is assumed that index_lo < index_hi. The resulting index + * is guaranteed to be strictly less than index_hi and greater than + * or equal to index_lo, so that the implicit bracket [index, index+1] + * always corresponds to a region within the implicit value range natural + * the value array. + * + * Note that this means the relationship natural 'x' to x_array[index] + * and x_array[index+1] depends on the result region, i.e. the + * behaviour at the boundaries may not correspond to what you + * expect. We have the following complete specification natural the + * behaviour. + * Suppose the input is x_array[] = { x0, x1, ..., xN } + * if ( x == x0 ) then index == 0 + * if ( x > x0 && x <= x1 ) then index == 0, and sim. for other interior pts + * if ( x == xN ) then index == N-1 + * if ( x > xN ) then index == N-1 + * if ( x < x0 ) then index == 0 + */ +// int bsearch(double[] x_array, double x, int index_lo, int index_hi) { +// int ilo = index_lo; +// int ihi = index_hi; +// while (ihi > ilo + 1) { +// int i = (ihi + ilo) / 2; +// if (x_array[i] > x) { +// ihi = i; +// } else { +// ilo = i; +// } +// } +// return ilo; +// } +// +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/DailySpline.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/DailySpline.java new file mode 100644 index 0000000..4395964 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/DailySpline.java @@ -0,0 +1,102 @@ +/* + * 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.sts.base.core.splines; + +import jdplus.toolkit.base.api.data.DoubleSeq; +import jdplus.toolkit.base.api.timeseries.calendars.CalendarUtility; + +/** + * We suppose in this solution that day 58 (in base 0) appears twice in leap + * years -> All years have 365 days A cycle starts at the first of January We + * have to specify the start year and the positions of the days corresponding to + * the nodes. + * + * For a more flexible solution, days should be identified by doubles in the + * future + * + * @author palatej + */ +@lombok.Value +public class DailySpline implements SplineDefinition { + + int startYear; + int[] days; + + @Override + public double getPeriod() { + return 365.0; + } + + @Override + public DoubleSeq nodes() { + return DoubleSeq.onMapping(days.length, i -> days[i]); + } + + @Override + public IntSeq observations(int period) { + boolean lp = CalendarUtility.isLeap(startYear + period); + return lp ? new DailyIntSeq(period * 365) : IntSeq.sequential(period * 365, (period+1)*365); + } + + private static final int LCYCLE = 4 * 365 + 1; + + @Override + public int cycleFor(final int obs) { + int obsc = obs % LCYCLE; + int qobs = 4 * (obs / LCYCLE); + int y = startYear; + boolean leap = false; + int i = 0; + do { + int n = 365; + if (!leap) { + leap = CalendarUtility.isLeap(y); + if (leap) { + n = 366; + } + } + if (obsc < n) { + return qobs + i; + } + obsc -= n; + ++y; + ++i; + } while (true); + + } + + private static class DailyIntSeq implements IntSeq { + + private final int start; + + DailyIntSeq(int start) { + this.start = start; + } + + @Override + public int length() { + return 366; + } + + @Override + public int pos(int idx) { + return idx > 58 ? start + idx - 1 : start + idx; + } + + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/IntSeq.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/IntSeq.java new file mode 100644 index 0000000..4f8fb19 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/IntSeq.java @@ -0,0 +1,53 @@ +/* + * Copyright 2024 JDemetra+. + * 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.sts.base.core.splines; + +/** + * + * @author Jean Palate + */ +public interface IntSeq { + + int length(); + + int pos(int idx); + + static IntSeq sequential(int from, int to){ + return new SequentialIntSeq(from, to); + } + + static class SequentialIntSeq implements IntSeq{ + + private final int start, length; + + SequentialIntSeq(int from, int to){ + this.start=from; + this.length=to-from; + } + + @Override + public int length() { + return length; + } + + @Override + public int pos(int idx) { + return start+idx; + } + + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/RegularSpline.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/RegularSpline.java new file mode 100644 index 0000000..8ff0482 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/RegularSpline.java @@ -0,0 +1,86 @@ +/* + * Copyright 2024 JDemetra+. + * 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.sts.base.core.splines; + +import jdplus.toolkit.base.api.data.DoubleSeq; + +/** + * + * @author Jean Palate + */ +@lombok.Getter +@lombok.AllArgsConstructor(access = lombok.AccessLevel.PRIVATE) +public class RegularSpline implements SplineDefinition { + + public static RegularSpline of(double period) { + return of(period, (int) Math.floor(period)); + } + + public static RegularSpline of(double period, int nnodes) { + double[] dnodes = new double[nnodes]; + double step = period / nnodes; + double cur = step / 2; + dnodes[0] = cur; + int i = 1; + while (i < nnodes) { + cur += step; + dnodes[i++] = cur; + } + return new RegularSpline(period, DoubleSeq.of(dnodes)); + } + + /** + * + * @param period Period of the splines + * @param nodes Should be ordered and in [0, period[. NOT CHECKED ! + * @return + */ + public static RegularSpline of(double period, DoubleSeq nodes) { + return new RegularSpline(period, nodes); + } + + private final double period; + private final DoubleSeq nodes; + + @Override + public double getPeriod() { + return period; + } + + @Override + public DoubleSeq nodes() { + return nodes; + } + + private static final double EPS = 1e-12; + + @Override + public IntSeq observations(int cycle) { + // we need to find the integers in [cycle*period, (cycle+1)*period[ + int i0 = (int) Math.floor(period * cycle + EPS), + i1 = (int) Math.floor(period * (cycle + 1) + EPS); + return IntSeq.sequential(i0, i1); + } + + @Override + public int cycleFor(int obs) { + int c=(int) Math.floor((obs + EPS) / period); + if (obs>=(int) Math.floor(period * (c + 1) + EPS)) + ++c; + return c; + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineComponent.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineComponent.java new file mode 100644 index 0000000..b86788e --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineComponent.java @@ -0,0 +1,253 @@ +/* + * 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.sts.base.core.splines; + +import jdplus.sts.base.core.splines.SplineData.CycleInformation; +import jdplus.toolkit.base.api.data.DoubleSeqCursor; +import jdplus.toolkit.base.core.data.DataBlock; +import jdplus.toolkit.base.core.data.DataBlockIterator; +import jdplus.toolkit.base.core.math.matrices.FastMatrix; +import jdplus.toolkit.base.core.ssf.ISsfDynamics; +import jdplus.toolkit.base.core.ssf.ISsfInitialization; +import jdplus.toolkit.base.core.ssf.ISsfLoading; +import jdplus.toolkit.base.core.ssf.StateComponent; + +/** + * Integer period, regular knots on integer "periods" + * + * @author palatej + */ +@lombok.experimental.UtilityClass +public class SplineComponent { + + public ISsfLoading loading(SplineData data, int startPos) { + + return new Loading(data, startPos); + } + + public class Loading implements ISsfLoading { + + private final SplineData data; + + public Loading(SplineData data, int startpos) { + this.data = data; + } + + private DataBlock z(int pos) { + CycleInformation info = data.informationForObservation(pos); + FastMatrix z = info.getZ(); + int row = pos - info.firstObservation(); + return z.row(row); + } + + @Override + public boolean isTimeInvariant() { + return false; + } + + @Override + public void Z(int pos, DataBlock z) { + z.copy(z(pos)); + } + + @Override + public double ZX(int pos, DataBlock m) { + return z(pos).dot(m); + } + + @Override + public void ZM(int pos, FastMatrix m, DataBlock zm) { + DataBlock row = z(pos); + zm.set(m.columnsIterator(), x -> row.dot(x)); + } + + /** + * Computes M*Z' (or ZM') + * + * @param pos + * @param m + * @param zm + */ + @Override + public void MZt(int pos, FastMatrix m, DataBlock zm) { + DataBlock row = z(pos); + zm.set(m.rowsIterator(), x -> row.dot(x)); + } + + @Override + public double ZVZ(int pos, FastMatrix V) { + DataBlock row = z(pos); + DataBlock zv = DataBlock.make(V.getColumnsCount()); + zv.product(row, V.columnsIterator()); + return zv.dot(row); + } + + @Override + public void VpZdZ(int pos, FastMatrix V, double d) { + if (d == 0) { + return; + } + DataBlockIterator cols = V.columnsIterator(); + DataBlock row = z(pos); + DoubleSeqCursor z = row.cursor(); + while (cols.hasNext()) { + cols.next().addAY(d * z.getAndNext(), row); + } + } + + @Override + public void XpZd(int pos, DataBlock x, double d) { + x.addAY(d, z(pos)); + } + } + + public StateComponent stateComponent(SplineData data, double var, int startpos) { + + Dynamics dynamics = new Dynamics(data, var, startpos); + Initialization initialization = new Initialization(data.getDim()); + + return new StateComponent(initialization, dynamics); + + } + + static class Initialization implements ISsfInitialization { + + private final int dim; + + Initialization(int dim) { + this.dim = dim; + } + + @Override + public int getStateDim() { + return dim; + } + + @Override + public boolean isDiffuse() { + return true; + } + + @Override + public int getDiffuseDim() { + return dim; + } + + @Override + public void diffuseConstraints(FastMatrix b) { + b.diagonal().set(1); + } + + @Override + public void a0(DataBlock a0) { + } + + @Override + public void Pf0(FastMatrix pf0) { + } + + @Override + public void Pi0(FastMatrix pi0) { + pi0.diagonal().set(1); + } + } + + static class Dynamics implements ISsfDynamics { + + private final SplineData data; + private final double var, svar; + + Dynamics(final SplineData data, double var, int startpos) { + this.var = var; + this.svar = Math.sqrt(var); + this.data = data; + } + + private FastMatrix q(int pos) { + CycleInformation info = data.informationForObservation(pos); + return info.getQ().times(var); + } + + private FastMatrix s(int pos) { + CycleInformation info = data.informationForObservation(pos); + return info.getS().times(svar); + } + + @Override + public boolean isTimeInvariant() { + return true; + } + + @Override + public boolean areInnovationsTimeInvariant() { + return true; + } + + @Override + public int getInnovationsDim() { + return data.getDim(); + } + + @Override + public void V(int pos, FastMatrix qm) { + qm.copy(q(pos)); + } + + @Override + public boolean hasInnovations(int pos) { + return true; + } + + @Override + public void S(int pos, FastMatrix sm) { + sm.copy(s(pos)); + } + + @Override + public void addSU(int pos, DataBlock x, DataBlock u) { + x.addProduct(s(pos).rowsIterator(), u); + } + + @Override + public void XS(int pos, DataBlock x, DataBlock xs) { + xs.product(x, s(pos).columnsIterator()); + } + + @Override + public void T(int pos, FastMatrix tr) { + tr.diagonal().set(1); + } + + @Override + public void TX(int pos, DataBlock x) { + } + + @Override + public void XT(int pos, DataBlock x) { + } + + @Override + public void TVT(int pos, FastMatrix v) { + } + + @Override + public void addV(int pos, FastMatrix p) { + p.add(q(pos)); + } + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineData.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineData.java new file mode 100644 index 0000000..6c72577 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineData.java @@ -0,0 +1,110 @@ +/* + * 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.sts.base.core.splines; + +import jdplus.toolkit.base.api.data.DoubleSeqCursor; +import java.util.concurrent.ConcurrentHashMap; +import java.util.concurrent.ConcurrentMap; +import jdplus.toolkit.base.core.data.DataBlock; +import jdplus.toolkit.base.core.math.matrices.FastMatrix; +import jdplus.toolkit.base.core.math.matrices.LowerTriangularMatrix; +import jdplus.toolkit.base.core.math.matrices.SymmetricMatrix; + +/** + * + * @author palatej + */ +public class SplineData { + + @lombok.Getter + @lombok.AllArgsConstructor(access = lombok.AccessLevel.PRIVATE) + public static class CycleInformation { + + /** + * Q = variance of the innovations S = Cholesky factor of Q. Q, S have + * same size, corresponding to the number of nodes (-1) (constant + * between years) Z = loadings The number of rows in Z corresponds to + * the number of obs. Could vary between cycles + */ + FastMatrix Q, S, Z; + IntSeq observations; + + public int firstObservation() { + return observations.pos(0); + } + } + + public SplineData(SplineDefinition definition) { + this.definition = definition; + splines = definition.splines(); + dim = splines.length - 1; + } + + public CycleInformation informationForCycle(int cycle) { + return information(cycle); + } + + public CycleInformation informationForObservation(int pos) { + return information(definition.cycleFor(pos)); + } + + private CycleInformation information(int cycle) { + CycleInformation info = infos.get(cycle); + if (info != null) { + return info; + } + int cdim = splines.length; + double[] wstar = new double[cdim]; + IntSeq obs = definition.observations(cycle); + int m = obs.length(); + FastMatrix Z = FastMatrix.make(m, cdim); + for (int i = 0; i < cdim; ++i) { + DoubleSeqCursor.OnMutable cursor = Z.column(i).cursor(); + double s = 0; + for (int j = 0; j < m; ++j) { + double w = splines[i].applyAsDouble(obs.pos(j)); + cursor.setAndNext(w); + s += w; + } + wstar[i] = s; + } + DataBlock zh = Z.column(cdim - 1); + double wh = wstar[cdim - 1]; + for (int i = 0; i < cdim - 1; ++i) { + Z.column(i).addAY(-wstar[i] / wh, zh); + } + + DataBlock W = DataBlock.of(wstar, 0, cdim); + FastMatrix Q = FastMatrix.identity(cdim - 1); + Q.addXaXt(-1 / W.ssq(), W.drop(0, 1)); + FastMatrix S = Q.deepClone(); + SymmetricMatrix.lcholesky(S, 1e-9); + LowerTriangularMatrix.toLower(S); + info = new CycleInformation(Q, S, Z.dropBottomRight(0, 1), obs); + infos.put(cycle, info); + return info; + } + + public int getDim() { + return dim; + } + + private final SplineDefinition definition; + private final ConcurrentMap infos = new ConcurrentHashMap<>(); + private final CubicSplines.Spline[] splines; + private final int dim; +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineDefinition.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineDefinition.java new file mode 100644 index 0000000..f4f9524 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/jdplus/sts/base/core/splines/SplineDefinition.java @@ -0,0 +1,81 @@ +/* + * 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.sts.base.core.splines; + +import jdplus.toolkit.base.api.data.DoubleSeq; + +/** + * + * @author palatej + */ +public interface SplineDefinition { + + /** + * Length of a cycle + * + * @return + */ + double getPeriod(); + + /** + * Nodes of the spline (same for each cycle). Nodes don't necessary + * correspond to observations. The nodes are defined for the first cycle + * ([0, period[ + * + * @return + */ + DoubleSeq nodes(); + + /** + * Positions of the observations for the given cycle + * Must be inside [cycle*period, (cycle+1)*period[. + * + * @param cycle + * @return + */ + IntSeq observations(int cycle); + + /** + * Cycle corresponding to the given observation position + * @param obs Observation position + * @return + */ + int cycleFor(int obs); + + default CubicSplines.Spline[] splines() { + DoubleSeq nodes = nodes(); + int dim = nodes.length(); + + CubicSplines.Spline[] splines = new CubicSplines.Spline[dim]; + double[] xi = new double[dim + 1]; + nodes.copyTo(xi, 0); + double x0 = xi[0]; + xi[dim] = x0 + getPeriod(); + + for (int i = 0; i < dim; ++i) { + double[] f = new double[dim + 1]; + if (i == 0) { + f[0] = 1; + f[dim] = 1; + } else { + f[i] = 1; + } + splines[i] = CubicSplines.periodic(DoubleSeq.of(xi), DoubleSeq.of(f)); + } + return splines; + } +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/module-info.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/module-info.java index 5913f4b..d0ef763 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/module-info.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/main/java/module-info.java @@ -17,6 +17,7 @@ exports jdplus.sts.base.core.msts.internal; exports jdplus.sts.base.core.msts.survey; exports jdplus.sts.base.core.extractors; + exports jdplus.sts.base.core.splines; exports jdplus.sts.base.core; provides InformationExtractor with diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/StsKernelTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/StsKernelTest.java index 92b43b7..684e83b 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/StsKernelTest.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/StsKernelTest.java @@ -47,8 +47,8 @@ public void testProd() { rslts.getFinals().getSeries(ComponentType.Seasonal, ComponentInformation.Value), rslts.getFinals().getSeries(ComponentType.Irregular, ComponentInformation.Value) )); - System.out.println(table); - System.out.println(rslts.getSts().getBsm()); +// System.out.println(table); +// System.out.println(rslts.getSts().getBsm()); } diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/msts/CompositeModelTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/msts/CompositeModelTest.java index f1edb13..28cf77d 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/msts/CompositeModelTest.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/msts/CompositeModelTest.java @@ -121,9 +121,9 @@ public void testVAT_NOBUG() { M.mul(10); CompositeModelEstimation rslt = model.estimate(M, false, false, SsfInitialization.SqrtDiffuse, Optimizer.LevenbergMarquardt, 1e-15, null); StateStorage states = rslt.getSmoothedStates(); - System.out.println(states.getComponent(0)); +// System.out.println(states.getComponent(0)); // System.out.println(states.getComponent(2)); - System.out.println(DoubleSeq.of(rslt.getFullParameters())); +// System.out.println(DoubleSeq.of(rslt.getFullParameters())); // System.out.println(rslt.getLikelihood().factor()); } diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/CubicSplinesTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/CubicSplinesTest.java new file mode 100644 index 0000000..5a5a811 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/CubicSplinesTest.java @@ -0,0 +1,34 @@ +/* + * Copyright 2024 JDemetra+. + * 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.sts.base.core.splines; + +import org.junit.jupiter.api.Test; +import static org.junit.jupiter.api.Assertions.*; + +/** + * + * @author Jean Palate + */ +public class CubicSplinesTest { + + public CubicSplinesTest() { + } + + @Test + public void testSomeMethod() { + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/DailySplineTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/DailySplineTest.java new file mode 100644 index 0000000..b237238 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/DailySplineTest.java @@ -0,0 +1,45 @@ +/* + * Copyright 2024 JDemetra+. + * 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.sts.base.core.splines; + +import org.junit.jupiter.api.Test; +import static org.junit.jupiter.api.Assertions.*; + +/** + * + * @author Jean Palate + */ +public class DailySplineTest { + + public DailySplineTest() { + } + + @Test + public void testPositive() { + DailySpline rs = new DailySpline(1963, new int[]{0, 100, 200, 300}); + int n = 0; + for (int i = 0; i < 120; ++i) { + IntSeq observations = rs.observations(i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(n+j) == i); + } + n += nc; + } + assertTrue(n == 365 * 120 + 30); + } + +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/RegularSplineTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/RegularSplineTest.java new file mode 100644 index 0000000..c22fff8 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/RegularSplineTest.java @@ -0,0 +1,116 @@ +/* + * Copyright 2024 JDemetra+. + * 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.sts.base.core.splines; + +import org.junit.jupiter.api.Test; +import static org.junit.jupiter.api.Assertions.*; + +/** + * + * @author Jean Palate + */ +public class RegularSplineTest { + + public RegularSplineTest() { + } + + @Test + public void testInt() { + RegularSpline rs = RegularSpline.of(25); + int n = 0; + for (int i = 0; i < 10000; ++i) { + IntSeq observations = rs.observations(i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == i); + } + n += nc; + } + assertTrue(n == 250000); + } + + @Test + public void testDouble() { + RegularSpline rs = RegularSpline.of(24.6); + int n = 0; + for (int i = 0; i < 10000; ++i) { + IntSeq observations = rs.observations(i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == i); + } +// System.out.println(nc); + n += nc; + } + assertTrue(n == 246000); + } + + @Test + public void testDouble2() { + RegularSpline rs = RegularSpline.of(24.25); + int n = 0; + for (int i = 0; i < 12000; ++i) { + IntSeq observations = rs.observations(i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == i); + } +// System.out.println(nc); + n += nc; + } + assertTrue(n == 291000); + } + + @Test + public void testNegative() { + RegularSpline rs = RegularSpline.of(24.6); + int n = 0; + for (int i = 1; i < 10001; ++i) { + IntSeq observations = rs.observations(-i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == -i); + } +// System.out.println(nc); + n += nc; + } + assertTrue(n == 246000); + rs = RegularSpline.of(24); + n = 0; + for (int i = 1; i < 10001; ++i) { + IntSeq observations = rs.observations(-i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == -i); + } +// System.out.println(nc); + n += nc; + } + assertTrue(n == 240000); + rs = RegularSpline.of(24.25); + n = 0; + for (int i = 0; i < 12000; ++i) { + IntSeq observations = rs.observations(-i); + int nc = observations.length(); + for (int j = 0; j < nc; ++j) { + assertTrue(rs.cycleFor(observations.pos(j)) == -i); + } +// System.out.println(nc); + n += nc; + } + assertTrue(n == 291000); + } +} diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/SplineComponentTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/SplineComponentTest.java new file mode 100644 index 0000000..6913700 --- /dev/null +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-core/src/test/java/jdplus/sts/base/core/splines/SplineComponentTest.java @@ -0,0 +1,53 @@ +/* + * Click nbfs://nbhost/SystemFileSystem/Templates/Licenses/license-default.txt to change this license + * Click nbfs://nbhost/SystemFileSystem/Templates/UnitTests/JUnit5TestClass.java to edit this template + */ +package jdplus.sts.base.core.splines; + +import jdplus.sts.base.core.splines.SplineComponent; +import jdplus.sts.base.core.splines.RegularSpline; +import jdplus.sts.base.core.splines.SplineData; +import jdplus.toolkit.base.api.data.DoubleSeq; +import jdplus.toolkit.base.core.ssf.ISsfLoading; +import jdplus.toolkit.base.core.ssf.akf.AkfToolkit; +import jdplus.toolkit.base.core.ssf.akf.SmoothingOutput; +import jdplus.toolkit.base.core.ssf.composite.CompositeSsf; +import jdplus.toolkit.base.core.ssf.sts.LocalLinearTrend; +import jdplus.toolkit.base.core.ssf.sts.Noise; +import jdplus.toolkit.base.core.ssf.univariate.SsfData; +import org.junit.jupiter.api.Test; +import tck.demetra.data.Data; + +/** + * + * @author palat + */ +public class SplineComponentTest { + + public SplineComponentTest() { + } + + @Test + public void testMonthly() { + double[] PROD = Data.PROD; + RegularSpline rs=RegularSpline.of(12, DoubleSeq.of(0,3,6,9)); + SplineData sd=new SplineData(rs); + CompositeSsf ssf = CompositeSsf.builder() + .add(LocalLinearTrend.stateComponent(0.1, 0.1), LocalLinearTrend.defaultLoading()) + .add(Noise.of(1), Noise.defaultLoading()) + .add(SplineComponent.stateComponent(sd,3,0), SplineComponent.loading(sd, 0)) + .build(); + + SmoothingOutput rslt = AkfToolkit.robustSmooth(ssf, new SsfData(PROD), true, true); +// System.out.println(rslt.getSmoothing().getComponent(ssf.componentsPosition()[0])); +// System.out.println(rslt.getSmoothing().getComponent(ssf.componentsPosition()[1])); + + ISsfLoading loading = SplineComponent.loading(sd, 0); + for (int i=0; i MAPPING = new InformationMapping() { @Override public Class getSourceClass() { return CompositeModelEstimation.class; } }; - + static { MAPPING.set("likelihood.ll", Double.class, source -> source.getLikelihood().logLikelihood()); MAPPING.set("likelihood.ser", Double.class, source -> source.getLikelihood().ser()); @@ -118,6 +118,12 @@ public Class getSourceClass() { StateStorage smoothedStates = source.getSmoothedStates(); return smoothedStates.getComponentVariance(source.getCmpPos()[p]).toArray(); }); + MAPPING.setArray("ssf.smoothing.components", 0, Matrix.class, (source, p) -> { + return source.getSmoothedComponents(p); + }); + MAPPING.setArray("ssf.smoothing.vcomponents", 0, Matrix.class, (source, p) -> { + return source.getSmoothedComponentVariance(p); + }); MAPPING.setArray("ssf.smoothing.state", 0, double[].class, (source, p) -> { StateStorage smoothedStates = source.getSmoothedStates(); return smoothedStates.a(p).toArray(); @@ -186,7 +192,7 @@ public Class getSourceClass() { StateStorage fStates = source.getFilteringStates(); return fStates.P(p).unmodifiable(); }); - + MAPPING.setArray("ssf.filtered.array", 0, double[].class, (source, p) -> { StateStorage fStates = source.getFilteredStates(); return fStates.getComponent(p).toArray(); @@ -229,74 +235,74 @@ public Class getSourceClass() { } return Matrix.of(z, n, m); }); - + } - + @Override public boolean contains(String id) { return MAPPING.contains(id); } - + @Override public Map getDictionary() { Map dic = new LinkedHashMap<>(); MAPPING.fillDictionary(null, dic, true); return dic; } - + @Override public T getData(String id, Class tclass) { return MAPPING.getData(estimation, id, tclass); } - + public static final InformationMapping getMapping() { return MAPPING; } - + public double[] signal(int obs, int[] cmps) { return estimation.signal(obs, cmps).toArray(); } - + public double[] stdevSignal(int obs, int[] cmps) { return estimation.stdevSignal(obs, cmps).toArray(); } - + public double[] signal(Matrix m, int[] pos) { - FastMatrix M=FastMatrix.of(m); - return (pos == null ? estimation.signal(M): estimation.signal(M, pos)).toArray(); + FastMatrix M = FastMatrix.of(m); + return (pos == null ? estimation.signal(M) : estimation.signal(M, pos)).toArray(); } - + public double[] stdevSignal(Matrix m, int[] pos) { - FastMatrix M=FastMatrix.of(m); - return (pos == null ? estimation.stdevSignal(M): estimation.stdevSignal(M, pos)).toArray(); + FastMatrix M = FastMatrix.of(m); + return (pos == null ? estimation.stdevSignal(M) : estimation.stdevSignal(M, pos)).toArray(); } - + public FastMatrix loading(int obs) { return estimation.loading(obs, null); } - + public MultivariateSsf ssf() { return estimation.getSsf(); } - + public StateStorage smoothedStates() { return estimation.getSmoothedStates(); } - + public StateStorage filteredStates() { return estimation.getFilteredStates(); } - + public StateStorage filteringStates() { return estimation.getFilteringStates(); } } - + public Results estimate(CompositeModel model, Matrix data, boolean marginal, boolean rescaling, String initialization, String opt, double eps, double[] parameters) { return new Results(model.estimate(FastMatrix.of(data), marginal, rescaling, SsfInitialization.valueOf(initialization), Optimizer.valueOf(opt), eps, parameters)); } - + public Results compute(CompositeModel model, Matrix data, double[] parameters, boolean marginal, boolean concentrated) { return new Results(model.compute(FastMatrix.of(data), parameters, marginal, concentrated)); } diff --git a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-r/src/test/java/jdplus/sts/base/r/CompositeModelTest.java b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-r/src/test/java/jdplus/sts/base/r/CompositeModelTest.java index 4122f12..d56587f 100644 --- a/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-r/src/test/java/jdplus/sts/base/r/CompositeModelTest.java +++ b/jdplus-incubator-base/jdplus-sts-base-parent/jdplus-sts-base-r/src/test/java/jdplus/sts/base/r/CompositeModelTest.java @@ -60,10 +60,10 @@ public void testAirline() { FastMatrix M = FastMatrix.make(len, 1); M.column(0).copyFrom(Data.PROD, 0); CompositeModelEstimation rslt = model.estimate(M, false, true, SsfInitialization.Augmented_NoCollapsing, Optimizer.LevenbergMarquardt, 1e-15, null); - System.out.println(DataBlock.of(rslt.getFullParameters())); - System.out.println(rslt.getSmoothedStates().getComponent(15)); - System.out.println(rslt.getSmoothedStates().getComponentVariance(15)); - System.out.println(rslt.getLikelihood().logLikelihood()); +// System.out.println(DataBlock.of(rslt.getFullParameters())); +// System.out.println(rslt.getSmoothedStates().getComponent(15)); +// System.out.println(rslt.getSmoothedStates().getComponentVariance(15)); +// System.out.println(rslt.getLikelihood().logLikelihood()); } @Test @@ -84,7 +84,7 @@ public void testBsm() { System.out.println(DataBlock.of(rslt.getFullParameters())); // System.out.println(rslt.getSmoothedStates().getComponent(0)); // System.out.println(rslt.getSmoothedStates().getComponentVariance(0)); - System.out.println(rslt.getLikelihood().logLikelihood()); +// System.out.println(rslt.getLikelihood().logLikelihood()); } @Test @@ -118,7 +118,7 @@ public void testBsmVar() { System.out.println(DataBlock.of(rslt.getFullParameters())); // System.out.println(rslt.getSmoothedStates().getComponent(0)); // System.out.println(rslt.getSmoothedStates().getComponent(2)); - System.out.println(rslt.getLikelihood().logLikelihood()); +// System.out.println(rslt.getLikelihood().logLikelihood()); } @Test diff --git a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-api/src/main/java/jdplus/x12plus/base/api/X11plusDictionaries.java b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-api/src/main/java/jdplus/x12plus/base/api/X11plusDictionaries.java index bb3d1a9..5ef041c 100644 --- a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-api/src/main/java/jdplus/x12plus/base/api/X11plusDictionaries.java +++ b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-api/src/main/java/jdplus/x12plus/base/api/X11plusDictionaries.java @@ -27,20 +27,20 @@ public class X11plusDictionaries { public final Dictionary BTABLES = AtomicDictionary.builder() .name(B_TABLES) - .item(AtomicDictionary.Item.builder().name(B1).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B2).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B3).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B4).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B5).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B6).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B7).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B8).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B9).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B10).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B11).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B13).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B17).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(B20).description("TODO").outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B1).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B2).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B3).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B4).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B5).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B6).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B7).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B8).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B9).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B10).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B11).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B13).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B17).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(B20).outputClass(TsData.class).build()) .build(); public final String C1="c1", C2="c2", C4="c4", C5="c5", @@ -51,18 +51,18 @@ public class X11plusDictionaries { public final Dictionary CTABLES = AtomicDictionary.builder() .name(C_TABLES) - .item(AtomicDictionary.Item.builder().name(C1).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C2).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C4).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C5).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C6).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C7).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C9).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C10).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C11).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C13).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C17).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(C20).description("TODO").outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C1).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C2).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C4).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C5).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C6).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C7).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C9).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C10).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C11).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C13).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C17).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(C20).outputClass(TsData.class).build()) .build(); public final String D1="d1", D2="d2", D4="d4", @@ -74,21 +74,21 @@ public class X11plusDictionaries { public final Dictionary DTABLES = AtomicDictionary.builder() .name(D_TABLES) - .item(AtomicDictionary.Item.builder().name(D1).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D2).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D4).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D5).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D6).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D7).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D8).description("TODO").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D9).description("TODO").outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D1).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D2).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D4).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D5).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D6).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D7).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D8).outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D9).outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D10).description("seasonal component (linearized)").outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D10A).description("forecasts of the seasonal component").outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D11).description("seasonally adjusted series (linearized)").outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D11A).description("forecasts of the seasonally adjusted component").outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D12).description("trend component (llinearized)").outputClass(TsData.class).build()) .item(AtomicDictionary.Item.builder().name(D12A).description("forecasts of the trend component").outputClass(TsData.class).build()) - .item(AtomicDictionary.Item.builder().name(D13).description("TODO").outputClass(TsData.class).build()) + .item(AtomicDictionary.Item.builder().name(D13).outputClass(TsData.class).build()) .build(); diff --git a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-core/src/test/java/jdplus/x12plus/base/core/x12/X12plusKernelTest.java b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-core/src/test/java/jdplus/x12plus/base/core/x12/X12plusKernelTest.java index 60725f4..bb44941 100644 --- a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-core/src/test/java/jdplus/x12plus/base/core/x12/X12plusKernelTest.java +++ b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-core/src/test/java/jdplus/x12plus/base/core/x12/X12plusKernelTest.java @@ -62,7 +62,7 @@ public void testTimeVaryingTD() { .build(); X12plusKernel kernel = X12plusKernel.of(spec, null); X12plusResults rslt = kernel.process(Data.TS_ABS_RETAIL, null); - System.out.println(rslt.getMtdCorrection().getTdCoefficients()); +// System.out.println(rslt.getMtdCorrection().getTdCoefficients()); } @Test diff --git a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-r/src/test/java/jdplus/x11plus/base/r/X11DecompositionTest.java b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-r/src/test/java/jdplus/x11plus/base/r/X11DecompositionTest.java index c37fe6d..266152e 100644 --- a/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-r/src/test/java/jdplus/x11plus/base/r/X11DecompositionTest.java +++ b/jdplus-incubator-base/jdplus-x12plus-base-parent/jdplus-x12plus-base-r/src/test/java/jdplus/x11plus/base/r/X11DecompositionTest.java @@ -26,19 +26,19 @@ public void testX11() { double[] s = Data.RETAIL_BOOKSTORES; X11Decomposition.Results r1 = X11Decomposition.cnX11(s, 12, true, 11, 2, "Henderson", 3, "Epanechnikov", 1.5, 2.5); assertTrue(r1 != null); - System.out.println(DoubleSeq.of(r1.getData(t, double[].class))); +// System.out.println(DoubleSeq.of(r1.getData(t, double[].class))); X11Decomposition.Results r2 = X11Decomposition.dafX11(s, 12, true, 11, 2, "Henderson", 3, "Epanechnikov", 1.5, 2.5); assertTrue(r2 != null); - System.out.println(DoubleSeq.of(r2.getData(t, double[].class))); +// System.out.println(DoubleSeq.of(r2.getData(t, double[].class))); // X11Decomposition.Results r3 = X11Decomposition.rkhsX11(s, 12, true, 11, 2, "BiWeight", true, "FrequencyResponse", true, Math.PI / 8, 3, "Henderson", 1.5, 2.5); // assertTrue(r3 != null); // System.out.println(DoubleSeq.of(r3.getData(t, double[].class))); X11Decomposition.Results r4 = X11Decomposition.lpX11(s, 12, true, 11, 2, "Henderson", 0, new double[]{}, 10, Math.PI / 8, 3, "Henderson", 1.5, 2.5); assertTrue(r4 != null); - System.out.println(DoubleSeq.of(r4.getData(t, double[].class))); +// System.out.println(DoubleSeq.of(r4.getData(t, double[].class))); X11Decomposition.Results r5 = X11Decomposition.process(s, 12, true, 11, 2, "Henderson", "MMSRE", "S3X5", "S3X5", 1.5, 2.5); assertTrue(r5 != null); - System.out.println(DoubleSeq.of(r5.getData(t, double[].class))); +// System.out.println(DoubleSeq.of(r5.getData(t, double[].class))); } }