text:
using BioMedQuery.Entrez
using SQLite
using XMLconvert

"""
pubmed_search_and_save(email, search_term, article_max::Int64,
db_path, verbose=false)

###Arguments

* email: valid email address (otherwise pubmed will block you)
* search_term : search string to submit to PubMed
e.g (asthma[MeSH Terms]) AND ("2001/01/29"[Date - Publication] : "2010"[Date - Publication])
see http://www.ncbi.nlm.nih.gov/pubmed/advanced for help constructing the string
* article_max : maximum number of articles to return. Defaults to 600,000
* db_path: path to output database
* verbose: if true, the NCBI xml response files are saved to current directory
"""
function pubmed_search_and_save(email, search_term, article_max,
    save_efetch_func, db_config, verbose=false)

    retstart = 0
    retmax = 10000  #e-utils only allows 10,000 at a time
    article_max = article_max

    if article_max < retmax
        retmax = article_max
    end

    db = nothing
    article_total = 0

    for rs=retstart:retmax:(article_max- 1)

        rm = rs + retmax
        if rm > article_max
            retmax = article_max - rs
        end

        println("Getting ", retmax, " articles, starting at index ", rs)

        #1. Formulate PubMed/MEDLINE search for articles between 2000 and 201
        #with obesity indicated as the major MeSH descriptor.
        println("------Searching Entrez--------")
        search_dic = Dict("db"=>"pubmed","term" => search_term,
        "retstart" => rs, "retmax"=>retmax, "tool" =>"BioJulia",
        "email" => email)
        esearch_response = esearch(search_dic)

        if verbose
            xmlASCII2file(esearch_response, "./esearch.xml")
        end

        #convert xml to dictionary
        esearch_dict = eparse(esearch_response)

        #2. Obtain PubMed/MEDLINE records (in MEDLINE or XML format) for
        #formulated search using NCBI E-Utilities.
        println("------Fetching Entrez--------")
        fetch_dic = Dict("db"=>"pubmed","tool" =>"BioJulia", "email" => email,
        "retmode" => "xml", "rettype"=>"null")
        #get the list of ids and perfom a fetch
        if !haskey(esearch_dict, "IdList")
            println("Error with esearch_dict:")
            println(esearch_dict)
            error("Response esearch_dict does not contain IdList")
        end

        ids = []
        for id_node in esearch_dict["IdList"][1]["Id"]
            push!(ids, id_node)
        end

        efetch_response = efetch(fetch_dic, ids)

        if verbose
            xmlASCII2file(efetch_response, "./efetch.xml")
        end

        efetch_dict = eparse(efetch_response)

        #save the results of an entrez fetch to a sqlite database
        println("------Saving to database--------")
        db = save_efetch_func(efetch_dict, db_config, verbose)

        #after the first pass - make sure the database is not deleted
        db_config[:overwrite] = false

        article_total+=length(ids)

        if (length(ids) < retmax)
            break
        end

    end

    println("Finished searching, total number of articles: ", article_total)
    return db
end


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
matplotlib.use("agg") # or try "macosx"
using SparseArrays

# geometry setup in domain [0,1]^2
solid_left = 0.45
solid_right = 0.55
solid_top = 0.5
solid_bottom = 0.52

chip_left = 0.48
chip_right = 0.52
chip_top = 0.5
chip_bottom = 0.505

delta=1e-5
delta2=1e-5

k_mold = 0.014531
k_chip = 2.60475
k_air = 0.64357
nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

filename = "mesh/CHT_2D.stl"
file_format = "stl"
mesh = Mesh(filename, file_format = file_format, degree=2)
mesh = Mesh(mesh.nodes ./ 0.030500000342726707, mesh.elems, -1, 2)
h = 0.1

# m = 100
# n = 100
# h = 1/n
# mesh = Mesh(m, n, h, degree=2)

nnode = mesh.nnode
nedge = mesh.nedge
ndof = mesh.ndof
nelem = mesh.nelem
ngauss = get_ngauss(mesh)

NT = 1    # number of iterations for Newton's method, 8 is good for m=400

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
chip_fem_idx = Array{Int64, 1}([])
# chip_fvm_idx = Array{Int64, 1}([])
chip_fem_top_idx = Array{Int64, 1}([])
bd = Array{Int64, 1}([])
fvm_bd = Array{Int64, 1}([])

for j = 1:nnode
    nodex, nodey = mesh.nodes[j, 1], mesh.nodes[j, 2]
    if nodex >= solid_left-delta2 && nodex <= solid_right+delta2 && nodey >= solid_top-delta2 && nodey <= solid_bottom+delta2
        # print(i, j)
        global solid_fem_idx = [solid_fem_idx; j]
        if nodex >= chip_left-delta2 && nodex <= chip_right+delta2 && nodey >= chip_top-delta2 && nodey <= chip_bottom+delta2
            global chip_fem_idx = [chip_fem_idx; j]
            if nodey <= chip_top+delta2
                global chip_fem_top_idx = [chip_fem_top_idx; j]
            end
        end
    end
    if abs(nodex-0.0) <= delta || abs(nodex-1.0) <= delta || abs(nodey-0.0) <= delta || abs(nodey-1.0) <= delta
        global bd = [bd; j]
    end
end

# fix chip_fem_top_idx
if size(chip_fem_top_idx, 1) == 0
    chip_fem_top_idx = chip_fem_idx
end

for j = 1:nedge
    edgex, edgey = (mesh.nodes[mesh.edges[j, 1], :] .+ mesh.nodes[mesh.edges[j, 2], :]) ./ 2
    if edgex >= solid_left-delta2 && edgex <= solid_right+delta2 && edgey >= solid_top-delta2 && edgey <= solid_bottom+delta2
        # print(i, j)
        global solid_fem_idx = [solid_fem_idx; nnode + j]
        if edgex >= chip_left-delta2 && edgex <= chip_right+delta2 && edgey >= chip_top-delta2 && edgey <= chip_bottom+delta2
            global chip_fem_idx = [chip_fem_idx; nnode + j]
            if edgex >= chip_left-delta2 && edgex <= chip_right+delta2 && edgey >= chip_top-delta2 && edgey <= chip_top+delta2
                global chip_fem_top_idx = [chip_fem_top_idx; nnode + j]
            end
        end
    end
    if abs(edgex-0.0) <= delta || abs(edgex-1.0) <= delta || abs(edgey-0.0) <= delta || abs(edgey-1.0) <= delta
        global bd = [bd; nnode + j]
    end
end

# bd_vec = zeros(ndof)
# for i in 1:length(bd)
#     bd_vec[ bd[i] ] = 1.0
# end
# figure();visualize_scalar_on_fem_points(bd_vec[1:nnode],mesh,with_mesh=true);savefig("bd.pdf")

# solid_vec = zeros(ndof)
# for i in 1:length(solid_fem_idx)
#     solid_vec[ solid_fem_idx[i] ] = 1.0
# end
# for i in 1:length(chip_fem_idx)
#     solid_vec[ chip_fem_idx[i] ] = 2.0
# end
# figure();visualize_scalar_on_fem_points(solid_vec[1:nnode],mesh,with_mesh=true);savefig("solid.pdf")

gaussxy = gauss_nodes(mesh)

for i in 1:nelem
    gaussx, gaussy = mean(gaussxy[6*i-5: 6*i, 1]), mean(gaussxy[6*i-5: 6*i, 2])
    # print("\nx=",gaussx, " y=",gaussy)
    if gaussx >= solid_left-delta2 && gaussx <= solid_right+delta2 && gaussy >= solid_top-delta2 && gaussy <= solid_bottom+delta2
        global solid_fvm_idx = [solid_fvm_idx; i]
    end
    if abs(gaussy - 0.0) <= h/2 + delta && ( abs(gaussx - 0.0) <= h + delta || abs(gaussx - 1.0) <= h + delta )
        global fvm_bd = [fvm_bd; i]
    end
end

k_fem = k_air * constant(ones(ndof))
k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
k_fem = scatter_update(k_fem, chip_fem_idx, k_chip * ones(length(chip_fem_idx)))
kgauss = dof_to_gauss_points(k_fem, mesh)  # try change to include k on edges

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
# heat_source_fem[chip_fem_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)  # try change to include k on edges

B = constant(compute_interaction_matrix(mesh)) # function not exist
# TODO: add unstructured mesh version of compute_interaction_matrix

# compute F
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, mesh))
# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = dof_to_gauss_points(k, mesh)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, mesh)

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
     bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

function compute_residual(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]
    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    interaction = compute_interaction_term(p, mesh) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, mesh)
    f2 = compute_fem_source_term1(vgauss.*uy, mesh)
    f3 = -interaction[1:ndof]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, mesh)
    g2 = compute_fem_source_term1(vgauss.*vy, mesh)
    g3 = -interaction[ndof+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = dof_to_gauss_points(T, mesh)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, mesh)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, mesh) * T - heat_source
    R = [F;G;H0;T0]
    # return R  # original
    # return[F;G;H0;constant(zeros(ndof))]  # fix T=0
    # return[F;G;constant(zeros(nelem));constant(zeros(ndof))]  # fix p=0 and T=0
end

function compute_jacobian(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]

    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    M1 = constant(compute_fem_mass_matrix1(ux, mesh))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, mesh))

    N1 = constant(compute_fem_mass_matrix1(vy, mesh))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, mesh))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, mesh))

    gradT = eval_grad_on_gauss_pts1(T, mesh)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, mesh))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, mesh))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = constant(compute_fem_mass_matrix1(-buoyance_coef * constant(ones(ndof)), mesh))
    T_mat = [SparseTensor(spzeros(ndof, ndof)); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(nelem, ndof))]

    # N = 2 * ndof + nelem
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros(ndof, nelem)) M]]

    # return [J[1:2*ndof+nelem, 1:2*ndof+nelem]  SparseTensor(spzeros(2*ndof+nelem, ndof))
    # SparseTensor(spzeros(ndof,2*ndof+nelem)) spdiag(ones(ndof))] # fix T=0

    # return [J[1:2*ndof, 1:2*ndof]  SparseTensor(spzeros(2*ndof, nelem+ndof))
    # SparseTensor(spzeros(nelem+ndof,2*ndof)) spdiag(ones(nelem+ndof))]  # fix p=0 and T=0
end

function solve_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, residual = impose_Dirichlet_boundary_conditions(J, residual, bd, zeros(length(bd)))

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(nelem+3*ndof))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, nelem+3*ndof))

sess = Session(); init(sess)
output = run_profile(sess, S)
save_profile("test.json")
error()
# J = run(sess, J); rank(J)
# J1 = run(sess, J1); rank(J1)
# J0 = run(sess, J0); rank(J0)

matwrite("xz_chip_unstructured_data.mat",
    Dict(
        "V"=>output[end, :]
    ))

# TODO: fix plot
u_out, v_out, p_out, T_out = output[NT+1,1:nnode], output[NT+1,ndof+1:ndof+nnode],
                             output[NT+1,2*ndof+1:2*ndof+nelem],output[NT+1,2*ndof+nelem+1:2*ndof+nelem+nnode]

figure(figsize=(10,10))
subplot(221)
title("u velocity")
visualize_scalar_on_fem_points(u_out .* u_std, mesh);#gca().invert_yaxis()
subplot(222)
title("v velocity")
visualize_scalar_on_fem_points(v_out .* u_std, mesh);#gca().invert_yaxis()
subplot(223)
visualize_scalar_on_fvm_points(p_out .* p_std, mesh);#gca().invert_yaxis()
title("pressure")
subplot(224)
title("temperature")
visualize_scalar_on_fem_points(T_out.* T_infty .+ T_infty, mesh);#gca().invert_yaxis()
tight_layout()
savefig("forward_solution_unstructured.pdf")

print("Solution range:",
    "\n [u velocity] \t min:", minimum(u_out .* u_std), ",\t max:", maximum(u_out .* u_std),
    "\n [v velocity] \t min:", minimum(v_out .* u_std), ",\t max:", maximum(v_out .* u_std),
    "\n [pressure]   \t min:", minimum(p_out .* p_std), ",\t max:", maximum(p_out .* p_std),
    "\n [temperature]\t min:", minimum(T_out.* T_infty .+ T_infty), ",\t\t\t max:", maximum(T_out.* T_infty .+ T_infty))

ext: jl


text:
module MultiAgentPOMDPs

using POMDPs

abstract type JointMDP{S, A<:AbstractVector} <: MDP{S, A} end

"""
    function n_agents(m::JointMDP)
Return the number of agents in the (PO)MDP
"""
function n_agents end

"""
    function agent_actions(m::JointMDP, idx::Int64, s::S) where S
Returns the discrete actions for the given agent index.

!!! note
    This will be called a LOT so it should not allocate each time....
"""
function agent_actions end

"""
    function agent_states(m::JointMDP, idx::Int64)
Returns the discrete states for the given agent index
"""
function agent_states end

"""
    function agent_actionindex(m::JointMDP, idx::Int64, a::A) where A
Returns the integer index of action `a` for agent `idx`. Used for discrete models only.
"""
function agent_actionindex end

"""
    function agent_actionindex(m::JointMDP, idx::Int64, s::S) where S
Returns the integer index of state `s` for agent `idx`. Used for discrete models only.
"""
function agent_stateindex end

"""
    function coordination_graph(m::JointMDP)
    function coordination_graph(m::JointMDP, s::S) where S
Returns the LightGraphs.SimpleGraph (or any appropriate structure) for the coordination graph.
"""
function coordination_graph end

export JointMDP, n_agents, agent_actions, agent_states, agent_actionindex, agent_stateindex, agent_reward, coordination_graph
end


text:
###############################################################################
## Least squares with nuclear norm regularization
#    min_x  0.5 Σᴹ (⟨Aᵢ, X⟩-yᵢ)² + λ ||X||_*
struct TracenormPb{Tr, Tm, Txman} <: CompositeProblem
    A::Vector{Matrix{Float64}}
    y::Vector{Float64}
    n1::Int
    n2::Int                         # matrix variable, dimension (n1,n2)
    m::Int
    regularizer::Tr
    x0::Txman
    M_x0::Tm
end

problem_dimension(pb::TracenormPb) = pb.n1*pb.n2

# 0th order
function f(pb::TracenormPb, X::AbstractMatrix)
    return 0.5 * sum( (dot(pb.A[i], X) - pb.y[i])^2 for i in 1:pb.m)
end
## TODO: improve this to use point structure.
f(pb::TracenormPb, X::SVDMPoint) = f(pb, X.U*Diagonal(X.S)*X.Vt)

# function f(pb::TracenormPb, x::AbstractVector)
#     X = reshape(pb, x)
#     return 0.5 * sum( (dot(pb.A[i], X) - pb.y[i])^2 for i in 1:pb.m)
# end

# 1st order
∇f(pb::TracenormPb, X::AbstractMatrix) = sum( pb.A[i] * (dot(pb.A[i], X) - pb.y[i]) for i in 1:pb.m)
function ∇f!(pb::TracenormPb, ∇f::AbstractMatrix, X::AbstractMatrix)
    ∇f .= 0.0
    for i in 1:pb.m
        ∇f .+= pb.A[i] * (dot(pb.A[i], X) - pb.y[i])
    end
    return ∇f
end

## TODO: improve this to account for point structure
∇f(pb::TracenormPb, X::SVDMPoint) = ∇f(pb, X.U*Diagonal(X.S)*X.Vt)
∇f!(pb::TracenormPb, res, X::SVDMPoint) = ∇f!(pb, res, X.U*Diagonal(X.S)*X.Vt)

# function ∇f(pb::TracenormPb, x::AbstractVector)
#     X = reshape(pb, x)
#     return vec(sum( pb.A[i] * (dot(pb.A[i], X) - pb.y[i]) for i in 1:pb.m))
# end

# 2nd order
∇²f_h(pb::TracenormPb, X::AbstractMatrix, H::AbstractMatrix) = sum( pb.A[i] * dot(pb.A[i], H) for i in 1:pb.m)
function ∇²f_h!(pb::TracenormPb, ∇²f::AbstractMatrix, X::AbstractMatrix, H::AbstractMatrix)
    ∇²f .= 0
    for i in 1:pb.m
        ∇²f .+= pb.A[i] * dot(pb.A[i], H)
    end
    return ∇²f
end

## TODO: improve this to account for point structure
∇²f_h(pb::TracenormPb, X::SVDMPoint, H) = ∇²f_h(pb, X.U*Diagonal(X.S)*X.Vt, X.U * H.M * X.Vt + X.U * H.Vt + H.U * X.Vt)
∇²f_h!(pb::TracenormPb, res, X::SVDMPoint, H) = ∇²f_h!(pb, res, X.U*Diagonal(X.S)*X.Vt, X.U * H.M * X.Vt + X.U * H.Vt + H.U * X.Vt)

# function ∇²f_h(pb::TracenormPb, x::AbstractVector, h::AbstractVector)
#     X = reshape(pb, x)
#     H = reshape(pb, h)
#     return vec(sum( pb.A[i] * dot(pb.A[i], H) for i in 1:pb.m))
# end

# conditioning
# get_gradlips(pb::TracenormPb) = opnorm(pb.A)^2
# get_μ_cvx(pb::TracenormPb) = (svdvals(pb.A)[end])^2

## Instance generation
function get_tracenorm_MLE(;n1=10, n2=12, m=5, sparsity=0.6, seed=1234, δ=0.01)
    # A is drawn from the standard normal distribution
    # b = Ax0 + e where e is taken from the normal distribution with standard deviation 0.001.
    # We set λ1 so that the original sparsity is ultimately recovered.
    @assert 0 ≤ sparsity ≤ 1

    A = Vector{Matrix{Float64}}(undef, m)
    for i in 1:m
        Random.seed!(seed+i)
        A[i] = rand(Normal(), n1, n2)
    end

    ## Generating structured signal
    nsingvals = min(n1, n2)
    rank = Int(ceil((1-sparsity) * nsingvals))
    optstructure = FixedRankMatrices(n1, n2, rank)

    Random.seed!(seed-1)
    x0 = project(optstructure, rand(Normal(), n1, n2))

    ## Noised measurements
    Random.seed!(seed)
    e = rand(Normal(0, δ^2), m)

    y = [dot(A[i], x0) for i in 1:m] .+ e

    return TracenormPb(A, y, n1, n2, m, regularizer_lnuclear(δ), x0, optstructure)
end


text:
cd(@__DIR__); include("setups/path16.jl");
gr(dpi = 250)

plot(size = (800, 250), framestyle = :none, xlim = [-10.5, N+13], ylim = [-5, 1])
path_pc_plot!(W, X; c = :yellow)

k = 1
path_pc_plot!(W, X .- [12 4]; c = :yellow,
        ind = GP_primal.inds[GP_primal.rs[k, 2]:GP_primal.rs[k+1, 2] - 1, 2])

k = 2
path_pc_plot!(W, X .- [-12 4]; c = :yellow,
        ind = GP_primal.inds[GP_primal.rs[k, 2]:GP_primal.rs[k+1, 2] - 1, 2])

plot!([6, -1], [-0.8, -3.3], line = :solid, c = :black, lw = 0.65)
plot!([11, 18], [-0.8, -3.3], line = :solid, c = :black, lw = 0.65)
annotate!(mean(X[:, 1]), 0.5,
    text(latexstring("V= \\{ \\delta_{1}, ..., \\delta_{16} \\}"), :bottom, 11))
annotate!(mean(X[:, 1]) - 12, -4 - 0.7,
    text(latexstring("V_{0}= \\{ \\delta_{1}, \\delta_{3}, \\delta_{5}, \\delta_{7}, \\delta_{10}, \\delta_{12}, \\delta_{14}, \\delta_{16} \\}"), :top, 11))
annotate!(mean(X[:, 1]) + 12, -4 - 0.7,
    text(latexstring("V_{1}= \\{ \\delta_{2}, \\delta_{4}, \\delta_{6}, \\delta_{8}, \\delta_{9}, \\delta_{11}, \\delta_{13}, \\delta_{15} \\}"), :top, 11))

plt = current()
savefig(plt, "../figs/Path16_1lev_PC.png")
display(plt)


text:
#JulieTest config file

include(Pkg.dir("JulieTest/src/watch.jl"))

# Configure the watcher

watch = ["src", "test"]
interval = 0.25 # in seconds

# Configure the test runner

# reporter = "dot"
reporter = "spec"

# start JulieTest
julieTest()


text:
module TestMult_GetDataframe

using Mimi
using NamedArrays
using Test

#####################################
#  LARGER MULTIREGIONAL TEST (2/3)  #
#####################################

#GROSS ECONOMY COMPONENT
@defcomp grosseconomy begin
    regions = Index()                           #Note that a regional index is defined here

    YGROSS  = Variable(index=[time, regions])   #Gross output
    K       = Variable(index=[time, regions])   #Capital
    l       = Parameter(index=[time, regions])  #Labor
    tfp     = Parameter(index=[time, regions])  #Total factor productivity
    s       = Parameter(index=[time, regions])  #Savings rate
    depk    = Parameter(index=[regions])        #Depreciation rate on capital - Note that it only has a region index
    k0      = Parameter(index=[regions])        #Initial level of capital
    share   = Parameter()                       #Capital share

    function run_timestep(p, v, d, t)
        #Define an equation for K
        for r in d.regions
            if is_first(t)
                v.K[t,r] = p.k0[r]
            else
                v.K[t,r] = (1 - p.depk[r])^5 * v.K[t-1,r] + v.YGROSS[t-1,r] * p.s[t-1,r] * 5
            end
        end

        #Define an equation for YGROSS
        for r in d.regions
            v.YGROSS[t,r] = p.tfp[t,r] * v.K[t,r]^p.share * p.l[t,r]^(1-p.share)
        end
    end
end

#EMISSIONS COMPONENT
@defcomp emissions begin
    regions     = Index()                           #The regions index must be specified for each component

    E           = Variable(index=[time, regions])   #Total greenhouse gas emissions
    E_Global    = Variable(index=[time])            #Global emissions (sum of regional emissions)
    sigma       = Parameter(index=[time, regions])  #Emissions output ratio
    YGROSS      = Parameter(index=[time, regions])  #Gross output - Note that YGROSS is now a parameter

    function run_timestep(p, v, d, t)
        #Define an eqation for E
        for r in d.regions
            v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
        end

        #Define an equation for E_Global
        for r in d.regions
            v.E_Global[t] = sum(v.E[t,:])
        end
    end
end

#DEFINE ALL THE PARAMETERS
l = Array{Float64}(undef,20,3)
for t in 1:20
    l[t,1] = (1. + 0.015)^t *2000
    l[t,2] = (1. + 0.02)^t * 1250
    l[t,3] = (1. + 0.03)^t * 1700
end

tfp = Array{Float64}(undef, 20,3)
for t in 1:20
    tfp[t,1] = (1 + 0.06)^t * 3.2
    tfp[t,2] = (1 + 0.03)^t * 1.8
    tfp[t,3] = (1 + 0.05)^t * 2.5
end

s = Array{Float64}(undef, 20,3)
for t in 1:20
    s[t,1] = 0.21
    s[t,2] = 0.15
    s[t,3] = 0.28
end

depk = [0.11, 0.135 ,0.15]
k0   = [50.5, 22., 33.5]

sigma = Array{Float64}(undef, 20,3)
for t in 1:20
    sigma[t,1] = (1. - 0.05)^t * 0.58
    sigma[t,2] = (1. - 0.04)^t * 0.5
    sigma[t,3] = (1. - 0.045)^t * 0.6
end

#FUNCTION TO RUN MY MODEL
function run_my_model()

    my_model = Model()

    set_dimension!(my_model, :time, collect(2015:5:2110))
    set_dimension!(my_model, :regions, ["Region1", "Region2", "Region3"])  #Note that the regions of your model must be specified here

    add_comp!(my_model, grosseconomy)
    add_comp!(my_model, emissions)

    set_param!(my_model, :grosseconomy, :l, l)
    set_param!(my_model, :grosseconomy, :tfp, tfp)
    set_param!(my_model, :grosseconomy, :s, s)
    set_param!(my_model, :grosseconomy, :depk,depk)
    set_param!(my_model, :grosseconomy, :k0, k0)
    set_param!(my_model, :grosseconomy, :share, 0.3)

    #set parameters for emissions component
    set_param!(my_model, :emissions, :sigma, sigma2)
    connect_param!(my_model, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    run(my_model)
    return(my_model)
end

# DEFINE ALL THE PARAMETERS using NAMEDARRAYS
region_labels = ["Region1", "Region2", "Region3"]
time_labels = collect(2015:5:2110)
num_time_labels = length(time_labels)

l2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:num_time_labels
    l2[t, 1] = (1. + 0.015)^t *2000
    l2[t, 2] = (1. + 0.02)^t * 1250
    l2[t, 3] = (1. + 0.03)^t * 1700
end

tfp2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:num_time_labels
    tfp2[t, 1] = (1 + 0.06)^t * 3.2
    tfp2[t, 2] = (1 + 0.03)^t * 1.8
    tfp2[t, 3] = (1 + 0.05)^t * 2.5
end

s2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:num_time_labels
    s2[t, 1] = 0.21
    s2[t, 2] = 0.15
    s2[t, 3] = 0.28
end

depk2 = NamedArray([0.11, 0.135 ,0.15], (region_labels,), (:regions,))
k02   = NamedArray([50.5, 22., 33.5], (region_labels,), (:regions,))

sigma2 = NamedArray(Array{Float64}(undef, 20, 3), (time_labels, region_labels), (:time, :regions))
for t in 1:num_time_labels
    sigma2[t, 1] = (1. - 0.05)^t * 0.58
    sigma2[t, 2] = (1. - 0.04)^t * 0.5
    sigma2[t, 3] = (1. - 0.045)^t * 0.6
end

run1 = run_my_model()

expanded_get_dataframe = getdataframe(run1, :grosseconomy =>(:YGROSS, :K))
curr_get_dataframe = getdataframe(run1, :grosseconomy, :YGROSS)
@test(expanded_get_dataframe[!, 3] == curr_get_dataframe[!, 3])

expanded_get_dataframe = getdataframe(run1, :grosseconomy => :YGROSS, :emissions => :E)
curr_get_dataframe = getdataframe(run1, :emissions, :E)
@test(expanded_get_dataframe[!, 4] == curr_get_dataframe[!, 3])

#################
## MORE TESTS ###
#################
my_model = Model()
set_dimension!(my_model, :time, collect(2015:5:2110))

#Testing that you cannot add two components of the same name
@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

par = collect(2015:5:2110)
add_comp!(my_model, testcomp1)

set_param!(my_model, :testcomp1, :par1, par)
run(my_model)

#Regular getdataframe
dataframe = getdataframe(my_model, :testcomp1 => :var1)
@test(dataframe[!, 2] == par)

#Test trying to getdataframe from component that does not exist
@test_throws ErrorException getdataframe(my_model, :testcomp1 => :var2)

end #module


text:
const TPQ = 960

const THIRTY_SECOND = (TPQ//8) // TPQ
const SIXTEENTH = (TPQ//4) // TPQ
const TRIPLET = (TPQ//3) // TPQ
const EIGHTH = (TPQ//2) // TPQ
const QUARTER = TPQ // TPQ
const HALF = (TPQ*2) // TPQ
const WHOLE = (TPQ*4) // TPQ


text:
using VortexDistributions

# simple test field
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi1 = Torus(copy(psi0), x, y)
psi2 = Torus(copy(psi0), x, y)

# WEIRD BEHAVIOUR WHEN VORTICES ARE THIS CLOSE TO EACH OTHER, NEED TO INVESTIGATE
# vp = PointVortex(.1,.3,1)
# vn = PointVortex(-.2,.7,-1)

xp = 1.1; yp = .3;
xn = -2.23; yn = .7;
vp = PointVortex(xp, yp, 1)
vn = PointVortex(xn, yn, -1)
dip = Dipole(vp,vn)
@test typeof(dip) <: VortexGroup

sp = ScalarVortex([vp, vn])
vortex!(psi1, sp)
vfound1 = findvortices(psi1)

periodic_dipole!(psi2, sp)
vfound2 = findvortices(psi2, periodic=true)

vpf1 = vfound1[2];
vnf1 = vfound1[1];
vpf2 = vfound2[2];
vnf2 = vfound2[1];

dx = Δ(psi1.x)
@test isapprox(vpf1.xv, xp, atol=dx/4)
@test isapprox(vpf1.yv, yp, atol=dx/4)
@test vpf1.qv == 1

@test isapprox(vnf1.xv, xn, atol=dx/4)
@test isapprox(vnf1.yv, yn, atol=dx/4)
@test vnf1.qv == -1

@test isapprox(vpf2.xv, xp, atol=dx/4)
@test isapprox(vpf2.yv, yp, atol=dx/4)
@test vpf2.qv == 1

@test isapprox(vnf2.xv, xn, atol=dx/4)
@test isapprox(vnf2.yv, yn, atol=dx/4)
@test vnf2.qv == -1


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Edward J. Xu
# 2019.5.28
########################################################################################################################

function doOptim(vecModelSub, modMas, vecConsRef, vecConsConvex, vecLambda, epsilon)
    ## 4,  Begin Iteration
    println("#### Begin Iteration ###########################################################") ########################
    extremePoints = [[]]
    extremePointForSub = [-1]
    done = false
    iter = 1
    while !done
        (vec_pi, vec_kappa, obj_master) = solveMas(modMas, vecConsRef, vecConsConvex)
        done = true
        ## Print the result
        # println("vec_pi = $(vec_pi)\n",
        #         "vec_kappa = $(vec_kappa)")
        # println("--------------------------------------------------------------------------------\n",
        #         "$(iter)-th iteration. obj_master = $(obj_master).\n",
        #         "--------------------------------------------------------------------------------")
        ## Column Generation
        costReduceBest = -1
        for k = 1: length(vecModelSub)
            (costReduce, vec_x_result) = solveSub(vecModelSub[k], vec_pi, vec_kappa[k])
            # println("Reduced cost of $(k)-th sub model is $costReduce.")
            if costReduce > costReduceBest
                costReduceBest = costReduce
            end
            # remember to look for negative reduced cost if we are minimizing.
            if costReduce > epsilon
                lambdaNew = addColToMaster(modMas, vecModelSub[k], vec_x_result, vecConsRef, vecConsConvex[k])
                push!(vecLambda, lambdaNew)
                push!(extremePoints, vec_x_result)
                push!(extremePointForSub, k)
                done = false
            end
        end
        ##
        iter += 1
        # println("best reduced cost is $costReduceBest")
    end
    ## 5,  Print the Result
    println("################################################################################\n", ######################
            "Optimization Done After $(iter-1) Iterations.\n",
            "obj_master = $(getobjectivevalue(modMas))\n",
            "################################################################################")
    vecLambdaResult = getvalue(vecLambda)
    return (vecLambdaResult, extremePointForSub, extremePoints)
end

function getVecStrNameVar(numSub, numXInSub)
    vecStrNameVar = Vector{String}(undef, numSub * numXInSub)
    idx = 1
    for i = 1: numSub
        for j = 1: numXInSub
            vecStrNameVar[idx] = "x_{$(i),$(j)}"
            idx += 1
        end
    end
    return vecStrNameVar
end

function printResult(vecLambdaResult, extremePointForSub, extremePoints, vecStrNameVar, matIndexSub)
    # compute values of original variables
    origVarValSub = []
    numSub = size(matIndexSub)[1]
    for s = 1:numSub
        push!(origVarValSub, zeros(length(matIndexSub[s,:])))
    end
    for p = 1: length(vecLambdaResult)
        if vecLambdaResult[p] > 0.000001
            println("lambda_$p = ", vecLambdaResult[p], ", sub = $(extremePointForSub[p]), \n",
                    "extreme point = $(extremePoints[p]).")
            origVarValSub[extremePointForSub[p]] += vecLambdaResult[p] * extremePoints[p]
        end
    end
    for s = 1:numSub
        #println("var val for sub problem $s: $(origVarValSub[s])")
        for t = 1: length(origVarValSub[s])
            if abs(origVarValSub[s][t]) > 0.000001
                println("$(vecStrNameVar[matIndexSub[s, t]]) = $(origVarValSub[s][t])")
            end
        end
    end
end


text:
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function pcarray(f::F, ::Type{A}, sz) where {F,A}
    a = f(A(undef, sz))
    fill!(a, zero(eltype(a)))
    return first(a)
end
function pcm(a1, a2; fillvalue = zero(eltype(a1)))
    v = mosaic(a1, a2; fillvalue=fillvalue)
    return first(v)
end
function pcmv(a; fillvalue = zero(eltype(a)))
    v = mosaicview(a; fillvalue=fillvalue)
    return first(v)
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    eltypes = (N0f8, N0f16, Float32, Float64)        # eltypes of parametric colors
    pctypes = (Gray, RGB, AGray, GrayA, ARGB, RGBA)  # parametric colors
    cctypes = (Gray24, AGray32, RGB24, ARGB32)       # non-parametric colors
    dims  = (1, 2, 3, 4)
    szs   = ((2,), (2, 2), (2, 2, 2), (2, 2, 2, 2))

    for T in eltypes
        @warnpcfail precompile(clamp01, (T,))
        @warnpcfail precompile(clamp01nan, (T,))
        @warnpcfail precompile(scaleminmax, (T, T))
        @warnpcfail precompile(scalesigned, (T,))
        @warnpcfail precompile(scalesigned, (T,T,T))
        for C in pctypes
            @warnpcfail precompile(clamp01, (C{T},))
            @warnpcfail precompile(clamp01nan, (C{T},))
            @warnpcfail precompile(colorsigned, (C{T},C{T}))
        end
    end
    for C in cctypes
        @warnpcfail precompile(clamp01, (C,))
        @warnpcfail precompile(clamp01nan, (C,))
        @warnpcfail precompile(colorsigned, (C,C))
    end
    # For the arrays, it's better to make them and exercise them so we get the getindex/setindex!
    # methods precompiled too.
    for sz in szs
        for T in eltypes
            T <: FixedPoint || continue
            pcarray(rawview, Array{T,length(sz)}, sz)
        end
        pcarray(normedview, Array{UInt8,length(sz)}, sz)
        pcarray(a->normedview(N0f8, a),  Array{UInt8, length(sz)}, sz)
        pcarray(a->normedview(N0f16, a), Array{UInt16,length(sz)}, sz)
    end
    for sz in szs
        for T in eltypes
            for C in pctypes
                nc = sizeof(C{T}) ÷ sizeof(T)
                if nc == 1
                    pcarray(a->colorview(C,    a), Array{T,length(sz)}, sz)
                    pcarray(a->colorview(C{T}, a), Array{T,length(sz)}, sz)
                else
                    pcarray(a->colorview(C,    a), Array{T,length(sz)+1}, (nc, sz...))
                    pcarray(a->colorview(C{T}, a), Array{T,length(sz)+1}, (nc, sz...))
                end
                pcarray(channelview, Array{C{T},length(sz)}, sz)
                if T<:FixedPoint
                    R = FixedPointNumbers.rawtype(T)
                    if nc == 1
                        pcarray(a->colorview(C,    normedview(T, a)), Array{R,length(sz)}, sz)
                        pcarray(a->colorview(C{T}, normedview(T, a)), Array{R,length(sz)}, sz)
                    else
                        pcarray(a->colorview(C,    normedview(T, a)), Array{R,length(sz)+1}, (nc, sz...))
                        pcarray(a->colorview(C{T}, normedview(T, a)), Array{R,length(sz)+1}, (nc, sz...))
                    end
                end
            end
        end
        T, C = Bool, Gray
        pcarray(a->colorview(C,    a), Array{T,length(sz)}, sz)
        pcarray(a->colorview(C{T}, a), Array{T,length(sz)}, sz)
        pcarray(channelview, Array{C{T},length(sz)}, sz)
    end
    for T in eltypes
        a = zeros(T, (2, 2))
        a3 = zeros(T, (2, 2, 2))
        pcm(a, a)
        pcmv(a3)
        for C in (Gray, RGB, GrayA, RGBA)
            a = zeros(C{T}, (2, 2))
            a3 = zeros(C{T}, (2, 2, 3))
            pcm(a, a)
            pcmv(a3)
            if C === RGB
                pcm(a, a; fillvalue=zero(Gray{T}))
            elseif C === RGBA
                pcm(a, a; fillvalue=zero(GrayA{T}))
            end
        end
    end
    pcm(zeros(Float32, 2, 2), zeros(Float64, 2, 2))  # heterogeneous
end


text:
export load_result
export replacement_summary

function load_result(save_file::String)
    isfile(save_file) || error("File $save_file does not exist.")

    result = load(save_file)["result"]

    return result
end

function replacement_summary(
    result::HLAModelResult;
    fdr::Bool = false,
    fisher_p::Bool = false
)
    sf = stanfit(result)
    si = stan_input(result)
    data = hla_data(result)
    posterior = extract(sf)

    alleles = hla_alleles(result)
    rs = replacements(result)

    df = DataFrame(
        allele = HLAAllele[],
        position = Int[],
        replacement = Char[],
        posterior_p = Float64[],
        log_odds_lower_95 = Float64[],
        log_odds_upper_95 = Float64[],
        fisher_p = Float64[],
        in_predicted_epitope = Int[],
        n_replacement_total = Int[],
        n_replacement_with_hla = Int[]
    )

    R = si["R"]
    D = si["D"]
    Z = si["Z"]

    if fisher_p
        fisher = Escape.run(Escape.FisherTest(), data; fdr = fdr)
    end

	for r in 1:length(rs)
        fisher_idx = fisher_p ? findfirst(x -> x == rs[r], fisher.replacements) : 0

        for d in 1:D
            beta_hla = posterior["beta_hla.$r.$d"]
            lower, upper = quantile(beta_hla, (0.025, 0.975))
            p_value = fisher_p ? fisher.p_values[fisher_idx][alleles[d]] : 0
            n_total = fisher_p ? sum(fisher.counts[fisher_idx][alleles[d]][1, :]) : 0
            n_with_hla = fisher_p ? fisher.counts[fisher_idx][alleles[d]][1, 1] : 0

            z = Z[r, d]

            push!(df,
                [
                    alleles[d],
                    position(rs[r]),
                    replacement(rs[r]),
                    mean(beta_hla .> 0),
                    lower,
                    upper,
                    p_value,
                    z,
                    n_total,
                    n_with_hla
                ]
            )
        end
    end

    sort!(df, :posterior_p, rev = true)

    return df
end

function reduce_size!(result::HLAModelResult)
    reduce_size!(result.sf)

    return result
end

function diagnostics(result::HLAModelResult)
    sf = stanfit(result)

    return sf.diagnostics
end

hla_model(result::HLAModelResult) = result.model
hla_data(result::HLAModelResult) = result.data
stanfit(result::HLAModelResult) = result.sf
stan_input(result::HLAModelResult) = stanfit(result).data
hla_alleles(result::HLAModelResult) = result.alleles
replacements(result::HLAModelResult) = result.replacements


text:
@testset "NumField/Elem" begin
  k = wildanger_field(3, 13)[1]
  c = cyclotomic_field(3)[1]
  q = quadratic_field(17)[1]

  @testset for K in [k, c, q]
    R, x = PolynomialRing(K, 3)

    f = sum(x)^2
    g = sum(x) * sum(t+2 for t = x)
    h = sum(x)

    @test gcd(f,g) == h

    @test isone(gcd(f+1,g))

  end
end


text:
"""
```
 set_regime_val!(p::Parameter{S},
    i::Int, v::S; override_bounds::Bool = false) where S <: Real
 set_regime_val!(p::Parameter{S},
    model_regime::Int, v::S, d::AbstractDict{Int, Int}; override_bounds::Bool = false) where S <: Real
```

sets the value in regime `i` of `p` to be `v`. By default, we enforce
the bounds that are currently in `p`, but the bounds can be ignoerd by
setting `override_bounds = true`.

The second method allows the user to pass a dictionary to permit the case where
there may be differences between the regimes of a regime-switching model and
the regimes for the parameters. For example, aside from regime-switching in parameters,
the model may also include other forms of regime-switching. To allow
estimation of regime-switching parameters in such a model, the dictionary `d`
maps each "model" regime to a "parameter" regime. In this way,
the second method specifies which "parameter" regime should be used at a given
"model" regime.
"""
function set_regime_val!(p::Parameter{S},
                         i::Int, v::S; override_bounds::Bool = false) where S <: Real
    if !haskey(p.regimes, :value)
        p.regimes[:value] = OrderedDict{Int,S}()
    end

    # First check if fixed and enforce valuebounds so
    # `set_regime_val!` mirrors `parameter(p::Parameter, newvalue::S)` functionality
    if haskey(p.regimes, :fixed) && haskey(p.regimes[:fixed], i) ? regime_fixed(p, i) : false
        if haskey(p.regimes[:value], i)
            return p.regimes[:value][i]
        else # If it doesn't exist yet, then we want to set the value for this regime
            p.regimes[:value][i] = v
            p.regimes[:valuebounds][i] = (v, v)
        end
    elseif (haskey(p.regimes, :valuebounds) && haskey(p.regimes[:valuebounds], i) ?
            ((regime_valuebounds(p, i)[1] <= v <= regime_valuebounds(p, i)[2]) || override_bounds) : false)
        p.regimes[:value][i] = v
    elseif (p.valuebounds[1] <= v <= p.valuebounds[2]) || override_bounds
        p.regimes[:value][i] = v
    elseif p.fixed
        # When parameters are initialized as non-regime-switching and fixed,
        # the valuebounds is automatically set to (p.value, p.value).
        # Unless valuebounds are set to regime-switching, it is not possible
        # to add extra regimes without using `override_bounds = true` (unless using the same value).
        # Note that the rest of ModelConstructors.jl assumes that if `p.fixed = true`
        # and haskey(p.regimes, :fixed) is false, then all regimes (if any) are also fixed.
        throw(ParamBoundsError("Parameter $(p.key) is fixed. Regimes cannot be added unless " *
                               "keyword `override_bounds` is set to `true`."))
    else
        throw(ParamBoundsError("New value of $(string(p.key)) ($(v)) is out of bounds ($(p.valuebounds))"))
    end
    return v
end

function set_regime_val!(p::Parameter{S}, model_regime::Int,
                         v::S, d::AbstractDict{Int, Int}; override_bounds::Bool = false) where S <: Real
    return set_regime_val!(p, d[model_regime], v; override_bounds = override_bounds)
end

"""
```
regime_val(p::Parameter{S}, i::Int) where S <: Real
regime_val(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
```

returns the value of `p` in regime `i` for the first method
and the value of `p` in regime `d[model_regime` for the second.
"""
function regime_val(p::Parameter{S}, i::Int) where S <: Real
    if !haskey(p.regimes, :value) || !haskey(p.regimes[:value], i)
        @error "regime_val(), Input Error: No regime $(i)"
    end
    return p.regimes[:value][i]
end

regime_val(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real = regime_val(p, d[model_regime])

"""
```
set_regime_prior!(p::Parameter{S}, i::Int, v)
set_regime_prior!(p::Parameter{S}, model_regime::Int, v, d::AbstractDict{Int, Int})
```

sets the prior in regime `i` of `p` to be `v`. The type of `v`
can be a `NullablePriorUnivariate`, `NullablePriorMultivariate`,
`ContinuousUnivariateDistribution`, or `ContinuousMultivariateDistribution'.

The second method allows the user to pass a dictionary to permit the case where
there may be differences between the regimes of a regime-switching model and
the regimes for the parameters. For example, aside from regime-switching in parameters,
the model may also include other forms of regime-switching. To allow
estimation of regime-switching parameters in such a model, the dictionary `d`
maps each "model" regime to a "parameter" regime. In this way,
the second method specifies which "parameter" regime should be used at a given
"model" regime.
"""
function set_regime_prior!(p::Parameter, i::Int, v::S) where {S <: Union{NullablePriorUnivariate, NullablePriorMultivariate}}
    if !haskey(p.regimes, :prior)
        p.regimes[:prior] = OrderedDict{Int, Union{NullablePriorUnivariate, NullablePriorMultivariate}}()
    end
    p.regimes[:prior][i] = v

    return v
end

function set_regime_prior!(p::Parameter, i::Int, v::S) where S <: ContinuousUnivariateDistribution
    return set_regime_prior!(p, i, NullablePriorUnivariate(v))
end

function set_regime_prior!(p::Parameter, i::Int, v::S) where S <: ContinuousMultivariateDistribution
    return set_regime_prior!(p, i, NullablePriorMultivariate(v))
end

function set_regime_prior!(p::Parameter, model_regime::Int, v::S,
                           d::AbstractDict{Int, Int}) where {S <: Union{NullablePriorUnivariate, NullablePriorMultivariate}}
    return set_regime_prior!(p, d[model_regime], v)
end

function set_regime_prior!(p::Parameter, model_regime::Int,
                           v::S, d::AbstractDict{Int, Int}) where S <: ContinuousUnivariateDistribution
    return set_regime_prior!(p, model_regime, NullablePriorUnivariate(v), d)
end

function set_regime_prior!(p::Parameter, model_regime::Int, v::S,
                           d::AbstractDict{Int, Int}) where S <: ContinuousMultivariateDistribution
    return set_regime_prior!(p, model_regime, NullablePriorMultivariate(v), d)
end

"""
```
regime_prior(p::Parameter{S}, i::Int) where S <: Real
regime_prior(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
```

returns the prior of `p` in regime `i` for the first method
and the prior of `p` in regime `d[model_regime]` for the second.
"""
function regime_prior(p::Parameter{S}, i::Int) where S <: Real
    if !haskey(p.regimes, :prior) || !haskey(p.regimes[:prior], i)
        @error "regime_prior(), Input Error: No regime $(i)"
    end
    return p.regimes[:prior][i]
end

regime_prior(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real = regime_prior(p, d[model_regime])

"""
```
set_regime_fixed!(p::Parameter{S}, i::Int, v::S; update_valuebounds::Interval = (NaN, NaN))
```

sets whether `p` is fixed in regime `i` of `p`. Set update_valuebounds to true to
set the valuebounds to match the fixed value.

The second method allows the user to pass a dictionary to permit the case where
there may be differences between the regimes of a regime-switching model and
the regimes for the parameters. For example, aside from regime-switching in parameters,
the model may also include other forms of regime-switching. To allow
estimation of regime-switching parameters in such a model, the dictionary `d`
maps each "model" regime to a "parameter" regime. In this way,
the second method specifies which "parameter" regime should be used at a given
"model" regime.
"""
function set_regime_fixed!(p::Parameter{S1}, i::Int, v::S; update_valuebounds::Interval{S1} = (NaN, NaN)) where {S <: Bool, S1 <: Real}
    if !haskey(p.regimes, :fixed)
        p.regimes[:fixed] = OrderedDict{Int, Bool}()
    end
    p.regimes[:fixed][i] = v

    if !haskey(p.regimes, :valuebounds)
        p.regimes[:valuebounds] = OrderedDict{Int, typeof(p.value)}()
    end

    if isnan(update_valuebounds[1]) || isnan(update_valuebounds[2])
        # Do nothing unless `v` is true (fixed in regime `i`) or
        # !haskey(p.regimes[:valuebounds], i) (no value for regime `i` for valuebounds dict)
        if v
            if (haskey(p.regimes, :value) ? haskey(p.regimes[:value], i) : false)
                # Regime is fixed, so valuebounds should be set to (value, value), if it exists
                p.regimes[:valuebounds][i] = (regime_val(p, i), regime_val(p, i))
            else
                error("Regime $(i) for parameter $(p.key) does not have a value. Set the value " *
                      "using `set_regime_val!` before making it a fixed value with `set_regime_fixed!`.")
            end
        elseif !haskey(p.regimes[:valuebounds], i)
            # If valuebounds is nonexistent in regime i, initialize to p.valuebounds
            p.regimes[:valuebounds][i] = p.valuebounds
        end
    else # Otherwise update the valuebounds. Note that there is no check that p.regimes[:value][i] lies inside the valuebounds
        p.regimes[:valuebounds][i] = update_valuebounds
    end

    return v
end

function set_regime_fixed!(p::Parameter, model_regime::Int, v::S,
                           d::AbstractDict{Int, Int}; update_valuebounds::Interval{S1} = (NaN, NaN)) where {S <: Bool, S1 <: Real}
    set_regime_fixed!(p, d[model_regime], v; update_valuebounds = update_valuebounds)
end

"""
```
regime_fixed(p::Parameter{S}, i::Int) where S <: Real
regime_fixed(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
```

returns whether `p` is fixed in regime `i` for the first method
and whether true `p` is fixed in regime `d[model_regime]` for the second method.
"""
function regime_fixed(p::Parameter{S}, i::Int) where S <: Real
    if !haskey(p.regimes, :fixed) || !haskey(p.regimes[:fixed], i)
        @error "regime_fixed(), Input Error: No regime $(i)"
    end
    return p.regimes[:fixed][i]
end

regime_fixed(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real = regime_fixed(p, d[model_regime])

"""
```
set_regime_valuebounds!(p::Parameter{S}, i::Int, v::S)
```

sets valuebounds for `p` in regime `i` to `v`.

The second method allows the user to pass a dictionary to permit the case where
there may be differences between the regimes of a regime-switching model and
the regimes for the parameters. For example, aside from regime-switching in parameters,
the model may also include other forms of regime-switching. To allow
estimation of regime-switching parameters in such a model, the dictionary `d`
maps each "model" regime to a "parameter" regime. In this way,
the second method specifies which "parameter" regime should be used at a given
"model" regime.
"""
function set_regime_valuebounds!(p::Parameter, i::Int, v::Interval{S}) where {S <: Real}
    if !haskey(p.regimes, :valuebounds)
        p.regimes[:valuebounds] = OrderedDict{Int, typeof(p.value)}()
    end
    p.regimes[:valuebounds][i] = v

    return v
end

function set_regime_valuebounds!(p::Parameter, model_regime::Int, v::Interval{S},
                                 d::AbstractDict{Int, Int}) where {S <: Real}
    set_regime_valuebounds!(p, d[model_regime], v)
end

"""
```
regime_valuebounds(p::Parameter{S}, i::Int) where S <: Real
regime_valuebounds(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
```

returns the `valuebounds` of `p`  in regime `i` for the first method
and the `valuebounds` of `p` in regime `d[model_regime]` for the second method.
"""
function regime_valuebounds(p::Parameter{S}, i::Int) where S <: Real
    if !haskey(p.regimes, :valuebounds) || !haskey(p.regimes[:valuebounds], i)
        @error "regime_valuebounds(), Input Error: No regime $(i)"
    end
    return p.regimes[:valuebounds][i]
end

regime_valuebounds(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real = regime_valuebounds(p, d[model_regime])

"""
```
toggle_regime!(p::Parameter{S}, i::Int) where S <: Real
toggle_regime!(pvec::ParameterVector{S}, i::Int) where S <: Real
toggle_regime!(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
toggle_regime!(pvec::ParameterVector{S}, model_regime::Int, d::AbstractDict{Symbol, <: AbstractDict{Int, Int}}) where S <: Real
toggle_regime!(pvec::ParameterVector{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real
```

changes the fields of `p` to regime `i`.

For example, if

```
p.regimes[:value] = OrderedDict{Int, Any}(1 => 1, 2 => 3)
```

then `toggle_regime!(p, 1)` will cause `p.value = 1` and `toggle_regime!(p, 2)`
will cause `p.value = 3`.

The third method allows the user to pass a dictionary to permit the case where
there may be differences between the regimes of a regime-switching model and
the regimes for the parameters. For example, aside from regime-switching in parameters,
the model may also include other forms of regime-switching. To allow
estimation of regime-switching parameters in such a model, the dictionary `d`
maps each "model" regime to a "parameter" regime. In this way,
the second method specifies which "parameter" regime should be used at a given
"model" regime.

The fourth method extends the third to a ParameterVector, with the possibility
that each parameter may have different mappings to the model regimes. Each key
of `d` corresponds to the key of a parameter, and each value of `d` is
the mapping for model regimes to the parameter regimes of `p.key`.
The fifth method is similar to the fourth but assumes
any regime-switching parameter has the same mapping from model regimes
to parameter regimes, hence the use of a common dictionary.
"""
function toggle_regime!(p::Parameter{S}, i::Int) where S <: Real
    if !isempty(p.regimes)
        for field in [:value, :valuebounds, :transform_parameterization,
                      :transform, :prior, :fixed]
            if haskey(p.regimes, field) && haskey(p.regimes[field], i)
                if field == :value
                    p.value = p.regimes[field][i]
                elseif field == :valuebounds
                    p.valuebounds = p.regimes[field][i]
                elseif field == :transform_parameterization
                    p.transform_parameterization = p.regimes[field][i]
                elseif field == :transform
                    p.transform = p.regimes[field][i]
                elseif field == :prior
                    p.prior = p.regimes[field][i]
                elseif field == :fixed
                    p.fixed = p.regimes[field][i]
                end
            elseif haskey(p.regimes, field) && !haskey(p.regimes[field], i)
                error("Regime $i for field $field not found")
            end
        end
    end
end

function toggle_regime!(pvec::ParameterVector{S}, i::Int) where {S <: Real}
    for p in pvec
        toggle_regime!(p, i)
    end
end

toggle_regime!(p::Parameter{S}, model_regime::Int, d::AbstractDict{Int, Int}) where S <: Real = toggle_regime!(p, d[model_regime])

function toggle_regime!(pvec::ParameterVector{S}, model_regime::Int, d::AbstractDict{Symbol, <: AbstractDict{Int, Int}}) where S <: Real
    for p in pvec
        if haskey(d, p.key)
            toggle_regime!(p, model_regime, d[p.key])
        end
    end
end

toggle_regime!(pvec::ParameterVector{S}, model_regime::Int, d::AbstractDict{Int, Int}) where {S <: Real} = toggle_regime!(pvec, d[model_regime])

"""
```
get_values(pvec::ParameterVector{S}; regime_switching::Bool = true) where {S <: Real}
```

constructs a vector of the underlying values in a `ParameterVector`, including
if there are regime-switching values.
"""
function get_values(pvec::ParameterVector{S}; regime_switching::Bool = true) where {S <: Real}

    if regime_switching  # Check if regime switching occurs
        np_reg = n_parameters_regime_switching(pvec)
        np     = length(pvec)
        if np == np_reg # No regime-switching
            vals = [x.value for x in pvec]
        else
            vals = Vector{S}(undef, np_reg)

            # An initial pass to find regime 1 values
            for i in 1:np
                if isempty(pvec[i].regimes)
                    vals[i] = pvec[i].value
                elseif haskey(pvec[i].regimes, :value)
                    vals[i] = pvec[i].regimes[:value][1]
                end
            end

            # A second loop to add in the extra regimes
            ct = np # a counter to assist us
            for i in 1:np
                if !isempty(pvec[i].regimes)
                    if haskey(pvec[i].regimes, :value)
                        for j in 1:length(pvec[i].regimes[:value])
                            if j > 1
                                ct += 1
                                vals[ct] = pvec[i].regimes[:value][j]
                            end
                        end
                    end
                end
            end
        end
    else # Regime switching doesn't occur, so just directly map
        vals = [x.value for x in pvec]
    end

    return vals
end

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


text:
using MAT
using PyPlot
using JOLI
using GenSPGL
using SeisJOLI
using Arpack
using LinearAlgebra

###load the weighted algorithm function
include("../script/Weighted_LR.jl")

###load the true data
datafile = MAT.matopen("../data/Full.mat")
data_full = read(datafile,"Full_Data")
datafull = copy(data_full[:,1:355,1:355])
nff,ns,nr = size(datafull)

###loaded the jittered index for subsampling
sub_index = MAT.matopen("../data/ind.mat")
ind = read(sub_index,"ind")
close(sub_index)

###define operators subsampling and joMH(source-receiver => midpoint-offset)
RM = joKron(joRestriction(ns,dropdims(trunc.(Int,ind),dims=1);DDT = Complex{Float64}), joDirac(nr;DDT = Complex{Float64}))
MH = joSRtoCMO(nr,ns)

###define rank information if rank=r_sub(convertional weighted method)
rank = 85;
r_sub = 25;

###define the final result to save
Final_result = zeros(Complex{Float64},nr,ns,nff)
Final_SNR = zeros(Float64,nff,1)

###weighted LR to interpolate the data
for sliceNum = 31:304  #(round([7Hz,75Hz]/df+1=[30,304]))

    if sliceNum == 31
	# load unweighted result as prior information
	pre_int = MAT.matopen(string("../data/Idx_", string(sliceNum-1),".mat"))
	pre_L = read(pre_int,"L1")
	pre_R = read(pre_int,"R1")
	close(pre_int)

	# using limited-subspace weighted method
    	SNR, Result, L1, R1 = Weighted_LR(pre_L, pre_R, datafull[sliceNum,:,:], RM, MH, rank, r_sub)
    	println("Num=",sliceNum)
    else
	# load prior information(neighbor low frequency results)
	pre_int = MAT.matopen(string("../julia_result/Idx_", string(sliceNum-1),".mat"))
	pre_L = read(pre_int,"L1")
	pre_R = read(pre_int,"R1")
	close(pre_int)

	# using limited-subspace weighted method
	SNR, Result, L1, R1 = Weighted_LR(pre_L, pre_R, datafull[sliceNum,:,:], RM, MH, rank, r_sub)
    	println("Num=",sliceNum)
    end
    Final_SNR[sliceNum,1] = SNR;
    Final_result[:,:,sliceNum] = Result;

    # save each frequency result
    file = matopen(string("../julia_result/Idx_",string(sliceNum),".mat"), "w")
    write(file,string("Final_result",string(sliceNum)),Result)
    write(file,"SNR",SNR)
    write(file,"L1",L1)
    write(file,"R1",R1)
    close(file);
end

### save all limited-subspace weighted result
file = matopen("../julia_result/recursive_data.mat", "w")
write(file,"Final_result",Final_result)
write(file,"Final_SNR",Final_SNR)
close(file)


text: Dict{String, Any}("best-day-ever" => Dict{String, Any}("value" => "1987-07-05T17:45:00Z", "type" => "datetime"), "numtheory" => Dict{String, Any}("perfection" => Any[Dict{String, Any}("value" => "6", "type" => "integer"), Dict{String, Any}("value" => "28", "type" => "integer"), Dict{String, Any}("value" => "496", "type" => "integer")], "boring" => Dict{String, Any}("value" => "false", "type" => "bool")))


text:
import TOML
build_artifact = joinpath(@__DIR__, "artifact")
isfile(build_artifact) && rm(build_artifact)
project = TOML.parsefile(Base.active_project())
@assert get(project["deps"], "JSON", nothing) === nothing
manifest = Base.get_deps(TOML.parsefile(joinpath(dirname(Base.active_project()), "Manifest.toml")))
json = manifest["JSON"][1]
@assert json["uuid"] == "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@assert json["version"] == "0.19.0"
touch(build_artifact)


text:
#=ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля

   РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так:
   нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке,
   следующий - за исключением двух последних клеток на Востоке, и т.д.
=#
function mountain!(r::Robot)
    to_Sud = moves!(r,Sud)
    to_West = moves!(r,West)
#по стандарту в левый нижний угол

   steps=putmarkers!(r,Ost) #нижняя линия в маркерах(робот в нижнем правом углу)
   back_moves!(r,West) #робот возвращаяется в левый угол
    while isborder(r,Nord)==false
       move!(r,Nord)
       putmarker!(r)
       steps-=1
       mark_part_line!(r,Ost,steps)
       back_moves!(r,West)
    end
    if isborder(r,Nord)==true
        back_moves!(r,Sud)
    end
#сначала маркеруется самая нижняя строчка,потом по while-end маркируются последующие строчки(на 1 клетку меньше)
#и под конец,если сверху перегородка-вернуться вниз

   moves!(r,Nord,to_Sud)
   moves!(r,Ost,to_West)
#возвращаемся в исходную точку

end
function back_moves!(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
    end
end
function mark_part_line!(r::Robot,side::HorizonSide,steps::Int)
    for _ in 1:steps
        move!(r,side)
        putmarker!(r)
    end
end
#маркеруем линию;на 1 меньше нижней

function putmarkers!(r::Robot, side::HorizonSide)
    steps=0
    putmarker!(r)
    while isborder(r,side)==false
        move!(r,side)
        steps+=1
        putmarker!(r)
    end
    return steps
end
#функция при которой вся линии помечается маркерами(включая первоначальную точку)и выводится кол-во ячеек в линии

function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end
#функция,при которой считаются шаги до границы

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end
#в исходную точку

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
using GR
#lines
x = LinRange(0, 10, 100)
y1 = sin.(x)
y2 = cos.(x)

figure()
subplot(2, 2, 1)
plot(x, y1, "r")

subplot(2, 2, 2)
plot(x, y2, "b")

subplot(2, 2, (3,4))
plot(x, y1, "r", x, y2, "b", labels = ["sin", "cos"])

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
module util

using Switch

export @enum

include("util.string.jl")
include("util.json.jl")
include("util.zip.jl")
include("util.uuid.jl")

macro enum(T,syms...)
    blk = quote
        immutable $(esc(T))
            n::Int32
            $(esc(T))(n::Integer) = new(n)
        end
        Base.show(io::IO, x::$(esc(T))) = print(io, $syms[x.n+1])
        Base.show(io::IO, x::Type{$(esc(T))}) = print(io, $(string("enum ", T, ' ', '(', join(syms, ", "), ')')))
    end
    for (i,sym) in enumerate(syms)
        push!(blk.args, :(const $(esc(sym)) = $(esc(T))($(i-1))))
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end

export @switch, @case, @default, break

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# Autogenerated wrapper script for ALPS_jll for aarch64-apple-darwin
export libalps

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ALPS")
JLLWrappers.@declare_library_product(libalps, "@rpath/libAlps.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libalps,
        "lib/libAlps.0.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#================================

    All kinds of functions related to pruning

=================================#

"""
    get_next_prune_constraint(com::CS.CoM, constraint_idxs_vec)

Check which function will be called for pruning next. This is based on `constraint_idxs_vec`. The constraint with the lowest
value is chosen and if two have the same lowest value the first one is chosen.
Return the best value and the constraint index. Return a constraint index of 0 if there is no constraint with a less than maximal value
"""
function get_next_prune_constraint(com::CS.CoM, constraint_idxs_vec)
    best_ci = 0
    best_open = typemax(Int)
    for ci in 1:length(constraint_idxs_vec)
        if constraint_idxs_vec[ci] < best_open
            best_ci = ci
            best_open = constraint_idxs_vec[ci]
        end
    end
    return best_open, best_ci
end

"""
    open_possibilities(search_space, indices)

Return the sum of possible values for the given list of indices. Does not count 1 if the value is fixed
"""
function open_possibilities(search_space, indices)
    open = 0
    for vi in indices
        if !isfixed(search_space[vi])
            open += nvalues(search_space[vi])
        end
    end
    return open
end

"""
    prune!(com::CS.CoM; pre_backtrack=false, all=false, only_once=false, initial_check=false)

Prune based on changes by initial solve or backtracking. The information for it is stored in each variable.
There are several parameters:
`pre_backtrack` when set to true `com.info.in_backtrack_calls` is incremented
`all` instead of only looking at changes each constraint is check at least once (if there are open values)
`only_once` Just run on the changed constraints or `all` once instead of repeatedly until nothing can be pruned
`initial_check` Checks on `all` constraints and also checks if variables are set for the whole constraint
    whether the constraint is fulfilled or the problem is infeasible.
Return whether it's still feasible
"""
function prune!(
    com::CS.CoM;
    pre_backtrack = false,
    all = false,
    only_once = false,
    initial_check = false,
)
    feasible = true
    N = typemax(Int)
    search_space = com.search_space
    prev_var_length = zeros(Int, length(search_space))
    constraint_idxs_vec = fill(N, length(com.constraints))
    # get all constraints which need to be called (only once)
    for var in search_space
        new_var_length = num_changes(var, com.c_step_nr)
        if new_var_length > 0 || all || initial_check
            prev_var_length[var.idx] = new_var_length
            for ci in com.subscription[var.idx]
                com.constraints[ci].is_deactivated && continue
                inner_constraint = com.constraints[ci]
                constraint_idxs_vec[inner_constraint.idx] =
                    open_possibilities(search_space, inner_constraint.indices)
            end
        end
    end

    # while we haven't called every constraint
    while true
        open_pos, ci = get_next_prune_constraint(com, constraint_idxs_vec)
        # no open values => don't need to call again
        if open_pos == 0 && !initial_check
            constraint_idxs_vec[ci] = N
            continue
        end
        # checked all
        if open_pos == N
            break
        end
        constraint_idxs_vec[ci] = N
        constraint = com.constraints[ci]

        changed!(com, constraint, constraint.fct, constraint.set)
        feasible =
            prune_constraint!(com, constraint, constraint.fct, constraint.set; logs = false)
        if !pre_backtrack
            com.info.in_backtrack_calls += 1
        else
            com.info.pre_backtrack_calls += 1
        end
        if !feasible
            break
        end

        # if we changed another variable increase the level of the constraints to call them later
        for vidx in constraint.indices
            var = search_space[vidx]
            new_var_length = num_changes(var, com.c_step_nr)
            if new_var_length > prev_var_length[var.idx]
                prev_var_length[var.idx] = new_var_length
                for ci in com.subscription[var.idx]
                    com.constraints[ci].is_deactivated && continue
                    # don't call the same constraint again.
                    # Each constraint should prune as much as possible
                    if ci != constraint.idx
                        inner_constraint = com.constraints[ci]
                        # if initial check or don't add constraints => update only those which already have open possibilities
                        if (only_once || initial_check) &&
                            constraint_idxs_vec[inner_constraint.idx] == N
                            continue
                        end
                        constraint_idxs_vec[inner_constraint.idx] =
                        open_possibilities(search_space, inner_constraint.indices)
                    end
                end
            end
        end
    end
    return feasible
end

"""
    restore_prune!(com::CS.CoM, prune_steps)

Prune the search space based on a list of backtracking indices `prune_steps`.
"""
function restore_prune!(com::CS.CoM, prune_steps)
    search_space = com.search_space
    for backtrack_idx in prune_steps
        step_nr = com.backtrack_vec[backtrack_idx].step_nr
        for var in search_space
            for change in view_changes(var, step_nr)
                fct_symbol = change[1]
                val = change[2]
                if fct_symbol == :fix
                    fix!(com, var, val; changes = false, check_feasibility = false)
                elseif fct_symbol == :rm
                    rm!(com, var, val; changes = false, check_feasibility = false)
                elseif fct_symbol == :remove_above
                    remove_above!(com, var, val; changes = false, check_feasibility = false)
                elseif fct_symbol == :remove_below
                    remove_below!(com, var, val; changes = false, check_feasibility = false)
                else
                    throw(ErrorException("There is no pruning function for $fct_symbol"))
                end
            end
        end
        com.c_backtrack_idx = backtrack_idx
    end
    call_restore_pruning!(com, prune_steps)
end

"""
    single_reverse_pruning!(search_space, vidx::Int, prune_int::Int, prune_fix::Int)

Reverse a single variable using `prune_int` (number of value removals) and `prune_fix` (new last_ptr if not 0).
"""
function single_reverse_pruning!(search_space, vidx::Int, prune_int::Int, prune_fix::Int)
    if prune_int > 0
        var = search_space[vidx]
        l_ptr = max(1, var.last_ptr)

        new_l_ptr = var.last_ptr + prune_int
        min_val, max_val = extrema(var.values[l_ptr:new_l_ptr])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.last_ptr = new_l_ptr
    end
    if prune_fix > 0
        var = search_space[vidx]
        var.last_ptr = prune_fix

        min_val, max_val = extrema(var.values[1:prune_fix])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.first_ptr = 1
    end
end

"""
    reverse_pruning!(com, backtrack_idx)

Reverse the changes made by a specific backtrack object
"""
function reverse_pruning!(com::CS.CoM, backtrack_idx::Int)
    com.c_backtrack_idx = backtrack_idx
    step_nr = com.backtrack_vec[backtrack_idx].step_nr
    search_space = com.search_space
    for var in search_space
        v_idx = var.idx
        changes = var.changes.changes
        ch_indices = var.changes.indices
        for change_id in ch_indices[step_nr+1]-1:-1:ch_indices[step_nr]
            single_reverse_pruning!(search_space, v_idx, changes[change_id][4], changes[change_id][3])
        end
    end
    subscriptions = com.subscription
    constraints = com.constraints
    for var in search_space
        num_changes(var, step_nr) == 0 && continue
        var.idx > length(subscriptions) && continue
        @inbounds for ci in subscriptions[var.idx]
            constraint = constraints[ci]
            single_reverse_pruning_constraint!(
                com,
                constraint,
                constraint.fct,
                constraint.set,
                var,
                backtrack_idx,
            )
        end
    end
    for constraint in constraints
        reverse_pruning_constraint!(
            com,
            constraint,
            constraint.fct,
            constraint.set,
            backtrack_idx,
        )
    end
    com.c_backtrack_idx = com.backtrack_vec[backtrack_idx].parent_idx
end


text:
using CSV, Distributions, DataFrames

# Load once and make global
arests = CSV.read(joinpath(@__DIR__, "../../../data/other/arestimates.csv"), DataFrame)

function calc_temp(p, v, d, tt, annual_year)
    # for every year, do the same calculations, but then with the new annual_year variables
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    # linear interpolation between the years of analysis, forward-stepping, up to gettime(tt)
    ## interpolation here is also linear for 2015-2020
    if is_first(tt)
        frac = annual_year - 2015
        fraction_timestep = frac / ((gettime(tt)) - 2015)

        v.pt_g_preliminarygmst_ann[yr] = (v.pt_g_preliminarygmst[tt]) * (fraction_timestep) + p.rt_g0_baseglobaltemp * (1 - fraction_timestep)  # difference_to_next_tt * ratio_to_there
    else
        frac = annual_year - gettime(tt - 1)
        fraction_timestep = frac / ((gettime(tt)) - (gettime(tt - 1))) # check if +1 might also need to feature here.

        v.pt_g_preliminarygmst_ann[yr] = (v.pt_g_preliminarygmst[tt] ) * (fraction_timestep) + v.pt_g_preliminarygmst[tt - 1] * (1 - fraction_timestep) # difference_to_next_tt * ratio_to_there
    end
    # Without surface albedo, just equal
    v.rt_g_globaltemperature_ann[yr] = v.pt_g_preliminarygmst_ann[yr]

    # Adding variability to global temperature
    if use_variability
        g_variationtemp = rand(Normal(0.0, sqrt(v.tvarerr_g_globaltemperatureerrorvariance)))
    else
        g_variationtemp = 0
    end
    if yr == 1
        v.rt_g_globaltemperature_ann[yr] = v.pt_g_preliminarygmst_ann[yr] + g_variationtemp
    else
        ptdiff = v.pt_g_preliminarygmst_ann[yr] - v.pt_g_preliminarygmst_ann[yr - 1]
        v.rt_g_globaltemperature_ann[yr] = v.tvarconst_g_globaltemperatureintercept + v.tvarar_g_globaltemperatureautoreg * (ptdiff + v.rt_g_globaltemperature_ann[yr - 1]) + v.tvargmst_g_globaltemperaturesmoothdep * v.rt_g_globaltemperature_ann[yr] + g_variationtemp
    end

    # Setting regional temperature (uses globally shocked temperature)
    for r in d.region
        if use_variability
            r_variationtemp = rand(Normal(0.0, sqrt(v.tvarerr_regionaltemperatureerrorvariance[r]))) # NEW - with variation
            # r_variationtemp = r_variationtemp * p.var_multiplier # for sensitivity analysis for variability
        else
            r_variationtemp = 0 # NEW - no variation
        end
        if yr == 1
            v.rtl_realizedtemperature_ann[yr, r] = v.rt_g_globaltemperature_ann[yr] * p.ampf_amplification[r] + r_variationtemp
        else
	    # Determine correction to conform in expectation to ampf * T_g
            v.rtl_realizedtemperature_ann[yr, r] = (p.ampf_amplification[r] / (v.tvarar_regionaltemperatureautoreg[r] + v.tvargmst_regionaltemperatureglobaldep[r])) * (v.tvarconst_regionaltemperatureintercept[r] + v.tvarar_regionaltemperatureautoreg[r] * (p.ampf_amplification[r] * ptdiff + v.rtl_realizedtemperature_ann[yr - 1, r]) + v.tvargmst_regionaltemperatureglobaldep[r] * v.rt_g_globaltemperature_ann[yr]) + r_variationtemp
        end
    end

    # Land average temperature
    v.rtl_g_landtemperature_ann[yr] = sum(v.rtl_realizedtemperature_ann[yr, :]' .* p.area') / sum(p.area)

    # Ocean average temperature
    v.rto_g_oceantemperature_ann[yr] = (p.area_e_eartharea * v.rt_g_globaltemperature_ann[yr] - sum(p.area) * v.rtl_g_landtemperature_ann[yr]) / (p.area_e_eartharea - sum(p.area))

end

@defcomp ClimateTemperature begin

    region = Index()
    year = Index()

    # Basic parameters
    area = Parameter(index=[region], unit="km^2")
    y_year_0 = Parameter(unit="year")
    y_year = Parameter(index=[time], unit="year")
    y_year_ann = Parameter(index=[year], unit="year")
    area_e_eartharea = Parameter(unit="km2", default=5.1e8)
    use_seaice = Parameter{Bool}()

    # Initial temperature outputs
    rt_g0_baseglobaltemp = Parameter(unit="degreeC", default=0.9461666666666667) # needed for feedback in CO2 cycle component
    rtl_0_baselandtemp = Variable(index=[region], unit="degreeC")
    rtl_g0_baselandtemp = Variable(unit="degreeC") # needed for feedback in CH4 and N2O cycles

    # variability parameters
    tvarorder_arestimatesrows = Parameter(index=[region], unit="none")
    tvarseed_coefficientsrandomseed = Parameter(unit="none")
    tvarerr_g_globaltemperatureerrorvariance = Variable(unit="degreeC^2")
    tvarconst_g_globaltemperatureintercept = Variable(unit="degreeC")
    tvarar_g_globaltemperatureautoreg = Variable(unit="none")
    tvargmst_g_globaltemperaturesmoothdep = Variable(unit="none")
    tvarerr_regionaltemperatureerrorvariance = Variable(index=[region], unit="degreeC^2")
    tvarconst_regionaltemperatureintercept = Variable(index=[region], unit="degreeC")
    tvarar_regionaltemperatureautoreg = Variable(index=[region], unit="none")
    tvargmst_regionaltemperatureglobaldep = Variable(index=[region], unit="none")

    # Rate of change of forcing
    fant_anthroforcing = Parameter(index=[time], unit="W/m2")
    ft_0_totalforcing0 = Parameter(unit="W/m2", default=3.202)
    fsd_g_0_directsulphate0 = Parameter(unit="W/m2", default=-0.46666666666666673)
    fsi_g_0_indirectsulphate0 = Parameter(unit="W/m2", default=-0.17529148701061475)

    # Climate sensitivity calculations
    tcr_transientresponse = Parameter(unit="degreeC", default=1.7666666666666668)
    frt_warminghalflife = Parameter(unit="year", default=28.333333333333332)

    ecs_climatesensitivity = Variable(unit="degreeC")

    # Unadjusted temperature calculations
    fslope_CO2forcingslope = Parameter(unit="W/m2", default=5.5)
    pt_g_preliminarygmst = Variable(index=[time], unit="degreeC")
    pt_g_preliminarygmst_ann = Variable(index=[year], unit="degreeC")

    # Surface albedo parameters
    alb_t_switch = Parameter(default=10.0)
    alb_saf_quadr_mean_t2_coeff = Parameter(default=-0.001955963673212411)
    alb_saf_quadr_mean_t1_coeff = Parameter(default=-0.006407351323331157)
    alb_saf_quadr_mean_t0_coeff = Parameter(default=0.36372124028281877)
    alb_saf_quadr_std = Parameter(default=0.10908771841232413)
    alb_saf_lin_mean = Parameter(default=0.0636051293985324)
    alb_saf_lin_std = Parameter(default=0.02907749205894582)
    alb_emulator_rand = Parameter(default=0)

    # Surface albedo internal variables
    alb_fsaf_y0 = Variable(unit="W/m2")
    alb_saf_y0 = Variable(unit="W/m2/degC")
    alb_fsaf_ecs = Variable(unit="W/m2")
    alb_saf_ecs = Variable(unit="W/m2/degC")
    alb_fsaf_t_switch = Variable(unit="W/m2")

    # Global outputs
    rt_g_globaltemperature = Variable(index=[time], unit="degreeC")
    rt_g_globaltemperature_ann = Variable(index=[year], unit="degreeC")
    rto_g_oceantemperature = Variable(index=[time], unit="degreeC")
    rto_g_oceantemperature_ann = Variable(index=[year], unit="degreeC")
    rtl_g_landtemperature = Variable(index=[time], unit="degreeC")
    rtl_g_landtemperature_ann = Variable(index=[year], unit="degreeC")

    # Regional outputs
    ampf_amplification = Parameter(index=[region])

    rtl_realizedtemperature = Variable(index=[time, region], unit="degreeC")
    rtl_realizedtemperature_ann = Variable(index=[year, region], unit="degreeC")

    function init(p, v, d)
        for rr in d.region
            v.rtl_0_baselandtemp[rr] = p.rt_g0_baseglobaltemp * p.ampf_amplification[rr]
        end

        # Equation 21 from Hope (2006): initial global land temperature
        v.rtl_g0_baselandtemp = sum(v.rtl_0_baselandtemp .* p.area) / sum(p.area)

        # Inclusion of transient climate response from Hope (2009)
        v.ecs_climatesensitivity = p.tcr_transientresponse / (1. - (p.frt_warminghalflife / 70.) * (1. - exp(-70. / p.frt_warminghalflife)))

        ## Surface albedo internal variables
        v.alb_fsaf_y0 = ((p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.rt_g0_baseglobaltemp) + (p.alb_saf_quadr_std * p.rt_g0_baseglobaltemp) * p.alb_emulator_rand)
        v.alb_saf_y0 = (p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t1_coeff * p.rt_g0_baseglobaltemp + p.alb_saf_quadr_mean_t0_coeff) + p.alb_saf_quadr_std * p.alb_emulator_rand

        if v.ecs_climatesensitivity <= p.alb_t_switch
            v.alb_fsaf_ecs = (p.alb_saf_quadr_mean_t2_coeff * (v.ecs_climatesensitivity^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.ecs_climatesensitivity^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.ecs_climatesensitivity) + (p.alb_saf_quadr_std * v.ecs_climatesensitivity) * p.alb_emulator_rand
        else
            v.alb_fsaf_ecs = (p.alb_saf_quadr_mean_t2_coeff * ( p.alb_t_switch^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * ( p.alb_t_switch^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.alb_t_switch + p.alb_saf_lin_mean * (v.ecs_climatesensitivity - p.alb_t_switch)) + (p.alb_saf_quadr_std * p.alb_t_switch + p.alb_saf_lin_std * (v.ecs_climatesensitivity - p.alb_t_switch)) * p.alb_emulator_rand
        end
        v.alb_saf_ecs = v.alb_fsaf_ecs / v.ecs_climatesensitivity
        v.alb_fsaf_t_switch = (p.alb_saf_quadr_mean_t2_coeff * (p.alb_t_switch^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.alb_t_switch^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.alb_t_switch) + (p.alb_saf_quadr_std * p.alb_t_switch) * p.alb_emulator_rand

        tvarconst, tvarar, tvargmst = tvar_getcoeffs(findfirst(arests[!, :region] .== "global"))
        v.tvarconst_g_globaltemperatureintercept = tvarconst
        v.tvarar_g_globaltemperatureautoreg = tvarar
        v.tvargmst_g_globaltemperaturesmoothdep = tvargmst
        v.tvarerr_g_globaltemperatureerrorvariance = tvar_geterror(findfirst(arests[!, :region] .== "global"))

        tvarconst, tvarar, tvargmst, tvarerr = tvar_getregions(p.tvarorder_arestimatesrows, tvar_getcoeffs)
        v.tvarconst_regionaltemperatureintercept[:] = tvarconst
        v.tvarar_regionaltemperatureautoreg[:] = tvarar
        v.tvargmst_regionaltemperatureglobaldep[:] = tvargmst
        v.tvarerr_regionaltemperatureerrorvariance[:] = tvarerr
        if p.tvarseed_coefficientsrandomseed != 0
            rng = MersenneTwister(trunc(Int, p.tvarseed_coefficientsrandomseed))
            tvarconst, tvarar, tvargmst = rand(rng, tvar_getmvnormal(findfirst(arests[!, :region] .== "global")))
            v.tvarconst_g_globaltemperatureintercept = tvarconst
            v.tvarar_g_globaltemperatureautoreg = tvarar
            v.tvargmst_g_globaltemperaturesmoothdep = tvargmst

            tvarconst, tvarar, tvargmst, tvarerr = tvar_getregions(p.tvarorder_arestimatesrows, rr -> rand(rng, tvar_getmvnormal(rr)))
            v.tvarconst_regionaltemperatureintercept[:] = tvarconst
            v.tvarar_regionaltemperatureautoreg[:] = tvarar
            v.tvargmst_regionaltemperatureglobaldep[:] = tvargmst
        end
    end

    function run_timestep(p, v, d, tt)
        # Rate of change of grand total forcing
        if is_first(tt)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            rate_fant = 0 # inferred from spreadsheet
            deltat = p.y_year[tt] - p.y_year_0
        elseif tt == TimestepIndex(2)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            rate_fant = (p.fant_anthroforcing[tt - 1] - fant0) / (p.y_year[tt - 1] - p.y_year_0)
            deltat = p.y_year[tt] - p.y_year[tt - 1]
        else
            rate_fant = (p.fant_anthroforcing[tt - 1] - p.fant_anthroforcing[tt - 2]) / (p.y_year[tt - 1] - p.y_year[tt - 2])
            deltat = p.y_year[tt] - p.y_year[tt - 1]
        end

        BB = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * rate_fant
        EXPT = exp(-deltat / p.frt_warminghalflife)

        if is_first(tt)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            AA = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * fant0
            v.pt_g_preliminarygmst[tt] = p.rt_g0_baseglobaltemp + (AA - p.rt_g0_baseglobaltemp) * (1 - EXPT) # Drop BB because 0
        else
            AA = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * p.fant_anthroforcing[tt - 1]
            v.pt_g_preliminarygmst[tt] = v.pt_g_preliminarygmst[tt - 1] + (AA - p.frt_warminghalflife * BB - v.pt_g_preliminarygmst[tt - 1]) * (1 - EXPT) + deltat * BB
        end

        if !p.use_seaice
            # Without surface albedo, just equal
            v.rt_g_globaltemperature[tt] = v.pt_g_preliminarygmst[tt]
        else
            alb_saf_approx = (((p.alb_saf_quadr_mean_t2_coeff * (v.pt_g_preliminarygmst[tt]^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.pt_g_preliminarygmst[tt]^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.pt_g_preliminarygmst[tt]) + (p.alb_saf_quadr_std * v.pt_g_preliminarygmst[tt]) * p.alb_emulator_rand) - v.alb_fsaf_y0) / (v.pt_g_preliminarygmst[tt] - p.rt_g0_baseglobaltemp)
            alb_saf_adjust = alb_saf_approx - v.alb_saf_ecs
            ecs_alb_adj = v.ecs_climatesensitivity / (1 - v.ecs_climatesensitivity * alb_saf_adjust / (log(2) * p.fslope_CO2forcingslope))
            frt_alb_adj = p.frt_warminghalflife / (1 - v.ecs_climatesensitivity * alb_saf_adjust / (log(2) * p.fslope_CO2forcingslope))
            if is_first(tt)
                alb_fsaf_adjust = ((p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.rt_g0_baseglobaltemp) + (p.alb_saf_quadr_std * p.rt_g0_baseglobaltemp) * p.alb_emulator_rand) - alb_saf_approx * p.rt_g0_baseglobaltemp
                v.rt_g_globaltemperature[tt] = p.rt_g0_baseglobaltemp + (ecs_alb_adj * (fant0 + alb_fsaf_adjust) / (log(2) * p.fslope_CO2forcingslope) - p.rt_g0_baseglobaltemp) * (1 - exp(-deltat / frt_alb_adj))
            else
                if v.rt_g_globaltemperature[tt - 1] <= p.alb_t_switch
                    alb_fsaf_adjust = ((p.alb_saf_quadr_mean_t2_coeff * (v.rt_g_globaltemperature[tt - 1]^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.rt_g_globaltemperature[tt - 1]^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.rt_g_globaltemperature[tt - 1]) + (p.alb_saf_quadr_std * v.rt_g_globaltemperature[tt - 1]) * p.alb_emulator_rand) - alb_saf_approx * v.rt_g_globaltemperature[tt - 1]
                else
                    alb_fsaf_adjust = v.alb_fsaf_t_switch + p.alb_saf_lin_mean * (v.rt_g_globaltemperature[tt - 1] - p.alb_t_switch) + (p.alb_saf_lin_std * (v.rt_g_globaltemperature[tt - 1] - p.alb_t_switch)) * p.alb_emulator_rand - alb_saf_approx * v.rt_g_globaltemperature[tt - 1]
                end
                v.rt_g_globaltemperature[tt] = v.rt_g_globaltemperature[tt - 1] + (ecs_alb_adj * (p.fant_anthroforcing[tt - 1] - rate_fant * frt_alb_adj + alb_fsaf_adjust) / (log(2) * p.fslope_CO2forcingslope) - v.rt_g_globaltemperature[tt - 1]) * (1 - exp(-deltat / frt_alb_adj)) + ecs_alb_adj * (rate_fant * deltat) / (log(2) * p.fslope_CO2forcingslope)
            end
        end

        # Adding adjustment, from Hope (2009)
        for rr in d.region
            v.rtl_realizedtemperature[tt, rr] = v.rt_g_globaltemperature[tt] * p.ampf_amplification[rr]
        end

        # Land average temperature
        v.rtl_g_landtemperature[tt] = sum(v.rtl_realizedtemperature[tt, :]' .* p.area') / sum(p.area)

        # Ocean average temperature
        v.rto_g_oceantemperature[tt] = (p.area_e_eartharea * v.rt_g_globaltemperature[tt] - sum(p.area) * v.rtl_g_landtemperature[tt]) / (p.area_e_eartharea - sum(p.area))

        # annual interpolation of the timestep-calculated variables.
        if is_first(tt)
            for annual_year = 2015:gettime(tt) # because p["y_year_0"] does not exist (globally)
                calc_temp(p, v, d, tt, annual_year) # NEW -- introduce function for calulating temperature
            end
        else
            for annual_year = (gettime(tt - 1) + 1):gettime(tt)
                calc_temp(p, v, d, tt, annual_year) # NEW -- use  function for calulating temperature
            end
        end
    end
end

function tvar_getcoeffs(rr::Int64)
    tvarconst = arests[rr, :intercept][1]
    tvarar = arests[rr, :ar][1]
    tvargmst = arests[rr, :gmst][1]

    tvarconst, tvarar, tvargmst
end

function tvar_geterror(rr::Int64)
    arests[rr, :varerror][1]
end

function tvar_getmvnormal(rr::Int64)
    tvarcoeffs = tvar_getcoeffs(rr)

    names = ["intercept", "ar", "gmst"]
    tvarsigma = zeros(3, 3)
    for ii in 1:3
        for jj in 1:3
            tvarsigma[ii, jj] = arests[rr, Symbol("var$(names[ii])*$(names[jj])")][1]
        end
    end

    MvNormal([tvarcoeffs...], tvarsigma)
end

function tvar_getregions(regionorder::Vector{Union{Missing,Float64}}, coeffproc::Function=tvar_getcoeffs)
    tvar_getregions(convert(Vector{Float64}, regionorder), coeffproc)
end

function tvar_getregions(regionorder::Vector{Float64}, coeffproc::Function=tvar_getcoeffs)
    regionorder = convert(Vector{Int64}, regionorder)
    tvar_getregions(regionorder, coeffproc)
end

function tvar_getregions(regionorder::Vector{Int64}, coeffproc::Function=tvar_getcoeffs)
    regionorder = convert(Vector{Int64}, regionorder)
    tvarconst = zeros(length(regionorder))
    tvarar = zeros(length(regionorder))
    tvargmst = zeros(length(regionorder))
    tvarerr = zeros(length(regionorder))
    for rr in 1:length(regionorder)
        tvarconst[rr], tvarar[rr], tvargmst[rr] = coeffproc(regionorder[rr])
        tvarerr[rr] = tvar_geterror(regionorder[rr])
    end

    tvarconst, tvarar, tvargmst, tvarerr
end

function addclimatetemperature(model::Model, use_seaice::Bool)
    climtemp = add_comp!(model, ClimateTemperature)

    climtemp[:use_seaice] = use_seaice

    region_keys = Mimi.dim_keys(model.md, :region)
    regionmapping = Dict{String,String}("EU" => "eu", "USA" => "usa", "OECD" => "oth", "Africa" => "afr", "China" => "chi+", "SEAsia" => "ind+", "LatAmerica" => "lat", "USSR" => "rus+")

    tvarorder = zeros(Int64, length(region_keys))
    for rr in 1:length(region_keys)
        tvarorder[rr] = findfirst(arests[!, :region] .== regionmapping[region_keys[rr]])
    end
    climtemp[:tvarorder_arestimatesrows] = tvarorder
    climtemp[:tvarseed_coefficientsrandomseed] = 0 # random draw not used

    return climtemp
end


text:
"""
    nlsolve!(nlsolver::NLSolver, nlcache::Union{NLNewtonCache,NLNewtonConstantCache}, integrator)

Perform numerically stable modified Newton iteration that is specialized for implicit
methods (see [^HS96] and [^HW96]).

It solves

```math
G(z) = dt⋅f(tmp + γ⋅z, p, t + c⋅h) - z = 0
```

by iterating

```math
W Δᵏ = dt⋅f(tmp + γ⋅zᵏ, p, t + c⋅h) - zᵏ
zᵏ⁺¹ = zᵏ + Δᵏ
```

where `W = M - dt⋅γJ`, `M` is the mass matrix, `dt` is the step size, `γ` is a
constant, and `J` is the Jacobian matrix.

It returns the tuple `(z, η, iter, fail_convergence)`, where `z` is the solution, `η` is
used to measure the iteration error (see [^HW96]), `iter` is the number of iterations, and
`fail_convergence` reports whether the algorithm succeeded.

[^HS96]: M.E.Hoseaa and L.F.Shampine, "Analysis and implementation of TR-BDF2",
Applied Numerical Mathematics, Volume 20, Issues 1–2, February 1996, Pages
21-37.
[doi:10.1016/0168-9274(95)00115-8](https://doi.org/10.1016/0168-9274(95)00115-8)

[^HW96]: Ernst Hairer and Gerhard Wanner, "Solving Ordinary Differential
Equations II, Springer Series in Computational Mathematics. ISBN
978-3-642-05221-7. Section IV.8.
[doi:10.1007/978-3-642-05221-7](https://doi.org/10.1007/978-3-642-05221-7)
"""
@muladd function nlsolve!(nlsolver::NLSolver, nlcache::NLNewtonConstantCache, integrator)
  @unpack t,dt,uprev,u,p = integrator
  @unpack z,tmp,κtol,c,γ,max_iter = nlsolver
  W = nlcache.W

  # precalculations
  mass_matrix = integrator.f.mass_matrix
  f = nlsolve_f(integrator)
  tstep = t + c*dt
  η = max(nlsolver.ηold,eps(eltype(integrator.opts.reltol)))^(0.8)

  # Newton iteration
  local ndz
  fail_convergence = true
  iter = 0
  while iter < max_iter
    iter += 1
    integrator.destats.nnonliniter += 1

    # evaluate function
    u = @. tmp + γ * z
    if mass_matrix === I
      ztmp = dt .* f(u, p, tstep) .- z
    else
      ztmp = dt .* f(u, p, tstep) .- mass_matrix * z
    end
    integrator.destats.nf += 1
    if DiffEqBase.has_invW(f)
      dz = _reshape(W * -_vec(ztmp), axes(ztmp)) # Here W is actually invW
    else
      dz = _reshape(W \ _vec(ztmp), axes(ztmp))
    end
    integrator.destats.nsolve += 1

    # compute norm of residuals
    iter > 1 && (ndzprev = ndz)
    ndz = integrator.opts.internalnorm(dz, tstep)

    # check divergence (not in initial step)
    if iter > 1
      θ = ndz / ndzprev
      if θ > 1 || ndz * θ^(max_iter - iter) > κtol * (1 - θ)
        # Newton method diverges
        break
      end
    end

    # update solution
    z = z .- dz

    # check stopping criterion
    iter > 1 && (η = θ / (1 - θ))
    if η * ndz < κtol && (iter > 1 || iszero(ndz) || !iszero(integrator.success_iter))
      # Newton method converges
      fail_convergence = false
      break
    end
  end

  if fail_convergence
    integrator.destats.nnonlinconvfail += 1
  end
  integrator.force_stepfail = fail_convergence
  z, η, iter, fail_convergence
end

@muladd function nlsolve!(nlsolver::NLSolver, nlcache::NLNewtonCache, integrator)
  @unpack t,dt,uprev,u,p,cache = integrator
  @unpack z,dz,tmp,ztmp,k,κtol,c,γ,max_iter = nlsolver
  @unpack W, new_W = nlcache
  cache = unwrap_cache(integrator, true)

  # precalculations
  mass_matrix = integrator.f.mass_matrix
  f = nlsolve_f(integrator)
  vecztmp = vec(ztmp); vecz = vec(z); vecdz = vec(dz)
  tstep = t + c*dt
  η = max(nlsolver.ηold,eps(eltype(integrator.opts.reltol)))^(0.8)

  # Newton iteration
  local ndz
  fail_convergence = true
  iter = 0
  while iter < max_iter
    iter += 1
    integrator.destats.nnonliniter += 1

    # evaluate function
    @. u = tmp + γ*z
    f(k, u, p, tstep)
    integrator.destats.nf += 1
    if mass_matrix === I
      @. ztmp = dt*k - z
    else
      mul!(vecztmp,mass_matrix,vecz)
      @. ztmp = dt*k - ztmp
    end
    if DiffEqBase.has_invW(f)
      mul!(vecdz,W,vecztmp) # Here W is actually invW
      @. vecdz = -vecdz
    else
      cache.linsolve(vecdz,W,vecztmp,iter == 1 && new_W)
    end
    integrator.destats.nsolve += 1

    # compute norm of residuals
    iter > 1 && (ndzprev = ndz)
    ndz = integrator.opts.internalnorm(dz, tstep)

    # check divergence (not in initial step)
    if iter > 1
      θ = ndz / ndzprev
      if θ > 1 || ndz * θ^(max_iter - iter) > κtol * (1 - θ)
        # Newton method diverges
        break
      end
    end

    # update solution
    @. z = z - dz

    # check stopping criterion
    iter > 1 && (η = θ / (1 - θ))
    if η * ndz < κtol && (iter > 1 || iszero(ndz) || !iszero(integrator.success_iter))
      # Newton method converges
      fail_convergence = false
      break
    end
  end

  if fail_convergence
    integrator.destats.nnonlinconvfail += 1
  end
  integrator.force_stepfail = fail_convergence
  z, η, iter, fail_convergence
end


text:
# [SOREN] TRUST REGION MODIFICATION OF NEWTON'S METHOD
# [N&W] Numerical optimization
# [Yuan] A review of trust region algorithms for optimization
abstract type TRSPSolver end
abstract type NearlyExactTRSP <: TRSPSolver end

include("solvers/NWI.jl")
include("solvers/Dogleg.jl")
include("solvers/NTR.jl")
include("solvers/TCG.jl")
#include("subproblemsolvers/TRS.jl") just make an example instead of relying onTRS.jl

function tr_return(;λ, ∇f, H, s, interior, solved, hard_case, Δ, m=nothing)
	m = m isa Nothing ? dot(∇f, s) + dot(s, H * s)/2 : m
	(p=s, mz=m, interior=interior, λ=λ, hard_case=hard_case, solved=solved, Δ=Δ)
end

function update_H!(H, h, λ=nothing)
  T = eltype(h)
  n = length(h)
  if !(λ ==T(0))
    for i = 1:n
      @inbounds H[i, i] = λ isa Nothing ? h[i] : h[i] + λ
    end
  end
  H
end


text:
using LazyGraphs
using Test

@testset "LazyGraphs.jl" begin
    # Write your tests here.
end


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

bernoulli(p) = rand() <= p

function poisson(lambda)
    k = 0
    p = exp(-lambda)
    s = p
    x = rand()
    while x > s
        k += 1
        p = p * lambda / k
        s += p
    end
    k
end

struct Householder{T} <: AbstractMatrix{T}
    v::Vector{T}
end

function Base.size(h::Householder)
    l = length(h.v)
    l,l
end

function Base.getindex(h::Householder, i, j)
    I[i,j] - 2*h.v[i]*h.v[j]
end

function random_householder(dim)
    v = randn(dim)
    normalize!(v)
    Householder(v)
end

function random_orthogonal(dim)
    ret = if dim == 1
        [(-1.)^bernoulli(0.5)]
    else
        dim1 = dim - 1
        m = random_householder(dim)
        rot1 = random_orthogonal(dim1)
        m * [rot1 zeros(dim1, 1); zeros(1, dim1) 1]
    end
    @assert ret * ret' ≈ Matrix(I,dim,dim)
    ret
end

function random_embedding(dim_target, dim_src)
    inc = Matrix(I,dim_target, dim_src)
    rot = random_orthogonal(dim_target)
    M = SMatrix{dim_target, dim_src}
    M(rot * inc)
end

function create_ball_points(npoints, dim_src;
                            codim = poisson(0.3),
                            p_boundary=1,
                            p_rep=1/sqrt(npoints),
                            p_shuffle=0.5,
                           )
    dim_target = dim_src + codim
    @assert dim_src <= dim_target
    @assert 0 <= p_boundary <= 1
    @assert 0 <= p_rep <= 1
    @assert 0 <= p_shuffle <= 1
    F = Float64
    P = SVector{dim_src, F}
    pts = P[]
    center = randn(dim_src)
    R = 10*rand()
    ball = MB.SqBall(center, R^2)
    while length(pts) < npoints
        dir = normalize!(randn(dim_src))
        r = if bernoulli(p_boundary)
            R
        else
            R*rand()
        end
        pt = center + r*dir
        while true
            push!(pts, pt)
            bernoulli(p_rep) || break
            (length(pts) == npoints) && break
        end
    end
    bernoulli(p_shuffle) && shuffle!(pts)
    @assert length(pts) == npoints
    @assert eltype(pts) == P
    if dim_target > dim_src
        embedding = random_embedding(dim_target, dim_src)
        center = embedding * center
        pts = map(pt -> embedding*pt, pts)
        ball = MB.SqBall(center, R^2)
    end
    make_ball_containing_pts(ball, pts)
end

function make_ball_containing_pts(ball, pts)
    R2 = ball.sqradius
    center = ball.center
    for pt in pts
        R2pt = MB.sqdist(pt, center)
        if R2pt > R2
            @assert R2pt ≈ R2
            R2 = R2pt
        end
    end
    MB.SqBall(center, R2), pts
end

function random_test(alg, npoints, dim;
                     rtol_inside=nothing,
                     atol_inside=nothing,
                     rtol_radius=nothing,
                     allow_broken::Bool=false,
                     kw...)
    ball_ref, pts = create_ball_points(npoints, dim; kw...)
    @assert MB.allinside(pts, ball_ref)
    P = eltype(pts)
    F = eltype(P)

    c, r = boundingsphere(pts, alg)
    ball = MB.SqBall(c, r^2)

    r_ref = MB.radius(ball_ref)
    if rtol_radius == nothing
        rtol_radius = sqrt(eps(max(r, r_ref)))
    end
    issmall = r <= r_ref || isapprox(r, r_ref; rtol=rtol_radius)
    if allow_broken && !issmall
        @test_broken issmall
    else
        @test issmall
    end

    if rtol_inside == nothing
        rtol_inside = 1e-6
    end
    if atol_inside == nothing
        atol_inside = 100eps(F)
    end
    contains_all_points = MB.allinside(pts, ball, rtol=rtol_inside, atol=atol_inside)
    if allow_broken && !contains_all_points
        @test_broken contains_all_points
        # pts = map(pt -> map(BigFloat,pt), pts)
        # @show length(pts)
        # @show length(first(pts))
        # @show ball_ref
        # @show pts
        # @show ball
    else
        @test contains_all_points
    end
end