Moved Mathematics class into tools.math package.

src/eva2/gui/GraphPointSet.java

@@ -20,8 +20,8 @@ import java.io.Serializable;
 import java.util.ArrayList;
 import java.util.Arrays;
 
-import eva2.tools.Mathematics;
 import eva2.tools.chart2d.*;
+import eva2.tools.math.Mathematics;
  /*==========================================================================*
  * CLASS DECLARATION
  *==========================================================================*/

src/eva2/gui/TopoPlot.java

@@ -15,9 +15,9 @@ package eva2.gui;
 import java.awt.*;
 
 import eva2.server.go.problems.Interface2DBorderProblem;
-import eva2.tools.Mathematics;
 import eva2.tools.chart2d.*;
 import eva2.tools.diagram.ColorBarCalculator;
+import eva2.tools.math.Mathematics;
 
 
 /*==========================================================================*

src/eva2/server/go/individuals/ESIndividualDoubleData.java

@@ -6,7 +6,7 @@ import eva2.server.go.operators.mutation.InterfaceMutation;
 import eva2.server.go.operators.mutation.MutateESGlobal;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /** This individual uses a real-valued genotype to code for double values.

src/eva2/server/go/individuals/codings/gp/AbstractGPNode.java

@@ -9,9 +9,9 @@ import eva2.gui.BeanInspector;
 import eva2.gui.GenericObjectEditor;
 import eva2.server.go.problems.GPFunctionProblem;
 import eva2.server.go.problems.InterfaceProgramProblem;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
 import eva2.tools.ReflectPackage;
+import eva2.tools.math.Mathematics;
 
 
 /** This gives an abstract node, with default functionality for get and set methods.

src/eva2/server/go/individuals/codings/gp/GPNodeProd.java

@@ -2,7 +2,7 @@ package eva2.server.go.individuals.codings.gp;
 
 
 import eva2.server.go.problems.InterfaceProgramProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /** 
  * A simple product node with a single, possibly vectorial (array), argument.

src/eva2/server/go/individuals/codings/gp/GPNodeSum.java

@@ -2,7 +2,7 @@ package eva2.server.go.individuals.codings.gp;
 
 
 import eva2.server.go.problems.InterfaceProgramProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /** 
  * A simple sum node with a single, possibly vectorial (array), argument.

src/eva2/server/go/operators/constraint/IntervalConstraint.java

@@ -5,7 +5,7 @@ import java.io.Serializable;
 import eva2.gui.GenericObjectEditor;
 import eva2.server.go.individuals.AbstractEAIndividual;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * A constraint for a parameter or a generic function to lie within certain bounds.

src/eva2/server/go/operators/crossover/CrossoverESPCX.java

@@ -10,7 +10,7 @@ import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.F1Problem;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /**

src/eva2/server/go/operators/crossover/CrossoverESSPX.java

@@ -10,7 +10,7 @@ import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.F1Problem;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /**

src/eva2/server/go/operators/crossover/CrossoverESUNDX.java

@@ -10,7 +10,7 @@ import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.F1Problem;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /**

src/eva2/server/go/operators/mutation/MutateESCorrVector.java

@@ -7,8 +7,8 @@ import eva2.server.go.individuals.AbstractEAIndividual;
 import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
 import eva2.tools.SelectedTag;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 

src/eva2/server/go/operators/mutation/MutateESCovarianceMatrixAdaption.java

@@ -4,7 +4,7 @@ import eva2.server.go.individuals.AbstractEAIndividual;
 import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 import eva2.tools.math.Jama.EigenvalueDecomposition;
 import eva2.tools.math.Jama.Matrix;

src/eva2/server/go/operators/mutation/MutateESMainVectorAdaption.java

@@ -4,7 +4,7 @@ import eva2.server.go.individuals.AbstractEAIndividual;
 import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /**

src/eva2/server/go/operators/mutation/MutateESPathLengthAdaption.java

@@ -4,7 +4,7 @@ import eva2.server.go.individuals.AbstractEAIndividual;
 import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /**

src/eva2/server/go/operators/mutation/MutateESRankMuCMA.java

@@ -13,8 +13,8 @@ import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
 import eva2.server.go.strategies.EvolutionStrategies;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 import eva2.tools.math.Jama.EigenvalueDecomposition;
 import eva2.tools.math.Jama.Matrix;

src/eva2/server/go/operators/paramcontrol/AbstractLinearParamAdaption.java

@@ -4,7 +4,7 @@ import java.io.Serializable;
 
 import eva2.gui.BeanInspector;
 import eva2.server.go.populations.Population;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * Simple linear adaption of a String property.

src/eva2/server/go/operators/paramcontrol/LinearParamAdaption.java

@@ -2,7 +2,7 @@ package eva2.server.go.operators.paramcontrol;
 
 import java.io.Serializable;
 
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * Linearly adapt a specific target String parameter.

src/eva2/server/go/operators/paramcontrol/PSOActivityFeedbackControl.java

@@ -4,7 +4,7 @@ import java.io.Serializable;
 
 import eva2.server.go.populations.Population;
 import eva2.server.go.strategies.ParticleSwarmOptimization;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * After the ANTS 08 paper by Yasuda et al., this implements an activity feedback control mechanism.

src/eva2/server/go/operators/paramcontrol/PSOInertnessAdaption.java

@@ -3,7 +3,7 @@ package eva2.server.go.operators.paramcontrol;
 import java.io.Serializable;
 
 import eva2.gui.GenericObjectEditor;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * Adapt PSO inertness linearly by time, from given start to end value.

src/eva2/server/go/operators/postprocess/PostProcess.java

@@ -41,8 +41,8 @@ import eva2.server.go.strategies.NelderMeadSimplex;
 import eva2.server.modules.GOParameters;
 import eva2.server.stat.InterfaceTextListener;
 import eva2.server.stat.StatsParameter;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
+import eva2.tools.math.Mathematics;
 
 
 /**

src/eva2/server/go/populations/Population.java

@@ -20,8 +20,8 @@ import eva2.server.go.operators.distancemetric.InterfaceDistanceMetric;
 import eva2.server.go.operators.distancemetric.PhenotypeMetric;
 import eva2.server.go.operators.selection.probability.AbstractSelProb;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 import eva2.tools.math.Jama.Matrix;
 import eva2.tools.tool.StatisticUtils;

src/eva2/server/go/problems/AbstractOptimizationProblem.java

@@ -28,7 +28,7 @@ import eva2.server.go.operators.terminators.EvaluationTerminator;
 import eva2.server.go.operators.terminators.PhenotypeConvergenceTerminator;
 import eva2.server.go.populations.Population;
 import eva2.server.go.strategies.InterfaceOptimizer;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * Created by IntelliJ IDEA.

src/eva2/server/go/problems/F6Problem.java

@@ -1,7 +1,7 @@
 package eva2.server.go.problems;
 
 import eva2.server.go.individuals.ESIndividualDoubleData;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.Jama.Matrix;
 
 /**

src/eva2/server/go/strategies/ClusterBasedNichingEA.java

@@ -27,10 +27,10 @@ import eva2.server.go.problems.Interface2DBorderProblem;
 import eva2.server.go.problems.InterfaceAdditionalPopulationInformer;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
 import eva2.server.go.problems.TF1Problem;
-import eva2.tools.Mathematics;
 import eva2.tools.chart2d.DPoint;
 import eva2.tools.chart2d.DPointIcon;
 import eva2.tools.chart2d.DPointSet;
+import eva2.tools.math.Mathematics;
 
 /** The infamous clustering based niching EA, still under construction.
  * It should be able to identify and track multiple global/local optima

src/eva2/server/go/strategies/DifferentialEvolution.java

@@ -17,7 +17,7 @@ import eva2.server.go.problems.AbstractOptimizationProblem;
 import eva2.server.go.problems.F1Problem;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /** 

src/eva2/server/go/strategies/DynamicParticleSwarmOptimization.java

@@ -8,8 +8,8 @@ import eva2.server.go.individuals.InterfaceESIndividual;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.AbstractOptimizationProblem;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
 import eva2.tools.SelectedTag;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 
 /** 

src/eva2/server/go/strategies/NelderMeadSimplex.java

@@ -13,7 +13,7 @@ import eva2.server.go.populations.SolutionSet;
 import eva2.server.go.problems.AbstractOptimizationProblem;
 import eva2.server.go.problems.AbstractProblemDouble;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /**
  * Nelder-Mead-Simplex does not guarantee an equal number of evaluations within each optimize call

src/eva2/server/go/strategies/ParticleSwarmOptimization.java

@@ -27,10 +27,10 @@ import eva2.server.go.problems.InterfaceAdditionalPopulationInformer;
 import eva2.server.go.problems.InterfaceOptimizationProblem;
 import eva2.server.go.problems.InterfaceProblemDouble;
 import eva2.tools.EVAERROR;
-import eva2.tools.Mathematics;
 import eva2.tools.SelectedTag;
 import eva2.tools.chart2d.DPoint;
 import eva2.tools.chart2d.DPointSet;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.RNG;
 import eva2.tools.math.Jama.Matrix;
 

src/eva2/server/stat/AbstractStatistics.java

@@ -18,8 +18,8 @@ import eva2.server.go.populations.InterfaceSolutionSet;
 import eva2.server.go.populations.Population;
 import eva2.server.go.problems.InterfaceAdditionalPopulationInformer;
 import eva2.server.go.strategies.InterfaceOptimizer;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
+import eva2.tools.math.Mathematics;
 
 /**
  * An abstract class handling statistics. Most important stuff happens in startOptPerformed, stopOptPerformed

src/eva2/server/stat/StatisticsStandalone.java

@@ -23,7 +23,7 @@ import java.util.List;
 
 import eva2.server.go.PopulationInterface;
 import eva2.server.go.problems.InterfaceAdditionalPopulationInformer;
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 
 /*

src/eva2/tools/chart2d/DGrid.java

@@ -19,7 +19,7 @@ package eva2.tools.chart2d;
 import java.awt.Color ;
 import java.awt.Graphics ;
 
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 /*==========================================================================*
  * CLASS DECLARATION

src/eva2/tools/chart2d/ScaledBorder.java

@@ -4,7 +4,7 @@ import java.awt.* ;
 import javax.swing.BorderFactory;
 import javax.swing.border.* ;
 
-import eva2.tools.Mathematics;
+import eva2.tools.math.Mathematics;
 
 import java.awt.geom.AffineTransform ;
 import java.text.NumberFormat;

src/eva2/tools/des/DESSolver.java

@@ -11,44 +11,53 @@ package eva2.tools.des;
 import java.io.Serializable;
 
 /**
- *
+ * 
  * TODO: comment missing
- *
+ * 
  * @since 2.0
- * @version
- * Copyright (c) ZBiT, University of Tübingen, Germany
- * Compiler: JDK 1.6.0
+ * @version Copyright (c) ZBiT, University of Tübingen, Germany Compiler:
+ *          JDK 1.6.0
  * @date Sep 10, 2007
- *
+ * 
  */
 public interface DESSolver extends Serializable {
+	
+	/**
+	 * 
+	 * @return
+	 */
+	public boolean isUnstable();
 
-  /**
-   * put your documentation comment here
-   *
-   * @param DES
-   * @param initalvalue
-   * @param x
-   * @param h
-   * @param steps
-   * @return
-   */
-  public double[][] solve(DESystem DES, double[] initalvalue, double x,
-      double h, int steps);
+	/**
+	 * put your documentation comment here
+	 * 
+	 * @param DES
+	 * @param initalvalue
+	 * @param x
+	 * @param h
+	 * @param steps
+	 * @return
+	 */
+	public double[][] solve(DESystem DES, double[] initalvalue, double x,
+			double h, int steps);
 
-  public double[][] solveAtTimePoints(DESystem DES, double[] initialvalue,
-      double[] timepoints);
-
-  public double[][] solveAtTimePointsWithInitialConditions(DESystem DES,
-      double[][] initconditions, double[] timepoints);
-
-  public boolean isUnstable();
-
-  /*
-   * public double[][] solveatTimepointsSSystem (double[] param, int order,
-   * double[] initialvalue, double[] timepoints); public boolean
-   * lastDESystemInvalid(); public void plotY(); public double[][]
-   * solveatTimepointsSSystemSeparated (double[] param, int order, double[]
-   * initialvalue, double[] timepoints, GEdata gedata, int tooptimize);
-   */
+	/**
+	 * 
+	 * @param DES
+	 * @param initialvalue
+	 * @param timepoints
+	 * @return
+	 */
+	public double[][] solveAtTimePoints(DESystem DES, double[] initialvalue,
+			double[] timepoints);
+	
+	/**
+	 * 
+	 * @param DES
+	 * @param initconditions
+	 * @param timepoints
+	 * @return
+	 */
+	public double[][] solveAtTimePointsWithInitialConditions(DESystem DES,
+			double[][] initconditions, double[] timepoints);
 }

src/eva2/tools/des/RKSolver.java

@@ -1,419 +1,469 @@
 package eva2.tools.des;
 
-import eva2.tools.Mathematics;
+import java.io.Serializable;
+
+import eva2.tools.math.Mathematics;
 
 /**
  * Title: JAVA-EVA Description: Runge-Kutta Method Copyright: Copyright (c) 2002
  * Company: University of Tübingen, Computer Architecture
- *
+ * 
  * @author Hannes Planatscher
  * @author Andreas Dräger
  * @author Marcel Kronfeld
  * @version 1.0 Status: works, but numerical inaccurate
  */
-public class RKSolver implements DESSolver, java.io.Serializable {
-  double                         stepSize    = 0.01;
-  boolean                        nonnegative = true;
-  boolean                        unstableFlag;
-  transient protected double[][] kVals       = null;
-  transient protected double[]   k0tmp, k1tmp, k2tmp;
-  private static boolean useLinearCalc = true;
+public class RKSolver implements DESSolver, Serializable {
+	/**
+	 * 
+	 */
+	private static final long serialVersionUID = -3383457963743552159L;
+	/**
+	 * 
+	 */
+	private static boolean useLinearCalc = true;
 
-  /**
+	/**
+	 * 
+	 * @param args
+	 */
+	public static void main(String args[]) {
+		new RKSolver(0.01);
+	}
+
+	/**
+	 * 
+	 */
+	transient protected double[] k0tmp, k1tmp, k2tmp;
+	/**
+	 * 
+	 */
+	transient protected double[][] kVals = null;
+	/**
+	 * 
+	 */
+	boolean nonnegative = true;
+	/**
+	 * 
+	 */
+	double stepSize = 0.01;
+	/**
+	 * 
+	 */
+	boolean unstableFlag;
+
+	/**
+	 * 
+	 */
+	public RKSolver() {
+	}
+
+	/**
+	 * A constructor.
+	 * 
+	 * @param withLinearCalc
+	 *            set whether the linear or old calculation method will be used.
+	 */
+	public RKSolver(boolean withLinearCalc) {
+		useLinearCalc = withLinearCalc;
+	}
+
+	/**
+	 * put your documentation comment here
+	 */
+	public RKSolver(double stepSize) {
+		this.stepSize = stepSize;
+	}
+
+	/**
+	 * @return
+	 */
+	public double getStepSize() {
+		return stepSize;
+	}
+
+	/**
+	 * @return
+	 */
+	public boolean isUnstable() {
+		return unstableFlag;
+	}
+
+	/**
+	 * @param DES
+	 * @param h
+	 * @param x
+	 * @param Ytemp
+	 * @return
+	 */
+	public double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) {
+		double[][] K = new double[4][];
+		K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp));
+		K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
+				Ytemp, Mathematics.svMult(0.5, K[0]))));
+		K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
+				Ytemp, Mathematics.svMult(0.5, K[1]))));
+		K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd(
+				Ytemp, K[2])));
+
+		double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0],
+				Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics
+						.vvAdd(Mathematics.svMult(2, K[2]), K[3]))));
+		for (int k = 0; k < change.length; k++) {
+			if (Double.isNaN(change[k])) {
+				unstableFlag = true;
+				change[k] = 0;
+				// return result;
+			}
+		}
+		return change;
+	}
+
+	/**
+	 * Linearized code for speed-up (no allocations).
+	 * 
+	 * @param DES
+	 * @param h
+	 * @param x
+	 * @param Ytemp
+	 * @return
+	 */
+	public void rkTerm2(DESystem DES, double h, double x, double[] Ytemp,
+			double[] res) {
+		if (kVals == null) { // "static" vectors which are allocated only once
+			k0tmp = new double[DES.getDESystemDimension()];
+			k1tmp = new double[DES.getDESystemDimension()];
+			k2tmp = new double[DES.getDESystemDimension()];
+			kVals = new double[4][DES.getDESystemDimension()];
+		}
+
+		// double[][] K = new double[4][];
+		DES.getValue(x, Ytemp, kVals[0]);
+		Mathematics.svMult(h, kVals[0], kVals[0]);
+
+		// K[0] = svMult(h, DES.getValue(x, Ytemp));
+
+		Mathematics.svMult(0.5, kVals[0], k0tmp);
+		Mathematics.vvAdd(Ytemp, k0tmp, k0tmp);
+		DES.getValue(x + h / 2, k0tmp, kVals[1]);
+		Mathematics.svMult(h, kVals[1], kVals[1]);
+
+		// K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
+		// K[0]))));
+
+		Mathematics.svMult(0.5, kVals[1], k1tmp);
+		Mathematics.vvAdd(Ytemp, k1tmp, k1tmp);
+		DES.getValue(x + h / 2, k1tmp, kVals[2]);
+		Mathematics.svMult(h, kVals[2], kVals[2]);
+
+		// K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
+		// K[1]))));
+
+		Mathematics.vvAdd(Ytemp, kVals[2], k2tmp);
+		DES.getValue(x + h, k2tmp, k1tmp);
+		Mathematics.svMult(h, k1tmp, kVals[3]);
+
+		// K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2])));
+
+		Mathematics.svMult(2, kVals[2], k0tmp);
+		Mathematics.vvAdd(k0tmp, kVals[3], k0tmp);
+
+		Mathematics.svMult(2, kVals[1], k1tmp);
+		Mathematics.vvAdd(k1tmp, k0tmp, k2tmp);
+
+		Mathematics.vvAdd(kVals[0], k2tmp, k1tmp);
+		Mathematics.svDiv(6, k1tmp, res);
+
+		// double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]),
+		// vvAdd(svMult(2, K[2]), K[3]))));
+		// for (int i=0; i<res.length; i++) {
+		// double diff = Math.abs(res[i]-change[i]);
+		// if (diff > 0.00000001) System.out.println("!!! ");
+		// }
+
+		// double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]),
+		// vvadd(svmult(2, kVals[2]), kVals[3]))));
+		for (int k = 0; k < res.length; k++) {
+			if (Double.isNaN(res[k])) {
+				unstableFlag = true;
+				res[k] = 0;
+				// return result;
+			}
+		}
+	}
+
+	/**
+	 * @param stepSize
+	 */
+	public void setStepSize(double stepSize) {
+		this.stepSize = stepSize;
+	}
+
+	/**
+	 * @param unstableFlag
+	 */
+	public void setUnstableFlag(boolean unstableFlag) {
+		this.unstableFlag = unstableFlag;
+	}
+
+	/**
+	 * Set whether the linear or old calculation method will be used.
+	 * 
+	 * @param withLinearCalc
+	 */
+	public void setWithLinearCalc(boolean withLinearCalc) {
+		useLinearCalc = withLinearCalc;
+	}
+
+	/**
+	 * put your documentation comment here
+	 * 
+	 * @param DES
+	 * @param initialValues
+	 * @param x
+	 * @param h
+	 * @param steps
+	 * @return
+	 */
+	public double[][] solve(DESystem DES, double[] initialValues, double x,
+			double h, int steps) {
+		double[] timeVector = new double[steps];
+		for (int i = 0; i < steps; i++)
+			timeVector[i] = x + i * h;
+		return solveAtTimePoints(DES, initialValues, timeVector);
+	}
+
+	public double[][] solveAtTimePoints(DESystem DES, double[] initialValues,
+			double[] timePoints) {
+		return solveAtTimePoints(DES, initialValues, timePoints, false);
+	}
+
+	/**
+	 * This method returns a matrix in which the first column includes all time
+	 * points. Every row is composed as time and all values at this time point.
+	 * It uses the same integration method than the regular
+	 * <code>solveatTimepoints</code> method.
+	 * 
+	 * @param DES
+	 * @param initialValues
+	 * @param timePoints
+	 * @return
+	 */
+	public double[][] solveAtTimePointsIncludingTime(DESystem DES,
+			double[] initialValues, double[] timePoints) {
+		return solveAtTimePoints(DES, initialValues, timePoints, true);
+	}
+
+	/**
    *
    */
-  public RKSolver() {
-  }
+	public double[][] solveAtTimePointsWithInitialConditions(DESystem DES,
+			double[][] initConditions, double[] timePoints) {
+		int order = DES.getDESystemDimension();
+		double[][] result = new double[timePoints.length][order];
+		result[0] = initConditions[0];
+		double x = timePoints[0];
+		for (int i = 1; i < timePoints.length; i++) {
+			double h = stepSize;
+			double[] Ytemp = new double[order];
+			int inbetweensteps = (int) Math
+					.floor((timePoints[i] - timePoints[i - 1]) / h);
+			Ytemp = (double[]) initConditions[i - 1].clone();
+			for (int j = 0; j < inbetweensteps; j++) {
+				double change[] = rkTerm(DES, h, x, Ytemp);
+				Ytemp = Mathematics.vvAdd(Ytemp, change);
+				x += h;
+			}
 
-  /**
-   * A constructor.
-   *
-   * @param withLinearCalc set whether the linear or old calculation method will be used.
-   */
-  public RKSolver(boolean withLinearCalc) {
-	  useLinearCalc = withLinearCalc;
-  }
+			h = timePoints[i] - x;
+			double change[] = rkTerm(DES, h, x, Ytemp);
 
-  
-  /**
-   * put your documentation comment here
-   */
-  public RKSolver(double stepSize) {
-    this.stepSize = stepSize;
-  }
-  
-  /**
-   * Set whether the linear or old calculation method will be used.
-   *
-   * @param withLinearCalc
-   */
-  public void setWithLinearCalc(boolean withLinearCalc) {
-	  useLinearCalc = withLinearCalc;
-  }
+			Ytemp = Mathematics.vvAdd(Ytemp, change);
 
-  /**
-   * put your documentation comment here
-   *
-   * @param DES
-   * @param initialValues
-   * @param x
-   * @param h
-   * @param steps
-   * @return
-   */
-  public double[][] solve(DESystem DES, double[] initialValues, double x,
-      double h, int steps) {
-    double[] timeVector = new double[steps];
-    for (int i = 0; i < steps; i++)
-      timeVector[i] = x + i * h;
-    return solveAtTimePoints(DES, initialValues, timeVector);
-  }
+			if (this.nonnegative) {
+				for (int k = 0; k < Ytemp.length; k++) {
+					if (Ytemp[k] < 0) {
+						Ytemp[k] = 0;
+					}
+				}
+			}
+			result[i] = Ytemp;
+			x += h;
+		}
 
-  /**
-   * @param DES
-   * @param initialValues
-   * @param timeBegin
-   * @param timeEnd
-   * @return
-   */
-  public double[][] solveByStepSize(DESystem DES, double[] initialValues,
-      double timeBegin, double timeEnd) {
-    return solveByStepSize(DES, initialValues, timeBegin, timeEnd, false);
-  }
+		return result;
 
-  /**
-   * @param DES
-   * @param initialValues
-   * @param timeBegin
-   * @param timeEnd
-   * @return
-   */
-  public double[][] solveByStepSizeIncludingTime(DESystem DES,
-      double[] initialValues, double timeBegin, double timeEnd) {
-    return solveByStepSize(DES, initialValues, timeBegin, timeEnd, true);
-  }
+	}
 
-  /**
-   * @param DES
-   * @param initialValues
-   * @param timeBegin
-   * @param timeEnd
-   * @return
-   */
-  private double[][] solveByStepSize(DESystem DES, double[] initialValues,
-      double timeBegin, double timeEnd, boolean time) {
-    int numsteps = (int) Math.round(((timeEnd - timeBegin) / stepSize) + 1);
-    unstableFlag = false;
-    // System.out.println(numsteps);
-    int order = DES.getDESystemDimension(), i;
-    double[][] result;
-    if (time) {
-      result = new double[numsteps][order + 1];
-      result[0][0] = timeBegin;
-      for (i = 0; i < order; i++)
-        result[0][i + 1] = initialValues[i];
-    } else {
-      result = new double[numsteps][order];
-      for (i = 0; i < order; i++)
-        result[0][i] = initialValues[i];
-    }
-    double x = timeBegin;
-    for (i = 1; i < numsteps; i++) {
-      double h = stepSize, change[] = null, Ytemp[] = null;
-      if (time) {
-        double tmp[] = new double[result[i - 1].length - 1];
-        System.arraycopy(result[i - 1], 1, tmp, 0, result[i - 1].length - 1);
-        change = rkTerm(DES, h, x, tmp);
-        Ytemp = Mathematics.vvAdd(tmp, change);
-      } else {
-        change = rkTerm(DES, h, x, result[i - 1]);
-        Ytemp = Mathematics.vvAdd(result[i - 1], change);
-      }
-      if (this.nonnegative) {
-        for (int k = 0; k < Ytemp.length; k++) {
-          if (Ytemp[k] < 0) Ytemp[k] = 0;
-        }
-      }
-      x += h;
-      if (time) {
-        System.arraycopy(Ytemp, 0, result[i], 1, Ytemp.length);
-        result[i][0] = x;
-      } else result[i] = Ytemp;
-    }
+	/**
+	 * @param DES
+	 * @param initialValues
+	 * @param timeBegin
+	 * @param timeEnd
+	 * @return
+	 */
+	public double[][] solveByStepSize(DESystem DES, double[] initialValues,
+			double timeBegin, double timeEnd) {
+		return solveByStepSize(DES, initialValues, timeBegin, timeEnd, false);
+	}
 
-    return result;
-  }
+	/**
+	 * @param DES
+	 * @param initialValues
+	 * @param timeBegin
+	 * @param timeEnd
+	 * @return
+	 */
+	public double[][] solveByStepSizeIncludingTime(DESystem DES,
+			double[] initialValues, double timeBegin, double timeEnd) {
+		return solveByStepSize(DES, initialValues, timeBegin, timeEnd, true);
+	}
 
-  public double[][] solveAtTimePoints(DESystem DES, double[] initialValues,
-      double[] timePoints) {
-    return solveAtTimePoints(DES, initialValues, timePoints, false);
-  }
+	/**
+	 * When set to <code>TRUE</code>, <code>includeTimes</code> will make the
+	 * solver to return a matrix with the first column containing the times. By
+	 * default the result of the ODE solver just returns the values for Y.
+	 * 
+	 * @param includeTimes
+	 */
+	private double[][] solveAtTimePoints(DESystem DES, double[] initialValues,
+			double[] timePoints, boolean includeTimes) {
+		// sorted timepoints!!!!!!!!!!!!!!!!!!!!!
+		int order = DES.getDESystemDimension();
+		double result[][], x = timePoints[0];
+		if (includeTimes) {
+			result = new double[timePoints.length][order + 1];
+			result[0][0] = timePoints[0];
+			for (int i = 1; i <= order; i++)
+				result[0][i] = initialValues[i - 1];
+		} else {
+			result = new double[timePoints.length][order];
+			for (int i = 0; i < order; i++)
+				result[0][i] = initialValues[i];
+		}
+		// System.out.println("JavaCalled");
+		unstableFlag = false;
 
-  /**
-   * This method returns a matrix in which the first column includes all time
-   * points. Every row is composed as time and all values at this time point. It
-   * uses the same integration method than the regular
-   * <code>solveatTimepoints</code> method.
-   *
-   * @param DES
-   * @param initialValues
-   * @param timePoints
-   * @return
-   */
-  public double[][] solveAtTimePointsIncludingTime(DESystem DES,
-      double[] initialValues, double[] timePoints) {
-    return solveAtTimePoints(DES, initialValues, timePoints, true);
-  }
+		double h = stepSize;
+		double change[] = new double[order];
+		double[] Ytemp = new double[order];
 
-  /**
-   * When set to <code>TRUE</code>, <code>includeTimes</code> will make the
-   * solver to return a matrix with the first column containing the times. By
-   * default the result of the ODE solver just returns the values for Y.
-   *
-   * @param includeTimes
-   */
-  private double[][] solveAtTimePoints(DESystem DES, double[] initialValues,
-      double[] timePoints, boolean includeTimes) {
-    // sorted timepoints!!!!!!!!!!!!!!!!!!!!!
-    int order = DES.getDESystemDimension();
-    double result[][], x = timePoints[0];
-    if (includeTimes) {
-      result = new double[timePoints.length][order + 1];
-      result[0][0] = timePoints[0];
-      for (int i = 1; i <= order; i++)
-        result[0][i] = initialValues[i - 1];
-    } else {
-      result = new double[timePoints.length][order];
-      for (int i = 0; i < order; i++)
-        result[0][i] = initialValues[i];
-    }
-    // System.out.println("JavaCalled");
-    unstableFlag = false;
+		for (int i = 1; i < timePoints.length; i++) {
+			h = stepSize;
 
-    double h = stepSize;
-    double change[] = new double[order];
-    double[] Ytemp = new double[order];
+			// int inbetweensteps = (int) Math.round((timePoints[i] -
+			// timePoints[i -
+			// 1]) / h + 1);
+			int inbetweensteps = (int) Math
+					.floor((timePoints[i] - timePoints[i - 1]) / h);
 
-    for (int i = 1; i < timePoints.length; i++) {
-      h = stepSize;
+			// System.out.println("inbetweensteps at " + i + ": " +
+			// inbetweensteps);
+			if (includeTimes)
+				System.arraycopy(result[i - 1], 1, Ytemp, 0,
+						result[i - 1].length - 1);
+			else
+				Ytemp = result[i - 1].clone();
 
-//       int inbetweensteps = (int) Math.round((timePoints[i] - timePoints[i -
-//       1]) / h + 1);
-      int inbetweensteps = (int) Math.floor((timePoints[i] - timePoints[i - 1])
-          / h);
+			for (int j = 0; j < inbetweensteps; j++) {
+				if (useLinearCalc)
+					rkTerm2(DES, h, x, Ytemp, change);
+				else
+					change = rkTerm(DES, h, x, Ytemp);
+				// System.out.println("aft change 0 " + change[0]);
 
-      //System.out.println("inbetweensteps at " + i + ": " + inbetweensteps);
-      if (includeTimes)
-        System.arraycopy(result[i - 1], 1, Ytemp, 0, result[i - 1].length - 1);
-      else Ytemp = result[i - 1].clone();
+				Mathematics.vvAdd(Ytemp, change, Ytemp);
 
-      for (int j = 0; j < inbetweensteps; j++) {
-          if (useLinearCalc) rkTerm2(DES, h, x, Ytemp, change);
-          else change = rkTerm(DES, h, x, Ytemp);
-//      	System.out.println("aft change 0 " + change[0]);
+				if (this.nonnegative) {
+					for (int k = 0; k < Ytemp.length; k++) {
+						if (Ytemp[k] < 0)
+							Ytemp[k] = 0;
+					}
+				}
 
-          Mathematics.vvAdd(Ytemp, change, Ytemp);
+				x += h;
+			}
 
-        if (this.nonnegative) {
-          for (int k = 0; k < Ytemp.length; k++) {
-            if (Ytemp[k] < 0) Ytemp[k] = 0;
-          }
-        }
+			h = timePoints[i] - x;
 
-        x += h;
-      }
+			if (useLinearCalc)
+				rkTerm2(DES, h, x, Ytemp, change);
+			else
+				change = rkTerm(DES, h, x, Ytemp);
 
-      h = timePoints[i] - x;
+			Mathematics.vvAdd(Ytemp, change, Ytemp);
 
-      if (useLinearCalc) rkTerm2(DES, h, x, Ytemp, change);
-      else change = rkTerm(DES, h, x, Ytemp);
+			if (this.nonnegative) {
+				for (int k = 0; k < Ytemp.length; k++) {
+					if (Ytemp[k] < 0)
+						Ytemp[k] = 0;
+				}
+			}
 
-      Mathematics.vvAdd(Ytemp, change, Ytemp);
+			if (includeTimes) {
+				result[i][0] = timePoints[i];
+				System.arraycopy(Ytemp, 0, result[i], 1, Ytemp.length);
+			} else
+				result[i] = Ytemp;
+			x += h;
 
-      if (this.nonnegative) {
-        for (int k = 0; k < Ytemp.length; k++) {
-          if (Ytemp[k] < 0) Ytemp[k] = 0;
-        }
-      }
+		}
 
-      if (includeTimes) {
-        result[i][0] = timePoints[i];
-        System.arraycopy(Ytemp, 0, result[i], 1, Ytemp.length);
-      } else result[i] = Ytemp;
-      x += h;
+		return result;
+	}
 
-    }
+	/**
+	 * @param DES
+	 * @param initialValues
+	 * @param timeBegin
+	 * @param timeEnd
+	 * @return
+	 */
+	private double[][] solveByStepSize(DESystem DES, double[] initialValues,
+			double timeBegin, double timeEnd, boolean time) {
+		int numsteps = (int) Math.round(((timeEnd - timeBegin) / stepSize) + 1);
+		unstableFlag = false;
+		// System.out.println(numsteps);
+		int order = DES.getDESystemDimension(), i;
+		double[][] result;
+		if (time) {
+			result = new double[numsteps][order + 1];
+			result[0][0] = timeBegin;
+			for (i = 0; i < order; i++)
+				result[0][i + 1] = initialValues[i];
+		} else {
+			result = new double[numsteps][order];
+			for (i = 0; i < order; i++)
+				result[0][i] = initialValues[i];
+		}
+		double x = timeBegin;
+		for (i = 1; i < numsteps; i++) {
+			double h = stepSize, change[] = null, Ytemp[] = null;
+			if (time) {
+				double tmp[] = new double[result[i - 1].length - 1];
+				System.arraycopy(result[i - 1], 1, tmp, 0,
+						result[i - 1].length - 1);
+				change = rkTerm(DES, h, x, tmp);
+				Ytemp = Mathematics.vvAdd(tmp, change);
+			} else {
+				change = rkTerm(DES, h, x, result[i - 1]);
+				Ytemp = Mathematics.vvAdd(result[i - 1], change);
+			}
+			if (this.nonnegative) {
+				for (int k = 0; k < Ytemp.length; k++) {
+					if (Ytemp[k] < 0)
+						Ytemp[k] = 0;
+				}
+			}
+			x += h;
+			if (time) {
+				System.arraycopy(Ytemp, 0, result[i], 1, Ytemp.length);
+				result[i][0] = x;
+			} else
+				result[i] = Ytemp;
+		}
 
-    return result;
-  }
-
-  /**
-   *
-   */
-  public double[][] solveAtTimePointsWithInitialConditions(DESystem DES,
-      double[][] initConditions, double[] timePoints) {
-    int order = DES.getDESystemDimension();
-    double[][] result = new double[timePoints.length][order];
-    result[0] = initConditions[0];
-    double x = timePoints[0];
-    for (int i = 1; i < timePoints.length; i++) {
-      double h = stepSize;
-      double[] Ytemp = new double[order];
-      int inbetweensteps = (int) Math.floor((timePoints[i] - timePoints[i - 1])
-          / h);
-      Ytemp = (double[]) initConditions[i - 1].clone();
-      for (int j = 0; j < inbetweensteps; j++) {
-        double change[] = rkTerm(DES, h, x, Ytemp);
-        Ytemp = Mathematics.vvAdd(Ytemp, change);
-        x += h;
-      }
-
-      h = timePoints[i] - x;
-      double change[] = rkTerm(DES, h, x, Ytemp);
-
-      Ytemp = Mathematics.vvAdd(Ytemp, change);
-
-      if (this.nonnegative) {
-        for (int k = 0; k < Ytemp.length; k++) {
-          if (Ytemp[k] < 0) {
-            Ytemp[k] = 0;
-          }
-        }
-      }
-      result[i] = Ytemp;
-      x += h;
-    }
-
-    return result;
-
-  }
-
-  /**
-   * @param DES
-   * @param h
-   * @param x
-   * @param Ytemp
-   * @return
-   */
-  public double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) {
-    double[][] K = new double[4][];
-    K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp));
-    K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(Ytemp, Mathematics.svMult(0.5, K[0]))));
-    K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(Ytemp, Mathematics.svMult(0.5, K[1]))));
-    K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd(Ytemp, K[2])));
-    
-    double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0], Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics.vvAdd(Mathematics.svMult(
-        2, K[2]), K[3]))));
-    for (int k = 0; k < change.length; k++) {
-      if (Double.isNaN(change[k])) {
-        unstableFlag = true;
-        change[k] = 0;
-        // return result;
-      }
-    }
-    return change;
-  }
-
-  /**
-   * Linearized code for speed-up (no allocations).
-   *
-   * @param DES
-   * @param h
-   * @param x
-   * @param Ytemp
-   * @return
-   */
-  public void rkTerm2(DESystem DES, double h, double x, double[] Ytemp,
-      double[] res) {
-    if (kVals == null) { // "static" vectors which are allocated only once
-      k0tmp = new double[DES.getDESystemDimension()];
-      k1tmp = new double[DES.getDESystemDimension()];
-      k2tmp = new double[DES.getDESystemDimension()];
-      kVals = new double[4][DES.getDESystemDimension()];
-    }
-
-//    double[][] K = new double[4][];
-    DES.getValue(x, Ytemp, kVals[0]);
-    Mathematics.svMult(h, kVals[0], kVals[0]);
-    
-//    K[0] = svMult(h, DES.getValue(x, Ytemp));
-
-    Mathematics.svMult(0.5, kVals[0], k0tmp);
-    Mathematics.vvAdd(Ytemp, k0tmp, k0tmp);
-    DES.getValue(x + h / 2, k0tmp, kVals[1]);
-    Mathematics.svMult(h, kVals[1], kVals[1]);
-    
-//    K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, K[0]))));
-
-    Mathematics.svMult(0.5, kVals[1], k1tmp);
-    Mathematics.vvAdd(Ytemp, k1tmp, k1tmp);
-    DES.getValue(x + h / 2, k1tmp, kVals[2]);
-    Mathematics.svMult(h, kVals[2], kVals[2]);
-    
-//    K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5, K[1]))));
-
-    Mathematics.vvAdd(Ytemp, kVals[2], k2tmp);
-    DES.getValue(x + h, k2tmp, k1tmp);
-    Mathematics.svMult(h, k1tmp, kVals[3]);
-    
-//    K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2])));
-    
-    Mathematics.svMult(2, kVals[2], k0tmp);
-    Mathematics.vvAdd(k0tmp, kVals[3], k0tmp);
-
-    Mathematics.svMult(2, kVals[1], k1tmp);
-    Mathematics.vvAdd(k1tmp, k0tmp, k2tmp);
-
-    Mathematics.vvAdd(kVals[0], k2tmp, k1tmp);
-    Mathematics.svDiv(6, k1tmp, res);
-
-//    double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]), vvAdd(svMult(2, K[2]), K[3]))));
-//    for (int i=0; i<res.length; i++) {
-//    	double diff = Math.abs(res[i]-change[i]);
-//    	if (diff > 0.00000001) System.out.println("!!! ");
-//    }
-
-    // double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]),
-    // vvadd(svmult(2, kVals[2]), kVals[3]))));
-    for (int k = 0; k < res.length; k++) {
-      if (Double.isNaN(res[k])) {
-        unstableFlag = true;
-        res[k] = 0;
-        // return result;
-      }
-    }
-  }
-
-  public static void main(String args[]) {
-    new RKSolver(0.01);
-  }
-
-  /**
-   * @return
-   */
-  public boolean isUnstable() {
-    return unstableFlag;
-  }
-
-  /**
-   * @param unstableFlag
-   */
-  public void setUnstableFlag(boolean unstableFlag) {
-    this.unstableFlag = unstableFlag;
-  }
-
-  /**
-   * @return
-   */
-  public double getStepSize() {
-    return stepSize;
-  }
-
-  /**
-   * @param stepSize
-   */
-  public void setStepSize(double stepSize) {
-    this.stepSize = stepSize;
-  }
+		return result;
+	}
 
 }

src/eva2/tools/math/Jama/Matrix.java

@@ -11,8 +11,8 @@ import java.util.Arrays;
 import java.util.Locale;
 
 import eva2.gui.BeanInspector;
-import eva2.tools.Mathematics;
 import eva2.tools.Pair;
+import eva2.tools.math.Mathematics;
 import eva2.tools.math.Jama.util.Maths;
 
 

src/eva2/tools/math/Mathematics.java

@@ -1,10 +1,10 @@
-package eva2.tools;
+package eva2.tools.math;
 
 import java.util.Arrays;
 import java.util.List;
 
 import eva2.server.go.tools.DoubleArrayComparator;
-import eva2.tools.math.RNG;
+import eva2.tools.EVAERROR;
 import eva2.tools.math.Jama.Matrix;
 import eva2.tools.math.interpolation.BasicDataSet;
 import eva2.tools.math.interpolation.InterpolationException;

src/eva2/tools/math/RNG.java

@@ -3,7 +3,6 @@ package eva2.tools.math;
 import java.util.ArrayList;
 import java.util.Random;
 
-import eva2.tools.Mathematics;
 
 /**
  *