List of usage examples for org.apache.commons.math3.ode FirstOrderIntegrator integrate
double integrate(FirstOrderDifferentialEquations equations, double t0, double[] y0, double t, double[] y) throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException;
From source file:beast.structuredCoalescent.distribution.Masco.java
public double calculateLogP() { // newly calculate tree intervals treeIntervalsInput.get().calculateIntervals(); // correctly calculate the daughter nodes at coalescent intervals in the case of // bifurcation or in case two nodes are at the same height treeIntervalsInput.get().swap();//from ww w.jav a 2 s .c o m // Set up ArrayLists for the indices of active lineages and the lineage state probabilities activeLineages = new ArrayList<Integer>(); lineStateProbs = new ArrayList<Double>(); // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards logP = 0; // set the current time double currTime = 0.0; // total number of intervals final int intervalCount = treeIntervalsInput.get().getIntervalCount(); // interval time counter int t = 0; // initialize the number of lineages nr_lineages = 0; // Captures the probabilities of lineages being in a state double[] p; // Initialize the migration rates matrix double[][] migration_rates = new double[states][states]; int c = 0; for (int k = 0; k < states; k++) { for (int l = 0; l < states; l++) { if (k != l) { migration_rates[k][l] = migrationRatesInput.get().getArrayValue(c); c++; } else { // diagonal migration_rates[k][l] = 0.0; } } } // Initialize the coalescent rates double[] coalescent_rates = new double[states]; for (int k = 0; k < states; k++) { coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2;//(epiModelInput.get().getF(t,k,k) / (Y.get(k)*Y.get(k))); } // integrate until there are no more tree intervals do { double nextIntervalTime = treeIntervalsInput.get().getInterval(t); // Length of the current interval final double duration = nextIntervalTime;// - currTime; // if the current interval has a length greater than 0, integrate if (duration > 0) { if (dependentHistory) p = new double[lineStateProbs.size()]; // Captures the probabilities of lineages being in a state else p = new double[lineStateProbs.size() + 1]; // Captures the probabilities of lineages being in a state, last one keeps track of the probability // convert the array list to double[] for (int i = 0; i < lineStateProbs.size(); i++) p[i] = lineStateProbs.get(i); // not needed if (!dependentHistory) p[lineStateProbs.size()] = 1; double[] p_for_ode = new double[p.length]; double ts = 0.0; // If proportial time step is true, set the integration time for the given interval // inverse proportional to the number of lineages if (propTimeStep) ts = timeStep / lineStateProbs.size(); else ts = timeStep; // Never choose a longer time step than the integration window if (duration < (ts / 2)) ts = duration / 2; FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts); // set the odes FirstOrderDifferentialEquations ode = new ode_masco(migration_rates, coalescent_rates, nr_lineages, states); // integrate integrator.integrate(ode, 0, p, duration, p_for_ode); // If the Dimension is larger than the maximum integer, at least one state prob is below 0 and the step is rejected if (ode.getDimension() == Integer.MAX_VALUE) { System.out.println(lineStateProbs.size()); System.out.println("lalalallal"); return Double.NEGATIVE_INFINITY; } for (int i = 0; i < lineStateProbs.size(); i++) lineStateProbs.set(i, p_for_ode[i]); } // update the time currTime = nextIntervalTime; // event is coalescent event if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) { logP += coalesce(t); nr_lineages--; } // event is sampling event if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) { logP += normalizeLineages(); addLineages(t); nr_lineages++; } // update the interval number t++; } while (t < intervalCount); //Compute likelihood of remaining tree intervals (coal events occuring before origin) if (Double.isInfinite(logP)) logP = Double.NEGATIVE_INFINITY; if (max_posterior < logP && logP < 0) { max_posterior = logP; max_mig = new double[states * (states - 1)]; max_coal = new double[states]; for (int i = 0; i < 1; i++) max_mig[i] = migrationRatesInput.get().getArrayValue(i); for (int i = 0; i < 1; i++) max_coal[i] = coalescentRatesInput.get().getArrayValue(i); } // System.exit(0); return logP; }
From source file:beast.structuredCoalescent.distribution.ExactStructuredCoalescent.java
public double calculateLogP() { // Calculate the tree intervals (time between events, which nodes participate at a event etc.) treeIntervalsInput.get().calculateIntervals(); treeIntervalsInput.get().swap();//from w ww .ja v a 2 s .c o m // Set up for lineage state probabilities activeLineages = new ArrayList<Integer>(); lineStateProbs = new ArrayList<Double>(); // Compute likelihood at each integration time and tree event starting at final sampling time and moving backwards logP = 0; // Initialize the line state probabilities // total number of intervals final int intervalCount = treeIntervalsInput.get().getIntervalCount(); // counts in which interval we are in int t = 0; nr_lineages = 0; // Captures the probabilities of lineages being in a state double[] p; // Initialize the migration rates matrix int c = 0; for (int k = 0; k < states; k++) { for (int l = 0; l < states; l++) { if (k != l) { migration_rates[c] = migrationRatesInput.get().getArrayValue(c); migration_map[k][l] = c; c++; } else { coalescent_rates[k] = coalescentRatesInput.get().getArrayValue(k) / 2; } } } boolean first = true; // integrate until there are no more tree intervals do { double nextIntervalTime = treeIntervalsInput.get().getInterval(t); // Length of the current interval final double duration = nextIntervalTime;// - currTime; // if the current interval has a length greater than 0, integrate if (duration > 0) { p = new double[jointStateProbabilities.size()]; // Captures the probabilities of lineages being in a state // convert the array list to double[] for (int i = 0; i < jointStateProbabilities.size(); i++) p[i] = jointStateProbabilities.get(i); double[] p_for_ode = new double[p.length]; double ts = timeStep; if (duration < timeStep) ts = duration / 2; // initialize integrator FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(ts); // set the odes FirstOrderDifferentialEquations ode = new ode_integrator(migration_rates, coalescent_rates, nr_lineages, states, connectivity, sums); // integrate integrator.integrate(ode, 0, p, duration, p_for_ode); // if the dimension is equal to the max integer, this means that a calculation // of a probability of a configuration resulted in a value below 0 and the // run will be stopped if (ode.getDimension() == Integer.MAX_VALUE) { System.out.println("lalalallal"); return Double.NEGATIVE_INFINITY; } // set the probabilities of the system being in a configuration again for (int i = 0; i < p_for_ode.length; i++) jointStateProbabilities.set(i, p_for_ode[i]); } /* * compute contribution of event to likelihood */ if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.COALESCENT) { nr_lineages--; logP += coalesce(t); } /* * add new lineage */ if (treeIntervalsInput.get().getIntervalType(t) == IntervalType.SAMPLE) { nr_lineages++; addLineages(t, first); first = false; } t++; } while (t < intervalCount); //Compute likelihood of remaining tree intervals (coal events occuring before origin) if (Double.isInfinite(logP)) logP = Double.NEGATIVE_INFINITY; if (max_posterior < logP && logP < 0) { max_posterior = logP; max_mig = new double[states * (states - 1)]; max_coal = new double[states]; for (int i = 0; i < 1; i++) max_mig[i] = migrationRatesInput.get().getArrayValue(i); for (int i = 0; i < 1; i++) max_coal[i] = coalescentRatesInput.get().getArrayValue(i); } return logP; }
From source file:eu.mihosoft.fx.tutorials.gravity.SolarSystem.java
/** * Initializes the ui frame listener./*from ww w.j a va 2 s . c om*/ * * @param ode ode to solve * @param y state vector (location and velocity) * @param m particle masses * @param ignoreFlags ignore flags (used for collisions) * @param dt time step size */ private void initFrameListener(FirstOrderDifferentialEquations ode, double[] y, final double[] m, final boolean[] ignoreFlags, double dt) { // local time scale (sec per h * h per day * days per year) final double localTimeScale = 3600 * 24 * 365; // integrator FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(dt * localTimeScale); double[] yPrev = new double[y.length]; double[] interpolatedY = new double[y.length]; // create frame listener frameListener = new AnimationTimer() { @Override public void handle(long now) { // thanks to http://gafferongames.com/game-physics/fix-your-timestep/ // measure elapsed time between last and current pulse (frame) double frameDuration = (now - lastTimeStamp) / 1e9; lastTimeStamp = now; // we don't allow frame durations above 2*dt if (frameDuration > 2 * dt) { frameDuration = 2 * dt; } // add elapsed time to remaining simulation interval remainingSimulationTime += frameDuration; // copy current state to prev state System.arraycopy(y, 0, yPrev, 0, yPrev.length); // simulate remaining interval while (remainingSimulationTime >= dt) { double scaledT = time * localTimeScale; double scaledDT = dt * localTimeScale * timeScale; // integrate one step try { integrator.integrate(ode, scaledT, y, scaledT + scaledDT, y); } catch (Exception ex) { ex.printStackTrace(System.err); } // remove integrated interval from remaining simulation time remainingSimulationTime -= dt; // update t time = time + dt; } // interpolate state double alpha = remainingSimulationTime / dt; // set interpolated state for (int i = 0; i < y.length; i++) { interpolatedY[i] = y[i] * alpha + yPrev[i] * (1.0 - alpha); } // update properties for visualization updateView(nodes, interpolatedY, m, ignoreFlags); } }; // finally, start the framle listener frameListener.start(); }
From source file:de.uni_erlangen.lstm.modelaccess.Model.java
/** * Run the model using set parameters/*from w w w . j ava 2 s . c o m*/ */ public void simulate() { finished = false; /* * Integrator selection */ //FirstOrderIntegrator integrator = new HighamHall54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); //FirstOrderIntegrator integrator = new DormandPrince54Integrator(1.0e-12, 100.0, 1.0e-12, 1.0e-12); //FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); //FirstOrderIntegrator integrator = new GraggBulirschStoerIntegrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); FirstOrderIntegrator integrator = new AdamsBashforthIntegrator(2, 1.0e-14, 100.0, 1.0e-10, 1.0e-10); //FirstOrderIntegrator integrator = new AdamsMoultonIntegrator(2, 1.0e-8, 100.0, 1.0e-10, 1.0e-10); // influent values, digester parameters, S_H_ion, dae system final DAEModel ode = new DAEModel(u, param, S_H_ion, dae, fix_pH); //FirstOrderDifferentialEquations ode = model; // Records progress StepHandler progHandler = new StepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(StepInterpolator interpolator, boolean isLast) { progress = interpolator.getCurrentTime(); } }; integrator.addStepHandler(progHandler); /* * Continuous model recorded in CSV */ if (onlineRecord) { final CSVWriter writer = new CSVWriter(); StepHandler stepHandler = new StepHandler() { double prevT = 0.0; public void init(double t0, double[] y0, double t) { } public void handleStep(StepInterpolator interpolator, boolean isLast) { double t = interpolator.getCurrentTime(); if (t - prevT > resolution) { // Add time to the beginning of the array double[] timemodel = new double[ode.getDimensions().length + 1]; timemodel[0] = t; // We need to pull variables (S_h2 and acid-base) directly from the model if using DAE for (int i = 1; i < timemodel.length; i++) { timemodel[i] = ode.getDimensions()[i - 1]; } writer.WriteArray(output_file, timemodel, true); prevT = t; } } }; integrator.addStepHandler(stepHandler); } /* * Add event handlers for discrete events * maxCheck - maximal time interval between switching function checks (this interval prevents missing sign changes in case the integration steps becomes very large) * conv - convergence threshold in the event time search * maxIt - upper limit of the iteration count in the event time search */ if (events.size() > 0) { for (DiscreteEvent event : events) { double maxCheck = Double.POSITIVE_INFINITY; double conv = 1.0e-20; int maxIt = 100; integrator.addEventHandler(event, maxCheck, conv, maxIt); } } integrator.integrate(ode, start, x, end, x); /* * Return the time that the discrete event occurred */ if (events.size() > 0) { for (DiscreteEvent event : events) { if (event.getTime() < end) { end = event.getTime(); } } } // We need to pull variables (S_h2 and acid-base) directly from the model x = ode.getDimensions(); finished = true; }
From source file:nl.rivm.cib.episim.model.disease.infection.MSEIRSTest.java
public static Observable<Map.Entry<Double, double[]>> deterministic(final SIRConfig config, final Supplier<FirstOrderIntegrator> integrators) { return Observable.create(sub -> { final double gamma = 1. / config.recovery(); final double beta = gamma * config.reproduction(); final double[] y0 = Arrays.stream(config.population()).mapToDouble(n -> n).toArray(); final double[] t = config.t(); try {//from w w w . j ava2 s . c o m final FirstOrderIntegrator integrator = integrators.get(); integrator.addStepHandler(new StepHandler() { @Override public void init(final double t0, final double[] y0, final double t) { publishCopy(sub, t0, y0); } @Override public void handleStep(final StepInterpolator interpolator, final boolean isLast) throws MaxCountExceededException { publishCopy(sub, interpolator.getInterpolatedTime(), interpolator.getInterpolatedState()); if (isLast) sub.onComplete(); } }); integrator.integrate(new FirstOrderDifferentialEquations() { @Override public int getDimension() { return y0.length; } @Override public void computeDerivatives(final double t, final double[] y, final double[] yp) { // SIR terms (flow rates) final double n = y[0] + y[1] + y[2], flow_si = beta * y[0] * y[1] / n, flow_ir = gamma * y[1]; yp[0] = -flow_si; yp[1] = flow_si - flow_ir; yp[2] = flow_ir; } }, t[0], y0, t[1], y0); } catch (final Exception e) { sub.onError(e); } }); }
From source file:org.orekit.utils.AngularCoordinatesTest.java
@Test public void testShiftWithAcceleration() throws OrekitException { double rate = 2 * FastMath.PI / (12 * 60); double acc = 0.001; double dt = 1.0; int n = 2000; final AngularCoordinates quadratic = new AngularCoordinates(Rotation.IDENTITY, new Vector3D(rate, Vector3D.PLUS_K), new Vector3D(acc, Vector3D.PLUS_J)); final AngularCoordinates linear = new AngularCoordinates(quadratic.getRotation(), quadratic.getRotationRate(), Vector3D.ZERO); final FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { public int getDimension() { return 4; }/* w w w. java 2 s.c o m*/ public void computeDerivatives(final double t, final double[] q, final double[] qDot) { final double omegaX = quadratic.getRotationRate().getX() + t * quadratic.getRotationAcceleration().getX(); final double omegaY = quadratic.getRotationRate().getY() + t * quadratic.getRotationAcceleration().getY(); final double omegaZ = quadratic.getRotationRate().getZ() + t * quadratic.getRotationAcceleration().getZ(); qDot[0] = 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ); qDot[1] = 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ); qDot[2] = 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ); qDot[3] = 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ); } }; FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12); integrator.addStepHandler(new StepNormalizer(dt / n, new FixedStepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(double t, double[] y, double[] yDot, boolean isLast) { Rotation reference = new Rotation(y[0], y[1], y[2], y[3], true); // the error in shiftedBy taking acceleration into account is cubic double expectedCubicError = 1.4544e-6 * t * t * t; Assert.assertEquals(expectedCubicError, Rotation.distance(reference, quadratic.shiftedBy(t).getRotation()), 0.0001 * expectedCubicError); // the error in shiftedBy not taking acceleration into account is quadratic double expectedQuadraticError = 5.0e-4 * t * t; Assert.assertEquals(expectedQuadraticError, Rotation.distance(reference, linear.shiftedBy(t).getRotation()), 0.00001 * expectedQuadraticError); } })); double[] y = new double[] { quadratic.getRotation().getQ0(), quadratic.getRotation().getQ1(), quadratic.getRotation().getQ2(), quadratic.getRotation().getQ3() }; integrator.integrate(ode, 0, y, dt, y); }
From source file:org.orekit.utils.TimeStampedAngularCoordinatesTest.java
private double[] interpolationErrors(final TimeStampedAngularCoordinates reference, double dt) throws OrekitException { final FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { public int getDimension() { return 4; }/*from w w w. ja va 2 s. c om*/ public void computeDerivatives(final double t, final double[] q, final double[] qDot) { final double omegaX = reference.getRotationRate().getX() + t * reference.getRotationAcceleration().getX(); final double omegaY = reference.getRotationRate().getY() + t * reference.getRotationAcceleration().getY(); final double omegaZ = reference.getRotationRate().getZ() + t * reference.getRotationAcceleration().getZ(); qDot[0] = 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ); qDot[1] = 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ); qDot[2] = 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ); qDot[3] = 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ); } }; final List<TimeStampedAngularCoordinates> complete = new ArrayList<TimeStampedAngularCoordinates>(); FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12); integrator.addStepHandler(new StepNormalizer(dt / 2000, new FixedStepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(double t, double[] y, double[] yDot, boolean isLast) { complete.add( new TimeStampedAngularCoordinates(reference.getDate().shiftedBy(t), new Rotation(y[0], y[1], y[2], y[3], true), new Vector3D(1, reference.getRotationRate(), t, reference.getRotationAcceleration()), reference.getRotationAcceleration())); } })); double[] y = new double[] { reference.getRotation().getQ0(), reference.getRotation().getQ1(), reference.getRotation().getQ2(), reference.getRotation().getQ3() }; integrator.integrate(ode, 0, y, dt, y); List<TimeStampedAngularCoordinates> sample = new ArrayList<TimeStampedAngularCoordinates>(); sample.add(complete.get(0)); sample.add(complete.get(complete.size() / 2)); sample.add(complete.get(complete.size() - 1)); double maxRotationError = 0; double maxRateError = 0; double maxAccelerationError = 0; for (TimeStampedAngularCoordinates acRef : complete) { TimeStampedAngularCoordinates interpolated = TimeStampedAngularCoordinates.interpolate(acRef.getDate(), AngularDerivativesFilter.USE_RRA, sample); double rotationError = Rotation.distance(acRef.getRotation(), interpolated.getRotation()); double rateError = Vector3D.distance(acRef.getRotationRate(), interpolated.getRotationRate()); double accelerationError = Vector3D.distance(acRef.getRotationAcceleration(), interpolated.getRotationAcceleration()); maxRotationError = FastMath.max(maxRotationError, rotationError); maxRateError = FastMath.max(maxRateError, rateError); maxAccelerationError = FastMath.max(maxAccelerationError, accelerationError); } return new double[] { maxRotationError, maxRateError, maxAccelerationError }; }
From source file:reactor.semibatchreactor.Simulation.java
public static void main(String[] args) { double[] initialConditions = { 10, 0, 0, 0, 300 }; double[] tempParameters = { 2500, 280, 2, 0.004, 75240, 300, -7.9076 * Math.pow(10, 7) }; double[] odeParameters = { 0.2, 6.14, 0, 2.37989273, 8.94266 * Math.pow(10, 12), 803373.6 }; SemiBatchReactor reactor = new SemiBatchReactor(initialConditions, tempParameters, odeParameters); FirstOrderIntegrator integrator = new DormandPrince853Integrator(0.1, 0.1, 0.001, 0.001); FirstOrderDifferentialEquations ode = reactor; integrator.addStepHandler(reactor.stepHandler); integrator.integrate(ode, 0.0, initialConditions, 10, initialConditions); }