From e03091e42a526943754943ee24cd84badac99ac7 Mon Sep 17 00:00:00 2001 From: apete Date: Sun, 22 Oct 2023 13:27:32 +0200 Subject: [PATCH] FFT --- CHANGELOG.md | 18 +- .../java/org/ojalgo/data/DataProcessors.java | 23 +- .../java/org/ojalgo/data/batch/BatchNode.java | 7 + .../java/org/ojalgo/data/image/ImageData.java | 167 +++- .../ojalgo/data/transform/DataTransform.java | 33 + .../transform/DiscreteFourierTransform.java | 803 ++++++++++++++++++ .../org/ojalgo/data/transform/ZTransform.java | 153 ++++ .../ojalgo/function/PrimitiveFunction.java | 62 +- .../org/ojalgo/function/UnaryFunction.java | 9 +- .../polynomial/AbstractPolynomial.java | 230 ++++- .../function/polynomial/BigPolynomial.java | 4 +- .../polynomial/ComplexPolynomial.java | 4 +- .../function/polynomial/PolynomialC128.java | 55 +- .../polynomial/PolynomialFunction.java | 21 +- .../function/polynomial/PolynomialQ128.java | 55 +- .../function/polynomial/PolynomialR032.java | 90 +- .../function/polynomial/PolynomialR064.java | 89 +- .../function/polynomial/PolynomialR128.java | 55 +- .../function/polynomial/PolynomialR256.java | 69 +- .../polynomial/PrimitivePolynomial.java | 4 +- .../polynomial/RationalPolynomial.java | 4 +- .../function/polynomial/ScalarPolynomial.java | 103 +++ .../ojalgo/function/series/FourierSeries.java | 162 ++++ .../function/series/PeriodicFunction.java | 149 ++++ .../java/org/ojalgo/matrix/BasicMatrix.java | 17 +- .../ojalgo/matrix/store/PhysicalStore.java | 66 +- .../org/ojalgo/matrix/store/SparseStore.java | 24 +- .../org/ojalgo/matrix/store/Subregion2D.java | 69 +- .../matrix/store/TransformableRegion.java | 12 + .../org/ojalgo/matrix/store/WrapperStore.java | 64 +- .../java/org/ojalgo/scalar/ComplexNumber.java | 96 ++- .../java/org/ojalgo/structure/Access1D.java | 2 +- src/main/java/org/ojalgo/type/TypeUtils.java | 32 +- src/test/java/org/ojalgo/TestUtils.java | 17 + .../ojalgo/data/transform/DFTBenchmark.java | 74 ++ .../data/transform/DataTransformTests.java | 28 + .../DiscreteFourierTransformTest.java | 469 ++++++++++ .../ojalgo/data/transform/ZtransformTest.java | 132 +++ .../polynomial/PolynomialImplTest.java | 41 +- .../function/series/FourierSeriesTest.java | 203 +++++ .../function/series/PeriodicFunctionTest.java | 132 +++ .../function/series/SeriesFunctionTests.java | 29 + .../org/ojalgo/scalar/ComplexNumberTest.java | 27 +- 43 files changed, 3550 insertions(+), 353 deletions(-) create mode 100644 src/main/java/org/ojalgo/data/transform/DataTransform.java create mode 100644 src/main/java/org/ojalgo/data/transform/DiscreteFourierTransform.java create mode 100644 src/main/java/org/ojalgo/data/transform/ZTransform.java create mode 100644 src/main/java/org/ojalgo/function/polynomial/ScalarPolynomial.java create mode 100644 src/main/java/org/ojalgo/function/series/FourierSeries.java create mode 100644 src/main/java/org/ojalgo/function/series/PeriodicFunction.java create mode 100644 src/test/java/org/ojalgo/data/transform/DFTBenchmark.java create mode 100644 src/test/java/org/ojalgo/data/transform/DataTransformTests.java create mode 100644 src/test/java/org/ojalgo/data/transform/DiscreteFourierTransformTest.java create mode 100644 src/test/java/org/ojalgo/data/transform/ZtransformTest.java create mode 100644 src/test/java/org/ojalgo/function/series/FourierSeriesTest.java create mode 100644 src/test/java/org/ojalgo/function/series/PeriodicFunctionTest.java create mode 100644 src/test/java/org/ojalgo/function/series/SeriesFunctionTests.java diff --git a/CHANGELOG.md b/CHANGELOG.md index 1aeb89fb1..638079589 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,13 +11,29 @@ Added / Changed / Deprecated / Fixed / Removed / Security > Corresponds to changes in the `develop` branch since the last release +### Added + +#### org.ojalgo.data + +- There is a new package `org.ojalgo.data.transform` that (so far) contain implementations of the `ZTransform` and the `DiscreteFourierTransform`. Most notably ojAlgo now contains an implementation of the FFT algorithm. +- `ImageData` has been extended with features to work with Fourier transforms of images, as well as a Gaussian blur function. + +#### org.ojalgo.function + +- The `PolynomialFunction` interface now extends `org.ojalgo.algebra.Ring` which means it is now possible to `add` and `multiply` polynomials. There's also been some refactoring of the internals of the polynomial implementations. +- Two new `PrimitiveFunction.Unary` implementations `PeriodicFunction` and `FourierSeries`. It is now possible to take any (periodic) function and, via FFT, get a Fourier series approximation of it (1 method call). + +#### org.ojalgo.scalar + +- Some minor additions to `ComplexNumber`. It is now possible to do complex number multiplication without creating new instances (calculate the real and imaginary parts separately). Added utilities to get the unit roots. + ## [53.1.1] – 2023-10-16 ### Added #### org.ojalgo.data -- New `ImageData` class that wraps a `java.awt.image.BufferedImage` and implements `MatrixStore`. Further it adds a few utility methods to simplify working image data - convert to grey scale, re-sample (change size), separate the colour channels... +- New `ImageData` class that wraps a `java.awt.image.BufferedImage` and implements `MatrixStore`. Further it adds a few utility methods to simplify working with image data - convert to grey scale, re-sample (change size), separate the colour channels... #### org.ojalgo.matrix diff --git a/src/main/java/org/ojalgo/data/DataProcessors.java b/src/main/java/org/ojalgo/data/DataProcessors.java index ab7471ac7..f67d040e6 100644 --- a/src/main/java/org/ojalgo/data/DataProcessors.java +++ b/src/main/java/org/ojalgo/data/DataProcessors.java @@ -52,27 +52,27 @@ public class DataProcessors { /** * Variables centered so that their average will be 0.0 */ - public static final Transformation2D CENTER = DataProcessors.newTransformation2D(ss -> SUBTRACT.by(ss.getMean())); + public static final Transformation2D CENTER = DataProcessors.newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean())); /** * Variables will be centered around 0.0 AND scaled to be [-1.0,1.0]. The minimum value will be * transformed to -1.0 and the maximum to +1.0. */ public static final Transformation2D CENTER_AND_SCALE = DataProcessors - .newTransformation2D(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by((ss.getMaximum() - ss.getMinimum()) / TWO))); + .newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by((ss.getMaximum() - ss.getMinimum()) / TWO))); /** * Variables scaled to be within [-1.0,1.0] (divide by largest magnitude regardless of sign). If all * values are positive the range will within [0.0,1.0]. If all are negative the range will be within * [-1.0,0.0] */ - public static final Transformation2D SCALE = DataProcessors.newTransformation2D(ss -> DIVIDE.by(ss.getLargest())); + public static final Transformation2D SCALE = DataProcessors.newColumnsTransformer(ss -> DIVIDE.by(ss.getLargest())); /** * Will normalise each variable - replace each value with its standard score. */ public static final Transformation2D STANDARD_SCORE = DataProcessors - .newTransformation2D(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by(ss.getStandardDeviation()))); + .newColumnsTransformer(ss -> SUBTRACT.by(ss.getMean()).andThen(DIVIDE.by(ss.getStandardDeviation()))); /** * Calculate the correlation matrix from a set of variables' samples. Each {@link Access1D} instance @@ -237,7 +237,20 @@ public static > M covariances(final Factory2D return retVal; } - public static Transformation2D newTransformation2D(final Function> definition) { + /** + * Creates a {@link Transformation2D} that will apply a {@link UnaryFunction} to each column. That unary + * function will be created by the provided {@link Function} using a {@link SampleSet} (of that column) as + * input. + *

+ * The constants {@link #CENTER}, {@link #SCALE}, {@link #CENTER_AND_SCALE} and {@link #STANDARD_SCORE} + * are predefined {@link Transformation2D} instances created by calling this method. + * + * @param definition A {@link Function} that will create a {@link UnaryFunction} from a {@link SampleSet} + * to be applied to each column + * @return A {@link Transformation2D} that will apply a {@link UnaryFunction} to each column + */ + public static Transformation2D newColumnsTransformer(final Function> definition) { + return new Transformation2D<>() { public > void transform(final T transformable) { diff --git a/src/main/java/org/ojalgo/data/batch/BatchNode.java b/src/main/java/org/ojalgo/data/batch/BatchNode.java index 332f8edf6..b87a9265c 100644 --- a/src/main/java/org/ojalgo/data/batch/BatchNode.java +++ b/src/main/java/org/ojalgo/data/batch/BatchNode.java @@ -51,6 +51,9 @@ * A batch processing data node for when there's no way to fit the data in memory. *

* Data is stored in sharded files, and data is written/consumed and processed concurrently. + *

+ * The data is processed in batches. Each batch is processed in a single thread. The number of threads is + * controlled by {@link #parallelism(IntSupplier)}. */ public final class BatchNode { @@ -182,18 +185,22 @@ private static final class TwoStepWrapper implements TwoStepMapper */ -public class ImageData implements MatrixStore, Mutate2D.Fillable { +public class ImageData implements MatrixStore, Mutate2D.Receiver { + + @FunctionalInterface + public static interface FrequencyDomainUpdater { + + /** + * Used as a callback from {@link #newTransformation(FrequencyDomainUpdater)}. + * + * @param distance The distance from the centre. The centre is the zero frequency. The distance is + * scaled so that the radius of the largest circle that can fit in the image rectangle is 100. + * @param value The current value + * @return The new value + */ + ComplexNumber invoke(double distance, ComplexNumber value); + + } static final class SingleChannel extends ImageData { @@ -102,6 +127,20 @@ public static ImageData copy(final Access2D values) { return retVal; } + /** + * Creates a new image, transforming the input (back) from the frequency domain to the spatial domain + * using the inverse discrete Fourier transform. The transformation also reverts the shift of the + * frequency representation – the input should be in a format as returned by {@link #toFrequencyDomain()}. + * The new image will be of the same dimensions as the input matrix. + * + * @see #toFrequencyDomain() + */ + public static ImageData fromFrequencyDomain(final MatrixStore transformed) { + ImageData retVal = ImageData.newGreyScale(transformed); + retVal.fillMatching(DiscreteFourierTransform.inverse2D(DiscreteFourierTransform.shift(transformed))); + return retVal; + } + public static ImageData newColour(final int nbRows, final int nbCols) { return new ImageData(new BufferedImage(nbCols, nbRows, BufferedImage.TYPE_INT_ARGB)); } @@ -118,6 +157,55 @@ public static ImageData newGreyScale(final Structure2D shape) { return ImageData.newGreyScale(shape.getRowDim(), shape.getColDim()); } + /** + * Creates a new transformation that can be used to transform a matrix of complex numbers in the frequency + * domain. The transformation will invoke the updater for each element in the input matrix, passing the + * distance from the centre of the matrix as the first argument. The distance is scaled so that the radius + * of the largest circle that can fit in the image rectangle is 100. + */ + public static Transformation2D newTransformation(final FrequencyDomainUpdater updater) { + + return new Transformation2D<>() { + + public > void transform(final T transformable) { + + int nbRows = transformable.getRowDim(); + int nbCols = transformable.getColDim(); + + double midRow = (nbRows - 1) / 2D; + double midCol = (nbCols - 1) / 2D; + + double scale = Math.min(midRow, midCol) / 100D; + + for (int j = 0; j < nbCols; j++) { + for (int i = 0; i < nbRows; i++) { + transformable.set(i, j, updater.invoke(MissingMath.hypot(i - midRow, j - midCol) / scale, transformable.get(i, j))); + } + } + } + }; + } + + /** + * Converts a matrix of complex numbers to an image of its power spectrum (log10 of the squared norms). In + * addition there is scaling applied to adapt the values towards the range [0,255]. + *

+ * This is meant to be used with the output from {@link #toFrequencyDomain()} to get a visual + * representation of the frequency content of an image. + */ + public static ImageData ofPowerSpectrum(final Access2D transformed) { + + PhysicalStore magnitudes = Primitive64Store.getComplexModulus(transformed); + + magnitudes.modifyAll(PrimitiveMath.POWER.parameter(2).andThen(PrimitiveMath.LOG10)); + + double average = magnitudes.aggregateAll(Aggregator.AVERAGE).doubleValue(); + + magnitudes.modifyAll(PrimitiveMath.MULTIPLY.by(127.5D / average)); + + return ImageData.copy(magnitudes); + } + public static ImageData read(final File file) { try { return new ImageData(ImageIO.read(file)); @@ -126,10 +214,37 @@ public static ImageData read(final File file) { } } + public static ImageData read(final File directory, final String fileName) { + return ImageData.read(new File(directory, fileName)); + } + public static ImageData wrap(final BufferedImage image) { return new ImageData(image); } + private static Kernel getGaussianKernel(final double sigma, final int radius) { + int size = radius * 2 + 1; + float[] data = new float[size]; + + double twoSigmaSquared = 2.0 * sigma * sigma; + double sigmaRoot = Math.sqrt(twoSigmaSquared * Math.PI); + double total = 0.0; + + for (int i = -radius; i <= radius; i++) { + double distance = i * i; + int index = i + radius; + data[index] = (float) (Math.exp(-distance / twoSigmaSquared) / sigmaRoot); + total += data[index]; + } + + // Normalize the kernel + for (int i = 0; i < data.length; i++) { + data[i] /= total; + } + + return new Kernel(size, 1, data); + } + static int toRanged(final int value) { return Math.max(0, Math.min(value, 255)); } @@ -141,12 +256,37 @@ static int toRanged(final int value) { myImage = image; } + /** + * Creates a new image (of the same type) blurring the input using a Gaussian blur kernel. + * + * @param sigma The standard deviation of the Gaussian blur kernel + */ + public ImageData applyGaussianBlur(final double sigma) { + + int radius = (int) (3.0 * sigma); + if (radius % 2 == 0) { + radius++; // Ensure radius is odd + } + + Kernel kernel = ImageData.getGaussianKernel(sigma, radius); + ConvolveOp op = new ConvolveOp(kernel, ConvolveOp.EDGE_NO_OP, null); + + BufferedImage blurredImage = new BufferedImage(myImage.getWidth(), myImage.getHeight(), myImage.getType()); + Graphics2D g2d = blurredImage.createGraphics(); + g2d.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY); + g2d.drawImage(myImage, op, 0, 0); + g2d.dispose(); + + return new ImageData(blurredImage); + } + @Override public byte byteValue(final int row, final int col) { return (byte) this.intValue(row, col); } /** + * Creates a new image, converting the input to the specified image type in the process. * {@link BufferedImage#TYPE_BYTE_GRAY}, {@link BufferedImage#TYPE_INT_ARGB} or any other valid type... */ public ImageData convertTo(final int imageType) { @@ -207,6 +347,13 @@ public int getColDim() { return myImage.getWidth(); } + /** + * @return The underlying {@link BufferedImage} + */ + public BufferedImage getImage() { + return myImage; + } + @Override public int getRowDim() { return myImage.getHeight(); @@ -337,6 +484,20 @@ public ImageData sliceRedChannel() { return new SingleChannel(myImage, MASK_RED, SHIFT_RED); } + /** + * Transforms the spatial representation of the image to its frequency representation using the discrete + * Fourier transform. + *

+ * In addition the frequency representation is shifted so that the zero frequency is in the centre of the + * image. + *

+ * When you're done processing the frequency representation, you can transform it back to a spatial + * representation using {@link #fromFrequencyDomain(MatrixStore)}. + */ + public PhysicalStore toFrequencyDomain() { + return DiscreteFourierTransform.shift(DiscreteFourierTransform.transform2D(this)).copy(); + } + /** * The file format is derived from the file name ending (png, jpg...) */ @@ -350,6 +511,10 @@ public void writeTo(final File file) { } } + public void writeTo(final File directory, final String fileName) { + this.writeTo(new File(directory, fileName)); + } + private double doubleValue(final double row, final double col) { return this.doubleValue(MissingMath.roundToInt(row), MissingMath.roundToInt(col)); } diff --git a/src/main/java/org/ojalgo/data/transform/DataTransform.java b/src/main/java/org/ojalgo/data/transform/DataTransform.java new file mode 100644 index 000000000..20efd17c1 --- /dev/null +++ b/src/main/java/org/ojalgo/data/transform/DataTransform.java @@ -0,0 +1,33 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +/** + * DataTransform is a functional interface that can be used to transform data. The transform method takes an + * input of type IN and returns an output of type OUT. + */ +@FunctionalInterface +public interface DataTransform { + + OUT transform(IN input); + +} diff --git a/src/main/java/org/ojalgo/data/transform/DiscreteFourierTransform.java b/src/main/java/org/ojalgo/data/transform/DiscreteFourierTransform.java new file mode 100644 index 000000000..7fe8e5b5f --- /dev/null +++ b/src/main/java/org/ojalgo/data/transform/DiscreteFourierTransform.java @@ -0,0 +1,803 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +import java.util.function.DoubleUnaryOperator; + +import org.ojalgo.array.ArrayR064; +import org.ojalgo.function.PrimitiveFunction; +import org.ojalgo.function.constant.ComplexMath; +import org.ojalgo.function.constant.PrimitiveMath; +import org.ojalgo.function.series.PeriodicFunction; +import org.ojalgo.function.special.PowerOf2; +import org.ojalgo.matrix.MatrixC128; +import org.ojalgo.matrix.MatrixC128.DenseReceiver; +import org.ojalgo.matrix.store.GenericStore; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.matrix.store.PhysicalStore; +import org.ojalgo.matrix.store.TransformableRegion; +import org.ojalgo.scalar.ComplexNumber; +import org.ojalgo.structure.Access1D; +import org.ojalgo.structure.Access2D.ColumnView; +import org.ojalgo.structure.Access2D.RowView; +import org.ojalgo.structure.Factory2D; +import org.ojalgo.structure.Mutate2D; +import org.ojalgo.structure.Structure1D; + +/** + * The discrete Fourier transform (DFT) converts a finite sequence of equally-spaced samples of a function + * into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which + * is a complex-valued function of frequency. The interval at which the DTFT is sampled is the reciprocal of + * the duration of the input sequence. An inverse DFT is a Fourier series, using the DTFT samples as + * coefficients of complex sinusoids at the corresponding DTFT frequencies. The DFT is therefore said to be a + * frequency domain representation of the original input sequence. + *

+ * The fast Fourier transform (FFT) is an algorithm for computing the DFT; it achieves its high speed by + * storing and reusing results of computations as it progresses. + *

+ * Calling the factory method {@linkplain #newInstance(int)} will return an FFT implementation if the size is + * a power of 2, and a full matrix implementation otherwise. + */ +public abstract class DiscreteFourierTransform implements DataTransform, MatrixStore> { + + public static final class Directive { + + /** + * Assume input complex (if false, will assume NOT complex and run simpler code). + */ + public final boolean complex; + /** + * Conjugate input before transforming and the output after (before returning). This is how the + * inverse transform is performed. + */ + public final boolean conjugate; + /** + * Scale the output by the number of elements. Either the transform or the inverse transform needs to + * be scaled. Usually it's done on the inverse. + */ + public final boolean scale; + + public Directive(final boolean complex, final boolean conjugate, final boolean scale) { + super(); + this.complex = complex; + this.conjugate = conjugate; + this.scale = scale; + } + + Directive withComplex(final boolean complex) { + return new Directive(complex, conjugate, scale); + } + + } + + static final class FFT extends DiscreteFourierTransform { + + /** + * Perform the remaining stage calculations (stage>=1). This is essentially "the algorithm". + */ + private static void doStages(final int nbStages, final ComplexNumber[] roots, final double[] workRe, final double[] workIm) { + + int size = roots.length; + + int nbHalf = 2; // Half the number of lanes per set + int nbLanesPerSet = 4; + int nbSets = size / 4; + + int index1, index2; + + for (int stage = 2; stage < nbStages; stage++) { + + nbHalf = nbLanesPerSet; + nbLanesPerSet *= 2; + nbSets /= 2; + + for (int set = 0; set < nbSets; set++) { + + index1 = set * nbLanesPerSet; + index2 = index1 + nbHalf; + + FFT.update(workRe, workIm, index1, index2); + + for (int lane = 1; lane < nbHalf; lane++) { + + index1++; + index2++; + + FFT.update(workRe, workIm, index1, index2, roots[lane * nbSets]); + } + } + } + } + + private static void setup0(final Access1D input, final boolean complex, final boolean conjugate, final double[] workRe, final double[] workIm) { + + if (complex) { + ComplexNumber value = ComplexNumber.valueOf(input.get(0)); + if (conjugate) { + workRe[0] = value.doubleValue(); + workIm[0] = -value.i; + } else { + workRe[0] = value.doubleValue(); + workIm[0] = value.i; + } + } else { + workRe[0] = input.doubleValue(0); + workIm[0] = PrimitiveMath.ZERO; + } + } + + private static void setup0(final double[] input, final double[] workRe, final double[] workIm) { + workRe[0] = input[0]; + workIm[0] = PrimitiveMath.ZERO; + } + + private static void setup1(final Access1D input, final boolean complex, final boolean conjugate, final double[] workRe, final double[] workIm) { + + if (complex) { + + ComplexNumber value1 = ComplexNumber.valueOf(input.get(0)); + ComplexNumber value2 = ComplexNumber.valueOf(input.get(1)); + + FFT.toWork(0, 1, value1.doubleValue(), conjugate ? -value1.i : value1.i, value2.doubleValue(), conjugate ? -value2.i : value2.i, workRe, + workIm); + + } else { + + double re1 = input.doubleValue(0); + double re2 = input.doubleValue(1); + + workRe[0] = re1 + re2; + workRe[1] = re1 - re2; + + workIm[0] = PrimitiveMath.ZERO; + workIm[1] = PrimitiveMath.ZERO; + } + } + + private static void setup1(final double[] input, final double[] workRe, final double[] workIm) { + + double re1 = input[0]; + double re2 = input[1]; + + workRe[0] = re1 + re2; + workRe[1] = re1 - re2; + + workIm[0] = PrimitiveMath.ZERO; + workIm[1] = PrimitiveMath.ZERO; + } + + private static void setup2(final Access1D input, final boolean complex, final boolean conjugate, final int[] reversed, final double[] workRe, + final double[] workIm) { + + double re1, re2, re3, re4; + double re01, re02, re03, re04; + + if (complex) { + + double im1, im2, im3, im4; + double im01, im02, im03, im04; + + for (int index1 = 0, index2, index3, index4, size = reversed.length; index1 < size; index1 += 4) { + index2 = index1 + 1; + index3 = index1 + 2; + index4 = index1 + 3; + + ComplexNumber value1 = ComplexNumber.valueOf(input.get(reversed[index1])); + ComplexNumber value2 = ComplexNumber.valueOf(input.get(reversed[index2])); + ComplexNumber value3 = ComplexNumber.valueOf(input.get(reversed[index3])); + ComplexNumber value4 = ComplexNumber.valueOf(input.get(reversed[index4])); + + if (conjugate) { + + re1 = value1.doubleValue(); + im1 = -value1.i; + + re2 = value2.doubleValue(); + im2 = -value2.i; + + re3 = value3.doubleValue(); + im3 = -value3.i; + + re4 = value4.doubleValue(); + im4 = -value4.i; + + } else { + + re1 = value1.doubleValue(); + im1 = value1.i; + + re2 = value2.doubleValue(); + im2 = value2.i; + + re3 = value3.doubleValue(); + im3 = value3.i; + + re4 = value4.doubleValue(); + im4 = value4.i; + } + + re01 = re1 + re2; + im01 = im1 + im2; + + re02 = re1 - re2; + im02 = im1 - im2; + + re03 = re3 + re4; + im03 = im3 + im4; + + re04 = re3 - re4; + im04 = im3 - im4; + + workRe[index1] = re01 + re03; + workRe[index2] = re02 + im04; + workRe[index3] = re01 - re03; + workRe[index4] = re02 - im04; + + workIm[index1] = im01 + im03; + workIm[index2] = im02 - re04; + workIm[index3] = im01 - im03; + workIm[index4] = im02 + re04; + } + + } else { + + for (int index1 = 0, index2, index3, index4, size = reversed.length; index1 < size; index1 += 4) { + index2 = index1 + 1; + index3 = index1 + 2; + index4 = index1 + 3; + + re1 = input.doubleValue(reversed[index1]); + re2 = input.doubleValue(reversed[index2]); + re3 = input.doubleValue(reversed[index3]); + re4 = input.doubleValue(reversed[index4]); + + re01 = re1 + re2; + re02 = re1 - re2; + re03 = re3 + re4; + re04 = re3 - re4; + + workRe[index1] = re01 + re03; + workRe[index2] = re02; + workRe[index3] = re01 - re03; + workRe[index4] = re02; + + workIm[index1] = PrimitiveMath.ZERO; + workIm[index2] = -re04; + workIm[index3] = PrimitiveMath.ZERO; + workIm[index4] = re04; + } + } + } + + private static void setup2(final double[] input, final int[] reversed, final double[] workRe, final double[] workIm) { + + for (int index1 = 0, index2, index3, index4, size = reversed.length; index1 < size; index1 += 4) { + index2 = index1 + 1; + index3 = index1 + 2; + index4 = index1 + 3; + + double re1 = input[reversed[index1]]; + double re2 = input[reversed[index2]]; + double re3 = input[reversed[index3]]; + double re4 = input[reversed[index4]]; + + double re01 = re1 + re2; + double re02 = re1 - re2; + double re03 = re3 + re4; + double re04 = re3 - re4; + + workRe[index1] = re01 + re03; + workRe[index2] = re02; + workRe[index3] = re01 - re03; + workRe[index4] = re02; + + workIm[index1] = PrimitiveMath.ZERO; + workIm[index2] = -re04; + workIm[index3] = PrimitiveMath.ZERO; + workIm[index4] = re04; + } + } + + /** + * Copy the results to the output data structure. In the copy-process transformations are performed. + */ + private static void toOutput(final double[] workRe, final double[] workIm, final boolean conjugate, final boolean scale, + final Mutate2D.ModifiableReceiver output) { + + int size = workRe.length; + + if (conjugate) { + if (scale) { + double divisor = size; + for (int i = 0; i < size; i++) { + output.set(i, ComplexNumber.of(workRe[i] / divisor, -workIm[i] / divisor)); + } + } else { + for (int i = 0; i < size; i++) { + output.set(i, ComplexNumber.of(workRe[i], -workIm[i])); + } + } + } else { + if (scale) { + double divisor = size; + for (int i = 0; i < size; i++) { + output.set(i, ComplexNumber.of(workRe[i] / divisor, workIm[i] / divisor)); + } + } else { + for (int i = 0; i < size; i++) { + output.set(i, ComplexNumber.of(workRe[i], workIm[i])); + } + } + } + } + + private static void toWork(final int index1, final int index2, final double re1, final double im1, final double re2, final double im2, + final double[] workRe, final double[] workIm) { + + workRe[index1] = re1 + re2; + workRe[index2] = re1 - re2; + + workIm[index1] = im1 + im2; + workIm[index2] = im1 - im2; + } + + private static void update(final double[] workRe, final double[] workIm, final int index1, final int index2) { + FFT.toWork(index1, index2, workRe[index1], workIm[index1], workRe[index2], workIm[index2], workRe, workIm); + } + + private static void update(final double[] workRe, final double[] workIm, final int index1, final int index2, final ComplexNumber scalar) { + + double re2 = workRe[index2]; + double im2 = workIm[index2]; + + FFT.toWork(index1, index2, workRe[index1], workIm[index1], scalar.multiplyRe(re2, im2), scalar.multiplyIm(re2, im2), workRe, workIm); + } + + private final int[] myBitReversedIndices; + private final int myStages; + private final ComplexNumber[] myUnitRoots; + private final double[] myWorkIm; + private final double[] myWorkRe; + + FFT(final int size) { + + super(size); + + myStages = DiscreteFourierTransform.toPowerOf2Exponent(size); + if (myStages < 0) { + throw new IllegalArgumentException(); + } + + myBitReversedIndices = DiscreteFourierTransform.lookupIndices(size); + myUnitRoots = DiscreteFourierTransform.lookupRoots(size); + + myWorkRe = new double[size]; + myWorkIm = new double[size]; + } + + @Override + public void transform(final Access1D input, final Directive directive, final Mutate2D.ModifiableReceiver output) { + + if (myStages == 0) { + FFT.setup0(input, directive.complex, directive.conjugate, myWorkRe, myWorkIm); + } else if (myStages == 1) { + FFT.setup1(input, directive.complex, directive.conjugate, myWorkRe, myWorkIm); + } else if (myStages == 2) { + FFT.setup2(input, directive.complex, directive.conjugate, myBitReversedIndices, myWorkRe, myWorkIm); + } else { + FFT.setup2(input, directive.complex, directive.conjugate, myBitReversedIndices, myWorkRe, myWorkIm); + FFT.doStages(myStages, myUnitRoots, myWorkRe, myWorkIm); + } + + FFT.toOutput(myWorkRe, myWorkIm, directive.conjugate, directive.scale, output); + } + + @Override + public MatrixStore transform(final double... input) { + + if (myStages == 0) { + FFT.setup0(input, myWorkRe, myWorkIm); + } else if (myStages == 1) { + FFT.setup1(input, myWorkRe, myWorkIm); + } else if (myStages == 2) { + FFT.setup2(input, myBitReversedIndices, myWorkRe, myWorkIm); + } else { + FFT.setup2(input, myBitReversedIndices, myWorkRe, myWorkIm); + FFT.doStages(myStages, myUnitRoots, myWorkRe, myWorkIm); + } + + PhysicalStore output = GenericStore.C128.makeDense(input.length, 1); + + FFT.toOutput(myWorkRe, myWorkIm, DEFAULT.conjugate, DEFAULT.scale, output); + + return output; + } + + } + + static final class FullMatrix extends DiscreteFourierTransform { + + private final ComplexNumber myDivisor; + private final PhysicalStore myVandermondeMatrix; + + FullMatrix(final int size) { + + super(size); + + myVandermondeMatrix = GenericStore.C128.make(size, size); + DiscreteFourierTransform.generate(myVandermondeMatrix); + + myDivisor = ComplexNumber.valueOf(size); + } + + @Override + public void transform(final Access1D input, final Directive directive, final Mutate2D.ModifiableReceiver output) { + + MatrixStore column = GenericStore.C128.makeWrapperColumn(input); + + if (directive.conjugate) { + column = column.onAll(ComplexMath.CONJUGATE); + } + + myVandermondeMatrix.multiply(column, TransformableRegion.cast(output)); + + if (directive.conjugate) { + output.modifyAll(ComplexMath.CONJUGATE); + } + + if (directive.scale) { + output.modifyAll(ComplexMath.DIVIDE.by(myDivisor)); + } + } + + } + + static final class Single extends DiscreteFourierTransform { + + Single() { + super(1); + } + + @Override + public void transform(final Access1D input, final Directive directive, final Mutate2D.ModifiableReceiver output) { + output.fillMatching(input); + } + + } + + private static final int[][] BIT_REVERSED_INDICES = new int[31][]; + private static final ComplexNumber[][] UNIT_ROOTS = new ComplexNumber[31][]; + + static final Directive DEFAULT = new Directive(false, false, false); + static final Directive INVERSE = new Directive(true, true, true); + + public static int[] getBitReversedIndices(final int size) { + return DiscreteFourierTransform.lookupIndices(size).clone(); + } + + public static Access1D getUnitRoots(final int size) { + return Access1D.wrap(DiscreteFourierTransform.lookupRoots(size)); + } + + public static MatrixStore inverse2D(final MatrixStore input) { + + PhysicalStore retVal = GenericStore.C128.makeDense(input.getRowDim(), input.getColDim()); + + DiscreteFourierTransform.transform2D(input, INVERSE, retVal); + + return retVal; + } + + /** + * Will return a functioning instance for any size, but if you want fast transformation (FFT) it needs to + * be a power of 2. + */ + public static DiscreteFourierTransform newInstance(final int size) { + + if (size < 1) { + throw new IllegalArgumentException(); + } else if (size == 1) { + return new Single(); + } else if (PowerOf2.isPowerOf2(size)) { + return new FFT(size); + } else { + return new FullMatrix(size); + } + } + + public static M newVandermonde(final Factory2D factory, final int size) { + + M matrix = factory.make(size, size); + + DiscreteFourierTransform.generate(matrix); + + return matrix; + } + + public static MatrixC128 newVandermondeMatrix(final int size) { + + DenseReceiver receiver = MatrixC128.FACTORY.makeDense(size, size); + + DiscreteFourierTransform.generate(receiver); + + return receiver.get(); + } + + /** + * Sample, and transform, a function using the Discrete Fourier Transform. + */ + public static MatrixStore sample(final DoubleUnaryOperator function, final PrimitiveFunction.SampleDomain sampleDomain) { + + double[] input = sampleDomain.arguments(); + for (int i = 0; i < input.length; i++) { + input[i] = function.applyAsDouble(input[i]); + } + + return DiscreteFourierTransform.newInstance(input.length).transform(input); + } + + /** + * @see #sample(DoubleUnaryOperator, PrimitiveFunction.SampleDomain) + */ + public static MatrixStore sample(final PeriodicFunction function, final int nbSamples) { + return DiscreteFourierTransform.sample(function, function.getSampleDomain(nbSamples)); + } + + /** + * There is a symmetry in the DFT matrix. The first half of the rows (and columns) are the complex + * conjugates of the second half. Furthermore, the first half correspond to positive frequencies, and the + * second half to negative frequencies. + *

+ * Re-arranging the elements of a matrix, shifting the first and second halves of the rows (and columns), + * puts the zero-frequency term in the middle, and the conjugate pairs at equal distances from the centre. + *

+ * This is useful when displaying the 2D DFT matrix as an image. + *

+ * To revert the shift, simply call {@linkplain #shift(MatrixStore)} again. However, this will not work + * correctly if there's an odd number of rows or columns. In that case the second call will not correctly + * revert the position of the zero-frequency term – it will end up in the last row/column instead of in + * the first. + */ + public static > MatrixStore shift(final MatrixStore matrix) { + + MatrixStore retVal = matrix; + + int nbRows = matrix.getRowDim(); + int nbCols = matrix.getColDim(); + + if (nbRows > 1) { + + int half = (nbRows + 1) / 2; + + MatrixStore first = retVal.limits(half, -1); + MatrixStore second = retVal.offsets(half, 0); + + retVal = second.below(first); + } + + if (nbCols > 1) { + + int half = (nbCols + 1) / 2; + + MatrixStore first = retVal.limits(-1, half); + MatrixStore second = retVal.offsets(0, half); + + retVal = second.right(first); + } + + return retVal; + } + + /** + * Perform a 2D Discrete Fourier Transform on the input matrix. The output will be a complex matrix of the + * same size. + */ + public static MatrixStore transform2D(final MatrixStore input) { + + PhysicalStore retVal = GenericStore.C128.makeDense(input.getRowDim(), input.getColDim()); + + DiscreteFourierTransform.transform2D(input, DEFAULT, retVal); + + return retVal; + } + + public static void transform2D(final MatrixStore input, final Directive directive, final TransformableRegion output) { + + int nbRows = input.getRowDim(); + int nbCols = input.getColDim(); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(nbRows); + PhysicalStore workOutput = GenericStore.C128.makeDense(nbRows, 1); + + Directive directive1 = directive; + for (ColumnView view : input.columns()) { + transformer.transform(view, directive1, workOutput); + output.fillColumn(view.column(), workOutput); + } + + if (nbRows != nbCols) { + transformer = DiscreteFourierTransform.newInstance(nbCols); + workOutput = GenericStore.C128.makeDense(nbCols, 1); + } + + Directive directive2 = directive.withComplex(true); + for (RowView view : output.rows()) { + transformer.transform(view, directive2, workOutput); + output.fillRow(view.row(), workOutput); + } + } + + private static ComplexNumber[] lookupRootsExponent(final int exponent) { + + if (exponent >= UNIT_ROOTS.length) { + + throw new IllegalArgumentException(); + + } else if (UNIT_ROOTS[exponent] != null) { + + return UNIT_ROOTS[exponent]; + + } else { + + if (exponent == 0) { + + return UNIT_ROOTS[0] = new ComplexNumber[] { ComplexNumber.ONE }; + + } else if (exponent == 1) { + + return UNIT_ROOTS[1] = new ComplexNumber[] { ComplexNumber.ONE, ComplexNumber.NEG }; + + } else if (exponent == 2) { + + return UNIT_ROOTS[2] = new ComplexNumber[] { ComplexNumber.ONE, ComplexNumber.N, ComplexNumber.NEG, ComplexNumber.I }; + + } else { + + ComplexNumber[] half = DiscreteFourierTransform.lookupRootsExponent(exponent - 1); + + ComplexNumber[] full = UNIT_ROOTS[exponent] = new ComplexNumber[half.length + half.length]; + + ComplexNumber[] unitRoots = ComplexNumber.newUnitRoots(full.length); + + for (int i = 0; i < half.length; i++) { + int ii = 2 * i; + int ii1 = ii + 1; + + full[ii] = half[i]; + full[ii1] = unitRoots[ii1]; + } + + return full; + } + } + } + + /** + * Function to perform bit-reversal on indices + */ + private static void reverseBits(final int[] indices) { + + int n = indices.length; + int bits = 31 - Integer.numberOfLeadingZeros(n); + + for (int i = 0; i < n; i++) { + int reversedIndex = 0; + for (int j = 0; j < bits; j++) { + reversedIndex |= (i >> j & 1) << bits - 1 - j; + } + // Swap indices[i] with indices[reversedIndex] + if (reversedIndex > i) { + int temp = indices[i]; + indices[i] = indices[reversedIndex]; + indices[reversedIndex] = temp; + } + } + } + + static void generate(final Mutate2D matrix) { + + int size = matrix.getMinDim(); + + double unitRootPhase = ComplexNumber.newUnitRoot(size).phase(); + + for (int j = 0; j < size; j++) { + for (int i = 0; i < size; i++) { + matrix.set(i, j, ComplexNumber.makePolar(PrimitiveMath.ONE, i * j * unitRootPhase)); + } + } + } + + static int[] lookupIndices(final int size) { + + int exponent = DiscreteFourierTransform.toPowerOf2Exponent(size); + + int[] retVal = BIT_REVERSED_INDICES[exponent]; + + if (retVal == null) { + + BIT_REVERSED_INDICES[exponent] = retVal = Structure1D.newIncreasingRange(0, size); + + DiscreteFourierTransform.reverseBits(retVal); + } + + return retVal; + } + + static ComplexNumber[] lookupRoots(final int size) { + + int exponent = DiscreteFourierTransform.toPowerOf2Exponent(size); + + return DiscreteFourierTransform.lookupRootsExponent(exponent); + } + + static int toPowerOf2Exponent(final int size) { + + int exponent = PowerOf2.find(size); + + if (exponent < 0) { + throw new IllegalArgumentException("Needs to be power of 2!"); + } + + return exponent; + } + + private final int mySize; + + DiscreteFourierTransform(final int size) { + super(); + mySize = size; + } + + public final void inverse(final Access1D input, final Mutate2D.ModifiableReceiver output) { + this.transform(input, INVERSE, output); + } + + public final MatrixStore inverse(final Access1D input) { + PhysicalStore output = GenericStore.C128.makeDense(input.size(), 1); + this.transform(input, INVERSE, output); + return output; + } + + @Override + public final MatrixStore transform(final Access1D input) { + PhysicalStore output = GenericStore.C128.makeDense(input.size(), 1); + this.transform(input, DEFAULT, output); + return output; + } + + public abstract void transform(Access1D input, Directive directive, Mutate2D.ModifiableReceiver output); + + public final void transform(final Access1D input, final Mutate2D.ModifiableReceiver output) { + this.transform(input, DEFAULT, output); + } + + public MatrixStore transform(final double... input) { + PhysicalStore output = GenericStore.C128.makeDense(input.length, 1); + this.transform(ArrayR064.wrap(input), DEFAULT, output); + return output; + } + + int size() { + return mySize; + } + +} diff --git a/src/main/java/org/ojalgo/data/transform/ZTransform.java b/src/main/java/org/ojalgo/data/transform/ZTransform.java new file mode 100644 index 000000000..ce7ebed6d --- /dev/null +++ b/src/main/java/org/ojalgo/data/transform/ZTransform.java @@ -0,0 +1,153 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +import java.util.function.UnaryOperator; + +import org.ojalgo.function.constant.PrimitiveMath; +import org.ojalgo.matrix.store.GenericStore; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.scalar.ComplexNumber; +import org.ojalgo.structure.Access1D; + +/** + * This class implements the Z-transform for a given complex number z. The transform is applied to a sequence + * of real numbers, and the output is a complex number. The Z-transform is defined as: Z{x[n]} = ∑x[n]z⁻ⁿ + * where n is the index of the sequence, x[n] is the value of the sequence at index n, and z is the complex + * number that the transform is applied to. + *

+ * Choosing the complex number z in the context of the Z-transform is often application-specific and depends + * on the goals of the analysis or the properties you want to study. There are some common choices and + * guidelines for specific scenarios: + *

    + *
  1. Unit Circle (|z| = 1): In many applications, especially in stability analysis and frequency response, z + * is chosen to lie on the unit circle (|z| = 1). This is because the unit circle in the complex plane + * corresponds to the frequency response of a system, and properties of a system can be analyzed by studying + * the behavior of the transfer function along the unit circle. + *
  2. Stability Analysis: For stability analysis of discrete-time systems, z is often chosen such that all + * poles of the system lie inside the unit circle. This ensures that the system is stable. + *
  3. Frequency Analysis: If you are interested in the frequency domain behavior of a system, you might + * choose z to be a complex exponential, e.g., z = exp(jω) where j is the imaginary unit and ω is the angular + * frequency. + *
  4. Arbitrary Complex Numbers: In some cases, you might be interested in the Z-transform at an arbitrary + * complex number z for general analysis. This could be done to study the behavior of the system at a specific + * point in the complex plane. + *
+ */ +public final class ZTransform implements DataTransform, ComplexNumber> { + + static final class ZOperator implements UnaryOperator { + + private final Access1D mySequence; + + ZOperator(final Access1D sequence) { + super(); + mySequence = sequence; + } + + @Override + public ComplexNumber apply(final ComplexNumber z) { + return ZTransform.doTransform(mySequence, z); + } + + } + + /** + * This method computes the discrete Fourier transform (DFT) of a sequence of real numbers. The DFT is + * computed using the Z-transform with z = exp(jω) where j is the imaginary unit and ω is the angular + * frequency. The DFT is defined as: X[k] = ∑x[n]exp(-jωkn) where n is the index of the sequence, x[n] is + * the value of the sequence at index n, and X[k] is the value of the DFT at index k. + *

+ * This method exists primarily for testing purposes. It is most likely better to use + * {@link DiscreteFourierTransform} to compute the DFT using an FFT algorithm. + */ + public static MatrixStore doDFT(final Access1D input) { + + UnaryOperator operator = ZTransform.newZOperator(input); + + int size = input.size(); + + GenericStore retVal = GenericStore.C128.make(size, 1); + + double angularFrequencyFactor = PrimitiveMath.TWO_PI / size; + + for (int i = 0; i < size; i++) { + retVal.set(i, 0, operator.apply(ComplexNumber.makePolar(PrimitiveMath.ONE, i * angularFrequencyFactor))); + } + + return retVal; + } + + public static UnaryOperator newZOperator(final Access1D sequence) { + return new ZOperator(sequence); + } + + public static ZTransform of(final ComplexNumber z) { + return new ZTransform(z); + } + + public static ZTransform of(final double angularFrequency) { + return new ZTransform(ComplexNumber.makePolar(PrimitiveMath.ONE, angularFrequency)); + } + + static ComplexNumber doTransform(final Access1D input, final ComplexNumber z) { + + ComplexNumber retVal = ComplexNumber.ZERO; + + for (int n = 0, N = input.size(); n < N; n++) { + retVal = retVal.add(z.power(-n).multiply(input.doubleValue(n))); + } + + return retVal; + } + + private final double myNormZ; + private final double myPhaseZ; + private final ComplexNumber myZ; + + ZTransform(final ComplexNumber z) { + + super(); + + myZ = z; + myNormZ = z.norm(); + myPhaseZ = z.phase(); + } + + /** + * Input is a sequence of real numbers. Output is a complex number. + */ + @Override + public ComplexNumber transform(final Access1D input) { + return ZTransform.doTransform(input, myZ); + + } + + double getNormZ() { + return myNormZ; + } + + double getPhaseZ() { + return myPhaseZ; + } + +} diff --git a/src/main/java/org/ojalgo/function/PrimitiveFunction.java b/src/main/java/org/ojalgo/function/PrimitiveFunction.java index 126f74dc7..b1aa21d02 100644 --- a/src/main/java/org/ojalgo/function/PrimitiveFunction.java +++ b/src/main/java/org/ojalgo/function/PrimitiveFunction.java @@ -24,6 +24,7 @@ import org.ojalgo.function.aggregator.AggregatorSet; import org.ojalgo.function.aggregator.PrimitiveAggregator; import org.ojalgo.function.constant.PrimitiveMath; +import org.ojalgo.function.special.PowerOf2; import org.ojalgo.type.context.NumberContext; /** @@ -39,19 +40,22 @@ public final class PrimitiveFunction extends FunctionSet { @FunctionalInterface public interface Binary extends BinaryFunction { - default Double invoke(final Double arg1, final Double arg2) { - return Double.valueOf(this.invoke(arg1.doubleValue(), arg2.doubleValue())); - } - + @Override default Double invoke(final Double arg1, final double arg2) { return Double.valueOf(this.invoke(arg1.doubleValue(), arg2)); } + @Override + default Double invoke(final Double arg1, final Double arg2) { + return Double.valueOf(this.invoke(arg1.doubleValue(), arg2.doubleValue())); + } + } @FunctionalInterface public interface Consumer extends VoidFunction { + @Override default void invoke(final Double arg) { this.invoke(arg.doubleValue()); } @@ -61,6 +65,7 @@ default void invoke(final Double arg) { @FunctionalInterface public interface Parameter extends ParameterFunction { + @Override default Double invoke(final Double arg, final int param) { return Double.valueOf(this.invoke(arg.doubleValue(), param)); } @@ -70,15 +75,64 @@ default Double invoke(final Double arg, final int param) { @FunctionalInterface public interface Predicate extends PredicateFunction { + @Override default boolean invoke(final Double arg) { return this.invoke(arg.doubleValue()); } } + public static final class SampleDomain { + + private final double myIncrement; + private final double myPeriod; + private final int myNumberOfSamples; + + public SampleDomain(final double period, final int nbSamples) { + super(); + myPeriod = period; + myNumberOfSamples = nbSamples; + myIncrement = period / nbSamples; + } + + /** + * Adjusts the number of samples to the smallest power of 2 that is not less than the current number + * of samples. + */ + public SampleDomain adjustToPowerOf2() { + return new SampleDomain(myPeriod, PowerOf2.smallestNotLessThan(myNumberOfSamples)); + } + + public double argumant(final int index) { + return index * myIncrement; + } + + public double[] arguments() { + double[] retVal = new double[myNumberOfSamples]; + for (int i = 0; i < myNumberOfSamples; i++) { + retVal[i] = i * myIncrement; + } + return retVal; + } + + public double increment() { + return myIncrement; + } + + public double period() { + return myPeriod; + } + + public int size() { + return myNumberOfSamples; + } + + } + @FunctionalInterface public interface Unary extends UnaryFunction { + @Override default Double invoke(final Double arg) { return Double.valueOf(this.invoke(arg.doubleValue())); } diff --git a/src/main/java/org/ojalgo/function/UnaryFunction.java b/src/main/java/org/ojalgo/function/UnaryFunction.java index 1fc87b059..ce9185f0a 100644 --- a/src/main/java/org/ojalgo/function/UnaryFunction.java +++ b/src/main/java/org/ojalgo/function/UnaryFunction.java @@ -37,7 +37,7 @@ static > boolean isZeroModified(final UnaryFunction f default UnaryFunction andThen(final UnaryFunction after) { ProgrammingError.throwIfNull(after); - return new UnaryFunction() { + return new UnaryFunction<>() { public double invoke(final double arg) { return after.invoke(UnaryFunction.this.invoke(arg)); @@ -54,26 +54,31 @@ public N invoke(final N arg) { }; } + @Override default N apply(final N arg) { return this.invoke(arg); } + @Override default double applyAsDouble(final double arg) { return this.invoke(arg); } default UnaryFunction compose(final UnaryFunction before) { ProgrammingError.throwIfNull(before); - return new UnaryFunction() { + return new UnaryFunction<>() { + @Override public double invoke(final double arg) { return UnaryFunction.this.invoke(before.invoke(arg)); } + @Override public float invoke(final float arg) { return UnaryFunction.this.invoke(before.invoke(arg)); } + @Override public N invoke(final N arg) { return UnaryFunction.this.invoke(before.invoke(arg)); } diff --git a/src/main/java/org/ojalgo/function/polynomial/AbstractPolynomial.java b/src/main/java/org/ojalgo/function/polynomial/AbstractPolynomial.java index c536a41da..0844c5b33 100644 --- a/src/main/java/org/ojalgo/function/polynomial/AbstractPolynomial.java +++ b/src/main/java/org/ojalgo/function/polynomial/AbstractPolynomial.java @@ -23,27 +23,31 @@ import java.math.BigDecimal; import java.util.List; +import java.util.Objects; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.BasicArray; +import org.ojalgo.data.transform.DiscreteFourierTransform; import org.ojalgo.function.constant.BigMath; +import org.ojalgo.function.constant.ComplexMath; +import org.ojalgo.function.special.PowerOf2; import org.ojalgo.matrix.decomposition.QR; +import org.ojalgo.matrix.store.GenericStore; import org.ojalgo.matrix.store.PhysicalStore; +import org.ojalgo.scalar.ComplexNumber; import org.ojalgo.series.NumberSeries; import org.ojalgo.structure.Access1D; import org.ojalgo.type.TypeUtils; +import org.ojalgo.type.context.NumberContext; -abstract class AbstractPolynomial> implements PolynomialFunction { +abstract class AbstractPolynomial, P extends AbstractPolynomial> implements PolynomialFunction { - private final Array1D myCoefficients; - private transient AbstractPolynomial myDerivative = null; - private transient AbstractPolynomial myPrimitive = null; + public static final NumberContext DEGREE_ACCURACY = NumberContext.of(16); - @SuppressWarnings("unused") - private AbstractPolynomial() { - this(null); - } + private final BasicArray myCoefficients; + private transient P myDerivative = null; + private transient P myPrimitive = null; - protected AbstractPolynomial(final Array1D coefficients) { + AbstractPolynomial(final BasicArray coefficients) { super(); @@ -51,13 +55,35 @@ protected AbstractPolynomial(final Array1D coefficients) { } @Override - public PolynomialFunction buildDerivative() { + public P add(final PolynomialFunction addend) { + + int leftDeg = this.degree(); + int righDeg = addend.degree(); + + int retSize = 1 + Math.max(leftDeg, righDeg); + + P retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + coefficients.add(l, this.get(l)); + } + + for (int r = 0; r <= righDeg; r++) { + coefficients.add(r, addend.get(r)); + } + + return retVal; + } + + @Override + public P buildDerivative() { if (myDerivative == null) { int tmpSize = Math.max(1, myCoefficients.size() - 1); - myDerivative = this.makeInstance(tmpSize); + myDerivative = this.newInstance(tmpSize); for (int i = 0; i < tmpSize; i++) { myDerivative.set(i, this.getDerivativeFactor(i)); @@ -68,13 +94,13 @@ public PolynomialFunction buildDerivative() { } @Override - public PolynomialFunction buildPrimitive() { + public P buildPrimitive() { if (myPrimitive == null) { int tmpSize = myCoefficients.size() + 1; - myPrimitive = this.makeInstance(tmpSize); + myPrimitive = this.newInstance(tmpSize); for (int i = 0; i < tmpSize; i++) { myPrimitive.set(i, this.getPrimitiveFactor(i)); @@ -86,12 +112,19 @@ public PolynomialFunction buildPrimitive() { @Override public long count() { - return this.size(); + return myCoefficients.count(); } @Override - public int degree() { - return myCoefficients.size() - 1; + public int degree(final NumberContext accuracy) { + + int retVal = myCoefficients.size() - 1; + + while (retVal > 0 && accuracy.isZero(this.norm(retVal))) { + retVal--; + } + + return retVal; } @Override @@ -100,34 +133,20 @@ public double doubleValue(final int power) { } @Override - public void estimate(final List x, final List y) { - this.estimate(Access1D.wrap(x), Access1D.wrap(y)); - } - - void estimate(final Access1D x, final Access1D y, final PhysicalStore.Factory store, final QR.Factory qr) { - - int tmpRowDim = Math.min(x.size(), y.size()); - int tmpColDim = this.size(); - - PhysicalStore tmpBody = store.make(tmpRowDim, tmpColDim); - PhysicalStore tmpRHS = store.make(tmpRowDim, 1); - - for (int i = 0; i < tmpRowDim; i++) { - - BigDecimal tmpX = BigMath.ONE; - BigDecimal tmpXfactor = TypeUtils.toBigDecimal(x.get(i)); - BigDecimal tmpY = TypeUtils.toBigDecimal(y.get(i)); - - for (int j = 0; j < tmpColDim; j++) { - tmpBody.set(i, j, tmpX); - tmpX = tmpX.multiply(tmpXfactor); - } - tmpRHS.set(i, 0, tmpY); + public boolean equals(final Object obj) { + if (this == obj) { + return true; + } + if (obj == null || this.getClass() != obj.getClass()) { + return false; } + AbstractPolynomial other = (AbstractPolynomial) obj; + return Objects.equals(myCoefficients, other.myCoefficients); + } - QR tmpQR = qr.make(tmpBody); - tmpQR.decompose(tmpBody); - this.set(tmpQR.getSolution(tmpRHS)); + @Override + public void estimate(final List x, final List y) { + this.estimate(Access1D.wrap(x), Access1D.wrap(y)); } @Override @@ -140,6 +159,11 @@ public N get(final long power) { return myCoefficients.get(power); } + @Override + public int hashCode() { + return Objects.hash(myCoefficients); + } + @Override public double invoke(final double arg) { @@ -168,6 +192,78 @@ public float invoke(final float arg) { return retVal; } + @Override + public abstract P multiply(final PolynomialFunction multiplicand); + + @Override + public final P power(final int power) { + + if (power < 0) { + + throw new IllegalArgumentException(); + + } else if (power == 0) { + + return this.one(); + + } else if (power == 1) { + + return (P) this; + + } else if (power == 2) { + + return this.multiply(this); + + } else if (power == 3) { + + P temp = this.multiply(this); + + return temp.multiply(this); + + } else if (power == 4) { + + P temp = this.multiply(this); + + return temp.multiply(temp); + + } else { + + int degree = this.degree(); + + int newSize = 1 + power * degree; + + int powerOf2 = PowerOf2.smallestNotLessThan(newSize); + + GenericStore work = GenericStore.C128.make(powerOf2, 1); + work.fillMatching(myCoefficients); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(powerOf2); + + transformer.transform(work, work); + + work.modifyAll(ComplexMath.POWER.parameter(power)); + + transformer.inverse(work, work); + + P retVal = this.newInstance(newSize); + + retVal.coefficients().fillMatching(work); + + return retVal; + } + } + + @Override + public final void set(final Access1D coefficients) { + if (coefficients.size() > myCoefficients.size()) { + throw new IllegalArgumentException(); + } + myCoefficients.reset(); + myCoefficients.fillMatching(coefficients); + myDerivative = null; + myPrimitive = null; + } + @Override public void set(final int power, final double coefficient) { myCoefficients.set(power, coefficient); @@ -182,15 +278,61 @@ public void set(final int power, final N coefficient) { myPrimitive = null; } + @Override + public final void set(final long power, final Comparable value) { + myCoefficients.set(power, value); + } + @Override public int size() { return myCoefficients.size(); } + @Override + public String toString() { + return myCoefficients.toString(); + } + protected abstract N getDerivativeFactor(int power); protected abstract N getPrimitiveFactor(int power); - protected abstract AbstractPolynomial makeInstance(int size); + protected abstract P newInstance(int size); + + BasicArray coefficients() { + return myCoefficients; + } + + void estimate(final Access1D x, final Access1D y, final PhysicalStore.Factory store, final QR.Factory qr) { + + int tmpRowDim = Math.min(x.size(), y.size()); + int tmpColDim = this.size(); + + PhysicalStore tmpBody = store.make(tmpRowDim, tmpColDim); + PhysicalStore tmpRHS = store.make(tmpRowDim, 1); + + for (int i = 0; i < tmpRowDim; i++) { + + BigDecimal tmpX = BigMath.ONE; + BigDecimal tmpXfactor = TypeUtils.toBigDecimal(x.get(i)); + BigDecimal tmpY = TypeUtils.toBigDecimal(y.get(i)); + + for (int j = 0; j < tmpColDim; j++) { + tmpBody.set(i, j, tmpX); + tmpX = tmpX.multiply(tmpXfactor); + } + tmpRHS.set(i, 0, tmpY); + } + + QR tmpQR = qr.make(tmpBody); + tmpQR.decompose(tmpBody); + this.set(tmpQR.getSolution(tmpRHS)); + } + + double norm(final int power) { + return Math.abs(this.doubleValue(power)); + } + + abstract P one(); } diff --git a/src/main/java/org/ojalgo/function/polynomial/BigPolynomial.java b/src/main/java/org/ojalgo/function/polynomial/BigPolynomial.java index ba54f6355..78d5b0478 100644 --- a/src/main/java/org/ojalgo/function/polynomial/BigPolynomial.java +++ b/src/main/java/org/ojalgo/function/polynomial/BigPolynomial.java @@ -23,7 +23,7 @@ import java.math.BigDecimal; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.BasicArray; /** * @deprecated v53 use {@link PolynomialR256} instead @@ -35,7 +35,7 @@ public BigPolynomial(final int degree) { super(degree); } - public BigPolynomial(final Array1D coefficients) { + public BigPolynomial(final BasicArray coefficients) { super(coefficients); } diff --git a/src/main/java/org/ojalgo/function/polynomial/ComplexPolynomial.java b/src/main/java/org/ojalgo/function/polynomial/ComplexPolynomial.java index a05514a86..76cb31774 100644 --- a/src/main/java/org/ojalgo/function/polynomial/ComplexPolynomial.java +++ b/src/main/java/org/ojalgo/function/polynomial/ComplexPolynomial.java @@ -21,7 +21,7 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.BasicArray; import org.ojalgo.scalar.ComplexNumber; /** @@ -34,7 +34,7 @@ public ComplexPolynomial(final int degree) { super(degree); } - public ComplexPolynomial(final Array1D coefficients) { + public ComplexPolynomial(final BasicArray coefficients) { super(coefficients); } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialC128.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialC128.java index 15dbe48ea..36da5d12b 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialC128.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialC128.java @@ -21,56 +21,34 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayC128; +import org.ojalgo.array.BasicArray; import org.ojalgo.matrix.decomposition.QR; import org.ojalgo.matrix.store.GenericStore; import org.ojalgo.scalar.ComplexNumber; import org.ojalgo.structure.Access1D; -public class PolynomialC128 extends AbstractPolynomial { +public class PolynomialC128 extends ScalarPolynomial { + + public static final PolynomialC128 ONE = PolynomialC128.wrap(ComplexNumber.ONE); + + public static PolynomialC128 wrap(final ComplexNumber... coefficients) { + return new PolynomialC128(ArrayC128.wrap(coefficients)); + } public PolynomialC128(final int degree) { - super(Array1D.C128.make(degree + 1)); + super(ArrayC128.make(degree + 1)); } - PolynomialC128(final Array1D coefficients) { + PolynomialC128(final BasicArray coefficients) { super(coefficients); } + @Override public void estimate(final Access1D x, final Access1D y) { this.estimate(x, y, GenericStore.C128, QR.C128); } - public ComplexNumber integrate(final ComplexNumber fromPoint, final ComplexNumber toPoint) { - - PolynomialFunction primitive = this.buildPrimitive(); - - ComplexNumber fromVal = primitive.invoke(fromPoint); - ComplexNumber toVal = primitive.invoke(toPoint); - - return toVal.subtract(fromVal); - } - - public ComplexNumber invoke(final ComplexNumber arg) { - - int power = this.degree(); - - ComplexNumber retVal = this.get(power); - - while (--power >= 0) { - retVal = this.get(power).add(arg.multiply(retVal)); - } - - return retVal; - } - - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, ComplexNumber.valueOf(coefficients.get(p))); - } - } - @Override protected ComplexNumber getDerivativeFactor(final int power) { int nextIndex = power + 1; @@ -86,8 +64,13 @@ protected ComplexNumber getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialC128(Array1D.C128.make(size)); + protected PolynomialC128 newInstance(final int size) { + return new PolynomialC128(ArrayC128.make(size)); + } + + @Override + PolynomialC128 one() { + return ONE; } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialFunction.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialFunction.java index f9e6023fd..65330d0e0 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialFunction.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialFunction.java @@ -23,16 +23,29 @@ import java.util.List; +import org.ojalgo.algebra.Ring; import org.ojalgo.function.BasicFunction.Differentiable; import org.ojalgo.function.BasicFunction.Integratable; import org.ojalgo.function.UnaryFunction; import org.ojalgo.series.NumberSeries; import org.ojalgo.structure.Access1D; +import org.ojalgo.structure.Mutate1D; +import org.ojalgo.type.context.NumberContext; -public interface PolynomialFunction> - extends UnaryFunction, Access1D, Differentiable>, Integratable> { +public interface PolynomialFunction> extends UnaryFunction, Access1D, Mutate1D, Differentiable>, + Integratable>, Ring> { - int degree(); + /** + * The largest exponent/power of the non-zero coefficients. + */ + default int degree() { + return this.degree(AbstractPolynomial.DEGREE_ACCURACY); + } + + /** + * The largest exponent/power of the non-zero (to the given accuracy) coefficients. + */ + int degree(NumberContext accuracy); void estimate(Access1D x, Access1D y); @@ -42,8 +55,6 @@ public interface PolynomialFunction> void set(Access1D coefficients); - void set(final int power, final double coefficient); - void set(final int power, final N coefficient); } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialQ128.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialQ128.java index d740118fc..461280ec4 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialQ128.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialQ128.java @@ -21,56 +21,34 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayQ128; +import org.ojalgo.array.BasicArray; import org.ojalgo.matrix.decomposition.QR; import org.ojalgo.matrix.store.GenericStore; import org.ojalgo.scalar.RationalNumber; import org.ojalgo.structure.Access1D; -public class PolynomialQ128 extends AbstractPolynomial { +public class PolynomialQ128 extends ScalarPolynomial { + + public static final PolynomialQ128 ONE = PolynomialQ128.wrap(RationalNumber.ONE); + + public static PolynomialQ128 wrap(final RationalNumber... coefficients) { + return new PolynomialQ128(ArrayQ128.wrap(coefficients)); + } public PolynomialQ128(final int degree) { - super(Array1D.Q128.make(degree + 1)); + super(ArrayQ128.make(degree + 1)); } - PolynomialQ128(final Array1D coefficients) { + PolynomialQ128(final BasicArray coefficients) { super(coefficients); } + @Override public void estimate(final Access1D x, final Access1D y) { this.estimate(x, y, GenericStore.Q128, QR.Q128); } - public RationalNumber integrate(final RationalNumber fromPoint, final RationalNumber toPoint) { - - PolynomialFunction primitive = this.buildPrimitive(); - - RationalNumber fromVal = primitive.invoke(fromPoint); - RationalNumber toVal = primitive.invoke(toPoint); - - return toVal.subtract(fromVal); - } - - public RationalNumber invoke(final RationalNumber arg) { - - int power = this.degree(); - - RationalNumber retVal = this.get(power); - - while (--power >= 0) { - retVal = this.get(power).add(arg.multiply(retVal)); - } - - return retVal; - } - - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, RationalNumber.valueOf(coefficients.get(p))); - } - } - @Override protected RationalNumber getDerivativeFactor(final int power) { int nextIndex = power + 1; @@ -86,8 +64,13 @@ protected RationalNumber getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialQ128(Array1D.Q128.make(size)); + protected PolynomialQ128 newInstance(final int size) { + return new PolynomialQ128(ArrayQ128.make(size)); + } + + @Override + PolynomialQ128 one() { + return ONE; } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialR032.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialR032.java index 9545d9ec7..ebfb9e060 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialR032.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialR032.java @@ -21,26 +21,57 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayR032; +import org.ojalgo.array.BasicArray; import org.ojalgo.function.constant.PrimitiveMath; import org.ojalgo.matrix.decomposition.QR; import org.ojalgo.matrix.store.Primitive32Store; import org.ojalgo.structure.Access1D; -public final class PolynomialR032 extends AbstractPolynomial { +public final class PolynomialR032 extends AbstractPolynomial { + + public static final PolynomialR032 ONE = PolynomialR032.wrap(1F); + + public static PolynomialR032 wrap(final float... coefficients) { + return new PolynomialR032(ArrayR032.wrap(coefficients)); + } public PolynomialR032(final int degree) { - super(Array1D.R032.make(degree + 1)); + super(ArrayR032.make(degree + 1)); } - PolynomialR032(final Array1D coefficients) { + PolynomialR032(final BasicArray coefficients) { super(coefficients); } + @Override + public PolynomialR032 add(final PolynomialFunction addend) { + + int leftDeg = this.degree(); + int righDeg = addend.degree(); + + int retSize = 1 + Math.max(leftDeg, righDeg); + + PolynomialR032 retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + coefficients.add(l, this.floatValue(l)); + } + + for (int r = 0; r <= righDeg; r++) { + coefficients.add(r, addend.floatValue(r)); + } + + return retVal; + } + + @Override public void estimate(final Access1D x, final Access1D y) { this.estimate(x, y, Primitive32Store.FACTORY, QR.R064); } + @Override public Double integrate(final Double fromPoint, final Double toPoint) { PolynomialFunction primitive = this.buildPrimitive(); @@ -51,15 +82,49 @@ public Double integrate(final Double fromPoint, final Double toPoint) { return Double.valueOf(toVal - fromVal); } + @Override public Double invoke(final Double arg) { return Double.valueOf(this.invoke(arg.doubleValue())); } - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, coefficients.doubleValue(p)); + @Override + public PolynomialR032 multiply(final PolynomialFunction multiplicand) { + + int leftDeg = this.degree(); + int righDeg = multiplicand.degree(); + + int retSize = 1 + leftDeg + righDeg; + + PolynomialR032 retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + + float left = this.floatValue(l); + + for (int r = 0; r <= righDeg; r++) { + + float right = multiplicand.floatValue(r); + + coefficients.add(l + r, left * right); + } + } + + return retVal; + } + + @Override + public PolynomialR032 negate() { + + int size = 1 + this.degree(); + + PolynomialR032 retVal = this.newInstance(size); + + for (int p = 0; p < size; p++) { + retVal.set(p, -this.floatValue(p)); } + + return retVal; } @Override @@ -77,8 +142,13 @@ protected Double getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialR032(Array1D.R032.make(size)); + protected PolynomialR032 newInstance(final int size) { + return new PolynomialR032(ArrayR032.make(size)); + } + + @Override + PolynomialR032 one() { + return ONE; } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialR064.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialR064.java index c8aa87426..ef504721f 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialR064.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialR064.java @@ -21,26 +21,57 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayR064; +import org.ojalgo.array.BasicArray; import org.ojalgo.function.constant.PrimitiveMath; import org.ojalgo.matrix.decomposition.QR; import org.ojalgo.matrix.store.Primitive64Store; import org.ojalgo.structure.Access1D; -public class PolynomialR064 extends AbstractPolynomial { +public class PolynomialR064 extends AbstractPolynomial { + + public static final PolynomialR064 ONE = PolynomialR064.wrap(PrimitiveMath.ONE); + + public static PolynomialR064 wrap(final double... coefficients) { + return new PolynomialR064(ArrayR064.wrap(coefficients)); + } public PolynomialR064(final int degree) { - super(Array1D.R064.make(degree + 1)); + super(ArrayR064.make(degree + 1)); } - PolynomialR064(final Array1D coefficients) { + PolynomialR064(final BasicArray coefficients) { super(coefficients); } + @Override + public PolynomialR064 add(final PolynomialFunction addend) { + + int leftDeg = this.degree(); + int righDeg = addend.degree(); + + int retSize = 1 + Math.max(leftDeg, righDeg); + + PolynomialR064 retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + coefficients.add(l, this.doubleValue(l)); + } + + for (int r = 0; r <= righDeg; r++) { + coefficients.add(r, addend.doubleValue(r)); + } + + return retVal; + } + + @Override public void estimate(final Access1D x, final Access1D y) { this.estimate(x, y, Primitive64Store.FACTORY, QR.R064); } + @Override public Double integrate(final Double fromPoint, final Double toPoint) { PolynomialFunction primitive = this.buildPrimitive(); @@ -51,15 +82,49 @@ public Double integrate(final Double fromPoint, final Double toPoint) { return Double.valueOf(toVal - fromVal); } + @Override public Double invoke(final Double arg) { return Double.valueOf(this.invoke(arg.doubleValue())); } - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, coefficients.doubleValue(p)); + @Override + public PolynomialR064 multiply(final PolynomialFunction multiplicand) { + + int leftDeg = this.degree(); + int righDeg = multiplicand.degree(); + + int retSize = 1 + leftDeg + righDeg; + + PolynomialR064 retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + + double left = this.doubleValue(l); + + for (int r = 0; r <= righDeg; r++) { + + double right = multiplicand.doubleValue(r); + + coefficients.add(l + r, left * right); + } + } + + return retVal; + } + + @Override + public PolynomialR064 negate() { + + int size = 1 + this.degree(); + + PolynomialR064 retVal = this.newInstance(size); + + for (int p = 0; p < size; p++) { + retVal.set(p, -this.doubleValue(p)); } + + return retVal; } @Override @@ -77,8 +142,12 @@ protected Double getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialR064(Array1D.R064.make(size)); + protected PolynomialR064 newInstance(final int size) { + return new PolynomialR064(ArrayR064.make(size)); } + @Override + PolynomialR064 one() { + return ONE; + } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialR128.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialR128.java index cebb423f7..c92ad7df5 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialR128.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialR128.java @@ -21,56 +21,34 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayR128; +import org.ojalgo.array.BasicArray; import org.ojalgo.matrix.decomposition.QR; import org.ojalgo.matrix.store.GenericStore; import org.ojalgo.scalar.Quadruple; import org.ojalgo.structure.Access1D; -public final class PolynomialR128 extends AbstractPolynomial { +public final class PolynomialR128 extends ScalarPolynomial { + + public static final PolynomialR128 ONE = PolynomialR128.wrap(Quadruple.ONE); + + public static PolynomialR128 wrap(final Quadruple... coefficients) { + return new PolynomialR128(ArrayR128.wrap(coefficients)); + } public PolynomialR128(final int degree) { - super(Array1D.R128.make(degree + 1)); + super(ArrayR128.make(degree + 1)); } - PolynomialR128(final Array1D coefficients) { + PolynomialR128(final BasicArray coefficients) { super(coefficients); } + @Override public void estimate(final Access1D x, final Access1D y) { this.estimate(x, y, GenericStore.R128, QR.R128); } - public Quadruple integrate(final Quadruple fromPoint, final Quadruple toPoint) { - - PolynomialFunction primitive = this.buildPrimitive(); - - Quadruple fromVal = primitive.invoke(fromPoint); - Quadruple toVal = primitive.invoke(toPoint); - - return toVal.subtract(fromVal); - } - - public Quadruple invoke(final Quadruple arg) { - - int power = this.degree(); - - Quadruple retVal = this.get(power); - - while (--power >= 0) { - retVal = this.get(power).add(arg.multiply(retVal)); - } - - return retVal; - } - - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, Quadruple.valueOf(coefficients.get(p))); - } - } - @Override protected Quadruple getDerivativeFactor(final int power) { int nextIndex = power + 1; @@ -86,8 +64,13 @@ protected Quadruple getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialR128(Array1D.R128.make(size)); + protected PolynomialR128 newInstance(final int size) { + return new PolynomialR128(ArrayR128.make(size)); + } + + @Override + PolynomialR128 one() { + return ONE; } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PolynomialR256.java b/src/main/java/org/ojalgo/function/polynomial/PolynomialR256.java index 5c7d7ad14..c7bf691d0 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PolynomialR256.java +++ b/src/main/java/org/ojalgo/function/polynomial/PolynomialR256.java @@ -23,26 +23,33 @@ import java.math.BigDecimal; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.ArrayR256; +import org.ojalgo.array.BasicArray; import org.ojalgo.function.constant.BigMath; import org.ojalgo.structure.Access1D; -import org.ojalgo.type.TypeUtils; /** * BigPolynomial * * @author apete */ -public class PolynomialR256 extends AbstractPolynomial { +public class PolynomialR256 extends AbstractPolynomial { + + public static final PolynomialR256 ONE = PolynomialR256.wrap(BigDecimal.ONE); + + public static PolynomialR256 wrap(final BigDecimal... coefficients) { + return new PolynomialR256(ArrayR256.wrap(coefficients)); + } public PolynomialR256(final int degree) { - super(Array1D.R256.make(degree + 1)); + super(ArrayR256.make(degree + 1)); } - PolynomialR256(final Array1D coefficients) { + PolynomialR256(final BasicArray coefficients) { super(coefficients); } + @Override public void estimate(final Access1D x, final Access1D y) { PolynomialQ128 delegate = new PolynomialQ128(this.degree()); @@ -52,6 +59,7 @@ public void estimate(final Access1D x, final Access1D y) { this.set(delegate); } + @Override public BigDecimal integrate(final BigDecimal fromPoint, final BigDecimal toPoint) { PolynomialFunction primitive = this.buildPrimitive(); @@ -62,6 +70,7 @@ public BigDecimal integrate(final BigDecimal fromPoint, final BigDecimal toPoint return toVal.subtract(fromVal); } + @Override public BigDecimal invoke(final BigDecimal arg) { int power = this.degree(); @@ -75,11 +84,44 @@ public BigDecimal invoke(final BigDecimal arg) { return retVal; } - public void set(final Access1D coefficients) { - int limit = Math.min(this.size(), coefficients.size()); - for (int p = 0; p < limit; p++) { - this.set(p, TypeUtils.toBigDecimal(coefficients.get(p))); + @Override + public PolynomialR256 multiply(final PolynomialFunction multiplicand) { + + int leftDeg = this.degree(); + int righDeg = multiplicand.degree(); + + int retSize = 1 + leftDeg + righDeg; + + PolynomialR256 retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + + BigDecimal left = this.get(l); + + for (int r = 0; r <= righDeg; r++) { + + BigDecimal right = multiplicand.get(r); + + coefficients.add(l + r, left.multiply(right)); + } } + + return retVal; + } + + @Override + public PolynomialR256 negate() { + + int size = 1 + this.degree(); + + PolynomialR256 retVal = this.newInstance(size); + + for (int p = 0; p < size; p++) { + retVal.set(p, this.get(p).negate()); + } + + return retVal; } @Override @@ -97,8 +139,13 @@ protected BigDecimal getPrimitiveFactor(final int power) { } @Override - protected AbstractPolynomial makeInstance(final int size) { - return new PolynomialR256(Array1D.R256.make(size)); + protected PolynomialR256 newInstance(final int size) { + return new PolynomialR256(ArrayR256.make(size)); + } + + @Override + PolynomialR256 one() { + return ONE; } } diff --git a/src/main/java/org/ojalgo/function/polynomial/PrimitivePolynomial.java b/src/main/java/org/ojalgo/function/polynomial/PrimitivePolynomial.java index 0602c9402..7b2346fa1 100644 --- a/src/main/java/org/ojalgo/function/polynomial/PrimitivePolynomial.java +++ b/src/main/java/org/ojalgo/function/polynomial/PrimitivePolynomial.java @@ -21,7 +21,7 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.BasicArray; /** * @deprecated v53 use {@link PolynomialR064} instead @@ -33,7 +33,7 @@ public PrimitivePolynomial(final int degree) { super(degree); } - public PrimitivePolynomial(final Array1D coefficients) { + public PrimitivePolynomial(final BasicArray coefficients) { super(coefficients); } diff --git a/src/main/java/org/ojalgo/function/polynomial/RationalPolynomial.java b/src/main/java/org/ojalgo/function/polynomial/RationalPolynomial.java index b7951deb4..5f7a466c5 100644 --- a/src/main/java/org/ojalgo/function/polynomial/RationalPolynomial.java +++ b/src/main/java/org/ojalgo/function/polynomial/RationalPolynomial.java @@ -21,7 +21,7 @@ */ package org.ojalgo.function.polynomial; -import org.ojalgo.array.Array1D; +import org.ojalgo.array.BasicArray; import org.ojalgo.scalar.RationalNumber; /** @@ -30,7 +30,7 @@ @Deprecated public final class RationalPolynomial extends PolynomialQ128 { - public RationalPolynomial(final Array1D coefficients) { + public RationalPolynomial(final BasicArray coefficients) { super(coefficients); } diff --git a/src/main/java/org/ojalgo/function/polynomial/ScalarPolynomial.java b/src/main/java/org/ojalgo/function/polynomial/ScalarPolynomial.java new file mode 100644 index 000000000..d94efe544 --- /dev/null +++ b/src/main/java/org/ojalgo/function/polynomial/ScalarPolynomial.java @@ -0,0 +1,103 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.polynomial; + +import org.ojalgo.array.BasicArray; +import org.ojalgo.scalar.Scalar; + +abstract class ScalarPolynomial, P extends ScalarPolynomial> extends AbstractPolynomial { + + ScalarPolynomial(final BasicArray coefficients) { + super(coefficients); + } + + @Override + public final N integrate(final N fromPoint, final N toPoint) { + + PolynomialFunction primitive = this.buildPrimitive(); + + N fromVal = primitive.invoke(fromPoint); + N toVal = primitive.invoke(toPoint); + + return toVal.subtract(fromVal).get(); + } + + @Override + public final N invoke(final N arg) { + + int power = this.degree(); + + Scalar retVal = this.get(power); + + while (--power >= 0) { + retVal = this.get(power).add(arg.multiply(retVal)); + } + + return retVal.get(); + } + + @Override + public P multiply(final PolynomialFunction multiplicand) { + + int leftDeg = this.degree(); + int righDeg = multiplicand.degree(); + + int retSize = 1 + leftDeg + righDeg; + + P retVal = this.newInstance(retSize); + BasicArray coefficients = retVal.coefficients(); + + for (int l = 0; l <= leftDeg; l++) { + + N left = this.get(l); + + for (int r = 0; r <= righDeg; r++) { + + N right = multiplicand.get(r); + + coefficients.add(l + r, left.multiply(right)); + } + } + + return retVal; + } + + @Override + public PolynomialFunction negate() { + + int size = 1 + this.degree(); + + P retVal = this.newInstance(size); + + for (int p = 0; p < size; p++) { + retVal.set(p, this.get(p).negate()); + } + + return retVal; + } + + @Override + double norm(final int power) { + return this.get(power).norm(); + } + +} diff --git a/src/main/java/org/ojalgo/function/series/FourierSeries.java b/src/main/java/org/ojalgo/function/series/FourierSeries.java new file mode 100644 index 000000000..e64865c98 --- /dev/null +++ b/src/main/java/org/ojalgo/function/series/FourierSeries.java @@ -0,0 +1,162 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.series; + +import java.util.function.DoubleUnaryOperator; + +import org.ojalgo.data.transform.DiscreteFourierTransform; +import org.ojalgo.function.PrimitiveFunction; +import org.ojalgo.function.PrimitiveFunction.SampleDomain; +import org.ojalgo.function.constant.PrimitiveMath; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.scalar.ComplexNumber; +import org.ojalgo.type.context.NumberContext; + +/** + * This is the real coefficient trigonometric form of the Fourier series. + */ +public class FourierSeries implements PrimitiveFunction.Unary { + + private static final NumberContext ACCURACY = NumberContext.of(8); + + public static FourierSeries estimate(final DoubleUnaryOperator function, final PrimitiveFunction.SampleDomain sampleDomain) { + + SampleDomain adjustedDomain = sampleDomain.adjustToPowerOf2(); + + MatrixStore transform = DiscreteFourierTransform.sample(function, adjustedDomain); + + int size = transform.size(); + int length = (size - 1) / 2; + double half = size / 2D; + + double[] cosCoefficients = new double[length]; + double[] sinCoefficients = new double[length]; + + double magnitude = PrimitiveMath.ONE; + + ComplexNumber base = transform.get(0).divide(size); + double baseNorm = base.norm(); + + if (!ACCURACY.isSmall(magnitude, baseNorm)) { + cosCoefficients[0] = base.getReal(); + sinCoefficients[0] = base.getImaginary(); + magnitude = baseNorm; + } + + for (int i = 1; i < length; i++) { + + ComplexNumber coefficient = transform.get(size - i).divide(half); + double real = coefficient.getReal(); + double imaginary = coefficient.getImaginary(); + + if (!ACCURACY.isSmall(magnitude, real)) { + cosCoefficients[i] = real; + } + if (!ACCURACY.isSmall(magnitude, imaginary)) { + sinCoefficients[i] = imaginary; + } + } + + return new FourierSeries(adjustedDomain.period(), cosCoefficients, sinCoefficients); + } + + public static FourierSeries estimate(final PeriodicFunction function, final int nbSamples) { + return FourierSeries.estimate(function, function.getSampleDomain(nbSamples)); + } + + private final double myAngularFrequency; + private final double[] myCosCoefficients; + private final double[] mySinCoefficients; + + /** + * @param period The period of the function + * @param coefficients The Fourier coefficients. The first coefficient is the constant term (the real part + * of that complex number), then the following coefficients are the coefficients of the cos and sin + * terms with increasing frequency. + */ + public FourierSeries(final double period, final ComplexNumber... coefficients) { + + super(); + + if (coefficients.length == 0) { + throw new IllegalArgumentException(); + } + + myAngularFrequency = PrimitiveMath.TWO_PI / period; + + myCosCoefficients = new double[coefficients.length]; + mySinCoefficients = new double[coefficients.length]; + + for (int i = 0; i < coefficients.length; i++) { + + ComplexNumber coefficient = coefficients[i]; + + if (coefficient != null) { + myCosCoefficients[i] = coefficient.getReal(); + mySinCoefficients[i] = coefficient.getImaginary(); + } + } + } + + FourierSeries(final double period, final double[] cosCoefficients, final double[] sinCoefficients) { + + super(); + + myAngularFrequency = PrimitiveMath.TWO_PI / period; + + myCosCoefficients = cosCoefficients; + mySinCoefficients = sinCoefficients; + + if (myCosCoefficients.length == 0) { + throw new IllegalArgumentException(); + } + } + + @Override + public double invoke(final double arg) { + + double retVal = PrimitiveMath.ZERO; + + for (int i = 0; i < myCosCoefficients.length; i++) { + + double coefficient = myCosCoefficients[i]; + + if (Double.isFinite(retVal) && coefficient != PrimitiveMath.ZERO) { + double a = arg * i * myAngularFrequency; + retVal += coefficient * Math.cos(a); + } + } + + for (int i = 1; i < mySinCoefficients.length; i++) { + + double coefficient = mySinCoefficients[i]; + + if (Double.isFinite(retVal) && coefficient != PrimitiveMath.ZERO) { + double a = arg * i * myAngularFrequency; + retVal += coefficient * Math.sin(a); + } + } + + return retVal; + } + +} diff --git a/src/main/java/org/ojalgo/function/series/PeriodicFunction.java b/src/main/java/org/ojalgo/function/series/PeriodicFunction.java new file mode 100644 index 000000000..ea736469d --- /dev/null +++ b/src/main/java/org/ojalgo/function/series/PeriodicFunction.java @@ -0,0 +1,149 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.series; + +import static org.ojalgo.function.constant.PrimitiveMath.*; + +import java.util.function.DoubleUnaryOperator; + +import org.ojalgo.function.PrimitiveFunction; +import org.ojalgo.function.PrimitiveFunction.SampleDomain; + + +/** + * A periodic function is a function that repeats its values in regular intervals or periods. The most + * important examples are the trigonometric functions, which repeat over intervals of length 2π radians. + *

+ * This class allows you to create a periodic function from any other function. The base function only needs + * to be defined over the interval specified in this class, The resulting function will repeat the base + * function in regular periods. The default period is 2π, but it's possible to specify any other period. The + * default origin is zero, but it's possible to specify any other origin. + */ +public final class PeriodicFunction implements PrimitiveFunction.Unary { + + /** + * https://en.wikipedia.org/wiki/Sawtooth_wave + *

+ * Domain: (-∞, ∞) by repeating/shifting the definition interval [0, 2π).
+ * Range: [-1, 1]. + */ + public static final PeriodicFunction SAWTOOTH = PeriodicFunction.of(arg -> Math.atan(Math.tan(arg / TWO)) / HALF_PI); + /** + * https://en.wikipedia.org/wiki/Sine_wave + *

+ * Domain: (-∞, ∞) by repeating/shifting the definition interval [0, 2π).
+ * Range: [-1, 1]. + */ + public static final PeriodicFunction SINE = PeriodicFunction.of(SIN); + /** + * https://en.wikipedia.org/wiki/Square_wave + *

+ * Domain: (-∞, ∞) by repeating/shifting the definition interval [0, 2π).
+ * Range: [-1, 1]. + */ + public static final PeriodicFunction SQUARE = PeriodicFunction.of(arg -> Math.signum(Math.sin(arg))); + /** + * https://en.wikipedia.org/wiki/Triangle_wave + *

+ * Domain: (-∞, ∞) by repeating/shifting the definition interval [0, 2π).
+ * Range: [-1, 1]. + */ + public static final PeriodicFunction TRIANGLE = PeriodicFunction.of(arg -> Math.asin(Math.sin(arg)) / HALF_PI); + + public static PeriodicFunction of(final double origin, final DoubleUnaryOperator shape) { + return new PeriodicFunction(origin, shape, TWO_PI); + } + + /** + * Origin at zero, and period 2π [0, 2π). + */ + public static PeriodicFunction of(final DoubleUnaryOperator shape) { + return new PeriodicFunction(ZERO, shape, TWO_PI); + } + + public static DoubleUnaryOperator of(final DoubleUnaryOperator shape, final double period) { + return new PeriodicFunction(ZERO, shape, period); + } + + /** + * Origin at -π, and period 2π [-π, π). + */ + public static PeriodicFunction ofCentered(final DoubleUnaryOperator shape) { + return new PeriodicFunction(-PI, shape, TWO_PI); + } + + private final double myOrigin; + private final double myPeriod; + private final DoubleUnaryOperator myShape; + + PeriodicFunction(final double origin, final DoubleUnaryOperator shape, final double period) { + super(); + myOrigin = origin; + myPeriod = period; + myShape = shape; + } + + public SampleDomain getSampleDomain(final int nbSamples) { + return new SampleDomain(myPeriod, nbSamples); + } + + @Override + public double invoke(final double arg) { + double local = arg; + while (local < myOrigin) { + local += myPeriod; + } + while (local >= myOrigin + myPeriod) { + local -= myPeriod; + } + return myShape.applyAsDouble(local); + } + + /** + * interval = [origin, origin+period) + */ + public PeriodicFunction withInterval(final double origin, final double period) { + return new PeriodicFunction(origin, myShape, period); + } + + /** + * interval = [origin, origin+period) + */ + public PeriodicFunction withOrigin(final double origin) { + return new PeriodicFunction(origin, myShape, myPeriod); + } + + /** + * interval = [origin, origin+period) + */ + public PeriodicFunction withPeriod(final double period) { + return new PeriodicFunction(myOrigin, myShape, period); + } + + /** + * The shape is the function that is repeated in regular intervals. + */ + public PeriodicFunction withShape(final DoubleUnaryOperator shape) { + return new PeriodicFunction(myOrigin, shape, myPeriod); + } + +} diff --git a/src/main/java/org/ojalgo/matrix/BasicMatrix.java b/src/main/java/org/ojalgo/matrix/BasicMatrix.java index caaa518bc..9e7f4e67f 100644 --- a/src/main/java/org/ojalgo/matrix/BasicMatrix.java +++ b/src/main/java/org/ojalgo/matrix/BasicMatrix.java @@ -231,7 +231,7 @@ public M conjugate() { } /** - * The returned instance can be have its elements mutated in various ways, while the size/shape is fixed. + * The returned instance can have its elements mutated in various ways, while the size/shape is fixed. * * @return A fully mutable matrix builder with the elements initially set to a copy of this matrix – * always creates a full dense copy. @@ -431,21 +431,20 @@ public long indexOfLargest() { *

*
    *
  • "right inverse": [this][right inverse]=[I]. You may calculate it using - * {@linkplain #solve(Access2D)}.
    • + * {@linkplain #solve(Access2D)}. *
  • "left inverse": [left inverse][this]=[I]. You may calculate it using {@linkplain #solve(Access2D)} - * and transposing.
    • + * and transposing. *
  • "generalised inverse": [this][generalised inverse][this]=[this]. Note that if [this] is singular or - * non-square, then [generalised inverse] is not unique.
    • + * non-square, then [generalised inverse] is not unique. *
  • "pseudoinverse": The generalised inverse (there are typically/possibly many) with the smallest - * frobenius norm is called the pseudoinverse. You may calculate it using the {@linkplain QR} or - * {@linkplain SingularValue} decompositions.
    • + * frobenius norm is called the (Moore-Penrose) pseudoinverse. You may calculate it using the + * {@linkplain QR} or {@linkplain SingularValue} decompositions. *
  • "inverse": *
      - *
    • If [left inverse]=[right inverse] then it is also [inverse].
      • + *
    • If [left inverse]=[right inverse] then it is also [inverse]. *
    • If [this] is square and has full rank then the [generalised inverse] is unique, with the - * [pseudoinverse] given, and equal to [inverse].
      • + * [pseudoinverse] given, and equal to [inverse]. *
    - *
    • *
* * @return The "best possible" inverse.... diff --git a/src/main/java/org/ojalgo/matrix/store/PhysicalStore.java b/src/main/java/org/ojalgo/matrix/store/PhysicalStore.java index 6b2608929..295240dd1 100644 --- a/src/main/java/org/ojalgo/matrix/store/PhysicalStore.java +++ b/src/main/java/org/ojalgo/matrix/store/PhysicalStore.java @@ -127,14 +127,14 @@ default RowsSupplier makeRowsSupplier(final int numberOfColumns) { return new RowsSupplier<>(this, numberOfColumns); } - default MatrixStore makeSingle(final N element) { - return new SingleStore<>(this, element); - } - default MatrixStore makeSingle(final double element) { return this.makeSingle(this.scalar().cast(element)); } + default MatrixStore makeSingle(final N element) { + return new SingleStore<>(this, element); + } + @Override default SparseStore makeSparse(final long rowsCount, final long columnsCount) { return SparseStore.makeSparse(this, rowsCount, columnsCount); @@ -164,6 +164,10 @@ default MatrixStore makeWrapper(final Access2D access) { return new WrapperStore<>(this, access); } + default MatrixStore makeWrapperColumn(final Access1D access) { + return new WrapperStore<>(access, this); + } + default MatrixStore makeZero(final long rowsCount, final long columnsCount) { return new ZeroStore<>(this, rowsCount, columnsCount); } @@ -193,6 +197,33 @@ default TensorFactory2D tensor2D() { */ List asList(); + default int indexOfLargestInColumn(final int row, final int col) { + long structure = this.countRows(); + long first = Structure2D.index(structure, row, col); + long limit = Structure2D.index(structure, 0L, col + 1L); + long step = 1L; + long largest = AMAX.invoke(this, first, limit, step); + return Math.toIntExact(largest % structure); + } + + default int indexOfLargestInRow(final int row, final int col) { + long structure = this.countRows(); + long first = Structure2D.index(structure, row, col); + long limit = Structure2D.index(structure, 0L, this.countColumns()); + long step = structure; + long largest = AMAX.invoke(this, first, limit, step); + return Math.toIntExact(largest / structure); + } + + default int indexOfLargestOnDiagonal(final int row, final int col) { + long structure = this.countRows(); + long first = Structure2D.index(structure, row, col); + long limit = Structure2D.index(structure, 0L, this.countColumns()); + long step = structure + 1L; + long largest = AMAX.invoke(this, first, limit, step); + return Math.toIntExact(largest / structure); + } + @Override default void modifyAny(final Transformation2D modifier) { modifier.transform(this); @@ -267,31 +298,4 @@ default void supplyTo(final TransformableRegion receiver) { */ void transformRight(Rotation transformation); - default int indexOfLargestInColumn(final int row, final int col) { - long structure = this.countRows(); - long first = Structure2D.index(structure, row, col); - long limit = Structure2D.index(structure, 0L, col + 1L); - long step = 1L; - long largest = AMAX.invoke(this, first, limit, step); - return Math.toIntExact(largest % structure); - } - - default int indexOfLargestInRow(final int row, final int col) { - long structure = this.countRows(); - long first = Structure2D.index(structure, row, col); - long limit = Structure2D.index(structure, 0L, this.countColumns()); - long step = structure; - long largest = AMAX.invoke(this, first, limit, step); - return Math.toIntExact(largest / structure); - } - - default int indexOfLargestOnDiagonal(final int row, final int col) { - long structure = this.countRows(); - long first = Structure2D.index(structure, row, col); - long limit = Structure2D.index(structure, 0L, this.countColumns()); - long step = structure + 1L; - long largest = AMAX.invoke(this, first, limit, step); - return Math.toIntExact(largest / structure); - } - } diff --git a/src/main/java/org/ojalgo/matrix/store/SparseStore.java b/src/main/java/org/ojalgo/matrix/store/SparseStore.java index 21d20380d..69356abfa 100644 --- a/src/main/java/org/ojalgo/matrix/store/SparseStore.java +++ b/src/main/java/org/ojalgo/matrix/store/SparseStore.java @@ -21,7 +21,8 @@ */ package org.ojalgo.matrix.store; -import static org.ojalgo.function.constant.PrimitiveMath.*; +import static org.ojalgo.function.constant.PrimitiveMath.E; +import static org.ojalgo.function.constant.PrimitiveMath.ZERO; import java.util.Arrays; @@ -34,7 +35,6 @@ import org.ojalgo.function.UnaryFunction; import org.ojalgo.function.VoidFunction; import org.ojalgo.function.aggregator.Aggregator; -import org.ojalgo.matrix.operation.MultiplyBoth; import org.ojalgo.scalar.ComplexNumber; import org.ojalgo.scalar.Quadruple; import org.ojalgo.scalar.Quaternion; @@ -244,7 +244,7 @@ static > void multiply(final SparseStore left, final private final SparseArray myElements; private final int[] myFirsts; private final int[] myLimits; - private TransformableRegion.FillByMultiplying myMultiplyer; + private final TransformableRegion.FillByMultiplying myMultiplier; SparseStore(final PhysicalStore.Factory factory, final int rowsCount, final int columnsCount) { @@ -257,11 +257,7 @@ static > void multiply(final SparseStore left, final // Arrays.fill(myLimits, 0); // Behövs inte, redan 0 Class tmpType = factory.scalar().zero().get().getClass(); - if (tmpType.equals(Double.class)) { - myMultiplyer = (TransformableRegion.FillByMultiplying) MultiplyBoth.newPrimitive64(rowsCount, columnsCount); - } else { - myMultiplyer = (TransformableRegion.FillByMultiplying) MultiplyBoth.newGeneric(rowsCount, columnsCount); - } + myMultiplier = Subregion2D.findMultiplier(tmpType, rowsCount, columnsCount); } @Override @@ -315,7 +311,7 @@ public void fillByMultiplying(final Access1D left, final Access1D right) { ProgrammingError.throwForMultiplicationNotPossible(); } - myMultiplyer.invoke(this, left, complexity, right); + myMultiplier.invoke(this, left, complexity, right); } @Override @@ -667,27 +663,27 @@ public void reduceRows(final Aggregator aggregator, final Mutate1D receiver) { @Override public TransformableRegion regionByColumns(final int... columns) { - return new Subregion2D.ColumnsRegion<>(this, myMultiplyer, columns); + return new Subregion2D.ColumnsRegion<>(this, myMultiplier, columns); } @Override public TransformableRegion regionByLimits(final int rowLimit, final int columnLimit) { - return new Subregion2D.LimitRegion<>(this, myMultiplyer, rowLimit, columnLimit); + return new Subregion2D.LimitRegion<>(this, myMultiplier, rowLimit, columnLimit); } @Override public TransformableRegion regionByOffsets(final int rowOffset, final int columnOffset) { - return new Subregion2D.OffsetRegion<>(this, myMultiplyer, rowOffset, columnOffset); + return new Subregion2D.OffsetRegion<>(this, myMultiplier, rowOffset, columnOffset); } @Override public TransformableRegion regionByRows(final int... rows) { - return new Subregion2D.RowsRegion<>(this, myMultiplyer, rows); + return new Subregion2D.RowsRegion<>(this, myMultiplier, rows); } @Override public TransformableRegion regionByTransposing() { - return new Subregion2D.TransposedRegion<>(this, myMultiplyer); + return new Subregion2D.TransposedRegion<>(this, myMultiplier); } @Override diff --git a/src/main/java/org/ojalgo/matrix/store/Subregion2D.java b/src/main/java/org/ojalgo/matrix/store/Subregion2D.java index eeb224127..d43bbb125 100644 --- a/src/main/java/org/ojalgo/matrix/store/Subregion2D.java +++ b/src/main/java/org/ojalgo/matrix/store/Subregion2D.java @@ -28,6 +28,7 @@ import org.ojalgo.function.UnaryFunction; import org.ojalgo.matrix.operation.MultiplyBoth; import org.ojalgo.structure.Access1D; +import org.ojalgo.structure.Mutate2D; abstract class Subregion2D> implements TransformableRegion { @@ -713,6 +714,72 @@ public void set(final long row, final long col, final Comparable value) { } + static final class WrapperRegion> extends Subregion2D { + + private final Mutate2D.ModifiableReceiver myBase; + + WrapperRegion(final Mutate2D.ModifiableReceiver base) { + super(Subregion2D.findMultiplier(base.get(0, 0).getClass(), base.getRowDim(), base.getColDim()), base.getRowDim(), + base.getColDim()); + myBase = base; + } + + @Override + public void add(final long row, final long col, final Comparable addend) { + myBase.add(row, col, addend); + } + + @Override + public void add(final long row, final long col, final double addend) { + myBase.add(row, col, addend); + } + + @Override + public long countColumns() { + return myBase.countColumns(); + } + + @Override + public long countRows() { + return myBase.countRows(); + } + + @Override + public double doubleValue(final int row, final int col) { + return myBase.doubleValue(row, col); + } + + @Override + public N get(final long row, final long col) { + return myBase.get(row, col); + } + + @Override + public void modifyOne(final long row, final long col, final UnaryFunction modifier) { + myBase.modifyOne(row, col, modifier); + } + + @Override + public void set(final int row, final int col, final double value) { + myBase.set(row, col, value); + } + + @Override + public void set(final long row, final long col, final Comparable value) { + myBase.set(row, col, value); + } + + } + + static > TransformableRegion.FillByMultiplying findMultiplier(final Class tmpType, final int rowsCount, + final int columnsCount) { + if (tmpType.equals(Double.class)) { + return (TransformableRegion.FillByMultiplying) MultiplyBoth.newPrimitive64(rowsCount, columnsCount); + } else { + return (TransformableRegion.FillByMultiplying) MultiplyBoth.newGeneric(rowsCount, columnsCount); + } + } + private final TransformableRegion.FillByMultiplying myMultiplier; @SuppressWarnings("unused") @@ -737,7 +804,7 @@ private Subregion2D() { @Override public final void fillByMultiplying(final Access1D left, final Access1D right) { - final int complexity = Math.toIntExact(left.count() / this.countRows()); + int complexity = Math.toIntExact(left.count() / this.countRows()); if (complexity != Math.toIntExact(right.count() / this.countColumns())) { ProgrammingError.throwForMultiplicationNotPossible(); } diff --git a/src/main/java/org/ojalgo/matrix/store/TransformableRegion.java b/src/main/java/org/ojalgo/matrix/store/TransformableRegion.java index f5ef71670..a08c38032 100644 --- a/src/main/java/org/ojalgo/matrix/store/TransformableRegion.java +++ b/src/main/java/org/ojalgo/matrix/store/TransformableRegion.java @@ -22,6 +22,7 @@ package org.ojalgo.matrix.store; import org.ojalgo.structure.Access1D; +import org.ojalgo.structure.Mutate2D; import org.ojalgo.structure.Mutate2D.ModifiableReceiver; import org.ojalgo.structure.Transformation2D; @@ -43,6 +44,15 @@ default void invoke(final TransformableRegion product, final Access1D left } + static > TransformableRegion cast(final Mutate2D.ModifiableReceiver target) { + if (target instanceof TransformableRegion) { + return (TransformableRegion) target; + } else { + return new Subregion2D.WrapperRegion<>(target); + } + } + + @Override default void exchangeColumns(final long colA, final long colB) { N valA, valB; for (long i = 0L, limit = this.countRows(); i < limit; i++) { @@ -53,6 +63,7 @@ default void exchangeColumns(final long colA, final long colB) { } } + @Override default void exchangeRows(final long rowA, final long rowB) { N valA, valB; for (long j = 0L, limit = this.countColumns(); j < limit; j++) { @@ -65,6 +76,7 @@ default void exchangeRows(final long rowA, final long rowB) { void fillByMultiplying(final Access1D left, final Access1D right); + @Override default void modifyAny(final Transformation2D modifier) { modifier.transform(this); } diff --git a/src/main/java/org/ojalgo/matrix/store/WrapperStore.java b/src/main/java/org/ojalgo/matrix/store/WrapperStore.java index b8aa6e173..8dabc0557 100644 --- a/src/main/java/org/ojalgo/matrix/store/WrapperStore.java +++ b/src/main/java/org/ojalgo/matrix/store/WrapperStore.java @@ -21,74 +21,56 @@ */ package org.ojalgo.matrix.store; -import org.ojalgo.ProgrammingError; +import org.ojalgo.scalar.Scalar.Factory; import org.ojalgo.structure.Access1D; import org.ojalgo.structure.Access2D; +import org.ojalgo.structure.Structure2D; /** * @author apete */ final class WrapperStore> extends FactoryStore { - private final Access2D myAccess; + private final Access1D myAccess; + private final Factory myScalarFactory; + private final int myStructure; - private WrapperStore(final PhysicalStore.Factory factory, final int rowsCount, final int columnsCount) { - super(factory, rowsCount, columnsCount); - myAccess = null; - ProgrammingError.throwForIllegalInvocation(); - } - - WrapperStore(final PhysicalStore.Factory factory, final Access2D access) { + WrapperStore(final Access1D access1D, final PhysicalStore.Factory factory) { - super(factory, (int) access.countRows(), (int) access.countColumns()); + super(factory, access1D.size(), 1); - myAccess = access; - } - - @Override - public double doubleValue(final int aRow, final int aCol) { - return myAccess.doubleValue(aRow, aCol); + myAccess = access1D; + myScalarFactory = factory.scalar(); + myStructure = access1D.size(); } - @Override - public N get(final int aRow, final int aCol) { - return this.physical().scalar().cast(myAccess.get(aRow, aCol)); - } + WrapperStore(final PhysicalStore.Factory factory, final Access2D access2D) { - @Override - public void multiply(final Access1D right, final TransformableRegion target) { - // TODO Auto-generated method stub - super.multiply(right, target); - } + super(factory, access2D.getRowDim(), access2D.getColDim()); - @Override - public MatrixStore multiply(final double scalar) { - // TODO Auto-generated method stub - return super.multiply(scalar); + myAccess = access2D; + myScalarFactory = factory.scalar(); + myStructure = access2D.getRowDim(); } @Override - public MatrixStore multiply(final MatrixStore right) { - // TODO Auto-generated method stub - return super.multiply(right); + public double doubleValue(final int index) { + return myAccess.doubleValue(index); } @Override - public MatrixStore multiply(final N scalar) { - // TODO Auto-generated method stub - return super.multiply(scalar); + public double doubleValue(final int row, final int col) { + return myAccess.doubleValue(Structure2D.index(myStructure, row, col)); } @Override - public N multiplyBoth(final Access1D leftAndRight) { - // TODO Auto-generated method stub - return super.multiplyBoth(leftAndRight); + public N get(final int row, final int col) { + return myScalarFactory.cast(myAccess.get(Structure2D.index(myStructure, row, col))); } @Override - public ElementsSupplier premultiply(final Access1D left) { - // TODO Auto-generated method stub - return super.premultiply(left); + public N get(final long index) { + return myScalarFactory.cast(myAccess.get(index)); } } diff --git a/src/main/java/org/ojalgo/scalar/ComplexNumber.java b/src/main/java/org/ojalgo/scalar/ComplexNumber.java index 1bbe71ea0..10e2d6bf7 100644 --- a/src/main/java/org/ojalgo/scalar/ComplexNumber.java +++ b/src/main/java/org/ojalgo/scalar/ComplexNumber.java @@ -80,13 +80,18 @@ public ComplexNumber zero() { }; /** - * Complex number {@code i}, satisfies i2 = -1; + * Complex number {@code i}, Z = (0.0 + 1.0i), satisfies i2 = -1; */ public static final ComplexNumber I = new ComplexNumber(PrimitiveMath.ZERO, PrimitiveMath.ONE); + /** * Complex number Z = (+∞ + 0.0i) */ public static final ComplexNumber INFINITY = ComplexNumber.makePolar(Double.POSITIVE_INFINITY, PrimitiveMath.ZERO); + /** + * Complex number {@code -i}, Z = (0.0 - 1.0i) + */ + public static final ComplexNumber N = new ComplexNumber(PrimitiveMath.ZERO, PrimitiveMath.NEG); /** * Complex number Z = (NaN + NaNi) */ @@ -163,44 +168,61 @@ public static boolean isSmall(final double comparedTo, final ComplexNumber value */ public static ComplexNumber makePolar(final double norm, final double phase) { - double tmpStdPhase = phase % PrimitiveMath.TWO_PI; - if (tmpStdPhase < PrimitiveMath.ZERO) { - tmpStdPhase += PrimitiveMath.TWO_PI; + double stdPhase = phase % PrimitiveMath.TWO_PI; + if (stdPhase < PrimitiveMath.ZERO) { + stdPhase += PrimitiveMath.TWO_PI; } - if (tmpStdPhase <= ARGUMENT_TOLERANCE) { + if (stdPhase <= ARGUMENT_TOLERANCE) { return new ComplexNumber(norm); - } - if (PrimitiveMath.ABS.invoke(tmpStdPhase - PrimitiveMath.PI) <= ARGUMENT_TOLERANCE) { + } else if (PrimitiveMath.ABS.invoke(stdPhase - PrimitiveMath.PI) <= ARGUMENT_TOLERANCE) { return new ComplexNumber(-norm); - } - double tmpRe = PrimitiveMath.ZERO; - if (norm != PrimitiveMath.ZERO) { - final double tmpCos = PrimitiveMath.COS.invoke(tmpStdPhase); - if (tmpCos != PrimitiveMath.ZERO) { - tmpRe = norm * tmpCos; + } else { + + double tmpRe = PrimitiveMath.ZERO; + if (norm != PrimitiveMath.ZERO) { + double tmpCos = PrimitiveMath.COS.invoke(stdPhase); + if (tmpCos != PrimitiveMath.ZERO) { + tmpRe = norm * tmpCos; + } } - } - double tmpIm = PrimitiveMath.ZERO; - if (norm != PrimitiveMath.ZERO) { - final double tmpSin = PrimitiveMath.SIN.invoke(tmpStdPhase); - if (tmpSin != PrimitiveMath.ZERO) { - tmpIm = norm * tmpSin; + double tmpIm = PrimitiveMath.ZERO; + if (norm != PrimitiveMath.ZERO) { + double tmpSin = PrimitiveMath.SIN.invoke(stdPhase); + if (tmpSin != PrimitiveMath.ZERO) { + tmpIm = norm * tmpSin; + } } - } - return new ComplexNumber(tmpRe, tmpIm); + return new ComplexNumber(tmpRe, tmpIm); + } } public static ComplexNumber makeRotation(final double angle) { return new ComplexNumber(PrimitiveMath.COS.invoke(angle), PrimitiveMath.SIN.invoke(angle)); } + public static ComplexNumber newUnitRoot(final int nbRoots) { + return ComplexNumber.makePolar(PrimitiveMath.ONE, -PrimitiveMath.TWO_PI / nbRoots); + } + + public static ComplexNumber[] newUnitRoots(final int nbRoots) { + + ComplexNumber[] retVal = new ComplexNumber[nbRoots]; + + for (int i = 0; i < retVal.length; i++) { + retVal[i] = ComplexNumber.makePolar(PrimitiveMath.ONE, i * -PrimitiveMath.TWO_PI / nbRoots); + + } + + return retVal; + } + /** * Static factory method returning a complex number from cartesian coordinates. * @@ -361,11 +383,13 @@ public ComplexNumber divide(final ComplexNumber arg) { return new ComplexNumber((myRealValue + r * i) / d, (i - r * myRealValue) / d); - } - final double r = tmpRe / tmpIm; - final double d = tmpIm + r * tmpRe; + } else { - return new ComplexNumber((r * myRealValue + i) / d, (r * i - myRealValue) / d); + final double r = tmpRe / tmpIm; + final double d = tmpIm + r * tmpRe; + + return new ComplexNumber((r * myRealValue + i) / d, (r * i - myRealValue) / d); + } } /** @@ -534,10 +558,10 @@ public long longValue() { @Override public ComplexNumber multiply(final ComplexNumber arg) { - final double tmpRe = arg.doubleValue(); - final double tmpIm = arg.i; + double argRe = arg.doubleValue(); + double argIm = arg.i; - return new ComplexNumber(myRealValue * tmpRe - i * tmpIm, myRealValue * tmpIm + i * tmpRe); + return new ComplexNumber(myRealValue * argRe - i * argIm, myRealValue * argIm + i * argRe); } /** @@ -551,6 +575,22 @@ public ComplexNumber multiply(final double arg) { return new ComplexNumber(myRealValue * arg, i * arg); } + /** + * The imaginary part of the complex number resulting when multiplying this with a complex number whose + * real and imaginary parts are argRe and argIm. + */ + public double multiplyIm(final double argRe, final double argIm) { + return myRealValue * argIm + i * argRe; + } + + /** + * The real part of the complex number resulting when multiplying this with a complex number whose real + * and imaginary parts are argRe and argIm. + */ + public double multiplyRe(final double argRe, final double argIm) { + return myRealValue * argRe - i * argIm; + } + /** * Performs the unary operation '-'. * diff --git a/src/main/java/org/ojalgo/structure/Access1D.java b/src/main/java/org/ojalgo/structure/Access1D.java index 5f0e45cd7..f09add4cb 100644 --- a/src/main/java/org/ojalgo/structure/Access1D.java +++ b/src/main/java/org/ojalgo/structure/Access1D.java @@ -251,7 +251,7 @@ public String toString() { } /** - * Tests if the two data strauctures are numerically equal to the given accuracy. (Only works with real + * Tests if the two data structures are numerically equal to the given accuracy. (Only works with real * numbers, and can't handle more than "double precision".) You have to implement your own version to * handle other cases. */ diff --git a/src/main/java/org/ojalgo/type/TypeUtils.java b/src/main/java/org/ojalgo/type/TypeUtils.java index ce960629b..048eb9700 100644 --- a/src/main/java/org/ojalgo/type/TypeUtils.java +++ b/src/main/java/org/ojalgo/type/TypeUtils.java @@ -56,7 +56,7 @@ public static String format(final String messagePattern, final Object... args) { int first = 0; int limit = patternLength; - final StringBuilder retVal = new StringBuilder(patternLength + (argsCount * 20)); + final StringBuilder retVal = new StringBuilder(patternLength + argsCount * 20); for (int a = 0; a < argsCount; a++) { @@ -64,32 +64,30 @@ public static String format(final String messagePattern, final Object... args) { if (limit == -1) { retVal.append(ASCII.SP); - if ((args[a] instanceof Access1D)) { + if (args[a] instanceof Access1D) { retVal.append(Access1D.toString((Access1D) args[a])); - } else if ((args[a] instanceof double[])) { + } else if (args[a] instanceof double[]) { retVal.append(Arrays.toString((double[]) args[a])); } else { retVal.append(args[a]); } } else { retVal.append(messagePattern.substring(first, limit)); - if ((args[a] instanceof Access1D)) { - retVal.append(Access1D.toString((Access1D) args[a])); - } else if ((args[a] instanceof double[])) { + if (args[a] instanceof double[]) { retVal.append(Arrays.toString((double[]) args[a])); - } else if ((args[a] instanceof float[])) { + } else if (args[a] instanceof float[]) { retVal.append(Arrays.toString((float[]) args[a])); - } else if ((args[a] instanceof long[])) { + } else if (args[a] instanceof long[]) { retVal.append(Arrays.toString((long[]) args[a])); - } else if ((args[a] instanceof int[])) { + } else if (args[a] instanceof int[]) { retVal.append(Arrays.toString((int[]) args[a])); - } else if ((args[a] instanceof short[])) { + } else if (args[a] instanceof short[]) { retVal.append(Arrays.toString((short[]) args[a])); - } else if ((args[a] instanceof byte[])) { + } else if (args[a] instanceof byte[]) { retVal.append(Arrays.toString((byte[]) args[a])); - } else if ((args[a] instanceof boolean[])) { + } else if (args[a] instanceof boolean[]) { retVal.append(Arrays.toString((boolean[]) args[a])); - } else if ((args[a] instanceof char[])) { + } else if (args[a] instanceof char[]) { retVal.append(Arrays.toString((char[]) args[a])); } else { retVal.append(args[a]); @@ -196,18 +194,18 @@ static boolean isSameDate(final Calendar aCal1, final Calendar aCal2) { boolean retVal = aCal1.get(Calendar.YEAR) == aCal2.get(Calendar.YEAR); - retVal = retVal && (aCal1.get(Calendar.MONTH) == aCal2.get(Calendar.MONTH)); + retVal = retVal && aCal1.get(Calendar.MONTH) == aCal2.get(Calendar.MONTH); - return retVal && (aCal1.get(Calendar.DAY_OF_MONTH) == aCal2.get(Calendar.DAY_OF_MONTH)); + return retVal && aCal1.get(Calendar.DAY_OF_MONTH) == aCal2.get(Calendar.DAY_OF_MONTH); } static boolean isSameTime(final Calendar aCal1, final Calendar aCal2) { boolean retVal = aCal1.get(Calendar.HOUR_OF_DAY) == aCal2.get(Calendar.HOUR_OF_DAY); - retVal = retVal && (aCal1.get(Calendar.MINUTE) == aCal2.get(Calendar.MINUTE)); + retVal = retVal && aCal1.get(Calendar.MINUTE) == aCal2.get(Calendar.MINUTE); - return retVal && (aCal1.get(Calendar.SECOND) == aCal2.get(Calendar.SECOND)); + return retVal && aCal1.get(Calendar.SECOND) == aCal2.get(Calendar.SECOND); } protected TypeUtils() { diff --git a/src/test/java/org/ojalgo/TestUtils.java b/src/test/java/org/ojalgo/TestUtils.java index 48047b829..ed0d72c1b 100644 --- a/src/test/java/org/ojalgo/TestUtils.java +++ b/src/test/java/org/ojalgo/TestUtils.java @@ -89,6 +89,17 @@ public static void assertBounds(final Comparable lower, final Comparable v } } + public static void assertComplexEquals(final Access1D expected, final Access1D actual) { + TestUtils.assertComplexEquals(expected, actual, EQUALS); + } + + public static void assertComplexEquals(final Access1D expected, final Access1D actual, final NumberContext context) { + TestUtils.assertEquals(expected, actual, context); + for (int i = 0; i < expected.size(); i++) { + TestUtils.assertEquals(expected.get(i), actual.get(i), context); + } + } + public static void assertEquals(final Access1D expected, final Access1D actual) { TestUtils.assertEquals(expected, actual, EQUALS); } @@ -395,6 +406,12 @@ public static void assertInRange(final int first, final int limit, final int act } } + public static void assertInRange(final double first, final double limit, final double actual) { + if (first > actual || actual >= limit) { + TestUtils.fail("Not in range!"); + } + } + public static void assertLessThan(final double reference, final double actual) { if (actual >= reference) { Assertions.fail(actual + " !< " + reference); diff --git a/src/test/java/org/ojalgo/data/transform/DFTBenchmark.java b/src/test/java/org/ojalgo/data/transform/DFTBenchmark.java new file mode 100644 index 000000000..a035d1921 --- /dev/null +++ b/src/test/java/org/ojalgo/data/transform/DFTBenchmark.java @@ -0,0 +1,74 @@ +package org.ojalgo.data.transform; + +import org.ojalgo.BenchmarkUtils; +import org.ojalgo.function.special.PowerOf2; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.matrix.store.Primitive64Store; +import org.ojalgo.random.Uniform; +import org.ojalgo.scalar.ComplexNumber; +import org.openjdk.jmh.annotations.Benchmark; +import org.openjdk.jmh.annotations.Param; +import org.openjdk.jmh.annotations.Scope; +import org.openjdk.jmh.annotations.Setup; +import org.openjdk.jmh.annotations.State; +import org.openjdk.jmh.runner.RunnerException; + +/** + * 
+Benchmark            (power)   Mode  Cnt           Score           Error    Units
+DFTBenchmark.fft           0  thrpt    3  3124513790.555 ± 305627548.431  ops/min
+DFTBenchmark.fft           1  thrpt    3  2398675511.894 ± 249205811.045  ops/min
+DFTBenchmark.fft           2  thrpt    3  1718574817.845 ±   4230212.102  ops/min
+DFTBenchmark.fft           3  thrpt    3   960108856.886 ±   1458412.598  ops/min
+DFTBenchmark.fft           4  thrpt    3   491959971.868 ±   5637718.754  ops/min
+DFTBenchmark.fft           5  thrpt    3   240706104.224 ±   1040731.406  ops/min
+DFTBenchmark.fft           6  thrpt    3   111525835.267 ±    421865.235  ops/min
+DFTBenchmark.fft           7  thrpt    3    51483058.691 ±    182468.393  ops/min
+DFTBenchmark.fft           8  thrpt    3    23777894.972 ±     73809.368  ops/min
+DFTBenchmark.matrix        0  thrpt    3  1576229925.659 ±  12786183.061  ops/min
+DFTBenchmark.matrix        1  thrpt    3   982580395.880 ±   2184563.826  ops/min
+DFTBenchmark.matrix        2  thrpt    3   535794233.758 ±   5384707.106  ops/min
+DFTBenchmark.matrix        3  thrpt    3   185259013.686 ±   2036997.500  ops/min
+DFTBenchmark.matrix        4  thrpt    3    58275950.781 ±    585827.402  ops/min
+DFTBenchmark.matrix        5  thrpt    3    15247998.933 ±    142923.255  ops/min
+DFTBenchmark.matrix        6  thrpt    3     4214174.031 ±    173260.531  ops/min
+DFTBenchmark.matrix        7  thrpt    3     1071304.704 ±     26403.251  ops/min
+DFTBenchmark.matrix        8  thrpt    3      262046.213 ±      4266.517  ops/min
+ *
+ */ +@State(Scope.Benchmark) +public class DFTBenchmark { + + public static void main(final String[] args) throws RunnerException { + BenchmarkUtils.run(DFTBenchmark.class); + } + + @Param({ "0", "1", "2", "3", "4", "5", "6", "7", "8" }) + public int power; + + DiscreteFourierTransform.FFT fft; + Primitive64Store input; + DiscreteFourierTransform.FullMatrix matrix; + + + + @Benchmark + public MatrixStore fft() { + return fft.transform(input.data); + } + + @Benchmark + public MatrixStore matrix() { + return matrix.transform(input.data); + } + + @Setup + public void setup() { + int size = PowerOf2.powerOfInt2(power); + input = Primitive64Store.FACTORY.makeFilled(size, 1, Uniform.of(-2, 4)); + fft = new DiscreteFourierTransform.FFT(input.size()); + matrix = new DiscreteFourierTransform.FullMatrix(input.size()); + + } + +} diff --git a/src/test/java/org/ojalgo/data/transform/DataTransformTests.java b/src/test/java/org/ojalgo/data/transform/DataTransformTests.java new file mode 100644 index 000000000..63eccc5e3 --- /dev/null +++ b/src/test/java/org/ojalgo/data/transform/DataTransformTests.java @@ -0,0 +1,28 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +abstract class DataTransformTests { + + static final boolean DEBUG = false; + +} diff --git a/src/test/java/org/ojalgo/data/transform/DiscreteFourierTransformTest.java b/src/test/java/org/ojalgo/data/transform/DiscreteFourierTransformTest.java new file mode 100644 index 000000000..6d7f44c2f --- /dev/null +++ b/src/test/java/org/ojalgo/data/transform/DiscreteFourierTransformTest.java @@ -0,0 +1,469 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +import org.junit.jupiter.api.Test; +import org.ojalgo.TestUtils; +import org.ojalgo.function.constant.ComplexMath; +import org.ojalgo.function.constant.PrimitiveMath; +import org.ojalgo.function.polynomial.PolynomialR064; +import org.ojalgo.function.special.PowerOf2; +import org.ojalgo.matrix.MatrixC128; +import org.ojalgo.matrix.MatrixC128.DenseReceiver; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.matrix.store.PhysicalStore; +import org.ojalgo.matrix.store.Primitive64Store; +import org.ojalgo.netio.BasicLogger; +import org.ojalgo.random.Uniform; +import org.ojalgo.scalar.ComplexNumber; +import org.ojalgo.structure.Access1D; +import org.ojalgo.type.context.NumberContext; + +final class DiscreteFourierTransformTest extends DataTransformTests { + + @Test + public void testBitReversal() { + + TestUtils.assertEquals(new int[] { 0 }, DiscreteFourierTransform.getBitReversedIndices(1)); + + TestUtils.assertEquals(new int[] { 0, 1 }, DiscreteFourierTransform.getBitReversedIndices(2)); + + TestUtils.assertEquals(new int[] { 0, 2, 1, 3 }, DiscreteFourierTransform.getBitReversedIndices(4)); + + TestUtils.assertEquals(new int[] { 0, 4, 2, 6, 1, 5, 3, 7 }, DiscreteFourierTransform.getBitReversedIndices(8)); + + TestUtils.assertEquals(new int[] { 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 }, DiscreteFourierTransform.getBitReversedIndices(16)); + } + + @Test + public void testCompareImplementationsUsingRandomInput() { + + NumberContext accuracy = NumberContext.of(8); + + for (int power = 0; power <= 8; power++) { + + int dim = PowerOf2.powerOfInt2(power); + + PhysicalStore input = Primitive64Store.FACTORY.makeFilled(dim, 1, Uniform.of(-2, 4)); + + DiscreteFourierTransform full = new DiscreteFourierTransform.FullMatrix(dim); + DiscreteFourierTransform fast = new DiscreteFourierTransform.FFT(dim); + + MatrixStore expected = full.transform(input); + MatrixStore actual = fast.transform(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Full", expected); + BasicLogger.debugMatrix("Fast", actual); + } + TestUtils.assertComplexEquals(expected, actual, accuracy); + + expected = full.inverse(expected); + actual = fast.inverse(actual); + + if (DEBUG) { + BasicLogger.debugMatrix("Input", input); + BasicLogger.debugMatrix("Full", expected); + BasicLogger.debugMatrix("Fast", actual); + } + TestUtils.assertEquals(input, expected, accuracy); + TestUtils.assertEquals(input, actual, accuracy); + } + } + + /** + * https://www.youtube.com/watch?v=x3QxJnI9jNI + */ + @Test + public void testDrUnderwoodsPhysicsYouTubePage() { + + // y = f(x) + double[] x = { 0, 2, 4, 6 }; + double[] y = { 1, 4, 3, 2 }; + + double dx = 2; + int N = 4; + + ComplexNumber[] amplitudes = new ComplexNumber[N]; + double[] frequencies = new double[N]; // angular frequencies, p + + double angularIncrement = PrimitiveMath.TWO_PI / N; + + for (int n = 0; n < frequencies.length; n++) { + frequencies[n] = n * PrimitiveMath.TWO_PI / (N * dx); + } + + ComplexNumber[] expectedAmplitudes = new ComplexNumber[N]; + double[] expectedAmplitudeNorms = new double[N]; + expectedAmplitudes[0] = ComplexNumber.of(10, 0); + expectedAmplitudes[1] = ComplexNumber.of(-2, -2); + expectedAmplitudes[2] = ComplexNumber.of(-2, 0); + expectedAmplitudes[3] = ComplexNumber.of(-2, 2); + expectedAmplitudeNorms[0] = 10.0; + expectedAmplitudeNorms[1] = 2.0 * PrimitiveMath.SQRT_TWO; + expectedAmplitudeNorms[2] = 2.0; + expectedAmplitudeNorms[3] = 2.0 * PrimitiveMath.SQRT_TWO; + + // Directly + + for (int n = 0; n < amplitudes.length; n++) { + + ComplexNumber amplitude = ComplexNumber.ZERO; + + for (int r = 0; r < N; r++) { + amplitude = amplitude.add(ComplexNumber.makePolar(y[r], -angularIncrement * n * r)); + } + + amplitudes[n] = amplitude; + } + + for (int n = 0; n < amplitudes.length; n++) { + if (DEBUG) { + BasicLogger.debug("norm({})={}", amplitudes[n], amplitudes[n].norm()); + } + TestUtils.assertEquals(expectedAmplitudes[n], amplitudes[n]); + TestUtils.assertEquals(expectedAmplitudeNorms[n], amplitudes[n].norm()); + } + + // Using ComplexNumber.newUnitRoot(N) + + ComplexNumber w = ComplexNumber.newUnitRoot(N); + + for (int n = 0; n < amplitudes.length; n++) { + + ComplexNumber amplitude = ComplexNumber.ZERO; + + for (int r = 0; r < N; r++) { + amplitude = amplitude.add(w.power(n * r).multiply(y[r])); + } + + amplitudes[n] = amplitude; + } + + for (int n = 0; n < amplitudes.length; n++) { + if (DEBUG) { + BasicLogger.debug("norm({})={}", amplitudes[n], amplitudes[n].norm()); + } + TestUtils.assertEquals(expectedAmplitudes[n], amplitudes[n]); + TestUtils.assertEquals(expectedAmplitudeNorms[n], amplitudes[n].norm()); + } + + // Using DiscreteFourierTransform + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(y.length); + + MatrixStore actual = transformer.transform(y); + + for (int n = 0; n < amplitudes.length; n++) { + ComplexNumber ampl = actual.get(n); + if (DEBUG) { + BasicLogger.debug("norm({})={}", ampl, ampl.norm()); + } + TestUtils.assertEquals(expectedAmplitudes[n], ampl); + TestUtils.assertEquals(expectedAmplitudeNorms[n], ampl.norm()); + } + } + + /** + * https://www.youtube.com/watch?v=hsi9-x1Ts4Y + */ + @Test + public void testECAcademy() { + + MatrixC128 input = MatrixC128.FACTORY.column(0.0, 1.0, 2.0, 3.0); + + DenseReceiver receiver = MatrixC128.FACTORY.makeDense(4, 1); + receiver.set(0, ComplexNumber.of(6, 0)); + receiver.set(1, ComplexNumber.of(-2, 2)); + receiver.set(2, ComplexNumber.of(-2, 0)); + receiver.set(3, ComplexNumber.of(-2, -2)); + MatrixC128 expected = receiver.get(); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(input.size()); + + MatrixStore actual = transformer.transform(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Expected", expected); + BasicLogger.debugMatrix("Actual", actual); + } + TestUtils.assertComplexEquals(expected, actual); + + MatrixStore inverse = transformer.inverse(actual); + if (DEBUG) { + BasicLogger.debugMatrix("Input", input); + BasicLogger.debugMatrix("Inverse", inverse); + } + TestUtils.assertComplexEquals(input, inverse); + + } + + /** + * https://www.youtube.com/watch?v=RHXvrzH0XvA + */ + @Test + public void testEnggCourseMadeEasy() { + + MatrixC128 input = MatrixC128.FACTORY.column(1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0); + + DenseReceiver receiver = MatrixC128.FACTORY.makeDense(8, 1); + receiver.set(0, ComplexNumber.of(4, 0)); + receiver.set(1, ComplexNumber.of(1, -2.4142135623730945)); + receiver.set(2, ComplexNumber.of(0, 0)); + receiver.set(3, ComplexNumber.of(1, -0.4142135623730945)); + receiver.set(4, ComplexNumber.of(0, 0)); + receiver.set(5, ComplexNumber.of(1, 0.4142135623730945)); + receiver.set(6, ComplexNumber.of(0, 0)); + receiver.set(7, ComplexNumber.of(1, 2.4142135623730945)); + MatrixC128 expected = receiver.get(); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(input.size()); + + MatrixStore actual = transformer.transform(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Expected", expected); + BasicLogger.debugMatrix("Actual", actual); + } + TestUtils.assertComplexEquals(expected, actual); + + MatrixStore inverse = transformer.inverse(actual); + if (DEBUG) { + BasicLogger.debugMatrix("Input", input); + BasicLogger.debugMatrix("Inverse", inverse); + } + TestUtils.assertComplexEquals(input, inverse); + } + + @Test + public void testGetUnitRoots() { + + for (int e = 9; e >= 0; e--) { + + int powerOf2 = PowerOf2.powerOfInt2(e); + + ComplexNumber[] expected = ComplexNumber.newUnitRoots(powerOf2); + Access1D actual = DiscreteFourierTransform.getUnitRoots(powerOf2); + + for (int i = 0; i < powerOf2; i++) { + TestUtils.assertEquals(expected[i], actual.get(i)); + } + } + } + + /** + * Primarily tests there are no obvious problems using dimensions not power of 2. + */ + @Test + public void testNonPowerOf2() { + + NumberContext accuracy = NumberContext.of(8); + + for (int dimension = 2; dimension <= 16; dimension++) { + + PhysicalStore input = Primitive64Store.FACTORY.makeFilled(dimension, 1, Uniform.of(-2, 4)); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(dimension); + + MatrixStore transformed = transformer.transform(input); + MatrixStore reverted = transformer.inverse(transformed); + + if (DEBUG) { + BasicLogger.debugMatrix("Input", input); + BasicLogger.debugMatrix("Transformed", transformed); + BasicLogger.debugMatrix("Reverted", reverted); + } + + TestUtils.assertEquals(input, reverted, accuracy); + } + } + + @Test + public void testShiftAndRevertEven() { + + for (int c = 2; c < 10; c += 2) { + for (int r = 2; r < 10; r += 2) { + + Primitive64Store original = Primitive64Store.FACTORY.make(r, c); + + for (int j = 0; j < c; j++) { + for (int i = 0; i < r; i++) { + original.set(i, j, i + j); + } + } + + if (DEBUG) { + BasicLogger.debugMatrix("Original " + r + "x" + c, original); + } + + MatrixStore shifted = DiscreteFourierTransform.shift(original); + + if (DEBUG) { + BasicLogger.debugMatrix("Shifted " + r + "x" + c, shifted); + } + + int cr = r / 2; + int cc = c / 2; + + TestUtils.assertEquals("(" + cr + "," + cc + ")", PrimitiveMath.ZERO, shifted.doubleValue(cr, cc)); + + MatrixStore reverted = DiscreteFourierTransform.shift(shifted); + + if (DEBUG) { + BasicLogger.debugMatrix("Reverted " + r + "x" + c, reverted); + } + + TestUtils.assertEquals(original, reverted); + } + } + } + + @Test + public void testShiftCenterEvenAndOdd() { + + for (int c = 3; c < 10; c++) { + for (int r = 3; r < 10; r++) { + + Primitive64Store original = Primitive64Store.FACTORY.make(r, c); + + for (int j = 0; j < c; j++) { + for (int i = 0; i < r; i++) { + original.set(i, j, i + j); + } + } + + if (DEBUG) { + BasicLogger.debugMatrix("Original " + r + "x" + c, original); + } + + MatrixStore shifted = DiscreteFourierTransform.shift(original); + + if (DEBUG) { + BasicLogger.debugMatrix("Shifted " + r + "x" + c, shifted); + } + + int cr = r / 2; + int cc = c / 2; + + TestUtils.assertEquals("(" + cr + "," + cc + ")", PrimitiveMath.ZERO, shifted.doubleValue(cr, cc)); + } + } + } + + /** + * https://www.youtube.com/watch?v=FaWSGmkboOs + */ + @Test + public void testSmartEngineer() { + + MatrixC128 input = MatrixC128.FACTORY.column(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0); + + DenseReceiver receiver = MatrixC128.FACTORY.makeDense(8, 1); + receiver.set(0, ComplexNumber.of(28, 0)); + receiver.set(1, ComplexNumber.of(-4, 9.656854249492376)); + receiver.set(2, ComplexNumber.of(-4, 4)); + receiver.set(3, ComplexNumber.of(-4, 1.656854249492376)); + receiver.set(4, ComplexNumber.of(-4, 0)); + receiver.set(5, ComplexNumber.of(-4, -1.656854249492376)); + receiver.set(6, ComplexNumber.of(-4, -4)); + receiver.set(7, ComplexNumber.of(-4, -9.656854249492376)); + MatrixC128 expected = receiver.get(); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(input.size()); + + MatrixStore actual = transformer.transform(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Expected", expected); + BasicLogger.debugMatrix("Actual", actual); + } + TestUtils.assertComplexEquals(expected, actual); + + MatrixStore inverse = transformer.inverse(actual); + if (DEBUG) { + BasicLogger.debugMatrix("Input", input); + BasicLogger.debugMatrix("Inverse", inverse); + } + TestUtils.assertComplexEquals(input, inverse); + } + + /** + * https://cs.stackexchange.com/questions/16266/show-how-to-do-fft-by-hand + */ + @Test + public void testStackExchangePolynomialMultiplication() { + + PolynomialR064 polynomialA = PolynomialR064.wrap(3, 1, 0, 0); + PolynomialR064 polynomialB = PolynomialR064.wrap(2, 0, 2, 0); + PolynomialR064 polynomialC = PolynomialR064.wrap(6, 2, 6, 2); + + TestUtils.assertEquals(polynomialC, polynomialA.multiply(polynomialB)); + + ComplexNumber unitRoot = ComplexNumber.newUnitRoot(4); + double unitRootPhase = unitRoot.phase(); + + DenseReceiver receiver = MatrixC128.FACTORY.makeDense(4, 4); + for (int j = 0; j < 4; j++) { + for (int i = 0; i < 4; i++) { + receiver.set(i, j, ComplexNumber.makePolar(PrimitiveMath.ONE, i * j * unitRootPhase)); + } + } + MatrixC128 matrix = receiver.get(); + + if (DEBUG) { + BasicLogger.debugMatrix("Transformation matrix", matrix); + } + TestUtils.assertComplexEquals(matrix, DiscreteFourierTransform.newVandermondeMatrix(4)); + + DiscreteFourierTransform transform = DiscreteFourierTransform.newInstance(4); + + MatrixC128 columnA = MatrixC128.FACTORY.columns(polynomialA); + MatrixC128 transfA = matrix.multiply(columnA); + MatrixC128 columnB = MatrixC128.FACTORY.columns(polynomialB); + MatrixC128 transfB = matrix.multiply(columnB); + + MatrixC128 transfC = transfA.onMatching(ComplexMath.MULTIPLY, transfB); + + if (DEBUG) { + BasicLogger.debugMatrix("Transformed A", transfA); + BasicLogger.debugMatrix("Transformed B", transfB); + BasicLogger.debugMatrix("Transformed C", transfC); + } + TestUtils.assertComplexEquals(transfA, transform.transform(columnA)); + TestUtils.assertComplexEquals(transfB, transform.transform(columnB)); + + MatrixC128 inverse = matrix.invert(); + + if (DEBUG) { + BasicLogger.debugMatrix("Inverse transformation", inverse); + } + + MatrixC128 columnC = inverse.multiply(transfC); + + if (DEBUG) { + BasicLogger.debugMatrix("Column C", columnC); + } + TestUtils.assertEquals(polynomialC, columnC); + TestUtils.assertEquals(polynomialC, transform.inverse(transfC)); + } + +} diff --git a/src/test/java/org/ojalgo/data/transform/ZtransformTest.java b/src/test/java/org/ojalgo/data/transform/ZtransformTest.java new file mode 100644 index 000000000..4b0043d0c --- /dev/null +++ b/src/test/java/org/ojalgo/data/transform/ZtransformTest.java @@ -0,0 +1,132 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.data.transform; + +import org.junit.jupiter.api.Test; +import org.ojalgo.TestUtils; +import org.ojalgo.array.ArrayR064; +import org.ojalgo.function.special.PowerOf2; +import org.ojalgo.matrix.MatrixC128; +import org.ojalgo.matrix.store.MatrixStore; +import org.ojalgo.matrix.store.PhysicalStore; +import org.ojalgo.matrix.store.Primitive64Store; +import org.ojalgo.netio.BasicLogger; +import org.ojalgo.random.Uniform; +import org.ojalgo.scalar.ComplexNumber; +import org.ojalgo.type.context.NumberContext; + +public class ZtransformTest extends DataTransformTests { + + @Test + public void testCompareDFTUsingRandomInput() { + + NumberContext accuracy = NumberContext.of(8); + + for (int exp = 1; exp <= 8; exp++) { + + int dim = PowerOf2.powerOfInt2(exp); + + PhysicalStore input = Primitive64Store.FACTORY.makeFilled(dim, 1, Uniform.of(-2, 4)); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(input.size()); + MatrixStore expected = transformer.transform(input); + + MatrixStore actual = ZTransform.doDFT(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Expected", expected); + BasicLogger.debugMatrix("Actual", actual); + } + TestUtils.assertComplexEquals(expected, actual, accuracy); + } + } + + /** + * https://www.youtube.com/watch?v=B4IyRw1zvvA&list=PLJ8LTUMGG9U6FGNmnLscY1fhSA7FiFrfO + */ + @Test + public void testDavidDorran() { + + // + + ArrayR064 sequence = ArrayR064.wrap(2.0, -3.0, -1.0, 4.0); + + ComplexNumber z = ComplexNumber.valueOf(1.0); + + ComplexNumber expected = ComplexNumber.valueOf(2.0); + + TestUtils.assertEquals(expected, ZTransform.doTransform(sequence, z)); + + // + + sequence = ArrayR064.wrap(3.0, -3.0, 3.0, -3.0); + + z = ComplexNumber.valueOf(1.0); + + expected = ComplexNumber.valueOf(0.0); + + TestUtils.assertEquals(expected, ZTransform.doTransform(sequence, z)); + + // + + sequence = ArrayR064.wrap(0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); + + z = ComplexNumber.of(0.3, -0.3); + + expected = ComplexNumber.of(1.333333333333333, 0.8333333333333333); + + TestUtils.assertEquals(expected, ZTransform.doTransform(sequence, z)); + + // + + sequence = ArrayR064.wrap(0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0); + + z = ComplexNumber.valueOf(-1.0); + + expected = ComplexNumber.of(0.0, 0.0); + + TestUtils.assertEquals(expected, ZTransform.doTransform(sequence, z)); + } + + /** + * https://www.youtube.com/watch?v=FaWSGmkboOs + * + * @see DiscreteFourierTransformTest#testSmartEngineer() + */ + @Test + public void testSmartEngineer() { + + MatrixC128 input = MatrixC128.FACTORY.column(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0); + + DiscreteFourierTransform transformer = DiscreteFourierTransform.newInstance(input.size()); + MatrixStore expected = transformer.transform(input); + + MatrixStore actual = ZTransform.doDFT(input); + + if (DEBUG) { + BasicLogger.debugMatrix("Expected", expected); + BasicLogger.debugMatrix("Actual", actual); + } + TestUtils.assertComplexEquals(expected, actual); + } + +} diff --git a/src/test/java/org/ojalgo/function/polynomial/PolynomialImplTest.java b/src/test/java/org/ojalgo/function/polynomial/PolynomialImplTest.java index b37fe51b7..73915ed35 100755 --- a/src/test/java/org/ojalgo/function/polynomial/PolynomialImplTest.java +++ b/src/test/java/org/ojalgo/function/polynomial/PolynomialImplTest.java @@ -30,6 +30,20 @@ public class PolynomialImplTest { + @Test + public void testAddPolynomials() { + + PolynomialR064 degree3 = PolynomialR064.wrap(1D, 1D, 1D, 1D); + PolynomialR064 degree2 = PolynomialR064.wrap(1D, 1D, 1D); + PolynomialR064 degree1 = PolynomialR064.wrap(1D, 1D); + PolynomialR064 degree0 = PolynomialR064.wrap(1D); + + PolynomialFunction expected = PolynomialR064.wrap(4D, 3D, 2D, 1D); + PolynomialFunction actual = degree3.add(degree2).add(degree1).add(degree0); + + TestUtils.assertEquals(expected, actual); + } + @Test public void testEstimation() { @@ -81,8 +95,33 @@ public void testEvaluation() { tmpPoly.set(2, 10.0); for (double i = -100.0; i <= 100; i = i + 10.0) { - TestUtils.assertEquals(5.0 + i + (10.0 * (i * i)), tmpPoly.invoke(i), 1E-14 / PrimitiveMath.THREE); + TestUtils.assertEquals(5.0 + i + 10.0 * (i * i), tmpPoly.invoke(i), 1E-14 / PrimitiveMath.THREE); } } + @Test + public void testPower() { + + PolynomialR064 degree3 = PolynomialR064.wrap(1D, 1D, 1D, 1D); + + TestUtils.assertEquals(3, degree3.degree()); + + TestUtils.assertEquals(0, degree3.power(0).degree()); + TestUtils.assertEquals(3, degree3.power(1).degree()); + TestUtils.assertEquals(6, degree3.power(2).degree()); + TestUtils.assertEquals(9, degree3.power(3).degree()); + TestUtils.assertEquals(12, degree3.power(4).degree()); + TestUtils.assertEquals(15, degree3.power(5).degree()); + + degree3 = PolynomialR064.wrap(2D, 1D, 4D, 3D); + + PolynomialFunction expected = degree3.multiply(degree3); + expected = expected.multiply(expected); + expected = expected.multiply(expected); + + PolynomialFunction actual = degree3.power(8); + + TestUtils.assertEquals(expected, actual); + } + } diff --git a/src/test/java/org/ojalgo/function/series/FourierSeriesTest.java b/src/test/java/org/ojalgo/function/series/FourierSeriesTest.java new file mode 100644 index 000000000..e62fedbdc --- /dev/null +++ b/src/test/java/org/ojalgo/function/series/FourierSeriesTest.java @@ -0,0 +1,203 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.series; + +import static org.ojalgo.function.constant.PrimitiveMath.*; + +import java.util.function.DoubleUnaryOperator; + +import org.junit.jupiter.api.Test; +import org.ojalgo.TestUtils; +import org.ojalgo.function.PrimitiveFunction.SampleDomain; +import org.ojalgo.netio.BasicLogger; +import org.ojalgo.type.context.NumberContext; + +public class FourierSeriesTest extends SeriesFunctionTests { + + private static final NumberContext ACCURACY = NumberContext.of(2); + + private static final double INITIAL = ZERO - PI - E; + + private static final double NUDGE = TEN * MACHINE_EPSILON; + + private static final SampleDomain SAMPLE = new SampleDomain(TWO * TWO_PI, 16); + + private static boolean isCloseTooAnyMultipleOfPI(final double arg) { + double cantidate = arg; + while (cantidate < PI) { + cantidate += PI; + } + while (cantidate > TWO_PI) { + cantidate -= PI; + } + cantidate = Math.min(cantidate - PI, TWO_PI - cantidate); + return cantidate < NINTH; + } + + static void doTest(final DoubleUnaryOperator function, final SampleDomain period) { + + FourierSeries series = FourierSeries.estimate(function, period); + + for (double arg = INITIAL; arg < TEN; arg += NINTH) { + double expected = function.applyAsDouble(arg); + double actual = series.invoke(arg); + if (DEBUG) { + BasicLogger.debug("f({}) = {} & {}", arg, expected, actual); + } + if (!FourierSeriesTest.isCloseTooAnyMultipleOfPI(arg)) { + // This is to avoid the discontinuities some functions have at multiples of PI. + TestUtils.assertEquals("" + arg, expected, actual, ACCURACY); + } + } + } + + static void doTest(final PeriodicFunction function, final int nbSamples) { + FourierSeriesTest.doTest(function, function.getSampleDomain(nbSamples)); + } + + @Test + public void testCos() { + FourierSeriesTest.doTest(Math::cos, SAMPLE); + } + + /** + * https://www.math.kit.edu/iana3/lehre/fourierana2014w/media/fstable141127.pdf + */ + @Test + public void testTableItem01() { + + PeriodicFunction function = PeriodicFunction.ofCentered(Math::abs); + + TestUtils.assertEquals(PI, function.applyAsDouble(-PI)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(-HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(ZERO)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(HALF_PI)); + TestUtils.assertEquals(PI, function.applyAsDouble(PI)); + + // Next period + TestUtils.assertEquals(PI, function.applyAsDouble(PI)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(TWO_PI - HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(TWO_PI)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(TWO_PI + HALF_PI)); + TestUtils.assertEquals(PI, function.applyAsDouble(TWO_PI + PI)); + + // Previous period + TestUtils.assertEquals(PI, function.applyAsDouble(-TWO_PI - PI)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(-TWO_PI - HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(-TWO_PI)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(-TWO_PI + HALF_PI)); + TestUtils.assertEquals(PI, function.applyAsDouble(-PI)); + + FourierSeriesTest.doTest(function, 64); + } + + /** + * https://www.math.kit.edu/iana3/lehre/fourierana2014w/media/fstable141127.pdf + */ + @Test + public void testTableItem02() { + + PeriodicFunction function = PeriodicFunction.ofCentered(DoubleUnaryOperator.identity()); + + TestUtils.assertEquals(-PI, function.applyAsDouble(-PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(ZERO)); + TestUtils.assertEquals(PI, function.applyAsDouble(PI - NUDGE)); + + FourierSeriesTest.doTest(function, 256); + } + + /** + * https://www.math.kit.edu/iana3/lehre/fourierana2014w/media/fstable141127.pdf + */ + @Test + public void testTableItem03() { + + PeriodicFunction function = PeriodicFunction.of(arg -> PI - arg); + + TestUtils.assertEquals(PI, function.applyAsDouble(ZERO)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(PI)); + TestUtils.assertEquals(-PI, function.applyAsDouble(TWO_PI - NUDGE)); + + FourierSeriesTest.doTest(function, 256); + } + + /** + * https://www.math.kit.edu/iana3/lehre/fourierana2014w/media/fstable141127.pdf + */ + @Test + public void testTableItem04() { + + PeriodicFunction function = PeriodicFunction.ofCentered(arg -> arg <= ZERO ? ZERO : arg); + + TestUtils.assertEquals(ZERO, function.applyAsDouble(-PI + NUDGE)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(-HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(ZERO)); + TestUtils.assertEquals(HALF_PI, function.applyAsDouble(HALF_PI)); + TestUtils.assertEquals(PI, function.applyAsDouble(PI - NUDGE)); + + FourierSeriesTest.doTest(function, 4096); + } + + /** + * https://www.math.kit.edu/iana3/lehre/fourierana2014w/media/fstable141127.pdf + */ + @Test + public void testTableItem05() { + + PeriodicFunction function = PeriodicFunction.ofCentered(arg -> Math.sin(arg) * Math.sin(arg)); + DoubleUnaryOperator series = arg -> 0.5 - 0.5 * Math.cos(2 * arg); + + TestUtils.assertEquals(ZERO, function.applyAsDouble(-PI)); + TestUtils.assertEquals(ONE, function.applyAsDouble(-HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(ZERO)); + TestUtils.assertEquals(ONE, function.applyAsDouble(HALF_PI)); + TestUtils.assertEquals(ZERO, function.applyAsDouble(PI - NUDGE)); + + FourierSeriesTest.doTest(function, 8); + } + + @Test + public void testSawtooth() { + FourierSeriesTest.doTest(PeriodicFunction.SAWTOOTH, 256); + } + + @Test + public void testSin() { + FourierSeriesTest.doTest(Math::sin, SAMPLE); + } + + @Test + public void testSine() { + FourierSeriesTest.doTest(PeriodicFunction.SINE, 8); + } + + @Test + public void testSquare() { + FourierSeriesTest.doTest(PeriodicFunction.SQUARE, 128); + } + + @Test + public void testTriangle() { + FourierSeriesTest.doTest(PeriodicFunction.TRIANGLE, 32); + } + +} diff --git a/src/test/java/org/ojalgo/function/series/PeriodicFunctionTest.java b/src/test/java/org/ojalgo/function/series/PeriodicFunctionTest.java new file mode 100644 index 000000000..75b0504cd --- /dev/null +++ b/src/test/java/org/ojalgo/function/series/PeriodicFunctionTest.java @@ -0,0 +1,132 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.series; + +import static org.ojalgo.function.constant.PrimitiveMath.*; + +import java.util.function.DoubleUnaryOperator; + +import org.junit.jupiter.api.Test; +import org.ojalgo.TestUtils; + +public class PeriodicFunctionTest extends SeriesFunctionTests { + + private static final double NUDGE = TEN * MACHINE_EPSILON; + + @Test + public void testOriginAndPeriod() { + + DoubleUnaryOperator shape = arg -> arg >= FOUR ? ONE : ZERO; + + PeriodicFunction function = new PeriodicFunction(THREE, shape, TWO); + + TestUtils.assertEquals(ZERO, function.invoke(3.5)); + TestUtils.assertEquals(ONE, function.invoke(4.5)); + + TestUtils.assertEquals(ZERO, function.invoke(1.5)); + TestUtils.assertEquals(ONE, function.invoke(2.5)); + + TestUtils.assertEquals(ZERO, function.invoke(5.5)); + TestUtils.assertEquals(ONE, function.invoke(6.5)); + + TestUtils.assertEquals(ZERO, function.invoke(-0.5)); + TestUtils.assertEquals(ONE, function.invoke(0.5)); + + TestUtils.assertEquals(ZERO, function.invoke(-2.5)); + TestUtils.assertEquals(ONE, function.invoke(-1.5)); + + TestUtils.assertEquals(ZERO, function.invoke(7.5)); + TestUtils.assertEquals(ONE, function.invoke(8.5)); + } + + @Test + public void testSawtooth() { + + PeriodicFunction function = PeriodicFunction.SAWTOOTH; + + for (int i = -3; i <= 3; i++) { + + double shift = TWO_PI * i; + + TestUtils.assertEquals(ZERO, function.invoke(shift + ZERO)); + TestUtils.assertEquals(HALF, function.invoke(shift + HALF_PI)); + TestUtils.assertEquals(ONE, function.invoke(shift + PI - NUDGE)); + TestUtils.assertEquals(NEG, function.invoke(shift + PI + NUDGE)); + TestUtils.assertEquals(-HALF, function.invoke(shift + PI + HALF_PI)); + TestUtils.assertEquals(ZERO, function.invoke(shift + TWO_PI)); + } + } + + @Test + public void testSine() { + + PeriodicFunction function = PeriodicFunction.SINE; + + for (int i = -3; i <= 3; i++) { + + double shift = TWO_PI * i; + + TestUtils.assertEquals(Math.sin(ZERO), function.invoke(shift + ZERO)); + TestUtils.assertEquals(Math.sin(HALF_PI), function.invoke(shift + HALF_PI)); + TestUtils.assertEquals(Math.sin(PI), function.invoke(shift + PI)); + TestUtils.assertEquals(Math.sin(PI + HALF_PI), function.invoke(shift + PI + HALF_PI)); + TestUtils.assertEquals(Math.sin(TWO_PI), function.invoke(shift + TWO_PI)); + } + } + + @Test + public void testSquare() { + + PeriodicFunction function = PeriodicFunction.SQUARE; + + for (int i = -3; i <= 3; i++) { + + double shift = TWO_PI * i; + + TestUtils.assertEquals(ONE, function.invoke(shift + ZERO + NUDGE)); + TestUtils.assertEquals(ONE, function.invoke(shift + HALF_PI)); + TestUtils.assertEquals(ONE, function.invoke(shift + PI - NUDGE)); + TestUtils.assertEquals(NEG, function.invoke(shift + PI + NUDGE)); + TestUtils.assertEquals(NEG, function.invoke(shift + PI + HALF_PI)); + TestUtils.assertEquals(NEG, function.invoke(shift + TWO_PI - NUDGE)); + } + } + + + @Test + public void testTriangle() { + + PeriodicFunction function = PeriodicFunction.TRIANGLE; + + for (int i = -3; i <= 3; i++) { + + double shift = TWO_PI * i; + + TestUtils.assertEquals(ZERO, function.invoke(shift + ZERO)); + TestUtils.assertEquals(ONE, function.invoke(shift + HALF_PI)); + TestUtils.assertEquals(ZERO, function.invoke(shift + PI)); + TestUtils.assertEquals(NEG, function.invoke(shift + PI + HALF_PI)); + TestUtils.assertEquals(ZERO, function.invoke(shift + TWO_PI)); + } + } + +} diff --git a/src/test/java/org/ojalgo/function/series/SeriesFunctionTests.java b/src/test/java/org/ojalgo/function/series/SeriesFunctionTests.java new file mode 100644 index 000000000..8285151f9 --- /dev/null +++ b/src/test/java/org/ojalgo/function/series/SeriesFunctionTests.java @@ -0,0 +1,29 @@ +/* + * Copyright 1997-2023 Optimatika + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +package org.ojalgo.function.series; + + +public abstract class SeriesFunctionTests { + + static final boolean DEBUG = false; + +} diff --git a/src/test/java/org/ojalgo/scalar/ComplexNumberTest.java b/src/test/java/org/ojalgo/scalar/ComplexNumberTest.java index 637309f45..1c377f486 100644 --- a/src/test/java/org/ojalgo/scalar/ComplexNumberTest.java +++ b/src/test/java/org/ojalgo/scalar/ComplexNumberTest.java @@ -25,6 +25,7 @@ import org.junit.jupiter.api.Test; import org.ojalgo.TestUtils; +import org.ojalgo.function.constant.PrimitiveMath; public class ComplexNumberTest extends ScalarTests { @@ -34,7 +35,6 @@ public void testModulus() { TestUtils.assertEquals(ONE, ComplexNumber.valueOf(NEG).getModulus()); TestUtils.assertEquals(NaN, ComplexNumber.NaN.getModulus()); TestUtils.assertEquals(FIVE, ComplexNumber.of(FOUR, THREE).getModulus()); - } @Test @@ -51,4 +51,29 @@ public void testPower() { TestUtils.assertEquals(64L, base.power(6)); } + @Test + public void testUnitRoots() { + + for (int p = 1; p <= 9; p++) { + + ComplexNumber[] unitRoots = ComplexNumber.newUnitRoots(p); + + for (int r = 0; r < unitRoots.length; r++) { + + ComplexNumber root = unitRoots[r]; + + TestUtils.assertEquals(PrimitiveMath.ONE, root.getModulus()); + // The range can be either [0, 2π] or [-π,π] + TestUtils.assertInRange(-PrimitiveMath.PI, PrimitiveMath.TWO_PI, root.getArgument()); + + ComplexNumber unit = root.power(p); + + TestUtils.assertEquals(PrimitiveMath.ONE, unit.getModulus()); + TestUtils.assertEquals(PrimitiveMath.ZERO, unit.getArgument()); + TestUtils.assertEquals(PrimitiveMath.ONE, unit.getReal()); + TestUtils.assertEquals(PrimitiveMath.ZERO, unit.getArgument()); + } + } + } + }