There are 5 methods in RationalNumber that take primitive double as input:

RationalNumber#valueOf(double) RationalNumber#add(double) RationalNumber#subtract(double) RationalNumber#multiply(double) RationalNumber#divide(double)

Their implementations all create various Objects before arriving at the long numerator and denominator of the "new" RationalNumber instance. It should be possible to go directly from the input double to corresponding long numerator and denominator without creating any intermediate Object.

double arg => long numerator, denominator

Deep neural network is the best tool to have image classification. Right now I have a repository named Deeplearning2C and it's generating a neural network in C and MATLAB. The problem is that Deeplearning2C is using Deeplearning4J. Heavy library that updates often. Deeplearning2C is made for Android as well, but the problem is the application is 500 Mb or more...

I'm planning to replace the library with your simple neural network library inside ojAlgo. So I have a question.

Is it possible to get all the matrices from the neural network?

ojAlgo usually has very competitive performance, especially considering it is pure Java. When comparing it to native LPSolve, it often performs almost to the same level.

However, I am occasionally observing strange performance degradation occasionally when the same model gets a few extra variables and constraints added to them (the additions wouldn't be constraining the original model). LPsolve would solve the model in approximately the same time, while ojAlgo slows down significantly. I don't think it is a Java vs native issue, although I could be wrong.

The attached 2 LPs demonstrate this (they are in the simple LPSolve text format). They are machine generated and might seem a bit redundant in some places, but they show the problem in isolation well.

The first file t_206.txt is solved by LPSolve (using the LPSolve IDE) in 281ms. The same LP is solved by ojAlgo in 323ms. Pretty close. It generates 261 variables and 474 constraints in ojAlgo's ExpressionBasedModel.

The second file t_235.txt is solved by LPSolve in 346ms. The same LP however is solved by ojAlgo in 2358ms. It generates 292 variables and 529 constraints on ojAlgo's ExpressionBasedModel.

I am parsing exactly the same files to both. I am using ojAlgo 47.0.0.

In case of ojAlgo what I am doing programmatically is that if the constraint has only one variable, the constant is divided by the coefficient, and the result is added as a bound to the variable itself, unless it already has a tighter bound (in case another constraint on the variable was defined). This way I don't add extra redundant constraints.

switch (inequality) {
           case LESS_THAN_EQUAL:
               //check if the existent limit is already tighter
               if ((variable.getUpperLimit() == null)
                       || (Double.compare(variable.getUpperLimit().doubleValue(), limit) > 0)) {
                   variable.upper(limit);
               }
               break;

           case GREATER_THAN_EQUAL:
               //check if the existent limit is already tighter
               if ((variable.getLowerLimit() == null)
                       || (Double.compare(variable.getLowerLimit().doubleValue(), limit) < 0)) {
                   variable.lower(limit);
               }
               break;

           case EQUAL:
               variable.level(limit);
               break;
}

If you diff the 2 files, you will see that t_235 just has some extra variables added, and constraints on them to get their value from other variables.

What could be the reason for this sudden performance drop? As more constraints get added I notice similar drops ... going up to 34 seconds. Is there anything I can do to understand a bit more what is going on?

t_206.txt t_235.txt

There are still quite a few corners cut; but all the tests are now passing, and I will be filing separate issues for the rest of the problems.

Above all, I need the initial feedback :)

Hi @apete,

I'm at the point where I'm ready to start deploying okAlgo 0.0.1 to Maven Central so myself and others can test it as a dependency.

I wanted to pass it by you first and offer to have it put under the optimatika domain rather than mine org.nield. I think the core features at the moment are DSL's and PuLP-like expression sytnax.

https://github.com/thomasnield/okAlgo

Thoughts?

This feature request is not related to a problem. I want to achieve something along the lines of the following example:

import org.ojalgo.netio.BasicLogger;
import org.ojalgo.optimisation.Expression;
import org.ojalgo.optimisation.ExpressionsBasedModel;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.Variable;

public class Example {
  public static Optimisation.State f() {
    ExpressionsBasedModel model = new ExpressionsBasedModel();

    Variable x = model.addVariable("x").weight(1.0);

    Expression c0 = model.addExpression("c0").upper(1);
    c0.set(x, 1);

    Expression c1 = model.addExpression("c1").lower(2);
    c1.set(x, 1);

    ExpressionsBasedModel model_simplified = model.simplify();

    if (model_simplified.isInfeasible()) {  // error: cannot be accessed from outside package
      return Optimisation.State.INFEASIBLE;
    } else {
      Optimisation.Result result = model_simplified.maximise();
      return result.getState();
    }
  }

  public static void main(final String[] args) {
    BasicLogger.debug(f());
  }
}

I have observed that model.simplify() returns an equivalent model when the model is either feasible or cannot be detected to be infeasible by the presolver. In this case calling model_simplified.maximise() appears to be reasonable. However, if the presolver detects that the problem is infeasible, conflicting constraints are removed and the returned problem is no longer equivalent to the original one. In this case the result of model_simplified.maximise() is not useful for me.

So I assume that the presolver is reliable and I want to return INFEASIBLE. There are two problems:

• The function isInfeasible() cannot be accessed as it is not public
• I am not familiar with ojAlgo and I am not sure that what I want to do is reasonable to begin with.

Alternatives: Something we do currently is:

1. Call model_simplified.maximise()
2. If the result is infeasible exit
3. If the problem is not solved within a given time we assume that the problem admits a solution and call model_simplified.maximise(). The first solve, however, seems redundant and all this feels quite hacky.

The reason I want to simplify() is that using the presolver really makes a difference on the problem instances we are dealing with.

Probably there is a standard/simpler/better/more reliable way to achieve this?

Thanks, Dimitar

I need the Lagrangian multipliers for a Quadratic Programming solution but they are not returned by result.getMultipliers() in this example using version 48.2.0 of ojAlgo:

package com.mycompany.testojalgo;

import java.util.Optional;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.Primitive64Store;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.convex.ConvexSolver;
import org.ojalgo.structure.Access1D;

public class ojAlgoQP {

   public static void main(String[] args) {
      testOjAlgoQuadraticProgramming2();
   }

   public static void testOjAlgoQuadraticProgramming2() {
//  QP Example 16.2 p453 in 'Numerical Optimization', 2ed, (2006), Jorge Nocedal and Stephen J. Wright.
//  minimize function F(x1,x2,x3) = 3*x1*x1 + 2*x1*x2 + x1*x3 + 2.5*x2*x2 + 2*x2*x3 + 2*x3*x3 - 8*x1 - 3*x2 - 3*x3
//  x = [x1, x2, x3]'
//  F(x) = 1/2*x'*H*x + x'*g
//  constraints x1 + x3 = 3, x2 + x3 = 0
//  A*x = b

//objectiveGradient
      Primitive64Store g = Primitive64Store.FACTORY.rows(new double[][]{
         {-8}, {-3}, {-3}
      });
//objectiveHessian
      Primitive64Store H = Primitive64Store.FACTORY.rows(new double[][]{
         {6, 2, 1},
         {2, 5, 2},
         {1, 2, 4}
      });
// constraint equations
      Primitive64Store A = Primitive64Store.FACTORY.rows(new double[][]{
         {1, 0, 1},
         {0, 1, 1}
      });
// required constraint values
      Primitive64Store b = Primitive64Store.FACTORY.rows(new double[][]{
         {3}, {0}
      });
      ConvexSolver.Builder builder = ConvexSolver.getBuilder();
      builder.equalities(A, b);
      builder.objective(H, g.negate());
      ConvexSolver solver = builder.build();
      Optimisation.Result result = solver.solve();

// How do I get the multipliers?  multipliers = Optional.empty
      Optional<Access1D<?>> multipliers = result.getMultipliers();
// value1 = -3.5
      double value1 = result.getValue();

// Verify result:
// x= [2.0, -0.9999999999999996, 0.9999999999999997]';
// value = -3.5
// residual =[-4.440892098500626E-16, 1.1102230246251565E-16]'
      Primitive64Store x = Primitive64Store.FACTORY.column(result.toRawCopy1D());
      MatrixStore<Double> value = x.transpose().multiply(H.multiply(0.5)).multiply(x).add(x.transpose().multiply(g));
      MatrixStore<Double> residual = A.multiply(x).subtract(b);

   }
}

Here is a minimal example of the SVD code

public class Svd {
	
	static Logger logger = LoggerFactory.getLogger(Svd.class);
	
	static public void svd(double A[][], double U[][], double S[], int rows, int columns) {
		
		// Create data
		Primitive64Matrix data = Primitive64Matrix.FACTORY.rows(A);
		
		// Print data
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < columns; j++) {
				System.out.print("\t" + data.get(i, j));
			}
			System.out.println("");
		}
		
		// Perform SVD
		SingularValue<Double> svd = SingularValue.PRIMITIVE.make();
		svd.decompose(data);
		MatrixStore<Double> svdU = svd.getU();
		MatrixStore<Double> svdS = svd.getD();
		MatrixStore<Double> svdV = svd.getV();
		logger.info("Done with Singular Value Decomposition");
		
		// Copy over to U, S
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < columns; j++) {
				U[i][j] = svdU.get(i, j);
				System.out.print("\t" + U[i][j]);
			}
			System.out.println("");
		}
		for(int i = 0; i < columns; i++) {
			S[i] = svdS.get(i, i);
			System.out.println(S[i]);
		}	
	}
}

The output looks like this:

	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
	-3.0	-2.0	-1.0	0.0	1.0	2.0	3.0
[main] INFO se.danielmartensson.fisherfaces.matlab.ojalgo.Svd - Done with Singular Value Decomposition
	-0.316227766016838	0.9486832980505135	-1.1102230246251565E-16	-1.3877787807814457E-17	-4.1633363423443376E-17	-1.3877787807814457E-17	2.7755575615628914E-17
	-0.31622776601683794	-0.10540925533894596	0.9428090415820636	-3.334181172154923E-18	-1.0002543516464769E-17	-1.721196897996938E-17	-2.1087213271319068E-17
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775786	0.9354143466934853	1.7040869798512316E-16	-6.803628124108537E-18	-2.7053155959738237E-19
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621228	0.9258200997725512	2.1524097680092275E-16	2.7485044056031528E-17
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209183	0.9128709291752767	2.0163162690960059E-16
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209188	-0.18257418583505536	0.8944271909999159
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209188	-0.18257418583505536	-0.223606797749979
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209188	-0.18257418583505536	-0.223606797749979
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209188	-0.18257418583505536	-0.223606797749979
	-0.31622776601683794	-0.10540925533894602	-0.11785113019775792	-0.13363062095621217	-0.15430334996209188	-0.18257418583505536	-0.223606797749979
16.733200530681515
0.0
0.0
0.0
0.0
0.0
0.0

But according to Octave, the output looks like this:

[U, D, V] = svd(X, 'econ');
X
U
D
------------------------------------
X =

  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3

U =

  -0.31623   0.94868  -0.00000   0.00000  -0.00000   0.00000   0.00000
  -0.31623  -0.10541  -0.33333  -0.33333  -0.33333  -0.33333  -0.33333
  -0.31623  -0.10541   0.91667  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333   0.91667  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333   0.91667  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333   0.91667  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333   0.91667
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333

D =

Diagonal Matrix

   16.73320          0          0          0          0          0          0
          0    0.00000          0          0          0          0          0
          0          0   -0.00000          0          0          0          0
          0          0          0    0.00000          0          0          0
          0          0          0          0    0.00000          0          0
          0          0          0          0          0    0.00000          0
          0          0          0          0          0          0    0.00000

As you can see. the U matrix is not correct compared with Octave. Why? Is it because I'm using PRIMITIVE? Don't know what it means.

I have also a question about eigendecomposition. Here is my example code:

public class Eig {
	
	static Logger logger = LoggerFactory.getLogger(Eig.class);
	
	// A*D = D*B*V - Find D and V
	static public void eig(double A[][], double B[][], double D[], double V[][], int rows) {
		
		// Create eigA and eigB
		Primitive64Matrix eigA = Primitive64Matrix.FACTORY.rows(A);
		Primitive64Matrix eigB = Primitive64Matrix.FACTORY.rows(B);
		
		// Find inverse of eigB - Or pseudo inverse if B is singular
		eigB = eigB.invert();
		
		// Multiply eigB*eigA
		Primitive64Matrix data = eigB.multiply(eigA);
		
		// Perform eigendecomposition
		Eigenvalue<Double> eig = Eigenvalue.PRIMITIVE.make(data,false);
		boolean success = eig.decompose(data);
		if(success == false)
			logger.error("Could not perform eigenvalue decomposition!");
		MatrixStore<Double> Deig = eig.getD();
		MatrixStore<Double> Veig = eig.getV();
		
		// Copy over to D, V
		for(int i = 0; i < rows; i++) {
			D[i] = Deig.get(i, i);
			for(int j = 0; j < rows; j++) {
				V[i][j] = Veig.get(i, j);
			}
		}	
	}
}

I call that code with this data and get the output:

// Create random data
		int rows = 5;
		double[][] A = new double[rows][rows];
		double[][] B = {{0.7868535918842593,	0.25206982118540255,	0.8502514736661777,	0.957057138020801,	0.057816696367313236},
				{0.13879101503702707,	0.7128202999301166,	0.1301621659119453,	0.31855375514373074,	0.5217887100001347},
				{0.32528138064626044,	0.375732472964022,	0.8722279148638662,	0.2623839802297663,	0.6554277853039092},
				{0.29653713361927037,	0.5156537331156116,	0.0769602762384829,	0.5727731028190134,	0.22446173846280937},
				{0.48182976689909873,	0.19256889992091686,	0.6541198727551285,	0.6268943810659438,	0.42040774186487684}};
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < rows; j++) {
				A[i][j] = i+j;
			}
		}
		System.out.println("Print A");
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < rows; j++) {
				System.out.print("\t" + A[i][j]);
			}
			System.out.println("");
		}
		System.out.println("Print B");
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < rows; j++) {
				System.out.print("\t" + B[i][j]);
			}
			System.out.println("");
		}
		
		// Perform eig
		double[][] V = new double[rows][rows];
		double[] D = new double[rows];
		Eig.eig(A, B, D, V, rows);
		
		System.out.println("Print eigenvalues");
		for(int i = 0; i < rows; i++) {
			System.out.println(D[i]);
		}	
		System.out.println("Print eigenvectors");
		for(int i = 0; i < rows; i++) {
			for(int j = 0; j < rows; j++) {
				System.out.print("\t" + V[i][j]);
			}
			System.out.println("");
		}

Print A
	0.0	1.0	2.0	3.0	4.0
	1.0	2.0	3.0	4.0	5.0
	2.0	3.0	4.0	5.0	6.0
	3.0	4.0	5.0	6.0	7.0
	4.0	5.0	6.0	7.0	8.0
Print B
	0.7868535918842593	0.25206982118540255	0.8502514736661777	0.957057138020801	0.057816696367313236
	0.13879101503702707	0.7128202999301166	0.1301621659119453	0.31855375514373074	0.5217887100001347
	0.32528138064626044	0.375732472964022	0.8722279148638662	0.2623839802297663	0.6554277853039092
	0.29653713361927037	0.5156537331156116	0.0769602762384829	0.5727731028190134	0.22446173846280937
	0.48182976689909873	0.19256889992091686	0.6541198727551285	0.6268943810659438	0.42040774186487684
Print eigenvalues
-72.24303468512798
19.68243960526661
-1.263846900438926E-13
2.1790180928485528E-14
-2.598102440286286E-14
Print eigenvectors
	-0.843534098386463	6.888766304572221	-0.4103808547114451	-0.15064677005171725	0.5300697841954638
	-0.021650343361678242	-0.7125750738079105	0.8421978260532967	0.5003740977706913	0.5722830367666923
	0.2626516887639604	-2.49185111743384	-0.18408075846905578	-0.7342555017208549	-1.4339375479081005
	0.45288497844928716	-4.075205813120237	-0.5169085423758242	0.5699757903365146	-0.96925315126561
	-0.11785912852251558	0.9110743367448255	0.26917232950303754	-0.18544761633463291	1.3008378782115613

But GNU Octave shows

V =

  -0.8435341  -0.8140592   0.5179533   0.4271032   0.0036528
  -0.0216503   0.0842064  -0.7341037  -0.0259704  -0.1009847
   0.2626517   0.2944670  -0.1408616  -0.7128079   0.5296932
   0.4528850   0.4815752   0.4122211  -0.2048857  -0.7710436
  -0.1178591  -0.1076635  -0.0552091   0.5165609   0.3386822

D =

Diagonal Matrix

  -7.2243e+01            0            0            0            0
            0   1.9682e+01            0            0            0
            0            0   1.2203e-12            0            0
            0            0            0   1.0295e-14            0
            0            0            0            0  -6.6871e-14

From the same data

A = [0.0 1.0	2.0	3.0	4.0
    1.0	2.0	3.0	4.0	5.0
    2.0	3.0	4.0	5.0	6.0
    3.0	4.0	5.0	6.0	7.0
    4.0	5.0	6.0	7.0	8.0];
  
B = [0.7868535918842593	0.25206982118540255	0.8502514736661777	0.957057138020801	0.057816696367313236
	0.13879101503702707	0.7128202999301166	0.1301621659119453	0.31855375514373074	0.5217887100001347
	0.32528138064626044	0.375732472964022	0.8722279148638662	0.2623839802297663	0.6554277853039092
	0.29653713361927037	0.5156537331156116	0.0769602762384829	0.5727731028190134	0.22446173846280937
	0.48182976689909873	0.19256889992091686	0.6541198727551285	0.6268943810659438	0.42040774186487684];
  

% A*D = D*B*V - Find V and D
[V, D] = eig(pinv(B)*A)

Notice that I'm using random matrices.

I was editing https://github.com/optimatika/ojAlgo/wiki/Getting-Started page to fix https://github.com/optimatika/ojAlgo/issues/71, then I wanted to update my revision comment and I think I somehow deleted the whole page. I also do not see any history there, hope you can rollback the deletion.

Just a thought. I've been messing around with ojAlgo to back a neutral network. Considering there is already an MIP solver, a neural network API would also be cool and this library seems to have all the pieces needed for it.

Here's what I've been doing so far.

https://github.com/thomasnield/kotlin_simple_neural_network/blob/master/src/main/kotlin/NeuralNetwork.kt

It seems to work, mostly. We have 9 tests failing, most of them due to small discrepancies, and we do need to investigate those - but overall it's surprisingly robust for a half an hour job.