List of usage examples for org.apache.commons.math3.ode FirstOrderDifferentialEquations getDimension
int getDimension();
From source file:beast.structuredCoalescent.distribution.IndependentStructuredCoalescent.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 w w w . j a v a 2 s.co 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 independent_ode_integrator(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) { 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); } return logP; }
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();/* w ww . j a va 2 s .co 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();//ww w. j ava2s . 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; }