/*
 * InfluenzaDgl.java created on 12.12.2005
 * de.explosys.sbt.core
 *
 * Copyright 2005 by Markus Schwehm, University of Tübingen
 */
package de.explosys.influenza.model;

import Jama.Matrix;

/**
 * InfluenzaDgl
 *
 * @author schwehm
 */
public class InfluenzaDgl implements DiffEqu {

    static final int AGE_CHILD = 0;
    static final int AGE_ADULT = 1;
    static final int AGE_ELDERLY = 2;
    static final int AGEGROUPS = 3;

    static final int RISK_LOW = 0;
    static final int RISK_HIGH = 1;
    static final int RISKGROUPS = 2;

    static final int MED_NO = 0;
    static final int MED_YES = 1;
    static final int MEDGROUPS = 2;

    static final int ESTAGE_1 = 0;
    static final int ESTAGE_2 = 1;
    static final int ESTAGE_LAST = 6;
    static final int ESTAGEGROUPS = 7;

    static final int ISTAGE_1 = 0;
    static final int ISTAGE_2 = 1;
    static final int ISTAGE_LAST = 18;
    static final int ISTAGEGROUPS = 19;

    static final int RSTAGE_1 = 0;
    static final int RSTAGE_2 = 1;
    static final int RSTAGE_LAST = 8;
    static final int RSTAGEGROUPS = 9;

    static final int ITYPE_A = 0;
    static final int ITYPE_M = 1;
    static final int ITYPE_V = 2;
    static final int ITYPE_W = 3;
    static final int ITYPE_X = 4;
    static final int ITYPE_H = 5;
    static final int ITYPEGROUPS = 6;

    static final int Ssize = AGEGROUPS * RISKGROUPS;
    static final int Esize = AGEGROUPS * RISKGROUPS * ESTAGEGROUPS;
    static final int Asize = AGEGROUPS * ISTAGEGROUPS;
    static final int Msize = AGEGROUPS * ISTAGEGROUPS;
    static final int Vsize = AGEGROUPS * ISTAGEGROUPS;
    static final int Wsize = AGEGROUPS * ISTAGEGROUPS * MEDGROUPS;
    static final int Xsize = AGEGROUPS * ISTAGEGROUPS;
    static final int Hsize = AGEGROUPS * ISTAGEGROUPS * MEDGROUPS;
    static final int Rsize = AGEGROUPS * RSTAGEGROUPS;
    static final int Dsize = AGEGROUPS;
    static final int Isize = AGEGROUPS;

    static final int Soffset = 0;
    static final int Eoffset = Soffset + Ssize;
    static final int Aoffset = Eoffset + Esize;
    static final int Moffset = Aoffset + Asize;
    static final int Voffset = Moffset + Msize;
    static final int Woffset = Voffset + Vsize;
    static final int Xoffset = Woffset + Wsize;
    static final int Hoffset = Xoffset + Xsize;
    static final int Roffset = Hoffset + Hsize;
    static final int Doffset = Roffset + Rsize;
    static final int Ioffset = Doffset + Dsize;
    static final int WorkReduction = Ioffset + Isize;
    static final int Outpatients = WorkReduction + 1;
    static final int Antivirals = Outpatients + 1;
    static final int Hospitalisation = Antivirals + 1;
    static final int WorkLossC = Hospitalisation + 1;
    
    /**
     * Number of individuals in each age group (per 100,000 total). 
     */
    double[] individuals =  { 16900, 60600, 22500 };
    
    /**
     * Fraction of individuals with high risk in each age group.
     */
    double[] ageRisk = { 0.06, 0.14, 0.47 };

    /**
     * Upper right half of the age mixing matrix.
     */
    double[] childrenMix = { 132.01, 27.69, 11.50};
    double[] workingMix = { 60.98, 26.06};
    double[] retiredMix = { 54.22 };

    /**
     * Average duration of the latent period [days].
     */
    double expDur = 1.9;

    /**
     * Average duration of the infectious period for young, adult and elderly moderately sick cases, respectively [days].
     */
    double[] infDurMild = { 7.0, 4.1, 4.1 };
    
    /**
     * Average duration of the infectious period of young, adult and elderly severely sick cases, respectively [days].
     */
    double[] infDurSevere = { 7.0, 7.0, 7.0 };
    
    /**
     * Average duration of the reconvalescence period [days].
     */
    double recDur = 5.0;

    /**
     * Fraction of asymptomatic cases out of all infected individuals.
     */
    double asymptomaticFraction = 0.33;
    
    /**
     * Fraction of severe cases out of all symptomatic cases.
     */
    double severeFraction = 0.5;

    /**
     * Fraction of young, adult and elderly extremely severe cases, respectively, out of all untreated severe cases in the low risk group.
     */
    double[] hospitalLow = {0.00187, 0.02339, 0.03560 };
    
    /**
     * Fraction of young, adult and elderly extremely severe cases, respectively, out of all untreated severe cases in the high risk group.
     */
    double[] hospitalHigh = { 0.01333, 0.02762, 0.07768 };

    /**
     * Fraction of young, adult and elderly untreated extremely severe cases, respectively, who die of influenza.
     */
    double[] deadFractionUntreated =  { 0.05541, 0.16531, 0.39505 };
    
    /**
     * Basic reproduction number.
     */
    double r0 = 2.0;

    /**
     * Contagiousness during prodromal period relative to moderately sick cases.
     */
    double eInfFact = 0.5;
    
    /**
     * Contagiousness of asymptomatic cases relative to moderately sick cases.
     */
    double infFactAsymptomatic = 0.5;
    
    /**
     * Contagiousness of severe and extremely severe cases relative to moderately sick cases.
     */
    double infFactSevere = 1.0;
    
    /**
     * Fraction of the population for which antiviral treatment is available.
     */
    double antiviralRessource = 1.0;

    /**
     * Average time after onset of symptoms when severe and extremely severe cases seek medical help [days].
     */
    double consultationDelay = 1.0;
    
    /**
     * Maximum time after onset of symptoms until which antiviral treatment can be given to cases [days].
     */
    double maxTreatmentDelay = 2;

    /**
     * Fraction of severe and extremely severe cases receiving antiviral treament.
     */
    double severeTreatFract = 0.0;
    
    /**
     * First day of antiviral treatment for severe and extremely severe cases.
     */
    double severeTreatBegin = 0.0;
    
    /**
     * Day at which antiviral treatment of severe cases is terminated.
     */
    double severeTreatEnd = 0.0;

    /**
     * Fraction of extremely severe cases receiving antiviral treament.
     */
    double extremeTreatFract = 0.0;

    /**
     * First day of antiviral treatment of extremely severe cases.
     */
    double extremeTreatBegin = 0.0;
    
    /**
     * Day at which antiviral treatment of extremely severe cases is terminated.
     */
    double extremeTreatEnd = 0.0;

    /**
     * Factor by which the contagiousness is reduced through antiviral treatment.
     */
    double treatEfficacyInfectiousness = 0.8;

    /**
     * Factor by which the infectious duration is reduced through antiviral treatment.
     */
    double infDurReduction = 0.25;
    
    /**
     * Fraction of hospitalisations which can be prevented by antiviral treatment.
     */
    double treatHospPrev = 0.5;

    /**
     * Fraction of deaths which can be prevented by antiviral treatment.
     */
    double deadPrevTreat = 0.8;

    /**
     * Factor by which the contacts of moderately sick cases can be reduced through partial isolation.
     */
    double isolationModerately = 0.0;

    /**
     * Factor by which the contacts of severe cases can be reduced through partial isolation.
     * The same reduction factor is used for extremely severe cases before they have seen a doctor and for extremely severe cases whose hospitalization could be prevented by antiviral treatment.
     */
    double isolationHome = 0.0;

    /**
     * Factor by which the contacts of extremely severe hospitalized cases can be reduced through partial isolation.
     */
    double isolationHospital = 0.0;

    /**
     * Factor by which the contacts in the general population are reduced.
     */
    double generalReductionFraction = 0.0;

    /**
     * Day at which the contactl reduction in the general population starts. 
     */
    double generalReductionBegin = 0.0;

    /**
     * Dayat which the general reduction of contacts is terminated.
     */
    double generalReductionEnd = 0.0;

    /**
     * Factor by which contacts among children are reduced by closing schools and day care centers.
     */
    double schoolClosingFraction = 0.0;

    /**
     * Day at which schools and day care centers are closed.
     */
    double schoolClosingBegin = 0.0;

    /**
     * Day at which schools and day care centers are opened again.
     */
    double schoolClosingEnd = 0.0;

    /**
     * Fraction of those contacts among children that are avoided through school closing which are actually re-distributed to adults.
     */
    double reDistribute;
    
    /**
     * Factor by which contacts among adult and elderly people are reduced by cancelling events of mass gathering.
     */
    double eventCancelingFraction = 0.0;

    /**
     * Day at which events of mass gathering are cancelled.
     */
    double eventCancelingBegin = 0.0;

    /**
     * Day at which events of mass gathering are permitted again.
     */
    double eventCancelingEnd = 0.0;

    // ------------------------------------------------------End of parameter section ----------
    
    /**
     * Distribution of the individuals over the age groups.
     */
    double[] ageDistribution = new double[AGEGROUPS];

    /**
     * Contact matrix depending on the age of the infectious case and the age of the susceptible individual.
     */
    double[][] ageMatrix = new double[AGEGROUPS][AGEGROUPS];
    
    /**
     * Transition rate for the latent period.
     */
    double delta; 

    /**
     * Average duration of the infectious period (by age, medication and course of disease).
     */
    double[][][] infDur = new double[AGEGROUPS][MEDGROUPS][ITYPEGROUPS];
    
    /**
     * Transition rate for reconvalescence.
     */
    double rho;
    
    /**
     * Fraction of hospitalized cases (by age and risk group).
     */
    double[][] hospFract = new double[AGEGROUPS][RISKGROUPS];

    /**
     * Case fatality among extremely severe cases.
     */
    double[][] deadFraction;

    /**
     * Contagiousness (by course of disease).
     */
    double[] infFact = new double[ITYPEGROUPS];

    /**
     * Rate at which severe and extremely severe cases seek medical help.
     */    
    double alpha;
    
    /**
     * Fraction of hospitalisations prevented by treatment.
     */
    double hospPrevTreat;

    /**
     * Course of disease (by age and risk group).
     */
    double[][][] destiny = new double[AGEGROUPS][RISKGROUPS][ITYPEGROUPS];

    /**
     * Transition rate during the contagious period (by age, medication and course of disease).
     */
    static double[][][] gamma = new double[AGEGROUPS][MEDGROUPS][ITYPEGROUPS];

    /**
     * Maximum stage of contagiousness during which antiviral treatment can be given (by age).
     */
    int[] maxTreatmentStage = new int[AGEGROUPS];

    /**
     * Effective contact matrix during the prodromal phase of the infection.
     */
    double[][] eBeta = new double[AGEGROUPS][AGEGROUPS];
    
    /**
     *  Effective contact matrix during the late (contagious) phase of the infection (by medication and course of disease).
     */
    double[][][][] beta = new double[AGEGROUPS][AGEGROUPS][MEDGROUPS][ITYPEGROUPS];

    /**
     * Rate at which extremely severe cases die (by age and medication).
     */
    double[][] tau = new double[AGEGROUPS][MEDGROUPS];
    
    /**
     * Initial number of susceptible individuals (by age and risk group). 
     */
    double susceptibles[][] = new double[AGEGROUPS][RISKGROUPS];

    /**
     * Initial fraction of newly infected individuals.
     */
    double initialExposed = 0.00001; //

    /**
     * Force of infection.
     */
    double[] lambda = new double[AGEGROUPS];
    
    double[] out;


    protected static String parameterName;

    protected static double parameter;

    /**
     * Differential equations 
     */
    public InfluenzaDgl() {
        initialize();
    }

    /**
     * Calculation of parameters.
     */
    protected void initialize() {
        
        // Distribution of individuals (by age).
        double total = 0.0;
        for (int age = 0; age < AGEGROUPS; age++) {
            total += individuals[age]; 
        }        
        for (int age = 0; age < AGEGROUPS; age++) {
            ageDistribution[age] = individuals[age] / total;
        }   
        
        // Distribution of individuals (by age and risk group).
        for (int age = 0; age < AGEGROUPS; age++) {
            susceptibles[age][RISK_LOW] = ageDistribution[age] * (1.0 - ageRisk[age]);
            susceptibles[age][RISK_HIGH] = ageDistribution[age] * ageRisk[age];
        }
        
        // Initialize the contact matrix.
        ageMatrix[AGE_CHILD][AGE_CHILD] = childrenMix[0];
        ageMatrix[AGE_CHILD][AGE_ADULT] = childrenMix[1];
        ageMatrix[AGE_CHILD][AGE_ELDERLY] = childrenMix[2];
        ageMatrix[AGE_ADULT][AGE_CHILD] = childrenMix[1];
        ageMatrix[AGE_ADULT][AGE_ADULT] = workingMix[0];
        ageMatrix[AGE_ADULT][AGE_ELDERLY] = workingMix[1];
        ageMatrix[AGE_ELDERLY][AGE_CHILD] = childrenMix[2];
        ageMatrix[AGE_ELDERLY][AGE_ADULT] = workingMix[1];
        ageMatrix[AGE_ELDERLY][AGE_ELDERLY] = retiredMix[0];
        
        // Initialize transition rate delta for the latent period.
        delta = ESTAGEGROUPS / expDur;
        
        // Initialize the average duration of the infectious periods (by age, medication and course of disease).
        for (int age = 0; age < AGEGROUPS; age++) {
            infDur[age][MED_NO][ITYPE_A] = infDurMild[age];
            infDur[age][MED_NO][ITYPE_M] = infDurMild[age];
            infDur[age][MED_NO][ITYPE_V] = infDurSevere[age];
            infDur[age][MED_NO][ITYPE_W] = infDurSevere[age];
            infDur[age][MED_NO][ITYPE_X] = infDurSevere[age];
            infDur[age][MED_NO][ITYPE_H] = infDurSevere[age];
            for (int type = 0; type < ITYPEGROUPS; type++) {
                infDur[age][MED_YES][type] = (1.0 - infDurReduction) * infDur[age][MED_NO][type];
            }
        }
        
        // Initialize the transition rate rho for the reconvalescent period.
        rho = RSTAGEGROUPS / recDur;
        
        // Initialize the fraction of extremely severe cases among all severe cases (by age and risk group).
        for (int age = 0; age < AGEGROUPS; age++) {
            hospFract[age][RISK_LOW] = hospitalLow[age];
            hospFract[age][RISK_HIGH] = hospitalHigh[age];
        }
        
        // Initialize what fraction of extremely severe cases dies (by age and medication).
        double xxx = (1.0 - deadPrevTreat) / (1.0 - hospPrevTreat);
        deadFraction = new double[AGEGROUPS][MEDGROUPS];
        for (int age = 0; age < AGEGROUPS; age++) {
            deadFraction[age][MED_NO] = deadFractionUntreated[age];
            deadFraction[age][MED_YES] = xxx * deadFractionUntreated[age];
        }

        // Initialize the contagiousness of cases (by course of disease).
        infFact[ITYPE_A] = infFactAsymptomatic;
        infFact[ITYPE_M] = 1.0;
        infFact[ITYPE_V] = infFactSevere;
        infFact[ITYPE_W] = infFactSevere;
        infFact[ITYPE_X] = infFactSevere;
        infFact[ITYPE_H] = infFactSevere;

        // Initialize the rate alpha at which severe and extremely severe cases seek medical help.
        alpha = 1.0 / consultationDelay;
                
        // Initialize the fraction of hospitalisations prevented by treatment.
        hospPrevTreat = treatHospPrev;
        
        // Initialize what courses of disease cases of different age and risk group will take
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int risk = 0; risk < RISKGROUPS; risk++) {
                destiny[age][risk][ITYPE_A] = asymptomaticFraction;
                destiny[age][risk][ITYPE_M] = (1.0 - asymptomaticFraction) * (1.0 - severeFraction);
                destiny[age][risk][ITYPE_W] = 0.0;
                destiny[age][risk][ITYPE_H] = 0.0;
                destiny[age][risk][ITYPE_V] = (1.0 - asymptomaticFraction) * severeFraction * (1.0 - hospFract[age][risk]);
                destiny[age][risk][ITYPE_X] = (1.0 - asymptomaticFraction) * severeFraction * hospFract[age][risk];
            }
        }
        
        // Initialize the transition rate gamma during the contagious period and the highest stage of contagiousness maxTreatmentStage during which antiviral treatment can be given.
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int med = 0; med < MEDGROUPS; med++) {
                for (int type = 0; type < ITYPEGROUPS; type++) {
                    gamma[age][med][type] = ISTAGEGROUPS / infDur[age][med][type];
                }
            }
            maxTreatmentStage[age] = (int) Math.ceil(maxTreatmentDelay * gamma(age, MED_NO, ITYPE_V));
        }

        // Find a factor x such that the resulting basic reproduction number is as desired.
        // This will implicitly initialize the effective contact matrix beta.
        new Root() {
            double eval(double x) {
                return computeR0(x);
            }
        }.findRoot(0, 1);

        // Modify the effective contact matrix beta to consider the effect of partial isolation of cases.
        for (int ageSus = 0; ageSus < AGEGROUPS; ageSus++) {
            for (int ageInf = 0; ageInf < AGEGROUPS; ageInf++) {
                for (int med = 0; med < MEDGROUPS; med++) {
                    beta[ageSus][ageInf][med][ITYPE_M] *= (1.0 - isolationModerately);
                    beta[ageSus][ageInf][med][ITYPE_V] *= (1.0 - isolationHome);
                    beta[ageSus][ageInf][med][ITYPE_W] *= (1.0 - isolationHome);
                    beta[ageSus][ageInf][med][ITYPE_X] *= (1.0 - isolationHome);
                    beta[ageSus][ageInf][med][ITYPE_H] *= (1.0 - isolationHospital);
                }
            }
        }

        // Initialize the rate tau at which extremely severe cases die (by age and medication).
        // tau is chosen such that the chosen fraction dies during the period of contagiousness.
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int med = 0; med < MEDGROUPS; med++) {
                final double gam = gamma(age, med, ITYPE_H);
                final double dead = deadFraction[age][med];
                tau[age][med] = new Root() {
                    double eval(double t) {
                        return computeDeadFraction(gam, t, dead);
                    }
                }.findRoot(0, 1);
            }
        }        
    }

    // Compute the basic reproduction number as the dominant eigenvalue of the next generation matrix.
    // The effective contact matrices eBeta (during the prodromal period) and beta (during the final contagious period) are also calculated in this iterative process.
    double computeR0(double b) {
        for (int ageSus = 0; ageSus < AGEGROUPS; ageSus++) {
            for (int ageInf = 0; ageInf < AGEGROUPS; ageInf++) {
                eBeta[ageSus][ageInf] = eInfFact *  b * ageMatrix[ageSus][ageInf];
                for (int type = 0; type < ITYPEGROUPS; type++) {
                    beta[ageSus][ageInf][MED_NO][type] = infFact[type]  * b * ageMatrix[ageSus][ageInf];
                    beta[ageSus][ageInf][MED_YES][type] = infFact[type] * b * ageMatrix[ageSus][ageInf] * (1.0-treatEfficacyInfectiousness);
                }
            }
        }
        double[][] nextGenerationMatrix = new double[AGEGROUPS][AGEGROUPS];
        for (int ageSus = 0; ageSus < AGEGROUPS; ageSus++) {
            for (int ageInf = 0; ageInf < AGEGROUPS; ageInf++) {
                nextGenerationMatrix[ageSus][ageInf] =  (susceptibles[ageInf][RISK_LOW]+susceptibles[ageInf][RISK_HIGH])  * eBeta[ageSus][ageInf] * expDur * 2.0 / ESTAGEGROUPS ;
                for (int risk = 0; risk < RISKGROUPS; risk++) {
                    for (int type = 0; type < ITYPEGROUPS; type++) {
                        double temp = susceptibles(ageInf, risk) * destiny(ageInf, risk, type) * beta(ageSus, ageInf, MED_NO, type) * infDur[ageInf][MED_NO][type];
                        if ((type == ITYPE_X) || (type == ITYPE_H))
                            temp *= (1.0 - deadFraction[ageInf][MED_NO]);
                        nextGenerationMatrix[ageSus][ageInf] += temp;
                    }
                }
            }
        }
        double[] realEigenvalues = new Matrix(nextGenerationMatrix).eig().getRealEigenvalues();
        double[] imagEigenvalues = new Matrix(nextGenerationMatrix).eig().getImagEigenvalues();
        double maxEigenvalue = Double.NEGATIVE_INFINITY;
        for (int i = 0; i < realEigenvalues.length; i++) {
            if(imagEigenvalues[i]!=0) System.err.println("Imaginary eigenvalue (" + realEigenvalues[i] + ", " + imagEigenvalues[i] + " ) ignored."); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
            else if (realEigenvalues[i] > maxEigenvalue) maxEigenvalue = realEigenvalues[i];
        }    
        return maxEigenvalue-r0;
    }

    // Calculate the difference between the chosen fraction who dies from influenza and the fraction that would die when the current parameters are used.
    // This procedure will be used to iteratively chose parameter values which will eventually lead to the correct case fatality.
    double computeDeadFraction(double gam, double t, double dead) {
        double result = 0.0;
        double product = 1.0;
        double temp = gam / (gam + t);
        for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
            result += product;
            product *= temp;
        }
        result *= t / (t + gam);
        return result - dead;
    }

     /**
     * @return initial values
     */
    public double[] init() {
        int  dimension = Ssize + Esize + Asize + Msize + Vsize + Wsize + Xsize + Hsize + Rsize + Dsize + Isize + 5;
        out = new double[dimension];

        double[] inVector = out.clone();
        for (int i = 0; i < inVector.length; i++) {
            inVector[i] = 0;
        }
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int risk = 0; risk < RISKGROUPS; risk++) {
                inVector[S(age, risk)] = susceptibles(age, risk);
            }
        }
        inVector[S(AGE_ADULT, RISK_LOW)] -= initialExposed;
        inVector[E(AGE_ADULT, RISK_LOW, ESTAGE_1)] = initialExposed;
        return inVector;
    }


    public void eval(double time, double[] y) {
        
        // Calculate the current force of infection.
        for (int ageSus = 0; ageSus < AGEGROUPS; ageSus++) {
            lambda[ageSus] = 0.0;
            for (int ageInf = 0; ageInf < AGEGROUPS; ageInf++) {
                
                double factor = 1.0;
                if ((ageSus == AGE_CHILD) && (ageInf==AGE_CHILD)) {
                     if ((time>=schoolClosingBegin) && (time<schoolClosingEnd)) {
                         factor = (1.0 - schoolClosingFraction);
                     }
                } else if ((ageSus != AGE_CHILD) && (ageInf!=AGE_CHILD)) {
                    if ((time>=eventCancelingBegin) && (time<eventCancelingEnd)) {
                        factor = (1.0 - eventCancelingFraction);
                    }
                } else {
                    if ((time>=schoolClosingBegin) && (time<schoolClosingEnd)) {
                        double temp;
                        if (ageSus < ageInf) {
                            temp = ageMatrix(AGE_CHILD, AGE_ADULT)*ageDistribution(AGE_ADULT)
                                     + ageMatrix(AGE_CHILD, AGE_ELDERLY)*ageDistribution(AGE_ELDERLY);
                        } else {
                            temp = ageMatrix(AGE_ADULT, AGE_CHILD)*ageDistribution(AGE_ADULT) 
                                     + ageMatrix(AGE_ELDERLY, AGE_CHILD)*ageDistribution(AGE_ELDERLY);
                        }                                               
                        factor = (schoolClosingFraction*reDistribute*ageMatrix(AGE_CHILD, AGE_CHILD)*ageDistribution(AGE_CHILD)/temp + 1);         
                    }
                }
                if ((time >= generalReductionBegin) && (time <  generalReductionEnd)) factor *= (1.0 - generalReductionFraction);
                double temp = 0.0;
                for (int eStage = ESTAGE_LAST - 1; eStage < ESTAGEGROUPS; eStage++) {
                    temp += (y[E(ageInf, RISK_LOW, eStage)] + y[E(ageInf, RISK_HIGH, eStage)])  * eBeta(ageSus, ageInf);
                }
                for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
                    temp += y[A(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_A)
                             + y[M(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_M)
                             + y[V(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_V)
                             + y[X(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_X)
                             + y[W(ageInf, iStage, MED_NO)] * beta(ageSus, ageInf, MED_NO, ITYPE_W)
                             + y[H(ageInf, iStage, MED_NO)] * beta(ageSus, ageInf, MED_NO, ITYPE_H)
                             + y[W(ageInf, iStage, MED_YES)] * beta(ageSus, ageInf, MED_YES, ITYPE_W)
                             + y[H(ageInf, iStage, MED_YES)] * beta(ageSus, ageInf, MED_YES, ITYPE_H);
                }                   
                lambda[ageSus] += temp * factor;
            }
        }
        
        // Calculate what fraction of severe cases are treated.
        double todaySevereTreatFract = 0.0;       
        if ((time>=severeTreatBegin) && (time<severeTreatEnd)) {
            todaySevereTreatFract = severeTreatFract;
            if(antiviralRessource<y[Antivirals]) todaySevereTreatFract = 0.0;
        }
        
        // Calculate what fraction of extremely severe cases are treated.
        double todayExtremeTreatFract = 0.0;       
        if ((time>=extremeTreatBegin) && (time<extremeTreatEnd)) {
            todayExtremeTreatFract = extremeTreatFract;
            if(antiviralRessource<y[Antivirals]) todayExtremeTreatFract = 0.0;
        }
        
        // Susceptible individuals are infected.
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int risk= 0; risk < RISKGROUPS; risk++) {
                     out[S(age, risk)] = -lambda[age] * y[S(age, risk)];
            }
        }
        
        for (int age = 0; age < AGEGROUPS; age++) {
            for (int risk= 0; risk < RISKGROUPS; risk++) {                    
                // The newly infected susceptible individuals start incubating the infection.
                out[E(age, risk, ESTAGE_1)] = lambda[age] * y[S(age, risk)] - delta * y[E(age, risk, ESTAGE_1)];
                // Incubating individuals proceed in the incubation period.
                for (int eStage = ESTAGE_2; eStage < ESTAGEGROUPS; eStage++) {
                    out[E(age, risk, eStage)] = delta * (y[E(age, risk, eStage-1)] - y[E(age, risk, eStage)]);
                }
            }
              
            // Some of the individuals who have just passed through the incubation period become asymptomatic cases now.
            out[A(age, ISTAGE_1)] 
                = delta * (destiny(age, RISK_LOW, ITYPE_A) * y[E(age, RISK_LOW, ESTAGE_LAST)] 
                             + destiny(age, RISK_HIGH, ITYPE_A) * y[E(age, RISK_HIGH, ESTAGE_LAST)])
                 - gamma(age, MED_NO, ITYPE_A) * y[A(age, ISTAGE_1)];
            // Some of the individuals who have just passed through the incubation period become moderately sick cases now.
           out[M(age, ISTAGE_1)] 
                = delta * (destiny(age, RISK_LOW, ITYPE_M) * y[E(age, RISK_LOW, ESTAGE_LAST)] 
                                + destiny(age, RISK_HIGH, ITYPE_M) * y[E(age, RISK_HIGH, ESTAGE_LAST)])
                 - gamma(age, MED_NO, ITYPE_M) * y[M(age, ISTAGE_1)];
           // Some of the individuals who have just passed through the incubation period become severe cases now.
           out[V(age, ISTAGE_1)]
                = delta * (destiny(age, RISK_LOW, ITYPE_V) * y[E(age, RISK_LOW, ESTAGE_LAST)] 
                             + destiny(age, RISK_HIGH, ITYPE_V) * y[E(age, RISK_HIGH, ESTAGE_LAST)])
                 - (gamma(age, MED_NO, ITYPE_V) + alpha) * y[V(age, ISTAGE_1)];
            // Some of the individuals who have just passed through the incubation period become extremely severe cases now.
           out[X(age, ISTAGE_1)] 
                = delta * (destiny(age, RISK_LOW, ITYPE_X) * y[E(age, RISK_LOW, ESTAGE_LAST)] 
                             + destiny(age, RISK_HIGH, ITYPE_X) * y[E(age, RISK_HIGH, ESTAGE_LAST)])
                 - (gamma(age, MED_NO, ITYPE_X) + (alpha + tau(age, MED_NO))) * y[X(age, ISTAGE_1)];

           // Some severe cases in contagiousness stage 1 have received treatment.
           // Some extremely severe cases in contagiousness stage 1 whose hospitalization could be avoided through antiviral treatment, join this group.
           out[W(age, ISTAGE_1, MED_YES)] 
               = alpha * (y[V(age, ISTAGE_1)] * todaySevereTreatFract + hospPrevTreat * y[X(age, ISTAGE_1)] * todayExtremeTreatFract) 
                - gamma(age, MED_YES, ITYPE_W) * y[W(age, ISTAGE_1, MED_YES)];
           // Some extremely severe cases in contagiousness stage 1 have received treatment and are hospitalized (hospitalization could not be prevented).
           out[H(age, ISTAGE_1, MED_YES)] 
               = alpha * todayExtremeTreatFract * (1.0 - hospPrevTreat) * y[X(age, ISTAGE_1)] 
                - (gamma(age, MED_YES, ITYPE_H) + tau(age, MED_YES)) * y[H(age, ISTAGE_1, MED_YES)];
           // Some severe cases in contagiousness stage 1 have seen the doctor but have not received treatment.
           out[W(age, ISTAGE_1, MED_NO)] 
               = alpha * (1.0 - todaySevereTreatFract) * y[V(age, ISTAGE_1)] 
                - gamma(age, MED_NO, ITYPE_W) * y[W(age, ISTAGE_1, MED_NO)];
           // Some extremely severe cases in contagiousness stage 1 have seen the doctor and have been hospitalized, but have not received treatment.
           out[H(age, ISTAGE_1, MED_NO)] 
               = alpha * (1.0 - todayExtremeTreatFract) * y[X(age, ISTAGE_1)]
                - (gamma(age, MED_NO, ITYPE_H) + tau(age, MED_NO)) * y[H(age, ISTAGE_1, MED_NO)];
           
            for (int iStage = ISTAGE_2; iStage < ISTAGEGROUPS; iStage++) {
                int previousStage = iStage - 1;
                // Progress through the contagious period for asymptomatic cases.
                out[A(age, iStage)] 
                    = gamma(age, MED_NO, ITYPE_A) * (y[A(age, previousStage)] - y[A(age, iStage)]);
                // Progress through the contagious period for moderately sick cases.
                out[M(age, iStage)] 
                    = gamma(age, MED_NO, ITYPE_M) * (y[M(age, previousStage)] - y[M(age, iStage)]);
                // Progress through the contagious period for severe cases (some seek medical help and possibly obtain treatment).
                out[V(age, iStage)] 
                    = gamma(age, MED_NO, ITYPE_V) * (y[V(age, previousStage)] - y[V(age, iStage)]) 
                    - alpha * y[V(age, iStage)];
                // Progress through the contagious period for extremely severe cases (some seek medical help, may be hospitalized and/or obtain antiviral treatment).
                out[X(age, iStage)] 
                    = gamma(age, MED_NO, ITYPE_X) * (y[X(age, previousStage)] - y[X(age, iStage)]) 
                    - (alpha + tau(age, MED_NO))  * y[X(age, iStage)];
                if (iStage + 1 < maxTreatmentStage(age)) {
                    // Some severe cases have just been treated, others simply progress in their contagious period.
                    // Some newly treated extremely severe cases whose hospitalization could be prevented by treatment, join this group.
                    out[W(age, iStage, MED_YES)] 
                        = gamma(age, MED_YES, ITYPE_W) * (y[W(age, previousStage, MED_YES)] 
                         - y[W(age, iStage, MED_YES)]) + alpha * (y[V(age, iStage)] * todaySevereTreatFract 
                        + hospPrevTreat * y[X(age, iStage)] * todayExtremeTreatFract);
                    // Some new severe cases have just seen the doctor without receiving antiviral treatment, others simply progress in their contagious period.
                    out[W(age, iStage, MED_NO)] 
                        = gamma(age, MED_NO, ITYPE_W) * (y[W(age, previousStage, MED_NO)] 
                         - y[W(age, iStage, MED_NO)]) + alpha * (1.0 - todaySevereTreatFract) * y[V(age, iStage)];
                    // Some new extremely severe cases have just been treated and hospitalized, other hospitalized casess simply progress in their contagious period.
                    out[H(age, iStage, MED_YES)] 
                        = gamma(age, MED_YES, ITYPE_H) * (y[H(age, previousStage, MED_YES)] 
                         - y[H(age, iStage, MED_YES)]) + (1.0 - hospPrevTreat) * todayExtremeTreatFract * alpha * y[X(age, iStage)] 
                         - tau(age, MED_YES) * y[H(age, iStage, MED_YES)];
                    // Some new extremely severe cases have just been hospitalized without receiving antiviral treatment; other hospitalized cases simply progress in their contagious period.
                    out[H(age, iStage, MED_NO)]  
                        = gamma(age, MED_NO, ITYPE_H) * (y[H(age, previousStage, MED_NO)] 
                         - y[H(age, iStage, MED_NO)]) + (1.0 - todayExtremeTreatFract) * alpha * y[X(age, iStage)] 
                         - tau(age, MED_NO) * y[H(age, iStage, MED_NO)];
                } else {
                    // Formerly treated severe cases progress in their contagious period.
                    out[W(age, iStage, MED_YES)] 
                        = gamma(age, MED_YES, ITYPE_W) * (y[W(age, previousStage, MED_YES)] 
                         - y[W(age, iStage, MED_YES)]);
                    // Some new severe cases have just seen the doctor without receiving antiviral treatment, others simply progress in their contagious period.
                    out[W(age, iStage, MED_NO)] 
                        = gamma(age, MED_NO, ITYPE_W) * (y[W(age, previousStage, MED_NO)] 
                         - y[W(age, iStage, MED_NO)]) + alpha * y[V(age, iStage)];
                    // Formerly treated hospitalized cases progress in their contagious period.
                   out[H(age, iStage, MED_YES)] 
                        = gamma(age, MED_YES, ITYPE_H) * (y[H(age, previousStage, MED_YES)]  
                         - y[H(age, iStage, MED_YES)]) - tau(age, MED_YES) * y[H(age, iStage, MED_YES)];
                   // Some new extremely severe cases have just seen the doctor and have been hospitalized without receiving antiviral treatment, other formerly untreated hospitalized cases simply progress in their contagious period.
                    out[H(age, iStage, MED_NO)] 
                        = gamma(age, MED_NO, ITYPE_H) * (y[H(age, previousStage, MED_NO)]  
                         - y[H(age, iStage, MED_NO)]) + alpha * y[X(age, iStage)] - tau(age, MED_NO) * y[H(age, iStage, MED_NO)];
                }
            }
        }
        
        // Initialize counters. 
        out[Outpatients] = 0.0;
        out[Antivirals] = 0.0;
        out[Hospitalisation] = 0.0;
        for (int age = 0; age < AGEGROUPS; age++) {
            out[R(age, RSTAGE_1)] = 0.0;
            for (int rStage = RSTAGE_2; rStage < RSTAGEGROUPS; rStage++) {
                out[R(age, rStage)] = 0.0;
            }
            out[D(age)] = 0.0;
            out[I(age)] = 0.0;
            
           // All severe and extremely severe cases have to go through a period of reconvalcescence after passing through their contagious period.
            out[R(age, RSTAGE_1)] += gamma(age, MED_NO, ITYPE_V)
            * y[V(age, ISTAGE_LAST)] + gamma(age, MED_NO, ITYPE_X)
            * y[X(age, ISTAGE_LAST)] + gamma(age, MED_NO, ITYPE_W)
            * y[W(age, ISTAGE_LAST, MED_NO)] + gamma(age, MED_NO, ITYPE_H)
            * y[H(age, ISTAGE_LAST, MED_NO)] + gamma(age, MED_YES, ITYPE_W)
            * y[W(age, ISTAGE_LAST, MED_YES)] + gamma(age, MED_YES, ITYPE_H)
            * y[H(age, ISTAGE_LAST, MED_YES)] - rho
            * y[R(age, RSTAGE_1)];

            // Progress during the reconvalescent period.
            for (int rStage = RSTAGE_2; rStage < RSTAGEGROUPS; rStage++) {
                out[R(age, rStage)] += rho * (y[R(age, rStage - 1)] - y[R(age, rStage)]);
            }

            // Update the number of cases who died.
            double dead = 0;
            for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
                dead += tau(age, MED_NO) * y[X(age, iStage)];
                for (int med = 0; med < MEDGROUPS; med++) {
                    dead += tau(age, med) * y[H(age, iStage, med)];
                }
            }
            out[D(age)] += dead;

            // Update the number of fully recovered immune people.
            out[I(age)] += gamma(age, MED_NO, ITYPE_A)
                    * y[A(age, ISTAGE_LAST)] + gamma(age, MED_NO, ITYPE_M)
                    * y[M(age, ISTAGE_LAST)] + rho * y[R(age, RSTAGE_LAST)];
            
            // Update the number of outpatient visits.
            for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
                out[Outpatients] += alpha  * (y[V(age, iStage)] + y[X(age, iStage)]);
            }
            
            // Update the number of antiviral doses used.
            for (int iStage = 0; iStage < maxTreatmentStage[age]; iStage++) {
                out[Antivirals] += alpha 
                        * (y[V(age, iStage)]* todaySevereTreatFract + y[X(age, iStage)]* todayExtremeTreatFract);
            }
            
            // Update the number of hospitalizations.
            for (int iStage = 0; iStage < maxTreatmentStage(age); iStage++) {
                out[Hospitalisation] += alpha
                        * (1.0 - todayExtremeTreatFract + todayExtremeTreatFract * (1.0 - hospPrevTreat)) * y[X(age, iStage)];
            }
            for (int iStage = maxTreatmentStage(age); iStage < ISTAGEGROUPS; iStage++) {
                out[Hospitalisation] += alpha * y[X(age, iStage)];
            }
        }   
        
        // Initialize counters.
        out[WorkReduction] = 0.0;
        out[WorkLossC] = 0.0;
        double dead = 0.0;
        for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
            dead += tau(AGE_ADULT, MED_NO) * y[X(AGE_ADULT, iStage)];
            for (int med = 0; med < MEDGROUPS; med++) {
                dead += tau(AGE_ADULT, med) * y[H(AGE_ADULT, iStage, med)];
            }
        }
        
        // Update the number of work days lost due to influenza.
        out[WorkReduction] += delta
        * (destiny(AGE_ADULT, RISK_LOW, ITYPE_V) * y[E(AGE_ADULT, RISK_LOW, ESTAGE_LAST)] 
        + destiny(AGE_ADULT, RISK_HIGH, ITYPE_V) * y[E(AGE_ADULT, RISK_HIGH, ESTAGE_LAST)]
        + destiny(AGE_ADULT, RISK_LOW, ITYPE_X) * y[E(AGE_ADULT, RISK_LOW, ESTAGE_LAST)]
        + destiny(AGE_ADULT, RISK_HIGH, ITYPE_X) * y[E(AGE_ADULT, RISK_HIGH, ESTAGE_LAST)])
        - rho * y[R(AGE_ADULT, RSTAGE_LAST)] - dead;
        
        // Update the cumulative number of work days lost due to influenza.
        for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
            out[WorkLossC] += y[V(AGE_ADULT, iStage)] + y[X(AGE_ADULT, iStage)];
            for (int med = 0; med < MEDGROUPS; med++) {
                out[WorkLossC] += y[W(AGE_ADULT, iStage, med)] + y[H(AGE_ADULT, iStage, med)];
            }
        }
        for (int rStage = 0; rStage < RSTAGEGROUPS; rStage++) {
            out[WorkLossC] += y[R(AGE_ADULT, rStage)];
        }
    }

    double susceptibles(int age, int risk) {
        return susceptibles[age][risk];
    }
    
    double ageDistribution(int age) {
        return ageDistribution[age];
    }
    
    int maxTreatmentStage(int age) {
        return maxTreatmentStage[age];
    }
    
    double ageMatrix(int ageSus, int ageInf) {
        return ageMatrix[ageSus][ageInf];
    }
    
    double destiny(int age, int risk, int type) {
        return destiny[age][risk][type];
    }
    
    double gamma(int age, int med, int type) {
        return gamma[age][med][type];
    }
    
    double eBeta(int ageSus, int ageInf) {
        return eBeta[ageSus][ageInf];
    }
    
    double beta(int ageSus, int ageInf, int med, int type) {
        return beta[ageSus][ageInf][med][type];
    }
    
    double tau(int age, int med) {
        return tau[age][med];
    }
    
    static int S(int age, int risk) {
        return Soffset + age * RISKGROUPS + risk;
    }

    static int E(int age, int risk, int eStages) {
        return Eoffset + (age * RISKGROUPS + risk) * ESTAGEGROUPS + eStages;
    }

    static int A(int age, int iStage) {
        return Aoffset + age * ISTAGEGROUPS + iStage;
    }

    static int M(int age, int iStage) {
        return Moffset + age * ISTAGEGROUPS + iStage;
    }

    static int V(int age, int iStage) {
        return Voffset + age * ISTAGEGROUPS + iStage;
    }

    static int W(int age, int iStage, int med) {
        return Woffset + (age * ISTAGEGROUPS + iStage) * MEDGROUPS + med;
    }

    static int X(int age, int iStage) {
        return Xoffset + age * ISTAGEGROUPS + iStage;
    }

    static int H(int age, int iStage, int med) {
        return Hoffset + (age * ISTAGEGROUPS + iStage) * MEDGROUPS + med;
    }

    static int R(int age, int rStage) {
        return Roffset + age * RSTAGEGROUPS + rStage;
    }

    static int D(int age) {
        return Doffset + age;
    }

    static int I(int age) {
        return Ioffset + age;
    }

    static abstract class Root {
        abstract double eval(double x);

        double findRoot(double a, double b) throws RuntimeException {
            if ((a <= 0) && (b >= 0))
                return find(a, b);
            else if ((b <= 0) && (a >= 0))
                return find(b, a);
            else
                throw new RuntimeException("Range error (" + a + ", " + b + ")");  //$NON-NLS-1$ //$NON-NLS-2$//$NON-NLS-3$
        }

        private double find(double a, double b) {
            double c = 0.5 * (a + b);
            if (b - a <= 1.0e-15)
                return c;
            double cval = eval(c);
            if (cval < 0)
                return find(c, b);
            if (cval > 0)
                return find(a, c);
            return c;
        }
    }
    
    public double[] getOutVector () {
        return out;
     }

    /**
     * @param args
     */
    public static void main(String[] args) {
//        parameterName = "r0";
//        double parameterMin = 1;
//        double parameterMax = 3;
       double parameterN = 0;
//        double parameterDelta = (parameterMax - parameterMin) / parameterN;
//
        for (int i = 0; i <= parameterN; i++) {
//            parameter = parameterMin + i * parameterDelta;

            // System.out.print("\nStarting with parameter " + parameterName + "
            // = " + parameter + "...");
            long time1 = System.currentTimeMillis();

            InfluenzaDgl influenzaDgl = new InfluenzaDgl();
            double[] startVector =  influenzaDgl.init();
            RungeKutta2 rk4 = new RungeKutta2(influenzaDgl);
            rk4.setError(0.1);
            rk4.setStart(0,startVector);

            int days = 400;
            int steps = 0;
            int oldsteps = 0;
            double[] S = new double[days + 1];
            double[] E = new double[days + 1];
            double[] A = new double[days + 1];
            double[] M = new double[days + 1];
            double[] S0 = new double[days + 1];
            double[] SA = new double[days + 1];
            double[] H0 = new double[days + 1];
            double[] HA = new double[days + 1];
            double[] R = new double[days + 1];
            double[] D = new double[days + 1];
            double[] I = new double[days + 1];
            double[] WR = new double[days + 1];
            double[] OPC = new double[days + 1];
            double[] OP = new double[days + 1];
            double[] HB = new double[days + 1];
            double[] Proof = new double[days + 1];

            double[] outVector = rk4.getY();
            for (int t = 0; t <= days; t++) {
                rk4.run(t);

                steps = rk4.getN() - oldsteps;
                oldsteps = rk4.getN();

                outVector = rk4.getY();

                S[t] = 0;
                E[t] = 0;
                A[t] = 0;
                M[t] = 0;
                S0[t] = 0;
                SA[t] = 0;
                H0[t] = 0;
                HA[t] = 0;
                R[t] = 0;
                D[t] = 0;
                I[t] = 0;
                Proof[t] = 0;
                for (int age = 0; age < AGEGROUPS; age++) {
                    for (int risk = 0; risk < RISKGROUPS; risk++) {
                        S[t] += outVector[S(age, risk)];
                        for (int eStage = 0; eStage < ESTAGEGROUPS; eStage++) {
                            E[t] += outVector[E(age, risk, eStage)];
                        }
                    }
                    for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
                        A[t] += outVector[A(age, iStage)];
                        M[t] += outVector[M(age, iStage)];
                        S0[t] += outVector[V(age, iStage)];
                        H0[t] += outVector[X(age, iStage)];
                        for (int med = 0; med < MEDGROUPS; med++) {
                            SA[t] += outVector[W(age, iStage, med)];
                            HA[t] += outVector[H(age, iStage, med)];
                        }
                    }

                    for (int rStage = 0; rStage < RSTAGEGROUPS; rStage++) {
                        R[t] += outVector[R(age, rStage)];
                    }
                    D[t] += outVector[D(age)];
                    I[t] += outVector[I(age)];
                }
                WR[t] = outVector[WorkReduction];
                OPC[t] = outVector[Outpatients];
                if (t>0) OP[t] = OPC[t] - OPC[t-1];
                
                Proof[t] = S[t] + E[t] + A[t] + M[t] + S0[t] + SA[t] + H0[t]
                        + HA[t]+ R[t] + D[t] + I[t];
//                 System.out.print(String.format("\nDay=%3d",t));
//                 System.out.print(String.format(", S=%7.4f", S[t ]*100));
//                 System.out.print(String.format(", E=%6.4f",E[t]*100));
//                 System.out.print(String.format(", A=%6.4f",A[t]*100));
//                 System.out.print(String.format(", M=%6.4f",M[t]*100));
//                 System.out.print(String.format(", V=%6.4f",S0[t]*100));
//                 System.out.print(String.format(", W=%6.4f",SA[t]*100));
//                 System.out.print(String.format(", X=%6.4f",H0[t]*100));
//                 System.out.print(String.format(", H=%6.4f",HA[t]*100));
//                 System.out.print(String.format(", R=%6.4f",R[t]*100));
//                 System.out.print(String.format(", D=%6.4f",D[t]*100));
//                 System.out.print(String.format(", I=%6.4f",I[t]*100));
//                 System.out.print(String.format(", WR=%6.4f",(WR[t]*100)));
//                 System.out.print(String.format(", steps=%3d", steps)); //$NON-NLS-1$
//                 System.out.print(String.format(", Proof=%6.4f",Proof[t]*100));
//                 System.out.print(", Error=" + (Proof[t]-1.0));
                
                
//                 System.out.print(String.format("\n%3d",t));
//                 System.out.print(String.format(", %6.4f",(WR[t]*100)));
//                 System.out.print(String.format(", %6.4f",(OP[t]*100)));
//                 System.out.print(String.format(", %6.4f",(HB[t]*100)));
//                 System.out.print(String.format(", %6.4f",(TD[t]*100)));
            }
//            System.out.print(String.format("\n%s=%7.4f", parameterName,
//                    parameter));
            System.out.print(String.format(", Dead=%7.4f", D[days] * 100));
            System.out.print(String.format(", Hospitalisation=%7.4f", influenzaDgl.out[Hospitalisation] * 100));
            System.out.print(String.format(", Outpatients=%7.4f", influenzaDgl.out[Outpatients] * 100));
            System.out.print(String.format(", Antivirals=%7.4f", influenzaDgl.out[Antivirals] * 100));
            System.out.print(String.format(", WorkLossC=%7.4f", influenzaDgl.out[WorkLossC] * 100));
            long time2 = System.currentTimeMillis();
             System.out.println("\nTime = " + (time2 - time1)); //$NON-NLS-1$
        }
    }

}