diff --git a/doap_math.rdf b/doap_math.rdf index ea3368f3f..85c4a08d2 100644 --- a/doap_math.rdf +++ b/doap_math.rdf @@ -30,10 +30,10 @@ Apache Commons Math The Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang. - - - - + + + + diff --git a/src/changes/changes.xml b/src/changes/changes.xml index dc6539cd2..132509b4f 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -59,6 +59,9 @@ If the output is not quite correct, check for invisible trailing spaces! when calling "optimize(...)" with linear constraints whose dimension does not match the dimension of the objective function. + + Fixed wrong event detection in case of close events pairs. + Reimplemented pow(double, double) in FastMath, for better accuracy in integral power cases and trying to fix erroneous JIT optimization again. diff --git a/src/main/java/org/apache/commons/math4/ode/AbstractIntegrator.java b/src/main/java/org/apache/commons/math4/ode/AbstractIntegrator.java index b9127ad86..e9326d9c6 100644 --- a/src/main/java/org/apache/commons/math4/ode/AbstractIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/AbstractIntegrator.java @@ -292,9 +292,12 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { * @param yDot placeholder array where to put the time derivative of the state vector * @exception MaxCountExceededException if the number of functions evaluations is exceeded * @exception DimensionMismatchException if arrays dimensions do not match equations settings + * @exception NullPointerException if the ODE equations have not been set (i.e. if this method + * is called outside of a call to {@link #integrate(ExpandableStatefulODE, double)} or {@link + * #integrate(FirstOrderDifferentialEquations, double, double[], double, double[])}) */ public void computeDerivatives(final double t, final double[] y, final double[] yDot) - throws MaxCountExceededException, DimensionMismatchException { + throws MaxCountExceededException, DimensionMismatchException, NullPointerException { evaluations.incrementCount(); expandable.computeDerivatives(t, y, yDot); } diff --git a/src/main/java/org/apache/commons/math4/ode/events/EventState.java b/src/main/java/org/apache/commons/math4/ode/events/EventState.java index fe3039a4c..1908440ce 100644 --- a/src/main/java/org/apache/commons/math4/ode/events/EventState.java +++ b/src/main/java/org/apache/commons/math4/ode/events/EventState.java @@ -296,7 +296,18 @@ public class EventState { ta = forward ? ta + convergence : ta - convergence; ga = f.value(ta); } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb))); - --i; + + if (forward ^ (ta >= tb)) { + // we were able to skip this spurious root + --i; + } else { + // we can't avoid this root before the end of the step, + // we have to handle it despite it is close to the former one + // maybe we have two very close roots + pendingEventTime = root; + pendingEvent = true; + return true; + } } else if (Double.isNaN(previousEventTime) || (FastMath.abs(previousEventTime - root) > convergence)) { pendingEventTime = root; diff --git a/src/main/java/org/apache/commons/math4/util/ResizableDoubleArray.java b/src/main/java/org/apache/commons/math4/util/ResizableDoubleArray.java index 09fd7482b..6377141a0 100644 --- a/src/main/java/org/apache/commons/math4/util/ResizableDoubleArray.java +++ b/src/main/java/org/apache/commons/math4/util/ResizableDoubleArray.java @@ -179,7 +179,6 @@ public class ResizableDoubleArray implements DoubleArray, Serializable { * The input array is copied, not referenced. * Other properties take default values: *
    - *
  • {@code initialCapacity = 16}
    • *
  • {@code expansionMode = MULTIPLICATIVE}
    • *
  • {@code expansionFactor = 2.0}
    • *
  • {@code contractionCriterion = 2.5}
    • @@ -189,7 +188,9 @@ public class ResizableDoubleArray implements DoubleArray, Serializable { * @since 2.2 */ public ResizableDoubleArray(double[] initialArray) { - this(DEFAULT_INITIAL_CAPACITY, + this(initialArray == null || initialArray.length == 0 ? + DEFAULT_INITIAL_CAPACITY : + initialArray.length, DEFAULT_EXPANSION_FACTOR, DEFAULT_CONTRACTION_DELTA + DEFAULT_EXPANSION_FACTOR, DEFAULT_EXPANSION_MODE, diff --git a/src/site/design/differentiation.puml b/src/site/design/differentiation.puml index ce24df740..8649345de 100644 --- a/src/site/design/differentiation.puml +++ b/src/site/design/differentiation.puml @@ -24,7 +24,7 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4 #ECEBD8 + package org.apache.commons.math4 #ECEBD8 { interface "FieldElement" as FieldElement_T_ { T add(T a) @@ -37,12 +37,12 @@ Field getField() } - package analysis #DDEBD8 + package analysis #DDEBD8 { interface UnivariateFunction { double value(double x) } - package differentiation #DDDBD8 + package differentiation #DDDBD8 { class DerivativeStructure { -DerivativeStructure(DSCompiler compiler) @@ -98,8 +98,8 @@ UnivariateDifferentiable <-- UnivariateDifferentiator : creates UnivariateDifferentiable --> DerivativeStructure : uses - end package - end package - end package + } + } + } @enduml diff --git a/src/site/design/oneD.puml b/src/site/design/oneD.puml index 368e9ab58..3f02d6e54 100644 --- a/src/site/design/oneD.puml +++ b/src/site/design/oneD.puml @@ -24,7 +24,7 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4.geometry #ECEBD8 + package org.apache.commons.math4.geometry #ECEBD8 { interface Space { +int getDimension() @@ -53,15 +53,15 @@ Space <-- Vector_S_ - package partitioning #DDEBD8 + package partitioning #DDEBD8 { interface "Region" as Region_S_ interface "Hyperplane" as Hyperplane_S_ interface "SubHyperplane" as SubHyperplane_S_ - end package + } - package euclidean #DDEBD8 + package euclidean #DDEBD8 { - package oned #DDDBD8 + package oned #DDDBD8 { class Euclidean1D class OrientedPoint @@ -72,10 +72,10 @@ Vector_S_ <|.. OrientedPoint Region_S_ <|.. IntervalSet - end package + } - end package + } - end package + } @enduml diff --git a/src/site/design/partitioning.puml b/src/site/design/partitioning.puml index dedfe24b2..bdf27dc83 100644 --- a/src/site/design/partitioning.puml +++ b/src/site/design/partitioning.puml @@ -24,9 +24,9 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4.geometry #ECEBD8 + package org.apache.commons.math4.geometry #ECEBD8 { - package partitioning #DDEBD8 + package partitioning #DDEBD8 { abstract "AbstractSubHyperplane" as AbstractSubHyperplane_S_T_ note left @@ -81,7 +81,7 @@ a region is a part of the space it may be empty, it may cover the whole space, - it may have infinite non-closed boudaries, + it may have infinite non-closed boundaries, it may be in several disconnected sub-parts, it may contain holes end note @@ -112,8 +112,7 @@ Region_S_ --> Side BSPTree_S_ <-- BSPTreeVisitor_S_ : visits - end package - - end package + } + } @enduml diff --git a/src/site/design/solver.puml b/src/site/design/solver.puml index 929a46d05..535a440f6 100644 --- a/src/site/design/solver.puml +++ b/src/site/design/solver.puml @@ -24,7 +24,7 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4.differentiation.solvers #ECEBD8 + package org.apache.commons.math4.differentiation.solvers #ECEBD8 { enum AllowedSolution { ANY_SIDE @@ -97,6 +97,6 @@ abstract class BaseSecantSolver AbstractPolynomialSolver <|-- LaguerreSolver AbstractUnivariateSolver <|-- BisectionSolver - end package + } @enduml diff --git a/src/site/design/threeD.puml b/src/site/design/threeD.puml index 81daf05a5..43175f632 100644 --- a/src/site/design/threeD.puml +++ b/src/site/design/threeD.puml @@ -24,7 +24,7 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4.geometry #ECEBD8 + package org.apache.commons.math4.geometry #ECEBD8 { interface Space { +int getDimension() @@ -53,15 +53,15 @@ Space <-- Vector_S_ - package partitioning #DDEBD8 + package partitioning #DDEBD8 { interface "Region" as Region_S_ interface "Hyperplane" as Hyperplane_S_ interface "SubHyperplane" as SubHyperplane_S_ - end package + } - package euclidean #DDEBD8 + package euclidean #DDEBD8 { - package threed #DDDBD8 + package threed #DDDBD8 { class Euclidean3D class Vector3D @@ -77,10 +77,10 @@ SubHyperplane_S_ <|.. SubPlane Region_S_ <|.. PolyhedronsSet - end package + } - end package + } - end package + } @enduml diff --git a/src/site/design/twoD.puml b/src/site/design/twoD.puml index 3de2827a6..b7f02d728 100644 --- a/src/site/design/twoD.puml +++ b/src/site/design/twoD.puml @@ -24,7 +24,7 @@ skinparam NoteFontColor #691616 skinparam ClassFontSize 11 - package org.apache.commons.math4.geometry #ECEBD8 + package org.apache.commons.math4.geometry #ECEBD8 { interface Space { +int getDimension() @@ -53,15 +53,15 @@ Space <-- Vector_S_ - package partitioning #DDEBD8 + package partitioning #DDEBD8 { interface "Region" as Region_S_ interface "Hyperplane" as Hyperplane_S_ interface "SubHyperplane" as SubHyperplane_S_ - end package + } - package euclidean #DDEBD8 + package euclidean #DDEBD8 { - package twod #DDDBD8 + package twod #DDDBD8 { class Euclidean2D class Vector2D @@ -75,10 +75,10 @@ SubHyperplane_S_ <|.. SubLine Region_S_ <|.. PolygonsSet - end package + } - end package + } - end package + } @enduml diff --git a/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java b/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java index b2041dcd3..bbb2bc51a 100644 --- a/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java +++ b/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java @@ -29,6 +29,7 @@ import org.apache.commons.math4.ode.SecondaryEquations; import org.apache.commons.math4.ode.events.EventHandler; import org.apache.commons.math4.ode.events.EventState; import org.apache.commons.math4.ode.nonstiff.DormandPrince853Integrator; +import org.apache.commons.math4.ode.nonstiff.LutherIntegrator; import org.apache.commons.math4.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math4.ode.sampling.DummyStepInterpolator; import org.junit.Assert; @@ -43,21 +44,9 @@ public class EventStateTest { final double r1 = 90.0; final double r2 = 135.0; final double gap = r2 - r1; - EventHandler closeEventsGenerator = new EventHandler() { - public void init(double t0, double[] y0, double t) { - } - public void resetState(double t, double[] y) { - } - public double g(double t, double[] y) { - return (t - r1) * (r2 - t); - } - public Action eventOccurred(double t, double[] y, boolean increasing) { - return Action.CONTINUE; - } - }; final double tolerance = 0.1; - EventState es = new EventState(closeEventsGenerator, 1.5 * gap, + EventState es = new EventState(new CloseEventsGenerator(r1, r2), 1.5 * gap, tolerance, 100, new BrentSolver(tolerance)); es.setExpandable(new ExpandableStatefulODE(new FirstOrderDifferentialEquations() { @@ -226,4 +215,63 @@ public class EventStateTest { } + @Test + public void testEventsCloserThanThreshold() + throws DimensionMismatchException, NumberIsTooSmallException, + MaxCountExceededException, NoBracketingException { + + FirstOrderDifferentialEquations equation = new FirstOrderDifferentialEquations() { + + public int getDimension() { + return 1; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = 1.0; + } + }; + + LutherIntegrator integrator = new LutherIntegrator(20.0); + CloseEventsGenerator eventsGenerator = + new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128); + integrator.addEventHandler(eventsGenerator, 1.0, 0.02, 1000); + double[] y = new double[1]; + double tEnd = integrator.integrate(equation, 0.0, y, 100.0, y); + Assert.assertEquals( 2, eventsGenerator.getCount()); + Assert.assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0); + + } + + private class CloseEventsGenerator implements EventHandler { + + final double r1; + final double r2; + int count; + + public CloseEventsGenerator(final double r1, final double r2) { + this.r1 = r1; + this.r2 = r2; + this.count = 0; + } + + public void init(double t0, double[] y0, double t) { + } + + public void resetState(double t, double[] y) { + } + + public double g(double t, double[] y) { + return (t - r1) * (r2 - t); + } + + public Action eventOccurred(double t, double[] y, boolean increasing) { + return ++count < 2 ? Action.CONTINUE : Action.STOP; + } + + public int getCount() { + return count; + } + + } + }