text:
using LinearAlgebra, SparseArrays,Statistics,DataFrames,DataFramesMeta

#To install Laplacians you need to put in Julia's Prompt "add Laplacians#master"
#For some reason Pkg.add() doesn't work with this.
#Pkg.add("MATLAB")
#using Laplacians
#using Pkg
#include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
using Test
using Parameters
using Random

run_full_test = get(ENV, "JLA_TESTCSV", "false")  == "true" ? true : false
pkg_dir = pkgdir(VarianceComponentsHDFE)

@testset "PruneNetwork" begin
    #Full Network
    y_full = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208; 0.372084 ;  0.958089]
    id_full = [1; 1; 2;2 ; 3; 4;4 ;5;5 ;6;6]
    firmid_full = [1;2;1;1;1;1;2;2;2;2;3]
    obs_id_full = collect(1:11)

    #Pruned Network
    y_p = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
    id_p = [1; 1; 2;2 ; 3; 4;4 ;5;5]
    firmid_p = [1;2;1;1;1;1;2;2;2]
    obs_id_p = collect(1:9)

    #Drop-Single Network
    y_lo = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
    id_lo = [1; 1; 2;2 ; 3;3 ;4;4]
    firmid_lo = [1;2;1;1;1;2;2;2]
    obs_id_lo1 = [1;2;3;4;6;7;8;9]
    obs_id_lo = collect(1:8)
    # @test find_connected_set(y_full,id_full,firmid_full; verbose=true) == (obs_id = obs_id_full , y = y_full , id = id_full, firmid = firmid_full )
    @test find_connected_set(y_full,id_full,firmid_full,  VCHDFESettings(print_level = 1)) == (obs_id = obs_id_full , y = y_full , first_id = id_full, second_id = firmid_full )

    @test prunning_connected_set(y_full,id_full,firmid_full, obs_id_full,  VCHDFESettings(print_level = 1)) == (obs_id = obs_id_p , y = y_p , first_id = id_p, second_id = firmid_p )

    @test drop_single_obs(y_p,id_p,firmid_p,obs_id_p) == (obs_id = obs_id_lo1 , y = y_lo , first_id = id_lo, second_id = firmid_lo )

end

@testset "Movers & Matches" begin
    #Drop-Single Network
    y_lo = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
    id_lo = [1; 1; 2;2 ; 3;3 ;4;4]
    firmid_lo = [1;2;1;1;1;2;2;2]
    obs_id_lo1 = [1;2;3;4;6;7;8;9]
    obs_id_lo = collect(1:8)

    @test  compute_matchid(firmid_lo,id_lo) == [1; 2 ; 3;3 ; 4; 5; 6; 6]
    @test  compute_movers(id_lo, firmid_lo) == (movers = Bool[1,1,0,0,1,1,0,0], T= [2, 2, 2, 2, 2, 2, 2, 2])
end

# @testset "Clustering" begin
#     id_lo = [1; 1; 2;2 ; 3;3 ;4;4]
#     obs_id_lo = collect(1:8)

#     #Leave-Out Network Matchid
#     match_id = [1; 2 ; 3;3 ; 4; 5; 6; 6]
# # Xtest = sparse([1;2;3;4;5;6;7;8;1;3;4;5],[1;1;2;2;3;3;4;4;5;5;5;5],[1.0;1.0;1.0;1.0;1.0;1.0;1.0;1.0;-1.0;-1.0;-1.0;-1.0])
# # xxtest = Xtest'*Xtest
# # compute_sol = approxcholSolver(xxtest;verbose=true)
# # settings_exact = VCHDFESettings(leverage_algorithm = ExactAlgorithm(), first_id_effects=true, cov_effects=true)
# # settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=3), first_id_effects=true, cov_effects=true)

#     @test check_clustering(id_lo).nnz_2 == 12

#     @test check_clustering(match_id).nnz_2 ==10

#     @test index_constr(obs_id_lo,id_lo, obs_id_lo) == [obs_id_lo  obs_id_lo  obs_id_lo   obs_id_lo]

#     @test index_constr(obs_id_lo,id_lo, match_id) == [obs_id_lo  obs_id_lo  [1;2;3;3;4;5;6;6]   obs_id_lo]

# end

@testset "Solvers" begin
    Random.seed!(1234)
    Xtest = sparse([1;2;3;4;5;6;7;8;1;3;4;5],[1;1;2;2;3;3;4;4;5;5;5;5],[1.0;1.0;1.0;1.0;1.0;1.0;1.0;1.0;-1.0;-1.0;-1.0;-1.0])
    xxtest = Xtest'*Xtest
    compute_sol = approxcholSolver(xxtest;verbose=true)

    rademacher1 = rand(1,8)  .> 0.5
    rademacher1 = rademacher1 - (rademacher1 .== 0)

    @test compute_sol( [Xtest[1,:]...] ;verbose=true) ≈ [ 0.25, -0.5, -0.25, 0.0, -0.5]
    @test compute_sol([rademacher1*Xtest...];verbose=true) ≈  [-0.5, 0.0, -1.5, 0.0, -1.0]
end

const K=0
@testset "LambdaMatrices" begin
    Xtest = sparse([1;2;3;4;5;6;7;8;1;3;4;5],[1;1;2;2;3;3;4;4;5;5;5;5],[1.0;1.0;1.0;1.0;1.0;1.0;1.0;1.0;-1.0;-1.0;-1.0;-1.0])
    y_lo = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.730622 ; 0.619239; 0.753429; 0.0863208]
    id_lo = [1; 1; 2;2 ; 3;3 ;4;4]
    firmid_lo = [1;2;1;1;1;2;2;2]
    obs_id_lo1 = [1;2;3;4;6;7;8;9]
    obs_id_lo = collect(1:8)
    match_id = [1; 2 ; 3;3 ; 4; 5; 6; 6]

    settings_exact = VCHDFESettings(leverage_algorithm = ExactAlgorithm(), first_id_effects=true, cov_effects=true)
    settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=3), first_id_effects=true, cov_effects=true)

    # @test    do_Pii(Xtest,obs_id_lo) == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.75,0.75,0.5,0.5,0.75,0.75,0.5,0.5] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_P == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.75,0.75,0.5,0.5,0.75,0.75,0.5,0.5] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_second == sparse( [1,2,5,6], [1,2,5,6], [0.5,0.5,0.5,0.5],8,8)
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_first == sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.625,0.625,0.375,0.375,0.625,0.625,0.375,0.375] )
    @test    eff_res(settings_exact.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_exact).Lambda_B_cov == sparse( [1,2,5,6], [1,2,5,6], [-0.25,-0.25,-0.25,-0.25],8,8)

    # For now, just test that this executes
    # @test eff_res(settings_JLA.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_JLA).Lambda_P != nothing

    # Original test
    # @test    isapprox(eff_res(settings_JLA.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_JLA).Lambda_P, sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.666667,0.666667,0.5,0.5,0.666667,0.666667,0.5,0.5] ), atol = 1e-4)

    # @test    isapprox(eff_res(settings_JLA.leverage_algorithm, Xtest,id_lo,firmid_lo,match_id, K, settings_JLA).Lambda_P, sparse( [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8], [0.99,0.75,0.5,0.5,0.75,0.75,0.5,0.5] ), atol = 1e-4)
end

@testset "A throughout Test" begin
    y_full = [0.146297; 0.29686 ;  0.54344; 0.432677 ; 0.464866 ; 0.730622 ; 0.619239; 0.753429; 0.0863208; 0.372084 ;  0.958089]
    id_full = [1; 1; 2;2 ; 3; 4;4 ;5;5 ;6;6]
    firmid_full = [1;2;1;1;1;1;2;2;2;2;3]
    obs_id_full = collect(1:11)
    @testset "Exact Algorithm, person_effects = true, cov_effects = true" begin
        settings_exact = VCHDFESettings(leverage_algorithm = ExactAlgorithm(), first_id_effects=true, cov_effects=true)

        @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y_full, id_full, firmid_full, settings_exact, nothing)

    # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y,first_id,second_id,controls,settings)
        @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_exact)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y_full,id_full,firmid_full,nothing,settings_exact)

        # println("θ_second: $θ_second \n θ_first:$θ_first \n, θCOV: $θCOV \n obs: $obs \n β: $β \n Dalpha: $Dalpha \n Fpsi: $Fpsi \n Pii: $Pii \n Bii_first: $Bii_first \n Bii_second: $Bii_second \n Bii_cov: $Bii_cov \n")
        @test β ≈ [0.2313735, 0.5076485000000001, 0.6847255, 0.4198749, 0.019589999999999996]
        @test [θ_second, θ_first, θCOV] ≈ [-0.004178833653678571, 0.0036744476824373262, 0.0018986001519821427]
        @test Dalpha ≈  [0.2313735, 0.2313735, 0.5076485000000001, 0.5076485000000001, 0.6847255, 0.6847255, 0.4198749, 0.4198749]
        @test Fpsi ≈  [0.019589999999999996, 0.0, 0.019589999999999996, 0.019589999999999996, 0.019589999999999996, 0.0, 0.0, 0.0]
        @test Pii ≈ diag(diagm([0.75, 0.75, 0.5, 0.5, 0.75, 0.75, 0.5, 0.5]))
        @test Bii_first ≈ diag(diagm([0.625, 0.625, 0.375, 0.375, 0.625, 0.625, 0.375, 0.375]))
        @test Bii_second ≈ diag(diagm([0.5, 0.5, 0, 0, 0.5, 0.5, 0, 0]))
        @test Bii_cov ≈ diag(diagm([-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0]))
    end
    @testset "Exact Algorithm, person_effects = false, cov_effects = false" begin
        settings_exact = VCHDFESettings(leverage_algorithm = ExactAlgorithm(), first_id_effects=false, cov_effects=false)

        @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y_full, id_full, firmid_full, settings_exact, nothing)

    # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y,first_id,second_id,controls,settings)
        @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_exact)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y_full,id_full,firmid_full,nothing,settings_exact)

        # println("θ_second: $θ_second \n θ_first:$θ_first \n, θCOV: $θCOV \n obs: $obs \n β: $β \n Dalpha: $Dalpha \n Fpsi: $Fpsi \n Pii: $Pii \n Bii_first: $Bii_first \n Bii_second: $Bii_second \n Bii_cov: $Bii_cov \n")
        @test β ≈ [0.2313735, 0.5076485000000001, 0.6847255, 0.4198749, 0.019589999999999996]
        @test θ_second ≈ -0.004178833653678571
        @test [θ_first, θCOV] == [nothing, nothing]
        @test Dalpha ≈  [0.2313735, 0.2313735, 0.5076485000000001, 0.5076485000000001, 0.6847255, 0.6847255, 0.4198749, 0.4198749]
        @test Fpsi ≈  [0.019589999999999996, 0.0, 0.019589999999999996, 0.019589999999999996, 0.019589999999999996, 0.0, 0.0, 0.0]
        @test Pii ≈ diag(diagm([0.75, 0.75, 0.5, 0.5, 0.75, 0.75, 0.5, 0.5]))
        @test Bii_first == nothing
        @test Bii_second ≈ diag(diagm([0.5, 0.5, 0, 0, 0.5, 0.5, 0, 0]))
        @test isequal(Bii_cov, nothing)
    end
    @testset "JLA Algorithm with three simulations, person and cov effects = true" begin
        Random.seed!(1234)
        settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=3), first_id_effects=true, cov_effects=true)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y_full,id_full,firmid_full,nothing,settings_JLA)

        @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y_full, id_full, firmid_full, settings_JLA, nothing)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y,first_id,second_id,controls,settings)
        @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_JLA)

        println("θ_second: $θ_second \n θ_first:$θ_first \n, θCOV: $θCOV \n obs: $obs \n β: $β \n Dalpha: $Dalpha \n Fpsi: $Fpsi \n Pii: $Pii \n Bii_first: $Bii_first \n Bii_second: $Bii_second \n Bii_cov: $Bii_cov \n")
        # @test β ≈ [0.2313735, 0.5076485000000001, 0.6847255, 0.4198749, 0.019589999999999996]
        @test isapprox([θ_second, θ_first, θCOV] ,  [-0.17707507536052222, -0.17145102440488802, 0.10101283333196641], atol=1e-1) #Rewrite this

        @test Pii ≈ diag(diagm([0.99, 0.99, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]))
        @test isapprox(Bii_first, diag(diagm([0.770833, 1.10417, 0.375, 0.375, 0.604167, 0.270833, 0.375, 0.375])), atol = 1e-4)
        @test isapprox(Bii_second, diag(diagm([1.41667, 1.41667, 0, 0, 1.41667, 1.41667, 0, 0])), atol = 1e-4)
        @test isapprox(Bii_cov, diag(diagm([-0.625, -0.708333, 0, 0, -0.791667, -0.541667, 0, 0])), atol = 1e-4)
    end
    @testset "JLA Algorithm with three simulations, person and cov effects = false" begin
        Random.seed!(1234)
        settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=3), first_id_effects=false, cov_effects=false)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y_full,id_full,firmid_full,nothing,settings_JLA)

        @unpack obs,  y  , first_id , second_id, controls = get_leave_one_out_set(y_full, id_full, firmid_full, settings_JLA, nothing)

        # @unpack θ_first, θ_second, θCOV, obs, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = compute_whole(y,first_id,second_id,controls,settings)
        @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,controls,settings_JLA)

        println("θ_second: $θ_second \n θ_first:$θ_first \n, θCOV: $θCOV \n obs: $obs \n β: $β \n Dalpha: $Dalpha \n Fpsi: $Fpsi \n Pii: $Pii \n Bii_first: $Bii_first \n Bii_second: $Bii_second \n Bii_cov: $Bii_cov \n")
        # @test β ≈ [0.2313735, 0.5076485000000001, 0.6847255, 0.4198749, 0.019589999999999996]
        @test isapprox(θ_second ,  -0.17707507536052222, atol=1e-2) #Rewrite this

        @test Pii ≈ diag(diagm([0.99, 0.99, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]))
        @test isequal(Bii_first, nothing)
        @test isapprox(Bii_second, diag(diagm([1.41667, 1.41667, 0, 0, 1.41667, 1.41667, 0, 0])), atol = 1e-4)
        @test isequal(Bii_cov,nothing)

    end
end

if run_full_test
    using CSV
    @testset "JLA Algorithm on Test.csv" begin
        Random.seed!(1234)
        settings_JLA = VCHDFESettings(leverage_algorithm = JLAAlgorithm(num_simulations=3155), first_id_effects=true, cov_effects=true)
        data = CSV.read(joinpath(pkg_dir,"test.csv"); header=false)
        id = data[:,1]
        firmid = data[:,2]
        y = data[:,4]

        @unpack  obs,  y  , first_id , second_id, controls  = get_leave_one_out_set(y,id, firmid,settings_JLA, nothing)

        @unpack θ_first, θ_second, θCOV, β, Dalpha, Fpsi, Pii, Bii_first, Bii_second, Bii_cov = leave_out_estimation(y,first_id,second_id,nothing,settings_JLA)

        @test isapprox([θ_second, θ_first, θCOV],  [0.010505964725588434 , 0.08534310058596604, 0.004550851841563358],atol=1e-2)

    end

end


text:
using Distributed
if length(workers())==1
   addprocs(14)
end
using Random, Test
@everywhere using Gurobi, JuMP, JuDGE

include("solvers/setup_gurobi.jl")

function knapsack_parallel(parallel,check)
   Random.seed!(50)
   # how many investments?
   numinvest = 5;

   # number of items to pick from in the knapsack?
   numitems = 1000

   # size of tree?
   degree = 6
   height = 4

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = ([1,1.8,3.5,6.8,13.5])*(1-((i-1)/(totalnodes*1.2)))
   end

   # investvol = [40,45,50,70]
   investvol = [1,2,4,8,16]
   initialcap = 0

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems))*2) + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*2# + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      #println("Formulating subproblem for node "*node.name)
      model = Model()
      @expansion(model, bag[1:numinvest], Bin)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver, parallel=parallel, check=check, sp_solver=JuDGE_SP_Solver)

   return num_variables(judy.master_problem)
end

function knapsack_advanced(seed::Int64, tree_size::Tuple{Int64,Int64}, numitems::Int64, custom_solve_settings::Symbol)
   Random.seed!(seed)
   # how many investments?
   numinvest = 6;

   # size of tree?
   degree = tree_size[2]
   height = tree_size[1]

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = ([1,1.8,3.5,6.8,13.5,25])*(1-((i-1)/(totalnodes*1.2)))
   end

   # investvol = [40,45,50,70]
   investvol = [1,2,4,8,16,32]
   initialcap = 0

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems))*2) + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*2# + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      print("Formulating subproblem for node ")
      println(node.name)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag[1:numinvest], Bin)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      set_optimizer_attribute(model,"MIPGap",0.0005)
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)

   if custom_solve_settings==:callbacks
      function oa(judge::JuDGEModel,stalled::Bool,initialize::Bool)
         function slowed()
      	  if judge.log[end].iter>3
      		  if judge.log[end].rlx_abs/judge.log[end-3].rlx_abs>0.8
      			  return true
      		  end
      	  end
      	  false
         end

         function relgap_factor(probability)
            if judge.ext[:optimizer_settings][:level]==4
               0.0001
            else
               10^-4/probability*20^(4-judge.ext[:optimizer_settings][:level])
            end
         end

         if initialize
            if !haskey(judge.ext[:optimizer_settings],:level)
               judge.ext[:optimizer_settings][:level]=1
            end

            println("\u1b[1F\u1b[130GInitializing optimizer settings at level "*string(judge.ext[:optimizer_settings][:level]))

            for (node,sp) in judge.sub_problems
              set_optimizer_attribute(sp,"MIPGap",relgap_factor(judge.probabilities[node]))
            end
            return false
         elseif judge.ext[:optimizer_settings][:level]==4
   	   	if stalled
   	   		return true
   		   else
   			   return false
   		   end
         elseif stalled || slowed()
      	  judge.ext[:optimizer_settings][:level]+=1

      	  for (node,sp) in judge.sub_problems
      		 set_optimizer_attribute(sp,"MIPGap",relgap_factor(judge.probabilities[node]))
      	  end

      	  if judge.ext[:optimizer_settings][:level]==4
      		 println("\u1b[1F\u1b[130GFinal optimizer settings")
      	  else
      		 println("\u1b[1F\u1b[130GIncremented optimizer settings to level "*string(judge.ext[:optimizer_settings][:level]))
      	  end
         end

         return false
      end

      function mp_cb(judge::JuDGEModel,abs_tol::Float64,rel_tol::Float64)
         function earlytermination(cb_data, cb_where)
         	if cb_where == GRB_CB_MIP
         		objbst = Ref{Cint}()
         		GRBcbget(cb_data,cb_where,GRB_CB_MIP_OBJBST,objbst)
         		objbnd = Ref{Cint}()
         		GRBcbget(cb_data,cb_where,GRB_CB_MIP_OBJBND,objbnd)
         		runtime = Ref{Cint}()
         		GRBcbget(cb_data,cb_where,GRB_CB_RUNTIME,runtime)
         		JuDGE.printright("Incumbent: "*string(objbst[])*" Bound: "*string(objbnd[])*" Runtime: "*string(runtime[]))
         		if (objbnd[]>judge.bounds.LB+abs_tol && objbnd[]>judge.bounds.LB+rel_tol*abs(judge.bounds.LB)) || objbst[]<judge.bounds.LB+abs_tol || objbst[]<judge.bounds.LB+rel_tol*abs(judge.bounds.LB)
         			GRBterminate(backend(judge.master_problem).optimizer.model.inner)
         		end
         	end
         	return
         end
      	MOI.set(judge.master_problem, Gurobi.CallbackFunction(), earlytermination)
      end

      judy=JuDGE.branch_and_price(judy,termination=Termination(rlx_reltol=10^-6,reltol=5*10^-3,inttol=10^-6),
                  optimizer_attributes=oa,mp_callback=nothing,max_no_int=5)
   elseif custom_solve_settings==:partialpricing
      blocks = JuDGE.get_groups(mytree,combine=height-1)
      judy=JuDGE.branch_and_price(judy,termination=Termination(rlx_reltol=10^-6,reltol=5*10^-3,inttol=10^-6),max_no_int=5,blocks=blocks)
   else
      judy=JuDGE.branch_and_price(judy,termination=Termination(rlx_reltol=10^-6,reltol=5*10^-3,inttol=10^-6),max_no_int=5)
   end

   println("Objective: "*string(objective_value(judy.master_problem)))
   JuDGE.print_expansions(judy,format=format_output,inttol=10^-6)

   return objective_value(judy.master_problem)
end

# Solve without callbacks for medium tree and 100 items
@time @test knapsack_advanced(50,(4,3),100,:standard) ≈ -14.05 atol = 2e-2

# Solve with callbacks for medium tree and 100 items
@time @test knapsack_advanced(50,(4,3),100,:callbacks) ≈ -14.05 atol = 2e-2

# Solve with parial pricing for medium tree and 100 items
@time @test knapsack_advanced(50,(4,3),100,:partialpricing) ≈ -14.05 atol = 2e-2

# Formulate with checks in serial mode
@time @test knapsack_parallel(false,true) ≈ 9071 atol = 1e-3

# Formulate with checks in parallel mode (much slower, so probably shouldn't use)
@time @test knapsack_parallel(true,true) ≈ 9071 atol = 1e-3

# Formulate without checks in parallel mode
@time @test knapsack_parallel(true,false) ≈ 9071 atol = 1e-3

ext: jl


text:
using SnoopCompile, AssociatedLegendrePolynomials
inf_timing = @snoopi tmin=0.005 begin
    # Invoking aliased wrappers forces inference of underlying implmentation
    for T in (Float64, Float32), I in (Int,)
        s = T(0.5)
        v = collect(range(T(-0.5), T(0.5), length=5))
        Ptab = LegendreUnitCoeff{T}(I(2))
        Λtab = LegendreSphereCoeff{T}(I(2))
        # Scalar argument, single (ℓ,m)
        Plm(I(2), I(2), s)
        Ptab(I(2), I(2), s)
        λlm(I(2), I(2), s)
        Λtab(I(2), I(2), s)
        # Scalar argument, vector ℓ, fixed m
        Plm(I(0):I(2), I(2), s)
        Ptab(I(0):I(2), I(2), s)
        λlm(I(0):I(2), I(2), s)
        Λtab(I(0):I(2), I(2), s)
        # Scalar argument, matrix over all ℓ, m
        Plm(I(0):I(2), I(0):I(2), s)
        Ptab(I(0):I(2), I(0):I(2), s)
        λlm(I(0):I(2), I(0):I(2), s)
        Λtab(I(0):I(2), I(0):I(2), s)

        # Broadcasted operations:
        #  Scalar argument, vector output
        Plm.(I(0):I(2), I(2), s)
        Ptab.(I(0):I(2), I(2), s)
        λlm.(I(0):I(2), I(2), s)
        Λtab.(I(0):I(2), I(2), s)
        #  Scalar argument, matrix output
        Plm.(I(0):I(2), I(0):I(2), s)
        Ptab.(I(0):I(2), I(0):I(2), s)
        λlm.(I(0):I(2), I(0):I(2), s)
        Λtab.(I(0):I(2), I(0):I(2), s)
        #  Vector argument, vector output
        Plm.(I(0):I(2), I(2), v)
        Ptab.(I(0):I(2), I(2), v)
        λlm.(I(0):I(2), I(2), v)
        Λtab.(I(0):I(2), I(2), v)
        #  Vector argument, matrix output
        Plm.(I(0):I(2), I(0):I(2), v)
        Ptab.(I(0):I(2), I(0):I(2), v)
        λlm.(I(0):I(2), I(0):I(2), v)
        Ptab.(I(0):I(2), I(0):I(2), v)

        # Force precompilation of displaying functions
        display(Ptab)
        display(Λtab)
    end
end

pc = SnoopCompile.parcel(inf_timing)
pc = filter!(p -> p.first === :AssociatedLegendrePolynomials, pc)
sort!(pc[:AssociatedLegendrePolynomials])

outfile = joinpath(tempdir(), "precompile")
println("Saving precompile scripts to $outfile")
SnoopCompile.write(outfile, pc)


text:
""
function constraint_thermal_limit_from_soft(pm::_PM.AbstractPowerModel, n::Int, f_idx, rate_a)
    l,i,j = f_idx
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_fr^2 + q_fr^2 <= (rate_a + sm_slack)^2)
end

""
function constraint_thermal_limit_to_soft(pm::_PM.AbstractPowerModel, n::Int, t_idx, rate_a)
    l,i,j = t_idx
    p_to = var(pm, n, :p, t_idx)
    q_to = var(pm, n, :q, t_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_to^2 + q_to^2 <= (rate_a + sm_slack)^2)
end

""
function constraint_gen_power_real_link(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int)
    pg_1 = var(pm, n_1, :pg, i)
    pg_2 = var(pm, n_2, :pg, i)

    JuMP.@constraint(pm.model, pg_1 == pg_2)
end

""
function constraint_gen_power_real_response(pm::_PM.AbstractPowerModel, nw_1::Int, nw_2::Int, i::Int, alpha)
    pg_base = var(pm, :pg, i, nw=nw_1)
    pg = var(pm, :pg, i, nw=nw_2)
    delta = var(pm, :delta, nw=nw_2)

    @constraint(pm.model, pg == pg_base + alpha*delta)
end

""
function constraint_gen_power_real_deviation(pm::_PM.AbstractPowerModel, n::Int, i, pg)
    pg_var = var(pm, n, :pg, i)
    pg_delta = var(pm, n, :pg_delta, i)

    JuMP.@constraint(pm.model, pg_var == pg + pg_delta)
end

""
function constraint_branch_contingency_ptdf_thermal_limit_from(pm::_PM.AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate)
end

""
function constraint_branch_contingency_ptdf_thermal_limit_to(pm::_PM.AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, -sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate)
end

""
function constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm::_PM.AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate + var(pm, :branch_cont_flow_vio, i))
end

""
function constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm::_PM.AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, -sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate + var(pm, :branch_cont_flow_vio, i))
end


text:
###############################################################################
#
# Calculate annual profit
#
###############################################################################
"""
    annual_profit(
                node::SPACSimple{FT},
                weather::Array{FT,2},
                Δt::FT = FT(1)
    ) where {FT<:AbstractFloat}

Calculate the profit in the growing season so as to optimize leaf investment,
    given
- `node` [`SPACSimple`] type struct
- `weather` Weather profile in a growing season
- `Δt` Time period in `[h]`
"""
function annual_profit(
            node::SPACSimple{FT},
            weather::Array{FT,2},
            Δt::FT = FT(1)
) where {FT<:AbstractFloat}
    # 0. unpack required values
    #    make sure construction cost must be postive
    @unpack c_cons, c_vmax, elevation, gaba, laba, latitude, vtoj = node;
    cons = laba * (c_cons + node.ps.PSM.v_cmax25 * c_vmax);

    # if cons > 0
    if cons > 0
        # 1. update the environmental constants based on node geographycal info
        ratio            = atmospheric_pressure_ratio(elevation);
        node.envir.P_AIR = ratio * P_ATM(FT);
        node.envir.P_O₂  = node.envir.P_AIR * FT(0.209);

        # 2. calculate the growing season canopy profit
        day  ::FT = FT(0);
        hour ::FT = FT(0);
        _tair::FT = FT(0);
        _dair::FT = FT(0);
        p_co2::FT = FT(0);
        r_all::FT = FT(0);
        wind ::FT = FT(0);
        rain ::FT = FT(0);
        gscp ::FT = FT(0);
        anet ::FT = FT(0);

        N = size(weather)[1]
        for i in 1:N
            # 2.1 read and update the hourly data
            day   = weather[i,2 ]
            hour  = weather[i,3 ]
            _tair = weather[i,7 ]
            _dair = weather[i,9 ]
            p_co2 = weather[i,10] * ratio
            r_all = weather[i,4 ]
            wind  = weather[i,6 ]
            rain  = weather[i,5 ]

            node.envir.t         = _tair + T_0(FT);
            node.envir.p_H₂O_sat = saturation_vapor_pressure( node.envir.t );
            node.envir.p_CO₂     = p_co2;
            node.envir.p_H₂O     = node.envir.p_H₂O_sat - _dair * 1000;
            node.envir.rh        = node.envir.p_H₂O / node.envir.p_H₂O_sat;
            node.envir.wind      = wind;

            # 2.2 if day time
            zenith = zenith_angle(latitude, day, hour, FT(0));
            if (r_all>0) & (zenith<=85)
                # 2.2.1 update the leaf partitioning
                big_leaf_partition!(node, zenith, r_all);
                @unpack frac_sh, frac_sl = node.container2L;

                # 2.2.2 optimize flows in each layer
                optimize_flows!(node);
                leaf_gas_exchange!(node, node.opt_f_sl, node.opt_f_sh);
                flow = node.opt_f_sl + node.opt_f_sh;
                anet = frac_sl * node.container2L.cont_sl.an + frac_sh * node.container2L.cont_sh.an;

                # 2.2.3 update drought history
                pressure_profile!(node.hs, node.p_soil, node.opt_f_sl, node.opt_f_sh, frac_sl);

            # 2.3 if night time
            else
                node.ps.t = max(200, leaf_temperature(node, r_all, FT(0)));
                photosystem_temperature_dependence!(node.ps.PSM, node.envir, node.ps.t);
                node.ps.p_H₂O_sat = saturation_vapor_pressure(node.ps.t);
                flow = FT(0);
                anet = -node.ps.PSM.r_d;
            end

            # 2.4 update soil moisture by converting flow to Kg per hour
            flow /= KG_H_2_MOL_S(FT);
            rain *= gaba * ρ_H₂O(FT) / 1000;
            soil_moisture!(node, flow-rain);

            # 2.5 add up the gscp
            gscp += anet
        end

        # 3. scale gscp to mol per year
        gscp *= laba * FT(1e-6) * 3600 * Δt;

        # 4. substract the construction costs
        gscp -= cons;

    else
        gscp = FT(-Inf);
    end

    return gscp
end


text:
using Images: otsu_threshold, imadjustintensity
using Unitful
using Unitful: μm
using OffsetArrays
using ProgressMeter

"""
    segment_cells(img, seed_channel, segment_channel, maskfunc)
"""
function segment_cells(img::AxisArray{T, 4},
                       seed_channel::Symbol,
                       segment_channel::Symbol,
                       maskfunc::Function) where {T}
    if !all(map(x->x in axisnames(img), [:time, :channel]))
        error("image requires :channel and :time axes")
    end
    tax = timeaxis(img)
    segments = Array{SegmentedImage{Array{Int64,2},Float32}}(undef, length(tax))
    ylen, xlen = size(img, Axis{:y}), size(img, Axis{:x})
    mask = Array{Bool}(undef, ylen, xlen)
    segslice = Array{T}(undef, ylen, xlen)
    seedslice = Array{T}(undef, ylen, xlen)
    seeds = Array{Bool}(undef, ylen, xlen)

    for (idx, tp) in enumerate(tax.val)

        # find seeds
        seedslice .= view(img, Axis{:channel}(seed_channel), Axis{:time}(tp))
        seeds .= seedslice .> otsu_threshold(seedslice)*1.1

        # find segmentation slice
        segslice .= view(img, Axis{:channel}(segment_channel), Axis{:time}(tp))

        # throw out seeds that are outside of segment areas
        seeds .&= (segslice .> 0.0)
        ImageMorphology.closing!(seeds)
        ImageMorphology.erode!(seeds)
        markers = label_components(seeds)

        # compute mask
        mask .= maskfunc(segslice)
        ImageMorphology.closing!(mask)

        segments[idx] = watershed(1 .- imadjustintensity(segslice), markers, mask=mask)
    end
    segments
end

"""
    build_tp_df(img, components; dist)

Build a `DataFrames.DataFrame` that is compatible with trackpys `link_df`
function. Needs to be converted to a `Pandas.DataFrame` before passing to
trackpy. `dist` is a 2-tuple of integers indicating the minimum and maximum
distance away in pixels from each cell to include in its local background
calculation.
"""
function build_tp_df(img::AxisArray{T1, 4},
                     components::AxisArray{T2, 3}; dist=(2, 10)) where {T1, T2 <: Integer}

    xstep = step(AxisArrays.axes(img, Axis{:x}).val)
    ystep = step(AxisArrays.axes(img, Axis{:y}).val)
    (xstep != ystep) && @warn "Different scaling for x and y axes is not supported"
    pixelarea = xstep * ystep
    particles = DataFrames.DataFrame[]

    @showprogress 1 "Computing..." for (idx, timepoint) in enumerate(timeaxis(img))
        component_slice = view(components, :, :, idx)

        # get all ids present in the current slice
        slice_ids = sort(unique(component_slice))
        # get the number of pixels for each id in the slice
        lengths = Images.component_lengths(component_slice)[slice_ids .+ 1]

        if eltype(pixelarea) <: Unitful.Area
            # convert boolean area to microns with the assumption that the pixel
            # space is the same in both x and y
            areas = round.(μm^2, lengths .* pixelarea, sigdigits=4)
            # filter out too large and too small particles
            correct_size = (10μm^2 .< areas .< 500μm^2)
        else
            areas = lengths
            correct_size = trues(length(lengths))
            # remove background
            correct_size[1] = false
        end

        ids = slice_ids[correct_size]
        # get the centroids for all ids and then select only the ids that appear
        # in the current slice and *also* have the correct size
        centroids = Images.component_centroids(component_slice)[slice_ids[correct_size] .+ 1]

        n = length(centroids)
        ys = map(f->f[1], centroids) # y corresponds to rows
        xs = map(f->f[2], centroids) # x corresponds to columns
        frames = fill(idx-1, n)

        # select the current time slice and enforce storage order to match
        # components
        slice = view(img, Axis{:y}(:), Axis{:x}(:), Axis{:channel}(:), Axis{:time}(timepoint))
        localities = get_localities(component_slice, ids, dist=dist)
        # dictionary of ids to areas
        data = OrderedDict(:x=>xs,
                           :y=>ys,
                           :frame=>frames,
                           :id=>ids,
                           :area=>areas[correct_size])

        cax = AxisArrays.axes(img, Axis{:channel})
        for c in cax
            channelslice = view(slice, Axis{:channel}(c))
            tfs = Float64[]
            medbkgs = Float64[]

            for (id, indices) in localities
                # total fluorescence is the sum of all signal in the actual
                # footprint of the object. We have to do a copy operation here
                # due to https://github.com/JuliaArrays/AxisArrays.jl/issues/179
                push!(tfs, sum(channelslice[component_slice .== id]))

                # median background is the median of background signal in the
                # locality of object
                push!(medbkgs, median(channelslice[indices]))
            end
            data[Symbol("tf_", c)] = tfs
            data[Symbol("medbkg_", c)] = medbkgs
        end
        push!(particles, DataFrames.DataFrame(data))
    end
    vcat(particles...)
end

function build_tp_df(img::AxisArray{T1, 3},
                     thresholds::AxisArray{T2, 3}; dist=(2, 10)) where {T1, T2 <: Bool}
    build_tp_df(AxisArray(reshape(img, size(img)..., 1),
                          AxisArrays.axes(img)..., Axis{:channel}([:slice])),
                thresholds;
                dist=dist
               )
end

function build_tp_df(img::AxisArray{T1, 4},
                     thresholds::AxisArray{T2, 3}; dist=(2, 10)) where {T1, T2 <: Bool}

    # we have to pass the underlying array due to
    # https://github.com/JuliaImages/ImageMorphology.jl/issues/21
    components = Images.label_components(thresholds.data, [axisdim(thresholds, Axis{:y}), axisdim(thresholds, Axis{:x})])
    build_tp_df(img, AxisArray(components, AxisArrays.axes(thresholds)); dist=dist)
end

function build_tp_df(img::AxisArray{T1, 3},
                     thresholds::AxisArray{T2, 3}; dist=(2, 10)) where {T1, T2 <: Integer}

    build_tp_df(AxisArray(reshape(img, size(img)..., 1),
                          AxisArrays.axes(img)..., Axis{:channel}([:slice])),
                thresholds;
                dist=dist
               )
end

"""
Given an `img` with at least `y`, `x`, and `t` axes and a 3 dimensional boolean
array, `thresholds`, in yxt order.
"""
function build_tp_df(img::AxisArray{T1, 4},
                     thresholds::BitArray{3}; dist=(2, 10)) where {T1}
    _axes = Tuple(AxisArrays.axes(img, Axis{ax}) for ax in (:y, :x, :time))
    build_tp_df(img, AxisArray(Bool.(thresholds), _axes...), dist=dist)
end

"""
    get_locality(foreground, cell_mask; dist=(mindist, maxdist))

Given a boolean matrix of the pixels belonging to a single cell,
`cell_mask`, and a boolean matrix of foreground pixels, `foreground`, this
function  identifies the local background ring around the cell that is at
minimum `mindist` away from every object (in pixels) and a maximum of `maxdist`
away from the target cell.
"""
function get_locality(foreground::AbstractArray{Bool, 2}, cellmask::AbstractArray{Bool, 2}; dist=(2, 10))
    all_localities = dist[1] .< distance_transform(feature_transform(foreground)) .< dist[2]
    locality = dist[1] .< distance_transform(feature_transform(cellmask)) .< dist[2]
    # only return true where this cells locality overlaps with other localities, this way we
    # insure that only foreground areas are counted in the locality
    locality .&= all_localities
end

"""
    get_localities(labels, ids; dist) -> Dict

Gets the local background areas given the output of `label_components` and a
list of objects whose areas should be found defined by `ids`. The local
background is defined as a ring around the object that is at
minimum `mindist` away from every object (in pixels) and a maximum of `maxdist`
away from the target object. The result is a dictionary mapping the object id to
all indices in its local background.

!!! warning
    It's really important that `labels` has all foreground objects labeled
    because otherwise they might inadvertently be included in an object's
    locality.
"""
function get_localities(labels::AbstractArray{Int, 2}, ids::Vector{Int}; dist=(2, 10))
    boxes = component_boxes(labels)
    # boxes always contains the background label as well, while there's no guarantee
    # that ids does as well
    (!(0 in ids)) && deleteat!(boxes, 1)

    # if we detect gaps in the numbering it's probably because some labels were
    # removed, which can make our foreground detection incorrect
    computed_ids = sort(unique(labels))[2:end]
    if computed_ids != collect(minimum(computed_ids):maximum(computed_ids))
        @warn "Do not remove values from `labels`, only from `ids`. Calc might be wrong!"
    end

    nx, ny = size(labels)
    localities = OrderedDict{Int, Vector{CartesianIndex{2}}}()

    allobjects = labels .> 0

    δ = dist[1] + dist[2]

    for id in ids
        # unpack bounding box
        (minx, miny), (maxx, maxy) = boxes[id]

        # extend window around bounding box by the max distance
        xrange = max(minx-δ, 1):min(maxx+δ, nx)
        yrange = max(miny-δ, 1):min(maxy+δ, ny)

        local_allobjects = view(allobjects, xrange, yrange)
        local_cellmask = view(labels, xrange, yrange) .== id
        @assert sum(local_cellmask) > 0 "No object of $id found"

        locality = get_locality(local_allobjects, local_cellmask; dist=dist)

        # use offset arrays to return CartesianIndices in the original image coordinates
        localities[id] = findall(OffsetArray(locality, xrange, yrange))
    end

    localities
end


text:
function p_bar_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_G_T /= K_n(cur_reactor)

  cur_G_T /= K_nT(cur_reactor)

  cur_G_T /= cur_reactor.T_bar

  cur_equation = Equation(-2, 0, 1, cur_G_T)

  cur_equation
end


text:
# TODO: Document this.
macro aggregate_vector(src, exprs...)
    Expr(
        :call,
        :AggregateVector,
        esc(src),
        Expr(
            :tuple,
            expression_expr.(exprs)...,
        ),
    )
end


text:
"""
# Stipple.Typography

Typography utilities for Stipple apps.
"""
module Typography

import Genie

"""
    `function header(args...; size::Int = 1, kwargs...)`

Generates code for HTML headers (H1 to H6) based on `size` which include styling for Twitter Bootstrap, the CSS
  framework used by Stipple.

### Example

```julia
julia> Typography.header("Hello", size = 2)
"<h2 class=\"text-h2\">Hello</h2>"
```
"""
function header(args...; size::Int = 1, kwargs...)
  1 <= size <= 6 || error("Invalid header size - expected 1:6")
  func = getproperty(Genie.Renderer.Html, Symbol("h$size"))
  func(class="text-h$size", args...; kwargs...) |> ParsedHTMLString
end

end


text:
using Revise, TopOpt

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force

problems = Any[PointLoadCantilever(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
            HalfMBB(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
            LBeam(Val{:Linear}, Float64),
            TieBeam(Val{:Quadratic}, Float64)]
problem_names = ["cantilever beam", "half MBB beam", "L-beam", "tie-beam"]
approx_objvals = [330.0, 175.0, 65.0, 1413.0]

i = 1

problem = problems[i]
# Parameter settings
V = 0.5 # volume fraction
xmin = 0.001 # minimum density
steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
reuse = false # adaptive penalty flag
convcriteria = MMA.KKTCriteria()
#penalty = TopOpt.PowerPenalty(1.0)
penalty = TopOpt.PowerPenalty(1.0)
pcont = Continuation(penalty, steps = steps, xmin = xmin, pmax = 5.0)

mma_options = options = MMA.Options(maxiter=1000, tol = MMA.Tolerances(kkttol = 1e-4))
if convcriteria isa MMA.KKTCriteria
    maxtol = 0.01 # maximum tolerance
    mintol = 0.001 # minimum tolerance

    b = log(mintol / maxtol) / steps
    a = maxtol / exp(b)
    mma_options_gen = TopOpt.MMAOptionsGen(steps = steps, initial_options = mma_options, kkttol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol))
else
    maxtol = 0.01 # maximum tolerance
    mintol = 0.0001 # minimum tolerance

    b = log(mintol / maxtol) / steps
    a = maxtol / exp(b)
    mma_options_gen = TopOpt.MMAOptionsGen(steps = steps, initial_options = mma_options, ftol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol))
end
csimp_options = TopOpt.CSIMPOptions(steps = steps,
                                    options_gen = mma_options_gen,
                                    p_gen = pcont,
                                    reuse = reuse
                                    )

# Define a finite element solver
solver = FEASolver(Displacement, Direct, problem, xmin = xmin, penalty = penalty)
# Define compliance objective
filtering = problem isa TopOptProblems.TieBeam ? false : true
obj = Objective(VolumeFunction(problem, solver))
constr = Constraint(TopOpt.GlobalStress(solver, 1.0), 1.0)
#cu_obj = TopOpt.cu(obj)
# Define volume constraint
# Define subproblem optimizer
#optimizer = MMAOptimizer{CPU}(cu_obj, constr, MMA.MMA87(),
#    ConjugateGradient(), maxiter=1000); optimizer.obj.fevals = 0
optimizer = MMAOptimizer{CPU}(obj, constr, MMA.MMA02(),
    ConjugateGradient(), options = mma_options, convcriteria = convcriteria); optimizer.obj.f.fevals = 0

# Define continuation SIMP optimizer
simp = SIMP(optimizer, penalty.p)

cont_simp = ContinuationSIMP(simp, steps, csimp_options)

# Solve
x0 = fill(1.0, length(solver.vars))
result = cont_simp(x0)


text:
abstract type MoistureBC end

"""
    Impermeable() :: MoistureBC

No moisture flux.
"""
struct Impermeable <: MoistureBC end
function atmos_moisture_boundary_state!(
    nf,
    bc_moisture::Impermeable,
    atmos,
    args...,
) end
function atmos_moisture_normal_boundary_flux_second_order!(
    nf,
    bc_moisture::Impermeable,
    atmos,
    args...,
) end

"""
    PrescribedMoistureFlux(fn) :: MoistureBC

Prescribe the net inward moisture flux across the boundary by `fn`, a function
with signature `fn(state, aux, t)`, returning the flux (in kg/m^2).
"""
struct PrescribedMoistureFlux{FN} <: MoistureBC
    fn::FN
end
function atmos_moisture_boundary_state!(
    nf,
    bc_moisture::PrescribedMoistureFlux,
    atmos,
    args...,
) end
function atmos_moisture_normal_boundary_flux_second_order!(
    nf,
    bc_moisture::PrescribedMoistureFlux,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diffusive⁻,
    hyperdiffusive⁻,
    aux⁻,
    state⁺,
    diffusive⁺,
    hyperdiffusive⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    nρd_q_tot = -bc_moisture.fn(state⁻, aux⁻, t)
    fluxᵀn.ρ += nρd_q_tot
    fluxᵀn.ρu += nρd_q_tot / state⁻.ρ .* state⁻.ρu
    # assumes EquilMoist
    fluxᵀn.moisture.ρq_tot += nρd_q_tot
end

"""
    BulkFormulaMoisture(fn) :: MoistureBC

Calculate the net inward moisture flux across the boundary using
the bulk formula. The drag coefficient is `C_q = fn_C_q(state, aux,
t, normu_int_tan)`. The surface q_tot at the boundary is `q_tot =
fn_q_tot(state, aux, t)`.

Return the flux (in kg m^-2 s^-1).
"""
struct BulkFormulaMoisture{FNX, FNM} <: MoistureBC
    fn_C_q::FNX
    fn_q_tot::FNM
end
function atmos_moisture_boundary_state!(
    nf,
    bc_moisture::BulkFormulaMoisture,
    atmos,
    args...,
) end
function atmos_moisture_normal_boundary_flux_second_order!(
    nf,
    bc_moisture::BulkFormulaMoisture,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diffusive⁻,
    hyperdiffusive⁻,
    aux⁻,
    state⁺,
    diffusive⁺,
    hyperdiffusive⁺,
    aux⁺,
    bctype,
    t,
    state_int⁻,
    diffusive_int⁻,
    aux_int⁻,
)

    u_int⁻ = state_int⁻.ρu / state_int⁻.ρ
    u_int⁻_tan = projection_tangential(atmos, aux_int⁻, u_int⁻)
    normu_int⁻_tan = norm(u_int⁻_tan)
    C_q = bc_moisture.fn_C_q(state⁻, aux⁻, t, normu_int⁻_tan)
    q_tot = bc_moisture.fn_q_tot(state⁻, aux⁻, t)
    q_tot_int = state_int⁻.moisture.ρq_tot / state_int⁻.ρ

    # TODO: use the correct density at the surface
    ρ_avg = average_density(state⁻.ρ, state_int⁻.ρ)
    # NOTE: difference from design docs since normal points outwards
    fluxᵀn.moisture.ρq_tot -= C_q * ρ_avg * normu_int⁻_tan * (q_tot - q_tot_int)
    fluxᵀn.ρ -= C_q * ρ_avg * normu_int⁻_tan * (q_tot - q_tot_int)
    fluxᵀn.ρu -= C_q * normu_int⁻_tan * (q_tot - q_tot_int) .* state_int⁻.ρu
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
abstract type PredictiveControlProblem end

struct PrimalNoCondensing <: PredictiveControlProblem end

struct PrimalFullCondensing <: PredictiveControlProblem end

struct DualFullCondensing <: PredictiveControlProblem end


text:
@system TemperatureDependence begin
    T: leaf_temperature ~ hold
    Tk(T): absolute_leaf_temperature ~ track(u"K")

    Tb: base_temperature => 25 ~ preserve(u"°C", parameter)
    Tbk(Tb): absolute_base_temperature ~ preserve(u"K")

    kT(T, Tk, Tb, Tbk; Ea(u"kJ/mol")): arrhenius_equation => begin
        exp(Ea * (T - Tb) / (u"R" * Tk * Tbk))
    end ~ call

    kTpeak(Tk, Tbk, kT; Ha(u"kJ/mol"), Hd(u"kJ/mol"), To(u"K")): peaked_function => begin
        R = u"R"
        ΔS = Hd / To + R * log(Ha / (Hd - Ha))
        kT(Ha) * (1 + exp((ΔS*Tbk - Hd) / (R*Tbk))) / (1 + exp((ΔS*Tk - Hd) / (R*Tk)))
    end ~ call

    Q10 => 2 ~ preserve(parameter)
    kTQ10(T, Tb, Q10): q10_rate => begin
        Q10^((T - Tb) / 10u"K")
    end ~ track
end

@system NitrogenDependence begin
    N: leaf_nitrogen_content ~ hold

    s => 2.9 ~ preserve(u"m^2/g", parameter)
    N0 => 0.25 ~ preserve(u"g/m^2", parameter)

    kN(N, s, N0): nitrogen_limited_rate => begin
        2 / (1 + exp(-s * (max(N0, N) - N0))) - 1
    end ~ track
end

@system CBase(TemperatureDependence, NitrogenDependence) begin
    Ci: intercellular_co2 ~ hold
    I2: effective_irradiance ~ hold
end


text:
using SpecialFunctions, Documenter

makedocs(
    modules = [SpecialFunctions],
    clean = false,
    format = :html,
    sitename = "SpecialFunctions.jl",
    authors = "Jeff Bezanson, Stefan Karpinski, Viral B. Shah, et al.",
    pages = [
        "Home" => "index.md",
        "Functions" => "special.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaMath/SpecialFunctions.jl.git",
    target = "build",
)


text:
module TestSubset

using Test, DataFrames, Statistics

const ≅ = isequal

@testset "subset and subset!" begin
    refdf = DataFrame(x = repeat(Any[true, false], 4),
                      y = repeat([true, false, missing, missing], 2),
                      z = repeat([1, 2, 3, 3], 2),
                      id = 1:8)

    for df in (copy(refdf), @view copy(refdf)[1:end-1, :])
        df2 = copy(df)
        @test subset(df, :x) ≅ filter(:x => identity, df)
        @test df ≅ df2
        @test subset(df, :x) isa DataFrame
        @test subset(df, :x, view=true) ≅ filter(:x => identity, df)
        @test subset(df, :x, view=true) isa SubDataFrame
        @test_throws ArgumentError subset(df, :y)
        @test_throws ArgumentError subset(df, :y, :x)
        @test subset(df, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :y, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(df, :x, :y, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :y, :x, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :x, :y, skipmissing=true, view=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(df, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(df, :x, :y, :id => ByRow(<(4)), skipmissing=true, view=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(df, :x, :id => ByRow(<(4))) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test subset(df, :x, :id => ByRow(<(4)), view=true) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test_throws ArgumentError subset(df)
        @test isempty(subset(df, :x, :x => ByRow(!)))
        @test_throws ArgumentError subset(df, :x => x -> false, :x => x -> missing)
        @test_throws ArgumentError subset(df, :x => x -> true, :x => x -> missing)
        @test_throws ArgumentError subset(df, :x => x -> true, :x => x -> 2)
    end

    for df in (copy(refdf), @view copy(refdf)[1:end-1, :]),
        gdf in (groupby(df, :z), groupby(df, :z)[[3, 2, 1]])
        df2 = copy(df)
        @test subset(gdf, :x) ≅ filter(:x => identity, df)
        @test df ≅ df2
        @test subset(gdf, :x) isa DataFrame
        @test subset(gdf, :x, ungroup=false) ≅
              groupby(filter(:x => identity, df), :z)
        @test subset(gdf, :x, ungroup=false) isa GroupedDataFrame{DataFrame}
        @test subset(gdf, :x, view=true) ≅ filter(:x => identity, df)
        @test subset(gdf, :x, view=true) isa SubDataFrame
        @test subset(gdf, :x, view=true, ungroup=false) ≅
              groupby(filter(:x => identity, df), :z)
        @test subset(gdf, :x, view=true, ungroup=false) isa GroupedDataFrame{<:SubDataFrame}
        @test_throws ArgumentError subset(gdf, :y)
        @test_throws ArgumentError subset(gdf, :y, :x)
        @test subset(gdf, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, :y, skipmissing=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :y, :y, skipmissing=true, view=true) ≅ filter(:y => x -> x === true, df)
        @test subset(gdf, :x, :y, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :y, :x, skipmissing=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :x, :y, skipmissing=true, view=true) ≅
              filter([:x, :y] => (x, y) -> x && y === true, df)
        @test subset(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true, view=true) ≅
              filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, df)
        @test subset(gdf, :x, :id => ByRow(<(4))) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test subset(gdf, :x, :id => ByRow(<(4)), view=true) ≅
              filter([:x, :id] => (x, id) -> x && id < 4, df)
        @test_throws ArgumentError subset(gdf)
        @test isempty(subset(gdf, :x, :x => ByRow(!)))
        @test_throws ArgumentError subset(gdf, :x => x -> false, :x => x -> missing)
        @test_throws ArgumentError subset(gdf, :x => x -> true, :x => x -> missing)
        @test_throws ArgumentError subset(gdf, :x => x -> true, :x => x -> 2)
    end

    df = copy(refdf)
    @test subset!(df, :x) === df
    @test subset!(df, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    @test_throws ArgumentError subset!(df, :y)
    @test df ≅ refdf
    df = copy(refdf)
    @test subset!(df, :y, skipmissing=true) === df
    @test subset!(df, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :y, skipmissing=true) === df
    @test subset!(df, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    @test subset!(df, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    @test_throws ArgumentError subset!(df)
    df = copy(refdf)
    @test isempty(subset!(df, :x, :x => ByRow(!)))
    @test isempty(df)

    df = copy(refdf)
    @test_throws ArgumentError subset!(df, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(df, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(df, :x => x -> true, :x => x -> 2)

    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x) === df

    df = copy(refdf)
    gdf = groupby(df, :z)
    gdf2 = subset!(gdf, :x, ungroup=false)
    @test gdf2 isa GroupedDataFrame{DataFrame}
    @test parent(gdf2) === df
    @test gdf2 ≅ groupby(df, :z) ≅ groupby(filter(:x => identity, refdf), :z)

    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf, :y)
    @test df ≅ refdf
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :y, skipmissing=true) === df
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, skipmissing=true) === df
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test subset!(gdf, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test isempty(subset!(gdf, :x, :x => ByRow(!)))
    @test isempty(df)
    df = copy(refdf)
    gdf = groupby(df, :z)
    @test_throws ArgumentError subset!(gdf, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> 2)

    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x) ≅ df ≅ filter(:x => identity, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf, :y)
    @test df ≅ refdf
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :y, skipmissing=true) ≅ df ≅ filter(:y => x -> x === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :y, skipmissing=true) ≅ df ≅
          filter([:x, :y] => (x, y) -> x && y === true, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :y, :id => ByRow(<(4)), skipmissing=true) ≅ df ≅
          filter([:x, :y, :id] => (x, y, id) -> x && y === true && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test subset!(gdf, :x, :id => ByRow(<(4))) ≅ df ≅
          filter([:x, :id] => (x, id) -> x && id < 4, refdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf)
    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test isempty(subset!(gdf, :x, :x => ByRow(!)))
    @test isempty(df)

    df = copy(refdf)
    gdf = groupby(df, :z)[[3, 2, 1]]
    @test_throws ArgumentError subset!(gdf, :x => x -> false, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> missing)
    @test_throws ArgumentError subset!(gdf, :x => x -> true, :x => x -> 2)

    @test_throws ArgumentError subset!(view(refdf, :, :), :x)
    @test_throws ArgumentError subset!(groupby(view(refdf, :, :), :z), :x)

    df = DataFrame(g=[2, 2, 1, 1, 1, 1, 3, 3, 3], x = 1:9)
    @test subset(df, :x => x -> x .< mean(x)) == DataFrame(g=[2, 2, 1, 1], x = 1:4)
    @test subset(groupby(df, :g), :x => x -> x .< mean(x)) ==
          DataFrame(g=[2, 1, 1, 3], x=[1, 3, 4, 7])

    @test_throws ArgumentError subset(df, :x => x -> missing)
    @test isempty(subset(df, :x => ByRow(x -> missing), skipmissing=true))
    @test_throws ArgumentError isempty(subset(df, :x => x -> missing, skipmissing=true))
    @test isempty(subset(df, :x => ByRow(x -> false)))
    @test_throws ArgumentError isempty(subset(df, :x => x -> false))
    @test subset(df, :x => ByRow(x -> true)) ≅ df
    @test_throws ArgumentError subset(df, :x => x -> true) ≅ df
    @test_throws ArgumentError subset(df, :x => x -> (a=x,))
    @test_throws ArgumentError subset(df, :x => (x -> (a=x,)) => AsTable)

    @test_throws ArgumentError subset(DataFrame(x=false, y=missing), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=missing, y=false), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=missing, y=false), :x)
    @test_throws ArgumentError subset(DataFrame(x=false, y=missing), :y)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :x, :y)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :y, :x)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :y)

    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :x, :y, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :x, :y, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1, y=false), :y, :x, skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=false, y=1), :y, skipmissing=true)

    @test_throws ArgumentError DataFrames._and()
    @test_throws ArgumentError DataFrames._and_missing()
end

@testset "subsetting requires passing vector" begin
    @test_throws ArgumentError subset(DataFrame(x=[]), :x => x -> 1)
    @test_throws AssertionError subset(DataFrame(x=[]), :x => ByRow(x -> 1))
    @test_throws ArgumentError subset(DataFrame(x=1:3), [] => () -> true)
    @test subset(DataFrame(x=1:3), [] => ByRow(() -> true)) == DataFrame(x=1:3)
    @test_throws ArgumentError subset(DataFrame(x = [0, 1]), :x => ==(0))
    @test_throws ArgumentError subset(DataFrame(x = [0, missing]), :x => ismissing)
    @test_throws ArgumentError subset(groupby(DataFrame(id = [0, 0, 1, 1],
                                                        x = [-1, 1, 3, 4]), :id),
                                      :x => (x -> sum(x) > 0))
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> fill(missing, length(x)))
    @test isempty(subset(DataFrame(x=1:3), :x => x -> fill(missing, length(x)),
                         skipmissing=true))
    @test subset(DataFrame(x=1:3), :x => x -> fill(true, length(x))) == DataFrame(x=1:3)
    @test subset(DataFrame(x=1:3), :x => x -> trues(length(x))) == DataFrame(x=1:3)
    @test isempty(subset(DataFrame(x=1:3), :x => x -> fill(false, length(x))))
    @test isempty(subset(DataFrame(x=1:3), :x => x -> falses(length(x))))
    @test_throws AssertionError subset(DataFrame(), [] => () -> Union{}[])
    @test_throws AssertionError subset(DataFrame(), [] => () -> Union{}[], skipmissing=true)
    @test subset(DataFrame(x=1:3:15, y=1:5), [:x, :y] => (x, y) -> iseven.(x) .& iseven.(y)) ==
          DataFrame(x=[4, 10], y=[2, 4])
    @test subset(DataFrame(x=1:3:15, y=1:5),
                 AsTable([:x, :y]) => v -> iseven.(v.x) .& iseven.(v.y)) ==
          DataFrame(x=[4, 10], y=[2, 4])

    @test subset(DataFrame(x=1:3), :x => x -> Any[true, false, true]) ==
          DataFrame(x=[1, 3])
    @test subset(DataFrame(x=1:3), :x => x -> view(Any[true, false, true], :)) ==
          DataFrame(x=[1, 3])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, missing])
    @test subset(DataFrame(x=1:3), :x => x -> Any[true, false, missing], skipmissing=true) ==
          DataFrame(x=[1])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, 1])
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> Any[true, false, 1],
                                      skipmissing=true)
    @test_throws ArgumentError subset(DataFrame(x=1:3), :x => x -> (true for i in 1:3))
end

end


text:
using ShaderToy
shadertoy("seascape.frag")


text:
# Different ways to simulate molecules

export
    VelocityVerlet,
    simulate!,
    Verlet,
    StormerVerlet,
    Langevin

# Forces are often expressed per mol but this dimension needs removing for use in the integrator
function remove_molar(x)
    fx = first(x)
    if dimension(fx) == u"𝐋 * 𝐍^-1 * 𝐓^-2"
        T = typeof(ustrip(fx))
        return x / T(Unitful.Na)
    else
        return x
    end
end

"""
    VelocityVerlet(; <keyword arguments>)

The velocity Verlet integrator.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct VelocityVerlet{T, C}
    dt::T
    coupling::C
    remove_CM_motion::Bool
end

function VelocityVerlet(; dt, coupling=NoCoupling(), remove_CM_motion=true)
    return VelocityVerlet(dt, coupling, remove_CM_motion)
end

"""
    simulate!(system, simulator, n_steps; parallel=true)

Run a simulation on a system according to the rules of the given simulator.
Custom simulators should implement this function.
"""
function simulate!(sys,
                    sim::VelocityVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    accels_t = accelerations(sys, neighbors; parallel=parallel)
    accels_t_dt = zero(accels_t)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        run_loggers!(sys, neighbors, step_n)

        sys.coords += sys.velocities .* sim.dt .+ (remove_molar.(accels_t) .* sim.dt ^ 2) ./ 2
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))
        accels_t_dt = accelerations(sys, neighbors; parallel=parallel)
        sys.velocities += remove_molar.(accels_t .+ accels_t_dt) .* sim.dt / 2

        sim.remove_CM_motion && remove_CM_motion!(sys)
        apply_coupling!(sys, sim, sim.coupling)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
            accels_t = accels_t_dt
        end
    end
    return sys
end

"""
    Verlet(; <keyword arguments>)

The leapfrog Verlet integrator.
This is a leapfrog integrator, so the velocities are offset by half a time step
behind the positions.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct Verlet{T, C}
    dt::T
    coupling::C
    remove_CM_motion::Bool
end

function Verlet(; dt, coupling=NoCoupling(), remove_CM_motion=true)
    return Verlet(dt, coupling, remove_CM_motion)
end

function simulate!(sys,
                    sim::Verlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        run_loggers!(sys, neighbors, step_n)

        accels_t = accelerations(sys, neighbors; parallel=parallel)
        sys.velocities += remove_molar.(accels_t) .* sim.dt

        sys.coords += sys.velocities .* sim.dt
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        sim.remove_CM_motion && remove_CM_motion!(sys)
        apply_coupling!(sys, sim, sim.coupling)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
        end
    end
    return sys
end

"""
    StormerVerlet(; <keyword arguments>)

The Störmer-Verlet integrator.
In this case the `velocities` given to the simulator act as the previous step
coordinates for the first step.
Does not currently work with units or thermostats.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
"""
struct StormerVerlet{T, C}
    dt::T
    coupling::C
end

StormerVerlet(; dt, coupling=NoCoupling()) = StormerVerlet(dt, coupling)

function simulate!(sys,
                    sim::StormerVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)

    for step_n in 1:n_steps
        run_loggers!(sys, neighbors, step_n)

        accels_t = accelerations(sys, neighbors; parallel=parallel)

        # Update coordinates
        coords_copy = sys.coords
        for i in 1:length(sys)
            sys.coords[i] = sys.coords[i] + vector(sys.velocities[i], sys.coords[i], sys.box_size) + remove_molar(accels_t[i]) * sim.dt ^ 2
            sys.coords[i] = wrap_coords.(sys.coords[i], sys.box_size)
        end
        sys.velocities = coords_copy

        apply_coupling!(sys, sim, sim.coupling)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n; parallel=parallel)
        end
    end
    return sys
end

"""
    Langevin(; <keyword arguments>)

The Langevin integrator, based on the Langevin Middle Integrator in OpenMM.
This is a leapfrog integrator, so the velocities are offset by half a time step
behind the positions.

# Arguments
- `dt::T`: the time step of the simulation.
- `temperature::K`: the temperature of the simulation.
- `friction::F`: the friction coefficient of the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct Langevin{S, K, F, T}
    dt::S
    temperature::K
    friction::F
    remove_CM_motion::Bool
    vel_scale::T
    noise_scale::T
end

function Langevin(; dt, temperature, friction, remove_CM_motion=true)
    vel_scale = exp(-dt * friction)
    noise_scale = sqrt(1 - vel_scale^2)
    return Langevin(dt, temperature, friction, remove_CM_motion, vel_scale, noise_scale)
end

function simulate!(sys,
                    sim::Langevin,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        run_loggers!(sys, neighbors, step_n)

        accels_t = accelerations(sys, neighbors; parallel=parallel)
        sys.velocities += remove_molar.(accels_t) .* sim.dt

        sys.coords += sys.velocities .* sim.dt / 2
        if isa(sys.coords, CuArray)
            noise = cu(velocity.(Array(mass.(sys.atoms)), (sim.temperature,);
                                    dims=n_dimensions(sys)))
        else
            noise = velocity.(mass.(sys.atoms), (sim.temperature,); dims=n_dimensions(sys))
        end
        sys.velocities = sys.velocities .* sim.vel_scale .+ noise .* sim.noise_scale

        sys.coords += sys.velocities .* sim.dt / 2
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        sim.remove_CM_motion && remove_CM_motion!(sys)
        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
        end
    end
    return sys
end


text:
const ADD = 0
const MUL = 1

function evaluate(line::String, i::Int, linelength::Int)::Tuple{Int,Int}
    res::Vector{Int} = [0]
    op::Int = 0
    while i <= linelength
        if line[i] == ' '
            i += 1
        elseif line[i] == '*'
            op = MUL
            i += 2
        elseif line[i] == '+'
            op = ADD
            i += 2
        elseif line[i] == '('
            (tmp, i) = evaluate(line, i + 1, linelength)
            if op == ADD
                res[end] += tmp
            elseif op == MUL
                push!(res, tmp)
            end
        elseif line[i] == ')'
            return (prod(res), i + 1)
        else
            tmp = parse(Int, line[i])
            if op == ADD
                res[end] += tmp
            elseif op == MUL
                push!(res, tmp)
            end
            i += 1
        end
    end
    return (prod(res), i)
end

function run(s)
    counter::Int = 0
    for line in eachline(IOBuffer(s))
        (res, _) = evaluate(line, 1, lastindex(line))
        counter += res
    end
    return counter
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
using Unitful,LaTeXDatax
using Test

@testset "LaTeXDatax.jl" begin
    io = IOBuffer();

    # Basic data printing
    LaTeXDatax.printdata(io,"String")
    @test String(take!(io)) == "String"

    LaTeXDatax.printdata(io,1.25)
    @test String(take!(io)) == "\$1.25\$"

    LaTeXDatax.printdata(io,3.141592; fmt="%.2f", unitformat=:siunitx)
    @test String(take!(io)) == "\\num{3.14}"

    # keyval printing
    LaTeXDatax.printkeyval(io, :a, 612.2u"nm")
    @test String(take!(io)) == "\\pgfkeyssetvalue{/datax/a}{\$612.2\\;\\mathrm{nm}\$}\n"

    # complete macro
    a = 2;
    b = 3.2u"m";
    @datax a b c=3*a d=27 unitformat:=:siunitx io:=io
    @test String(take!(io)) == """
    \\pgfkeyssetvalue{/datax/a}{\\num{2}}
    \\pgfkeyssetvalue{/datax/b}{\\SI{3.2}{\\meter}}
    \\pgfkeyssetvalue{/datax/c}{\\num{6}}
    \\pgfkeyssetvalue{/datax/d}{\\num{27}}
    """
end


text:
"""

## spabac!

Function spabac! performs Cholesky forward and back-substitution on
a symmetric skyline global matrix. The loads vector is updated.

###Arguments
```julia
spabac!(kv, loads, kdiag)
```

### Arguments
```julia
* kv::Vector{Float64}       : Skyline vector of global stiffness matrix
* loads::Vector{Float64}    : Load vector (Updated)
* kdiag::Vector{Int}      : Diagonal elemnt index vector
```
"""
function spabac!(kv::Vector{Float64}, loads::Vector{Float64}, kdiag::Vector{Int})
  local x::Float64
  n = size(kdiag, 1)
  loads[1]=loads[1]/kv[1]
  for i in 2:n
    ki = kdiag[i]-i
    l = kdiag[i-1] - ki + 1
    x = loads[i]
    if l !== i
      m = i - 1
      for j in l:m
        x -= kv[ki+j]*loads[j]
      end
    end
    loads[i] = x/kv[ki+i]
  end
  for it in 2:n
    i = n + 2 - it
    ki = kdiag[i] - i
    x = loads[i]/kv[ki+i]
    loads[i] = x
    l = kdiag[i-1] - ki + 1
    if l !== i
      m = i - 1
      for k in l:m
        loads[k] -= x*kv[ki+k]
      end
    end
  end
  loads[1] /= kv[1]
  loads
end


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.LinearGaussian
using Test

setRNGs(0)

model, theta, ys, ko = LinearGaussian.defaultLGModel(5)

numParticles = 2^16

smcio = SMCIO{model.particle, model.pScratch}(numParticles, model.maxn, 1,
  false)
smc!(model, smcio)

@test smcio.logZhats ≈ ko.logZhats atol=0.1

nthreads = Threads.nthreads()
if nthreads > 1
  smcio = SMCIO{model.particle, model.pScratch}(numParticles, model.maxn,
    nthreads, false)
  smc!(model, smcio)

  @test smcio.logZhats ≈ ko.logZhats atol=0.1
end

lM = LinearGaussian.makelM(theta)
smcio = SMCIO{model.particle, model.pScratch}(numParticles, model.maxn,
  nthreads, true)
smc!(model, smcio)
path = Vector{model.particle}(undef, smcio.n)
for i in 1:smcio.n
  path[i] = model.particle()
end
SequentialMonteCarlo.pickParticleBS!(path, smcio, lM)


text:
versionof(pkg::Module) = Pkg.dependencies()[Base.PkgId(pkg).uuid].version

"""
Return the type's fields as a tuple
"""
@generated fieldvalues(x) = Expr(:tuple, (:(x.$f) for f=fieldnames(x))...)
@generated fields(x) = Expr(:tuple, (:($f=x.$f) for f=fieldnames(x))...)

"""
Rewrites `@! x = f(args...)` to `x = f!(x,args...)`

Special cases for `*` and `\\` forward to `mul!` and `ldiv!`, respectively.
"""
macro !(ex)
    if @capture(ex, x_ = f_(args__; kwargs_...))
        esc(:($(Symbol(string(f,"!")))($x,$(args...); $kwargs...)))
    elseif @capture(ex, x_ = f_(args__))
        if f == :*
            f = :mul
        elseif f==:\
            f = :ldiv
        end
        esc(:($x = $(Symbol(string(f,"!")))($x,$(args...))::typeof($x))) # ::typeof part helps inference sometimes
    else
        error("Usage: @! x = f(...)")
    end
end

nan2zero(x::T) where {T} = isfinite(x) ? x : zero(T)
@adjoint nan2zero(x::T) where {T} = nan2zero(x), Δ -> (isfinite(x) ? Δ : zero(T),)

""" Return a tuple with the expression repeated n times """
macro repeated(ex,n)
    :(tuple($(repeated(esc(ex),n)...)))
end

"""
Pack some variables in a dictionary

```julia
> x = 3
> y = 4
> @dict x y z=>5
Dict(:x=>3,:y=>4,:z=>5)
```
"""
macro dict(exs...)
    kv(ex::Symbol) = :($(QuoteNode(ex))=>$(esc(ex)))
    kv(ex) = isexpr(ex,:call) && ex.args[1]==:(=>) ? :($(QuoteNode(ex.args[2]))=>$(esc(ex.args[3]))) : error()
    :(Dict($((kv(ex) for ex=exs)...)))
end

macro namedtuple(exs...)
    Base.depwarn("@namedtuple(x,y) is deprecated and will be removed soon, just use the built-in Julia (;x,y) now.", nothing)
    if length(exs)==1 && isexpr(exs[1],:tuple)
        exs = exs[1].args
    end
    kv(ex::Symbol) = :($(esc(ex))=$(esc(ex)))
    kv(ex) = isexpr(ex,:(=)) ? :($(esc(ex.args[1]))=$(esc(ex.args[2]))) : error()
    Expr(:tuple, (kv(ex) for ex=exs)...)
end

# these allow pinv and sqrt of SMatrices of Diagonals to work correctly, which
# we use for the T-E block of the covariance. hopefully some of this can be cut
# down on in the futue with some PRs into StaticArrays.
permutedims(A::SMatrix{2,2}) = @SMatrix[A[1] A[3]; A[2] A[4]]
function sqrt(A::SMatrix{2,2,<:Diagonal})
    a,b,c,d = A
    s = sqrt(a*d-b*c)
    t = pinv(sqrt(a+d+2s))
    @SMatrix[t*(a+s) t*b; t*c t*(d+s)]
end
function pinv(A::SMatrix{2,2,<:Diagonal})
    a,b,c,d = A
    idet = pinv(a*d-b*c)
    @SMatrix[d*idet -(b*idet); -(c*idet) a*idet]
end

# some usefule tuple manipulation functions:

# see: https://discourse.julialang.org/t/efficient-tuple-concatenation/5398/10
# and https://github.com/JuliaLang/julia/issues/27988
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

# see https://discourse.julialang.org/t/any-way-to-make-this-one-liner-type-stable/10636/2
using Base: tuple_type_cons, tuple_type_head, tuple_type_tail, first, tail
map_tupleargs(f,::Type{T}) where {T<:Tuple} =
    (f(tuple_type_head(T)), map_tupleargs(f,tuple_type_tail(T))...)
map_tupleargs(f,::Type{T},::Type{S}) where {T<:Tuple,S<:Tuple} =
    (f(tuple_type_head(T),tuple_type_head(S)), map_tupleargs(f,tuple_type_tail(T),tuple_type_tail(S))...)
map_tupleargs(f,::Type{T},s::Tuple) where {T<:Tuple} =
    (f(tuple_type_head(T),first(s)), map_tupleargs(f,tuple_type_tail(T),tail(s))...)
map_tupleargs(f,::Type{<:Tuple{}}...) = ()
map_tupleargs(f,::Type{<:Tuple{}},::Tuple) = ()

# returns the base parametric type with all type parameters stripped out
basetype(::Type{T}) where {T} = T.name.wrapper
@generated function basetype(t::UnionAll)
    unwrap_expr(s::UnionAll, t=:t) = unwrap_expr(s.body, :($t.body))
    unwrap_expr(::DataType, t) = t
    :($(unwrap_expr(t.parameters[1])).name.wrapper)
end

function ensuresame(args...)
    @assert all(args .== Ref(args[1]))
    args[1]
end

tuple_type_len(::Type{<:NTuple{N,Any}}) where {N} = N

ensure1d(x::Union{Tuple,AbstractArray}) = x
ensure1d(x) = (x,)

# https://discourse.julialang.org/t/dispatching-on-the-result-of-unwrap-unionall-seems-weird/25677
# has some background related to this function. we can simplify this in 1.6
function typealias(t)
    if isconcretetype(t)
        ta = typealias_def(t)
        if !isnothing(ta)
            return ta
        end
    end
    sprint(io -> invoke(Base.show_datatype, Tuple{IO,DataType}, io, t))
end
typealias_def(t) = nothing

@doc doc"""
```
@subst sum(x*$(y+1) for x=1:2)
```

becomes

```
let tmp=(y+1)
    sum(x*tmp for x=1:2)
end
```

to aid in writing clear/succinct code that doesn't recompute things
unnecessarily.
"""
macro subst(ex)

    subs = []
    ex = postwalk(ex) do x
        if isexpr(x, Symbol(raw"$"))
            var = gensym()
            push!(subs, :($(esc(var))=$(esc(x.args[1]))))
            var
        else
            x
        end
    end

    quote
        let $(subs...)
            $(esc(ex))
        end
    end

end

"""
    @ondemand(Package.function)(args...; kwargs...)
    @ondemand(Package.Submodule.function)(args...; kwargs...)

Just like calling `Package.function` or `Package.Submodule.function`, but
`Package` will be loaded on-demand if it is not already loaded. The call is no
longer inferrable.
"""
macro ondemand(ex)
    get_root_package(x) = @capture(x, a_.b_) ? get_root_package(a) : x
    quote
        @eval import $(get_root_package(ex))
        (args...; kwargs...) -> Base.invokelatest($(esc(ex)), args...; kwargs...)
    end
end

get_kwarg_names(func::Function) = Vector{Symbol}(Base.kwarg_decl(first(methods(func))))

# https://discourse.julialang.org/t/is-there-a-way-to-modify-captured-variables-in-a-closure/31213/16
@static if versionof(Adapt) < v"3.1.0"
    @generated function adapt_structure(to, f::F) where {F<:Function}
        if fieldcount(F) == 0
            :f
        else
            quote
                captured_vars = $(Expr(:tuple, (:(adapt(to, f.$x)) for x=fieldnames(F))...))
                $(Expr(:new, :($(F.name.wrapper){map(typeof,captured_vars)...}), (:(captured_vars[$i]) for i=1:fieldcount(F))...))
            end
        end
    end
end

adapt_structure(to, d::Dict) = Dict(k => adapt(to, v) for (k,v) in d)

@doc doc"""

    cpu(x)

Recursively move an object to CPU memory. See also [`gpu`](@ref).
"""
cpu(x) = adapt_structure(Array, x)

@doc doc"""

    @cpu! x y

Equivalent to `x = cpu(x)`, `y = cpu(y)`, etc... for any number of
listed variables. See [`cpu`](@ref).
"""
macro cpu!(vars...)
    :(begin; $((:($(esc(var)) = cpu($(esc(var)))) for var in vars)...); nothing; end)
end

@doc doc"""

    gpu(x)

Recursively move an object to GPU memory. Note that, unlike `cu(x)`,
this does not change the `eltype` of any underlying arrays. See also
[`cpu`](@ref).
"""
function gpu end # defined in gpu.jl only when CUDA.jl is loaded

function corrify(H)
    H = copy(H)
    σ = sqrt.(abs.(diag(H)))
    for i=1:checksquare(H)
        H[i,:] ./= σ
        H[:,i] ./= σ
    end
    H
end

struct LazyPyImport
    pkg
end
getproperty(p::LazyPyImport, s::Symbol) = getproperty(@ondemand(PyCall.pyimport)(getfield(p,:pkg)), s)

@doc doc"""

    lazy_pyimport(s)

Like `pyimport(s)`, but doesn't actually load anything (not even
PyCall) until a property of the returned module is accessed, allowing
this to go in `__init__` and still delay loading PyCall, as well as
preventing a Julia module load error if a Python module failed to load.
"""
lazy_pyimport(s) = LazyPyImport(s)

@doc doc"""
    @ismain()

Return true if the current file is being run as a script.
"""
macro ismain()
    (__source__ != nothing) && (String(__source__.file) == abspath(PROGRAM_FILE))
end

firsthalf(x) = x[1:end÷2]
lasthalf(x) = x[end÷2:end]

USE_SUM_KBN = true
use_sum_kbn!(flag) = USE_SUM_KBN=flag

# type-stable combination of summing and dropping dims, which uses
# either sum or sum_kbn (to reduce roundoff error), depending on
# CMBLensing.USE_SUM_KBN constant
function sum_dropdims(A::AbstractArray{T,N}; dims=:) where {T,N}
    if (dims == (:)) || (N == length(dims))
        if USE_SUM_KBN
            sum_kbn(cpu(A))
        else
            sum(A)
        end :: T
    else
        if USE_SUM_KBN
            dropdims(mapslices(sum_kbn, cpu(A), dims=dims), dims=dims)
        else
            dropdims(mapslices(sum,         A,  dims=dims), dims=dims)
        end :: Array{T,N-length(dims)}
    end
end
@adjoint sum_dropdims(A) = sum_dropdims(A), Δ -> (fill!(similar(A),Δ),)

# for mixed eltype, which Loess stupidly does not support
Loess.loess(x::AbstractVector, y::AbstractVector; kwargs...) =
    loess(collect.(zip(promote.(x,y)...))...; kwargs...)

expnorm(x) = exp.(x .- maximum(x))

# MacroTool's is broken https://github.com/FluxML/MacroTools.jl/issues/154
_isdef(ex) = @capture(ex, function f_(arg__) body_ end)

"""

    @⌛ [label] code ...
    @⌛ [label] function_definition() = ....

Label a section of code to be timed. If a label string is not
provided, the first form uses the code itselfs as a label, the second
uses the function name, and its the body of the function which is
timed.

To run the timer and print output, returning the result of the
calculation, use

    @show⌛ run_code()

Timing uses `TimerOutputs.get_defaulttimer()`.
"""
macro ⌛(args...)
    if length(args)==1
        label, ex = nothing, args[1]
    else
        label, ex = esc(args[1]), args[2]
    end
    source_str = last(splitpath(string(__source__.file)))*":"*string(__source__.line)
    if _isdef(ex)
        sdef = splitdef(ex)
        if isnothing(label)
            label = "$(string(sdef[:name]))(…)  ($source_str)"
        end
        sdef[:body] = quote
            CMBLensing.@timeit $label $(sdef[:body])
        end
        esc(combinedef(sdef))
    else
        if isnothing(label)
            label = "$(Base._truncate_at_width_or_chars(string(prewalk(rmlines,ex)),26))  ($source_str)"
        end
        :(@timeit $label $(esc(ex)))
    end
end

"""
See [`@⌛`](@ref)
"""
macro show⌛(ex)
    quote
        reset_timer!(get_defaulttimer())
        result = $(esc(ex))
        show(get_defaulttimer())
        result
    end
end

# used in a couple of places to create a Base.promote_rule-like system
# where you can specify a set of rules for promotion via dispatch but
# don't need to write a method for both orders
select_known_rule(rule, x, y) = select_known_rule(rule, x, y, rule(x,y), rule(y,x))
select_known_rule(rule, x, y, R₁::Any,       R₂::Unknown) = R₁
select_known_rule(rule, x, y, R₁::Unknown,   R₂::Any)     = R₂
select_known_rule(rule, x, y, R₁::Any,       R₂::Any)     = (R₁ == R₂) ? R₁ : error("Conflicting rules.")
select_known_rule(rule, x, y, R₁::Unknown,   R₂::Unknown) = unknown_rule_error(rule, x, y)

"""
    @auto_adjoint foo(args...; kwargs...) = body

is equivalent to

    _foo(args...; kwargs...) = body
    foo(args...; kwargs...) = _foo(args...; kwargs...)
    @adjoint foo(args...; kwargs...) = Zygote.pullback(_foo, args...; kwargs...)

That is, it defines the function as well as a Zygote adjoint which
takes a gradient explicitly through the body of the function, rather
than relying on rules which may be defined for `foo`. Mainly useful in
the case that `foo` is a common function with existing rules, but
which you do *not* want to be used.
"""
macro auto_adjoint(funcdef)
    sdef = splitdef(funcdef)
    name = sdef[:name]
    sdef[:name] = symname = gensym(string(name))
    defs = []
    push!(defs, combinedef(sdef))
    sdef[:name] = name
    sdef[:body] = :($symname($(sdef[:args]...); $(sdef[:kwargs]...)))
    push!(defs, :(Core.@__doc__ $(combinedef(sdef))))
    sdef[:body] = :($Zygote.pullback($symname, $(sdef[:args]...); $(sdef[:kwargs]...)))
    push!(defs, :($Zygote.@adjoint $(combinedef(sdef))))
    esc(Expr(:block, defs...))
end

string_trunc(x) = Base._truncate_at_width_or_chars(string(x), displaysize(stdout)[2]-14)

import NamedTupleTools
NamedTupleTools.select(d::Dict, keys) = (;(k=>d[k] for k in keys)...)

@init @require ComponentArrays="b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" begin
    using ComponentArrays
    # a Zygote-compatible conversion of ComponentVector to a NamedTuple
    Base.convert(::Type{NamedTuple}, x::ComponentVector) = NamedTuple{keys(x)}([x[k] for k in keys(x)])
    @adjoint Base.convert(::Type{NamedTuple}, x::ComponentVector) = convert(NamedTuple, x), Δ -> (nothing, ComponentArray(Δ))
end


text:
function get_thermal_standard_uc_template()
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    return template
end

function get_thermal_dispatch_template_network(network=CopperPlatePowerModel)
    template = ProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, MonitoredLine, StaticBranchBounds)
    set_device_model!(template, Line, StaticBranch)
    set_device_model!(template, Transformer2W, StaticBranch)
    set_device_model!(template, TapTransformer, StaticBranch)
    set_device_model!(template, HVDCLine, HVDCLossless)
    return template
end

function get_template_basic_uc_simulation()
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, InterruptibleLoad, StaticPowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    return template
end

function get_template_standard_uc_simulation()
    template = get_template_basic_uc_simulation()
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    return template
end

function get_template_nomin_ed_simulation(network=CopperPlatePowerModel)
    template = ProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    return template
end

function get_template_hydro_st_uc(network=CopperPlatePowerModel)
    template = ProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment),
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch),
    set_device_model!(template, PowerLoad, StaticPowerLoad),
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad),
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage),
    return template
end

function get_template_hydro_st_ed(network=CopperPlatePowerModel, duals=[])
    template = ProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    return template
end

function get_template_dispatch_with_network(network=StandardPTDFModel)
    template = ProblemTemplate(network)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalBasicDispatch)
    set_device_model!(template, Line, StaticBranch)
    set_device_model!(template, Transformer2W, StaticBranch)
    set_device_model!(template, TapTransformer, StaticBranch)
    set_device_model!(template, HVDCLine, HVDCLossless)
    return template
end