text:
module ProblemSet

using DifferentialEquations
using Memoize
export DEProblem, get_ode_problem,
        get_problem_key, problem_info, problem_plot

struct DEProblem
    fun::Function
    phi::Vector
    bounds::Vector
    data::Vector
    t::AbstractArray
end

Base.copy(p::DEProblem) = DEProblem(p.fun, p.phi, p.bounds,
                                    p.data, p.t)

function problem_info(p::DEProblem)::String
    num_states = length(p.data[1])
    num_params = length(p.phi)
    sep = "\n--- // ---\n"

    rep = sep
    res *= "Number of State Variables:\n\t $(num_states)"
    res *= "\n\t Initial Condition: \n\t $(p.data[1])"
    res *= "\n\t Number of Parameters: \n\t $(num_params) "
    res *= "\n\t Parameter Values: \n\t $(p.phi)"
    res *= "\n\t Bounds: \n\t $(p.bounds)"
    res *= sep
    return res
end

function get_problem_key(i::Int=0)::String
    keys = ["fic","frc",
            "ccgo","bep",
            "mhp","lvp",
            "bbg","fhn",
            "mapk",
            "gosc",
            "tsmp","cho","exponential"]
    if i > 0
        return keys[i]
    else
        return keys
    end
end

@memoize function get_problem(p::String)::DEProblem
    if p == "exponential"
        # ----- Problem 0 -----

        function f_exp(dz_dt, z, phi, t)
            dz_dt[1] = z[1]*phi[1]
        end

        k1 = 1.0
        phi = [k1]
        bounds = [0., 100.]
        ini_cond = [1.]

        t = range(0.0, stop=10.0, length=10)
        de_prob = ODEProblem(f_exp, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u

        return DEProblem(f_exp, phi, bounds, ode_data, t)

    elseif p == "fic"
        # ----- Problem 1 -----

        function fic(dz_dt, z, phi, t)
            r_1 = phi[1]*z[1]
            r_2 = phi[2]*z[2]

            dz_dt[1] = -r_1
            dz_dt[2] = r_1 - r_2
        end

        k1 = 5.0035
        k2 = 1.
        phi = [k1, k2]
        ini_cond = [1., 0.]
        bounds = [[0. for i in 1:length(phi)],
                [10. for i in 1:length(phi)]]
        #=>
        floudas_samples = [
                            1. 0.606 0.368 0.223 0.135 0.082 0.050 0.030 0.018 0.011 0.007;
                            0. 0.373 0.564 0.647 0.669 0.656 0.642 0.583 0.539 0.494 0.451
                            ]
        <=#
        t = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
        de_prob = ODEProblem(fic, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u

        return DEProblem(fic, phi, bounds, ode_data, t)

    elseif p == "frc"
        # ----- Problem 2 -----

        function frc(dz_dt, z, phi, t)
            r_1 = phi[1]*z[1]
            r_2 = phi[2]*z[2]
            r_3 = phi[3]*z[2]
            r_4 = phi[4]*z[3]

            dz_dt[1] = - r_1 + r_2
            dz_dt[2] = r_1 - r_2 - r_3 + r_4
            dz_dt[3] = - r_4 + r_3
        end

        k1 = 4.
        k2 = 2.
        k3 = 40.013
        k4 = 20.007
        phi = [k1, k2, k3, k4]
        ini_cond = [1., 0., 0.]
        bounds = [[0., 0., 10., 10.],
                [10., 10., 50., 50.]]
        #=>
        floudas_samples = [
                            1. 0.8241 0.6852 0.5747 0.4867 0.4166 0.3608 0.3164 0.2810 0.2529 0.2304 0.2126 0.1984 0.1870 0.1870 0.1709 0.1651 0.1606 0.1570 0.1541 0.1518;
                            0. 0.0937 0.1345 0.1654 0.1899 0.2094 0.2249 0.2373 0.2472 0.2550 0.2613 0.2662 0.2702 0.2733 0.2759 0.2779 0.2794 0.2807 0.2817 0.2825 0.2832;
                            0. 0.0821 0.1802 0.2598 0.3233 0.3738 0.4141 0.4461 0.4717 0.4920 0.5082 0.5210 0.5313 0.5395 0.5460 0.5511 0.5553 0.5585 0.5612 0.5632 0.5649
                           ]
        <=#

        t = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0]
        de_prob = ODEProblem(frc, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(frc, phi, bounds, ode_data, t)

    elseif p == "ccgo"
        # ----- Problem 3 -----

        function ccgo(dz_dt, z, phi, t)
            r_1 = phi[1]*z[1]^2
            r_2 = phi[2]*z[2]
            r_3 = phi[3]*z[1]^2

            dz_dt[1] = - r_1 - r_3
            dz_dt[2] = r_1 - r_2
        end

        k1 = 12.214
        k2 = 7.9798
        k3 = 2.2216
        phi = [k1, k2, k3]
        ini_cond = [1., 0.]
        bounds = [[0. for i in 1:length(phi)],
                [20. for i in 1:length(phi)]]
        #=>
        floudas_samples = [
                            1. 0.7307 0.5982 0.4678 0.4267 0.3436 0.3126 0.2808 0.2692 0.2210 0.2122 0.1903 0.1735 0.1615 0.1240 0.1190 0.1109 0.0890 0.0820 0.0745 0.0639;
                            0. 0.1954 0.2808 0.3175 0.3047 0.2991 0.2619 0.2391 0.2210 0.1898 0.1801 0.1503 0.1030 0.0964 0.0581 0.0471 0.0413 0.0367 0.0219 0.0124 0.0089
                            ]
        <=#
        t =[0.0,.025,.05,.075,.1,.125,.150,.175,.2,.225,.25,.3,.35,.4,.45,.5,.55,.65,.75,.85,.95]
        de_prob = ODEProblem(ccgo, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(ccgo, phi, bounds, ode_data, t)

    elseif p == "bep"
        # ----- Problem 4 -----

        function bep(dz_dt, z, phi, t)
            dz_dt[1] = phi[1]*(126.2 - z[1])*(91.9 - z[1])^2 - phi[2]*z[1]^2
        end

        k1 = 4.5704*10^(-6)
        k2 = 2.7845*10^(-4)
        phi = [k1, k2]
        ini_cond = [0.]
        bounds = [[0. for i in 1:length(phi)],
                [.1 for i in 1:length(phi)]]
        #=>
        floudas_samples = [
                            0. 1.4 6.3 10.4 14.2 17.6 21.4 23.0 27.0 30.4 34.4 38.8 41.6 43.5 45.3
                            ]
        <=#
        t = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,9.0,11.0,14.0,19.0,24.0,29.0,39.0]
        de_prob = ODEProblem(bep, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(bep, phi, bounds, ode_data, t)

    elseif p == "mhp"
        # ----- Problem 5 -----

        function mhp(dz_dt, z, phi, t)
            dz_dt[1] = - (2*phi[1] - ((phi[1]*z[2])/((phi[2] + phi[5])*z[1] + z[2])) + phi[3] + phi[4])*z[1]
            dz_dt[2] = ((phi[1]*z[1])*(phi[2]*z[1] - z[2]))/((phi[2] + phi[5])*z[1] + z[2]) + phi[3]*z[1]
            dz_dt[3] = ((phi[1]*z[1])*(z[2] + phi[5]*z[1]))/((phi[2] + phi[5])*z[1] + z[2]) + phi[4]*z[1]
        end

        k1 = 5.2407
        k2 = 1.2176
        k3 = 0.
        k4 = 0.
        k5 = 0.
        phi = [k1, k2, k3, k4, k5]
        ini_cond = [1., 0., 0.]
        bounds = [[0. for i in 1:length(phi)],
                [100. for i in 1:length(phi)]]
        #=>
        floudas_samples = [
                            1. 0.461 0.426 0.383 0.305 0.195 0.170 0.139 0.112  0.112 0.090 0.082 0.066 0.053 0.043 0.041 0.029;
                            0. 0.114 0.135 0.157 0.194 0.231 0.234 0.228 0.228  0.226 0.220 0.214 0.178 0.188 0.183 0.184 0.166;
                            0. 0.018 0.035 0.045 0.047 0.084 0.095 0.111 0.134 0.168 0.148 0.157 0.206 0.206 0.214 0.213 0.230
                            ]
        <=#
        t = [0.,0.050,0.065,0.080,0.123,0.233,0.273,0.354,0.397,0.418,0.502,0.553,0.681,0.750,0.916,0.937,1.122]
        de_prob = ODEProblem(mhp, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(mhp, phi, bounds, ode_data, t)

    elseif p == "lvp"
        # ----- Problem 6 -----

        function lvp(dz_dt, z, phi, t)
            dz_dt[1] = phi[1]*z[1]*(1-z[2])
            dz_dt[2] = phi[2]*z[2]*(z[1]-1)
        end

        k1 = 3.2434
        k2 = 0.9209
        phi = [k1, k2]
        ini_cond = [1.2, 1.1]
        bounds = [[0. for i in 1:length(phi)],
                [10. for i in 1:length(phi)]]
        #=>
        floudas_samples = [
                            1.2 0.7990 0.8731 1.2487 1.0362 0.7483 1.0024 1.2816 0.8944 0.7852 1.1527;
                            1.1 1.0758 0.8711 0.9393 1.1468 1.0027 0.8577 1.0274 1.1369 0.9325 0.9074
                            ]
        <=#
        t = [0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.]
        de_prob = ODEProblem(lvp, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(lvp, phi, bounds, ode_data, t)

    elseif p == "bbg"
        # ----- BBG -----
        """
        Biomass Batch Growth
        Robust and Efficient Parameter Estimation in Dynamic
        Models of Biological Systems
        [ROEF]
        """

        phi = [0.4, 5, 0.05, 0.5]

        bounds = [[10^(-5) for i in 1:length(phi)],
                [10^(2) for i in 1:length(phi)]]

        ini_cond = [2.0, 30.0]
        t = range(0.0, stop=12.0, length=7)

        function bbg(dz_dt, z, phi, t)
            mi, Ks, Kd, yield = phi
            Cb, Cs = z
            dz_dt[1] = mi*(Cs*Cb/(Ks+Cs))-Kd*Cb
            dz_dt[2] = -(mi/yield)*(Cs*Cb/(Ks+Cs))
        end

        de_prob = ODEProblem(bbg, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(bbg, phi, bounds, ode_data, t)

    elseif p == "fhn"
        # ----- FHN -----
        """
        FitzHugh-Nagumo neural spike
        Robust and Efficient Parameter Estimation in Dynamic
        Models of Biological Systems
        [ROEF]
        """

        phi = [0.2, 0.2, 3.0]

        bounds = [[10^(-5) for i in 1:length(phi)],
                [10^(5) for i in 1:length(phi)]]

        ini_cond = [-1.0, 1.0]
        t = range(0.0, stop=20.0, length=7)

        function fhn(dz_dt, z, phi, t)
            a, b, c = phi
            V, R = z
            u = 0
            dz_dt[1] = c*(V - (V^3)/3 + R) + u
            dz_dt[2] = -(1/c)*(V - a + b*R)
        end

        de_prob = ODEProblem(fhn, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(fhn, phi, bounds, ode_data, t)

    elseif p == "mapk"
        # ----- MPK -----
        """
        Kholodenko MAPK signalling pathway (MAPK)
        ROEF
        """
        phi = [2.5, 0.25, 0.75, 0.75, 0.5, 0.5]

        bounds = [[0.01 for i in 1:length(phi)],
                [50.0 for i in 1:length(phi)]]

        ini_cond = [90.0, 10.0, 280.0, 10.0, 10.0, 280.0, 10.0, 10.0]
        #t = range(0.0, stop=20.0, length=7)
        t = [50.0, 100.0, 150.0, 200.0, 300.0,
            400.0, 500.0, 600.0, 800.0, 1000.0]

        function mapk(dz_dt, z, phi, t)
            J0V1, J1V2, J4V5, J5V6, J8V9, J9V10 = phi
            J0Ki, J0n, J0k1, J1KK2, J2k3, J2KK3, J3k4,
            J3KK4, J4KK5, J5KK6, J6k7, J6KK7, J7k8, J7KK8,
            J8KK9, J9KK10 = convert(Array{eltype(phi)},
                [9.0, 1.0, 10.0, 8.0, 0.025, 15.0, 0.025, 15.0, 15.0,
                15.0, 0.025, 15.0, 0.025, 15.0, 15.0, 15.0])

            RJ0 = J0V1*(z[1]/(1+((z[8]/J0Ki)^J0n)*(J0k1+z[1])))
            RJ1 = J1V2*(z[2]/(J1KK2+z[2]))
            RJ2 = J2k3*(z[2]*z[3]/(J2KK3+z[3]))
            RJ3 = J3k4*(z[2]*z[4]/(J3KK4+z[4]))
            RJ4 = J4V5*(z[5]/(J4KK5+z[5]))
            RJ5 = J5V6*(z[4]/(J5KK6+z[4]))
            RJ6 = J6k7*(z[5]*z[6]/(J6KK7+z[6]))
            RJ7 = J7k8*(z[5]*z[7]/(J7KK8+z[7]))
            RJ8 = J8V9*(z[8]/(J8KK9+z[8]))
            RJ9 = J9V10*(z[7]/(J9KK10+z[7]))

            dz_dt[1] = - RJ0 + RJ1
            dz_dt[2] = RJ0 - RJ1
            dz_dt[3] = - RJ2 + RJ5
            dz_dt[4] = RJ2 - RJ3 + RJ4 - RJ5
            dz_dt[5] = RJ3 - RJ4
            dz_dt[6] = - RJ6 + RJ9
            dz_dt[7] = RJ6 - RJ7 + RJ8 - RJ9
            dz_dt[8] = RJ7 - RJ8
        end

        de_prob = ODEProblem(mapk, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(mapk, phi, bounds, ode_data, t)

    elseif p == "gosc"
        # ----- GOsc -----
        """
        Goodwin Oscillator
        Robust and Efficient Parameter Estimation in Dynamic
        Models of Biological Systems
        [ROEF]
        """
        phi = [1.0, 0.1, 1.0, 0.1, 1.0, 0.1, 1.0, 10.0]

        bounds = [[10^(-2) for i in 1:length(phi)],
                [10^(2) for i in 1:length(phi)]]
        bounds[1][end] = 1.0
        bounds[2][end] = 12.0

        ini_cond = [0.1, 0.2, 2.5]
        t = range(0.0, stop=240.0, length=10)

        function gosc(dz_dt, z, phi, t)
            k1, k2, k3, k4, k5, k6, Ki, n = phi

            dz_dt[1] = k1*Ki^n/(Ki^n + z[3]^n) - k2*z[1]
            dz_dt[2] = k3*z[1] - k4*z[2]
            dz_dt[3] = k5*z[2] - k6*z[3]
        end

        de_prob = ODEProblem(gosc, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        new_prob = DEProblem(gosc, phi, bounds, ode_data, t)

        #=>
        # ----- TGFB -----
        """
        TGF - β signalling pathway model
        Robust and Efficient Parameter Estimation in Dynamic
        Models of Biological Systems
        [ROEF]
        """
        phi = [0.00015, 0.023, 0.01, 0.01, 0.01, 0.1, 0.000403, 0.0026, 0.0056, 0.002, 0.016, 5.7, 0.00657, 0.0017, 1.0, 0.0008, 0.0001, 0.0021, 0.001, 9000.0, 1800.0]

        bounds = [[10^(-3) for i in 1:length(phi)],
                [10^(3) for i in 1:length(phi)]]
        bounds[:,end] = [1, 12]

        ini_cond = [0.1, 0.2, 2.5]
        t = range(0.0, stop=240.0, length=10)

        function f_gosc(dz_dt, z, phi, t)
            CTGFb, CTGFbR, CTGFbTGFbR, CTGFbTGFbRP, CISmadTGFbTGFbRP, CSmad, CSmadP, CCoSmad, CSmadPSmaP, CSmadPCoSmad, CSmadN, CSmadPSmadPN, CSmadPN, CSmadPCoSmadN, CCoSmadN, CISmadmRNA1, CISmadmRNA2, CISmad = z
            k1, k2, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19 = phi
            k3 = 0.01
            k20 = 9000.0
            k21 = 1800.0
            r1 = k1*CTGFbTGFbR
            r2 = k2*CTGFbR*CTGFb
            r3 = k3*CTGFbTGFbR*(1-e^(-((t-k20)/(k21)))^10)
            r4 = k4*CTGFbTGFbRP
            r5 = k5*CTGFbTGFbRP*CISmad
            r6 = k6*CISmadTGFbTGFbRP
            r7 = k7*C Smad*CTGFbTGFbRP
            r8 = k8*CSmad
            r9 = k9*CSmadN
            r10 = k10*2*CSmadP*CSmadP
            r11 = k11*CSmadPCoSmad
            r12 = k12*CSmadP*CCoSmad
            r13 = k11*CSmadP*CoSmad
            r14 = k8*CCoSmad
            r15 = k9*CCoSmadN
            r16 = k12*k8*CSmadPSmaP
            r17 = k8*CSmadP
            r18 = k9*CSmadPN
            r19 = k12*k8*CSmadPCoSmad
            r20 = k13*CSmadPN
            r21 = k10*2*CSmadPNCSmadPN
            r22 = k11*CSmadPSmadPN
            r23 = k10*CSmadPN*CCoSmadN
            r24 = k11*CSmadPCoSmadN
            r25 = k14*((CSmadPCoSmadN^2)\(CSmadPCoSmadN^2+k15^2))
            r26 = k16*CISmadmRNA1
            r27 = k17*CISmadmRNA2
            r28 = k18*CISmadmRNA2
            r29 = k19*CISmad

            dz_dt[1] = r1 - r2
            dz_dt[1] = r1 - r2
            dz_dt[1] = - r1 + r2 - r3 + r4 + r6
            dz_dt[1] = r3 - r4 - r5
            dz_dt[1] = r5 - r6
            dz_dt[1] = - r7 - r8 + r9
            dz_dt[1] = r7 - r10 + r11 - r12 + r13 - r17 + r18
            dz_dt[1] = - r12 + r13 - r14 + r15
            dz_dt[1] = r10 - r11 - r16
            dz_dt[1] = r12 - r13 - r19
            dz_dt[1] = r8 - r9 + r20
            dz_dt[1] = r16 + r21 - r22
            dz_dt[1] = r17 - r18 - r20 - r21 + r22 - r23 + r24
            dz_dt[1] = r19 + r23 - r24
            dz_dt[1] = r14 - r15 - r23 + r24
            dz_dt[1] = r25 - r26
            dz_dt[1] = r26 - r27
            dz_dt[1] = r28 - r29 - r5 + r6
        end

        de_prob = ODEProblem(f_fhn, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, OwrenZen3(), saveat=t)
        ode_data = de_sol.u

        new_prob = DEProblem(f_gosc, phi, bounds, ode_data, t)
        push!(problem_set, new_prob)
        <=#

    elseif p == "tsmp"
        # ----- TSP -----
        """
        Three Step Metabolic Pathway from:
        Parameter Estimation in Biochemical Pathways:
        A Comparison of Global Optimization Methods
        [BioParamEst-GO]
        [ROEF]
        """

        s = [.1, .46416, 2.1544, 10]
        p = [.05, .13572, .36840, .1]

        phi = [1.0, 1.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0,
                1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 2.0, 1.0,
                .1, 1.0, .1, .1, 1.0, .1, .1, 1.0, .1, 1.0,
                1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

        bounds = [[10^(-5) for i in 1:length(phi)],
                [10^(5) for i in 1:length(phi)]]
        for i in [2,4,8,10,14,16]
            bounds[1][i] = 0.1
            bounds[2][i] = 10
        end

        ini_cond = [.6667, .5725, .4176, .4, .3641, .2946, 1.419, .9346]
        t = range(0, stop=120, length=21)

        function tsmp(dz_dt, z, phi, t)
            G1 = z[1]
            G2 = z[2]
            G3 = z[3]

            E1 = z[4]
            E2 = z[5]
            E3 = z[6]

            M1 = z[7]
            M2 = z[8]

            S = s[3]
            P = p[4]

            V1 = phi[1]
            Ki1 = phi[2]
            ni1 = phi[3]
            Ka1 = phi[4]
            na1 = phi[5]
            k1 = phi[6]
            V2 = phi[7]
            Ki2 = phi[8]
            ni2 = phi[9]
            Ka2 = phi[10]
            na2 = phi[11]
            k2 = phi[12]
            V3 = phi[13]
            Ki3 = phi[14]
            ni3 = phi[15]
            Ka3 = phi[16]
            na3 = phi[17]
            k3 = phi[18]
            V4 = phi[19]
            K4 = phi[20]
            k4 = phi[21]
            V5 = phi[22]
            K5 = phi[23]
            k5 = phi[24]
            V6 = phi[25]
            K6 = phi[26]
            k6 = phi[27]
            kcat1 = phi[28]
            Km1 = phi[29]
            Km2 = phi[30]
            kcat2 = phi[31]
            Km3 = phi[32]
            Km4 = phi[33]
            kcat3 = phi[34]
            Km5 = phi[35]
            Km6 = phi[36]

            dz_dt[1] = (V1/(1+(P/Ki1)^ni1+(Ka1/S)^na1))-k1*G1
            dz_dt[2] = (V2/(1+(P/Ki2)^ni2+(Ka2/M1)^na2))-k2*G2
            dz_dt[3] = (V3/(1+(P/Ki3)^ni3+(Ka3/M2)^na3))-k3*G3

            dz_dt[4] = (V4*G1/(K4+G1))-k4*E1
            dz_dt[5] = (V5*G2/(K5+G2))-k5*E2
            dz_dt[6] = (V6*G3/(K6+G3))-k6*E3

            dz_dt[7] = ((kcat1*E1*(1/Km1)*(S-M1)))/(1+(S/Km1)+(M1/Km2)) - ((kcat2*E2*(1/Km3)*(M1-M2)))/(1+(M1/Km3)+(M2/Km4))
            dz_dt[8] = ((kcat2*E2*(1/Km3)*(M1-M2)))/(1+(M1/Km3)+(M2/Km4)) - ((kcat3*E3*(1/Km5)*(M2-P)))/(1+(M2/Km5)+(P/Km6))
        end

        de_prob = ODEProblem(tsmp, ini_cond, (t[1],t[end]), phi)
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(tsmp, phi, bounds, ode_data, t)

    elseif p == "cho"
    # ----- CHO -----
    # [BioPreDyn]

        x0 = zeros(Float64, 35)
        x0[1] = 5
        x0[2] = 1
        x0[3] = 5
        x0[4] = 5
        x0[5] = 1
        x0[6] = 1
        x0[7] = 1
        x0[8] = 1
        x0[9] = 1
        x0[10] = 1
        x0[11] = 1
        x0[12] = 1
        x0[13] = 1
        x0[14] = 1
        x0[15] = 1
        x0[16] = 1
        x0[17] = 1
        x0[18] = 1
        x0[19] = 1
        x0[20] = 1
        x0[21] = 1
        x0[22] = 1
        x0[23] = 1
        x0[24] = 1
        x0[25] = 1
        x0[26] = 1
        x0[27] = 1
        x0[28] = 1
        x0[29] = 1
        x0[30] = 1
        x0[31] = 1
        x0[32] = 1
        x0[33] = 1
        x0[34] = 1
        x0[35] = 6

        p = zeros(Float64, 117)
        p[1] = 1000
        p[2] = 1000
        p[3] = 1000
        p[4] = 1000
        p[5] = 1000
        p[6] = 1000
        p[7] = 1000
        p[8] = 1000
        p[9] = 1
        p[10] = 1
        p[11] = 0.7
        p[12] = 0.7
        p[13] = 0.7
        p[14] = 0.7
        p[15] = 0.2
        p[16] = 0.2
        p[17] = 1
        p[18] = 1
        p[19] = 0.7
        p[20] = 0.7
        p[21] = 0.05
        p[22] = 1
        p[23] = 0.5
        p[24] = 0.7
        p[25] = 0.5
        p[26] = 2
        p[27] = 0.7
        p[28] = 0.7
        p[29] = 2
        p[30] = 0.2
        p[31] = 0.2
        p[32] = 1
        p[33] = 0.7
        p[34] = 0.2
        p[35] = 0.21
        p[36] = 1
        p[37] = 0.7
        p[38] = 1
        p[39] = 0.2
        p[40] = 0.2
        p[41] = 0.05
        p[42] = 0.01
        p[43] = 0.7
        p[44] = 1
        p[45] = 0.5
        p[46] = 1
        p[47] = 0.04
        p[48] = 0.01
        p[49] = 0.01
        p[50] = 0.01
        p[51] = 1
        p[52] = 1
        p[53] = 0.5
        p[54] = 1
        p[55] = 1
        p[56] = 0.5
        p[57] = 0.1
        p[58] = 0.1
        p[59] = 0.5
        p[60] = 1
        p[61] = 1
        p[62] = 1
        p[63] = 0.3
        p[64] = 1
        p[65] = 1
        p[66] = 0.7
        p[67] = 0.4
        p[68] = 0.7
        p[69] = 1
        p[70] = 0.7
        p[71] = 0.7
        p[72] = 2
        p[73] = 0.2
        p[74] = 2
        p[75] = 1
        p[76] = 0.7
        p[77] = 0.7
        p[78] = 0.2
        p[79] = 0.7
        p[80] = 0.5
        p[81] = 0.21
        p[82] = 0.2
        p[83] = 0.2
        p[84] = 0.01
        p[85] = 2
        p[86] = 2
        p[87] = 2
        p[88] = 2
        p[89] = 2
        p[90] = 2
        p[91] = 2
        p[92] = 2
        p[93] = 1
        p[94] = 1
        p[95] = 1
        p[96] = 1
        p[97] = 1
        p[98] = 0.5
        p[99] = 0.7
        p[100] = 0.7
        p[101] = 2
        p[102] = 0.2
        p[103] = 2
        p[104] = 1
        p[105] = 1
        p[106] = 0.7
        p[107] = 1
        p[108] = 0.7
        p[109] = 1
        p[110] = 0.7
        p[111] = 1
        p[112] = 0.7
        p[113] = 2
        p[114] = 2
        p[115] = 0.7
        p[116] = 0.2
        p[117] = 0.2

        t = collect(1.0:20.0)
        bounds = [[10^-5 for i in 1:length(p)],
                [10^5 for i in 1:length(p)]]

        function cho(dx_dt, x, par, t)
            p = zeros(347)
            p[118] = 1
            p[119] = 1
            p[120] = 1
            p[121] = 1
            p[122] = 1
            p[123] = 1
            p[124] = 1
            p[125] = 1
            p[126] = 1
            p[127] = 1
            p[128] = 1
            p[129] = 1
            p[130] = 1
            p[131] = 1
            p[132] = 1
            p[133] = 1
            p[134] = 1
            p[135] = 1
            p[136] = 1
            p[137] = 1
            p[138] = 1
            p[139] = 1
            p[140] = 1
            p[141] = 1
            p[142] = 1
            p[143] = 1
            p[144] = 1
            p[145] = 1
            p[146] = 1
            p[147] = 1
            p[148] = 1
            p[149] = 1
            p[150] = 0
            p[151] = 141.471
            p[152] = 0.9
            p[153] = 0.1
            p[154] = 83402
            p[155] = 127877
            p[156] = 603.414
            p[157] = 603.414
            p[158] = 5.02845
            p[159] = 21722.9
            p[160] = 603.414
            p[161] = 361200
            p[162] = 673524
            p[163] = 1.62679e+06
            p[164] = 164390
            p[165] = 36512.9
            p[166] = 166804
            p[167] = 39530
            p[168] = 1810.24
            p[169] = 127877
            p[170] = 166201
            p[171] = 1.94454e+06
            p[172] = 5.02845
            p[173] = 317755
            p[174] = 0
            p[175] = 532531
            p[176] = 38926.6
            p[177] = 40133.4
            p[178] = 689982
            p[179] = 4.45873e+06
            p[180] = 2413.66
            p[181] = 83402
            p[182] = 603.414
            p[183] = 603.414
            p[184] = 2413.66
            p[185] = 533135
            p[186] = 100000
            p[187] = 1
            p[188] = 1000
            p[189] = 1000
            p[190] = 30
            p[191] = 1000
            p[192] = 1000
            p[193] = 1000
            p[194] = 1000
            p[195] = 1000
            p[196] = 1000
            p[197] = 1000
            p[198] = 3000
            p[199] = 1000
            p[200] = 1000
            p[201] = 1000
            p[202] = 100
            p[203] = 100
            p[204] = 100
            p[205] = 100
            p[206] = 1000
            p[207] = 1000
            p[208] = 1000
            p[209] = 1000
            p[210] = 1000
            p[211] = 1000
            p[212] = 1000
            p[213] = 1000
            p[214] = 1000
            p[215] = 3000
            p[216] = 100000
            p[217] = 1000
            p[218] = 1000
            p[219] = 1000
            p[220] = 1000
            p[221] = 1000
            p[222] = 1000
            p[223] = 1000
            p[224] = 1000
            p[225] = 1000
            p[226] = 1000
            p[227] = 1000
            p[228] = 0
            p[229] = 1
            p[230] = 0
            p[231] = 0
            p[232] = 1
            p[233] = 1
            p[234] = 1
            p[235] = -0.7
            p[236] = -0.7
            p[237] = 0.7
            p[238] = 0.7
            p[239] = -0.2
            p[240] = -0.2
            p[241] = 1
            p[242] = 1
            p[243] = -0.7
            p[244] = -0.7
            p[245] = 0.05
            p[246] = 1
            p[247] = 0.5
            p[248] = -0.7
            p[249] = -0.5
            p[250] = 2
            p[251] = 0.7
            p[252] = 0.7
            p[253] = -2
            p[254] = -0.2
            p[255] = -0.2
            p[256] = 1
            p[257] = 0.7
            p[258] = -0.2
            p[259] = -0.21
            p[260] = 0
            p[261] = 1
            p[262] = 0.7
            p[263] = 1
            p[264] = -0.2
            p[265] = -0.2
            p[266] = 0.05
            p[267] = 0
            p[268] = -0.01
            p[269] = 0.7
            p[270] = 1
            p[271] = -0.5
            p[272] = -1
            p[273] = 0.04
            p[274] = -0.01
            p[275] = -0.01
            p[276] = -0.01
            p[277] = 1
            p[278] = 1
            p[279] = 0.5
            p[280] = -1
            p[281] = -1
            p[282] = -0.5
            p[283] = 0.1
            p[284] = -0.1
            p[285] = 0.5
            p[286] = 1
            p[287] = 1
            p[288] = -1
            p[289] = -0.3
            p[290] = -1
            p[291] = 1
            p[292] = 0.7
            p[293] = -0.4
            p[294] = -0.7
            p[295] = 1
            p[296] = -0.7
            p[297] = 0.7
            p[298] = 2
            p[299] = -0.2
            p[300] = -2
            p[301] = 1
            p[302] = 0.7
            p[303] = 0.7
            p[304] = 0.2
            p[305] = 0.7
            p[306] = -0.5
            p[307] = -0.21
            p[308] = -0.2
            p[309] = -0.2
            p[310] = 0
            p[311] = 0
            p[312] = -0.01
            p[313] = 0.7
            p[314] = -0.2
            p[315] = 2
            p[316] = 2
            p[317] = -2
            p[318] = -2
            p[319] = 2
            p[320] = 2
            p[321] = -2
            p[322] = -2
            p[323] = 1
            p[324] = 1
            p[325] = -1
            p[326] = -1
            p[327] = 1
            p[328] = -0.5
            p[329] = 0.7
            p[330] = 0.7
            p[331] = 2
            p[332] = -0.2
            p[333] = -2
            p[334] = 1
            p[335] = 1
            p[336] = -0.7
            p[337] = 1
            p[338] = -0.7
            p[339] = 1
            p[340] = -0.7
            p[341] = 1
            p[342] = -0.7
            p[343] = 2
            p[344] = -2
            p[345] = 0.7
            p[346] = -0.2
            p[347] = -0.2

            for i in 1:length(par)
                p[i] = par[i]
            end

            y = zeros(promote_type(eltype(x),eltype(t)),102)
            y[1] = log(x[1])
            y[2] = log(x[2])
            y[3] = log(x[3])
            y[4] = log(x[4])
            y[5] = log(x[5])
            y[6] = log(x[6])
            y[7] = log(x[7])
            y[8] = log(x[8])
            y[9] = log(x[9])
            y[10] = log(x[10])
            y[11] = log(x[11])
            y[12] = log(x[12])
            y[13] = log(x[13])
            y[14] = log(x[14])
            y[15] = log(x[15])
            y[16] = log(x[16])
            y[17] = log(x[17])
            y[18] = log(x[18])
            y[19] = log(x[19])
            y[20] = log(x[20])
            y[21] = log(x[21])
            y[22] = log(x[22])
            y[23] = log(x[23])
            y[24] = log(x[24])
            y[25] = log(x[25])
            y[26] = log(x[26])
            y[27] = log(x[27])
            y[28] = log(x[28])
            y[29] = log(x[29])
            y[30] = log(x[30])
            y[31] = log(x[31])
            y[32] = log(x[32])
            y[33] = log(x[33])
            y[34] = log(x[34])
            y[35] = log(x[35])
            y[36] = x[1]*p[186]
            y[37] = x[2]*p[187]
            y[38] = x[3]*p[188]
            y[39] = x[4]*p[189]
            y[40] = x[5]*p[190]
            y[41] = x[6]*p[191]
            y[42] = x[7]*p[192]
            y[43] = x[8]*p[193]
            y[44] = x[9]*p[194]
            y[45] = x[10]*p[195]
            y[46] = x[11]*p[196]
            y[47] = x[12]*p[197]
            y[48] = x[13]*p[198]
            y[49] = x[14]*p[199]
            y[50] = x[15]*p[200]
            y[51] = x[16]*p[201]
            y[52] = x[17]*p[202]
            y[53] = x[18]*p[203]
            y[54] = x[19]*p[204]
            y[55] = x[20]*p[205]
            y[56] = x[21]*p[206]
            y[57] = x[22]*p[207]
            y[58] = x[23]*p[208]
            y[59] = x[24]*p[209]
            y[60] = x[25]*p[210]
            y[61] = x[26]*p[211]
            y[62] = x[27]*p[212]
            y[63] = x[28]*p[213]
            y[64] = x[29]*p[214]
            y[65] = x[30]*p[215]
            y[66] = x[31]*p[216]
            y[67] = x[32]*p[217]
            y[68] = x[33]*p[218]
            y[69] = x[34]*p[219]
            if t>5 && t<10
                y[70] = p[150]
            else
                y[70] = 0.00000000000000000
            end
            y[71] = 0.00000000000000000
            y[72] = 0.00000000000000000
            y[73] = 0.00000000000000000
            y[74] = 0.00000000000000000
            y[75] = 0.00000000000000000
            y[76] = p[123]*p[159]*(1.00000000000000000+p[233]*y[8]+p[234]*y[9]+p[235]*y[6]+p[236]*y[7])
            y[77] = p[124]*p[160]*(1.00000000000000000+p[237]*y[12]+p[238]*y[13]+p[239]*y[10]+p[240]*y[11])
            y[78] = p[125]*p[161]*(1.00000000000000000+p[241]*y[15]+p[242]*y[6]+p[243]*y[8]+p[244]*y[14]+p[245]*y[16])
            y[79] = p[126]*p[162]*(1.00000000000000000+p[246]*y[16]+p[247]*y[7]+p[248]*y[15]+p[249]*y[9])
            y[80] = p[127]*p[163]*(1.00000000000000000+p[250]*y[9]+p[251]*y[19]+p[252]*y[20]+p[253]*y[7]+p[254]*y[17]+p[255]*y[18])
            y[81] = p[128]*p[164]*(1.00000000000000000+p[256]*y[22]+p[257]*y[11]+p[258]*y[21]+p[259]*y[13]+p[260]*y[13])
            y[82] = p[129]*p[165]*(1.00000000000000000+p[261]*y[15]+p[262]*y[7]+p[263]*y[21]+p[264]*y[8]+p[265]*y[9]+p[266]*y[32]+p[267]*y[9]+p[268]*y[30])
            y[83] = p[130]*p[166]*(1.00000000000000000+p[269]*y[25]+p[270]*y[26]+p[271]*y[23]+p[272]*y[24]+p[273]*y[11]+p[274]*y[11]+p[275]*y[13]+p[276]*y[22])
            y[84] = p[131]*p[167]*(1.00000000000000000+p[277]*y[28]+p[278]*y[29]+p[279]*y[23]+p[280]*y[12]+p[281]*y[27]+p[282]*y[25]+p[283]*y[27]+p[284]*y[10])
            y[85] = p[132]*p[168]*(1.00000000000000000+p[285]*y[32]+p[286]*y[22]+p[287]*y[16]+p[288]*y[15]+p[289]*y[30]+p[290]*y[27])
            y[86] = p[133]*p[169]*(1.00000000000000000+p[291]*y[21]+p[292]*y[23]+p[293]*y[2]+p[294]*y[25])
            y[87] = p[134]*p[170]*(1.00000000000000000+p[295]*y[24]+p[296]*y[22])
            y[88] = p[135]*p[171]*(1.00000000000000000+p[297]*y[20]+p[298]*y[17]+p[299]*y[18]+p[300]*y[19])
            y[89] = p[136]*(p[172]/(p[209]/p[220]*(p[210]/p[221])*(p[197]/p[222])*(p[218]/p[223])*(p[219]/p[224])*(p[214]/p[225])*(p[195]/p[226])*(p[198]/p[227])/((1.00000000000000000+p[209]/p[220])*(1.00000000000000000+p[210]/p[221])*(1.00000000000000000+p[197]/p[222])*(1.00000000000000000+p[218]/p[223])*(1.00000000000000000+p[219]/p[224])*(1.00000000000000000+p[214]/p[225])*(1.00000000000000000+p[195]/p[226])*(1.00000000000000000+p[198]/p[227]))))*(x[24]*p[209]/p[220]*(x[25]*p[210]/p[221])*(x[12]*p[197]/p[222])*(x[33]*p[218]/p[223])*(x[34]*p[219]/p[224])*(x[29]*p[214]/p[225])*(x[10]*p[195]/p[226])*(x[13]*p[198]/p[227])/((1.00000000000000000+x[24]*p[209]/p[220])*(1.00000000000000000+x[25]*p[210]/p[221])*(1.00000000000000000+x[12]*p[197]/p[222])*(1.00000000000000000+x[33]*p[218]/p[223])*(1.00000000000000000+x[34]*p[219]/p[224])*(1.00000000000000000+x[29]*p[214]/p[225])*(1.00000000000000000+x[10]*p[195]/p[226])*(1.00000000000000000+x[13]*p[198]/p[227])))
            y[90] = p[137]*p[173]*(1.00000000000000000+p[301]*y[8]+p[302]*y[7]+p[303]*y[19]+p[304]*y[31]+p[305]*y[32]+p[306]*y[16]+p[307]*y[9]+p[308]*y[17]+p[309]*y[30]+p[310]*y[31]+p[311]*y[9]+p[312]*y[15])
            y[91] = p[138]*p[174]*(1.00000000000000000+p[313]*y[18]+p[314]*y[20])
            y[92] = p[139]*p[175]*(1.00000000000000000+p[315]*y[11]+p[316]*y[30]+p[317]*y[32]+p[318]*y[13])
            y[93] = p[140]*p[176]*(1.00000000000000000+p[319]*y[8]+p[320]*y[27]+p[321]*y[28]+p[322]*y[16])
            y[94] = p[141]*p[177]*(1.00000000000000000+p[323]*y[14]+p[324]*y[12]+p[325]*y[29]+p[326]*y[6])
            y[95] = p[142]*p[178]*(1.00000000000000000+p[327]*y[13]+p[328]*y[11])
            y[96] = p[143]*p[179]*(1.00000000000000000+p[329]*y[32]+p[330]*y[31]+p[331]*y[18]+p[332]*y[30]+p[333]*y[20])
            y[97] = p[144]*p[180]*(1.00000000000000000+p[334]*y[27]+p[335]*y[31]+p[336]*y[16])
            y[98] = p[145]*p[181]*(1.00000000000000000+p[337]*y[1]+p[338]*y[26])
            y[99] = p[146]*p[182]*(1.00000000000000000+p[339]*y[3]+p[340]*y[33])
            y[100] = p[147]*p[183]*(1.00000000000000000+p[341]*y[4]+p[342]*y[34])
            y[101] = p[148]*p[184]*(1.00000000000000000+p[343]*y[6]+p[344]*y[12])
            y[102] = p[149]*p[185]*(1.00000000000000000+p[345]*y[18]+p[346]*y[31]+p[347]*y[20])

            dx_dt[1] = (y[71]-y[98]*1.00000000000000000/p[151]-x[1]*p[186]/x[35]*y[70])/p[186]
            dx_dt[2] = ((-y[72])+y[86]*1.00000000000000000/p[151]-x[2]*p[187]/x[35]*y[70])/p[187]
            dx_dt[3] = (y[73]-y[99]*1.00000000000000000/p[151]-x[3]*p[188]/x[35]*y[70])/p[188]
            dx_dt[4] = (y[74]-y[100]*1.00000000000000000/p[151]-x[4]*p[189]/x[35]*y[70])/p[189]
            dx_dt[5] = ((-y[75])+y[89]*1.00000000000000000/p[151]-x[5]*p[190]/x[35]*y[70])/p[190]
            dx_dt[6] = (y[76]-y[78]+y[94]*p[152]/p[153]-y[101]*p[152]/p[153])/p[191]
            dx_dt[7] = (y[76]-y[79]+y[80]-2.00000000000000000*y[82]*p[152]/p[153]-y[90])/p[192]
            dx_dt[8] = ((-y[76])+y[78]+y[82]*p[152]/p[153]-y[90]-y[93]*p[152]/p[153])/p[193]
            dx_dt[9] = ((-y[76])+y[79]-y[80]+2.00000000000000000*y[82]*p[152]/p[153]+y[90])/p[194]
            dx_dt[10] = (y[77]-120.00000000000000000*y[89])/p[195]
            dx_dt[11] = (y[77]-y[81]+1260.00000000000000000*y[89]-y[92]+y[95])/p[196]
            dx_dt[12] = ((-y[77])+y[84]-240.00000000000000000*y[89]-y[94]+y[101])/p[197]
            dx_dt[13] = ((-y[77])+y[81]-1260.00000000000000000*y[89]+y[92]-y[95])/p[198]
            dx_dt[14] = (y[78]-y[94]*p[152]/p[153])/p[199]
            dx_dt[15] = ((-y[78])+y[79]-y[82]*p[152]/p[153]+y[85]*p[152]/p[153])/p[200]
            dx_dt[16] = ((-y[79])-y[85]*p[152]/p[153]+y[90]+y[93]*p[152]/p[153]+y[97]*p[152]/p[153])/p[201]
            dx_dt[17] = (2.00000000000000000*y[80]-2.00000000000000000*y[88]+2.00000000000000000*y[90])/p[202]
            dx_dt[18] = (4.00000000000000000*y[80]+6.00000000000000000*y[88]-y[91]-3.00000000000000000*y[96]-y[102]*p[152]/p[153])/p[203]
            dx_dt[19] = ((-2.00000000000000000)*y[80]+2.00000000000000000*y[88]-2.00000000000000000*y[90])/p[204]
            dx_dt[20] = ((-4.00000000000000000)*y[80]-6.00000000000000000*y[88]+y[91]+3.00000000000000000*y[96]+y[102]*p[152]/p[153])/p[205]
            dx_dt[21] = (y[81]-y[82]-y[86])/p[206]
            dx_dt[22] = ((-y[81])-y[85]+y[87])/p[207]
            dx_dt[23] = (y[83]-y[84]-y[86]+120.00000000000000000*y[89])/p[208]
            dx_dt[24] = (y[83]-y[87]-120.00000000000000000*y[89])/p[209]
            dx_dt[25] = ((-y[83])+y[84]+y[86]-120.00000000000000000*y[89])/p[210]
            dx_dt[26] = ((-0.50000000000000000)*y[83]+y[98])/p[211]
            dx_dt[27] = (y[84]+y[85]-y[93]-y[97])/p[212]
            dx_dt[28] = ((-y[84])+120.00000000000000000*y[89]+y[93])/p[213]
            dx_dt[29] = ((-y[84])-120.00000000000000000*y[89]+y[94])/p[214]
            dx_dt[30] = (y[85]*p[152]/p[153]+y[90]-y[92]*p[152]/p[153]+y[96])/p[215]
            dx_dt[31] = ((-y[90])-y[96]-y[97]*p[152]/p[153]+y[102]*p[152]/p[153])/p[216]
            dx_dt[32] = ((-y[85])*p[152]/p[153]-y[90]+y[92]*p[152]/p[153]-y[96])/p[217]
            dx_dt[33] = ((-120.00000000000000000)*y[89]+y[99])/p[218]
            dx_dt[34] = ((-120.00000000000000000)*y[89]+y[100])/p[219]
            dx_dt[35] = y[70]
        end

        de_prob = ODEProblem(cho, x0, (t[1],t[end]), p,
                       isoutofdomain=(u,p,t) -> any(x -> (x < 0), vcat(u,p)))
        de_sol = solve(de_prob, AutoVern7(Rodas5()), saveat=t)
        ode_data = de_sol.u
        return DEProblem(cho, p, bounds, ode_data, t)

    end
end

end

max_stars_repo_name: atgmello/DEParamEst

max_issues_repo_name: atgmello/DEParamEst

max_forks_repo_name: atgmello/DEParamEst


text:
function solve(
    prob::Union{DiffEqBase.AbstractODEProblem{uType,tType,isinplace},
                DiffEqBase.AbstractSDEProblem{uType,tType,isinplace}},
    alg::AlgType,
    timeseries=[],ts=[],ks=[];
    verbose=true,
    dt = nothing,
    timeseries_errors=true,
    callback=nothing,kwargs...) where {uType,tType,isinplace,AlgType<:BridgeAlgorithm}

    if dt == nothing
        error("dt required for fixed timestep methods.")
    end

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        warned && warn_compat()
    end

    if callback != nothing || :callback in keys(prob.kwargs)
        error("Bridge is not compatible with callbacks.")
    end

    u0 = prob.u0
    p = prob.p

    if isinplace
        f = (t,u) -> (du = similar(u); prob.f(du,u,p,t); Diagonal(du))
        if typeof(prob) <: DiffEqBase.AbstractSDEProblem
            g = (t,u) -> (du = similar(u); prob.g(du,u,p,t,u); Diagonal(du))
        end
    else
        f = (t,u) -> prob.f(u,p,t)
        if typeof(prob) <: DiffEqBase.AbstractSDEProblem
            if typeof(u0) <: Number
                g = (t,u) -> prob.g(u,p,t)
            else
                g = (t,u) -> Diagonal(prob.g(u,p,t))
            end
        end
    end

    t = prob.tspan[1]:dt:prob.tspan[2]
    if typeof(prob) <: DiffEqBase.AbstractSDEProblem
        W = Bridge.samplepath(t, zero(u0))
        samp = Bridge.sample(t, Bridge.Wiener{typeof(u0)}())
        if typeof(alg) <: BridgeEuler
            u = Bridge.solve!(Bridge.Euler(), W, u0, samp, (f,g))
        elseif typeof(alg) <: BridgeHeun
            u = Bridge.solve!(Bridge.StochasticHeun(), W, u0, samp, (f,g))
        elseif typeof(alg) <: BridgeSRK && typeof(u0) <: Number
            u = Bridge.solve!(Bridge.StochasticRungeKutta(), W, u0, samp, (f,g))
        else
            error("BridgeSRK is not compatible with non Number types")
        end
    else # ODE
        samp = Bridge.SamplePath(t,Vector{typeof(u0)}(undef,length(t)))
        W = nothing
        if typeof(alg) <: BridgeR3
            u = Bridge.solve!(Bridge.R3(),samp,u0,f)
        elseif typeof(alg) <: BridgeBS3
            u,_ = Bridge.solve!(Bridge.BS3(),samp,u0,f)
        end
    end

    DiffEqBase.build_solution(prob,  alg, u.tt, u.yy,
                   W = W,
                   timeseries_errors = timeseries_errors,
                   retcode = :Success)
end

ext: jl


text:
module Percolation

using Random,Distributions,Luxor

function set_seed(seed::Int64)
  Random.seed!(seed)
end

function poisson_clock(t0::Float64,t1::Float64,rate::Float64)::Vector{Float64}
  dt = Exponential(1/rate)
  t = t0 + rand(dt)
  ts = []
  while t < t1
    push!(ts,t)
    t += rand(dt)
  end
  ts
end

@enum EventType death infectl infectr

struct Event
  type::EventType
  site::Int64
  time::Float64
end

struct Perc
  events::Vector{Event}

  function Perc(es::Vector{Event})
    new(es)
  end

  function Perc(n::Int,T::Float64,λ::Float64)
    deaths = [[Event(death, i, clk) for clk in poisson_clock(0.,T,1.)] for i in 1:n]
    infectls = [[Event(infectl, i, clk) for clk in poisson_clock(0.,T,λ)] for i in 2:n]
    infectrs = [[Event(infectr, i, clk) for clk in poisson_clock(0.,T,λ)] for i in 1:(n-1)]
    events = [vcat(deaths...); vcat(infectls...); vcat(infectrs...)]
    sort!(events,by=e -> e.time)
    new(events)
  end
end

function simulate(init::BitVector, t1::Float64, perc::Perc)
  state = copy(init)
  for e in perc.events
    if e.time > t1
      break
    end
    if e.type == death
      state[e.site] = 0
    elseif e.type == infectl
      state[e.site-1] |= state[e.site]
    elseif e.type == infectr
      state[e.site+1] |= state[s.site]
    end
  end
  return state
end

function reverse_perc(perc::Perc, T::Float64)
  reverse_time_events = map(perc.events) do e
    if e.type == death
      Event(death,e.site,T - e.time)
    elseif e.type == infectl
      Event(infectr,e.site - 1,T - e.time)
    elseif e.type == infectr
      Event(infectl,e.site + 1,T - e.time)
    end
  end
  Perc(reverse(reverse_time_events))
end

const AliveInterval = Tuple{Float64, Union{Float64,Nothing}}

const AliveIntervals = Vector{AliveInterval}

struct PercResult
  ints :: Vector{AliveIntervals}
  linfects :: Vector{Vector{Float64}}
  rinfects :: Vector{Vector{Float64}}
  function PercResult(init::BitVector)
    n = length(init)
    ints = [AliveIntervals(init[i] ? [(0,nothing)] : []) for i in 1:n]
    linfects = [[] for i in 1:n]
    rinfects = [[] for i in 1:n]
    new(ints,linfects,rinfects)
  end
  function PercResult(ints::Vector{AliveIntervals},linfects::Vector{Vector{Float64}},rinfects::Vector{Vector{Float64}})
    new(ints,linfects,rinfects)
  end
end

alive(pr::PercResult, site::Int64) = length(pr.ints[site]) > 0 && pr.ints[site][end][2] == nothing

function infect(pr::PercResult,site::Int64, t::Float64)
  if !alive(pr,site)
    push!(pr.ints[site],(t,nothing))
  end
end

function kill(pr::PercResult,site::Int64,t::Float64)
  if alive(pr,site)
    pr.ints[site][end] = (pr.ints[site][end][1],t)
  end
end

function reverse_int(int::AliveInterval, T::Float64)
  if int[2] == nothing
    (0,T - int[1])
  else
    (T - int[2],T - int[1])
  end
end

function reverse_pr(pr::PercResult,T::Float64)
  PercResult(
    Vector{AliveIntervals}(map(is -> reverse(reverse_int.(is,T)), pr.ints)),
    Vector{Vector{Float64}}([[[]] ; map(ri -> reverse(T .- ri), pr.rinfects)[1:end-1]]),
    Vector{Vector{Float64}}([map(li -> reverse(T .- li), pr.linfects)[2:end] ; [[]]])
  )
end

function simulate_pr(init::BitVector, t1::Float64, perc::Perc)
  n = length(init)
  # makes the logic a bit easier if the dead ones died before we started
  pr = PercResult(init)
  for e in perc.events
    if e.time > t1
      break
    end
    if !alive(pr, e.site)
      continue
    end
    if e.type == death
      kill(pr,e.site,e.time)
    elseif e.type == infectl
      infect(pr,e.site-1,e.time)
      push!(pr.linfects[e.site],e.time)
    elseif e.type == infectr
      infect(pr,e.site+1,e.time)
      push!(pr.rinfects[e.site],e.time)
    end
  end
  pr
end

struct AnimParams
  pr::PercResult
  perc::Perc
  T::Float64
  width::Float64
  height::Float64
  frames::Int64
  function AnimParams(pr::PercResult,perc::Perc,T::Float64,w::Float64,h::Float64,fr::Int64)
    new(pr,perc,T,w,h,fr)
  end
  function AnimParams(init::BitVector, T::Float64, λ::Float64, width::Float64, height::Float64, frames::Int64)
    n = length(init)
    perc = Perc(n,T,λ)
    pr = simulate_pr(init, T, perc)
    new(pr,perc,T,width,height,frames)
  end
end

function reverse_pr(ap::AnimParams)
  n = length(ap.pr.ints)
  rev_perc = reverse_perc(ap.perc, ap.T)
  rev_init = BitVector(map(i -> alive(ap.pr,i), 1:n))
  rev_pr = reverse_pr(simulate_pr(rev_init, ap.T, rev_perc),ap.T)
  AnimParams(rev_pr,ap.perc,ap.T,ap.width,ap.height,ap.frames)
end

function draw_intervals(ap::AnimParams,cur_t::Float64)
  pr = ap.pr
  n = length(pr.ints)
  perc = ap.perc
  width = ap.width
  height = ap.height
  T = ap.T
  # height scale and width scale
  ws = width / n
  hs = height / T
  transform([0.9,0,0,-0.9,0,0])
  sethue("grey")
  pt(site,t) = Point((site - ((n+1) / 2)) * ws, (t - T/2) * hs)
  for i in 1:n
    line(pt(i,0),pt(i,T),:stroke)
    for e in perc.events
      if e.type == death
        circle(pt(e.site,e.time),4,:fill)
      elseif e.type == infectl
        arrow(pt(e.site,e.time),pt(e.site-1,e.time))
      elseif e.type == infectr
        arrow(pt(e.site,e.time),pt(e.site+1,e.time))
      end
    end
  end
  sethue("red")
  setline(4)
  for i in 1:n
    for int in pr.ints[i]
      if int[1] <= cur_t
        t2 = min(cur_t,int[2] == nothing ? T : int[2])
        line(pt(i,int[1]),pt(i,t2),:stroke)
        if t2 < min(cur_t, T)
          circle(pt(i,t2),6,:fill)
        end
      end
    end
    for t in pr.linfects[i]
      if t < cur_t
        arrow(pt(i,t),pt(i-1,t))
      end
    end
    for t in pr.rinfects[i]
      if t < cur_t
        arrow(pt(i,t),pt(i+1,t))
      end
    end
  end
end

function anim_contact_process(ap::AnimParams, path) where {T,S}
  t_scale = 1.5 * ap.T / ap.frames
  function scene(s,framenumber)
    draw_intervals(ap, t_scale * framenumber)
  end
  m = Movie(ap.width,ap.height,"Contact Process")
  animate(m,[Scene(m,scene,1:ap.frames)],creategif=true,pathname=path)
end

end


text:
# generates a random positive definite matrix.
function generaterandomposdefmat(n)
    A = rand(n,n)
    out = A'*A
    return forcesymmetric(out)
end

function generaterandomsym(n)
    return extractsym(rand(n,n))
end

function generaterandomskew(n)
    return extractskew(rand(n,n))
end

function generaterandomortho(n)
    A = randn(n,n)
    B = A*A'
    s,Q = LinearAlgebra.eigen(B)
    return Q
end

function generateStiefel(n,p)
    Q = generaterandomortho(p)
    out = zeros(Float64, n,p)
    out[1:p,1:end] = Q
    return out
end

function generaterandomhankel(n)
    a = rand(2*n-1)
    out = Matrix{Float64}(undef,n,n)
    for j = 1:n
        for i = 1:n
            out[i,j] = a[i+j-1]
        end
    end
    return out
end

# function generateisotopysubgroup(n,p)
#     O = generaterandomortho(n-p)
#     out = zeros(Float64, n,n)
#     out[1:p,1:p] = eye(p)
#     out[(n-p):end, (n-p):end] = O
#     return out
# end

# generate an invertible nxn indefinite matrix, with k positive eigenvalues.
function generateindefmatrix(n,k=rand(2:n-1))
    @assert 2 <= k <= n

    U,S,V = LinearAlgebra.svd(rand(n,n))
    for i = 1:k
        S[i] = rand() # rand() is between [0,1]
    end
    for i = k+1:length(S)
        S[i] = -rand()
    end
    return U*LinearAlgebra.diagm(0=>S)*U'
end

function generateposdefmatrix(n)
    U,S,V = LinearAlgebra.svd(rand(n,n))
    S = rand(n)
    out = U*LinearAlgebra.diagm(0=>S)*U'

    return forcesymmetric(out)
end

function generatenegdefmatrix(n)
    U,S,V = LinearAlgebra.svd(rand(n,n))
    S = -rand(n)
    out = U*LinearAlgebra.diagm(0=>S)*U'

    return forcesymmetric(out)
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
real_checked(A::AbstractArray) = A
function real_checked(A::AbstractArray{Complex{T}}) where T
    epsT = eps(real(eltype(A)))
    rtol = 1000epsT
    atol = 100epsT

    if any(abs(imag(x)) > rtol * abs(x) && abs(imag(x)) > atol for x in A)
        relerror = imag(A) ./ abs.(A)
        relerror[map(x -> abs(imag(x)) < atol, A)] .= 0
        @warn "Large imaginary part" norm(relerror)
    end

    real(A)
end

"""
Zero all elements of a Fourier space array which have no complex-conjugate partner
and may thus lead to an imaginary component in real space (after an iFFT).
"""
function clear_without_conjugate!(A::AbstractArray{T, 3}) where {T <: Complex}
    fft_size = size(A)
    # For even-length grids the element with largest negative integer coordinate
    # in G_vectors has no complex-conjugated partner. In order to avoid having
    # a spurious imaginary component in real-space, we set these elements to zero
    iseven(fft_size[1]) && (A[div(end, 2) + 1, :, :] .= 0)
    iseven(fft_size[2]) && (A[:, div(end, 2) + 1, :] .= 0)
    iseven(fft_size[3]) && (A[:, :, div(end, 2) + 1] .= 0)

    A
end


text:
"""
	Structure representing context-sensitive grammar
	Extends ExprRules.Grammar with constraints
"""
struct ContextGrammar <: GrammarType
	rules::Vector{Any}
	types::Vector{Symbol}
	isterminal::BitVector
	iseval::BitVector
	bytype::Dict{Symbol, Vector{Int}}
	childtypes::Vector{Vector{Symbol}}
	constraints::Vector{Constraint}
end

"""
    get_childtypes(rule::Any, types::AbstractVector{Symbol})
Returns the child types of a production rule.
"""
function get_childtypes(rule::Any, types::AbstractVector{Symbol})
    retval = Symbol[]
    if isa(rule, Expr)
        for arg in rule.args
            append!(retval, get_childtypes(arg, types))
        end
    elseif rule ∈ types
        push!(retval, rule)
    end
    return retval
end

"""
	Reduces the set of possible children of a node using the grammar's constraints
"""
function propagate_contraints(grammar::ContextGrammar, context::GrammarContext, child_rules::Vector{Int})
	domain = child_rules

	for propagator in grammar.constraints
		domain = propagate(propagator, context, domain)
	end

	return domain
end

"""
	reimplementation of the ExprRules._next_state!
	Change: child expressions are filtered so that the constraints are not violated

"""
function _next_state!(node::RuleNode, grammar::ContextGrammar, max_depth::Int, context::GrammarContext)

	if max_depth < 1
	    return (node, false) # did not work
	elseif isterminal(grammar, node)
	    # do nothing
	    if iseval(grammar, node.ind) && (node._val === nothing)  # evaluate the rule
		node._val = eval(grammar.rules[node.ind].args[2])
	    end
	    return (node, false) # cannot change leaves
	else # !isterminal
	    # if node is not terminal and doesn't have children, expand every child
	    if isempty(node.children)
		if max_depth ≤ 1
		    return (node,false) # cannot expand
		end

		child_index = 1  # keepa track of which child we are processing now (needed for context)

		# build out the node
		for c in child_types(grammar, node)
		    worked = false
		    i = 0
		    child = RuleNode(0)

		    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
		    push!(new_context.nodeLocation, child_index)

		    child_rules = [x for x in grammar[c]]  # select all applicable rules
		    child_rules = propagate_contraints(grammar, new_context, child_rules)  # filter out those that violate constraints

		    while !worked && i < length(child_rules)
			i += 1
			child = RuleNode(child_rules[i])

			if iseval(grammar, child.ind) # if rule needs to be evaluated (_())
			    child._val = eval(grammar.rules[child.ind].args[2])
			end

			worked = true
			if !isterminal(grammar, child)
			    child, worked = _next_state!(child, grammar, max_depth-1, new_context)
			end
		    end
		    if !worked
			return (node, false) # did not work
		    end
		    push!(node.children, child)

		    child_index += 1
		end

		return (node, true)
	    else # not empty
		# make one change, starting with rightmost child
		worked = false
		child_index = length(node.children) + 1
		while !worked && child_index > 1
		    child_index -= 1
		    child = node.children[child_index]

		    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
		    push!(new_context.nodeLocation, child_index)

		    child, child_worked = _next_state!(child, grammar, max_depth-1, new_context)
		    while !child_worked
			child_type = return_type(grammar, child)

			child_rules = [x for x in grammar[child_type]]  # get all applicable rules
			child_rules = propagate_contraints(grammar, new_context, child_rules)  # filter ones that violate constraints

			i = something(findfirst(isequal(child.ind), child_rules), 0)
			if i < length(child_rules)
			    child_worked = true
			    child = RuleNode(child_rules[i+1])

			    # node needs to be evaluated
			    if iseval(grammar, child.ind)
				child._val = eval(grammar.rules[child.ind].args[2])
			    end

			    if !isterminal(grammar, child)
				child, child_worked = _next_state!(child, grammar, max_depth-1, new_context)
			    end
			    node.children[child_index] = child
			else
			    break
			end
		    end

		    if child_worked
			worked = true

			# reset remaining children
			for child_index2 in child_index+1 : length(node.children)
			    c = child_types(grammar, node)[child_index2]
			    worked = false
			    i = 0
			    child = RuleNode(0)

			    new_context = GrammarContext(context.originalExpr, deepcopy(context.nodeLocation))
			    push!(new_context.nodeLocation, child_index2)

			    child_rules = [x for x in grammar[c]]  # take all applicable rules
			    child_rules = propagate_contraints(grammar, new_context, child_rules)  # reomove ones that violate constraints

			    while !worked && i < length(child_rules)
				i += 1
				child = RuleNode(child_rules[i])

				if iseval(grammar, child.ind)
				    child._val = eval(grammar.rules[child.ind].args[2])
				end

				worked = true
				if !isterminal(grammar, child)
				    child, worked = _next_state!(child, grammar, max_depth-1, new_context)
				end
			    end
			    if !worked
				break
			    end
			    node.children[child_index2] = child
			end
		    end
		end

		return (node, worked)
	    end
	end
    end

mutable struct ContextExpressionIterator <: ExprIter
	grammar::ContextGrammar
	max_depth::Int
	sym::Symbol
end

function Base.iterate(iter::ContextExpressionIterator)
	init_node = RuleNode(0)  # needed for propagating constraints on the root node
	init_context = GrammarContext(init_node)

	grammar, sym, max_depth = iter.grammar, iter.sym, iter.max_depth

	# propagate constraints on the root node
	sym_rules = [x for x in grammar[sym]]
	sym_rules = propagate_contraints(grammar, init_context, sym_rules)
	#node = RuleNode(grammar[sym][1])
	node = RuleNode(sym_rules[1])

	if isterminal(grammar, node)
	    return (deepcopy(node), node)
	else
	    context = GrammarContext(node)
	    node, worked =  _next_state!(node, grammar, max_depth, context)
	    while !worked
		# increment root's rule
		rules = [x for x in grammar[sym]]
		rules = propagate_contraints(grammar, init_context, rules) # propagate constraints on the root node

		i = something(findfirst(isequal(node.ind), rules), 0)
		if i < length(rules)
		    node, worked = RuleNode(rules[i+1]), true
		    if !isterminal(grammar, node)
			node, worked = _next_state!(node, grammar, max_depth, context)
		    end
		else
		    break
		end
	    end
	    return worked ? (deepcopy(node), node) : nothing
	end
end

function Base.iterate(iter::ContextExpressionIterator, state::RuleNode)
	grammar, max_depth = iter.grammar, iter.max_depth
	context = GrammarContext(state)
	node, worked = _next_state!(state, grammar, max_depth, context)

	while !worked
	    # increment root's rule
	    init_node = RuleNode(0)  # needed for propagating constraints on the root node
    	    init_context = GrammarContext(init_node)

	    rules = [x for x in grammar[iter.sym]]
	    rules = propagate_contraints(grammar, init_context, rules)

	    i = something(findfirst(isequal(node.ind), rules), 0)
	    if i < length(rules)
		node, worked = RuleNode(rules[i+1]), true
		if !isterminal(grammar, node)
		    context = GrammarContext(node)
		    node, worked = _next_state!(node, grammar, max_depth, context)
		end
	    else
		break
	    end
	end
	return worked ? (deepcopy(node), node) : nothing
end

macro cgrammar(ex)
	rules = Any[]
	types = Symbol[]
	bytype = Dict{Symbol,Vector{Int}}()
	for e in ex.args
	    if isa(e, Expr)
		if e.head == :(=)
		    s = e.args[1] # name of return type
		    rule = e.args[2] # expression?
		    rvec = Any[]
		    _parse_rule!(rvec, rule)
		    for r in rvec
			push!(rules, r)
			push!(types, s)
			bytype[s] = push!(get(bytype, s, Int[]), length(rules))
		    end
		end
	    end
	end
	alltypes = collect(keys(bytype))
	is_terminal = [isterminal(rule, alltypes) for rule in rules]
	is_eval = [iseval(rule) for rule in rules]
	childtypes = [get_childtypes(rule, alltypes) for rule in rules]
	return ContextGrammar(rules, types, is_terminal, is_eval, bytype, childtypes, [])
end

_parse_rule!(v::Vector{Any}, r) = push!(v, r)

function _parse_rule!(v::Vector{Any}, ex::Expr)
	if ex.head == :call && ex.args[1] == :|
	     terms = length(ex.args) == 2 ?
		collect(interpret(ex.args[2])) :    #|(a:c) case
		ex.args[2:end]                      #a|b|c case
	    for t in terms
		_parse_rule!(v, t)
	    end
	else
	    push!(v, ex)
	end
end

"""
    Add constraint to the grammar
"""
addconstraint!(grammar::ContextGrammar, cons::Constraint) = push!(grammar.constraints, cons)


text:
# Test classified collections

using DataStructures
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

# classified lists

c = classified_lists(Compat.ASCIIString, Int)

push!(c, "low", 1)
push!(c, "low", 2)
push!(c, "low", 3)
push!(c, "high", 4)
push!(c, "high", 5)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

@test c["low"] == [1, 2, 3]
@test c["high"] == [4, 5]

@test length(c) == 2
@test collect(keys(c)) == Compat.ASCIIString["high","low"]

pop!(c,"low")
@test !haskey(c,"low")

# classified sets

c = classified_sets(Compat.ASCIIString, Int)

push!(c, "low", 1)
push!(c, "low", 2)
push!(c, "low", 3)
push!(c, "low", 1)
push!(c, "low", 2)

push!(c, "high", 4)
push!(c, "high", 5)
push!(c, "high", 5)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

@test isa(c["low"], Set{Int})
@test isa(c["high"], Set{Int})

@test sort(collect(c["low"])) == [1, 2, 3]
@test sort(collect(c["high"])) == [4, 5]

# classified counters

c = classified_counters(Compat.ASCIIString, Float64)

push!(c, "low", 1.)
push!(c, "low", 2.)
push!(c, "low", 3.)
push!(c, "low", 1.)
push!(c, "low", 2.)
push!(c, "low", 2.)

push!(c, "high", 4.)
push!(c, "high", 5.)
push!(c, "high", 5.)

@test haskey(c, "low")
@test haskey(c, "high")
@test !haskey(c, "mid")

cl = c["low"]
ch = c["high"]

@test isa(cl, Accumulator{Float64, Int})
@test isa(ch, Accumulator{Float64, Int})

@test cl[1.] == 2
@test cl[2.] == 3
@test cl[3.] == 1
@test ch[4.] == 1
@test ch[5.] == 2


text:
#=
diffl.jl
Left finite differences "in-place" (pre-allocated outpus)

Could use StaticKernels.jl for this; see timing test in ../../time.

Inspired by:
https://docs.julialang.org/en/latest/manual/performance-tips/#Pre-allocating-outputs-1

2019-06-22 Jeff Fessler, University of Michigan
=#

export diffl, diffl!, diffl_adj, diffl_adj!, diffl_map

using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO

"""
    diffl!(g::AbstractArray, x::AbstractArray, dim::Int ; ...)

Apply left finite difference operator to input array `x`,
storing the result "in-place" in pre-allocated output array `g`.

Arrays `g` and `x` must have the same size, and cannot alias.
The "first" elements of `g` are zero for dimension `dim`.
The default is `dim=1`.

Option:
- `add::Bool = false` use `x[i] + x[i-1]` instead of `x[i] - x[i-1]`
- `edge::Symbol = :zero` set the first elements of dimension `dim` to 0

Choose `edge=:circ` to use circulant (aka periodic) boundary conditions.
Choose `edge=:none` to leave the first elements untouched.

In 1D, if `x = [2, 6, 7]` then `g = [0, 4, 1]` with default options.
"""
function diffl!(
    g::AbstractArray{Tg,N},
    x::AbstractArray{Tx,N},
    dim::Int ;
    edge::Symbol=:zero,
    add::Bool=false,
) where {Tg,Tx,N}

    Base.require_one_based_indexing(g) && Base.require_one_based_indexing(x)
    1 <= dim <= N || throw(ArgumentError("dimension $dim out of range (1:$N)"))
    size(g) != size(x) && throw(DimensionMismatch("sizes g=>$(size(g)) vs x=>$(size(x))"))

    Nd = size(x, dim)
    if add
        @inbounds selectdim(g, dim, 2:Nd) .=
        selectdim(x, dim, 2:Nd) .+ selectdim(x, dim, 1:(Nd-1))
    else
        @inbounds selectdim(g, dim, 2:Nd) .=
        selectdim(x, dim, 2:Nd) .- selectdim(x, dim, 1:(Nd-1))
    end

    # handle edge conditions
    g1 = selectdim(g, dim, 1)
    if edge === :zero
        @inbounds g1 .= zero(Tg)
    elseif edge === :circ
        if add
            @inbounds g1 .= selectdim(x, dim, 1) .+ selectdim(x, dim, Nd)
        else
            @inbounds g1 .= selectdim(x, dim, 1) .- selectdim(x, dim, Nd)
        end
    else
        edge != :none && throw(ArgumentError("edge $edge"))
        # caution: in this case g1 is untouched, possibly undef
    end

    return g
end

# for default dim=1 case
diffl!(g::AbstractArray, x::AbstractArray ; kwargs...) = diffl!(g, x, 1 ; kwargs...)

"""
    diffl!(g::AbstractArray, x::AbstractArray, dims::AbstractVector{Int} ; ...)

When `x` is a `N`-dimensional array, the `i`th slice of the `g` array
(along its last dimension) is the `diffl!` of `x` along `dims[i]`.
This is useful for total variation (TV) and other regularizers
that need finite differences along multiple dimensions.
"""
function diffl!(
    g::AbstractArray{Tg,Ng},
    x::AbstractArray{Tx,Nx},
    dims::AbstractVector{Int} ;
    kwargs...,
) where {Tg,Tx,Ng,Nx}

    Ng != Nx+1 && throw(DimensionMismatch("Ng=$Ng Nx=$Nx"))
    size(g) != (size(x)..., length(dims)) &&
        throw(DimensionMismatch("sizes g=>$(size(g)) vs x=>$(size(x))"))

    for (i,d) in enumerate(dims)
        diffl!(selectdim(g, Ng, i), x, d ; kwargs...)
    end
    return g
end

"""
    g = diffl(x::AbstractArray ; ...)
Allocating version of `diffl!` along `dim=1`
"""
diffl(x::AbstractArray ; kwargs...) = diffl(x, 1 ; kwargs...)

"""
    g = diffl(x::AbstractArray, dim::Int ; ...)
Allocating version of `diffl!` along `dim`
"""
diffl(x::AbstractArray, dim::Int ; kwargs...) = diffl!(similar(x), x, dim ; kwargs...)

"""
    g = diffl(x::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of `diffl!` for `dims`
"""
diffl(x::AbstractArray, dims::AbstractVector{Int} ; kwargs...) =
    diffl!(similar(x, size(x)..., length(dims)), x, dims ; kwargs...)

"""
    diffl_adj!(z, g, dim::Int ; ...)

Adjoint of left finite difference `diffl!`, in-place.
Arrays `z` and `g` must be same size.
See `diffl` for details.
"""
function diffl_adj!(
    z::AbstractArray{Tz,N},
    g::AbstractArray{Tg,N},
    dim::Int ;
    reset0::Bool=true,
    edge::Symbol=:zero,
    add::Bool=false,
) where {Tz,Tg,N}

    1 <= dim <= N || throw(ArgumentError("dimension $dim out of range (1:$N)"))
    size(z) != size(g) && throw(DimensionMismatch("sizes z=>$(size(z)) vs g=>$(size(g))"))

    # todo: handle reset0 better
    if reset0
        z .= zero(Tz)
    end

    Nd = size(g, dim)

    if edge === :zero
        @inbounds selectdim(z, dim, 2:Nd) .+= selectdim(g, dim, 2:Nd)
    elseif edge === :circ
    #    selectdim(z, dim, 1:Nd) .+= selectdim(g, dim, 1:Nd)
        @inbounds z .+= g
    else
        edge != :none && throw(ArgumentError("edge $edge"))
        # in this case g1 is unused, even if undef
    end

    if add
        @inbounds selectdim(z, dim, 1:(Nd-1)) .+= selectdim(g, dim, 2:Nd)
        if edge === :circ
            @inbounds selectdim(z, dim, Nd) .+= selectdim(g, dim, 1)
        end
    else
        @inbounds selectdim(z, dim, 1:(Nd-1)) .-= selectdim(g, dim, 2:Nd)
        if edge === :circ
            @inbounds selectdim(z, dim, Nd) .-= selectdim(g, dim, 1)
        end
    end

    return z
end

"""
    diffl_adj!(z::AbstractArray, g::AbstractArray, dims::AbstractVector{Int} ; ...)

Adjoint of `diffl!` for multiple dimensions `dims`.
Here `g` must have one more dimension than `z`.
"""
function diffl_adj!(
    z::AbstractArray{Tz,Nz},
    g::AbstractArray{Tg,Ng},
    dims::AbstractVector{Int} ;
    kwargs...,
) where {Tz,Tg,Nz,Ng}

    Ng != Nz+1 && throw(DimensionMismatch("Ng=$Ng Nz=$Nz"))
    size(g) != (size(z)..., length(dims)) &&
        throw(DimensionMismatch("sizes g=>$(size(g)) vs z=>$(size(z))"))

    for (i,d) in enumerate(dims)
        diffl_adj!(z, selectdim(g, Ng, i), d ; reset0 = (i==1), kwargs...)
    end
    return z
end

"""
    z = diffl_adj(g::AbstractArray ; ...)
Allocating version of `diffl!` along `dim=1`
"""
diffl_adj(g::AbstractArray ; kwargs...) = diffl_adj(g, 1 ; kwargs...)

"""
    z = diffl(g::AbstractArray, dim::Int ; ...)
Allocating version of `diffl!` along `dim`
"""
diffl_adj(g::AbstractArray, dim::Int ; kwargs...) =
    diffl_adj!(similar(g), g, dim ; reset0=true, kwargs...)

"""
    z = diffl_adj(g::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of `diffl!` for `dims`
"""
function diffl_adj(
    g::AbstractArray{T,N},
    dims::AbstractVector{Int} ;
    kwargs...,
) where {T,N}

    size(g)[end] != length(dims) &&
        throw(DimensionMismatch("sizes g=>$(size(g)) vs dims=$dims"))
    return diffl_adj!(similar(g, size(g)[1:(N-1)]...), g, dims ; kwargs...)
end

"""
    T = diffl_map(N::Dims{D}, dims::AbstractVector{Int} ; kwargs...)
    T = diffl_map(N::Dims{D}, dim::Int ; kwargs...)

in
- `N::Dims` image size

options: see `diffl!`
- `T::Type` for `LinearMapAA`, default `Float32`
- `operator::Bool = true` use `false` for `LinearMapAM`

out
- `T` `LinearMapAA` object for computing finite differences via `T*x`
using `diffl!` and `diffl_adj!`
"""
function diffl_map(
    N::Dims{D},
    dims::AbstractVector{Int} ;
    T::Type=Float32,
    edge::Symbol = :zero,
    operator::Bool = true, # !
    kwargs...,
) where {D}

    !all(1 .<= dims .<= D) && throw(ArgumentError("dims $dims"))
    (edge == :none) && throw("edge=$edge unsupported")

    forw!(g,x) = diffl!(g, x, dims ; edge=edge, kwargs...)
    back!(z,g) = diffl_adj!(z, g, dims ; edge=edge, kwargs...)

    if operator
        return LinearMapAA(forw!, back!,
            (length(dims), 1) .* prod(N) ;
            prop = (name="diffl_map", N=N, dims=dims),
            T=T, operator=true,
            idim=N, odim=(N...,length(dims)),
        )
    end

    gshape = g -> reshape(g, N..., length(dims))
    return LinearMapAA(
        (g,x) -> forw!(gshape(g), reshape(x, N)),
        (z,g) -> vec(back!(reshape(z,N), gshape(g))),
        (length(dims), 1) .* prod(N),
        (name="diffl_map", N=N, dims=dims),
        T=T,
    )
end

# for single dimension case
function diffl_map(
    N::Dims{D},
    dim::Int ;
    T::Type = Float32,
    edge::Symbol = :zero,
    operator::Bool = true, # !
    kwargs...,
) where {D}

    (1 .<= dim .<= D) || throw(ArgumentError("dim $dim"))
    (edge == :none) && throw("edge=$edge unsupported")

    forw!(g,x) = diffl!(g, x, dim ; edge=edge, kwargs...)
    back!(z,g) = diffl_adj!(z, g, dim ; edge=edge, kwargs...)

    if operator
        return LinearMapAA(forw!, back!,
            (1, 1) .* prod(N) ;
            prop = (name="diffl_map", N=N, dim=dim),
            T=T, operator=true,
            idim=N, odim=N,
        )
    end

    return LinearMapAA(
        (g,x) -> forw!(reshape(g, N), reshape(x, N)),
        (z,g) -> vec(back!(reshape(z,N), reshape(g, N))),
        (1, 1) .* prod(N) ;
        prop = (name="diffl_map", N=N, dim=dim),
        T=T,
    )
end

# default dim=1
diffl_map(N::Dims ; kwargs...) = diffl_map(N, 1, ; kwargs...)


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module ReverseAD

import ForwardDiff
import MathOptInterface
import ..Nonlinear
import SparseArrays

const MOI = MathOptInterface

# Override basic math functions to return NaN instead of throwing errors.
# This is what NLP solvers expect, and sometimes the results aren't needed
# anyway, because the code may compute derivatives wrt constants.
import NaNMath:
    sin,
    cos,
    tan,
    asin,
    acos,
    acosh,
    atanh,
    log,
    log2,
    log10,
    lgamma,
    log1p,
    pow,
    sqrt

include("Coloring/Coloring.jl")
include("graph_tools.jl")
include("types.jl")
include("utils.jl")

include("reverse_mode.jl")
include("forward_over_reverse.jl")
include("mathoptinterface_api.jl")

include("precompile.jl")

end  # module

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# A simple style processor
#using PrettyPrint  # for debug
using OrderedCollections
# To find the position and size of an element
#
# 1. get its x, y, w, h properties.
# 2. check element `position` property in the style context to see if it is `relative`.
#    If it is, create relative lengths, otherwise, get the unit from the global
#    setting and create absolute length.
# 3. check its alignment in the style context, if there is alignment property,
#    create an attribute.
# 4. will do margin and padding later...
abstract type VGElement end

PVDict= OrderedDict{String, String}  # Dict of (property=>value)
SPVDict = OrderedDict{String, OrderedDict{String, String}} # selector => PVDict
CssRules = OrderedDict{String, String}  # selector => one CSS rule in string

#TODO: complete this dictionary
const default_pvd = OrderedDict([
    ("fill", "none"),
    ("stroke", "black"),
    ("stroke-width", "1"),
    ("--vg-unit", "px"),
    ("--vg-showlabel", "true"),
    ("--vg-suppress", "false"),
    ("--vg-shape", "bracket"),
    ("--vg-offset", "0.2"),
    ("--vg-distance", "5"),
    ("--vg-side", "above"),
    ("--vg-brace-height", "5"),
    ("--vg-brace-stroke", "blue"),
    ("--vg-brace-stroke-width", "1"),
    ("--vg-label-side", "above"),
    ("--vg-label-anchor", "middle"),
    ("--vg-label-distance", "8"),
    ("--vg-label-font-size", "18"),
    ("--vg-inner-center-vshift", "-2")            # hack
   ])

# minimum set of properties needed for drawing rect, circle, polyline, etc
const supp_properties = OrderedDict([
    ("fill", default_pvd["fill"]),
    ("stroke", default_pvd["stroke"]),
    ("stroke-width", default_pvd["stroke-width"])
    ])

# we strip leading and trailing whitespaces for a property name and values, but does not
# consolidate inner whitespaces...
function cssrules_to_spvdict(css::CssRules)
    styles = SPVDict()

    for k in keys(css)
        rule_list = split(strip(css[k]), ';')
        # NOTE: empty string could appear in rule_list, e.g., the splitting
        # result of "color:red;"
        pvd = OrderedDict{String, String}()
        for r in rule_list
            isempty(r) && continue
            keyval_pair = split(r, ':')
            length(keyval_pair) == 1 && error("Invalid CSS format")
            #if property keyval_pair[1] already exists in rdict, the last one will be the
            #finally accepted one
            pvd[strip(keyval_pair[1])] = strip(keyval_pair[2])
            # keyval_pair[2] is a string containing space-separated values
            # for example "solid black" (for a property such as `border`)
            # that can be further splitted by `split(str, r" +")`
        end
        # FIXME: what if k already exits in styles? we should merge instead of
        # replace
        styles[k] = pvd
    end
    return styles
end

# TODO: support combinator selectors
function calculate_specificity(e::VGElement, selector::String)
    sels = split(strip(selector), r"\s+")
    if length(sels) > 1
        error("We can only process simple selectors so far:), ", @__FILE__, " ", @__LINE__)
    end

    # we only compare simple selectors with `sels[1]` for now
    #(string(typeof(e)) == string(@__MODULE__)*"."*sels[1]) && return 1
    (string(typeof(e)) == sels[1]) && return 1
    e.class != "" && "."*e.class == sels[1] && return 10
    e.id != "" && "#"*e.id == sels[1] && return 100
    return 0
end

# Not too ambitious at this time. Just extract all **custom** CSS rules for `el`
# and collect them in a property-value dictionary
function get_style_PVD(e::VGElement, spv_dict::SPVDict)
    fin_dict = PVDict()
    specificity = OrderedDict{String, Int}()
    for selector in keys(spv_dict)
        s = calculate_specificity(e, selector)
        s == 0 && continue
        pv_dict = spv_dict[selector]
        for p in keys(pv_dict)
            if !(p in keys(fin_dict)) || specificity[p] <= s
                fin_dict[p] = pv_dict[p]
                specificity[p] = s
            end
        end
    end
    return fin_dict
end

function get_property(pvd::PVDict, p::String, v::String="")
    global default_pvd

    haskey(pvd, p) && return pvd[p]
    haskey(default_pvd, p) && return default_pvd[p]
    !isempty(v) && return v
    return ""
end

function isnumber(s::String)
    dots, others= 0, 0, 0
    for c in s
        if isdigit(c); continue; end
        if c == '.'; dots += 1; continue; end
        others += 1
        return false
    end
    return dots == 1
end

# return a numerical value, with unit converted to px
# + a pure number input is converted to a Float64 and returned
# + a number with unit `px` is stripped off of the `px`, converted to a Float64, and
#   returned
# TODO: actually units other than `px` are not supported for now.
function get_numeric_px(val::String)
    val = strip(val)
    length(val) <= 2 && error("invalid numerical format")
    # a pure number is converted to a Float64 and returned
    isnumber(val) && return parse(Float64, val)
    unit = val[end-1:end]
    num = val[1:end-2]
    !isnumber(num) && error("invalid numerical format")
    unit == "px" && return parse(Float64, num)
    error("length units other than `px` are not supported!")
end

const tag_of = Dict([
    (:Box, "rect"),
    (:Circ, "circle"),
    (:Polyline, "polyline")
    ])

# convert a single rule dictionary of (propert=>value)s into a CSS rule string
function stringify(pvd::PVDict)
    css = ""
    for property in keys(pvd)
        css *= property*':'*pvd[property]*';'
    end
    return css
end

# Remove all **custom** properties and convert one PVDict `pvd` into a CSS rule string.
function nstringify(pvd::PVDict)
    css = ""
    for property in keys(pvd)
        (length(property) > 5 && property[1:5] == "--vg-") && continue
        css *= property*':'*pvd[property]*';'
    end
    return css
end

# Remove all custom properties and convert a set of PVDicts into a CSS string
function nstringify_css(spvd::SPVDict)
    css = ""
    for selector in keys(spvd)
        css *= selector*'{'*nstringify(spvd[selector])*'}'
    end
    return css
end


text:
_hasbg(pl::Plot) = pl.theme.background_color != nothing
_addbg!(pl::Plot) = pl.theme.background_color=colorant"white"
_removebg!(pl::Plot) = pl.theme.background_color=nothing

# works for both plots and subplots arranged by hstack/vstack
function _saveplot(pl,format::Symbol,outPath::AbstractString; width=29.7cm, height=21cm, kwargs...)
    filename = string(outPath,'.', format==:jssvg ? "js.svg" : format)
    if format==:png
        draw(PNG(filename, width, height; kwargs...), pl)
    elseif format==:pdf
        draw(PDF(filename, width, height; kwargs...), pl)
    elseif format==:svg
        draw(SVG(filename, width, height; kwargs...), pl)
    elseif format==:jssvg
        draw(SVGJS(filename, width, height; kwargs...), pl)
    else
        error("Unknown image format $format")
    end
    filename
end

function saveplot(pl::Plot,format::Symbol,outPath::AbstractString; solidBG=format==:png, kwargs...)
    changeBG = solidBG && !_hasbg(pl)
    changeBG && _addbg!(pl)
    filename = _saveplot(pl,format,outPath; kwargs...)
    changeBG && _removebg!(pl)
    filename
end

function saveplot(pl::AbstractArray{Plot},format::Symbol,outPath::AbstractString; solidBG=format==:png, kwargs...)
    changeBG = .!map(_hasbg, pl) .& solidBG
    for (p,c) in zip(pl,changeBG)
        c && _addbg!(p)
    end

    combined = hstack(mapslices(vstack,pl,2)...) # stack vertically and then horizontally
    filename = _saveplot(combined,format,outPath; kwargs...)

    for (p,c) in zip(pl,changeBG)
        c && _removebg!(p)
    end

    filename
end

saveplot(pl,formats::AbstractArray{Symbol},outPath::AbstractString; kwargs...) = map(fmt->saveplot(pl,fmt,outPath;kwargs...), formats)


text:
using FEniCS
using Test

FEniCS.set_log_level(FEniCS.WARNING)

examples_dir = joinpath(@__DIR__, "..", "examples")
@assert ispath(examples_dir)
example_filenames = readdir(examples_dir)
@assert "tutorial1.jl" in example_filenames

@testset "Example $(filename)" for filename in example_filenames
    path = joinpath(examples_dir, filename)
    @test (include(path) ;true)
    # Test that PyCall is not used in the example.
    # Use of PyCall indicates, that we did not wrap enough
    # funcionality
    s = read(path, String)
    @test !occursin("PyCall", s)
end

@testset "Creation" begin
   @test include("test_create.jl")
   @test include("test_pycreate.jl")
   @test include("test_jfem.jl")
end;

@testset "assign" begin
    mesh = UnitIntervalMesh(10)
    V = FunctionSpace(mesh, "P", 1)
    u = interpolate(Expression("x[0]", degree=1), V)
    arr1 = get_array(u)
    arr2 = randn(size(arr1))
    assign(u, arr2)
    @test get_array(u) == arr2
    assign(u, arr1)
    @test get_array(u) == arr1
end

#tests relating to interface.jl file
@testset "Interface" begin
   global mesh = UnitSquareMesh(2,2)
   dbc(x,y) = 1
   global u = feMesh(mesh,dbc)
   @test (u.n_nodes==9)
   @test (u.n_elements==8)
   @test (u.n_internal_nodes==1)
   @test (u.n_boundary_nodes==8)
   @test (u.node_vals == [1,1,1,1,0,1,1,1,1])
end;


text:
struct Fallback end

struct PatchStore{K,V}
    dct::Dict{K,V}
    prev::Dict{K,V}
end

"Data type for object returned by `Expr.splitdef` function."
const FunctionDef = Dict{Symbol, Any}


text:
module MyProject

using Expectations, Distributions

function foo(μ = 1., σ = 2.)
    println("Modified foo definition again")
    d = Normal(μ, σ)
    E = expectation(d)
    return E(x -> sin(x))
end

function bar()
    return 10
end

export foo

end # module


text:
@time include("src/required.jl")

outcome = "bart"

@load joinpath("data", "jld2", "$(outcome)-distributions.jld2") distributions

@load joinpath("data", "jld2", "$(outcome)-distributions_xtras.jld2") prob_distrib lower_distrib upper_distrib

prob_distrib[1].grid
filter(!isnothing, prob_distrib[1].grid)

plot_layer(prob_distrib[1])
plot_layer(prob_distrib[1]; c=:BuPu)

# Create matrix Y
Y = Ymatrix(prob_distrib; transform=false)

# Create matrix of observed sites only
Yobs = _Yobs(Y)

# Sort Yobs by rows & columns sums
Ysort = sortslices(Yobs; dims=1, by=sum) |> y -> sortslices(y; dims=2, by=sum, rev=true)
# Heatmap of Yobs
Yobs_plot = heatmap(
    Ysort;
    # title = "$(titlecase(outcome)) matrix Y (sorted by row and column sum)",
    ylabel="Sites",
    yticks=:none,
    xlabel="Species number",
)

## Richness

# Get richness
richness = calculate_richness(Y, prob_distrib[1])
richness.grid

# Plot richness
richness_plot = plot_layer(
    richness;
    c=:viridis,
    # title = "Richness ($outcome distributions)",
    clim=(0.0, maximum(richness)),
    colorbar_title="Number of species per site",
)

## LCBD

# Get LCBD values
lcbd = calculate_lcbd(Y, prob_distrib[1]; transform=true, relative=true)

# Plot relative values
lcbdtr_plot = plot_layer(
    lcbd;
    c=:viridis,
    # title = "LCBD values per site ($(outcome) distributions, hellinger transformed)",
    colorbar_title="Relative LCBD value",
)

## Relationship

# Plot relationship as histogram2d
rel2d_plot = histogram2d(
    richness,
    lcbd;
    c=:viridis,
    bins=40, # title = "Relationship",
    xlabel="Richness",
    ylabel="LCBD",
    colorbar_title="Number of sites",
    xlim=(0.0, Inf),
    ylim=(0.0, 1.0),
)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
# This file was generated, do not modify it. # hide
train, test = partition(collect(eachindex(y)), 0.7, shuffle=true, rng=StableRNG(612));

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
Dynamically import the top-level module given by `sym`, and return it if
possible.
"""
function dynamic_import_toplevel_module(sym)::Nullable{Module}
    info("dynamic import: $sym")
    try
        eval(Main, :(import $sym))
        Nullable(getfield(Main, sym))
    catch
        Nullable()
    end
end


text:
include("test_exports.jl")
include("test_intersection.jl")
include("test_vector_cache.jl")
include("test_util.jl")
include("test_box_types.jl")
include("test_extensions.jl")


text:
function PrintArray(_array)
    output = ""
    for i in _array
        output = output .. i .. ","
    end
    print(output)
end

# Linear Array Search is O(n)
# Non-recursive
function Array_Search_Linear(_array, _needle)
    for i in _array
        if _array[i] == _needle then
            return i - 1
        end
    end
    return -1
end

# Binary Array Search is O(n), if unsorted but while sorted it is O(Log n).
# Recursive
function Array_Search_Binary(_array, _needle, _index, _endex)
    if _index == nil
        _index = 1
    end
    if _endex == nil
        _endex = #_array
    end

    if _array[_index] == _needle then
        return _index - 1
    end
    if _array[_endex] == _needle then
        return _endex - 1
    end

    if _index >= _endex then
        return -1
    end

    mid = math.floor((_index + (_endex - 1)) / 2)

    if _array[mid] == _needle then
        return mid - 1
    end

    if _array[mid] > _needle then
        return Array_Search_Binary(_array, _needle, _index, mid - 1)
    end

    return Array_Search_Binary(_array, _needle, mid + 1, _endex)
end

function Array_Search_Jump(_array, _needle)
    arraySize = #_array
    step = math.floor(math.sqrt(arraySize))

    prev = 1

    while(_array[math.min(step, arraySize) - 1] < _needle)
        prev = step
        step = step + math.floor(math.sqrt(arraySize))
        if prev >= arraySize then
            return -1
        end
    end

    while(_array[prev] < _needle)
        prev = prev + 1
        if prev == math.min(step, arraySize) then
            return -1
        end
    end

    if (_array[prev] == _needle) then
        return prev
    end
    return -1
end

function Array_Interpolation_Search(_array, _needle)
    low = 1
    high = #_array - 1

    while low <= high and _needle >= _array[low] and _needle <= _array[high]
        if low == high then
            if _array[low] == _needle then
                return low
            end
            return -1
        end

        pos = low + math.floor(((high - low + 0.0) / (_array[high] - _array[low])) * (_needle - _array[low]))
        if _array[pos] == _needle then
            return pos
        end
        if _array[pos] < _needle then
            low = pos + 1
        else
            high = pos - 1
        end
    end

    return -1
end

function Array_Exponential_Search(_array, _needle)
    if _array[1] == _needle then
        return 1
    end

    i = 1
    while (i < #_array and _array[i] <= _needle)
        i = i * 2
    end

    return Array_Search_Binary(_array, _needle, math.floor(i / 2), math.min(i, #_array - 1))
end

function Array_selectionSort(_array)
    min_idx = 1
    for i = 1, #_array - 1
        min_idx = i
        for j = i+1, #_array
            if _array[j] < _array[min_idx] then
                min_idx = j
            end
        end

        _array[min_idx], _array[i] = _array[i], _array[min_idx]
    end

    PrintArray(_array)
    return -1 # Ignore this, just to hide false-positive errors.
end

function Array_bubbleSort(_array)
    for i = 1, #_array
        for j = 1, #_array - i - 1
            if _array[j] > _array[j + 1] then
                _array[j], _array[j + 1] = _array[j + 1], _array[j]
            end
        end
    end

    PrintArray(_array)
    return -1 # Ignore this, just to hide false-positive errors.
end

function Array_insertionSort(_array)
    for i = 1, #_array
        val = _array[i]
        j = i - 1
        while (j >= 1 and _array[j] > val)
            _array[j + 1] = _array[j]
            j = j - 1
        end

        _array[j + 1] = val
    end

    PrintArray(_array)
    return -1 # Ignore this, just to hide false-positive errors.
end

function array_partition(_array, _low, _high)
    pivot = _array[_high]
    i = (_low - 1)
    for j = _low, _high
        if _array[j] < pivot then
            i = i + 1
            _array[i], _array[j] = _array[j], _array[i]
        end
    end
    _array[i + 1], _array[_high] = _array[_high], _array[i + 1]
    return (i + 1)
end
function Array_quickSort(_array, _low, _high, _init)
    if _low == nil then _low = 1 end
    if _high == nil then _high = #_array end
    if _init == nil then _init = true end

    if _low < _high then
        part = array_partition(_array, _low, _high)
        Array_quickSort(_array, _low, part - 1, false)
        Array_quickSort(_array, part + 1, _high, false)
    end
    if _init then
        PrintArray(_array)
    end
    return -1 # Ignore this, just to hide false-positive errors.
end

function array_getMax(_array)
    mx = _array[1]
    for i = 2, #_array
        if _array[i] > mx then
            mx = _array[i]
        end
    end
    return mx
end
function array_countsort(_array, _exp)
    output = table.clone(_array)
    count = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
    for i = 1, #_array
        cIndx = (math.floor(_array[i] / _exp) % 10) + 1
        count[cIndx] = count[cIndx] + 1
    end
    for i = 2, #count
        count[i] = count[i] + count[i - 1]
    end

    for i = #_array, 1, -1
        cIndx = (math.floor(_array[i] / _exp) % 10) + 1

        output[count[cIndx] - 1] = _array[i]
        count[cIndx] = count[cIndx] - 1
    end

    for i = 1, #_array
        _array[i] = output[i]
    end
end
function Array_radixSort(_array)
    m = array_getMax(_array)
    exp = 1
    while m / exp > 0
        array_countsort(_array, exp)
        exp = exp * 10
    end
    PrintArray(_array)
    return -1 # Ignore this, just to hide false-positive errors.
end


text:
# run on the RCI cluster
using DrWatson
@quickactivate
using PrettyTables
using Statistics
using FileIO, BSON, DataFrames

#DrWatson.datadir(path) = joinpath("/home/skvarvit/generativead/GenerativeAD.jl/data", path)
df_tabular = load(datadir("evaluation/tabular_eval.bson"))[:df];

scores = ["reconstruction-sampled", "reconstruction-mean", "jacodeco"]
models = ["vae", "vae_full", "vae_simple", "aae", "aae_full", "wae", "wae_full"]

mean_scores = Dict()
for model in models
	mdf = filter(r->r.modelname == model, df_tabular)
	mean_scores[model] = []
	for score in scores
		subdf = filter(r-> occursin(score, r.parameters), mdf)
		push!(mean_scores[model], mean(subdf.tst_eval_t))
	end
end
mean_scores
df_out = DataFrame()
df_out[!,:model] = models
for (i,score) in enumerate(scores)
	df_out[!,Symbol(score)] = map(m->mean_scores[m][i],models)
end

println(df_out)


text:
using GRAPE
using Pkg
using Documenter

# Generate examples
include("generate.jl")

DocMeta.setdocmeta!(GRAPE, :DocTestSetup, :(using GRAPE); recursive = true)

PROJECT_TOML = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))
VERSION = PROJECT_TOML["version"]
NAME = PROJECT_TOML["name"]
AUTHORS = join(PROJECT_TOML["authors"], ", ") * " and contributors"
GITHUB = "https://github.com/JuliaQuantumControl/GRAPE.jl"

println("Starting makedocs")

makedocs(;
    modules = [GRAPE],
    authors = AUTHORS,
    repo = "$GITHUB/blob/{commit}{path}#{line}",
    sitename = "GRAPE.jl",
    format = Documenter.HTML(;
        prettyurls = true,
        canonical = "https://juliaquantumcontrol.github.io/GRAPE.jl",
        assets = String[],
        mathengine = KaTeX(),
        footer="[$NAME.jl]($GITHUB) v$VERSION docs powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl)."
    ),
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Examples" => [
            "List of Examples" => "examples/index.md",
            "Example 1 (TLS)" => "examples/simple_state_to_state.md",
        ],
        "API" => "api.md",
        "History" => "history.md",
    ],
)

println("Finished makedocs")

rm(joinpath(@__DIR__, "build", "examples", ".gitignore"))

deploydocs(; repo = "github.com/JuliaQuantumControl/GRAPE.jl", devbranch = "master")


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

@testset "Test Array getindex" begin
    AG.read("ospy/data4/cropped_aster.img") do ds
        @testset "Dataset indexing" begin
            @testset "dims dropped correctly" begin
                @test typeof(ds[:, :, :]) <: Array{UInt8,3}
                @test typeof(ds[:, :, 1]) <: Array{UInt8,2}
                @test typeof(ds[:, 1, 1]) <: Array{UInt8,1}
                @test typeof(ds[1, 1, 1]) <: UInt8
                @test typeof(ds[:, :]) <: Array{UInt8,2}
                @test typeof(ds[:, 1]) <: Array{UInt8,1}
                @test typeof(ds[1, 1]) <: UInt8
            end
            @testset "range indexing" begin
                buffer = ds[1:AG.width(ds), 1:AG.height(ds), 1:1]

                @test typeof(buffer) <: Array{UInt8,3}
                total = sum(buffer)
                count = sum(buffer .> 0)

                # [python] total/count = 74985/863 ≈ 86.88876013904982
                @test total / count ≈ 86.88876013904982

                # [python] total / (band.shape[0] * band.shape[1])  -- in python
                #  = 27.741398446170923
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 27.741398446170923
            end
            @testset "colon indexing" begin
                buffer = ds[:, :, 1]

                @test buffer == ds[:, :]
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 86.88876013904982
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 27.741398446170923
            end
            @testset "int indexing" begin
                #   [python] band_1[35, 32] -> 130 (only one instance of 130)
                #       35, 32 in Python = 33, 36 in Julia
                #       130 = 0x82 in hex
                @test ds[33, 36, 1] == 0x82
            end
        end

        @testset "RasterBand indexing" begin
            band = AG.getband(ds, 1)
            @testset "dims dropped correctly" begin
                @test typeof(band[:, :]) <: Array{UInt8,2}
                @test typeof(band[:, 1]) <: Array{UInt8,1}
                @test typeof(band[1, 1]) <: UInt8
            end
            @testset "range indexing" begin
                buffer = band[1:AG.width(band), 1:AG.height(band)]
                @test typeof(buffer) <: Array{UInt8,2}
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 86.88876013904982
                @test total / (AG.height(band) * AG.width(band)) ≈ 27.741398446170923
            end
            @testset "colon indexing" begin
                buffer = band[:, :]
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 86.88876013904982
                @test total / (AG.height(band) * AG.width(band)) ≈ 27.741398446170923
            end
            @testset "int indexing" begin
                @test band[33, 36] == 0x82
            end
        end
    end
end

cp("ospy/data4/cropped_aster.img", "ospy/data4/cropped_aster_write.img"; force=true)

@testset "Test Array setindex" begin
    AG.read("ospy/data4/cropped_aster_write.img"; flags=AG.OF_Update) do ds
        @testset "Dataset setindex" begin
            @test ds[33, 36, 1] == 0x82
            ds[33, 36, 1] = 0x00
            @test ds[33, 36, 1] == 0x00
            ds[33:33, 36:36, 1:1] = reshape([0x01], 1, 1, 1)
            @test ds[33, 36, 1] == 0x01
            ds[33:33, 36:36, 1] = reshape([0x02], 1, 1)
            @test ds[33, 36, 1] == 0x02
            ds[33:33, 36, 1] = [0x03]
            @test ds[33, 36, 1] == 0x03
            ds[33, 36] = 0x04
            @test ds[33, 36] == 0x04
            ds[33:33, 36:36] = reshape([0x82], 1, 1)
            @test ds[33, 36] == 0x82
            buffer = ds[:, :]
            ds[:, :] = buffer .* 0x00
            @test sum(ds[:, :]) == 0x00
            ds[:, 1:35] = buffer[:, 1:35]
            ds[:, 36:end] = buffer[:, 36:end]
            @test sum(buffer) / sum(buffer .> 0) ≈ 86.88876013904982
            @test_throws DimensionMismatch ds[:, 501:end] = [1, 2, 3]
        end

        @testset "RasterBand setindex" begin
            band = AG.getband(ds, 1)
            @test band[33, 36] == 0x82
            band[33:33, 36] = [0x00]
            @test band[33, 36] == 0x00
            band[33, 36] = 0x01
            @test band[33, 36] == 0x01
            band[33:33, 36:36] = reshape([0x82], 1, 1)
            @test band[33, 36] == 0x82
            buffer = band[:, :]
            band[:, :] = buffer .* 0x00
            @test sum(band[:, :]) == 0x00
                ds[:, 1:35] = buffer[:, 1:35]
                ds[:, 36:end] = buffer[:, 36:end]
            @test sum(buffer) / sum(buffer .> 0) ≈ 86.88876013904982
            @test_throws DimensionMismatch band[:, 501:end] = [1, 2, 3]
        end
    end
end

@testset "Test Array constructor" begin
    AG.read("ospy/data4/cropped_aster_write.img"; flags=AG.OF_Update) do ds
        buffer = Array(ds)
        @test typeof(buffer) <: Array{UInt8,3}
        total = sum(buffer[:, :, 1:1])
        count = sum(buffer[:, :, 1:1] .> 0)
        @test total / count ≈ 86.88876013904982
        @test total / (AG.height(ds) * AG.width(ds)) ≈ 27.741398446170923
    end
end


text:
function Median()
    return Median(())
end

function copy(obj::Median)
    return jcall(obj, "copy", Median, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function evaluate(obj::AbstractUnivariateStatistic)
    return jcall(obj, "evaluate", jdouble, ())
end

function evaluate(obj::Median, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_data(obj::AbstractUnivariateStatistic)
    return jcall(obj, "getData", Vector{jdouble}, ())
end

function get_estimation_type(obj::Median)
    return jcall(obj, "getEstimationType", Percentile_EstimationType, ())
end

function get_kth_selector(obj::Median)
    return jcall(obj, "getKthSelector", KthSelector, ())
end

function get_na_n_strategy(obj::Median)
    return jcall(obj, "getNaNStrategy", NaNStrategy, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function set_data(obj::AbstractUnivariateStatistic, arg0::Vector{jdouble})
    return jcall(obj, "setData", void, (Vector{jdouble},), arg0)
end

function set_data(obj::AbstractUnivariateStatistic, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "setData", void, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function with_estimation_type(obj::Median, arg0::Percentile_EstimationType)
    return jcall(obj, "withEstimationType", Median, (Percentile_EstimationType,), arg0)
end

function with_kth_selector(obj::Median, arg0::KthSelector)
    return jcall(obj, "withKthSelector", Median, (KthSelector,), arg0)
end

function with_na_n_strategy(obj::Median, arg0::NaNStrategy)
    return jcall(obj, "withNaNStrategy", Median, (NaNStrategy,), arg0)
end


text:
@testset "factorial_ij" begin
    @test all(
        QSB.factorial_div(i, j) == factorial(big(i)) / factorial(big(j))
        for i in 0:40, j in 0:40
    )
end

@testset "z" begin
    x_range = -10:0.01:10
    p_range = -6:0.01:6
    @test all(
        QSB.z(x, p) == sqrt(2.)*(x + p*im)
        for x in x_range, p in p_range
    )

    xs = collect(x_range)
    ps = collect(p_range)
    @test QSB.z(xs, ps) == sqrt(2.).*(xs .+ ps' .* im)
end

@testset "gaussian_function" begin
    tol = 1e-14
    x_range = -10:0.01:10
    p_range = -6:0.01:6
    @test all(
        isapprox(QSB.gaussian_function(x, p), exp(-(x^2 + p^2)) / π, atol=tol)
        for x in x_range, p in p_range
    )

    xs = collect(x_range)
    ps = collect(p_range)
    @test isapprox(
        QSB.gaussian_function(xs, ps),
        exp.(-(xs.^2 .+ (ps').^2)) ./ π,
        atol=tol
    )
end

@testset "(-1)^i" begin
    @test all(
        QSB.neg_one_to_power_of(i) == (-1)^i
        for i in 1:50
    )
end

@testset "coefficient_of_wave_function" begin
    tol = 1e-14
    @test all(
        isapprox(
            QSB.coefficient_of_wave_function(m, n),
            ((n ≥ m) ? (-1)^(m-1) : (-1)^(n-1)) *
            sqrt(factorial(big(min(m, n))) / factorial(big(max(m, n)))),
            atol=tol
        )
        for m in 1:70, n in 1:70
    )
end

@testset "z_to_power" begin
    tol = 1e-14

    m_range = n_range = 1:70
    x_range = -1:0.1:1
    p_range = -0.6:0.1:0.6
    @test all(
        QSB.z_to_power(m, n, x, p) == (
            (n >= m) ? sqrt(2.)*(x - p*im) : sqrt(2.)*(x + p*im)
        )^(max(m, n) - min(m, n))
        for m in m_range, n in n_range, x in x_range, p in p_range
    )

    xs = collect(x_range)
    ps = collect(p_range)
    @test all(
        QSB.z_to_power(m, n, xs, ps) == QSB.z_to_power.(m, n, xs, ps')
        for m in m_range, n in n_range
    )
end

@testset "laguerre" begin
    m_range = n_range = 1:70
    x_range = -1:0.1:1
    p_range = -0.6:0.1:0.6

    result = Bool[]
    # n >= m
    for m in m_range, n in m:n_range.stop, x in x_range, p in p_range
        push!(
            result,
            abs(laguerre(m, n, x, p) - laguerre(m-1, n-m, abs2(sqrt(2.)*(x+p*im)))) < 1e-7
        )
    end
    # n < m
    for n in n_range, m in n:m_range.stop, x in x_range, p in p_range
        push!(
            result,
            abs(laguerre(m, n, x, p) - laguerre(n-1, m-n, abs2(sqrt(2.)*(x+p*im)))) < 1e-7
        )
    end
    @test all(result)

    xs = collect(x_range)
    ps = collect(p_range)
    @test all(
        laguerre(m, n, xs, ps) == laguerre.(m, n, xs, ps')
        for m in m_range, n in n_range
    )
end