SIRModel.bundle.js

(function (f) {
    if (typeof exports === "object" && typeof module !== "undefined") {
        module.exports = f()
    } else if (typeof define === "function" && define.amd) {
        define([], f)
    } else {
        var g;
        if (typeof window !== "undefined") {
            g = window
        } else if (typeof global !== "undefined") {
            g = global
        } else if (typeof self !== "undefined") {
            g = self
        } else {
            g = this
        }
        g.sirModel = f()
    }
})(function () {
    var define, module, exports;
    return (function () {
        function r(e, n, t) {
            function o(i, f) {
                if (!n[i]) {
                    if (!e[i]) {
                        var c = "function" == typeof require && require;
                        if (!f && c) return c(i, !0);
                        if (u) return u(i, !0);
                        var a = new Error("Cannot find module '" + i + "'");
                        throw a.code = "MODULE_NOT_FOUND", a
                    }
                    var p = n[i] = {
                        exports: {}
                    };
                    e[i][0].call(p.exports, function (r) {
                        var n = e[i][1][r];
                        return o(n || r)
                    }, p, p.exports, r, e, n, t)
                }
                return n[i].exports
            }
            for (var u = "function" == typeof require && require, i = 0; i < t.length; i++) o(t[i]);
            return o
        }
        return r
    })()({
        1: [function (require, module, exports) {
            "use strict";
            /**
             * An implementation of ODEX, by E. Hairer and G. Wanner, ported from the Fortran ODEX.F.
             * The original work carries the BSD 2-clause license, and so does this.
             *
             * Copyright (c) 2016 Colin Smith.
             * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following
             * disclaimer.
             * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the
             * following disclaimer in the documentation and/or other materials provided with the distribution.
             *
             * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
             * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
             * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
             * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
             * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
             * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
             * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
             */
            Object.defineProperty(exports, "__esModule", {
                value: true
            });
            var Outcome;
            (function (Outcome) {
                Outcome[Outcome["Converged"] = 0] = "Converged";
                Outcome[Outcome["MaxStepsExceeded"] = 1] = "MaxStepsExceeded";
                Outcome[Outcome["EarlyReturn"] = 2] = "EarlyReturn";
            })(Outcome = exports.Outcome || (exports.Outcome = {}));
            var Solver = (function () {
                function Solver(n) {
                    this.n = n;
                    this.uRound = 2.3e-16;
                    this.maxSteps = 10000;
                    this.initialStepSize = 1e-4;
                    this.maxStepSize = 0;
                    this.maxExtrapolationColumns = 9;
                    this.stepSizeSequence = 0;
                    this.stabilityCheckCount = 1;
                    this.stabilityCheckTableLines = 2;
                    this.denseOutput = false;
                    this.denseOutputErrorEstimator = true;
                    this.denseComponents = undefined;
                    this.interpolationFormulaDegree = 4;
                    this.stepSizeReductionFactor = 0.5;
                    this.stepSizeFac1 = 0.02;
                    this.stepSizeFac2 = 4.0;
                    this.stepSizeFac3 = 0.8;
                    this.stepSizeFac4 = 0.9;
                    this.stepSafetyFactor1 = 0.65;
                    this.stepSafetyFactor2 = 0.94;
                    this.relativeTolerance = 1e-5;
                    this.absoluteTolerance = 1e-5;
                    this.debug = false;
                }
                Solver.prototype.grid = function (dt, out) {
                    if (!this.denseOutput)
                        throw new Error('Must set .denseOutput to true when using grid');
                    var components = this.denseComponents;
                    if (!components) {
                        components = [];
                        for (var i = 0; i < this.n; ++i)
                            components.push(i);
                    }
                    var t;
                    return function (n, xOld, x, y, interpolate) {
                        if (n === 1) {
                            var v = out(x, y);
                            t = x + dt;
                            return v;
                        }
                        while (t <= x) {
                            var yf = [];
                            for (var _i = 0, components_1 = components; _i < components_1.length; _i++) {
                                var i = components_1[_i];
                                yf.push(interpolate(i, t));
                            }
                            var v = out(t, yf);
                            if (v === false)
                                return false;
                            t += dt;
                        }
                    };
                };
                // Make a 1-based 2D array, with r rows and c columns. The initial values are undefined.
                Solver.dim2 = function (r, c) {
                    var a = new Array(r + 1);
                    for (var i = 1; i <= r; ++i)
                        a[i] = Solver.dim(c);
                    return a;
                };
                // Generate step size sequence and return as a 1-based array of length n.
                Solver.stepSizeSequence = function (nSeq, n) {
                    var a = new Array(n + 1);
                    a[0] = 0;
                    switch (nSeq) {
                        case 1:
                            for (var i = 1; i <= n; ++i)
                                a[i] = 2 * i;
                            break;
                        case 2:
                            a[1] = 2;
                            for (var i = 2; i <= n; ++i)
                                a[i] = 4 * i - 4;
                            break;
                        case 3:
                            a[1] = 2;
                            a[2] = 4;
                            a[3] = 6;
                            for (var i = 4; i <= n; ++i)
                                a[i] = 2 * a[i - 2];
                            break;
                        case 4:
                            for (var i = 1; i <= n; ++i)
                                a[i] = 4 * i - 2;
                            break;
                        case 5:
                            for (var i = 1; i <= n; ++i)
                                a[i] = 4 * i;
                            break;
                        default:
                            throw new Error('invalid stepSizeSequence selected');
                    }
                    return a;
                };
                // Integrate the differential system represented by f, from x to xEnd, with initial data y.
                // solOut, if provided, is called at each integration step.
                Solver.prototype.solve = function (f, x, y0, xEnd, solOut) {
                    var _this = this;
                    // Make a copy of y0, 1-based. We leave the user's parameters alone so that they may be reused if desired.
                    var y = [0].concat(y0);
                    var dz = Solver.dim(this.n);
                    var yh1 = Solver.dim(this.n);
                    var yh2 = Solver.dim(this.n);
                    if (this.maxSteps <= 0)
                        throw new Error('maxSteps must be positive');
                    var km = this.maxExtrapolationColumns;
                    if (km <= 2)
                        throw new Error('maxExtrapolationColumns must be > 2');
                    var nSeq = this.stepSizeSequence || (this.denseOutput ? 4 : 1);
                    if (nSeq <= 3 && this.denseOutput)
                        throw new Error('stepSizeSequence incompatible with denseOutput');
                    if (this.denseOutput && !solOut)
                        throw new Error('denseOutput requires a solution observer function');
                    if (this.interpolationFormulaDegree <= 0 || this.interpolationFormulaDegree >= 7)
                        throw new Error('bad interpolationFormulaDegree');
                    var icom = [0]; // icom will be 1-based, so start with a pad entry.
                    var nrdens = 0;
                    if (this.denseOutput) {
                        if (this.denseComponents) {
                            for (var _i = 0, _a = this.denseComponents; _i < _a.length; _i++) {
                                var c = _a[_i];
                                // convert dense components requested into one-based indexing.
                                if (c < 0 || c > this.n)
                                    throw new Error('bad dense component: ' + c);
                                icom.push(c + 1);
                                ++nrdens;
                            }
                        } else {
                            // if user asked for dense output but did not specify any denseComponents,
                            // request all of them.
                            for (var i = 1; i <= this.n; ++i) {
                                icom.push(i);
                            }
                            nrdens = this.n;
                        }
                    }
                    if (this.uRound <= 1e-35 || this.uRound > 1)
                        throw new Error('suspicious value of uRound');
                    var hMax = Math.abs(this.maxStepSize || xEnd - x);
                    var lfSafe = 2 * km * km + km;

                    function expandToArray(x, n) {
                        // If x is an array, return a 1-based copy of it. If x is a number, return a new 1-based array
                        // consisting of n copies of the number.
                        var tolArray = [0];
                        if (Array.isArray(x)) {
                            return tolArray.concat(x);
                        } else {
                            for (var i = 0; i < n; ++i)
                                tolArray.push(x);
                            return tolArray;
                        }
                    }
                    var aTol = expandToArray(this.absoluteTolerance, this.n);
                    var rTol = expandToArray(this.relativeTolerance, this.n);
                    var _b = [0, 0, 0, 0],
                        nEval = _b[0],
                        nStep = _b[1],
                        nAccept = _b[2],
                        nReject = _b[3];
                    // call to core integrator
                    var nrd = Math.max(1, nrdens);
                    var ncom = Math.max(1, (2 * km + 5) * nrdens);
                    var dens = Solver.dim(ncom);
                    var fSafe = Solver.dim2(lfSafe, nrd);
                    // Wrap f in a function F which hides the one-based indexing from the customers.
                    var F = function (x, y, yp) {
                        var ret = f(x, y.slice(1));
                        for (var i = 0; i < ret.length; ++i)
                            yp[i + 1] = ret[i];
                    };
                    var odxcor = function () {
                        // The following three variables are COMMON/CONTEX/
                        var xOldd;
                        var hhh;
                        var kmit;
                        var acceptStep = function (n) {
                            // Returns true if we should continue the integration. The only time false
                            // is returned is when the user's solution observation function has returned false,
                            // indicating that she does not wish to continue the computation.
                            xOld = x;
                            x += h;
                            if (_this.denseOutput) {
                                // kmit = mu of the paper
                                kmit = 2 * kc - _this.interpolationFormulaDegree + 1;
                                for (var i = 1; i <= nrd; ++i)
                                    dens[i] = y[icom[i]];
                                xOldd = xOld;
                                hhh = h; // note: xOldd and hhh are part of /CONODX/
                                for (var i = 1; i <= nrd; ++i)
                                    dens[nrd + i] = h * dz[icom[i]];
                                var kln = 2 * nrd;
                                for (var i = 1; i <= nrd; ++i)
                                    dens[kln + i] = t[1][icom[i]];
                                // compute solution at mid-point
                                for (var j = 2; j <= kc; ++j) {
                                    var dblenj = nj[j];
                                    for (var l = j; l >= 2; --l) {
                                        var factor = Math.pow((dblenj / nj[l - 1]), 2) - 1;
                                        for (var i = 1; i <= nrd; ++i) {
                                            ySafe[l - 1][i] = ySafe[l][i] + (ySafe[l][i] - ySafe[l - 1][i]) / factor;
                                        }
                                    }
                                }
                                var krn = 4 * nrd;
                                for (var i = 1; i <= nrd; ++i)
                                    dens[krn + i] = ySafe[1][i];
                                // compute first derivative at right end
                                for (var i = 1; i <= n; ++i)
                                    yh1[i] = t[1][i];
                                F(x, yh1, yh2);
                                krn = 3 * nrd;
                                for (var i = 1; i <= nrd; ++i)
                                    dens[krn + i] = yh2[icom[i]] * h;
                                // THE LOOP
                                for (var kmi = 1; kmi <= kmit; ++kmi) {
                                    // compute kmi-th derivative at mid-point
                                    var kbeg = (kmi + 1) / 2 | 0;
                                    for (var kk = kbeg; kk <= kc; ++kk) {
                                        var facnj = Math.pow((nj[kk] / 2), (kmi - 1));
                                        iPt = iPoint[kk + 1] - 2 * kk + kmi;
                                        for (var i = 1; i <= nrd; ++i) {
                                            ySafe[kk][i] = fSafe[iPt][i] * facnj;
                                        }
                                    }
                                    for (var j = kbeg + 1; j <= kc; ++j) {
                                        var dblenj = nj[j];
                                        for (var l = j; l >= kbeg + 1; --l) {
                                            var factor = Math.pow((dblenj / nj[l - 1]), 2) - 1;
                                            for (var i = 1; i <= nrd; ++i) {
                                                ySafe[l - 1][i] = ySafe[l][i] + (ySafe[l][i] - ySafe[l - 1][i]) / factor;
                                            }
                                        }
                                    }
                                    krn = (kmi + 4) * nrd;
                                    for (var i = 1; i <= nrd; ++i)
                                        dens[krn + i] = ySafe[kbeg][i] * h;
                                    if (kmi === kmit)
                                        continue;
                                    // compute differences
                                    for (var kk = (kmi + 2) / 2 | 0; kk <= kc; ++kk) {
                                        var lbeg = iPoint[kk + 1];
                                        var lend = iPoint[kk] + kmi + 1;
                                        if (kmi === 1 && nSeq === 4)
                                            lend += 2;
                                        var l = void 0;
                                        for (l = lbeg; l >= lend; l -= 2) {
                                            for (var i = 1; i <= nrd; ++i) {
                                                fSafe[l][i] -= fSafe[l - 2][i];
                                            }
                                        }
                                        if (kmi === 1 && nSeq === 4) {
                                            l = lend - 2;
                                            for (var i = 1; i <= nrd; ++i)
                                                fSafe[l][i] -= dz[icom[i]];
                                        }
                                    }
                                    // compute differences
                                    for (var kk = (kmi + 2) / 2 | 0; kk <= kc; ++kk) {
                                        var lbeg = iPoint[kk + 1] - 1;
                                        var lend = iPoint[kk] + kmi + 2;
                                        for (var l = lbeg; l >= lend; l -= 2) {
                                            for (var i = 1; i <= nrd; ++i) {
                                                fSafe[l][i] -= fSafe[l - 2][i];
                                            }
                                        }
                                    }
                                }
                                interp(nrd, dens, kmit);
                                // estimation of interpolation error
                                if (_this.denseOutputErrorEstimator && kmit >= 1) {
                                    var errint = 0;
                                    for (var i = 1; i <= nrd; ++i)
                                        errint += Math.pow((dens[(kmit + 4) * nrd + i] / scal[icom[i]]), 2);
                                    errint = Math.sqrt(errint / nrd) * errfac[kmit];
                                    hoptde = h / Math.max(Math.pow(errint, (1 / (kmit + 4))), 0.01);
                                    if (errint > 10) {
                                        h = hoptde;
                                        x = xOld;
                                        ++nReject;
                                        reject = true;
                                        return true;
                                    }
                                }
                                for (var i = 1; i <= n; ++i)
                                    dz[i] = yh2[i];
                            }
                            for (var i = 1; i <= n; ++i)
                                y[i] = t[1][i];
                            ++nAccept;
                            if (solOut) {
                                // If denseOutput, we also want to supply the dense closure.
                                if (solOut(nAccept + 1, xOld, x, y.slice(1), _this.denseOutput && contex(xOldd, hhh, kmit, dens, icom)) === false)
                                    return false;
                            }
                            // compute optimal order
                            var kopt;
                            if (kc === 2) {
                                kopt = Math.min(3, km - 1);
                                if (reject)
                                    kopt = 2;
                            } else {
                                if (kc <= k) {
                                    kopt = kc;
                                    if (w[kc - 1] < w[kc] * _this.stepSizeFac3)
                                        kopt = kc - 1;
                                    if (w[kc] < w[kc - 1] * _this.stepSizeFac4)
                                        kopt = Math.min(kc + 1, km - 1);
                                } else {
                                    kopt = kc - 1;
                                    if (kc > 3 && w[kc - 2] < w[kc - 1] * _this.stepSizeFac3)
                                        kopt = kc - 2;
                                    if (w[kc] < w[kopt] * _this.stepSizeFac4)
                                        kopt = Math.min(kc, km - 1);
                                }
                            }
                            // after a rejected step
                            if (reject) {
                                k = Math.min(kopt, kc);
                                h = posneg * Math.min(Math.abs(h), Math.abs(hh[k]));
                                reject = false;
                                return true; // goto 10
                            }
                            if (kopt <= kc) {
                                h = hh[kopt];
                            } else {
                                if (kc < k && w[kc] < w[kc - 1] * _this.stepSizeFac4) {
                                    h = hh[kc] * a[kopt + 1] / a[kc];
                                } else {
                                    h = hh[kc] * a[kopt] / a[kc];
                                }
                            }
                            // compute stepsize for next step
                            k = kopt;
                            h = posneg * Math.abs(h);
                            return true;
                        };
                        var midex = function (j) {
                            var dy = Solver.dim(_this.n);
                            // Computes the jth line of the extrapolation table and
                            // provides an estimation of the optional stepsize
                            var hj = h / nj[j];
                            // Euler starting step
                            for (var i = 1; i <= _this.n; ++i) {
                                yh1[i] = y[i];
                                yh2[i] = y[i] + hj * dz[i];
                            }
                            // Explicit midpoint rule
                            var m = nj[j] - 1;
                            var njMid = (nj[j] / 2) | 0;
                            for (var mm = 1; mm <= m; ++mm) {
                                if (_this.denseOutput && mm === njMid) {
                                    for (var i = 1; i <= nrd; ++i) {
                                        ySafe[j][i] = yh2[icom[i]];
                                    }
                                }
                                F(x + hj * mm, yh2, dy);
                                if (_this.denseOutput && Math.abs(mm - njMid) <= 2 * j - 1) {
                                    ++iPt;
                                    for (var i = 1; i <= nrd; ++i) {
                                        fSafe[iPt][i] = dy[icom[i]];
                                    }
                                }
                                for (var i = 1; i <= _this.n; ++i) {
                                    var ys = yh1[i];
                                    yh1[i] = yh2[i];
                                    yh2[i] = ys + 2 * hj * dy[i];
                                }
                                if (mm <= _this.stabilityCheckCount && j <= _this.stabilityCheckTableLines) {
                                    // stability check
                                    var del1 = 0;
                                    for (var i = 1; i <= _this.n; ++i) {
                                        del1 += Math.pow((dz[i] / scal[i]), 2);
                                    }
                                    var del2 = 0;
                                    for (var i = 1; i <= _this.n; ++i) {
                                        del2 += Math.pow(((dy[i] - dz[i]) / scal[i]), 2);
                                    }
                                    var quot = del2 / Math.max(_this.uRound, del1);
                                    if (quot > 4) {
                                        ++nEval;
                                        atov = true;
                                        h *= _this.stepSizeReductionFactor;
                                        reject = true;
                                        return;
                                    }
                                }
                            }
                            // final smoothing step
                            F(x + h, yh2, dy);
                            if (_this.denseOutput && njMid <= 2 * j - 1) {
                                ++iPt;
                                for (var i = 1; i <= nrd; ++i) {
                                    fSafe[iPt][i] = dy[icom[i]];
                                }
                            }
                            for (var i = 1; i <= _this.n; ++i) {
                                t[j][i] = (yh1[i] + yh2[i] + hj * dy[i]) / 2;
                            }
                            nEval += nj[j];
                            // polynomial extrapolation
                            if (j === 1)
                                return; // was j.eq.1
                            var dblenj = nj[j];
                            var fac;
                            for (var l = j; l > 1; --l) {
                                fac = Math.pow((dblenj / nj[l - 1]), 2) - 1;
                                for (var i = 1; i <= _this.n; ++i) {
                                    t[l - 1][i] = t[l][i] + (t[l][i] - t[l - 1][i]) / fac;
                                }
                            }
                            err = 0;
                            // scaling
                            for (var i = 1; i <= _this.n; ++i) {
                                var t1i = Math.max(Math.abs(y[i]), Math.abs(t[1][i]));
                                scal[i] = aTol[i] + rTol[i] * t1i;
                                err += Math.pow(((t[1][i] - t[2][i]) / scal[i]), 2);
                            }
                            err = Math.sqrt(err / _this.n);
                            if (err * _this.uRound >= 1 || (j > 2 && err >= errOld)) {
                                atov = true;
                                h *= _this.stepSizeReductionFactor;
                                reject = true;
                                return;
                            }
                            errOld = Math.max(4 * err, 1);
                            // compute optimal stepsizes
                            var exp0 = 1 / (2 * j - 1);
                            var facMin = Math.pow(_this.stepSizeFac1, exp0);
                            fac = Math.min(_this.stepSizeFac2 / facMin, Math.max(facMin, Math.pow((err / _this.stepSafetyFactor1), exp0) / _this.stepSafetyFactor2));
                            fac = 1 / fac;
                            hh[j] = Math.min(Math.abs(h) * fac, hMax);
                            w[j] = a[j] / hh[j];
                        };
                        var interp = function (n, y, imit) {
                            // computes the coefficients of the interpolation formula
                            var a = new Array(31); // zero-based: 0:30
                            // begin with Hermite interpolation
                            for (var i = 1; i <= n; ++i) {
                                var y0_1 = y[i];
                                var y1 = y[2 * n + i];
                                var yp0 = y[n + i];
                                var yp1 = y[3 * n + i];
                                var yDiff = y1 - y0_1;
                                var aspl = -yp1 + yDiff;
                                var bspl = yp0 - yDiff;
                                y[n + i] = yDiff;
                                y[2 * n + i] = aspl;
                                y[3 * n + i] = bspl;
                                if (imit < 0)
                                    continue;
                                // compute the derivatives of Hermite at midpoint
                                var ph0 = (y0_1 + y1) * 0.5 + 0.125 * (aspl + bspl);
                                var ph1 = yDiff + (aspl - bspl) * 0.25;
                                var ph2 = -(yp0 - yp1);
                                var ph3 = 6 * (bspl - aspl);
                                // compute the further coefficients
                                if (imit >= 1) {
                                    a[1] = 16 * (y[5 * n + i] - ph1);
                                    if (imit >= 3) {
                                        a[3] = 16 * (y[7 * n + i] - ph3 + 3 * a[1]);
                                        if (imit >= 5) {
                                            for (var im = 5; im <= imit; im += 2) {
                                                var fac1 = im * (im - 1) / 2;
                                                var fac2 = fac1 * (im - 2) * (im - 3) * 2;
                                                a[im] = 16 * (y[(im + 4) * n + i] + fac1 * a[im - 2] - fac2 * a[im - 4]);
                                            }
                                        }
                                    }
                                }
                                a[0] = (y[4 * n + i] - ph0) * 16;
                                if (imit >= 2) {
                                    a[2] = (y[n * 6 + i] - ph2 + a[0]) * 16;
                                    if (imit >= 4) {
                                        for (var im = 4; im <= imit; im += 2) {
                                            var fac1 = im * (im - 1) / 2;
                                            var fac2 = im * (im - 1) * (im - 2) * (im - 3);
                                            a[im] = (y[n * (im + 4) + i] + a[im - 2] * fac1 - a[im - 4] * fac2) * 16;
                                        }
                                    }
                                }
                                for (var im = 0; im <= imit; ++im)
                                    y[n * (im + 4) + i] = a[im];
                            }
                        };
                        var contex = function (xOld, h, imit, y, icom) {
                            return function (c, x) {
                                var i = 0;
                                for (var j = 1; j <= nrd; ++j) {
                                    // careful: customers describe components 0-based. We record indices 1-based.
                                    if (icom[j] === c + 1)
                                        i = j;
                                }
                                if (i === 0)
                                    throw new Error('no dense output available for component ' + c);
                                var theta = (x - xOld) / h;
                                var theta1 = 1 - theta;
                                var phthet = y[i] + theta * (y[nrd + i] + theta1 * (y[2 * nrd + i] * theta + y[3 * nrd + i] * theta1));
                                if (imit < 0)
                                    return phthet;
                                var thetah = theta - 0.5;
                                var ret = y[nrd * (imit + 4) + i];
                                for (var im = imit; im >= 1; --im) {
                                    ret = y[nrd * (im + 3) + i] + ret * thetah / im;
                                }
                                return phthet + Math.pow((theta * theta1), 2) * ret;
                            };
                        };
                        // preparation
                        var ySafe = Solver.dim2(km, nrd);
                        var hh = Solver.dim(km);
                        var t = Solver.dim2(km, _this.n);
                        // Define the step size sequence
                        var nj = Solver.stepSizeSequence(nSeq, km);
                        // Define the a[i] for order selection
                        var a = Solver.dim(km);
                        a[1] = 1 + nj[1];
                        for (var i = 2; i <= km; ++i) {
                            a[i] = a[i - 1] + nj[i];
                        }
                        // Initial Scaling
                        var scal = Solver.dim(_this.n);
                        for (var i = 1; i <= _this.n; ++i) {
                            scal[i] = aTol[i] + rTol[i] + Math.abs(y[i]);
                        }
                        // Initial preparations
                        var posneg = xEnd - x >= 0 ? 1 : -1;
                        var k = Math.max(2, Math.min(km - 1, Math.floor(-Solver.log10(rTol[1] + 1e-40) * 0.6 + 1.5)));
                        var h = Math.max(Math.abs(_this.initialStepSize), 1e-4);
                        h = posneg * Math.min(h, hMax, Math.abs(xEnd - x) / 2);
                        var iPoint = Solver.dim(km + 1);
                        var errfac = Solver.dim(2 * km);
                        var xOld = x;
                        var iPt = 0;
                        if (solOut) {
                            if (_this.denseOutput) {
                                iPoint[1] = 0;
                                for (var i = 1; i <= km; ++i) {
                                    var njAdd = 4 * i - 2;
                                    if (nj[i] > njAdd)
                                        ++njAdd;
                                    iPoint[i + 1] = iPoint[i] + njAdd;
                                }
                                for (var mu = 1; mu <= 2 * km; ++mu) {
                                    var errx = Math.sqrt(mu / (mu + 4)) * 0.5;
                                    var prod = Math.pow((1 / (mu + 4)), 2);
                                    for (var j = 1; j <= mu; ++j)
                                        prod *= errx / j;
                                    errfac[mu] = prod;
                                }
                                iPt = 0;
                            }
                            // check return value and abandon integration if called for
                            if (false === solOut(nAccept + 1, xOld, x, y.slice(1))) {
                                return Outcome.EarlyReturn;
                            }
                        }
                        var err = 0;
                        var errOld = 1e10;
                        var hoptde = posneg * hMax;
                        var w = Solver.dim(km);
                        w[1] = 0;
                        var reject = false;
                        var last = false;
                        var atov;
                        var kc = 0;
                        var STATE;
                        (function (STATE) {
                            STATE[STATE["Start"] = 0] = "Start";
                            STATE[STATE["BasicIntegrationStep"] = 1] = "BasicIntegrationStep";
                            STATE[STATE["ConvergenceStep"] = 2] = "ConvergenceStep";
                            STATE[STATE["HopeForConvergence"] = 3] = "HopeForConvergence";
                            STATE[STATE["Accept"] = 4] = "Accept";
                            STATE[STATE["Reject"] = 5] = "Reject";
                        })(STATE || (STATE = {}));
                        var state = STATE.Start;
                        loop: while (true) {
                            _this.debug && console.log('STATE', STATE[state], nStep, xOld, x, h, k, kc, hoptde);
                            switch (state) {
                                case STATE.Start:
                                    atov = false;
                                    // Is xEnd reached in the next step?
                                    if (0.1 * Math.abs(xEnd - x) <= Math.abs(x) * _this.uRound)
                                        break loop;
                                    h = posneg * Math.min(Math.abs(h), Math.abs(xEnd - x), hMax, Math.abs(hoptde));
                                    if ((x + 1.01 * h - xEnd) * posneg > 0) {
                                        h = xEnd - x;
                                        last = true;
                                    }
                                    if (nStep === 0 || !_this.denseOutput) {
                                        F(x, y, dz);
                                        ++nEval;
                                    }
                                    // The first and last step
                                    if (nStep === 0 || last) {
                                        iPt = 0;
                                        ++nStep;
                                        for (var j = 1; j <= k; ++j) {
                                            kc = j;
                                            midex(j);
                                            if (atov)
                                                continue loop;
                                            if (j > 1 && err <= 1) {
                                                state = STATE.Accept;
                                                continue loop;
                                            }
                                        }
                                        state = STATE.HopeForConvergence;
                                        continue;
                                    }
                                    state = STATE.BasicIntegrationStep;
                                    continue;
                                case STATE.BasicIntegrationStep:
                                    // basic integration step
                                    iPt = 0;
                                    ++nStep;
                                    if (nStep >= _this.maxSteps) {
                                        return Outcome.MaxStepsExceeded;
                                    }
                                    kc = k - 1;
                                    for (var j = 1; j <= kc; ++j) {
                                        midex(j);
                                        if (atov) {
                                            state = STATE.Start;
                                            continue loop;
                                        }
                                    }
                                    // convergence monitor
                                    if (k === 2 || reject) {
                                        state = STATE.ConvergenceStep;
                                    } else {
                                        if (err <= 1) {
                                            state = STATE.Accept;
                                        } else if (err > Math.pow(((nj[k + 1] * nj[k]) / 4), 2)) {
                                            state = STATE.Reject;
                                        } else
                                            state = STATE.ConvergenceStep;
                                    }
                                    continue;
                                case STATE.ConvergenceStep:
                                    midex(k);
                                    if (atov) {
                                        state = STATE.Start;
                                        continue;
                                    }
                                    kc = k;
                                    if (err <= 1) {
                                        state = STATE.Accept;
                                        continue;
                                    }
                                    state = STATE.HopeForConvergence;
                                    continue;
                                case STATE.HopeForConvergence:
                                    // hope for convergence in line k + 1
                                    if (err > Math.pow((nj[k + 1] / 2), 2)) {
                                        state = STATE.Reject;
                                        continue;
                                    }
                                    kc = k + 1;
                                    midex(kc);
                                    if (atov)
                                        state = STATE.Start;
                                    else if (err > 1)
                                        state = STATE.Reject;
                                    else
                                        state = STATE.Accept;
                                    continue;
                                case STATE.Accept:
                                    if (!acceptStep(_this.n))
                                        return Outcome.EarlyReturn;
                                    state = STATE.Start;
                                    continue;
                                case STATE.Reject:
                                    k = Math.min(k, kc, km - 1);
                                    if (k > 2 && w[k - 1] < w[k] * _this.stepSizeFac3)
                                        k -= 1;
                                    ++nReject;
                                    h = posneg * hh[k];
                                    reject = true;
                                    state = STATE.BasicIntegrationStep;
                            }
                        }
                        return Outcome.Converged;
                    };
                    var outcome = odxcor();
                    return {
                        y: y.slice(1),
                        outcome: outcome,
                        nStep: nStep,
                        xEnd: xEnd,
                        nAccept: nAccept,
                        nReject: nReject,
                        nEval: nEval
                    };
                };
                return Solver;
            }());
            // return a 1-based array of length n. Initial values undefined.
            Solver.dim = function (n) {
                return Array(n + 1);
            };
            Solver.log10 = function (x) {
                return Math.log(x) / Math.LN10;
            };
            exports.Solver = Solver;

        }, {}],
        2: [function (require, module, exports) {
            
            
            
            
            
            
            
            
            
            
            
            
            
            //Set up Initial Conditions
            function CityParams(beta, gamma, S0, I0) {

                this.beta = beta; // infection rate
                this.gamma = gamma; // Recovery rate (1/days)
                this.S0 = S0; //Initial susceptible fraction
                this.I0 = I0; //Initial infected fraction
                this.R0 = 1 - this.S0 - this.I0; //Initial recovered fraction
                this.y0 = [this.S0, this.I0, this.R0]; //Combining into array
            }
            
            /**Main SIR function containing the key differential equations*/
            var SIR_function = function (cityParams) {
                return function (t, y) {

                    beta = cityParams.beta; //infection rate
                    gamma = cityParams.gamma; //recover rate

                    var S = y[0];
                    var I = y[1];
                    //        var R = y[2];
                    var dSdt = -beta * S * I;
                    var dIdt = beta * S * I - gamma * I;
                    var dRdt = gamma * I;

                    var dydt = [dSdt, dIdt, dRdt];

                    return (dydt);
                }
            }


            /**Differential Equation Solver returns the S,I,R values per shift for a given city and total days and number of shifts per day*/
            function SIRMODEL(cityParams, days, shifts) {

                var odex = require('odex'); //npm module for solving ODEs
                var s = new odex.Solver(3); //initializing the main solver 

                var solutionArray = [] //array to store the S,I,R values (in that order) for each shift

                var init_conditions = cityParams.y0; //Initial S, I, R, values when simulation begins

                for (t = 0; t < days; t += 1 / shifts) {

                    var solution = s.solve(SIR_function(cityParams), 0, init_conditions, t); //Function, initial time, initial conditions, final time
                    solutionArray[t * 2] = solution.y;

                }
                //    console.log(solutionArray);
                return solutionArray; //array with the S,I,R values (in that order) for each shift
            }


            /**Returns the probability of infection for a given location using: 
                a given city's SIR model (SIRMODEL), 
                the day since start of game starting with first day as 1 (e.g., 18)
                the shift on that day (1 or 2), 
                total people at the location (e.g., 30), 
                probability of coming into contact with any person at that location (e.g., 0.3)
                probability of contact with an infected person transmitting infection (usually 0.05)
                */
            function Infection_prob_city(cityDiseaseModel, dayNumber, shiftNumber, totalPeople, contactProb, infectiousContactProb) {

                totalShifts = (dayNumber) * 2 + shiftNumber; //total shifts passed so far, used to calculate index on the Disease Model array
                curr_I = cityDiseaseModel[totalShifts - 1][1] * totalPeople; //Number of infected people in current space
                contacts_I = contactProb * curr_I; //Expected contacts with infected person 
                infected_prob = infectiousContactProb * contacts_I; //Probability of getting infected with those contacts
                return infected_prob;
            }

            /** If you already know how many sick people are in a given location, 
             * then you can use the Infection_prob_local function to calculate the probability of infection at a given location. 
             * This function only needs three inputs:
             * the number of sick people in the location (e.g., 30)
             * probability of coming into contact with any person at that location (e.g., 0.8)
             * probability that contact with an infected person will spread infection (usually 0.05)
             */
            function Infection_prob_local(numSickPeople, contactProb, infectiousContactProb) {
                contacts_I = contactProb * (numSickPeople); //Expected contacts with infected NPCs
                infected_prob = infectiousContactProb * contacts_I; //Probability of getting infected
                
                if(infected_prob>=1){
                    infected_prob = 1;
                }
                if(infected_prob<=0){
                    infected_prob = 0;
                }
                
                return infected_prob;
            }


            /** Returns true or false depending on whether player gets infected */
            function IsInfected(infected_prob) {
                scaledProb = infected_prob * 100;
                var randomNumber = Math.floor(Math.random() * 101);

                if (randomNumber < scaledProb)
                    return true;
                else
                    return false;
            }


            /** Returns the S,I,R values for a given day and shift */
            function CurrentSIRData(cityDiseaseModel, dayNumber, shiftNumber) {
                totalShifts = (dayNumber) * 2 + shiftNumber; //total shifts passed so far, used to calculate index on the Disease Model array
                return [cityDiseaseModel[totalShifts - 1][0], cityDiseaseModel[totalShifts - 1][1], cityDiseaseModel[totalShifts - 1][2]];
            }



            module.exports.CityParams = CityParams;
            module.exports.SIRMODEL = SIRMODEL;
            // module.exports.Infection_prob = Infection_prob;
            module.exports.Infection_prob_city = Infection_prob_city;
            module.exports.Infection_prob_local = Infection_prob_local;
            module.exports.IsInfected = IsInfected;
            module.exports.CurrentSIRData = CurrentSIRData;

            //browserify SIRModel.js --standalone sirModel > /Users/pampam/OneDrive\ -\ Georgia\ Institute\ of\ Technology/DILAC/DILAC-Town/public/schedule/SIRModel.bundle.js 

        }, {
            "odex": 1
        }]
    }, {}, [2])(2)
});