| /** |
| * <copyright> |
| * </copyright> |
| * |
| * $Id$ |
| */ |
| package org.eclipse.stem.solvers.rk.impl; |
| |
| import java.util.ArrayList; |
| import java.util.HashMap; |
| import java.util.Iterator; |
| |
| import org.apache.commons.math.ode.DerivativeException; |
| import org.apache.commons.math.ode.FirstOrderDifferentialEquations; |
| import org.apache.commons.math.ode.FirstOrderIntegrator; |
| import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator; |
| import org.eclipse.core.runtime.Preferences; |
| import org.eclipse.emf.common.notify.Notification; |
| import org.eclipse.emf.common.util.BasicEList; |
| import org.eclipse.emf.common.util.EList; |
| |
| import org.eclipse.emf.ecore.EAttribute; |
| import org.eclipse.emf.ecore.EClass; |
| |
| import org.eclipse.emf.ecore.impl.ENotificationImpl; |
| |
| import org.eclipse.stem.core.graph.DynamicLabel; |
| import org.eclipse.stem.core.graph.IntegrationLabel; |
| import org.eclipse.stem.core.graph.IntegrationLabelValue; |
| import org.eclipse.stem.core.graph.LabelValue; |
| import org.eclipse.stem.core.model.Decorator; |
| import org.eclipse.stem.core.model.IntegrationDecorator; |
| import org.eclipse.stem.core.model.STEMTime; |
| import org.eclipse.stem.core.model.TransformationDecorator; |
| import org.eclipse.stem.core.solver.impl.SolverImpl; |
| import org.eclipse.stem.core.trigger.Trigger; |
| |
| import org.eclipse.stem.solvers.rk.DormandPrince853; |
| import org.eclipse.stem.solvers.rk.RkPackage; |
| import org.eclipse.stem.ui.Activator; |
| |
| /** |
| * <!-- begin-user-doc --> |
| * An implementation of the model object '<em><b>Dormand Prince853</b></em>'. |
| * <!-- end-user-doc --> |
| * <p> |
| * The following features are implemented: |
| * <ul> |
| * <li>{@link org.eclipse.stem.solvers.rk.impl.DormandPrince853Impl#getRelativeTolerance <em>Relative Tolerance</em>}</li> |
| * <li>{@link org.eclipse.stem.solvers.rk.impl.DormandPrince853Impl#getAbsoluteTolerance <em>Absolute Tolerance</em>}</li> |
| * <li>{@link org.eclipse.stem.solvers.rk.impl.DormandPrince853Impl#getMinStep <em>Min Step</em>}</li> |
| * <li>{@link org.eclipse.stem.solvers.rk.impl.DormandPrince853Impl#getMaxStep <em>Max Step</em>}</li> |
| * </ul> |
| * </p> |
| * |
| * @generated |
| */ |
| public class DormandPrince853Impl extends SolverImpl implements DormandPrince853 { |
| /** |
| * The default value of the '{@link #getRelativeTolerance() <em>Relative Tolerance</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getRelativeTolerance() |
| * @generated |
| * @ordered |
| */ |
| protected static final double RELATIVE_TOLERANCE_EDEFAULT = 1.0E-9; |
| |
| /** |
| * The cached value of the '{@link #getRelativeTolerance() <em>Relative Tolerance</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getRelativeTolerance() |
| * @generated |
| * @ordered |
| */ |
| protected double relativeTolerance = RELATIVE_TOLERANCE_EDEFAULT; |
| |
| /** |
| * The default value of the '{@link #getAbsoluteTolerance() <em>Absolute Tolerance</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getAbsoluteTolerance() |
| * @generated |
| * @ordered |
| */ |
| protected static final double ABSOLUTE_TOLERANCE_EDEFAULT = 1.0E-5; |
| |
| /** |
| * The cached value of the '{@link #getAbsoluteTolerance() <em>Absolute Tolerance</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getAbsoluteTolerance() |
| * @generated |
| * @ordered |
| */ |
| protected double absoluteTolerance = ABSOLUTE_TOLERANCE_EDEFAULT; |
| |
| /** |
| * The default value of the '{@link #getMinStep() <em>Min Step</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getMinStep() |
| * @generated |
| * @ordered |
| */ |
| protected static final double MIN_STEP_EDEFAULT = 1.0E-8; |
| |
| /** |
| * The cached value of the '{@link #getMinStep() <em>Min Step</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getMinStep() |
| * @generated |
| * @ordered |
| */ |
| protected double minStep = MIN_STEP_EDEFAULT; |
| |
| /** |
| * The default value of the '{@link #getMaxStep() <em>Max Step</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getMaxStep() |
| * @generated |
| * @ordered |
| */ |
| protected static final double MAX_STEP_EDEFAULT = 1.0; |
| |
| /** |
| * The cached value of the '{@link #getMaxStep() <em>Max Step</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getMaxStep() |
| * @generated |
| * @ordered |
| */ |
| protected double maxStep = MAX_STEP_EDEFAULT; |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated NOT |
| */ |
| public DormandPrince853Impl() { |
| super(); |
| } |
| |
| // Cache to improve performance |
| private HashMap<Integer, double []> arraysBySize = new HashMap<Integer, double[]>(); |
| |
| /** |
| * Step one |
| */ |
| @Override |
| public boolean step(STEMTime time, long timeDelta, int cycle) { |
| |
| EList<TransformationDecorator> transDecorators = new BasicEList<TransformationDecorator>(); |
| // Validate all decorators that return deltas to make sure |
| // they are of deterministic nature. |
| |
| for(Decorator decorator:this.getCanonicalGraph().getDecorators()){ |
| if(decorator instanceof IntegrationDecorator) { |
| IntegrationDecorator idec = (IntegrationDecorator)decorator; |
| if(!idec.isDeterministic()) { |
| Activator.logError("Error, decorator: "+idec+" is not deterministic. The Runge Kutta Integrator can only handle deterministic models.", new Exception()); |
| return false; |
| } |
| } |
| if(decorator instanceof TransformationDecorator) transDecorators.add((TransformationDecorator)decorator); |
| } |
| |
| short num_threads; |
| Activator act = org.eclipse.stem.ui.Activator.getDefault(); |
| if(act != null) { |
| final Preferences preferences = act.getPluginPreferences(); |
| num_threads = (short)preferences.getInt(org.eclipse.stem.ui.preferences.PreferenceConstants.SIMULATION_THREADS); |
| } else num_threads = 2; // Just so we can run inside junit test |
| |
| partitioner.setNumProcesses(num_threads); |
| |
| // Find triggers and make sure they are invoked |
| for(Decorator decorator:this.getCanonicalGraph().getDecorators()) { |
| if(decorator instanceof Trigger || decorator instanceof TransformationDecorator) { |
| decorator.updateLabels(time, timeDelta, cycle); |
| } |
| } |
| |
| // First initialize the Y and temp label values from the current |
| // label values. |
| |
| for(Decorator decorator:this.getCanonicalGraph().getDecorators()) { |
| if(decorator instanceof IntegrationDecorator) { |
| EList<DynamicLabel>allLabels = decorator.getLabelsToUpdate(); |
| for (final Iterator<DynamicLabel> currentStateLabelIter = allLabels.iterator(); currentStateLabelIter.hasNext();) { |
| // It's a standard disease model with a standard disease model label |
| final IntegrationLabel iLabel = (IntegrationLabel) currentStateLabelIter.next(); |
| ((IntegrationLabelValue)iLabel.getProbeValue()).set((IntegrationLabelValue)iLabel.getCurrentValue()); |
| ((IntegrationLabelValue)iLabel.getTempValue()).set((IntegrationLabelValue)iLabel.getCurrentValue()); |
| ((IntegrationLabelValue)iLabel.getTempValue()).prepareCycle(); |
| ((IntegrationLabelValue)iLabel.getProbeValue()).prepareCycle(); |
| } |
| } |
| } |
| |
| FirstOrderIntegrator dp853 = new DormandPrince853Integrator(this.getMinStep(), this.getMaxStep(), this.getAbsoluteTolerance(), this.getRelativeTolerance()); |
| |
| ArrayList<IntegrationDecorator> iDecorators = new ArrayList<IntegrationDecorator>(); |
| |
| for(Decorator d:getCanonicalGraph().getDecorators()) { |
| if(d instanceof IntegrationDecorator) |
| iDecorators.add((IntegrationDecorator)d); |
| } |
| |
| // The initial value y0 will be a (potentially very large) array containing the |
| // values for each decorator and each patch. Figure out how large the array needs to be |
| |
| double [] y0; |
| int size=0; |
| for(IntegrationDecorator sdm:iDecorators) |
| for(int i=0;i<num_threads;++i) { |
| // We need to use the partitioner to make sure this works in distributed STEM |
| EList<DynamicLabel> list = partitioner.partitionDecoratorLabels(sdm, i); |
| int sz=0; |
| // We need to check all since it's not guaranteed that the labels a decorator |
| // is updating all are the same size |
| for(int j=0;j<list.size();++j) |
| size += list.get(j).getCurrentValue().eClass().getEAllAttributes().size(); |
| } |
| |
| if(arraysBySize.containsKey(size)) |
| y0 = arraysBySize.get(size); |
| else { |
| y0 = new double[size]; |
| arraysBySize.put(size, y0); |
| } |
| |
| |
| int index = 0; |
| for(IntegrationDecorator sdm:iDecorators) |
| for(int i=0;i<num_threads;++i) { |
| Iterator<DynamicLabel> iter = partitioner.partitionDecoratorLabels(sdm, i) |
| .iterator(); |
| while(iter.hasNext()) { |
| DynamicLabel lab = iter.next(); |
| int sz = lab.getCurrentValue().eClass().getEAllAttributes().size(); |
| double []ytmp = null; |
| if(arraysBySize.containsKey(sz)) // Could be same array as y0, which is okay |
| ytmp = arraysBySize.get(sz); |
| else { |
| ytmp = new double[sz]; |
| arraysBySize.put(sz, ytmp); |
| } |
| |
| // We must use the temp value. It was copied from the current value above |
| // but it's been prepared ensuring that variables such as incidence has been reset. |
| getValues(((IntegrationLabel)lab).getTempValue(), ytmp); |
| for(int j=0;j<ytmp.length;++j) |
| y0[index++]=ytmp[j]; |
| } |
| } |
| double t0 = (double)cycle; |
| STEMDiffEquation sde = new STEMDiffEquation(time, timeDelta, iDecorators, num_threads, size); |
| try { |
| dp853.integrate(sde, t0, y0, t0+1.0, y0); |
| } catch(Exception e) { |
| Activator.logError(e.getMessage(), e); |
| return false; // failed, will force simulation to stop |
| } |
| |
| |
| |
| index = 0; |
| for(IntegrationDecorator sdm:iDecorators) |
| for(int i=0;i<num_threads;++i) { |
| EList<DynamicLabel> labels = partitioner.partitionDecoratorLabels(sdm, i); |
| for(DynamicLabel lab:labels) { |
| int sz = lab.getCurrentValue().eClass().getEAllAttributes().size(); |
| double []ytmp =null; |
| if(arraysBySize.containsKey(sz)) |
| ytmp = arraysBySize.get(sz); |
| else { |
| ytmp = new double[sz]; |
| arraysBySize.put(sz, ytmp); |
| } |
| for(int j=0;j<ytmp.length;++j) |
| ytmp[j] = y0[index++]; |
| setValues((IntegrationLabelValue)lab.getNextValue(), ytmp); |
| |
| } |
| } |
| |
| return true; |
| } |
| |
| void getValues(IntegrationLabelValue lv, double []result) { |
| // Assume order is always same here |
| int i=0; |
| EList<EAttribute>attrList = lv.eClass().getEAllAttributes(); |
| for(int j=0;j<attrList.size();++j) { |
| EAttribute ea = attrList.get(j); |
| result[i++] = lv.get(ea); |
| } |
| } |
| |
| private void setValues(IntegrationLabelValue ilv, double [] d) { |
| int i=0; |
| EList<EAttribute>attrList = ilv.eClass().getEAllAttributes(); |
| for(int j=0;j<attrList.size();++j) { |
| EAttribute ea = attrList.get(j); |
| if(ea.isChangeable()) ilv.set(ea, d[i++]); |
| else i++; |
| } |
| } |
| |
| public class STEMDiffEquation implements FirstOrderDifferentialEquations { |
| STEMTime time; |
| long timeDelta; |
| EList<DynamicLabel> labelList = new BasicEList<DynamicLabel>(); |
| ArrayList<IntegrationDecorator> decorators; |
| int num_threads; |
| int dimensions; |
| |
| public STEMDiffEquation(STEMTime t, long timeDelta, ArrayList<IntegrationDecorator>decs, int numThreads, int dimensions) { |
| this.time = t; |
| this.timeDelta = timeDelta; |
| decorators = decs; |
| num_threads = numThreads; |
| this.dimensions = dimensions; |
| } |
| |
| public int getDimension() { |
| return dimensions; |
| } |
| |
| public void computeDerivatives(double t, double[] y, double[] yDot) |
| throws DerivativeException { |
| |
| int index = 0; |
| double [] ytemp = null; |
| for(int i=0;i<decorators.size();++i) { |
| IntegrationDecorator sdm = decorators.get(i); |
| for(int j=0;j<num_threads;++j) { |
| EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, j); |
| for(int k=0;k<llist.size();++k) { |
| DynamicLabel l = llist.get(k); |
| // Set both probe and temp. When (if) we decide to handle multi-threaded code, |
| // we need to revisit this and figure out how to set the two values. Temp |
| // is supposed to be safe to read by another concurrent thread. |
| IntegrationLabelValue probeV = ((IntegrationLabel)l).getProbeValue(); |
| IntegrationLabelValue tempV = ((IntegrationLabel)l).getTempValue(); |
| int dim = probeV.eClass().getEAllAttributes().size(); |
| if(arraysBySize.containsKey(dim)) |
| ytemp = arraysBySize.get(dim); |
| else { |
| ytemp = new double[dim]; |
| arraysBySize.put(dim, ytemp); |
| } |
| |
| for(int m=0;m<dim;++m) |
| ytemp[m] = y[index++]; |
| |
| DormandPrince853Impl.this.setValues(probeV, ytemp); |
| DormandPrince853Impl.this.setValues(tempV, ytemp); |
| } |
| } |
| } |
| |
| // Do the magic |
| for(int i=0;i<decorators.size();++i) { |
| IntegrationDecorator sdm = decorators.get(i); |
| for(int j=0;j<num_threads;++j) { |
| EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, j); |
| sdm.calculateDelta(time, t, timeDelta, labels); |
| } |
| } |
| for(int i=0;i<decorators.size();++i) { |
| IntegrationDecorator sdm = decorators.get(i); |
| for(int j=0;j<num_threads;++j) { |
| EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, j); |
| sdm.applyExternalDeltas(time, t, timeDelta, labels); |
| } |
| } |
| |
| // Copy over the deltas to the output |
| |
| index = 0; |
| for(int i=0;i<decorators.size();++i) { |
| IntegrationDecorator sdm = decorators.get(i); |
| for(int j=0;j<num_threads;++j) { |
| EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, j); |
| for(DynamicLabel l:llist) { |
| IntegrationLabelValue deltaV = ((IntegrationLabel)l).getDeltaValue(); |
| int dim = deltaV.eClass().getEAllAttributes().size(); |
| if(arraysBySize.containsKey(dim)) |
| ytemp = arraysBySize.get(dim); |
| else { |
| ytemp = new double[dim]; |
| arraysBySize.put(dim, ytemp); |
| } |
| |
| DormandPrince853Impl.this.getValues(((IntegrationLabel)l).getDeltaValue(), ytemp); |
| |
| for(int k=0;k<dim;++k) |
| yDot[index++] = ytemp[k]; |
| } |
| } |
| } |
| } |
| } |
| |
| /** |
| * Reset the solver |
| * @generated NOT |
| */ |
| @Override |
| public void reset() { |
| this.setInitialized(false); |
| } |
| |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| protected EClass eStaticClass() { |
| return RkPackage.Literals.DORMAND_PRINCE853; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getRelativeTolerance() { |
| return relativeTolerance; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setRelativeTolerance(double newRelativeTolerance) { |
| double oldRelativeTolerance = relativeTolerance; |
| relativeTolerance = newRelativeTolerance; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, RkPackage.DORMAND_PRINCE853__RELATIVE_TOLERANCE, oldRelativeTolerance, relativeTolerance)); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getAbsoluteTolerance() { |
| return absoluteTolerance; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setAbsoluteTolerance(double newAbsoluteTolerance) { |
| double oldAbsoluteTolerance = absoluteTolerance; |
| absoluteTolerance = newAbsoluteTolerance; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, RkPackage.DORMAND_PRINCE853__ABSOLUTE_TOLERANCE, oldAbsoluteTolerance, absoluteTolerance)); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getMinStep() { |
| return minStep; |
| } |
| |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setMinStep(double newMinStep) { |
| double oldMinStep = minStep; |
| minStep = newMinStep; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, RkPackage.DORMAND_PRINCE853__MIN_STEP, oldMinStep, minStep)); |
| } |
| |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getMaxStep() { |
| return maxStep; |
| } |
| |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setMaxStep(double newMaxStep) { |
| double oldMaxStep = maxStep; |
| maxStep = newMaxStep; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, RkPackage.DORMAND_PRINCE853__MAX_STEP, oldMaxStep, maxStep)); |
| } |
| |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public Object eGet(int featureID, boolean resolve, boolean coreType) { |
| switch (featureID) { |
| case RkPackage.DORMAND_PRINCE853__RELATIVE_TOLERANCE: |
| return getRelativeTolerance(); |
| case RkPackage.DORMAND_PRINCE853__ABSOLUTE_TOLERANCE: |
| return getAbsoluteTolerance(); |
| case RkPackage.DORMAND_PRINCE853__MIN_STEP: |
| return getMinStep(); |
| case RkPackage.DORMAND_PRINCE853__MAX_STEP: |
| return getMaxStep(); |
| } |
| return super.eGet(featureID, resolve, coreType); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public void eSet(int featureID, Object newValue) { |
| switch (featureID) { |
| case RkPackage.DORMAND_PRINCE853__RELATIVE_TOLERANCE: |
| setRelativeTolerance((Double)newValue); |
| return; |
| case RkPackage.DORMAND_PRINCE853__ABSOLUTE_TOLERANCE: |
| setAbsoluteTolerance((Double)newValue); |
| return; |
| case RkPackage.DORMAND_PRINCE853__MIN_STEP: |
| setMinStep((Double)newValue); |
| return; |
| case RkPackage.DORMAND_PRINCE853__MAX_STEP: |
| setMaxStep((Double)newValue); |
| return; |
| } |
| super.eSet(featureID, newValue); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public void eUnset(int featureID) { |
| switch (featureID) { |
| case RkPackage.DORMAND_PRINCE853__RELATIVE_TOLERANCE: |
| setRelativeTolerance(RELATIVE_TOLERANCE_EDEFAULT); |
| return; |
| case RkPackage.DORMAND_PRINCE853__ABSOLUTE_TOLERANCE: |
| setAbsoluteTolerance(ABSOLUTE_TOLERANCE_EDEFAULT); |
| return; |
| case RkPackage.DORMAND_PRINCE853__MIN_STEP: |
| setMinStep(MIN_STEP_EDEFAULT); |
| return; |
| case RkPackage.DORMAND_PRINCE853__MAX_STEP: |
| setMaxStep(MAX_STEP_EDEFAULT); |
| return; |
| } |
| super.eUnset(featureID); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public boolean eIsSet(int featureID) { |
| switch (featureID) { |
| case RkPackage.DORMAND_PRINCE853__RELATIVE_TOLERANCE: |
| return relativeTolerance != RELATIVE_TOLERANCE_EDEFAULT; |
| case RkPackage.DORMAND_PRINCE853__ABSOLUTE_TOLERANCE: |
| return absoluteTolerance != ABSOLUTE_TOLERANCE_EDEFAULT; |
| case RkPackage.DORMAND_PRINCE853__MIN_STEP: |
| return minStep != MIN_STEP_EDEFAULT; |
| case RkPackage.DORMAND_PRINCE853__MAX_STEP: |
| return maxStep != MAX_STEP_EDEFAULT; |
| } |
| return super.eIsSet(featureID); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public String toString() { |
| if (eIsProxy()) return super.toString(); |
| |
| StringBuffer result = new StringBuffer(super.toString()); |
| result.append(" (relativeTolerance: "); |
| result.append(relativeTolerance); |
| result.append(", absoluteTolerance: "); |
| result.append(absoluteTolerance); |
| result.append(", minStep: "); |
| result.append(minStep); |
| result.append(", maxStep: "); |
| result.append(maxStep); |
| result.append(')'); |
| return result.toString(); |
| } |
| |
| } //DormandPrince853Impl |