core/org.eclipse.stem.solvers.rk/src/org/eclipse/stem/solvers/rk/impl/DormandPrince853Impl.java - stem/org.eclipse.stem - Git at Google

 /**
  * <copyright>
  * </copyright>
  *
  * $Id$
  */
 package org.eclipse.stem.solvers.rk.impl;

 import java.util.ArrayList;
 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.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();
 	}


 	/**
 	 * Step one
 	 */
 	@Override
 	public boolean step(STEMTime time, long timeDelta, int cycle) {

 		// Validate all decorators that return deltas to make sure
 		// they are of deterministic nature.

 		for(Decorator decorator:this.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;
 				}
 			}

 		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.getDecorators()) {
 			if(decorator instanceof Trigger) {
 				decorator.updateLabels(time, timeDelta, cycle);
 			}
 		}

 		// First initialize the Y and temp label values from the current
 		// label values.

 		for(Decorator decorator:this.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: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();
 			}

 		y0 = new double[size];

 		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();
 					double []ytmp = new double[lab.getCurrentValue().eClass().getEAllAttributes().size()];
 					// 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) {
 					double []ytmp = new double[lab.getCurrentValue().eClass().getEAllAttributes().size()];
 					for(int j=0;j<ytmp.length;++j)
 						ytmp[j] = y0[index++];
 					setValues((IntegrationLabelValue)lab.getNextValue(), ytmp);

 				}
 			}

 		return true;
 	}

 	void getValues(LabelValue lv, double []result) {
 		// Assume order is always same here
 		int i=0;
 		for(EAttribute ea:lv.eClass().getEAllAttributes()) {
 			result[i++] = (Double)lv.eGet(ea);
 		}
 	}

 	private void setValues(IntegrationLabelValue ilv, double [] d) {
 		int i=0;
 		for(EAttribute ea:ilv.eClass().getEAllAttributes())
 			if(ea.isChangeable()) ilv.eSet(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(IntegrationDecorator sdm:decorators)
 				for(int i=0;i<num_threads;++i) {
 					EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
 					for(DynamicLabel l:llist) {
 						// 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(ytemp == null || ytemp.length != dim)
 							ytemp = new double[dim];

 						for(int j=0;j<dim;++j)
 							ytemp[j] = y[index++];

 						DormandPrince853Impl.this.setValues(probeV, ytemp);
 						DormandPrince853Impl.this.setValues(tempV, ytemp);
 					}
 				}

 			// Do the magic
 			for(IntegrationDecorator sdm:decorators)
 				for(int i=0;i<num_threads;++i) {
 					EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
 					sdm.calculateDelta(time, timeDelta, labels);
 				}
 			for(IntegrationDecorator sdm:decorators)
 				for(int i=0;i<num_threads;++i) {
 					EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
 					sdm.applyExternalDeltas(time, timeDelta, labels);
 				}

 			// Copy over the deltas to the output

 			index = 0;
 			for(IntegrationDecorator sdm:decorators)
 				for(int i=0;i<num_threads;++i) {
 					EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
 					for(DynamicLabel l:llist) {
 						IntegrationLabelValue deltaV = ((IntegrationLabel)l).getDeltaValue();
 						int dim =  deltaV.eClass().getEAllAttributes().size();
 						if(ytemp == null || ytemp.length != dim)
 							ytemp = new double[dim];

 						DormandPrince853Impl.this.getValues(((IntegrationLabel)l).getDeltaValue(), ytemp);

 						for(int j=0;j<dim;++j)
 							yDot[index++] = ytemp[j];
 					}
 				}
 		}
 	}

 	/**
 	 * 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
	/**
	* <copyright>
	* </copyright>
	*
	* $Id$
	*/
	package org.eclipse.stem.solvers.rk.impl;

	import java.util.ArrayList;
	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.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();
	}


	/**
	* Step one
	*/
	@Override
	public boolean step(STEMTime time, long timeDelta, int cycle) {

	// Validate all decorators that return deltas to make sure
	// they are of deterministic nature.

	for(Decorator decorator:this.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;
	}
	}

	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.getDecorators()) {
	if(decorator instanceof Trigger) {
	decorator.updateLabels(time, timeDelta, cycle);
	}
	}

	// First initialize the Y and temp label values from the current
	// label values.

	for(Decorator decorator:this.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: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();
	}

	y0 = new double[size];

	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();
	double []ytmp = new double[lab.getCurrentValue().eClass().getEAllAttributes().size()];
	// 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) {
	double []ytmp = new double[lab.getCurrentValue().eClass().getEAllAttributes().size()];
	for(int j=0;j<ytmp.length;++j)
	ytmp[j] = y0[index++];
	setValues((IntegrationLabelValue)lab.getNextValue(), ytmp);

	}
	}

	return true;
	}

	void getValues(LabelValue lv, double []result) {
	// Assume order is always same here
	int i=0;
	for(EAttribute ea:lv.eClass().getEAllAttributes()) {
	result[i++] = (Double)lv.eGet(ea);
	}
	}

	private void setValues(IntegrationLabelValue ilv, double [] d) {
	int i=0;
	for(EAttribute ea:ilv.eClass().getEAllAttributes())
	if(ea.isChangeable()) ilv.eSet(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(IntegrationDecorator sdm:decorators)
	for(int i=0;i<num_threads;++i) {
	EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
	for(DynamicLabel l:llist) {
	// 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(ytemp == null \|\| ytemp.length != dim)
	ytemp = new double[dim];

	for(int j=0;j<dim;++j)
	ytemp[j] = y[index++];

	DormandPrince853Impl.this.setValues(probeV, ytemp);
	DormandPrince853Impl.this.setValues(tempV, ytemp);
	}
	}

	// Do the magic
	for(IntegrationDecorator sdm:decorators)
	for(int i=0;i<num_threads;++i) {
	EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
	sdm.calculateDelta(time, timeDelta, labels);
	}
	for(IntegrationDecorator sdm:decorators)
	for(int i=0;i<num_threads;++i) {
	EList<DynamicLabel> labels = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
	sdm.applyExternalDeltas(time, timeDelta, labels);
	}

	// Copy over the deltas to the output

	index = 0;
	for(IntegrationDecorator sdm:decorators)
	for(int i=0;i<num_threads;++i) {
	EList<DynamicLabel> llist = DormandPrince853Impl.this.partitioner.partitionDecoratorLabels(sdm, i);
	for(DynamicLabel l:llist) {
	IntegrationLabelValue deltaV = ((IntegrationLabel)l).getDeltaValue();
	int dim = deltaV.eClass().getEAllAttributes().size();
	if(ytemp == null \|\| ytemp.length != dim)
	ytemp = new double[dim];

	DormandPrince853Impl.this.getValues(((IntegrationLabel)l).getDeltaValue(), ytemp);

	for(int j=0;j<dim;++j)
	yDot[index++] = ytemp[j];
	}
	}
	}
	}

	/**
	* 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