/*
 * InfluenzaDgl.java created on 12.12.2005
 * de.explosys.influenza.model
 *
 *  Version 1.2 (18.05.2006)
 *
 * Copyright 2006 by Markus Schwehm, University of Tübingen
 */
package de.explosys.influenza.model;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;

/**
 * A deterministic model for pandemic influenza. 
 *
 * @author schwehm
 */
public class InfluenzaDgl implements DifferentialEquation {

    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;
    
    double[][] m;
    double bb;
    
    /**
     * Number of individuals in each age group (per 100,000 total). 
     */
    double[] individuals =  { 16792, 61731, 21477 };
    
    /**
     * 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.5;

    /**
     * Fraction of infectiosity during first half of infectious period.
     */
    double halfInfectiosity = 0.9;
    
    /**
     * Number of prodromal stages
     */
    int prodromalStages = 2;
    
    /**
     * 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 = 500.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 = 500.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;

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

    /**
     * Fraction of contacts between children that take place at school.
     */
    double schoolContactsFraction = 0.5;

    /**
     * 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 = 0.0;
    
    /**
     * 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 ----------
    
    /**
     * Total number of individuals in the population.
     */
    double total;
    
    /**
     * 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];
    
    /**
     * Contagiousness distribution over the infectious stages.
     */
    double[] contagiousness = new double[ISTAGEGROUPS];
    
    /**
     * 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];
    
    /**
     * Initial number of susceptible individuals (by age and risk group). 
     */
    double susceptibles[][] = new double[AGEGROUPS][RISKGROUPS];

    /**
     * The eigenvecor for the largest eigenvalue (used for initialisation of 
     * the population).
     */
    double[] eigenvector; 
    
    /**
     * Initial fraction of newly infected individuals.
     */
    double initialExposed; 

    /**
     * Force of infection.
     */
    double[] lambda = new double[AGEGROUPS];
    
    /**
     * Differential equations 
     */
    public InfluenzaDgl() {
        initialize();
    }

    /**
     * Calculation of parameters.
     */
    protected void initialize() {        
        
        // Distribution of individuals (by age).
        total = 0.0;
        for (int age = 0; age < AGEGROUPS; age++) {
            total += individuals[age]; 
        } 
        
        if (total > 0) {
            initialExposed = 1.0 / total;  
            for (int age = 0; age < AGEGROUPS; age++) {
                ageDistribution[age] = individuals[age] / total;
            }   
        } else {
            initialExposed = 0.0;
            for (int age = 0; age < AGEGROUPS; age++) {
                ageDistribution[age] = 0.0;
            }               
        }
        
        // 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; // potential division by zero!
        
        // 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; // potential division by zero!
        
        // 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 the fraction of hospitalisations prevented by treatment.
        hospPrevTreat = treatHospPrev;
        
        // Initialize what fraction of extremely severe cases dies 
        // (by age and medication).
        deadFraction = new double[AGEGROUPS];
        for (int age = 0; age < AGEGROUPS; age++) {
            deadFraction[age] = 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 contagiousness distribution.
        new Root() {
            double eval(double mu) {
                return computeF50(mu);
            }
        }.findRoot(0, 1);
        
//        System.out.print("\nHalfInfectiosity=" + halfInfectiosity);
//        for (int i = 0; i < contagiousness.length; i++) {
//            System.out.print(String.format("%8.5f ", contagiousness[i]));
//        }        

        // Initialize the rate alpha at which severe and extremely severe cases seek medical help.
        alpha = 1.0 / consultationDelay; // potential division by zero!
                
        // 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]; // potential division by zero!
                }
            }
            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++) {
            final double dead = deadFraction[age];
            final double gam = gamma(age, MED_NO, ITYPE_H);
            tau[age] = 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 = ITYPE_A; type <= ITYPE_M; type++) {
                    beta[ageSus][ageInf][MED_NO][type] = infFact[type]  * b * ageMatrix[ageSus][ageInf];
                }
            }
        }
        
        for (int type = ITYPE_V; type <= ITYPE_H; type++) {
            // severely ill children do not go to school
            beta[AGE_CHILD][AGE_CHILD][MED_NO][type] = infFact[type]  * b * ageMatrix[AGE_CHILD][AGE_CHILD] * (1.0 - schoolContactsFraction);
            
            // severely ill children redistribute contacts to adult and elderly
            double factor = ageMatrix[AGE_CHILD][AGE_CHILD] * ageDistribution(AGE_CHILD)
                                 / (ageMatrix[AGE_ADULT][AGE_CHILD] * ageDistribution(AGE_ADULT) // potential division by zero!
                                  +ageMatrix[AGE_ELDERLY][AGE_CHILD] * ageDistribution(AGE_ELDERLY));
            
            beta[AGE_ADULT][AGE_CHILD][MED_NO][type] = infFact[type] * b * ageMatrix[AGE_ADULT][AGE_CHILD] * (1.0 + factor * schoolContactsFraction);
            
            beta[AGE_ELDERLY][AGE_CHILD][MED_NO][type] = infFact[type]  * b * ageMatrix[AGE_ELDERLY][AGE_CHILD] * (1.0 + factor * schoolContactsFraction);
            
            // all other contacts remain untouched
            for (int ageSus = AGE_CHILD; ageSus < AGEGROUPS; ageSus++) {
                for (int ageInf = AGE_ADULT; ageInf < AGEGROUPS; ageInf++) {
                    beta[ageSus][ageInf][MED_NO][type] = infFact[type]  * b * ageMatrix[ageSus][ageInf];
                }
            }       
        }
        
        for (int ageSus = 0; ageSus < AGEGROUPS; ageSus++) {
            for (int ageInf = 0; ageInf < AGEGROUPS; ageInf++) {
                for (int type = ITYPE_A; type < ITYPEGROUPS; type++) {
                    beta[ageSus][ageInf][MED_YES][type] = beta[ageSus][ageInf][MED_NO][type] * (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[ageSus][RISK_LOW]+susceptibles[ageSus][RISK_HIGH]) 
                                                                       * eBeta[ageSus][ageInf] * expDur * prodromalStages / ESTAGEGROUPS ;
                for (int risk = 0; risk < RISKGROUPS; risk++) {
                    for (int type = 0; type < ITYPEGROUPS; type++) {
                        double temp = susceptibles(ageSus, 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]);
                        nextGenerationMatrix[ageSus][ageInf] += temp;
                    }
                }
            }
        }
        EigenvalueDecomposition eigen = new Matrix(nextGenerationMatrix).eig();
        
        double[] realEigenvalues = eigen.getRealEigenvalues();
        double[] imagEigenvalues = eigen.getImagEigenvalues();
        double maxEigenvalue = Double.NEGATIVE_INFINITY;
        int maxI = -1; 
        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) {
                maxI = i;
                maxEigenvalue = realEigenvalues[i];
            }
        }    
        eigenvector = new double[AGEGROUPS];
        double sum = 0.0;
        for (int age = 0; age < AGEGROUPS; age++) {
            eigenvector[age] = eigen.getV().get(age, maxI);
            sum+=eigenvector[age];
        }
        for (int age = 0; age < AGEGROUPS; age++) {
            eigenvector[age] /= sum; // potential division by zero!
        }
        m = nextGenerationMatrix;
        bb = b;
        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;
    }

    double computeF50(double mu) {
        double product = 1.0;
        double sum = ISTAGEGROUPS;
        if (mu < 1.0) sum =  (1.0 - Math.pow(mu, ISTAGEGROUPS)) / (1.0 - mu); 
        for (int i = 0; i < ISTAGEGROUPS; i++) {
            contagiousness[i] = ISTAGEGROUPS  * product / sum; // potential division by zero!
            product *= mu;
        }
        double allStagesInfectiosity = 0.0;
        for (int i = 0; i < ISTAGEGROUPS; i++) {
            allStagesInfectiosity += contagiousness[i];
        }
        int halfInfectiousStages = ISTAGEGROUPS / 2;
        double halfStagesInfectiosity = 0.0;
        for (int i = 0; i < halfInfectiousStages; i++) {
            halfStagesInfectiosity += contagiousness[i];
        }
        if (halfInfectiousStages+halfInfectiousStages < ISTAGEGROUPS) {
            halfStagesInfectiosity += 0.5 * contagiousness[halfInfectiousStages];
        }

//        double mu1 = 1.0 - mu;
//        int halfInfectiousStages = ISTAGEGROUPS / 2;
//        double halfStagesInfectiosity = 0.0;
//        double term = 1.0;
//        for (int i = 1; i <= halfInfectiousStages; i++) {
//            halfStagesInfectiosity += term;
//            term *= mu1;
//        }
//        double allStagesInfectiosity = halfStagesInfectiosity;
//        for (int i = halfInfectiousStages+1; i <= ISTAGEGROUPS; i++) {
//            allStagesInfectiosity += term;
//            term *= mu1;
//        }
        return halfInfectiosity - (halfStagesInfectiosity / allStagesInfectiosity);
    }
    
     /**
     * @return initial values
     */
    public double[] getInitialY() {
        int  dimension = Ssize + Esize + Asize + Msize + Vsize + Wsize + Xsize + Hsize + Rsize + Dsize + Isize + 5;

        double[] inVector = new double[dimension];
        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);
            }
        }
        for (int age = 0; age < AGEGROUPS; age++) {
            inVector[S(age, RISK_LOW)] -= (1.0 - ageRisk[age]) * eigenvector[age] / total; // potential division by zero!
            inVector[E(age, RISK_LOW, ESTAGE_1)] = (1.0 - ageRisk[age]) *  eigenvector[age] / total; // potential division by zero!
            inVector[S(age, RISK_HIGH)] -= ageRisk[age] * eigenvector[age] / total; // potential division by zero!
            inVector[E(age, RISK_HIGH, ESTAGE_1)] = ageRisk[age] * eigenvector[age] / total; // potential division by zero!
        } 
        
        return inVector;
    }

    public double getInitialX() {
        return 0.0;
    }

    public void eval(double time, double[] y, double[] out) {
        
        // 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 cancellingFactor = 1.0;
                double redistributionFactor1 = 0.0;
                double redistributionFactor2 = 0.0;
                if ((ageSus == AGE_CHILD) && (ageInf==AGE_CHILD)) {
                     if ((time>=schoolClosingBegin) && (time<schoolClosingEnd)) {
                         cancellingFactor = (1.0 - schoolContactsFraction);
                     }
                } else if ((ageSus != AGE_CHILD) && (ageInf!=AGE_CHILD)) {
                    if ((time>=eventCancelingBegin) && (time<eventCancelingEnd)) {
                        cancellingFactor = (1.0 - eventCancelingFraction);
                    }
                } else {
                    if ((time>=schoolClosingBegin) && (time<schoolClosingEnd)) {
                        double factor = ageMatrix(AGE_CHILD, AGE_CHILD)*ageDistribution(AGE_CHILD)
                                             / (ageMatrix(AGE_ADULT, AGE_CHILD)*ageDistribution(AGE_ADULT) // potential division by zero!
                                             + ageMatrix(AGE_ELDERLY, AGE_CHILD)*ageDistribution(AGE_ELDERLY));      
                        
                        redistributionFactor1 = schoolContactsFraction * reDistribute * factor;   
                        if(ageSus==AGE_CHILD) redistributionFactor2 = redistributionFactor1;
                    }
                }
                double generalFactor = 1.0;
                if ((time >= generalReductionBegin) && (time <  generalReductionEnd)) generalFactor *= (1.0 - generalReductionFraction);
                double temp = 0.0;
                for (int eStage = ESTAGEGROUPS - prodromalStages; eStage < ESTAGEGROUPS; eStage++) {
                    temp += (y[E(ageInf, RISK_LOW, eStage)] + y[E(ageInf, RISK_HIGH, eStage)])  * eBeta(ageSus, ageInf) * (cancellingFactor + redistributionFactor1);
                }
                for (int iStage = 0; iStage < ISTAGEGROUPS; iStage++) {
                    temp += (y[A(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_A) * (cancellingFactor + redistributionFactor1) 
                             + y[M(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_M) * (cancellingFactor + redistributionFactor1)
                             
                             + y[V(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_V) * (1.0 + redistributionFactor2)
                             + y[X(ageInf, iStage)] * beta(ageSus, ageInf, MED_NO, ITYPE_X) * (1.0 + redistributionFactor2)
                             + y[W(ageInf, iStage, MED_NO)] * beta(ageSus, ageInf, MED_NO, ITYPE_W) * (1.0 + redistributionFactor2)
                             + y[H(ageInf, iStage, MED_NO)] * beta(ageSus, ageInf, MED_NO, ITYPE_H) * (1.0 + redistributionFactor2)
                             + y[W(ageInf, iStage, MED_YES)] * beta(ageSus, ageInf, MED_YES, ITYPE_W) * (1.0 + redistributionFactor2)
                             + y[H(ageInf, iStage, MED_YES)] * beta(ageSus, ageInf, MED_YES, ITYPE_H) * (1.0 + redistributionFactor2)) * contagiousness[iStage];
                }                   
                lambda[ageSus] += temp * generalFactor;
            }
        }
        // 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))) * 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)) * 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)) * 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))  * 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) * 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) * 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) * 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) * 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) * y[X(age, iStage)];
                for (int med = 0; med < MEDGROUPS; med++) {
                    dead += tau(age) * 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) * y[X(AGE_ADULT, iStage)];
            for (int med = 0; med < MEDGROUPS; med++) {
                dead += tau(AGE_ADULT) * 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) {
        return tau[age];
    }
    
    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;
        }
    }
    
    /**
     * @param args
     */
    public static void main(String[] args) {
        long time1 = System.currentTimeMillis();

        InfluenzaDgl influenzaDgl = new InfluenzaDgl();
        double[] startVector =  influenzaDgl.getInitialY();
        double startX = influenzaDgl.getInitialX();
        
        RungeKutta4 rk4 = new RungeKutta4(influenzaDgl);
        
        rk4.setStart(startX ,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[] 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.printf("Day=%3d",t); //$NON-NLS-1$
            System.out.printf(", S=%7.4f", S[t ]*100); //$NON-NLS-1$
            System.out.printf(", E=%6.4f",E[t]*100); //$NON-NLS-1$
            System.out.printf(", A=%6.4f",A[t]*100); //$NON-NLS-1$
            System.out.printf(", M=%6.4f",M[t]*100); //$NON-NLS-1$
            System.out.printf(", V=%6.4f",S0[t]*100); //$NON-NLS-1$
            System.out.printf(", W=%6.4f",SA[t]*100); //$NON-NLS-1$
            System.out.printf(", X=%6.4f",H0[t]*100); //$NON-NLS-1$
            System.out.printf(", H=%6.4f",HA[t]*100); //$NON-NLS-1$
            System.out.printf(", R=%6.4f",R[t]*100); //$NON-NLS-1$
            System.out.printf(", D=%6.4f",D[t]*100); //$NON-NLS-1$
            System.out.printf(", I=%6.4f",I[t]*100); //$NON-NLS-1$
            System.out.printf(", WR=%6.4f",(WR[t]*100)); //$NON-NLS-1$
            System.out.printf(", steps=%3d", steps);  //$NON-NLS-1$
            System.out.printf(", Proof=%6.4f",Proof[t]*100); //$NON-NLS-1$
            System.out.println(", Error=" + (Proof[t]-1.0)); //$NON-NLS-1$
        }
        long time2 = System.currentTimeMillis();
         System.out.println("\nTime = " + (time2 - time1)); //$NON-NLS-1$
    }
}