text:
using CoupledSystems
using ForwardDiff
using Test

@testset "Named Variables" begin

    # default values for each variable
    x1 = rand()
    y1 = rand(10)
    z1 = rand(10, 10)

    # define each named variable
    @var x = x1
    @var y = y1
    @var z = z1

    # create tuple of named variables
    vars = (x, y, z)

    # combine values into a single vector
    v = combine(vars)

    # separate variables
    x2, y2, z2 = separate(vars, v)

    # check that variable values didn't change
    @test x1 == x2
    @test y1 == y2
    @test z1 == z2
end

@testset "ExplicitComponent" begin

    # This uses the paraboloid example from OpenMDAO

    # define variables and set defaults
    @var x = 0.0
    @var y = 0.0
    @var fxy = 0.0

    # construct paraboloid function (and define inputs and outputs)
    fin = (x, y)
    fout = (fxy,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct optional component output function
    f! = function(y, x)
        y[1] = (x[1]-3)^2 + x[1]*x[2] + (x[2]+4)^2 - 3
        return y
    end

    # construct optional component jacobian function
    df! = function(dydx, x)
        dydx[1,1] = 2*(x[1]-3) + x[2]
        dydx[1,2] = x[1] + 2*(x[2]+4)
        return dydx
    end

    # construct optional combined component output and jacobian function
    fdf! = function(y, dydx, x)
        f!(y, x)
        df!(dydx, x)
        return y, dydx
    end

    # define optional jacobian matrix storage
    dydx = zeros(1,2)

    # test a bunch of different constructors
    comp1 = ExplicitComponent(func, fin, fout, foutin; f=f!, df=df!, fdf=fdf!, dydx=dydx)
    comp2 = ExplicitComponent(func, fin, fout, foutin; df=df!, fdf=fdf!, dydx=dydx)
    comp3 = ExplicitComponent(func, fin, fout, foutin; df=df!, fdf=fdf!)
    comp4 = ExplicitComponent(func, fin, fout, foutin; df=df!)
    comp5 = ExplicitComponent(func, fin, fout, foutin; fdf=fdf!)
    comp6 = ExplicitComponent(fin, fout, foutin; f=f!, df=df!, fdf=fdf!, dydx=dydx)
    comp7 = ExplicitComponent(fin, fout, foutin; f=f!, df=df!, fdf=fdf!)
    comp8 = ExplicitComponent(fin, fout, foutin; f=f!, df=df!)
    comp9 = ExplicitComponent(fin, fout, foutin; fdf=fdf!)

    # make sure all the constructors yield the same results
    for comp in [comp2, comp3, comp4, comp5, comp6, comp7, comp8, comp9]

        # storage for outputs and jacobian
        y = zeros(1)
        dydx = zeros(1,2)

        x = rand(2)
        @test outputs(comp1, x) == outputs(comp, x)
        @test outputs!(comp1, y, x) == outputs!(comp, y, x)
        @test outputs!(comp1, x) == outputs!(comp, x)
        @test outputs!!(comp1, x) == outputs!!(comp, x)
        @test outputs(comp1) == outputs(comp)

        x = rand(2)
        @test jacobian(comp1, x) == jacobian(comp, x)
        @test jacobian!(comp1, dydx, x) == jacobian!(comp, dydx, x)
        @test jacobian!(comp1, x) == jacobian!(comp, x)
        @test jacobian!!(comp1, x) == jacobian!!(comp, x)
        @test jacobian(comp1) == jacobian(comp)

        x = rand(2)
        @test outputs_and_jacobian(comp1, x) == outputs_and_jacobian(comp, x)
        @test outputs_and_jacobian!(comp1, y, dydx, x) == outputs_and_jacobian!(comp, y, dydx, x)
        @test outputs_and_jacobian!(comp1, x) == outputs_and_jacobian!(comp, x)
        @test outputs_and_jacobian!!(comp1, x) == outputs_and_jacobian!!(comp, x)
        @test outputs_and_jacobian(comp1) == outputs_and_jacobian(comp)
    end
end

@testset "ExplicitComponent - input/output arguments" begin

    # define variables and set defaults
    @var x1 = [7.0; 8.0]
    @var x2 = 0.0
    @var y = [2.0;]
    @var fxy = 3.0

    # Main constraint is that all unmutable outputs must appear first (in the order specified in foutin),
    #  and then all the input parameters (specified in fin)

    # construct dummy function with all kinds of arguments
    fin = (x1,x2)
    fout = (fxy,)
    foutin = (y,)
    func = function(y,x1,x3)
        y[1] = 2*x1[1]+x1[2]
        fxy = 10. *x3[1]
        return fxy
    end

    comp1 = ExplicitComponent(func, fin, fout, foutin)

    X = [1.0; 4.0; 2.0] #Caution: this must be large enough to group all inputs
    Y = zeros(2)        #Caution: this must be large enough to group all outpus (unmutable and mutable)
    Y_res = [20.,6.]    #expected output

    @test outputs(comp1, X) == Y_res
    @test outputs!(comp1, Y, X) == Y_res
    @test Y == Y_res

    # ---------------------------------------------

    # construct dummy function with NO unmutable arguments, at leat none that we use
    fin = (x1,x2)
    fout = ()
    foutin = (y,)
    #let's use the same func

    comp2 = ExplicitComponent(func, fin, fout, foutin)

    X = [1.0; 4.0; 2.0] #Caution: this must be large enough to group all inputs
    Y = zeros(1)        #Caution: this must be large enough to group all outpus (unmutable and mutable)
    Y_res = [6.,]    #expected output

       @test outputs(comp2, X) == Y_res
    @test outputs!(comp2, Y, X) == Y_res
    @test Y == Y_res

    # ---------------------------------------------
    # define variables and set defaults
    @var x1 = [7.0; 8.0]
    @var x2 = 6.0
    @var y1 = [0.0;]
    @var y2 = [1.0; 2.0]
    @var f1 = 3.0
    @var f2 = [4.0;5.0]

    # construct dummy function with many in/out arguments, both vectors and scalars, in mixed orders
    fin = (x2,x1)
    fout = (f1,f2)
    foutin = (y1,y2)
    func = function(y1,y2,a1,a2)
        y1[1] = sum(a2)
        y2 .= 2 .* a2[1:2]
        f1 = 100 .* a1
        f2 = [a2[1];a1]
        return f1, f2
    end

    comp3 = ExplicitComponent(func, fin, fout, foutin)

    A1 = 2.
    A2 = [1.0; 4.0]
    Y1 = zeros(1)
    Y2 = zeros(2)
    F1,F2 = func(Y1,Y2,A1,A2)

    X = [2.0;  1.0; 4.0] #Caution: this must be large enough to group all inputs
    Y = zeros(6)        #Caution: this must be large enough to group all outpus (unmutable and mutable)
    Y_res = vcat(F1,F2,Y1,Y2)    #expected output

    @test outputs(comp3, X) == Y_res
    @test outputs!(comp3, Y, X) == Y_res
    @test Y == Y_res
end

@testset "ExplicitComponent - Derivatives" begin

    # This uses the paraboloid example from OpenMDAO

    # define variables and set defaults
    @var x = 0.0
    @var y = 0.0
    @var fxy = 0.0

    # construct template function
    fin = (x, y)
    fout = (fxy,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct optional component output function
    f! = function(y, x)
        y[1] = (x[1]-3)^2 + x[1]*x[2] + (x[2]+4)^2 - 3
        return y
    end

    # construct optional component jacobian function
    df! = function(dydx, x)
        dydx[1,1] = 2*(x[1]-3) + x[2]
        dydx[1,2] = x[1] + 2*(x[2]+4)
        return dydx
    end

    # construct optional combined component output and jacobian function
    fdf! = function(y, dydx, x)
        f!(y, x)
        df!(dydx, x)
        return y, dydx
    end

    x = zeros(2)
    y = zeros(1)
    dydx = zeros(1,2)
    comp = ExplicitComponent(func, fin, fout, foutin; f=f!, df=df!, fdf=fdf!)

    comp_FAD = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())
    comp_RAD = ExplicitComponent(func, fin, fout, foutin; deriv=ReverseAD())
    comp_FFD = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardFD())
    comp_CFD = ExplicitComponent(func, fin, fout, foutin; deriv=CentralFD())
    comp_XFD = ExplicitComponent(func, fin, fout, foutin; deriv=ComplexFD())

    for dcomp in [comp_FAD, comp_RAD, comp_FFD, comp_CFD, comp_XFD]

        x = rand(2)
        @test isapprox(outputs(comp, x), outputs(dcomp, x), atol=1e-6)
        @test isapprox(outputs!(comp, y, x), outputs!(dcomp, y, x), atol=1e-6)
        @test isapprox(outputs!(comp, x), outputs!(dcomp, x), atol=1e-6)
        @test isapprox(outputs!!(comp, x), outputs!!(dcomp, x), atol=1e-6)
        @test isapprox(outputs(comp), outputs(dcomp), atol=1e-6)

        x = rand(2)
        @test isapprox(jacobian(comp, x), jacobian(dcomp, x), atol=1e-6)
        @test isapprox(jacobian!(comp, dydx, x), jacobian!(dcomp, dydx, x), atol=1e-6)
        @test isapprox(jacobian!(comp, x), jacobian!(dcomp, x), atol=1e-6)
        @test isapprox(jacobian!!(comp, x), jacobian!!(dcomp, x), atol=1e-6)
        @test isapprox(jacobian(comp), jacobian(dcomp), atol=1e-6)

        x = rand(2)
        @test all(isapprox.(outputs_and_jacobian(comp, x), outputs_and_jacobian(dcomp, x), atol=1e-6))
        @test all(isapprox.(outputs_and_jacobian!(comp, y, dydx, x), outputs_and_jacobian!(dcomp, y, dydx, x), atol=1e-6))
        @test all(isapprox.(outputs_and_jacobian!(comp, x), outputs_and_jacobian!(dcomp, x), atol=1e-6))
        @test all(isapprox.(outputs_and_jacobian!!(comp, x), outputs_and_jacobian!!(dcomp, x), atol=1e-6))
        @test all(isapprox.(outputs_and_jacobian(comp), outputs_and_jacobian(dcomp), atol=1e-6))
    end
end

@testset "Explicit to Implicit" begin

    # This uses the paraboloid example from OpenMDAO

    # define variables and set defaults
    @var x = 0.0
    @var y = 0.0
    @var fxy = 0.0

    # construct template function
    fin = (x, y)
    fout = (fxy,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct explicit component
    xcomp = ExplicitComponent(func, fin, fout, foutin)

    # construct implicit component from explicit component
    icomp = ImplicitComponent(xcomp)

    r = zeros(1)
    drdx = zeros(1,2)
    drdy = zeros(1,1)

    x = rand(2)
    y = rand(1)
    r1 = residuals(icomp, x, y)
    r2 = residuals!(icomp, r, x, y)
    r3 = residuals!(icomp, x, y)
    r4 = residuals!!(icomp, x, y)
    r5 = residuals(icomp)
    @test r1 == r2 == r3 == r4 == r5

    x = rand(2)
    y = rand(1)
    drdx1 = residual_input_jacobian(icomp, x, y)
    drdx2 = residual_input_jacobian!(icomp, drdx, x, y)
    drdx3 = residual_input_jacobian!(icomp, x, y)
    drdx4 = residual_input_jacobian!!(icomp, x, y)
    drdx5 = residual_input_jacobian(icomp)
    @test drdx1 == drdx2 == drdx3 == drdx4 == drdx5

    x = rand(2)
    y = rand(1)
    drdy1 = residual_output_jacobian(icomp, x, y)
    drdy2 = residual_output_jacobian!(icomp, drdy, x, y)
    drdy3 = residual_output_jacobian!(icomp, x, y)
    drdy4 = residual_output_jacobian!!(icomp, x, y)
    drdy5 = residual_output_jacobian(icomp)
    @test drdy1 == drdy2 == drdy3 == drdy4 == drdy5

    x = rand(2)
    y = rand(1)
    r1, drdx1 = residuals_and_input_jacobian(icomp, x, y)
    r2, drdx2 = residuals_and_input_jacobian!(icomp, r, drdx, x, y)
    r3, drdx3 = residuals_and_input_jacobian!(icomp, x, y)
    r4, drdx4 = residuals_and_input_jacobian!!(icomp, x, y)
    r5, drdx5 = residuals_and_input_jacobian(icomp)
    @test r1 == r2 == r3 == r4 == r5
    @test drdx1 == drdx2 == drdx3 == drdx4 == drdx5

    x = rand(2)
    y = rand(1)
    r1, drdy1 = residuals_and_output_jacobian(icomp, x, y)
    r2, drdy2 = residuals_and_output_jacobian!(icomp, r, drdy, x, y)
    r3, drdy3 = residuals_and_output_jacobian!(icomp, x, y)
    r4, drdy4 = residuals_and_output_jacobian!!(icomp, x, y)
    r5, drdy5 = residuals_and_output_jacobian(icomp)
    @test r1 == r2 == r3 == r4 == r5
    @test drdy1 == drdy2 == drdy3 == drdy4 == drdy5

    x = rand(2)
    y = rand(1)
    r1, drdx1, drdy1 = residuals_and_jacobians(icomp, x, y)
    r2, drdx2, drdy2 = residuals_and_jacobians!(icomp, r, drdx, drdy, x, y)
    r3, drdx3, drdy3 = residuals_and_jacobians!(icomp, x, y)
    r4, drdx4, drdy4 = residuals_and_jacobians!!(icomp, x, y)
    r5, drdx5, drdy5 = residuals_and_jacobians(icomp)
    @test r1 == r2 == r3 == r4 == r5
    @test drdx1 == drdx2 == drdx3 == drdx4 == drdx5
    @test drdy1 == drdy2 == drdy3 == drdy4 == drdy5
end

@testset "ExplicitSystem" begin

    # define system variables and defaults
    @var x = 0.0
    @var y = 0.0
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var fp = 0.0
    @var fq = 0.0
    @var ft = zeros(2)

    # construct template function for paraboloid component
    fin = (x, y)
    fout = (fp,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct paraboloid component
    paraboloid = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Second Component: Quadratic

    # construct template function for quadratic component
    fin = (a, b, c, fp)
    fout = (fq,)
    foutin = ()
    func = function(a, b, c, fp)
        fq = a*fp^2 + b*fp + c*fp + 1
        return fq
    end

    # construct quadratic component
    quadratic = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Third Component: Trigometric Functions

    # construct template function for trigometric component
    fin = (fp, fq)
    fout = ()
    foutin = (ft,)
    func = function(ft, fp, fq)
        ft[1] = sin(fp)
        ft[2] = cos(fq)
        return ft
    end

    # construct trigometric component
    trigometric = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Combined System

    argin = (x, y, a, b, c)
    argout = (ft,)
    components = (paraboloid, quadratic, trigometric)

    sys = ExplicitSystem(components, argin, argout)

    y = zeros(2)
    dydx = zeros(2,5)

    x = rand(5)
    y1 = outputs!(sys, y, x)
    y2 = outputs!(sys, x)
    y3 = outputs!!(sys, x)
    y4 = outputs!!!(sys, x)
    y5 = outputs(sys, x)
    @test y1 == y2 == y3 == y4 == y5

    x = rand(5)
    dydx1 = jacobian!(sys, dydx, x)
    dydx2 = jacobian!(sys, x)
    dydx3 = jacobian!!(sys, x)
    dydx4 = jacobian!!!(sys, x)
    dydx5 = jacobian(sys, x)
    @test dydx1 == dydx2 == dydx3 == dydx4 == dydx5

    x = rand(5)
    y1, dydx1 = outputs_and_jacobian!(sys, y, dydx, x)
    y2, dydx2 = outputs_and_jacobian!(sys, x)
    y3, dydx3 = outputs_and_jacobian!!(sys, x)
    y4, dydx4 = outputs_and_jacobian!!!(sys, x)
    y5, dydx5 = outputs_and_jacobian(sys, x)
    @test y1 == y2 == y3 == y4 == y5
    @test dydx1 == dydx2 == dydx3 == dydx4 == dydx5
end

@testset "ExplicitSystem - Derivatives" begin

    # define system variables and defaults
    @var x = 0.0
    @var y = 0.0
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var fp = 0.0
    @var fq = 0.0
    @var ft = zeros(2)

    # construct template function for paraboloid component
    fin = (x, y)
    fout = (fp,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct paraboloid component
    paraboloid = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Second Component: Quadratic

    # construct template function for quadratic component
    fin = (a, b, c, fp)
    fout = (fq,)
    foutin = ()
    func = function(a, b, c, fp)
        fq = a*fp^2 + b*fp + c*fp + 1
        return fq
    end

    # construct quadratic component
    quadratic = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Third Component: Trigometric Functions

    # construct template function for trigometric component
    fin = (fp, fq)
    fout = ()
    foutin = (ft,)
    func = function(ft, fp, fq)
        ft[1] = sin(fp)
        ft[2] = cos(fq)
        return ft
    end

    # construct trigometric component
    trigometric = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Combined System
    argin = (x, y, a, b, c)
    argout = (ft,)
    components = (paraboloid, quadratic, trigometric)

    sys = ExplicitSystem(components, argin, argout)

    # initialize storage for outputs and jacobian
    y = zeros(2)
    dydx = zeros(2,5)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_f = jacobian!(sys, dydx, x, Forward())
    @test isapprox(dydx_ad, dydx_f)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_r = jacobian!(sys, dydx, x, Reverse())
    @test isapprox(dydx_ad, dydx_r)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_f = jacobian!(sys, x, Forward())
    @test isapprox(dydx_ad, dydx_f)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_r = jacobian!(sys, x, Reverse())
    @test isapprox(dydx_ad, dydx_r)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_f = jacobian!!(sys, x, Forward())
    @test isapprox(dydx_ad, dydx_f)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_r = jacobian!!(sys, x, Reverse())
    @test isapprox(dydx_ad, dydx_r)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_f = jacobian!!!(sys, x, Forward())
    @test isapprox(dydx_ad, dydx_f)

    x = rand(5)
    dydx_ad = ForwardDiff.jacobian(x -> outputs(sys, x), x)
    dydx_r = jacobian!!!(sys, x, Reverse())
    @test isapprox(dydx_ad, dydx_r)
end

@testset "ImplicitComponent" begin

    # This uses the quadratic example from OpenMDAO

    # variables
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var x = 0.0

    # template residual function
    argin = (a, b, c)
    argout = (x,)
    func = function(r, a, b, c, x)
        r[1] = a*x^2 + b*x + c
        return r
    end

    # vector form of residual function
    f! = function(r, x, y)
        r[1] = x[1]*y[1]^2 + x[2]*y[1] + x[3]
        return r
    end

    # jacobian wrt inputs
    dfdx! = function(drdx, x, y)
        drdx[1,1] = y[1]^2
        drdx[1,2] = y[1]
        drdx[1,3] = 1
        return drdx
    end

    # jacobian wrt outputs
    dfdy! = function(drdy, x, y)
        drdy[1,1] = 2*x[1]*y[1] + x[2]
        return drdy
    end

    # residual and jacobian wrt inputs
    fdfdx! = function(r, drdx, x, y)
        f!(r, x, y)
        dfdx!(drdx, x, y)
        return r, drdx
    end

    # residual and jacobian wrt outputs
    fdfdy! = function(r, drdy, x, y)
        f!(r, x, y)
        dfdy!(drdy, x, y)
        return r, drdy
    end

    # residual and jacobians
    fdf! = function(r, drdx, drdy, x, y)
        f!(r, x, y)
        dfdx!(drdx, x, y)
        dfdy!(drdy, x, y)
        return r, drdx, drdy
    end

    # storage for results
    r = zeros(1)
    drdx = zeros(1,3)
    drdy = zeros(1,1)

    comp1 = ImplicitComponent(func, argin, argout, r; f=f!, dfdx=dfdx!, dfdy=dfdy!, fdfdx=fdfdx!,
        fdfdy=fdfdy!, fdf=fdf!, drdx=drdx, drdy=drdy)
    comp2 = ImplicitComponent(func, argin, argout, r; dfdx=dfdx!, dfdy=dfdy!, fdfdx=fdfdx!,
        fdfdy=fdfdy!, fdf=fdf!, drdx=drdx, drdy=drdy)
    comp3 = ImplicitComponent(func, argin, argout, r; dfdx=dfdx!, dfdy=dfdy!, fdfdx=fdfdx!,
        fdfdy=fdfdy!, fdf=fdf!)
    comp4 = ImplicitComponent(func, argin, argout, r; dfdx=dfdx!, dfdy=dfdy!)
    comp5 = ImplicitComponent(func, argin, argout; fdfdx=fdfdx!, fdfdy=fdfdy!)
    comp6 = ImplicitComponent(func, argin, argout; fdf=fdf!)

    comp7 = ImplicitComponent(argin, argout; f=f!, dfdx=dfdx!, dfdy=dfdy!, fdfdx=fdfdx!,
        fdfdy=fdfdy!, fdf=fdf!, drdx=drdx, drdy=drdy)
    comp8 = ImplicitComponent(argin, argout; f=f!, dfdx=dfdx!, dfdy=dfdy!, fdfdx=fdfdx!,
        fdfdy=fdfdy!, fdf=fdf!)
    comp9 = ImplicitComponent(argin, argout; f=f!, dfdx=dfdx!, dfdy=dfdy!)
    comp10 = ImplicitComponent(argin, argout; fdfdx=fdfdx!, fdfdy=fdfdy!)
    comp11 = ImplicitComponent(argin, argout; fdf=fdf!)

    for comp in [comp2, comp3, comp4, comp5, comp6, comp7, comp8, comp9, comp10, comp11]
        x = rand(3)
        y = rand(1)
        @test residuals(comp1, x, y) == residuals(comp, x, y)
        @test residuals!(comp1, r, x, y) == residuals!(comp, r, x, y)
        @test residuals!(comp1, x, y) == residuals!(comp, x, y)
        @test residuals!!(comp1, x, y) == residuals!!(comp, x, y)
        @test residuals(comp1) == residuals(comp)

        x = rand(3)
        y = rand(1)
        @test residual_input_jacobian(comp1, x, y) == residual_input_jacobian(comp, x, y)
        @test residual_input_jacobian!(comp1, drdx, x, y) == residual_input_jacobian!(comp, drdx, x, y)
        @test residual_input_jacobian!(comp1, x, y) == residual_input_jacobian!(comp, x, y)
        @test residual_input_jacobian!!(comp1, x, y) == residual_input_jacobian!!(comp, x, y)
        @test residual_input_jacobian(comp1) == residual_input_jacobian(comp)

        x = rand(3)
        y = rand(1)
        @test residual_output_jacobian(comp1, x, y) == residual_output_jacobian(comp, x, y)
        @test residual_output_jacobian!(comp1, drdy, x, y) == residual_output_jacobian!(comp, drdy, x, y)
        @test residual_output_jacobian!(comp1, x, y) == residual_output_jacobian!(comp, x, y)
        @test residual_output_jacobian!!(comp1, x, y) == residual_output_jacobian!!(comp, x, y)
        @test residual_output_jacobian(comp1) == residual_output_jacobian(comp)

        x = rand(3)
        y = rand(1)
        @test residuals_and_input_jacobian(comp1, x, y) == residuals_and_input_jacobian(comp, x, y)
        @test residuals_and_input_jacobian!(comp1, r, drdx, x, y) == residuals_and_input_jacobian!(comp, r, drdx, x, y)
        @test residuals_and_input_jacobian!(comp1, x, y) == residuals_and_input_jacobian!(comp, x, y)
        @test residuals_and_input_jacobian!!(comp1, x, y) == residuals_and_input_jacobian!!(comp, x, y)
        @test residuals_and_input_jacobian(comp1) == residuals_and_input_jacobian(comp)

        x = rand(3)
        y = rand(1)
        @test residuals_and_output_jacobian(comp1, x, y) == residuals_and_output_jacobian(comp, x, y)
        @test residuals_and_output_jacobian!(comp1, r, drdy, x, y) == residuals_and_output_jacobian!(comp, r, drdy, x, y)
        @test residuals_and_output_jacobian!(comp1, x, y) == residuals_and_output_jacobian!(comp, x, y)
        @test residuals_and_output_jacobian!!(comp1, x, y) == residuals_and_output_jacobian!!(comp, x, y)
        @test residuals_and_output_jacobian(comp1) == residuals_and_output_jacobian(comp)

        x = rand(3)
        y = rand(1)
        @test residuals_and_jacobians(comp1, x, y) == residuals_and_jacobians(comp, x, y)
        @test residuals_and_jacobians!(comp1, r, drdx, drdy, x, y) == residuals_and_jacobians!(comp, r, drdx, drdy, x, y)
        @test residuals_and_jacobians!(comp1, x, y) == residuals_and_jacobians!(comp, x, y)
        @test residuals_and_jacobians!!(comp1, x, y) == residuals_and_jacobians!!(comp, x, y)
        @test residuals_and_jacobians(comp1) == residuals_and_jacobians(comp)
    end
end

@testset "ImplicitComponent - Derivatives" begin

    # This uses the quadratic example from OpenMDAO

    # variables
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var x = 0.0

    # template residual function
    fin = (a, b, c)
    fout = (x,)
    func = function(r, a, b, c, x)
        r[1] = a*x^2 + b*x + c
        return r
    end

    # jacobian wrt inputs
    dfdx! = function(drdx, x, y)
        drdx[1,1] = y[1]^2
        drdx[1,2] = y[1]
        drdx[1,3] = 1
        return drdx
    end

    # jacobian wrt outputs
    dfdy! = function(drdy, x, y)
        drdy[1,1] = 2*x[1]*y[1] + x[2]
        return drdy
    end

    comp = ImplicitComponent(func, fin, fout; dfdx=dfdx!, dfdy=dfdy!)

    comp_FAD = ImplicitComponent(func, fin, fout; xderiv=ForwardAD(), yderiv=ForwardAD())
    comp_RAD = ImplicitComponent(func, fin, fout; xderiv=ReverseAD(), yderiv=ReverseAD())
    comp_FFD = ImplicitComponent(func, fin, fout; xderiv=ForwardFD(), yderiv=ForwardFD())
    comp_CFD = ImplicitComponent(func, fin, fout; xderiv=CentralFD(), yderiv=CentralFD())
    comp_XFD = ImplicitComponent(func, fin, fout; xderiv=ComplexFD(), yderiv=ComplexFD())

    # storage for outputs
    r = zeros(1)
    drdx = zeros(1, 3)
    drdy = zeros(1, 1)

    for dcomp in [comp_FAD, comp_RAD, comp_FFD, comp_CFD, comp_XFD]

        x = rand(3)
        y = rand(1)
        @test isapprox(residuals(comp, x, y), residuals(dcomp, x, y), atol=1e-6)
        @test isapprox(residuals!(comp, r, x, y), residuals!(dcomp, r, x, y), atol=1e-6)
        @test isapprox(residuals!(comp, x, y), residuals!(dcomp, x, y), atol=1e-6)
        @test isapprox(residuals!!(comp, x, y), residuals!!(dcomp, x, y), atol=1e-6)
        @test isapprox(residuals(comp), residuals(dcomp), atol=1e-6)

        x = rand(3)
        y = rand(1)
        @test isapprox(residual_input_jacobian(comp, x, y), residual_input_jacobian(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_input_jacobian!(comp, drdx, x, y), residual_input_jacobian!(dcomp, drdx, x, y), atol=1e-6)
        @test isapprox(residual_input_jacobian!(comp, x, y), residual_input_jacobian!(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_input_jacobian!!(comp, x, y), residual_input_jacobian!!(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_input_jacobian(comp), residual_input_jacobian(dcomp), atol=1e-6)

        x = rand(3)
        y = rand(1)
        @test isapprox(residual_output_jacobian(comp, x, y), residual_output_jacobian(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_output_jacobian!(comp, drdy, x, y), residual_output_jacobian!(dcomp, drdy, x, y), atol=1e-6)
        @test isapprox(residual_output_jacobian!(comp, x, y), residual_output_jacobian!(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_output_jacobian!!(comp, x, y), residual_output_jacobian!!(dcomp, x, y), atol=1e-6)
        @test isapprox(residual_output_jacobian(comp), residual_output_jacobian(dcomp), atol=1e-6)

        x = rand(3)
        y = rand(1)
        @test all(isapprox.(residuals_and_input_jacobian(comp, x, y), residuals_and_input_jacobian(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_input_jacobian!(comp, r, drdx, x, y), residuals_and_input_jacobian!(dcomp, r, drdx, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_input_jacobian!(comp, x, y), residuals_and_input_jacobian!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_input_jacobian!!(comp, x, y), residuals_and_input_jacobian!!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_input_jacobian(comp), residuals_and_input_jacobian(dcomp), atol=1e-6))

        x = rand(3)
        y = rand(1)
        @test all(isapprox.(residuals_and_output_jacobian(comp, x, y), residuals_and_output_jacobian(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_output_jacobian!(comp, r, drdy, x, y), residuals_and_output_jacobian!(dcomp, r, drdy, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_output_jacobian!(comp, x, y), residuals_and_output_jacobian!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_output_jacobian!!(comp, x, y), residuals_and_output_jacobian!!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_output_jacobian(comp), residuals_and_output_jacobian(dcomp), atol=1e-6))

        x = rand(3)
        y = rand(1)
        @test all(isapprox.(residuals_and_jacobians(comp, x, y), residuals_and_jacobians(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_jacobians!(comp, r, drdx, drdy, x, y), residuals_and_jacobians!(dcomp, r, drdx, drdy, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_jacobians!(comp, x, y), residuals_and_jacobians!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_jacobians!!(comp, x, y), residuals_and_jacobians!!(dcomp, x, y), atol=1e-6))
        @test all(isapprox.(residuals_and_jacobians(comp), residuals_and_jacobians(dcomp), atol=1e-6))

    end
end

@testset "Implicit to Explicit" begin

    # This uses the quadratic example from OpenMDAO

    # variables
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var x = 0.0

    # template residual function
    fin = (a, b, c)
    fout = (x,)
    func = function(r, a, b, c, x)
        r[1] = a*x^2 + b*x + c
        return r
    end

    # create implicit component
    icomp = ImplicitComponent(func, fin, fout)

    # create explicit component
    xcomp = ExplicitComponent(icomp)

    # inputs
    x = [2, 0, -1]

    # storage for outputs
    y = zeros(1)
    dydx = zeros(1, 3)

    y1 = outputs!(xcomp, y, x)
    y2 = outputs!(xcomp, x)
    y3 = outputs!!(xcomp, x)
    y4 = outputs!!!(xcomp, x)
    y5 = outputs(xcomp, x)
    @test isapprox(y1, y2, atol=1e-6)
    @test isapprox(y1, y3, atol=1e-6)
    @test isapprox(y1, y4, atol=1e-6)
    @test isapprox(y1, y5, atol=1e-6)

    dydx1 = jacobian!(xcomp, dydx, x)
    dydx2 = jacobian!(xcomp, x)
    dydx3 = jacobian!!(xcomp, x)
    dydx4 = jacobian!!!(xcomp, x)
    dydx5 = jacobian(xcomp, x)
    @test isapprox(dydx1, dydx2, atol=1e-6)
    @test isapprox(dydx1, dydx3, atol=1e-6)
    @test isapprox(dydx1, dydx4, atol=1e-6)
    @test isapprox(dydx1, dydx5, atol=1e-6)

    y1, dydx1 = outputs_and_jacobian!(xcomp, y, dydx, x)
    y2, dydx2 = outputs_and_jacobian!(xcomp, x)
    y3, dydx3 = outputs_and_jacobian!!(xcomp, x)
    y4, dydx4 = outputs_and_jacobian!!!(xcomp, x)
    y5, dydx5 = outputs_and_jacobian(xcomp, x)
    @test isapprox(y1, y2, atol=1e-6)
    @test isapprox(y1, y3, atol=1e-6)
    @test isapprox(y1, y4, atol=1e-6)
    @test isapprox(y1, y5, atol=1e-6)
    @test isapprox(dydx1, dydx2, atol=1e-6)
    @test isapprox(dydx1, dydx3, atol=1e-6)
    @test isapprox(dydx1, dydx4, atol=1e-6)
    @test isapprox(dydx1, dydx5, atol=1e-6)

end

@testset "ImplicitSystem" begin

    # We use an explicit system here, modified to become implicit
    # define system variables and defaults
    @var x = 0.0
    @var y = 0.0
    @var a = 0.0
    @var b = 0.0
    @var c = 0.0
    @var fp = 0.0
    @var fq = 0.0
    @var ft = zeros(2)

    # construct template function for paraboloid component
    fin = (x, y)
    fout = (fp,)
    foutin = ()
    func = (x,y) -> (x - 3)^2 + x*y + (y + 4)^2 - 3

    # construct paraboloid component
    paraboloid = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Second Component: Quadratic

    # construct template function for quadratic component
    fin = (a, b, c, fp)
    fout = (fq,)
    foutin = ()
    func = function(a, b, c, fp)
        fq = a*fp^2 + b*fp + c*fp + 1
        return fq
    end

    # construct quadratic component
    quadratic = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # Third Component: Trigometric Functions

    # construct template function for trigometric component
    fin = (fp, fq)
    fout = ()
    foutin = (ft,)
    func = function(ft, fp, fq)
        ft[1] = sin(fp)
        ft[2] = cos(fq)
        return ft
    end

    # construct trigometric component
    trigometric = ExplicitComponent(func, fin, fout, foutin; deriv=ForwardAD())

    # components in implicit system
    components = (paraboloid, quadratic, trigometric)

    # arguments to implicit system
    argin = (x,y,a,b,c)

    # construct implicit system
    sys = ImplicitSystem(components, argin)

    r = zeros(4)
    drdx = zeros(4, 5)
    drdy = zeros(4, 4)

    x = rand(5)
    y = rand(4)
    r1 = residuals!(sys, r, x, y)
    r2 = residuals!(sys, x, y)
    r3 = residuals!!(sys, x, y)
    r4 = residuals!!!(sys, x, y)
    r5 = residuals(sys, x, y)
    @test isapprox(r1, r2)
    @test isapprox(r1, r3)
    @test isapprox(r1, r4)
    @test isapprox(r1, r5)

    x = rand(5)
    y = rand(4)
    drdx1 = residual_input_jacobian!(sys, drdx, x, y)
    drdx2 = residual_input_jacobian!(sys, x, y)
    drdx3 = residual_input_jacobian!!(sys, x, y)
    drdx4 = residual_input_jacobian!!!(sys, x, y)
    drdx5 = residual_input_jacobian(sys, x, y)
    @test isapprox(drdx1, drdx2)
    @test isapprox(drdx1, drdx3)
    @test isapprox(drdx1, drdx4)
    @test isapprox(drdx1, drdx5)

    x = rand(5)
    y = rand(4)
    drdy1 = residual_output_jacobian!(sys, drdy, x, y)
    drdy2 = residual_output_jacobian!(sys, x, y)
    drdy3 = residual_output_jacobian!!(sys, x, y)
    drdy4 = residual_output_jacobian!!!(sys, x, y)
    drdy5 = residual_output_jacobian(sys, x, y)
    @test isapprox(drdy1, drdy2)
    @test isapprox(drdy1, drdy3)
    @test isapprox(drdy1, drdy4)
    @test isapprox(drdy1, drdy5)

    x = rand(5)
    y = rand(4)
    r1, drdx1 = residuals_and_input_jacobian!(sys, r, drdx, x, y)
    r2, drdx2 = residuals_and_input_jacobian!(sys, x, y)
    r3, drdx3 = residuals_and_input_jacobian!!(sys, x, y)
    r4, drdx4 = residuals_and_input_jacobian!!!(sys, x, y)
    r5, drdx5 = residuals_and_input_jacobian(sys, x, y)
    @test isapprox(r1, r2)
    @test isapprox(r1, r3)
    @test isapprox(r1, r4)
    @test isapprox(r1, r5)
    @test isapprox(drdx1, drdx2)
    @test isapprox(drdx1, drdx3)
    @test isapprox(drdx1, drdx4)
    @test isapprox(drdx1, drdx5)

    x = rand(5)
    y = rand(4)
    r1, drdy1 = residuals_and_output_jacobian!(sys, r, drdy, x, y)
    r2, drdy2 = residuals_and_output_jacobian!(sys, x, y)
    r3, drdy3 = residuals_and_output_jacobian!!(sys, x, y)
    r4, drdy4 = residuals_and_output_jacobian!!!(sys, x, y)
    r5, drdy5 = residuals_and_output_jacobian(sys, x, y)
    @test isapprox(r1, r2)
    @test isapprox(r1, r3)
    @test isapprox(r1, r4)
    @test isapprox(r1, r5)
    @test isapprox(drdy1, drdy2)
    @test isapprox(drdy1, drdy3)
    @test isapprox(drdy1, drdy4)
    @test isapprox(drdy1, drdy5)

    x = rand(5)
    y = rand(4)
    r1, drdx1, drdy1 = residuals_and_jacobians!(sys, r, drdx, drdy, x, y)
    r2, drdx2, drdy2 = residuals_and_jacobians!(sys, x, y)
    r3, drdx3, drdy3 = residuals_and_jacobians!!(sys, x, y)
    r4, drdx4, drdy4 = residuals_and_jacobians!!!(sys, x, y)
    r5, drdx5, drdy5 = residuals_and_jacobians(sys, x, y)
    @test isapprox(r1, r2)
    @test isapprox(r1, r3)
    @test isapprox(r1, r4)
    @test isapprox(r1, r5)
    @test isapprox(drdx1, drdx2)
    @test isapprox(drdx1, drdx3)
    @test isapprox(drdx1, drdx4)
    @test isapprox(drdx1, drdx5)
    @test isapprox(drdy1, drdy2)
    @test isapprox(drdy1, drdy3)
    @test isapprox(drdy1, drdy4)
    @test isapprox(drdy1, drdy5)

end

@testset "Sellar" begin

    # --- System Variables --- #
    @var x = 0.0
    @var y1 = 0.0
    @var y2 = 0.0
    @var z1 = 0.0
    @var z2 = 0.0
    @var f = 0.0
    @var g1 = 0.0
    @var g2 = 0.0

    # --- Define Discipline 1 --- #

    # describe function inputs/outputs
    fin = (x, y2, z1, z2) # inputs
    fout = (y1,) # normal outputs
    fout_m = () # in-place (mutating) outputs

    # create function for discipline 1
    function f_d1(x, y2, z1, z2)
        y1 = z1^2 + z2 + x - 0.2*y2
        return y1
    end

    # construct explicit component for discipline 1
    d1 = ExplicitComponent(f_d1, fin, fout, fout_m; deriv=ForwardFD())

    # --- Define Discipline 2 --- #

    # describe function inputs/outputs
    fin = (y1, z1, z2) # inputs
    fout = (y2,) # normal outputs
    fout_m = () # in-place (mutating) outputs

    # create function for discipline 2
    function f_d2(y1, z1, z2)
        y2 = sqrt(y1) + z1 + z2
        return y2
    end

    # construct explicit component for discipline 2
    d2 = ExplicitComponent(f_d2, fin, fout, fout_m; deriv=ForwardFD())

    # --- Define Objective --- #

    # describe function inputs/outputs
    fin = (x, y1, y2, z1) # inputs
    fout = (f,) # normal outputs
    fout_m = () # in-place (mutating) outputs

    # create objective function
    function f_obj(x, y1, y2, z1)
        f = x^2 + z1 + y1 + exp(-y2) # objective
        return f
    end

    # construct explicit component for objective
    obj = ExplicitComponent(f_obj, fin, fout, fout_m; deriv=ForwardFD())

    # --- Define Constraint 1 --- #

    # describe function inputs/outputs
    fin = (y1,) # inputs
    fout = (g1,) # normal outputs
    fout_m = () # in-place (mutating) outputs

    # define first constraint function
    function f_c1(y1)
        g1 = 3.16 - y1
        return g1
    end

    # construct explicit component for constraint 1
    c1 = ExplicitComponent(f_c1, fin, fout, fout_m; deriv=ForwardFD())

    # --- Define Constraint 2 --- #

    # describe function inputs/outputs
    fin = (y2,) # inputs
    fout = (g2,) # normal outputs
    fout_m = () # in-place (mutating) outputs

    # define second constraint function
    function f_c2(y2)
        g2 = y2 - 24.0
        return g2
    end

    # construct explicit component for constraint 2
    c2 = ExplicitComponent(f_c2, fin, fout, fout_m; deriv=ForwardFD())

    # --- Construct Implicit System --- #

    # components in the implicit system
    components_isys = (d1, d2)

    # inputs to the implicit system
    inputs_isys = (x, z1, z2)

    # implicit system construction
    isys = ImplicitSystem(components_isys, inputs_isys)

    # --- Construct Explicit Component from Implicit System --- #
    mda = ExplicitComponent(isys; solver=Newton())

    # --- Construct Explicit System --- #

    # components in the sellar problem
    components_sellar = (mda, obj, c1, c2)

    # inputs to the sellar problem
    inputs_sellar = (x, z1, z2)

    # outputs from the sellar problem
    outputs_sellar = (f, g1, g2)

    sellar = ExplicitSystem(components_sellar, inputs_sellar, outputs_sellar)

    # --- Query Results --- #

    # input arguments to the Sellar problem, expressed as a single vector
    X = [0.29, 0.78, 0.60]

    # outputs from the Sellar problem, expressed as a single vector
    Y = outputs!(sellar, X)

    # jacobian of the outputs with respect to the inputs, expressed as a matrix
    dYdX = jacobian!(sellar, X)

    # combined evaluation of outputs and jacobian
    Y, dYdX = outputs_and_jacobian!(sellar, X)

    # --- Test Results against AD --- #

    # # Verify using forward mode automatic differentiation
    # dYdX_ad = ForwardDiff.jacobian((X) -> outputs(sellar, X), X)
    dYdX_ad = [1.4486568466809213 2.0897560103596686 0.6033081762203991; -0.9099208777535168 -1.2374923948590968 -0.7279367033191703; 0.4503956112356196 1.6125380257032502 1.3603164834112624]

    @test isapprox(dYdX, dYdX_ad, atol=1e-6)
end


text:
We might want to document in Bytes type that it always applies big endian approach. If you want little endian, do bswap.
doesn't this solve the problem? Or do you think the transformation approach will yield better results?

ext: jl


text:
function calc_discontinuity(p, v, d, t, annual_year, r)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    # global variables only have to be done once. [r goes from 1 to 8]
    if isequal(r, 1)

        v.idis_lossfromdisc_ann[yr] = max(0, p.rt_g_globaltemperature_ann[yr] - p.tdis_tolerabilitydisc) # global temperature minus 1.5 (default)

        if is_first(t)
            if v.idis_lossfromdisc_ann[yr]*(p.pdis_probability/100) > p.rand_discontinuity # (temp - 1.5)*0.2 >? .5 (default)
                v.occurdis_occurrencedummy_ann[yr] = 1
            else
                v.occurdis_occurrencedummy_ann[yr] = 0
            end
            v.expfdis_discdecay_ann[yr]=exp(-(p.y_year[t] - p.y_year_0)/p.distau_discontinuityexponent) # left same as timestep version
        else
            if v.idis_lossfromdisc_ann[yr]*(p.pdis_probability/100) > p.rand_discontinuity
                v.occurdis_occurrencedummy_ann[yr] = 1
            elseif v.occurdis_occurrencedummy_ann[yr-1] == 1
                v.occurdis_occurrencedummy_ann[yr] = 1
            else
                v.occurdis_occurrencedummy_ann[yr] = 0
            end
            v.expfdis_discdecay_ann[yr]=exp(-(p.y_year[t] - p.y_year[t-1])/p.distau_discontinuityexponent) # left same as timestep version
        end
    end

    v.igdpeqdis_eqdiscimpact_ann[yr,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP_ann[yr,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

    if is_first(t)
        v.igdp_realizeddiscimpact_ann[yr,r]=v.occurdis_occurrencedummy_ann[yr]*(1-v.expfdis_discdecay_ann[yr])*v.igdpeqdis_eqdiscimpact_ann[yr,r]
    else
        v.igdp_realizeddiscimpact_ann[yr,r]=v.igdp_realizeddiscimpact_ann[yr-1,r]+v.occurdis_occurrencedummy_ann[yr]*(1-v.expfdis_discdecay_ann[yr])*(v.igdpeqdis_eqdiscimpact_ann[yr,r]-v.igdp_realizeddiscimpact_ann[yr-1,r])
    end

    if v.igdp_realizeddiscimpact_ann[yr,r] < p.isatg_saturationmodification    #### analyse this in dissertation?
        v.isat_satdiscimpact_ann[yr,r] = v.igdp_realizeddiscimpact_ann[yr,r]
    else
        v.isat_satdiscimpact_ann[yr,r] = p.isatg_saturationmodification + (100-p.isatg_saturationmodification)*((v.igdp_realizeddiscimpact_ann[yr,r]-p.isatg_saturationmodification)/((100-p.isatg_saturationmodification)+(v.igdp_realizeddiscimpact_ann[yr,r] - p.isatg_saturationmodification)))
    end

    v.isat_per_cap_DiscImpactperCapinclSaturation_ann[yr,r] = (v.isat_satdiscimpact_ann[yr,r]/100)*p.rgdp_per_cap_NonMarketRemainGDP_ann[yr,r]
    v.rcons_per_cap_DiscRemainConsumption_ann[yr,r] = p.rcons_per_cap_NonMarketRemainConsumption_ann[yr,r] - v.isat_per_cap_DiscImpactperCapinclSaturation_ann[yr,r]

    # # for analysis / debugging / code analysis
    # if isequal(r, 8)
    #     v.rcons_per_cap_DiscRemainConsumption_ann_sum =  sum(v.rcons_per_cap_DiscRemainConsumption_ann[yr,:])
    # end

end

@defcomp Discontinuity begin

  region = Index()
  year = Index() # NEW - for annual version

  y_year=Parameter(index=[time], unit="year")
  y_year_0 = Parameter(unit="year")
  y_year_ann = Parameter(index=[year], unit="year") # NEW - for annual version
  yagg_periodspan = Parameter(index=[time], unit="year") # added for doing in-component summation

  rand_discontinuity = Parameter(unit="unitless", default=.5)

  irefeqdis_eqdiscimpact=Variable(index=[region], unit="%")
  wincf_weightsfactor_sea = Parameter(index=[region], unit="")
  wdis_gdplostdisc=Parameter(unit="%", default=3.)

  igdpeqdis_eqdiscimpact=Variable(index=[time,region], unit="%")
  rgdp_per_cap_NonMarketRemainGDP=Parameter(index=[time,region], unit="\$/person")
  GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="unitless", default=34298.93698672955)
  ipow_incomeexponent=Parameter(unit="unitless", default=-0.13333333333333333)
    # origin in this component, no need for in-component interpolation, as they are calculated in run_timestep
  igdpeqdis_eqdiscimpact_ann=Variable(index=[year,region], unit="%")
    # from MarketDamagesBurke
  rgdp_per_cap_NonMarketRemainGDP_ann=Parameter(index=[year,region], unit="\$/person")

  igdp_realizeddiscimpact=Variable(index=[time,region], unit="%")
  occurdis_occurrencedummy=Variable(index=[time], unit="unitless")
  expfdis_discdecay=Variable(index=[time], unit="unitless")
    # origin in this component, no need for in-component interpolation, as they are calculated in run_timestep
  igdp_realizeddiscimpact_ann = Variable(index=[year,region], unit="%")
  occurdis_occurrencedummy_ann = Variable(index=[year], unit="unitless") # NEEDS review
  expfdis_discdecay_ann = Variable(index=[year], unit="unitless") # NEEDS review
  occurdis_occurrencedummy_ann_sum = Variable(unit="unitless") # for analysis / debugging / code analysis
  occurdis_occurrencedummy_sum = Variable(unit="unitless") # for analysis / debugging / code analysis

  distau_discontinuityexponent=Parameter(unit="unitless", default=20.)

  idis_lossfromdisc=Variable(index=[time], unit="degreeC")
  tdis_tolerabilitydisc=Parameter(unit="degreeC", default=1.5)
  rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
  pdis_probability=Parameter(unit="%/degreeC", default=20.)
    # origin in this component, no need for in-component interpolation, as they are calculated in run_timestep
  idis_lossfromdisc_ann = Variable(index=[year], unit="degreeC")
    # from ClimateTemperature
  rt_g_globaltemperature_ann = Parameter(index=[year], unit="degreeC")

  isatg_saturationmodification = Parameter(unit="unitless")
  isat_satdiscimpact=Variable(index=[time,region], unit="%")
    # origin in this component, no need for in-component interpolation, as they are calculated in run_timestep
  isat_satdiscimpact_ann = Variable(index=[year,region], unit="%")

  isat_per_cap_DiscImpactperCapinclSaturation=Variable(index=[time,region], unit="\$/person")
  isat_per_cap_DiscImpactperCapinclSaturation_sum=Variable(unit="\$/person") # NEW - for annual version MC
  isat_per_cap_DiscImpactperCapinclSaturation_ann_sum=Variable(unit="\$/person") # NEW - for annual version MC
  rcons_per_cap_DiscRemainConsumption=Variable(index=[time, region], unit = "\$/person")
  rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit = "\$/person")
    # origin in this component, no need for in-component interpolation, as they are calculated in run_timestep
  isat_per_cap_DiscImpactperCapinclSaturation_ann = Variable(index=[year,region], unit="\$/person")
  rcons_per_cap_DiscRemainConsumption_sum =  Variable(unit = "\$/person") # for analysis / debugging / code analysis
  rcons_per_cap_DiscRemainConsumption_ann = Variable(index=[year, region], unit = "\$/person")  # for analysis / debugging / code analysis
  rcons_per_cap_DiscRemainConsumption_ann_sum = Variable(unit = "\$/person") # for analysis / debugging / code analysis
    # from MarketDamagesBurke
  rcons_per_cap_NonMarketRemainConsumption_ann = Parameter(index=[year, region], unit = "\$/person")

    function run_timestep(p, v, d, t)

        v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - p.tdis_tolerabilitydisc)

        if is_first(t)
            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year_0)/p.distau_discontinuityexponent)
        else
            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            elseif v.occurdis_occurrencedummy[t-1] == 1
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year[t-1])/p.distau_discontinuityexponent)
        end

        for r in d.region
            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor_sea[r]*p.wdis_gdplostdisc

            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

            if is_first(t)
                v.igdp_realizeddiscimpact[t,r]=v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*v.igdpeqdis_eqdiscimpact[t,r]
            else
                v.igdp_realizeddiscimpact[t,r]=v.igdp_realizeddiscimpact[t-1,r]+v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*(v.igdpeqdis_eqdiscimpact[t,r]-v.igdp_realizeddiscimpact[t-1,r])
            end

            if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
            else
                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100-p.isatg_saturationmodification)*((v.igdp_realizeddiscimpact[t,r]-p.isatg_saturationmodification)/((100-p.isatg_saturationmodification)+(v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
            end
            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r]/100)*p.rgdp_per_cap_NonMarketRemainGDP[t,r]
            v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]

            # calculate  for this specific year
            if is_first(t)
                # v.rcons_per_cap_DiscRemainConsumption_ann_sum = 0
                for annual_year = 2015:(gettime(t))
                    calc_discontinuity(p, v, d, t, annual_year, r)
                end
                if isequal(r,8)
                    v.occurdis_occurrencedummy_sum = v.occurdis_occurrencedummy[t]*p.yagg_periodspan[t]
                    v.rcons_per_cap_DiscRemainConsumption_sum = sum(v.rcons_per_cap_DiscRemainConsumption[t,:])*p.yagg_periodspan[t]
                    v.isat_per_cap_DiscImpactperCapinclSaturation_sum = sum(v.isat_per_cap_DiscImpactperCapinclSaturation[t,:])*p.yagg_periodspan[t]
                end
            else
                for annual_year = (gettime(t-1)+1):(gettime(t))
                     calc_discontinuity(p, v, d, t, annual_year, r)

                     if isequal(annual_year, 2300)
                         if isequal(r,8)
                             v.occurdis_occurrencedummy_ann_sum = sum(v.occurdis_occurrencedummy_ann[:])
                             v.rcons_per_cap_DiscRemainConsumption_ann_sum =  sum(v.rcons_per_cap_DiscRemainConsumption_ann[:,:])
                             v.isat_per_cap_DiscImpactperCapinclSaturation_ann_sum = sum(v.isat_per_cap_DiscImpactperCapinclSaturation_ann[:,:])
                         end
                     end
                end
                if isequal(r,8)
                    # for analysis / debugging / code analysis
                    v.occurdis_occurrencedummy_sum = v.occurdis_occurrencedummy_sum + v.occurdis_occurrencedummy[t]*p.yagg_periodspan[t]
                    v.rcons_per_cap_DiscRemainConsumption_sum = v.rcons_per_cap_DiscRemainConsumption_sum + sum(v.rcons_per_cap_DiscRemainConsumption[t,:])*p.yagg_periodspan[t]
                    v.isat_per_cap_DiscImpactperCapinclSaturation_sum = v.isat_per_cap_DiscImpactperCapinclSaturation_sum + sum(v.isat_per_cap_DiscImpactperCapinclSaturation[t,:])*p.yagg_periodspan[t]
                end

            end

        end
    end
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
"""
# MAXFETCHSIZE sets the default rowset fetch size
# used in retrieving resultset blocks from queries
"""
const MAXFETCHSIZE = 65535

# success codes
const SQL_SUCCESS           = Int16(0)
const SQL_SUCCESS_WITH_INFO = Int16(1)

# error codes
const SQL_ERROR             = Int16(-1)
const SQL_INVALID_HANDLE    = Int16(-2)

# status codes
const SQL_STILL_EXECUTING   = Int16(2)
const SQL_NO_DATA           = Int16(100)

const RETURN_VALUES = Dict(SQL_ERROR   => "SQL_ERROR",
                           SQL_NO_DATA => "SQL_NO_DATA",
                           SQL_INVALID_HANDLE  => "SQL_INVALID_HANDLE",
                           SQL_STILL_EXECUTING => "SQL_STILL_EXECUTING")

# SQLAllocHandle
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712455(v=vs.85).aspx"
# Description: allocates an environment, connection, statement, or descriptor handle
# Valid handle types
const SQL_HANDLE_ENV  = Int16(1)
const SQL_HANDLE_DBC  = Int16(2)
const SQL_HANDLE_STMT = Int16(3)
const SQL_HANDLE_DESC = Int16(4)
const SQL_NULL_HANDLE = C_NULL

# SQLSetEnvAttr
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709285(v=vs.85).aspx"
# Description: sets attributes that govern aspects of environments
# Valid attributes; valid values for attribute are indented
const SQL_ATTR_CONNECTION_POOLING = 201
const SQL_CP_OFF = UInt(0)
const SQL_CP_ONE_PER_DRIVER = UInt(1)
const SQL_CP_ONE_PER_HENV = UInt(2)
const SQL_CP_DEFAULT = SQL_CP_OFF
const SQL_ATTR_CP_MATCH = 202
const SQL_CP_RELAXED_MATCH = UInt(1)
const SQL_CP_STRICT_MATCH = UInt(0)
const SQL_ATTR_ODBC_VERSION = 200
const SQL_OV_ODBC2 = 2
const SQL_OV_ODBC3 = 3
const SQL_ATTR_OUTPUT_NTS = 10001
const SQL_TRUE = 1
const SQL_FALSE = 0

# Valid attributes
const SQL_ATTR_ACCESS_MODE = 101
const SQL_MODE_READ_ONLY = UInt(1)
const SQL_MODE_READ_WRITE = UInt(0)
#const SQL_ATTR_ASYNC_DBC_EVENT
#pointer
#const SQL_ATTR_ASYNC_DBC_FUNCTIONS_ENABLE
#const SQL_ASYNC_DBC_ENABLE_ON = UInt()
#const SQL_ASYNC_DBC_ENABLE_OFF = UInt()
#const SQL_ATTR_ASYNC_DBC_PCALLBACK
#pointer
#const SQL_ATTR_ASYNC_DBC_PCONTEXT
#pointer
const SQL_ATTR_ASYNC_ENABLE = 4
const SQL_ASYNC_ENABLE_OFF = UInt(0)
const SQL_ASYNC_ENABLE_ON = UInt(1)
const SQL_ATTR_AUTOCOMMIT = 102
const SQL_AUTOCOMMIT_OFF = UInt(0)
const SQL_AUTOCOMMIT_ON = UInt(1)
const SQL_ATTR_CONNECTION_TIMEOUT = 113
#uint of how long you want the connection timeout
const SQL_ATTR_CURRENT_CATALOG = 109
#string/Ptr{UInt8} of default database to use
#const SQL_ATTR_DBC_INFO_TOKEN
#pointer
const SQL_ATTR_ENLIST_IN_DTC = 1207
#pointer: Pass a DTC OLE transaction object that specifies the transaction to export to
# SQL Server, or SQL_DTC_DONE to end the connection's DTC association.
const SQL_ATTR_LOGIN_TIMEOUT = 103
#uint of how long you want the login timeout
const SQL_ATTR_METADATA_ID = 10014
#SQL_TRUE, SQL_FALSE
const SQL_ATTR_ODBC_CURSORS = 110
const SQL_CUR_USE_IF_NEEDED = UInt(0)
const SQL_CUR_USE_ODBC = UInt(1)
const SQL_CUR_USE_DRIVER = UInt(2)
const SQL_ATTR_PACKET_SIZE = 112
#uint for network packet size
const SQL_ATTR_QUIET_MODE = 111
#window handle pointer
const SQL_ATTR_TRACE = 104
const SQL_OPT_TRACE_OFF = UInt(0)
const SQL_OPT_TRACE_ON = UInt(1)
const SQL_ATTR_TRACEFILE = 105
#A null-terminated character string containing the name of the trace file.
const SQL_ATTR_TRANSLATE_LIB = 106
# A null-terminated character string containing the name of a library containing the functions SQLDriverToDataSource and
# SQLDataSourceToDriver that the driver accesses to perform tasks such as character set translation.
const SQL_ATTR_TRANSLATE_OPTION = 107
#A 32-bit flag value that is passed to the translation DLL.
const SQL_ATTR_TXN_ISOLATION = 108
#A 32-bit bitmask that sets the transaction isolation level for the current connection.

#Valid value_length
const SQL_IS_POINTER = -4
const SQL_IS_INTEGER = -6
const SQL_IS_UINTEGER = -5
const SQL_NTS = -3

#Valid driver_prompt
const SQL_DRIVER_COMPLETE = UInt16(1)
const SQL_DRIVER_COMPLETE_REQUIRED = UInt16(3)
const SQL_DRIVER_NOPROMPT = UInt16(0)
const SQL_DRIVER_PROMPT = UInt16(2)

#valid iotype
const SQL_PARAM_INPUT = Int16(1)
const SQL_PARAM_OUTPUT = Int16(4)
const SQL_PARAM_INPUT_OUTPUT = Int16(2)
#const SQL_PARAM_INPUT_OUTPUT_STREAM = Int16()
#const SQL_PARAM_OUTPUT_STREAM = Int16()

#valid fetch_orientation
const SQL_FETCH_NEXT = Int16(1)
const SQL_FETCH_PRIOR = Int16(4)
const SQL_FETCH_FIRST = Int16(2)
const SQL_FETCH_LAST = Int16(3)
const SQL_FETCH_ABSOLUTE = Int16(5)
const SQL_FETCH_RELATIVE = Int16(6)
const SQL_FETCH_BOOKMARK = Int16(8)

const SQL_POSITION = UInt16(0) #SQLSetPos
const SQL_REFRESH = UInt16(1) #SQLSetPos
const SQL_UPDATE = UInt16(2) #SQLSetPos
const SQL_DELETE = UInt16(3) #SQLSetPos
#valid lock_type
const SQL_LOCK_NO_CHANGE = UInt16(0) #SQLSetPos
const SQL_LOCK_EXCLUSIVE = UInt16(1) #SQLSetPos
const SQL_LOCK_UNLOCK = UInt16(2) #SQLSetPos

#valid completion_type
const SQL_COMMIT = Int16(0) #SQLEndTran
const SQL_ROLLBACK = Int16(1) #SQLEndTran

#valid unique
const SQL_INDEX_ALL = UInt16(1)
const SQL_INDEX_CLUSTERED = UInt16(1)
const SQL_INDEX_HASHED = UInt16(2)
const SQL_INDEX_OTHER = UInt16(3)
const SQL_INDEX_UNIQUE = UInt16(0)
#valid reserved
const SQL_ENSURE = UInt16(1)
const SQL_QUICK = UInt16(0)

#valid id_type
const SQL_BEST_ROWID        = Int16(1) #SQLSpecialColumns
const SQL_ROWVER            = Int16(2) #SQLSpecialColumns
#valid scope
const SQL_SCOPE_CURROW      = Int16(0) #SQLSpecialColumns
const SQL_SCOPE_SESSION     = Int16(2) #SQLSpecialColumns
const SQL_SCOPE_TRANSACTION = Int16(1) #SQLSpecialColumns
#valid nullable
const SQL_NO_NULLS          = Int16(0) #SQLSpecialColumns
const SQL_NULLABLE          = Int16(1) #SQLSpecialColumns
#const SQL_NULLABLE_UNKNOWN = Int16() #SQLSpecialColumns

const SQL_ATTR_AUTO_IPD = 10001
#SQL_TRUE, SQL_FALSE
const SQL_ATTR_CONNECTION_DEAD = 1209
const SQL_CD_TRUE = 1
const SQL_CD_FALSE = 0

const SQL_ATTR_ROW_STATUS_PTR = 25
const SQL_ATTR_ROWS_FETCHED_PTR  = 26
const SQL_ATTR_ROW_ARRAY_SIZE = 27

const SQL_CLOSE = UInt16(0)
const SQL_RESET_PARAMS = UInt16(3)
const SQL_UNBIND = UInt16(2)


text:
import GeometryBasics
import Observables
import Makie

export Arrow
export update

"""
    Represents an arrow.

fields: width, color, tail and head.
"""
mutable struct Arrow <: Sprite
    width::Float64
    color::Any
    tail::Observables.Observable
    head::Observables.Observable
end

"""
    Arrow(tail, head, scene, [width, [color]])

Construct an Arrow with the given `tail`, `head`, `scene`, and the optional arguments:
`width` and `color`.
"""
function Arrow(tail::ℝ³,
               head::ℝ³,
               scene::Makie.Scene;
               width::Float64 = 0.05,
               color::Makie.RGBAf = Makie.RGBAf(1.0, 0.0, 0.0, 1.0),
               transparency::Bool = false)
    tailobservable = Observables.Observable([GeometryBasics.Point3f(vec(tail)...)])
    headobservable = Observables.Observable([GeometryBasics.Point3f(vec(head)...)])
    colorobservable = Observables.Observable(color)
    Makie.arrows!(scene,
                  tailobservable,
                  headobservable,
                  arrowsize = Makie.Vec3f(width / 2, width / 2, 2width),
                  linecolor = colorobservable,
                  arrowcolor = colorobservable,
                  linewidth = width,
                  transparency = transparency)
    Arrow(width, colorobservable, tailobservable, headobservable)
end

"""
    update(arrow, tail, head)

Update an Arrow by changing its observable with the given `arrow`, `tail` and `head`.
"""
function update(arrow::Arrow, tail::ℝ³, head::ℝ³)
    arrow.tail[] = [GeometryBasics.Point3f(vec(tail)...)]
    arrow.head[] = [GeometryBasics.Point3f(vec(head)...)];
end

"""
    update(arrow, color)

Update an Arrow by changing its observable with the given `color`.
"""
function update(arrow::Arrow, color::Makie.RGBAf)
    arrow.color[] = color
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
"""
CORDIC (COordinate Rotation DIgital Computer), a.k.a. Volder's algorithm, for computing `sin` and
`cos` functions efficiently.
"""
module CORDIC

export cordic

# Plagiarizing @njuffa's answer from https://stackoverflow.com/a/23097989/.
#
# These coefficients are obtained by Remez algorithm rather than from Taylor series.
function arctan(x::Float64)::Float64
    # argument reduction:
    #
    # arctan (-x) = -arctan(x)
    # arctan (1/x) = 1/2 * pi - arctan (x), when x > 0
    z = abs(x)
    a = (z > 1.0) ? 1.0 / z : z

    # evaluate minimax polynomial approximation
    s = a * a  # a ^ 2
    q = s * s  # a ^ 4
    o = q * q  # a ^ 8

    # use Estrin's scheme for low-order terms
    p = fma(fma(fma(-0x1.53e1d2a25ff34p-16, s, 0x1.d3b63dbb65af4p-13),
                q,
                fma(-0x1.312788dde0801p-10, s, 0x1.f9690c82492dbp-9)),
            o,
            fma(fma(-0x1.2cf5aabc7cef3p-7, s, 0x1.162b0b2a3bfcep-6),
                q,
                fma(-0x1.a7256feb6fc5cp-6, s, 0x1.171560ce4a483p-5)))

    # use Horner's scheme for high-order terms
    p = fma(fma(fma(fma(fma(fma(fma(fma(fma(fma(fma(fma(p, s,
            -0x1.4f44d841450e1p-5), s,
            +0x1.7ee3d3f36bb94p-5), s,
            -0x1.ad32ae04a9fd1p-5), s,
            +0x1.e17813d66954fp-5), s,
            -0x1.11089ca9a5bcdp-4), s,
            +0x1.3b12b2db51738p-4), s,
            -0x1.745d022f8dc5cp-4), s,
            +0x1.c71c709dfe927p-4), s,
            -0x1.2492491fa1744p-3), s,
            +0x1.99999999840d2p-3), s,
            -0x1.555555555544cp-2) * s, a, a)
    # back substitution based on argument reduction
    r = (z > 1.0) ? (0x1.921fb54442d18p+0 - p) : p
    copysign(r, x)
end

ANGLES = @. arctan(exp2(-(0:27)))
K_VALUES = cumprod(@. 1.0 / sqrt(1.0 + exp2(-(0:2:50))))

"Returns a pair of `(sin(beta), cos(beta))` using CORDIC."
function cordic(beta::Float64, n_iterations::Integer=25)::Tuple{Float64, Float64}
    if beta < -pi / 2.0
        return .-cordic(beta + pi, n_iterations)
    elseif beta > pi / 2.0
        return .-cordic(beta - pi, n_iterations)
    end

    v = [1.0, 0.0]
    angle = ANGLES[1]

    i = 0
    while i < n_iterations
        σ = sign(beta)
        if σ == 0.0
            break
        end

        x = v[1] - σ * ldexp(v[2], -i)
        y = σ * ldexp(v[1], -i) + v[2]
        v = [x, y]
        beta -= σ * angle
        i += 1

        # +1 because the iterating variable starts from 0 (2 ^ 0, 2 ^ (-1), ...).
        if i + 1 > length(ANGLES)
            angle /= 2.0
        else
            angle = ANGLES[i + 1]
        end
    end

    if i > 0
        v *= K_VALUES[min(i, length(K_VALUES))]
    end
    (v[2], v[1])
end

const flatten = Iterators.flatten
function smoke()
    X = range(-pi*4, stop=pi*4, length=800)
    println("arctan error: ", max(@. abs(arctan(X) - atan(X))...))
    println("cordic error: ", max(abs.(flatten(cordic.(X)) .- flatten(sincos.(X)))...))
end

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
module Petals

using ..Ahorn, Maple

placements = Maple.Petals

function Ahorn.canFgBg(effect::Maple.Petals)
    return true, true
end

end


text:
# Explicit 1D heat equation
using Checkpointing
using Plots
using Zygote

mutable struct Heat
    Tnext::Vector{Float64}
    Tlast::Vector{Float64}
    n::Int
    λ::Float64
    tsteps::Int
end

function advance(heat)
    next = heat.Tnext
    last = heat.Tlast
    λ = heat.λ
    n = heat.n
    for i in 2:(n-1)
        next[i] = last[i] + λ*(last[i-1]-2*last[i]+last[i+1])
    end
    return nothing
end

function sumheat(heat::Heat)
    # AD: Create shadow copy for derivatives
    shadowheat = Heat(zeros(n), zeros(n), 0, 0.0, 0)
    @checkpoint_struct revolve heat shadowheat for i in 1:tsteps
        heat.Tlast .= heat.Tnext
        advance(heat)
    end
    return reduce(+, heat.Tnext)
end

n = 100
Δx=0.1
Δt=0.001
# Select μ such that λ ≤ 0.5 for stability with μ = (λ*Δt)/Δx^2
λ = 0.5
# time steps
tsteps = 500

# Create object from struct
heat = Heat(zeros(n), zeros(n), n, λ, tsteps)

# Boundary conditions
heat.Tnext[1]   = 20.0
heat.Tnext[end] = 0

# Set up AD
# Number of available snapshots
snaps = 4
verbose = 0
revolve = Revolve(tsteps, snaps; verbose=verbose)

# Compute gradient
g = Zygote.gradient(sumheat,heat)

# Plot function values
plot(heat.Tnext)
# Plot gradient with respect to sum(T).
plot(g[1].Tnext[2:end-1])


text:
################################################################################
#VERSION 1
#linesearch is an intermediary between NewtonStopLS and a 1D solver.
function preparelinesearch(
  stp::NLPStopping,
  x₀::AbstractVector,
  d::AbstractVector,
  f₀::AbstractFloat,
  ∇f₀::AbstractVector,
  τ₀::Real,
  τ₁::Real,
)

  # initialize the LSModel
  g₀ = dot(∇f₀, d)
  ϕ = LSModel(stp.pb, x₀, d, τ₀, f₀, g₀)

  # convert the Armijo and Wolfe criteria to an asymetric interval [α,β]
  α = (τ₁ - τ₀) * g₀
  β = -(τ₁ + τ₀) * g₀

  # reuse the stopping
  ϕstp = LS_Stopping(
    ϕ,
    LSAtT(0.0, ht = f₀, gt = ∇f₀),
    optimality_check = (p, s) -> optim_check_LS(p, s, α, β),
    main_stp = stp,
    max_iter = 40,
    atol = 0.0, # to rely only on the Armijo-Wolfe conditions
    rtol = 0.0, # otherwise, tolerance may prohibit convergence
    unbounded_threshold = 1e100,
  )

  # optimize in the interval [0.0,Inf]
  ϕstp = min_1D(ϕ, 0.0, Inf, α, β, stp = ϕstp)

  # unpack the results
  if ϕstp.meta.optimal
    t, ft, gt = ϕstp.current_state.x, ϕstp.current_state.ht, ϕstp.current_state.gt
  else
    t, ft, gt = 0.0, f₀, ∇f₀
    stp.meta.fail_sub_pb = true
  end

  return t, ft, gt
end

function optim_check_LS(p, s::LSAtT, α::Float64, β::Float64)
  return max(α - s.gt, s.gt - β, 0.0)
end

################################################################################
#VERSION 2
#Not an equivalent way would be with tol_check:
function preparelinesearch(
  stp::NLPStopping,
  x₀::AbstractVector,
  d::AbstractVector,
  f₀::AbstractFloat,
  ∇f₀::AbstractVector,
  τ₀::Real,
  τ₁::Real,
)

  # initialize the LSModel
  g₀ = dot(∇f₀, d)
  ϕ = LSModel(stp.pb, x₀, d, τ₀, f₀, g₀)
  #Instead of redefining LSModel we can use the known ADNLPModel ? (not optimal though)
  #ϕ = ADNLPModel(t -> (obj(stp.pb, x₀ + t*d) - f₀ - τ₀ * t * g₀), [0.0], lvar = [0.0], uvar = [Inf])

  # convert the Armijo and Wolfe criteria to an asymetric interval [α,β]
  α = (τ₁ - τ₀) * g₀
  β = -(τ₁ + τ₀) * g₀

  # reuse the stopping
  ϕstp = NLPStopping(
    ϕ,
    NLPAtX([0.0], fx = f₀, gx = ∇f₀),
    optimality_check = unconstrained_check,
    main_stp = stp,
    max_iter = 40,
    unbounded_threshold = 1e100,
    tol_check = (a, b, c) -> β,
    tol_check_neg = (a, b, c) -> α,
  )

  # optimize in the interval [0.0,Inf]
  ϕstp = min_1D(ϕ, 0.0, Inf, stp = ϕstp)

  # unpack the results
  if ϕstp.meta.optimal
    t, ft, gt = ϕstp.current_state.x, ϕstp.current_state.ht, ϕstp.current_state.gt
  else
    t, ft, gt = 0.0, f₀, ∇f₀
    stp.meta.fail_sub_pb = true
  end

  return t, ft, gt
end


text:
module ImplicitGraphs

using DataStructures

export ImplicitGraph

import Base: getindex, show, eltype
import SimpleGraphs: deg, find_path, dist, has
export deg, find_path, dist, has

"""
`ImplicitGraph{T}(has_vertex, out_neighbors)` constructs a new `ImplicitGraph`
whose vertices have type `T`.
* `has_vertex(v::T)::Bool` is a function that checks if `v` is in the graph.
* `out_neighbors(v::T)::Vector{T}` is a function that takes an object `v` of
type `T` and returns a list of `v`'s out neighbors.
"""
struct ImplicitGraph{T}
    has_vertex::Function
    out_neighbors::Function
end

function getindex(G::ImplicitGraph{T}, v::T) where {T}
    if !has(G, v)
        error("This graph does not have a vertex $v")
    end
    G.out_neighbors(v)
end

eltype(G::ImplicitGraph{T}) where {T} = T

function getindex(G::ImplicitGraph{T}, v::T, w::T) where {T}
    if !has(G, v) || !has(G, w)
        error("One or both of vertices $v and $w are not in this graph")
    end
    return in(w, G[v])
end

has(G::ImplicitGraph{T}, v::T, w::T) where {T} = G[v, w]

"""
`deg(G::ImplicitGraph,v)` returns the degree of vertex `v`
in the graph `G`.
"""
deg(G::ImplicitGraph{T}, v::T) where {T} = length(G[v])

"""
`has(G::ImplicitGraph,v)` checks if `v` is a vertex of `G`.

`has(G,v,w)` checks if `(v,w)` is an edge of `G`.
"""
function has(G::ImplicitGraph{T}, v::T) where {T}
    return G.has_vertex(v)
end

function show(io::IO, G::ImplicitGraph{T}) where {T}
    print(io, "ImplicitGraph{$T}")
end

include("iGraphs.jl")
include("find_path.jl")
include("guided_path_finder.jl")

end # module


text:
# test for whole optimization on the space of topologies
# for h<= hmax
# Claudia March 2015
# based on test_optTopLevel.jl, but now we use a more identifiable Case G, H to begin with
# and now we use afterOptBLALL

# -------------------5taxon tree------------------

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("Tree_output.csv")
d = readTableCF(df)

# starting tree:
tree = "((6,4),(7,8),10);"
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
currT = readTopologyUpdate("prueba_tree.txt");
printEdges(currT)

epsilon = eps()
N = 100

@time newT = optTopLevel!(currT,epsilon,d,0);
#got 5.34957 at [0.2,0.1] after 28 iterations (returned FTOL_reached)
#loglik_1 = 5.34957
#found minimizer topology at step 1 with -loglik=5.34957 and ht_min=[0.2,0.1]
#elapsed time: 8.910952599 seconds (91065584 bytes allocated, 0.59% gc time)
printEdges(newT)
# forgot to copy, but true tree!

# ------------------5taxon network 1 hybridization: Case H-----------------
# starting topology: Case G

## include("../examples/case_h_example2.jl");
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseH_output2.csv",df)

include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output2.csv")
d = readTableCF(df)

@time optTopLevel!(currT,d,1)

tree = string("(1,2,((7,(8)#H5:5.814544267883624):0.9977876663423212,(#H5:1.0,4):0.0):1.9430580774498776);")
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
newT = readTopologyUpdate("prueba_tree.txt");

printEdges(newT)

newT.ht

q1 = Quartet(1,["1","4","2","7"],[0.5,0.4,0.1]);
q2 = Quartet(2,["1","4","8","7"],[0.5,0.4,0.1]);
q3 = Quartet(3,["8","4","2","7"],[0.5,0.4,0.1]);
q4 = Quartet(4,["1","8","2","7"],[0.5,0.4,0.1]);
q5 = Quartet(5,["1","4","2","8"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(newT,d)

wrongdf = writeExpCF(d.quartet)
writetable("CaseH_output_wrong_optTop_startCaseG.csv",wrongdf)

# ------------------
# starting topology: Case F
include("../examples/case_f_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output2.csv")
d = readTableCF(df)

@time optTopLevel!(currT,d,1)

tree = string("(4,6,(#H2:0.7392085405544356::0.046179825120885414,(7,(10,(8)#H2:0.0::0.9538201748791146):0.9803511144374873):2.212878358589699):0.00000073506);")
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
newT = readTopologyUpdate("prueba_tree.txt");

printEdges(newT)

#include("../examples/case_h_example2.jl");
newT.ht

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

dd = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(newT,dd)

wrongdf = writeExpCF(d.quartet)
writetable("CaseH2_output_wrong_optTop_startCaseF.csv",wrongdf)

#elapsed time: 4.236802192 seconds (274215268 bytes allocated, 4.26% gc time)
#WARNING: newT.loglik 1.4887847814495139 not really close to 0.0, you might need to redo with another starting point

# ------------------5taxon network 1 hybridization: Case F-----------------
# starting topology: Case G
include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

## include("../examples/case_f_example2.jl");
## parameters!(net)
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseF_output2.csv",df)

# real network: Case F
df = readtable("CaseF_output.csv")
df = readtable("CaseF_output2.csv")
d = readTableCF(df)

@time optTopLevel!(currT,d,1)

printEdges(newT)

# starting topology: Case H
include("../examples/case_h_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case F
df = readtable("CaseF_output.csv")
df = readtable("CaseF_output2.csv")
d = readTableCF(df)

@time optTopLevel!(currT,d,1)

printEdges(newT)


text:
t1 = time_ns()
#using BrainwaveExperimental
#using Base.Test

include("test_examples.jl")
# write your own tests here
#@test 1 == 1
t2 = time_ns()
println(float(t2-t1)*1e9)


text:
module CellQuadraturesTests

using Test
using FillArrays
using Gridap.Helpers
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Integration
using Gridap.CellData
using Gridap.TensorValues

using Gridap.Fields: MockField, MockBasis, OtherMockBasis

degree = 3
quad = Quadrature(HEX,degree)
test_quadrature(quad)
@test sum(get_weights(quad)) ≈ 1

degree = (1,2,3)
quad = Quadrature(HEX,degree)
test_quadrature(quad)
@test sum(get_weights(quad)) ≈ 1

degree = 3
quad = Quadrature(TET,degree)
test_quadrature(quad)
@test sum(get_weights(quad)) ≈ 0.5*1/3

l = 10
quad = CellQuadrature(degree,[QUAD,],Fill(1,l))
q = get_coordinates(quad)
@test isa(q,Fill)
w = get_weights(quad)
@test isa(w,Fill)

quad = CellQuadrature(degree,[QUAD,],ones(Int,l))
q = get_coordinates(quad)
@test isa(q,CompressedArray)
w = get_weights(quad)
@test isa(w,CompressedArray)

ndofs = 4

z = 2.0
v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
a = MockBasis{2}(v,ndofs)
b = MockBasis{2}(w,ndofs)
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)
g = MockField{2}(v)

l = 10
zl = [ z for  i in 1:l]
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)
gl = fill(g,l)
al = Fill(a,l)
bl = fill(b,l)

gf = GenericCellField(gl,ϕl,Val(true))
af = GenericCellField(al,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
bf = GenericCellField(bl,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
zf = convert_to_cell_field(zl,ϕl)
df = af*zf
dft = trialize_cell_basis(df)

@test sum(integrate(1,ϕl,quad)) ≈ 213.33333333333323
@test sum(integrate(zl,ϕl,quad)) ≈  426.66666666666646

cm = integrate(df⋅dft,ϕl,quad)
@test cm[1] ≈ fill(560.0,4,4)

cv = integrate(df⋅v,ϕl,quad)
@test cv[1] ≈ fill(359.9999999999998,4)

quad1 = CellQuadrature(degree,[QUAD,],Fill(1,4))
quad2 = CellQuadrature(degree,[TRI,],Fill(1,6))
quad = lazy_append(quad1,quad2)
@test isa(get_coordinates(quad),AppendedArray)
@test isa(get_weights(quad),AppendedArray)
@test isa(get_array(quad),AbstractArray{<:Quadrature})
@test isa(get_array(quad),AppendedArray)

end # module


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

using ModelBlocks
using Test
using LinearAlgebra
using Plots

variables = Variables([
    Variable("X", 1, (-100, 100), "", ""),
    Variable("Y", 1, (-100, 100), "", ""),
]);

parameters = Variables([
    Variable("x", 0.25, (-2, 2), "", ""),
    Variable("y", 0.375, (-2, 2), "", ""),
]);

reactions = [
    GeneralReaction("Sines", (dvdt, t, v, p, e) -> begin
        dvdt.X += (p.x - v.X);
        dvdt.Y += (p.y - v.Y);
    end),
];

block = Block(variables, parameters, reactions);
setTimeRange!(block, 0:10);

outputs = Variables([
    Variable("o1", 0.0, (0, 100), "?", ""),
    Variable("o2", 0.0, (0, 100), "?", ""),
]);

setOutputDefinition!(block, outputs, (variables, parameters, timeRange, solution, outputs) -> begin
        outputs.o1 = solution.X[end] ^ 2 + solution.Y[end] ^ 2;
        outputs.o2 = 0;
        return outputs;
    end);

@time ppop = generatePPopFarthest([block], [
    ("x", -2., 2.),
    ("y", -2., 2.)
], Dict(
    "o1" => (0., 1.),
    "o2" => (0., 1.)
), 100; MaxTime = 120, DistanceFactor = 1.);

Plots.plot(ppop[1,:], ppop[2,:], seriestype = :scatter)

squareSums = ones(size(ppop, 2), 1) * sum(ppop .^ 2, dims=1);
distances = squareSums + squareSums' - 2 * ppop' * ppop;
distances = reshape(distances, 1, :)[:];

Plots.histogram(log.(distances[distances .> 1e-5]), xlims=(-8, 2), ylims=(0,1000))

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
module Logs

using ProgressLogging: asprogress
using InteractiveUtils
using Dates: Dates
using Logging

import Term:
    Theme,
    textlen,
    escape_brackets,
    unescape_brackets,
    reshape_text,
    has_markup,
    int,
    highlight,
    term_theme,
    truncate,
    ltrim_str

import ..Consoles: console_width, console_height, change_scroll_region, move_to_line
import ..Renderables: AbstractRenderable
import ..Style: apply_style
import ..Tprint: tprintln
import ..Boxes: ROUNDED
import ..Progress:
    ProgressBar,
    ProgressJob,
    DescriptionColumn,
    ProgressColumn,
    render,
    start!,
    stop!,
    addjob!,
    removejob!,
    update!,
    getjob,
    SeparatorColumn,
    PercentageColumn

export TermLogger, install_term_logger

DEFAULT_LOGGER = global_logger()

# ---------------------------- create term logger ---------------------------- #
"""
    TermLogger

Custom logger type.
"""
struct TermLogger <: AbstractLogger
    io::IO
    theme::Theme
    pbar::ProgressBar
end
TermLogger(theme::Theme) = TermLogger(stderr, theme)

function TermLogger(io::IO, theme::Theme)
    return TermLogger(
        io,
        theme,
        ProgressBar(;
            transient = true,
            columns = [
                DescriptionColumn,
                SeparatorColumn,
                ProgressColumn,
                SeparatorColumn,
                PercentageColumn,
            ],
        ),
    )
end

# set logger beavior
Logging.min_enabled_level(logger::TermLogger) = Logging.Info

Logging.shouldlog(logger::TermLogger, level, _module, group, id) = true
Logging.catch_exceptions(logger::TermLogger) = true

# --------------------------- handle logger message -------------------------- #

"""
    print_closing_line(color::String, width::Int)

Print the final line of a log message with style and date info
"""
function print_closing_line(color::String, width::Int = 48)
    tprintln(
        "  {$color bold dim}$(ROUNDED.bottom.left)" *
        "$(ROUNDED.row.mid)"^(width) *
        "{/$color bold dim}",
    )
    _date = Dates.format(Dates.now(), "e, dd u yyyy")
    _time = Dates.format(Dates.now(), "HH:MM:SS")
    pad = width - textlen(_date * _time) - 2
    return tprintln(
        " "^pad * "{dim}$_date{/dim} {bold dim underline}$_time{/bold dim underline}";
        highlight = false,
    )
end

"""
    handle_progress(logger::TermLogger, prog)

Handle progress information passed by `ProgressLogging`

It creates/adds/removes `ProgressJob`s to the logger's
`ProgressBar` to create progress visualizations.
"""
function handle_progress(logger::TermLogger, prog)
    pbar = logger.pbar
    pbar.running || start!(pbar)

    # if progress not started yet, ignore
    if isnothing(prog.fraction) || prog.fraction == 0
        return nothing
    end

    # check if a job exists for this progress, or add one
    job = getjob(pbar, prog.id)
    job = if isnothing(job)
        addjob!(
            pbar;
            description = length(prog.name) > 0 ? prog.name : "running...",
            N = 100,
            transient = false,
            id = prog.id,
        )
    else
        job
    end

    # if done, remove job.
    if prog.done || prog.fraction == 1.0
        job.i = 100
        stop!(job)
    else
        update!(job; i = (Int ∘ floor)(prog.fraction * 100))
    end

    # render
    if all(map(j -> j.finished, pbar.jobs))
        map(j -> removejob!(pbar, j), pbar.jobs)
        stop!(pbar)
    else
        render(pbar)
    end
end

"""
    handle_message(logger::TermLogger, lvl, msg, _mod, group, id, file, line; kwargs...)

Handle printing of log messages, with style!.

In addition to the log message and info such as file/line and time of log,
it prints kwargs styled by their type.
"""
function Logging.handle_message(
    logger::TermLogger,
    lvl,
    msg,
    _mod,
    group,
    id,
    file,
    line;
    kwargs...,
)
    _progress = asprogress(lvl, msg, _mod, group, id, file, line; kwargs...)
    isnothing(_progress) || return handle_progress(logger, _progress)

    # get name of function where logging message was called from
    fname = ""
    for frame in stacktrace()
        if "$(frame.file)" == file && frame.line == line
            fname = ".{underline}$(frame.func){/underline}"
        end
    end

    # prepare styles
    color = if lvl == Logging.Info
        logger.theme.info
    elseif lvl == Logging.Debug
        logger.theme.debug
    elseif lvl == Logging.Warn
        logger.theme.warn
    elseif lvl == Logging.Error
        logger.theme.error
    else
        "#90CAF9"
    end

    outline_markup = "$color dim"
    vert = "{$outline_markup}" * ROUNDED.mid.left * "{/$outline_markup}"

    # style message
    msg = if msg isa AbstractString
        reshape_text(has_markup(msg) ? msg : "{#8abeff}$msg{/#8abeff}", 64)
    else
        string(msg)
    end

    # get the first line of information
    content = "{$color underline bold}@$(string(lvl)){/$color underline bold} {#edb46f }($(_mod)$fname):{/#edb46f }"

    # for multi-lines message, print each line separately.
    msg_lines::Vector{AbstractString} = split(msg, "\n")
    for n in 2:length(msg_lines)
        msg_lines[n] = "  $vert   " * " "^textlen(content) * "{#8abeff}" * msg_lines[n]
    end

    length(msg_lines) > 0 && (content *= "  " * msg_lines[1])

    tprintln(content; highlight = false)
    tprintln.(msg_lines[2:end]; highlight = false)

    # if no kwargs we can just quit
    if length(kwargs) == 0 || length(msg_lines) == 0
        print_closing_line(color)
        return nothing
    end

    # get padding width
    _types = string.(typeof.(collect(values(kwargs))))
    _types = map(t -> ltrim_str(t, 22), _types)
    for (i, _type) in enumerate(values(kwargs))
        typeof(_type) <: Function && (_types[i] = string(Function))
    end

    wpad = max(textlen.((_types))...) + 2
    ks = truncate.(string.(keys(kwargs)), 28)
    namepad = max(textlen.(ks)...)

    # print all kwargs
    tprintln("  $vert"; highlight = false)
    for (k, v, _type) in zip(ks, values(kwargs), _types)
        # get line stub
        pad = wpad - textlen(_type) - 1
        line =
            "  $vert" *
            " "^pad *
            " {$(logger.theme.type) dim}($(_type)){/$(logger.theme.type) dim}" *
            " {bold #e0e0e0}$k{/bold #e0e0e0}"

        epad = namepad - textlen(string(k))
        line *= " "^epad * " {bold red}={/bold red} "
        lpad = textlen(line) - 4

        # get value style
        if v isa Number
            _style = logger.theme.number
        elseif v isa Symbol
            _style = logger.theme.symbol
        elseif v isa AbstractString
            _style = logger.theme.string
        elseif v isa AbstractVector
            _style = logger.theme.number
            _size = length(v)
            v = escape_brackets(string(v))
            v = textlen(v) > 33 ? v[1:30] * "..." : v
            v *= "\n {white}$(_size) {/white}{dim}items{/dim}"
        elseif v isa AbstractArray || v isa AbstractMatrix
            _style = logger.theme.number
            _size = size(v)
            v = ltrim_str("$(typeof(v)) {dim}<: $(supertypes(typeof(v))[end-1]){/dim}", 30)
            v *=
                "\n {dim}shape: {default white}" *
                join(string.(_size), " × ") *
                "{/default white}{/dim}"

        elseif v isa Function
            _style = logger.theme.func
        elseif v isa AbstractRenderable
            _style = "default"
            v = "$(typeof(v)) \e[2m(size: $(v.measure))\e[0m"
        else
            _style = nothing
        end

        # print value lines
        try
            v = reshape_text(ltrim_str(string(v), 177), 60)
        catch
            v = ltrim_str(string(v), 60)
        end
        vlines = split(v, "\n")

        vlines = if !isnothing(_style)
            map(ln -> "{" * _style * "}" * ln * "{/" * _style * "}", vlines)
        else
            highlight.(vlines)
        end

        tprintln(line * vlines[1]; highlight = false)
        if length(vlines) ≥ 1
            for ln in vlines[2:end]
                tprintln("  $vert " * " "^lpad * ln; highlight = false)
            end
        end
    end
    return print_closing_line(color)
end

# ---------------------------- install term logger --------------------------- #

"""
    install_term_logger(theme::Theme=theme)

Install `TermLogger` as the global logging system.

`theme::Theme` can be passed to specify the theme to use for styling objects.
"""
function install_term_logger(theme::Theme = term_theme[])
    _logger = TermLogger(theme)
    return global_logger(_logger)
end

function uninstall_term_logger()
    _lg = global_logger(DEFAULT_LOGGER)
    logger = global_logger(DEFAULT_LOGGER)
    return logger
end

end


text:
### Control the ego car in the occluded_crosswalk environment ###
# Each policy should implement its own version of propagate corresponding to the returned action type
function AutomotiveDrivingModels.propagate(veh::Vehicle, action::Float64, roadway::Roadway, Δt::Float64)
    x_ = veh.state.posG.x + veh.state.v*Δt + action*Δt^2/2
    if x_ <= veh.state.posG.x
        x_ = veh.state.posG.x
    end
    v_ = veh.state.v + action*Δt
    if v_ <= 0.
        v_ = 0.
    elseif v_ >= 8.0 #TODO parameterized
        v_ = 8.0
    end
    return VehicleState(VecSE2(x_, veh.state.posG.y, veh.state.posG.θ), roadway, v_)
end

"""
    The driver model to be evaluated, contains the policy and the sensor model. The action is an
    acceleration updated by the observe! function at each time step
"""
mutable struct AVDriver{A, B, O, P, U, S} <: DriverModel{Float64}
    update_freq::Int64
    tick::Int64
    env::CrosswalkEnv
    a::A# acceleration
    b::B # belief
    o::O
    ego::Vehicle
    sensor::S
    policy::P # replace by policy type
    updater::U
end

AutomotiveDrivingModels.get_name(model::AVDriver) = "AVDriver"
Base.rand(model::AVDriver) = model.a # deterministic policy

function AutomotiveDrivingModels.observe!(model::AVDriver,
                                          scene::EntityFrame{VehicleState, VehicleDef, Int64},
                                          roadway::Roadway,
                                          egoid::Int)
    # @assert scene[1].id == egoid
    model.o = measure(scene, model.sensor, model.env)
    b_ = update(model.updater, model.b, model.a, model.o)
    # f = open("timing.csv", "a")
    # write(f,"belief, $btime, $(scene.n - 1) \n")
    # close(f)
    model.b = deepcopy(b_)
    if model.tick % model.update_freq == 0
        model.a = action(model.policy, model.b) # policy
        # f = open("timing.csv", "a")
        # write(f, "action, $atime, $(scene.n - 1) \n")
        # close(f)
    end
    model.tick += 1
end

"""
    initial_ego(env::CrosswalkEnv, rng::AbstractRNG)
generate an initial state for the ego car, returns a Vehicle object
"""
function initial_ego(env::CrosswalkEnv, rng::AbstractRNG)
    x0 = env.params.start_pos
    y0 = 0. # parameterize
    v0 = rand(rng, Distributions.Uniform(6., env.params.speed_limit))
    # v0 = 8.

    return Vehicle(VehicleState(VecSE2(x0, y0, 0.0), env.roadway.segments[1].lanes[1], env.roadway, v0),
                   VehicleDef(), 1)
end

function reset_model!(model::AVDriver, ego::Vehicle)
    model.tick = 0.
    model.a = 0.
    model.o = Dict{Int64, SingleOCState}()
    # model.b = initialstate(model.policy.pomdp, model.sensor.n_ped)
    b0 = Dict{Int64, SingleOCDistribution}()
    b0[-1] = initialstate(policy.pomdp, ego.state)
    model.b = b0
    model.ego = ego
    reset_policy!(model.policy)
end