text:
using CellListMap
using StaticArrays
using Test

# Loads Unitful and ForwardDiff
include("../src/examples/generic_types.jl")

@testset "CellListMap.jl" begin

    if Threads.nthreads() == 1
        println("""

             WARNING: Ideally, run a multi-threaded test to check the parallel versions.

        """)
    end

    # Number of particles, sides and cutoff
    sides = @SVector [250.,250.,250.]
    cutoff = 10.
    box = Box(sides,cutoff)

    # Particle positions
    N = 2000
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    # Add some pathological coordinates
    x[1] = -sides/2
    x[2] = -sides/2 + @SVector [ nextfloat(0.), nextfloat(0.), sides[3]*rand() ]
    x[3] = -sides/2 + @SVector [ prevfloat(0.), prevfloat(0.), sides[3]*rand() ]
    x[4] = sides/2 + @SVector [ nextfloat(0.), nextfloat(0.), sides[3]*rand() ]
    x[5] = sides/2 + @SVector [ prevfloat(0.), prevfloat(0.), sides[3]*rand() ]
    x[10] = sides
    x[11] = sides + @SVector [ nextfloat(0.), nextfloat(0.), sides[3]*rand() ]
    x[12] = sides + @SVector [ prevfloat(0.), prevfloat(0.), sides[3]*rand() ]
    x[13] = @SVector [ nextfloat(0.), nextfloat(0.), sides[3]*rand() ]
    x[14] = @SVector [ prevfloat(0.), prevfloat(0.), sides[3]*rand() ]
    x[100] = @SVector [sides[1]/2, -sides[2]/2, 2*sides[3]]
    y = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]

    # Initialize auxiliary linked lists
    cl = CellList(x,box)

    # Function to be evalulated for each pair: sum of displacements on x
    f(x,y,avg_dx) = avg_dx + abs(x[1] - y[1])

    naive = CellListMap.map_naive!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,x,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=false) ≈ naive

    # Test updating the cell lists
    new_cl = deepcopy(cl)
    new_x = copy(x) .+ [rand(SVector{3,Float64}) for _ in 1:N ]
    new_sides = sides + rand(SVector{3,Float64})
    new_cutoff = cutoff + rand()
    new_box = Box(sides,cutoff)
    new_cl = CellListMap.UpdateCellList!(new_x,new_box,new_cl)
    new_naive = CellListMap.map_naive!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_x,new_box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_box,new_cl,parallel=true) ≈ new_naive
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_box,new_cl,parallel=false) ≈ new_naive
    new_box = Box([ 250   0  10
                     10 250   0
                      0   0 250 ],cutoff)
    new_cl = CellListMap.UpdateCellList!(new_x,new_box,new_cl)
    new_naive = CellListMap.map_naive!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_x,new_box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_box,new_cl,parallel=true) ≈ new_naive
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,new_box,new_cl,parallel=false) ≈ new_naive

    # Test the input as a matrix
    xmat = zeros(3,N)
    for i in 1:N
        for j in 1:3
            xmat[j,i] = x[i][j]
        end
    end
    cl_mat = CellList(xmat,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl_mat,parallel=true) ≈ naive
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl_mat,parallel=false) ≈ naive

    # Check if changing lcell breaks something
    box = Box(sides,cutoff,lcell=1); cl = CellList(x,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    box = Box(sides,cutoff,lcell=2); cl = CellList(x,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    box = Box(sides,cutoff,lcell=3); cl = CellList(x,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    box = Box(sides,cutoff,lcell=5); cl = CellList(x,box)
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive

    # Test if changing the number of batches breaks anything
    cl = CellList(x,box,nbatches=(3,5))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    cl = CellList(x,box,nbatches=(1,1))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    cl = CellList(x,box,nbatches=(1,7))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    cl = CellList(x,box,nbatches=(7,1))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    cl = CellList(x,box,nbatches=(13,17))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
    cl = CellList(x,box,nbatches=(4,16))
    @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive

    # Test random cells of all possible types
    for N in 2:3,
        M in rand(10:20),
        UnitCellType in [ TriclinicCell, OrthorhombicCell ],
        parallel in [ false, true ],
        lcell in 1:3
        @test CellListMap.check_random_cells(
            N,M,
            UnitCellType=UnitCellType,
            parallel=parallel,
            lcell=lcell,
            show_progress=false
        )[1]
    end

    # Function to be evalulated for each pair: build distance histogram
    function build_histogram!(d2,hist)
        d = sqrt(d2)
        ibin = floor(Int,d) + 1
        hist[ibin] += 1
        return hist
    end

    naive = CellListMap.map_naive!( (x,y,i,j,d2,hist) -> build_histogram!(d2,hist), zeros(Int,10),x,box)
    @test map_pairwise!( (x,y,i,j,d2,hist) -> build_histogram!(d2,hist), zeros(Int,10),box,cl,parallel=true) ≈ naive
    @test map_pairwise!( (x,y,i,j,d2,hist) -> build_histogram!(d2,hist), zeros(Int,10),box,cl,parallel=false) ≈ naive

    # Function to be evalulated for each pair: build distance histogram
    function potential(i,j,d2,u,mass)
        d = sqrt(d2)
        u = u - 9.8*mass[i]*mass[j]/d
        return u
    end

    # Run pairwise computation
    mass = rand(N)
    naive = CellListMap.map_naive!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,x,box)
    @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,box,cl,parallel=true) ≈ naive
    @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,box,cl,parallel=false) ≈ naive

    # Function to be evalulated for each pair: build distance histogram
    function calc_forces!(x,y,i,j,d2,mass,forces)
        G = 9.8*mass[i]*mass[j]/d2
        d = sqrt(d2)
        df = (G/d)*(x - y)
        forces[i] = forces[i] - df
        forces[j] = forces[j] + df
        return forces
    end

    # forces
    forces = [ zeros(SVector{3,Float64}) for i in 1:N ]

    # Run pairwise computation
    naive = CellListMap.map_naive!( (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces), copy(forces),x,box)
    @test map_pairwise!( (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces), copy(forces),box,cl,parallel=true) ≈ naive
    @test map_pairwise!( (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces), copy(forces),box,cl,parallel=false) ≈ naive

    #
    # Compute some properties of disjoint sets
    #
    box = Box(sides,cutoff,lcell=1)
    cl = CellList(x,y,box)

    naive = CellListMap.map_naive!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,x,y,box)
    @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=true) ≈ naive

    # Test disjoint sets, with matrices
    xmat = zeros(3,length(x))
    for i in 1:length(x)
        for j in 1:3
            xmat[j,i] = x[i][j]
        end
    end
    ymat = zeros(3,length(y))
    for i in 1:length(y)
        for j in 1:3
            ymat[j,i] = y[i][j]
        end
    end
    cl_mat = CellList(xmat,ymat,box)
    @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0, box, cl_mat, parallel=false) ≈ naive
    @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0, box, cl_mat, parallel=true) ≈ naive

    # Check different lcell
    box = Box(sides,cutoff,lcell=3)
    cl = CellList(x,y,box)
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=true) ≈ naive

    # Test if changing the number of batches breaks anything
    cl = CellList(x,y,box,nbatches=(1,1))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    cl = CellList(x,y,box,nbatches=(3,5))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    cl = CellList(x,y,box,nbatches=(7,1))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    cl = CellList(x,y,box,nbatches=(1,7))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    cl = CellList(x,y,box,nbatches=(4,16))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive
    cl = CellList(x,y,box,nbatches=(13,17))
    @test map_pairwise!( (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass), 0.0,box,cl,parallel=false) ≈ naive

    # Test the examples, to check further if the parallelization didn't break something
    N = 100_000
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    y = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    @test CellListMap.Examples.average_displacement(parallel=true) ≈ CellListMap.Examples.average_displacement(parallel=false)[1]
    @test CellListMap.Examples.distance_histogram(parallel=true,x=x) ≈ CellListMap.Examples.distance_histogram(parallel=false,x=x)
    @test CellListMap.Examples.gravitational_potential(parallel=true,x=x) ≈ CellListMap.Examples.gravitational_potential(parallel=false,x=x)
    @test CellListMap.Examples.gravitational_force(parallel=true,x=x) ≈ CellListMap.Examples.gravitational_force(parallel=false,x=x)
    @test count(CellListMap.Examples.nearest_neighbor(parallel=true,x=x,y=y) .≈ CellListMap.Examples.nearest_neighbor(parallel=false,x=x,y=y)) == 3

    function pair_match(p1,p2)
        p1[3] ≈ p2[3] || return false
        p1[1] == p2[1] && p1[2] == p2[2] && return true
        p1[1] == p2[2] && p1[2] == p2[1] && return true
    end
    pairs1 = sort!(CellListMap.Examples.neighborlist(parallel=true,x=x),by=x->x[3])
    pairs2 = sort!(CellListMap.Examples.neighborlist(parallel=false,x=x),by=x->x[3])
    @test length(pairs1) == length(pairs2)
    @test count(pair_match.(pairs1,pairs2)) == length(pairs1)

    x = [ sides .* rand(SVector{3,Float64}) for i in 1:1_500 ]
    y = [ sides .* rand(SVector{3,Float64}) for i in 1:1_500_000 ]
    @test count(CellListMap.Examples.nearest_neighbor_nopbc(parallel=true,x=x,y=y) .≈ CellListMap.Examples.nearest_neighbor_nopbc(parallel=false,x=x,y=y)) == 3

    # invert x and y to test swap
    ixy = CellListMap.Examples.nearest_neighbor_nopbc(parallel=false,x=x,y=y)
    iyx = CellListMap.Examples.nearest_neighbor_nopbc(parallel=false,x=y,y=x)
    @test ( ixy[1] == iyx[2] && ixy[2] == iyx[1] && ixy[3] ≈ iyx[3] )
    ixy = CellListMap.Examples.nearest_neighbor_nopbc(parallel=true,x=x,y=y)
    iyx = CellListMap.Examples.nearest_neighbor_nopbc(parallel=true,x=y,y=x)
    @test ( ixy[1] == iyx[2] && ixy[2] == iyx[1] && ixy[3] ≈ iyx[3] )

    # Test some fractional box lengths with the packmol test
    @test CellListMap.Examples.packmol(parallel=false, sides=[46.4,32.1,44.7], tol=3.14, UnitCellType=TriclinicCell)[1]
    @test CellListMap.Examples.packmol(parallel=true, sides=[18.4,30.1,44], tol=2, UnitCellType=TriclinicCell)[1]
    @test CellListMap.Examples.packmol(parallel=false, sides=[46.4,32.1,44.7], tol=3.14, UnitCellType=OrthorhombicCell)[1]
    @test CellListMap.Examples.packmol(parallel=true, sides=[18.4,30.1,44], tol=2, UnitCellType=OrthorhombicCell)[1]

    # Testing the propagation of types in automatic differentiation
    @test generic_types(false) == (true,u"nm^2",Measurement{Float64})

    # Test when we have pathologically few number of particles
    x = [ Float64[1,1,1] ]
    y = [ Float64[1.05,1,1], Float64[0,0,0]  ]
    @test CellListMap.neighborlist(x,y,0.1)[1] == (1, 1, 0.050000000000000044)
    z = [ Float64[1,1,1], Float64[1.05,1,1], Float64[0,0,0]  ]
    @test CellListMap.neighborlist(z,0.1)[1] == (1, 2, 0.050000000000000044)

end

include("./namd/compare_with_namd.jl")
include("./neighborlist.jl")


text:
function solve_counterfactual_tvf(ap, mp)

    ################################
    ################################
    # required data inputs
    ################################
    ################################

    ### BASELINE ECONOMY TIME PATHS
    # μ_baseline: baseline migration shares ((J+1)Rx(J+1)RxT)
    # xbilat_baseline: baseline bilateral expenditures (JNxNxT)
    # wages_baseline: baseline change in wages (JxNxT)
    # λ_baseline: baseline change in expenditure shares (JNxNxT)
    # labor_baseline: baseline change in labor for mobile regions ((J+1)xRxT)
    # L0: baseline initial labor for mobile regions [actual data] ((J+1)Rx1)

    ### INITIAL VALUES AND PARAMETERS
    # VARjn00: structures value added (JxN)
    # VALjn00: labor value added (JxN)
    # λ00: initial expenditure shares (JNxN)
    # xbilat00: initial bilateral expenditures (JNxN)
    # Sn00: trade imbalances not explained by global portfolio (Nx1)
    # α: final consumption shares (JxN)
    # io: share of each region in global portfolio (Nx1)
    # Ψ: share of structures in structure-labor aggregate (JxN)
    # G: share of materials in gross output (JNxN)
    # γ: share of value added in gross output (JxN)
    # θinv: inverse of trade elasticity (Jx1)

    ### INITIAL GUESSES
    # H: initial guesses of udots ((J+1)xRxT)
    # ωs: initial guesses of factor prices (JxNxT)

    ### COUNTERFACTUAL SHOCKS
    # Zs: counterfactual Zdot / baseline Zdot (JxNxT)

    ################################
    ################################
    # counterfactual shocks
    ################################
    ################################

    # counterfactual Zdot / baseline Zdot
    Zs = [ones(mp.J, mp.N) for i = 1:mp.T]
    china=[25.6097
            7.6093
            5.2773
            9.2464
            7.6897
            45.7781
            5.9298
            5.6533
            64.5661
            38.9064
            196.6509
            3.5324]
    china = china.^(1/28)
    for t = 1:28
        Zs[t][1:12, 57] = china
    end

    ################################
    ################################
    # load baseline and initial data
    ################################
    ################################
    #
    baseline_economy = load("data/Baseline_economy.jld")
    extract(baseline_economy)
    μ00 = baseline_economy["μ_baseline"][1]

    #baseline_economy = matread("../data/test-cdp/Baseline_economy.mat")
    counterfactual_economy = matread("data/Counterfactual_economy.mat")

    # Usual initial conditions
    base_year_full = matread("data/Base_year.mat")
    extract(base_year_full)

    # # initial guess for counterfactual trajectory of changes in utility udot = exp(Vt+1 - V)^1/ν
    Hc = [ones((mp.J + 1), mp.R) for i = 1:mp.T]
    Hc = [reshape(counterfactual_economy["V"][:,t], mp.J + 1, mp.R) for t = 1:mp.T]

    # set initial conditions
    VALjn0 = VALjn00 # labor value added
    VARjn0 = VARjn00 # structures value added
    λ = Din00 # bilateral trade shares
    Sn = Sn00 #trade imbalances not explained by global portfolio (they are zero by construction)
    Sn_prime = zeros(size(Sn)) # trade balance [exogenous] 2nd term on RHS of eq (14)
    xbilat = xbilat00 # bilateral trade flows
    κ_hat = ones(mp.J * mp.N, mp.N) # relative changes in trade costs [exogenous]

    io = base_year_full["io"] # region shares
    α = alphas # sector shares in final demand
    γ = gamma # share of value added in output
    Ψ = B # share of structures in structure-labor aggregate

    G = base_year_full["G"]' # share of materials (TRANSPOSE SO WE PULL COLUMNS NOT ROWS)
    G_vec = [G[:,n:n+mp.J-1] for n = 1:mp.J:size(G,2)] # vector form for quicker algebra

    θinv = T # inverse of theta
    θinv_stack = similar(λ[:,1]) # stack θinv JxJ submatrices for easier multiplication later
    for j in 1:mp.J
        index = 1+(j-1)*mp.N:mp.N*j
        θinv_stack[index,:] = ones(mp.N)*θinv[j]
    end
    θinv_stack = repeat(θinv_stack, 1, mp.N)

    # storage vectors
    wages_out = [ones(mp.J, mp.N) for i = 1:mp.T]
    λ_out = [similar(λ) for i = 1:mp.T]
    xbilat_out = [similar(VARjn0) for i = 1:mp.T]
    VARjn0_out = [similar(VARjn0) for i = 1:mp.T]
    VALjn0_out = [similar(VARjn0) for i = 1:mp.T]
    real_wages = [ones(mp.J, mp.N) for i = 1:mp.T]
    Ldyn = similar(labor_baseline)
    μ = [zeros(mp.R * (mp.J+1), mp.R * (mp.J + 1)) for t in 1:mp.T]

    wages_out[1] = ones(mp.J, mp.N)
    λ_out[1] = Din00
    xbilat_out[1] = xbilat_baseline[1]
    ωs = [wages_baseline[t+1].*labor_baseline[t+1].^Ψ for t in 1:mp.T-1]

    # initial phat guess
    phats = [ones(mp.J, mp.N) for i = 1:mp.T] # initial factor price guesses

    ################################
    ################################
    # outer loop
    ################################
    ################################

    iter = 1

    println("Starting algorithm.")

    while (iter <= ap.maxit) & (ap.Hmax > ap.tol_out)

        # reassign initial values
        VALjn0 = VALjn00 # labor value added
        VARjn0 = VARjn00 # structures value added
        λ = Din00 # bilateral trade shares

        Hcaux = reshape.(Hc, 1, (mp.J + 1) * mp.R)
        Hcaux = repeat.(Hcaux, (mp.J + 1) * mp.R, 1)

        # Step 2 in algorithm described in Appendix 3 Part II (create path of migration matrices)
        # Calculate mu1_tilde that refers the mobility matrix vartheta in the
        # appendix that has information on the jump in values
        μ = [zeros(mp.R * (mp.J+1), mp.R * (mp.J + 1)) for t in 1:mp.T]
        μ[1] = μ_baseline[2].*(Hcaux[2].^mp.beta);

        # Computing the mobility matrix at t=1
        num = μ[1].*(Hcaux[3].^mp.beta)
        num[isnan.(num)] .= 0
        den = sum(num', dims = 1)'
        μ[2] = num./den

        # And now for all the other time periods
        for t = 2:mp.T-2
            num = (μ_baseline[t+1]./μ_baseline[t]).*μ[t].*(Hcaux[t+2].^mp.beta)
            num[isnan.(num)] .= 0
           den = sum(num')'
           den = sum(num', dims = 1)'
           μ[t+1] = num./den
       end

        # using the mobility flows to construct the evolution of labor allocations in the baseline economy
        L00 = reshape(L0, mp.J+1, mp.R)
        L00_aux0 = reshape(L00, (mp.J+1)*mp.R, 1)
        L00_aux1 = repeat(L00_aux0, 1, (mp.J+1)*mp.R)
        L00_aux1 = μ00.*L00_aux1
        L1 = sum(L00_aux1, dims = 1)'
        L1 = reshape(L1, mp.J+1, mp.R)

        # compute full labor path
        Ldyn = similar(labor_baseline)
        Ldyn[2] = L1
        Ldyn[1] = reshape(L0, mp.J+1, mp.R)

        for t = 2:(mp.T-1)
            aux = μ[t] .* repeat(reshape(Ldyn[t], (mp.J+1) * mp.R, 1), 1, (mp.J+1) * mp.R)
            aux = sum(aux, dims = 1)'
            Ldyn[t+1] = reshape(aux, (mp.J+1), mp.R)
        end

        Ldyn[mp.T] .= 0
        Ltemp = [Ldyn[t][2:mp.J+1,:] for t = 1:mp.T]

        ################################
        # PART 2: given labor supply, solve for temporary
        # equilibrium (static production problem)
        ################################

        # Solve static subproblem at time t
        for t = 1:(mp.T-2)

            println("Counterfactual economy: solving period $(t) of $(mp.T).")

            Z_hat = 1 ./Zs[t] # change in productivity relative to baseline
            phat = phats[t]
            labor = ones(mp.J, mp.N) # [exogenous] change in employment in other countries
            labor[:,1:mp.R] = Ltemp[t+1] ./ Ltemp[t]

            # lambdas
            λ_1 = λ_baseline[t+1]
            λ_0 = λ_baseline[t]
            # set initial guess to be baseline economy value
            ω0 = wages_baseline[t+1].*labor_baseline[t+1].^Ψ
            ω = ωs[t]
            ω, wages, Phat, rental_rate, phat, VARjn_prime, VALjn_prime, X_prime, λ_prime, xbilat_prime =
                solve_temp_equilibrium(
                ω, ω0, phat, labor, VARjn0, VALjn0, λ, λ_0, λ_1, Sn_prime,
                κ_hat, Z_hat, α, io, θinv, θinv_stack, Ψ, G, G_vec, γ,
                mp, ap
            )

            # Get next period values
            λ = λ_prime
            VARjn0 = VARjn_prime
            VALjn0 = VALjn_prime

            ωs[t] = ω
            wages_out[t+1] = wages
            λ_out[t+1] = λ_prime
            VARjn0_out = VARjn_prime
            VALjn0_out = VALjn_prime
            xbilat_out[t+1] = xbilat_prime
            real_wages[t+1] = (wages./wages_baseline[t+1])./Phat
            phats[t] = phat

        end

        # From real wages, compute H
        real_wages_aux = [[1; (real_wages[t][:,1:mp.R]).^(1/mp.ν)] for t in 1:mp.T]
        real_wages_aux = reshape.(real_wages_aux, mp.R*(mp.J+1), 1) # 1st term on rhs of eq (17)

        # update changes in value at t = 1
        H0 = repeat(reshape(Hc[2], mp.R * (mp.J + 1), 1)', mp.R * (mp.J + 1), 1)
        μ1 = μ_baseline[2].*H0.^mp.beta
        aux = μ1 .* Hcaux[3].^mp.beta
        aux = sum(aux', dims = 1)'
        H0 = aux.*real_wages_aux[2]

        # update changes in value at other times t
        real_wages_aux = [[1; (real_wages[t][:,1:mp.R]).^(1/mp.ν)] for t in 1:mp.T]
        real_wages_aux = repeat.(reshape.(real_wages_aux, mp.R*(mp.J+1), 1), 1, mp.R*(mp.J+1)) # 1st term on rhs of eq (17)

        H_temp = [zeros(mp.R * (mp.J + 1), mp.R * (mp.J + 1)) for t in 1:mp.T]
        for t in 1:mp.T-1
            H_temp[t+1] = μ_baseline[t+1] ./ μ_baseline[t] .* μ[t] .* real_wages_aux[t+1]
            H_temp[t+1][isnan.(H_temp[t+1])] .= 0
        end

        num = [H_temp[t] .* Hcaux[t+1].^mp.beta for t in 1:mp.T-1] # RHS of eq (17) before summation
        push!(num, zeros(size(num[1])))

        Hc_new = [vec(sum(num[t]', dims = 1)) for t = 3:mp.T] # new path of y's
        pushfirst!(Hc_new, vec(H0))
        pushfirst!(Hc_new, vec(zeros(size(H0))))
        Hc_new[mp.T] .= 1.
        Hc[mp.T] .= 1
        Hc_old = reshape.(Hc, (mp.J+1)*mp.R, 1)

        # CDP error criterion
        last_error = ap.Hmax
        ap.Hmax = maximum([maximum(abs.(Hc_new[t] - Hc_old[t])) for t in 2:mp.T])
        println("Outer loop iteration $iter error: $(ap.Hmax), previous iteration error: $(last_error).")

        Hc = [0.5 * reshape(Hc_new[t], (mp.J+1), mp.R) + 0.5 * Hc[t] for t in 1:mp.T]

        # Tighten up tolerance once close, switch back to solving
        if ap.Hmax < ap.tol_out .* 2.0
            if ap.tol_in < 1e-7
                println("Decreasing tolerance.")
            end
            ap.tol_in = max(ap.tol_in ./ 10, 1e-7)
            howard_iter = ap.howard_its + 1
        end

        iter += 1

    end

    save("data/counterfactual_economy.jld",
    "wages", wages_out,
    "real_wages", real_wages,
    "λ", λ_out,
    "xbilat", xbilat_out,
    "Ldyn", Ldyn,
    "μ", μ,
    "Hc", Hc
    )

end

ext: jl

max_stars_repo_name: rafserqui/class-repo

max_issues_repo_name: rafserqui/class-repo

max_forks_repo_name: rafserqui/class-repo


text:
using CirrusCI
using Test

@test iscirrus()


text:
@testset "CDS_MVp_mat2CDS_tests" begin
  TI=Int32
  TF=Float32;
  n1=30;
  n2=20;
  N=n1*n2;
  x=randn(TF,N);
  y=zeros(TF,N);

  #test structured matrix
  #comp_grid = compgrid( (25, 25),( n1,n2 ) );
  comp_grid = compgrid( (TF(25), TF(25)),( n1,n2 ) );
  (TV_OP, AtA_diag, dense, TD_n, banded)=get_TD_operator(comp_grid,"TV",TF);
  A=TV_OP'*TV_OP;
  A=convert(SparseMatrixCSC{TF,TI},A);

  Ax_native = A*x;

  (R,offset)=mat2CDS(A)
  fill!(y,0)
  Ax_CDS = CDS_MVp(size(R,1),size(R,2),R,offset,x,y)
  @test isapprox(Ax_native,Ax_CDS,rtol=10*eps())

  #test random matrix
  x=randn(1000);
  y=zeros(1000);
  A=sprandn(1000,1000,0.1)
  Ax_native = A*x;

  (R,offset)=mat2CDS(A)
  fill!(y,0)
  Ax_CDS = CDS_MVp(size(R,1),size(R,2),R,offset,x,y)
  @test isapprox(Ax_native,Ax_CDS,rtol=10*eps())

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# This file was generated, do not modify it. # hide
sch = schema(X)
println(rpad(" Name", 28), "| Scitype")
println("-"^45)
for (name, scitype) in zip(sch.names, sch.scitypes)
    println(rpad("$name", 30), scitype)
end


text:
@testset "utils tests" begin
    @testset "knn_search" begin
        @inferred knn_search(rand(5, 100), 5, Euclidean())
        @inferred knn_search(rand(5, 100), 5, Euclidean(), Val(:approximate))
        @testset "pairwise tests" begin
            data = [0. 0. 0.; 0. 1.5 2.]
            true_knns = [2 3 2; 3 1 1]
            true_dists = [1.5 .5 .5; 2. 1.5 2.]
            knns, dists = knn_search(data, 2, Euclidean(), Val(:pairwise))
            @test knns == true_knns
            @test dists == true_dists
        end

        @testset "precomputed tests" begin
            dist_mat = [0. 2. 1.;
                        2. 0. 3.;
                        1. 3. 0.]
            true_knns = [3 1 1;
                         2 3 2]
            true_dists = [1 2 1;
                          2 3 3]
            knns, dists = knn_search(dist_mat, 2, :precomputed)
            @test knns == true_knns
            @test dists == true_dists

        end

    end

    @testset "combine_fuzzy_sets tests" begin
        A = [1.0 0.1; 0.4 1.0]
        union_res = [1.0 0.46; 0.46 1.0]
        res = combine_fuzzy_sets(A, 1.0)
        @test isapprox(res, union_res)
        inter_res = [1.0 0.04; 0.04 1.0]
        res = combine_fuzzy_sets(A, 0.0)
        @test isapprox(res, inter_res)

    end

    @testset "fit_ab tests" begin
        @test all((1, 2) .== fit_ab(-1, 0, 1, 2))
    end
end


text:
# no longer used +++++

function Base.show(io :: IO,  ws :: WorkerSet{F1}) where F1
    print(io,  "WorkerSet:",  ws.switches);
end

## ----------  Properties

ModelParams.has_pvector(wss :: WorkerSet) = false;

n_school(ws :: WorkerSet) = n_school(ws.switches);
n_school(wss :: WorkerSetSwitches) = wss.nSchool;

# Fixed utility while working
# Expressed as utility from fixed leisure
# Bad input argument +++
function u_fixed(ws :: WorkerSet{F1},  utilS,  s :: Integer) where F1
    if utilS.hasWorkLeisure
        u = UtilityFunctionsLH.CRRA(utilS.lCurvature, true);
        util = utilS.wtLeisure .* UtilityFunctionsLH.utility(u, utilS.leisure_sV[s]);
    else
        util = zero(F1);
    end
    return util
end

function wage_by_school(ws :: WorkerSet)
	return exp.(ModelParams.values(ws.wages))
end

function wage(ws :: WorkerSet, s :: Integer)
    return wage_by_school(ws)[s]
end

## ---------------  Constructors

"""
	$(SIGNATURES)

Constructs the worker set from switches and experience profile.
"""
function make_worker_set(objId :: ObjectId, wss :: WorkerSetSwitches, xpProfileV :: Vector{F1}) where F1 <: AbstractFloat
    # objId = make_child_id(model_id(),  :workerSet);
    # Assumes that all experience profiles are the same +++
    # xp = exper_profile(ds, HSG, T = wss.retireAge);
    wages = init_wages(make_child_id(objId, :wages), n_school(wss));

    # pvec = worker_set_pvector(objId, wss, uS);

    ws = WorkerSet(objId = objId, switches = wss,
        wages = wages, xp = xpProfileV[1 : wss.retireAge]);
    return ws
end

function make_test_worker_set(wss :: WorkerSetSwitches)
    return make_worker_set(ObjectId(:workerSet), wss,
        collect(range(0.0, 0.6, length = 70)))
end

# Log Wages by schooling (for workers)
    # change for wage object +++++
function init_wages(objId :: ObjectId, nSchool :: Integer)
    # This is a rough guess. The scale of human capital is on the order of 1
    # Starting wages of HSGs are only on the order of $10k; $20k for CGs (in Oksana's data).
    dataWageV = collect(range(10_000.0, 20_000.0, length = nSchool));
    # modelWageV = dollars_data_to_model(dataWageV, :perYear);
    modelWageV = dataWageV ./ 1000.0;  # +++
    logWageV = log.(modelWageV);
    ns = length(dataWageV);

    # Intercept and increments for `IncreasingVector`
    wageHSG = logWageV[1];
    wageInter = Param(:x0, "Log wage HSG", "wHSG",
        wageHSG, wageHSG, wageHSG - 1.0, wageHSG + 0.5,  true);

    # dWageV = diff(logWageV);
    dWageV = fill(0.05, ns - 1);
    dWage = Param(:dxV, "Log wage gradient", "dw",
        dWageV, dWageV, fill(0.01, ns-1), fill(0.3, ns-1), true);

    # ownId = make_child_id(objId, :wages);
    pvec = ParamVector(objId = objId,  pv = [wageInter, dWage]);
    return IncreasingVector(objId, pvec,  wageHSG, dWageV)
end

"""
	$(SIGNATURES)

Make a worker. `s` is the school level.
"""
function make_worker(ws :: WorkerSet{F1},  utilS,  R :: F1,  s :: Integer) where F1
    # Log utility is handled by setting `uCurvature = 1.0`
    u = UtilityFunctionsLH.CRRA(utilS.cCurvature, true);
    uFct = WorkerUtility(u,  utilS.discFactor,  u_fixed(ws, utilS, s));
    retireIncome = zero(F1); # give workers retirement income +++
    kMin, kMax = k_range(ws);
    hMin, hMax = h_range(ws);
    return Worker(uFct, ws.xp, wage(ws, s), R, ws.switches.retireAge,
        ws.switches.retireDuration, retireIncome,
        kMin, kMax, hMin, hMax)
end

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


text:
using CHCoordinateTransformations
using Test

@testset "ch_coordinates" begin

#testing the parameters according to page 19 of 20 of this PDF  https://www.swisstopo.admin.ch/content/swisstopo-internet/de/online/calculation-services/_jcr_content/contentPar/tabs/items/dokumente_und_publik/tabPar/downloadlist/downloadItems/8_1467103085694.download/refsys_d.pdf

#LV03, LN02
have=[602030.680 191775.030 897.915;
617306.300 268507.300 456.064;
776668.105 265372.681 1042.624;
497313.292 145625.438 1207.434;
722758.810 87649.670 1636.600]

#altitudes from step 3
altitudes2=[897.361;
457.138;
1043.616;
1206.367;
1634.472]

have_mod=hcat(have[:,1:2],altitudes2)

want=[    7 27 54.983506 46 52 37.540562 947.149;
7 40 6.983077 47 34 1.385301 504.935;
9 47 3.697723 47 30 55.172797 1089.372;
6 6 7.326361 46 27 14.690021 1258.274;
9 1 16.389053 45 55 45.438020 1685.027]

for i=1:size(want,1)
    thisDatum=have_mod[i,:]
    calculated_result=LV03toETRS89(thisDatum...)

    #check against expected values
    deltah=calculated_result[3]-want[i,7]
    deltacoords=vcat(compound(calculated_result[1]).-want[i,4:6],compound(calculated_result[2]).-want[i,1:3])
    deltas=vcat(deltah,deltacoords)
    @show deltas
    @test maximum(abs.(deltas))<0.05
    #1 second is roughly 30 meters (for Switzerland)
    #0.05 is roughly 1.5 meters
end

end


text:
using LazyArrays, Test

@testset "BroadcastArray" begin
    a = randn(6)
    b = BroadcastArray(exp, a)
    @test BroadcastArray(b) == BroadcastVector(b) == b

    @test b ==  Vector(b) == exp.(a)
    @test b[2:5] isa BroadcastVector
    @test b[2:5] == exp.(a[2:5])

    @test exp.(b) isa BroadcastVector
    @test b .+ SVector(1,2,3,4,5,6) isa BroadcastVector
    @test SVector(1,2,3,4,5,6) .+ b isa BroadcastVector

    A = randn(6,6)
    B = BroadcastArray(exp, A)

    @test Matrix(B) == exp.(A)

    C = BroadcastArray(+, A, 2)
    @test C == A .+ 2
    D = BroadcastArray(+, A, C)
    @test D == A + C

    @test sum(B) ≈ sum(exp, A)
    @test sum(C) ≈ sum(A .+ 2)
    @test prod(B) ≈ prod(exp, A)
    @test prod(C) ≈ prod(A .+ 2)

    x = Vcat([3,4], [1,1,1,1,1], 1:3)
    @test x .+ (1:10) isa Vcat
    @test (1:10) .+ x isa Vcat
    @test x + (1:10) isa Vcat
    @test (1:10) + x isa Vcat
    @test x .+ (1:10) == (1:10) .+ x == (1:10) + x == x + (1:10) == Vector(x) + (1:10)

    @test exp.(x) isa Vcat
    @test exp.(x) == exp.(Vector(x))
    @test x .+ 2 isa Vcat
    @test (x .+ 2).args[end] ≡ x.args[end] .+ 2 ≡ 3:5
    @test x .* 2 isa Vcat
    @test 2 .+ x isa Vcat
    @test 2 .* x isa Vcat

    A = Vcat([[1 2; 3 4]], [[4 5; 6 7]])
    @test A .+ Ref(I) == Ref(I) .+ A == Vcat([[2 2; 3 5]], [[5 5; 6 8]])

    @test_broken BroadcastArray(*,1.1,[1 2])[1] == 1.1

    B = BroadcastArray(*, Diagonal(randn(5)), randn(5,5))
    @test B == broadcast(*,B.args...)
    @test colsupport(B,1) == rowsupport(B,1) == 1:1
    @test colsupport(B,3) == rowsupport(B,3) == 3:3
    @test colsupport(B,5) == rowsupport(B,5) == 5:5
    B = BroadcastArray(*, Diagonal(randn(5)), 2)
    @test B == broadcast(*,B.args...)
    @test colsupport(B,1) == rowsupport(B,1) == 1:1
    @test colsupport(B,3) == rowsupport(B,3) == 3:3
    @test colsupport(B,5) == rowsupport(B,5) == 5:5
    B = BroadcastArray(*, Diagonal(randn(5)), randn(5))
    @test B == broadcast(*,B.args...)
    @test colsupport(B,1) == rowsupport(B,1) == 1:1
    @test colsupport(B,3) == rowsupport(B,3) == 3:3
    @test colsupport(B,5) == rowsupport(B,5) == 5:5

    B = BroadcastArray(+, Diagonal(randn(5)), 2)
    @test colsupport(B,1) == rowsupport(B,1) == 1:5
    @test colsupport(B,3) == rowsupport(B,3) == 1:5
    @test colsupport(B,5) == rowsupport(B,5) == 1:5
end

@testset "vector*matrix broadcasting #27" begin
    H = [1., 0.]
    @test Mul(H, H') .+ 1 == H*H' .+ 1
    B =  randn(2,2)
    @test Mul(H, H') .+ B == H*H' .+ B
end

@testset "BroadcastArray +" begin
    a = BroadcastArray(+, randn(400), randn(400))
    b = similar(a)
    copyto!(b, a)
    if VERSION ≥ v"1.1"
        @test @allocated(copyto!(b, a)) == 0
    end
    @test b == a
end

@testset "Lazy range" begin
    @test broadcasted(LazyArrayStyle{1}(), +, 1:5) ≡ 1:5
    @test broadcasted(LazyArrayStyle{1}(), +, 1, 1:5) ≡ 2:6
    @test broadcasted(LazyArrayStyle{1}(), +, 1:5, 1) ≡ 2:6

    @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5)) ≡ Fill(2,5)
    @test broadcasted(LazyArrayStyle{1}(), +, 1, Fill(2,5)) ≡ Fill(3,5)
    @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), 1) ≡ Fill(3,5)
    @test broadcasted(LazyArrayStyle{1}(), +, Ref(1), Fill(2,5)) ≡ Fill(3,5)
    @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), Ref(1)) ≡ Fill(3,5)
    @test broadcasted(LazyArrayStyle{1}(), +, 1, Fill(2,5)) ≡ Fill(3,5)
    @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), Fill(3,5)) ≡ Fill(5,5)

    @test broadcasted(LazyArrayStyle{1}(), *, Zeros(5), Zeros(5)) ≡ Zeros(5)
    b = BroadcastArray(exp, randn(5))
    @test b .* Zeros(5) ≡ Zeros(5)
    @test Zeros(5) .* b ≡ Zeros(5)
end


text:
module LazyGroupBy

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end LazyGroupBy

export grouped

using BangBang: merge!!
using BangBang.NoBang: SingletonVector, SingletonDict
using Base: add_sum, mul_prod
using Base.Broadcast: Broadcasted
using InitialValues: InitialValues
using OnlineStats: Mean, Variance, value, nobs
using Statistics: Statistics, mean, std, var
using Transducers
using Transducers: IdentityTransducer, ZipSource, extract_transducer

if !@isdefined only
    using Compat: only
end

const FoldFunction0 = Union{typeof(foldl),typeof(reduce),typeof(dreduce)}
if @isdefined foldxl
    const FoldFunction1 = Union{typeof(foldxl),typeof(foldxt),typeof(foldxd)}
else
    const FoldFunction1 = Union{}
    const foldxl = foldl
    const foldxt = reduce
    const foldxd = dreduce
end
const FoldFunction = Union{FoldFunction0,FoldFunction1}

# Used only for `KWBroadcastedInner` below:
function air end
struct Aired{T}
    value::T
end
Broadcast.broadcasted(::typeof(air), x) = Aired(x)
Broadcast.materialize(x::Aired) = x.value

# A hack to support broadcasting with keyword arguments:
const KWBroadcastedInner = let T = typeof(air.(identity.(; dummy = 1)).f)
    getfield(parentmodule(T), nameof(T))
end
# TODO: Use `FixKwargs` after:
# https://github.com/JuliaLang/julia/pull/36093

"""
    grouped(key, collection)

Create a lazy associative (dict-like) object grouped by a function
`key`.  Actual per-group reduction can be initiated by the dot-call
(broadcasting) of the "reducers" like `foldl` and `reduce`.

# Examples
```jldoctest
julia> using LazyGroupBy

julia> length.(grouped(isodd, 1:7))
Transducers.GroupByViewDict{Bool,Int64,…} with 2 entries:
  false => 3
  true  => 4
```
"""
grouped(key, collection) = GroupedBy(key, collection)

struct GroupedBy{K,C}
    key::K
    collection::C
end

struct GroupedByStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:GroupedBy}) = GroupedByStyle()
Base.BroadcastStyle(::GroupedByStyle, ::Broadcast.BroadcastStyle) = GroupedByStyle()

Base.broadcastable(x::GroupedBy) = x
Broadcast.instantiate(bc::Broadcasted{GroupedByStyle}) = bc

unwrap0d(x::Ref) = x[]
unwrap0d(x) = Base.IteratorSize(x) isa Base.HasShape{0} ? only(x) : x

_namedtuple(kwargs) = kwargs.data::NamedTuple

function Base.copy(bc::Broadcasted{GroupedByStyle})
    if bc.f isa KWBroadcastedInner
        return impl(bc.f.f, _namedtuple(bc.f.kwargs), map(unwrap0d, bc.args)...)
    else
        return impl(bc.f, NamedTuple(), map(unwrap0d, bc.args)...)
    end
end

check_kwargs(_, kwargs::Union{NamedTuple{()},NamedTuple{(:init,)}}) = kwargs
check_kwargs(f, kwargs::NamedTuple) =
    throw(ArgumentError("invalid keyword arguments for `$f`: $kwargs"))

_groupby_kwargs(; init = Transducers.DefaultInit, kwargs...) = ((init,), kwargs)

function _impl_fold(fold, kwargs, rf, xf, group)
    (init,), kwargs = _groupby_kwargs(; kwargs...)
    prexf, collection = extract_transducer(group.collection)
    gxf = opcompose(prexf, GroupBy(group.key, opcompose(Map(last), xf), rf, init))
    d = fold(right, gxf, collection; init = nothing, kwargs...)
    if d === nothing
        error("input collection is empty or all filtered out")
    end
    return d
end

impl(fold::FoldFunction, kwargs, rf, xf::Transducer, group::GroupedBy) =
    _impl_fold(fold, kwargs, rf, xf, group)
impl(fold::FoldFunction, kwargs, rf, xs) = impl(fold, kwargs, rf, Map(identity), xs)

impl(::typeof(mapfoldl), kwargs, f, rf, xs) = impl(foldl, kwargs, rf, Map(f), xs)
impl(::typeof(mapreduce), kwargs, f, rf, xs) = impl(reduce, kwargs, rf, Map(f), xs)

impl(::typeof(collect), kwargs, xf::Transducer, group::GroupedBy) = foldl(
    right,
    GroupBy(group.key, opcompose(Map(last), xf, Map(SingletonVector)), append!!),
    group.collection,
)
impl(::typeof(collect), kwargs, group::GroupedBy) =
    impl(collect, kwargs, IdentityTransducer(), group)

impl(::typeof(tcollect), kwargs, xf::Transducer, group::GroupedBy) = foldxt(
    right,
    GroupBy(group.key, opcompose(Map(last), xf, Map(SingletonVector)), append!!),
    group.collection,
)
impl(::typeof(tcollect), kwargs, group::GroupedBy) =
    impl(tcollect, kwargs, IdentityTransducer(), group)

impl(::typeof(dcollect), kwargs, xf::Transducer, group::GroupedBy) = foldxd(
    right,
    GroupBy(group.key, opcompose(Map(last), xf, Map(SingletonVector)), append!!),
    group.collection;
    kwargs...,
)
impl(::typeof(dcollect), kwargs, group::GroupedBy) =
    impl(dcollect, kwargs, IdentityTransducer(), group)

impl(::typeof(map), kwargs, f, group::GroupedBy) = impl(collect, kwargs, Map(f), group)

impl(::typeof(sum), kwargs::NamedTuple, group::GroupedBy) =
    impl(foldl, kwargs, add_sum, IdentityTransducer(), group)
impl(::typeof(sum), kwargs::NamedTuple, f, group::GroupedBy) =
    impl(foldl, check_kwargs(sum, kwargs), add_sum, Map(f), group)

impl(::typeof(prod), kwargs::NamedTuple, group::GroupedBy) =
    impl(foldl, kwargs, mul_prod, IdentityTransducer(), group)
impl(::typeof(prod), kwargs::NamedTuple, f, group::GroupedBy) =
    impl(foldl, check_kwargs(prod, kwargs), mul_prod, Map(f), group)

struct AsBool{F}
    f::F
end
AsBool(::Type{T}) where {T} = AsBool{Type{T}}(T)
@inline (f::AsBool)(x) = f.f(x)::Bool

impl(::typeof(count), ::NamedTuple{()}, group::GroupedBy) =
    impl(foldl, (init = 0,), +, Map(AsBool(identity)), group)
impl(::typeof(count), ::NamedTuple{()}, f, group::GroupedBy) =
    impl(foldl, (init = 0,), +, Map(AsBool(f)), group)
impl(::typeof(length), ::NamedTuple{()}, group::GroupedBy) =
    impl(foldl, (init = 0,), +, Map(_ -> 1), group)

# impl(::typeof(any), ::NamedTuple{()}, group::GroupedBy) =
#     impl(foldl, NamedTuple(), |, IdentityTransducer(), group)
impl(::typeof(any), ::NamedTuple{()}, f, group::GroupedBy) =
    impl(foldl, NamedTuple(), |, Map(f), group)

# impl(::typeof(all), ::NamedTuple{()}, group::GroupedBy) =
#     impl(foldl, NamedTuple(), &, IdentityTransducer(), group)
impl(::typeof(all), ::NamedTuple{()}, f, group::GroupedBy) =
    impl(foldl, NamedTuple(), &, Map(f), group)

impl(::typeof(minimum), kwargs::NamedTuple, group::GroupedBy) =
    impl(minimum, kwargs, identity, group)
impl(::typeof(minimum), kwargs::NamedTuple, f, group::GroupedBy) =
    impl(mapfoldl, check_kwargs(minimum, kwargs), f, min, group)
impl(::typeof(maximum), kwargs::NamedTuple, group::GroupedBy) =
    impl(maximum, kwargs, identity, group)
impl(::typeof(maximum), kwargs::NamedTuple, f, group::GroupedBy) =
    impl(mapfoldl, check_kwargs(maximum, kwargs), f, max, group)

impl(::typeof(extrema), kwargs::NamedTuple, group::GroupedBy) =
    impl(extrema, kwargs, identity, group)
impl(::typeof(extrema), kwargs::NamedTuple, f, group::GroupedBy) =
    impl(mapfoldl, check_kwargs(extrema, kwargs), DupY(f), extrema_rf, group)

struct DupY{F} <: Function
    f::F
end
DupY(f::Type{T}) where {T} = DupY{Type{T}}(f)
@inline (f::DupY)(x) = (y = f.f(x); (y, y))

extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))

impl(::typeof(keys), ::NamedTuple{()}, group::GroupedBy) = foldl(
    right,
    GroupBy(
        group.key ∘ last,
        Map(((_key, (idx, _val)),) -> SingletonVector(idx)),
        append!!,
    ),
    pairs(group.collection),
)

impl(::typeof(pairs), ::NamedTuple{()}, group::GroupedBy) = foldl(
    right,
    GroupBy(
        group.key ∘ last,
        Map(((_key, (idx, val)),) -> SingletonDict(idx => val)),
        merge!!,
    ),
    pairs(group.collection),
)

impl(::typeof(view), ::NamedTuple{()}, group::GroupedBy) =
    Dict(k => view(group.collection, idx) for (k, idx) in impl(keys, NamedTuple(), group))

left(x, _) = x
left(x) = x
InitialValues.@def_monoid left

# `identity` here is pretty useless. Maybe `_ -> true`?
# impl(::typeof(findfirst), ::NamedTuple{()}, group::GroupedBy) =
#     impl(findfirst, NamedTuple(), identity, group)
impl(::typeof(findfirst), ::NamedTuple{()}, f, group::GroupedBy) = foldl(
    right,
    GroupBy(
        group.key ∘ last,
        opcompose(
            Filter(((_key, (_idx, val)),) -> f(val)),
            Map(((_key, (idx, _val)),) -> idx),
        ),
        left,
    ),
    pairs(group.collection),
)

# impl(::typeof(findlast), ::NamedTuple{()}, group::GroupedBy) =
#     impl(findlast, NamedTuple(), identity, group)
impl(::typeof(findlast), ::NamedTuple{()}, f, group::GroupedBy) = foldl(
    right,
    GroupBy(
        group.key ∘ last,
        opcompose(Filter(((_key, (_idx, val)),) -> f(val)), Map(((_key, (idx, _val)),) -> idx)),
        right,
    ),
    pairs(group.collection),
)

impl(::typeof(findall), kwargs, f, group::GroupedBy) = foldl(
    right,
    GroupBy(
        group.key ∘ last,
        opcompose(
            Filter(((_key, (_idx, val)),) -> f(val)),
            Map(((_key, (idx, _val)),) -> SingletonVector(idx)),
        ),
        append!!,
    ),
    pairs(group.collection),
)

# Statistics
function _group_foldl_os(::NamedTuple{()}, os, xf, group)
    d = impl(foldl, NamedTuple(), os, xf, group)
    return Dict(k => value(v) for (k, v) in pairs(d))
end

impl(::typeof(mean), kwargs, group::GroupedBy) =
    impl(mean, kwargs, identity, group::GroupedBy)
impl(::typeof(mean), kwargs, f, group::GroupedBy) = impl(mean, kwargs, Map(f), group)
impl(::typeof(mean), kwargs, xf::Transducer, group::GroupedBy) =
    _group_foldl_os(kwargs, Mean(), xf, group)

impl(::typeof(var), kwargs, group::GroupedBy) =
    impl(var, kwargs, identity, group::GroupedBy)
impl(::typeof(var), kwargs, f, group::GroupedBy) = impl(var, kwargs, Map(f), group)
function impl(::typeof(var), kwargs, xf::Transducer, group::GroupedBy)
    d = impl(foldl, NamedTuple(), Variance(), xf, group)
    return Dict(k => nobs(v) < 2 ? NaN : value(v) for (k, v) in pairs(d))
end

impl(::typeof(std), kwargs, group::GroupedBy) =
    impl(std, kwargs, identity, group::GroupedBy)
impl(::typeof(std), kwargs, f, group::GroupedBy) = impl(std, kwargs, Map(f), group)
function impl(::typeof(std), ::NamedTuple{()}, xf::Transducer, group::GroupedBy)
    d = impl(foldl, NamedTuple(), Variance(), xf, group)
    return Dict(k => nobs(v) < 2 ? NaN : sqrt(value(v)) for (k, v) in pairs(d))
end

# Helper function for Documenter
_is_public(x) =
    x === (@__MODULE__) ||
    x === grouped ||
    parentmodule(x) in (Base, Statistics, Transducers)

function _import_docs()
    processed = Symbol[]
    # Use markdown files as docstring:
    for file in readdir(joinpath(@__DIR__, "docs"))
        endswith(file, ".md") || continue
        name = Symbol(file[1:end-length(".md")])
        path = joinpath(@__DIR__, "docs", file)
        # Don't fail when somehow importing docstrings doesn't work (we
        # don't loose any functionalities for that).
        try
            include_dependency(path)
            str = read(path, String)
            str = replace(str, r"^```julia"m => "```jldoctest $name")
            @eval @doc $str $name
        catch err
            @error(
                "Failed to import docstring for $name",
                exception = (err, catch_backtrace()),
            )
        end
        push!(processed, name)
    end
    return processed
end
_import_docs()

end # module


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ a11a5250-c5a3-11eb-2cb2-65a0eb1d220d
md"""
### Test Code for Module LogReg
"""

# ╔═╡ cff452a8-1123-4228-bdba-a6613e02c5bd
module script include("script.jl") end

# ╔═╡ a54d8027-2db7-4b0e-8794-f614e6b5ad50
LogReg = script.LogReg

# ╔═╡ 65be2db3-70fb-4cc3-ae2d-d53d9c94740f
import Plots

# ╔═╡ bfefa9b3-ecfb-496c-9176-9882672541d8
import Random

# ╔═╡ 394a6ebc-7544-499b-ad97-7b4452b1555e
begin
# generate linearly separated random data
X_data = Random.rand(0.0:0.5:100.0, (800, 2)) # first column X-axis, second column Y-axis
Y_data = Array{Int64}(undef, size(X_data)[1])
linear_func(x) = 2 * x - 3
for i = 1:size(X_data)[1]
    if X_data[i, 2] > linear_func(X_data[i, 1]) + randn()*20
        Y_data[i] = 0
    else
        Y_data[i] = 1
    end
end
size(X_data), size(Y_data)
end

# ╔═╡ 389887b7-2f1e-4dda-90c0-40f42f1eb5e7
# scale data
X_data_scaled = LogReg.scale(X_data);

# ╔═╡ 649f29fe-9f96-4be3-bff7-dfbbf0c1700f
let
# plot the data
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
Plots.gr()
Plots.scatter(plot_X1, plot_Y1, leg=false, c="red", background_color=Plots.RGB(0.2, 0.2, 0.2))
Plots.scatter!(plot_X2, plot_Y2, leg=false, c="blue")
Plots.plot!(plot_X3, plot_Y3, leg=false, c="orange", linewidth=5)
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ c7df3c2c-54d8-4177-ba02-4d0098087cf6
let
# run logistic regression
weights = LogReg.train(X_data_scaled, Y_data, learning_rate=0.1, max_iter=10, early_stop=false)
predictions = LogReg.predict(X_data_scaled, weights)
string("Accuracy: ", LogReg.accuracy(predictions, Y_data))
end

# ╔═╡ c4cdf49e-3929-4393-a9fc-5db66cd2a812
let
# increase iteractions
weights = LogReg.train(X_data_scaled, Y_data, learning_rate=0.1, max_iter=100, early_stop=false)
predictions = LogReg.predict(X_data_scaled, weights)
string("Accuracy: ", LogReg.accuracy(predictions, Y_data))
end

# ╔═╡ 5e698890-39fb-47dc-a975-53d6c8e716a6
let
# increase iteractions
weights = LogReg.train(X_data_scaled, Y_data, learning_rate=0.1, max_iter=1000, early_stop=false)
predictions = LogReg.predict(X_data_scaled, weights)
string("Accuracy: ", LogReg.accuracy(predictions, Y_data))
end

# ╔═╡ 6bb2b4a7-2f56-4d42-a5ac-29333a69abba
let
# increase iteractions
weights = LogReg.train(X_data_scaled, Y_data, learning_rate=0.1, max_iter=10000, early_stop=false)
predictions = LogReg.predict(X_data_scaled, weights)
string("Accuracy: ", LogReg.accuracy(predictions, Y_data))
end

# ╔═╡ 50f6b449-6e24-48c1-8d33-dc13c1fb5b9c
let
# plot comparison graph
plot_X1 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_Y1 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 0]
plot_X2 = [X_data[i, 1] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_Y2 = [X_data[i, 2] for i in 1:(size(X_data)[1]) if Y_data[i] == 1]
plot_X3 = [i for i in 0:100]
plot_Y3 = [linear_func(m) for m in plot_X3]
Plots.gr()
Plots.scatter(plot_X1, plot_Y1, c="red", background_color=Plots.RGB(0.2, 0.2, 0.2), label="")
Plots.scatter!(plot_X2, plot_Y2, c="blue", label="")
Plots.plot!(plot_X3, plot_Y3, c="orange", linewidth=5, label="Actual")
accuracy = Array{Float64}(undef, 4)
for (i, iternum) in enumerate([10, 100, 1000, 10000])
    weights = LogReg.train(X_data_scaled, Y_data, learning_rate=0.1, max_iter=iternum)
    accuracy[i] = LogReg.accuracy(LogReg.predict(X_data_scaled, weights), Y_data)
    @assert size(weights)[1] == 3
    plot_Y3_ = -(weights[3] .+ (weights[1] .* plot_X3)) ./ weights[2]
    Plots.plot!(plot_X3, plot_Y3_, palette=:lightrainbow, linewidth=3, label="iter=$iternum", alpha=0.5)
end
Plots.xlims!((0.0, 100.0))
Plots.ylims!((0.0, 100.0))
end

# ╔═╡ Cell order:
# ╟─a11a5250-c5a3-11eb-2cb2-65a0eb1d220d
# ╠═cff452a8-1123-4228-bdba-a6613e02c5bd
# ╠═a54d8027-2db7-4b0e-8794-f614e6b5ad50
# ╠═65be2db3-70fb-4cc3-ae2d-d53d9c94740f
# ╠═bfefa9b3-ecfb-496c-9176-9882672541d8
# ╠═394a6ebc-7544-499b-ad97-7b4452b1555e
# ╠═389887b7-2f1e-4dda-90c0-40f42f1eb5e7
# ╠═649f29fe-9f96-4be3-bff7-dfbbf0c1700f
# ╠═c7df3c2c-54d8-4177-ba02-4d0098087cf6
# ╠═c4cdf49e-3929-4393-a9fc-5db66cd2a812
# ╠═5e698890-39fb-47dc-a975-53d6c8e716a6
# ╠═6bb2b4a7-2f56-4d42-a5ac-29333a69abba
# ╠═50f6b449-6e24-48c1-8d33-dc13c1fb5b9c

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using MOTProject
using Test

@testset "MOTProject.jl" begin

end


text:
const J0_2N = (
    big"3.133239376997663645548490085151484674892E16", big"-5.479944965767990821079467311839107722107E14",
    big"6.290828903904724265980249871997551894090E12", big"-3.633750176832769659849028554429106299915E10",
    big"1.207743757532429576399485415069244807022E8", big"-2.107485999925074577174305650549367415465E5",
    big"1.562826808020631846245296572935547005859E2"
)
const J0_2D = (
    big"2.005273201278504733151033654496928968261E18", big"2.063038558793221244373123294054149790864E16",
    big"1.053350447931127971406896594022010524994E14",big"3.496556557558702583143527876385508882310E11",
    big"8.249114511878616075860654484367133976306E8",big"1.402965782449571800199759247964242790589E6",
    big"1.619910762853439600957801751815074787351E3", big"1.00000000000000000000000000000000000000000"
)
const P16_IN = (
    big"-1.901689868258117463979611259731176301065E-16", big"-1.798743043824071514483008340803573980931E-13",
    big"-6.481746687115262291873324132944647438959E-11", big"-1.150651553745409037257197798528294248012E-8",
    big"-1.088408467297401082271185599507222695995E-6", big"-5.551996725183495852661022587879817546508E-5",
    big"-1.477286941214245433866838787454880214736E-3", big"-1.882877976157714592017345347609200402472E-2",
    big"-9.620983176855405325086530374317855880515E-2", big"-1.271468546258855781530458854476627766233E-1"
)

const P16_ID = (
    big"2.704625590411544837659891569420764475007E-15", big"2.562526347676857624104306349421985403573E-12",
    big"9.259137589952741054108665570122085036246E-10", big"1.651044705794378365237454962653430805272E-7",
    big"1.573561544138733044977714063100859136660E-5", big"8.134482112334882274688298469629884804056E-4",
    big"2.219259239404080863919375103673593571689E-2", big"2.976990606226596289580242451096393862792E-1",
    big"1.713895630454693931742734911930937246254E0", big"3.231552290717904041465898249160757368855E0",
    big"1.000000000000000000000000000000000000000E0"
)

const Q16_IN = (
    big"2.343640834407975740545326632205999437469E-18", big"2.667978112927811452221176781536278257448E-15",
    big"1.178415018484555397390098879501969116536E-12", big"2.622049767502719728905924701288614016597E-10",
    big"3.196908059607618864801313380896308968673E-8", big"2.179466154171673958770030655199434798494E-6",
    big"8.139959091628545225221976413795645177291E-5", big"1.563900725721039825236927137885747138654E-3",
    big"1.355172364265825167113562519307194840307E-2", big"3.928058355906967977269780046844768588532E-2",
    big"1.107891967702173292405380993183694932208E-2",
)

const Q16_ID = (
    big"3.199850952578356211091219295199301766718E-17", big"3.652601488020654842194486058637953363918E-14",
    big"1.620179741394865258354608590461839031281E-11", big"3.629359209474609630056463248923684371426E-9",
    big"4.473680923894354600193264347733477363305E-7", big"3.106368086644715743265603656011050476736E-5",
    big"1.198239259946770604954664925153424252622E-3", big"2.446041004004283102372887804475767568272E-2",
    big"2.403235525011860603014707768815113698768E-1", big"9.491006790682158612266270665136910927149E-1",
    big"1.000000000000000000000000000000000000000E0"
)

const P8_16N = (
    big"-2.335166846111159458466553806683579003632E-15", big"-1.382763674252402720401020004169367089975E-12",
    big"-3.192160804534716696058987967592784857907E-10", big"-3.744199606283752333686144670572632116899E-8",
    big"-2.439161236879511162078619292571922772224E-6", big"-9.068436986859420951664151060267045346549E-5",
    big"-1.905407090637058116299757292660002697359E-3", big"-2.164456143936718388053842376884252978872E-2",
    big"-1.212178415116411222341491717748696499966E-1", big"-2.782433626588541494473277445959593334494E-1",
    big"-1.670703190068873186016102289227646035035E-1"
)

const P8_16D = (
    big"3.321126181135871232648331450082662856743E-14", big"1.971894594837650840586859228510007703641E-11",
    big"4.571144364787008285981633719513897281690E-9", big"5.396419143536287457142904742849052402103E-7",
    big"3.551548222385845912370226756036899901549E-5", big"1.342353874566932014705609788054598013516E-3",
    big"2.899133293006771317589357444614157734385E-2", big"3.455374978185770197704507681491574261545E-1",
    big"2.116616964297512311314454834712634820514E0", big"5.850768316827915470087758636881584174432E0",
    big"5.655273858938766830855753983631132928968E0",
    big"1.000000000000000000000000000000000000000E0"
)

const P5_8N = (
    big"-1.270478335089770355749591358934012019596E-12", big"-4.007588712145412921057254992155810347245E-10",
    big"-4.815187822989597568124520080486652009281E-8", big"-2.867070063972764880024598300408284868021E-6",
    big"-9.218742195161302204046454768106063638006E-5", big"-1.635746821447052827526320629828043529997E-3",
    big"-1.570376886640308408247709616497261011707E-2", big"-7.656484795303305596941813361786219477807E-2",
    big"-1.659371030767513274944805479908858628053E-1", big"-1.185340550030955660015841796219919804915E-1",
    big"-8.920026499909994671248893388013790366712E-3",
)

const P5_8D = (
    big"1.806902521016705225778045904631543990314E-11", big"5.728502760243502431663549179135868966031E-9",
    big"6.938168504826004255287618819550667978450E-7", big"4.183769964807453250763325026573037785902E-5",
    big"1.372660678476925468014882230851637878587E-3", big"2.516452105242920335873286419212708961771E-2",
    big"2.550502712902647803796267951846557316182E-1", big"1.365861559418983216913629123778747617072E0",
    big"3.523825618308783966723472468855042541407E0", big"3.656365803506136165615111349150536282434E0",
    big"1.000000000000000000000000000000000000000E0"
)
const Q5_8N = (
    big"1.750399094021293722243426623211733898747E-13", big"6.483426211748008735242909236490115050294E-11",
    big"9.279430665656575457141747875716899958373E-9", big"6.696634968526907231258534757736576340266E-7",
    big"2.666560823798895649685231292142838188061E-5", big"6.025087697259436271271562769707550594540E-4",
    big"7.652807734168613251901945778921336353485E-3", big"5.226269002589406461622551452343519078905E-2",
    big"1.748390159751117658969324896330142895079E-1", big"2.378188719097006494782174902213083589660E-1",
    big"8.383984859679804095463699702165659216831E-2",
)

const Q5_8D = (
    big"2.389878229704327939008104855942987615715E-12", big"8.926142817142546018703814194987786425099E-10",
    big"1.294065862406745901206588525833274399038E-7", big"9.524139899457666250828752185212769682191E-6",
    big"3.908332488377770886091936221573123353489E-4", big"9.250427033957236609624199884089916836748E-3",
    big"1.263420066165922645975830877751588421451E-1", big"9.692527053860420229711317379861733180654E-1",
    big"3.937813834630430172221329298841520707954E0", big"7.603126427436356534498908111445191312181E0",
    big"5.670677653334105479259958485084550934305E0", big"1.000000000000000000000000000000000000000E0"
)

const P4_5N = (
    big"-9.791405771694098960254468859195175708252E-10", big"-1.917193059944531970421626610188102836352E-7",
    big"-1.393597539508855262243816152893982002084E-5", big"-4.881863490846771259880606911667479860077E-4",
    big"-8.946571245022470127331892085881699269853E-3", big"-8.707474232568097513415336886103899434251E-2",
    big"-4.362042697474650737898551272505525973766E-1", big"-1.032712171267523975431451359962375617386E0",
    big"-9.630502683169895107062182070514713702346E-1", big"-2.251804386252969656586810309252357233320E-1",
)

const P4_5D = (
    big"1.392555487577717669739688337895791213139E-8", big"2.748886559120659027172816051276451376854E-6",
    big"2.024717710644378047477189849678576659290E-4", big"7.244868609350416002930624752604670292469E-3",
    big"1.373631762292244371102989739300382152416E-1", big"1.412298581400224267910294815260613240668E0",
    big"7.742495637843445079276397723849017617210E0", big"2.138429269198406512028307045259503811861E1",
    big"2.651547684548423476506826951831712762610E1", big"1.167499382465291931571685222882909166935E1",
    big"1.000000000000000000000000000000000000000E0"
)

const P3r2_4N = (
    big"-2.589155123706348361249809342508270121788E-8", big"-3.746254369796115441118148490849195516593E-6",
    big"-1.985595497390808544622893738135529701062E-4", big"-5.008253705202932091290132760394976551426E-3",
    big"-6.529469780539591572179155511840853077232E-2", big"-4.468736064761814602927408833818990271514E-1",
    big"-1.556391252586395038089729428444444823380E0", big"-2.533135309840530224072920725976994981638E0",
    big"-1.605509621731068453869408718565392869560E0", big"-2.518966692256192789269859830255724429375E-1"
)

const Q8_16N = (
    big"1.001954266485599464105669390693597125904E-17", big"7.545499865295034556206475956620160007849E-15",
    big"2.267838684785673931024792538193202559922E-12", big"3.561909705814420373609574999542459912419E-10",
    big"3.216201422768092505214730633842924944671E-8", big"1.731194793857907454569364622452058554314E-6",
    big"5.576944613034537050396518509871004586039E-5", big"1.051787760316848982655967052985391418146E-3",
    big"1.102852974036687441600678598019883746959E-2", big"5.834647019292460494254225988766702933571E-2",
    big"1.290281921604364618912425380717127576529E-1", big"7.598886310387075708640370806458926458301E-2",
)
const Q8_16D = (
    big"1.368001558508338469503329967729951830843E-16", big"1.034454121857542147020549303317348297289E-13",
    big"3.128109209247090744354764050629381674436E-11", big"4.957795214328501986562102573522064468671E-9",
    big"4.537872468606711261992676606899273588899E-7", big"2.493639207101727713192687060517509774182E-5",
    big"8.294957278145328349785532236663051405805E-4", big"1.646471258966713577374948205279380115839E-2",
    big"1.878910092770966718491814497982191447073E-1", big"1.152641605706170353727903052525652504075E0",
    big"3.383550240669773485412333679367792932235E0", big"3.823875252882035706910024716609908473970E0",
    big"1.000000000000000000000000000000000000000E0"
)
const Q4_5N = (
    big"2.233870042925895644234072357400122854086E-11", big"5.146223225761993222808463878999151699792E-9",
    big"4.459114531468296461688753521109797474523E-7", big"1.891397692931537975547242165291668056276E-5",
    big"4.279519145911541776938964806470674565504E-4", big"5.275239415656560634702073291768904783989E-3",
    big"3.468698403240744801278238473898432608887E-2", big"1.138773146337708415188856882915457888274E-1",
    big"1.622717518946443013587108598334636458955E-1", big"7.249040006390586123760992346453034628227E-2",
    big"1.941595365256460232175236758506411486667E-3",
)

const Q4_5D = (
    big"3.049977232266999249626430127217988047453E-10", big"7.120883230531035857746096928889676144099E-8",
    big"6.301786064753734446784637919554359588859E-6", big"2.762010530095069598480766869426308077192E-4",
    big"6.572163250572867859316828886203406361251E-3", big"8.752566114841221958200215255461843397776E-2",
    big"6.487654992874805093499285311075289932664E-1", big"2.576550017826654579451615283022812801435E0",
    big"5.056392229924022835364779562707348096036E0", big"4.179770081068251464907531367859072157773E0",
    big"1.000000000000000000000000000000000000000E0"
)

const P3r2_4D = (
    big"3.682353957237979993646169732962573930237E-7",
    big"5.386741661883067824698973455566332102029E-5",
    big"2.906881154171822780345134853794241037053E-3",
    big"7.545832595801289519475806339863492074126E-2",
    big"1.029405357245594877344360389469584526654E0",
    big"7.565706120589873131187989560509757626725E0",
    big"2.951172890699569545357692207898667665796E1",
    big"5.785723537170311456298467310529815457536E1",
    big"5.095621464598267889126015412522773474467E1",
    big"1.602958484169953109437547474953308401442E1",
    big"1.000000000000000000000000000000000000000E0"
)

const P2r7_3r2N = (
    big"-1.917322340814391131073820537027234322550E-7",
    big"-1.966595744473227183846019639723259011906E-5",
    big"-7.177081163619679403212623526632690465290E-4",
    big"-1.206467373860974695661544653741899755695E-2",
    big"-1.008656452188539812154551482286328107316E-1",
    big"-4.216016116408810856620947307438823892707E-1",
    big"-8.378631013025721741744285026537009814161E-1",
    big"-6.973895635309960850033762745957946272579E-1",
    big"-1.797864718878320770670740413285763554812E-1",
    big"-4.098025357743657347681137871388402849581E-3"
)

const P2r7_3r2D = (
    big"2.726858489303036441686496086962545034018E-6",
    big"2.840430827557109238386808968234848081424E-4",
    big"1.063826772041781947891481054529454088832E-2",
    big"1.864775537138364773178044431045514405468E-1",
    big"1.665660052857205170440952607701728254211E0",
    big"7.723745889544331153080842168958348568395E0",
    big"1.810726427571829798856428548102077799835E1",
    big"1.986460672157794440666187503833545388527E1",
    big"8.645503204552282306364296517220055815488E0",
    big"1.000000000000000000000000000000000000000E0"
)

const P2r3_2r7N = (
    big"-1.594642785584856746358609622003310312622E-6",
    big"-1.323238196302221554194031733595194539794E-4",
    big"-3.856087818696874802689922536987100372345E-3",
    big"-5.113241710697777193011470733601522047399E-2",
    big"-3.334229537209911914449990372942022350558E-1",
    big"-1.075703518198127096179198549659283422832E0",
    big"-1.634174803414062725476343124267110981807E0",
    big"-1.030133247434119595616826842367268304880E0",
    big"-1.989811539080358501229347481000707289391E-1",
    big"-3.246859189246653459359775001466924610236E-3"
)

const P2r3_2r7D = (
    big"2.267936634217251403663034189684284173018E-5",
    big"1.918112982168673386858072491437971732237E-3",
    big"5.771704085468423159125856786653868219522E-2",
    big"8.056124451167969333717642810661498890507E-1",
    big"5.687897967531010276788680634413789328776E0",
    big"2.072596760717695491085444438270778394421E1",
    big"3.801722099819929988585197088613160496684E1",
    big"3.254620235902912339534998592085115836829E1",
    big"1.104847772130720331801884344645060675036E1",
    big"1.000000000000000000000000000000000000000E0"
)

const P2_2r3N = (
    big"-1.001042324337684297465071506097365389123E-4",
    big"-6.289034524673365824853547252689991418981E-3",
    big"-1.346527918018624234373664526930736205806E-1",
    big"-1.268808313614288355444506172560463315102E0",
    big"-5.654126123607146048354132115649177406163E0",
    big"-1.186649511267312652171775803270911971693E1",
    big"-1.094032424931998612551588246779200724257E1",
    big"-3.728792136814520055025256353193674625267E0",
    big"-3.000348318524471807839934764596331810608E-1",
)

const P2_2r3D = (
    big"1.423705538269770974803901422532055612980E-3",
    big"9.171476630091439978533535167485230575894E-2",
    big"2.049776318166637248868444600215942828537E0",
    big"2.068970329743769804547326701946144899583E1",
    big"1.025103500560831035592731539565060347709E2",
    big"2.528088049697570728252145557167066708284E2",
    big"2.992160327587558573740271294804830114205E2",
    big"1.540193761146551025832707739468679973036E2",
    big"2.779516701986912132637672140709452502650E1",
    big"1.000000000000000000000000000000000000000E0"
)

const Q2r7_3r2N = (
    big"7.584861620402450302063691901886141875454E-8",
    big"9.300939338814216296064659459966041794591E-6",
    big"4.112108906197521696032158235392604947895E-4",
    big"8.515168851578898791897038357239630654431E-3",
    big"8.971286321017307400142720556749573229058E-2",
    big"4.885856732902956303343015636331874194498E-1",
    big"1.334506268733103291656253500506406045846E0",
    big"1.681207956863028164179042145803851824654E0",
    big"8.165042692571721959157677701625853772271E-1",
    big"9.805848115375053300608712721986235900715E-2"
)
const Q2r7_3r2D = (
    big"1.035586492113036586458163971239438078160E-6",
    big"1.301999337731768381683593636500979713689E-4",
    big"5.993695702564527062553071126719088859654E-3",
    big"1.321184892887881883489141186815457808785E-1",
    big"1.528766555485015021144963194165165083312E0",
    big"9.561463309176490874525827051566494939295E0",
    big"3.203719484883967351729513662089163356911E1",
    big"5.497294687660930446641539152123568668447E1",
    big"4.391158169390578768508675452986948391118E1",
    big"1.347836630730048077907818943625789418378E1",
    big"1.000000000000000000000000000000000000000E0"
)

const Q3r2_4N = (
    big"6.126167301024815034423262653066023684411E-10",
    big"1.043969327113173261820028225053598975128E-7",
    big"6.592927270288697027757438170153763220190E-6",
    big"2.009103660938497963095652951912071336730E-4",
    big"3.220543385492643525985862356352195896964E-3",
    big"2.774405975730545157543417650436941650990E-2",
    big"1.258114008023826384487378016636555041129E-1",
    big"2.811724258266902502344701449984698323860E-1",
    big"2.691837665193548059322831687432415014067E-1",
    big"7.949087384900985370683770525312735605034E-2",
    big"1.229509543620976530030153018986910810747E-3",
)
const Q3r2_4D = (
    big"8.364260446128475461539941389210166156568E-9",
    big"1.451301850638956578622154585560759862764E-6",
    big"9.431830010924603664244578867057141839463E-5",
    big"3.004105101667433434196388593004526182741E-3",
    big"5.148157397848271739710011717102773780221E-2",
    big"4.901089301726939576055285374953887874895E-1",
    big"2.581760991981709901216967665934142240346E0",
    big"7.257105880775059281391729708630912791847E0",
    big"1.006014717326362868007913423810737369312E1",
    big"5.879416600465399514404064187445293212470E0",
    big"1.000000000000000000000000000000000000000E0"
)

const Q2r3_2r7N = (
    big"4.455027774980750211349941766420190722088E-7",
    big"4.031998274578520170631601850866780366466E-5",
    big"1.273987274325947007856695677491340636339E-3",
    big"1.818754543377448509897226554179659122873E-2",
    big"1.266748858326568264126353051352269875352E-1",
    big"4.327578594728723821137731555139472880414E-1",
    big"6.892532471436503074928194969154192615359E-1",
    big"4.490775818438716873422163588640262036506E-1",
    big"8.649615949297322440032000346117031581572E-2",
    big"7.261345286655345047417257611469066147561E-4",
)
const Q2r3_2r7D = (
    big"6.082600739680555266312417978064954793142E-6",
    big"5.693622538165494742945717226571441747567E-4",
    big"1.901625907009092204458328768129666975975E-2",
    big"2.958689532697857335456896889409923371570E-1",
    big"2.343124711045660081603809437993368799568E0",
    big"9.665894032187458293568704885528192804376E0",
    big"2.035273104990617136065743426322454881353E1",
    big"2.044102010478792896815088858740075165531E1",
    big"8.445937177863155827844146643468706599304E0",
    big"1.000000000000000000000000000000000000000E0"
)

const Q2_2r3N = (
    big"2.817566786579768804844367382809101929314E-6",
    big"2.122772176396691634147024348373539744935E-4",
    big"5.501378031780457828919593905395747517585E-3",
    big"6.355374424341762686099147452020466524659E-2",
    big"3.539652320122661637429658698954748337223E-1",
    big"9.571721066119617436343740541777014319695E-1",
    big"1.196258777828426399432550698612171955305E0",
  big"6.069388659458926158392384709893753793967E-1",
  big"9.026746127269713176512359976978248763621E-2",
  big"5.317668723070450235320878117210807236375E-4",
)

const Q2_2r3D = (
    big"3.846924354014260866793741072933159380158E-5",
    big"3.017562820057704325510067178327449946763E-3",
    big"8.356305620686867949798885808540444210935E-2",
    big"1.068314930499906838814019619594424586273E0",
    big"6.900279623894821067017966573640732685233E0",
    big"2.307667390886377924509090271780839563141E1",
    big"3.921043465412723970791036825401273528513E1",
    big"3.167569478939719383241775717095729233436E1",
    big"1.051023841699200920276198346301543665909E1",
    big"1.000000000000000000000000000000000000000E0"
)

const Y0_2N = (
    big"-1.062023609591350692692296993537002558155E19",
    big"2.542000883190248639104127452714966858866E19",
    big"-1.984190771278515324281415820316054696545E18",
    big"4.982586044371592942465373274440222033891E16",
    big"-5.529326354780295177243773419090123407550E14",
    big"3.013431465522152289279088265336861140391E12",
    big"-7.959436160727126750732203098982718347785E9",
    big"8.230845651379566339707130644134372793322E6",
)
const Y0_2D = (
    big"1.438972634353286978700329883122253752192E20",
    big"1.856409101981569254247700169486907405500E18",
    big"1.219693352678218589553725579802986255614E16",
    big"5.389428943282838648918475915779958097958E13",
    big"1.774125762108874864433872173544743051653E11",
    big"4.522104832545149534808218252434693007036E8",
    big"8.872187401232943927082914504125234454930E5",
    big"1.251945613186787532055610876304669413955E3",
    big"1.000000000000000000000000000000000000000E0"
)


text:
export Move

export MOVE_NULL, from, ispromotion, movefromstring, promotion, to, tostring

"""
    Move

Type representing a chess move.

A `Move` value is usually obtained by asking a chess board for moves, or by
parsing a move string in UCI or SAN format (in the latter case, we also need
a board).
"""
struct Move
    val::Int
end

Base.show(io::IO, m::Move) = print(io, "Move($(tostring(m)))")

"""
    Move(from::Square, to::Square)

Low-level constructor for creating a move with the given `from` and `to`
squares.
"""
function Move(from::Square, to::Square)
    Move((to.val - 1) | ((from.val - 1) << 6))
end

"""
    Move(from::Square, to::Square, promotion::PieceType)

Low-level constructor for creating a move with the given `from` and `to`
squares and promotion piece type.
"""
function Move(from::Square, to::Square, promotion::PieceType)
    Move((to.val - 1) | ((from.val - 1) << 6) | (promotion.val << 12))
end

const MOVE_NULL = Move(SQ_A1, SQ_A1)

"""
    from(m::Move)

The source square of a move.

# Examples

```julia-repl
julia> Move(SQ_D2, SQ_D4)
Move(d2d4)

julia> from(Move(SQ_G1, SQ_F3))
SQ_G1

julia> from(Move(SQ_C7, SQ_C8, QUEEN))
SQ_C7
```
"""
function from(m::Move)::Square
    Square(((m.val >> 6) & 63) + 1)
end

"""
    to(m::Move)

The destination square of a move.

# Examples

```julia-repl
julia> to(Move(SQ_G1, SQ_F3))
SQ_F3

julia> to(Move(SQ_C7, SQ_C8, QUEEN))
SQ_C8
true
```
"""
function to(m::Move)::Square
    Square((m.val & 63) + 1)
end

"""
    ispromotion(m::Move)

Determine whether a move is a promotion move.

# Examples

```julia-repl
julia> ispromotion(Move(SQ_G1, SQ_F3))
false

julia> ispromotion(Move(SQ_C7, SQ_C8, QUEEN))
true
```
"""
function ispromotion(m::Move)::Bool
    m.val & (7 << 12) ≠ 0
end

"""
    promotion(m::Move)

Find the promotion piece type of a move.

Use this function only after first using `ispromotion` to determine whether
the move is a promotion move at all.

# Examples

```julia-repl
julia> promotion(Move(SQ_C7, SQ_C8, QUEEN))
QUEEN

julia> promotion(Move(SQ_B2, SQ_B1, KNIGHT))
KNIGHT
```
"""
function promotion(m::Move)::PieceType
    PieceType((m.val >> 12) & 7)
end

"""
    tostring(m::Move)

Convert a move to a string in UCI notation.

# Examples

```julia-repl
julia> tostring(Move(SQ_G1, SQ_F3))
"g1f3"

julia> tostring(Move(SQ_E2, SQ_E1, KNIGHT))
"e2e1n"
```
"""
function tostring(m::Move)::String
    if m == MOVE_NULL
        "0000"
    else
        tostring(from(m)) * tostring(to(m)) * (ispromotion(m) ? tochar(promotion(m)) : "")
    end
end

"""
    movefromstring(s::String)

Convert a UCI move string to a move.

Returns `nothing` if the input string is not a valid UCI move.

# Examples

```julia-repl
julia> movefromstring("d2d4") == Move(SQ_D2, SQ_D4)
true

julia> movefromstring("h7h8q") == Move(SQ_H7, SQ_H8, QUEEN)
true

julia> movefromstring("f7f9") == nothing
true

julia> movefromstring("") == nothing
true
```
"""
function movefromstring(s::String)::Union{Move,Nothing}
    if s == "0000"
        return MOVE_NULL
    elseif length(s) < 4
        nothing
    else
        f = squarefromstring(s[1:2])
        t = squarefromstring(s[3:4])
        if f ≠ nothing && t ≠ nothing
            p = piecetypefromchar(get(s, 5, '?'))
            if isnothing(p)
                Move(f, t)
            else
                Move(f, t, p)
            end
        else
            nothing
        end
    end
end


text:
"""
    ObjectSetting

The current settings of an [`Object`](@ref) which are saved in `object.current_setting`.

# Fields
- `line_width::Float64`: the current line width
- `mul_line_width::Float64`: the current multiplier for line width.
    The actual line width is then: `mul_line_width * line_width`
- `opacity::Float64`: the current opacity
- `mul_opacity::Float64`: the current multiplier for opacity.
    The actual opacity is then: `mul_opacity * opacity`
- `fontsize::Float64` the current font size
- `show_object::Bool` is set to false if scale would be 0.0 which is forbidden by Cairo
- `current_scale::Tuple{Float64, Float64}`: the current scale
- `desired_scale::Tuple{Float64, Float64}`: the new desired scale
- `mul_scale::Float64`: the multiplier for the new desired scale.
    The actual new scale is then: `mul_scale * desired_scale`
"""
mutable struct ObjectSetting
    line_width::Float64
    mul_line_width::Float64 # the multiplier of line width is between 0 and 1
    opacity::Float64
    mul_opacity::Float64 # the multiplier of opacity is between 0 and 1
    fontsize::Float64
    # scale has three fields instead of just the normal two
    # current scale
    # desired scale and scale multiplier => `desired_scale*mul_scale` is the new desired scale
    # the scale change needs to be computed using `current_scale` and the desired scale
    # current_scale should never be 0 as this breaks scaleto has various other bad effects
    # see: https://github.com/JuliaGraphics/Luxor.jl/issues/114
    # in this case show will be set to false and the object will not be called
    show_object::Bool
    current_scale::Scale
    desired_scale::Scale
    mul_scale::Float64 # the multiplier of scale is between 0 and 1

    ObjectSetting() =
        new(2.0, 1.0, 1.0, 1.0, 10.0, true, Scale(1.0, 1.0), Scale(1.0, 1.0), 1.0)
end

"""
    update_ObjectSetting!(as::ObjectSetting, by::ObjectSetting)

Set the fields of `as` to the same as `by`. Basically copying them over.
"""
function update_ObjectSetting!(as::ObjectSetting, by::ObjectSetting)
    as.line_width = by.line_width
    as.mul_line_width = by.mul_line_width
    as.opacity = by.opacity
    as.mul_opacity = by.mul_opacity
    as.fontsize = by.fontsize
    as.show_object = by.show_object
    as.current_scale = by.current_scale
    as.desired_scale = by.desired_scale
    as.mul_scale = by.mul_scale
end

function update_background_settings!(setting::ObjectSetting, object::AbstractObject)
    in_global_layer = get(object.opts, :in_global_layer, false)
    if in_global_layer
        update_ObjectSetting!(setting, object.current_setting)
    end
end

function update_object_settings!(object::AbstractObject, setting::ObjectSetting)
    update_ObjectSetting!(object.current_setting, setting)
end


text:
# cairo interface (cairo_xcb_x)

import Cairo

# cairo_surface_t *
# cairo_xcb_surface_create (xcb_connection_t *connection,
#                           xcb_drawable_t drawable,
#                           xcb_visualtype_t *visual,
#                           int width,
#                           int height);

function cairo_xcb_surface_create(connection::xcb_connection_t, drawable::xcb_drawable_t,
	visual::xcb_visualtype_t,width::Int32,height::Int32)

    ptr = ccall((:cairo_xcb_surface_create,Cairo._jl_libcairo), Ptr{Void},
                (Ptr{xcb_connection_t}, xcb_drawable_t, Ptr{xcb_visualtype_t}, Int32, Int32),
                &connection, drawable, &visual, width, height)
    Cairo.CairoSurface(ptr, width, height)
end


text:
required_fields = Dict(
# byr (Birth Year) - four digits; at least 1920 and at most 2002.
"byr" => function (s::SubString)
    x = parse(Int64, s)
    return 1920 <= x && x <= 2002
end,
# # iyr (Issue Year) - four digits; at least 2010 and at most 2020.
"iyr" => function (s::SubString)
    x = parse(Int64, s)
    return 2010 <= x && x <= 2020
end,
# # eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
"eyr" => function (s::SubString)
    x = parse(Int64, s)
    return 2020 <= x && x <= 2030
end,
# # hgt (Height) - a number followed by either cm or in:
#     # If cm, the number must be at least 150 and at most 193.
#     # If in, the number must be at least 59 and at most 76.
"hgt" => function (s::SubString)
    x = parse(Int64, s[1:length(s) - 2])
    unit = s[length(s) - 1:length(s)]
    if unit == "cm"
        # cm case
        return 150 <= x && x <= 193
    else
        # in case
        return 59 <= x && x <= 76
    end
end,
# hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
"hcl" => function (s::SubString)
    # TODO: check for alphanumeric
    return length(s) == 7 && s[1] == '#'
end,
# # ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
"ecl" => function (s::SubString)
    return s in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
end,
# # pid (Passport ID) - a nine-digit number, including leading zeroes.
"pid" => function (s::SubString)
    # has 9 digits?
    if ! (length(s) == 9)
        return false
    end
    # is numeric?
    try
        parse(Int64, s)
        return true
    catch
        return false
    end
end,
# # "cid", # confusing: this is described as a "required" field at first but then "optional" later
)

function is_valid(passport::Dict)
    result = true
    for k in keys(required_fields)
        if ! haskey(passport, k)
            result = false
            break
        end

        rule = required_fields[k]
        if ! rule(passport[k])
            println("Failed rule ", k, " for passport ", passport)
            result = false
            break
        end
    end

    # println("IsValid? => ", result)
    return result
end

function compute(rows)
    valid_count = 0
    passport = Dict()
    for r in rows
        # end current passport
        if r == ""
            if is_valid(passport)
                valid_count += 1
            end
            passport = Dict()
        else
            # continue reading details of current passport
            items = split(r, " ")
            for i in items
                key_val = split(i, ":")
                key = key_val[1]
                passport[key_val[1]] = key_val[2]
            end
        end
    end

    # make sure to include last passport in file
    if is_valid(passport)
        valid_count += 1
    end

    return valid_count
end

function run(fname)
    println("Running", fname, " ...")
    f = open(fname, "r")
    rows = readlines(f)
    close(f)

    out = compute(rows)
    println("Final result for ", fname, " was: ", out)
    return out
end

@assert run("4ex2.txt") == 0
@assert run("4ex3.txt") == 4

run("4.txt")


text:
module MultivariateFunctions

using Dates
using SchumakerSpline
using GLM
using DataFrames: DataFrame
using DataStructures: OrderedDict
using Combinatorics: permutations
using LinearAlgebra: Symmetric, inv, det
using Optim

# This includes functions for use of MultivariateFunctions with dates.
include("date_conversions.jl")
export years_between, years_from_global_base, period_length
# The abstract type MultivariateFunction and all structs are implemented here.
# In addition we have operator reversals and some supporting functions.
include("0_structs_and_generic_reversals.jl")
export MultivariateFunction, PE_Function, Sum_Of_Functions, Piecewise_Function
export Sum_Of_Piecewise_Functions, PE_Unit
export change_base, trim_piecewise_function, sort, convert_to_linearly_rescale_inputs
export evaluate, rebadge, underlying_dimensions, convert, ≂
include("1_algebra.jl")
export +, -, *, /, ^
include("2_calculus.jl")
export derivative, all_derivatives, Hessian, jacobian, uniroot, find_local_optima
export integral
include("3_1D_splines_and_interpolation.jl")
export create_quadratic_spline, create_constant_interpolation_to_right
export create_constant_interpolation_to_left, create_linear_interpolation
include("chebyshevs.jl")
export get_chevyshevs_up_to
include("4_chebyshev_approximation.jl")
export create_chebyshev_approximation
include("5_ols_regression.jl")
export create_ols_approximation, create_saturated_ols_approximation
include("6_HighDimensionalApproximation.jl")
export create_recursive_partitioning, create_mars_spline, trim_mars_spline
end


text:
module InfrastructureSystemsPlots

export get_stacked_plot_data
export get_bar_plot_data
export get_stacked_generation_data
export get_bar_gen_data
export bar_plot
export stack_plot
export report
export make_fuel_dictionary
export fuel_plot

import Dates
import TimeSeries
import RecipesBase

#I/O Imports
import Colors
import DataFrames
import Feather
import Weave
import PowerSystems
import InfrastructureSystems

const PSY = PowerSystems
const IS = InfrastructureSystems
include("plot_results.jl")
include("fuel_results.jl")
include("plot_recipes.jl")
include("call_plots.jl")
include("make_report.jl")

end


text:
module Lagrangian
using Compat, JuMP

export
# Solve for the Lagrangian duals
lagrangian_method!,
# Relax, solve, then store
lagrangiansolve!,
# Structs to implement methods
LevelMethod,
SubgradientMethod,
KelleyMethod,
# Structs to solve the primal
LinearProgram,
# Each method needs to know about the problem being solved
LinearProgramData,
# Utilities
getslack,
isclose, closetozero,
Absolute, Relative, Unit

abstract type AbstractLagrangianMethod end
abstract type AbstractProblemClass end
struct LinearProgram <: AbstractProblemClass end

mutable struct LinearProgramData{M<:AbstractLagrangianMethod, C<:AbstractProblemClass}
    obj::QuadExpr                                # objective
    constraints::Vector{<:ConstraintRef} # constraints being relaxed
    relaxed_bounds::Vector{Float64}             # bounds on constraints when we relax them
    senses::Vector{Symbol}                      # we will cache the sense of constraints
    slacks::Vector{AffExpr}                     # also cache Ax-b
    old_bound::Vector{Float64}                  # cache before relaxing constraints
    method::M                                   # method parameters
    pc::C                                       # problem class
end

include("utils.jl")
include("levelmethod.jl")
include("subgradient.jl")
include("kelley.jl")

"""
    LinearProgramData(obj::QuadExpr, constraints::Vector{<:ConstraintRef}, relaxed_bounds::Vector{Float64}; method=LevelMethod(), problem_class=LinearProgram())

Creates a `LinearProgramData` object for calling `lagrangiansolve!`.

# Arguments
* obj:               Objective of the linear program without relaxation.
* constraints:       A vector of type `JuMP.ConstraintRef` with the contraints to be relaxed.
* relaxed_bounds:    A RHS that will replace the RHS of the constraint being relaxed (choose a big number for ≤ or = constraints, and something small for ≥).
* method:            Solving parameters, of type `AbstractLagrangianMethod`.
* problem_class:     To overload how the primal problem is solved.
"""
function LinearProgramData(obj::QuadExpr, constraints::Vector{<:ConstraintRef}, relaxed_bounds::Vector{Float64}; method=LevelMethod(), problem_class=LinearProgram())
    @assert length(constraints) == length(relaxed_bounds)
    LinearProgramData(obj, constraints,
        relaxed_bounds,
        sense.(constraints),
        getslack.(constraints),
        zeros(Float64, length(constraints)),
        method,
        problem_class
    )
end

function relaxandcache!(l::LinearProgramData, m::JuMP.Model)
    for (i, c) in enumerate(l.constraints)
        con = m.linconstr[c.idx]
        if l.senses[i] == :le
            l.old_bound[i] = con.ub
            con.ub = l.relaxed_bounds[i]
        elseif l.senses[i] == :ge
            l.old_bound[i] = con.lb
            con.lb = l.relaxed_bounds[i]
        else
            l.old_bound[i] = con.lb
            con.ub =  l.relaxed_bounds[i]
            con.lb = -Inf
        end
    end
end
function recover!(l::LinearProgramData, m::JuMP.Model, π::Vector{Float64})
    for (i, c) in enumerate(l.constraints)
        con = m.linconstr[c.idx]
        m.linconstrDuals[c.idx] = -π[i]
        if l.senses[i] == :le
            con.ub = l.old_bound[i]
        elseif l.senses[i] == :ge
            con.lb = l.old_bound[i]
        else
            con.lb = con.ub = l.old_bound[i]
        end
    end
end

"""
    lagrangian_method!{M<:AbstractLagrangianMethod}(l::LinearProgramData{M}, m::JuMP.Model, π::Vector{Float64})

Solve the Lagrangian dual problem.
"""
lagrangian_method!(l::LinearProgramData{M}, m::JuMP.Model, π::Vector{Float64}) where {M<:AbstractLagrangianMethod} = error("No solve method defined.")

"""
    lagrangiansolve!{M<:AbstractLagrangeinMethod}(l::LinearProgramData{M}, m::JuMP.Model, π::Vector{Float64})

Relax the primal, solve the Lagrangian dual, reset the primal.
"""
function lagrangiansolve!(l::LinearProgramData, m::JuMP.Model, π::Vector{Float64})

    @assert length(π) == length(l.constraints)
    @assert !isempty(l.constraints)

    # Relax bounds, cache old bounds
    relaxandcache!(l, m)

    # Solve
    status, obj = lagrangian_method!(l, m, π)

    # Set objective
    m.objVal = obj
    # JuMP.setobjective(m, getobjectivesense(m), l.obj) # bypass quad terms prep
    m.obj = l.obj

    # Set duals and reset constraint bounds
    recover!(l, m, π)

    return status, obj
end

end # module


text:
# Determine the index of the first Pkg REPL command
const first_cmd_idx = something(findfirst(x -> !startswith(x, "--"), ARGS), length(ARGS)+1)
const JLPKG_ARGS = ARGS[1:first_cmd_idx-1]
const PKG_REPL_ARGS = ARGS[first_cmd_idx:end]

# Parse --julia option (not supported on Windows)
if !Sys.iswindows()
    r = r"^--julia=(.+)$"
    idx = findlast(x -> match(r, x) !== nothing, JLPKG_ARGS)
    if idx !== nothing
        julia = match(r, JLPKG_ARGS[idx]).captures[1]
        deleteat!(ARGS, findall(x -> match(r, x) !== nothing, JLPKG_ARGS))
        f = @__FILE__ # JuliaLang/julia #28188
        cmd = Base.julia_cmd()
        cmd.exec[1] = julia # swap out the executable
        filter!(x -> !startswith(x, "-J"), cmd.exec) # filter out incompatible sysimg
        push!(cmd.exec, "--color=$(Base.have_color ? "yes" : "no")")
        pipe = pipeline(`$(cmd) $(f) $(ARGS)`; stdout=stdout, stderr=stderr)
        try
            exit(!success(pipe))
        catch err
            printstyled(stderr, "Error: "; bold=true, color=:red)
            io = IOBuffer()
            showerror(io, err)
            printstyled(stderr, String(take!(io)), '\n'; color=:red)
            exit(1)
        end
    end
end

# Parse --version option
if "--version" in JLPKG_ARGS
    println(stdout, "jlpkg version 1.1.1, julia version $(VERSION)")
    exit(0)
end

function isvalid(arg)
    return arg == "--update" || arg == "--version" ||
           arg == "--help" || startswith(arg, "--project")
end

# Check input and parse --help option
if isempty(ARGS) || isempty(PKG_REPL_ARGS) || "--help" in JLPKG_ARGS ||
        findfirst(!isvalid, JLPKG_ARGS) !== nothing
    exit_code = 0
    io = stdout
    if !("--help" in JLPKG_ARGS)
        io = stderr
        if isempty(ARGS)
            println(io, "No input arguments, showing help:\n")
            exit_code = 1
        elseif (idx = findfirst(!isvalid, JLPKG_ARGS); idx !== nothing)
            println(io, "Invalid argument `$(JLPKG_ARGS[idx])`, showing help:\n")
            exit_code = 1
        elseif isempty(PKG_REPL_ARGS)
            println(io, "No Pkg REPL arguments, showing help:\n")
            exit_code = 1
        end
    end
    # Print help
    printstyled(io, "NAME\n"; bold=true)
    println(io, """
           jlpkg - command line interface (CLI) to Pkg, Julia's package manager
    """)
    printstyled(io, "SYNOPSIS\n"; bold=true)
    println(io, """
           jlpkg [--options] <pkg-cmds>...
    """)
    printstyled(io, "OPTIONS\n"; bold=true)
    println(io, """
           <pkg-cmds>...
               Commands to the Pkg REPL mode. Execute the `help` command
               (e.g. `jlpkg help`) to list all available commands, and execute
               `jlpkg ?cmd` (or `jlpkg help cmd`) to show detailed help for a
               specific command. See https://julialang.github.io/Pkg.jl/v1/repl/
               for documentation of the syntax and the available commands.

           --project[=<path>]
               Set the home project/environment.
               Equivalent to Julia's `--project` switch.

           --julia=<path>
               Specify path to, or command for, the Julia executable.
               Overrides the executable set when installing the CLI.

           --update
               Allow the subsequent commands to update package registries.

           --version
               Show jlpkg and julia version numbers.

           --help
               Show this message.
    """)
    printstyled(io, "EXAMPLES\n"; bold=true)
    print(io, """
           · Add the Example package to the package environment located at `path`:
               \$ jlpkg --project=path add Example

           · Update package registries and add the JSON package:
               \$ jlpkg --update add JSON

           · Show the help for the `add` Pkg REPL command:
               \$ jlpkg ?add
    """)
    exit(exit_code)
end

# Parse --project option
let
    r = r"^--project(=(.+))?$"
    idx = findlast(x->match(r, x) !== nothing, JLPKG_ARGS)
    if idx === nothing # --project not given
        project = get(ENV, "JULIA_PROJECT", nothing)
    else # --project given
        m = match(r, JLPKG_ARGS[idx[end]])
        if m.captures[2] === nothing
            project = "@."
        else # m.captures !== nothing
            project = m.captures[2]
        end
    end
    Base.HOME_PROJECT[] =
        project === nothing ? nothing :
        project == "" ? nothing :
        project == "@." ? Base.current_project() :
        abspath(expanduser(project))
end

# Load Pkg; circumvent user-modified LOAD_PATH
const LOAD_PATH = copy(Base.LOAD_PATH)
try
    push!(empty!(Base.LOAD_PATH), joinpath(Sys.STDLIB, "Pkg"))
    using Pkg
catch
    printstyled(stderr, "Error: "; bold=true, color=:red)
    printstyled(stderr, "could not load Pkg.\n"; color=:red)
    rethrow()
finally
    append!(empty!(Base.LOAD_PATH), LOAD_PATH)
end

# Parse --update option
let
    idx = findlast(==("--update"), JLPKG_ARGS)
    if idx === nothing
        Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
    end
end

# Swap out --compile=min and --optimize=0 if we are running tests
# since we don't want that to propagate to the test-subprocess,
# and packages expect, and should be, tested with default options.
if "test" in PKG_REPL_ARGS
    o = Base.JLOptions()
    o′ = Base.JLOptions((x === :compile_enabled ? Int8(1) :
        x === :opt_level ? Int8(2) :
        getfield(o, x) for x in fieldnames(Base.JLOptions))...)
    unsafe_store!(cglobal(:jl_options, Base.JLOptions), o′)
end

# Run Pkg REPL mode with PKG_REPL_ARGS
try
    Pkg.REPLMode.pkgstr(join(PKG_REPL_ARGS, " "))
    if !isempty(PKG_REPL_ARGS) && (PKG_REPL_ARGS[1] == "help" || startswith(PKG_REPL_ARGS[1], '?'))
        # The help command uses `display` which does not add a trailing \n
        println(stdout)
    end
    exit(0)
catch err
    if err isa Pkg.Types.PkgError
        printstyled(stderr, "PkgError: "; bold=true, color=:red)
    elseif err isa Pkg.Types.ResolverError
        printstyled(stderr, "ResolverError: "; bold=true, color=:red)
    else
        rethrow()
    end
    io = IOBuffer()
    showerror(io, err)
    printstyled(stderr, String(take!(io)), '\n'; color=:red)
    exit(1)
end