| package org.eclipse.stem.diseasemodels.forcing.impl; |
| /******************************************************************************* |
| * Copyright (c) 2010 IBM Corporation and others. |
| * All rights reserved. This program and the accompanying materials |
| * are made available under the terms of the Eclipse Public License v1.0 |
| * which accompanies this distribution, and is available at |
| * http://www.eclipse.org/legal/epl-v10.html |
| * |
| * Contributors: |
| * IBM Corporation - initial API and implementation |
| *******************************************************************************/ |
| |
| import java.io.FileWriter; |
| import java.util.ArrayList; |
| import java.util.Calendar; |
| |
| import org.eclipse.emf.common.notify.Notification; |
| import org.eclipse.emf.ecore.EClass; |
| import org.eclipse.emf.ecore.impl.ENotificationImpl; |
| import org.eclipse.stem.core.model.STEMTime; |
| import org.eclipse.stem.diseasemodels.forcing.ForcingPackage; |
| import org.eclipse.stem.diseasemodels.forcing.GaussianForcingDiseaseModel; |
| import org.eclipse.stem.diseasemodels.standard.DiseaseModelLabelValue; |
| import org.eclipse.stem.diseasemodels.standard.SILabelValue; |
| import org.eclipse.stem.diseasemodels.standard.SIRLabelValue; |
| import org.eclipse.stem.diseasemodels.standard.StandardDiseaseModelLabel; |
| import org.eclipse.stem.diseasemodels.standard.StandardDiseaseModelLabelValue; |
| import org.eclipse.stem.diseasemodels.standard.impl.SIRLabelValueImpl; |
| import org.eclipse.stem.diseasemodels.standard.impl.StochasticSIRDiseaseModelImpl; |
| |
| /** |
| * <!-- begin-user-doc --> |
| * An implementation of the model object '<em><b>Gaussian Forcing Disease Model</b></em>'. |
| * <!-- end-user-doc --> |
| * <p> |
| * The following features are implemented: |
| * <ul> |
| * <li>{@link org.eclipse.stem.diseasemodels.forcing.impl.GaussianForcingDiseaseModelImpl#getSigma2 <em>Sigma2</em>}</li> |
| * <li>{@link org.eclipse.stem.diseasemodels.forcing.impl.GaussianForcingDiseaseModelImpl#getModulationPeriod <em>Modulation Period</em>}</li> |
| * <li>{@link org.eclipse.stem.diseasemodels.forcing.impl.GaussianForcingDiseaseModelImpl#getModulationPhaseShift <em>Modulation Phase Shift</em>}</li> |
| * <li>{@link org.eclipse.stem.diseasemodels.forcing.impl.GaussianForcingDiseaseModelImpl#getModulationFloor <em>Modulation Floor</em>}</li> |
| * </ul> |
| * </p> |
| * |
| * @generated |
| */ |
| public class GaussianForcingDiseaseModelImpl extends StochasticSIRDiseaseModelImpl implements GaussianForcingDiseaseModel { |
| /** |
| * The default value of the '{@link #getSigma2() <em>Sigma2</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getSigma2() |
| * @generated |
| * @ordered |
| */ |
| protected static final double SIGMA2_EDEFAULT = 100.0; |
| /** |
| * The cached value of the '{@link #getSigma2() <em>Sigma2</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getSigma2() |
| * @generated |
| * @ordered |
| */ |
| protected double sigma2 = SIGMA2_EDEFAULT; |
| /** |
| * The default value of the '{@link #getModulationPeriod() <em>Modulation Period</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationPeriod() |
| * @generated |
| * @ordered |
| */ |
| protected static final double MODULATION_PERIOD_EDEFAULT = 365.25; |
| /** |
| * The cached value of the '{@link #getModulationPeriod() <em>Modulation Period</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationPeriod() |
| * @generated |
| * @ordered |
| */ |
| protected double modulationPeriod = MODULATION_PERIOD_EDEFAULT; |
| /** |
| * The default value of the '{@link #getModulationPhaseShift() <em>Modulation Phase Shift</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationPhaseShift() |
| * @generated |
| * @ordered |
| */ |
| protected static final double MODULATION_PHASE_SHIFT_EDEFAULT = 0.0; |
| /** |
| * The cached value of the '{@link #getModulationPhaseShift() <em>Modulation Phase Shift</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationPhaseShift() |
| * @generated |
| * @ordered |
| */ |
| protected double modulationPhaseShift = MODULATION_PHASE_SHIFT_EDEFAULT; |
| /** |
| * The default value of the '{@link #getModulationFloor() <em>Modulation Floor</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationFloor() |
| * @generated |
| * @ordered |
| */ |
| protected static final double MODULATION_FLOOR_EDEFAULT = 0.6; |
| /** |
| * The cached value of the '{@link #getModulationFloor() <em>Modulation Floor</em>}' attribute. |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @see #getModulationFloor() |
| * @generated |
| * @ordered |
| */ |
| protected double modulationFloor = MODULATION_FLOOR_EDEFAULT; |
| |
| private Calendar calendar = Calendar.getInstance(); |
| private Calendar calendar2 = Calendar.getInstance(); |
| private static final double MILLIS_PER_DAY = 1000.0*60.0*60.0*24.0; |
| |
| private static final int WINDOWSIZE = 30; |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated NOT |
| */ |
| public GaussianForcingDiseaseModelImpl() { |
| super(); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| protected EClass eStaticClass() { |
| return ForcingPackage.Literals.GAUSSIAN_FORCING_DISEASE_MODEL; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getSigma2() { |
| return sigma2; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setSigma2(double newSigma2) { |
| double oldSigma2 = sigma2; |
| sigma2 = newSigma2; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__SIGMA2, oldSigma2, sigma2)); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getModulationPeriod() { |
| return modulationPeriod; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setModulationPeriod(double newModulationPeriod) { |
| double oldModulationPeriod = modulationPeriod; |
| modulationPeriod = newModulationPeriod; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PERIOD, oldModulationPeriod, modulationPeriod)); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getModulationPhaseShift() { |
| return modulationPhaseShift; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setModulationPhaseShift(double newModulationPhaseShift) { |
| double oldModulationPhaseShift = modulationPhaseShift; |
| modulationPhaseShift = newModulationPhaseShift; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PHASE_SHIFT, oldModulationPhaseShift, modulationPhaseShift)); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public double getModulationFloor() { |
| return modulationFloor; |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| public void setModulationFloor(double newModulationFloor) { |
| double oldModulationFloor = modulationFloor; |
| modulationFloor = newModulationFloor; |
| if (eNotificationRequired()) |
| eNotify(new ENotificationImpl(this, Notification.SET, ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_FLOOR, oldModulationFloor, modulationFloor)); |
| } |
| |
| /** |
| * To create a user defined (experimental) disease mode, the user may wish |
| * to extend the existing models but simply override two methods. |
| * computeDiseaseDeltas() sets the work flow for the model. The following example |
| * is derived from the standard code for all build in models but in the line below |
| * where we define final "double transmisionRate" we show and example modification |
| * where a periodic forcing factor { sin(freq*t) } is added to the code as an example |
| * modification. |
| * |
| * @param time |
| * current time |
| * @param currentState |
| * the current state of the population |
| * @param diseaseLabel |
| * the disease label for which the state transitions are being |
| * computed. |
| * @param timeDelta |
| * the time period (milliseconds) over which the population |
| * members transition to new states |
| * @return a disease state label value that contains the number of |
| * delta changes in each disease state |
| * |
| * Users can modify the method below to create their own model. |
| * |
| * @see org.eclipse.stem.diseasemodels.standard.impl.SIImpl#computeDiseaseDeltas(StandardDiseaseModelLabelValue, |
| * StandardDiseaseModelLabel, long) |
| */ |
| private ArrayList<Long> writtedTimes = new ArrayList<Long>(); |
| private FileWriter fw; |
| private Long firstTime = new Long(Long.MAX_VALUE); |
| @Override |
| public StandardDiseaseModelLabelValue computeDiseaseDeltas( |
| final STEMTime time, |
| final StandardDiseaseModelLabelValue currentState, |
| final StandardDiseaseModelLabel diseaseLabel, final long timeDelta, DiseaseModelLabelValue returnValue) { |
| final SIRLabelValue currentSIR = (SIRLabelValue) currentState; |
| |
| long currentMillis = time.getTime().getTime(); |
| double modulationPeriod = getModulationPeriod(); |
| double phase = getModulationPhaseShift(); |
| double sigma2 = getSigma2(); |
| |
| synchronized(this) { |
| if(firstTime.longValue() == Long.MAX_VALUE) |
| firstTime = Long.valueOf(time.getTime().getTime()); |
| } |
| |
| // Get the day from time and adjust for the phase |
| double day=0; |
| // Shared calendar object not thread safe |
| synchronized(calendar) { |
| calendar.setTimeInMillis(currentMillis+(long)(phase*MILLIS_PER_DAY)); |
| calendar2.setTimeInMillis(firstTime.longValue()); |
| day = (calendar.getTimeInMillis() - calendar2.getTimeInMillis())/MILLIS_PER_DAY; |
| |
| } |
| double intDay = day; |
| int year = (int)(day / modulationPeriod); |
| int nextDayYear = (int)((day+1) / modulationPeriod); |
| |
| double modulatedTransmissionRate = 0.0; |
| |
| // Smoothing |
| if(day % modulationPeriod < WINDOWSIZE || day % modulationPeriod > modulationPeriod-WINDOWSIZE) { |
| double pos; |
| double fday = ((day % modulationPeriod)-modulationPeriod/2.0)/modulationPeriod; |
| double f1 = (getAdjustedTransmissionRate(timeDelta)) * (getModulationFloor() + (1-getModulationFloor())*Math.exp(-(Math.pow(fday,2))/(2*sigma2))); |
| double f2 = (getAdjustedTransmissionRate(timeDelta)) * (getModulationFloor() + (1-getModulationFloor())*Math.exp(-(Math.pow(fday-1,2))/(2*sigma2))); |
| double f3 = (getAdjustedTransmissionRate(timeDelta)) * (getModulationFloor() + (1-getModulationFloor())*Math.exp(-(Math.pow(fday+1,2))/(2*sigma2))); |
| |
| |
| if((day % modulationPeriod) > modulationPeriod - WINDOWSIZE) |
| pos = WINDOWSIZE - (modulationPeriod - (day % modulationPeriod)); |
| else pos = WINDOWSIZE+(day % modulationPeriod); |
| |
| pos = Math.round(pos); |
| double smooth = 0; |
| if(pos == WINDOWSIZE && nextDayYear == year) |
| smooth = 0.5*f1 + 0.5*f3; |
| else if(pos == WINDOWSIZE && nextDayYear == year+1) |
| smooth = 0.5*f1 +0.5*f2; |
| else if (pos < WINDOWSIZE) |
| smooth = ((2*WINDOWSIZE-pos)/(2*WINDOWSIZE))*f1 + (pos/(2*WINDOWSIZE))*f2; |
| else |
| smooth = (pos/(2*WINDOWSIZE))*f1 + ((2*WINDOWSIZE-pos)/(2*WINDOWSIZE))*f3; |
| |
| modulatedTransmissionRate = smooth; |
| } else { |
| double fday = ((day % modulationPeriod)-modulationPeriod/2.0)/modulationPeriod; |
| modulatedTransmissionRate = (getAdjustedTransmissionRate(timeDelta)) * (getModulationFloor() + (1-getModulationFloor())*Math.exp(-(Math.pow(fday,2))/(2*sigma2))); |
| } |
| |
| |
| /* if(diseaseLabel.getNode().getURI().toString().hashCode() % 2 == 0 |
| && ((int)intDay) == 3285) |
| transmissionRate *=10; |
| else |
| if(diseaseLabel.getNode().getURI().toString().hashCode() % 2 == 1 |
| && ((int)intDay) == 3285) |
| transmissionRate *=0.1; |
| */ |
| |
| |
| if(!this.isFrequencyDependent()) modulatedTransmissionRate *= getTransmissionRateScaleFactor(diseaseLabel); |
| |
| /* synchronized(writtedTimes) { |
| if(!writtedTimes.contains(time.getTime().getTime())) { |
| try { |
| if(fw == null) fw = new FileWriter("beta.csv"); |
| fw.write((int)intDay+","+modulatedTransmissionRate+"\n"); |
| fw.flush(); |
| } catch(Exception e) { |
| e.printStackTrace(); |
| } |
| writtedTimes.add(time.getTime().getTime()); |
| } |
| } |
| */ |
| // The effective Infectious population is a dimensionles number normalize by total |
| // population used in teh computation of bets*S*i where i = Ieffective/Pop. |
| // This includes a correction to the current |
| // infectious population (Ieffective) based on the conserved exchange of people (circulation) |
| // between regions. Note that this is no the "arrivals" and "departures" which are |
| // a different process. |
| final double effectiveInfectious = getNormalizedEffectiveInfectious(diseaseLabel.getNode(), diseaseLabel, currentSIR.getI()); |
| |
| /* |
| * Compute state transitions |
| * |
| * Regarding computing the number of transitions from Susceptible to Exposed: |
| * In a linear model the "effective" number of infectious people is just |
| * the number of infectious people In a nonlinear model we have a |
| * nonLinearity exponent that is > 1 this models the effect of immune |
| * system saturation when Susceptible people are exposed to large |
| * numbers of infectious people. then the "effective" number of |
| * infectious people is I^nonLinearity exponent to allow for either |
| * linear or nonlinear models we always calculate I^nonLinearity |
| * exponent and allow nonLinearity exponent >= 1.0 |
| */ |
| double numberOfInfectedToRecovered = getAdjustedRecoveryRate(timeDelta) |
| * currentSIR.getI(); |
| double numberOfRecoveredToSusceptible = getAdjustedImmunityLossRate(timeDelta) |
| * currentSIR.getR(); |
| // Need to checked what do do here. If non linear coefficient is not 1 and |
| // the effective infectious is negative (which is possible), what do do? |
| // Let's fall back on the linear method for now. |
| double numberOfSusceptibleToInfected = 0.0; |
| if(getNonLinearityCoefficient() != 1.0 && effectiveInfectious >= 0.0) |
| numberOfSusceptibleToInfected = modulatedTransmissionRate |
| * currentSIR.getS()* Math.pow(effectiveInfectious, getNonLinearityCoefficient()); |
| else |
| numberOfSusceptibleToInfected = modulatedTransmissionRate |
| * currentSIR.getS()* effectiveInfectious; |
| |
| |
| // Determine delta S |
| final double deltaS = numberOfRecoveredToSusceptible - numberOfSusceptibleToInfected; |
| // Determine delta I |
| final double deltaI = numberOfSusceptibleToInfected- numberOfInfectedToRecovered; |
| // Determine delta R |
| final double deltaR = numberOfInfectedToRecovered - numberOfRecoveredToSusceptible; |
| |
| SIRLabelValueImpl ret = (SIRLabelValueImpl)returnValue; |
| ret.setS(deltaS); |
| ret.setI(deltaI); |
| ret.setIncidence(numberOfSusceptibleToInfected); |
| ret.setR(deltaR); |
| ret.setDiseaseDeaths(0); |
| return ret; |
| } // computeTransitions |
| |
| /** |
| * ModelSpecificAdjustments for a Stochastic model adds noise to or adjusts |
| * the disease state transition values by multiplying |
| * the additions by a random variable r ~ (1+/-x) with x small. |
| * The requirements that no more individuals can be moved from a state than are |
| * already in that state is still enforced. |
| * |
| */ |
| @Override |
| public void doModelSpecificAdjustments( |
| final StandardDiseaseModelLabelValue state) { |
| final SILabelValue currentSI = (SILabelValue) state; |
| double oldI = currentSI.getI(); |
| double Inoisy = currentSI.getI()* computeNoise(); |
| double change = oldI-Inoisy; |
| currentSI.setI(Inoisy); |
| double newS = currentSI.getS() + change; |
| if(newS < 0.0) { |
| // Need to rescale |
| double scale = (currentSI.getS() + newS) / currentSI.getS(); |
| currentSI.setI(Inoisy*scale); |
| } else currentSI.setS(newS); |
| return; |
| } // doModelSpecificAdjustments |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public Object eGet(int featureID, boolean resolve, boolean coreType) { |
| switch (featureID) { |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__SIGMA2: |
| return getSigma2(); |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PERIOD: |
| return getModulationPeriod(); |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PHASE_SHIFT: |
| return getModulationPhaseShift(); |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_FLOOR: |
| return getModulationFloor(); |
| } |
| return super.eGet(featureID, resolve, coreType); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public void eSet(int featureID, Object newValue) { |
| switch (featureID) { |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__SIGMA2: |
| setSigma2((Double)newValue); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PERIOD: |
| setModulationPeriod((Double)newValue); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PHASE_SHIFT: |
| setModulationPhaseShift((Double)newValue); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_FLOOR: |
| setModulationFloor((Double)newValue); |
| return; |
| } |
| super.eSet(featureID, newValue); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public void eUnset(int featureID) { |
| switch (featureID) { |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__SIGMA2: |
| setSigma2(SIGMA2_EDEFAULT); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PERIOD: |
| setModulationPeriod(MODULATION_PERIOD_EDEFAULT); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PHASE_SHIFT: |
| setModulationPhaseShift(MODULATION_PHASE_SHIFT_EDEFAULT); |
| return; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_FLOOR: |
| setModulationFloor(MODULATION_FLOOR_EDEFAULT); |
| return; |
| } |
| super.eUnset(featureID); |
| } |
| |
| /** |
| * <!-- begin-user-doc --> |
| * <!-- end-user-doc --> |
| * @generated |
| */ |
| @Override |
| public boolean eIsSet(int featureID) { |
| switch (featureID) { |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__SIGMA2: |
| return sigma2 != SIGMA2_EDEFAULT; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PERIOD: |
| return modulationPeriod != MODULATION_PERIOD_EDEFAULT; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_PHASE_SHIFT: |
| return modulationPhaseShift != MODULATION_PHASE_SHIFT_EDEFAULT; |
| case ForcingPackage.GAUSSIAN_FORCING_DISEASE_MODEL__MODULATION_FLOOR: |
| return modulationFloor != MODULATION_FLOOR_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(" (sigma2: "); //$NON-NLS-1$ |
| result.append(sigma2); |
| result.append(", modulationPeriod: "); //$NON-NLS-1$ |
| result.append(modulationPeriod); |
| result.append(", modulationPhaseShift: "); //$NON-NLS-1$ |
| result.append(modulationPhaseShift); |
| result.append(", modulationFloor: "); //$NON-NLS-1$ |
| result.append(modulationFloor); |
| result.append(')'); |
| return result.toString(); |
| } |
| |
| } //GaussianForcingDiseaseModelImpl |