text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -7.107954588851326 <= q <= 8.439913276188518)

                     add_NL_constraint(m, :(sigmoid(-0.4953615570310643 + -0.5197659840811282*sigmoid(0.7941458217392698 + 0.49653685598476827*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.07122185861106178*sigmoid(-0.06420985369672083 + 0.9131512287471262*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.004105284291951339*sigmoid(-0.8124442669202567 + 0.4437121866105662*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + sigmoid(0.7344840166166842 + -0.5219551557665656*sigmoid(0.7941458217392698 + 0.49653685598476827*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.12237894168399821*sigmoid(-0.06420985369672083 + 0.9131512287471262*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.5945370485807278*sigmoid(-0.8124442669202567 + 0.4437121866105662*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + sigmoid(-0.7546740103129919 + -0.5743271753296177*sigmoid(0.7941458217392698 + 0.49653685598476827*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.959961363428667*sigmoid(-0.06420985369672083 + 0.9131512287471262*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.0685664132974031*sigmoid(-0.8124442669202567 + 0.4437121866105662*sigmoid(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*sigmoid(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*sigmoid(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using CairoMakie
using Images, FileIO
using ImageTransformations
using WindMixing
using JLD2
using WindMixing: animate_profiles_fluxes_final

# convective_adjustment = parse(Bool, ARGS[1])
convective_adjustment = false
# params_type = ARGS[2]
params_type = "18simBFGST0.8nograd"
# NN_type = ARGS[3]
NN_type = "leakyrelu"
# num = parse(Int, ARGS[1])
# NN_version = ARGS[2]

num=1
NN_version = "old"

if convective_adjustment
    CA_str = "CA"
else
    CA_str = "noCA"
end

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot", "z.png")))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
)

EXTRACTED_FILE_DIR = joinpath(pwd(), "extracted_training_output")
if NN_version == "old"
    EXTRACTED_FILE_PATH = joinpath(EXTRACTED_FILE_DIR, "NDE_18sim_windcooling_windheating_$(params_type)_divide1f5_gradient_smallNN_$(NN_type)_rate_2e-4_T0.8_extracted.jld2")
    FILE_DIR = joinpath(pwd(), "final_results", "18sim_old")
else
    EXTRACTED_FILE_PATH = joinpath(EXTRACTED_FILE_DIR, "NDE_18sim_windcooling_windheating_$(params_type)_divide1f5_gradient_smallNN_$(NN_type)_rate_2e-4_T0.8_1e-4_extracted.jld2")
    FILE_DIR = joinpath(pwd(), "final_data_2", "NDE_18sim_windcooling_windheating_$(params_type)_divide1f5_gradient_smallNN_$(NN_type)_rate_2e-4_T0.8_1e-4_$(CA_str)")
end

if !ispath(FILE_DIR)
    mkdir(FILE_DIR)
end

file = jldopen(EXTRACTED_FILE_PATH)
train_files = file["training_info/train_files"]
loss_scalings = file["training_info/loss_scalings"]
train_parameters = file["training_info/parameters"]
close(file)

test_files_all = [
    # "wind_-5e-4_cooling_3e-8_new",
    # "wind_-3.5e-4_cooling_3e-8_new",
    # "wind_-2e-4_cooling_3e-8_new",

    # "wind_-5e-4_heating_-3e-8_new",
    # "wind_-3.5e-4_heating_-3e-8_new",
    # "wind_-2e-4_heating_-3e-8_new",

    # "wind_-5e-4_cooling_2e-8_new",
    # "wind_-5e-4_cooling_1e-8_new",
    # "wind_-3.5e-4_cooling_2e-8_new",
    # "wind_-3.5e-4_cooling_1e-8_new",
    # "wind_-2e-4_cooling_2e-8_new",
    # "wind_-2e-4_cooling_1e-8_new",
    # "wind_-5e-4_heating_-2e-8_new",
    # "wind_-5e-4_heating_-1e-8_new",
    # "wind_-3.5e-4_heating_-2e-8_new",
    # "wind_-3.5e-4_heating_-1e-8_new",
    # "wind_-2e-4_heating_-2e-8_new",
    # "wind_-2e-4_heating_-1e-8_new",

    # "wind_-4.5e-4_cooling_2.5e-8",
    # "wind_-2.5e-4_cooling_1.5e-8",
    # "wind_-4.5e-4_cooling_1.5e-8",
    # "wind_-2.5e-4_cooling_2.5e-8",

    # "wind_-4.5e-4_heating_-2.5e-8",
    # "wind_-2.5e-4_heating_-1.5e-8",
    # "wind_-4.5e-4_heating_-1.5e-8",
    # "wind_-2.5e-4_heating_-2.5e-8",

    # "wind_-5.5e-4_diurnal_5.5e-8",
    # "wind_-1.5e-4_diurnal_5.5e-8",

    # "wind_-5.5e-4_new",

    # "wind_-5.5e-4_heating_-3.5e-8",
    # "wind_-1.5e-4_heating_-3.5e-8",
    # "wind_-5.5e-4_cooling_3.5e-8",
    "wind_-1.5e-4_cooling_3.5e-8",
]

test_files = [test_files_all[num]]

ν₀ = train_parameters["ν₀"]
ν₋ = train_parameters["ν₋"]
ΔRi = train_parameters["ΔRi"]
Riᶜ = train_parameters["Riᶜ"]
Pr = train_parameters["Pr"]

# solve_oceananigans_modified_pacanowski_philander_nn(test_files, EXTRACTED_FILE_PATH, FILE_DIR,
#                                                         timestep=1, convective_adjustment=convective_adjustment)

diurnal = occursin("diurnal", test_files[1])

if num <= 18
    DIR_NAME = "train_$(test_files[1])"
    animation_type = "Training"
elseif num >= 19 && num <= 26
    DIR_NAME = "test_$(test_files[1])"
    animation_type = "Interpolating"
else
    DIR_NAME = "test_$(test_files[1])"
    animation_type = "Extrapolating"
end

DIR_NAME = "test_$(test_files[1])"
animation_type = "Extrapolating"

n_trainings = length(train_files)

plot_data = NDE_profile_oceananigans(joinpath(FILE_DIR, DIR_NAME), train_files, test_files,
                                  ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                  loss_scalings=loss_scalings,
                                  OUTPUT_PATH=joinpath("C:\\Users\\xinle\\Documents\\OceanParameterizations.jl\\final_results\\18sim_old", DIR_NAME, "profiles_fluxes_oceananigans.jld2"))

# animate_profiles_fluxes_final(plot_data, axis_images,
#                     joinpath(FILE_DIR, DIR_NAME, "$(test_files[1])_60fps"),
#                     animation_type=animation_type, n_trainings=n_trainings, training_types="Wind + Cooling, Wind + Heating", fps=60)

ext: jl


text:
__precompile__()

module TSAnalysis

    # Libraries
    using Distributed;
    using Dates, Logging;
    using LinearAlgebra, Distributions, Statistics;

    # Custom dependencies
    const local_path = dirname(@__FILE__);
    include("$local_path/types.jl");
    include("$local_path/methods.jl");
    include("$local_path/kalman.jl");

    # Export
    export JVector, JArray, KalmanSettings, ImmutableKalmanSettings, MutableKalmanSettings, KalmanStatus;
    export mean_skipmissing, std_skipmissing, is_vector_in_matrix, demean, lag, companion_form, ext_companion_form;
    export kfilter!, kforecast, ksmoother;
end


text:
module ThreadManager

export Thread

type Thread
	id::Integer
	mutex::Threads.Mutex
	Thread(id, m) = new(id, m)
end

Mutex = Threads.Mutex()
List = Thread[]

function safe_call(t::Thread, f::Function)
	Threads.lock(t.mutex)
	f(t)
	Threads.unlock(t.mutex)
end

function reset()
	global Mutex
	global List
	Mutex = Threads.Mutex()
	List = Thread[]
end

function run(a = Function[])
	reset()
	max=Threads.nthreads()
	i=0; Threads.@threads for f in a
		t=Thread(Threads.threadid(), Mutex)
		push!(List, t)
		f(t)
		if i >= max break end
		i+=1
  end
end

end # ThreadManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
function ManifoldTests.test_reversediff(M::Manifold, pts, tv)
    return for (p, X) in zip(pts, tv)
        exp_f(t) = distance(M, p, exp(M, p, t[1] * X))
        d12 = norm(M, p, X)
        for t in 0.1:0.1:0.9
            Test.@test d12 ≈ ReverseDiff.gradient(exp_f, [t])[1]
        end

        retract_f(t) = distance(M, p, retract(M, p, t[1] * X))
        for t in 0.1:0.1:0.9
            Test.@test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0
        end
    end
end


text:
@inline ngens(R::Nemo.Generic.MPolyRing) = R.num_vars

function _prod(A, b)
  for a in A
    b = b * a
  end
  return b
end

#aparently, should be called evaluate, talk to Bill...
#definitely non-optimal, in particular for automorphisms
function msubst(f::Generic.MPoly{T}, v::Array{NfRelElem{T}, 1}) where T
  k = base_ring(parent(f))
  n = length(v)
  @assert n == ngens(parent(f))
  r = zero(parent(v[1]))
  L = parent(v[1])
  for i=1:length(f)
    #@show prod(v[j]^f.exps[j, i] for j=1:n)
    s = _prod((v[j]^f.exps[n - j + 1, i] for j=1:n), one(L))
    r += f.coeffs[i]* s
  end
  return r
end
function msubst(f::Generic.MPoly{T}, v::Array{NfRelNSElem{T}, 1}) where T
  k = base_ring(parent(f))
  n = length(v)
  @assert n == ngens(parent(f))
  r = zero(k)
  for i=1:length(f)
    r += f.coeffs[i]*prod(v[j]^f.exps[n - j + 1, i] for j=1:n)
  end
  return r
end

################################################################################
#
#  Permutation group from NfRelNSToNfRelNSMor
#
################################################################################

function _get_poly_from_elem(a::NfRelNSElem{nf_elem}, Qxy)
  K = base_field(parent(a))
  Qx = parent(K.pol)
  p = change_base_ring(a.data, x -> evaluate(Qx(x), gen(Qxy, nvars(Qxy))))
  p1 = evaluate(p, [i for i in 1:ngens(parent(a))], gens(Qxy)[1:end-1])
  res = coeff(p1, [0 for i = 1:nvars(parent(p1))])
  return res
end

function multivariate_from_tower(f::Generic.MPoly{nf_elem}, Qxy)
  M = MPolyBuildCtx(Qxy)
  K = base_ring(f)
  Qx = parent(K.pol)
  cvzip = zip(coefficients(f), exponent_vectors(f))
  for (c, v) in cvzip
    pc = Qx(c)
    for i = degree(pc):-1:0
      cpc = coeff(pc, i)
      if iszero(cpc)
        continue
      end
      vn = vcat(v, i)
      push_term!(M, cpc, vn)
    end
  end
  return finish(M)::fmpq_mpoly
end

function (Rxy::NmodMPolyRing)(f::fmpq_mpoly)
  R = base_ring(Rxy)
  res = change_base_ring(f, x -> divexact(R(numerator(x)), R(denominator(x))), Rxy)
  return res
end

function _get_polys_from_auto(f::NfRelNSToNfRelNSMor_nf_elem, Qxy)
  res = Vector{elem_type(Qxy)}(undef, nvars(Qxy))
  if isdefined(f.image_data.base_field_map_data, :prim_image)
    # ap is a constant, but there is no easy way to coerce to the base field
    ap = f.image_data.base_field_map_data.prim_image.data.coeffs[1]
  else
    ap = gen(base_field(codomain(f)))
  end
  K = base_field(codomain(f))
  res[nvars(Qxy)] = evaluate(parent(K.pol)(ap), gen(Qxy, nvars(Qxy)))
  for i = 1:nvars(Qxy)-1
    res[i] = multivariate_from_tower(image_generators(f)[i].data, Qxy)
  end
  return res
end

function Base.hash(f::nmod_mpoly, h::UInt)
  return UInt(1)
end

function permutation_group1(G::Vector{NfRelNSToNfRelNSMor_nf_elem})
  L = domain(G[1])
  K = base_field(L)
  dK = absolute_degree(L)
  d1 = numerator(discriminant(L, FlintQQ))
  p = 2
  while divisible(d1, p)
    p = next_prime(p)
  end
  R = ResidueRing(FlintZZ, p, cached = false)
  Rm, gRm = PolynomialRing(R, ngens(L)+1, cached = false)
  fmod = Vector{nmod_mpoly}(undef, ngens(L)+1)
  RQm, gRQm = PolynomialRing(FlintQQ, ngens(L)+1, cached = false)
  p1 = K.pol
  p1Q = evaluate(p1, gRQm[end])
  fmod[1] = Rm(p1Q)
  for i = 1:ngens(L)
    pp = L.pol[i]
    pp1 = multivariate_from_tower(pp, RQm)
    fmod[i+1] = Rm(pp1)
  end
  permutations = Array{Array{Int, 1},1}(undef, length(G))
  for i = 1:length(G)
    permutations[i] = Vector{Int}(undef, dK)
  end
  pols = Vector{Vector{nmod_mpoly}}(undef, dK)
  pols[1] = gRm
  ind = 2
  gpols = nmod_mpoly[Rm(elel) for elel in _get_polys_from_auto(G[1], RQm)]
  if gpols != gRm
    pols[ind] = gpols
    ind += 1
    gpol = nmod_mpoly[compose_mod(gpols[i], [j for j = 1:nvars(Rm)], pols[2], fmod) for i = 1:length(gpols)]
    while gRm != gpol
      pols[ind] = gpol
      ind += 1
      gpol = nmod_mpoly[compose_mod(gpol[i], [j for j in 1:nvars(Rm)], pols[2], fmod) for i = 1:length(gpols)]
    end
  end
  for i in 2:length(G)
    pi = nmod_mpoly[Rm(x) for x in _get_polys_from_auto(G[i], RQm)]
    if !(pi in pols[1:ind-1])
      previous_order = ind - 1
      pols[ind] = pi
      ind += 1
      for j in 2:previous_order
        pols[ind] = nmod_mpoly[compose_mod(pols[j][s], Int[z for z in 1:nvars(Rm)], pi, fmod) for s = 1:length(pi)]
        ind += 1
      end
      if ind - 1 == dK
        break
      end
      rep_pos = previous_order + 1
      while rep_pos <= ind - 1
        for k in 1:i
          po = nmod_mpoly[Rm(elel) for elel in _get_polys_from_auto(G[k], RQm)]
          att = nmod_mpoly[compose_mod(pols[rep_pos][s], Int[i for i in 1:nvars(Rm)], po, fmod) for s = 1:length(pols[rep_pos])]
          if !(att in pols[1:ind-1])
            pols[ind] = att
            ind += 1
            for j in 2:previous_order
              pols[ind] = nmod_mpoly[compose_mod(pols[j][s], Int[z for z in 1:nvars(Rm)], att, fmod) for s = 1:length(pols[j])]
              ind += 1
            end
            if ind - 1 == dK
              break
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  #Now, I have the images mod p
  Dcreation = Vector{Tuple{Vector{nmod_mpoly}, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end

  gen_pols = Vector{nmod_mpoly}[nmod_mpoly[Rm(y) for y in _get_polys_from_auto(x, RQm)] for x in G]
  D = Dict(Dcreation)
  for s = 1:length(G)
    permutations[s][1] = D[gen_pols[s]]
    for i = 2:length(pols)
      permutations[s][i] = D[nmod_mpoly[compose_mod(gen_pols[s][t], Int[i for i in 1:nvars(Rm)], pols[i], fmod) for t = 1:length(gen_pols[s])]]
    end
  end
  return permutations
end

@doc Markdown.doc"""
    compose_mod(a::AbstractAlgebra.MPolyElem{T}, vars::Vector{Int}, vals::Vector{MPolyElem{T}}, mod::Vector{MPolyElem{T}}) where T <: FieldElement
Evaluate the polynomial by substituting in the supplied values in the array `vals` for
the corresponding variables with indices given by the array `vars`. The evaluation will
succeed if multiplication is defined between elements of the coefficient ring of $a$ and
elements of `vals`. The result will be reduced modulo "mod". If "mod" is a Groebner basis for the ideal
the elements generate.
"""
function compose_mod(a::S, vars::Vector{Int}, vals::Vector{S}, mod::Vector{S}) where S <:MPolyElem{T} where T <: RingElem
  unique(vars) != vars && error("Variables not unique")
  length(vars) != length(vals) && error("Number of variables does not match number of values")
  for i = 1:length(vars)
    if vars[i] < 1 || vars[i] > nvars(parent(a))
      error("Variable index not in range")
    end
  end
  if length(vars) == 0
    return a
  end
  powers = Dict{Int, S}[Dict{Int, S}() for i in 1:length(vals)]
  return _compose_mod(a, vars, vals, powers, mod)::S
end

function powermod(a::S, i::Union{Int, fmpz}, modu::Vector{S}) where S <:MPolyElem{T} where T <: RingElem
  if i == 0
    return one(parent(a))
  end
  if i == 1
    b = divrem(a, modu)[2]
    return b
  end
  if mod(i, 2) == 0
    j = div(i, 2)
    b = powermod(a, j, modu)
    b = b*b
    b = divrem(b, modu)[2]
    return b
  end
  b = divrem(a * powermod(a, i - 1, modu), modu)[2]
  return b
end

function mulmod(a::S, b::S, mod::Vector{S}) where S <:MPolyElem{T} where T <: RingElem
  return divrem(a*b, mod)[2]
end

function _compose_mod(a, vars, vals, powers, modu)
  S = parent(a)
  r = AbstractAlgebra.Generic.geobucket(S)
  cvzip = zip(coefficients(a), exponent_vectors(a))
  for (c, v) in cvzip
    t = one(S)
    for j = 1:length(vars)
      varnum = vars[j]
      exp = v[varnum]
      if !haskey(powers[j], exp)
        powers[j][exp] = powermod(vals[j], exp, modu)
      end
      t = mulmod(t, powers[j][exp], modu)
      v[varnum] = 0
    end
    M = MPolyBuildCtx(S)
    push_term!(M, c, v)
    push!(r, mulmod(t, finish(M), modu))
  end
  return finish(r)
end

function change_base_ring(p::MPolyElem{T}, g, new_polynomial_ring) where {T <: RingElement}
  cvzip = zip(coefficients(p), exponent_vectors(p))
  M = MPolyBuildCtx(new_polynomial_ring)
  for (c, v) in cvzip
    res = g(c)
    if !iszero(res)
      push_term!(M, g(c), v)
    end
  end
  return finish(M)::elem_type(new_polynomial_ring)
end

function Base.:(*)(f::Hecke.NfToNfRel, g::Hecke.NfRelToNfRelNSMor_nf_elem)
  @assert codomain(f) === domain(g)
  return hom(domain(f), codomain(g), g(f(gen(domain(f)))))
end
#
#function hom(K::AnticNumberField, L::NfRelNS{nf_elem}, img_gen::NfRelNSElem{nf_elem})
#  return Hecke.NfToNfRelNSMor(K, L, img_gen)
#end
#
#function image(f::Hecke.NfToNfRelNSMor, a::nf_elem)
#  K = parent(a)
#  Qx = parent(K.pol)
#  return evaluate(Qx(a), f.img_gen)
#end
#
#function preimage(phi::Hecke.NfToNfRelNSMor, a::NfRelNSElem{nf_elem})
#  @assert isdefined(phi, :preimg_base_field) && isdefined(phi, :preimgs)
#  f = data(a)
#  K = codomain(phi)
#  k = base_field(K)
#  R = parent(k.pol)
#  g = map_coefficients(x -> evaluate(R(x), phi.preimg_base_field), f)
#  return evaluate(g, phi.preimgs)
#end
#
#
function degrees(L::NfRelNS)
  return Int[total_degree(x) for x in L.pol]
end

function automorphisms(L::NfRelNS{T}) where T
  c = get_special(L, :automorphisms)
  if c !== nothing
    return c
  end
  auts = _automorphisms(L)
  Hecke.set_special(L, :automorphisms => auts)
  return auts
end

function _automorphisms(L::NfRelNS{T}) where T
  Kx, _ = PolynomialRing(base_field(L), "x", cached = false)
  rts = Vector{elem_type(L)}[roots(to_univariate(Kx, x), L) for x in L.pol]
  auts = Vector{morphism_type(L)}(undef, prod(length(x) for x in rts))
  ind = 1
  it = cartesian_product_iterator([1:length(rts[i]) for i in 1:length(rts)], inplace = true)
  for i in it
    embs = NfRelNSElem{T}[rts[j][i[j]] for j = 1:length(rts)]
    auts[ind] = hom(L, L, embs)
    ind += 1
  end
  return auts
end


text:
using Distributions
using DataFrames

include("utils/mctools.jl")

function getsim()
    mcs = @defsim begin

        ############################################################################
        # Define random variables (RVs)
        ############################################################################

        # The folllowing RVs are in more than one component.  For clarity they are
        # set here as opposed to below within the blocks of RVs separated by component
        # so that they are not set more than once.

        save_savingsrate = TriangularDist(10, 20, 15) # components: MarketDamages, MarketDamagesBurke, NonMarketDamages. GDP, SLRDamages
        tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.) # components: MarketDamages, NonMarketDamages, SLRDamages, Discountinuity

        wincf_weightsfactor_sea["USA"] = TriangularDist(.6, 1, .8) # components: SLRDamages, Discountinuity (weights for market and nonmarket are non-stochastic and uniformly 1)
        wincf_weightsfactor_sea["OECD"] = TriangularDist(.4, 1.2, .8)
        wincf_weightsfactor_sea["USSR"] = TriangularDist(.2, .6, .4)
        wincf_weightsfactor_sea["China"] = TriangularDist(.4, 1.2, .8)
        wincf_weightsfactor_sea["SEAsia"] = TriangularDist(.4, 1.2, .8)
        wincf_weightsfactor_sea["Africa"] = TriangularDist(.4, .8, .6)
        wincf_weightsfactor_sea["LatAmerica"] = TriangularDist(.4, .8, .6)

        automult_autonomoustechchange = TriangularDist(0.5, 0.8, 0.65)  # components: AdaptationCosts, AbatementCosts

        # The following RVs are divided into blocks by component

        # CO2cycle
        air_CO2fractioninatm = TriangularDist(57, 67, 62)
        res_CO2atmlifetime = TriangularDist(50, 100, 70)
        # ccf_CO2feedback = TriangularDist(0, 0, 0) # only usable if lb <> ub
        ccfmax_maxCO2feedback = TriangularDist(10, 30, 20)
        stay_fractionCO2emissionsinatm = TriangularDist(0.25, 0.35, 0.3)
        ce_0_basecumCO2emissions = TriangularDist(1830000, 2240000, 2040000)
        a1_percentco2oceanlong = TriangularDist(4.3,	41.6, 23.0)
        a2_percentco2oceanshort = TriangularDist(23.1, 30.1, 26.6)
        a3_percentco2land = TriangularDist(11.4, 42.5, 27.0)
        t1_timeco2oceanlong = TriangularDist(248.9, 376.2, 312.5)
        t2_timeco2oceanshort = TriangularDist(25.9, 43.9, 34.9)
        t3_timeco2land = TriangularDist(2.8, 5.7, 4.3)
        rt_g0_baseglobaltemp = TriangularDist(0.903, 0.989, 0.946)

        # SiBCASA Permafrost
        perm_sib_af = TriangularDist(1.42609149897258, 2.32504747848815, 1.87556948873036)
        perm_sib_sens_c_co2 = TriangularDist(28191.1555428869, 35688.3253432574, 31939.7404430722)
        perm_sib_lag_c_co2 = TriangularDist(35.4926669856915, 87.8949041341782, 61.6937855599349)
        perm_sib_pow_c_co2 = TriangularDist(0.107020247715729, 0.410961185142816, 0.258990716429273)
        perm_sib_sens_c_ch4 = TriangularDist(1240.3553299183, 3348.11995329232, 2294.23764160531)
        perm_sib_lag_c_ch4 = TriangularDist(75.1943160023131, 337.382510123922, 206.288413063117)
        perm_sib_pow_c_ch4 = TriangularDist(-0.108779283732708, 0.610889007954489, 0.25105486211089)

        # JULES Permafrost
        perm_jul_af = TriangularDist(1.70960411816136, 2.16221162526313, 1.93590787171224)
        perm_jul_sens_c_co2 = TriangularDist(24726.8035695649, 99008.7553497378, 61867.7794596514)
        perm_jul_lag_c_co2 = TriangularDist(252.558368389676, 834.674343162273, 543.616355775975)
        perm_jul_pow_c_co2 = TriangularDist(-0.226045987062471, 1.14010750072118, 0.457030756829357)
        perm_jul_ch4_co2_c_ratio = TriangularDist(2.77492291880781, 9.52902519167579, 6.04453870625663)

        # SulphateForcing
        d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
        ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)

        # ClimateTemperature
        frt_warminghalflife = TriangularDist(10, 55, 20)        # from PAGE-ICE v6.2 documentation
        tcr_transientresponse = TriangularDist(0.8, 2.7, 1.8)   # from PAGE-ICE v6.2 documentation
        alb_emulator_rand = TriangularDist(-1., 1., 0.)
        ampf_amplification["EU"] = TriangularDist(1.05, 1.53, 1.23)
        ampf_amplification["USA"] = TriangularDist(1.16, 1.54, 1.32)
        ampf_amplification["OECD"] = TriangularDist(1.14, 1.31, 1.21)
        ampf_amplification["USSR"] = TriangularDist(1.41, 1.9, 1.64)
        ampf_amplification["China"] = TriangularDist(1, 1.3, 1.21)
        ampf_amplification["SEAsia"] = TriangularDist(0.84, 1.15, 1.04)
        ampf_amplification["Africa"] = TriangularDist(0.99, 1.42, 1.22)
        ampf_amplification["LatAmerica"] = TriangularDist(0.9, 1.18, 1.04)

        # SeaLevelRise
        s0_initialSL = TriangularDist(0.17, 0.21, 0.19)                             # taken from PAGE-ICE v6.20 default
        sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)                            # median sensitivity to GMST changes
        sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)                           # asymptote for pre-industrial
        sltau_SLresponsetime = Gamma(16.0833333333333333, 24.)                      # fat-tailed distribution of time constant T_sl, sea level response time, from mode=362, mean = 386

        # GDP
        isat0_initialimpactfxnsaturation = TriangularDist(15, 25, 20)

        # MarketDamages
        iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
        W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
        pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
        ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)

        # MarketDamagesBurke
        impf_coeff_lin = TriangularDist(-0.0139791885347898, -0.0026206307945989, -0.00829990966469437)
        impf_coeff_quadr = TriangularDist(-0.000599999506482576, -0.000400007300924579, -0.000500003403703578)

        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["EU"] = TriangularDist(6.76231496767033, 13.482086163781, 10.1222005657257)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["USA"] = TriangularDist(9.54210085883826, 17.3151395362191, 13.4286201975287)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["OECD"] = TriangularDist(9.07596053028087, 15.0507477943984, 12.0633541623396)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["USSR"] = TriangularDist(3.01320548016903, 11.2132204366259, 7.11321295839747)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["China"] = TriangularDist(12.2330402806912, 17.7928749427573, 15.0129576117242)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["SEAsia"] = TriangularDist(23.3863348263352, 26.5136231383473, 24.9499789823412)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["Africa"] = TriangularDist(20.1866940491107, 23.5978086497453, 21.892251349428)
        MarketDamagesBurke_rtl_abs_0_realizedabstemperature["LatAmerica"] = TriangularDist(19.4846849750102, 22.7561130637973, 21.1203990194037)

        # NonMarketDamages
        tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.)
        iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
        w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
        pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
        ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)

        # SLRDamages
        scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
        # iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
        W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
        pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
        ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)

        # Discontinuity
        rand_discontinuity = Uniform(0, 1)
        tdis_tolerabilitydisc = TriangularDist(1, 2, 1.5)
        pdis_probability = TriangularDist(10, 30, 20)
        wdis_gdplostdisc = TriangularDist(1, 5, 3)
        ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
        distau_discontinuityexponent = TriangularDist(10, 30, 20)

        # EquityWeighting
        civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
        ptp_timepreference = TriangularDist(0.1, 2, 1)
        emuc_utilityconvexity = TriangularDist(0.5, 2, 1)

        # AbatementCosts
        AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50, 6.0, -22)
        AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-67, 6.0, -30)
        AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-20, 6.0, -7.0)
        AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50, 50, 0)

        AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 40, 20)
        AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)
        AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)
        AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)

        AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400, -100, -200)
        AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000, -1000, -4000)
        AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000, 0, -7000)
        AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400, -100, -200)

        AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60, 80, 70)
        AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35, 70, 50)
        AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35, 70, 50)
        AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60, 80, 70)

        AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100, 700, 400)
        AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000, 10000, 6000)
        AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000, 60000, 20000)
        AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100, 600, 300)

        AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000, 200000, 150000)
        AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500, 2500, 2000)
        AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30, 80, 50)
        AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500, 2500, 2000)

        # the following variables need to be set, but set the same in all 4 abatement cost components
        # note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables

        emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8, 1.2, 1.0)
        emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8, 1.2, 1.0)
        emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65, 1.35, 1.0)
        emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5, 1.5, 1.0)
        emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5, 1.5, 1.0)
        emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5, 1.5, 1.0)
        emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5, 1.5, 1.0)

        q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75, 1.5, 1.0)
        q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75, 1.25, 1.0)
        q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4, 1.0, 0.7)
        q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4, 1.0, 0.7)
        q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4, 1.0, 0.7)
        q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4, 1.0, 0.7)
        q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4, 1.0, 0.7)

        cmaxf_maxcostfactor["USA"] = TriangularDist(0.8, 1.2, 1.0)
        cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0, 1.5, 1.2)
        cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4, 1.0, 0.7)
        cmaxf_maxcostfactor["China"] = TriangularDist(0.8, 1.2, 1.0)
        cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1, 1.5, 1.2)
        cmaxf_maxcostfactor["Africa"] = TriangularDist(1, 1.5, 1.2)
        cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4, 1.0, 0.7)

        q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3, 1.2, 0.7)
        qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1, 1.5, 1.3)
        c0mult_mostnegativecostinfinalyear = TriangularDist(0.5, 1.2, 0.8)
        curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25, 0.8, 0.45)
        curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1, 0.7, 0.4)
        cross_experiencecrossoverratio = TriangularDist(0.1, 0.3, 0.2)
        learn_learningrate = TriangularDist(0.05, 0.35, 0.2)

        # AdaptationCosts
        AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
        AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
        AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
        AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
        AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
        AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)

        cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
        cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
        cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
        cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
        cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
        cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
        cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)

        ############################################################################
        # Indicate which parameters to save for each model run
        ############################################################################

        save(EquityWeighting.td_totaldiscountedimpacts,
             EquityWeighting.tpc_totalaggregatedcosts,
             EquityWeighting.tac_totaladaptationcosts,
             EquityWeighting.te_totaleffect,
             CO2Cycle.c_CO2concentration,
             TotalForcing.ft_totalforcing,
             ClimateTemperature.rt_g_globaltemperature,
             SeaLevelRise.s_sealevel,
             SLRDamages.rgdp_per_cap_SLRRemainGDP,
             MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP,
             NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
             Discontinuity.rgdp_per_cap_NonMarketRemainGDP)

    end # de
    return mcs
end

# Reformat the RV results into the format used for testing
function reformat_RV_outputs(samplesize::Int; output_path::String=joinpath(@__DIR__, "../output"))

    # create vectors to hold results of Monte Carlo runs
    td = zeros(samplesize);
    tpc = zeros(samplesize);
    tac = zeros(samplesize);
    te = zeros(samplesize);
    ft = zeros(samplesize);
    rt_g = zeros(samplesize);
    s = zeros(samplesize);
    c_co2concentration = zeros(samplesize);
    rgdppercap_slr = zeros(samplesize);
    rgdppercap_market = zeros(samplesize);
    rgdppercap_nonmarket = zeros(samplesize);
    rgdppercap_disc = zeros(samplesize);

    # load raw data
    # no filter
    td      = load_RV("EquityWeighting_td_totaldiscountedimpacts", "td_totaldiscountedimpacts"; output_path=output_path)
    tpc     = load_RV("EquityWeighting_tpc_totalaggregatedcosts", "tpc_totalaggregatedcosts"; output_path=output_path)
    tac     = load_RV("EquityWeighting_tac_totaladaptationcosts", "tac_totaladaptationcosts"; output_path=output_path)
    te      = load_RV("EquityWeighting_te_totaleffect", "te_totaleffect"; output_path=output_path)

    # time index
    c_co2concentration = load_RV("CO2Cycle_c_CO2concentration", "c_CO2concentration"; output_path=output_path)
    ft      = load_RV("TotalForcing_ft_totalforcing", "ft_totalforcing"; output_path=output_path)
    rt_g    = load_RV("ClimateTemperature_rt_g_globaltemperature", "rt_g_globaltemperature"; output_path=output_path)
    s       = load_RV("SeaLevelRise_s_sealevel", "s_sealevel"; output_path=output_path)

    # region index
    rgdppercap_slr          = load_RV("SLRDamages_rgdp_per_cap_SLRRemainGDP", "rgdp_per_cap_SLRRemainGDP"; output_path=output_path)
    rgdppercap_market       = load_RV("MarketDamagesBurke_rgdp_per_cap_MarketRemainGDP", "rgdp_per_cap_MarketRemainGDP"; output_path=output_path)
    rgdppercap_nonmarket    = load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path=output_path)
    rgdppercap_disc         = load_RV("Discontinuity_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path=output_path)

    # resave data
    df = DataFrame(td=td, tpc=tpc, tac=tac, te=te, c_co2concentration=c_co2concentration, ft=ft, rt_g=rt_g, sealevel=s, rgdppercap_slr=rgdppercap_slr, rgdppercap_market=rgdppercap_market, rgdppercap_nonmarket=rgdppercap_nonmarket, rgdppercap_di=rgdppercap_disc)
    save(joinpath(output_path, "mimipagemontecarlooutput.csv"), df)
end

function do_monte_carlo_runs(samplesize::Int, scenario::String="RCP4.5 & SSP2", output_path::String=joinpath(@__DIR__, "../output"))
    # get simulation
    mcs = getsim()

    # get a model
    m = getpage(scenario)
    run(m)

    # Run
    res = run(mcs, m, samplesize; trials_output_filename=joinpath(output_path, "trialdata.csv"), results_output_dir=output_path)

    # reformat outputs for testing and analysis
    reformat_RV_outputs(samplesize, output_path=output_path)
end


text:
module h2o
  #Right now, development is progressing assuming H2O cluster is started externally

  using JSON
  import Requests: get, put, post, readall

  #global vars
  export H2Oip, H2Oport

  #types
  export ClusterInfo

  #functions
  #not exporting functions by convention at this point

  #Type files
  include("types/clusterinfo.jl")

  #Function files
  include("cluster.jl")
  include("helpers.jl")
  include("ingestion.jl")
  include("info.jl")

end # module


text:
module BaseTestMulti

export MultiTestSet
export @mtestset
export @mtest_values_vary, @mtest_values_are, @mtest_values_include, @mtest_that_sometimes, @mtest_distributed_as

import Base.Test: record, finish
using Base.Test: AbstractTestSet, DefaultTestSet, Result, Pass, Fail, Error, ExecutionResult, Returned, Threw, get_testset
using HypothesisTests
using Distributions

# wrap DefaultTestSet with additional fields for samples taken at mtest macros
# we need DefaulTestSet since it handles pretty printing and accummulation of results
type MultiTestSet <: AbstractTestSet
	defaultts::DefaultTestSet
	mtests::Dict{Symbol, Function}
	samples::Dict{Symbol, Vector}
	repetition_count::Int
	repetitions::Int
	alpha::Float64
	function MultiTestSet(desc; reps::Int=30, alpha::Float64=0.01)
		reps >= 1 || error("reps must be at least 1")
		0. <= alpha <= 1. || error("alpha must be between 0 and 1")
		new(DefaultTestSet(desc), Dict{Symbol, Function}(), Dict{Symbol, Vector}(), 0, reps, alpha)
	end
end

# for normal tests, record against the wrapped default test set
record(ts::MultiTestSet, child::AbstractTestSet) = record(ts.defaultts, child)
record(ts::MultiTestSet, res::Result) = record(ts.defaultts, res)

# at end of multitest set, evaluate samples against corresponding mtest closures and record results in the wrapped default test set,
# then finish the wrapped test set so that results accumulate appropriately
function finish(ts::MultiTestSet)
	for (id, mtestclosure) in ts.mtests
		s = get(ts.samples, id, [])
		res = try
				result, macrosym, evaluand, params = mtestclosure(s, ts.alpha)
				testdesc = string(macrosym) * (isempty(params) ? "" : " " * params) * " " * evaluand * "\n      Sample: " * string(s)
				# this description (macrosym evaluand & any parameters) will be reported by DefaultTestSet after the label "Expression:"
				# the spaces after the \n aligns the label "Sample:" on the second line to this
				result ? Pass(macrosym, testdesc, nothing, nothing) : Fail(macrosym, testdesc, nothing, nothing)
			catch _e
				Error(:mtest, nothing, _e, catch_backtrace())
			end
		record(ts.defaultts, res)
	end
	finish(ts.defaultts)
end

# determine whether to continue with more repetitions
# this could be extend with accumulator-type functionality from AutoTest so that number of repeats (sample size) is adjusted automatically
function continue_sampling()
	ts = get_testset()
	ts.repetition_count += 1
	ts.repetition_count <= ts.repetitions
end

# syntatic sugar for @testset MultiTestSet, also insert a while loop around test set block that is controlled automatically
macro mtestset(args...)
	testsetex = args[end]
	@assert isa(testsetex, Expr) && (testsetex.head in (:block, :for))
	testsetblock = testsetex.head == :block ? testsetex : testsetex.args[2]
	@assert isa(testsetblock, Expr) && (testsetblock.head == :block)
	mtestloopex = quote
		while BaseTestMulti.continue_sampling()
			$(testsetblock)
		end
	end
	if testsetex.head == :block
		testsetex = mtestloopex
	else
		testsetex.args[2] = mtestloopex
	end
	newargs = (:MultiTestSet, args[1:end-1]..., testsetex)
	mc = :( @testset($(newargs...)) )
	esc(mc)
end

# register a mtest macro first time it is executed, to record both id and the test function (as a closure including any parameters
# such as expected results)
register_mtest(ts::MultiTestSet, id::Symbol, mtestclosure::Function) = ts.mtests[id] = mtestclosure

is_mtest_registered(ts::MultiTestSet, id::Symbol) = haskey(ts.mtests, id)

# add result to sample (or report error during test)
function add_to_mtest_sample(ts::MultiTestSet, id::Symbol, result::ExecutionResult, origex)
    if isa(result, Returned)
		s = get!(Vector{Any}, ts.samples, id)
		push!(s, result.value)
    else
        # The predicate couldn't be evaluated without throwing an
        # exception, so that is an Error and not a Fail
        @assert isa(result, Threw)
        testres = Error(:test_error, origex, result.exception, result.backtrace)
	    record(ts, testres)
    end
end

# simplified form of Base.Test.get_test_result: special handling for comparisons is removed as we don't need it (but see mtest_that_sometimes)
function get_mtest_result_expr(ex)
    testret = :(Returned($(esc(ex)), nothing))
    resultex = quote
        try
            $testret
        catch _e
            Threw(_e, catch_backtrace())
        end
    end
    Base.remove_linenums!(resultex)
    resultex
end

# common code for all mtest macros: output is code to register mtest instance and its associated closure (if first time executed),
# and then add value to the sample
function mtest_macro(ex, paramex, mtestclosureex, macrosym)
    origex = Expr(:inert, ex)
    resultex = get_mtest_result_expr(ex)
	id = gensym(:mtest)
	idstr = string(id)
	tsvar = gensym(:ts)
    quote
		$tsvar = get_testset()
		isa($tsvar, MultiTestSet) || error("To apply the @" * $(string(macrosym)) * " macro, replace @testset with @mtestset")
		if !is_mtest_registered($tsvar, Symbol($idstr))
			$paramex
			register_mtest($tsvar, Symbol($idstr), $mtestclosureex)
		end
		add_to_mtest_sample($tsvar, Symbol($idstr), $resultex, $origex)
	end
end

# mtest_values_vary
macro mtest_values_vary(ex)
	macrosym = :mtest_values_vary
	paramex = :( nothing )
	mtestclosureex = :( (_s,_a)->(length(unique(_s))>1, Symbol($(string(macrosym))), $(string(ex)), "") )
	mtest_macro(ex, paramex, mtestclosureex, macrosym)
end

# mtest_values_are
macro mtest_values_are(expex, ex)
	macrosym = :mtest_values_are
	expvar = gensym(:exp)
	paramex = :( $expvar = $(esc(expex)) )
	mtestclosureex = :( (_s,_a)->(sort(unique(_s))==sort(unique($expvar)), Symbol($(string(macrosym))), $(string(ex)), string($expvar)) )
	mtest_macro(ex, paramex, mtestclosureex, macrosym)
end

# mtest_values_include
macro mtest_values_include(expex, ex)
	macrosym = :mtest_values_include
	expvar = gensym(:exp)
	paramex = :( $expvar = $(esc(expex)) )
	mtestclosureex = :( (_s,_a)->(issubset($expvar, _s), Symbol($(string(macrosym))), $(string(ex)), string($expvar)) )
	mtest_macro(ex, paramex, mtestclosureex, macrosym)
end

# mtest_that_sometimes
macro mtest_that_sometimes(ex)
	macrosym = :mtest_that_sometimes
	paramex = :( nothing )
	mtestclosureex = :( (_s,_a)->(any(_s), Symbol($(string(macrosym))), $(string(ex)), "") )
	mtest_macro(ex, paramex, mtestclosureex, macrosym)
end

# mtest_distributed_as
# comparison need not be a Distribution - needs only to permit rand( ,n); a vector would work
# test is Mann Whitney Ranksum against a sample of the same length with significance level alpha
# NOTE: THIS IS NOT IDEAL:  Ranksum is NOT a test of same distributions, only of stochastic equivalence, i.e. if a sample is drawn from each,
# the chance that the first is greater than the second is 0.5.
# Ideally would like to use a test such as Kolmogorov-Smirnov, but KS applies only to continuous distribution and is unrealiable when there are
# many tied values (this is seen in practice with ExactOneSampleKSTest and ApproximateTwoSampleKSTest from Distributions.jl)
macro mtest_distributed_as(distex, ex)
	macrosym = :mtest_distributed_as
	distvar = gensym(:dist)
	paramex = :( $distvar = $(esc(distex)) )
	mtestclosureex = :( (_s,_a)->(pvalue(MannWhitneyUTest(convert(Vector{Real},_s),rand($distvar, length(_s)))) > _a,
		 					Symbol($(string(macrosym))), $(string(ex)), string($distvar) ) )
	# mtestclosureex = :( (_s,_a)->(pvalue(isa($distvar, UnivariateDistribution) ?
	# 								ExactOneSampleKSTest(convert(Vector{Real},_s), $distvar) :
	# 								ApproximateTwoSampleKSTest(convert(Vector{Real},_s),rand($distvar, length(_s)))) > _a,
	# 	 					Symbol($(string(macrosym))), $(string(ex)), string($distvar) ) )
	mtest_macro(ex, paramex, mtestclosureex, macrosym)
end

end # module


text:
# Multidimensional arrays
z = zeros(Float64, 2, 3)
println(typeof(z))

# Declare array of dimension n x m
n, m = 2, 4
arr = Array{Int}(undef, 2, 4)
println(arr)
println(size(arr))
arr2 = Array{Int}(undef, 3, 2, 2)
println(arr2)
println(size(arr2))

s = ones(String, 1, 3)
println(s)
# Note that s is considered a "row matrix / vector".
# The typical array is considered a "column matrix / vector".
println(s ==  ["", "", ""])

a = [1 2 3; 4 5 6]
println(a) # 2 rows, 3 columns
println(a[1, 2])
a[2, 3] = 10
println(a)

u = a[:,2:end]
println(ℯ.^(im*u))


text:
using Nextion
using Test

@testset "checkbox" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexCheckbox0 = NexCheckbox(nexSerial, Name("c0"), cid=ComponentID(4))
    nexCheckbox1 = NexCheckbox(nexSerial, Name("c1"), cid=ComponentID(5))
    nexCheckbox2 = NexCheckbox(nexSerial, Name("c2"), cid=ComponentID(7))

    send(nexSerial, "page pg_dsb_chk_rad")

    nexCheckbox0.value = false
    nexCheckbox1.value = false
    nexCheckbox2.value = false

    @test !nexCheckbox0.value

    sleep(1)

    nexCheckbox0.value = true

    @test nexCheckbox0.value

    sleep(1)

    nexCheckbox0.value = false

    @test !nexCheckbox0.value

    sleep(1)

    close(nexSerial)

    @test true

end


text:
using NamedDims
using NamedDims: names
using Test
using Tracker

@testset "Binary broadcasting operations (.+)" begin
    nda = NamedDimsArray{(:a,)}(ones(3))

    @testset "standard case" begin
        @test nda .+ nda == 2ones(3)
        @test dimnames(nda .+ nda) == (:a,)

        @test nda .+ nda .+ nda == 3ones(3)
        @test dimnames(nda .+ nda .+ nda) == (:a,)

        # in-place
        @test dimnames(nda .= 0 .* nda .+ 7) == (:a,)
        @test unname(nda .= 0 .* nda .+ 7) == 7*ones(3)
    end

    @testset "partially named dims" begin
        ndx = NamedDimsArray{(:x, :_)}(ones(3, 5))
        ndy = NamedDimsArray{(:_, :y)}(ones(3, 5))

        lhs = ndx .+ ndy
        rhs = ndy .+ ndx
        @test dimnames(lhs) == (:x, :y) == dimnames(rhs)
        @test lhs == 2ones(3, 5) == rhs
    end

    @testset "Dimension disagreement" begin
        @test_throws DimensionMismatch .+(
            NamedDimsArray{(:a, :b, :c, :d)}(zeros(3, 3, 3, 3)),
            NamedDimsArray{(:w, :x, :y, :z)}(ones(3, 3, 3, 3))
        )
    end

    @testset "named and unnamed" begin
        lhs_sum = .+(
            NamedDimsArray{(:a,:b,:c,:d)}(zeros(3, 3, 3, 3)),
            ones(3, 3, 3, 3)
        )
        @test lhs_sum == ones(3, 3, 3, 3)
        @test dimnames(lhs_sum) == (:a, :b, :c, :d)

        rhs_sum = .+(
            zeros(3, 3, 3, 3),
            NamedDimsArray{(:w, :x, :y, :z)}(ones(3, 3, 3, 3))
        )
        @test rhs_sum == ones(3, 3, 3, 3)
        @test dimnames(rhs_sum) == (:w, :x, :y, :z)
    end

    @testset "broadcasting" begin
        v = NamedDimsArray{(:time,)}(zeros(3,))
        m = NamedDimsArray{(:time, :value)}(ones(3, 3))
        s = 0

        @test v .+ m == ones(3, 3) == m .+ v
        @test s .+ m == ones(3, 3) == m .+ s
        @test s .+ v .+ m == ones(3, 3) == m .+ s .+ v

        @test dimnames(v .+ m) == (:time, :value) == dimnames(m .+ v)
        @test dimnames(s .+ m) == (:time, :value) == dimnames(m .+ s)
        @test dimnames(s .+ v .+ m) == (:time, :value) == dimnames(m .+ s .+ v)
    end

    @testset "Mixed array types" begin
        casts = (
            NamedDimsArray{(:foo, :bar)},  # Named Matrix
            x->NamedDimsArray{(:foo,)}(x[:, 1]),  # Named Vector
            x->NamedDimsArray{(:foo, :bar)}(x[:, 1:1]),  # Named Single Column Matrix
            identity, # Matrix
            x->x[:, 1], # Vector
            x->x[:, 1:1], # Single Column Matrix
            first, # Scalar
         )
        for (T1, T2, T3) in Iterators.product(casts, casts, casts)
            all(isequal(identity), (T1, T2, T3)) && continue
            !any(isequal(NamedDimsArray{(:foo, :bar)}), (T1, T2, T3)) && continue

            total = T1(ones(3, 6)) .+ T2(2ones(3, 6)) .+ T3(3ones(3, 6))
            @test total == 6ones(3, 6)
            @test dimnames(total) == (:foo, :bar)
        end
    end

    @testset "Regression test again #8b" begin
        # https://github.com/invenia/NamedDims.jl/issues/8#issuecomment-490124369
        nda = NamedDimsArray{(:x,:y,:z)}(ones(10,20,30))
        @test nda .+ ones(1,20) == 2ones(10,20,30)
        @test dimnames(nda .+ ones(1,20)) == (:x, :y, :z)
    end

    @testset "in-place assignment .=" begin
        ab = NamedDimsArray(rand(2,2), (:a, :b))
        a_ = NamedDimsArray(rand(2,2), (:a, :_))
        ba = NamedDimsArray(rand(2,2), (:b, :a))
        ac = NamedDimsArray(rand(2,2), (:a, :c))
        z = zeros(2,2)

        # https://github.com/invenia/NamedDims.jl/issues/71
        @test_throws DimensionMismatch z .= ab .+ ba
        @test_throws DimensionMismatch z .= ab .+ ac
        @test_throws DimensionMismatch a_ .= ab .+ ac
        @test_throws DimensionMismatch ab .= a_ .+ ac
        @test_throws DimensionMismatch ac .= ab .+ ba

        # check that dest is written into:
        @test dimnames(z .= ab .+ ba') == (:a, :b)
        @test z == (ab.data .+ ba.data')
        @test z isa Array  # has not itself magically gained names

        @test dimnames(z .= ab .+ a_) == (:a, :b)
        @test dimnames(a_ .= ba' .+ ab) == (:a, :b)

        # https://github.com/invenia/NamedDims.jl/issues/175
        nda1 = NamedDimsArray(zeros(2, 2), (:a, :b))
        nda2 = NamedDimsArray([1.0 2.0;], (:a, :b))
        @test (nda1 .= nda2) == (parent(nda1) .= parent(nda2))
    end
end

@testset "Competing Wrappers" begin
    nda = NamedDimsArray(ones(4), :foo)
    ta = TrackedArray(5*ones(4))
    ndt = NamedDimsArray(TrackedArray(5*ones(4)), :foo)

    arrays = (nda, ta, ndt)
    @testset "$a .- $b" for (a, b) in Iterators.product(arrays, arrays)
        a === b && continue
        @test typeof(nda .- ta) <: NamedDimsArray
        @test typeof(parent(nda .- ta)) <: TrackedArray
    end
end

@testset "Broadcasting against Tuples" begin
    # https://github.com/invenia/NamedDims.jl/issues/140
    nda = NamedDimsArray(zeros(4), :foo)
    @test (nda .+ (1,2,3,4)) isa NamedDimsArray{(:foo,)}
    @test (nda .+ (1,2,3,4)) == [1, 2, 3, 4]

    @test ((1,2,3,4) .+ nda) isa NamedDimsArray{(:foo,)}
    @test ((1,2,3,4) .+ nda) == [1, 2, 3, 4]
end

function foo(a, b, c, d)
    @. a -= 0.5 * d * b * c
    return a
end
@testset "Regression test against allocations in broadcasting #187" begin
    # https://github.com/invenia/NamedDims.jl/issues/187
    a = NamedDimsArray{(:z,)}(rand(5))
    b = NamedDimsArray{(:z,)}(rand(5))
    c = rand()
    d = rand()
    @test_modern 0 == @ballocated foo($a, $b, $c, $d)
end


text:
module ImageInTerminal

using Requires
using Crayons
using ImageCore
using ImageBase: restrict

export

    colorant2ansi,
    imshow,
    imshow256,
    imshow24bit

include("colorant2ansi.jl")
include("encodeimg.jl")
include("imshow.jl")

# -------------------------------------------------------------------
# overload default show in the REPL for colorant (arrays)

const colormode = TermColorDepth[TermColor256()]
const should_render_image = Bool[true]
const encoder_backend = [:ImageInTerminal]

"""
    use_256()

Triggers `imshow256` automatically if an array of colorants is to
be displayed in the julia REPL. (This is the default)
"""
use_256() = (colormode[1] = TermColor256(); should_render_image[1] = true)

"""
    use_24bit()

Triggers `imshow24bit` automatically if an array of colorants is to
be displayed in the julia REPL.
Call `ImageInTerminal.use_256()` to restore default behaviour.
"""
use_24bit() = (colormode[1] = TermColor24bit(); should_render_image[1] = true)

"""
    disable_encoding()

Disable the image encoding feature and show images as if they are normal arrays.

This can be restored by calling `ImageInTerminal.enable_encoding()`.
"""
disable_encoding() = (should_render_image[1] = false)

"""
    enable_encoding()

Enable the image encoding feature and show images in terminal.

This can be disabled by calling `ImageInTerminal.disable_encoding()`. To choose between
different encoding method, call `ImageInTerminal.use_256()` or `ImageInTerminal.use_24bit()`.
"""
enable_encoding() = (should_render_image[1] = true)

# colorant arrays
function Base.show(
        io::IO, mime::MIME"text/plain",
        img::AbstractArray{<:Colorant})
    if should_render_image[1]
        println(io, summary(img), ":")
        ImageInTerminal.imshow(io, img, colormode[1])
    else
        invoke(Base.show, Tuple{typeof(io), typeof(mime), AbstractArray}, io, mime, img)
    end
end

# colorant
function Base.show(io::IO, mime::MIME"text/plain", color::Colorant)
    if should_render_image[1]
        fgcol = _colorant2ansi(color, colormode[1])
        chr = _charof(alpha(color))
        print(io, Crayon(foreground = fgcol), chr, chr, " ")
        print(io, Crayon(foreground = :white), color)
        print(io, Crayon(reset = true))
    else
        invoke(Base.show, Tuple{typeof(io), typeof(mime), Any}, io, mime, color)
    end
end

function __init__()
    # use 24bit if the terminal supports it
    lowercase(get(ENV, "COLORTERM", "")) in ("24bit", "truecolor") && use_24bit()
    enable_encoding()

    if VERSION < v"1.6.0-DEV.888" && Sys.iswindows()
        # https://discourse.julialang.org/t/image-in-repl-does-not-correct/46359
        @warn "ImageInTerminal is not supported for Windows platform: Julia at least v1.6.0 is required."
        disable_encoding()
    end

    # Sixel requires Julia at least v1.6. We don't want to maintain an ImageInTerminal branch
    # for old Julia versions so here we use Requires to conditionally load Sixel as an advanced
    # image encoding choice. All ImageInTerminal functionality is still there even without Sixel
    # -- well, basically.
    @require Sixel="45858cf5-a6b0-47a3-bbea-62219f50df47" begin
        if Sixel.is_sixel_supported()
            encoder_backend[1] = :Sixel
        end
        sixel_encode(args...; kwargs...) = Sixel.sixel_encode(args...; kwargs...)
    end
end

end # module


text:
@testset "Merge" begin

    left = GorFile("left.gor")
    right = GorFile("right.gor")

    m = left |> GOR.merge(right)
    @test m |> collect |> length == 200

    @test eltype(m) == NamedTuple{(:Chrom, :Pos, :Val), Tuple{String, Int64, String}}

    using DataFrames

    df1 = Tables.namedtupleiterator(DataFrame(Chrom = [1,2,3], Pos = [1,2,3],
                                              val1 = [1,2,3]))
    df2 = Tables.namedtupleiterator(DataFrame(Chrom = [1,2,3], Pos = [1,2,3],
                                              val1 = [1.0, 2, missing],
                                              val2 = [1,2,3]))

    m2 = df1 |> GOR.merge(df2)

    @test eltype(m2) == NamedTuple{(:Chrom, :Pos, :val1, :val2),Tuple{Int64,Int64,Union{Missing, Float64},Union{Missing, Int64}}}

    @test m2 |> collect |> length == 6

    using Tables

    @test Tables.istable(m)
    @test Tables.rowaccess(m)
    @test Tables.schema(m).names == (:Chrom, :Pos, :Val)
    @test Tables.schema(m).types == (String, Int64, String)
    @test Tables.schema(m) == Tables.Schema{(:Chrom, :Pos, :Val), Tuple{String, Int64, String}}()

end


text:
"""
Welcome to EnvironmentMigrators.jl.

Use `EnvironmentMigrators.wizard()` for interactive use.

# High level interface (unexported)
* `list_shared_environments([depot])`
* `select_shared_environments([depot])`
* `backup_current_environment(; fileaction = cp)`
* `migrate_selected_environment(selected_environment; backup = true)`
* `migrate_current_environment_to(path)`
* `wizard([depot])`

# Exported types
* `SimpleEnvironmentMigrator`
* `SimpleBackupOnlyEnvironmentMigrator`

# Exported methods
* `migrate(m::EnvironmentMigrators.AbstractEnvironmentMigrator)`
* `backup(m::EnvironmentMigrators.AbstractEnvironmentMigrator)`
"""
module EnvironmentMigrators

# https://github.com/JuliaLang/julia/pull/34896
@static if isdefined(Base, :Experimental) &&
           isdefined(Base.Experimental, Symbol("@optlevel"))
    Base.Experimental.@optlevel 1
end

using Pkg
using REPL.TerminalMenus
using Dates

export SimpleEnvironmentMigrator, SimpleBackupOnlyEnvironmentMigrator
export migrate, backup

include("types.jl")

"""
    list_shared_environments([depot])

List the top level shared environments in $(joinpath(DEPOT_PATH[1], "environments")).
`depot` defaults to `first(DEPOT_PATH)`.
"""
function list_shared_environments(depot = first(DEPOT_PATH))
    shared_environments = joinpath(depot, "environments")
    if !isdir(shared_environments)
        return String[]
    else
        return readdir(shared_environments)
    end
end

"""
    select_shared_environments([depot])

Present an interactive menu to select a shared environment to copy from.
Also provide options to select the current enviornment for backup or quit.
`depot` defaults to `first(DEPOT_PATH)`.
"""
function select_shared_environments(depot = first(DEPOT_PATH))
    envs = list_shared_environments()
    options = copy(envs)
    for idx in eachindex(options)
        if options[idx] == "v$(VERSION.major).$(VERSION.minor)"
            options[idx] = options[idx] * "(current version)"
        end
    end
    push!(options, "Backup the current environment: $(Base.active_project()).")
    push!(options, "Quit. Do Nothing.")
    menu = RadioMenu(options)
    println("Use the arrow keys to move the cursor. Press enter to select.")
    println("Please select a shared environment to copy to $(Base.active_project()):")
    menu_idx = request(menu)
    if menu_idx <= length(envs)
        return joinpath(depot, "environments", envs[menu_idx])
    elseif menu_idx == length(envs) + 1
        return Base.active_project()
    elseif menu_idx == length(options)
        @info "Quiting. No action taken."
        return ""
    end
end

"""
    backup_current_environments([depot]; fileaction=cp)

Backup the Project.toml and Manifest.toml to `backups/timestamp` where
`timestamp = Dates.format(Dates.now(), "yyyy-mm-dd_HH_MM_SS")``.

`fileaction` is a keyword parameter for a function that accepts two parameters.
`fileaction` defaults to `cp` for copy. An alternative is `mv` for move.
`depot` defaults to `first(DEPOT_PATH)`.
"""
function backup_current_environment(; fileaction=cp)
    @info "Current environment:"
    m = SimpleBackupOnlyEnvironmentMigrator(fileaction)
    Pkg.activate(m.target_project_toml)
    Pkg.status()
    backup(m)
end

"""
    migrate_selected_environment(selected_env)

Migrate the selected environment specified as a directory in selected_env to
the current active project.

`selected_env` should be the absolute path to an environment _directory_.
"""
function migrate_selected_environment(selected_env; backup = true)
    if !isdir(selected_env) && isfile(selected_env) && endswith(selected_env, "Project.toml")
        selected_env = dirname(selected_env)
    end
    selected_project_toml = joinpath(selected_env, "Project.toml")
    selected_julia_project_toml = joinpath(selected_env, "JuliaProject.toml")
    if isfile(selected_julia_project_toml)
        selected_project_toml = selected_julia_project_toml
    end
    selected_manifest_toml = joinpath(selected_env, "Manifest.toml")
    selected_julia_manifest_toml = joinpath(selected_env, "JuliaManifest.toml")
    if isfile(selected_julia_manifest_toml)
        selected_manifest_toml = selected_julia_manifest_toml
    end
    m = SimpleEnvironmentMigrator(selected_project_toml, selected_manifest_toml)
    migrate(m; backup = backup)
    return nothing
end

"""
    wizard([depot])

Run an interactive terminal based wizard.

`depot` defaults to `first(DEPOT_PATH)`.
"""
function wizard(depot = first(DEPOT_PATH))
    selected_env = select_shared_environments(depot)
    if isempty(selected_env)
        return nothing
    end
    fileaction = if selected_env == Base.active_project()
        @info "Backing up the current environment" Base.active_project()
        cp
    else
        mv
    end
    backup_current_environment(; fileaction=fileaction)
    if selected_env != Base.active_project()
        migrate_selected_environment(selected_env; backup = false)
    end
    return nothing
end

"""
    migrate_current_environment_to(path::AbstractString)

Migrate the current environment to the indicated path.
"""
function migrate_current_environment_to(path::AbstractString; update = false)
    mkpath(path)
    source_project_toml, source_manifest_toml = get_current_project_and_manifest()
    target_project_toml, target_manifest_toml = joinpath.(path, basename.((source_project_toml, source_manifest_toml)))
    m = SimpleEnvironmentMigrator(source_project_toml, source_manifest_toml,
                                  target_project_toml, target_manifest_toml)
    migrate(m; update = update)
end

function __init__()
    @info "Thank you for using EnvironmentMigrators.jl. Please run `EnvironmentMigrators.wizard()` to begin."
end

end


text:
module AMRVW

using LinearAlgebra
import LinearAlgebra: eigvals, diagm
import Base: adjoint, *

# package code goes here

include("rotators.jl")
include("transformations.jl")
include("chains.jl")
include("diagonal.jl")

include("rfactorization.jl")
include("qfactorization.jl")
include("qrfactorization.jl")

include("diagonal-block.jl")
include("bulge-step.jl")
include("create-bulge.jl")

include("RDS.jl")
include("CSS.jl")
include("pencil.jl")
include("twisted.jl")

include("AMRVW_algorithm.jl")
include("companion.jl")

include("utils.jl")
include("diagnostics.jl")

end # module


text:
jltype2nctype={Int16=>NC_SHORT,
               Int32=>NC_INT,
               Int64=>NC_LONG,
               Float32=>NC_FLOAT,
               Float64=>NC_DOUBLE}

function _cchartostring(in)
  iname=1
  name=""
  #println(in)
  if (length(in)>0)
    while (char(in[iname])!='\0' && iname<length(in))
      name=string(name,char(in[iname]))
      iname=iname+1
    end
  end
  return name
end

function _nc_op(fname::String,omode::Uint16)
  # Open netcdf file
  ida=Array(Int32,1)
  NetCDF._nc_open_c(fname,omode,ida)
  id=ida[1]
  NC_VERBOSE ? println("Successfully opened ",fname," dimid=",id) : nothing
  return id
end

function _nc_inq_dim(id::Integer,idim::Integer)
  namea=Array(Uint8,NC_MAX_NAME+1);lengtha=Array(Int32,1)
  NetCDF._nc_inq_dim_c(id,idim,namea,lengtha)
  name=_cchartostring(namea)
  dimlen=lengtha[1]
  NC_VERBOSE ? println("Successfully read from file") : nothing
  NC_VERBOSE ? println("name=",name," dimlen=",dimlen) : nothing
  return (name,dimlen)
end

function _nc_inq_dimid(id::Integer,name::String)
  dimida=Array(Int32,1)
  try
    NetCDF._nc_inq_dimid_c(id,name,dimida)
  catch
    dimida[1]=-1
  end
  NC_VERBOSE ? println("Successfully read from file") : nothing
  return dimida[1]
end

function _ncf_inq(id::Integer)
  # Inquire number of codes in netCDF file
  ndima=Array(Int32,1);nvara=Array(Int32,1);ngatta=Array(Int32,1);nunlimdimida=Array(Int32,1)
  _nc_inq_c(id,ndima,nvara,ngatta,nunlimdimida)
  ndim=ndima[1]
  nvar=nvara[1]
  ngatt=ngatta[1]
  nunlimdimid=nunlimdimida[1]
  NetCDF.NC_VERBOSE ? println("Successfully read from file") : nothing
  NetCDF.NC_VERBOSE ? println("ndim=",ndim," nvar=",nvar," ngatt=",ngatt," numlimdimid=",nunlimdimid) : nothing
  return (ndim,nvar,ngatt,nunlimdimid)
end

function _nc_inq_attname(ncid::Integer,varid::Integer,attnum::Integer)
  # Get attribute name from attribute number
  namea=Array(Uint8,NC_MAX_NAME+1)
  _nc_inq_attname_c(ncid,varid,attnum,namea)
  name=_cchartostring(namea)
  NC_VERBOSE ? println("Successfully read attribute name") : nothing
  NC_VERBOSE ? println("name=",name) : nothing
  return name
end

function _nc_inq_att(ncid::Integer,varid::Integer,attnum::Integer)
  # First get attribute name
  name=_nc_inq_attname(ncid,varid,attnum)
  NC_VERBOSE ? println(name) : nothing
  #Then find out about attribute
  typea=Array(Int32,1);nvals=Array(Int32,1)
  _nc_inq_att_c(ncid,varid,name,typea,nvals)
  attype=typea[1]
  NC_VERBOSE ? println("Successfully read attribute type and number of vals") : nothing
  NC_VERBOSE ? println("atttype=",attype," nvals=",nvals[1]) : nothing
  text=_nc_get_att(ncid,varid,name,attype,nvals[1])
  return (name,text)
end

function _nc_put_att(ncid::Integer,varid::Integer,name,val)
  if (start(val)>0)
    if (typeof(val[1])<:Char)
      attlen=length(val)+1
      attype=NC_CHAR
    else
      attlen=length(val)
      attype=jltype2nctype[typeof(val[1])]
    end
  else
    if (typeof(val)<:Char)
      attlen=2
      attype=NC_CHAR
    else
      attlen=1
      attype=jltype2nctype[typeof(val)]
      val=[val]
    end
  end
  NC_VERBOSE ? println(ncid, varid, name, attype, attlen, val) : nothing
  if (attype==NC_CHAR)
    _nc_put_att_text_c(ncid,varid,name,attlen,val)
  elseif (attype==NC_SHORT)
    _nc_put_att_short_c(ncid,varid,name,attype,attlen,val)
  elseif (attype==NC_INT)
    _nc_put_att_int_c(ncid,varid,name,attype,attlen,int32(val))
  elseif (attype==NC_FLOAT)
    _nc_put_att_float_c(ncid,varid,name,attype,attlen,val)
  elseif (attype==NC_DOUBLE)
    _nc_put_att_float_c(ncid,varid,name,NC_FLOAT,attlen,float32(val))
  elseif (attype==NC_BYTE)
    _nc_put_att_ubyte_c(ncid,varid,name,attype,attlen,val)
  else
    valsa="Could not read attribute, currently unsupported datatype by the netcdf package"
  end
end

function _nc_get_att(ncid::Integer,varid::Integer,name,attype::Integer,attlen::Integer)
  if (attype==NC_CHAR)
    valsa=Array(Uint8,attlen+5)
    _nc_get_att_text_c(ncid,varid,name,valsa)
    valsa=string(_cchartostring(valsa))
  elseif (attype==NC_SHORT)
    valsa=Array(Int16,attlen)
    _nc_get_att_short_c(ncid,varid,name,valsa)
  elseif (attype==NC_INT)
    valsa=Array(Int32,attlen)
    _nc_get_att_int_c(ncid,varid,name,valsa)
  elseif (attype==NC_LONG)
    valsa=Array(Int64,attlen)
    _nc_get_att_long_c(ncid,varid,name,valsa)
  elseif (attype==NC_FLOAT)
    valsa=Array(Float32,attlen)
    _nc_get_att_float_c(ncid,varid,name,valsa)
  elseif (attype==NC_DOUBLE)
    valsa=Array(Float64,attlen)
    _nc_get_att_double_c(ncid,varid,name,valsa)
  elseif (attype==NC_BYTE)
    valsa=Array(Uint8,attlen)
    _nc_get_att_ubyte_c(ncid,varid,name,valsa)
  else
    valsa="Could not read attribute, currently unsupported datatype by the netcdf package"
  end
  return valsa
end

function _ncv_inq(nc::NcFile,varid::Integer)
  id=nc.ncid
  ndim=length(nc.dim)
  # Inquire variables in the file
  namea=Array(Uint8,NC_MAX_NAME+1);xtypea=Array(Int32,1);ndimsa=Array(Int32,1);dimida=Array(Int32,ndim);natta=Array(Int32,1)
  _nc_inq_var_c(id,varid,namea,xtypea,ndimsa,dimida,natta)
  NC_VERBOSE ? println("dimida=",dimida," ndimsa=",ndimsa) : nothing
  nctype=xtypea[1]
  vndim=ndimsa[1]
  dimids=vndim>0 ? dimida[1:vndim] : []
  natts=natta[1]
  NC_VERBOSE ? println("Successfully read from file") : nothing
  name=_cchartostring(namea)
  isdimvar=false
  for n in nc.dim
    if (n[2].name==name)
      isdimvar=true
      break
    end
  end
  NC_VERBOSE ? println("name=",name," nctype=",nctype," dimids=",dimids," natts=",natts," vndim=",vndim) : nothing
  return (name,nctype,dimids,natts,vndim,isdimvar)
end

function _getvarindexbyname(nc::NcFile,varname::String)
  va=nothing
  for v in nc.vars
    va = v[2].name==varname ? v[2] : va
  end
  return va
end

function getvarbyid(nc::NcFile,varid::Integer)
  va=nothing
  for v in nc.vars
    va = v[2].varid==varid ? v[2] : va
  end
  return va
end

function getdimidbyname(nc::NcFile,dimname::String)
  da=nothing
  for d in nc.dim
    da = d.name==dimname ? d : da
  end
  return da.dimid
end

function getdimnamebyid(nc::NcFile,dimid::Integer)
  da=nothing
  for d in nc.dim
    da = d[2].dimid==dimid ? d[2] : da
  end
  return da.name
end

function _readdimdvar(ncid::Integer,dim::NcDim)
  start=0
  p=dim.dimlen
  #Determine size of Array
  retvalsa=Array(Float64,p)
  _nc_get_vara_double_c(ncid,varid,start,count,retvalsa)
  NC_VERBOSE ? println("Successfully read dimension from file ",dim.name) : nothing
  dim.vals=retvalsa
  end

function _nc_getatts_all(ncid::Integer,varid::Integer,natts::Integer)
  atts=Dict{Any,Any}()
  for attnum=0:natts-1
    gatt=_nc_inq_att(ncid,varid,attnum)
    v=gatt[2]
    if ((length(v)==1) & !(typeof(v)<:String))
      v=v[1]
    end
    atts[gatt[1]]=v
  end
  NC_VERBOSE ? println(atts) : nothing
  return atts
end

function _readdimvars(nc::NcFile)
  for d in nc.dim
    for v in nc.vars
      if (d[2].name==v[2].name)
        NC_VERBOSE ?println(d[2].name," ",v[2].name) : nothing
        d[2].vals=readvar(nc,v[2].varid,[1],[-1])
        d[2].atts=v[2].atts
      end
    end
  end
end

function parsedimargs(dim)
  idim=0
  dimlen=nothing
  dimvals=nothing
  dimatts=nothing
  name=nothing
  #Determine number of dimensions
  ndim=0
  for a in dim
    if (typeof(a)<:String)
      ndim=ndim+1
    end
    NC_VERBOSE ? println(a) : nothing
  end
  d=Array(NcDim,ndim)
  idim=1
  for a in dim
    NC_VERBOSE ? println(a,idim) : nothing
    if (typeof(a)<:String)
      #Assume a name is given
      #first create an NcDim object from the last dim
      if (name!=nothing)
        d[idim]=finalizedim(dimlen,dimvals,dimatts,name)
        idim=idim+1
        dimlen=nothing
        dimvals=nothing
        dimatts=nothing
        name=nothing
      end
      name=a
    elseif (typeof(a)<:Integer)
      #Assume a dimension length is given
      dimlen=a
    elseif (typeof(a)<:AbstractArray)
      #Assume dimension values are given
      if dimvals==nothing
        dimvals=float64(a)
        dimlen=length(dimvals)
      else
        error ("Dimension values of $name defined more than once")
      end
    elseif (typeof(a)<:Dict)
      #Assume attributes are given
      dimatts= dimatts==nothing ? a : error("Dimension attributes of $name defined more than once")
    end
  end
  d[idim]=finalizedim(dimlen,dimvals,dimatts,name)
  return(d)
end

function finalizedim(dimlen,dimvals,dimatts,name)
  if ((dimlen==nothing) & (dimvals==nothing))
    dimlen=1
  end
  if ((dimlen!=nothing) & (dimvals==nothing))
    dimvals=float64([1:dimlen])
  end
  if (dimatts==nothing)
    dimatts={"missval"=>-9999}
  end
  return(NcDim(name,dimvals,dimatts))
end


text:
# TODO: use threshold to optimize and avoid infinity

export Sigmoid

"""
    Sigmoid{T, OutputType}

Sigmoid activation.
"""
mutable struct Sigmoid{T, O} <: AbstractBlock
    output::O
    Sigmoid{T, O}() where {T, O} = new{T, O}()
end

Sigmoid(::Type{T}, ::Type{O}) where {T, O} =
    Sigmoid{T, O}()

Sigmoid(::Type{T}; nbatch=1) where T = Sigmoid(T, nbatch > 1 ? Matrix{T} : Vector{T})
Sigmoid(;nbatch=1) = Sigmoid(Float64; nbatch=nbatch)

function forward(op::Sigmoid{T, O}, input::O) where {T, O}
    op.output = 1 ./ (1 + exp(-input))
    return op.output
end

function backward(op::Sigmoid{T, O}, grad::O) where {T <: Real, O}
    grad .* (1.0 .- op.output) .* op.output
end


text:
function solve_diag(λ, Δ, g̃, seuil, α, ϵ; M = [0.0])
    # printstyled("on est ici ↓ \n", color = :red)
    T = eltype(α)
    # @show T
    if M == [0.0]
        # M = ones(Δ) ;
        # M = fill(T(1.0), size(Δ))
        M = ones(T, size(Δ))
    end
    λin = λ
    λtry = seuil
    dtry = -(Δ .+ λtry * M) .\ g̃
    normdtry = sqrt(dtry⋅dtry)
    seuiltry = normdtry / α
    d̃ = dtry

    #println(" Biss $biss_iter : λ = $λtry  seuil = $seuiltry ")

    # secant like heuristic to help Newton process by providing a good starting point
    global iter_bis = 0
    if λin < seuil
        while (λtry > seuiltry && (iter_bis < 10))
            λtry = λ + (λtry - λ) / max((seuil - seuiltry), T(2.0))
            dtry = -(Δ .+ λtry * M) .\ g̃
            normdtry = sqrt(dtry⋅dtry)
            seuiltry = normdtry / α
            global iter_bis += 1
            #println(" λ computation, iter_bis : $iter_bis  λtry  $λtry")
        end
        λ = λtry
        d̃ = dtry
        normd̃ = normdtry
        seuil = seuiltry
    end

    # Newton iterations
    global iter_nwt = 0
    while (abs(seuil - λ)/max(seuil,λ) > ϵ) && (iter_nwt < 40)
        dotd̃ = (Δ .+ λ * M) .\ (M .* d̃)
        Δλ = λ * (seuil - λ) / (seuil + λ * (λ * (d̃ ⋅ dotd̃) / (d̃ ⋅ d̃)))
        λ = max(λ + Δλ, λin + T(1.0e-10))
        λ = max(λ + Δλ, λin + (sqrt(eps(T) ./ 100.0)))
        d̃ = -(Δ .+ λ * M) .\ g̃
        normd̃ = sqrt(d̃⋅d̃)
        seuil = normd̃ / α

        @assert (λ >= 0.0)
        global iter_nwt += 1
    end
    #println(" λ computation, iter_bis : $iter_bis  iter_nwt : $iter_nwt  λ  $λ")
    #try assert((g̃ + 0.5*Δ.*d̃)⋅d̃ <= 0.0)  catch  @bp  end
    #println("*******SolveDiag:  g̃⋅d̃ = $(g̃⋅d̃), 0.5 d̃'Δd̃ = $(0.5*(Δ .* d̃)⋅d̃)")

    return d̃,λ
end


text:
real_mod(x) = real(x)
function call_to_put(C, underlying, zero_rate, opt::EuropeanOption)
    S = underlying.S0
    T = opt.T
    d = FinancialMonteCarlo.integral(underlying.d, T) / T
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    K = opt.K
    return opt.isCall ? C : (C - S * exp(-d * T) + K * exp(-r * T))
end

function call_to_put(C, underlying, zero_rate, opt::BinaryEuropeanOption)
    T = opt.T
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    return opt.isCall ? C : (exp(-r * T) - C)
end

function midpoint_definite_integral(f, xmin, xmax, N)
    x = range(xmin, length = N, stop = xmax)
    dx = (xmax - xmin) / (N - 1)
    sum_ = zero(typeof(f(xmin))) * dx * 0
    weights_ = @. 3 + (-1)^((0:(N-1)) + 1)
    @views weights_[1] = 1
    @views weights_[N] = 1
    for (x_, w) in zip(x, weights_)
        x_ != 0 && !isnan(f(x_)) && (sum_ += f(x_) * w * dx / 3)
    end
    return sum_
end


text:
@testset "TFIDF" begin

    doc1 = "a a a sample text text"
    doc2 = "another example example text text"
    doc3 = ""
    doc4 = "another another text text text text"

    crps = Corpus(StringDocument.([doc1, doc2, doc3, doc4]))

    update_lexicon!(crps)
    m = DocumentTermMatrix(crps)

    # Terms are in alphabetical ordering
    correctweights =[0.5  0.0  0.0  1/6  1/3
                     0.0  0.2  0.4  0.0  0.4
                     0.0  0.0  0.0  0.0  0.0
                     0.0  1/3  0.0  0.0  2/3]

    myweights = tf(m)
    @test myweights == correctweights

    myweights = tf(dtm(m))
    @test myweights ≈ sparse(correctweights)
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf(dtm(m, :dense))
    @test isnan(sum(myweights)) == 0
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(m));
    tf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = float(dtm(m, :dense));
    tf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    # Terms are in alphabetical ordering
    correctweights = [0.6931471805599453 0.0 0.0 0.23104906018664842 0.09589402415059362
	              0.0 0.13862943611198905 0.5545177444479562 0.0 0.11507282898071235
	              0.0 0.0 0.0 0.0 0.0
	              0.0 0.23104906018664842  0.0 0.0 0.19178804830118723]

    myweights = tf_idf(m)
    @test myweights ≈ correctweights

    myweights = tf_idf(dtm(m))
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf_idf(dtm(m, :dense))
    @test isnan(sum(myweights)) == 0
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(m));
    tf_idf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = float(dtm(m, :dense));
    tf_idf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    @testset "Okapi BM25" begin
        doc1 = "a a a sample text text"
        doc2 = "another example example text text"
        doc3 = ""
        doc4 = "another another text text text text"

        crps = Corpus(StringDocument.([doc1, doc2, doc3, doc4]))

        update_lexicon!(crps)
        m = DocumentTermMatrix(crps)

        max_tol = 1e-5
        approx_eq(m1::AbstractMatrix{T}, m2::AbstractMatrix{T}; tol=1e-6
                ) where T = begin
            t = all(abs.(m1-m2) .<= tol)
            return t
        end

        correctweights = [1.29959 0.0 0.0 1.89031 0.405067;
                          0.0 0.882404 1.54025 0.0 0.405067;
                          0.0 0.0 0.0 0.0 0.0;
                          0.0 1.40179 0.0 0.0 0.676646]

        myweights = bm_25(m)

        @test approx_eq(myweights, correctweights, tol=max_tol)

        myweights = bm_25(dtm(m))
        @test approx_eq(myweights, correctweights, tol=max_tol)
        @test typeof(myweights) <: SparseMatrixCSC

        myweights = bm_25(Matrix(dtm(m)))
        @test isnan(sum(myweights)) == 0
        @test approx_eq(myweights, correctweights, tol=max_tol)
        @test typeof(myweights) <: Matrix

        myweights = float(dtm(crps))
        myweights = bm_25(myweights)
        @test approx_eq(Float64.(myweights), correctweights, tol=max_tol)
        @test typeof(myweights) <: SparseMatrixCSC
        @test eltype(myweights) == typeof(1.0)

        myweights = float(Matrix(dtm(crps)));
        myweights = bm_25(myweights)
        @test approx_eq(Float64.(myweights), correctweights, tol=max_tol)
        @test typeof(myweights) <: Matrix
        @test eltype(myweights) == typeof(1.0)

        @test_throws MethodError bm_25!(DocumentTermMatrix(crps))
    end

	@testset "cosine similarity `cos_similarity()`" begin
		crps = Corpus( StringDocument.([
		    "to be or not to be",
		    "to sing or not to sing",
		    "to talk or to silence"]) )
		update_lexicon!(crps)
		d = dtm(crps)
		tfm = tf_idf(d)
		cs = cos_similarity(tfm)
		@test cs==[
			1.0                  0.03293177886245518  0.0;
			0.03293177886245518  1.0                  0.0;
			0.0                  0.0                  1.0
		]
    end
end


text:
using WriteVTK
using Printf

let

    dir = string("../../erickson/output_files/", ARGS[1])
    u_file = open(string(dir, "/us.dat"))

    @printf "Opened u files...\n"

    us = eachline(u_file)

    @printf "Read u and v data\n..."

    x = map(x->parse(Float64, x), split.(Iterators.take(us, 1))[1])
    y = map(x->parse(Float64, x), split.(Iterators.take(us, 2))[1])

    count = 1
    @printf "Converting to vtk...\n"
    for line in us
        @printf "\r%d" count
        vtk = vtk_grid(string(dir, "/vtk_filesU/cycles_", count, ".vtr"), x, y)
        u = map(x-> parse(Float64, x), split.(line))
        vtk["u"] = u
        saved_files = vtk_save(vtk)
        count += 1
    end

    close(u_file)

    @printf ("Done.\n")

end