Skip to content

Commit

Permalink
Merge pull request #19 from StaticBeagle/adding-sparse-matrix
Browse files Browse the repository at this point in the history 
Adding Sparse Matrix
  • Loading branch information
@StaticBeagle
StaticBeagle authored Nov 11, 2024
2 parents 393c622 + 31b3a03 commit 525f7c6
Show file tree
Hide file tree
Showing 50 changed files with 4,216 additions and 2,743 deletions.
1 change: 0 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@ Please see [SciPy](https://github.com/StaticBeagle/ETK4J/blob/master/SciPy).
Last but not least, this project includes code that was translated from [numal](https://github.com/JeffBezanson/numal),
and also from [Math.NET](https://www.mathdotnet.com/) please see [Math.NET](https://github.com/StaticBeagle/ETK4J/blob/master/Math.NET.txt).

Some parts of the library do not follow Java naming conventions and the main reason for this is to align with the more familiar Matlab and SciPy syntax.
Arrays are mainly used throughout the library in order to use native doubles but the use of `List`s is encouraged.

There's a set of examples that show how to use some classes contained in the library. The examples
Expand Down
2 changes: 1 addition & 1 deletion pom.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
<modelVersion>4.0.0</modelVersion>
<groupId>com.wildbitsfoundry</groupId>
<artifactId>ETK4J</artifactId>
<version>1.1.0</version>
<version>2.0.0</version>
<name>Engineering Toolkit for Java</name>
<description>Tools and implementation of mathematical methods or backing math for engineering problems and
applications.
Expand Down
34 changes: 17 additions & 17 deletions src/main/java/com/wildbitsfoundry/etk4j/control/LinearTimeInvariantSystem.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@

import com.wildbitsfoundry.etk4j.math.MathETK;
import com.wildbitsfoundry.etk4j.math.complex.Complex;
import com.wildbitsfoundry.etk4j.math.linearalgebra.EigenvalueDecomposition;
import com.wildbitsfoundry.etk4j.math.linearalgebra.Matrix;
import com.wildbitsfoundry.etk4j.math.linearalgebra.EigenvalueDecompositionDense;
import com.wildbitsfoundry.etk4j.math.linearalgebra.MatrixDense;
import com.wildbitsfoundry.etk4j.util.ComplexArrays;
import com.wildbitsfoundry.etk4j.util.DoubleArrays;

Expand Down Expand Up @@ -48,7 +48,7 @@ public enum IntegrationMethod {
* @throws IllegalArgumentException If the length of the initial conditions is different from the number of states.
* @see <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lsim.html">lsim</a>
*/
protected TimeResponse lsim(double[][] input, double[] time, double[] initialConditions,
protected TimeResponse lSim(double[][] input, double[] time, double[] initialConditions,
StateSpace ss, IntegrationMethod integrationMethod) {
double[][] U = DoubleArrays.transpose(input);

Expand All @@ -59,10 +59,10 @@ protected TimeResponse lsim(double[][] input, double[] time, double[] initialCon
throw new IllegalArgumentException("The time array must have at least one element.");
}

Matrix A = ss.getA();
Matrix B = ss.getB();
Matrix C = ss.getC();
Matrix D = ss.getD();
MatrixDense A = ss.getA();
MatrixDense B = ss.getB();
MatrixDense C = ss.getC();
MatrixDense D = ss.getD();

final int noStates = A.getRowCount();
final int noInputs = B.getColumnCount();
Expand Down Expand Up @@ -105,7 +105,7 @@ protected TimeResponse lsim(double[][] input, double[] time, double[] initialCon
M[i] = new double[noStates + noInputs];
}

Matrix expMT = new Matrix(M).transpose().expm();
MatrixDense expMT = new MatrixDense(M).transpose().expm();
double[][] Ad = expMT.subMatrix(0, noStates - 1, 0, noStates - 1).getAs2dArray();
double[][] Bd = expMT.subMatrix(noStates, expMT.getRowCount() - 1, 0, noStates - 1).getAs2dArray();
for (int i = 1; i < noSteps; ++i) {
Expand All @@ -120,15 +120,15 @@ protected TimeResponse lsim(double[][] input, double[] time, double[] initialCon
for (int i = 0; i < noStates; ++i) {
M[i] = DoubleArrays.concatenateAll(A.getRow(i), B.getRow(i), new double[noInputs]);
}
double[][] identity = Matrix.identity(noInputs).getAs2dArray();
double[][] identity = MatrixDense.identity(noInputs).getAs2dArray();
for (int i = noStates, j = 0; i < noStates + noInputs; ++i, ++j) {
M[i] = DoubleArrays.concatenate(new double[noStates + noInputs], identity[j]);
}
for (int i = noStates + noInputs; i < noStates + 2 * noInputs; ++i) {
M[i] = new double[noStates + 2 * noInputs];
}

Matrix expMT = new Matrix(M).transpose().expm();
MatrixDense expMT = new MatrixDense(M).transpose().expm();
double[][] Ad = expMT.subMatrix(0, noStates - 1, 0, noStates - 1).getAs2dArray();
double[][] Bd1 = expMT.subMatrix(noStates + noInputs, expMT.getRowCount() - 1, 0, noStates - 1).getAs2dArray();
double[][] Bd0 = expMT.subMatrix(noStates, noStates + noInputs - 1, 0, noStates - 1).getAs2dArray();
Expand Down Expand Up @@ -159,8 +159,8 @@ private static double[] dot(double[] a, double[][] b) {
if (a.length != b.length) {
throw new IllegalArgumentException("The number of elements in a must match the number of rows in b.");
}
Matrix A = new Matrix(a, 1);
A.multiplyEquals(new Matrix(b));
MatrixDense A = new MatrixDense(a, 1);
A.multiplyEquals(new MatrixDense(b));
return A.getArray();
}

Expand Down Expand Up @@ -241,7 +241,7 @@ public StepResponse step(double[] time, double[] initialConditions) {
* Helper function to support all the step overloads.
*
* @param time The time vector. Can be null then
* {@link #generateDefaultResponseTimes(Matrix, int)} will be used
* {@link #generateDefaultResponseTimes(MatrixDense, int)} will be used
* @param initialConditions Initial conditions of the system.
* @param numberOfPoints Number of points to be used in case the time vector is null.
* @return The StepReponse of the system.
Expand All @@ -251,7 +251,7 @@ protected StepResponse step(double[] time, double[] initialConditions, int numbe
time = time == null ? generateDefaultResponseTimes(ss.getA(), numberOfPoints) : time;
double[][] U = new double[1][];
U[0] = DoubleArrays.ones(time.length);
TimeResponse lSim = lsim(U, time, initialConditions, ss, IntegrationMethod.ZERO_ORDER_HOLD);
TimeResponse lSim = lSim(U, time, initialConditions, ss, IntegrationMethod.ZERO_ORDER_HOLD);
return new StepResponse(lSim.getTime(), lSim.getResponse()[0]);
}

Expand All @@ -267,8 +267,8 @@ protected StepResponse step(double[] time, double[] initialConditions, int numbe
* @param numberOfPoints The number of points to generate.
* @return The default response times.
*/
protected static double[] generateDefaultResponseTimes(Matrix A, int numberOfPoints) {
EigenvalueDecomposition eig = A.eig();
protected static double[] generateDefaultResponseTimes(MatrixDense A, int numberOfPoints) {
EigenvalueDecompositionDense eig = A.eig();
double[] realEig = eig.getRealEigenvalues();
for (int i = 0; i < realEig.length; ++i) {
realEig[i] = Math.abs(realEig[i]);
Expand Down Expand Up @@ -535,7 +535,7 @@ public SISOTimeResponse simulateTimeResponse(double[] input, double[] time,
IntegrationMethod integrationMethod) {
double[][] U = new double[1][time.length];
U[0] = input;
TimeResponse tr = lsim(U, time, initialConditions, this.toStateSpace(), integrationMethod);
TimeResponse tr = lSim(U, time, initialConditions, this.toStateSpace(), integrationMethod);
return new SISOTimeResponse(tr.getTime(), tr.getResponse()[0], tr.getEvolutionOfStateVector());
}

Expand Down
94 changes: 47 additions & 47 deletions src/main/java/com/wildbitsfoundry/etk4j/control/StateSpace.java
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
package com.wildbitsfoundry.etk4j.control;

import com.wildbitsfoundry.etk4j.math.complex.Complex;
import com.wildbitsfoundry.etk4j.math.linearalgebra.ComplexMatrix;
import com.wildbitsfoundry.etk4j.math.linearalgebra.Matrix;
import com.wildbitsfoundry.etk4j.math.linearalgebra.ComplexMatrixDense;
import com.wildbitsfoundry.etk4j.math.linearalgebra.MatrixDense;
import com.wildbitsfoundry.etk4j.math.linearalgebra.NonSquareMatrixException;
import com.wildbitsfoundry.etk4j.util.DoubleArrays;

Expand All @@ -13,19 +13,19 @@
*/
public class StateSpace extends LinearTimeInvariantSystem {

private Matrix A;
private Matrix B;
private Matrix C;
private Matrix D;
private MatrixDense A;
private MatrixDense B;
private MatrixDense C;
private MatrixDense D;

/**
* Constructs a null system where all the matrices are empty matrices.
*/
public StateSpace() {
this.A = Matrix.empty();
this.B = Matrix.empty();
this.C = Matrix.empty();
this.D = Matrix.empty();
this.A = MatrixDense.empty();
this.B = MatrixDense.empty();
this.C = MatrixDense.empty();
this.D = MatrixDense.empty();
}

/**
Expand All @@ -36,10 +36,10 @@ public StateSpace() {
* @param D The feed through matrix in array form.
*/
public StateSpace(double[][] A, double[][] B, double[][] C, double[][] D) {
this.A = A == null ? null : new Matrix(A);
this.B = B == null ? null : new Matrix(B);
this.C = C == null ? null : new Matrix(C);
this.D = D == null ? null : new Matrix(D);
this.A = A == null ? null : new MatrixDense(A);
this.B = B == null ? null : new MatrixDense(B);
this.C = C == null ? null : new MatrixDense(C);
this.D = D == null ? null : new MatrixDense(D);
}

/**
Expand All @@ -49,43 +49,43 @@ public StateSpace(double[][] A, double[][] B, double[][] C, double[][] D) {
* @param C The state to output matrix.
* @param D The feed through matrix.
*/
public StateSpace(Matrix A, Matrix B, Matrix C, Matrix D) {
this.A = A == null ? null : new Matrix(A);
this.B = B == null ? null : new Matrix(B);
this.C = C == null ? null : new Matrix(C);
this.D = D == null ? null : new Matrix(D);
public StateSpace(MatrixDense A, MatrixDense B, MatrixDense C, MatrixDense D) {
this.A = A == null ? null : new MatrixDense(A);
this.B = B == null ? null : new MatrixDense(B);
this.C = C == null ? null : new MatrixDense(C);
this.D = D == null ? null : new MatrixDense(D);
}

/**
* Get the state {@link Matrix}.
* Get the state {@link MatrixDense}.
* @return A copy of the state matrix A.
*/
public Matrix getA() {
return new Matrix(A);
public MatrixDense getA() {
return new MatrixDense(A);
}

/**
* Get input to state {@link Matrix}.
* Get input to state {@link MatrixDense}.
* @return A copy of the input to state matrix B.
*/
public Matrix getB() {
return new Matrix(B);
public MatrixDense getB() {
return new MatrixDense(B);
}

/**
* Get the state to output {@link Matrix}.
* Get the state to output {@link MatrixDense}.
* @return A copy of the state to output matrix C.
*/
public Matrix getC() {
return new Matrix(C);
public MatrixDense getC() {
return new MatrixDense(C);
}

/**
* Get the feed through {@link Matrix}.
* Get the feed through {@link MatrixDense}.
* @return A copy of the feed through matrix D.
*/
public Matrix getD() {
return new Matrix(D);
public MatrixDense getD() {
return new MatrixDense(D);
}

@Override
Expand All @@ -110,10 +110,10 @@ public StateSpace toStateSpace() {
*/
public TransferFunction[] toTransferFunction(int input) {
StateSpace ss = normalize(this);
Matrix A = ss.A;
Matrix B = ss.B;
Matrix C = ss.C;
Matrix D = ss.D;
MatrixDense A = ss.A;
MatrixDense B = ss.B;
MatrixDense C = ss.C;
MatrixDense D = ss.D;
int nin = D.getColumnCount();
int nout = D.getRowCount();
if(input >= nin) {
Expand Down Expand Up @@ -144,7 +144,7 @@ public TransferFunction[] toTransferFunction(int input) {
for(int k = 0; k < nout; ++k) {
double[] Ck = C.getRow(k);
double Dk = D.get(k, 0);
num[k] = A.subtract(new Matrix(dot(B.getArray(), Ck))).poly();
num[k] = A.subtract(new MatrixDense(dot(B.getArray(), Ck))).poly();
DoubleArrays.addElementWiseInPlace(num[k], DoubleArrays.multiplyElementWise(den, Dk - 1));
tfs[k] = new TransferFunction(num[k], den);
}
Expand Down Expand Up @@ -181,16 +181,16 @@ private StateSpace normalize(StateSpace ss) {
}

if(result.A == null) {
result.A = Matrix.empty();
result.A = MatrixDense.empty();
}
if(result.B == null) {
result.B = Matrix.empty();
result.B = MatrixDense.empty();
}
if(result.C == null) {
result.C = Matrix.empty();
result.C = MatrixDense.empty();
}
if(result.D == null) {
result.D = Matrix.empty();
result.D = MatrixDense.empty();
}
result.A = restore(result.A, p, p);
result.B = restore(result.B, p, q);
Expand All @@ -204,7 +204,7 @@ private StateSpace normalize(StateSpace ss) {
Copyright (c) 2001-2002 Enthought, Inc. 2003-2022, SciPy Developers.
All rights reserved. See https://github.com/StaticBeagle/ETK4J/blob/master/SciPy.
*/
private Integer[] shapeOrNull(Matrix m) {
private Integer[] shapeOrNull(MatrixDense m) {
return m != null ? new Integer[] {m.getRowCount(), m.getColumnCount()} : new Integer[2];
}

Expand All @@ -225,9 +225,9 @@ private Integer choiceNotNull(Integer ...params) {
Copyright (c) 2001-2002 Enthought, Inc. 2003-2022, SciPy Developers.
All rights reserved. See https://github.com/StaticBeagle/ETK4J/blob/master/SciPy.
*/
private Matrix restore(Matrix m, int rows, int cols) {
private MatrixDense restore(MatrixDense m, int rows, int cols) {
if(m.isEmpty()) {
return new Matrix(rows, cols);
return new MatrixDense(rows, cols);
} else {
if(m.getRowCount() != rows && m.getColumnCount() != cols) {
throw new RuntimeException("The input arrays have incompatible shapes.");
Expand Down Expand Up @@ -324,7 +324,7 @@ public TimeResponse simulateTimeResponse(double[][] input, double[] time, Integr
*/
public TimeResponse simulateTimeResponse(double[][] input, double[] time, double[] initialConditions,
IntegrationMethod integrationMethod) {
return lsim(input, time, initialConditions, this, integrationMethod);
return lSim(input, time, initialConditions, this, integrationMethod);
}

@Override
Expand All @@ -338,9 +338,9 @@ public Complex evaluateAt(double w) {
* @return The complex frequency response of the system.
*/
public Complex[] evaluateMIMOAt(double w) {
ComplexMatrix inner = Matrix.identity(A.getRowCount()).multiply(Complex.fromImaginary(w)).subtract(A).inv();
ComplexMatrix outer = C.multiply(inner).multiply(B);
return ComplexMatrix.fromRealMatrix(D).add(outer).transpose().getArray();
ComplexMatrixDense inner = MatrixDense.identity(A.getRowCount()).multiply(Complex.fromImaginary(w)).subtract(A).inv();
ComplexMatrixDense outer = C.multiply(inner).multiply(B);
return ComplexMatrixDense.fromRealMatrix(D).add(outer).transpose().getArray();
}

/**
Expand Down
6 changes: 3 additions & 3 deletions src/main/java/com/wildbitsfoundry/etk4j/control/TransferFunction.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import com.wildbitsfoundry.etk4j.constants.ConstantsETK;
import com.wildbitsfoundry.etk4j.math.MathETK;
import com.wildbitsfoundry.etk4j.math.complex.Complex;
import com.wildbitsfoundry.etk4j.math.linearalgebra.Matrix;
import com.wildbitsfoundry.etk4j.math.linearalgebra.MatrixDense;
import com.wildbitsfoundry.etk4j.math.polynomials.Polynomial;
import com.wildbitsfoundry.etk4j.math.polynomials.RationalFunction;
import com.wildbitsfoundry.etk4j.util.ComplexArrays;
Expand Down Expand Up @@ -645,13 +645,13 @@ public StateSpace toStateSpace() {
System.arraycopy(den, 1, fRow, 0, fRow.length);
DoubleArrays.multiplyElementWiseInPlace(fRow, -1.0);

double[][] eye = Matrix.identity(k - 2, k - 1).getAs2dArray();
double[][] eye = MatrixDense.identity(k - 2, k - 1).getAs2dArray();
double[][] A = new double[eye.length + 1][];
A[0] = fRow;
for (int i = 0; i < eye.length; ++i) {
A[i + 1] = Arrays.copyOf(eye[i], eye[0].length);
}
double[][] B = Matrix.identity(k - 1, 1).getAs2dArray();
double[][] B = MatrixDense.identity(k - 1, 1).getAs2dArray();
double[][] C = new double[1][];
double[][] outer = DoubleArrays.outer(new double[]{numPadded[0]}, Arrays.copyOfRange(den, 1, den.length));
C[0] = DoubleArrays.subtractElementWise(Arrays.copyOfRange(numPadded, 1, numPadded.length), outer[0]);
Expand Down
Loading

0 comments on commit 525f7c6

Please sign in to comment.