text:
##############################################################################
##
## Solve hamilton-jacobi-bellman
##
##############################################################################

##############################
##
## FDSpec
##
##############################
"""
Use finite difference method to solve for HJB equation.

### INPUTS
- `fd`    : Finite Difference scheme
- `sol`   : holds solution information
- `twoap` : household problem description
"""
function solve_hjb!(twoap::TwoAssetsProb, fd::FDSpec, sol::SolutionExp;
                    maxit::Int = 400, crit::Float64 = 1e-8, monotonicity = true, verbose::Bool = true)

    bn, an, zn = length(fd.b), length(fd.a), length(fd.z)

    #== Initialize Vold ==#
    Vold = zeros(bn,an,zn)

    #== Construct functions ==#
    compute_pol!(solut, V_old) = comp_pol!(twoap, fd, solut, V_old)
    update_V!(solut, V_old)  = updateV!(twoap, fd, solut, V_old, test_mon = monotonicity)

    for it in 1:maxit
        copy!(Vold, vcat(sol.V...))
        # vcat is important in the FDExp case where fd.V is vector of vectors

        #== Compute the Optimal Policy ==#
        compute_pol!(sol, Vold)

        #== Update vⁿ ==#
        update_V!(sol, Vold)

        # check convergence
        distance = chebyshev(vcat(sol.V...), vec(Vold))
        # distance = norm(vec(fd.V)-vec(Vold), Inf)
        if distance < crit
            println("hjb solved : $(it) iterations")
            # println(out_text,"hjb solved : $(it) iterations")
            break
        else
            # update V using newV
            verbose && (it % 10 == 0) && @printf("  value function iteration %d, distance %.4f \n", it,distance)
            # verbose && (it % 10 == 0) && @printf(out_text,"  value function iteration %d, distance %.4f \n", it,distance)
        end
    end
end

function optconsump(twoap::TwoAssetsProb, ∂Vb::Float64, liqinc::Float64, ℓutil::Float64)

    invγ   = 1.0/twoap.γ

    # ∂Vb>=0 ? c = ∂Vb^(-invγ) : error("negative value not allowed")
    # REVIEW:110 trick to not get stuck in case of non-monotonicity
    ∂Vb>=0 ? c = ∂Vb^(-invγ) + ℓutil : c = liqinc
    sc = liqinc - c
    dV = utilfn(twoap, c, ℓutil) + sc*∂Vb

    return c, sc, dV
end

function comp_pol!(twoap::TwoAssetsProb, fd::FDSpec, sol::SolutionExp, V::Array{Float64,3})

    #== Construct fns ==#
    dχinv1(∂ratio, a) = dχinv(twoap, ∂ratio, a)
    χ1(d,a) = χ(twoap, d, a)
    utilfn1(c, ℓutil) = utilfn(twoap, c, ℓutil)

    #== Check if function is changing ==#
    # println(dχinv1(1.1, 1.0))
    # println(χ1(1e-5, 100.0))

    #== Unpack from FD ==#
    z = fd.z; a = fd.a ; b = fd.b;
    Δbgrid = fd.Δbgrid
    Δagrid = fd.Δagrid
    netbinc = fd.netbinc

    ℓutilgrid = fd.ℓutilgrid

    # precompute the lengths
    bn = length(b); an = length(a); zn = length(z)

    invdb = 1.0./Δbgrid
    invda = 1.0./Δagrid

    for zi in 1:zn, ai in 1:an, bi in 1:bn

        # == Derivative w.r.t B == #
        bi<bn && ( ∂Vbᶠ = ( V[bi+1, ai, zi] - V[bi, ai, zi] ) * invdb[bi] )
        bi>1  && ( ∂Vbᴮ = ( V[bi, ai, zi] - V[bi-1, ai, zi] ) * invdb[bi-1] )

        # == Derivative w.r.t A == #
        ai<an && ( ∂Vaᶠ = ( V[bi, ai+1, zi] - V[bi, ai, zi] ) * invda[ai] )
        ai>1  && ( ∂Vaᴮ = ( V[bi, ai, zi] - V[bi, ai-1, zi] ) * invda[ai-1] )

        ### CONSUMPTION decision ###
        liqinc = netbinc[bi,zi]
        ℓutil  = ℓutilgrid[zi]

        # == Foward == #
        if bi<bn
            cᶠ, scᶠ, Hcᶠ  = optconsump(twoap, ∂Vbᶠ, liqinc, ℓutil)
        else
            scᶠ, dVᶠ = 0.0, -1.0e12
        end
        scᶠ>0.0 ? (validᶠ = 1 ) : (validᶠ = 0 )

        # == Backward == #
        if bi>1
            cᴮ, scᴮ, Hcᴮ  = optconsump(twoap, ∂Vbᴮ, liqinc, ℓutil)
        else
            scᴮ, Hcᴮ = 0.0, -1.0e12
        end
        scᴮ<0.0 ? (validᴮ = 1) : (validᴮ = 0)

        #== Catch Other cases ==#
        c⁰  = liqinc
        sc⁰ = 0.0
        Hc⁰ = utilfn1(c⁰, ℓutil)

        if validᶠ==1   && ( validᴮ==0 || Hcᶠ>=Hcᴮ ) && ( Hcᶠ>=Hc⁰ )
            sol.c[bi, ai, zi]  = cᶠ
            sol.sc[bi, ai, zi] = scᶠ
        elseif validᴮ==1 && ( validᶠ==0 || Hcᴮ>=Hcᶠ ) && ( Hcᴮ>=Hc⁰ )
            sol.c[bi, ai, zi]  = cᴮ
            sol.sc[bi, ai, zi] = scᴮ
        else
            sol.c[bi, ai, zi]  = c⁰
            sol.sc[bi, ai, zi] = sc⁰
        end

        #== DEPOSIT decision ==#
        #== b backward, a forward ==#
        if bi>1 && ai<an
            dᴮᶠ = dχinv1(∂Vaᶠ/∂Vbᴮ, a[ai])
            # change in utility
            Hdᴮᶠ = ∂Vaᶠ * dᴮᶠ - ∂Vbᴮ*( dᴮᶠ + χ1(dᴮᶠ,a[ai]) )

            (dᴮᶠ>0.0 && Hdᴮᶠ>=0.0) ? ( validᴮᶠ=1 ) : ( validᴮᶠ=0 )
        else
            validᴮᶠ = 0; Hdᴮᶠ = -1.0e-12
        end

        #== b forward, a backward ==#
        if bi < bn && ai > 1
            dᶠᴮ = dχinv1(∂Vaᴮ/∂Vbᶠ, a[ai])
            # change in utility
            Hdᶠᴮ = ∂Vaᴮ * dᶠᴮ - ∂Vbᶠ*( dᶠᴮ + χ1(dᶠᴮ,a[ai]) )

            (dᶠᴮ<-χ1(dᶠᴮ,a[ai]) && Hdᶠᴮ>=0.0) ? ( validᶠᴮ=1 ) : ( validᶠᴮ=0 )
        else
            validᶠᴮ = 0; Hdᶠᴮ = -1.0e-12
        end

        #== b backward, a backward ==#
        if bi>1 && ai>1
            dᴮᴮ = dχinv1(∂Vaᴮ/∂Vbᴮ, a[ai])
            # change in utility
            Hdᴮᴮ = ∂Vaᴮ * dᴮᴮ - ∂Vbᴮ*( dᴮᴮ + χ1(dᴮᴮ,a[ai]) )

            (dᴮᴮ>-χ1(dᴮᴮ, a[ai]) && dᴮᴮ<=0 && Hdᴮᴮ>=0.0) ? ( validᴮᴮ=1 ) : ( validᴮᴮ=0 )

            # REVIEW:80 BOUNDARY Adjustment 01
            #        force the use of ∂Vbᴮ in case bi==bn
            # bi == bn && (dᴮᴮ<=0 && Hdᴮᴮ>=0.0) && (validᴮᴮ = 1)
        else
            validᴮᴮ = 0; Hdᴮᴮ = -1.0e-12
        end

        #== Check and assign deposit ==#
        if    validᴮᶠ==1                 && (validᶠᴮ==0 || Hdᴮᶠ>=Hdᶠᴮ)   && (validᴮᴮ==0 || Hdᴮᶠ>=Hdᴮᴮ)

            sol.d[bi, ai, zi] = dᴮᶠ
        # ------------------------------------------------------------------------------------------ #
        elseif (validᴮᶠ==0 || Hdᶠᴮ>=Hdᴮᶠ) &&  validᶠᴮ==1                  && (validᴮᴮ==0 || Hdᶠᴮ>=Hdᴮᴮ)

            sol.d[bi, ai, zi] = dᶠᴮ
        # ------------------------------------------------------------------------------------------ #
        elseif (validᴮᶠ==0 || Hdᴮᴮ>=Hdᴮᶠ) && (validᶠᴮ==0 || Hdᴮᴮ>=Hdᶠᴮ)   &&  validᴮᴮ==1

            sol.d[bi, ai, zi] = dᴮᴮ
        # ------------------------------------------------------------------------------------------ #
        elseif validᴮᶠ==0 && validᶠᴮ==0 && validᴮᴮ==0

            sol.d[bi, ai, zi] = 0.0
        end

        # u[bi, ai, zi] = utilfn1( c[bi, ai, zi] )
        # bdot[bi, ai, zi] = sc[bi, ai, zi] - d[bi, ai, zi] - χ1(d, a[ai])
    end

    Void
end

function updateV!(twoap::TwoAssetsProb, fd::FDExp, sol::SolutionExp, V::Array{Float64,3}; test_mon = true)

    #== Construct function ==#
    χ1(d,a) = χ(twoap, d, a)
    utilfn1(c, ℓutil) = utilfn(twoap, c, ℓutil)
    #== Idiosyncratic component ==#
    λdiag = fd.λdiag
    λtoff = fd.λoff'

    #== HOUSEHOLD Parameters ==#
    γ, ρ, _, σ, ψ = _unpackparams(twoap)

    #== Optimal policies ==#
    c = sol.c; sc = sol.sc; d = sol.d

    #== Extract GRIDS ==#
    b = fd.b; a = fd.a; z = fd.z
    Δbgrid, Δagrid = fd.Δbgrid, fd.Δagrid
    netainc = fd.netainc;
    ℓutilgrid = fd.ℓutilgrid

    invΔ = fd.invΔ

    #== SOLUTION Matrices to fill ==#
    Vnew = sol.V
    A    = sol.A
    Au   = sol.Au

    #== Storage matrices ==#
    B  = fd.B
    b̃  = fd.b̃

    # precompute the lengths
    bn = length(b); an = length(a); zn = length(z); abn = bn*an

    invdb = 1.0./Δbgrid
    invda = 1.0./Δagrid
    inv1γ = 1.0/(1.0-γ)

    for zi in 1:zn
        #== set values to zero ==#
        fill!(nonzeros(A[zi]), zero(Float64))
        fill!(nonzeros(Au[zi]), zero(Float64))

        for abi in 1:abn
        # NOTE:90 Positions of A that I need to fill up
            # A[zi][abi-bn, abi] - v(i  ,j-1,k) --> - invda*[ d⁻ ]                  (SHOULD BE POSITIVE)
            # A[zi][abi-1 , abi] - v(i+1,j  ,k) --> - invdb*[ (scᴮ)⁻ + (sdᴮ)⁻ ]     (SHOULD BE POSITIVE)
            # A[zi][abi   , abi] - v(i  ,j  ,k) --> + invdb*[ (scᴮ)⁻ + (sdᴮ)⁻ - ( (scᶠ)⁺ - (sdᶠ)⁺ ) ]
            #                                       + invda*[ d⁻ - (d⁺ + ξwz + rᵃa) ]
            # A[zi][abi+1 , abi] - v(i-1,j  ,k) --> + invdb*[ (scᶠ)⁺ + (sdᶠ)⁺ ]     (SHOULD BE POSITIVE)
            # A[zi][abi+bn, abi] - v(i  .j+1,k) --> + invda*[ d⁺ + ξwx + rᵃa ]      (SHOULD BE POSITIVE)

            ai = afromab(abi)
            bi = bfromab(abi)
            ℓutil  = ℓutilgrid[zi]

            #== RHS of the Bellman equation==#
            b̃[zi][abi] = utilfn1(c[bi, ai, zi], ℓutil) + V[bi, ai, zi] * invΔ + dot(λtoff[:,zi] ,V[bi,ai,:][:] )

            ## COMPUTE THE ENTRIES ##
            # ==================================================================
            #== bdrift ==#
            bdriftᶠ = max( sc[bi, ai, zi], 0.0 ) + max( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 )
            bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) + min( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 )

            bi <  bn && ( budriftᶠ = max( sc[bi, ai, zi] - d[bi, ai, zi] - χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
            bi <  bn && ( budriftᴮ = min( sc[bi, ai, zi] - d[bi, ai, zi] - χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
            bi == bn && ( budriftᶠ = max( sc[bi, ai, zi] - d[bi-1, ai, zi] - χ1(d[bi-1, ai, zi], a[ai]) , 0.0 ) )
            bi == bn && ( budriftᴮ = min( sc[bi, ai, zi] - d[bi-1, ai, zi] - χ1(d[bi-1, ai, zi], a[ai]) , 0.0 ) )

            #== adrift ==#
            adriftᶠ = max( d[bi, ai, zi], 0.0 ) + netainc[ai,zi]
            adriftᴮ = min( d[bi, ai, zi], 0.0 )

            bi <  bn && ( audriftᶠ = max( d[bi, ai, zi] + netainc[ai,zi], 0.0 ) )
            bi <  bn && ( audriftᴮ = min( d[bi, ai, zi] + netainc[ai,zi], 0.0 ) )
            bi == bn && ( audriftᶠ = max( d[bi-1, ai, zi] + netainc[ai,zi], 0.0 ) )
            bi == bn && ( audriftᴮ = min( d[bi-1, ai, zi] + netainc[ai,zi], 0.0 ) )

            #== BOUNDARY Adjustment 02 ==#
            # bi <  bn && ( bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) + min( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
            # bi == bn && ( bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) )
            #== BOUNDARY Adjustment 03 ==#
            # ai <  an && ( adriftᴮ = min( d[bi, ai, zi], 0.0 ) )
            # ai == an && ( adriftᴮ = d[bi, ai, zi] + netainc[ai,zi] )

            # v[i, j- 1, k]
            if adriftᴮ != 0.0 && ai>1
                row = abfromab(ai-1,bi)
                A[zi][row, abi] = -adriftᴮ * invda[ai-1]
            end

            # v[i, j- 1, k] KF
            if audriftᴮ != 0.0 && ai>1
                row = abfromab(ai-1,bi)
                Au[zi][row, abi] = -audriftᴮ * invda[ai-1]
            end

            # v[i-1, j, k]
            if bdriftᴮ != 0.0 && bi>1
                row = abfromab(ai,bi-1)
                A[zi][row, abi] = -bdriftᴮ * invdb[bi-1]
            end

            # v[i-1, j, k] KF
            if budriftᴮ != 0.0 && bi>1
                row = abfromab(ai,bi-1)
                Au[zi][row, abi] = -budriftᴮ * invdb[bi-1]
            end

            # v[i, j, k]
            val = 0.0
            ( bdriftᴮ != 0.0 && bi > 1 ) && ( val += bdriftᴮ * invdb[bi-1] )
            ( bdriftᶠ != 0.0 && bi < bn) && ( val -= bdriftᶠ * invdb[bi]   )
            ( adriftᴮ != 0.0 && ai > 1)  && ( val += adriftᴮ * invda[ai-1] )
            ( adriftᶠ != 0.0 && ai < an) && ( val -= adriftᶠ * invda[ai]   )
            val != 0.0 && ( A[zi][abi, abi] = val )

            # v[i, j, k] KF
            val = 0.0
            ( budriftᴮ != 0.0 && bi > 1 ) && ( val += budriftᴮ * invdb[bi-1] )
            ( budriftᶠ != 0.0 && bi < bn) && ( val -= budriftᶠ * invdb[bi]   )
            ( audriftᴮ != 0.0 && ai > 1)  && ( val += audriftᴮ * invda[ai-1] )
            ( audriftᶠ != 0.0 && ai < an) && ( val -= audriftᶠ * invda[ai]   )
            val != 0.0 && ( Au[zi][abi, abi] = val )

            # v[i-1, j, k]
            if bdriftᶠ != 0.0 && bi<bn
                row = abfromab(ai,bi+1)
                A[zi][row, abi] = bdriftᶠ * invdb[bi]
            end

            # v[i-1, j, k] KF
            if budriftᶠ != 0.0 && bi<bn
                row = abfromab(ai,bi+1)
                Au[zi][row, abi] = budriftᶠ * invdb[bi]
            end

            # v[i, j+1, k]
            if adriftᶠ != 0.0 && ai<an
                row = abfromab(ai+1,bi)
                A[zi][row, abi] = adriftᶠ * invda[ai]
            end

            # v[i, j+1, k]
            if audriftᶠ != 0.0 && ai<an
                row = abfromab(ai+1,bi)
                Au[zi][row, abi] = audriftᶠ * invda[ai]
            end

        end
    end

    ## Set B = diag(invΔ + ρ) - A
    # TODO: think of a better way to do this part
    @inbounds for zi in 1:zn
        Avals = nonzeros(A[zi])
        Bvals = nonzeros(B[zi])
        #== clean the sotorage matrix before assigning values ==#
        fill!(Bvals, zero(Float64))
        Brows = rowvals(B[zi])
        for abi in 1:abn
            for k in nzrange(B[zi], abi)
                # loop over elements in the column ij
                row = Brows[k]
                Bvals[k] = - Avals[k]
                row == abi && (Bvals[k] += invΔ + ρ - λdiag[zi])
                #NOTE:120 `λdiag` is entering only in this ste
            end
        end
        Vnew[zi][:] = B[zi]' \ b̃[zi]
    end

    #== Check for monotonicity ==#
    if test_mon
        for zi in 1:zn, abi in 1:abn
            ai = afromab(abi)
            bi = bfromab(abi)
            #= Check mon in the b dimension =#
            if bi>1
                Vnew[zi][abi]<Vnew[zi][abfromab(ai,bi-1)] &&
                @printf("   non-mon on b dimension at = (%d,%d)  \n", bi,ai)
            end

            #= Check mon in the a dimension =#
            if ai>1
                Vnew[zi][abi]<Vnew[zi][abfromab(ai-1,bi)] &&
                @printf("   non-mon on a dimension at = (%d,%d)  \n", bi,ai)
            end
        end
    end
    Void
end

"""
Update value function based on Fortran code
"""
function updateV!(fd::FDImp, sol::SolutionImp, twoap::TwoAssetsProb, V::Array{Float64,3}; test_mon = true)

    #== HOUSEHOLD Parameters ==#
    γ, ρ, ξ, _, _ = _unpackparams(twoap)

    #== Optimal policies ==#
    c = sol.c; sc = sol.sc; d = sol.d

    #== Extract grids and matrices ==#
    b = fd.b; a = fd.a; z = fd.z
    Δbgrid, Δagrid = fd.Δbgrid, fd.Δagrid
    netainc = fd.netainc;
    ℓutilgrid = fd.ℓutilgrid
    invΔ = fd.invΔ

    #== Solution matrices to fill ==#
    V  = sol.V
    A  = sol.A
    Au = sol.Au

    #== Storage matrices to fill ==#
    b̃ = fd.b̃
    Λ = fd.Λ
    B = fd.B

    # precompute the lengths
    bn = length(b); an = length(a); zn = length(z); abn = bn*an

    invdb = 1.0./Δbgrid
    invda = 1.0./Δagrid
    inv1γ = 1.0/(1.0-γ)

    ijk = zero(Int)

    #== set past values to zero ==#
    fill!(nonzeros(Au), zero(Float64))
    fill!(nonzeros(A), zero(Float64))

    for zi in 1:zn, abi in 1:abn
    # NOTE:190 Positions of A that I need to fill
        # A[ijk-bn, ijk] - v(i  ,j-1,k) --> - invda*[ d⁻ ]
        # A[ijk- 1, ijk] - v(i+1,j  ,k) --> - invdb*[ (scᴮ)⁻ + (sdᴮ)⁻ ]
        # A[ijk   , ijk] - v(i  ,j  ,k) -->   invdb*[ (scᴮ)⁻ + (sdᴮ)⁻ -(scᶠ)⁺ -(sdᶠ)⁺] + invda*[ d⁻ -(d⁺ + ξwz + rᵃa) ]
        # A[ijk+ 1, ijk] - v(i-1,j  ,k) -->   invdb*[ (scᶠ)⁺ + (sdᶠ)⁺ ]
        # A[ijk+bn, ijk] - v(i  .j+1,k) -->   invda*[ d⁺ + ξwx + rᵃa ]

        ijk +=1

        ai = afromab(abi)
        bi = bfromab(abi)
        ℓutil  = ℓutilgrid[zi]

        #== RHS of the Bellman equation==#
        b̃[ijk] =  utilfn1(c[bi, ai, zi], ℓutil)  + V[bi, ai, zi] * invΔ
        # ==================================================================

        ### COMPUTE THE ENTRIES ###

        #== bdrift ==#
        bdriftᶠ = max( sc[bi, ai, zi], 0.0 ) + max( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 )
        bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) + min( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 )
        #== bdrift for KF ==#
        bi <  bn && ( budriftᶠ = max( sc[bi, ai, zi] - d[bi, ai, zi] - χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
        bi <  bn && ( budriftᴮ = min( sc[bi, ai, zi] - d[bi, ai, zi] - χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
        bi == bn && ( budriftᶠ = max( sc[bi, ai, zi] - d[bi-1, ai, zi] - χ1(d[bi-1, ai, zi], a[ai]) , 0.0 ) )
        bi == bn && ( budriftᴮ = min( sc[bi, ai, zi] - d[bi-1, ai, zi] - χ1(d[bi-1, ai, zi], a[ai]) , 0.0 ) )

        #== adrift ==#
        adriftᶠ = max( d[bi, ai, zi], 0.0 ) + netainc[ai,zi]
        adriftᴮ = min( d[bi, ai, zi], 0.0 )
        # for KF
        bi <  bn && ( audriftᶠ = max( d[bi, ai, zi] + netainc[ai,zi], 0.0 ) )
        bi <  bn && ( audriftᴮ = min( d[bi, ai, zi] + netainc[ai,zi], 0.0 ) )
        bi == bn && ( audriftᶠ = max( d[bi-1, ai, zi] + netainc[ai,zi], 0.0 ) )
        bi == bn && ( audriftᴮ = min( d[bi-1, ai, zi] + netainc[ai,zi], 0.0 ) )

        #== REVIEW:100 BOUNDARY ADJUSTMENT 02 ==#
        # bi <  bn && ( bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) + min( -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) , 0.0 ) )
        # bi == bn && ( bdriftᴮ = min( sc[bi, ai, zi], 0.0 ) -d[bi, ai, zi] -χ1(d[bi, ai, zi], a[ai]) )
        #== REVIEW:100 BOUNDARY ADJUSTMENT 03 ==#
        # ai <  an && ( adriftᴮ = min( d[bi, ai, zi], 0.0 ) )
        # ai == an && ( adriftᴮ = d[bi, ai, zi] + netainc[ai,zi] )

        # v[i, j- 1, k]
        if adriftᴮ != 0.0 && ai>1
            row = abzfromabz(ai-1,bi,zi)
            A[row, ijk] = -adriftᴮ * invda[ai-1]
        end

        # v[i, j- 1, k] KF
        if audriftᴮ != 0.0 && ai>1
            row = abzfromabz(ai-1,bi,zi)
            Au[row, ijk] = -audriftᴮ * invda[ai-1]
        end

        # v[i-1, j, k]
        if bdriftᴮ != 0.0 && bi>1
            row = abzfromabz(ai,bi-1,zi)
            A[row,ijk] = -bdriftᴮ*invdb[bi-1]
        end

        # v[i-1, j, k] KF
        if budriftᴮ != 0.0 && bi>1
            row = abzfromabz(ai,bi-1,zi)
            Au[row,ijk] = -budriftᴮ*invdb[bi-1]
        end

        # v[i, j, k]
        val = 0.0
        ( bdriftᴮ != 0.0 && bi > 1 ) && ( val += bdriftᴮ * invdb[bi-1] )
        ( bdriftᶠ != 0.0 && bi < bn) && ( val -= bdriftᶠ * invdb[bi]   )
        ( adriftᴮ != 0.0 && ai > 1)  && ( val += adriftᴮ * invda[ai-1] )
        ( adriftᶠ != 0.0 && ai < an) && ( val -= adriftᶠ * invda[ai]   )
        val != 0.0 && ( A[ijk, ijk] = val )

        # v[i, j, k] KF
        val = 0.0
        ( budriftᴮ != 0.0 && bi > 1 ) && ( val += budriftᴮ * invdb[bi-1] )
        ( budriftᶠ != 0.0 && bi < bn) && ( val -= budriftᶠ * invdb[bi]   )
        ( audriftᴮ != 0.0 && ai > 1)  && ( val += audriftᴮ * invda[ai-1] )
        ( audriftᶠ != 0.0 && ai < an) && ( val -= audriftᶠ * invda[ai]   )
        val != 0.0 && ( Au[ijk, ijk] = val )

        # v[i-1, j, k]
        if bdriftᶠ != 0.0 && bi<bn
            row = abzfromabz(ai,bi+1,zi)
            A[row, ijk] = bdriftᶠ*invdb[bi]
        end

        # v[i-1, j, k] KF
        if budriftᶠ != 0.0 && bi<bn
            row = abzfromabz(ai,bi+1,zi)
            Au[row, ijk] = budriftᶠ*invdb[bi]
        end

        # v[i, j+1, k]
        if adriftᶠ != 0.0 && ai<an
            row = abzfromabz(ai+1,bi,zi)
            A[row, ijk] = adriftᶠ * invda[ai]
        end

        # v[i, j+1, k] KF
        if audriftᶠ != 0.0 && ai<an
            row = abzfromabz(ai+1,bi,zi)
            Au[row, ijk] = audriftᶠ * invda[ai]
        end

    end
    #== set A = A + Λ ==#
    # NOTE:70 Notice that all matrices are initialized transposed
    #        they are going to be transposed at the END
    Λvals = nonzeros(Λ)
    Avals = nonzeros(A)
    copy!(Avals, Avals + Λvals)

    #== Set B = diag(invΔ + ρ) - A ==#
    # TODO: think of a better way to do this part
    Bvals = nonzeros(B)
    fill!(Bvals, zero(Float64))
    Brows = rowvals(B)
    ijk = zero(Int)
    @inbounds for zi in 1:zn, ai in 1:an, bi in 1:bn
        ijk += 1
        for k in nzrange(B, ijk)
            # loop over elements in the column ijk
            row = Brows[k]
            Bvals[k] = - Avals[k]
            if row == ijk
                Bvals[k] += invΔ + ρ
            end
        end
    end

    # solve for new V (NOTE:280 the B transposed)
    V[:] = B' \ b̃

    return Void
end

##############################
##
## TwoAssetsFD
##
# soldated
##############################
# """
# Use finite difference method to solve for HJB equation.
# """
# function solve_hjb!(fd::TwoAssetsFD;
#                     maxit::Int = 100,
#                     crit::Float64 = 1e-5,
#                     verbose::Bool = true)
#
#     for iter in 1:maxit
#
#         ## UPDATE newV
#         updateV!(fd)
#
#         # check convergence
#         # distance = chebyshev(vec(fd.newV), vec(fd.V))
#         distance = norm(vec(fd.newV)-vec(fd.V),Inf)
#         if distance < crit
#             if verbose
#                 println("hjb solved : $(iter) iterations")
#             end
#             break
#         else
#             # update V using newV
#             fd.V = fd.newV
#             @printf("Value function iteration %d, distance %.4f \n", iter,distance)
#         end
#     end
# end
#
# """
#
# Update value function based on MATLAB code
# """
# function updateV!(fd::TwoAssetsFD)
#
#     #== Parameters ==#
#     γ = fd.twoap.γ; ρ = fd.twoap.ρ; ξ = fd.twoap.ξ; rᴬ = fd.twoap.rᴬ;
#     rᴮ⁺ = fd.twoap.rᴮ; w = fd.twoap.w; χ₀ = fd.twoap.χ₀; χ₁ = fd.twoap.χ₁
#     rᴮ⁻ = fd.twoap.rᴮ + fd.twoap.wedge
#
#     #== FD method ==#
#     invΔ = fd.invΔ;
#     z = fd.z; a = fd.a ; b = fd.b
#     Δb = fd.Δb; Δa = fd.Δa
#     u = fd.u;
#     B = fd.B; A = fd.A; Λ = fd.Λ; V = fd.V
#
#     # precompute the lengths
#     bn = length(b)
#     an = length(a)
#     zn = length(z)
#     N = an*bn*zn
#
#     invdb = 1.0/(Δb)
#     invda = 1.0/(Δa)
#     invγ  = 1/γ
#     inv1γ = 1.0/(1.0-γ)
#
#     V = reshape(V, bn, an, zn)
#
#     # set A = Λ
#     # NOTE:20 Notice that all matrices are initialized transposed
#     #        they are going to be transposed at the END
#     Λvals = nonzeros(Λ)
#     Avals = nonzeros(A)
#     copy!(Avals, Λvals)
#
#     # update A and u
#     ij = zero(Int)
#     @inbounds for zi in 1:zn, ai in 1:an, bi in 1:bn
#         ij += 1
#
#         rᴮ = (b[bi]>0)*rᴮ⁺ + ( 1-(b[bi]>0) )*rᴮ⁻ # set the interest rate
#
#         ## == Derivative w.r.t B == ##
#         bi<bn ? ( ∂Vbᶠ = ( V[bi+1, ai, zi] - V[bi, ai, zi] ) * invdb ) : ( ∂Vbᶠ = ( ( 1-ξ ) * w * z[zi] + rᴮ⁺ * b[end] )^(-γ) )
#         bi>1  ? ( ∂Vbᴮ = ( V[bi, ai, zi] - V[bi-1, ai, zi] ) * invdb ) : ( ∂Vbᴮ = ( ( 1-ξ ) * w * z[zi] + rᴮ⁻ * b[1]   )^(-γ) )
#
#         ## == Derivative w.r.t A == ##
#         ai<an ? ( ∂Vaᶠ = ( V[bi, ai+1, zi] - V[bi, ai, zi] ) * invda ) : ∂Vaᶠ = zero(Float64) # make sure below that don't use ∂Vaᶠ
#         ai>1  ? ( ∂Vaᴮ = ( V[bi, ai, zi] - V[bi, ai-1, zi] ) * invda ) : ∂Vaᴮ = zero(Float64)
#
#         ## CONSUMPTION decision ##
#         cᶠ = ∂Vbᶠ^(-invγ); cᴮ = ∂Vbᴮ^(-invγ)
#         scᴮ = (1-ξ) * z[zi] * w + b[bi] * rᴮ - cᴮ
#         scᶠ = (1-ξ) * z[zi] * w + b[bi] * rᴮ - cᶠ
#
#         ## DEPOSIT decision ##
#         dᴮᴮ = foc_dep( ∂Vbᴮ, ∂Vaᴮ, a[ai] )
#         dᴮᶠ = foc_dep( ∂Vbᴮ, ∂Vaᶠ, a[ai] )
#         dᶠᴮ = foc_dep( ∂Vbᶠ, ∂Vaᴮ, a[ai] )
#         dᶠᶠ = foc_dep( ∂Vbᶠ, ∂Vaᶠ, a[ai] )
#
#     # DOUBT: having a hard time to understand these constraints...
#         #    Although the deposits are computed for bi=1,bn, later on
#         #    these are set to zero when doing upwind for b
#
#         #== eq (13),(14) + set the boundaries in a ==#
#         dᴮ = min( dᴮᴮ, 0.0 ) + max(dᴮᶠ, 0.0)     # eq (13)/ note dᴮᴮ > dᴮᶠ
#         ai==1  && ( dᴮ = (dᴮᶠ>0)*dᴮᶠ )              # Moll: make sure d>=0 at amin, don't use VaB(:,1,:)
#         ai==an && ( dᴮ = (dᴮᴮ<0)*dᴮᴮ )              # Moll: make sure d<=0 at amax, don't use VaF(:,J,:)
#
#         dᶠ = min( dᶠᴮ,0.0 ) + max(dᶠᶠ, 0.0)         # eq (14)
#         ai==1  && ( dᶠ = (dᶠᶠ>0)*dᶠᶠ )
#         ai==an && ( dᶠ = (dᶠᴮ<0)*dᶠᴮ )
#
#         sdᴮ = - dᴮ - χ(dᴮ,a[ai])
#         sdᶠ = - dᶠ - χ(dᶠ,a[ai])
#         bi==bn && ( sdᶠ = min(sdᶠ,0.0) )
#
#         # == consumption UPWIND SCHEME == #
#         c₀ = (1-ξ) * z[zi] * w + b[bi] * rᴮ
#         Icᴮ = (scᴮ<0) ; Icᶠ = (scᶠ>0 && scᴮ>=0); Ic⁰ = 1 - Icᴮ - Icᶠ
#
#         c  = Icᴮ*cᴮ + Icᶠ*cᶠ + Ic⁰*c₀
#         u[ij] = inv1γ * c^(1-γ) + invΔ * V[ij]  # notice that V[ij] still valid
#
#         # == dep UPWIND SCHEME for Vb == #
#         Idᶠ = (sdᶠ>0); Idᴮ = (sdᴮ<0 && sdᶠ<=0)
#
#         # NOTE:30 Boundary constraints of b
#         bi==1  && ( Idᴮ =0  )               # make sure don't use ∂Vbᴮ  if ib=1  for deposit decision
#         bi==bn && ((Idᴮ,Idᶠ)=(1,0))         # make sure don't use ∂Vbᶠ  if ib=bn for deposit decision
#                                             # and FORCE ∂Vbᴮ is used at bn for deposit decision
#         Id⁰ = 1 - Idᶠ - Idᴮ                 # INACTION Region
#
#         #== dep UPWINDE SCHEMEq for Va- eq (17) ==#
#         d⁻ = min(Idᴮ * dᴮᴮ + Idᶠ * dᶠᴮ,0.0)       # use backward for ∂Va
#         d⁺ = max(Idᴮ * dᴮᶠ + Idᶠ * dᶠᶠ,0.0)       # use forward for ∂Va
#         mf = d⁺ + ξ*w*z[zi] + rᴬ*a[ai]
#
#         # NOTE:40 Boundary constraints of a
#         #   Aggregate all terms in ∂Vaᴮ at an
#         ai==an && (mf = 0; d⁻ = ξ*w*z[zi] + rᴬ*a[an] + d⁻)
#
#         # NOTE:50 Positions of A that I need to fill
#             # A[ij-bn, ij] - v(i  ,j-1,k) --> - invda*[ d⁻ ]
#             # A[ij- 1, ij] - v(i+1,j  ,k) --> - invdb*[(scᴮ)⁻ + (sdᴮ)⁻ ]
#             # A[ij   , ij] - v(i  ,j  ,k) --> invdb*[ (scᴮ)⁻ - (scᶠ)⁺ + (sdᴮ)⁻ - (sdᶠ)⁺ ] + invda*[ d⁻ - (d⁺ + ξwz + rᵃa) ]
#             # A[ij+ 1, ij] - v(i-1,j  ,k) --> invdb*[ (scᶠ)⁺ + (sdᶠ)⁺ ]
#             # A[ij+bn, ij] - v(i  .j+1,k) --> invda*[ d⁺ + ξwx + rᵃa ]
#             #
#             # the matrix a is changing its entrances... Note that code from aigary only changed when
#             # it was nonzero...
#
#         current = Array(Float64,0)
#         zi == 2   && push!(current, 0.0)
#
#         ij>bn     && ( ai>1  ? push!(current, - invda * d⁻ ) : push!(current, 0.0) )
#         ij>1      && ( bi>1  ? push!(current, - invdb *( Icᴮ*scᴮ + Idᴮ*sdᴮ ) ) : push!(current, 0.0) )
#         push!(current, invdb *( Icᴮ*scᴮ + Idᴮ*sdᴮ - Icᶠ*scᶠ - Idᶠ*sdᶠ ) + invda *( d⁻ - mf ) )
#         ij<N      && ( bi<bn ? push!(current, invdb *( Icᶠ*scᶠ + Idᶠ*sdᶠ ) ) : push!(current,0.0) )
#         ij<=N-bn  && ( ai<an ? push!(current, invda * mf ) : push!(current,0.0) )
#
#         zi == 1   && push!(current, 0.0)
#
#         Avals[nzrange(A, ij)] += current
#     end
#
#     Bvals = nonzeros(B)
#     Brows = rowvals(B)
#     ij = zero(Int)
#     @inbounds for zi in 1:zn, ai in 1:an, bi in 1:bn
#         ij += 1
#         for k in nzrange(B, ij)
#             # loop over elements in the column ij
#             row = Brows[k]
#             Bvals[k] = - Avals[k]
#             if row == ij
#                 Bvals[k] += invΔ + ρ
#             end
#         end
#     end
#
#     #== Solve for new V ==#
#     #NOTE:60 the B transposed
#     fd.newV = B' \ u
#
#     return Void
# end


text:
module TestSelect

using Relations
using Base.Test

a = [1, 2, 3]
b = [4, 5, 6]
c = ["a", "b", "c"]
d = [:a, :b, :c]

for T in (Array, NullableArray)
    global src = Relation(a=T(a), b=T(b), c=T(c), d=T(d))
    global _src = copy(src)

    for (field, col) in zip(Relations.attributes(src), Relations.columns(src))
        @eval res = @with src select($field)
        @eval _res = Relation($field = $T($col))
        @test isequal(src, _src)
        @test isequal(res, _res)
    end

    f = a .* b
    _res1 = Relation(f=T(f))
    res1 = @with src select(f = a * b)
    @test isequal(src, _src)
    @test isequal(res1, _res1)
end

# non-standard lifting semantics

# isnull(x)
src = Relation(
    a = NullableArray(collect(1:5), [true, false, true, false, true]))
_res = Relation(
    b = [true, false, true, false, true])

res = @with src select(b = isnull(a))
@test isequal(res, _res)

end

ext: jl


text:
using Documenter, Fermi

makedocs(
    sitename="Fermi.jl",
    authors = "Gustavo Aroeira",
    format = Documenter.HTML(
        sidebar_sitename = false
    ),
    pages = [
        "Home" => "index.md",
        "Modules" => "modules.md",
        "Core" => "core.md",
        "Backend" => "backend.md",
        "Methods" => "methods.md",
        "Contributing" => "contributing.md",
        "Index" => "indice.md"
    ]
)

deploydocs(
    repo = "github.com/FermiQC/Fermi.jl.git",
)


text:
module Adwin

include("./Adwin2.jl")
include("./AdwinNaive.jl")

end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
## Banded LU decomposition

# this is currently a hack and not the nicest of the interfaces.
# if the LAPACK storage was built in the BandedMatrix we could get rid of
# the last three fields of this type, and have the first point to the parent
# BandedMatrix, rather then to a generic chunk of memory
struct BandedLU{T} <: Factorization{T}
    data::Matrix{T}        # banded matrix plus additional storage
    ipiv::Vector{BlasInt}  # vector of pivots
    l::Int                 # lower bandwidth
    u::Int                 # upper bandwidth
    m::Int                 # number of rows
end

# conversion
convert(::Type{BandedLU{T}}, B::BandedLU{S}) where {T<:Number, S<:Number} =
    BandedLU{T}(convert(Matrix{T}, B.data), B.ipiv, B.l, B.u, B.m)

# size of the parent array
size(A::BandedLU) = (A.m, size(A.data, 2))
size(A::BandedLU, k::Integer) = k <= 0 ? error("dimension out of range") :
                                k == 1 ? A.m :
                                k == 2 ? size(A.data, 2) : 1

# LU factorisation with pivoting. This makes a copy!
function lu(A::BandedMatrix{T}) where {T<:Number}
    # copy to a blas type that allows calculation of the factorisation in place.
    S = _promote_to_blas_type(T, T)
    # copy into larger array of size (2l+u*1)×n, i.e. l additional rows
    m, n = size(A)
    data = Array{S}(undef, 2*A.l+A.u+1, n)
    data[(A.l+1):end, :] = A.data
    data, ipiv = gbtrf!(m, A.l, A.u, data)
    BandedLU{S}(data, ipiv, A.l, A.u, m)
end
lu(F::BandedLU) = F # no op

lu(A::AbstractBandedMatrix) = lu(convert(BandedMatrix, A))
lu(A::Transpose{<:Any,<:AbstractBandedMatrix}) = lu(convert(BandedMatrix, A))
lu(A::Adjoint{<:Any,<:AbstractBandedMatrix}) = lu(convert(BandedMatrix, A))

adjoint(F::BandedLU) = Adjoint(F)
transpose(F::BandedLU) = Transpose(F)

# TODO: Finish. There was something weird about the pivoting
# function Base.getproperty(F::BandedLU{T}, d::Symbol) where T
#     if d == :L
#         m, n = size(F)
#         l,u = F.l, F.u
#         return _BandedMatrix([Ones{T}(1,n); view(F.data,l+u+2:2l+u+1, :)], m, l,0)
#     elseif d == :U
#         m, n = size(F)
#         l,u = F.l, F.u
#         return _BandedMatrix(F.data[l+1:l+u+1, :], n, 0, u)
#     elseif d == :p
#         return LinearAlgebra.ipiv2perm(getfield(F, :ipiv), m)
#     elseif d == :P
#         m, n = size(F)
#         return Matrix{T}(I, m, m)[:,invperm(F.p)]
#     else
#         getfield(F, d)
#     end
# end
#
# # iteration for destructuring into components
# Base.iterate(S::BandedLU) = (S.L, Val(:U))
# Base.iterate(S::BandedLU, ::Val{:U}) = (S.U, Val(:p))
# Base.iterate(S::BandedLU, ::Val{:p}) = (S.p, Val(:done))
# Base.iterate(S::BandedLU, ::Val{:done}) = nothing

## Utilities

checksquare(A::BandedLU) = (A.m == size(A.data, 2) ||
    throw(DimensionMismatch("matrix must be matrix is not square: dimensions are $(size(A))")))

## Conversion/Promotion

# Returns the narrowest blas type given the eltypes of A and b in A*x=b
function _promote_to_blas_type(::Type{T}, ::Type{S}) where {T<:Number, S<:Number}
    TS = Base.promote_op(/, T, S)
    # promote to narrowest type
    TS <: Complex       && return Base.promote_op(/, TS, ComplexF32)
    TS <: AbstractFloat && return Base.promote_op(/, TS, Float32)
    error("Cannot convert objects of element type $(T), $(S) to a `BlasFloat` type")
end

# Converts A and b to the narrowest blas type
for typ in (BandedMatrix, BandedLU, (Transpose{V,BandedLU{V}} where V), (Adjoint{V,BandedLU{V}} where V))
    @eval function _convert_to_blas_type(A::$typ, B::AbstractVecOrMat{T}) where {T<:Number}
        TS = _promote_to_blas_type(eltype(A), eltype(B))
        AA = convert($typ{TS}, A)
        BB = convert(Array{TS, ndims(B)}, B)
        AA, BB # one of these two might make a copy
    end
end

# basic interface
(\)(A::Union{BandedLU{T}, BandedMatrix{T}}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
    ldiv!(A, copy(B)) # makes a copy


text:
# author: Ziyi Yin
# This script trains a Fourier Neural Operator which maps Poisson coefficients to solutions
# This script re-implements the script in Zongyi Li's repository

using DrWatson
@quickactivate "FNO"
import Pkg; Pkg.instantiate()

using PyPlot
using BSON
using Flux, Random, FFTW, Zygote, NNlib
using MAT, Statistics, LinearAlgebra
using CUDA
using ProgressMeter

try
    CUDA.device()
    global gpu_flag=true
catch e
    println("CUDA.device() found no GPU device on this machine.")
    global gpu_flag=false
end

CUDA.culiteral_pow(::typeof(^), a::Complex{Float32}, b::Val{2}) = real(conj(a)*a)
CUDA.sqrt(a::Complex) = cu(sqrt(a))
Base.broadcasted(::typeof(sqrt), a::Base.Broadcast.Broadcasted) = Base.broadcast(sqrt, Base.materialize(a))

include("utils.jl")

Random.seed!(3)

mutable struct SpectralConv2d{T,N}
    weights1::AbstractArray{T,N}
    weights2::AbstractArray{T,N}
end

@Flux.functor SpectralConv2d

# Constructor
function SpectralConv2d(in_channels::Integer, out_channels::Integer, modes1::Integer, modes2::Integer)
    scale = (1f0 / (in_channels * out_channels))
    if gpu_flag
        weights1 = scale*rand(Complex{Float32}, modes1, modes2, in_channels, out_channels) |> gpu
        weights2 = scale*rand(Complex{Float32}, modes1, modes2, in_channels, out_channels) |> gpu
    else
        weights1 = scale*rand(Complex{Float32}, modes1, modes2, in_channels, out_channels)
        weights2 = scale*rand(Complex{Float32}, modes1, modes2, in_channels, out_channels)
    end
    return SpectralConv2d{Complex{Float32}, 4}(weights1, weights2)
end

function compl_mul2d(x::AbstractArray{Complex{Float32}}, y::AbstractArray{Complex{Float32}})
    # complex multiplication
    # x in (modes1, modes2, input channels, batchsize)
    # y in (modes1, modes2, input channels, output channels)
    # output in (modes1,modes2,output channles,batchsize)
    x_per = permutedims(x,[4,3,1,2]) # batchsize*in_channels*modes1*modes2
    y_per = permutedims(y,[3,4,1,2]) # in_channels*out_channels*modes1*modes2
    x_resh = reshape(x_per,size(x_per,1),size(x_per,2),:) # batchsize*in_channels*(modes1*modes2)
    y_resh = reshape(y_per,size(y_per,1),size(y_per,2),:) # in_channels*out_channels*(modes1*modes2)
    out_resh = batched_mul(x_resh,y_resh) # batchsize*out_channels*(modes1*modes2)
    out_per = reshape(out_resh,size(out_resh,1),size(out_resh,2),size(x,1),size(x,2)) # batchsize*out_channels*modes1*modes2
    out = permutedims(out_per,[3,4,2,1])
    return out
end

function (L::SpectralConv2d)(x::AbstractArray{Float32})
    # x in (size_x, size_y, channels, batchsize)
    x_ft = rfft(x,[1,2])
    modes1 = size(L.weights1,1)
    modes2 = size(L.weights1,2)
    out_ft = cat(cat(compl_mul2d(x_ft[1:modes1, 1:modes2,:,:], L.weights1),
                0f0im .* view(x_ft, 1:modes1, 1:size(x_ft,2)-2*modes2, :, :),
                compl_mul2d(x_ft[1:modes1, end-modes2+1:end,:,:], L.weights2),dims=2),
                0f0im .* view(x_ft, 1:size(x_ft,1)-modes1, :, :, :),dims=1)
    x = irfft(out_ft, size(x,1),[1,2])
end

mutable struct SimpleBlock2d
    fc0::Conv
    conv0::SpectralConv2d
    conv1::SpectralConv2d
    conv2::SpectralConv2d
    conv3::SpectralConv2d
    w0::Conv
    w1::Conv
    w2::Conv
    w3::Conv
    bn0::BatchNorm
    bn1::BatchNorm
    bn2::BatchNorm
    bn3::BatchNorm
    fc1::Conv
    fc2::Conv
end

@Flux.functor SimpleBlock2d

function SimpleBlock2d(modes1::Integer, modes2::Integer, width::Integer)
    block = SimpleBlock2d(
        Conv((1, 1), 3=>width),
        SpectralConv2d(width, width, modes1, modes2),
        SpectralConv2d(width, width, modes1, modes2),
        SpectralConv2d(width, width, modes1, modes2),
        SpectralConv2d(width, width, modes1, modes2),
        Conv((1, 1), width=>width),
        Conv((1, 1), width=>width),
        Conv((1, 1), width=>width),
        Conv((1, 1), width=>width),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        Conv((1, 1), width=>128),
        Conv((1, 1), 128=>1)
    )
    return block
end

function (B::SimpleBlock2d)(x::AbstractArray{Float32})
    x = B.fc0(x)
    x1 = B.conv0(x)
    x2 = B.w0(x)
    x = B.bn0(x1+x2)
    x = relu.(x)
    x1 = B.conv1(x)
    x2 = B.w1(x)
    x = B.bn1(x1+x2)
    x = relu.(x)
    x1 = B.conv2(x)
    x2 = B.w2(x)
    x = B.bn2(x1+x2)
    x = relu.(x)
    x1 = B.conv3(x)
    x2 = B.w3(x)
    x = B.bn3(x1+x2)
    x = B.fc1(x)
    x = relu.(x)
    x = B.fc2(x)
    return x
end

mutable struct Net2d
    conv1::SimpleBlock2d
end

@Flux.functor Net2d

function Net2d(modes::Integer, width::Integer)
    return Net2d(SimpleBlock2d(modes,modes,width))
end

function (NN::Net2d)(x::AbstractArray{Float32})
    x = NN.conv1(x)
    x = dropdims(x,dims=3)
end

ntrain = 1000
ntest = 100

batch_size = 20
learning_rate = 1f-3

epochs = 100
step_size = 100
gamma = 5f-1

modes = 12
width = 32

r = 5
h = Int(((421 - 1)/r) + 1)
s = h

### To load the train/test datasets, first download the Darcy_421.zip from https://drive.google.com/drive/folders/1UnbQh2WWc6knEHbLn-ZaXrKUZhp7pjt-

TRAIN = matread("data/piececonst_r421_N1024_smooth1.mat")
x_train_ = convert(Array{Float32},TRAIN["coeff"][1:ntrain,1:r:end,1:r:end][:,1:s,1:s])
y_train_ = convert(Array{Float32},TRAIN["sol"][1:ntrain,1:r:end,1:r:end][:,1:s,1:s])

TEST = matread("data/piececonst_r421_N1024_smooth2.mat")
x_test_ = convert(Array{Float32},TEST["coeff"][1:ntest,1:r:end,1:r:end][:,1:s,1:s])
y_test = convert(Array{Float32},TEST["sol"][1:ntest,1:r:end,1:r:end][:,1:s,1:s])

x_train_ = permutedims(x_train_,[2,3,1])
x_test_ = permutedims(x_test_,[2,3,1])
y_train_ = permutedims(y_train_,[2,3,1])
y_test = permutedims(y_test,[2,3,1])

x_normalizer = UnitGaussianNormalizer(x_train_)
x_train_ = encode(x_normalizer,x_train_)
x_test_ = encode(x_normalizer,x_test_)

y_normalizer = UnitGaussianNormalizer(y_train_)
y_train = encode(y_normalizer,y_train_)

x = reshape(collect(range(0f0,stop=1f0,length=s)), :, 1)
z = reshape(collect(range(0f0,stop=1f0,length=s)), 1, :)

grid = zeros(Float32,s,s,2)
grid[:,:,1] = repeat(z,s)
grid[:,:,2] = repeat(x',s)'

x_train = zeros(Float32,s,s,3,ntrain)
x_train[:,:,1,:] = x_train_

for i = 1:ntrain
    x_train[:,:,2,i] = grid[:,:,1]
    x_train[:,:,3,i] = grid[:,:,2]
end

x_test = zeros(Float32,s,s,3,ntest)
x_test[:,:,1,:] = x_test_

for i = 1:ntest
    x_test[:,:,2,i] = grid[:,:,1]
    x_test[:,:,3,i] = grid[:,:,2]
end

train_loader = Flux.Data.DataLoader((x_train, y_train); batchsize = batch_size, shuffle = true)
test_loader = Flux.Data.DataLoader((x_test, y_test); batchsize = batch_size, shuffle = false)

y_normalizer.mean_ = y_normalizer.mean_ |> gpu
y_normalizer.std_ = y_normalizer.std_   |> gpu
y_normalizer.eps_ = y_normalizer.eps_   |> gpu

NN = Net2d(modes, width) |> gpu

w = Flux.params(NN)
Flux.trainmode!(NN, true)
opt = Flux.Optimise.ADAMW(learning_rate, (0.9f0, 0.999f0), 1f-4)

prog = Progress(ntrain * epochs)

Loss = zeros(Float32,epochs*Int(ntrain/batch_size))

iter = 0
for ep = 1:epochs
    for (x,y) in train_loader
        global iter = iter + 1
        grads = gradient(w) do
            if gpu_flag
                x = x |> gpu
                y = y |> gpu
            end
            out = decode(y_normalizer,NN(x))
            y_n = decode(y_normalizer,y)
            global loss = Flux.mse(out,y_n;agg=sum)
            return loss
        end
        Loss[iter] = loss
        for p in w
            Flux.Optimise.update!(opt, p, grads[p])
        end
        ProgressMeter.next!(prog; showvalues = [(:loss, loss), (:epoch, ep)])
    end
end

if gpu_flag
    NN = NN |> cpu
    w = convert.(Array,w) |> cpu
end

# Define result directory
mkpath(datadir("TrainedNet"))
BSON.@save "data/TrainedNet/darcynet_$epochs.bson" NN w batch_size Loss modes width learning_rate epochs gamma step_size ntrain ntest s r h


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

==(w::WeakRef, v::WeakRef) = isequal(w.value, v.value)
==(w::WeakRef, v) = isequal(w.value, v)
==(w, v::WeakRef) = isequal(w, v.value)

"""
    finalizer(f, x)

Register a function `f(x)` to be called when there are no program-accessible references to
`x`, and return `x`. The type of `x` must be a `mutable struct`, otherwise the behavior of
this function is unpredictable.
"""
function finalizer(@nospecialize(f), @nospecialize(o))
    if isimmutable(o)
        error("objects of type ", typeof(o), " cannot be finalized")
    end
    ccall(:jl_gc_add_finalizer_th, Cvoid, (Ptr{Cvoid}, Any, Any),
          Core.getptls(), o, f)
    return o
end

function finalizer(f::Ptr{Cvoid}, o::T) where T
    @_inline_meta
    if isimmutable(T)
        error("objects of type ", T, " cannot be finalized")
    end
    ccall(:jl_gc_add_ptr_finalizer, Cvoid, (Ptr{Cvoid}, Any, Ptr{Cvoid}),
          Core.getptls(), o, f)
    return o
end

"""
    finalize(x)

Immediately run finalizers registered for object `x`.
"""
finalize(@nospecialize(o)) = ccall(:jl_finalize_th, Cvoid, (Ptr{Cvoid}, Any,),
                                   Core.getptls(), o)

module GC

"""
    GC.gc()

Perform garbage collection.

!!! warning
    Excessive use will likely lead to poor performance.
"""
gc(full::Bool=true) = ccall(:jl_gc_collect, Cvoid, (Int32,), full)

"""
    GC.enable(on::Bool)

Control whether garbage collection is enabled using a boolean argument (`true` for enabled,
`false` for disabled). Return previous GC state.

!!! warning
    Disabling garbage collection should be used only with caution, as it can cause memory
    use to grow without bound.
"""
enable(on::Bool) = ccall(:jl_gc_enable, Int32, (Int32,), on) != 0

"""
    GC.@preserve x1 x2 ... xn expr

Temporarily protect the given objects from being garbage collected, even if they would
otherwise be unreferenced.

The last argument is the expression during which the object(s) will be preserved.
The previous arguments are the objects to preserve.
"""
macro preserve(args...)
    syms = args[1:end-1]
    for x in syms
        isa(x, Symbol) || error("Preserved variable must be a symbol")
    end
    s, r = gensym(), gensym()
    esc(quote
        $s = $(Expr(:gc_preserve_begin, syms...))
        $r = $(args[end])
        $(Expr(:gc_preserve_end, s))
        $r
    end)
end

end # module GC

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
(==){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = ((a.lo == b.lo) & (a.hi == b.hi))
(!=){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = ((a.lo != b.lo) | (a.hi != b.hi))
(<=){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = ((a.lo <= b.lo) & (a.hi <= b.hi))
(>=){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = ((a.lo >= b.lo) & (a.hi >= b.hi))
(< ){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = (!(a >= b))
(> ){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = (!(a <= b))

(isequal){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = (a == b)
(isless ){G<:Grasp, R<:Real}(a::Rvl{G,R}, b::Rvl{G,R}) = (a <  b)


text:
#=
Provides a function to solve for asset prices under one set of beliefs in the
Harrison -- Kreps model.

Authors: Shunsuke Hori
=#
"""
Function to Solve Single Beliefs
"""
function price_singlebeliefs(transition, dividend_payoff; beta=.75)
    # First compute inverse piece
    imbq_inv = inv(eye(size(transition,1)) - beta*transition)

    # Next compute prices
    prices = beta * ((imbq_inv*transition)* dividend_payoff)

    return prices
end


text:
export
    Line,
    get_polar_angle,
    get_distance,
    get_side

immutable Line
    C::VecE2
    θ::Float64
end
Line(A::VecE2, B::VecE2) = Line((A+B)/2, atan2(B - A))

Base.:-(line::Line, V::VecE2) = Line(line.C - V, line.θ)
Base.:+(line::Line, V::VecE2) = Line(line.C + V, line.θ)

get_polar_angle(line::Line) = line.θ
rot(line::Line, Δθ::Float64) = Line(line.C, line.θ+Δθ)
Base.rot180(line::Line) = rot(line, 1π)
Base.rotl90(line::Line) = rot(line,  π/2)
Base.rotr90(line::Line) = rot(line, -π/2)

"""
The distance between the line and the point P
"""
function get_distance(line::Line, P::VecE2)

    ab = polar(1.0, line.θ)
    pb = P - line.C

    denom = abs2(ab)
    if denom ≈ 0.0
        return 0.0
    end

    r = dot(ab, pb)/denom
    return abs(P - (line.C + r*ab))
end

"""
What side of the line you are on
Is -1 if on the left, 1 if on the right, and 0 if on the line
"""
function get_side(line::Line, p::VecE2, ε::Float64=1e-10)
    ab = polar(1.0, line.θ)
    signed_dist = ab.x*(p.y-line.C.y) - ab.y*(p.x-line.C.x)
    if abs(signed_dist) < ε
        return 0
    else
        return sign(signed_dist)
    end
end


text:
"""
$(DocStringExtensions.README)
"""
module DataDrivenDiffEq

using DocStringExtensions
using LinearAlgebra
using DiffEqBase
using CommonSolve
using ModelingToolkit

using Distributions
using QuadGK
using Statistics
using DataInterpolations
using Parameters
using Random
using Measurements

using Requires
using ProgressMeter
using Reexport
using Compat
using DocStringExtensions
using RecipesBase

@reexport using DiffEqBase: solve
@reexport using ModelingToolkit: states, parameters, independent_variable, observed, controls, get_iv
@reexport using DataInterpolations: ConstantInterpolation, LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox, Curvefit
using Symbolics: scalarize, variable

using ModelingToolkit: AbstractSystem
using ModelingToolkit: value, operation, arguments, istree, get_observed
# Basis and Koopman
abstract type AbstractBasis <: AbstractSystem end
abstract type AbstractKoopman <: AbstractBasis end
# Collect the DataInterpolations Methods into an Interpolation Type
abstract type AbstractInterpolationMethod end
abstract type CollocationKernel end

# Algortihms for Koopman
abstract type AbstractKoopmanAlgorithm end

# Abstract symbolic_regression
abstract type AbstractSymbolicRegression end

# Problem and solution
abstract type AbstractDataDrivenProblem{dType, cType, probType} end
abstract type AbstractDataDrivenSolution end

# Optimizer
abstract type AbstractProximalOperator end;
abstract type AbstractOptimizer{T} end;
abstract type AbstractSubspaceOptimizer{T} <: AbstractOptimizer{T} end;

## Basis

include("./basis/build_function.jl")
include("./basis/utils.jl")
include("./basis/type.jl")
export Basis
export jacobian, dynamics
export free_parameters, implicit_variables

include("./utils/basis_generators.jl")
export chebyshev_basis, monomial_basis, polynomial_basis
export sin_basis, cos_basis, fourier_basis

include("./utils/collocation.jl")
export InterpolationMethod
export EpanechnikovKernel, UniformKernel, TriangularKernel,QuarticKernel
export TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel
export LogisticKernel, SigmoidKernel, SilvermanKernel
export collocate_data

include("./utils/utils.jl")
export AIC, AICC, BIC
export optimal_shrinkage, optimal_shrinkage!
export burst_sampling, subsample

## Sparse Regression

include("./optimizers/Optimize.jl")
export SoftThreshold, HardThreshold,ClippedAbsoluteDeviation
export sparse_regression!
export init, init!, set_threshold!, get_threshold
export STLSQ, ADMM, SR3
export ImplicitOptimizer

## Koopman

include("./koopman/type.jl")
export Koopman
export operator, generator
export is_stable, is_discrete, is_continuous
export modes, frequencies, outputmap, updatable
export update!

include("./koopman/algorithms.jl")
export DMDPINV, DMDSVD, fbDMD, TOTALDMD

## Problem and Solution
# Use to distinguish the problem types
@enum DDProbType begin
    Direct=1 # Direct problem without further information
    Discrete=2 # Time discrete problem
    Continuous=3 # Time continous problem
end

# Define some alias type for easier dispatch
const AbstractDirectProb{N,C} = AbstractDataDrivenProblem{N,C,DDProbType(1)}
const AbstractDiscreteProb{N,C} = AbstractDataDrivenProblem{N,C,DDProbType(2)}
const AbstracContProb{N,C} = AbstractDataDrivenProblem{N,C,DDProbType(3)}

include("./problem/type.jl")

export DataDrivenProblem
export DiscreteDataDrivenProblem, ContinuousDataDrivenProblem, DirectDataDrivenProblem
export is_autonomous, is_discrete, is_direct, is_continuous, is_parametrized, has_timepoints
export is_valid, @is_applicable, get_name

include("./problem/sample.jl")
export DataSampler, Split, Batcher

# Result selection
select_by(x, y::AbstractMatrix) = y
select_by(x, sol) = select_by(Val(x), sol)

include("./solution.jl")
export DataDrivenSolution
export result, parameters, parameter_map, algorithm
export output, metrics, error, aic, determination, get_problem

include("./solve/common.jl")
export DataDrivenCommonOptions
include("./solve/sparse_identification.jl")
#include("./solve/sindy.jl")
include("./solve/koopman.jl")
#export solve

include("./recipes/problem_result.jl")

# Optional
function __init__()
    # Load and export OccamNet
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin

        using .Flux
        include("./symbolic_regression/occamnet.jl")

        export OccamNet,ProbabilityLayer
        export set_temp!
        export probability, logprobability
        export probabilities, logprobabilities
        export OccamSR

        @info "DataDrivenDiffEq : OccamNet is available."
    end

    @require SymbolicRegression = "8254be44-1295-4e6a-a16d-46603ac705cb" begin

        using .SymbolicRegression
        include("./symbolic_regression/symbolic_regression.jl")
        export EQSearch

        @info "DataDrivenDiffEq : Symboolic Regression is available."
    end

end

end # module


text:
#=
format: ASCII
The structure of the file is as follows:
    1 x comment line "# header"
    1 x comment line with labels of header variables
    1 x header (exactly 1 line):
    1 x common data
    n x configuration(s) of the players
    1 x comment line "# body"
    1 x comment line with labels of data block variables
    m x data set (exactly 1 line):
    1 x common data of simulation step
    n x data of the players

The following representations shall be used for the different data types:
    strings must be quoted ("OwnPlayer")
    floating point numbers shall be in scientific notation with 16 digits
    (0.1234567890123456E+4)
    integer numbers shall contain the relevant digits
    state masks are written as (unsigned) integer numbers
    The separator between data entries is a comma (,).
=#

export
    read_vires_csv,
    export_to_bosch_csv,

    inbouds,
    carind_inbounds,
    frameind_inbounds,

    set!,
    get_nplayers,
    get_nframes,
    get_column_symbol,
    get_pos_inertial,
    get_orientation,
    get_speed_body,
    get_inertial_rate,
    get_acceleration_body

type PlayerHeader
    ptype           :: Int     # player type
    id              :: UInt    # player ID
    name            :: String  # player name
    file            :: String  # ???
    length          :: Float64 # player length [m]
    width           :: Float64 # player width [m]
    center_offset_x :: Float64 # x-offset between center of geometry and player reference point (in body frame) [m]
    center_offset_y :: Float64 # x-offset between center of geometry and player reference point (in body frame) [m]
    model           :: Int     # graphical model
    category        :: UInt    # category
    ctrl            :: Int     # controller type
    ctrl_arg        :: Int     # controller argument
    master          :: Int     # player acted as master in the simulation {0,1}

    PlayerHeader() = new()
end

type ViresCSV
    header::Dict{Symbol,String}
    player_headers::Vector{PlayerHeader}
    body::DataFrame
end

function read_vires_csv(pathname::String)
    open(pathname) do fin

        nbytes = filesize(pathname)

        line = readline(fin)
        nbytes -= length(line)
	    @assert(contains(line, "# header"))

        # extract header variables and put them in a dictionary Dict{Symbol,String}

        line = readline(fin)
        nbytes -= length(line)
        substrings = split(line, ",") # split_by_commas
        map!(str->replace(str, r"(^\s+|\s+\Z)", ""), substrings) # removed_leading_and_trailing_whitespace
        map!(str->replace(str, " ", "_"), substrings) # replace_spaces_with_underscores
        map!(str->replace(str, r"[.,]", ""), substrings) # remove periods and commas
        substrings = filter!(s->!isempty(s), substrings)
        dictkeys = map(str->Symbol(str), substrings)

        line = readline(fin)
        nbytes -= length(line)
        substrings = String[]
        i = 1
        while i < length(line)
            k = findnext(line, '\"', i)
            if k != 0
                l = findnext(line, '\"', k+1)
                @assert(l > 0)
                append!(substrings, split(line[i:k-1]))
                push!(substrings, line[k+1:l-1])
                i = l+1
            else
                append!(substrings, split(line[i:end]))
                break
            end
        end

        map!(str->replace(str, r"(^\s+|\s+\Z)", ""), substrings) # removed_leading_and_trailing_whitespace
        map!(str->replace(str, " ", "_"), substrings) # replace_spaces_with_underscores
        map!(str->replace(str, r"[.,]", ""), substrings) # remove periods and commas
        dictvalues = filter(s->!isempty(s), substrings)

        # println(dictvalues)
        # println(length(dictkeys))
        # println(length(dictvalues))
        @assert(length(dictkeys) == length(dictvalues))

        player_headers = PlayerHeader[]

        header_variables = Dict{Symbol,String}()
        for (key,value) in zip(dictkeys, dictvalues)
            if key == :player_type
                push!(player_headers, PlayerHeader())
                player_headers[end].ptype = parse(Int,value)
            elseif key == :player_ID
                player_headers[end].id = parse(Int,value)
            elseif key == :player_name
                player_headers[end].name = value
            elseif key == :file_name
                player_headers[end].file = value
            elseif key == :file_length
                player_headers[end].length = float(value)
            elseif key == :player_width
                player_headers[end].width = float(value)
            elseif key == :player_offset_x
                player_headers[end].offset_x = float(value)
            elseif key == :player_offset_y
                player_headers[end].ctr_offset_y = float(value)
            elseif key == :player_model
                player_headers[end].model = parse(Int,value)
            elseif key == :player_category
                player_headers[end].category = parse(Int,value)
            elseif key == :player_ctrl
                player_headers[end].ctrl = parse(Int,value)
            elseif key == :player_ctrl_arg
                player_headers[end].ctrl_arg = parse(Int,value)
            elseif key == :master
                player_headers[end].master = parse(Int,value)
            else
                header_variables[key] = value
            end
        end
        @assert(parse(Int,get(header_variables, :no_of_players, "0")) == length(player_headers))

        line = readline(fin)
        nbytes -= length(line)
        @assert(contains(line, "# body"))
        body = DataFrames.readtable(fin, nbytes; header=true)

        ViresCSV(header_variables, player_headers, body)
    end
end

inbounds(csv::ViresCSV, carind::Int, frameind::Int) =
    carind_inbounds(csv, carind) && frameind_inbounds(csv, frameind)
carind_inbounds(csv::ViresCSV, carind::Int) = 0 ≤ carind < get_nplayers(csv)
frameind_inbounds(csv::ViresCSV, frameind::Int) = 1 ≤ frameind ≤ nrow(csv.body)

Base.get(csv::ViresCSV, frameind::Int, sym::Symbol) = csv.body[frameind, sym]
function Base.get(csv::ViresCSV, carind::Int, frameind::Int, str::String)
    @assert(carind_inbounds(csv, carind))
    @assert(frameind_inbounds(csv, frameind))
    sym = get_column_symbol(carind, str)
    csv.body[frameind, sym]
end
function set!(csv::ViresCSV, carind::Int, frameind::Int, str::String, value::Any)
    @assert(inbounds(csv, carind, frameind))
    str2 = get_column_symbol(carind, str)
    csv.body[frameind, str2] = value
    csv
end

get_nplayers(csv::ViresCSV) = length(csv.player_headers)
get_nframes(csv::ViresCSV) = nrow(csv.body)
function get_column_symbol(carind::Int, str::String)
    if carind != 0 # not ego car
        str = @sprintf("%s_%d", str, carind)
    end
    symbol(str)
end
function get_pos_inertial(csv::ViresCSV, carind::Int, frameind::Int)
    # returns a tuple, (x::Float64,y::Float64,z::Float64) in inertial for the vehicle
    x = get(csv, carind, frameind, "inertialX")::Float64
    y = get(csv, carind, frameind, "inertalY")::Float64
    z = get(csv, carind, frameind, "inertialZ")::Float64
    (x,y,z)
end
function get_orientation(csv::ViresCSV, carind::Int, frameind::Int)
    # returns a tuple, (heading::Float64,pitch::Float64,roll::Float64) in inertial for the vehicle
    h = get(csv, carind, frameind, "hdg")::Float64
    p = get(csv, carind, frameind, "pitch")::Float64
    r = get(csv, carind, frameind, "roll")::Float64
    (h,p,r)
end
function get_speed_body(csv::ViresCSV, carind::Int, frameind::Int)
    # returns a tuple, (vx::Float64,vy::Float64,vz::Float64) containing
    # the speed in [m/s] inertial for the vehicle
    vx = get(csv, carind, frameind, "vPlayerX")::Float64
    vy = get(csv, carind, frameind, "vPlayerY")::Float64
    vz = get(csv, carind, frameind, "vPlayerZ")::Float64
    (vx,vy,vz)
end
function get_inertial_rate(csv::ViresCSV, carind::Int, frameind::Int)
    # returns a tuple, (hdot::Float64,pdot::Float64,rdot::Float64)
    # containing the rate of change of orientation in the inertial frame for the vehicle
    dh = get(csv, carind, frameind, "hdgDot")::Float64
    dp = get(csv, carind, frameind, "pitchDot")::Float64
    dr = get(csv, carind, frameind, "rollDot")::Float64
    (dh,dp,dr)
end
function get_acceleration_body(csv::ViresCSV, carind::Int, frameind::Int)
    # returns a tuple, (ax::Float64,ay::Float64,az::Float64) containing
    # the speed in [m/s] inertial for the vehicle
    ax = get(csv, carind, frameind, "aPlayerX")::Float64
    ay = get(csv, carind, frameind, "aPlayerY")::Float64
    az = get(csv, carind, frameind, "aPlayerZ")::Float64
    (ax,ay,az)
end

function _euler2quat( roll::Real, pitch::Real, yaw::Real )

    #=
    Converts roll, pitch, yaw in radians to the corresponding unit quaternion (w,x,y,z)
    =#

    c_r = cos(0.5roll)
    s_r = sin(0.5roll)
    c_p = cos(0.5pitch)
    s_p = sin(0.5pitch)
    c_y = cos(0.5yaw)
    s_y = sin(0.5yaw)

    w2 = c_r*c_p*c_y + s_r*s_p*s_y
    x2 = s_r*c_p*c_y - c_r*s_p*s_y
    y2 = c_r*s_p*c_y + s_r*c_p*s_y
    z2 = c_r*c_p*s_y - s_r*s_p*c_y

    return [w2, x2, y2, z2]
end
function _global2ego( egocarGx::Real, egocarGy::Real, egocarYaw::Real, posGx::Real, posGy::Real )

    #=
    Convert an (x,y) in inertial coordinates to (x,y) relative to the ego vehicle's body frame
    =#

    # translate to be relative to ego car
    pt = [posGx-egocarGx, posGy-egocarGy]

    # rotate to be in ego car frame
    cθ, sθ = cos(egocarYaw), sin(egocarYaw)
    R = [ cθ sθ;
         -sθ cθ]
    pt = R*pt

    return (pt[1], pt[2]) # posEx, posEy
end
function _ego2global(egocarGx::Real, egocarGy::Real, egocarYaw::Real, posEx::Real, posEy::Real )

    #=
    Convert an (x,y) in body coordinates to (x,y) relative to the global frame
    =#

    pt = [posEx, posEy]
    eo = [egocarGx, egocarGy]

    # rotate
    c, s = cos(egocarYaw), sin(egocarYaw)
    R = [ c -s;
          s  c]
    pt = R*pt

    # translate
    pt += eo

    (pt[1], pt[2]) # posGx, posGy
end
function export_to_bosch_csv(io::IO, csv::ViresCSV;
    start_frame::Int=5, # first frame to export
    end_frame::Int=get_nframes(csv) # last frame to export
    )
    #=
    Exports the ViresCSV to a csv file following the Bosch CSV format.
    The Bosch CSV format (made by Jorge), is a CSV file with the following columns:

    entry                          row number, starting at 0
    timings                        frame time, unix epoch
    status                         ego vehicle status (99=CONTROL_STATUS_VIRES_AUTO for ghostdriver)
    global position x/y/z          3 columns, one for each of the inertial coordinates, [m]
    global rotation w/x/y/z        quaternion defining the vehicle orientation in global coordinates
    odom velocity x/y/z            ego velocity in the body frame [m/s]
    odom acceleration x/y/z        ego acceleration in the body frame [m/s²]
    lane                           deprecated
    offsetwrtcenterline            deprecated
    closest centerline point x/y   deprecated
    right lane x/y                 deprecated
    # repeated for a number of extra vehicles, with increasing digit, starting at 0
    car id0                        vehicle unique id, integer
    ego x0/y0                      vehicle ego-relative position [m]
    v x0/yo                        vehicle ego-relative velocity [m/s]
    global x0/y0                   vehicle global position [m]
    global v_x0/v_y0               vehicle global velocity [m/s]
    car angle0                     vehicle heading w.r.t global
    lane0                          deprecated
    offset wrt centerline0         deprecated
    tangent theta0                 deprecated
    angle to lane0                 deprecated
    lane v_x0                      deprecated
    lane v_y0                      deprecated
    lane x0                        deprecated
    lane y0                        deprecated
    kappa0                         deprecated
    kappa_prime0                   deprecated
    reference lane0                deprecated
    reference lane x0              deprecated
    reference lane y0              deprecated
    =#

    CONTROL_STATUS_VIRES_AUTO = 99

    n_other_cars = get_nplayers(csv) - 1
    print(io, "entry, timings, status, global position x, global position y, global position z, global rotation w,global rotation x, global rotation y, global rotation z, odom velocity x, odom velocity y, odom velocity z, odom acceleration x, odom acceleration y, odom acceleration z, lane, offsetwrtcenterline, closest centerline point x, closest centerline point y, right lane x, right lane y")
    for i = 1 : n_other_cars
        id = i-1
        print(io, ",car id$(id),ego x$(id),ego y$(id),v x$(id),v y$(id),global x$(id),global y$(id),global v_x$(id),global v_y$(id),car angle$(id),lane$(id),offset wrt centerline$(id),tangent theta$(id),angle to lane$(id),lane v_x$(id),lane v_y$(id),lane x$(id),lane y$(id),kappa$(id),kappa_prime$(id),reference lane$(id),reference lane x$(id),reference lane y$(id)",)
    end
    print(io, "\n")

    export_frameind_count = 0
    for i = start_frame : end_frame
        x,y,z = get_pos_inertial(csv, 0, i)
        h,p,r = get_orientation(csv, 0, i)
        quat = _euler2quat(r, p, h)
        vx,vy,vz = get_speed_body(csv, 0, i)
        ax,ay,az = get_acceleration_body(csv, 0, i)
        # vxi,vyi = _ego2global(0,0,h,vx,vy)

        print(io, export_frameind_count, ",", get(csv, i, :_simTime), ",", CONTROL_STATUS_VIRES_AUTO, ",", x, ",", y,
            ",", z, ",", quat[1], ",", quat[2], ",", quat[3], ",", quat[4], ",", vx, ",", vy, ",",
            vz, ",", ax, ",", ay, ",", az, ",", 1, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0)

        for j = 1 : n_other_cars
            id = j-1
            cx,cy,cz = get_pos_inertial(csv, j, i)
            cvxb,cvyb,cvzb = get_speed_body(csv, j, i)
            ch = get_orientation(csv, j, i)[1]
            cvx,cvy = _ego2global(0,0,ch,cvxb,cvyb)
            ex,ey = _global2ego(x, y, h, cx, cy)
            velEx,velEy = _global2ego(0, 0, h, cvx, cvy)

            print(io, ",", id, ",", ex, ",", ey, ",", velEx, ",", velEy, ",", cx, ",", vy, ",",
                  cvx, ",", cvy, ",", get(csv, j, i, "hdg"), ",", 1, ",", 0.0, ",", 0.0, ",",
                  0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0, ",", 0.0)
        end
        print(io, "\n")

        export_frameind_count += 1
    end
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
function remove_extension(s)
    arr = split(s, ".")
    return (length(arr) == 1) ? s : join(arr[1:(end-1)])
end

function pdf_to_eps(abspath, fname, dir)
    eps_name = joinpath(dir, "$(fname).eps")
    run(`inkscape $abspath -E $eps_name --export-ignore-filters --export-ps-level=2`)
    return eps_name
end

function printfig(fname; kwargs...)
    dir, base = Base.Filesystem.dirname(fname), Base.Filesystem.basename(fname)
    filename = remove_extension(base)
    workdir = mktempdir()
    pdf_name = joinpath(workdir, "$(filename).pdf")
    savefig(pdf_name)
    eps_path = pdf_to_eps(pdf_name, filename, workdir)
    cp(eps_path, joinpath(dir, basename(eps_path)), remove_destination=true)
end


text:
# classification demo

# avoid unrelated libz error
# using ImageMagick
# using Makie

using KernelDensityEstimate
using Gadfly
using Distributions
using StatsBase, Random
using IteratedProbabilisticClassification
# don't mind syntax warnings, these are being removed from package dependencies over time

# modified EM classification algorithm has these major tuning parameter
global localparams = TuningParameters(
  Mfair=300,   # Cluster distribution approximation accuracy, more is better: 50 is bare minium, and computation slow above ~1000
  EMiters=40   # expectation maximization iterations to refine Cluster belief estimate classification assignments: 5 is minimum, and diminishing returns beyond ~70
)

# Example: simulated data with ground truth for reference
N1, N2 = 5000, 1000
# N1, N2 = 4000, 3000
data, groundtruth = simdata02_true(N1=N1,N2=N2)

# actual distribution of data
Gadfly.set_default_plot_size(15cm, 10cm)
# what does the data look like
Gadfly.plot(x=data.samples[1,:],y=data.samples[2,:], Geom.histogram2d)

# Ground truth distribution of underlying data clusters
# Gadfly.set_default_plot_size(15cm, 10cm)
# a generic 1D plotting tool for IteratedProbablisticClassification
plotUtil2D([2;1], groundtruth=groundtruth)

# belief from expert prediction, also used for initialization
c1_expt_pts = rand(MvNormal([5.0;5.0],[3.0 -1.0; -1.0 3.0]),100)
c1_expert = kde!(reshape([5.0,5.0],2,1), [3.0,3.0]) # Pretty much a normal distribution
c1_expert = kde!(c1_expt_pts)
c2_expert = kde!(reshape([-5.0,-6.0],2,1), [4.0,4.0]) # Pretty much a normal distribution

# cs is the main structure which classification algorithm will operate on, and modify during execution
cs = ClassificationSystem(
  2,                              # number of categories or classes available for assignment
  Dict(1=>"class1", 2=>"class2"), # names of classes
  Dict(1=>"blue", 2=>"red"),      # default plotting colors
  Dict(1=>c1_expert, 2=>c2_expert), # user forcing behaviour (expert guidance); required
  Dict(1=>deepcopy(c1_expert), 2=>deepcopy(c2_expert)), # initialize temporal prediction (forward-backward smoothing)
  Dict(1=>deepcopy(c1_expert), 2=>deepcopy(c2_expert)), # initialize current belief (0 iterations)
  rand(Categorical([0.5;0.5]),length(data.samples)) # initialize samples with random assignment, 50/50% in this 2 class case
);
println()

# expert may modify expertBelief using output from previous run of classification algorithm
# pts = getPoints(cs.currentBelief[1])
# cs.expertBelief[1] = kde!(dispersePoints(pts, MvNormal([-1.0],[1.0]) ))
# pts = getPoints(cs.currentBelief[2])
# cs.expertBelief[2] = kde!(dispersePoints(pts, MvNormal([2.0],[1.0]) ))
# println()

Gadfly.set_default_plot_size(20cm, 7cm)
# plotUtil1D(cs=cs) # if we didn't have gt
plotUtil2D(1,cs=cs)

# simulation data allows us to debug with absolute knowledge
dbg = defaultDebugResults()

# do the classification
stats = EMClassificationRun!(localparams, cs, data, debug=dbg, groundtruth=groundtruth);
println()

Gadfly.set_default_plot_size(20cm, 10cm)
plotUtil1D(cs=cs, groundtruth=groundtruth,
  drawcurrent=true,
  expertcolor=Dict(1=>"gray",2=>"gray")
)

sum(cs.assignment .== 1), sum(cs.assignment .== 2)

Gadfly.set_default_plot_size(20cm, 10cm)
plot(
layer(x=data.samples[1,cs.assignment .== 2], Geom.histogram, Theme(default_color=colorant"red")),
layer(x=data.samples[1,cs.assignment .== 1], Geom.histogram)
)

# proxy to convergence of classification algorithm, always available from stats structure
Gadfly.set_default_plot_size(20cm, 15cm)
plotEMStatus(localparams,stats)

# only available when ground truth data is available
Gadfly.set_default_plot_size(20cm, 15cm)
plotClassificationStats(localparams, dbg)

# pts = getPoints(cs.expertBelief[1])
# pts2 = dispersePoints(pts, MvNormal([0.0],[3.0]) )
# p2 = kde!(pts2)
# plotKDE([cs.expertBelief[1];p2])

#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

if !isdefined(@__MODULE__, :DefsCommentId)
    const DefsCommentId = String
else
    @warn("Skipping redefinition of DefsCommentId to String")
end


text:
module Day21

export day21

using CombinedParsers
using CombinedParsers.Regexp: word, words, whitespace_maybe, whitespace_horizontal

function parseinput(filename)
    p = Sequence(:ingredients => join(Repeat(word), whitespace_maybe), " (contains " , :allergens => join(Repeat(word), ", "), ")")
    possibilities = Dict{String,Union{Set{String}, Nothing}}()
    ingredients = Set{String}()
    for line in readlines(filename)
        t = line |> p
        for allergen in t[:allergens]
            tmp = get!(possibilities, allergen, nothing)
            possibilities[allergen] = isnothing(tmp) ? Set(t[:ingredients]) : intersect(tmp, Set(t[:ingredients]))
            union!(ingredients, Set(t[:ingredients]))
        end
    end
    (ingredients, possibilities)
end

function day21(::Val{:a}, filename)
    (all_ingredients, possibilities) = parseinput(filename)
    contains_allergen = Set{String}()
    changes = true
    while changes
        changes = false
        todo = Set{String}()
        for (allergen, ingredients) in possibilities
            if length(ingredients) == 1
                changes = true
                ingred = pop!(ingredients)
                push!(todo, ingred)
                push!(contains_allergen, ingred)
                delete!(possibilities, allergen)
            end
        end
        for ingred in todo
            for allergen in keys(possibilities)
                delete!(possibilities[allergen], ingred)
            end
        end
    end
    allergen_free = setdiff(all_ingredients,contains_allergen)
    counts = Dict{String,Int}()
    for line in readlines(filename)
        words = split(line[1:findfirst('(', line)-1], ' ')
        for ingred in allergen_free
            if in(ingred, words)
                counts[ingred] = get(counts, ingred, 0) + 1
            end
        end
    end
    sum(values(counts))
end

function day21(::Val{:b}, filename)
    (all_ingredients, possibilities) = parseinput(filename)
    contained = Dict{String,String}()
    changes = true
    while changes
        changes = false
        todo = Set{String}()
        for (allergen, ingredients) in possibilities
            if length(ingredients) == 1
                changes = true
                ingred = pop!(ingredients)
                push!(todo, ingred)
                contained[allergen] = ingred
                delete!(possibilities, allergen)
            end
        end
        for ingred in todo
            for allergen in keys(possibilities)
                delete!(possibilities[allergen], ingred)
            end
        end
    end
    allergens = sort(collect(keys(contained)))
    result = ""
    for allergen in allergens
        result *= "," * contained[allergen]
    end
    result[2:end]
end

end


text:
include("Ecl.jl")

using DataFrames

df = read_all_results("C:\\Test\\VOLVE_2016.SMSPEC", r"^W");


text:
@safetestset "ishexcolor" begin
    using Validator

    @test Validator.ishexcolor("#ff0000ff") === true
    @test Validator.ishexcolor("#ff0034") === true
    @test Validator.ishexcolor("#CCCCCC") === true
    @test Validator.ishexcolor("0f38") === true
    @test Validator.ishexcolor("fff") === true
    @test Validator.ishexcolor("#f00") === true

    @test Validator.ishexcolor("#ff") === false
    @test Validator.ishexcolor("fff0a") === false
    @test Validator.ishexcolor("#ff12FG") === false
end


text:
using Documenter

push!(LOAD_PATH, "../src/")
using LinearInterpolators

DEPLOYDOCS = (get(ENV, "CI", nothing) == "true")

makedocs(
    sitename = "Linear interpolators for Julia",
    format = Documenter.HTML(
        prettyurls = DEPLOYDOCS,
    ),
    authors = "Éric Thiébaut and contributors",
    pages = ["index.md",
             "install.md",
             "interpolation.md",
             "library.md",
             "notes.md"]
)

if DEPLOYDOCS
    deploydocs(
        repo = "github.com/emmt/LinearInterpolators.jl.git",
    )
end


text:
module DoublePendulum

using NeuralNetworkAnalysis, LaTeXStrings
using NeuralNetworkAnalysis: Specification
import Plots, DifferentialEquations
using Plots: plot, plot!, savefig, font, Measures.mm

# problem

m = 0.5
L = 0.5
c = 0.0
g = 1.0
gL = g/L
mL = m*L^2

@taylorize function double_pendulum!(dx, x, p, t)
    x₁, x₂, x₃, x₄, T₁, T₂ = x

    Δ12 = x₁ - x₂
    ★ = cos(Δ12)
    x3sin12 = x₃^2 * sin(Δ12)
    x4sin12 = x₄^2 * sin(Δ12) / 2
    gLsin1 = gL * sin(x₁)
    gLsin2 = gL * sin(x₂)
    T1_frac = (T₁ - c * x₃) / (2 * mL)
    T2_frac = (T₂ - c * x₄) / mL
    bignum = x3sin12 - ★ * (gLsin1 - x4sin12 + T1_frac) + gLsin2 + T2_frac
    denom = ★^2 / 2 - 1

    dx[1] = x₃
    dx[2] = x₄
    dx[3] = ★ * bignum / (2 * denom) - x4sin12 + gLsin1 + T1_frac
    dx[4] = - bignum / denom
end;

controller = read_nnet_sherlock(@modelpath "" "controller_Double-Pendulum_mr_sherlock")

X₀ = BallInf(fill(1.15, 4), 0.15)
X₀_small = BallInf(fill(1.000005, 4), 0.000005)
U₀ = ZeroSet(2)
vars_idx = Dict(:state_vars=>1:4, :control_vars=>5:6)
ivp = @ivp(x' = double_pendulum!(x), dim: 6, x(0) ∈ X₀ × U₀)
ivp_small = @ivp(x' = double_pendulum!(x), dim: 6, x(0) ∈ X₀_small × U₀)

period = 0.02
k = 20
T = k * period
T_warmup = 2 * period

prob = ControlledPlant(ivp, controller, vars_idx, period)
prob_small = ControlledPlant(ivp_small, controller, vars_idx, period)

# simulation

println("simulation")
res = @timed simulate(prob, T=T, trajectories=10, include_vertices=true)
sim = res.value
print_timed(res)

# reachability analysis

alg = TMJets20(abstol=1e-9, orderT=8, orderQ=1)
alg_nn = Ai2()

function benchmark(; T=T, silent::Bool=false)
    silent || println("flowpipe construction")
    res = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res.value
    silent || print_timed(res)

    silent || println("flowpipe construction (small initial states)")
    res = @timed solve(prob_small, T=T, alg_nn=alg_nn, alg=alg)
    sol_small = res.value
    silent || print_timed(res)

    return sol, sol_small
end

benchmark(T=T_warmup, silent=true)
res = @timed benchmark()
sol, sol_small = res.value

# plotting

vars = (1, 2)
fig = plot(xlab=L"\theta", ylab=L"\theta'", leg=:topleft,
		   legendfontsize=5,
		   tickfont=font(25, "Times"),
		   guidefontsize=25,
		   xguidefont=font(25, "Times"),
		   yguidefont=font(25, "Times"),
		   xticks=([1, 1.25, 1.5, 1.75, 2], ["1.0", "1.25", "1.5", "1.75", "2.0"]),
		   bottom_margin=0mm, left_margin=0mm, right_margin=4mm, top_margin=0mm, size=(500, 450))
plot!(fig, sol, vars=vars, color=:orange, lab="")
plot!(project(X₀, vars), c=:blue, lab=L"X_0", alpha=0.6)
plot!(fig, sol_small, vars=vars, color=:red, lab="")
plot!(project(X₀_small, vars), c=:cyan, lab=L"X_0'", alpha=1)
plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
savefig("doublependulum_more_robust_x1x2_juliareach.png")

end
nothing


text:
#=

prec=precision(midpt); # bits of significand (signficant radix-2 digits
prec, typeof(midpt)
# (57,ArbFloats.ArbFloat{57})

1 + round(Int,log2(ufp2(midpt))) -  round(Int,log2(ulp2(midpt)))
# 57

prec == (1 + round(Int,log2(ufp2(midpt))) -  round(Int,log2(ulp2(midpt))))
# true

=#

"""
   x.midpoint -> (significand, exponentOf2)
                  [0.5,1.0)     2^expo
   x.radius   -> (radial significand, radial exponentOf2)
"""
function frexp(x::ArfFloat{P}) where {P}
    exponentOf2 = x.exponentOf2
    significandOf = deepcopy(x)
    significandOf.exponentOf2 = 0
    return significandOf, exponentOf2
end

function ldexp(s::ArfFloat{P}, e::Int) where {P}
    z = deepcopy(s)
    z.exponentOf2 = e
    return z
end

function frexp(x::ArbFloat{P}) where {P}
    significandOf, exponentOf2 = frexp(ArfFloat{P}(x))
    return ArbFloat{P}(significandOf), exponentOf2
end

function ldexp(s::ArbFloat{P}, e::Int) where {P}
    z = deepcopy(s)
    z.exponentOf2 = e
    return z
end
ldexp(x::Tuple{ArbFloats.ArbFloat{P}, Int}) where {P} = ldexp(x[1],x[2])

#=

      gte_bits2digs(b_bits)    ~>  d_digits
      d_digits | d digits suffice to encode and recover b bits without error

      gte_digs2bits(d_digits)  ~>  b_bits
      b_bits |  b bits suffice to encode and recover d digits without error

      lte_bits2digs(b_bits)    ~>  c_digits
      c_digits | b bits are necessary encode and recover c digits without error

      lte_digs2bits(d_digits)  ~>  a_bits
      a_bits |  d digits are necessary encode and recover a bits without error

      lte_digs2bits( gte_bits2digs(b_bits  ) ) == b_bits
      lte_bits2digs( gte_digs2bits(d_digits) ) == d_digits

const log2_log10 = log(10,2)  # 0.3010299956639812 ~= log(2)/log(10)
const log10_log2 = log(2,10)  # 3.321928094887362  ~= log(10)/log(2)

lte_bits2digs(nbits::Int) = floor(Int, nbits * log2_log10)
lte_digs2bits(ndigs::Int) = floor(Int, ndigs * log10_log2)
=#

lte_bits2digs(nbits::Int) = floor( Int, nbits * 0.30102999566398125 ) # log(10,2) RoundUp
gte_bits2digs(nbits::Int)  = ceil( Int, nbits * 0.3010299956639812  ) # log(10,2) RoundDown
lte_digs2bits(ndigs::Int) = floor( Int, ndigs * 3.3219280948873626  ) # log(2,10) RoundUp
gte_digs2bits(ndigs::Int)  = ceil( Int, ndigs * 3.3219280948873617  ) # log(2,10) RoundDown

"""
logarithm_base(x)
"""
function log_base(x::Real, base::Int)
   z = if base == 2
           log2(x)
        elseif base == 10
           log10(x)
        else
           log(x) / log(base)
        end
   return z
end
log_base(x::ArbFloat{P}, base::Int) where {P} = ArbFloats.logbase(x,base)

#=
    position_first_place: the radix position of the most significant nonzero bit|digit

    pfp{T<:AbstractFloat}(x::T, base::Int=2)
    pfp{T<:AbstractFloat}(x::T)              == pfp(x,  2)  ==  pfp2{T<:AbstractFloat}(x::T)
    pfp{T<:AbstractFloat}(x::T, base=10)     == pfp(x, 10)  ==  pfp10{T<:AbstractFloat}(x::T)

    position_last_place: the radix position of the least significant nonzero bit|digit

    plp{T<:AbstractFloat}(x::T, base::Int=2)
    plp{T<:AbstractFloat}(x::T)              == p;p(x,  2)  ==  plp2{T<:AbstractFloat}(x::T)
    plp{T<:AbstractFloat}(x::T, base=10)     == plp(x, 10)  ==  p;p10{T<:AbstractFloat}(x::T)

    unit_first_place: the radix *value) of the most significant nonzero bit|digit

    pfp{T<:AbstractFloat}(x::T, base::Int=2)
    pfp{T<:AbstractFloat}(x::T)              == pfp(x,  2)  ==  pfp2{T<:AbstractFloat}(x::T)
    pfp{T<:AbstractFloat}(x::T, base=10)     == pfp(x, 10)  ==  pfp10{T<:AbstractFloat}(x::T)

    unit_last_place: the radix *value* of the least significant nonzero bit|digit

    plp{T<:AbstractFloat}(x::T, base::Int=2)
    plp{T<:AbstractFloat}(x::T)              == plp(x,  2)  ==  plp2{T<:AbstractFloat}(x::T)
    plp{T<:AbstractFloat}(x::T, base=10)     == plp(x, 10)  ==  plp10{T<:AbstractFloat}(x::T)

=#

"""
position_first_place
determine the position of the most significant nonzero bit|digit
"""
function pfp(x::T, base::Int=2) where {T <: Real}
   z = 0 # if x==0.0
   if notzero(x)
       z = floor( Int, log_base(abs(x), base) )
   end
   return z
end
pfp(x::ArbFloat{P}, base::Int=2) where {P} =
    return iszero(x) ? 0 : floor( Int, log_base(abs(smartvalue(x)), base) )
"""
binary position_first_place
determine the position of the most significant nonzero bit
"""
pfp2(x::T) where {T <: Real} = (x==zero(T) ? 0 : floor( Int, log2(abs(x)) ))
pfp2(x::ArbFloat{P}) where {P} = (iszero(x) ? 0 : floor( Int, log2(abs(smartvalue(x))) ))
"""
decimal position_first_place
determine the position of the most significant nonzero digit
"""
pfp10(x::T) where {T <: Real} = (x==zero(T) ? 0 : floor( Int, log10(abs(x)) ))
pfp10(x::ArbFloat{P}) where {P} = (iszero(x) ? 0 : floor( Int, log10(abs(smartvalue(x))) ))

"""
ufp is unit_first_place
the float value given by a 1 at the position of
  the most significant nonzero bit|digit in _x_
"""
function ufp(x::AbstractFloat, base::Int=2)
   z = pfp(x, base)
   b = convert(Float64, base)
   return b^z
end
function ufp(x::ArbFloat{P}, base::Int=2) where {P}
   z = pfp(x, base)
   return Float64(base)^z
end
ufp(x::Integer, base::Int=2) = ufp(Float64(x), base)
"ufp2 is unit_first_place in base 2"
ufp2(x::T) where {T <: Real} = 2.0^pfp2(x)
ufp2(x::ArbFloat{P}) where {P} = 2.0^pfp2(x)
ufp2(x::Integer) = ufp2(Float64(x))
"ufp10 is unit_first_place in base 10"
ufp10(x::T) where {T <: Real} = 10.0^pfp10(x)
ufp10(x::ArbFloat{P}) where {P} = 10.0^pfp10(x)
ufp10(x::Integer) = ufp10(Float64(x))
"""
ulp   is unit_last_place
the float value given by a 1 at the position of
  the least significant nonzero bit|digit in _x_
"""
function ulp(x::Real, precision::Int, base::Int)
   unitfp  = base==2 ? ufp2(x) : (base==10 ? ufp10(x) : throw(DomainError()))
   twice_u = 2.0^(1-precision)
   return twice_u * unitfp
end
ulp(x::T, base::Int=2) where {T <: AbstractFloat}  =
    ulp(x, 1+Base.significand_bits(T), base)
ulp(x::ArbFloat{P}, base::Int=2) where {P}  =
    ulp(x, P, base)
ulp(x::Integer, base::Int=2) = ulp(Float64(x), base)

"""ulp2  is unit_last_place base 2"""
function ulp2(x::Real, precision::Int)
   unitfp  = ufp2(x)
   twice_u = 2.0^(1-precision)
   return (*)(promote(twice_u, unitfp)...)
end
function ulp2(x::ArbFloat{P}) where {P}
   unitfp  = ufp2(x)
   twice_u = 2.0^(1-P)
   return (*)(promote(twice_u, unitfp)...)
end
ulp2(x::T) where {T <: AbstractFloat}  = ulp2(x, 1+Base.significand_bits(T))
ulp2(x::Integer) = ulp2(Float64(x))

"""ulp10 is unit_last_place base 10"""
function ulp10(x::Real, bitprecision::Int)
    unit10fp = ufp10(x)
    digitprecision = lte_bits2digs(bitprecision)
    twice_u10 = 10.0^(1-digitprecision)
    return twice_u10 * unit10fp
end
function ulp10(x::ArbFloat{P}) where {P}
    unit10fp = ufp10(x)
    digitprecision = lte_bits2digs(P)
    twice_u10 = 10.0^(1-digitprecision)
    return twice_u10 * unit10fp
end
ulp10(x::T) where {T <: AbstractFloat} = ulp10( x, (1+Base.significand_bits(T)) )
ulp10(x::Integer) = ulp10(Float64(x))

function eps(::Type{T}) where {T <: ArbFloat}
    P = precision(T)
    return ldexp(one(T), 1-P)
end
function eps(x::ArbFloat{P}) where {P}
    a = eps(ArbFloat{P})
    return midpoint(abs(x)*a)
end

# !!revisit!!
function nextfloat(x::ArbFloat{P}) where {P}
    m,r = midpoint_radius(x)
    e = eps(m)
    n = m
    i = 1
    while n == m
      n = midpoint(n + i*e)
      i += 1
    end
    return midpoint_radius(n, r)
end

# !!revisit!!
function prevfloat(x::ArbFloat{P}) where {P}
    m,r = midpoint_radius(x)
    e = eps(m)
    n = m
    i = 1
    while n == m
      n = midpoint(n - i*e)
      i += 1
    end
    return midpoint_radius(n, r)
end


text:
# 宣言方法がちょっとだけ違う
fruits = Set(["apple","banana","orange"])
println(fruits)

# 型を明示的に指定できる
fruits = Set{String}(["apple","banana","orange"])
println(fruits)

# 空集合
fruits_empty = Set()

# 要素の追加

fruits = Set{String}()
push!(fruits,"apple")
println(fruits)

# 型がおかしいからエラー
# push!(fruits,1)
# println(fruits)

# 存在があるかどうかをチェック

fruits = Set(["apple","orange"])
println(in("apple",fruits))
println(in("banana",fruits))

# 集合の演算

fruits = Set(["apple","banana","orange"])
another_fruits = Set(["peach","pineapple"])
drink = Set(["coffee","milk","cocoa"])
number = Set([1,2,3,4])

# Union
println(union(fruits,drink))

println(union(fruits,number))

# intersection
println(intersect(fruits, another_fruits) )

# difference
println(setdiff(fruits,another_fruits))
println(setdiff(another_fruits,fruits))


text:
using BinDeps
using Compat.Sys: BINDIR, isapple, iswindows

@BinDeps.setup

ENV["JULIA_ROOT"] = abspath(BINDIR, "../../")

# include alias for WinRPM library
libsndfile = library_dependency("libsndfile", aliases=["libsndfile-1"])

provides(AptGet, "libsndfile1-dev", libsndfile)
provides(Pacman, "libsndfile", libsndfile)

@static if isapple()
    using Homebrew
    provides(Homebrew.HB, "libsndfile", libsndfile)
end

@static if iswindows()
    using WinRPM
    provides(WinRPM.RPM, "libsndfile1", libsndfile, os = :Windows)
end

@BinDeps.install Dict(:libsndfile => :libsndfile)


text:
\#-----------------------------------------------------------------------# Partition
"""
    Partition(stat)
    Partition(stat, nparts)

Split a data stream into `nparts` (default 100) where each part is summarized by `stat`.

# Example

    o = Partition(Extrema())
    fit!(o, cumsum(randn(10^5)))

    using Plots
    plot(o)
"""
mutable struct Partition{T, I, O <: OnlineStat{T}} <: OnlineStat{T}
    parts::Vector{Pair{Tuple{Int,Int}, O}} # (a,b) => stat
    b::Int
    init::I
    n::Int
end
Partition(init::Function, b::Int=100) = Partition([(1,1) => init()], b, init, 0)
Partition(o::OnlineStat, b::Int=100) = Partition(() -> copy(o), b)

OnlineStatsBase.additional_info(o::Partition) = (;b = o.b)
value(o::Partition) = (o.parts[1][1][1], o.parts[end][1][2]) => value(reduce(merge!, last.(o.parts), init=o.init()))

function _fit!(o::Partition, y)
    n = (o.n += 1)
    parts = o.parts
    lastpart = parts[end]
    (a, b), stat = lastpart
    if a ≤ n ≤ b
        _fit!(stat, y)
    else
        lastpart = (n, n + nobs(stat) - 1) => fit!(o.init(), y)
        push!(parts, lastpart)
    end
    length(parts) > o.b && nobs(lastpart[2]) == (b - a + 1) && merge_next!(parts)
end

function merge_next!(parts)
    n = nobs(parts[1][2])
    i = 1
    for (j, p) in enumerate(parts)
        nobs(p[2]) < n && (i = j; break;)
    end
    # a, b = parts[i], parts[i+1]
    ((a1, b1), stat1), ((a2,b2), stat2) = parts[i], parts[i+1]
    parts[i] = (min(a1, a2), max(b1, b2)) => merge!(stat1, stat2)
    deleteat!(parts, i + 1)
end

# Assumes `a` goes before `b`
function _merge!(a::Partition, b::Partition)
    n = nobs(a)
    a.n += b.n
    for p in b.parts
        push!(a.parts, (p[1][1] + n, p[1][2] + n) => copy(p[2]))
    end
    while length(a.parts) > a.b
        merge_next!(a.parts)
    end
    a
end

#-----------------------------------------------------------------------# IndexedPartition
"""
    IndexedPartition(T, stat, b=100)

Summarize data with `stat` over a partition of size `b` where the data is indexed by a
variable of type `T`.

# Example

    x, y = randn(10^5), randn(10^6)
    o = IndexedPartition(Float64, KHist(10))
    fit!(o, zip(x, y))

    using Plots
    plot(o)
"""
mutable struct IndexedPartition{I, T, O <: OnlineStat{T}, P <: Pair{<:Tuple, O}} <: OnlineStat{TwoThings{I,T}}
    parts::Vector{P}
    b::Int
    init::O
    n::Int
end
function IndexedPartition(I::Type, o::O, b::Int=100) where {T, O<:OnlineStat{T}}
    IndexedPartition{I,T, O,Pair{Tuple{I,I}, O}}(Pair{Tuple{I,I}, O}[], b, o, 0)
end
value(o::IndexedPartition) = nobs(o) == 0 ?
    nothing => value(copy(o.init)) :
    (o.parts[1][1][1], o.parts[end][1][2]) => value(reduce(merge!, last.(o.parts), init=copy(o.init)))

function _fit!(o::IndexedPartition{I}, xy) where {I}
    x, y = xy
    n = o.n += 1
    isempty(o.parts) && push!(o.parts, (I(x), I(x)) => copy(o.init))
    i = findfirst(p -> p[1][1] ≤ x ≤ p[1][2], o.parts)
    if isnothing(i)
        push!(o.parts, (I(x), I(x)) => fit!(copy(o.init), y))
    else
        _fit!(o.parts[i][2], y)
    end
    length(o.parts) > o.b && indexed_merge_next!(sort!(o.parts, by = x -> x[1][1]))
end

function indexed_merge_next!(parts)
    diffs = map(neighbors(parts)) do (a, b)
        (b[1][1] - a[1][2]) * round(Int, middle(nobs(a[2]), nobs(b[2])))
    end
    _, i = findmin(diffs)
    parts[i] = (parts[i][1][1], parts[i+1][1][2]) => merge!(parts[i][2], parts[i + 1][2])
    deleteat!(parts, i + 1)
end

function _merge!(p1::IndexedPartition, p2::IndexedPartition)
    sort!(append!(p1.parts, p2.parts), by = x -> x[1][1])
    p1.n += p2.n
    while length(p1.parts) > p1.b
        indexed_merge_next!(p1.parts)
    end
end

#-----------------------------------------------------------------------------# KIndexedPartition
"""
    KIndexedPartition(T, stat_init, k=100)

Similar to [`IndexedPartition`](@ref), but indexes the first variable by centroids (like [`KHist`](@ref))
rather than intervals.

- Note: `stat_init` must be a function e.g. `() -> Mean()`

# Example
    using Plots

    o = KIndexedPartition(Float64, () -> KHist(10))

    fit!(o, zip(randn(10^6), randn(10^6)))

    plot(o)
"""
mutable struct KIndexedPartition{I,T,O<:OnlineStat{T},F} <: OnlineStat{TwoThings{I,T}}
    parts::Vector{Pair{I,O}}
    k::Int
    init::F
    n::Int
end
function KIndexedPartition(I::Type, init::Base.Callable, k::Int=100)
    o = init()
    T, O = OnlineStatsBase.input(o), typeof(o)
    KIndexedPartition{I,T,O,typeof(init)}(Pair{I,O}[], k, init, 0)
end
value(o::KIndexedPartition) =
    (o.parts[1][1][1], o.parts[end][1][2]) => value(reduce(merge!, last.(o.parts), init=o.init()))

function _fit!(o::KIndexedPartition, xy)
    x, y = xy
    parts = o.parts
    newpart = x => fit!(o.init(), y)
    insert!(parts, searchsortedfirst(parts, newpart; by=first), newpart)
    length(parts) > o.k && kindexed_merge_next!(parts)
end

function kindexed_merge_next!(parts)
    mindiff = Inf
    i = 0
    for (j, (a,b)) in enumerate(neighbors(parts))
        d = first(b) - first(a)
        if d < mindiff && 1 < j < (length(parts) - 1)
            mindiff = d
            i = j
        end
    end
    a = parts[i]
    b = parts[i + 1]
    n = nobs(last(a)) + nobs(last(b))
    parts[i] = smooth(first(a), first(b), nobs(last(b)) / n) => merge!(a[2], b[2])
    deleteat!(parts, i + 1)
end

function _merge!(a::KIndexedPartition, b::KIndexedPartition)
    sort!(append!(a.parts, b.parts), by=first)
    while length(a.parts) > a.k
        kindexed_merge_next!(a.parts)
    end
end


text:
to_dash(t::Any) = t

struct Component
    name ::String
    type ::String
    namespace ::String
    props ::Dict{Symbol, Any}
    available_props ::Set{Symbol}
    wildcard_regex ::Union{Nothing, Regex}
    function Component(name::String, type::String, namespace::String,
                            props::Vector{Symbol}, wildcard_props::Vector{Symbol}; kwargs...)

        available_props = Set{Symbol}(props)
        wildcard_regex::Union{Nothing, Regex} = nothing
        if !isempty(wildcard_props)
            wildcard_regex = Regex(join(string.(wildcard_props), "|"))
        end
        component = new(name, type, namespace, Dict{Symbol, Any}(), available_props, wildcard_regex)
        for (prop, value) in kwargs
            Base.setproperty!(component, prop, value)
        end
        return component
    end
end

get_name(comp::Component) = getfield(comp, :name)
get_type(comp::Component) = getfield(comp, :type)
get_namespace(comp::Component) = getfield(comp, :namespace)
get_available_props(comp::Component) = getfield(comp, :available_props)
get_wildcard_regex(comp::Component) = getfield(comp, :wildcard_regex)
get_props(comp::Component) = getfield(comp, :props)

function Base.getproperty(comp::Component, prop::Symbol)
    !Base.hasproperty(comp, prop) && error("Component $(get_name(comp)) has no property $(prop)")
    props = get_props(comp)
    return haskey(props, prop) ? props[prop] : nothing
end

function Base.setproperty!(comp::Component, prop::Symbol, value)
    !Base.hasproperty(comp, prop) && error("Component $(get_name(comp)) has no property $(prop)")
    props = get_props(comp)
    push!(props, prop=>to_dash(value))
end

function check_whild(wildcard_regex::Union{Nothing, Regex}, name::Symbol)
    isnothing(wildcard_regex) && return false
    return startswith(string(name), wildcard_regex)
end

function Base.hasproperty(comp::Component, prop::Symbol)
    return in(prop, get_available_props(comp)) || check_whild(get_wildcard_regex(comp), prop)
end

Base.propertynames(comp::Component) = collect(get_available_props(comp))

push_prop!(component::Component, prop::Symbol, value) = push!(component.props, prop=>to_dash(value))

JSON2.@format Component begin
    name => (exclude = true,)
    available_props => (exclude = true,)
    wildcard_regex => (exclude = true,)
end

max_stars_repo_name: waralex/DashBase

max_issues_repo_name: waralex/DashBase

max_forks_repo_name: waralex/DashBase


text: include("StandardForm.jl")

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
function ConstantThrustManeuver(arg0::AbsoluteDate, arg1::jdouble, arg2::AttitudeProvider, arg3::AbstractConstantThrustPropulsionModel)
    return ConstantThrustManeuver((AbsoluteDate, jdouble, AttitudeProvider, AbstractConstantThrustPropulsionModel), arg0, arg1, arg2, arg3)
end

function ConstantThrustManeuver(arg0::AbsoluteDate, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::AttitudeProvider, arg5::Vector3D)
    return ConstantThrustManeuver((AbsoluteDate, jdouble, jdouble, jdouble, AttitudeProvider, Vector3D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function ConstantThrustManeuver(arg0::AbsoluteDate, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::AttitudeProvider, arg5::Vector3D, arg6::JString)
    return ConstantThrustManeuver((AbsoluteDate, jdouble, jdouble, jdouble, AttitudeProvider, Vector3D, JString), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function ConstantThrustManeuver(arg0::AbsoluteDate, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::Vector3D)
    return ConstantThrustManeuver((AbsoluteDate, jdouble, jdouble, jdouble, Vector3D), arg0, arg1, arg2, arg3, arg4)
end

function ConstantThrustManeuver(arg0::AbsoluteDate, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::Vector3D, arg5::JString)
    return ConstantThrustManeuver((AbsoluteDate, jdouble, jdouble, jdouble, Vector3D, JString), arg0, arg1, arg2, arg3, arg4, arg5)
end

function ConstantThrustManeuver(arg0::AttitudeProvider, arg1::DateBasedManeuverTriggers, arg2::AbstractConstantThrustPropulsionModel)
    return ConstantThrustManeuver((AttitudeProvider, DateBasedManeuverTriggers, AbstractConstantThrustPropulsionModel), arg0, arg1, arg2)
end

function get_direction(obj::ConstantThrustManeuver)
    return jcall(obj, "getDirection", Vector3D, ())
end

function get_duration(obj::ConstantThrustManeuver)
    return jcall(obj, "getDuration", jdouble, ())
end

function get_end_date(obj::ConstantThrustManeuver)
    return jcall(obj, "getEndDate", AbsoluteDate, ())
end

function get_flow_rate(obj::ConstantThrustManeuver)
    return jcall(obj, "getFlowRate", jdouble, ())
end

function get_isp(obj::ConstantThrustManeuver)
    return jcall(obj, "getISP", jdouble, ())
end

function get_start_date(obj::ConstantThrustManeuver)
    return jcall(obj, "getStartDate", AbsoluteDate, ())
end

function get_thrust(obj::ConstantThrustManeuver)
    return jcall(obj, "getThrust", jdouble, ())
end

function get_thrust_vector(obj::ConstantThrustManeuver)
    return jcall(obj, "getThrustVector", Vector3D, ())
end

function is_firing(obj::ConstantThrustManeuver, arg0::AbsoluteDate)
    return jcall(obj, "isFiring", jboolean, (AbsoluteDate,), arg0)
end

function is_firing(obj::ConstantThrustManeuver, arg0::FieldSpacecraftState)
    return jcall(obj, "isFiring", jboolean, (FieldSpacecraftState,), arg0)
end

function is_firing(obj::ConstantThrustManeuver, arg0::SpacecraftState)
    return jcall(obj, "isFiring", jboolean, (SpacecraftState,), arg0)
end


text:
struct STRtree{T}
    ptr::Ptr{GEOSSTRtree}  # void pointer to the tree
    items::T  # any geometry for which an envelope can be derived
end

"""
    STRtree(items; nodecapacity=10, context::GEOSContext=_context)

Create a Sort Tile Recursive tree for fast intersection queries of objects for which an
envelope is defined. The tree cannot be changed after its creation.

## Arguments
- `items`: The items to store in the tree.
- `nodecapacity`: The maximum number of items that can be stored per tree node (default 10).
- `context`: The GEOS context in which the tree should be created in. Defaults to the global
  LibGEOS context.
"""
function STRtree(items; nodecapacity = 10, context::GEOSContext = _context)
    tree = LibGEOS.GEOSSTRtree_create_r(context.ptr, nodecapacity)
    for item in items
        envptr = envelope(item).ptr
        GEOSSTRtree_insert_r(context.ptr, tree, envptr, pointer_from_objref(item))
    end
    return STRtree(tree, items)
end

"""
    query(tree::STRtree, geometry; context::GEOSContext=_context)

Returns the objects within `tree`, whose envelope intersects the envelope of `geometry`.

## Arguments
- `tree`: The STRtree to query
- `geometry`: The LibGEOS geometry (e.g. Polygon) to run the query for
- `context`: The GEOS context. Defaults to the global LibGEOS context.
"""
function query(
    tree::STRtree{T},
    geometry::Geometry;
    context::GEOSContext = _context,
) where {T}
    matches = eltype(T)[]

    # called for each overlapping item in the tree, used to push the item to matches
    function callback(item, userdata)::Ptr{Cvoid}
        item_ = unsafe_pointer_to_objref(item)
        userdata_ = unsafe_pointer_to_objref(userdata)
        push!(userdata_, item_)
        return C_NULL
    end

    cfun_callback = @cfunction($callback, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}))
    GEOSSTRtree_query_r(
        context.ptr,
        tree.ptr,
        geometry.ptr,
        cfun_callback,
        pointer_from_objref(matches),
    )
    return matches
end