text:
# Copyright © 2021, UChicago Argonne, LLC
# All Rights Reserved.

# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

function get_vbounds(obs, operators)
    # Initialization
    operators = Set(operators)
    (v_lb, v_ub) = (-1e+09, 1e+09)
    factor1 = 10
    factor2 = 2
    eps = 1e-09

    # if !isempty(intersect(operators, Set("EL")))
    #     factor1 = 5
    #     factor2 = 2
    # end

    ## Compute min/max of indep/dep variables
    # (z_lb, z_ub) = (minimum(obs[:,end]), maximum(obs[:,end]))
    # (x_lb, x_ub) = (minimum(obs[:,1:end-1]), maximum(obs[:,1:end-1]))
    (obs_lb, obs_ub) = (minimum(obs), maximum(obs))

    ## (1) Shrink bounds
    (v_lb, v_ub) = (maximum([v_lb, factor1^(-sign(obs_lb)) * obs_lb]), minimum([v_ub, factor1 * obs_ub]))

    ## (2) Expand bounds if it is too tight
    (v_lb, v_ub) = (minimum([v_lb, factor2^(-sign(obs_lb)) * obs_lb]), maximum([v_ub, factor2 * obs_ub]))

    ## (3) make bounds symmetric and +- epsilon
    magnitude = maximum([abs(v_lb), abs(v_ub)]) + eps
    (v_lb, v_ub) = (-magnitude, magnitude)

    v_lb, v_ub
end

function get_bounds(obs, operators, params)
    bounds = Dict()
    v_lb, v_ub = get_vbounds(obs, operators)
    merge!(bounds, Dict("v_lb"=>v_lb, "v_ub"=>v_ub))
    merge!(bounds, Dict("c_lb"=>get(params, "DBL_CONSTANT_LB", -2.0), "c_ub"=>get(params, "DBL_CONSTANT_UB", 2.0)))

    bounds["cst_abslb"] = 1/(4*pi)
    bounds["mult_child_abslb"] = 0.01
    bounds["div_child_abslb"] = 0.01
    bounds["sqrt_rhs_lb"] = 0.01
    bounds["exp_rhs_lb"] = v_lb > 0 ? max(v_lb, log(v_lb)) : v_lb
    bounds["exp_rhs_ub"] = min(v_ub, log(v_ub))
    bounds["log_rhs_lb"] = max(v_lb, exp(v_lb))
    bounds["log_rhs_ub"] = min(v_ub, exp(v_ub))

    bounds
end

function get_df_vars(nodes, y_indexes, num_obs)
    df_vars = DataFrame(vartype=String[], nodeid=Int[], dataid=Int[], operator=Any[], sol=Float64[])
    for (n,o) in y_indexes
        append!(df_vars, Dict("vartype"=>"y", "nodeid"=>n, "dataid"=>0, "operator"=>o, "sol"=>0.0))
    end

    for n in nodes, i=1:num_obs
        append!(df_vars, Dict("vartype"=>"v", "nodeid"=>n, "dataid"=>i, "operator"=>0, "sol"=>0.0))
    end

    for n in nodes
        append!(df_vars, Dict("vartype"=>"c", "nodeid"=>n, "dataid"=>0, "operator"=>0, "sol"=>0.0))
    end

    for i=1:num_obs
        append!(df_vars, Dict("vartype"=>"e", "nodeid"=>0, "dataid"=>i, "operator"=>0, "sol"=>0.0))
    end

    df_vars.index = 1:nrow(df_vars)
    select!(df_vars, ["index", "vartype", "nodeid", "dataid", "operator", "sol"])
    df_vars
end

function create_decision_variables(model, nodes, operA, y_indexes, num_obs, bounds, params)
    @variable(model, y[n in nodes, o in operA; (n,o) in y_indexes], Bin)
    @variable(model, bounds["v_lb"] <= v[1:num_obs, nodes] <= bounds["v_ub"])
    if get(params, "INT_CONSTANT_ISINTEGER", 0) == 1
        @variable(model, bounds["c_lb"] <= c[nodes] <= bounds["c_ub"], Int)
    else
        @variable(model, bounds["c_lb"] <= c[nodes] <= bounds["c_ub"])
    end
    @variable(model, -1e+09 <= e[1:num_obs] <= 1e+09)
    df_vars = get_df_vars(nodes, y_indexes, num_obs)
    model, y, c, v, e, df_vars
end

function set_objective_function(model, e, num_obs)
    @objective(model, Min, 1 / num_obs * sum(e.^2))
    model
end

function limit_number_of_constant(model, y, y_indexes, numcst)
    @constraint(model, sum(y for (n,o) in y_indexes if o == 'c') <= numcst)
    model
end

function add_e_deining_constraints(model, v, e, obs, num_obs)
    @constraint(model, edef[i in 1:num_obs], e[i] == v[i,1] - obs[i,end])
    model
end

function add_tree_defining_constraints(model, formulation, y, y_indexes, nodes, nleaves, operBU, operB, operU, operL, indvars)
    if occursin("Coz", formulation)
        @constraints(model, begin
            # [Cozad]-(1d) at most one operator is active
            active_at_most[n1 in nodes], sum(y[n,o] for (n,o) in y_indexes if n == n1) <= 1
            # [Cozad]-(1i) at least one independent variable is active
            active_indvar, sum(y[n,o] for (n,o) in y_indexes if o in indvars) >= 1
            # [Cozad]-(1e) if binary/unary is active, then rch is active
            active_rch[n1 in nleaves],
                sum(y[n1,o] for o in operBU) <=
                sum(y[n,o] for (n,o) in y_indexes if n == 2*n1+1)
            # [Cozad]-(1f) if binary is active, then lch is active
            active_lch[n1 in nleaves],
                sum(y[n1,o] for o in operB) <=
                sum(y[n,o] for (n,o) in y_indexes if n == 2*n1)
            # [Cozad]-(1h)
            inactive_rch[n1 in nleaves],
                sum(y[n1,o] for o in operL) <=
                1 - sum(y[n,o] for (n,o) in y_indexes if n == 2*n1+1)
            # [Cozad]-(1g)
            inactive_lch[n1 in nleaves],
                sum(y[n1,o] for o in union(operU, operL)) <=
                1 - sum(y[n,o] for (n,o) in y_indexes if n == 2*n1)
        end)
    elseif occursin("Imp", formulation)
        @constraints(model, begin
            # [Cozad]-(1d) at most one operator is active
            active_at_most[n1 in nodes], sum(y[n,o] for (n,o) in y_indexes if n == n1) <= 1
            # [Cozad]-(1i) at least one independent variable is active
            active_indvar, sum(y[n,o] for (n,o) in y_indexes if o in indvars) >= 1
            # improvement of (1e-1h). Also, improvement of [Newman]'s (18-23)
            active_inactive_rch[n1 in nleaves],
                sum(y[n1,o] for o in operBU) ==
                sum(y[n,o] for (n,o) in y_indexes if n == 2*n1+1)
            active_inactive_lch[n1 in nleaves],
                sum(y[n1,o] for o in operB) ==
                sum(y[n,o] for (n,o) in y_indexes if n == 2*n1)
        end)
    else
        nothing
    end

    model
end

function add_value_defining_constraints(model, formulation, y, c, v, y_indexes, bounds, nodes, nleaves, obs, num_obs, operA, operB, indvars)
    ### Constrs for defining domain bounds
    @info "Constrs for defining domain bounds"

    v_lb = bounds["v_lb"]
    v_ub = bounds["v_ub"]
    c_lb = bounds["c_lb"]
    c_ub = bounds["c_ub"]

    if 'C' in operA && bounds["cst_abslb"] > 0
        @constraint(model, domain_cst[n in nodes], c[n]^2 >= bounds["cst_abslb"]^2 * y[n,'C'])
    end

    if ('*' in operA && bounds["mult_child_abslb"] > 0) || ('D' in operA && bounds["div_child_abslb"] > 0)
        @constraint(model, domain_lch_abslb[i in 1:num_obs, n in nleaves],
        v[i,2*n]^2 >=
            ( '*' in operA ? bounds["mult_child_abslb"]^2 * y[n,'*'] : 0 ) +
            ( 'D' in operA ? bounds["div_child_abslb"]^2 * y[n,'D'] : 0 ) )

        @constraint(model, domain_rch_abslb[i in 1:num_obs, n in nleaves],
            v[i,2*n+1]^2 >=
                ( '*' in operA ? bounds["mult_child_abslb"]^2 * y[n,'*'] : 0 ) +
                ( 'D' in operA ? bounds["div_child_abslb"]^2 * y[n,'D'] : 0 ) )
    end

    if ('R' in operA && bounds["sqrt_rhs_lb"] > 0) || ('E' in operA && bounds["exp_rhs_lb"] > v_lb) ||
        ('L' in operA && bounds["log_rhs_lb"] > v_lb)
        @constraint(model, domain_rhs_lb[i in 1:num_obs, n in nleaves],
            v[i,2*n+1] >= v_lb +
                ( ('R' in operA && bounds["sqrt_rhs_lb"] > 0) ? (bounds["sqrt_rhs_lb"] - v_lb) * y[n,'D'] : 0 ) +
                ( ('E' in operA && bounds["exp_rhs_lb"] > v_lb) ? (bounds["exp_rhs_lb"] - v_lb) * y[n,'E'] : 0 ) +
                ( ('L' in operA && bounds["log_rhs_lb"] > v_lb) ? (bounds["log_rhs_lb"] - v_lb) * y[n,'L'] : 0 ) )
    end

    if ('E' in operA && bounds["exp_rhs_ub"] < v_ub) || ('L' in operA && bounds["log_rhs_ub"] < v_ub)
        @constraint(model, domain_rhs_ub[i in 1:num_obs, n in nleaves],
            v[i,2*n+1] <= v_ub +
                ( ('E' in operA && bounds["exp_rhs_ub"] < v_ub) ? (bounds["exp_rhs_ub"] - v_ub) * y[n,'E'] : 0 ) +
                ( ('L' in operA && bounds["log_rhs_ub"] < v_ub) ? (bounds["log_rhs_ub"] - v_ub) * y[n,'L'] : 0 ) )
    end

    if 'S' in operA
        @constraints(model, begin
            domain_asincos_node_lb[i in 1:num_obs, n in nleaves], v[i,n] >= v_lb*(1 - y[n,'S'])
            domain_asincos_node_ub[i in 1:num_obs, n in nleaves], v[i,n] <= y[n,'S'] + v_ub*(1 - y[n,'S'])
            domain_asincos_rch_lb[i in 1:num_obs, n in nleaves], v[i,2*n+1] >= -y[n,'S'] + v_lb*(1 - y[n,'S'])
            domain_asincos_rch_ub[i in 1:num_obs, n in nleaves], v[i,2*n+1] <= y[n,'S'] + v_ub*(1 - y[n,'S'])
        end)
    end

    ### Constrs for defining v
    if occursin("Coz", formulation)
        @info "Constrs for defining v (Cozad) $(formulation)"

        @constraints(model, begin
            v_inactive_ub[i in 1:num_obs, n in nodes], v[i,n] <= v_ub * sum(y[n1,o] for (n1,o) in y_indexes if n1 == n)
            v_inactive_lb[i in 1:num_obs, n in nodes], v[i,n] >= v_lb * sum(y[n1,o] for (n1,o) in y_indexes if n1 == n)
        end)

        @constraints(model, begin
            v_indvars_ub[i in 1:num_obs, n in nodes, o in indvars], v[i,n] - obs[i,o] <= (v_ub - obs[i,o]) * (1 - y[n,o])
            v_indvars_lb[i in 1:num_obs, n in nodes, o in indvars], v[i,n] - obs[i,o] >= (v_lb - obs[i,o]) * (1 - y[n,o])
        end)

        if 'P' in operA
            @constraints(model, begin
                v_pi_ub[i in 1:num_obs, n in nodes], v[i,n] - pi <= (v_ub - v_lb) * (1 - y[n,'P'])
                v_pi_lb[i in 1:num_obs, n in nodes], v[i,n] - pi >= (v_lb - v_ub) * (1 - y[n,'P'])
            end)
        end

        if 'C' in operA
            @constraints(model, begin
                v_cst_ub[i in 1:num_obs, n in nodes], v[i,n] - c[n] <= (v_ub - c_lb) * (1 - y[n,'C'])
                v_cst_lb[i in 1:num_obs, n in nodes], v[i,n] - c[n] >= (v_lb - c_ub) * (1 - y[n,'C'])
            end)
        end
    elseif occursin("Imp", formulation)
        @info "Constrs for defining v (New) $(formulation)"

        @constraints(model, begin
            v_indvars_ub[i in 1:num_obs, n1 in nodes],
                v[i,n1] <=
                    v_ub * sum(y[n,o] for (n,o) in y_indexes if n==n1 && o in operB) +
                    min(v_ub, sqrt(v_ub)) * ((n1,'R') in y_indexes ? y[n1,'R'] : 0) +
                    1 * ((n1,'S') in y_indexes ? y[n1,'S'] : 0) +
                    min(v_ub, exp(v_ub)) * ((n1,'E') in y_indexes ? y[n1,'E'] : 0) +
                    min(v_ub, log(v_ub)) * ((n1,'L') in y_indexes ? y[n1,'L'] : 0) +
                    c_ub * ((n1,'C') in y_indexes ? y[n1,'C'] : 0) +
                    pi * ((n1,'P') in y_indexes ? y[n1,'P'] : 0) +
                    sum(obs[i,o] * y[n1,o] for o in indvars)
            v_indvars_lb[i in 1:num_obs, n1 in nodes],
                v[i,n1] >=
                    v_lb * sum(y[n,o] for (n,o) in y_indexes if n==n1 && o in operB) +
                    max(v_lb, sqrt(bounds["sqrt_rhs_lb"])) * ((n1,'R') in y_indexes ? y[n1,'R'] : 0) +
                    max(v_lb, exp(bounds["exp_rhs_lb"])) * ((n1,'E') in y_indexes ? y[n1,'E'] : 0) +
                    max(v_lb, log(bounds["log_rhs_lb"])) * ((n1,'L') in y_indexes ? y[n1,'L'] : 0) +
                    c_lb * ((n1,'C') in y_indexes ? y[n1,'C'] : 0) +
                    pi * ((n1,'P') in y_indexes ? y[n1,'P'] : 0) +
                    sum(obs[i,o] * y[n1,o] for o in indvars)
        end)

        if 'C' in operA
            @constraints(model, begin
                v_cst_ub[i in 1:num_obs, n in nodes],
                    v[i,n] - c[n] <= (v_ub - c_lb) * (1 - y[n,'C'])
                v_cst_lb[i in 1:num_obs, n in nodes],
                    v[i,n] - c[n] >= (v_lb - c_ub) * (1 - y[n,'C'])
            end)
        end
    end

    if '+' in operA
        @constraints(model, begin
            v_plus_ub[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] + v[i,2*n+1]) <= (v_ub - 2 * v_lb) * (1 - y[n,'+'])
            v_plus_lb[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] + v[i,2*n+1]) >= (v_lb - 2 * v_ub) * (1 - y[n,'+'])
        end)
    end

    if '-' in operA
        @constraints(model, begin
            v_minus_ub[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] - v[i,2*n+1]) <= (2 * v_ub - v_lb) * (1 - y[n,'-'])
            v_minus_lb[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] - v[i,2*n+1]) >= (2 * v_lb - v_ub) * (1 - y[n,'-'])
        end)
    end

    if '*' in operA
        @constraints(model, begin
            v_mult_ub[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] * v[i,2*n+1]) <=
                (v_ub - min(v_lb^2, v_ub^2, v_lb*v_ub)) * (1 - y[n,'*'])
            v_mult_lb[i in 1:num_obs, n in nleaves],
                v[i,n] - (v[i,2*n] * v[i,2*n+1]) >=
                (v_lb - max(v_lb^2, v_ub^2)) * (1 - y[n,'*'])
        end)
    end

    if 'D' in operA
        @constraints(model, begin
            v_div_ub[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,2*n+1] - v[i,2*n] <=
                (max(v_lb^2, v_ub^2) - v_lb) * (1 - y[n,'D'])
            v_div_lb[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,2*n+1] - v[i,2*n] >=
                (min(v_lb^2, v_ub^2, v_lb*v_ub) - v_ub) * (1 - y[n,'D'])
        end)
    end

    if 'R' in operA
        @constraints(model, begin
            v_sqrt_ub[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,n] - v[i,2*n+1] <=
                (max(v_lb^2, v_ub^2) - v_lb) * (1 - y[n,'R'])
            v_sqrt_lb[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,n] - v[i,2*n+1] >=
                (- v_ub) * (1 - y[n,'R'])
        end)
    end

    if 'S' in operA
        @constraints(model, begin
            v_asincos_ub[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,n] + v[i,2*n+1] * v[i,2*n+1] <= y[n,'S'] + 2 * max(v_lb^2, v_ub^2) * (1 - y[n,'S'])
            v_asincos_lb[i in 1:num_obs, n in nleaves],
                v[i,n] * v[i,n] + v[i,2*n+1] * v[i,2*n+1] >= y[n,'S']
        end)
    end

    if 'E' in operA
        @NLconstraints(model, begin
            v_exp_ub[i in 1:num_obs, n in nleaves],
                v[i,n] - exp(v[i,2*n+1]) <= v_ub * (1 - y[n,'E'])
            v_exp_lb[i in 1:num_obs, n in nleaves],
                v[i,n] - exp(v[i,2*n+1]) >= (v_lb - exp(v_ub)) * (1 - y[n,'E'])
        end)
    end

    if 'L' in operA
        @NLconstraints(model, begin
            v_log_ub[i in 1:num_obs, n in nleaves],
                exp(v[i,n]) - v[i,2*n+1] <= (exp(v_ub) - v_lb) * (1 - y[n,'L'])
            v_log_lb[i in 1:num_obs, n in nleaves],
                exp(v[i,n]) - v[i,2*n+1] >= (-v_ub) * (1 - y[n,'L'])
        end)
    end

    model
end

function add_redundancy_eliminating_constraints(model, formulation, y, c, y_indexes, nodes, nleaves, operA, operU)
    ### Nested Operations, e.g., no exp(log(x)) or log(exp(x))
    if issubset(Set("EL"), operA)
        @constraints(model, begin
            redun_explog[n in nleaves; 4*n+3 in nodes], y[n,'E'] + y[2*n+1,'L'] <= 1
            redun_logexp[n in nleaves; 4*n+3 in nodes], y[n,'L'] + y[2*n+1,'E'] <= 1
        end)
    end

    ### Constant Operations
    if 'C' in operA && !JuMP.is_integer(c[1])
        # no operation of two constants and no
        @constraint(model, redun_cst_oper[n in nleaves], y[2*n,'C'] + y[2*n+1,'C'] <= 1)

        # no x-cst, x/cst, and no unary operation with a constant
        if occursin("Coz", formulation)
            @info "Constrs to remove redundancy (Cozad version) $(formulation)"
            @constraints(model, begin
                redun_cst_unary[n in nleaves], y[2*n+1,'C'] + sum(y[nn,o] for (nn,o) in y_indexes if nn == n && o in operU) <= 1
                redun_cst_minus[n in nleaves], y[2*n+1,'C'] + sum(y[nn,o] for (nn,o) in y_indexes if nn == n && o == '-') <= 1
                redun_cst_div[n in nleaves], y[2*n+1,'C'] + sum(y[nn,o] for (nn,o) in y_indexes if nn == n && o == 'D') <= 1
            end)
        elseif occursin("Imp", formulation)
            @info "Constrs to remove redundancy (New version) $(formulation)"
            oper = intersect(operA, "+*")
            @constraint(model, redun_cst_rhs[n in nleaves], y[2*n+1,'C'] <= sum(y[n,o] for o in intersect(operA, "+*")))
        end
    end

    ### Association Properties
    if occursin("Imp", formulation)
        @info "Constrs to remove redundancy (New version 2) $(formulation)"
        nodes_constr = intersect(get_nodes_complete(nodes), get_nodes_grandparents(nodes))
        if issubset("+-", operA)
            @constraint(model, redun_pm[n in nodes_constr], y[n,'+'] + y[2*n+1,'-'] <= 1)
        end
        if issubset("*D", operA)
            @constraint(model, redun_pd[n in nodes_constr], y[n,'*'] + y[2*n+1,'D'] <= 1)
        end
    end

    model
end

function add_implication_cuts(model)
    model
end

function add_symmetry_breaking_constraints(model, y, v, nodes, nleaves, bounds)
    @info "Constrs to break symmetry"
    nodes_constr = intersect(get_nodes_complete(nodes), nleaves)
    @constraint(model, sym[n in nodes_constr], v[1,2*n] - v[1,2*n+1] >= (bounds["v_lb"] - bounds["v_ub"]) * (1 - y[n,'+'] - y[n,'*']))
    model
end

function get_MINLP_model(nodes, obs, operators, formulation, params)
    ### Declare useful parameters
    leaves		= get_leaves(nodes)
    nleaves		= get_nonleaves(nodes)
    num_obs 	= size(obs, 1)
    num_indvars = size(obs, 2) - 1
    indvars		= OrderedSet(1:num_indvars)
    operA		= union(operators, indvars)
    operB       = intersect(operators, Set("+-*D"))
    operU       = intersect(operators, Set("RSEL"))
    operBU      = union(operB, operU)
    operL       = union(intersect(operators, Set("CP")), indvars)
    y_indexes   = [(n,o) for n in nodes for o in operA if n in nleaves || o in operL]

    ### Create a model
    model = Model()
    bounds = get_bounds(obs, operators, params)
    model, y, c, v, e, df_vars = create_decision_variables(model, nodes, operA, y_indexes, num_obs, bounds, params)
    model = set_objective_function(model, e, num_obs)
    if get(params, "INT_NUM_CONSTANT", Inf) < length(nodes) / 2
        model = limit_number_of_constant(model, y, y_indexes, get(params, "INT_NUM_CONSTANT", Inf))
    end
    model = add_e_deining_constraints(model, v, e, obs, num_obs)
    model = add_tree_defining_constraints(model, formulation, y, y_indexes, nodes, nleaves, operBU, operB, operU, operL, indvars)
    model = add_value_defining_constraints(model, formulation, y, c, v, y_indexes, bounds, nodes, nleaves, obs, num_obs, operA, operB, indvars)
    if occursin("-R", formulation) || occursin("-F", formulation)
        model = add_redundancy_eliminating_constraints(model, formulation, y, c, y_indexes, nodes, nleaves, operA, operU)
    end
    if occursin("-S", formulation) || occursin("-F", formulation)
        model = add_symmetry_breaking_constraints(model, y, v, nodes, nleaves, bounds)
    end
    # if occursin("Imp", formulation)
    #     model = add_implication_cuts(model)
    # end

    model, df_vars, bounds
end

function save_solution(dir, model, df_vars, obs; save_df_sol=true, save_df_vars=false)
    num_results = JuMP.result_count(model)
    df_sol = DataFrame(index=Int[], formula=String[], objval=Float64[], treesize=Int[])
    if num_results > 0
        for i in 1:num_results
            if JuMP.has_values(model; result = i)
                df_vars.sol = JuMP.value.(JuMP.all_variables(model), result=i)
                ysol, csol = get_ysol_csol(df_vars)
                @info "ysol and csol", ysol, csol
                formula = get_formula(ysol, csol)[1]
                objval = compute_err_formula(obs, formula, "mse"; err2Inf=true)
                append!(df_sol, Dict("index"=>i, "formula"=>formula, "objval"=>objval, "treesize"=>length(ysol)))
            end
        end
        df_sol.time = repeat([JuMP.solve_time(model)], nrow(df_sol))
        df_sol.bnbnodes = repeat([JuMP.node_count(model)], nrow(df_sol))
    end

    ### Save the summary of the list of the solutions
    sort!(df_sol, ["objval"])
    if save_df_sol
        CSV.write(dir * "df_sol.csv", df_sol)
    end

    ### Save the best solution
    if num_results > 0
        df_vars.sol = JuMP.value.(JuMP.all_variables(model), result=df_sol[1,"index"])
        ysol, csol = get_ysol_csol(df_vars)
        if save_df_vars
            CSV.write(dir * "df_vars.csv", df_vars)
        end

        return df_sol[1,"objval"], ysol, csol, df_vars
    else
        return Inf, Dict(), Dict(), df_vars
    end
end

function MINLP(dir; datafile="data.txt", paramfile="params.txt")
    params = read_parameters(dir * paramfile)
    operators = params["OPERATORS"]
    maxdepth = params["INT_MAXDEPTH"]
    timelimit = params["DBL_TIMELIMIT"]
    formulation = params["FORMULATION"]
    loglevel = params["LOGLEVEL"]
    io = open_logger(dir * "log.log", loglevel)

    ### Parameters
    nodes = OrderedSet(1:(2^(maxdepth+1)-1))
    obs = DelimitedFiles.readdlm(dir * datafile)

    ### Build a MINLP model
    model, df_vars, bounds = get_MINLP_model(nodes, obs, operators, formulation, params)

    ### Set the MINLP solver
    set_optimizer(model, SCIP.Optimizer)
    set_optimizer_attribute(model, "limits/time", timelimit)

    ### Solve
    optimize!(model)

    ### Save the result
    save_solution(dir, model, df_vars, obs)

    close(io)
end

function change_c_variables(model, df_vars, params)

end

max_stars_repo_name: jongeunkim/STreCH

max_issues_repo_name: jongeunkim/STreCH

max_forks_repo_name: jongeunkim/STreCH


text:
using Test, QuantumPEPS, Yao
using LinearAlgebra, Random

@testset "pswap gate" begin
    pb = QuantumPEPS.pswap(6, 2, 4)
    reg = rand_state(6)
    @test copy(reg) |> pb ≈ invoke(apply!, Tuple{ArrayReg, PutBlock}, copy(reg), pb)
    @test copy(reg) |> pb ≈ reg

    dispatch!(pb, π)
    @test copy(reg) |> pb ≈ -im*(copy(reg) |> swap(6, 2, 4))
    @test copy(reg) |> pb |> isnormalized

    dispatch!(pb, :random)
    @test copy(reg) |> pb ≈ invoke(apply!, Tuple{ArrayReg, PutBlock}, copy(reg), pb)
end

@testset "bag" begin
    bag = Bag(X)
    reg = zero_state(1)
    @test reg |> bag ≈ ArrayReg([0im, 1.0])
    @test mat(bag) == mat(X)
    @test isreflexive(bag)
    @test ishermitian(bag)
    @test isunitary(bag)

    setcontent!(bag, Z)
    print(bag)
    @test reg |> bag ≈ ArrayReg([0im, -1.0])
    @test mat(bag) == mat(Z)
end

@testset "basis rotor" begin
    mblock = Measure(1)

    reg = repeat(ArrayReg([0im,1]/sqrt(2)), 10)
    reg |> basis_rotor(Z) |> mblock
    @test mblock.results == fill(1,10)

    reg = repeat(ArrayReg([1+0im,1]/sqrt(2)), 10)
    reg |> basis_rotor(X) |> mblock
    @test mblock.results == fill(0,10)

    reg = repeat(ArrayReg([im,1.0]/sqrt(2)), 10)
    reg |> basis_rotor(Y) |> mblock
    @test mblock.results == fill(1,10)
end

@testset "QPEPSConfig" begin
    c = QPEPSConfig(4, 2, 4, 1)
    @test nqubits(c) == 12
    @test mbit(c, 3) == 3
    @test vbit(c, 4) == 11:12
    @test vbit(c, 3, 2) == 10
    @test nspins(c) == 4*4+8
    @test nparameters(get_circuit(c)) == 4*20
    @test collect_blocks(Measure, get_circuit(c))|>length == 6
    @test collect_blocks(HGate, get_circuit(c))|>length == 12
end

@testset "measure" begin
    mr = MeasureResult(4, [1 2 3; 4 5 6])
    @test mr.nx == 3
    @test mr.nbatch == 2
    @test mr.nm == 4
    @test mr[1, 1] == [1, 0]
    @test mr[1, 2] == [0, 0]
    @test mr[1, 3] == [0, 1]
    @test mr[1, 4] == [0, 0]

    @test mr[1, 1] == [1, 0]
    @test mr[2, 1] == [0, 1]
    @test mr[3, 1] == [1, 0]
    @test mr[CartesianIndex(3, 1)] == [1, 0]
    @test mr[3] == [1, 0]
    @test mr[10] == [0, 0]
end

@testset "qpeps machine" begin
    Random.seed!(3)
    c = QPEPSConfig(4, 2, 4, 1)
    qpeps = QPEPSMachine(c, zero_state(12; nbatch=100))
    rt = qpeps.runtime
    @test length(rt.rotblocks) == 4*20
    @test length(rt.mblocks) == 6
    @test length(rt.mbasis) == 6

    @test collect_blocks(I2Gate, qpeps.runtime.circuit) |> length == 6
    chbasis!(qpeps.runtime, Y)
    @test collect_blocks(RotationGate{<:Any,<:Any,<:XGate}, qpeps.runtime.circuit) |> length == 6
    @test collect_blocks(I2Gate, qpeps.runtime.circuit) |> length == 0

    model = J1J2(6, 4; J2=0.5, periodic=false)
    @test nspins(qpeps.config) == nspins(model)
    mres = gensample(qpeps, Z)
    @test mres.nbatch == 100
    @test mres.nx == model.size[1]
    @test mres.nm == model.size[2]
    @test isapprox(energy(qpeps, model), -0.75*12, atol=0.1)
    @test collect_blocks(I2Gate, qpeps.runtime.circuit) |> length == 6
end

@testset "qmps machine" begin
    Random.seed!(3)
    c = QMPSConfig(5, 12, 3)
    qmps = QPEPSMachine(c, zero_state(6; nbatch=100))
    rt = qmps.runtime
    @test length(rt.rotblocks) == 180
    @test length(rt.mblocks) == 16
    @test length(rt.mbasis) == 16

    @test collect_blocks(I2Gate, qmps.runtime.circuit) |> length == 16
    chbasis!(qmps.runtime, Y)
    @test collect_blocks(RotationGate{<:Any,<:Any,<:XGate}, qmps.runtime.circuit) |> length == 16
    @test collect_blocks(I2Gate, qmps.runtime.circuit) |> length == 0

    model = J1J2(4, 4; J2=0.5, periodic=false)
    @test nspins(qmps.config) == nspins(model)
    mres = gensample(qmps, Z)
    @test mres.nbatch == 100
    @test mres.nx == nspins(model)
    @test mres.nm == 1
    @test isapprox(energy(qmps, model), -0.75*8, atol=0.1)
    @test collect_blocks(I2Gate, qmps.runtime.circuit) |> length == 16
end

@testset "j1j2" begin
    j1j2 = J1J2(4; periodic=false, J2=0.5)
    @test get_bonds(j1j2) == [(1, 2, 1.0),(2, 3, 1.0), (3, 4, 1.0), (1,3, 0.5), (2,4, 0.5)]
    j1j2 = J1J2(4; periodic=true, J2=0.5)
    @test get_bonds(j1j2) == [(1, 2, 1.0), (2, 3, 1.0), (3, 4, 1.0), (4, 1, 1.0), (1,3, 0.5), (2,4, 0.5), (3,1, 0.5), (4,2, 0.5)]

    j1j2 = J1J2(3, 3; periodic=false, J2=0.5)
    @test nspins(j1j2) == 9
    prs = [1=>2, 1=>4, 2=>3, 2=>5, 2=>4, 3=>6, 3=>5, 4=>5, 4=>7, 5=>6, 5=>8, 5=>1, 5=>7, 6=>9, 6=>2, 6=>8, 7=>8, 8=>9, 8=>4, 9=>5]
    vs =  [1.0,   1,    1,    1,   0.5,  1.0,   0.5,  1.0,  1.0,  1.0, 1.0,  0.5,   0.5,  1.0,  0.5, 0.5,   1.0, 1.0,  0.5,  0.5]
    tps = [(i,j,v) for ((i,j), v) in zip(prs, vs)]
    @test sort(get_bonds(j1j2)) == sort(tps)
    j1j2 = J1J2(3, 3; periodic=true, J2=0.5)
    @test sort(get_bonds(j1j2)) == sort([(1,2,1.0), (1,4,1.0), (1,9,0.5), (1,6,0.5), (2,3,1.0), (2,5,1.0), (2,7,0.5), (2,4,0.5), (3,1,1.0), (3,6,1.0), (3,8,0.5), (3,5,0.5),
        (4,5,1.0), (4,7,1.0), (4,3,0.5), (4,9,0.5), (5,6,1.0), (5,8,1.0), (5,1,0.5), (5,7,0.5), (6,4,1.0), (6,9,1.0), (6,2,0.5), (6,8,0.5), (7,8,1.0), (7,1,1.0), (7,6,0.5), (7,3,0.5),
        (8,9,1.0), (8,2,1.0), (8,4,0.5), (8,1,0.5), (9,7,1.0), (9,3,1.0), (9,5,0.5), (9,2,0.5)])
end

@testset "train" begin
    model = J1J2(2, 4; J2=0.5, periodic=false)
    h = hamiltonian(model)
    EG = eigen(mat(h) |> Matrix).values[1]
    @test EG/nspins(model) ≈ -0.45522873432070216
end

ext: jl


text:
#! /Users/malcolm/bin/julia -q
#
# Do some command processing here
# usage: ./echo.jl [ -q ] [ -p <port> ]  [-h <host> ]

epipe = "/tmp/epipe"

server = listen(epipe)
while true
  conn = accept(server)
  @async begin
    try
      while true
        s0 = readline(conn)
        s1 = chomp(s0)
        if length(chomp(s1)) > 0
          s2 = reverse(s1)
          if s2 == "."
            println("Done.")
            close(conn)
            exit(0)
          else
            write(conn,string(s2,"\r\n"))
          end
        end
      end
    catch err
      println("Connection lost:  $err")
      exit(1)
    end
  end
end

max_stars_repo_path: Module 3/Chapter09/epipe.jl

max_issues_repo_path: Module 3/Chapter09/epipe.jl

max_forks_repo_path: Module 3/Chapter09/epipe.jl


text:
using DifferentialEquations
using ApproxFun, Sundials
using PyPlot

#using Pkg
#Pkg.add("DiffEqOperators")
using DiffEqOperators, LinearAlgebra
begin
	S = Fourier()
	#S = Chebyshev()
	n = 2^8
	x = points(S, n)
	T = ApproxFun.plan_transform(S, n)
	Ti = ApproxFun.plan_itransform(S, n)

	d = 0.2
	nu = 2e-2
	dt = 1e-3
	Nmax =  (3000+1)
	t0 = Nmax * dt

	u₀ = T*(cos.(x).*(sin.(x).+1));
end

function d_dx(S, n, N)
    D = Derivative(S,n)
    return D[1:N,1:N]
end

Dx(S,m) = d_dx(S,m,n)
L1 = Dx(S,1)
L2 = Dx(S,2)
L4 = Dx(S,4)

A = DiffEqArrayOperator(-d*Diagonal(L2)-nu*Diagonal(L4))

function f(dû,û,tmp,t)
  # Transform u back to point-space
  mul!(tmp,Ti,û)
  tmp= 0.5.*tmp.*tmp
  mul!(tmp, T, tmp)
  mul!(dû,-L1,tmp)
  #mul!(dû,T,tmp) # Transform back to Fourier space
end

prob = SplitODEProblem(A, f, u₀, (0.0,t0), similar(u₀));

sol = solve(prob, ETDRK4(), dt=dt)

plot(Ti*sol(0.001))
plot(Ti*sol(3))
#plot(Ti*sol(3))

#
u = zeros(eltype(sol),(n, Nmax))
for it in 1:Nmax
   print(dt*it)
   u[:,it] = Ti*sol(dt*it)
   #println(u[:,it])
end
size(u)

IM = PyPlot.imshow(abs.(u))
cb = PyPlot.colorbar(IM, orientation="horizontal")
#gcf()

max_stars_repo_name: aetherks/ksequation

max_issues_repo_name: aetherks/ksequation

max_forks_repo_name: aetherks/ksequation


text:
@testset "Core features" begin

  @safetestset "MySQL configuration" begin
    using SearchLight

    conn_info = SearchLight.Configuration.load("mysql_connection.yml")

    @test conn_info["adapter"] == "MySQL"
    @test conn_info["host"] == "127.0.0.1"
    @test conn_info["password"] == "root"
    @test conn_info["config"]["log_level"] == ":debug"
    @test conn_info["config"]["log_queries"] == true
    @test conn_info["port"] == 3306
    @test conn_info["username"] == "root"
    @test conn_info["database"] == "searchlight_tests"
  end;


text:
using FiniteDiff, LinearAlgebra, Test
using SparsityDetection, SparseArrays

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const D = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 8
const X = reshape([i for i in 1:N for j in 1:N],N,N)
const Y = reshape([j for i in 1:N for j in 1:N],N,N)
const α₁ = 1.0.*(X.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
u0 = zeros(N,N,3)

const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N);
# Define the discretized PDE as an ODE function
function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = D*(MyA + AMx)
  @. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

input = rand(N,N,3)
output = similar(input)
sparsity_pattern = jacobian_sparsity(f,output,input,nothing,0.0)
jac_sparsity = Float64.(sparse(sparsity_pattern))

function testf2(u)
	_u = reshape(u,size(input)...)
	du = similar(_u)
	f(du,_u,nothing,0.0)
	vec(du)
end
findifjac = FiniteDiff.finite_difference_jacobian(testf2,input)
@test all((abs.(findifjac) .> 0) .== (jac_sparsity .> 0))


text:
#-------------------------------------------------------
#  set up FUND and FAIR models
#-------------------------------------------------------

using MimiFAIR
using MimiFUND
using Mimi

include(joinpath(@__DIR__, "helper.jl")) # path relative to this file

# set some constants
year = 2030
pulse_size = 1e7

plots_output_dir = joinpath(@__DIR__, "..", "output")
isdir(plots_output_dir) || mkpath(plots_output_dir)

# start with default MimiFUND model
m = get_fundfair(;rcp_scenario = "RCP85", FAIR_first = 1765, FAIR_last = 2500, FUND_first = 1950, FUND_last = 2500)
run(m)

#-------------------------------------------------------
# create marginal model and perturb
#-------------------------------------------------------

pulse_size = 1e7
mm = Mimi.create_marginal_model(m, pulse_size)

# A component for an emissions pulse to be used in social cost calculations. Computes the `output` vector by adding
#   `add` to `input`. This is similar to the Mimi.adder component, except that it allows missing values to be passed through.

@defcomp emissionspulse begin
    add    = Parameter(index=[time])
    input  = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.output[t] = Mimi.@allow_missing(p.input[t]) + p.add[t]
    end
end

# Add additional emissions to mm.modified
add_comp!(mm.modified, emissionspulse, after = :emissions)
nyears = length(Mimi.time_labels(m))
addem = zeros(nyears)
baseyear = 1765

# pulse is spread over ten years, and emissions components is in Mt so divide by 1e7, and convert from CO2 to C because emissions component is in MtC
addem[(year - baseyear + 1):(year - baseyear + 1) + 9] .= pulse_size / 1e7 * 12/44

set_param!(mm.modified, :emissionspulse, :add, addem)

# Reconnect the appropriate emissions in m
connect_param!(mm.modified, :emissionspulse, :input, :emissions, :mco2, zeros(nyears), backup_offset = 1)
connect_param!(mm.modified, :multiplier, :input, :emissionspulse, :output)
run(mm)

# check results
# mm.modified
# mm.modified[:emissions, :mco2] # 736 element array, missings for pre-1950 years
# mm.modified[:emissionspulse, :input] # 736 element array, 0s for pre-1950 years
# mm.modified[:emissionspulse, :input][266:275]
# mm.modified[:emissionspulse, :add] # 736 element array, all zeros except 10 x 0.2727 between years corresponding to 2030 and 2040 (mm.modified[:emissionspulse, :add][266:275])
# mm.modified[:emissionspulse, :add][266:275]
# mm.modified[:emissionspulse, :output] # 736 element array, 0s for pre-1950 years
# mm.modified[:emissionspulse, :output][266:275]

# mm.modified[:emissionspulse, :output][266:275] == mm.modified[:emissionspulse, :input][266:275] + mm.modified[:emissionspulse, :add][266:275] # true

# mm.modified[:multiplier, :input]
# mm.modified[:multiplier, :input] == mm.modified[:emissionspulse, :output] # true
# mm.modified[:multiplier, :multiply]
# mm.modified[:multiplier, :output] # missings for pre-1950 years?
# mm.modified[:multiplier, :output][187:end] == mm.modified[:multiplier, :input][187:end] .* 0.001 # true
# mm.modified[:multiplier, :input][187:end]

# mm.modified[:co2_cycle, :E_CO₂]
# mm.modified[:co2_cycle, :E_CO₂][187:end] == mm.modified[:multiplier, :output][187:end] # true
# mm.modified[:co2_cycle, :E_CO₂][1:186] == FAIR_CO₂_backup[1:186]

#-------------------------------------------------------
# calculate SCC
#-------------------------------------------------------

marginaldamage = mm[:impactaggregation, :loss]
# marginaldamage = mm.modified[:impactaggregation, :loss] - mm.base[:impactaggregation, :loss]

# calculate discount factors
prtp = 0.05
new_years = collect(1765:1:2500)
year_index = findfirst(isequal(year), new_years)

df = zeros(length(new_years), 16)
for i in 1:length(new_years)
    if i >= year_index
        df[i,:] .= 1/(1+prtp)^(i-year_index)
    end
end

# calculate SCC
scc = sum(skipmissing(marginaldamage .* df))

## compare with FUND SCC
mfund = MimiFUND.get_model()
run(mfund)
scc_fund = MimiFUND.compute_scco2(mfund, year = 2030, eta = 0., prtp = 0.05, equity_weights = false) #...wut

#-------------------------------------------------------
#  plot results
#-------------------------------------------------------

using VegaLite
using DataFrames

# get marginal FUND model
mmfund = MimiFUND.get_marginal_model(mfund; year = year, gas = :CO2, pulse_size = pulse_size)
run(mmfund)

# get marginal FAIR model
mfair = MimiFAIR.get_model()
run(mfair)

# plot marginal damages
fundfair_md = sum(mm[:impactaggregation, :loss], dims = 2)[:,1];
# fundfairvals_perturbed = sum(mm.modified[:impactaggregation, :loss], dims = 2)[:,1];
# fundfairvals_base = sum(mm.base[:impactaggregation, :loss], dims = 2)[:,1];

fund_md = vcat(
    fill(missing, length(FAIR_first:FUND_first)-1),
    sum(mmfund[:impactaggregation, :loss][1:length(FUND_first:FAIR_last),:], dims = 2)
)
fund_md = fund_md[:,1]; # make a vector

df = DataFrame(
    :Year => Mimi.time_labels(m),
    :FUND => fund_md,
    :FUNDFAIR => fundfair_md)
    # :FUNDFAIR_base => fundfairvals_base,
    # :FUNDFAIR_perturbed => fundfairvals_perturbed)

# stack(df, [:FUND, :FUNDFAIR_base, :FUNDFAIR_perturbed]) |>
stack(df, [:FUND, :FUNDFAIR]) |>

@vlplot(
    :line,
    x = :Year,
    y = :value,
    color = :variable
) |>
save("fundfair_marginaldamages.png")

# plot temperature delta
fundfair_tempdelta =  mm[:temperature, :T]
fund_tempdelta = vcat(
    fill(missing, length(FAIR_first:FUND_first)-1),
    mmfund[:climateco2cycle, :temp][1:length(FUND_first:FAIR_last),:]
)
fund_tempdelta = fund_tempdelta[:,1]; # make a vector

df = DataFrame(
    :Year => Mimi.time_labels(m),
    :FUND => fund_tempdelta,
    # :FAIR => fairvals,
    :FUNDFAIR => fundfair_tempdelta
)
stack(df, [:FUND, :FUNDFAIR]) |>
# stack(df, [:FUND, :FAIR, :FUNDFAIR]) |>

@vlplot(
    :line,
    x = :Year,
    y = :value,
    color = :variable
) |>
save("fundfair_tempdelta.png")

## using plots instead of vegalite

using Plots

# temp deltas for 1GtCO2 perturbation
fund_tempdelta[186:736]
fundfair_tempdelta[186:736]

# years = Mimi.time_labels(m)
years = collect(1950:1:2500)

Plots.plot(years, fund_tempdelta[186:736],
        title = "Temperature Deltas from 1tCO2 Emissions Pulse",
        xaxis = "Year",
        ylabel = "Degrees C",
        label = "FUND")
Plots.plot!(years, fundfair_tempdelta[186:736], label = "FUNDFAIR (RCP8.5)")

Plots.savefig(joinpath(plots_output_dir, "fundfair_tempdelta"))

# marginal damages

Plots.plot(years, fund_md[186:736],
        title = "Marginal Damages from 1tCO2 Emissions Pulse",
        xaxis = "Year",
        ylabel = "\$",
        label = "FUND")
Plots.plot!(years, fundfair_md[186:736], label = "FUNDFAIR (RCP8.5)")

Plots.savefig(joinpath(plots_output_dir, "fundfair_marginaldamages"))

# emissions
fundfair_base_emissions = mm.base[:co2_cycle, :E_CO₂]
fund_base_emissions = mfund[:emissions, :mco2][1:551]/1000
fair_emissions = mfair[:co2_cycle, :E_CO₂][186:736]

Plots.plot(years, fundfair_base_emissions[186:736],
        title = "CO2 Emissions",
        xaxis = "Year",
        ylabel = "GtC / year",
        label = "FUNDFAIR / FUND")
# Plots.plot!(years, fund_base_emissions, label = "FUND Base")
Plots.plot!(years, fair_emissions, label = "FAIR (RCP8.5)")

Plots.savefig("fundfair_emissions")

# atmospheric CO2
fund_acco2 = mfund[:climateco2cycle, :acco2][1:551]
fundfair_acco2 = mm.base[:co2_cycle, :C][186:736]
fair_acco2 = mfair[:co2_cycle, :C][186:736]

Plots.plot(years, fund_acco2,
        title = "Atmospheric CO2 Concentration",
        xaxis = "Year",
        ylabel = "ppm",
        label = "FUND")
Plots.plot!(years, fundfair_acco2, label = "FUNDFAIR")
Plots.plot!(years, fair_acco2, label = "FAIR (RCP8.5)")

Plots.savefig("fundfair_acco2")

# radiative forcing

mfair85 = MimiFAIR.get_model(rcp_scenario = "RCP85")
run(mfair85)

mfair45 = MimiFAIR.get_model(rcp_scenario = "RCP45")
run(mfair45)

mfair26 = MimiFAIR.get_model(rcp_scenario = "RCP26")
run(mfair26)

fund_rf = mfund[:climateforcing, :radforc][1:551]
fundfair_rf_base = mm.base[:temperature, :F][186:736]
fundfair_rf_perturbed = mm.modified[:temperature, :F][186:736]

fair85_rf = mfair85[:temperature, :F][186:736]
fair45_rf = mfair45[:temperature, :F][186:736]
fair26_rf = mfair26[:temperature, :F][186:736]

Plots.plot(years, fair85_rf,
        title = "Total Radiative Forcing",
        xaxis = "Year",
        ylabel = "",
        label = "FAIR (RCP8.5)",
        legend = :topleft)
Plots.plot!(years, fair45_rf, label = "FAIR (RCP4.5)")
Plots.plot!(years, fair26_rf, label = "FAIR (RCP2.6)")
Plots.plot!(years, fundfair_rf_base, label = "FUNDFAIR (Base)")
Plots.plot!(years, fundfair_rf_perturbed, label = "FUNDFAIR (Perturbed)")
Plots.plot!(years, fund_rf, label = "FUND")

Plots.savefig(joinpath(plots_output_dir, "fundfair_acco2"))

# change in rf
fundfair_rf_delta = mm[:temperature, :F][186:736]
fund_rf_delta = mmfund[:climateforcing, :radforc][1:551]

Plots.plot(years, fund_rf_delta,
        title = "Change in Radiative Forcing from 1tCO2 Pulse",
        xaxis = "Year",
        ylabel = "",
        label = "FUND")
        # legend = :topright)
Plots.plot!(years, fundfair_rf_delta, label = "FUNDFAIR")

Plots.savefig(joinpath(plots_output_dir, "fundfair_rfdelta"))

# change in atmospheric CO2
fund_acco2_delta = mmfund[:climateco2cycle, :acco2][1:551]
fundfair_acco2_delta = mm[:co2_cycle, :C][186:736]

Plots.plot(years, fund_acco2_delta,
        title = "Change in Atmospheric CO2 Concentration",
        xaxis = "Year",
        ylabel = "ppm",
        label = "FUND")
Plots.plot!(years, fundfair_acco2_delta, label = "FUNDFAIR")

Plots.savefig(joinpath(plots_output_dir, "fundfair_acco2_delta"))


text:
"""
    sphMapping( Pos, HSML, M,
                Rho, Bin_Quant,
                Weights=Rho;
                param::mappingParameters,
                kernel::SPHKernel,
                show_progress::Bool=true,
                parallel::Bool=false,
                reduce_image::Bool=true,
                filter_particles::Bool=true,
                dimensions::Int=2)

Maps the data in `Bin_Quant` to a grid. Parameters of mapping are supplied in
`param` and the kernel to be used in `kernel`.

# Arguments
- `Pos`: Matrix (3xNpart) with particle positions.
- `HSML`: Array with particle hsml.
- `M`: Array with particle masses.
- `Rho`: Array with particle densities.
- `Bin_Quant`: Array with particle quantity to be mapped.
- `Weights`: Array with weights. Defaults to density-weighted.
- `kernel::SPHKernel`: Kernel object to be used.
- `show_progress::Bool=true`: Show progress bar.
- `parallel::Bool=true`: Run on multiple processors.
- `reduce_image::Bool=true`: If weights need to be applied or not. Set to `false` for [`part_weight_one`](@ref) and [`part_weight_physical`](@ref).
- `filter_particles::Bool=true`: Find the particles that are actually contained in the image.
- `dimensions::Int=2`: Number of mapping dimensions (2 = to grid, 3 = to cube).
"""
function sphMapping(Pos, HSML, M,
                    Rho, Bin_Quant,
                    Weights=Rho;
                    param::mappingParameters,
                    kernel::SPHKernel,
                    show_progress::Bool=true,
                    parallel::Bool=false,
                    reduce_image::Bool=true,
                    filter_particles::Bool=true,
                    dimensions::Int=2)

    # store number of input particles
    N_in = size(Bin_Quant,1)

    # First check if all particles are centered around 0 and shift them if they are not
    if show_progress
        @info "Centering on [0.0, 0.0, 0.0]"
        t1 = time_ns()
    end
    Pos, par = check_center_and_move_particles(Pos, param)

    if show_progress
        t2 = time_ns()
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    # filter particles if they are contained in the image
    if show_progress
        @info "Filtering particles..."
        t1 = time_ns()
    end

    if filter_particles
        p_in_image = filter_particles_in_image(Pos, HSML, param)
    else
        p_in_image = trues(N_in)
    end

    if show_progress
        t2 = time_ns()
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    # if this is not a float it has units, which need to be stripped
    if !(typeof(Bin_Quant[1]) <: AbstractFloat)

        if show_progress
            @info "Stripping units..."
            t1 = time_ns()
        end

        # allocate reduced arrays
        x       = ustrip.(Pos[:,p_in_image])
        hsml    = ustrip.(HSML[p_in_image])
        m       = ustrip.(M[p_in_image])
        rho     = ustrip.(Rho[p_in_image])
        bin_q   = ustrip.(Bin_Quant[p_in_image])
        weights = ustrip.(Weights[p_in_image])

    else
        if show_progress
            @info "Assigning arrays..."
            t1 = time_ns()
        end

        # allocate reduced arrays
        x       = Pos[:,p_in_image]
        hsml    = HSML[p_in_image]
        m       = M[p_in_image]
        rho     = Rho[p_in_image]
        bin_q   = Bin_Quant[p_in_image]
        weights = Weights[p_in_image]
    end

    if show_progress
        t2 = time_ns()
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    N_map = length(m)

    @info "Particles in image: $N_map / $N_in"

    if show_progress
        @info "Mapping..."
        t1 = time_ns()
    end

    if (dimensions == 2)

        if !parallel

            image = sphMapping_2D(x, hsml, m, rho, bin_q, weights;
                                param=par, kernel=kernel,
                                show_progress=show_progress)

            if show_progress
                t2 = time_ns()
                @info "  elapsed: $(output_time(t1,t2)) s"
            end

            if !reduce_image
                image[:,2] .= 1.0
            end

            return reduce_image_2D( image,
                            param.Npixels[1], param.Npixels[2] )
        else
            @info "Running on $(nworkers()) cores."

            # Number of particles
            N = length(m)

            # allocate an array of Future objects
            futures = Array{Future}(undef, nworkers())

            # 'Domain decomposition':
            # calculate array slices for each worker
            batch = domain_decomposition(N, nworkers())

            # start remote processes
            for (i, id) in enumerate(workers())
                futures[i] = @spawnat id sphMapping_2D(x[:,batch[i]], hsml[batch[i]],
                                                        m[batch[i]], rho[batch[i]],
                                                        bin_q[batch[i]], weights[batch[i]];
                                                        param=par, kernel=kernel,
                                                        show_progress=false)
            end

            image = sum(fetch.(futures))

            if show_progress
                t2 = time_ns()
                @info "  elapsed: $(output_time(t1,t2)) s"
            end

            if !reduce_image
                image[:,2] .= 1.0
            end

            return reduce_image_2D( image,
                        param.Npixels[1], param.Npixels[2] )
        end

    elseif (dimensions == 3 )
        if !parallel
            image = sphMapping_3D(x, hsml, m, rho, bin_q, weights;
                                param=par, kernel=kernel,
                                show_progress=show_progress)

            if show_progress
                t2 = time_ns()
                @info "  elapsed: $(output_time(t1,t2)) s"
            end

            if !reduce_image
                image[:,2] .= 1.0
            end

            return reduce_image_3D( image,
                        param.Npixels[1], param.Npixels[2], param.Npixels[3] )
        else
            @info "Running on $(nworkers()) cores."

            N = length(m)
            futures = Array{Future}(undef, nworkers())

            # 'Domain decomposition':
            # calculate array slices for each worker
            batch = domain_decomposition(N, nworkers())

            # start remote processes
            for (i, id) in enumerate(workers())
                futures[i] = @spawnat id sphMapping_3D(x[:,batch[i]], hsml[batch[i]],
                                                    m[batch[i]], rho[batch[i]],
                                                    bin_q[batch[i]], weights[batch[i]];
                                                    param=par, kernel=kernel,
                                                    show_progress=false)
            end

            # get and reduce results
            image = sum(fetch.(futures))

            if show_progress
                t2 = time_ns()
                @info "  elapsed: $(output_time(t1,t2)) s"
            end

            if !reduce_image
                image[:,2] .= 1.0
            end

            return reduce_image_3D( image,
                        param.Npixels[1], param.Npixels[2], param.Npixels[3] )
        end
    end

end


text:
@testset "Random Strings" begin
  @test isa(rand(String), String)
  @test length(rand(String)) == 8
  @test length(rand(String, 10)) == 10

  @test_broken length(rand(String)) == 10
  @test_broken rand(String) == "abcdefgh"
end


text:
export GMRES, BiCGStabl, Exact

abstract type CorrectionSolver end

struct GMRES <: CorrectionSolver
    n::Int
    iterations::Int
    tolerance::Float64
    function GMRES(n::Int; iterations::Int = 5, tolerance::Float64 = 1e-6)
        new(n, iterations, tolerance)
    end
end

struct BiCGStabl <: CorrectionSolver
    max_mv_products::Int
    tolerance::Float64
    l::Int
    r_shadow
    rs
    us
    γ
    M

    function BiCGStabl(n::Int; l::Int = 2, max_mv_products::Int = 5, tolerance::Float64 = 1e-6)
        r_shadow = rand(ComplexF64, n)
        rs = Matrix{ComplexF64}(undef, n, l + 1)
        us = Matrix{ComplexF64}(undef, n, l + 1)
        γ = Vector{ComplexF64}(undef, l)
        M = Matrix{ComplexF64}(undef, l + 1, l + 1)

        new(max_mv_products, tolerance, l, r_shadow, rs, us, γ, M)
    end
end

struct Exact <: CorrectionSolver end

function solve_deflated_correction!(solver::Exact, A, x, Q, θ, r::AbstractVector, tol)
    # The exact solver is mostly useful for testing Jacobi-Davidson
    # method itself and should result in quadratic convergence.
    # However, in general the correction equation should be solved
    # only approximately in a fixed number of GMRES / BiCGStab
    # iterations to reduce computation costs.

    # Here we solve the augmented system
    # [A - θI, Q; Q' 0][t; y] = [-r; 0] for t,
    # which is equivalent to solving (I - QQ')(A - θI)(I - QQ')t = r
    # for t ⟂ Q.

    n = size(A, 1)
    m = size(Q, 2)
    Ã = [(A - θ * speye(n)) Q; Q' zeros(m, m)]
    rhs = [r; zeros(m, 1)]
    copyto!(x, (Ã \ rhs)[1 : n])
end

function solve_deflated_correction!(solver::BiCGStabl, A, x, Q, θ, r::AbstractVector, tol)
    n = size(A, 1)
    T = eltype(x)

    # y = (I - QQ')(A - θI)x
    matrix = LinearMap{T}((y, x) -> begin
        mul!(y, A, x)
        axpy!(-θ, x, y)
        just_orthogonalize!(Q, x, DGKS)
    end, n, ismutating = true)

    mv_products = 0

    residual = view(solver.rs, :, 1)

    # The initial residual rs[:, 1] = b - A * x
    copyto!(residual, r)
    fill!(view(solver.us, :, 1), zero(T))
    fill!(x, zero(T))

    ω = σ = one(T)

    just_orthogonalize!(Q, residual, DGKS)
    nrm = norm(residual)

    # Stopping condition based on relative tolerance.
    reltol = nrm * tol

    iterable = BiCGStabIterable(matrix, solver.l, x, solver.r_shadow, solver.rs, solver.us,
        solver.max_mv_products, mv_products, reltol, nrm,
        Identity(),
        solver.γ, ω, σ, solver.M
    )

    for res = iterable
        # println(res)
    end

    nothing
end

function solve_generalized_correction_equation!(solver::Exact, A, B, preconditioner, x, Q, Z, precZ, QZ, α, β, r, spare_vector, tol)
  n = size(A, 1)
  m = size(Q, 2)
  # Assuming both A and B are sparse while Q and Z are dense, let's try to avoid constructing a huge dense matrix.
  # Let C = β * A - α * B

  # We have to solve:
  # |C  Z| |t| = |-r|
  # |Q' O| |z|   |0 |

  # Use the Schur complement trick with S = -Q' * inv(C) * Z
  # |C  Z| = |I            O| |C Z|
  # |Q' O|   |Q' * inv(C)  I| |O S|

  # And solve two systems with inv(C) occurring multiple times.
  C = β * A - α * B # Completely sparse
  y = Q' * (C \ r)
  S = Q' * (C \ Z) # Schur complement
  z = -S \ y
  x .= C \ (-r - Z * z)
end

struct QZpreconditioner{T}
    Q
    Z
    QZ
    preconditioner
end

function ldiv!(QZ::QZpreconditioner{T}, x) where {T}
    # x ← preconditioner \ x
    ldiv!(QZ.preconditioner, x)

    # x ← (I - Z inv(Q'Z) Q')x
    h = QZ.Q' * x
    ldiv!(QZ.QZ.LU, h)
    gemv!('N', -one(T), QZ.Z, h, one(T), x)
end

"""
Solve the problem ``(I - ZZ^*)(βA - αB)(I - QQ^*)t = b``
this is simplified in Krylov subspaces as solving
``(βA - αB)t = b``
with the left preconditioner `Pl = (I - Z inv(Q'Z) Q')`
"""
function solve_generalized_correction_equation!(solver::GMRES, A, B, preconditioner, x, Q, Z, precZ, QZ, α, β, r, spare_vector, tol)
    n = size(A, 1)
    @assert n == solver.n
    T = eltype(x)

    # Preconditioner
    Pl = QZpreconditioner{T}(Q, precZ, QZ, preconditioner)

    # Linear operator
    Ax_minus_Bx = LinearMap{T}((y, x) -> begin
        # y = (βA - αB) * x
        mul!(y, B, x, -α, false)
        mul!(y, A, x, β, true)
    end, n, ismutating = true)

    # Start with a zero guess
    fill!(x, zero(T))

    iterable = IterativeSolvers.gmres_iterable!(x, Ax_minus_Bx, r, Pl=Pl, tol=tol,
                                                restart=solver.iterations,
                                                initially_zero=true)

    for res = iterable
        # println(res)
    end

    nothing
end

function solve_generalized_correction_equation!(solver::BiCGStabl, A, B, preconditioner, x, Q, Z, precZ, QZ, α, β, r, spare_vector, tol)
    n = size(A, 1)
    T = eltype(x)

    Pl = QZpreconditioner{T}(Q, precZ, QZ, preconditioner)

    Ax_minus_Bx = LinearMap{T}((y, x) -> begin
        # y = (βA - αB) * x
        mul!(y, B, x, -α, false)
        mul!(y, A, x, β, true)
    end, n, ismutating = true)

    mv_products = 0

    residual = view(solver.rs, :, 1)

    # The initial residual rs[:, 1] = b - A * x
    copyto!(residual, r)
    fill!(view(solver.us, :, 1), zero(T))
    fill!(x, zero(T))

    # Apply the left preconditioner
    ldiv!(Pl, residual)

    ω = σ = one(T)

    nrm = norm(residual)

    # Stopping condition based on relative tolerance.
    reltol = nrm * tol

    iterable = BiCGStabIterable(Ax_minus_Bx, solver.l, x, solver.r_shadow, solver.rs, solver.us,
        solver.max_mv_products, mv_products, reltol, nrm,
        Pl,
        solver.γ, ω, σ, solver.M
    )

    for res = iterable
        # println(res)
    end

    nothing
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to create the REPL mode `pager`.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   This code was adapted from the Pkg.jl project.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _create_pager_repl_mode(repl::REPL.AbstractREPL, main::LineEdit.Prompt)

Create the REPL mode `pager`. `repl` must be the active REPL, and `main` must be
the main REPL mode (julia prompt).

"""
function _create_pager_repl_mode(repl::REPL.AbstractREPL, main::LineEdit.Prompt)
    # Prompt of the pager mode
    # ==========================================================================

    # In this case, we will use the same completion provider as the julia
    # prompt.
    tp_mode = LineEdit.Prompt(
        _tp_mode_prompt;
        prompt_prefix = repl.options.hascolor ? Base.text_colors[:magenta] : "",
        prompt_suffix = "",
        complete = REPL.REPLCompletionProvider(),
        sticky = true
    )

    # This function is called when the user hits return after typing a command.
    tp_mode.on_done = (s, buf, ok) -> begin
        ok || return REPL.transition(s, :abort)

        # Take the input command.
        input = String(take!(buf))
        REPL.reset(repl)

        # Process the input command inside the pager mode.
        _tp_mode_do_cmd(repl, input)

        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main)
    end

    tp_mode.repl = repl

    # Check if the expression is incomplete, and, if so, request for another
    # line.
    tp_mode.on_enter = REPL.return_callback

    # Let the command history be equal to the julia prompt.
    hp = main.hist
    hp.mode_mapping[:pager] = tp_mode
    tp_mode.hist = hp

    # Create the sub mode: pager help.
    tp_help_mode = _create_pager_help_repl_mode(repl, main, tp_mode)

    # Key mappings
    # ==========================================================================

    # We want to support all the default keymap prefixes.
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, tp_mode)

    # We also want to support reverse searching.
    search_promt, skeymap = LineEdit.setup_search_keymap(hp)

    # Key mappings used in the pager mode:
    mk = REPL.mode_keymap(main)

    # Assign `?` as the key map to switch to pager help mode.
    help_mode_transition_keymap = Dict{Any, Any}(
        '?' => function(s, args...)
            # We must only switch to pager mode if `|` is typed at the beginning
            # of the line.
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, tp_help_mode) do
                    LineEdit.state(s, tp_help_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, '?')
            end
        end
    )

    tp_mode_keymaps = Dict{Any, Any}[
        mk,
        prefix_keymap,
        skeymap,
        help_mode_transition_keymap,
        LineEdit.history_keymap,
        LineEdit.default_keymap,
        LineEdit.escape_defaults,
    ]

    tp_mode.keymap_dict = LineEdit.keymap(tp_mode_keymaps)

    return tp_mode
end

"""
    _create_pager_help_repl_mode( repl::REPL.AbstractREPL, main::LineEdit.Prompt, tp_mode::LineEdit.Prompt )

Create the REPL mode `pager help`. `repl` must be the active REPL, `main` must
be the main REPL mode (julia prompt), and `tp_mode` must be the REPL mode
`pager`.

"""
function _create_pager_help_repl_mode(
    repl::REPL.AbstractREPL,
    main::LineEdit.Prompt,
    tp_mode::LineEdit.Prompt
)
    # Prompt of the pager help mode
    # ==========================================================================

    tp_help_mode = LineEdit.Prompt(
        _tp_help_mode_prompt;
        prompt_prefix = repl.options.hascolor ? Base.text_colors[:yellow] : "",
        prompt_suffix = "",
        complete = REPL.REPLCompletionProvider(),
        sticky = false
    )

    tp_help_mode.on_done = (s, buf, ok) -> begin
        # Take the input command.
        input = String(take!(buf))
        REPL.reset(repl)

        # Process the input command inside the pager mode.
        _tp_help_mode_do_cmd(repl, input)

        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, tp_mode)
    end

    # Key mappings
    # ==========================================================================

    hp = main.hist
    hp.mode_mapping[:pager_help] = tp_help_mode
    tp_help_mode.hist = hp

    # We want to support all the default keymap prefixes.
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, tp_mode)

    # We also want to support reverse searching.
    search_promt, skeymap = LineEdit.setup_search_keymap(hp)

    # Key mappings used in the pager mode:
    mk = REPL.mode_keymap(main)

    tp_help_mode.repl = repl

    tp_help_mode_keymaps = Dict{Any, Any}[
        mk,
        prefix_keymap,
        skeymap,
        LineEdit.history_keymap,
        LineEdit.default_keymap,
        LineEdit.escape_defaults,
    ]

    tp_help_mode.keymap_dict = LineEdit.keymap(tp_help_mode_keymaps)

    return tp_help_mode
end

"""
    _init_pager_repl_mode(repl::AbstractREPL)

Initialize the pager mode in the `repl`.

"""
function _init_pager_repl_mode(repl::AbstractREPL)
    # Get the main REPL mode (julia prompt).
    main_mode = repl.interface.modes[1]

    # Create the pager mode.
    tp_mode = _create_pager_repl_mode(repl, main_mode)

    # Add the new mode to the REPL interfaces.
    push!(repl.interface.modes, tp_mode)

    # Assign `|` as the key map to switch to pager mode.
    keymap = Dict{Any, Any}(
        '|' => function(s, args...)
            # We must only switch to pager mode if `|` is typed at the beginning
            # of the line.
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, tp_mode) do
                    LineEdit.state(s, tp_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, '|')
            end
        end
    )

    # Add the key map that initialize the pager mode to the default REPL key
    # mappings.
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, keymap)

    return nothing
end

################################################################################
#                     Command treatment for the REPL modes
################################################################################

"""
    _tp_mode_do_cmd(repl::REPL.AbstractREPL, input::String)

Execute the actions when a command has been received in the REPL mode `pager`.
`repl` must be the active REPL, and `input` is a string with the command.
"""
function _tp_mode_do_cmd(repl::REPL.AbstractREPL, input::String)
    if !isinteractive() && !PRINTED_REPL_WARNING[]
        @warn "The parger mode is intended for interaction use only, and should not be used from scripts."
        PRINTED_REPL_WARNING[] = true
    end

    # The `stdout` will be redirected inside the try/catch. Hence, we need to
    # store the old one to restore it if everything fails.
    old_stdout = stdout

    try
        # Create a buffer that will replace `stdout`.
        buf = IOBuffer()
        io = IOContext(
            IOContext(buf, stdout),
            :limit => false
        )

        # Redirect `stdout` to the new buffer.
        Base.eval(:(stdout = $io))

        # First, we need to split the buffer into lines.
        lines = split(input, '\n', keepempty = true)

        # Variable to assemble the command, which can have multiple lines.
        cmd = ""

        # Loop throught the lines.
        for i ∈ 1:length(lines)
            cmd *= lines[i] * "\n"
            ast = Base.parse_input_line(cmd)

            # If the command is incomplete, we need to wait for another line.
            ast !== nothing && ast.head == :incomplete && continue

            # Evaluate the input.
            result = Core.eval(Main, ast)

            # Make `ans` be the last evaluated expression.
            #
            # This code was copied from the function `eval_user_input` in REPL
            # stdlib.
            ccall(:jl_set_global, Cvoid, (Any, Any, Any), Main, :ans, result)

            # If the user added `;` at the end of the command, then we should
            # not show the output.
            if !REPL.ends_with_semicolon(cmd)
                # If the output is not `nothing`, call `show` with
                # `MIME("text/plain")` to render the object.
                if result !== nothing
                    Base.invokelatest(show, stdout, MIME("text/plain"), result)
                    write(stdout, '\n')
                end
            end

            # Clear the current command to receive the next one.
            cmd = ""
        end

        # Restore the old stdout.
        Base.eval(:(stdout = $old_stdout))

        # Take everything and display in the pager using `auto` mode. In this
        # case, the pager will only be called if there is not space in the
        # display to show everything.
        pager(String(take!(buf)); auto = true)

        close(io)
    catch err
        Base.display_error(repl.t.err_stream, err, Base.catch_backtrace())
    finally
        Base.eval(:(stdout = $old_stdout))
    end

    return nothing
end

"""
    _tp_help_mode_do_cmd(repl::REPL.AbstractREPL, input::String)

Execute the actions when a command has been received in the REPL mode `pager
help`. `repl` must be the active REPL, and `input` is a string with the command.
"""
function _tp_help_mode_do_cmd(repl::REPL.AbstractREPL, input::String)
    # We do not need to verify if we are in a interactive environment because
    # this mode is only accessible throught pager mode, which already checks it.
    try
        # Create a buffer that will replace `stdout`.
        buf = IOBuffer()
        io = IOContext(
            IOContext(buf, stdout),
            :limit => false
        )

        # Get the AST that generates the help.
        ast = Base.invokelatest(REPL.helpmode, io, input)

        # Evaluate the AST, which returns a Markdown object.
        response = Core.eval(Main, ast)

        # Render the output.
        show(io, MIME("text/plain"), response)
        write(io, '\n')

        # Take everything and display in the pager using `auto` mode. In this
        # case, the pager will only be called if there is not space in the
        # display to show everything.
        pager(String(take!(buf)); auto = true)

        close(io)
    catch err
        Base.display_error(repl.t.err_stream, err, Base.catch_backtrace())
    end

    return nothing
end

################################################################################
#                                   Prompts
################################################################################

_tp_mode_prompt() = "pager> "
_tp_help_mode_prompt() = "pager?> "


text:
# General case - to be used with bc::NoBC or m::PrescribedXModels
function soil_boundary_flux!(
    nf,
    bc::AbstractBoundaryConditions,
    m::AbstractSoilComponentModel,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
) end

function soil_boundary_state!(
    nf,
    bc::AbstractBoundaryConditions,
    m::AbstractSoilComponentModel,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)

end

# Dirichlet Methods for SoilHeat and SoilWater

"""
    function soil_boundary_state!(
        nf,
        bc::Dirichlet,
        water::SoilWaterModel,
        land::LandModel,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Dirichlet-type method for `soil_boundary_state!` for the
`SoilWaterModel`.
"""
function soil_boundary_state!(
    nf,
    bc::Dirichlet,
    water::SoilWaterModel,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    bc_function = bc.state_bc
    state⁺.soil.water.ϑ_l = bc_function(aux⁻, t)

end

"""
    function soil_boundary_state!(
        nf,
        bc::Dirichlet,
        heat::SoilHeatModel,
        land::LandModel,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Dirichlet-type method for `soil_boundary_state!` for the
`SoilHeatModel`.
"""
function soil_boundary_state!(
    nf,
    bc::Dirichlet,
    heat::SoilHeatModel,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    bc_function = bc.state_bc
    param_set = parameter_set(land)
    ϑ_l, θ_i = get_water_content(land.soil.water, aux⁻, state⁻, t)
    θ_l = volumetric_liquid_fraction(ϑ_l, land.soil.param_functions.porosity)
    ρc_s = volumetric_heat_capacity(
        θ_l,
        θ_i,
        land.soil.param_functions.ρc_ds,
        param_set,
    )

    ρe_int_bc =
        volumetric_internal_energy(θ_i, ρc_s, bc_function(aux⁻, t), param_set)

    state⁺.soil.heat.ρe_int = ρe_int_bc
end

# Neumann conditions for SoilHeat and SoilWater

"""
    function soil_boundary_flux!(
        nf,
        bc::Neumann,
        water::SoilWaterModel,
        land::LandModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Neumann method for `soil_boundary_flux!` for the
`SoilWaterModel`.
"""
function soil_boundary_flux!(
    nf,
    bc::Neumann,
    water::SoilWaterModel,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    bc_function = bc.scalar_flux_bc
    # on faces at the top and bottom of the domain,
    # the vector n̂ can point in either ± ẑ, depending on which side of the domain
    # we are on. The user supplies a scalar flux F, such that the flux at the boundary
    # is assumed to be F⃗ = F ẑ - i.e. the user doesn't need to worry about the normal vector.
    # so, we take the absolute value of n̂ here.
    # The minus sign is because the condition is applied on minus the flux.
    # the same argument applies for the other directions.
    diff⁺.soil.water.K∇h = abs.(n̂) * (-bc_function(aux⁻, t))
end

"""
    function soil_boundary_flux!(
        nf,
        bc::Neumann,
        heat::SoilHeatModel,
        land::LandModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Neumann method for `soil_boundary_flux!` for the
`SoilHeatModel`.
"""
function soil_boundary_flux!(
    nf,
    bc::Neumann,
    heat::SoilHeatModel,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    bc_function = bc.scalar_flux_bc
    # on faces at the top and bottom of the domain,
    # the vector n̂ can point in either ± ẑ, depending on which side of the domain
    # we are on. The user supplies a scalar flux F, such that the flux at the boundary
    # is assumed to be F⃗ = F ẑ - i.e. the user doesn't need to worry about the normal vector.
    # so, we take the absolute value of n̂ here.
    # The minus sign is because the condition is applied on minus the flux.
    # the same argument applies for the other directions.
    diff⁺.soil.heat.κ∇T = abs.(n̂) * (-bc_function(aux⁻, t))
end

# SurfaceDriven conditions for for SoilHeat and SoilWater

"""
    function soil_boundary_flux!(
        nf,
        bc::SurfaceDrivenWaterBoundaryConditions,
        water::SoilWaterModel,
        land::LandModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Surface Driven BC method for `soil_boundary_flux!` for the
`SoilWaterModel`.
"""
function soil_boundary_flux!(
    nf,
    bc::SurfaceDrivenWaterBoundaryConditions,
    water::SoilWaterModel,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)

    diff⁺.soil.water.K∇h = compute_surface_grad_bc(
        land.soil,
        bc.runoff_model,
        bc.precip_model,
        n̂,
        state⁻,
        diff⁻,
        aux⁻,
        t,
    )

end

"""
    soil_boundary_flux!(
        nf,
        bc::SurfaceDrivenHeatBoundaryConditions,
        heat::SoilHeatModel,
        land::LandModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Surface Driven BC method for `soil_boundary_flux!` for the
`SoilHeatModel`.
"""
function soil_boundary_flux!(
    nf,
    bc::SurfaceDrivenHeatBoundaryConditions,
    heat::SoilHeatModel,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)

    net_surface_flux = compute_net_radiative_energy_flux(bc.nswf_model, t)
    diff⁺.soil.heat.κ∇T = n̂ * (-net_surface_flux)
end


text:
struct Solution
    X::Array{Vector{Int}, 2}
end

function Solution(model::Model)
    @assert length(model.varsNotFixed) == 0 "Mais il se passe quoi là ..."

    X = Array{Vector{Int}, 2}(undef, model.g, model.w)

    for g in 1:model.g
        for w in 1:model.w
            X[g, w] = copy(model.X[g, w].lowerBound)
        end
    end

    return Solution(X)
end


text:
module TestTreiberStack

using Reagents
using Test

include("../../../examples/treiberstack.jl")

function test_1()
    stack = TreiberStack{Int}()
    push!(stack, 1)
    @test pop!(stack) == 1
end

function test_123()
    stack = TreiberStack{Int}()
    push!(stack, 1)
    push!(stack, 2)
    push!(stack, 3)
    @test pop!(stack) == 3
    @test pop!(stack) == 2
    @test pop!(stack) == 1
end

function test_spawn_pops()
    ns = [Threads.nthreads()]
    if Threads.nthreads() > 1
        push!(ns, Threads.nthreads() - 1)
    end
    @testset for nconsumers in ns
        test_spawn_pops(nconsumers)
    end
end

function test_spawn_pops(nconsumers)
    stack = TreiberStack{Int}()
    for _ in 1:nconsumers
        push!(stack, 0)
    end
    m = 100 * nconsumers
    for i in 1:m
        push!(stack, i)
    end

    vectors = [Int[] for _ in 1:nconsumers]
    Base.Experimental.@sync begin
        for xs in vectors
            Threads.@spawn begin
                while true
                    x = pop!(stack)
                    if x == 0
                        break
                    else
                        push!(xs, x)
                    end
                end
            end
        end
        for i in m+1:2m
            push!(stack, i)
        end
    end
    @debug "" map(length, vectors)
    popped = sort!(reduce(vcat, vectors))
    @test length(popped) >= m
    @test popped == 1:length(popped)
end

function test_paired_pop_pop()
    s1 = TreiberStack{Int}()
    s2 = TreiberStack{Int}()
    push!(s1, 1)
    push!(s2, 2)
    @test (trypopping(s1) & trypopping(s2))() == (Some(1), Some(2))
end

function test_paired_pop_push()
    s1 = TreiberStack{Int}()
    s2 = TreiberStack{Int}()
    push!(s1, 1)
    (trypopping(s1) ⨟ Reagents.Map(something) ⨟ pushing(s2))()
    @test pop!(s2) == 1
end

end  # module


text:
"""
    parent(taxon::NCBITaxon)

Returns the taxon from which the argument is descended.
"""
function Base.parent(tax::NCBITaxon)
    position = findfirst(isequal(tax.id), NCBITaxonomy.nodes_table.tax_id)
    parent_id = NCBITaxonomy.nodes_table.parent_tax_id[position]
    isnothing(parent_id) && return nothing
    return taxon(parent_id)
end
"""
    rank(taxon::NCBITaxon)

Returns the rank of a taxon.
"""
function rank(tax::NCBITaxon)
    position = findfirst(isequal(tax.id), NCBITaxonomy.nodes_table.tax_id)
    return NCBITaxonomy.nodes_table.rank[position]
end

"""
    lineage(tax::NCBITaxon; stop_at::NCBITaxon=ncbi"root")

Returns an array of `NCBITaxon` going up to the root of the taxonomy, or to the
optionally specified `stop_at` taxonomic node.
"""
function lineage(tax::NCBITaxon; stop_at::NCBITaxon=ncbi"root")
    full_lineage = [tax]
    while last(full_lineage) != stop_at
        push!(full_lineage, parent(last(full_lineage)))
    end
    return reverse(full_lineage)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PDEAssembler.jl/blob/master/LICENSE

using JuliaFEM
using JuliaFEM.Preprocess
add_elements! = JuliaFEM.add_elements!

#field_elements, boundary_elements = get_unit_square()

meshfile = Pkg.dir("FEMPlates", "examples", "dkt", "mesh.med")
mesh = aster_read_mesh(meshfile)

using FEMPlates

field_elements = create_elements(mesh, "PLATE")
update!(field_elements, "youngs modulus", 4880.0)
update!(field_elements, "poissons ratio", 1/3)
update!(field_elements, "thickness", 0.5)
update!(field_elements, "distributed load", 10.0)
plate = Problem(DKT, "test problem", 3)
add_elements!(plate, field_elements)

boundary_elements = create_elements(mesh, "BORDER")
update!(boundary_elements, "displacement 1", 0.0)
bc = Problem(Dirichlet, "fixed", 3, "displacement")
add_elements!(bc, boundary_elements)

analysis = Analysis(Linear)
add_problems!(analysis, [plate, bc])
run!(analysis)

X = plate("geometry", 0.0)
u = plate("displacement", 0.0)
N = length(X)
x = [X[i][1] for i=1:N]
y = [X[i][2] for i=1:N]
w = [u[i][1] for i=1:N]

using Plots
surface(x, y, w)
gui()
using PyPlot
elset_name = :PLATE
node_ids = sort(collect(keys(mesh.nodes)))
node_x = [mesh.nodes[id][1] for id in node_ids]
node_y = [mesh.nodes[id][2] for id in node_ids]
element_ids = mesh.element_sets[elset_name]
element_conn = [mesh.elements[id]-1 for id in element_ids]
triplot(node_x, node_y, element_conn, "-k")
axis("equal")
axis("off")


text:
# Initialization algorithms
#
#   Each algorithm is represented by a subtype of SeedingAlgorithm
#
#   Let alg be an instance of such an algorithm, then it should
#   support the following usage:
#
#       initseeds!(iseeds, alg, X)
#       initseeds_by_costs!(iseeds, alg, costs)
#
#   Here:
#       - iseeds:   a vector of resultant indexes of the chosen seeds
#       - alg:      the seeding algorithm instance
#       - X:        the data matrix, each column being a data point
#       - costs:    pre-computed pairwise cost matrix.
#
#   This function returns iseeds
#

"""
    SeedingAlgorithm

Base type for all seeding algorithms.

Each seeding algorithm should implement the two functions: [`initseeds!`](@ref)
and [`initseeds_by_costs!`](@ref).
"""
abstract type SeedingAlgorithm end

"""
    initseeds(alg::Union{SeedingAlgorithm, Symbol},
              X::AbstractMatrix, k::Integer) -> Vector{Int}

Select `k` seeds from a ``d×n`` data matrix `X` using the `alg`
algorithm.

`alg` could be either an instance of [`SeedingAlgorithm`](@ref) or a symbolic
name of the algorithm.

Returns the vector of `k` seed indices.
"""
initseeds(alg::SeedingAlgorithm, X::AbstractMatrix{<:Real}, k::Integer) =
    initseeds!(Vector{Int}(undef, k), alg, X)

"""
    initseeds_by_costs(alg::Union{SeedingAlgorithm, Symbol},
                       costs::AbstractMatrix, k::Integer) -> Vector{Int}

Select `k` seeds from the ``n×n`` `costs` matrix using algorithm `alg`.

Here, `costs[i, j]` is the cost of assigning points `i`` and ``j``
to the same cluster. One may, for example, use the squared Euclidean distance
between the points as the cost.

Returns the vector of `k` seed indices.
"""
initseeds_by_costs(alg::SeedingAlgorithm, costs::AbstractMatrix{<:Real}, k::Integer) =
    initseeds_by_costs!(Vector{Int}(undef, k), alg, costs)

seeding_algorithm(s::Symbol) =
    s == :rand ? RandSeedAlg() :
    s == :kmpp ? KmppAlg() :
    s == :kmcen ? KmCentralityAlg() :
    throw(ArgumentError("Unknown seeding algorithm $s"))

function check_seeding_args(n::Integer, k::Integer)
    k >= 1 || throw(ArgumentError("The number of seeds ($k) must be positive."))
    k <= n || throw(ArgumentError("Cannot select more seeds ($k) than data points ($n)."))
end

check_seeding_args(X::AbstractMatrix, iseeds::AbstractVector) =
    check_seeding_args(size(X, 2), length(iseeds))

initseeds(algname::Symbol, X::AbstractMatrix{<:Real}, k::Integer) =
    initseeds(seeding_algorithm(algname), X, k)::Vector{Int}

initseeds_by_costs(algname::Symbol, costs::AbstractMatrix{<:Real}, k::Integer) =
    initseeds_by_costs(seeding_algorithm(algname), costs, k)

# use specified vector of seeds
function initseeds(iseeds::AbstractVector{<:Integer}, X::AbstractMatrix{<:Real}, k::Integer)
    length(iseeds) == k ||
        throw(ArgumentError("The length of seeds vector ($(length(iseeds))) differs from the number of seeds requested ($k)"))
    check_seeding_args(X, iseeds)
    n = size(X, 2)
    # check that seed indices are fine
    for (i, seed) in enumerate(iseeds)
        (1 <= seed <= n) || throw(ArgumentError("Seed #$i refers to an incorrect data point ($seed)"))
    end
    # NOTE no duplicate checks are done, should we?
    convert(Vector{Int}, iseeds)
end
initseeds_by_costs(iseeds::AbstractVector{<:Integer}, costs::AbstractMatrix{<:Real}, k::Integer) =
    initseeds(iseeds, costs, k) # NOTE: passing costs as X, but should be fine since only size(X, 2) is used

function copyseeds!(S::AbstractMatrix{<:AbstractFloat},
                    X::AbstractMatrix{<:Real},
                    iseeds::AbstractVector)
    d, n = size(X)
    k = length(iseeds)
    size(S) == (d, k) ||
        throw(DimensionMismatch("Inconsistent seeds matrix dimensions: $((d, k)) expected, $(size(S)) given."))
    return copyto!(S, view(X, :, iseeds))
end

"""
    RandSeedAlg <: SeedingAlgorithm

Random seeding (`:rand`).

Chooses an arbitrary subset of ``k`` data points as cluster seeds.
"""
struct RandSeedAlg <: SeedingAlgorithm end

"""
    initseeds!(iseeds::AbstractVector{Int}, alg::SeedingAlgorithm,
               X::AbstractMatrix) -> iseeds

Initialize `iseeds` with the indices of cluster seeds for the `X` data matrix
using the `alg` seeding algorithm.
"""
function initseeds!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real})
    check_seeding_args(X, iseeds)
    sample!(1:size(X, 2), iseeds; replace=false)
end

"""
    initseeds_by_costs!(iseeds::AbstractVector{Int}, alg::SeedingAlgorithm,
                        costs::AbstractMatrix) -> iseeds

Initialize `iseeds` with the indices of cluster seeds for the `costs` matrix
using the `alg` seeding algorithm.

Here, `costs[i, j]` is the cost of assigning points ``i`` and ``j``
to the same cluster. One may, for example, use the squared Euclidean distance
between the points as the cost.
"""
function initseeds_by_costs!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real})
    check_seeding_args(X, iseeds)
    sample!(1:size(X,2), iseeds; replace=false)
end

"""
    KmppAlg <: SeedingAlgorithm

Kmeans++ seeding (`:kmpp`).

Chooses the seeds sequentially. The probability of a point to be chosen is
proportional to the minimum cost of assigning it to the existing seeds.

# References
> D. Arthur and S. Vassilvitskii (2007).
> *k-means++: the advantages of careful seeding.*
> 18th Annual ACM-SIAM symposium on Discrete algorithms, 2007.
"""
struct KmppAlg <: SeedingAlgorithm end

function initseeds!(iseeds::IntegerVector, alg::KmppAlg,
                    X::AbstractMatrix{<:Real},
                    metric::PreMetric = SqEuclidean())
    n = size(X, 2)
    k = length(iseeds)
    check_seeding_args(n, k)

    # randomly pick the first center
    p = rand(1:n)
    iseeds[1] = p

    if k > 1
        mincosts = colwise(metric, X, view(X, :, p))
        mincosts[p] = 0

        # pick remaining (with a chance proportional to mincosts)
        tmpcosts = zeros(n)
        for j = 2:k
            p = wsample(1:n, mincosts)
            iseeds[j] = p

            # update mincosts
            c = view(X, :, p)
            colwise!(tmpcosts, metric, X, view(X, :, p))
            updatemin!(mincosts, tmpcosts)
            mincosts[p] = 0
        end
    end

    return iseeds
end

function initseeds_by_costs!(iseeds::IntegerVector, alg::KmppAlg,
                             costs::AbstractMatrix{<:Real})
    n = size(costs, 1)
    k = length(iseeds)
    check_seeding_args(n, k)

    # randomly pick the first center
    p = rand(1:n)
    iseeds[1] = p

    if k > 1
        mincosts = costs[:, p]
        mincosts[p] = 0

        # pick remaining (with a chance proportional to mincosts)
        for j = 2:k
            p = wsample(1:n, mincosts)
            iseeds[j] = p

            # update mincosts
            updatemin!(mincosts, view(costs, :, p))
            mincosts[p] = 0
        end
    end

    return iseeds
end

"""
    KmCentralityAlg <: SeedingAlgorithm

K-medoids initialization based on centrality (`:kmcen`).

Choose the ``k`` points with the highest *centrality* as seeds.

# References
> Hae-Sang Park and Chi-Hyuck Jun.
> *A simple and fast algorithm for K-medoids clustering.*
> doi:10.1016/j.eswa.2008.01.039
"""
struct KmCentralityAlg <: SeedingAlgorithm end

function initseeds_by_costs!(iseeds::IntegerVector, alg::KmCentralityAlg,
                             costs::AbstractMatrix{<:Real})

    n = size(costs, 1)
    k = length(iseeds)
    check_seeding_args(n, k)

    # compute score for each item
    coefs = vec(sum(costs, dims=2))
    for i = 1:n
        @inbounds coefs[i] = inv(coefs[i])
    end

    # scores[j] = \sum_j costs[i,j] / (\sum_{j'} costs[i,j'])
    #           = costs[i,j] * coefs[i]
    scores = costs'coefs

    # lower score indicates better seeds
    sp = sortperm(scores)
    for i = 1:k
        @inbounds iseeds[i] = sp[i]
    end
    return iseeds
end

initseeds!(iseeds::IntegerVector, alg::KmCentralityAlg, X::AbstractMatrix{<:Real},
           metric::PreMetric = SqEuclidean()) =
    initseeds_by_costs!(iseeds, alg, Distances.pairwise(metric, X, dims=2))


text:
using JuMP, Compat.Test, JAMSDWriter

# Example testing basic use of NLExpr with JAMSDWriter.jl

if !isdef(:solver); solver = JAMSDWriter.JAMSDSolver(); end

@testset "example: jump_nlexpr" begin
    m = Model(solver=solver)

    n = 30
    l = -ones(n); l[1] = 0
    u = ones(n)
    @variable(m, l[i] <= x[i=1:n] <= u[i])
    @NLexpression(m, f1, x[1])
    @NLexpression(m, g, 1 + 9 * sum(x[j] ^ 2 for j = 2:n) / (n - 1))
    @NLexpression(m, h, 1 - (f1 / g) ^ 2)
    @NLexpression(m, f2, g * h)

    setvalue(x[1], 1)
    setvalue(x[2:n], zeros(n - 1))
    @NLobjective(m, :Min, f2)

    @NLconstraint(m, log(exp(x[2])) <= u[2])

    @test solve(m) == :Optimal
    @test isapprox(getvalue(x[1]), 1.0, atol=1e-5)
    @test isapprox(getvalue(x[2:end]), zeros(n - 1), atol=1e-5)
    @test isapprox(getvalue(f1), 1.0, atol=1e-5)
    @test isapprox(getvalue(f2), 0.0, atol=1e-5)
    @test isapprox(getvalue(g), 1.0, atol=1e-5)
    @test isapprox(getvalue(h), 0.0, atol=1e-5)
    @test isapprox(getobjectivevalue(m), 0.0, atol=1e-5)
end


text:
function minute(line::String)
    return parse(Int, line[16:17])
end

function foo()
    lines = sort(readlines(open("data/4.txt", "r")))
    x = Dict{Int, Array{Int, 1}}()
    guard = -1
    for line in lines
        if occursin("Guard", line)
            start = findfirst(isequal('#'), line)
            stop = findnext(isequal(' '), line, start)
            guard = parse(Int, line[start+1:stop-1])
            if !(guard in keys(x))
                x[guard] = zeros(Int, 60)
            end
        elseif occursin("falls", line)
            start = minute(line) + 1
            x[guard][start:end] .+= 1
        elseif occursin("wakes", line)
            start = minute(line) + 1
            x[guard][start:end] .-= 1
        end
    end

    sleepiest = -1
    time_asleep = -1
    for (guard, sleep) in x
        sum_sleep = sum(sleep)
        if sum_sleep > time_asleep
            sleepiest = guard
            time_asleep = sum_sleep
        end
    end

    max_minute = argmax(x[sleepiest]) - 1
    return sleepiest * max_minute
end

println(foo())


text:
"""
   initialized(petsclib)

Check if `petsclib` is initialized

# External Links
$(_doc_external("Sys/PetscInitialized"))
"""
function Initialized end
@for_petsc function initialized(::$UnionPetscLib)
    r_flag = Ref{PetscBool}()
    @chk ccall(
        (:PetscInitialized, $petsc_library),
        PetscErrorCode,
        (Ptr{PetscBool},),
        r_flag,
    )
    return r_flag[] == PETSC_TRUE
end

"""
   initialize([petsclib])

Initialized the `petsclib`, if no `petsclib` is given then all `PETSc.petsclibs`
will be initialized.

Additionally:
 - This will initialize MPI if it has not already been initialized.
 - It will disable the PETSc signal handler (via
   $(_petsc_link("Sys/PetscPopSignalHandler"))
 - Add an [`atexit`](https://docs.julialang.org/en/v1/base/base/#Base.atexit)
   hook to call [`PETSc.finalize`](@ref).

# External Links
$(_doc_external("Sys/PetscInitializeNoArguments"))
"""
function Initialize end

function initialize()
    map(initialize, petsclibs)
    return nothing
end

@for_petsc function initialize(::$UnionPetscLib)
    if !initialized($petsclib)
        MPI.Initialized() || MPI.Init()
        @chk ccall(
            (:PetscInitializeNoArguments, $petsc_library),
            PetscErrorCode,
            (),
        )

        # disable signal handler
        @chk ccall((:PetscPopSignalHandler, $petsc_library), PetscErrorCode, ())

        atexit(() -> finalize($petsclib))
    end
    return nothing
end

"""
   finalize(petsclib)

Finalize the `petsclib`, if no `petsclib` is given then all `PETSc.petsclibs`
will be finalized.

# External Links
$(_doc_external("Sys/PetscFinalize"))
"""
function finalize end

function finalize()
    map(finalize, petsclibs)
    return nothing
end

@for_petsc function finalize(::$UnionPetscLib)
    if !finalized($petsclib)
        @chk ccall((:PetscFinalize, $petsc_library), PetscErrorCode, ())
    end
    return nothing
end

"""
   finalized(petsclib)

Check if `petsclib` is finalized

# External Links
$(_doc_external("Sys/PetscFinalized"))
"""
function finalized end
@for_petsc function finalized(::$UnionPetscLib)
    r_flag = Ref{PetscBool}()
    @chk ccall(
        (:PetscFinalized, $petsc_library),
        PetscErrorCode,
        (Ptr{PetscBool},),
        r_flag,
    )
    return r_flag[] == PETSC_TRUE
end


text:
using Test
using ProductizeBayes

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

sm = VersionedModel("bernoulli", bernoulli_model)
rc = stan_sample(sm.stan_model; data, num_chains=2, seed=-1)

@test sm.hash == "5dce4091b38059e657f64ee7222f2b74778d852c051ad33bf0d646d1f0f66a05237510dd556214a343f5da8981faa5aef866a3cb109897f6074adbd38c8835bd"

if success(rc)
  sdf = read_summary(sm.stan_model)
  @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05

  (samples, parameters) = read_samples(sm.stan_model, :array;
    return_parameters=true)
  @test size(samples) == (1000, 1, 2)
  @test length(parameters) == 1

  (samples, parameters) = read_samples(sm.stan_model, :array;
    return_parameters=true, include_internals=true)
  @test size(samples) == (1000, 8, 2)
  @test length(parameters) == 8

  samples = read_samples(sm.stan_model, :array;
    include_internals=true)
  @test size(samples) == (1000, 8, 2)

  samples = read_samples(sm.stan_model, :array)
  @test size(samples) == (1000, 1, 2)
end


text:
module test_bukdu_routing

using Test

using Bukdu

struct WelcomeController <: ApplicationController
    conn::Conn
end

struct Route
    C::Type{<:ApplicationController}
    action
    param_types::Dict{Symbol,DataType}
    path_params::Vector{Pair{String,Any}}
    pipelines::Vector{Function}
end

index(::WelcomeController) = "hello"

C = WelcomeController
action = index

url = "/a/:b/:c"
segments = split(url, '/')
param_types = Dict{Symbol,DataType}(:b=>Int, :c=>Int)
(vals, path_params) = Bukdu.Routing.penetrate_segments(segments)
method = "GET"
pipelines = []
@eval route(::Val{Symbol($method)}, $(vals...)) = Route($C, $action, $param_types, Vector{Pair{String,Any}}($(path_params...)), $pipelines)

target_path = "/a/25/36"
segments = split(target_path, '/')
vals = [Val(Symbol(seg)) for seg in segments]
r = route(Val(Symbol(method)), vals...)

@test r.C === WelcomeController
@test r.action === index
@test r.path_params == Vector{Pair{String,Any}}(["b"=>"25", "c"=>"36"])

# https://discourse.julialang.org/t/write-a-rest-interface-like-flask/18538
function update(c::WelcomeController)
    c.path_params
end
get("/update/region/:region/site/:site_id/channel/:channel_id/", WelcomeController, update, :region=>String, :site_id=>Int, :channel_id=>Int)

Plug.Loggers.config[:path_pad] = 36
@test Router.call(get, "/update/region/west/site/1/channel/2").got == Assoc("region"=>"west", "site_id"=>1, "channel_id"=>2)

end # module test_bukdu_routing