Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Updated contributed sources - Mantissa 7 upgrade. JIRA: MATH-162 

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@488828 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz 2006-12-19 22:11:07 +00:00
parent 8242bc2644
commit cc73bfb42f
123 changed files with 2212 additions and 4719 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/mantissa/src/org/spaceroots/mantissa/MessagesResources.java
View File

@ -33,7 +33,7 @@ public class MessagesResources
}
public Object[][] getContents() {
return contents;
return (Object[][]) contents.clone();
}
static final Object[][] contents = {

2
src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
View File

@ -32,7 +32,7 @@ public class MessagesResources_fr
}
public Object[][] getContents() {
return contents;
return (Object[][]) contents.clone();
}
static final Object[][] contents = {

79
src/mantissa/src/org/spaceroots/mantissa/algebra/Chebyshev.java
View File

@ -17,18 +17,15 @@
package org.spaceroots.mantissa.algebra;
import java.util.ArrayList;
import java.util.List;
/**
* This class implements Chebyshev polynomials.
* <p>Chebyshev polynomials can be defined by the following recurrence
* relations:
* <pre>
* T0(X) = 1
* T1(X) = X
* Tk+1(X) = 2X Tk(X) - Tk-1(X)
* T<sub>0</sub>(X) = 1
* T<sub>1</sub>(X) = X
* T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
* </pre></p>
* @version $Id: Chebyshev.java 1705 2006-09-17 19:57:39Z luc $
@ -39,67 +36,35 @@ import java.util.List;
public class Chebyshev
extends OrthogonalPolynomial {
/** Generator for the Chebyshev polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// Tk+1(X) = 2X Tk(X) - Tk-1(X)
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2l),
new RationalNumber(1l));
}
};
/** Simple constructor.
* Build a degree 0 Chebyshev polynomial
*/
public Chebyshev() {
super(0, l, maxDegree);
super(0, generator);
}
/** Simple constructor.
* Build a degree d Chebyshev polynomial
* @param d degree of the polynomial
* @param degree degree of the polynomial
*/
public Chebyshev(int d) {
super(d, l, maxDegree);
public Chebyshev(int degree) {
super(degree, generator);
}
/** Initialize the recurrence coefficients.
* The recurrence relation is
* <pre>Tk+1(X) = 2X Tk(X) - Tk-1(X)</pre>
* @param k index of the current step
* @param b2k coefficient to initialize (b2k = a2k / a1k)
* @param b3k coefficient to initialize (b3k = a3k / a1k)
* @param b4k coefficient to initialize (b4k = a4k / a1k)
*/
protected void initRecurrenceCoefficients(int k,
RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k) {
b2k.reset(0l);
b3k.reset(2l);
b4k.reset(1l);
}
/** Set the maximal degree of already computed polynomials.
* @param d maximal degree of already computed polynomials
*/
protected void setMaxDegree(int d) {
maxDegree = d;
}
private static final long serialVersionUID = 8367010179599693222L;
/** List holding the coefficients of the polynomials computed so far. */
private static List l;
/** Maximal degree of the polynomials computed so far. */
private static int maxDegree;
/** Build the first two polynomials. */
static {
l = new ArrayList ();
// T0(X) = 1
l.add(new RationalNumber(1l));
// T1(X) = X
l.add(new RationalNumber(0l));
l.add(new RationalNumber(1l));
maxDegree = 1;
}
private static final long serialVersionUID = -893367988717182601L;
}

79
src/mantissa/src/org/spaceroots/mantissa/algebra/Hermite.java
View File

@ -17,18 +17,15 @@
package org.spaceroots.mantissa.algebra;
import java.util.ArrayList;
import java.util.List;
/**
* This class implements Hermite polynomials.
* <p>Hermite polynomials can be defined by the following recurrence
* relations:
* <pre>
* H0(X) = 1
* H1(X) = 2X
* Hk+1(X) = 2X Hk(X) - 2k Hk-1(X)
* H<sub>0</sub>(X) = 1
* H<sub>1</sub>(X) = 2X
* H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
* </pre></p>
* @version $Id: Hermite.java 1705 2006-09-17 19:57:39Z luc $
@ -38,67 +35,35 @@ import java.util.List;
public class Hermite
extends OrthogonalPolynomial {
/** Generator for the Hermite polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(2l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// Hk+1(X) = 2X Hk(X) - 2k Hk-1(X)
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2l),
new RationalNumber(k * 2l));
}
};
/** Simple constructor.
* Build a degree 0 Hermite polynomial
*/
public Hermite() {
super(0, l, maxDegree);
super(0, generator);
}
/** Simple constructor.
* Build a degree d Hermite polynomial
* @param d degree of the polynomial
* @param degree degree of the polynomial
*/
public Hermite(int d) {
super(d, l, maxDegree);
public Hermite(int degree) {
super(degree, generator);
}
/** Initialize the recurrence coefficients.
* The recurrence relation is
* <pre>Hk+1(X) = 2X Hk(X) - 2k Hk-1(X)</pre>
* @param k index of the current step
* @param b2k coefficient to initialize (b2k = a2k / a1k)
* @param b3k coefficient to initialize (b3k = a3k / a1k)
* @param b4k coefficient to initialize (b4k = a4k / a1k)
*/
protected void initRecurrenceCoefficients(int k,
RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k) {
b2k.reset(0l);
b3k.reset(2l);
b4k.reset(2l * k);
}
/** Set the maximal degree of already computed polynomials.
* @param d maximal degree of already computed polynomials
*/
protected void setMaxDegree(int d) {
maxDegree = d;
}
private static final long serialVersionUID = -4639726453485128770L;
/** Table holding the coefficients of the polynomials computed so far. */
private static List l;
/** Maximal degree of the polynomials computed so far. */
private static int maxDegree;
/** Build the first two polynomials. */
static {
l = new ArrayList ();
// H0(X) = 1
l.add(new RationalNumber(1l));
// H1(X) = 2X
l.add(new RationalNumber(0l));
l.add(new RationalNumber(2l));
maxDegree = 1;
}
private static final long serialVersionUID = 7910082423686662133L;
}

81
src/mantissa/src/org/spaceroots/mantissa/algebra/Laguerre.java
View File

@ -17,18 +17,15 @@
package org.spaceroots.mantissa.algebra;
import java.util.ArrayList;
import java.util.List;
/**
* This class implements Laguerre polynomials.
* <p>Laguerre polynomials can be defined by the following recurrence
* relations:
* <pre>
* L0(X) = 1
* L1(X) = 1 - X
* (k+1) Lk+1(X) = (2k + 1 - X) Lk(X) - k Lk-1(X)
* L<sub>0</sub>(X) = 1
* L<sub>1</sub>(X) = 1 - X
* (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
* </pre></p>
* @version $Id: Laguerre.java 1705 2006-09-17 19:57:39Z luc $
@ -38,68 +35,36 @@ import java.util.List;
public class Laguerre
extends OrthogonalPolynomial {
/** Generator for the Laguerre polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(1l),
new RationalNumber(-1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// (k+1) Lk+1(X) = (2k + 1 - X) Lk(X) - k Lk-1(X)
long kP1 = k + 1;
setRecurrenceCoefficients(new RationalNumber(2 * k + 1, kP1),
new RationalNumber(-1l, kP1),
new RationalNumber(k, kP1));
}
};
/** Simple constructor.
* Build a degree 0 Laguerre polynomial
*/
public Laguerre() {
super(0, l, maxDegree);
super(0, generator);
}
/** Simple constructor.
* Build a degree d Laguerre polynomial
* @param d degree of the polynomial
* @param degree degree of the polynomial
*/
public Laguerre(int d) {
super(d, l, maxDegree);
public Laguerre(int degree) {
super(degree, generator);
}
/** Initialize the recurrence coefficients.
* The recurrence relation is
* <pre>(k+1) Lk+1(X) = (2k + 1 - X) Lk(X) - k Lk-1(X)</pre>
* @param k index of the current step
* @param b2k coefficient to initialize (b2k = a2k / a1k)
* @param b3k coefficient to initialize (b3k = a3k / a1k)
* @param b4k coefficient to initialize (b4k = a4k / a1k)
*/
protected void initRecurrenceCoefficients(int k,
RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k) {
long kP1 = k + 1;
b2k.reset(2 * k + 1, kP1);
b3k.reset(-1l, kP1);
b4k.reset(k, kP1);
}
/** Set the maximal degree of already computed polynomials.
* @param d maximal degree of already computed polynomials
*/
protected void setMaxDegree(int d) {
maxDegree = d;
}
private static final long serialVersionUID = -750526984136835515L;
/** List holding the coefficients of the polynomials computed so far. */
private static List l;
/** Maximal degree of the polynomials computed so far. */
private static int maxDegree;
/** Build the first two polynomials. */
static {
l = new ArrayList ();
// L0(X) = 1
l.add(new RationalNumber(1l));
// L1(X) = 1 - X
l.add(new RationalNumber(1l));
l.add(new RationalNumber(-1l));
maxDegree = 1;
}
private static final long serialVersionUID = 3213856667479179710L;
}

81
src/mantissa/src/org/spaceroots/mantissa/algebra/Legendre.java
View File

@ -17,18 +17,15 @@
package org.spaceroots.mantissa.algebra;
import java.util.ArrayList;
import java.util.List;
/**
* This class implements Legendre polynomials.
* <p>Legendre polynomials can be defined by the following recurrence
* relations:
* <pre>
* P0(X) = 1
* P1(X) = X
* (k+1) Pk+1(X) = (2k+1) X Pk(X) - k Pk-1(X)
* P<sub>0</sub>(X) = 1
* P<sub>1</sub>(X) = X
* (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
* </pre></p>
* @version $Id: Legendre.java 1705 2006-09-17 19:57:39Z luc $
@ -39,68 +36,36 @@ import java.util.List;
public class Legendre
extends OrthogonalPolynomial {
/** Generator for the Legendre polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// (k+1) Pk+1(X) = (2k+1) X Pk(X) - k Pk-1(X)
long kP1 = k + 1;
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2 * k + 1, kP1),
new RationalNumber(k, kP1));
}
};
/** Simple constructor.
* Build a degree 0 Legendre polynomial
*/
public Legendre() {
super(0, l, maxDegree);
super(0, generator);
}
/** Simple constructor.
* Build a degree d Legendre polynomial
* @param d degree of the polynomial
* @param degree degree of the polynomial
*/
public Legendre(int d) {
super(d, l, maxDegree);
public Legendre(int degree) {
super(degree, generator);
}
/** Initialize the recurrence coefficients.
* The recurrence relation is
* <pre>(k+1) Pk+1(X) = (2k+1) X Pk(X) - k Ok-1(X)</pre>
* @param k index of the current step
* @param b2k coefficient to initialize (b2k = a2k / a1k)
* @param b3k coefficient to initialize (b3k = a3k / a1k)
* @param b4k coefficient to initialize (b4k = a4k / a1k)
*/
protected void initRecurrenceCoefficients(int k,
RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k) {
long kP1 = k + 1;
b2k.reset(0l);
b3k.reset(2 * k + 1, kP1);
b4k.reset(k, kP1);
}
/** Set the maximal degree of already computed polynomials.
* @param d maximal degree of already computed polynomials
*/
protected void setMaxDegree(int d) {
maxDegree = d;
}
private static final long serialVersionUID = 428266828791532209L;
/** List holding the coefficients of the polynomials computed so far. */
private static List l;
/** Maximal degree of the polynomials computed so far. */
private static int maxDegree;
/** Build the first two polynomials. */
static {
l = new ArrayList ();
// P0(X) = 1
l.add(new RationalNumber(1l));
// P1(X) = X
l.add(new RationalNumber(0l));
l.add(new RationalNumber(1l));
maxDegree = 1;
}
private static final long serialVersionUID = 4014485393845978429L;
}

120
src/mantissa/src/org/spaceroots/mantissa/algebra/OrthogonalPolynomial.java
View File

@ -17,18 +17,17 @@
package org.spaceroots.mantissa.algebra;
import java.util.List;
/**
* This class is the base class for orthogonal polynomials.
* <p>Orthogonal polynomials can be defined by recurrence relations like:
* <pre>
* O0(X) = some 0 degree polynomial
* O1(X) = some first degree polynomial
* a1k Ok+1(X) = (a2k + a3k X) Ok(X) - a4k Ok-1(X)
* O<sub>0</sub>(X) = some 0 degree polynomial
* O<sub>1</sub>(X) = some first degree polynomial
* a<sub>1,k</sub> O<sub>k+1</sub>(X) = (a<sub>2,k</sub> + a<sub>3,k</sub> X) O<sub>k</sub>(X) - a<sub>4,k</sub> O<sub>k-1</sub>(X)
* </pre>
* where a0k, a1k, a2k and a3k are simple expressions which either are
* where a<sub>1,k</sub>, a<sub>2,k</sub>, a<sub>3,k</sub> and
* a<sub>4,k</sub> are simple expressions which either are
* constants or depend on k.</p>
* @version $Id: OrthogonalPolynomial.java 1705 2006-09-17 19:57:39Z luc $
@ -40,112 +39,11 @@ public abstract class OrthogonalPolynomial
/** Simple constructor.
* Build a degree d orthogonal polynomial
* @param d degree of the polynomial
* @param l list containing all coefficients already computed
* @param maxDegree maximal degree of computed coefficients, this
* coefficient <em>must</em> be greater or equal to 1, i.e. the
* derived class <em>must</em> have initialized the first two
* polynomials of degree 0 and 1 before this constructor can be
* called.
* @param degree degree of the polynomial
* @param generator coefficients generator for the current type of polynomials
*/
protected OrthogonalPolynomial(int d, List l, int maxDegree) {
if (d > maxDegree) {
computeUpToDegree(d, l, maxDegree);
}
// coefficient for polynomial 0 is l [0]
// coefficient for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
// coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
// coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
// coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
// coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
// coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
// ...
int start = d * (d + 1) / 2;
a = new RationalNumber[d+1];
for (int i = 0; i <= d; ++i) {
a[i] = new RationalNumber((RationalNumber) l.get(start + i));
}
unknown = null;
}
/** Initialize the recurrence coefficients.
* The recurrence relation is
* <pre>a1k Ok+1(X) = (a2k + a3k X) Ok(X) - a4k Ok-1(X)</pre>
* @param k index of the current step
* @param b2k coefficient to initialize (b2k = a2k / a1k)
* @param b3k coefficient to initialize (b3k = a3k / a1k)
* @param b4k coefficient to initialize (b4k = a4k / a1k)
*/
protected abstract void initRecurrenceCoefficients(int k,
RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k);
/** Set the maximal degree of already computed polynomials.
* @param d maximal degree of already computed polynomials
*/
protected abstract void setMaxDegree(int d);
/** Compute all the polynomial coefficients up to a given degree.
* @param d maximal degree
* @param l list containing all coefficients already computed
* @param maxDegree maximal degree of computed coefficients
*/
protected void computeUpToDegree(int d, List l, int maxDegree) {
RationalNumber b2k = new RationalNumber();
RationalNumber b3k = new RationalNumber();
RationalNumber b4k = new RationalNumber();
int startK = (maxDegree - 1) * maxDegree / 2;
for (int k = maxDegree; k < d; ++k) {
// start indices of two previous polynomials Ok(X) and Ok-1(X)
int startKm1 = startK;
startK += k;
// a1k Ok+1(X) = (a2k + a3k X) Ok(X) - a4k Ok-1(X)
// we use bik = aik/a1k
initRecurrenceCoefficients(k, b2k, b3k, b4k);
RationalNumber ckPrev = null;
RationalNumber ck = (RationalNumber)l.get(startK);
RationalNumber ckm1 = (RationalNumber)l.get(startKm1);
// degree 0 coefficient
RationalNumber coeff = RationalNumber.multiply(ck, b2k);
coeff.multiplyAndSubtractFromSelf(ckm1, b4k);
l.add(coeff);
// degree 1 to degree k-1 coefficients
for (int i = 1; i < k; ++i) {
ckPrev = ck;
ck = (RationalNumber)l.get(startK + i);
ckm1 = (RationalNumber)l.get(startKm1 + i);
coeff = RationalNumber.multiply(ck, b2k);
coeff.multiplyAndAddToSelf(ckPrev, b3k);
coeff.multiplyAndSubtractFromSelf(ckm1, b4k);
l.add(coeff);
}
// degree k coefficient
ckPrev = ck;
ck = (RationalNumber)l.get(startK + k);
coeff = RationalNumber.multiply(ck, b2k);
coeff.multiplyAndAddToSelf(ckPrev, b3k);
l.add(coeff);
// degree k+1 coefficient
l.add(RationalNumber.multiply(ck, b3k));
}
setMaxDegree(d);
protected OrthogonalPolynomial(int degree, CoefficientsGenerator generator) {
a = generator.getCoefficients(degree);
}
}

856
src/mantissa/src/org/spaceroots/mantissa/algebra/Polynomial.java
View File

File diff suppressed because it is too large Load Diff

241
src/mantissa/src/org/spaceroots/mantissa/algebra/PolynomialFraction.java
View File

@ -17,26 +17,28 @@
package org.spaceroots.mantissa.algebra;
import java.io.Serializable;
import java.math.BigInteger;
/**
* This class implements fractions of polynomials with one unknown and
* rational coefficients.
* <p>Instances of this class are immutable.</p>
* @version $Id: PolynomialFraction.java 1705 2006-09-17 19:57:39Z luc $
* @author L. Maisonobe
*/
public class PolynomialFraction {
public class PolynomialFraction implements Serializable {
/**
* Simple constructor.
* Build a constant null fraction
*/
public PolynomialFraction() {
this(new Polynomial.Rational(new RationalNumber(0l)),
new Polynomial.Rational(new RationalNumber(1l)));
this(new Polynomial.Rational(RationalNumber.ZERO),
new Polynomial.Rational(RationalNumber.ONE));
}
/**
@ -47,8 +49,8 @@ public class PolynomialFraction {
* @exception ArithmeticException if the denominator is null
*/
public PolynomialFraction(long numerator, long denominator) {
this(new Polynomial.Rational(new RationalNumber(numerator)),
new Polynomial.Rational(new RationalNumber(denominator)));
this(new Polynomial.Rational(numerator),
new Polynomial.Rational(denominator));
}
/**
@ -90,21 +92,20 @@ public class PolynomialFraction {
throw new ArithmeticException("null denominator");
}
p = new Polynomial.Rational(numerator);
q = new Polynomial.Rational(denominator);
p = numerator;
q = denominator;
RationalNumber[] a = q.getCoefficients();
if (a[a.length - 1].isNegative()) {
p.negateSelf();
q.negateSelf();
p = (Polynomial.Rational) p.negate();
q = (Polynomial.Rational) q.negate();
}
simplify();
}
/**
* Simple constructor.
/** Simple constructor.
* Build a fraction from a single integer
* @param l value of the fraction
*/
@ -121,8 +122,7 @@ public class PolynomialFraction {
this(i, BigInteger.ONE);
}
/**
* Simple constructor.
/** Simple constructor.
* Build a fraction from a single rational number
* @param r value of the fraction
*/
@ -130,178 +130,95 @@ public class PolynomialFraction {
this(r.getNumerator(), r.getDenominator());
}
/**
* Simple constructor.
* Build a fraction from a single Polynom
/** Simple constructor.
* Build a fraction from a single Polynomial
* @param p value of the fraction
*/
public PolynomialFraction(Polynomial.Rational p) {
this(p, new Polynomial.Rational(new RationalNumber(1l)));
this(p, new Polynomial.Rational(1l));
}
/**
* Copy-constructor.
* @param f fraction to copy
/** Negate the instance.
* @return a new polynomial fraction opposite to the instance
*/
public PolynomialFraction(PolynomialFraction f) {
p = new Polynomial.Rational(f.p);
q = new Polynomial.Rational(f.q);
public PolynomialFraction negate() {
return new PolynomialFraction((Polynomial.Rational) p.negate(), q);
}
/**
* Negate the instance
/** Add a polynomial fraction to the instance.
* @param f polynomial fraction to add.
* @return a new polynomial fraction
*/
public void negateSelf() {
p.negateSelf();
public PolynomialFraction add(PolynomialFraction f) {
return new PolynomialFraction(p.multiply(f.q).add(f.p.multiply(q)),
q.multiply(f.q));
}
/**
* Negate a fraction.
* @param f fraction to negate
* @return a new fraction which is the opposite of f
/** Subtract a fraction from the instance.
* @param f polynomial fraction to subtract.
* @return a new polynomial fraction
*/
public static PolynomialFraction negate(PolynomialFraction f) {
PolynomialFraction copy = new PolynomialFraction(f);
copy.negateSelf();
return copy;
public PolynomialFraction subtract(PolynomialFraction f) {
return new PolynomialFraction(p.multiply(f.q).subtract(f.p.multiply(q)),
q.multiply(f.q));
}
/**
* Add a fraction to the instance.
* @param f fraction to add.
/** Multiply the instance by a polynomial fraction.
* @param f polynomial fraction to multiply by
* @return a new polynomial fraction
*/
public void addToSelf(PolynomialFraction f) {
PolynomialFraction sum = add(this, f);
p = sum.p;
q = sum.q;
public PolynomialFraction multiply(PolynomialFraction f) {
PolynomialFraction product =
new PolynomialFraction(p.multiply(f.p), q.multiply(f.q));
product.simplify();
return product;
}
/** Add two fractions.
* @param f1 first fraction
* @param f2 second fraction
* @return a new fraction which is the sum of f1 and f2
*/
public static PolynomialFraction add(PolynomialFraction f1,
PolynomialFraction f2) {
Polynomial.Rational num =
Polynomial.Rational.add(Polynomial.Rational.multiply(f1.p, f2.q),
Polynomial.Rational.multiply(f2.p, f1.q));
Polynomial.Rational den = Polynomial.Rational.multiply(f1.q, f2.q);
return new PolynomialFraction(num, den);
}
/**
* Subtract a fraction to the instance.
* @param f fraction to subtract.
*/
public void subtractFromSelf(PolynomialFraction f) {
PolynomialFraction diff = subtract(this, f);
p = diff.p;
q = diff.q;
}
/** Subtract two fractions.
* @param f1 first fraction
* @param f2 second fraction
* @return a new fraction which is the difference f1 minus f2
*/
public static PolynomialFraction subtract(PolynomialFraction f1,
PolynomialFraction f2) {
Polynomial.Rational num =
Polynomial.Rational.subtract(Polynomial.Rational.multiply(f1.p, f2.q),
Polynomial.Rational.multiply(f2.p, f1.q));
Polynomial.Rational den = Polynomial.Rational.multiply(f1.q, f2.q);
return new PolynomialFraction(num, den);
}
/** Multiply the instance by a fraction.
* @param f fraction to multiply by
*/
public void multiplySelf(PolynomialFraction f) {
p.multiplySelf(f.p);
q.multiplySelf(f.q);
simplify();
}
/** Multiply two fractions.
* @param f1 first fraction
* @param f2 second fraction
* @return a new fraction which is the product of f1 and f2
*/
public static PolynomialFraction multiply(PolynomialFraction f1,
PolynomialFraction f2) {
PolynomialFraction copy = new PolynomialFraction(f1);
copy.multiplySelf(f2);
return copy;
}
/** Divide the instance by a fraction.
* @param f fraction to divide by
/** Divide the instance by a polynomial fraction.
* @param f polynomial fraction to divide by
* @return a new polynomial fraction
* @exception ArithmeticException if f is null
*/
public void divideSelf(PolynomialFraction f) {
public PolynomialFraction divide(PolynomialFraction f) {
if (f.p.isZero()) {
throw new ArithmeticException("divide by zero");
}
p.multiplySelf(f.q);
q.multiplySelf(f.p);
Polynomial.Rational newP = p.multiply(f.q);
Polynomial.Rational newQ = q.multiply(f.p);
RationalNumber[] a = q.getCoefficients();
RationalNumber[] a = newQ.getCoefficients();
if (a[a.length - 1].isNegative()) {
p.negateSelf();
q.negateSelf();
newP = (Polynomial.Rational) newP.negate();
newQ = (Polynomial.Rational) newQ.negate();
}
simplify();
PolynomialFraction result = new PolynomialFraction(newP, newQ);
result.simplify();
return result;
}
/** Divide two fractions.
* @param f1 first fraction
* @param f2 second fraction
* @return a new fraction which is the quotient of f1 by f2
*/
public static PolynomialFraction divide(PolynomialFraction f1,
PolynomialFraction f2) {
PolynomialFraction copy = new PolynomialFraction(f1);
copy.divideSelf(f2);
return copy;
}
/** Invert the instance.
* Replace the instance by its inverse.
* @exception ArithmeticException if the instance is null
* @return the inverse of the instance
* @exception ArithmeticException if the instance is zero
*/
public void invertSelf() {
public PolynomialFraction invert() {
if (p.isZero()) {
throw new ArithmeticException("divide by zero");
}
Polynomial.Rational tmp = p;
p = q;
q = tmp;
RationalNumber[] a = p.getCoefficients();
PolynomialFraction inverse =
(a[a.length - 1].isNegative())
? new PolynomialFraction((Polynomial.Rational) q.negate(),
(Polynomial.Rational) p.negate())
: new PolynomialFraction(q, p);
inverse.simplify();
return inverse;
RationalNumber[] a = q.getCoefficients();
if (a[a.length - 1].isNegative()) {
p.negateSelf();
q.negateSelf();
}
simplify();
}
/** Invert a fraction.
* @param f fraction to invert
* @return a new fraction which is the inverse of f
*/
public static PolynomialFraction invert(PolynomialFraction f) {
PolynomialFraction copy = new PolynomialFraction(f);
copy.invertSelf();
return copy;
}
/** Simplify a fraction.
@ -318,8 +235,8 @@ public class PolynomialFraction {
*/
private void simplify() {
Polynomial.Rational a = new Polynomial.Rational(p);
Polynomial.Rational b = new Polynomial.Rational(q);
Polynomial.Rational a = p;
Polynomial.Rational b = q;
if (a.getDegree() < b.getDegree()) {
Polynomial.Rational tmp = a;
a = b;
@ -342,30 +259,28 @@ public class PolynomialFraction {
if (q.getDegree() == 0) {
if (! q.isOne()) {
RationalNumber f = q.getCoefficients()[0];
f.invertSelf();
p.multiplySelf(f);
p = (Polynomial.Rational) p.divide(q.getCoefficients()[0]);
q = new Polynomial.Rational(1l);
}
} else {
BigInteger lcm = p.getDenominatorsLCM();
if (lcm.compareTo(BigInteger.ONE) != 0) {
p.multiplySelf(lcm);
q.multiplySelf(lcm);
p = (Polynomial.Rational) p.multiply(lcm);
q = (Polynomial.Rational) q.multiply(lcm);
}
lcm = q.getDenominatorsLCM();
if (lcm.compareTo(BigInteger.ONE) != 0) {
p.multiplySelf(lcm);
q.multiplySelf(lcm);
p = (Polynomial.Rational) p.multiply(lcm);
q = (Polynomial.Rational) q.multiply(lcm);
}
}
if (q.getCoefficients()[q.getDegree()].isNegative()) {
p.negateSelf();
q.negateSelf();
p = (Polynomial.Rational) p.negate();
q = (Polynomial.Rational) q.negate();
}
}
@ -386,16 +301,6 @@ public class PolynomialFraction {
return q;
}
/** Set the name of the unknown (to appear during conversions to
* strings).
* @param name name to set (if null, the default 'x' value will be
* used)
*/
public void setUnknownName(String name) {
p.setUnknownName(name);
q.setUnknownName(name);
}
public String toString() {
if (p.isZero()) {
return "0";
@ -436,4 +341,6 @@ public class PolynomialFraction {
/** Denominator. */
private Polynomial.Rational q;
private static final long serialVersionUID = 6033909492898954748L;
}

349
src/mantissa/src/org/spaceroots/mantissa/algebra/RationalNumber.java
View File

@ -17,16 +17,25 @@
package org.spaceroots.mantissa.algebra;
import java.io.Serializable;
import java.math.BigInteger;
/**
* This class implements reduced rational numbers.
* <p>Instances of this class are immutable.</p>
* @version $Id: RationalNumber.java 1705 2006-09-17 19:57:39Z luc $
* @author L. Maisonobe
*/
public class RationalNumber {
public class RationalNumber implements Serializable {
/** Zero as a rational numer. */
public static final RationalNumber ZERO = new RationalNumber(0l);
/** One as a rational numer. */
public static final RationalNumber ONE = new RationalNumber(1l);
/**
* Simple constructor.
@ -37,63 +46,14 @@ public class RationalNumber {
q = BigInteger.ONE;
}
/**
* Simple constructor.
/** Simple constructor.
* Build a rational number from a numerator and a denominator.
* @param numerator numerator of the rational number
* @param denominator denominator of the rational number
* @exception ArithmeticException if the denominator is zero
*/
public RationalNumber(long numerator, long denominator) {
reset(numerator, denominator);
}
/**
* Simple constructor.
* Build a rational number from a numerator and a denominator.
* @param numerator numerator of the rational number
* @param denominator denominator of the rational number
* @exception ArithmeticException if the denominator is zero
*/
public RationalNumber(BigInteger numerator, BigInteger denominator) {
reset(numerator, denominator);
}
/**
* Simple constructor.
* Build a rational number from a single integer
* @param l value of the rational number
*/
public RationalNumber(long l) {
p = BigInteger.valueOf(l);
q = BigInteger.ONE;
}
/**
* Simple constructor.
* Build a rational number from a single integer
* @param i value of the rational number
*/
public RationalNumber(BigInteger i) {
p = i;
q = BigInteger.ONE;
}
/**
* Copy-constructor.
* @param r rational number to copy
*/
public RationalNumber(RationalNumber r) {
p = r.p;
q = r.q;
}
/** Reset the instance from a numerator and a denominator.
* @param numerator numerator of the rational number
* @param denominator denominator of the rational number
* @exception ArithmeticException if the denominator is zero
*/
public void reset(long numerator, long denominator) {
if (denominator == 0l) {
throw new ArithmeticException("divide by zero");
}
@ -110,12 +70,14 @@ public class RationalNumber {
}
/** Reset the instance from a numerator and a denominator.
/** Simple constructor.
* Build a rational number from a numerator and a denominator.
* @param numerator numerator of the rational number
* @param denominator denominator of the rational number
* @exception ArithmeticException if the denominator is zero
*/
public void reset(BigInteger numerator, BigInteger denominator) {
public RationalNumber(BigInteger numerator, BigInteger denominator) {
if (denominator.signum() == 0) {
throw new ArithmeticException("divide by zero");
}
@ -132,291 +94,166 @@ public class RationalNumber {
}
/** Reset the instance from a single integer
/** Simple constructor.
* Build a rational number from a single integer
* @param l value of the rational number
*/
public void reset(long l) {
public RationalNumber(long l) {
p = BigInteger.valueOf(l);
q = BigInteger.ONE;
}
/** Reset the instance from a single integer
/** Simple constructor.
* Build a rational number from a single integer
* @param i value of the rational number
*/
public void reset(BigInteger i) {
public RationalNumber(BigInteger i) {
p = i;
q = BigInteger.ONE;
}
/** Reset the instance from another rational number.
* @param r rational number to copy
/** Negate the instance.
* @return a new rational number, opposite to the isntance
*/
public void reset(RationalNumber r) {
p = r.p;
q = r.q;
public RationalNumber negate() {
return new RationalNumber(p.negate(), q);
}
/**
* Negate the instance
/** Add an integer to the instance.
* @param l integer to add
* @return a new rational number which is the sum of the instance and l
*/
public void negateSelf() {
p = p.negate();
public RationalNumber add(long l) {
return add(BigInteger.valueOf(l));
}
/**
* Negate a rational number.
* @param r rational number to negate
* @return a new rational number which is the opposite of r
/** Add an integer to the instance.
* @param l integer to add
* @return a new rational number which is the sum of the instance and l
*/
public static RationalNumber negate(RationalNumber r) {
RationalNumber copy = new RationalNumber(r);
copy.negateSelf();
return copy;
public RationalNumber add(BigInteger l) {
return new RationalNumber(p.add(q.multiply(l)), q);
}
/**
* Add a rational number to the instance.
* @param r rational number to add.
/** Add a rational number to the instance.
* @param r rational number to add
* @return a new rational number which is the sum of the instance and r
*/
public void addToSelf(RationalNumber r) {
p = p.multiply(r.q).add(r.p.multiply(q));
q = q.multiply(r.q);
simplify();
public RationalNumber add(RationalNumber r) {
return new RationalNumber(p.multiply(r.q).add(r.p.multiply(q)),
q.multiply(r.q));
}
/** Add two rational numbers.
* @param r1 first rational number
* @param r2 second rational number
* @return a new rational number which is the sum of r1 and r2
/** Subtract an integer from the instance.
* @param l integer to subtract
* @return a new rational number which is the difference the instance minus l
*/
public static RationalNumber add(RationalNumber r1, RationalNumber r2) {
return new RationalNumber(r1.p.multiply(r2.q).add(r2.p.multiply(r1.q)),
r1.q.multiply(r2.q));
public RationalNumber subtract(long l) {
return subtract(BigInteger.valueOf(l));
}
/**
* Subtract a rational number to the instance.
* @param r rational number to subtract.
/** Subtract an integer from the instance.
* @param l integer to subtract
* @return a new rational number which is the difference the instance minus l
*/
public void subtractFromSelf(RationalNumber r) {
p = p.multiply(r.q).subtract(r.p.multiply(q));
q = q.multiply(r.q);
simplify();
public RationalNumber subtract(BigInteger l) {
return new RationalNumber(p.subtract(q.multiply(l)), q);
}
/** Subtract two rational numbers.
* @param r1 first rational number
* @param r2 second rational number
* @return a new rational number which is the difference r1 minus r2
/** Subtract a rational number from the instance.
* @param r rational number to subtract
* @return a new rational number which is the difference the instance minus r
*/
public static RationalNumber subtract(RationalNumber r1, RationalNumber r2) {
return new RationalNumber(r1.p.multiply(r2.q).subtract(r2.p.multiply(r1.q)),
r1.q.multiply(r2.q));
public RationalNumber subtract(RationalNumber r) {
return new RationalNumber(p.multiply(r.q).subtract(r.p.multiply(q)),
q.multiply(r.q));
}
/** Multiply the instance by an integer.
* @param l integer to multiply by
* @return a new rational number which is the produc of the instance by l
*/
public void multiplySelf(long l) {
p = p.multiply(BigInteger.valueOf(l));
simplify();
public RationalNumber multiply(long l) {
return multiply(BigInteger.valueOf(l));
}
/** Multiply the instance by an integer.
* @param i integer to multiply by
*/
public void multiplySelf(BigInteger i) {
p = p.multiply(i);
simplify();
}
/** Multiply a rational number by an integer.
* @param l integer to multiply by
* @return a new rational number which is the produc of the instance by l
*/
public static RationalNumber multiply(RationalNumber r, long l) {
return new RationalNumber(r.p.multiply(BigInteger.valueOf(l)), r.q);
}
/** Multiply a rational number by an integer.
* @param i integer to multiply by
*/
public static RationalNumber multiply(RationalNumber r, BigInteger i) {
return new RationalNumber(r.p.multiply(i), r.q);
public RationalNumber multiply(BigInteger l) {
return new RationalNumber(p.multiply(l), q);
}
/** Multiply the instance by a rational number.
* @param r rational number to multiply by
* @param r rational number to multiply the instance with
* @return a new rational number which is the product of the instance and r
*/
public void multiplySelf(RationalNumber r) {
p = p.multiply(r.p);
q = q.multiply(r.q);
simplify();
}
/** Multiply two rational numbers.
* @param r1 first rational number
* @param r2 second rational number
* @return a new rational number which is the product of r1 and r2
*/
public static RationalNumber multiply(RationalNumber r1, RationalNumber r2) {
return new RationalNumber(r1.p.multiply(r2.p),
r1.q.multiply(r2.q));
public RationalNumber multiply(RationalNumber r) {
return new RationalNumber(p.multiply(r.p), q.multiply(r.q));
}
/** Divide the instance by an integer.
* @param l integer to divide by
* @return a new rational number which is the quotient of the instance by l
* @exception ArithmeticException if l is zero
*/
public void divideSelf(long l) {
if (l == 0l) {
throw new ArithmeticException("divide by zero");
} else if (l > 0l) {
q = q.multiply(BigInteger.valueOf(l));
} else {
p = p.negate();
q = q.multiply(BigInteger.valueOf(-l));
}
simplify();
public RationalNumber divide(long l) {
return divide(BigInteger.valueOf(l));
}
/** Divide the instance by an integer.
* @param i integer to divide by
* @param l integer to divide by
* @return a new rational number which is the quotient of the instance by l
* @exception ArithmeticException if l is zero
*/
public void divideSelf(BigInteger i) {
public RationalNumber divide(BigInteger l) {
if (i.signum() == 0) {
if (l.signum() == 0) {
throw new ArithmeticException("divide by zero");
} else if (i.signum() > 0) {
q = q.multiply(i);
} else {
p = p.negate();
q = q.multiply(i.negate());
}
simplify();
if (l.signum() > 0) {
return new RationalNumber(p, q.multiply(l));
}
}
return new RationalNumber(p.negate(), q.multiply(l.negate()));
/** Divide a rational number by an integer
* @param r rational number
* @param l integer
* @return a new rational number which is the quotient of r by l
* @exception ArithmeticException if l is zero
*/
public static RationalNumber divide(RationalNumber r, long l) {
RationalNumber copy = new RationalNumber(r);
copy.divideSelf(l);
return copy;
}
/** Divide a rational number by an integer
* @param r rational number
* @param i integer
* @return a new rational number which is the quotient of r by l
* @exception ArithmeticException if l is zero
*/
public static RationalNumber divide(RationalNumber r, BigInteger i) {
RationalNumber copy = new RationalNumber(r);
copy.divideSelf(i);
return copy;
}
/** Divide the instance by a rational number.
* @param r rational number to divide by
* @return a new rational number which is the quotient of the instance by r
* @exception ArithmeticException if r is zero
*/
public void divideSelf(RationalNumber r) {
public RationalNumber divide(RationalNumber r) {
if (r.p.signum() == 0) {
throw new ArithmeticException("divide by zero");
}
p = p.multiply(r.q);
q = q.multiply(r.p);
BigInteger newP = p.multiply(r.q);
BigInteger newQ = q.multiply(r.p);
if (q.signum() < 0) {
p = p.negate();
q = q.negate();
}
return (newQ.signum() < 0) ? new RationalNumber(newP.negate(),
newQ.negate())
: new RationalNumber(newP, newQ);
simplify();
}
/** Divide two rational numbers.
* @param r1 first rational number
* @param r2 second rational number
* @return a new rational number which is the quotient of r1 by r2
* @exception ArithmeticException if r2 is zero
*/
public static RationalNumber divide(RationalNumber r1, RationalNumber r2) {
RationalNumber copy = new RationalNumber(r1);
copy.divideSelf(r2);
return copy;
}
/** Invert the instance.
* Replace the instance by its inverse.
* @return the inverse of the instance
* @exception ArithmeticException if the instance is zero
*/
public void invertSelf() {
public RationalNumber invert() {
if (p.signum() == 0) {
throw new ArithmeticException("divide by zero");
}
BigInteger tmp = p;
p = q;
q = tmp;
return (q.signum() < 0) ? new RationalNumber(q.negate(), p.negate())
: new RationalNumber(q, p);
if (q.signum() < 0) {
p = p.negate();
q = q.negate();
}
}
/** Invert a rational number.
* @param r rational number to invert
* @return a new rational number which is the inverse of r
* @exception ArithmeticException if r is zero
*/
public static RationalNumber invert(RationalNumber r) {
return new RationalNumber(r.q, r.p);
}
/**
* Add the product of two rational numbers to the instance.
* This operation is equivalent to
* <code>addToSelf(RationalNumber.multiply(r1, r2))</code> except
* that no intermediate simplification is attempted.
* @param r1 first term of the product to add
* @param r2 second term of the product to add
*/
public void multiplyAndAddToSelf(RationalNumber r1, RationalNumber r2) {
BigInteger r1qr2q = r1.q.multiply(r2.q);
p = p.multiply(r1qr2q).add(r1.p.multiply(r2.p).multiply(q));
q = q.multiply(r1qr2q);
simplify();
}
/**
* Subtract the product of two rational numbers from the instance.
* This operation is equivalent to
* <code>subtractFromSelf(RationalNumber.multiply(r1, r2))</code>
* except that no intermediate simplification is attempted.
* @param r1 first term of the product to subtract
* @param r2 second term of the product to subtract
*/
public void multiplyAndSubtractFromSelf(RationalNumber r1, RationalNumber r2) {
BigInteger r1qr2q = r1.q.multiply(r2.q);
p = p.multiply(r1qr2q).subtract(r1.p.multiply(r2.p).multiply(q));
q = q.multiply(r1qr2q);
simplify();
}
/** Simplify a rational number by removing common factors.
@ -431,16 +268,14 @@ public class RationalNumber {
}
}
/**
* Get the numerator.
/** Get the numerator.
* @return the signed numerator
*/
public BigInteger getNumerator() {
return p;
}
/**
* Get the denominator.
/** Get the denominator.
* @return the denominator (always positive)
*/
public BigInteger getDenominator() {
@ -533,4 +368,6 @@ public class RationalNumber {
/** Denominator. */
private BigInteger q;
private static final long serialVersionUID = -324954393137577531L;
}

102
src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java
View File

@ -20,7 +20,6 @@ package org.spaceroots.mantissa.fitting;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import java.util.Iterator;
import org.spaceroots.mantissa.estimation.*;
@ -40,72 +39,30 @@ import org.spaceroots.mantissa.estimation.*;
public abstract class AbstractCurveFitter
implements EstimationProblem, Serializable {
/**
* Simple constructor.
/** Simple constructor.
* @param n number of coefficients in the underlying function
* @param maxIterations maximum number of iterations allowed
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
* @param steadyStateThreshold steady state detection threshold, the
* problem has reached a steady state (read converged) if
* <code>Math.abs (Jn - Jn-1) < Jn * convergence</code>, where
* <code>Jn</code> and <code>Jn-1</code> are the current and
* preceding criterion value (square sum of the weighted residuals
* of considered measurements).
* @param epsilon threshold under which the matrix of the linearized
* problem is considered singular (see {@link
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @param estimator estimator to use for the fitting
*/
protected AbstractCurveFitter(int n,
int maxIterations,
double convergence,
double steadyStateThreshold,
double epsilon) {
coefficients = new EstimatedParameter[n];
measurements = new ArrayList();
measurementsArray = null;
this.maxIterations = maxIterations;
this.steadyStateThreshold = steadyStateThreshold;
this.convergence = convergence;
this.epsilon = epsilon;
protected AbstractCurveFitter(int n, Estimator estimator) {
coefficients = new EstimatedParameter[n];
measurements = new ArrayList();
this.estimator = estimator;
}
/**
* Simple constructor.
/** Simple constructor.
* @param coefficients first estimate of the coefficients. A
* reference to this array is hold by the newly created object. Its
* elements will be adjusted during the fitting process and they will
* be set to the adjusted coefficients at the end.
* @param maxIterations maximum number of iterations allowed
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
* @param steadyStateThreshold steady state detection threshold, the
* problem has reached a steady state (read converged) if
* <code>Math.abs (Jn - Jn-1) < Jn * convergence</code>, where
* <code>Jn</code> and <code>Jn-1</code> are the current and
* preceding criterion value (square sum of the weighted residuals
* of considered measurements).
* @param epsilon threshold under which the matrix of the linearized
* problem is considered singular (see {@link
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @param estimator estimator to use for the fitting
*/
protected AbstractCurveFitter(EstimatedParameter[] coefficients,
int maxIterations,
double convergence,
double steadyStateThreshold,
double epsilon) {
Estimator estimator) {
this.coefficients = coefficients;
measurements = new ArrayList();
measurementsArray = null;
this.maxIterations = maxIterations;
this.steadyStateThreshold = steadyStateThreshold;
this.convergence = convergence;
this.epsilon = epsilon;
this.coefficients = coefficients;
measurements = new ArrayList();
this.estimator = estimator;
}
/** Add a weighted (x,y) pair to the sample.
@ -114,7 +71,6 @@ public abstract class AbstractCurveFitter
* @param y ordinate, we have <code>y = f (x)</code>
*/
public void addWeightedPair(double weight, double x, double y) {
measurementsArray = null;
measurements.add(new FitMeasurement(weight, x, y));
}
@ -131,9 +87,8 @@ public abstract class AbstractCurveFitter
*/
public double[] fit()
throws EstimationException {
// perform the fit using a linear least square estimator
new GaussNewtonEstimator(maxIterations, convergence,
steadyStateThreshold, epsilon).estimate(this);
// perform the fit
estimator.estimate(this);
// extract the coefficients
double[] fittedCoefficients = new double[coefficients.length];
@ -146,14 +101,7 @@ public abstract class AbstractCurveFitter
}
public WeightedMeasurement[] getMeasurements() {
if (measurementsArray == null) {
measurementsArray = new FitMeasurement[measurements.size()];
int i = 0;
for (Iterator iterator = measurements.iterator(); iterator.hasNext(); ++i) {
measurementsArray[i] = (FitMeasurement) iterator.next();
}
}
return measurementsArray;
return (WeightedMeasurement[]) measurements.toArray(new FitMeasurement[measurements.size()]);
}
/** Get the unbound parameters of the problem.
@ -161,14 +109,14 @@ public abstract class AbstractCurveFitter
* @return unbound parameters
*/
public EstimatedParameter[] getUnboundParameters() {
return coefficients;
return (EstimatedParameter[]) coefficients.clone();
}
/** Get all the parameters of the problem.
* @return parameters
*/
public EstimatedParameter[] getAllParameters() {
return coefficients;
return (EstimatedParameter[]) coefficients.clone();
}
/** Utility method to sort the measurements with respect to the abscissa.
@ -205,10 +153,6 @@ public abstract class AbstractCurveFitter
}
}
// make sure subsequent calls to getMeasurements
// will not use the unsorted array
measurementsArray = null;
}
/** Get the value of the function at x according to the current parameters value.
@ -272,15 +216,7 @@ public abstract class AbstractCurveFitter
/** Measurements vector */
protected List measurements;
/** Measurements array.
* This array contains the same entries as measurements_, but in a
* different structure.
*/
private FitMeasurement[] measurementsArray;
private int maxIterations;
private double convergence;
private double steadyStateThreshold;
private double epsilon;
/** Estimator for the fitting problem. */
private Estimator estimator;
}

11
src/mantissa/src/org/spaceroots/mantissa/fitting/F2FP2Iterator.java
View File

@ -58,12 +58,13 @@ class F2FP2Iterator
// get the raw values from the underlying FFPIterator
VectorialValuedPair point = ffpIterator.nextSamplePoint();
double[] y = point.y;
// hack the values (to avoid building a new object)
double[] y = point.getY();
y[0] *= y[0];
y[1] *= y[1];
return point;
// square the values
return new VectorialValuedPair(point.x,
new double[] {
y[0] * y[0], y[1] * y[1]
});
}

124
src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java
View File

@ -32,32 +32,32 @@ import org.spaceroots.mantissa.estimation.EstimationException;
* <p>The algorithm used to guess the coefficients is as follows:</p>
* <p>We know f (t) at some sampling points ti and want to find a,
* omega and phi such that f (t) = a cos (omega t + phi).
* <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
* &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
* </p>
*
* <p>From the analytical expression, we can compute two primitives :
* <pre>
* If2 (t) = int (f^2) = a^2 * [t + S (t)] / 2
* If'2 (t) = int (f'^2) = a^2 * omega^2 * [t - S (t)] / 2
* where S (t) = sin (2 * (omega * t + phi)) / (2 * omega)
* If2 (t) = &int; f<sup>2</sup> = a<sup>2</sup> &times; [t + S (t)] / 2
* If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
* where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
* </pre>
* </p>
*
* <p>We can remove S between these expressions :
* <pre>
* If'2 (t) = a^2 * omega ^ 2 * t - omega ^ 2 * If2 (t)
* If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
* </pre>
* </p>
*
* <p>The preceding expression shows that If'2 (t) is a linear
* combination of both t and If2 (t): If'2 (t) = A * t + B * If2 (t)
* combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
* </p>
*
* <p>From the primitive, we can deduce the same form for definite
* integrals between t1 and ti for each ti :
* integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
* <pre>
* If2 (ti) - If2 (t1) = A * (ti - t1) + B * (If2 (ti) - If2 (t1))
* If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
* </pre>
* </p>
*
@ -66,62 +66,60 @@ import org.spaceroots.mantissa.estimation.EstimationException;
* each sample points.
* </p>
*
* <p>For a bilinear expression z (xi, yi) = A * xi + B * yi, the
* coefficients a and b that minimize a least square criterion
* Sum ((zi - z (xi, yi))^2) are given by these expressions:</p>
* <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
* coefficients A and B that minimize a least square criterion
* &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
* <pre>
*
* Sum (yi^2) Sum (xi zi) - Sum (xi yi) Sum (yi zi)
* A = ------------------------------------------------
* Sum (xi^2) Sum (yi^2) - Sum (xi yi) Sum (xi yi)
* &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* A = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* Sum (xi^2) Sum (yi zi) - Sum (xi yi) Sum (xi zi)
* B = ------------------------------------------------
* Sum (xi^2) Sum (yi^2) - Sum (xi yi) Sum (xi yi)
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
* B = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
* </pre>
* </p>
*
*
* <p>In fact, we can assume both a and omega are positive and
* compute them directly, knowing that A = a^2 * omega^2 and that
* B = - omega^2. The complete algorithm is therefore:</p>
* <p>In fact, we can assume both a and &omega; are positive and
* compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
* B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
* <pre>
*
* for each ti from t1 to t(n-1), compute:
* f (ti)
* f' (ti) = (f (t(i+1)) - f(t(i-1))) / (t(i+1) - t(i-1))
* xi = ti - t1
* yi = int (f^2) from t1 to ti
* zi = int (f'^2) from t1 to ti
* update the sums Sum (xi^2), Sum (yi^2),
* Sum (xi yi), Sum (xi zi)
* and Sum (yi zi)
* for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
* f (t<sub>i</sub>)
* f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
* x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
* y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
* end for
*
* |-------------------------------------------------
* \ | Sum (yi^2) Sum (xi zi) - Sum (xi yi) Sum (yi zi)
* a = \ | ------------------------------------------------
* \| Sum (xi yi) Sum (xi zi) - Sum (xi^2) Sum (yi zi)
* |--------------------------
* \ | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* a = \ | ------------------------
* \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
*
*
* |-------------------------------------------------
* \ | Sum (xi yi) Sum (xi zi) - Sum (xi^2) Sum (yi zi)
* omega = \ | ------------------------------------------------
* \| Sum (xi^2) Sum (yi^2) - Sum (xi yi) Sum (xi yi)
* |--------------------------
* \ | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* &omega; = \ | ------------------------
* \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* </pre>
* </p>
* <p>Once we know omega, we can compute:
* <p>Once we know &omega;, we can compute:
* <pre>
* fc = omega * f (t) * cos (omega * t) - f' (t) * sin (omega * t)
* fs = omega * f (t) * sin (omega * t) + f' (t) * cos (omega * t)
* fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
* fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
* </pre>
* </p>
* <p>It appears that <code>fc = a * omega * cos (phi)</code> and
* <code>fs = -a * omega * sin (phi)</code>, so we can use these
* expressions to compute phi. The best estimate over the sample is
* <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
* <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
* expressions to compute &phi;. The best estimate over the sample is
* given by averaging these expressions.
* </p>
@ -138,7 +136,8 @@ public class HarmonicCoefficientsGuesser
implements Serializable{
public HarmonicCoefficientsGuesser(AbstractCurveFitter.FitMeasurement[] measurements) {
this.measurements = measurements;
this.measurements =
(AbstractCurveFitter.FitMeasurement[]) measurements.clone();
a = Double.NaN;
omega = Double.NaN;
}
@ -159,7 +158,7 @@ public class HarmonicCoefficientsGuesser
guessPhi();
}
/** Estimate a first guess of the a and omega coefficients.
/** Estimate a first guess of the a and &omega; coefficients.
* @exception ExhaustedSampleException if the sample is exhausted.
@ -185,25 +184,25 @@ public class HarmonicCoefficientsGuesser
SampledFunctionIterator sampler =
new EnhancedSimpsonIntegratorSampler(iter);
VectorialValuedPair p0 = sampler.nextSamplePoint();
double p0X = p0.getX();
double[] p0Y = p0.getY();
double p0X = p0.x;
double[] p0Y = p0.y;
// get the points for the linear model
while (sampler.hasNext()) {
VectorialValuedPair point = sampler.nextSamplePoint();
double pX = point.getX();
double[] pY = point.getY();
double pX = point.x;
double[] pY = point.y;
double dx = pX - p0X;
double dy0 = pY[0] - p0Y[0];
double dy1 = pY[1] - p0Y[1];
double x = pX - p0X;
double y = pY[0] - p0Y[0];
double z = pY[1] - p0Y[1];
sx2 += dx * dx;
sy2 += dy0 * dy0;
sxy += dx * dy0;
sxz += dx * dy1;
syz += dy0 * dy1;
sx2 += x * x;
sy2 += y * y;
sxy += x * y;
sxz += x * z;
syz += y * z;
}
@ -219,7 +218,7 @@ public class HarmonicCoefficientsGuesser
}
/** Estimate a first guess of the phi coefficient.
/** Estimate a first guess of the &phi; coefficient.
* @exception ExhaustedSampleException if the sample is exhausted.
@ -237,12 +236,11 @@ public class HarmonicCoefficientsGuesser
while (iter.hasNext()) {
VectorialValuedPair point = iter.nextSamplePoint();
double omegaX = omega * point.getX();
double[] pY = point.getY();
double omegaX = omega * point.x;
double cosine = Math.cos(omegaX);
double sine = Math.sin(omegaX);
fcMean += omega * pY[0] * cosine - pY[1] * sine;
fsMean += omega * pY[0] * sine + pY[1] * cosine;
fcMean += omega * point.y[0] * cosine - point.y[1] * sine;
fsMean += omega * point.y[0] * sine + point.y[1] * cosine;
}
phi = Math.atan2(-fsMean, fcMean);

58
src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java
View File

@ -17,9 +17,10 @@
package org.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.estimation.*;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.FunctionException;
@ -40,8 +41,37 @@ import org.spaceroots.mantissa.functions.FunctionException;
*/
public class HarmonicFitter
extends AbstractCurveFitter
implements EstimationProblem, Serializable {
extends AbstractCurveFitter {
/** Simple constructor.
* @param estimator estimator to use for the fitting
*/
public HarmonicFitter(Estimator estimator) {
super(3, estimator);
coefficients[0] = new EstimatedParameter("a", 2.0 * Math.PI);
coefficients[1] = new EstimatedParameter("omega", 0.0);
coefficients[2] = new EstimatedParameter("phi", 0.0);
firstGuessNeeded = true;
}
/**
* Simple constructor.
* <p>This constructor can be used when a first estimate of the
* coefficients is already known.</p>
* @param coefficients first estimate of the coefficients.
* A reference to this array is hold by the newly created
* object. Its elements will be adjusted during the fitting process
* and they will be set to the adjusted coefficients at the end.
* @param estimator estimator to use for the fitting
*/
public HarmonicFitter(EstimatedParameter[] coefficients,
Estimator estimator) {
super(coefficients, estimator);
firstGuessNeeded = false;
}
/**
* Simple constructor.
@ -58,14 +88,13 @@ public class HarmonicFitter
* problem is considered singular (see {@link
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @deprecated replaced by {@link #HarmonicFitter(Estimator)}
* as of version 7.0
*/
public HarmonicFitter(int maxIterations, double convergence,
double steadyStateThreshold, double epsilon) {
super(3, maxIterations, convergence, steadyStateThreshold, epsilon);
coefficients[0] = new EstimatedParameter("a", 2.0 * Math.PI);
coefficients[1] = new EstimatedParameter("omega", 0.0);
coefficients[2] = new EstimatedParameter("phi", 0.0);
firstGuessNeeded = true;
this(new GaussNewtonEstimator(maxIterations, convergence,
steadyStateThreshold, epsilon));
}
/**
@ -92,14 +121,15 @@ public class HarmonicFitter
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @deprecated replaced by {@link #HarmonicFitter(EstimatedParameter[],
* Estimator)} as of version 7.0
*/
public HarmonicFitter(EstimatedParameter[] coefficients,
int maxIterations, double convergence,
double steadyStateThreshold, double epsilon) {
super(coefficients,
maxIterations, convergence,
steadyStateThreshold, epsilon);
firstGuessNeeded = false;
this(coefficients,
new GaussNewtonEstimator(maxIterations, convergence,
steadyStateThreshold, epsilon));
}
public double[] fit()

82
src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java
View File

@ -17,9 +17,9 @@
package org.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.estimation.*;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
/** This class implements a curve fitting specialized for polynomials.
@ -38,11 +38,47 @@ import org.spaceroots.mantissa.estimation.*;
*/
public class PolynomialFitter
extends AbstractCurveFitter
implements EstimationProblem, Serializable {
extends AbstractCurveFitter {
/**
* Simple constructor.
/** Simple constructor.
* <p>The polynomial fitter built this way are complete polynoms,
* ie. a n-degree polynom has n+1 coefficients. In order to build
* fitter for sparse polynoms (for example <code>a x^20 - b
* x^30</code>, on should first build the coefficients array and
* provide it to {@link
* #PolynomialFitter(PolynomialCoefficient[], int, double, double,
* double)}.</p>
* @param degree maximal degree of the polynom
* @param estimator estimator to use for the fitting
*/
public PolynomialFitter(int degree, Estimator estimator) {
super(degree + 1, estimator);
for (int i = 0; i < coefficients.length; ++i) {
coefficients[i] = new PolynomialCoefficient(i);
}
}
/** Simple constructor.
* <p>This constructor can be used either when a first estimate of
* the coefficients is already known (which is of little interest
* because the fit cost is the same whether a first guess is known or
* not) or when one needs to handle sparse polynoms like <code>a
* x^20 - b x^30</code>.</p>
* @param coefficients first estimate of the coefficients.
* A reference to this array is hold by the newly created
* object. Its elements will be adjusted during the fitting process
* and they will be set to the adjusted coefficients at the end.
* @param estimator estimator to use for the fitting
*/
public PolynomialFitter(PolynomialCoefficient[] coefficients,
Estimator estimator) {
super(coefficients, estimator);
}
/** Simple constructor.
* <p>The polynomial fitter built this way are complete polynoms,
* ie. a n-degree polynom has n+1 coefficients. In order to build
@ -66,23 +102,18 @@ public class PolynomialFitter
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @deprecated replaced by {@link #PolynomialFitter(int,Estimator)}
* as of version 7.0
*/
public PolynomialFitter(int degree,
int maxIterations, double convergence,
double steadyStateThreshold, double epsilon) {
super(degree + 1,
maxIterations, steadyStateThreshold,
convergence, epsilon);
for (int i = 0; i < coefficients.length; ++i) {
coefficients[i] = new PolynomialCoefficient(i);
}
this(degree,
new GaussNewtonEstimator(maxIterations, steadyStateThreshold,
convergence, epsilon));
}
/**
* Simple constructor.
/** Simple constructor.
* <p>This constructor can be used either when a first estimate of
* the coefficients is already known (which is of little interest
@ -108,13 +139,15 @@ public class PolynomialFitter
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(
* org.spaceroots.mantissa.linalg.Matrix,double) SquareMatrix.solve}).
* @deprecated replaced by {@link #PolynomialFitter(PolynomialCoefficient[],
* Estimator)} as of version 7.0
*/
public PolynomialFitter(PolynomialCoefficient[] coefficients,
int maxIterations, double convergence,
double steadyStateThreshold, double epsilon) {
super(coefficients,
maxIterations, steadyStateThreshold,
convergence, epsilon);
this(coefficients,
new GaussNewtonEstimator(maxIterations, steadyStateThreshold,
convergence, epsilon));
}
/** Get the value of the function at x according to the current parameters value.
@ -135,9 +168,12 @@ public class PolynomialFitter
* @return partial derivative
*/
public double partial(double x, EstimatedParameter p) {
return Math.pow(x, ((PolynomialCoefficient) p).degree);
if (p instanceof PolynomialCoefficient) {
return Math.pow(x, ((PolynomialCoefficient) p).degree);
}
throw new RuntimeException("internal error");
}
private static final long serialVersionUID = -226724596015163603L;
private static final long serialVersionUID = -744904084649890769L;
}

4
src/mantissa/src/org/spaceroots/mantissa/functions/scalar/ComputableFunction.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.functions.scalar;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
/** This interface represents scalar functions of one real variable.
@ -42,7 +44,7 @@ import org.spaceroots.mantissa.functions.FunctionException;
* @author L. Maisonobe
*/
public interface ComputableFunction {
public interface ComputableFunction extends Serializable {
/** Get the value of the function at the specified abscissa.
* @param x current abscissa

4
src/mantissa/src/org/spaceroots/mantissa/functions/scalar/SampledFunction.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.functions.scalar;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
/** This interface represent sampled scalar functions.
@ -45,7 +47,7 @@ import org.spaceroots.mantissa.functions.FunctionException;
* @author L. Maisonobe
*/
public interface SampledFunction {
public interface SampledFunction extends Serializable {
/** Get the number of points in the sample.
* @return number of points in the sample

4
src/mantissa/src/org/spaceroots/mantissa/functions/vectorial/ComputableFunction.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.functions.vectorial;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
/** This interface represents vectorial functions of one real variable.
@ -42,7 +44,7 @@ import org.spaceroots.mantissa.functions.FunctionException;
* @author L. Maisonobe
*/
public interface ComputableFunction {
public interface ComputableFunction extends Serializable {
/** Get the dimension of the vectorial values of the function.
* @return dimension
*/

4
src/mantissa/src/org/spaceroots/mantissa/functions/vectorial/SampledFunction.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.functions.vectorial;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
/** This interface represent sampled vectorial functions.
@ -45,7 +47,7 @@ import org.spaceroots.mantissa.functions.FunctionException;
* @author L. Maisonobe
*/
public interface SampledFunction {
public interface SampledFunction extends Serializable {
/** Get the number of points in the sample.
* @return number of points in the sample

51
src/mantissa/src/org/spaceroots/mantissa/functions/vectorial/VectorialValuedPair.java
View File

@ -22,7 +22,7 @@ import java.io.Serializable;
/** This class represents an (x, f(x)) pair for vectorial functions.
* <p>A vectorial function is a function of one vectorial parameter x whose
* value is a vector. This class is used has a simple placeholder to
* value is a vector. This class is used has a simple immutable placeholder to
* contain both an abscissa and the value of the function at this
* abscissa.</p>
@ -44,56 +44,15 @@ public class VectorialValuedPair
*/
public VectorialValuedPair(double x, double[] y) {
this.x = x;
this.y = y;
}
/**
* Copy-constructor.
* @param p point to copy
*/
public VectorialValuedPair(VectorialValuedPair p) {
x = p.x;
y = p.y;
}
/**
* Getter for the abscissa.
* @return value of the abscissa
*/
public double getX() {
return x;
}
/**
* Getter for the ordinate.
* @return value of the ordinate
*/
public double[] getY() {
return y;
}
/**
* Setter for the abscissa.
* @param x new value for the abscissa
*/
public void setX(double x) {
this.x = x;
}
/**
* Setter for the ordinate.
* @param y new value for the ordinate
*/
public void setY(double[] y) {
this.y = y;
this.y = (double[]) y.clone();
}
/** Abscissa of the point. */
private double x;
public final double x;
/** Vectorial ordinate of the point, y = f (x). */
private double[] y;
public final double[] y;
private static final long serialVersionUID = -1336411215846160578L;
private static final long serialVersionUID = -7397116933564410103L;
}

101
src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
View File

@ -32,8 +32,8 @@ import java.io.Serializable;
* <code>Rotation</code> instance (see the various constructors and
* getters). In addition, a rotation can also be built implicitely
* from a set of vectors and their image.</p>
* <p>This implies that this class can be used to transform one
* representation into another one. For example, converting a rotation
* <p>This implies that this class can be used to convert from one
* representation to another one. For example, converting a rotation
* matrix into a set of Cardan angles from can be done using the
* followong single line of code:</p>
* <pre>
@ -463,7 +463,7 @@ public class Rotation implements Serializable {
* <p>Beware that many people routinely use the term Euler angles even
* for what really are Cardan angles (this confusion is especially
* widespread in the aerospace business where Roll, Pitch and Yaw angles
* are often tagged as Euler angles).</p>
* are often wrongly tagged as Euler angles).</p>
* @param order order of rotations to use
* @param alpha1 angle of the first elementary rotation
@ -472,94 +472,19 @@ public class Rotation implements Serializable {
*/
public Rotation(RotationOrder order,
double alpha1, double alpha2, double alpha3) {
this(computeRotation(order, alpha1, alpha2, alpha3));
}
/** Copy constructor.
* <p>This constructor is used only for the sake of Cardan/Euler
* angles handling.</p>
* @param r rotation to copy
*/
private Rotation(Rotation r) {
q0 = r.q0;
q1 = r.q1;
q2 = r.q2;
q3 = r.q3;
}
/** Build a rotation from three Cardan or Euler elementary rotations.
* @param order order of rotations to use
* @param alpha1 angle of the first elementary rotation
* @param alpha2 angle of the second elementary rotation
* @param alpha3 angle of the third elementary rotation
*/
public static Rotation computeRotation(RotationOrder order,
double alpha1,
double alpha2,
double alpha3) {
if (order == RotationOrder.XYZ) {
return compose(new Rotation(Vector3D.plusI, alpha1),
new Rotation(Vector3D.plusJ, alpha2),
new Rotation(Vector3D.plusK, alpha3));
} else if (order == RotationOrder.XZY) {
return compose(new Rotation(Vector3D.plusI, alpha1),
new Rotation(Vector3D.plusK, alpha2),
new Rotation(Vector3D.plusJ, alpha3));
} else if (order == RotationOrder.YXZ) {
return compose(new Rotation(Vector3D.plusJ, alpha1),
new Rotation(Vector3D.plusI, alpha2),
new Rotation(Vector3D.plusK, alpha3));
} else if (order == RotationOrder.YZX) {
return compose(new Rotation(Vector3D.plusJ, alpha1),
new Rotation(Vector3D.plusK, alpha2),
new Rotation(Vector3D.plusI, alpha3));
} else if (order == RotationOrder.ZXY) {
return compose(new Rotation(Vector3D.plusK, alpha1),
new Rotation(Vector3D.plusI, alpha2),
new Rotation(Vector3D.plusJ, alpha3));
} else if (order == RotationOrder.ZYX) {
return compose(new Rotation(Vector3D.plusK, alpha1),
new Rotation(Vector3D.plusJ, alpha2),
new Rotation(Vector3D.plusI, alpha3));
} else if (order == RotationOrder.XYX) {
return compose(new Rotation(Vector3D.plusI, alpha1),
new Rotation(Vector3D.plusJ, alpha2),
new Rotation(Vector3D.plusI, alpha3));
} else if (order == RotationOrder.XZX) {
return compose(new Rotation(Vector3D.plusI, alpha1),
new Rotation(Vector3D.plusK, alpha2),
new Rotation(Vector3D.plusI, alpha3));
} else if (order == RotationOrder.YXY) {
return compose(new Rotation(Vector3D.plusJ, alpha1),
new Rotation(Vector3D.plusI, alpha2),
new Rotation(Vector3D.plusJ, alpha3));
} else if (order == RotationOrder.YZY) {
return compose(new Rotation(Vector3D.plusJ, alpha1),
new Rotation(Vector3D.plusK, alpha2),
new Rotation(Vector3D.plusJ, alpha3));
} else if (order == RotationOrder.ZXZ) {
return compose(new Rotation(Vector3D.plusK, alpha1),
new Rotation(Vector3D.plusI, alpha2),
new Rotation(Vector3D.plusK, alpha3));
} else { // last possibility is ZYZ
return compose(new Rotation(Vector3D.plusK, alpha1),
new Rotation(Vector3D.plusJ, alpha2),
new Rotation(Vector3D.plusK, alpha3));
}
}
/** Override the instance by the composition of three rotations.
* @param r3 last (outermost) rotation to compose
* @param r2 intermediate rotation to compose
* @param r1 first (innermost) rotation to compose
*/
private static Rotation compose(Rotation r3, Rotation r2, Rotation r1) {
return r3.applyTo(r2.applyTo(r1));
Rotation r1 = new Rotation(order.getA1(), alpha1);
Rotation r2 = new Rotation(order.getA2(), alpha2);
Rotation r3 = new Rotation(order.getA3(), alpha3);
Rotation composed = r1.applyTo(r2.applyTo(r3));
q0 = composed.q0;
q1 = composed.q1;
q2 = composed.q2;
q3 = composed.q3;
}
/** Revert a rotation.
* Build a rotation which reverse the effect of another
* rotation. This means that is r(u) = v, then r.revert (v) = u. The
* rotation. This means that if r(u) = v, then r.revert(v) = u. The
* instance is not changed.
* @return a new rotation whose effect is the reverse of the effect
* of the instance
@ -1087,6 +1012,6 @@ public class Rotation implements Serializable {
/** Third coordinate of the vectorial part of the quaternion. */
private final double q3;
private static final long serialVersionUID = 7264384082212242475L;
private static final long serialVersionUID = 5127795878493115119L;
}

75
src/mantissa/src/org/spaceroots/mantissa/geometry/RotationOrder.java
View File

@ -36,9 +36,16 @@ public final class RotationOrder {
* This is a utility class that cannot be instantiated by the user,
* so its only constructor is private.
* @param name name of the rotation order
* @param a1 axis of the first rotation
* @param a2 axis of the second rotation
* @param a3 axis of the third rotation
*/
private RotationOrder(String name) {
private RotationOrder(String name,
Vector3D a1, Vector3D a2, Vector3D a3) {
this.name = name;
this.a1 = a1;
this.a2 = a2;
this.a3 = a3;
}
/** Get a string representation of the instance.
@ -48,79 +55,121 @@ public final class RotationOrder {
return name;
}
/** Get the axis of the first rotation.
* @return axis of the first rotation
*/
public Vector3D getA1() {
return a1;
}
/** Get the axis of the second rotation.
* @return axis of the second rotation
*/
public Vector3D getA2() {
return a2;
}
/** Get the axis of the second rotation.
* @return axis of the second rotation
*/
public Vector3D getA3() {
return a3;
}
/** Set of Cardan angles.
* this ordered set of rotations is around X, then around Y, then
* around Z
*/
public static final RotationOrder XYZ = new RotationOrder("XYZ");
public static final RotationOrder XYZ =
new RotationOrder("XYZ", Vector3D.plusI, Vector3D.plusJ, Vector3D.plusK);
/** Set of Cardan angles.
* this ordered set of rotations is around X, then around Z, then
* around Y
*/
public static final RotationOrder XZY = new RotationOrder("XZY");
public static final RotationOrder XZY =
new RotationOrder("XZY", Vector3D.plusI, Vector3D.plusK, Vector3D.plusJ);
/** Set of Cardan angles.
* this ordered set of rotations is around Y, then around X, then
* around Z
*/
public static final RotationOrder YXZ = new RotationOrder("YXZ");
public static final RotationOrder YXZ =
new RotationOrder("YXZ", Vector3D.plusJ, Vector3D.plusI, Vector3D.plusK);
/** Set of Cardan angles.
* this ordered set of rotations is around Y, then around Z, then
* around X
*/
public static final RotationOrder YZX = new RotationOrder("YZX");
public static final RotationOrder YZX =
new RotationOrder("YZX", Vector3D.plusJ, Vector3D.plusK, Vector3D.plusI);
/** Set of Cardan angles.
* this ordered set of rotations is around Z, then around X, then
* around Y
*/
public static final RotationOrder ZXY = new RotationOrder("ZXY");
public static final RotationOrder ZXY =
new RotationOrder("ZXY", Vector3D.plusK, Vector3D.plusI, Vector3D.plusJ);
/** Set of Cardan angles.
* this ordered set of rotations is around Z, then around Y, then
* around X
*/
public static final RotationOrder ZYX = new RotationOrder("ZYX");
public static final RotationOrder ZYX =
new RotationOrder("ZYX", Vector3D.plusK, Vector3D.plusJ, Vector3D.plusI);
/** Set of Euler angles.
* this ordered set of rotations is around X, then around Y, then
* around X
*/
public static final RotationOrder XYX = new RotationOrder("XYX");
public static final RotationOrder XYX =
new RotationOrder("XYX", Vector3D.plusI, Vector3D.plusJ, Vector3D.plusI);
/** Set of Euler angles.
* this ordered set of rotations is around X, then around Z, then
* around X
*/
public static final RotationOrder XZX = new RotationOrder("XZX");
public static final RotationOrder XZX =
new RotationOrder("XZX", Vector3D.plusI, Vector3D.plusK, Vector3D.plusI);
/** Set of Euler angles.
* this ordered set of rotations is around Y, then around X, then
* around Y
*/
public static final RotationOrder YXY = new RotationOrder("YXY");
public static final RotationOrder YXY =
new RotationOrder("YXY", Vector3D.plusJ, Vector3D.plusI, Vector3D.plusJ);
/** Set of Euler angles.
* this ordered set of rotations is around Y, then around Z, then
* around Y
*/
public static final RotationOrder YZY = new RotationOrder("YZY");
public static final RotationOrder YZY =
new RotationOrder("YZY", Vector3D.plusJ, Vector3D.plusK, Vector3D.plusJ);
/** Set of Euler angles.
* this ordered set of rotations is around Z, then around X, then
* around Z
*/
public static final RotationOrder ZXZ = new RotationOrder("ZXZ");
public static final RotationOrder ZXZ =
new RotationOrder("ZXZ", Vector3D.plusK, Vector3D.plusI, Vector3D.plusK);
/** Set of Euler angles.
* this ordered set of rotations is around Z, then around Y, then
* around Z
*/
public static final RotationOrder ZYZ = new RotationOrder("ZYZ");
public static final RotationOrder ZYZ =
new RotationOrder("ZYZ", Vector3D.plusK, Vector3D.plusJ, Vector3D.plusK);
/** Name of the rotations order. */
private final String name;
/** Axis of the first rotation. */
private final Vector3D a1;
/** Axis of the second rotation. */
private final Vector3D a2;
/** Axis of the third rotation. */
private final Vector3D a3;
}

177
src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java
View File

@ -19,13 +19,14 @@ package org.spaceroots.mantissa.geometry;
import java.io.Serializable;
/** This class implements vectors in a three-dimensional space.
* <p>Vector3D are guaranteed to be immutable objects.</p>
* @version $Id: Vector3D.java 1705 2006-09-17 19:57:39Z luc $
* <p>Instance of this class are guaranteed to be immutable.</p>
* @version $Id: Vector3D.java 1716 2006-12-13 22:56:35Z luc $
* @author L. Maisonobe
*/
public class Vector3D implements Serializable {
public class Vector3D
implements Serializable {
/** First canonical vector (coordinates : 1, 0, 0). */
public static final Vector3D plusI = new Vector3D(1, 0, 0);
@ -42,7 +43,7 @@ public class Vector3D implements Serializable {
/** Third canonical vector (coordinates : 0, 0, 1). */
public static final Vector3D plusK = new Vector3D(0, 0, 1);
/** Opposite of the third canonical vector (coordinates : 0, 0, -1). */
/** Opposite of the third canonical vector (coordinates : 0, 0, -1). */
public static final Vector3D minusK = new Vector3D(0, 0, -1);
/** Simple constructor.
@ -98,58 +99,57 @@ public class Vector3D implements Serializable {
/** Linear constructor
* Build a vector from two other ones and corresponding scale factors.
* The vector built will be a * u + b * v
* @param a first scale factor
* @param u first base (unscaled) vector
* @param b second scale factor
* @param v second base (unscaled) vector
* The vector built will be a1 * u1 + a2 * u2
* @param a1 first scale factor
* @param u1 first base (unscaled) vector
* @param a2 second scale factor
* @param u2 second base (unscaled) vector
*/
public Vector3D(double a, Vector3D u, double b, Vector3D v) {
this.x = a * u.x + b * v.x;
this.y = a * u.y + b * v.y;
this.z = a * u.z + b * v.z;
public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2) {
this.x = a1 * u1.x + a2 * u2.x;
this.y = a1 * u1.y + a2 * u2.y;
this.z = a1 * u1.z + a2 * u2.z;
}
/** Linear constructor
* Build a vector from three other ones and corresponding scale factors.
* The vector built will be a * u + b * v + c * w
* @param a first scale factor
* @param u first base (unscaled) vector
* @param b second scale factor
* @param v second base (unscaled) vector
* @param c third scale factor
* @param w third base (unscaled) vector
* The vector built will be a1 * u1 + a2 * u2 + a3 * u3
* @param a1 first scale factor
* @param u1 first base (unscaled) vector
* @param a2 second scale factor
* @param u2 second base (unscaled) vector
* @param a3 third scale factor
* @param u3 third base (unscaled) vector
*/
public Vector3D(double a, Vector3D u, double b, Vector3D v,
double c, Vector3D w) {
this.x = a * u.x + b * v.x + c * w.x;
this.y = a * u.y + b * v.y + c * w.y;
this.z = a * u.z + b * v.z + c * w.z;
public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2,
double a3, Vector3D u3) {
this.x = a1 * u1.x + a2 * u2.x + a3 * u3.x;
this.y = a1 * u1.y + a2 * u2.y + a3 * u3.y;
this.z = a1 * u1.z + a2 * u2.z + a3 * u3.z;
}
/** Linear constructor
* Build a vector from four other ones and corresponding scale factors.
* The vector built will be a * t + b * u + c * v + d * w
* @param a first scale factor
* @param t first base (unscaled) vector
* @param b second scale factor
* @param u second base (unscaled) vector
* @param c third scale factor
* @param v third base (unscaled) vector
* @param d third scale factor
* @param w third base (unscaled) vector
* The vector built will be a1 * u1 + a2 * u2 + a3 * u3 + a4 * u4
* @param a1 first scale factor
* @param u1 first base (unscaled) vector
* @param a2 second scale factor
* @param u2 second base (unscaled) vector
* @param a3 third scale factor
* @param u3 third base (unscaled) vector
* @param a4 fourth scale factor
* @param u4 fourth base (unscaled) vector
*/
public Vector3D(double a, Vector3D t, double b, Vector3D u,
double c, Vector3D v, double d, Vector3D w) {
this.x = a * t.x + b * u.x + c * v.x + d * w.x;
this.y = a * t.y + b * u.y + c * v.y + d * w.y;
this.z = a * t.z + b * u.z + c * v.z + d * w.z;
public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2,
double a3, Vector3D u3, double a4, Vector3D u4) {
this.x = a1 * u1.x + a2 * u2.x + a3 * u3.x + a4 * u4.x;
this.y = a1 * u1.y + a2 * u2.y + a3 * u3.y + a4 * u4.y;
this.z = a1 * u1.z + a2 * u2.z + a3 * u3.z + a4 * u4.z;
}
/** Get the abscissa of the vector.
* @return abscissa of the vector
* @see #Vector3D(double, double, double)
* @see #setX(double)
*/
public double getX() {
return x;
@ -158,7 +158,6 @@ public class Vector3D implements Serializable {
/** Get the ordinate of the vector.
* @return ordinate of the vector
* @see #Vector3D(double, double, double)
* @see #setY(double)
*/
public double getY() {
return y;
@ -167,7 +166,6 @@ public class Vector3D implements Serializable {
/** Get the height of the vector.
* @return height of the vector
* @see #Vector3D(double, double, double)
* @see #setZ(double)
*/
public double getZ() {
return z;
@ -196,38 +194,50 @@ public class Vector3D implements Serializable {
return Math.asin(z / getNorm());
}
/** Normalize a vector.
* @param v vector to normalize
* @return a new vector equal to v / ||v||
* @exception ArithmeticException if the norm of the instance is null
/** Add a vector to the instance.
* @param v vector to add
* @return a new vector
*/
public static Vector3D normalize(Vector3D v) {
double norm = v.getNorm();
if (norm == 0) {
public Vector3D add(Vector3D v) {
return new Vector3D(x + v.x, y + v.y, z + v.z);
}
/** Add a scaled vector to the instance.
* @param factor scale factor to apply to v before adding it
* @param v vector to add
* @return a new vector
*/
public Vector3D add(double factor, Vector3D v) {
return new Vector3D(x + factor * v.x, y + factor * v.y, z + factor * v.z);
}
/** Subtract a vector from the instance.
* @param v vector to subtract
* @return a new vector
*/
public Vector3D subtract(Vector3D v) {
return new Vector3D(x - v.x, y - v.y, z - v.z);
}
/** Subtract a scaled vector from the instance.
* @param factor scale factor to apply to v before subtracting it
* @param v vector to subtract
* @return a new vector
*/
public Vector3D subtract(double factor, Vector3D v) {
return new Vector3D(x - factor * v.x, y - factor * v.y, z - factor * v.z);
}
/** Normalize the instance.
* @return a new normalized vector
* @exception ArithmeticException if the norm is null
*/
public Vector3D normalize() {
double s = getNorm();
if (s == 0) {
throw new ArithmeticException("null norm");
}
double inv = 1.0 / norm;
return new Vector3D(inv * v.x, inv * v.y, inv * v.z);
}
/** Add two vectors.
* Add two vectors and return the sum as a new vector
* @param v1 first vector
* @param v2 second vector
* @return a new vector equal to v1 + v2
*/
public static Vector3D add(Vector3D v1, Vector3D v2) {
return new Vector3D(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
}
/** Subtract two vectors.
* Subtract two vectors and return the difference as a new vector
* @param v1 first vector
* @param v2 second vector
* @return a new vector equal to v1 - v2
*/
public static Vector3D subtract(Vector3D v1, Vector3D v2) {
return new Vector3D(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
return multiply(1 / s);
}
/** Get a vector orthogonal to the instance.
@ -238,7 +248,7 @@ public class Vector3D implements Serializable {
* following example shows how to build a frame having the k axis
* aligned with the known vector u :
* <pre><code>
* Vector3D k = Vector3D.normalize(u);
* Vector3D k = u.normalize();
* Vector3D i = k.orthogonal();
* Vector3D j = Vector3D.crossProduct(k, i);
* </code></pre></p>
@ -298,22 +308,19 @@ public class Vector3D implements Serializable {
}
/** Get the opposite of a vector.
* @param u vector to revert
* @return a new vector which is -u
/** Get the opposite of the instance.
* @return a new vector which is opposite to the instance
*/
public static Vector3D negate(Vector3D u) {
return new Vector3D(-u.x, -u.y, -u.z);
public Vector3D negate() {
return new Vector3D(-x, -y, -z);
}
/** Multiply a vector by a scalar
* Multiply a vectors by a scalar and return the product as a new vector
/** Multiply the instance by a scalar
* @param a scalar
* @param v vector
* @return a new vector equal to a * v
* @return a new vector
*/
public static Vector3D multiply(double a, Vector3D v) {
return new Vector3D(a * v.x, a * v.y, a * v.z);
public Vector3D multiply(double a) {
return new Vector3D(a * x, a * y, a * z);
}
/** Compute the dot-product of two vectors.
@ -345,6 +352,6 @@ public class Vector3D implements Serializable {
/** Height. */
private final double z;
private static final long serialVersionUID = 484345009325358136L;
private static final long serialVersionUID = 7318440192750283659L;
}

5
src/mantissa/src/org/spaceroots/mantissa/linalg/DiagonalMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class implements diagonal matrices of linear algebra.
* @version $Id: DiagonalMatrix.java 1705 2006-09-17 19:57:39Z luc $
@ -27,8 +25,7 @@ import java.io.Serializable;
*/
public class DiagonalMatrix
extends SquareMatrix
implements Serializable, Cloneable {
extends SquareMatrix {
/** Simple constructor.
* This constructor builds a diagonal matrix of specified order, all

5
src/mantissa/src/org/spaceroots/mantissa/linalg/GeneralMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class represents matrices of the most general type.
* <p>This class is the basic implementation of matrices to use when
@ -30,8 +28,7 @@ import java.io.Serializable;
*/
public class GeneralMatrix
extends Matrix
implements Serializable {
extends Matrix {
/** Simple constructor.
* Build a matrix with null elements.

8
src/mantissa/src/org/spaceroots/mantissa/linalg/GeneralSquareMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class implements general square matrices of linear algebra.
* @version $Id: GeneralSquareMatrix.java 1705 2006-09-17 19:57:39Z luc $
@ -27,8 +25,7 @@ import java.io.Serializable;
*/
public class GeneralSquareMatrix
extends SquareMatrix
implements Serializable, Cloneable {
extends SquareMatrix {
/** Simple constructor.
* This constructor builds a square matrix of specified order, all
@ -264,9 +261,6 @@ public class GeneralSquareMatrix
}
}
// release the memory as soon as possible
work = null;
lower = new LowerTriangularMatrix(rows, lowerData);
upper = new UpperTriangularMatrix(rows, upperData);

5
src/mantissa/src/org/spaceroots/mantissa/linalg/LowerTriangularMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class implements lower triangular matrices of linear algebra.
* @version $Id: LowerTriangularMatrix.java 1705 2006-09-17 19:57:39Z luc $
@ -27,8 +25,7 @@ import java.io.Serializable;
*/
public class LowerTriangularMatrix
extends SquareMatrix
implements Serializable, Cloneable {
extends SquareMatrix {
/** Simple constructor.
* This constructor builds a lower triangular matrix of specified order, all

2
src/mantissa/src/org/spaceroots/mantissa/linalg/Matrix.java
View File

@ -78,7 +78,7 @@ public abstract class Matrix
this.rows = rows;
this.columns = columns;
this.data = data;
this.data = (data == null) ? null : (double[]) data.clone();
}

2
src/mantissa/src/org/spaceroots/mantissa/linalg/NonNullRange.java
View File

@ -52,7 +52,7 @@ import java.io.Serializable;
*/
class NonNullRange
implements Serializable, Cloneable {
implements Serializable {
/** Index in row/column count of the first non-null element. */
public final int begin;

5
src/mantissa/src/org/spaceroots/mantissa/linalg/SquareMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class factor all services common to square matrices of linear algebra.
* <p>This class is the base class of all square matrix
@ -31,8 +29,7 @@ import java.io.Serializable;
*/
public abstract class SquareMatrix
extends Matrix
implements Serializable, Cloneable {
extends Matrix {
/** Simple constructor.
* Build a matrix with null elements.
* @param order order of the matrix

5
src/mantissa/src/org/spaceroots/mantissa/linalg/SymetricalMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class implements symetrical matrices of linear algebra.
* @version $Id: SymetricalMatrix.java 1705 2006-09-17 19:57:39Z luc $
@ -27,8 +25,7 @@ import java.io.Serializable;
*/
public class SymetricalMatrix
extends GeneralSquareMatrix
implements Serializable, Cloneable {
extends GeneralSquareMatrix {
/** Simple constructor.
* This constructor builds a symetrical matrix of specified order, all

5
src/mantissa/src/org/spaceroots/mantissa/linalg/UpperTriangularMatrix.java
View File

@ -17,8 +17,6 @@
package org.spaceroots.mantissa.linalg;
import java.io.Serializable;
/** This class implements upper triangular matrices of linear algebra.
* @version $Id: UpperTriangularMatrix.java 1705 2006-09-17 19:57:39Z luc $
@ -27,8 +25,7 @@ import java.io.Serializable;
*/
public class UpperTriangularMatrix
extends SquareMatrix
implements Serializable, Cloneable {
extends SquareMatrix {
/** Simple constructor.
* This constructor builds a upper triangular matrix of specified order, all

22
src/mantissa/src/org/spaceroots/mantissa/ode/AbstractStepInterpolator.java
View File

@ -131,15 +131,8 @@ public abstract class AbstractStepInterpolator
interpolatedTime = interpolator.interpolatedTime;
if (interpolator.currentState != null) {
int dimension = interpolator.currentState.length;
currentState = new double[dimension];
System.arraycopy(interpolator.currentState, 0, currentState, 0,
dimension);
interpolatedState = new double[dimension];
System.arraycopy(interpolator.interpolatedState, 0, interpolatedState, 0,
dimension);
currentState = (double[]) interpolator.currentState.clone();
interpolatedState = (double[]) interpolator.interpolatedState.clone();
} else {
currentState = null;
interpolatedState = null;
@ -189,7 +182,14 @@ public abstract class AbstractStepInterpolator
* @return a copy of the instance.
*/
public abstract Object clone();
public Object clone() {
try {
return super.clone();
} catch (CloneNotSupportedException cnse) {
// should never happen
return null;
}
}
/** Shift one step forward.
* Copy the current time into the previous time, hence preparing the
@ -291,7 +291,7 @@ public abstract class AbstractStepInterpolator
* @return state vector at time {@link #getInterpolatedTime}
*/
public double[] getInterpolatedState() {
return interpolatedState;
return (double[]) interpolatedState.clone();
}

2
src/mantissa/src/org/spaceroots/mantissa/ode/ClassicalRungeKuttaIntegrator.java
View File

@ -47,7 +47,7 @@ package org.spaceroots.mantissa.ode;
public class ClassicalRungeKuttaIntegrator
extends RungeKuttaIntegrator {
private static final String methodName = new String("classical Runge-Kutta");
private static final String methodName = "classical Runge-Kutta";
private static final double[] c = {
1.0 / 2.0, 1.0 / 2.0, 1.0

2
src/mantissa/src/org/spaceroots/mantissa/ode/DormandPrince54Integrator.java
View File

@ -47,7 +47,7 @@ package org.spaceroots.mantissa.ode;
public class DormandPrince54Integrator
extends RungeKuttaFehlbergIntegrator {
private static final String methodName = new String("Dormand-Prince 5(4)");
private static final String methodName = "Dormand-Prince 5(4)";
private static final double[] c = {
1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0

16
src/mantissa/src/org/spaceroots/mantissa/ode/DormandPrince54StepInterpolator.java
View File

@ -68,18 +68,10 @@ class DormandPrince54StepInterpolator
} else {
int dimension = interpolator.v1.length;
v1 = new double[dimension];
v2 = new double[dimension];
v3 = new double[dimension];
v4 = new double[dimension];
System.arraycopy(interpolator.v1, 0, v1, 0, dimension);
System.arraycopy(interpolator.v2, 0, v2, 0, dimension);
System.arraycopy(interpolator.v3, 0, v3, 0, dimension);
System.arraycopy(interpolator.v4, 0, v4, 0, dimension);
v1 = (double[]) interpolator.v1.clone();
v2 = (double[]) interpolator.v2.clone();
v3 = (double[]) interpolator.v3.clone();
v4 = (double[]) interpolator.v4.clone();
vectorsInitialized = interpolator.vectorsInitialized;
}

2
src/mantissa/src/org/spaceroots/mantissa/ode/DormandPrince853Integrator.java
View File

@ -55,7 +55,7 @@ package org.spaceroots.mantissa.ode;
public class DormandPrince853Integrator
extends RungeKuttaFehlbergIntegrator {
private static final String methodName = new String("Dormand-Prince 8 (5, 3)");
private static final String methodName = "Dormand-Prince 8 (5, 3)";
private static final double sqrt6 = Math.sqrt(6.0);

11
src/mantissa/src/org/spaceroots/mantissa/ode/DummyStepInterpolator.java
View File

@ -49,7 +49,7 @@ public class DummyStepInterpolator
* to create the step interpolators by cloning an uninitialized
* model and latter initializing the copy.
*/
protected DummyStepInterpolator() {
public DummyStepInterpolator() {
super();
}
@ -83,15 +83,6 @@ public class DummyStepInterpolator
super(interpolator);
}
/** Copy the instance.
* the copy is a deep copy: its arrays are separated from the
* original arrays of the instance
* @return a copy of the instance.
*/
public Object clone() {
return new DummyStepInterpolator(this);
}
/** Compute the state at the interpolated time.
* In this class, this method does nothing: the interpolated state
* is always the state at the end of the current step.

2
src/mantissa/src/org/spaceroots/mantissa/ode/EulerIntegrator.java
View File

@ -50,7 +50,7 @@ package org.spaceroots.mantissa.ode;
public class EulerIntegrator
extends RungeKuttaIntegrator {
private static final String methodName = new String("Euler");
private static final String methodName = "Euler";
private static final double[] c = {
};

2
src/mantissa/src/org/spaceroots/mantissa/ode/GillIntegrator.java
View File

@ -46,7 +46,7 @@ package org.spaceroots.mantissa.ode;
public class GillIntegrator
extends RungeKuttaIntegrator {
private static final String methodName = new String("Gill");
private static final String methodName = "Gill";
private static final double sqrt2 = Math.sqrt(2.0);

6
src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
View File

@ -88,7 +88,7 @@ package org.spaceroots.mantissa.ode;
public class GraggBulirschStoerIntegrator
extends AdaptiveStepsizeIntegrator {
private static final String methodName = new String("Gragg-Bulirsch-Stoer");
private static final String methodName = "Gragg-Bulirsch-Stoer";
/** Simple constructor.
* Build a Gragg-Bulirsch-Stoer integrator with the given step
@ -713,7 +713,7 @@ public class GraggBulirschStoerIntegrator
// estimate if there is a chance convergence will
// be reached on next iteration, using the
// asymptotic evolution of error
double ratio = sequence [k] * sequence[k+1]
double ratio = ((double) sequence [k] * sequence[k+1])
/ (sequence[0] * sequence[0]);
if (error > ratio * ratio) {
// we don't expect to converge on next iteration
@ -740,7 +740,7 @@ public class GraggBulirschStoerIntegrator
// estimate if there is a chance convergence will
// be reached on next iteration, using the
// asymptotic evolution of error
double ratio = sequence[k+1] / sequence[0];
double ratio = ((double) sequence[k+1]) / sequence[0];
if (error > ratio * ratio) {
// we don't expect to converge on next iteration
// we reject the step immediately

2
src/mantissa/src/org/spaceroots/mantissa/ode/HighamHall54Integrator.java
View File

@ -35,7 +35,7 @@ package org.spaceroots.mantissa.ode;
public class HighamHall54Integrator
extends RungeKuttaFehlbergIntegrator {
private static final String methodName = new String("Higham-Hall 5(4)");
private static final String methodName = "Higham-Hall 5(4)";
private static final double[] c = {
2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0

2
src/mantissa/src/org/spaceroots/mantissa/ode/MidpointIntegrator.java
View File

@ -43,7 +43,7 @@ package org.spaceroots.mantissa.ode;
public class MidpointIntegrator
extends RungeKuttaIntegrator {
private static final String methodName = new String("midpoint");
private static final String methodName = "midpoint";
private static final double[] c = {
1.0 / 2.0

3
src/mantissa/src/org/spaceroots/mantissa/ode/StepNormalizer.java
View File

@ -104,8 +104,7 @@ public class StepNormalizer
interpolator.setInterpolatedTime(lastTime);
double[] state = interpolator.getInterpolatedState();
lastState = new double[state.length];
System.arraycopy(state, 0, lastState, 0, lastState.length);
lastState = (double[]) state.clone();
// take the integration direction into account
forward = (interpolator.getCurrentTime() >= lastTime);

2
src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
View File

@ -42,6 +42,8 @@ import org.spaceroots.mantissa.roots.BrentSolver;
class SwitchState
implements ComputableFunction, ConvergenceChecker {
private static final long serialVersionUID = 6944466361876662425L;
/** Switching function. */
private SwitchingFunction function;

4
src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.ode;
import java.io.Serializable;
/** This interface represents a switching function.
*
* <p>A switching function allows to handle discrete events in
@ -45,7 +47,7 @@ package org.spaceroots.mantissa.ode;
*
*/
public interface SwitchingFunction {
public interface SwitchingFunction extends Serializable {
/** Stop indicator.
* <p>This value should be used as the return value of the {@link

2
src/mantissa/src/org/spaceroots/mantissa/ode/ThreeEighthesIntegrator.java
View File

@ -46,7 +46,7 @@ package org.spaceroots.mantissa.ode;
public class ThreeEighthesIntegrator
extends RungeKuttaIntegrator {
private static final String methodName = new String("3/8");
private static final String methodName = "3/8";
private static final double[] c = {
1.0 / 3.0, 2.0 / 3.0, 1.0

22
src/mantissa/src/org/spaceroots/mantissa/optimization/DirectSearchOptimizer.java
View File

@ -243,7 +243,7 @@ public abstract class DirectSearchOptimizer {
}
RandomVectorGenerator rvg =
new CorrelatedRandomVectorGenerator(statistics.getMean(null),
new CorrelatedRandomVectorGenerator(statistics.getMean(),
statistics.getCovarianceMatrix(null),
new UniformRandomGenerator(seed));
setMultiStart(starts, rvg);
@ -344,7 +344,7 @@ public abstract class DirectSearchOptimizer {
if (i < n) {
System.arraycopy(vertexA, i, vertex, i, n - i);
}
simplex[i] = new PointCostPair(vertex);
simplex[i] = new PointCostPair(vertex, Double.NaN);
}
}
@ -356,7 +356,7 @@ public abstract class DirectSearchOptimizer {
int n = vertices.length - 1;
simplex = new PointCostPair[n + 1];
for (int i = 0; i <= n; ++i) {
simplex[i] = new PointCostPair(vertices[i]);
simplex[i] = new PointCostPair(vertices[i], Double.NaN);
}
}
@ -369,11 +369,11 @@ public abstract class DirectSearchOptimizer {
double[] vertex = generator.nextVector();
int n = vertex.length;
simplex = new PointCostPair[n + 1];
simplex[0] = new PointCostPair(vertex);
simplex[0] = new PointCostPair(vertex, Double.NaN);
// fill up the vertex
for (int i = 1; i <= n; ++i) {
simplex[i] = new PointCostPair(generator.nextVector());
simplex[i] = new PointCostPair(generator.nextVector(), Double.NaN);
}
}
@ -428,7 +428,7 @@ public abstract class DirectSearchOptimizer {
* minimizes} has not been called
*/
public PointCostPair[] getMinima() {
return minima;
return (PointCostPair[]) minima.clone();
}
/** Minimizes a cost function.
@ -522,8 +522,8 @@ public abstract class DirectSearchOptimizer {
// evaluate the cost at all non-evaluated simplex points
for (int i = 0; i < simplex.length; ++i) {
PointCostPair pair = simplex[i];
if (! pair.isEvaluated()) {
pair.setCost(evaluateCost(pair.getPoint()));
if (Double.isNaN(pair.cost)) {
simplex[i] = new PointCostPair(pair.point, evaluateCost(pair.point));
}
}
@ -538,7 +538,7 @@ public abstract class DirectSearchOptimizer {
protected void replaceWorstPoint(PointCostPair pointCostPair) {
int n = simplex.length - 1;
for (int i = 0; i < n; ++i) {
if (simplex[i].getCost() > pointCostPair.getCost()) {
if (simplex[i].cost > pointCostPair.cost) {
PointCostPair tmp = simplex[i];
simplex[i] = pointCostPair;
pointCostPair = tmp;
@ -555,8 +555,8 @@ public abstract class DirectSearchOptimizer {
} else if (o2 == null) {
return -1;
} else {
double cost1 = ((PointCostPair) o1).getCost();
double cost2 = ((PointCostPair) o2).getCost();
double cost1 = ((PointCostPair) o1).cost;
double cost2 = ((PointCostPair) o2).cost;
return (cost1 < cost2) ? -1 : ((o1 == o2) ? 0 : +1);
}
}

10
src/mantissa/src/org/spaceroots/mantissa/optimization/MultiDirectional.java
View File

@ -54,7 +54,7 @@ public class MultiDirectional
// save the original vertex
PointCostPair[] original = simplex;
double originalCost = original[0].getCost();
double originalCost = original[0].cost;
// perform a reflection step
double reflectedCost = evaluateNewSimplex(original, 1.0);
@ -93,24 +93,24 @@ public class MultiDirectional
private double evaluateNewSimplex(PointCostPair[] original, double coeff)
throws CostException {
double[] xSmallest = original[0].getPoint();
double[] xSmallest = original[0].point;
int n = xSmallest.length;
// create the linearly transformed simplex
simplex = new PointCostPair[n + 1];
simplex[0] = original[0];
for (int i = 1; i <= n; ++i) {
double[] xOriginal = original[i].getPoint();
double[] xOriginal = original[i].point;
double[] xTransformed = new double[n];
for (int j = 0; j < n; ++j) {
xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]);
}
simplex[i] = new PointCostPair(xTransformed);
simplex[i] = new PointCostPair(xTransformed, Double.NaN);
}
// evaluate it
evaluateSimplex();
return simplex[0].getCost();
return simplex[0].cost;
}

37
src/mantissa/src/org/spaceroots/mantissa/optimization/NelderMead.java
View File

@ -61,16 +61,16 @@ public class NelderMead
int n = simplex.length - 1;
// interesting costs
double smallest = simplex[0].getCost();
double secondLargest = simplex[n-1].getCost();
double largest = simplex[n].getCost();
double[] xLargest = simplex[n].getPoint();
double smallest = simplex[0].cost;
double secondLargest = simplex[n-1].cost;
double largest = simplex[n].cost;
double[] xLargest = simplex[n].point;
// compute the centroid of the best vertices
// (dismissing the worst point at index n)
double[] centroid = new double[n];
for (int i = 0; i < n; ++i) {
double[] x = simplex[i].getPoint();
double[] x = simplex[i].point;
for (int j = 0; j < n; ++j) {
centroid[j] += x[j];
}
@ -90,9 +90,7 @@ public class NelderMead
if ((smallest <= costR) && (costR < secondLargest)) {
// accept the reflected point
PointCostPair r = new PointCostPair(xR);
r.setCost(costR);
replaceWorstPoint(r);
replaceWorstPoint(new PointCostPair(xR, costR));
} else if (costR < smallest) {
@ -105,14 +103,10 @@ public class NelderMead
if (costE < costR) {
// accept the expansion point
PointCostPair e = new PointCostPair(xE);
e.setCost(costE);
replaceWorstPoint(e);
replaceWorstPoint(new PointCostPair(xE, costE));
} else {
// accept the reflected point
PointCostPair r = new PointCostPair(xR);
r.setCost(costR);
replaceWorstPoint(r);
replaceWorstPoint(new PointCostPair(xR, costR));
}
} else {
@ -128,9 +122,7 @@ public class NelderMead
if (costC <= costR) {
// accept the contraction point
PointCostPair c = new PointCostPair(xC);
c.setCost(costC);
replaceWorstPoint(c);
replaceWorstPoint(new PointCostPair(xC, costC));
return;
}
@ -145,23 +137,20 @@ public class NelderMead
if (costC < largest) {
// accept the contraction point
PointCostPair c = new PointCostPair(xC);
c.setCost(costC);
replaceWorstPoint(c);
replaceWorstPoint(new PointCostPair(xC, costC));
return;
}
}
// perform a shrink
double[] xSmallest = simplex[0].getPoint();
double[] xSmallest = simplex[0].point;
for (int i = 1; i < simplex.length; ++i) {
PointCostPair pair = simplex[i];
double[] x = pair.getPoint();
double[] x = simplex[i].point;
for (int j = 0; j < n; ++j) {
x[j] = xSmallest[j] + sigma * (x[j] - xSmallest[j]);
}
pair.setCost(Double.NaN);
simplex[i] = new PointCostPair(x, Double.NaN);
}
evaluateSimplex();

62
src/mantissa/src/org/spaceroots/mantissa/optimization/PointCostPair.java
View File

@ -18,72 +18,26 @@
package org.spaceroots.mantissa.optimization;
/** This class holds a point and its associated cost.
* <p>A cost/point pair is not evaluated at build time. Its associated
* cost set to <code>Double.NaN</code> until it is evaluated.</p>
* <p>This is a simple immutable container.</p>
* @author Luc Maisonobe
* @version $Id: PointCostPair.java 1705 2006-09-17 19:57:39Z luc $
* @version $Id: PointCostPair.java 1709 2006-12-03 21:16:50Z luc $
* @see CostFunction
*/
public class PointCostPair {
/** Build a point/cost pair with non-evaluated cost.
/** Build a point/cost pair.
* @param point point coordinates
* @param cost point cost
*/
public PointCostPair(double[] point) {
this.point = point;
cost = Double.NaN;
}
/** Reset the point coordinates.
* <p>Resetting the points coordinates automatically reset the cost
* to non-evaluated</p>
* @param point new point coordinates
* @return old point coordinates (this can be re-used to put the
* coordinates of another point without re-allocating an array)
*/
public double[] setPoint(double[] point) {
double[] oldPoint = this.point;
this.point = point;
cost = Double.NaN;
return oldPoint;
}
/** Get the point coordinates.
* @return point coordinates
*/
public double[] getPoint() {
return point;
}
/** Set the cost.
* @param cost cost to store in the instance (can be
* <code>Double.NaN</code> to reset the instance to non-evaluated)
*/
public void setCost(double cost) {
public PointCostPair(double[] point, double cost) {
this.point = (double[]) point.clone();
this.cost = cost;
}
/** Get the cost.
* @return cost associated to the point (or <code>Double.NaN</code>
* if the instance is not evaluated)
*/
public double getCost() {
return cost;
}
/** Check if the cost has been evaluated.
* <p>The cost is considered to be non-evaluated if it is
* <code>Double.isNaN(pair.getCost())</code> would return true</p>
* @return true if the cost has been evaluated
*/
public boolean isEvaluated() {
return ! Double.isNaN(cost);
}
/** Point coordinates. */
private double[] point;
public final double[] point;
/** Cost associated to the point. */
private double cost;
public final double cost;
}

2
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/EnhancedSimpsonIntegrator.java
View File

@ -45,7 +45,7 @@ public class EnhancedSimpsonIntegrator
try {
while (true) {
sum = sampler.nextSamplePoint().getY();
sum = sampler.nextSamplePoint().y;
}
} catch(ExhaustedSampleException e) {
}

22
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/EnhancedSimpsonIntegratorSampler.java
View File

@ -85,15 +85,15 @@ public class EnhancedSimpsonIntegratorSampler
try {
next = iter.nextSamplePoint();
double h1 = current.getX() - previous.getX();
double h2 = next.getX() - current.getX();
double h1 = current.x - previous.x;
double h2 = next.x - current.x;
double cP = (h1 + h2) * (2 * h1 - h2) / (6 * h1);
double cC = (h1 + h2) * (h1 + h2) * (h1 + h2) / (6 * h1 * h2);
double cN = (h1 + h2) * (2 * h2 - h1) / (6 * h2);
double[] pY = previous.getY();
double[] cY = current.getY();
double[] nY = next.getY();
double[] pY = previous.y;
double[] cY = current.y;
double[] nY = next.y;
for (int i = 0; i < sum.length; ++i) {
sum [i] += cP * pY[i] + cC * cY[i] + cN * nY[i];
}
@ -101,18 +101,16 @@ public class EnhancedSimpsonIntegratorSampler
} catch(ExhaustedSampleException e) {
// we have an incomplete step at the end of the sample
// we use a trapezoid scheme for this last step
double halfDx = 0.5 * (current.getX() - previous.getX());
double[] pY = previous.getY();
double[] cY = current.getY();
double halfDx = 0.5 * (current.x - previous.x);
double[] pY = previous.y;
double[] cY = current.y;
for (int i = 0; i < sum.length; ++i) {
sum [i] += halfDx * (pY[i] + cY[i]);
}
return new VectorialValuedPair(current.getX(), sum);
return new VectorialValuedPair(current.x, sum);
}
double[] values = new double[sum.length];
System.arraycopy(sum, 0, values, 0, sum.length);
return new VectorialValuedPair(next.getX(), values);
return new VectorialValuedPair(next.x, (double[]) sum.clone());
}

2
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/RiemannIntegrator.java
View File

@ -51,7 +51,7 @@ public class RiemannIntegrator
try {
while (true) {
sum = sampler.nextSamplePoint().getY();
sum = sampler.nextSamplePoint().y;
}
} catch(ExhaustedSampleException e) {
}

8
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/RiemannIntegratorSampler.java
View File

@ -86,15 +86,13 @@ public class RiemannIntegratorSampler
// performs one step of a Riemann scheme
VectorialValuedPair previous = current;
current = iter.nextSamplePoint();
double step = (current.getX() - previous.getX());
double[] pY = previous.getY();
double step = (current.x - previous.x);
double[] pY = previous.y;
for (int i = 0; i < sum.length; ++i) {
sum[i] += step * pY[i];
}
double[] values = new double[sum.length];
System.arraycopy(sum, 0, values, 0, sum.length);
return new VectorialValuedPair (current.getX(), values);
return new VectorialValuedPair (current.x, (double[]) sum.clone());
}

2
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/TrapezoidIntegrator.java
View File

@ -42,7 +42,7 @@ public class TrapezoidIntegrator
try {
while (true) {
sum = sampler.nextSamplePoint().getY();
sum = sampler.nextSamplePoint().y;
}
} catch(ExhaustedSampleException e) {
}

10
src/mantissa/src/org/spaceroots/mantissa/quadrature/vectorial/TrapezoidIntegratorSampler.java
View File

@ -87,16 +87,14 @@ public class TrapezoidIntegratorSampler
VectorialValuedPair previous = current;
current = iter.nextSamplePoint();
double halfDx = 0.5 * (current.getX() - previous.getX());
double[] pY = previous.getY();
double[] cY = current.getY();
double halfDx = 0.5 * (current.x - previous.x);
double[] pY = previous.y;
double[] cY = current.y;
for (int i = 0; i < sum.length; ++i) {
sum[i] += halfDx * (pY[i] + cY[i]);
}
double[] values = new double[sum.length];
System.arraycopy(sum, 0, values, 0, sum.length);
return new VectorialValuedPair (current.getX(), values);
return new VectorialValuedPair (current.x, (double[]) sum.clone());
}

16
src/mantissa/src/org/spaceroots/mantissa/random/CorrelatedRandomVectorGenerator.java
View File

@ -81,13 +81,12 @@ public class CorrelatedRandomVectorGenerator
});
throw new IllegalArgumentException(message);
}
this.mean = mean;
this.mean = (double[]) mean.clone();
factorize(covariance);
this.generator = generator;
normalized = new double[rank];
correlated = new double[order];
}
@ -114,7 +113,6 @@ public class CorrelatedRandomVectorGenerator
this.generator = generator;
normalized = new double[rank];
correlated = new double[order];
}
@ -240,10 +238,8 @@ public class CorrelatedRandomVectorGenerator
}
/** Generate a correlated random vector.
* @return a random vector as an array of double. The generator
* <em>will</em> reuse the same array for each call, in order to
* save the allocation time, so the user should keep a copy by
* himself if he needs so.
* @return a random vector as an array of double. The returned array
* is created at each call, the caller can do what it wants with it.
*/
public double[] nextVector() {
@ -253,6 +249,7 @@ public class CorrelatedRandomVectorGenerator
}
// compute correlated vector
double[] correlated = new double[mean.length];
for (int i = 0; i < correlated.length; ++i) {
correlated[i] = mean[i];
for (int j = 0; j < rank; ++j) {
@ -279,9 +276,6 @@ public class CorrelatedRandomVectorGenerator
/** Storage for the normalized vector. */
private double[] normalized;
/** Storage for the random vector. */
private double[] correlated;
private static final long serialVersionUID = -4754497552287369719L;
private static final long serialVersionUID = -88563624902398453L;
}

8
src/mantissa/src/org/spaceroots/mantissa/random/GaussianRandomGenerator.java
View File

@ -33,9 +33,6 @@ import java.util.Random;
public class GaussianRandomGenerator
implements NormalizedRandomGenerator {
/** Underlying generator. */
Random generator;
/** Create a new generator.
* The seed of the generator is related to the current time.
*/
@ -64,4 +61,9 @@ public class GaussianRandomGenerator
return generator.nextGaussian();
}
/** Underlying generator. */
private Random generator;
private static final long serialVersionUID = 5504568059866195697L;
}

4
src/mantissa/src/org/spaceroots/mantissa/random/NormalizedRandomGenerator.java
View File

@ -17,6 +17,8 @@
package org.spaceroots.mantissa.random;
import java.io.Serializable;
/** This interface represent a normalized random generator for
* scalars.
* Normalized generator should provide null mean and unit standard
@ -24,7 +26,7 @@ package org.spaceroots.mantissa.random;
* @version $Id: NormalizedRandomGenerator.java 1705 2006-09-17 19:57:39Z luc $
* @author L. Maisonobe
*/
public interface NormalizedRandomGenerator {
public interface NormalizedRandomGenerator extends Serializable {
/** Generate a random scalar with null mean and unit standard deviation.
* <p>This method does <strong>not</strong> specify the shape of the

18
src/mantissa/src/org/spaceroots/mantissa/random/UncorrelatedRandomVectorGenerator.java
View File

@ -46,11 +46,10 @@ public class UncorrelatedRandomVectorGenerator
if (mean.length != standardDeviation.length) {
throw new IllegalArgumentException("dimension mismatch");
}
this.mean = mean;
this.standardDeviation = standardDeviation;
this.mean = (double[]) mean.clone();
this.standardDeviation = (double[]) standardDeviation.clone();
this.generator = generator;
random = new double[mean.length];
}
@ -72,7 +71,6 @@ public class UncorrelatedRandomVectorGenerator
}
this.generator = generator;
random = new double[dimension];
}
@ -84,13 +82,12 @@ public class UncorrelatedRandomVectorGenerator
}
/** Generate a correlated random vector.
* @return a random vector as an array of double. The generator
* <em>will</em> reuse the same array for each call, in order to
* save the allocation time, so the user should keep a copy by
* himself if he needs so.
* @return a random vector as an array of double. The returned array
* is created at each call, the caller can do what it wants with it.
*/
public double[] nextVector() {
double[] random = new double[mean.length];
for (int i = 0; i < random.length; ++i) {
random[i] = mean[i] + standardDeviation[i] * generator.nextDouble();
}
@ -108,9 +105,6 @@ public class UncorrelatedRandomVectorGenerator
/** Underlying scalar generator. */
NormalizedRandomGenerator generator;
/** Storage for the random vector. */
private double[] random;
private static final long serialVersionUID = -3323293740860311151L;
private static final long serialVersionUID = -9094322067568302961L;
}

16
src/mantissa/src/org/spaceroots/mantissa/random/UniformRandomGenerator.java
View File

@ -34,13 +34,6 @@ import java.util.Random;
public class UniformRandomGenerator
implements NormalizedRandomGenerator {
private static final double SQRT3 = Math.sqrt(3.0);
private static final double TWOSQRT3 = 2.0 * Math.sqrt(3.0);
/** Underlying generator. */
Random generator;
/** Create a new generator.
* The seed of the generator is related to the current time.
*/
@ -71,4 +64,13 @@ public class UniformRandomGenerator
return TWOSQRT3 * generator.nextDouble() - SQRT3;
}
/** Underlying generator. */
private Random generator;
private static final double SQRT3 = Math.sqrt(3.0);
private static final double TWOSQRT3 = 2.0 * Math.sqrt(3.0);
private static final long serialVersionUID = -6913329325753217654L;
}

90
src/mantissa/src/org/spaceroots/mantissa/random/VectorialSampleStatistics.java
View File

@ -19,8 +19,6 @@ package org.spaceroots.mantissa.random;
import org.spaceroots.mantissa.linalg.SymetricalMatrix;
import java.util.Arrays;
/** This class compute basic statistics on a scalar sample.
* @version $Id: VectorialSampleStatistics.java 1705 2006-09-17 19:57:39Z luc $
* @author L. Maisonobe
@ -65,16 +63,6 @@ public class VectorialSampleStatistics {
sum2 = null;
}
/** Allocate all the arrays. */
private void allocate() {
min = new double[dimension];
minIndices = new int[dimension];
max = new double[dimension];
maxIndices = new int[dimension];
sum = new double[dimension];
sum2 = new double[dimension * (dimension + 1) / 2];
}
/** Add one point to the instance.
* @param x value of the sample point
* @exception IllegalArgumentException if there is a dimension
@ -85,14 +73,13 @@ public class VectorialSampleStatistics {
if (n == 0) {
dimension = x.length;
allocate();
Arrays.fill(minIndices, 0);
Arrays.fill(maxIndices, 0);
System.arraycopy(x, 0, min, 0, dimension);
System.arraycopy(x, 0, max, 0, dimension);
System.arraycopy(x, 0, sum, 0, dimension);
dimension = x.length;
minIndices = new int[dimension];
maxIndices = new int[dimension];
min = (double[]) x.clone();
max = (double[]) x.clone();
sum = (double[]) x.clone();
sum2 = new double[dimension * (dimension + 1) / 2];
int k = 0;
for (int i = 0; i < dimension; ++i) {
@ -153,14 +140,12 @@ public class VectorialSampleStatistics {
if (n == 0) {
dimension = s.dimension;
allocate();
System.arraycopy(s.min, 0, min, 0, dimension);
System.arraycopy(s.minIndices, 0, minIndices, 0, dimension);
System.arraycopy(s.max, 0, max, 0, dimension);
System.arraycopy(s.maxIndices, 0, maxIndices, 0, dimension);
System.arraycopy(s.sum, 0, sum, 0, dimension);
System.arraycopy(s.sum2, 0, sum2, 0, sum2.length);
min = (double[]) s.min.clone();
minIndices = (int[]) s.minIndices.clone();
max = (double[]) s.max.clone();
maxIndices = (int[]) s.maxIndices.clone();
sum = (double[]) s.sum.clone();
sum2 = (double[]) s.sum2.clone();
} else {
int k = 0;
@ -203,26 +188,22 @@ public class VectorialSampleStatistics {
* of the sample at which the minimum was encountered can be
* retrieved with the {@link #getMinIndices getMinIndices}
* method.</p>
* @return minimal value in the sample (the array is a reference to
* an internal array that changes each time something is added to
* the instance, the caller should neither change it nor rely on its
* value in the long term)
* @return minimal value in the sample (a new array is created
* at each call, the caller may do what it wants to with it)
* @see #getMinIndices
*/
public double[] getMin() {
return min;
return (double[]) min.clone();
}
/** Get the indices at which the minimal value occurred in the sample.
* @return a vector reporting at which occurrence each component of
* the sample reached its minimal value (the array is a reference to
* an internal array that changes each time something is added to
* the instance, the caller should neither change it nor rely on its
* value in the long term)
* the sample reached its minimal value (a new array is created
* at each call, the caller may do what it wants to with it)
* @see #getMin
*/
public int[] getMinIndices() {
return minIndices;
return (int[]) minIndices.clone();
}
/** Get the maximal value in the sample.
@ -232,43 +213,34 @@ public class VectorialSampleStatistics {
* of the sample at which the maximum was encountered can be
* retrieved with the {@link #getMaxIndices getMaxIndices}
* method.</p>
* @return maximal value in the sample (the array is a reference to
* an internal array that changes each time something is added to
* the instance, the caller should neither change it nor rely on its
* value in the long term)
* @return maximal value in the sample (a new array is created
* at each call, the caller may do what it wants to with it)
* @see #getMaxIndices
*/
public double[] getMax() {
return max;
return (double[]) max.clone();
}
/** Get the indices at which the maximal value occurred in the sample.
* @return a vector reporting at which occurrence each component of
* the sample reached its maximal value (the array is a reference to
* an internal array that changes each time something is added to
* the instance, the caller should neither change it nor rely on its
* value in the long term)
* the sample reached its maximal value (a new array is created
* at each call, the caller may do what it wants to with it)
* @see #getMax
*/
public int[] getMaxIndices() {
return maxIndices;
return (int[]) maxIndices.clone();
}
/** Get the mean value of the sample.
* @param mean placeholder where to store the array, if null a new
* array will be allocated
* @return mean value of the sample or null if the sample is empty
* and hence the dimension of the vectors is still unknown
* (reference to mean if it was non-null, reference to a new array
* otherwise)
* @return mean value of the sample or an empty array
* if the sample is empty (a new array is created
* at each call, the caller may do what it wants to with it)
*/
public double[] getMean(double[] mean) {
public double[] getMean() {
if (n == 0) {
return null;
}
if (mean == null) {
mean = new double[dimension];
return new double[0];
}
double[] mean = new double[dimension];
for (int i = 0; i < dimension; ++i) {
mean[i] = sum[i] / n;
}

8
src/mantissa/src/org/spaceroots/mantissa/utilities/ArrayMapper.java
View File

@ -77,14 +77,14 @@ public class ArrayMapper {
}
/** Get the internal data array.
* @return internal data array
/** Get the data array.
* @return copy of the data array
*/
public double[] getInternalDataArray() {
public double[] getDataArray() {
if (internalData == null) {
internalData = new double [size];
}
return internalData;
return (double[]) internalData.clone();
}
/** Map data from the internal array to the domain objects.

14
src/mantissa/src/org/spaceroots/mantissa/utilities/MappableArray.java
View File

@ -43,24 +43,16 @@ public class MappableArray
/** Simple constructor.
* Build a mappable array from an existing array
* @param array array to use
* @param doReallocate true if a new array should be allocated and
* initialized using the other argument, false if the instance
* should reference the existing array throughout its lifetime
*/
public MappableArray(double[] array, boolean doReallocate) {
if (doReallocate) {
internalArray = new double[array.length];
System.arraycopy(array, 0, internalArray, 0, array.length);
} else {
internalArray = array;
}
public MappableArray(double[] array) {
internalArray = (double[]) array.clone();
}
/** Get the array stored in the instance.
* @return array stored in the instance
*/
public double[] getArray () {
return internalArray;
return (double[]) internalArray.clone();
}
/** Get the dimension of the internal array.

17
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/ChebyshevTest.java
View File

@ -26,7 +26,7 @@ public class ChebyshevTest
super(name);
}
public void aatestOne() {
public void testOne() {
assertTrue(new Chebyshev().isOne());
}
@ -44,16 +44,16 @@ public class ChebyshevTest
}
public void aatestBounds() {
public void testBounds() {
for (int k = 0; k < 12; ++k) {
Chebyshev Tk = new Chebyshev(k);
OrthogonalPolynomial Tk = new Chebyshev(k);
for (double x = -1.0; x <= 1.0; x += 0.02) {
assertTrue(Math.abs(Tk.valueAt(x)) < (1.0 + 1.0e-12));
}
}
}
public void aatestDifferentials() {
public void testDifferentials() {
for (int k = 0; k < 12; ++k) {
Polynomial.Rational Tk0 = new Chebyshev(k);
@ -64,12 +64,11 @@ public class ChebyshevTest
Polynomial.Rational g1 = new Polynomial.Rational(-1l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(-1l, 0l, 1l);
Polynomial.Rational Tk0g0 = Polynomial.Rational.multiply(Tk0, g0);
Polynomial.Rational Tk1g1 = Polynomial.Rational.multiply(Tk1, g1);
Polynomial.Rational Tk2g2 = Polynomial.Rational.multiply(Tk2, g2);
Polynomial.Rational Tk0g0 = Tk0.multiply(g0);
Polynomial.Rational Tk1g1 = Tk1.multiply(g1);
Polynomial.Rational Tk2g2 = Tk2.multiply(g2);
Polynomial.Rational d =
Polynomial.Rational.add(Tk0g0, Polynomial.Rational.add(Tk1g1, Tk2g2));
Polynomial.Rational d = Tk0g0.add(Tk1g1.add(Tk2g2));
assertTrue(d.isZero());
}

9
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/HermiteTest.java
View File

@ -55,12 +55,11 @@ public class HermiteTest
Polynomial.Rational g1 = new Polynomial.Rational(-2l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(1l);
Polynomial.Rational Hk0g0 = Polynomial.Rational.multiply(Hk0, g0);
Polynomial.Rational Hk1g1 = Polynomial.Rational.multiply(Hk1, g1);
Polynomial.Rational Hk2g2 = Polynomial.Rational.multiply(Hk2, g2);
Polynomial.Rational Hk0g0 = Hk0.multiply(g0);
Polynomial.Rational Hk1g1 = Hk1.multiply(g1);
Polynomial.Rational Hk2g2 = Hk2.multiply(g2);
Polynomial.Rational d =
Polynomial.Rational.add(Hk0g0, Polynomial.Rational.add(Hk1g1, Hk2g2));
Polynomial.Rational d = Hk0g0.add(Hk1g1.add(Hk2g2));
assertTrue(d.isZero());
}

11
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/LaguerreTest.java
View File

@ -61,19 +61,18 @@ public class LaguerreTest
Polynomial.Rational g1 = new Polynomial.Rational(-1l, 1l);
Polynomial.Rational g2 = new Polynomial.Rational(1l, 0l);
Polynomial.Rational Lk0g0 = Polynomial.Rational.multiply(Lk0, g0);
Polynomial.Rational Lk1g1 = Polynomial.Rational.multiply(Lk1, g1);
Polynomial.Rational Lk2g2 = Polynomial.Rational.multiply(Lk2, g2);
Polynomial.Rational Lk0g0 = Lk0.multiply(g0);
Polynomial.Rational Lk1g1 = Lk1.multiply(g1);
Polynomial.Rational Lk2g2 = Lk2.multiply(g2);
Polynomial.Rational d =
Polynomial.Rational.add(Lk0g0, Polynomial.Rational.add(Lk1g1, Lk2g2));
Polynomial.Rational d = Lk0g0.add(Lk1g1.add(Lk2g2));
assertTrue(d.isZero());
}
}
public void checkLaguerre(Laguerre p, long denominator, String reference) {
assertTrue(Laguerre.multiply(p, denominator).toString().equals(reference));
assertTrue(p.multiply(denominator).toString().equals(reference));
}
public static Test suite() {

11
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/LegendreTest.java
View File

@ -55,12 +55,11 @@ public class LegendreTest
Polynomial.Rational g1 = new Polynomial.Rational(-2l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(-1l, 0l, 1l);
Polynomial.Rational Pk0g0 = Polynomial.Rational.multiply(Pk0, g0);
Polynomial.Rational Pk1g1 = Polynomial.Rational.multiply(Pk1, g1);
Polynomial.Rational Pk2g2 = Polynomial.Rational.multiply(Pk2, g2);
Polynomial.Rational Pk0g0 = Pk0.multiply(g0);
Polynomial.Rational Pk1g1 = Pk1.multiply(g1);
Polynomial.Rational Pk2g2 = Pk2.multiply(g2);
Polynomial.Rational d =
Polynomial.Rational.add(Pk0g0, Polynomial.Rational.add(Pk1g1, Pk2g2));
Polynomial.Rational d = Pk0g0.add(Pk1g1.add(Pk2g2));
assertTrue(d.isZero());
}
@ -92,7 +91,7 @@ public class LegendreTest
}
public void checkLegendre(Legendre p, long denominator, String reference) {
assertTrue(Legendre.multiply(p, denominator).toString().equals(reference));
assertTrue(p.multiply(denominator).toString().equals(reference));
}
public static Test suite() {

28
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/PolynomialDoubleTest.java
View File

@ -49,20 +49,14 @@ public class PolynomialDoubleTest
public void testConversion() {
Polynomial.Rational r = new Polynomial.Rational(1l, 3l, -5l);
r.multiplySelf(new RationalNumber(1l, 2l));
r = (Polynomial.Rational) r.multiply(new RationalNumber(1l, 2l));
Polynomial.Double p = new Polynomial.Double(r);
checkPolynomial(p, "-2.5 + 1.5 x + 0.5 x^2");
}
public void testString() {
Polynomial.Double p = new Polynomial.Double(1.0, 3.0, -5.0);
checkPolynomial(p, "-5.0 + 3.0 x + x^2");
p.setUnknownName("alpha");
checkPolynomial(p, "-5.0 + 3.0 alpha + alpha^2");
p.setUnknownName(null);
checkPolynomial(p, "-5.0 + 3.0 x + x^2");
checkPolynomial(new Polynomial.Double(3.0, -2.0, 0.0),
"-2.0 x + 3.0 x^2");
checkPolynomial(new Polynomial.Double(3.0, -2.0, 1.0),
@ -75,22 +69,20 @@ public class PolynomialDoubleTest
"1.0 + 3.0 x^2");
checkPolynomial(new Polynomial.Double(0.0),
"0");
}
public void testAddition() {
Polynomial.Double p1 = new Polynomial.Double(1.0, -2.0);
Polynomial.Double p2 = new Polynomial.Double(0.0, -1.0, 2.0);
assertTrue(Polynomial.Double.add(p1, p2).isZero());
assertTrue(p1.add(p2).isZero());
p2 = new Polynomial.Double(p1);
p2.addToSelf(p2);
p2 = p1.add(p1);
checkPolynomial(p2, "-4.0 + 2.0 x");
p1 = new Polynomial.Double(2.0, -4.0, 1.0);
p2 = new Polynomial.Double(-2.0, 3.0, -1.0);
p1.addToSelf(p2);
p1 = p1.add(p2);
assertEquals(1, p1.getDegree());
checkPolynomial(p1, "-x");
@ -99,15 +91,15 @@ public class PolynomialDoubleTest
public void testSubtraction() {
Polynomial.Double p1 = new Polynomial.Double(1.0, -2.0);
assertTrue(Polynomial.Double.subtract(p1, p1).isZero());
assertTrue(p1.subtract(p1).isZero());
Polynomial.Double p2 = new Polynomial.Double(6.0, -2.0);
p2.subtractFromSelf(p1);
p2 = p2.subtract(p1);
checkPolynomial(p2, "5.0 x");
p1 = new Polynomial.Double(2.0, -4.0, 1.0);
p2 = new Polynomial.Double(2.0, 3.0, -1.0);
p1.subtractFromSelf(p2);
p1 = p1.subtract(p2);
assertEquals(1, p1.getDegree());
checkPolynomial(p1, "2.0 - 7.0 x");
@ -117,12 +109,12 @@ public class PolynomialDoubleTest
Polynomial.Double p1 = new Polynomial.Double(2.0, -3.0);
Polynomial.Double p2 = new Polynomial.Double(1.0, 2.0, 3.0);
checkPolynomial(Polynomial.Double.multiply(p1, p2), "-9.0 + x^2 + 2.0 x^3");
checkPolynomial(p1.multiply(p2), "-9.0 + x^2 + 2.0 x^3");
p1 = new Polynomial.Double(1.0, 0.0);
p2 = new Polynomial.Double(p1);
p2 = p1;
for (int i = 2; i < 10; ++i) {
p2.multiplySelf(p1);
p2 = p2.multiply(p1);
checkPolynomial(p2, "x^" + i);
}

100
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/PolynomialFractionTest.java
View File

@ -55,18 +55,18 @@ public class PolynomialFractionTest
public void testInvert() {
PolynomialFraction f = new PolynomialFraction(2l, 4l);
f.invertSelf();
f= f.invert();
checkValue(f, "2");
f.invertSelf();
f = f.invert();
checkValue(f, "1/2");
f = new PolynomialFraction(120l);
f.invertSelf();
f = f.invert();
checkValue(f, "1/120");
f = new PolynomialFraction(0l, 4l);
try {
f.invertSelf();
f = f.invert();
fail("an exception should have been thrown");
} catch (ArithmeticException e) {
} catch (Exception e) {
@ -74,7 +74,7 @@ public class PolynomialFractionTest
}
f = new PolynomialFraction(307692l, 999999l);
PolynomialFraction fInverse = PolynomialFraction.invert(f);
PolynomialFraction fInverse = f.invert();
checkValue(fInverse, "13/4");
checkValue(f, "4/13");
@ -83,23 +83,18 @@ public class PolynomialFractionTest
public void testAddition() {
PolynomialFraction f1 = new PolynomialFraction(4l, 6l);
f1.addToSelf(f1);
f1 = f1.add(f1);
checkValue(f1, "4/3");
checkValue(PolynomialFraction.add(new PolynomialFraction(17l, 3l),
new PolynomialFraction(-17l, 3l)),
checkValue(new PolynomialFraction(17l, 3l).add(new PolynomialFraction(-17l, 3l)),
"0");
checkValue(PolynomialFraction.add(new PolynomialFraction(2l, 3l),
new PolynomialFraction(3l, 4l)),
checkValue(new PolynomialFraction(2l, 3l).add(new PolynomialFraction(3l, 4l)),
"17/12");
checkValue(PolynomialFraction.add(new PolynomialFraction(1l, 6l),
new PolynomialFraction(2l, 6l)),
checkValue(new PolynomialFraction(1l, 6l).add(new PolynomialFraction(2l, 6l)),
"1/2");
checkValue(PolynomialFraction.add(new PolynomialFraction(4l, 5l),
new PolynomialFraction(-3l, 4l)),
checkValue(new PolynomialFraction(4l, 5l).add(new PolynomialFraction(-3l, 4l)),
"1/20");
checkValue(PolynomialFraction.add(new PolynomialFraction(-3l, 4l),
new PolynomialFraction(4l, 5l)),
checkValue(new PolynomialFraction(-3l, 4l).add(new PolynomialFraction(4l, 5l)),
"1/20");
}
@ -107,44 +102,32 @@ public class PolynomialFractionTest
public void testSubtraction() {
PolynomialFraction f1 = new PolynomialFraction(4l, 6l);
f1.subtractFromSelf(f1);
checkValue(f1, "0");
checkValue(f1.subtract(f1), "0");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(7l, 3l),
new PolynomialFraction(-7l, 3l)),
checkValue(new PolynomialFraction(7l, 3l).subtract(new PolynomialFraction(-7l, 3l)),
"14/3");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(3l, 4l),
new PolynomialFraction(2l, 3l)),
checkValue(new PolynomialFraction(3l, 4l).subtract(new PolynomialFraction(2l, 3l)),
"1/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(3l, 4l),
new PolynomialFraction(-2l, 3l)),
checkValue(new PolynomialFraction(3l, 4l).subtract(new PolynomialFraction(-2l, 3l)),
"17/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(-3l, 4l),
new PolynomialFraction(2l, 3l)),
checkValue(new PolynomialFraction(-3l, 4l).subtract(new PolynomialFraction(2l, 3l)),
"-17/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(-3l, 4l),
new PolynomialFraction(-2l, 3l)),
checkValue(new PolynomialFraction(-3l, 4l).subtract(new PolynomialFraction(-2l, 3l)),
"-1/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(2l, 3l),
new PolynomialFraction(3l, 4l)),
checkValue(new PolynomialFraction(2l, 3l).subtract(new PolynomialFraction(3l, 4l)),
"-1/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(-2l, 3l),
new PolynomialFraction(3l, 4l)),
checkValue(new PolynomialFraction(-2l, 3l).subtract(new PolynomialFraction(3l, 4l)),
"-17/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(2l, 3l),
new PolynomialFraction(-3l, 4l)),
checkValue(new PolynomialFraction(2l, 3l).subtract(new PolynomialFraction(-3l, 4l)),
"17/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(-2l, 3l),
new PolynomialFraction(-3l, 4l)),
checkValue(new PolynomialFraction(-2l, 3l).subtract(new PolynomialFraction(-3l, 4l)),
"1/12");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(1l, 6l),
new PolynomialFraction(2l, 6l)),
checkValue(new PolynomialFraction(1l, 6l).subtract(new PolynomialFraction(2l, 6l)),
"-1/6");
checkValue(PolynomialFraction.subtract(new PolynomialFraction(1l, 2l),
new PolynomialFraction(1l, 6l)),
checkValue(new PolynomialFraction(1l, 2l).subtract(new PolynomialFraction(1l, 6l)),
"1/3");
}
@ -152,23 +135,17 @@ public class PolynomialFractionTest
public void testMultiplication() {
PolynomialFraction f = new PolynomialFraction(2l, 3l);
f.multiplySelf(new PolynomialFraction(9l,4l));
checkValue(f, "3/2");
checkValue(f.multiply(new PolynomialFraction(9l,4l)), "3/2");
checkValue(PolynomialFraction.multiply(new PolynomialFraction(1l, 2l),
new PolynomialFraction(0l)),
checkValue(new PolynomialFraction(1l, 2l).multiply(new PolynomialFraction(0l)),
"0");
checkValue(PolynomialFraction.multiply(new PolynomialFraction(4l, 15l),
new PolynomialFraction(-5l, 2l)),
checkValue(new PolynomialFraction(4l, 15l).multiply(new PolynomialFraction(-5l, 2l)),
"-2/3");
checkValue(PolynomialFraction.multiply(new PolynomialFraction(-4l, 15l),
new PolynomialFraction(5l, 2l)),
checkValue(new PolynomialFraction(-4l, 15l).multiply(new PolynomialFraction(5l, 2l)),
"-2/3");
checkValue(PolynomialFraction.multiply(new PolynomialFraction(4l, 15l),
new PolynomialFraction(5l, 2l)),
checkValue(new PolynomialFraction(4l, 15l).multiply(new PolynomialFraction(5l, 2l)),
"2/3");
checkValue(PolynomialFraction.multiply(new PolynomialFraction(-4l, 15l),
new PolynomialFraction(-5l, 2l)),
checkValue(new PolynomialFraction(-4l, 15l).multiply(new PolynomialFraction(-5l, 2l)),
"2/3");
}
@ -176,29 +153,24 @@ public class PolynomialFractionTest
public void testDivision() {
PolynomialFraction f = new PolynomialFraction(2l, 3l);
f.divideSelf(new PolynomialFraction(4l,9l));
checkValue(f, "3/2");
;
checkValue(f.divide(new PolynomialFraction(4l,9l)), "3/2");
try {
PolynomialFraction.divide(new PolynomialFraction(1l, 2l),
new PolynomialFraction(0l));
new PolynomialFraction(1l, 2l).divide(new PolynomialFraction(0l));
fail("an exception should have been thrown");
} catch (ArithmeticException e) {
} catch (Exception e) {
fail("wrong exception caught");
}
checkValue(PolynomialFraction.divide(new PolynomialFraction(4l, 15l),
new PolynomialFraction(-2l, 5l)),
checkValue(new PolynomialFraction(4l, 15l).divide(new PolynomialFraction(-2l, 5l)),
"-2/3");
checkValue(PolynomialFraction.divide(new PolynomialFraction(-4l, 15l),
new PolynomialFraction(2l, 5l)),
checkValue(new PolynomialFraction(-4l, 15l).divide(new PolynomialFraction(2l, 5l)),
"-2/3");
checkValue(PolynomialFraction.divide(new PolynomialFraction(4l, 15l),
new PolynomialFraction(2l, 5l)),
checkValue(new PolynomialFraction(4l, 15l).divide(new PolynomialFraction(2l, 5l)),
"2/3");
checkValue(PolynomialFraction.divide(new PolynomialFraction(-4l, 15l),
new PolynomialFraction(-2l, 5l)),
checkValue(new PolynomialFraction(-4l, 15l).divide(new PolynomialFraction(-2l, 5l)),
"2/3");
}

25
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/PolynomialRationalTest.java
View File

@ -57,10 +57,6 @@ public class PolynomialRationalTest
Polynomial.Rational p = new Polynomial.Rational(1l, 3l, -5l);
checkPolynomial(p, "-5 + 3 x + x^2");
p.setUnknownName("alpha");
checkPolynomial(p, "-5 + 3 alpha + alpha^2");
p.setUnknownName(null);
checkPolynomial(p, "-5 + 3 x + x^2");
checkPolynomial(new Polynomial.Rational(3l, -2l, 0l), "-2 x + 3 x^2");
checkPolynomial(new Polynomial.Rational(3l, -2l, 1l), "1 - 2 x + 3 x^2");
@ -75,15 +71,14 @@ public class PolynomialRationalTest
Polynomial.Rational p1 = new Polynomial.Rational(1l, -2l);
Polynomial.Rational p2 = new Polynomial.Rational(0l, -1l, 2l);
assertTrue(Polynomial.Rational.add(p1, p2).isZero());
assertTrue(p1.add(p2).isZero());
p2 = new Polynomial.Rational(p1);
p2.addToSelf(p2);
p2 = p1.add(p1);
checkPolynomial(p2, "-4 + 2 x");
p1 = new Polynomial.Rational(2l, -4l, 1l);
p2 = new Polynomial.Rational(-2l, 3l, -1l);
p1.addToSelf(p2);
p1 = p1.add(p2);
assertEquals(1, p1.getDegree());
checkPolynomial(p1, "-x");
@ -92,15 +87,15 @@ public class PolynomialRationalTest
public void testSubtraction() {
Polynomial.Rational p1 = new Polynomial.Rational(1l, -2l);
assertTrue(Polynomial.Rational.subtract(p1, p1).isZero());
assertTrue(p1.subtract(p1).isZero());
Polynomial.Rational p2 = new Polynomial.Rational(6l, -2l);
p2.subtractFromSelf(p1);
p2 = p2.subtract(p1);
checkPolynomial(p2, "5 x");
p1 = new Polynomial.Rational(2l, -4l, 1l);
p2 = new Polynomial.Rational(2l, 3l, -1l);
p1.subtractFromSelf(p2);
p1 = p1.subtract(p2);
assertEquals(1, p1.getDegree());
checkPolynomial(p1, "2 - 7 x");
@ -110,12 +105,12 @@ public class PolynomialRationalTest
Polynomial.Rational p1 = new Polynomial.Rational(2l, -3l);
Polynomial.Rational p2 = new Polynomial.Rational(1l, 2l, 3l);
checkPolynomial(Polynomial.Rational.multiply(p1, p2), "-9 + x^2 + 2 x^3");
checkPolynomial(p1.multiply(p2), "-9 + x^2 + 2 x^3");
p1 = new Polynomial.Rational(1l, 0l);
p2 = new Polynomial.Rational(p1);
p2 = p1;
for (int i = 2; i < 10; ++i) {
p2.multiplySelf(p1);
p2 = p2.multiply(p1);
checkPolynomial(p2, "x^" + i);
}
@ -128,7 +123,7 @@ public class PolynomialRationalTest
checkPolynomial(p, "3/4 - 1/6 x + 2/5 x^2");
BigInteger lcm = p.getDenominatorsLCM();
assertEquals(BigInteger.valueOf(60l), lcm);
p.multiplySelf(lcm);
p = (Polynomial.Rational) p.multiply(lcm);
checkPolynomial(p, "45 - 10 x + 24 x^2");
}

98
src/mantissa/tests-src/org/spaceroots/mantissa/algebra/RationalNumberTest.java
View File

@ -54,19 +54,17 @@ public class RationalNumberTest
public void testInvert() {
RationalNumber f = new RationalNumber(2l, 4l);
f.invertSelf();
RationalNumber f = new RationalNumber(2l, 4l).invert();
checkValue(f, "2");
f.invertSelf();
f = f.invert();
checkValue(f, "1/2");
f = new RationalNumber(120l);
f.invertSelf();
f = new RationalNumber(120l).invert();
checkValue(f, "1/120");
f = new RationalNumber(0l, 4l);
try {
f.invertSelf();
f.invert();
fail("an exception should have been thrown");
} catch (ArithmeticException e) {
} catch (Exception e) {
@ -74,7 +72,7 @@ public class RationalNumberTest
}
f = new RationalNumber(307692l, 999999l);
RationalNumber fInverse = RationalNumber.invert(f);
RationalNumber fInverse = f.invert();
checkValue(fInverse, "13/4");
checkValue(f, "4/13");
@ -83,23 +81,18 @@ public class RationalNumberTest
public void testAddition() {
RationalNumber f1 = new RationalNumber(4l, 6l);
f1.addToSelf(f1);
f1 = f1.add(f1);
checkValue(f1, "4/3");
checkValue(RationalNumber.add(new RationalNumber(17l, 3l),
new RationalNumber(-17l, 3l)),
checkValue(new RationalNumber(17l, 3l).add(new RationalNumber(-17l, 3l)),
"0");
checkValue(RationalNumber.add(new RationalNumber(2l, 3l),
new RationalNumber(3l, 4l)),
checkValue(new RationalNumber(2l, 3l).add(new RationalNumber(3l, 4l)),
"17/12");
checkValue(RationalNumber.add(new RationalNumber(1l, 6l),
new RationalNumber(2l, 6l)),
checkValue(new RationalNumber(1l, 6l).add(new RationalNumber(2l, 6l)),
"1/2");
checkValue(RationalNumber.add(new RationalNumber(4l, 5l),
new RationalNumber(-3l, 4l)),
checkValue(new RationalNumber(4l, 5l).add(new RationalNumber(-3l, 4l)),
"1/20");
checkValue(RationalNumber.add(new RationalNumber(-3l, 4l),
new RationalNumber(4l, 5l)),
checkValue(new RationalNumber(-3l, 4l).add(new RationalNumber(4l, 5l)),
"1/20");
}
@ -107,44 +100,33 @@ public class RationalNumberTest
public void testSubtraction() {
RationalNumber f1 = new RationalNumber(4l, 6l);
f1.subtractFromSelf(f1);
f1 = f1.subtract(f1);
checkValue(f1, "0");
checkValue(RationalNumber.subtract(new RationalNumber(7l, 3l),
new RationalNumber(-7l, 3l)),
checkValue(new RationalNumber(7l, 3l).subtract(new RationalNumber(-7l, 3l)),
"14/3");
checkValue(RationalNumber.subtract(new RationalNumber(3l, 4l),
new RationalNumber(2l, 3l)),
checkValue(new RationalNumber(3l, 4l).subtract(new RationalNumber(2l, 3l)),
"1/12");
checkValue(RationalNumber.subtract(new RationalNumber(3l, 4l),
new RationalNumber(-2l, 3l)),
checkValue(new RationalNumber(3l, 4l).subtract(new RationalNumber(-2l, 3l)),
"17/12");
checkValue(RationalNumber.subtract(new RationalNumber(-3l, 4l),
new RationalNumber(2l, 3l)),
checkValue(new RationalNumber(-3l, 4l).subtract(new RationalNumber(2l, 3l)),
"-17/12");
checkValue(RationalNumber.subtract(new RationalNumber(-3l, 4l),
new RationalNumber(-2l, 3l)),
checkValue(new RationalNumber(-3l, 4l).subtract(new RationalNumber(-2l, 3l)),
"-1/12");
checkValue(RationalNumber.subtract(new RationalNumber(2l, 3l),
new RationalNumber(3l, 4l)),
checkValue(new RationalNumber(2l, 3l).subtract(new RationalNumber(3l, 4l)),
"-1/12");
checkValue(RationalNumber.subtract(new RationalNumber(-2l, 3l),
new RationalNumber(3l, 4l)),
checkValue(new RationalNumber(-2l, 3l).subtract(new RationalNumber(3l, 4l)),
"-17/12");
checkValue(RationalNumber.subtract(new RationalNumber(2l, 3l),
new RationalNumber(-3l, 4l)),
checkValue(new RationalNumber(2l, 3l).subtract(new RationalNumber(-3l, 4l)),
"17/12");
checkValue(RationalNumber.subtract(new RationalNumber(-2l, 3l),
new RationalNumber(-3l, 4l)),
checkValue(new RationalNumber(-2l, 3l).subtract(new RationalNumber(-3l, 4l)),
"1/12");
checkValue(RationalNumber.subtract(new RationalNumber(1l, 6l),
new RationalNumber(2l, 6l)),
checkValue(new RationalNumber(1l, 6l).subtract(new RationalNumber(2l, 6l)),
"-1/6");
checkValue(RationalNumber.subtract(new RationalNumber(1l, 2l),
new RationalNumber(1l, 6l)),
checkValue(new RationalNumber(1l, 2l).subtract(new RationalNumber(1l, 6l)),
"1/3");
}
@ -152,23 +134,18 @@ public class RationalNumberTest
public void testMultiplication() {
RationalNumber f = new RationalNumber(2l, 3l);
f.multiplySelf(new RationalNumber(9l,4l));
f = f.multiply(new RationalNumber(9l,4l));
checkValue(f, "3/2");
checkValue(RationalNumber.multiply(new RationalNumber(1l, 2l),
new RationalNumber(0l)),
checkValue(new RationalNumber(1l, 2l).multiply(new RationalNumber(0l)),
"0");
checkValue(RationalNumber.multiply(new RationalNumber(4l, 15l),
new RationalNumber(-5l, 2l)),
checkValue(new RationalNumber(4l, 15l).multiply(new RationalNumber(-5l, 2l)),
"-2/3");
checkValue(RationalNumber.multiply(new RationalNumber(-4l, 15l),
new RationalNumber(5l, 2l)),
checkValue(new RationalNumber(-4l, 15l).multiply(new RationalNumber(5l, 2l)),
"-2/3");
checkValue(RationalNumber.multiply(new RationalNumber(4l, 15l),
new RationalNumber(5l, 2l)),
checkValue(new RationalNumber(4l, 15l).multiply(new RationalNumber(5l, 2l)),
"2/3");
checkValue(RationalNumber.multiply(new RationalNumber(-4l, 15l),
new RationalNumber(-5l, 2l)),
checkValue(new RationalNumber(-4l, 15l).multiply(new RationalNumber(-5l, 2l)),
"2/3");
}
@ -176,29 +153,24 @@ public class RationalNumberTest
public void testDivision() {
RationalNumber f = new RationalNumber(2l, 3l);
f.divideSelf(new RationalNumber(4l,9l));
f = f.divide(new RationalNumber(4l,9l));
checkValue(f, "3/2");
try {
RationalNumber.divide(new RationalNumber(1l, 2l),
new RationalNumber(0l));
new RationalNumber(1l, 2l).divide(new RationalNumber(0l));
fail("an exception should have been thrown");
} catch (ArithmeticException e) {
} catch (Exception e) {
fail("wrong exception caught");
}
checkValue(RationalNumber.divide(new RationalNumber(4l, 15l),
new RationalNumber(-2l, 5l)),
checkValue(new RationalNumber(4l, 15l).divide(new RationalNumber(-2l, 5l)),
"-2/3");
checkValue(RationalNumber.divide(new RationalNumber(-4l, 15l),
new RationalNumber(2l, 5l)),
checkValue(new RationalNumber(-4l, 15l).divide(new RationalNumber(2l, 5l)),
"-2/3");
checkValue(RationalNumber.divide(new RationalNumber(4l, 15l),
new RationalNumber(2l, 5l)),
checkValue(new RationalNumber(4l, 15l).divide(new RationalNumber(2l, 5l)),
"2/3");
checkValue(RationalNumber.divide(new RationalNumber(-4l, 15l),
new RationalNumber(-2l, 5l)),
checkValue(new RationalNumber(-4l, 15l).divide(new RationalNumber(-2l, 5l)),
"2/3");
}

1
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/AllTests.java
View File

@ -30,6 +30,7 @@ public class AllTests {
suite.addTest(WeightedMeasurementTest.suite());
suite.addTest(GaussNewtonEstimatorTest.suite());
suite.addTest(LevenbergMarquardtEstimatorTest.suite());
suite.addTest(MinpackTest.suite());
return suite;

10
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/GaussNewtonEstimatorTest.java
View File

@ -106,7 +106,7 @@ public class GaussNewtonEstimatorTest
for (int i = 0; i < gridSize; ++i) {
for (int j = 0; j < gridSize; ++j) {
String name = new Integer(k).toString();
String name = Integer.toString(k);
perfectPars[2 * k] = new EstimatedParameter("x" + name, i);
perfectPars[2 * k + 1] = new EstimatedParameter("y" + name, j);
++k;
@ -194,7 +194,7 @@ public class GaussNewtonEstimatorTest
}
private class Distance extends WeightedMeasurement {
private static class Distance extends WeightedMeasurement {
public Distance(double weight, double measuredValue,
EstimatedParameter x1, EstimatedParameter y1,
@ -243,15 +243,15 @@ public class GaussNewtonEstimatorTest
}
public WeightedMeasurement[] getMeasurements() {
return measurements;
return (WeightedMeasurement[]) measurements.clone();
}
public EstimatedParameter[] getUnboundParameters() {
return unboundPars;
return (EstimatedParameter[]) unboundPars.clone();
}
public EstimatedParameter[] getAllParameters() {
return randomizedPars;
return (EstimatedParameter[]) randomizedPars.clone();
}
private EstimatedParameter[] perfectPars;

1415
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/LevenbergMarquardtEstimatorTest.java
View File

File diff suppressed because it is too large Load Diff

20
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/WeightedMeasurementTest.java
View File

@ -24,16 +24,18 @@ public class WeightedMeasurementTest
public WeightedMeasurementTest(String name) {
super(name);
p1 = null;
p2 = null;
}
public void testConstruction() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1);
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getWeight(), 3.0);
checkValue(m.getMeasuredValue(), theoretical() + 0.1);
}
public void testIgnored() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1);
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
assertTrue(!m.isIgnored());
m.setIgnored(true);
assertTrue(m.isIgnored());
@ -42,7 +44,7 @@ public class WeightedMeasurementTest
}
public void testTheory() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1);
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getTheoreticalValue(), theoretical());
checkValue(m.getResidual(), 0.1);
@ -92,21 +94,25 @@ public class WeightedMeasurementTest
}
}
private class MyMeasurement
private static class MyMeasurement
extends WeightedMeasurement {
public MyMeasurement(double weight, double measuredValue) {
public MyMeasurement(double weight, double measuredValue,
WeightedMeasurementTest testInstance) {
super(weight, measuredValue);
this.testInstance = testInstance;
}
public double getTheoreticalValue() {
return theoretical();
return testInstance.theoretical();
}
public double getPartial(EstimatedParameter p) {
return partial(p);
return testInstance.partial(p);
}
private transient WeightedMeasurementTest testInstance;
private static final long serialVersionUID = -246712922500792332L;
}

10
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/AbstractCurveFitterTest.java
View File

@ -21,6 +21,7 @@ import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.estimation.WeightedMeasurement;
public class AbstractCurveFitterTest
@ -28,6 +29,7 @@ public class AbstractCurveFitterTest
public AbstractCurveFitterTest(String name) {
super(name);
fitter = null;
}
public void testAlreadySorted() {
@ -74,15 +76,15 @@ public class AbstractCurveFitterTest
fitter = new DummyFitter();
}
public void tearOff() {
public void tearDown() {
fitter = null;
}
private class DummyFitter
private static class DummyFitter
extends AbstractCurveFitter {
public DummyFitter() {
super(10, 10, 0.0, 0.0, 0.0);
super(10, new LevenbergMarquardtEstimator());
}
public double valueAt(double x) {
@ -97,7 +99,7 @@ public class AbstractCurveFitterTest
sortMeasurements();
}
private static final long serialVersionUID = -5453139487565082528L;
private static final long serialVersionUID = 4016396219767783678L;
}

33
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java
View File

@ -21,6 +21,7 @@ import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.estimation.WeightedMeasurement;
public class HarmonicFitterTest
@ -34,8 +35,8 @@ public class HarmonicFitterTest
throws EstimationException {
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter = new HarmonicFitter(20, 1.0e-7,
1.0e-10, 1.0e-10);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
for (double x = 0.0; x < 1.3; x += 0.01) {
fitter.addWeightedPair(1.0, x, f.valueAt(x));
}
@ -45,12 +46,12 @@ public class HarmonicFitterTest
HarmonicFunction fitted = new HarmonicFunction(coeffs[0],
coeffs[1],
coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-12);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 1.0e-12);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.0e-12);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-13);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 1.0e-13);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.0e-13);
for (double x = -1.0; x < 1.0; x += 0.01) {
assertTrue(Math.abs(f.valueAt(x) - fitted.valueAt(x)) < 1.0e-12);
assertTrue(Math.abs(f.valueAt(x) - fitted.valueAt(x)) < 1.0e-13);
}
}
@ -60,8 +61,8 @@ public class HarmonicFitterTest
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter = new HarmonicFitter(20, 1.0e-7,
1.0e-10, 1.0e-10);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addWeightedPair(1.0, x,
f.valueAt(x) + 0.01 * randomizer.nextGaussian());
@ -70,9 +71,9 @@ public class HarmonicFitterTest
double[] coeffs = fitter.fit();
new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-3);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 3.5e-3);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 2.0e-2);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 2.7e-3);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.3e-2);
WeightedMeasurement[] measurements = fitter.getMeasurements();
for (int i = 0; i < measurements.length; ++i) {
@ -88,8 +89,8 @@ public class HarmonicFitterTest
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter = new HarmonicFitter(100, 1.0e-7,
1.0e-10, 1.0e-10);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
// build a regularly spaced array of measurements
int size = 100;
@ -120,9 +121,9 @@ public class HarmonicFitterTest
double[] coeffs = fitter.fit();
new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-3);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 3.5e-3);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 2.0e-2);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.5e-2);
WeightedMeasurement[] measurements = fitter.getMeasurements();
for (int i = 0; i < measurements.length; ++i) {
@ -142,7 +143,7 @@ public class HarmonicFitterTest
return a - twoPi * Math.floor((a + Math.PI - ref) / twoPi);
}
private class HarmonicFunction {
private static class HarmonicFunction {
public HarmonicFunction(double a, double omega, double phi) {
this.a = a;
this.omega = omega;

42
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java
View File

@ -21,6 +21,9 @@ import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
public class PolynomialFitterTest
extends TestCase {
@ -38,9 +41,8 @@ public class PolynomialFitterTest
p.initCoeff (i, randomizer.nextGaussian());
}
PolynomialFitter fitter = new PolynomialFitter(degree,
10, 1.0e-7,
1.0e-10, 1.0e-10);
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
for (int i = 0; i <= degree; ++i) {
fitter.addWeightedPair(1.0, i, p.valueAt(i));
}
@ -66,9 +68,8 @@ public class PolynomialFitterTest
p.initCoeff(i, randomizer.nextGaussian());
}
PolynomialFitter fitter = new PolynomialFitter(degree,
10, 1.0e-7,
1.0e-10, 1.0e-10);
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addWeightedPair(1.0, x,
p.valueAt(x) + 0.1 * randomizer.nextGaussian());
@ -85,18 +86,28 @@ public class PolynomialFitterTest
}
public void testUnsolvableProblem()
throws EstimationException {
public void testRedundantSolvable() {
// Levenberg-Marquardt should handle redundant information gracefully
checkUnsolvableProblem(new LevenbergMarquardtEstimator(), true);
}
public void testRedundantUnsolvable() {
// Gauss-Newton should not be able to solve redundant information
checkUnsolvableProblem(new GaussNewtonEstimator(10, 1.0e-7, 1.0e-7,
1.0e-10),
false);
}
private void checkUnsolvableProblem(Estimator estimator,
boolean solvable) {
Random randomizer = new Random(1248788532l);
for (int degree = 0; degree < 10; ++degree) {
Polynom p = new Polynom(degree);
for (int i = 1; i <= degree; ++i) {
for (int i = 0; i <= degree; ++i) {
p.initCoeff(i, randomizer.nextGaussian());
}
PolynomialFitter fitter = new PolynomialFitter(degree,
10, 1.0e-7,
1.0e-10, 1.0e-10);
PolynomialFitter fitter = new PolynomialFitter(degree, estimator);
// reusing the same point over and over again does not bring
// information, the problem cannot be solved in this case for
@ -106,13 +117,12 @@ public class PolynomialFitterTest
fitter.addWeightedPair(1.0, 0.0, p.valueAt(0.0));
}
boolean gotIt = false;
try {
fitter.fit();
assertTrue(solvable || (degree == 0));
} catch(EstimationException e) {
gotIt = true;
assertTrue((! solvable) && (degree > 0));
}
assertTrue((degree == 0 && ! gotIt) || (degree > 0 && gotIt));
}
@ -122,7 +132,7 @@ public class PolynomialFitterTest
return new TestSuite(PolynomialFitterTest.class);
}
private class Polynom {
private static class Polynom {
public Polynom(int degree) {
coeffs = new double[degree + 1];

42
src/mantissa/tests-src/org/spaceroots/mantissa/functions/scalar/BasicSampledFunctionIteratorTest.java
View File

@ -71,24 +71,7 @@ public class BasicSampledFunctionIteratorTest
throws ExhaustedSampleException, FunctionException {
BasicSampledFunctionIterator iter =
new BasicSampledFunctionIterator(new SampledFunction() {
private boolean fireException = false;
public int size() {
return 2;
}
public ScalarValuedPair samplePointAt(int i)
throws FunctionException {
if (fireException) {
throw new FunctionException("boom");
}
fireException = true;
return new ScalarValuedPair(0.0, 0.0);
}
});
new BasicSampledFunctionIterator(new ExceptionGeneratingFunction());
boolean exceptionOccurred = false;
try {
@ -112,9 +95,10 @@ public class BasicSampledFunctionIteratorTest
return new TestSuite(BasicSampledFunctionIteratorTest.class);
}
private class Function
private static class Function
implements SampledFunction {
private static final long serialVersionUID = -5071329620086891960L;
private double begin;
private double step;
private int n;
@ -142,4 +126,24 @@ public class BasicSampledFunctionIteratorTest
}
}
private static class ExceptionGeneratingFunction
implements SampledFunction {
private static final long serialVersionUID = 1417147976215668305L;
private boolean fireException = false;
public int size() {
return 2;
}
public ScalarValuedPair samplePointAt(int i)
throws FunctionException {
if (fireException) {
throw new FunctionException("boom");
}
fireException = true;
return new ScalarValuedPair(0.0, 0.0);
}
}
}

25
src/mantissa/tests-src/org/spaceroots/mantissa/functions/scalar/ComputableFunctionSamplerTest.java
View File

@ -131,15 +131,7 @@ public class ComputableFunctionSamplerTest
public void testUnderlyingException() {
ComputableFunctionSampler sampler =
new ComputableFunctionSampler(new ComputableFunction() {
public double valueAt(double x)
throws FunctionException {
if (x < 0.5) {
return -x;
}
throw new FunctionException("upper half range exception");
}
},
new ComputableFunctionSampler(new ExceptionGeneratingFunction(),
0.0, 0.1, 11);
boolean exceptionOccurred = false;
@ -164,9 +156,10 @@ public class ComputableFunctionSamplerTest
return new TestSuite(ComputableFunctionSamplerTest.class);
}
private class Function
private static class Function
implements ComputableFunction {
private static final long serialVersionUID = -7173012970400285826L;
private double min;
private double max;
@ -188,4 +181,16 @@ public class ComputableFunctionSamplerTest
}
private static class ExceptionGeneratingFunction
implements ComputableFunction {
private static final long serialVersionUID = 7853080602731012102L;
public double valueAt(double x)
throws FunctionException {
if (x < 0.5) {
return -x;
}
throw new FunctionException("upper half range exception");
}
}
}

55
src/mantissa/tests-src/org/spaceroots/mantissa/functions/vectorial/BasicSampledFunctionIteratorTest.java
View File

@ -38,9 +38,9 @@ public class BasicSampledFunctionIteratorTest
for (int i = 0; i < 10; ++i) {
assertTrue(iter.hasNext());
VectorialValuedPair pair = iter.nextSamplePoint();
assertTrue(Math.abs(pair.getX() - 0.1 * i) < 1.0e-10);
assertTrue(Math.abs(pair.getY()[0] + 0.1 * i) < 1.0e-10);
assertTrue(Math.abs(pair.getY()[1] + 0.2 * i) < 1.0e-10);
assertTrue(Math.abs(pair.x - 0.1 * i) < 1.0e-10);
assertTrue(Math.abs(pair.y[0] + 0.1 * i) < 1.0e-10);
assertTrue(Math.abs(pair.y[1] + 0.2 * i) < 1.0e-10);
}
}
@ -72,27 +72,7 @@ public class BasicSampledFunctionIteratorTest
throws ExhaustedSampleException, FunctionException {
BasicSampledFunctionIterator iter =
new BasicSampledFunctionIterator(new SampledFunction() {
private boolean fireException = false;
public int size() {
return 2;
}
public int getDimension() {
return 2;
}
public VectorialValuedPair samplePointAt(int i)
throws FunctionException {
if (fireException) {
throw new FunctionException("boom");
}
fireException = true;
return new VectorialValuedPair(0.0, null);
}
});
new BasicSampledFunctionIterator(new ExceptionGeneratingFunction());
boolean exceptionOccurred = false;
try {
@ -116,9 +96,10 @@ public class BasicSampledFunctionIteratorTest
return new TestSuite(BasicSampledFunctionIteratorTest.class);
}
private class Function
private static class Function
implements SampledFunction {
private static final long serialVersionUID = -6049535144225908344L;
private double begin;
private double step;
private int n;
@ -154,4 +135,28 @@ public class BasicSampledFunctionIteratorTest
}
}
private static class ExceptionGeneratingFunction
implements SampledFunction {
private static final long serialVersionUID = 3750401068561053681L;
private boolean fireException = false;
public int size() {
return 2;
}
public int getDimension() {
return 2;
}
public VectorialValuedPair samplePointAt(int i)
throws FunctionException {
if (fireException) {
throw new FunctionException("boom");
}
fireException = true;
return new VectorialValuedPair(0.0, new double[] { 0, 1 });
}
}
}

118
src/mantissa/tests-src/org/spaceroots/mantissa/functions/vectorial/ComputableFunctionSamplerTest.java
View File

@ -36,15 +36,15 @@ public class ComputableFunctionSamplerTest
assertTrue(sampler.size() == 11);
assertTrue(sampler.getDimension() == 2);
assertTrue(Math.abs(sampler.samplePointAt(0).getX() - 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[0] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[1] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getX() - 0.495) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[0] + 0.495) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[1] + 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getX() - 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getY()[0] + 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getY()[1] + 1.980) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).x - 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[0] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[1] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).x - 0.495) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[0] + 0.495) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[1] + 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).x - 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).y[0] + 0.990) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).y[1] + 1.980) < 1.0e-10);
}
@ -59,15 +59,15 @@ public class ComputableFunctionSamplerTest
assertTrue(sampler.size() == 11);
assertTrue(sampler.getDimension() == 2);
assertTrue(Math.abs(sampler.samplePointAt(0).getX() - 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[0] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[1] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getX() - 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[0] + 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[1] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getX() - 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getY()[0] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).getY()[1] + 2.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).x - 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[0] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[1] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).x - 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[0] + 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[1] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).x - 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).y[0] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(10).y[1] + 2.0) < 1.0e-10);
}
@ -83,15 +83,15 @@ public class ComputableFunctionSamplerTest
assertTrue(sampler.size() == 12);
assertTrue(sampler.getDimension() == 2);
assertTrue(Math.abs(sampler.samplePointAt(0).getX() - 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[0] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[1] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getX() - 0.415) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[0] + 0.415) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).getY()[1] + 0.830) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).getX() - 0.913) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).getY()[0] + 0.913) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).getY()[1] + 1.826) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).x - 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[0] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[1] + 0.000) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).x - 0.415) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[0] + 0.415) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(5).y[1] + 0.830) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).x - 0.913) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).y[0] + 0.913) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(11).y[1] + 1.826) < 1.0e-10);
}
@ -107,15 +107,15 @@ public class ComputableFunctionSamplerTest
assertTrue(sampler.size() == 13);
assertTrue(sampler.getDimension() == 2);
assertTrue(Math.abs(sampler.samplePointAt(0).getX() - 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[0] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).getY()[1] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).getX() - 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).getY()[0] + 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).getY()[1] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).getX() - 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).getY()[0] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).getY()[1] + 2.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).x - 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[0] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(0).y[1] + 0.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).x - 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).y[0] + 0.5) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(6).y[1] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).x - 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).y[0] + 1.0) < 1.0e-10);
assertTrue(Math.abs(sampler.samplePointAt(12).y[1] + 2.0) < 1.0e-10);
}
@ -146,24 +146,7 @@ public class ComputableFunctionSamplerTest
public void testUnderlyingException() {
ComputableFunctionSampler sampler =
new ComputableFunctionSampler(new ComputableFunction() {
public int getDimension() {
return 2;
}
public double[] valueAt(double x)
throws FunctionException {
if (x < 0.5) {
double[] res = new double[2];
res[0] = -x;
res[1] = -2.0 * x;
return res;
}
throw new FunctionException("upper half range exception");
}
},
new ComputableFunctionSampler(new ExceptionGeneratingFunction(),
0.0, 0.1, 11);
boolean exceptionOccurred = false;
@ -188,7 +171,7 @@ public class ComputableFunctionSamplerTest
return new TestSuite(ComputableFunctionSamplerTest.class);
}
private class Function
private static class Function
implements ComputableFunction {
private double min;
private double max;
@ -216,6 +199,31 @@ public class ComputableFunctionSamplerTest
return values;
}
private static final long serialVersionUID = -1859103913610458563L;
}
private static class ExceptionGeneratingFunction
implements ComputableFunction {
public int getDimension() {
return 2;
}
public double[] valueAt(double x)
throws FunctionException {
if (x < 0.5) {
double[] res = new double[2];
res[0] = -x;
res[1] = -2.0 * x;
return res;
}
throw new FunctionException("upper half range exception");
}
private static final long serialVersionUID = 2849780376767626912L;
}
}

21
src/mantissa/tests-src/org/spaceroots/mantissa/functions/vectorial/VectorialValuedPairTest.java
View File

@ -31,9 +31,9 @@ public class VectorialValuedPairTest
tab[0] = -8.4;
tab[1] = -3.2;
VectorialValuedPair pair = new VectorialValuedPair(1.2, tab);
assertTrue(Math.abs(pair.getX() - 1.2) < 1.0e-10);
assertTrue(Math.abs(pair.getY()[0] + 8.4) < 1.0e-10);
assertTrue(Math.abs(pair.getY()[1] + 3.2) < 1.0e-10);
assertTrue(Math.abs(pair.x - 1.2) < 1.0e-10);
assertTrue(Math.abs(pair.y[0] + 8.4) < 1.0e-10);
assertTrue(Math.abs(pair.y[1] + 3.2) < 1.0e-10);
}
public void testCopyConstructor() {
@ -41,13 +41,14 @@ public class VectorialValuedPairTest
tab[0] = -8.4;
tab[1] = -3.2;
VectorialValuedPair pair1 = new VectorialValuedPair(1.2, tab);
VectorialValuedPair pair2 = new VectorialValuedPair(pair1);
assertTrue(Math.abs(pair2.getX() - pair1.getX()) < 1.0e-10);
assertTrue(Math.abs(pair2.getY()[0] - pair1.getY()[0]) < 1.0e-10);
assertTrue(Math.abs(pair2.getY()[1] - pair1.getY()[1]) < 1.0e-10);
assertTrue(Math.abs(pair2.getX() - 1.2) < 1.0e-10);
assertTrue(Math.abs(pair2.getY()[0] + 8.4) < 1.0e-10);
assertTrue(Math.abs(pair2.getY()[1] + 3.2) < 1.0e-10);
VectorialValuedPair pair2 = new VectorialValuedPair(pair1.x,
pair1.y);
assertTrue(Math.abs(pair2.x - pair1.x) < 1.0e-10);
assertTrue(Math.abs(pair2.y[0] - pair1.y[0]) < 1.0e-10);
assertTrue(Math.abs(pair2.y[1] - pair1.y[1]) < 1.0e-10);
assertTrue(Math.abs(pair2.x - 1.2) < 1.0e-10);
assertTrue(Math.abs(pair2.y[0] + 8.4) < 1.0e-10);
assertTrue(Math.abs(pair2.y[1] + 3.2) < 1.0e-10);
}
public static Test suite() {

18
src/mantissa/tests-src/org/spaceroots/mantissa/geometry/RotationTest.java
View File

@ -59,7 +59,7 @@ public class RotationTest
checkAngle(r.getAngle(), 2 * Math.PI / 3);
try {
r = new Rotation(new Vector3D(0, 0, 0), 2 * Math.PI / 3);
new Rotation(new Vector3D(0, 0, 0), 2 * Math.PI / 3);
fail("an exception should have been thrown");
} catch (ArithmeticException e) {
} catch (Exception e) {
@ -81,10 +81,9 @@ public class RotationTest
Vector3D u = new Vector3D(3, 2, 1);
Vector3D v = new Vector3D(-4, 2, 2);
Rotation r = new Rotation(u, v);
checkVector(r.applyTo(Vector3D.multiply(v.getNorm(), u)),
Vector3D.multiply(u.getNorm(), v));
checkVector(r.applyTo(u.multiply(v.getNorm())), v.multiply(u.getNorm()));
checkAngle(new Rotation(u, Vector3D.negate(u)).getAngle(), Math.PI);
checkAngle(new Rotation(u, u.negate()).getAngle(), Math.PI);
}
@ -96,14 +95,14 @@ public class RotationTest
Vector3D v2 = new Vector3D(-2, 0, 2);
Rotation r = new Rotation(u1, u2, v1, v2);
checkVector(r.applyTo(Vector3D.plusI), Vector3D.plusK);
checkVector(r.applyTo(Vector3D.plusJ), Vector3D.negate(Vector3D.plusI));
checkVector(r.applyTo(Vector3D.plusJ), Vector3D.minusI);
r = new Rotation(u1, u2, Vector3D.negate(u1), Vector3D.negate(u2));
r = new Rotation(u1, u2, u1.negate(), u2.negate());
Vector3D axis = r.getAxis();
if (Vector3D.dotProduct(axis, Vector3D.plusK) > 0) {
checkVector(axis, Vector3D.plusK);
} else {
checkVector(axis, Vector3D.negate(Vector3D.plusK));
checkVector(axis, Vector3D.minusK);
}
checkAngle(r.getAngle(), Math.PI);
@ -187,8 +186,9 @@ public class RotationTest
{ 0.0, 1.0, 0.0 },
{ 1.0, 0.0, 0.0 } };
r = new Rotation(m5, 1.0e-7);
fail("an exception should have been thrown");
fail("got " + r + ", should have caught an exception");
} catch (NotARotationMatrixException e) {
// expected
} catch (Exception e) {
fail("wrong exception caught");
}
@ -329,7 +329,7 @@ public class RotationTest
}
private void checkVector(Vector3D v1, Vector3D v2) {
assertTrue(Vector3D.subtract(v1, v2).getNorm() < 1.0e-10);
assertTrue(v1.subtract(v2).getNorm() < 1.0e-10);
}
private void checkAngle(double a1, double a2) {

19
src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java
View File

@ -43,29 +43,29 @@ public class Vector3DTest
Vector3D v1 = new Vector3D(1, 2, 3);
Vector3D v2 = new Vector3D(-3, -2, -1);
v1 = Vector3D.subtract(v1, v2);
v1 = v1.subtract(v2);
checkVector(v1, new Vector3D(4, 4, 4));
checkVector(Vector3D.subtract(v2, v1), new Vector3D(-7, -6, -5));
checkVector(v2.subtract(v1), new Vector3D(-7, -6, -5));
}
public void testAdd() {
Vector3D v1 = new Vector3D(1, 2, 3);
Vector3D v2 = new Vector3D(-3, -2, -1);
v1 = Vector3D.add(v1, v2);
v1 = v1.add(v2);
checkVector(v1, new Vector3D(-2, 0, 2));
checkVector(Vector3D.add(v2, v1), new Vector3D(-5, -2, 1));
checkVector(v2.add(v1), new Vector3D(-5, -2, 1));
}
public void testScalarProduct() {
Vector3D v = new Vector3D(1, 2, 3);
v = Vector3D.multiply(3, v);
v = v.multiply(3);
checkVector(v, new Vector3D(3, 6, 9));
checkVector(Vector3D.multiply(0.5, v), new Vector3D(1.5, 3, 4.5));
checkVector(v.multiply(0.5), new Vector3D(1.5, 3, 4.5));
}
@ -101,17 +101,16 @@ public class Vector3DTest
public void testAngularSeparation() {
Vector3D v1 = new Vector3D(2, -1, 4);
Vector3D k = Vector3D.normalize(v1);
Vector3D k = v1.normalize();
Vector3D i = k.orthogonal();
Vector3D v2 = new Vector3D(Math.cos(1.2), k, Math.sin(1.2), i);
Vector3D v2 = k.multiply(Math.cos(1.2)).add(i.multiply(Math.sin(1.2)));
assertTrue(Math.abs(Vector3D.angle(v1, v2) - 1.2) < 1.0e-12);
}
private void checkVector(Vector3D v1, Vector3D v2) {
assertTrue(Vector3D.subtract(v1, v2).getNorm() < 1.0e-12);
assertTrue(v1.subtract(v2).getNorm() < 1.0e-12);
}
public static Test suite() {

2
src/mantissa/tests-src/org/spaceroots/mantissa/linalg/DiagonalMatrixTest.java
View File

@ -149,7 +149,7 @@ public class DiagonalMatrixTest
boolean gotIt = false;
try {
d.setElement(3, 3, 0.0);
result = d.solve(b, 1.0e-10);
d.solve(b, 1.0e-10);
} catch (SingularMatrixException e) {
gotIt = true;
}

193
src/mantissa/tests-src/org/spaceroots/mantissa/linalg/GeneralMatrixTest.java
View File

@ -54,26 +54,14 @@ public class GeneralMatrixTest
}
public void testElements() {
Matrix m = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
checkMatrix(m, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix(m, new BilinearPattern(1.0, 0.01));
}
public void testCopy() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
GeneralMatrix m2 = new GeneralMatrix(m1);
@ -86,20 +74,12 @@ public class GeneralMatrixTest
assertTrue(m2.getRows() == m1.getRows());
assertTrue(m2.getColumns() == m1.getColumns());
checkMatrix(m2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix(m2, new BilinearPattern(1.0, 0.01));
}
public void testDuplicate() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = m1.duplicate();
assertTrue(m2 instanceof GeneralMatrix);
@ -113,11 +93,7 @@ public class GeneralMatrixTest
assertTrue(m2.getRows() == m1.getRows());
assertTrue(m2.getColumns() == m1.getColumns());
checkMatrix (m2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix (m2, new BilinearPattern(1.0, 0.01));
}
@ -133,53 +109,29 @@ public class GeneralMatrixTest
public void testAddOK() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = buildMatrix(m1.getRows(),
m1.getColumns(),
new ElementPattern() {
public double value(int i, int j) {
return 100 * i - 0.01 * j;
}
});
new BilinearPattern(100, -0.01));
Matrix m3 = m1.add(m2);
checkMatrix(m3, new ElementPattern() {
public double value(int i, int j) {
return 101 * i;
}
});
checkMatrix(m3, new BilinearPattern(101, 0));
}
public void testSelfAdd() {
GeneralMatrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralMatrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = buildMatrix(m1.getRows(),
m1.getColumns(),
new ElementPattern() {
public double value(int i, int j) {
return 100 * i - 0.01 * j;
}
});
new BilinearPattern(100, -0.01));
m1.selfAdd(m2);
checkMatrix(m1, new ElementPattern() {
public double value(int i, int j) {
return 101 * i;
}
});
checkMatrix(m1, new BilinearPattern(101, 0));
}
@ -195,53 +147,29 @@ public class GeneralMatrixTest
public void testSubOK() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = buildMatrix(m1.getRows(),
m1.getColumns(),
new ElementPattern() {
public double value(int i, int j) {
return 100 * i - 0.01 * j;
}
});
new BilinearPattern(100, -0.01));
Matrix m3 = m1.sub(m2);
checkMatrix(m3, new ElementPattern() {
public double value(int i, int j) {
return 0.02 * j - 99 * i;
}
});
checkMatrix(m3, new BilinearPattern(-99, 0.02));
}
public void testSelfSub() {
GeneralMatrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralMatrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = buildMatrix(m1.getRows(),
m1.getColumns(),
new ElementPattern() {
public double value(int i, int j) {
return 100 * i - 0.01 * j;
}
});
new BilinearPattern(100, -0.01));
m1.selfSub(m2);
checkMatrix(m1, new ElementPattern() {
public double value(int i, int j) {
return 0.02 * j - 99 * i;
}
});
checkMatrix(m1, new BilinearPattern(-99, 0.02));
}
@ -257,85 +185,46 @@ public class GeneralMatrixTest
public void testMulMOK() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = buildMatrix(m1.getColumns(), 4, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - j;
}
});
Matrix m2 = buildMatrix(m1.getColumns(), 4, new BilinearPattern(2, -1));
Matrix m3 = m1.mul(m2);
checkMatrix(m3, new ElementPattern() {
public double value(int i, int j) {
int p = 10; // must be equal to m1.getColumns()
return p * ((2 * i - 0.01 *j) * (p - 1) / 2.0
- i* j
+ (p - 1) * (2 * p - 1) / 300.0);
}
});
checkMatrix(m3, new ComplexPattern(m1.getColumns()));
}
public void testMulD() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = m1.mul(2.5);
checkMatrix(m2, new ElementPattern() {
public double value(int i, int j) {
return 2.5 * (i + 0.01 * j);
}
});
checkMatrix(m2, new BilinearPattern(2.5, 0.025));
}
public void testSelfMul() {
Matrix m = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
m.selfMul(2.5);
checkMatrix(m, new ElementPattern() {
public double value(int i, int j) {
return 2.5 * (i + 0.01 * j);
}
});
checkMatrix(m, new BilinearPattern(2.5, 0.025));
}
public void testTranspose() {
Matrix m1 = buildMatrix(5, 10, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m1 = buildMatrix(5, 10, new BilinearPattern(1.0, 0.01));
Matrix m2 = m1.getTranspose();
assertTrue(m1.getRows() == m2.getColumns());
assertTrue(m1.getColumns() == m2.getRows());
checkMatrix(m2, new ElementPattern() {
public double value(int i, int j) {
return 0.01 * i + j;
}
});
checkMatrix(m2, new BilinearPattern(0.01, 1.0));
}
@ -343,12 +232,36 @@ public class GeneralMatrixTest
return new TestSuite(GeneralMatrixTest.class);
}
public interface ElementPattern {
private interface ElementPattern {
public double value(int i, int j);
}
private static class BilinearPattern implements ElementPattern {
public BilinearPattern(double coeffI, double coeffJ) {
this.coeffI = coeffI;
this.coeffJ = coeffJ;
}
public double value(int i, int j) {
return coeffI * i + coeffJ * j;
}
private final double coeffI;
private final double coeffJ;
}
private static class ComplexPattern implements ElementPattern {
public ComplexPattern(int p) {
this.p = p;
}
public double value(int i, int j) {
return p * ((2 * i - 0.01 *j) * (p - 1) / 2.0
- i* j
+ (p - 1) * (2 * p - 1) / 300.0);
}
private final int p;
}
public GeneralMatrix buildMatrix(int rows, int columns,
ElementPattern pattern) {
BilinearPattern pattern) {
GeneralMatrix m = new GeneralMatrix(rows, columns);
for (int i = 0; i < m.getRows(); ++i) {

122
src/mantissa/tests-src/org/spaceroots/mantissa/linalg/GeneralSquareMatrixTest.java
View File

@ -54,26 +54,14 @@ public class GeneralSquareMatrixTest
}
public void testElements() {
Matrix m = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
Matrix m = buildMatrix(5, new BilinearPattern(1.0, 0.01));
checkMatrix(m, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix(m, new BilinearPattern(1.0, 0.01));
}
public void testCopy() {
GeneralSquareMatrix m1 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralSquareMatrix m1 = buildMatrix(5, new BilinearPattern(1.0, 0.01));
GeneralSquareMatrix m2 = new GeneralSquareMatrix(m1);
@ -86,20 +74,12 @@ public class GeneralSquareMatrixTest
assertTrue(m2.getRows() == m1.getRows());
assertTrue(m2.getColumns() == m1.getColumns());
checkMatrix(m2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix(m2, new BilinearPattern(1.0, 0.01));
}
public void testDuplicate() {
GeneralSquareMatrix m1 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralSquareMatrix m1 = buildMatrix(5, new BilinearPattern(1.0, 0.01));
Matrix m2 = m1.duplicate();
assertTrue(m2 instanceof GeneralSquareMatrix);
@ -113,59 +93,31 @@ public class GeneralSquareMatrixTest
assertTrue(m2.getRows() == m1.getRows());
assertTrue(m2.getColumns() == m1.getColumns());
checkMatrix(m2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
checkMatrix(m2, new BilinearPattern(1.0, 0.01));
}
public void testSelfAdd() {
GeneralSquareMatrix m1 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralSquareMatrix m1 = buildMatrix(5, new BilinearPattern(1.0, 0.01));
GeneralSquareMatrix m2 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.03 * j;
}
});
GeneralSquareMatrix m2 = buildMatrix(5, new BilinearPattern(2, -0.03));
m1.selfAdd(m2);
checkMatrix(m1, new ElementPattern() {
public double value(int i, int j) {
return 3 * i - 0.02 * j;
}
});
checkMatrix(m1, new BilinearPattern(3, -0.02));
}
public void testSelfSub() {
GeneralSquareMatrix m1 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return i + 0.01 * j;
}
});
GeneralSquareMatrix m1 = buildMatrix(5, new BilinearPattern(1.0, 0.01));
GeneralSquareMatrix m2 = buildMatrix(5, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.03 * j;
}
});
GeneralSquareMatrix m2 = buildMatrix(5, new BilinearPattern(2, -0.03));
m1.selfSub(m2);
checkMatrix(m1, new ElementPattern() {
public double value(int i, int j) {
return 0.04 * j - i;
}
});
checkMatrix(m1, new BilinearPattern(-1, 0.04));
}
@ -196,15 +148,16 @@ public class GeneralSquareMatrixTest
result = p.a.solve(p.b, 1.0e-10);
checkSolve(p, result);
boolean gotIt = false;
try {
p = buildProblem3();
result = p.a.solve(p.b, 1.0e-10);
fail("got " + result + ", should have caught an exception");
} catch(SingularMatrixException e) {
gotIt = true;
// expected
} catch(Exception e) {
fail("wrong exception caught: " + e.getMessage());
}
assertTrue(gotIt);
}
public void testInverse()
@ -214,28 +167,21 @@ public class GeneralSquareMatrixTest
a = buildProblem1().a;
inverse = a.getInverse(1.0e-10);
checkMatrix(a.mul(inverse), new ElementPattern() {
public double value(int i, int j) {
return (i == j) ? 1.0 : 0.0;
}
});
checkMatrix(a.mul(inverse), new IdentityPattern());
a = buildProblem2().a;
inverse = a.getInverse(1.0e-10);
checkMatrix(a.mul(inverse), new ElementPattern() {
public double value(int i, int j) {
return (i == j) ? 1.0 : 0.0;
}
});
checkMatrix(a.mul(inverse), new IdentityPattern());
boolean gotIt = false;
try {
a = buildProblem3().a;
inverse = a.getInverse(1.0e-10);
fail("got " + inverse + ", should have caught an exception");
} catch(SingularMatrixException e) {
gotIt = true;
// expected
} catch(Exception e) {
fail("wrong exception caught: " + e.getMessage());
}
assertTrue(gotIt);
}
@ -243,10 +189,28 @@ public class GeneralSquareMatrixTest
return new TestSuite(GeneralSquareMatrixTest.class);
}
public interface ElementPattern {
private interface ElementPattern {
public double value(int i, int j);
}
private static class BilinearPattern implements ElementPattern {
public BilinearPattern(double coeffI, double coeffJ) {
this.coeffI = coeffI;
this.coeffJ = coeffJ;
}
public double value(int i, int j) {
return coeffI * i + coeffJ * j;
}
private final double coeffI;
private final double coeffJ;
}
private static class IdentityPattern implements ElementPattern {
public double value(int i, int j) {
return (i == j) ? 1.0 : 0.0;
}
}
public GeneralSquareMatrix buildMatrix(int order,
ElementPattern pattern) {
GeneralSquareMatrix m = new GeneralSquareMatrix(order);
@ -270,7 +234,7 @@ public class GeneralSquareMatrixTest
}
}
private class LinearProblem {
private static class LinearProblem {
public GeneralSquareMatrix a;
public Matrix x;
public Matrix b;

88
src/mantissa/tests-src/org/spaceroots/mantissa/linalg/LowerTriangularMatrixTest.java
View File

@ -38,8 +38,7 @@ public class LowerTriangularMatrixTest
} else {
boolean gotIt = false;
try {
l.setElement
(i, j, -1.3);
l.setElement(i, j, -1.3);
} catch(ArrayIndexOutOfBoundsException e) {
gotIt = true;
}
@ -48,58 +47,34 @@ public class LowerTriangularMatrixTest
}
}
checkMatrix(l, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix(l, new BilinearPattern(1.0, 0.1));
}
public void testCopy() {
LowerTriangularMatrix l1 = buildMatrix(4, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
LowerTriangularMatrix l1 = buildMatrix(4, new BilinearPattern(1.0, 0.01));
LowerTriangularMatrix l2 = new LowerTriangularMatrix (l1);
checkMatrix (l2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix (l2, new BilinearPattern(1.0, 0.01));
}
public void testDuplicate() {
LowerTriangularMatrix l1 = buildMatrix(4, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
LowerTriangularMatrix l1 = buildMatrix(4, new BilinearPattern(1.0, 0.01));
Matrix l2 = l1.duplicate();
assertTrue(l2 instanceof LowerTriangularMatrix);
checkMatrix(l2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix(l2, new BilinearPattern(1.0, 0.01));
}
public void testTranspose() {
LowerTriangularMatrix l = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
LowerTriangularMatrix l = buildMatrix(7, new BilinearPattern(1.0, 0.1));
Matrix transposed = l.getTranspose();
assertTrue(transposed instanceof UpperTriangularMatrix);
@ -114,46 +89,23 @@ public class LowerTriangularMatrixTest
}
public void testSelfAdd() {
LowerTriangularMatrix l1 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 3 * i - 0.2 * j;
}
});
LowerTriangularMatrix l1 = buildMatrix(7, new BilinearPattern(3, -0.2));
LowerTriangularMatrix l2 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.4 * j; }
});
LowerTriangularMatrix l2 = buildMatrix(7, new BilinearPattern(2, -0.4));
l1.selfAdd(l2);
checkMatrix(l1, new ElementPattern() {
public double value(int i, int j) {
return 5 * i - 0.6 * j;
}
});
checkMatrix(l1, new BilinearPattern(5, -0.6));
}
public void testSelfSub() {
LowerTriangularMatrix l1 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 3 * i - 0.2 * j;
}
});
LowerTriangularMatrix l1 = buildMatrix(7, new BilinearPattern(3, -0.2));
LowerTriangularMatrix l2 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.4 * j;
}
});
LowerTriangularMatrix l2 = buildMatrix(7, new BilinearPattern(2, -0.4));
l1.selfSub(l2);
checkMatrix(l1, new ElementPattern() {
public double value(int i, int j) {
return i + 0.2 * j;
}
});
checkMatrix(l1, new BilinearPattern(1, 0.2));
}
public void testDeterminant() {
@ -199,7 +151,7 @@ public class LowerTriangularMatrixTest
boolean gotIt = false;
try {
l.setElement(3, 3, 0.0);
result = l.solve(b, 1.0e-10);
l.solve(b, 1.0e-10);
} catch(SingularMatrixException e) {
gotIt = true;
}
@ -235,6 +187,18 @@ public class LowerTriangularMatrixTest
public double value(int i, int j);
}
private static class BilinearPattern implements ElementPattern {
public BilinearPattern(double coeffI, double coeffJ) {
this.coeffI = coeffI;
this.coeffJ = coeffJ;
}
public double value(int i, int j) {
return coeffI * i + coeffJ * j;
}
private final double coeffI;
private final double coeffJ;
}
public LowerTriangularMatrix buildMatrix(int order,
ElementPattern pattern) {
LowerTriangularMatrix m = new LowerTriangularMatrix(order);

86
src/mantissa/tests-src/org/spaceroots/mantissa/linalg/UpperTriangularMatrixTest.java
View File

@ -46,58 +46,34 @@ public class UpperTriangularMatrixTest
}
}
checkMatrix(u, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix(u, new BilinearPattern(1.0, 0.1));
}
public void testCopy() {
UpperTriangularMatrix u1 = buildMatrix(4, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
UpperTriangularMatrix u1 = buildMatrix(4, new BilinearPattern(1.0, 0.1));
UpperTriangularMatrix u2 = new UpperTriangularMatrix(u1);
checkMatrix(u2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix(u2, new BilinearPattern(1.0, 0.1));
}
public void testDuplicate() {
UpperTriangularMatrix u1 = buildMatrix(4, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
UpperTriangularMatrix u1 = buildMatrix(4, new BilinearPattern(1.0, 0.1));
Matrix u2 = u1.duplicate();
assertTrue(u2 instanceof UpperTriangularMatrix);
checkMatrix(u2, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
checkMatrix(u2, new BilinearPattern(1.0, 0.1));
}
public void testTranspose() {
UpperTriangularMatrix u = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return i + 0.1 * j;
}
});
UpperTriangularMatrix u = buildMatrix(7, new BilinearPattern(1.0, 0.1));
Matrix transposed = u.getTranspose();
assertTrue(transposed instanceof LowerTriangularMatrix);
@ -112,47 +88,23 @@ public class UpperTriangularMatrixTest
}
public void testSelfAdd() {
UpperTriangularMatrix u1 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 3 * i - 0.2 * j;
}
});
UpperTriangularMatrix u1 = buildMatrix(7, new BilinearPattern(3, -0.2));
UpperTriangularMatrix u2 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.4 * j;
}
});
UpperTriangularMatrix u2 = buildMatrix(7, new BilinearPattern(2, -0.4));
u1.selfAdd(u2);
checkMatrix(u1, new ElementPattern() {
public double value(int i, int j) {
return 5 * i - 0.6 * j;
}
});
checkMatrix(u1, new BilinearPattern(5, -0.6));
}
public void testSelfSub() {
UpperTriangularMatrix u1 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 3 * i - 0.2 * j;
}
});
UpperTriangularMatrix u1 = buildMatrix(7, new BilinearPattern(3, -0.2));
UpperTriangularMatrix u2 = buildMatrix(7, new ElementPattern() {
public double value(int i, int j) {
return 2 * i - 0.4 * j;
}
});
UpperTriangularMatrix u2 = buildMatrix(7, new BilinearPattern(2, -0.4));
u1.selfSub(u2);
checkMatrix(u1, new ElementPattern() {
public double value(int i, int j) {
return i + 0.2 * j;
}
});
checkMatrix(u1, new BilinearPattern(1, 0.2));
}
public void testDeterminant() {
@ -199,7 +151,7 @@ public class UpperTriangularMatrixTest
boolean gotIt = false;
try {
u.setElement(3, 3, 0.0);
result = u.solve(b, 1.0e-10);
u.solve(b, 1.0e-10);
} catch(SingularMatrixException e) {
gotIt = true;
}
@ -235,6 +187,18 @@ public class UpperTriangularMatrixTest
public double value(int i, int j);
}
private static class BilinearPattern implements ElementPattern {
public BilinearPattern(double coeffI, double coeffJ) {
this.coeffI = coeffI;
this.coeffJ = coeffJ;
}
public double value(int i, int j) {
return coeffI * i + coeffJ * j;
}
private final double coeffI;
private final double coeffJ;
}
public UpperTriangularMatrix buildMatrix(int order,
ElementPattern pattern) {
UpperTriangularMatrix m = new UpperTriangularMatrix (order);

83
src/mantissa/tests-src/org/spaceroots/mantissa/ode/ClassicalRungeKuttaIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.spaceroots.mantissa.ode;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.fitting.PolynomialFitter;
public class ClassicalRungeKuttaIntegratorTest
@ -72,11 +73,9 @@ public class ClassicalRungeKuttaIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -95,9 +94,8 @@ public class ClassicalRungeKuttaIntegratorTest
public void testOrder()
throws EstimationException, DerivativeException,
IntegratorException {
PolynomialFitter fitter = new PolynomialFitter(1,
10, 1.0e-7, 1.0e-10,
1.0e-10);
PolynomialFitter fitter =
new PolynomialFitter(1, new LevenbergMarquardtEstimator());
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
@ -112,11 +110,9 @@ public class ClassicalRungeKuttaIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -177,38 +173,45 @@ public class ClassicalRungeKuttaIntegratorTest
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
integ.setStepHandler(new StepHandler() {
private double maxError = 0;
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
if (isLast) {
// even with more than 1000 evaluations per period,
// RK4 is not able to integrate such an eccentric
// orbit with a good accuracy
assertTrue(maxError > 0.005);
}
}
});
integ.setStepHandler(new KeplerHandler(pb));
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
}
private static class KeplerHandler implements StepHandler {
public KeplerHandler(TestProblem3 pb) {
this.pb = pb;
reset();
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
if (isLast) {
// even with more than 1000 evaluations per period,
// RK4 is not able to integrate such an eccentric
// orbit with a good accuracy
assertTrue(maxError > 0.005);
}
}
private double maxError = 0;
private TestProblem3 pb;
}
public static Test suite() {
return new TestSuite(ClassicalRungeKuttaIntegratorTest.class);
}

2
src/mantissa/tests-src/org/spaceroots/mantissa/ode/ContinuousOutputModelTest.java
View File

@ -25,6 +25,8 @@ public class ContinuousOutputModelTest
public ContinuousOutputModelTest(String name) {
super(name);
pb = null;
integ = null;
}
public void testBoundaries()

170
src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince54IntegratorTest.java
View File

@ -91,7 +91,7 @@ public class DormandPrince54IntegratorTest
}
private class DP54SmallLastHandler implements StepHandler {
private static class DP54SmallLastHandler implements StepHandler {
public DP54SmallLastHandler(double minStep) {
lastSeen = false;
@ -170,11 +170,9 @@ public class DormandPrince54IntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -197,43 +195,7 @@ public class DormandPrince54IntegratorTest
FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
integ.setStepHandler(new StepHandler() {
private int nbSteps = 0;
private double maxError = 0;
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
nbSteps = 0;
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
++nbSteps;
for (int a = 1; a < 10; ++a) {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
double interp = ((10 - a) * prev + a * curr) / 10;
interpolator.setInterpolatedTime(interp);
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
}
if (isLast) {
assertTrue(maxError < 7.0e-10);
assertTrue(nbSteps < 400);
}
}
});
integ.setStepHandler(new KeplerHandler(pb));
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -254,47 +216,97 @@ public class DormandPrince54IntegratorTest
FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
integ.setStepHandler(new StepHandler() {
private boolean firstTime = true;
private double minStep = 0;
private double maxStep = 0;
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
firstTime = true;
minStep = 0;
maxStep = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double step = Math.abs(interpolator.getCurrentTime()
- interpolator.getPreviousTime());
if (firstTime) {
minStep = Math.abs(step);
maxStep = minStep;
firstTime = false;
} else {
if (step < minStep) {
minStep = step;
}
if (step > maxStep) {
maxStep = step;
}
}
if (isLast) {
assertTrue(minStep < (1.0 / 450.0));
assertTrue(maxStep > (1.0 / 4.2));
}
}
});
integ.setStepHandler(new VariableHandler());
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
}
private static class KeplerHandler implements StepHandler {
public KeplerHandler(TestProblem3 pb) {
this.pb = pb;
reset();
}
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
nbSteps = 0;
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
++nbSteps;
for (int a = 1; a < 10; ++a) {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
double interp = ((10 - a) * prev + a * curr) / 10;
interpolator.setInterpolatedTime(interp);
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
}
if (isLast) {
assertTrue(maxError < 7.0e-10);
assertTrue(nbSteps < 400);
}
}
private int nbSteps;
private double maxError;
private TestProblem3 pb;
}
private static class VariableHandler implements StepHandler {
public VariableHandler() {
firstTime = true;
minStep = 0;
maxStep = 0;
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
firstTime = true;
minStep = 0;
maxStep = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double step = Math.abs(interpolator.getCurrentTime()
- interpolator.getPreviousTime());
if (firstTime) {
minStep = Math.abs(step);
maxStep = minStep;
firstTime = false;
} else {
if (step < minStep) {
minStep = step;
}
if (step > maxStep) {
maxStep = step;
}
}
if (isLast) {
assertTrue(minStep < (1.0 / 450.0));
assertTrue(maxStep > (1.0 / 4.2));
}
}
private boolean firstTime;
private double minStep;
private double maxStep;
}
public static Test suite() {
return new TestSuite(DormandPrince54IntegratorTest.class);
}

207
src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
View File

@ -130,11 +130,9 @@ public class DormandPrince853IntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -157,43 +155,7 @@ public class DormandPrince853IntegratorTest
FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
integ.setStepHandler(new StepHandler() {
private int nbSteps = 0;
private double maxError = 0;
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
nbSteps = 0;
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
++nbSteps;
for (int a = 1; a < 10; ++a) {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
double interp = ((10 - a) * prev + a * curr) / 10;
interpolator.setInterpolatedTime(interp);
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
}
if (isLast) {
assertTrue(maxError < 2.4e-10);
assertTrue(nbSteps < 150);
}
}
});
integ.setStepHandler(new KeplerHandler(pb));
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -214,42 +176,7 @@ public class DormandPrince853IntegratorTest
FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
integ.setStepHandler(new StepHandler() {
private boolean firstTime = true;
private double minStep = 0;
private double maxStep = 0;
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
firstTime = true;
minStep = 0;
maxStep = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double step = Math.abs(interpolator.getCurrentTime()
- interpolator.getPreviousTime());
if (firstTime) {
minStep = Math.abs(step);
maxStep = minStep;
firstTime = false;
} else {
if (step < minStep) {
minStep = step;
}
if (step > maxStep) {
maxStep = step;
}
}
if (isLast) {
assertTrue(minStep < (1.0 / 100.0));
assertTrue(maxStep > (1.0 / 2.0));
}
}
});
integ.setStepHandler(new VariableHandler());
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -267,35 +194,13 @@ public class DormandPrince853IntegratorTest
FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
integ.setStepHandler(new StepHandler() {
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
}
});
integ.setStepHandler(DummyStepHandler.getInstance());
integ.integrate(pb1,
pb1.getInitialTime(), pb1.getInitialState(),
pb1.getFinalTime(), new double[pb1.getDimension()]);
int callsWithoutDenseOutput = pb1.getCalls();
integ.setStepHandler(new StepHandler() {
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
interpolator.setInterpolatedTime(0.5*(prev + curr));
}
});
integ.setStepHandler(new InterpolatingStepHandler());
integ.integrate(pb2,
pb2.getInitialTime(), pb2.getInitialState(),
pb2.getFinalTime(), new double[pb2.getDimension()]);
@ -316,6 +221,104 @@ public class DormandPrince853IntegratorTest
assertEquals(8.0, y[0], 1.0e-12);
}
private static class KeplerHandler implements StepHandler {
public KeplerHandler(TestProblem3 pb) {
this.pb = pb;
reset();
}
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
nbSteps = 0;
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
++nbSteps;
for (int a = 1; a < 10; ++a) {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
double interp = ((10 - a) * prev + a * curr) / 10;
interpolator.setInterpolatedTime(interp);
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
}
if (isLast) {
assertTrue(maxError < 2.4e-10);
assertTrue(nbSteps < 150);
}
}
private int nbSteps;
private double maxError;
private TestProblem3 pb;
}
private static class VariableHandler implements StepHandler {
public VariableHandler() {
reset();
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
firstTime = true;
minStep = 0;
maxStep = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double step = Math.abs(interpolator.getCurrentTime()
- interpolator.getPreviousTime());
if (firstTime) {
minStep = Math.abs(step);
maxStep = minStep;
firstTime = false;
} else {
if (step < minStep) {
minStep = step;
}
if (step > maxStep) {
maxStep = step;
}
}
if (isLast) {
assertTrue(minStep < (1.0 / 100.0));
assertTrue(maxStep > (1.0 / 2.0));
}
}
private boolean firstTime = true;
private double minStep = 0;
private double maxStep = 0;
}
private static class InterpolatingStepHandler implements StepHandler {
public boolean requiresDenseOutput() {
return true;
}
public void reset() {
}
public void handleStep(StepInterpolator interpolator,
boolean isLast)
throws DerivativeException {
double prev = interpolator.getPreviousTime();
double curr = interpolator.getCurrentTime();
interpolator.setInterpolatedTime(0.5*(prev + curr));
}
}
public static Test suite() {
return new TestSuite(DormandPrince853IntegratorTest.class);
}

22
src/mantissa/tests-src/org/spaceroots/mantissa/ode/EulerIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.spaceroots.mantissa.ode;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.fitting.PolynomialFitter;
public class EulerIntegratorTest
@ -59,11 +60,9 @@ public class EulerIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -84,9 +83,8 @@ public class EulerIntegratorTest
public void testOrder()
throws EstimationException, DerivativeException,
IntegratorException {
PolynomialFitter fitter = new PolynomialFitter(1,
10, 1.0e-7, 1.0e-10,
1.0e-10);
PolynomialFitter fitter =
new PolynomialFitter(1, new LevenbergMarquardtEstimator());
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
@ -101,11 +99,9 @@ public class EulerIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),

5
src/mantissa/tests-src/org/spaceroots/mantissa/ode/EulerStepInterpolatorTest.java
View File

@ -55,8 +55,7 @@ public class EulerStepInterpolatorTest
double t0 = 0;
double[] y0 = {0.0, 1.0, -2.0};
double[] y = new double[y0.length];
System.arraycopy(y0, 0, y, 0, y0.length);
double[] y = (double[]) y0.clone();
double[][] yDot = { new double[y0.length] };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
interpolator.reinitialize(new DummyEquations(), y, yDot, true);
@ -154,7 +153,7 @@ public class EulerStepInterpolatorTest
}
private class DummyEquations
private static class DummyEquations
implements FirstOrderDifferentialEquations {
public int getDimension() {
return 0;

2
src/mantissa/tests-src/org/spaceroots/mantissa/ode/FirstOrderConverterTest.java
View File

@ -69,7 +69,7 @@ public class FirstOrderConverterTest
return new TestSuite(FirstOrderConverterTest.class);
}
private class Equations
private static class Equations
implements SecondOrderDifferentialEquations {
private int n;

83
src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.spaceroots.mantissa.ode;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.fitting.PolynomialFitter;
public class GillIntegratorTest
@ -59,11 +60,9 @@ public class GillIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -82,9 +81,8 @@ public class GillIntegratorTest
public void testOrder()
throws EstimationException, DerivativeException,
IntegratorException {
PolynomialFitter fitter = new PolynomialFitter(1,
10, 1.0e-7, 1.0e-10,
1.0e-10);
PolynomialFitter fitter =
new PolynomialFitter(1, new LevenbergMarquardtEstimator());
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
@ -99,11 +97,9 @@ public class GillIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb);
integ.setStepHandler(handler);
SwitchingFunction[] functions = pb.getSwitchingFunctions();
if (functions != null) {
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -164,33 +160,7 @@ public class GillIntegratorTest
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
FirstOrderIntegrator integ = new GillIntegrator(step);
integ.setStepHandler(new StepHandler() {
private double maxError = 0;
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
if (isLast) {
// even with more than 1000 evaluations per period,
// RK4 is not able to integrate such an eccentric
// orbit with a good accuracy
assertTrue(maxError > 0.001);
}
}
});
integ.setStepHandler(new KeplerStepHandler(pb));
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -206,6 +176,39 @@ public class GillIntegratorTest
assertEquals(8.0, y[0], 1.0e-12);
}
private static class KeplerStepHandler implements StepHandler {
public KeplerStepHandler(TestProblem3 pb) {
this.pb = pb;
reset();
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) {
double[] interpolatedY = interpolator.getInterpolatedState ();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
if (isLast) {
// even with more than 1000 evaluations per period,
// RK4 is not able to integrate such an eccentric
// orbit with a good accuracy
assertTrue(maxError > 0.001);
}
}
private double maxError;
private TestProblem3 pb;
}
public static Test suite() {
return new TestSuite(GillIntegratorTest.class);
}

Some files were not shown because too many files have changed in this diff Show More