text:
function maybe_track(arg, has_argument_grad::Bool, tape)
    has_argument_grad ? track(arg, tape) : arg
end

@noinline function ReverseDiff.special_reverse_exec!(
        instruction::ReverseDiff.SpecialInstruction{D}) where {D <: Distribution}
    dist::D = instruction.func
    args_maybe_tracked = instruction.input
    score_tracked = instruction.output
    arg_grads = logpdf_grad(dist, map(value, args_maybe_tracked)...)
    value_tracked = args_maybe_tracked[1]
    value_grad = arg_grads[1]
    if istracked(value_tracked)
        if has_output_grad(dist)
            increment_deriv!(value_tracked, value_grad * deriv(score_tracked))
        else
            error("Gradient required but not available for return value of distribution $dist")
        end
    end
    for (i, (arg_maybe_tracked, grad, has_grad)) in enumerate(
            zip(args_maybe_tracked[2:end], arg_grads[2:end], has_argument_grads(dist)))
        if istracked(arg_maybe_tracked)
            if has_grad
                increment_deriv!(arg_maybe_tracked, grad * deriv(score_tracked))
            else
                error("Gradient required but not available for argument $i of $dist")
            end
        end
    end
    nothing
end

###################
# accumulate_param_gradients! #
###################

mutable struct GFBackpropParamsState
    trace::DynamicDSLTrace
    score::TrackedReal
    tape::InstructionTape
    visitor::AddressVisitor
    scaler::Float64

    # only those tracked parameters that are in scope (not including splice calls)
    tracked_params::Dict{Symbol,Any}

    # tracked parameters for all (nested) splice calls
    splice_tracked_params::Dict{Tuple{GenerativeFunction,Symbol},Any}
end

function track_params(tape, params)
    tracked_params = Dict{Symbol,Any}()
    for (name, value) in params
        tracked_params[name] = track(value, tape)
    end
    tracked_params
end

function GFBackpropParamsState(trace::DynamicDSLTrace, tape, params, scaler)
    tracked_params = track_params(tape, params)
    splice_tracked_params = Dict{Tuple{GenerativeFunction,Symbol},Any}()
    score = track(0., tape)
    GFBackpropParamsState(trace, score, tape, AddressVisitor(), scaler,
        tracked_params, splice_tracked_params)
end

function read_param(state::GFBackpropParamsState, name::Symbol)
    value = state.tracked_params[name]
    value
end

function traceat(state::GFBackpropParamsState, dist::Distribution{T},
              args_maybe_tracked, key) where {T}
    local retval::T
    visit!(state.visitor, key)
    retval = get_choice(state.trace, key).retval
    args = map(value, args_maybe_tracked)
    score_tracked = track(logpdf(dist, retval, args...), state.tape)
    record!(state.tape, ReverseDiff.SpecialInstruction, dist,
        (retval, args_maybe_tracked...,), score_tracked)
    state.score += score_tracked
    retval
end

struct BackpropParamsRecord
    gen_fn::GenerativeFunction
    subtrace::Any
    scaler::Float64
end

function traceat(state::GFBackpropParamsState, gen_fn::GenerativeFunction{T},
              args_maybe_tracked, key) where {T}
    local retval::T
    visit!(state.visitor, key)
    subtrace = get_call(state.trace, key).subtrace
    retval = get_retval(subtrace)
    if accepts_output_grad(gen_fn)
        retval_maybe_tracked = track(retval, state.tape)
        @assert istracked(retval_maybe_tracked)
    else
        retval_maybe_tracked = retval
        @assert !istracked(retval_maybe_tracked)
    end
    record!(state.tape, ReverseDiff.SpecialInstruction,
        BackpropParamsRecord(gen_fn, subtrace, state.scaler),
        (args_maybe_tracked...,), retval_maybe_tracked)
    retval_maybe_tracked
end

function splice(state::GFBackpropParamsState, gen_fn::DynamicDSLFunction,
                args_maybe_tracked::Tuple)

    # save previous tracked parameter scope
    prev_tracked_params = state.tracked_params

    # construct new tracked parameter scope
    state.tracked_params = Dict{Symbol,Any}()
    for name in keys(gen_fn.params)
        if haskey(state.splice_tracked_params, (gen_fn, name))
            # parameter was already tracked in another splice call
            state.tracked_params[name] = state.splice_tracked_params[(gen_fn, name)]
        else
            # parameter was not already tracked
            tracked = track(get_param(gen_fn, name), state.tape)
            state.tracked_params[name] = tracked
            state.splice_tracked_params[(gen_fn, name)] = tracked
        end
    end

    retval_maybe_tracked = exec(gen_fn, state, args_maybe_tracked)

    # restore previous tracked parameter scope
    state.tracked_params = prev_tracked_params

    retval_maybe_tracked
end

@noinline function ReverseDiff.special_reverse_exec!(
        instruction::ReverseDiff.SpecialInstruction{BackpropParamsRecord})
    record = instruction.func
    gen_fn = record.gen_fn
    args_maybe_tracked = instruction.input
    retval_maybe_tracked = instruction.output
    if accepts_output_grad(gen_fn)
        @assert istracked(retval_maybe_tracked)
        retval_grad = deriv(retval_maybe_tracked)
    else
        @assert !istracked(retval_maybe_tracked)
        retval_grad = nothing
    end
    arg_grads = accumulate_param_gradients!(record.subtrace, retval_grad, record.scaler)

    # if has_argument_grads(gen_fn) is true for a given argument, then that
    # argument may or may not be tracked. if has_argument_grads(gen_fn) is
    # false for a given argument, then that argument must not be tracked.
    # otherwise it is an error.
    # note: code duplication with accumulate_param_gradients!
    for (i, (arg, grad, has_grad)) in enumerate(
            zip(args_maybe_tracked, arg_grads, has_argument_grads(gen_fn)))
        if has_grad && istracked(arg)
            increment_deriv!(arg, grad)
        elseif !has_grad && istracked(arg)
            error("Gradient required but not available for argument $i of $gen_fn")
        end
    end
    nothing
end

function accumulate_param_gradients!(trace::DynamicDSLTrace, retval_grad, scaler=1.)
    gen_fn = trace.gen_fn
    tape = new_tape()
    state = GFBackpropParamsState(trace, tape, gen_fn.params, scaler)
    args = get_args(trace)
    args_maybe_tracked = (map(maybe_track,
        args, gen_fn.has_argument_grads, fill(tape, length(args)))...,)
    retval_maybe_tracked = exec(gen_fn, state, args_maybe_tracked)

    # note: code duplication with choice_gradients
    if accepts_output_grad(gen_fn)
        if retval_grad == nothing
            error("Return value gradient required but not provided")
        end
        if istracked(retval_maybe_tracked)
            deriv!(retval_maybe_tracked, retval_grad)
        end
        # note: if accepts_output_grad(gen_fn) and
        # !istracked(retval_maybe_tracked), this means the return value did not
        # depend on the gradient source elements for this trace. that is okay.
    else
        if retval_grad != nothing
            error("Return value gradient not supported, but got $retval_grad != nothing")
        end
    end
    seed!(state.score)
    reverse_pass!(tape)

    # increment the gradient accumulators for trainable parameters in scope
    for (name, tracked) in state.tracked_params
        gen_fn.params_grad[name] += deriv(tracked) * state.scaler
    end

    # increment the gradient accumulators for trainable parameters in splice calls
    for ((spliced_gen_fn, name), tracked) in state.splice_tracked_params
        spliced_gen_fn.params_grad[name] += deriv(tracked) * state.scaler
    end

    # return gradients with respect to arguments with gradients, or nothing
    input_grads::Tuple = (map((arg, has_grad) -> has_grad ? deriv(arg) : nothing,
                             args_maybe_tracked, gen_fn.has_argument_grads)...,)
    input_grads
end

##################
# choice_gradients #
##################

mutable struct GFBackpropTraceState
    trace::DynamicDSLTrace
    score::TrackedReal
    tape::InstructionTape
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
    selection::Selection
    tracked_choices::Trie{Any,TrackedReal}
    value_choices::DynamicChoiceMap
    gradient_choices::DynamicChoiceMap
end

function GFBackpropTraceState(trace, selection, params, tape)
    score = track(0., tape)
    visitor = AddressVisitor()
    tracked_choices = Trie{Any,TrackedReal}()
    value_choices = choicemap()
    gradient_choices = choicemap()
    GFBackpropTraceState(trace, score, tape, visitor, params,
        selection, tracked_choices, value_choices, gradient_choices)
end

function fill_gradient_map!(gradient_choices::DynamicChoiceMap,
                             tracked_trie::Trie{Any,TrackedReal})
    for (key, tracked) in get_leaf_nodes(tracked_trie)
        set_value!(gradient_choices, key, deriv(tracked))
    end
    # NOTE: there should be no address collision between these primitive
    # choices and the gen_fn invocations, as enforced by the visitor
    for (key, subtrie) in get_internal_nodes(tracked_trie)
        @assert !has_value(gradient_choices, key) && isempty(get_submap(gradient_choices, key))
        gradient_submap = choicemap()
        fill_gradient_map!(gradient_submap, subtrie)
        set_submap!(gradient_choices, key, gradient_submap)
    end
end

function fill_value_map!(value_choices::DynamicChoiceMap,
                          tracked_trie::Trie{Any,TrackedReal})
    for (key, tracked) in get_leaf_nodes(tracked_trie)
        set_value!(value_choices, key, value(tracked))
    end
    # NOTE: there should be no address collision between these primitive
    # choices and the gen_fn invocations, as enforced by the visitor
    for (key, subtrie) in get_internal_nodes(tracked_trie)
        @assert !has_value(value_choices, key) && !isempty(get_submap(value_choices, key))
        value_choices_submap = choicemap()
        fill_value_map!(value_choices_submap, subtrie)
        set_submap!(value_choices, key, value_choices_submap)
    end
end

function traceat(state::GFBackpropTraceState, dist::Distribution{T},
              args_maybe_tracked, key) where {T}
    local retval::T
    visit!(state.visitor, key)
    retval = get_choice(state.trace, key).retval
    args = map(value, args_maybe_tracked)
    score_tracked = track(logpdf(dist, retval, args...), state.tape)
    if key in state.selection
        tracked_retval = track(retval, state.tape)
        set_leaf_node!(state.tracked_choices, key, tracked_retval)
        record!(state.tape, ReverseDiff.SpecialInstruction, dist,
            (tracked_retval, args_maybe_tracked...,), score_tracked)
        state.score += score_tracked
        return tracked_retval
    else
        record!(state.tape, ReverseDiff.SpecialInstruction, dist,
            (retval, args_maybe_tracked...,), score_tracked)
        state.score += score_tracked
        return retval
    end
end

struct BackpropTraceRecord
    gen_fn::GenerativeFunction
    subtrace::Any
    selection::Selection
    value_choices::DynamicChoiceMap
    gradient_choices::DynamicChoiceMap
    key::Any
end

function traceat(state::GFBackpropTraceState, gen_fn::GenerativeFunction{T,U},
              args_maybe_tracked, key) where {T,U}
    local retval::T
    local subtrace::U
    visit!(state.visitor, key)
    subtrace = get_call(state.trace, key).subtrace
    get_gen_fn(subtrace) === gen_fn || gen_fn_changed_error(key)
    retval = get_retval(subtrace)
    if accepts_output_grad(gen_fn)
        retval_maybe_tracked = track(retval, state.tape)
        @assert istracked(retval_maybe_tracked)
    else
        retval_maybe_tracked = retval
        @assert !istracked(retval_maybe_tracked)
    end
    selection = state.selection[key]
    record = BackpropTraceRecord(gen_fn, subtrace, selection, state.value_choices,
        state.gradient_choices, key)
    record!(state.tape, ReverseDiff.SpecialInstruction, record, (args_maybe_tracked...,), retval_maybe_tracked)
    retval_maybe_tracked
end

function splice(state::GFBackpropTraceState, gen_fn::DynamicDSLFunction,
                args_maybe_tracked::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args_maybe_tracked)
    state.params = prev_params
    retval
end

@noinline function ReverseDiff.special_reverse_exec!(
        instruction::ReverseDiff.SpecialInstruction{BackpropTraceRecord})
    record = instruction.func
    gen_fn = record.gen_fn
    args_maybe_tracked = instruction.input
    retval_maybe_tracked = instruction.output
    if accepts_output_grad(gen_fn)
        @assert istracked(retval_maybe_tracked)
        retval_grad = deriv(retval_maybe_tracked)
    else
        @assert !istracked(retval_maybe_tracked)
        retval_grad = nothing
    end
    (arg_grads, value_choices, gradient_choices) = choice_gradients(
        record.subtrace, record.selection, retval_grad)
    @assert isempty(get_submap(record.gradient_choices, record.key))
    @assert !has_value(record.gradient_choices, record.key)
    set_submap!(record.gradient_choices, record.key, gradient_choices)
    set_submap!(record.value_choices, record.key, value_choices)

    # if has_argument_grads(gen_fn) is true for a given argument, then that
    # argument may or may not be tracked. if has_argument_grads(gen_fn) is
    # false for a given argument, then that argument must not be tracked.
    # otherwise it is an error.
    # note: code duplication with accumulate_param_gradients!
    for (i, (arg, grad, has_grad)) in enumerate(
            zip(args_maybe_tracked, arg_grads, has_argument_grads(gen_fn)))
        if has_grad && istracked(arg)
            increment_deriv!(arg, grad)
        elseif !has_grad && istracked(arg)
            error("Gradient required but not available for argument $i of $gen_fn")
        end
    end
    nothing
end

function choice_gradients(trace::DynamicDSLTrace, selection::Selection, retval_grad)
    gen_fn = trace.gen_fn
    tape = new_tape()
    state = GFBackpropTraceState(trace, selection, gen_fn.params, tape)
    args = get_args(trace)
    args_maybe_tracked = (map(maybe_track,
        args, gen_fn.has_argument_grads, fill(tape, length(args)))...,)
    retval_maybe_tracked = exec(gen_fn, state, args_maybe_tracked)

    # note: code duplication with accumulate_param_gradients!
    if accepts_output_grad(gen_fn)
        if retval_grad == nothing
            error("Return value gradient required but not provided")
        end
        if istracked(retval_maybe_tracked)
            deriv!(retval_maybe_tracked, retval_grad)
        end
        # note: if accepts_output_grad(gen_fn) and
        # !istracked(retval_maybe_tracked), this means the return value did not
        # depend on the gradient source elements for this trace. that is okay.
    else
        if retval_grad != nothing
            error("Return value gradient not supported, but got $retval_grad != nothing")
        end
    end

    seed!(state.score)
    reverse_pass!(tape)

    # fill trace gradient with gradients with respect to primitive random choices
    fill_gradient_map!(state.gradient_choices, state.tracked_choices)
    fill_value_map!(state.value_choices, state.tracked_choices)

    # return gradients with respect to arguments with gradients, or nothing
    # NOTE: if a value isn't tracked the gradient is nothing
    input_grads::Tuple = (map((arg, has_grad) -> has_grad ? deriv(arg) : nothing,
                             args_maybe_tracked, gen_fn.has_argument_grads)...,)

    (input_grads, state.value_choices, state.gradient_choices)
end

max_stars_repo_name: dpmerrell/Gen

max_issues_repo_name: rameshputalapattu/Gen

max_forks_repo_name: rameshputalapattu/Gen


text:
#=
# 3. Applying pulse forcing to a flow
Here, we will show how to apply a transient force to a flow.
=#

using ViscousFlow
#-
using Plots

#=
For this case, we first seek to introduce a pulse to an otherwise quiescent fluid,
with no boundaries. The pulse shape will be a smooth Gaussian-distributed force, directed upward.

We will start with the usual steps: specify the problem parameters and the discretization.
=#
my_params = Dict()

my_params["Re"] = 200
xlim = (-2.0,2.0)
ylim = (-2.0,4.0)
my_params["grid Re"] = 4.0
g = setup_grid(xlim,ylim,my_params)

#=
### Construct the forcing
There are a variety of forcing types available to us. Here, we will
use *area-type forcing*, which is distributed over a region of space.
We will create a short-lived force applied in the vertical direction in this region. For the
spatial distribution of this forcing, we make use of a Gaussian field. We will center it at (0,0), and
give it a strength of 10, using the `SpatialGaussian` constructor. The x component
of the force we set to zero, using the `EmptySpatialField` constructor.
=#
σx = 0.5
σy = 0.1
x0 = y0 = 0.0
amp = 10
force_dist = [EmptySpatialField(),SpatialGaussian(σx,σy,x0,y0,amp)];

#=

When we introduce forcing in any problem, it is essential to
supply a *model function* that supplies the instantaneous strength of the
forcing. This gives us the opportunity to specify the transient
aspect of the forcing. In this example, we modulate the force in time with a Gaussian,
with a small half-width in time equal to $\sigma_t = 0.1$, effectively creating a pulse.
We center it at time $t_0 = 0.1$ (using the `>>` shift operator). This is carried out in the function below,
in which we create the modulating function `modfcn`.
Note that `fr.generated_field()` provides the spatial distribution that
we specified with our spatial field `force_dist`. This is multiplied
by `modfcn`, our Gaussian pulse, evaluated at the current time.
=#
function forcing_model!(σ,T,t,fr::AreaRegionCache,phys_params)
    σt = 0.1
    t0 = 0.1
    modfcn = Gaussian(σt,sqrt(π*σt^2)) >> t0
    σ .= modfcn(t)*fr.generated_field()
end

#=
We pack the model function and the spatial field together
into an `AreaForcingModel`, which we then put into a forcing dictionary.
It is important to use the key "forcing models" for this.
=#
afm = AreaForcingModel(forcing_model!,spatialfield=force_dist)
forcing_dict = Dict("forcing models" => afm);

#=
### Construct the system structure
We supply these forcing characteristics with the keyword argument `forcing` as we
set up the usual system structure:
=#
sys = viscousflow_system(g,phys_params=my_params,forcing=forcing_dict);

#=
The remaining steps go just as they did for the previous examples. We will simulate
the pulse for 4 time units:
=#
u0 = init_sol(sys)
tspan = (0.0,10.0) # longer than we need, but just being safe
integrator = init(u0,tspan,sys)
#-
step!(integrator,4.0)

#=
### Examine
Let's look at a sequence of snapshots of the forced flow. As we see, this
short-lived force creates a dipole-like pair of counter-rotating vortices that propagates upward under its
own influence.
=#
sol = integrator.sol
plt = plot(layout = (1,4), size = (800, 300), legend=:false)
tsnap = 0.1:1.0:3.1
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,title="t = $(round(t,digits=2))")
end
plt
#=
And here is another look a the flow, this time with streamlines:
=#
plot(streamfunction(integrator),sys,title="Streamfunction at t = $(round(integrator.t,digits=2))")

#=
We can also supply more than one pulse. If we don't change the spatial
distribution, then we just redefine the model function and run it again:
Let's give 3 pulses, each separated by 1 time unit.
=#
function forcing_model!(σ,T,t,fr::AreaRegionCache,phys_params)
    σt = 0.1
    t0 = [0.1, 1.1, 2.1]
    modfcn = (Gaussian(σt,sqrt(π*σt^2)) >> t0[1]) +
             (Gaussian(σt,sqrt(π*σt^2)) >> t0[2]) +
             (Gaussian(σt,sqrt(π*σt^2)) >> t0[3])
    σ .= modfcn(t)*fr.generated_field()
end
#-
u0 = init_sol(sys)
integrator = init(u0,tspan,sys)
#-
step!(integrator,4.0)
#=
### Examine
In this case, the pulses coalesce with one another
=#
sol = integrator.sol
plt = plot(layout = (1,4), size = (800, 300), legend=:false)
tsnap = 0.1:1.0:3.1
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,title="t = $(round(t,digits=2))")
end
plt

ext: jl


text:
push!(LOAD_PATH, "../")
push!(LOAD_PATH, "./")

include("hopperfunctions.jl")

####################################### Experiment Parameters / HyperParameters

model_file = "$(modeldir)/challenge_hopper.xml"

myskip = 5
dtype = Float64

my_mjsys       = mjw.load_model(model_file, myskip, "normal")

function getNPG()
    T = 800
    K = 32 #64
    niter = 80
    gamma = 0.995
    gae = 0.98

    fullFIM = true #false

    norm_step_size = 0.05

    cg_iter= 18
    cg_reg = 1e-6
    cg_tol = 1e-10

    my_policy      = Policy.GLP{dtype}(my_mjsys.ns, my_mjsys.nu) # inputs: n, m
    #my_policy      = Policy.NN{dtype}(my_mjsys.ns, my_mjsys.nu, 32) # inputs: n, m

    pgmodel        = PolicyGradModel(copy(my_policy.theta),
                                     my_mjsys,
                                     my_policy,
                                     T, K, T, 10, niter)

    # NPG baseline aka value function approximation
    baseline    = Baseline.Quadratic{Float64}(1e-5, my_mjsys.ns, T, K)
    #baseline    = Baseline.NN{Float64}(my_mjsys.ns, T, K, 128, 0.001, Baseline.getNNfeatures)

    pg_specs   = NPGStrategy{dtype}(pgmodel,
                                    baseline,
                                    fullFIM,
                                    norm_step_size,
                                    gamma, gae,
                                    cg_iter, cg_reg, cg_tol, myfuncs)

    return pgmodel, pg_specs
end

function getOpt()
    opt_K = 120 #240 #256 #512
    H = div(128, myskip)
    T = div(1000, myskip)
    mpcT = 100

    opt   = TrajOptModel(my_mjsys, H, T, opt_K;
                         theta = zeros(my_mjsys.nu, H)#,
                         #s0 = [my_mjsys.d.qpos; my_mjsys.d.qvel]
                        )

    opt_specs = MPPIStrategy(opt,
                             Matrix(0.2*I, my_mjsys.nu, my_mjsys.nu), # sigma
                             0.8,
                             gamma,
                             myfuncs)

    return opt, opt_specs
end


text:
# TODO: Rename this file

using Plots
using DSP
using FFTW
include("NoiseWaveform.jl")

"""
write_binary_file(data,filename,dtype)

# Arguments
- `data`: Data to be saved
- `filename::String`: Name of the file to be saved (path relative to current directory)
- `dtype::DataType`: Data type of the data to be saved
                     Default: `Matrix{ComplexF32}`
"""
function write_binary_file(data, filename::String, dtype::DataType = Matrix{ComplexF32})
    f = open(filename, "w")
    data = convert(dtype, data)
    write(f, data)
    close(f)
end

# PCFM Waveform Parameters
# ---------------
nWaveforms = 1
# Oversampled code length
m = 300
# Oversampling factor
k = 3
# Log base for log-FTE computation
a = 10

# Window Parameters
# -----------------
# Window length
l = 2 * m - 1
# Full-width 3dB (normalized) bandwidth
bandwidth = 0.1
sigma = bandwidth / (2 * sqrt(2 * log(2)))
# Window standard deviation
# sigma = 0.1
u = DSP.Windows.gaussian((l, 1), sigma)
# In practice, we probably won't be able to optimize the PSD below -50 dB, so
# clip all window values below it
u[findall(<(-50), 10 * log10.(u))] .= 10^-5

# Simulation Parameters
# ---------------------
maxIter = 1000
ϵ = 1e-5
filename = "figures/gaussian.gif"
(x, s) = NoiseWaveform.optimize(u, nWaveforms, k, maxIter = maxIter, savePlots = true, filename = filename)

write_binary_file(s, "data/gaussian.bin")


text:
using ProfileView
using BenchmarkTools
using PolynomialBasis
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function plane_distance_function(coords, normal, x0)
    return (coords .- x0)' * normal
end

x0 = [0.0,0.0]
normal = [1.0,0.0]
polyorder = 3
numqp = 5
levelset = InterpolatingPolynomial(1,2,polyorder)
levelsetcoeffs = plane_distance_function(levelset.basis.points,normal,x0)
update!(levelset,levelsetcoeffs)

# @profiler run_quadrature_iterations(levelset,numqp,100)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
using Libdl

MPI_PATH = get(ENV, "JULIA_MPI_PATH", nothing)
MPI_LIBRARY_PATH = get(ENV, "JULIA_MPI_LIBRARY_PATH") do
    MPI_PATH !== nothing ? joinpath(MPI_PATH,"lib") : nothing
end
MPI_INCLUDE_PATH = get(ENV, "JULIA_MPI_INCLUDE_PATH") do
    MPI_PATH !== nothing ? joinpath(MPI_PATH,"include") : nothing
end
mpicc = get(ENV, "JULIA_MPICC") do
    if MPI_PATH !== nothing
        joinpath(MPI_PATH,"bin","mpicc")
    else
        "mpicc"
    end
end

const mpiexec_path = get(ENV, "JULIA_MPIEXEC") do
    if MPI_PATH !== nothing && Sys.isexecutable(joinpath(MPI_PATH,"bin","mpiexec"))
        joinpath(MPI_PATH,"bin","mpiexec")
    else
        Sys.which("mpiexec")
    end
end

const mpiexec_args = split(get(ENV, "JULIA_MPIEXEC_ARGS", ""))

const libmpi = get(ENV, "JULIA_MPI_LIBRARY") do
    libmpi = find_library(["libmpi", "msmpi", "libmpich"],
                          MPI_LIBRARY_PATH !== nothing ? [MPI_LIBRARY_PATH] : [])
    if libmpi == ""
        error("No MPI library found.\nEnsure an MPI implementation is loaded, or set the `JULIA_MPI_PATH` variable.")
    end
    # expand paths
    dlpath(libmpi)
end

if haskey(ENV, "JULIA_MPI_CFLAGS")
    CFLAGS = split(ENV["JULIA_MPI_CFLAGS"])
else
    lname, = split(basename(libmpi),'.')
    if startswith(lname, "lib")
        lname = lname[4:end]
    end
    CFLAGS = ["-l$lname"]
    if MPI_LIBRARY_PATH !== nothing
        push!(CFLAGS, "-L$(MPI_LIBRARY_PATH)")
    end
    if MPI_INCLUDE_PATH !== nothing
        push!(CFLAGS, "-I$(MPI_INCLUDE_PATH)")
    end
end

libsize = filesize(libmpi)

@info "Using MPI library $libmpi"

function Get_version()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    if Sys.iswindows()
        ccall((:MPI_Get_version, libmpi), stdcall, Cint,
              (Ptr{Cint}, Ptr{Cint}), major, minor)
    else
        ccall((:MPI_Get_version, libmpi), Cint,
              (Ptr{Cint}, Ptr{Cint}), major, minor)
    end
    VersionNumber(major[], minor[])
end

MPI_VERSION = Get_version()

function Get_library_version()
    # There is no way to query at runtime what the length of the buffer should be.
    # https://github.com/mpi-forum/mpi-issues/issues/159
    # 8192 is the maximum value of MPI_MAX_LIBRARY_VERSION_STRING across known
    # implementations.
    buf = Array{UInt8}(undef, 8192)
    buflen = Ref{Cint}()
    if Sys.iswindows()
        ccall((:MPI_Get_library_version, libmpi), stdcall, Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    else
        ccall((:MPI_Get_library_version, libmpi), Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    end
    resize!(buf, buflen[])
    return String(buf)
end

@info "MPI library version:\n$(Get_library_version())"

open("deps.jl","w") do f
    println(f, "# This file is automatically generated")
    println(f, "# Do not edit")
    println(f)
    println(f, :(const libmpi = $libmpi))
    println(f, :(const libmpi_size = $libsize))
    println(f, :(const MPI_VERSION = $MPI_VERSION))
    println(f, :(const mpiexec_path = $mpiexec_path))

    if Sys.iswindows()
        println(f, :(include("consts_msmpi.jl")))
    else
        include("gen_consts.jl")

        run(`$mpicc gen_consts.c -o gen_consts $CFLAGS`)
        run(`$mpiexec_path $mpiexec_args -n 1 ./gen_consts`)

        println(f, :(include("consts.jl")))
    end
end


text:
using DebuggerFramework
using ASTInterpreter2

# Simple evaling of function argument
function evalfoo1(x,y)
    x+y
end
frame = ASTInterpreter2.enter_call_expr(:($(evalfoo1)(1,2)))
res = DebuggerFramework.eval_code(nothing, frame, "x")
@assert res == 1

res = DebuggerFramework.eval_code(nothing, frame, "y")
@assert res == 2

# Evaling with sparams
function evalsparams(x::T) where T
    x
end
frame = ASTInterpreter2.enter_call_expr(:($(evalsparams)(1)))
res = DebuggerFramework.eval_code(nothing, frame, "x")
@assert res == 1

res = DebuggerFramework.eval_code(nothing, frame, "T")
@assert res == Int


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#Commonly used 3rd party packages
#Install all in one swoop.

#Update all packages -- NOTE this is dangerous, as Julia is still a young
#language. Updates frequently break things:
# Pkg.update()

Pkg.add("DifferentialEquations")
# Pkg.add("Gurobi") #Optional, for direct method, requires license
Pkg.add("Plots")
Pkg.add("PlotlyJS")
Pkg.add("Ipopt") #Open-source alternative to Gurobi
# Pkg.add("JPLEphemeris") #for interplanetary problems, not CRTBP
# Pkg.add("MATLAB") #Optional, used for plotting, requires license
Pkg.add("JuMP")

#Trigger pre-compilation:
using DifferentialEquations
# using Gurobi
using Plots
using PlotlyJS
using Ipopt
# using JPLEphemeris
# using MATLAB
using JuMP

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
export SCFconfig, runHF, runHFcore

using LinearAlgebra: dot, Hermitian, \, det, I
using PiecewiseQuadratics: indicator
using SeparableOptimization
using Combinatorics: powerset

const TelLB = Float64 # Union{} to loose constraint
const TelUB = Float64 # Any to loose constraint

getXcore1(S::Matrix{T}) where {TelLB<:T<:TelUB} = S^(-0.5) |> Array

const getXmethods = Dict{Int, Function}(1=>getXcore1)

getX(S::Matrix{T}; method::Int=1) where {TelLB<:T<:TelUB} = getXmethods[method](S)

function getC(X::Matrix{T1}, F::Matrix{T2};
              outputCx::Bool=false, outputEmo::Bool=false, stabilizeSign::Bool=true) where
             {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB}
    ϵ, Cₓ = eigen(X'*F*X |> Hermitian, sortby=x->x)
    outC = outputCx ? Cₓ : X*Cₓ
    # Stabilize the sign factor of each column.
    stabilizeSign && for j = 1:size(outC, 2)
       outC[:, j] *= (outC[1,j] < 0 ? -1 : 1)
    end
    outputEmo ? (outC, ϵ) : outC
end

function breakSymOfC(C::Matrix{T}) where {TelLB<:T<:TelUB}
    C2 = copy(C)
    l = min(size(C2)[1], 2)
    C2[1:l, 1:l] .= 0 # Breaking spin symmetry.
    # C2[l, :] .= 0 # Another way.
    (copy(C), C2)
end

function getCfromGWH(S::Matrix{T1}, Hcore::Matrix{T2}, K::Float64=1.75; X=getX(S),
                     forUHF::Bool=false) where
                    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB}
    l = size(Hcore)[1]
    H = zero(Hcore)
    for i in 1:l, j in 1:l
        H[i,j] = K * S[i,j] * (Hcore[i,i] + Hcore[j,j]) * 0.5
    end
    C = getC(X, H)
    forUHF ? breakSymOfC(C) : C
end

function getCfromHcore(X::Matrix{T1}, Hcore::Matrix{T2}; forUHF::Bool=false) where
                      {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB}
    C = getC(X, Hcore)
    forUHF ? breakSymOfC(C) : C
end

function getCfromSAD(S::Matrix{T1}, Hcore::Matrix{T2}, HeeI::Array{T3, 4},
                     bs::Vector{<:AbstractGTBasisFuncs},
                     nuc::Vector{String}, nucCoords::Vector{<:AbstractArray};
                     X=getX(S), forUHF::Bool=false,
                     scfConfig=SCFconfig([:ADIIS], [1e-10])) where
                    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    D₁ = zero(Hcore)
    D₂ = zero(Hcore)
    N₁tot = 0
    N₂tot = 0
    order = sortperm(nuc, by=x->AtomicNumberList[x])
    for (atm, coord) in zip(nuc[order], nucCoords[order])
        N = getCharge(atm)
        N₁ = N ÷ 2
        N₂ = N - N₁
        if N₂ > N₁ && N₂tot > N₁tot
            temp = N₁
            N₁ = N₂
            N₂ = temp
        end
        h1 = coreH(bs, [atm], [coord])
        res = runHFcore(scfConfig,
                        (N₁, N₂), h1, HeeI, S, X, getCfromHcore(X, h1, forUHF=true))
        D₁ += res.D[1]
        D₂ += res.D[2]
        N₁tot += N₁
        N₂tot += N₂
    end
    Dᵀ = D₁ + D₂
    if forUHF
        getC.(Ref(X), getF.(Ref(Hcore), Ref(HeeI), (D₁, D₂), Ref(Dᵀ)))
    else
        getC(X, getF(Hcore, HeeI, Dᵀ.*0.5, Dᵀ))
    end
end

const guessCmethods =
    Dict(  :GWH => (forUHF, S, X, Hcore, _...) -> getCfromGWH(S, Hcore; X, forUHF),
         :Hcore => (forUHF, S, X, Hcore, _...) -> getCfromHcore(X, Hcore; forUHF),
           :SAD => (forUHF, S, X, Hcore, HeeI, bs, nuc, nucCoords) ->
                   getCfromSAD(S, Hcore, HeeI, bs, nuc, nucCoords; X, forUHF))

guessC(method::Symbol, forUHF::Bool, S::Matrix{T1}, X::Matrix{T1}, Hcore::Matrix{T2},
       HeeI::Array{T3, 4}, bs::Vector{<:AbstractGTBasisFuncs},
       nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) where
      {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB} =
guessCmethods[method](forUHF, S, X, Hcore, HeeI, bs, nuc, nucCoords)

getD(C::Matrix{T}, Nˢ::Int) where {TelLB<:T<:TelUB} =
@views (C[:,1:Nˢ]*C[:,1:Nˢ]') |> Hermitian |> Array
# Nˢ: number of electrons with the same spin.

getD(X::Matrix{T1}, F::Matrix{T2}, Nˢ::Int) where {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB} =
getD(getC(X, F), Nˢ)

function getGcore(HeeI::Array{T1, 4}, DJ::Matrix{T2}, DK::Matrix{T3}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    G = zero(DJ)
    l = size(G)[1]
    for μ = 1:l, ν = 1:l # fastest
        G[μ, ν] = dot(transpose(DJ), @view HeeI[μ,ν,:,:]) - dot(DK, @view HeeI[μ,:,:,ν])
    end
    G |> Hermitian |> Array
end

# RHF
getG(HeeI::Array{T1, 4}, D::Matrix{T2}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB} =
getGcore(HeeI, 2D, D)

# UHF
getG(HeeI::Array{T1, 4}, D::Matrix{T2}, Dᵀ::Matrix{T3}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB} =
getGcore(HeeI, Dᵀ, D)

getF(Hcore::Matrix{T1}, G::Matrix{T2}) where
            {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB} =
        Hcore + G

# RHF
getF(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, D::Matrix{T3}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB} =
getF(Hcore, getG(HeeI, D))

# UHF
getF(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, D::Matrix{T3}, Dᵀ::Matrix{T4}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB} =
getF(Hcore, getG(HeeI, D, Dᵀ))

# RHF or UHF
getF(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, Ds::NTuple{N, Matrix{T3}}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, N} =
getF(Hcore, getG(HeeI, Ds...))

getE(Hcore::Matrix{T1}, F::Matrix{T2}, D::Matrix{T3}) where
    {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB} =
dot(transpose(D), 0.5*(Hcore + F))

# RHF
getEᵀcore(Hcore::Matrix{T1}, F::Matrix{T2}, D::Matrix{T3}) where
         {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB} =
2*getE(Hcore, F, D)

# UHF
getEᵀcore(Hcore::Matrix{T1}, Fᵅ::Matrix{T2},   Dᵅ::Matrix{T3},
          Fᵝ::Matrix{T4},    Dᵝ::Matrix{T5}) where
         {TelLB<:T1<:TelUB,  TelLB<:T2<:TelUB, TelLB<:T3<:TelUB,
          TelLB<:T4<:TelUB,  TelLB<:T5<:TelUB} =
getE(Hcore, Fᵅ, Dᵅ) + getE(Hcore, Fᵝ, Dᵝ)

# RHF
function getEᵀ(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, C::Matrix{T3}, N::Int) where
              {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    D = getD(C, N÷2)
    F = getF(Hcore, HeeI, D)
    getEᵀcore(Hcore, F, D)
end

# UHF
function getEᵀ(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, (Cᵅ,Cᵝ)::NTuple{2, Matrix{T3}},
               (Nᵅ,Nᵝ)::NTuple{2, Int}) where
              {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    Dᵅ = getD(Cᵅ, Nᵅ)
    Dᵝ = getD(Cᵝ, Nᵝ)
    Dᵀ = Dᵅ + Dᵝ
    Fᵅ = getF(Hcore, HeeI, Dᵅ, Dᵀ)
    Fᵝ = getF(Hcore, HeeI, Dᵝ, Dᵀ)
    getEᵀcore(Hcore, Fᵅ, Dᵅ, Fᵝ, Dᵝ)
end

function getCFDE(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, X::Matrix{T3},
                 Ds::Vararg{Matrix{T4}, N}) where
                {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB, N}
    Fnew = getF(Hcore, HeeI, Ds)
    Enew = getE(Hcore, Fnew, Ds[1])
    Cnew = getC(X, Fnew)
    [Cnew, Fnew, Ds[1], Enew] # Fnew is based on latest variables.
end

# RHF
function initializeSCF(Hcore::Matrix{T1}, HeeI::Array{T2, 4}, C::Matrix{T3}, N::Int) where
                      {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    Nˢ = N÷2
    D = getD(C, Nˢ)
    F = getF(Hcore, HeeI, D)
    E = getE(Hcore, F, D)
    HFtempVars(:RHF, Nˢ, C, F, D, E, 2D, 2E)
end

# UHF
function initializeSCF(Hcore::Matrix{T1}, HeeI::Array{T2, 4},
                       Cs::NTuple{2, Matrix{T3}},
                       Ns::NTuple{2, Int}) where
                       {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    Ds = getD.(Cs, Ns)
    Dᵀs = [Ds |> sum]
    Fs = getF.(Ref(Hcore), Ref(HeeI), Ds, Ref(Dᵀs[]))
    Es = getE.(Ref(Hcore), Fs, Ds)
    Eᵀs = [Es |> sum]
    res = HFtempVars.(:UHF, Ns, Cs, Fs, Ds, Es)
    res[1].shared.Dtots = res[2].shared.Dtots = Dᵀs
    res[1].shared.Etots = res[2].shared.Etots = Eᵀs
    res |> Tuple
end

const Doc_SCFconfig_OneRowTable = "| `:DIIS`, `:EDIIS`, `:ADIIS` | Subspace size (>1); "*
                                  "Coefficient solver(`:ADMM`-> ADMM solver,"*
                                  " `:LCM` -> Lagrange solver) | "*
                                  "`DIISsize::Int`; `solver::Symbol` | `15`; `:ADMM` |"

"""

    SCFconfig{N} <: ImmutableParameter{SCFconfig, Any}

The `struct` for SCF iteration configurations.

≡≡≡ Field(s) ≡≡≡

`methods::NTuple{N, Symbol}`: The applied methods. The available methods are their
configurations (in terms of keyword arguments):

| Methods | Configuration(s) | keyword argument(s) | Default value(s) |
| :----   | :---:            | :---:               | ----:            |
| `:DS` | Damping strength: [0,1] | `dampingStrength::Float64` | `0.0` |
$(Quiqbox.Doc_SCFconfig_OneRowTable)

`intervals`: The stopping (skipping) thresholds for the required methods.

`methodConfigs`: The additional keywords arguments for each method stored as `Tuple`s of
`Pair`s.

`oscillateThreshold`: The threshold for oscillating convergence.

≡≡≡ Initialization Method(s) ≡≡≡

    SCFconfig(methods::Vector{Symbol}, intervals::Vector{Float64},
              configs::Dict{Int, <:Vector{<:Pair}}=Dict(1=>Pair[]);
              oscillateThreshold::Float64=1e-5) ->
    SCFconfig{N}

`methods` and `intervals` are the methods to be applied and their stopping (skipping)
thresholds respectively; the length of those two `Vector`s should be the same. `configs`
specifies the additional keyword arguments for each methods by a `Pair` of which the `Int`
key `i` is for `i`th method and the pointed `Vector{<:Pair}` is the pairs of keyword
arguments and their values respectively.

≡≡≡ Example(s) ≡≡≡

julia> SCFconfig([:SD, :ADIIS, :DIIS], [1e-4, 1e-12, 1e-13],
                 Dict(2=>[:solver=>:LCM])
SCFconfig{3}((:SD, :ADIIS, :DIIS), (0.0001, 1.0e-12, 1.0e-13), ((), (:solver => :LCM,),
()), 1.0e-5)
"""
struct SCFconfig{N} <: ImmutableParameter{SCFconfig, Any}
    methods::NTuple{N, Symbol}
    intervals::NTuple{N, Float64}
    methodConfigs::NTuple{N, Tuple{Vararg{Pair}}}
    oscillateThreshold::Float64

    function SCFconfig(methods::Vector{Symbol}, intervals::Vector{Float64},
                       configs::Dict{Int, <:Vector{<:Pair}}=Dict(1=>Pair[]);
                       oscillateThreshold::Float64=1e-5)
        l = length(methods)
        kwPairs = [Pair[] for _=1:l]
        for i in keys(configs)
            kwPairs[i] = configs[i]
        end
        new{length(methods)}(Tuple(methods), Tuple(intervals), Tuple(kwPairs .|> Tuple),
                             oscillateThreshold)
    end
end

const defaultSCFconfig = SCFconfig([:ADIIS, :DIIS, :ADIIS], [1e-4, 1e-6, 1e-15])

mutable struct HFinterrelatedVars <: HartreeFockintermediateData
    Dtots::Vector{Matrix{T}} where {TelLB<:T<:TelUB}
    Etots::Vector{Float64}

    HFinterrelatedVars() = new()
    HFinterrelatedVars(Dtots, Etots) = new(Dtots, Etots)
end

"""
    HFtempVars{HFtype, N} <: HartreeFockintermediateData

The container to store the intermediate values (only of the same spin configuration) for
each iteration during the Hartree-Fock SCF procedure.

≡≡≡ Field(s) ≡≡≡

`Cs::Array{Array{T1, 2}, 1} where {$(TelLB)<:T1<:$(TelUB)}`: Coefficient matrices.

`Fs::Array{Array{T2, 2}, 1} where {$(TelLB)<:T2<:$(TelUB)}`: Fock matrices

`Ds::Array{Array{T3, 2}, 1} where {$(TelLB)<:T3<:$(TelUB)}`: Density matrices corresponding
to only spin configuration. For RHF each elements means (unconverged) 0.5*Dᵀ.

`Es::Array{Float64, 1}`: Part of Hartree-Fock energy corresponding to only spin
configuration. For RHF each element means (unconverged) 0.5*E0HF.

`shared.Dtots::Array{Array{T, 2}, 1} where {$(TelLB)<:T<:$(TelUB)}`: The total density
matrices.

`shared.Etots::Array{Float64, 1}`: The total Hartree-Fock energy.

**NOTE: For UHF, there are 2 `HFtempVars` being updated during the SCF iterations, and
change the field `shared.Dtots` or `shared.Etots` of one container will affect the other
one's.**
"""
struct HFtempVars{HFtype, N} <: HartreeFockintermediateData
    Cs::Vector{Matrix{T1}} where {TelLB<:T1<:TelUB}
    Fs::Vector{Matrix{T2}} where {TelLB<:T2<:TelUB}
    Ds::Vector{Matrix{T3}} where {TelLB<:T3<:TelUB}
    Es::Vector{Float64}
    shared::HFinterrelatedVars

    HFtempVars(HFtype::Symbol, N, C, F, D, E, vars...) =
    new{HFtype, N}([C], [F], [D], [E], HFinterrelatedVars([[x] for x in vars]...))

    HFtempVars(HFtype::Symbol, Nˢ::Int, Cs::Vector{Matrix{T1}}, Fs::Vector{Matrix{T2}},
               Ds::Vector{Matrix{T3}}, Es::Vector{<:Real}, Dtots::Vector{Matrix{T4}},
               Etots::Vector{<:Real}) where
              {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB} =
    new{HFtype, Nˢ}(Cs, Fs, Ds, Es, HFinterrelatedVars(Dtots, Etots))
end

"""

    HFfinalVars{T, N, Nb} <: HartreeFockFinalValue{T}

The container of the final values after a Hartree-Fock SCF procedure.

≡≡≡ Field(s) ≡≡≡

`E0HF::Float64`: Hartree-Fock energy of the electronic Hamiltonian.

`C::Union{Array{T1, 2}, NTuple{2, Array{T1, 2}}} where {$(TelLB)<:T1<:$(TelUB)}`:
Coefficient matrix(s) for one spin configuration.

`F::Union{Array{T2, 2}, NTuple{2, Array{T2, 2}}} where {$(TelLB)<:T2<:$(TelUB)}`: Fock
matrix(s) for one spin configuration.

`D::Union{Array{T3, 2}, NTuple{2, Array{T3, 2}}} where {$(TelLB)<:T3<:$(TelUB)}`: Density
matrix(s) for one spin configuration.

`Emo::Union{Array{Float64, 1}, NTuple{2, Array{Float64, 1}}}`: Energies of molecular
orbitals.

`occu::Union{Array{Int, 1}, NTuple{2, Array{Int, 1}}}`: occupation numbers of molecular
orbitals.

`temp::Union{HFtempVars{T}, NTuple{2, HFtempVars{T}}}` the intermediate values.

`isConverged::Bool`: Whether the SCF procedure is converged in the end.
"""
struct HFfinalVars{T, N, Nb} <: HartreeFockFinalValue{T}
    E0HF::Float64
    C::Union{Matrix{T1}, NTuple{2, Matrix{T1}}} where {TelLB<:T1<:TelUB}
    F::Union{Matrix{T2}, NTuple{2, Matrix{T2}}} where {TelLB<:T2<:TelUB}
    D::Union{Matrix{T3}, NTuple{2, Matrix{T3}}} where {TelLB<:T3<:TelUB}
    Emo::Union{Vector{Float64}, NTuple{2, Vector{Float64}}}
    occu::Union{Vector{Int}, NTuple{2, Vector{Int}}}
    temp::Union{HFtempVars{T}, NTuple{2, HFtempVars{T}}}
    isConverged::Bool

    function HFfinalVars(X::Matrix{T}, vars::HFtempVars{:RHF}, isConverged::Bool) where
             {TelLB<:T<:TelUB}
        C = vars.Cs[end]
        F = vars.Fs[end]
        D = vars.Ds[end]
        E0HF = vars.shared.Etots[end]
        _, Emo = getC(X, F, outputEmo=true)
        Nˢ = typeof(vars).parameters[2]
        occu = vcat(2*ones(Int, Nˢ), zeros(Int, size(X, 1) - Nˢ))
        temp = vars
        new{:RHF, 2Nˢ, length(Emo)}(E0HF, C, F, D, Emo, occu, temp, isConverged)
    end

    function HFfinalVars(X::Matrix{T}, (αVars, βVars)::NTuple{2, HFtempVars{:UHF}},
                         isConverged::Bool) where {TelLB<:T<:TelUB}
        C = (αVars.Cs[end], βVars.Cs[end])
        F = (αVars.Fs[end], βVars.Fs[end])
        D = (αVars.Ds[end], βVars.Ds[end])
        E0HF = αVars.shared.Etots[end]
        res = getC.(Ref(X), F, outputEmo=true)
        Emo = getindex.(res, 2)
        Nˢs = (typeof(αVars).parameters[2], typeof(βVars).parameters[2])
        occu = vcat.(ones.(Int, Nˢs), zeros.(Int, size(X, 1) .- Nˢs))
        temp = (αVars, βVars)
        new{:UHF, Nˢs |> sum, length(Emo[1])}(E0HF, C, F, D, Emo, occu, temp, isConverged)
    end
end

"""
    runHF(bs::Union{BasisSetData, Array{<:AbstractGTBasisFuncs, 1}},
          nuc::Array{String, 1},
          nucCoords::Array{<:AbstractArray, 1},
          N::Int=getCharge(nuc);
          initialC::Union{Array{T, 2}, NTuple{2, Array{T, 2}}, Symbol}=:SAD,
          HFtype::Symbol=:RHF,
          scfConfig::SCFconfig=defaultSCFconfig,
          earlyTermination::Bool=true,
          printInfo::Bool=true,
          maxSteps::Int=1000) where {$(TelLB)<:T<:$(TelUB)} -> HFfinalVars

Main function to run Hartree-Fock in Quiqbox.

=== Positional argument(s) ===

`bs::Union{BasisSetData, Array{<:AbstractGTBasisFuncs, 1}}`: Basis set.

`nuc::Array{String, 1}`: The element symbols of the nuclei for the Molecule.

`nucCoords::Array{<:AbstractArray, 1}`: The coordinates of the nuclei.

`N::Int`: The total number of electrons.

=== Keyword argument(s) ===

`initialC::Union{Array{T, 2}, NTuple{2, Array{T, 2}}, Symbol}`: Initial guess of the
coefficient matrix(s) C of the molecular orbitals.

`HFtype::Symbol`: Hartree-Fock type. Available values are `:RHF` and `:UHF`.

`scfConfig::SCFconfig`: SCF iteration configuration.

`earlyTermination::Bool`: Whether automatically early terminate (skip) a convergence method
when its performance becomes unstable or poor.

`printInfo::Bool`: Whether print out the information of each iteration step.

`maxSteps::Int`: Maximum allowed iteration steps regardless of whether the SCF converges.
"""
function runHF(bs::Vector{<:AbstractGTBasisFuncs},
               nuc::Vector{String},
               nucCoords::Vector{<:AbstractArray},
               N::Int=getCharge(nuc);
               initialC::Union{Matrix{T}, NTuple{2, Matrix{T}}, Symbol}=:SAD,
               HFtype::Symbol=:RHF,
               scfConfig::SCFconfig=defaultSCFconfig,
               earlyTermination::Bool=true,
               printInfo::Bool=true,
               maxSteps::Int=1000) where {TelLB<:T<:TelUB}
    gtb = GTBasis(bs, false)
    runHF(gtb, nuc, nucCoords, N; initialC, HFtype, scfConfig,
          earlyTermination, printInfo, maxSteps)
end

function runHF(gtb::BasisSetData,
               nuc::Vector{String},
               nucCoords::Vector{<:AbstractArray},
               N::Int=getCharge(nuc);
               initialC::Union{Matrix{T}, NTuple{2, Matrix{T}}, Symbol}=:SAD,
               HFtype::Symbol=:RHF,
               scfConfig::SCFconfig=defaultSCFconfig,
               earlyTermination::Bool=true,
               printInfo::Bool=true,
               maxSteps::Int=1000) where {TelLB<:T<:TelUB}
    @assert length(nuc) == length(nucCoords)
    @assert typeof(gtb).parameters[1] >= ceil(sum(N)/2)
    @assert N > (HFtype == :RHF) "$(HFtype) requires more than $(LBofN) electron to run."
    HFtype == :UHF && (N = (N÷2, N-N÷2))
    Hcore = gtb.getHcore(nuc, nucCoords)
    X = getX(gtb.S)
    initialC isa Symbol && (initialC = guessC(initialC, (HFtype == :UHF), gtb.S, X, Hcore,
                                              gtb.eeI, gtb.basis, nuc, nucCoords))
    runHFcore(scfConfig, N, Hcore, gtb.eeI, gtb.S, X, initialC;
              printInfo, maxSteps, earlyTermination)
end

"""

    runHFcore(scfConfig::SCFconfig{L},
              N::Union{NTuple{2, Int}, Int},
              Hcore::Array{T1, 2},
              HeeI::Array{T2, 4},
              S::Array{T3, 2},
              X::Array{T4, 2}=getX(S),
              C::Union{Array{T, 2}, NTuple{2, Array{T, 2}}}=
              getCfromGWH(S, Hcore; X, forUHF=(length(N)==2));
              earlyTermination::Bool=true,
              printInfo::Bool=false,
              maxSteps::Int=1000) where {$(TelLB)<:T1<:$(TelUB),
                                         $(TelLB)<:T2<:$(TelUB),
                                         $(TelLB)<:T3<:$(TelUB),
                                         $(TelLB)<:T4<:$(TelUB),
                                         $(TelLB)<:T5<:$(TelUB), L}

The core function of `runHF`.

=== Positional argument(s) ===

`scfConfig::SCFconfig`: SCF iteration configuration.

`N::Union{NTuple{2, Int}, Int}`: The total number of electrons or the numbers of electrons
with different spins respectively. When the latter is input, an UHF is performed.

`Hcore::Array{T1, 2}`: Core Hamiltonian of electronic Hamiltonian.

`HeeI::Array{T2, 4}`: The electron-electron interaction Hamiltonian which includes both the
Coulomb interactions and the Exchange Correlations.

`S::Array{T3, 2}`: Overlap matrix of the corresponding basis set.

`X::Array{T4, 2}`: Orthogonal transformation matrix of S. Default value is S^(-0.5).

`C::Union{Array{T, 2}, NTuple{2, Array{T, 2}}}`: Initial guess of the coefficient matrix(s)
C of the molecular orbitals.

=== Keyword argument(s) ===

`earlyTermination::Bool`: Whether automatically early terminate (skip) a convergence method
when its performance becomes unstable or poor.

`printInfo::Bool`: Whether print out the information of each iteration step.

`maxSteps::Int`: Maximum allowed iteration steps regardless of whether the SCF converges.
"""
function runHFcore(scfConfig::SCFconfig{L},
                   N::Union{NTuple{2, Int}, Int},
                   Hcore::Matrix{T1},
                   HeeI::Array{T2, 4},
                   S::Matrix{T3},
                   X::Matrix{T4}=getX(S),
                   C::Union{Matrix{T5}, NTuple{2, Matrix{T5}}}=
                   getCfromGWH(S, Hcore; X, forUHF=(length(N)==2));
                   earlyTermination::Bool=true,
                   printInfo::Bool=false,
                   maxSteps::Int=1000) where {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB,
                                              TelLB<:T3<:TelUB, TelLB<:T4<:TelUB,
                                              TelLB<:T5<:TelUB, L}
    @assert maxSteps > 0
    vars = initializeSCF(Hcore, HeeI, C, N)
    Etots = (vars isa Tuple) ? vars[1].shared.Etots : vars.shared.Etots
    HFtypeStr =  length(N) == 1 ? "RHF" : "UHF"
    printInfo && println(rpad(HFtypeStr*" | Initial Gauss", 22), "E = $(Etots[end])")
    isConverged = true
    EtotMin = Etots[]
    for (method, kws, breakPoint, i) in
        zip(scfConfig.methods, scfConfig.methodConfigs, scfConfig.intervals, 1:L)

        while true
            iStep = length(Etots)
            iStep <= maxSteps || (isConverged = false) || break
            HF!(method, N, Hcore, HeeI, S, X, vars; kws...)
            printInfo && (iStep % floor(log(4, iStep) + 1) == 0 || iStep == maxSteps) &&
            println(rpad("Step $(iStep)", 10),
                    rpad("#$(i) ($(method))", 12),
                    "E = $(Etots[end])")
            abs(Etots[end]-Etots[end-1]) > breakPoint || (isConverged = true) && break
            flag, Std = isOscillateConverged(Etots,
                                             10^(log(10, breakPoint)÷2),
                                             returnStd=true)
            flag && (isConverged = Std > scfConfig.oscillateThreshold ? false : true; break)
            earlyTermination && (Etots[end] - EtotMin) / abs(EtotMin) > 0.2 &&
            (printInfo && println("Early termination of ", method,
                                  " due to the poor performance.");
             isConverged = false; break)
        end

    end
    negStr = isConverged ? "is " : "has not "
    printInfo && println("The SCF procedure ", negStr, "converged.\n")
    HFfinalVars(X, vars, isConverged)
end

function SDcore(Nˢ::Int, Hcore::Matrix{T1}, HeeI::Array{T2, 4}, X::Matrix{T3},
                F::Matrix{T4}, D::Matrix{T5}; dampingStrength::Float64=0.0, _kws...) where
               {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB,
                TelLB<:T5<:TelUB}
    @assert 0 <= dampingStrength <= 1 "The range of `dampingStrength`::Float64 is [0,1]."
    Dnew = getD(X, F, Nˢ)
    (1-dampingStrength)*Dnew + dampingStrength*D
end

function xDIIScore(method::Symbol, Nˢ::Int, Hcore::Matrix{T1}, HeeI::Array{T2, 4},
                   S::Matrix{T3}, X::Matrix{T4}, Fs::Vector{Matrix{T5}},
                   Ds::Vector{Matrix{T6}},Es::Vector{Float64};
                   DIISsize::Int=15, solver::Symbol=:ADMM,
                   _kws...) where
                  {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB,
                   TelLB<:T5<:TelUB, TelLB<:T6<:TelUB}
    DIISmethod, convexConstraint, permuteData = DIISmethods[method]
    is = permuteData ? sortperm(Es) : (:)
    ∇s = (@view Fs[is])[1:end .> end-DIISsize]
    Ds = (@view Ds[is])[1:end .> end-DIISsize]
    Es = (@view Es[is])[1:end .> end-DIISsize]
    vec, B = DIISmethod(∇s, Ds, Es, S)
    c = constraintSolver(vec, B, solver; convexConstraint)
    grad = c.*∇s |> sum
    getD(X, grad |> Hermitian |> Array, Nˢ) # grad == F.
end

const DIISmethods = Dict( :DIIS => ((∇s, Ds,  _, S)->DIIScore(∇s, Ds, S),   false, true),
                         :EDIIS => ((∇s, Ds, Es, _)->EDIIScore(∇s, Ds, Es), true, false),
                         :ADIIS => ((∇s, Ds,  _, _)->ADIIScore(∇s, Ds),     true, false))

function EDIIScore(∇s::Vector{Matrix{T1}}, Ds::Vector{Matrix{T2}},
                   Es::Vector{Float64}) where {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB}
    len = length(Ds)
    B = ones(len, len)
    for i=1:len, j=1:len
        B[i,j] = -dot(Ds[i]-Ds[j], ∇s[i]-∇s[j])
    end
    Es, B
end

function ADIIScore(∇s::Vector{Matrix{T1}}, Ds::Vector{Matrix{T2}}) where
                  {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB}
    len = length(Ds)
    B = ones(len, len)
    vec = [dot(D - Ds[end], ∇s[end]) for D in Ds]
    for i=1:len, j=1:len
        B[i,j] = dot(Ds[i]-Ds[len], ∇s[j]-∇s[len])
    end
    vec, B
end

function DIIScore(∇s::Vector{Matrix{T1}}, Ds::Vector{Matrix{T2}}, S::Matrix{T3}) where
                 {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}
    len = length(Ds)
    B = ones(len, len)
    vec = zeros(len)
    for i=1:len, j=1:len
        B[i,j] = dot(∇s[i]*Ds[i]*S - Ds[i]*S*∇s[i], ∇s[j]*Ds[j]*S - Ds[j]*S*∇s[j])
    end
    vec, B
end

const SD = ((Nˢ::Int, Hcore::Matrix{T1}, HeeI::Array{T2, 4}, _dm::Any, X::Matrix{T3},
             tVars::HFtempVars; kws...) where
            {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB}) ->
      SDcore(Nˢ, Hcore, HeeI, X, tVars.Fs[end], tVars.Ds[end]; kws...)

const xDIIS = (method::Symbol) ->
              ((Nˢ::Int, Hcore::Matrix{T1}, HeeI::Array{T2, 4}, S::Matrix{T3},
                X::Matrix{T4}, tVars::HFtempVars; kws...) where
               {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB,
                TelLB<:T3<:TelUB, TelLB<:T4<:TelUB}) ->
              xDIIScore(method, Nˢ, Hcore, HeeI, S, X, tVars.Fs, tVars.Ds, tVars.Es; kws...)

const SCFmethods = [:SD, :DIIS, :ADIIS, :EDIIS]

const SCFmethodSelector = Dict(SCFmethods .=>
                               [SD, xDIIS(:DIIS), xDIIS(:ADIIS), xDIIS(:EDIIS)])

function HF!(SCFmethod::Symbol, N::Union{NTuple{2, Int}, Int},
             Hcore::Matrix{T1}, HeeI::Array{T2, 4}, S::Matrix{T3}, X::Matrix{T4},
             tVars::Union{HFtempVars{:RHF}, NTuple{2, HFtempVars{:UHF}}}; kws...) where
            {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB}
    res = HFcore(SCFmethod, N, Hcore, HeeI, S, X, tVars; kws...)
    pushHFtempVars!(tVars, res)
end

# RHF
function HFcore(SCFmethod::Symbol, N::Int,
                Hcore::Matrix{T1}, HeeI::Array{T2, 4}, S::Matrix{T3}, X::Matrix{T4},
                rVars::HFtempVars{:RHF}; kws...) where
               {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB}
    D = SCFmethodSelector[SCFmethod](N÷2, Hcore, HeeI, S, X, rVars; kws...)
    partRes = getCFDE(Hcore, HeeI, X, D)
    partRes..., 2D, 2partRes[end]
end

function pushHFtempVars!(rVars::HFtempVars,
                         res::Tuple{Matrix{T1}, Matrix{T2}, Matrix{T3}, Float64,
                                    Matrix{T4}, Float64}) where
                        {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB,
                         TelLB<:T4<:TelUB}
    push!(rVars.Cs, res[1])
    push!(rVars.Fs, res[2])
    push!(rVars.Ds, res[3])
    push!(rVars.Es, res[4])
    push!(rVars.shared.Dtots, res[5])
    push!(rVars.shared.Etots, res[6])
end

# UHF
function HFcore(SCFmethod::Symbol, Ns::NTuple{2, Int},
                Hcore::Matrix{T1}, HeeI::Array{T2, 4}, S::Matrix{T3}, X::Matrix{T4},
                uVars::NTuple{2, HFtempVars{:UHF}}; kws...) where
               {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB, TelLB<:T4<:TelUB}
    Ds = SCFmethodSelector[SCFmethod].(Ns, Ref(Hcore), Ref(HeeI),
                                Ref(S), Ref(X), uVars; kws...)
    Dᵀnew = Ds |> sum
    partRes = getCFDE.(Ref(Hcore), Ref(HeeI), Ref(X), Ds, Ref(Dᵀnew))
    Eᵀnew = partRes[1][end] + partRes[2][end]
    (partRes[1]..., Dᵀnew, Eᵀnew), (partRes[2]..., Dᵀnew, Eᵀnew)
end

function pushHFtempVars!(uVars::NTuple{2, HFtempVars{:UHF}},
                         res::NTuple{2, Tuple{Matrix{T1}, Matrix{T2}, Matrix{T3}, Float64,
                              Matrix{T4}, Float64}}) where
                        {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB, TelLB<:T3<:TelUB,
                         TelLB<:T4<:TelUB}
    pushHFtempVars!.(uVars, res)
    pop!(uVars[1].shared.Dtots)
    pop!(uVars[1].shared.Etots)
end

function popHFtempVars!(rVars::HFtempVars)
    pop!(rVars.Cs)
    pop!(rVars.Fs)
    pop!(rVars.Ds)
    pop!(rVars.Es)
    pop!(rVars.shared.Dtots)
    pop!(rVars.shared.Etots)
end

function popHFtempVars!(uVars::NTuple{2, HFtempVars{:UHF}})
    fields = [:Cs, :Fs, :Ds, :Es]
    for field in fields pop!.(getfield.(uVars, field)) end
    pop!(uVars[1].shared.Dtots)
    pop!(uVars[1].shared.Etots)
end

# Default
function ADMMSolver(vec::Vector{Float64}, B::Matrix{Float64}; convexConstraint::Bool=true)
    len = length(vec)
    A = ones(len) |> transpose |> Array
    b = Float64[1.0]
    g = convexConstraint ? fill(indicator(0, 1), len) : fill(indicator(-Inf, Inf), len)
    params = SeparableOptimization.AdmmParams(B, vec, A, b, g)
    settings = SeparableOptimization.Settings(; ρ=ones(1), σ=ones(len), compute_stats=true)
    vars, _ = SeparableOptimization.optimize(params, settings)
    vars.x
end

function CMSolver(vec::Vector, B::Matrix; convexConstraint=true, ϵ::Float64=1e-5)
    len = length(vec)
    getA = (B)->[B  ones(len); ones(1, len) 0]
    b = vcat(-vec, 1)
    local c
    while true
        A = getA(B)
        while det(A) == 0
            B += ϵ*I
            A = getA(B)
        end
        x = A \ b
        c = x[1:end-1]
        (findfirst(x->x<0, c) !== nothing && convexConstraint) || (return c)
        idx = (sortperm(abs.(c)) |> powerset |> collect)
        popfirst!(idx)

        for is in idx
            Atemp = A[1:end .∉ Ref(is), 1:end .∉ Ref(is)]
            btemp = b[begin:end .∉ Ref(is)]
            det(Atemp) == 0 && continue
            xtemp = Atemp \ btemp
            c = xtemp[1:end-1]
            for i in sort(is)
                insert!(c, i, 0.0)
            end
            (findfirst(x->x<0, c) !== nothing) || (return c)
        end

        B += ϵ*I
    end
    c
end

const ConstraintSolvers = Dict(:ADMM=>ADMMSolver, :LCM=>CMSolver)

constraintSolver(vec::Vector{T1}, B::Matrix{T2},
                 solver::Symbol=:ADMM; convexConstraint::Bool=true) where
                {TelLB<:T1<:TelUB, TelLB<:T2<:TelUB} =
ConstraintSolvers[solver](vec, B; convexConstraint)


text:
#####
##### This file provides functions that conditionally-evaluate interpolation operators
##### near boundaries in bounded directions.
#####
##### For example, the function _symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, c) either
#####
#####     1. Always returns symmetric_interpolate_xᶠᵃᵃ if the x-direction is Periodic; or
#####
#####     2. Returns symmetric_interpolate_xᶠᵃᵃ if the x-direction is Bounded and index i is not
#####        close to the boundary, or a second-order interpolation if i is close to a boundary.
#####

using Oceananigans.Grids: AbstractUnderlyingGrid, Bounded

const AUG = AbstractUnderlyingGrid

# Left-biased buffers are smaller by one grid point on the right side; vice versa for right-biased buffers
                                                                              # outside left | outside right buffer
@inline    outside_symmetric_buffer(i, N, ::AbstractAdvectionScheme{Nᴮ}) where Nᴮ = i > Nᴮ     && i < N + 1 - Nᴮ
@inline  outside_left_biased_buffer(i, N, ::AbstractAdvectionScheme{Nᴮ}) where Nᴮ = i > Nᴮ     && i < N + 1 - (Nᴮ - 1)
@inline outside_right_biased_buffer(i, N, ::AbstractAdvectionScheme{Nᴮ}) where Nᴮ = i > Nᴮ - 1 && i < N + 1 - Nᴮ

const ADV = AbstractAdvectionScheme
const WVI = WENOVectorInvariant

for bias in (:symmetric, :left_biased, :right_biased)

    for (d, ξ) in enumerate((:x, :y, :z))

        code = [:ᵃ, :ᵃ, :ᵃ]

        for loc in (:ᶜ, :ᶠ)
            code[d] = loc
            second_order_interp = Symbol(:ℑ, ξ, code...)
            interp = Symbol(bias, :_interpolate_, ξ, code...)
            alt_interp = Symbol(:_, interp)

            # Simple translation for Periodic directions (fallback)
            @eval $alt_interp(i, j, k, grid::AUG, scheme::ADV, args...) = $interp(i, j, k, grid, scheme, args...)

            outside_buffer = Symbol(:outside_, bias, :_buffer)

            # Conditional high-order interpolation in Bounded directions
            if ξ == :x
                @eval begin
                    @inline $alt_interp(i, j, k, grid::AUG{FT, <:Bounded}, scheme::ADV, ψ) where FT =
                        ifelse($outside_buffer(i, grid.Nx, scheme),
                               $interp(i, j, k, grid, scheme, ψ),
                               $second_order_interp(i, j, k, grid, ψ))

                    @inline $alt_interp(i, j, k, grid::AUG{FT, <:Bounded}, scheme::WVI, ζ, VI, u, v) where FT =
                        ifelse($outside_buffer(i, grid.Nx, scheme),
                            $interp(i, j, k, grid, scheme, ζ, VI, u, v),
                            $second_order_interp(i, j, k, grid, ζ, u, v))
                end
            elseif ξ == :y
                @eval begin
                    @inline $alt_interp(i, j, k, grid::AUG{FT, TX, <:Bounded}, scheme::ADV, ψ) where {FT, TX} =
                        ifelse($outside_buffer(j, grid.Ny, scheme),
                               $interp(i, j, k, grid, scheme, ψ),
                               $second_order_interp(i, j, k, grid, ψ))

                    @inline $alt_interp(i, j, k, grid::AUG{FT, TX, <:Bounded}, scheme::WVI, ζ, VI, u, v) where {FT, TX} =
                        ifelse($outside_buffer(j, grid.Ny, scheme),
                               $interp(i, j, k, grid, scheme, ζ, VI, u, v),
                               $second_order_interp(i, j, k, grid, ζ, u, v))
                end
            elseif ξ == :z
                @eval begin
                    @inline $alt_interp(i, j, k, grid::AUG{FT, TX, TY, <:Bounded}, scheme::ADV, ψ) where {FT, TX, TY} =
                        ifelse($outside_buffer(k, grid.Nz, scheme),
                               $interp(i, j, k, grid, scheme, ψ),
                               $second_order_interp(i, j, k, grid, ψ))

                    @inline $alt_interp(i, j, k, grid::AUG{FT, TX, TY, <:Bounded}, scheme::WVI, ∂z, VI, u) where {FT, TX, TY} =
                        ifelse($outside_buffer(k, grid.Nz, scheme),
                                $interp(i, j, k, grid, scheme, ∂z, VI, u),
                                $second_order_interp(i, j, k, grid, ∂z, u))
                end
            end
        end
    end
end


text:
## Parallel Test
## ==============
using Sigma
import Sigma: SimplexBenchmark, weighted_mid_split, rand_partial_split,
              weighted_partial_split
import Sigma: dreal, z3, dreal3
import Sigma: SigmaSMT, SigmaAI
import Sigma: cond_sample_tlmh
import Sigma: runbenchmarks

SplitBenchmarks = [SimplexBenchmark(i,[:sample_distribution,
                                        :accumulative_KL,
                                        :total_time,])
                   for i = 5:5]

mh_captures = [:start_loop, :refinement_depth]
all_splits = [weighted_partial_split]

# SMT algorithms
SMTAlgorithms = [SigmaSMT(mh_captures, solver, sampler, nprocs, split)
  for nprocs = [1:40],
      solver = [dreal],
      split = all_splits,
      sampler = [cond_sample_tlmh]][:]

AIAlgorithms = [SigmaAI(mh_captures, sampler, nprocs, split)
  for nprocs = [1:40],
      split = all_splits,
      sampler = [cond_sample_tlmh]][:]

dimbenchmarks = SplitBenchmarks
dimalgorithms = vcat(AIAlgorithms,SMTAlgorithms)

runbenchmarks(dimalgorithms,dimbenchmarks)

  # function extractdata()
  #   ys = Float64[]
  #   xs = Int64[]
  #   for (key,value) in data
  #     if isa(key[1],Sigma.SigmaSMT)
  #       push!(xs,key[1].ncores)
  #       push!(ys,value[:total_time][1])
  #     end
  #   end
  # end

  # [1,2,3,4][2,3]
  # x = [37,31,1,8,22,28,24,33,17,25,19,4,40,27,18,2,5,9,3,11,14,10,30,35,13,34,21,36,16,39,7,26,12,29,23,38,32,20,6,15]
  # y = [0.784782175,0.663763885,4.684491724,0.679153679,0.631475561,0.667785573,0.851447618,0.719184689,0.974614338,0.742905691,0.698064856,0.696794102,0.694817954,0.812678691,0.773986244,1.991063103,0.716975545,0.795551611,1.776576016,0.664803417,0.753942407,0.731583818,0.67289623,0.881060887,0.6393919,0.725649725,0.906088166,0.736668336,0.65141427,0.704332105,0.794086698,0.76449543,0.665638965,0.795386601,0.715931824,0.759590093,0.793091758,0.743893172,0.907585122,0.646210437]

  # Dict(x,y)[2]
  # Dict(x,y)[3]
  # Dict(x,y)[4]

  # plot(x = [1,2,3,4,5], y = [Dict(x,y)[1],Dict(x,y)[2],Dict(x,y)[3],Dict(x,y)[4],Dict(x,y)[5]], Geom.line(), Guide.xlabel("Number of Processors")

  #     ,Guide.ylabel("Runtime(s)"))

  # # SMT
  # x = [33,40,17,11,25,34,26,28,19,10,38,2,39,13,29,22,3,37,6,1,7,15,27,4,18,36,35,32,31,5,30,16,23,12,21,14,8,20,24,9]
  # y = [29.046722726,30.595505191,31.956475104,32.494068589,31.76427943,31.770079567,30.993137568,30.178511233,30.431204351,33.392272495,32.249786349,45.953667995,30.785560356,33.241083416,29.967980785,29.683638232,29.945061056,30.622136539,33.0408161,48.719286872,34.90705489,30.945352803,30.708056486,31.934607426,30.088652276,28.413354981,30.677160721,30.861685656,28.987676502,32.061471006,30.513410012,30.595890113,27.883390746,32.454835351,30.295537073,33.004421108,31.591543316,36.806507415,30.000664329,31.45082343]

  # using DataFrames
  # smtdf = DataFrame(x=x, y = y, label = "SMT")
  # aidf = DataFrame(x=x, y = y, label = "AI")

  # def = vcat(smtdf,aidf)
  # smtplot = plot(smtdf, x="x", y="y", color="label",Geom.line, Scale.discrete_color_manual("red"),Guide.xlabel("Number of processors"),Guide.ylabel("Runtime (s)"))

  # draw(PDF("/home/zenna/Dropbox/writing/conferences/icml2015/figures/smt1.pdf",4*1inch,3*1inch), smtplot)

  # x = [1,2,3,4]
  # y = [Dict(x,y)[1],Dict(x,y)[2],Dict(x,y)[3],Dict(x,y)[4]]
  # y
  # x = [37,31,1,8,22,28,24,33,17,25,19,4,40,27,18,2,5,9,3,11,14,10,30,35,13,34,21,36,16,39,7,26,12,29,23,38,32,20,6,15]
  # y = [0.784782175,0.663763885,4.684491724,0.679153679,0.631475561,0.667785573,0.851447618,0.719184689,0.974614338,0.742905691,0.698064856,0.696794102,0.694817954,0.812678691,0.773986244,1.991063103,0.716975545,0.795551611,1.776576016,0.664803417,0.753942407,0.731583818,0.67289623,0.881060887,0.6393919,0.725649725,0.906088166,0.736668336,0.65141427,0.704332105,0.794086698,0.76449543,0.665638965,0.795386601,0.715931824,0.759590093,0.793091758,0.743893172,0.907585122,0.646210437]
  # lsmt = layer(x = [1,2,3,4], y = [Dict(x,y)[1],Dict(x,y)[2],Dict(x,y)[3],Dict(x,y)[4]],Geom.line)

  # lai
  # plot(lai,lsmt,Guide.ylabel("Runtime(s)"),Guide.xlabel("Runtime(s)"),Guide.title("Parallel Performance over 4 cores"))

  # q =plot(x = [37,31,1,8,22,28,24,33,17,25,19,4,40,27,18,2,5,9,3,11,14,10,30,35,13,34,21,36,16,39,7,26,12,29,23,38,32,20,6,15],
  #      y = [0.784782175,0.663763885,4.684491724,0.679153679,0.631475561,0.667785573,0.851447618,0.719184689,0.974614338,0.742905691,0.698064856,0.696794102,0.694817954,0.812678691,0.773986244,1.991063103,0.716975545,0.795551611,1.776576016,0.664803417,0.753942407,0.731583818,0.67289623,0.881060887,0.6393919,0.725649725,0.906088166,0.736668336,0.65141427,0.704332105,0.794086698,0.76449543,0.665638965,0.795386601,0.715931824,0.759590093,0.793091758,0.743893172,0.907585122,0.646210437],Geom.line)

  # plot(x = [33,40,17,11,25,34,26,28,19,10,38,2,39,13,29,22,3,37,6,1,7,15,27,4,18,36,35,32,31,5,30,16,23,12,21,14,8,20,24,9], y = [29.046722726,30.595505191,31.956475104,32.494068589,31.76427943,31.770079567,30.993137568,30.178511233,30.431204351,33.392272495,32.249786349,45.953667995,30.785560356,33.241083416,29.967980785,29.683638232,29.945061056,30.622136539,33.0408161,48.719286872,34.90705489,30.945352803,30.708056486,31.934607426,30.088652276,28.413354981,30.677160721,30.861685656,28.987676502,32.061471006,30.513410012,30.595890113,27.883390746,32.454835351,30.295537073,33.004421108,31.591543316,36.806507415,30.000664329,31.45082343],Geom.line)

  # # #SMT
  # # plot(x=[18,12,17,37,33,21,35,40,28,6,13,3,31,23,20,11,39,34,29,32,30,15,24,22,26,7,19,16,2,36,14,8,25,38,10,4,27,5,9,1],
  # #      y= [30.088652276,32.454835351,31.956475104,30.622136539,29.046722726,30.295537073,30.677160721,30.595505191,30.178511233,33.0408161,33.241083416,29.945061056,28.987676502,27.883390746,36.806507415,32.494068589,30.785560356,31.770079567,29.967980785,30.861685656,30.513410012,30.945352803,30.000664329,29.683638232,30.993137568,34.90705489,30.431204351,30.595890113,45.953667995,28.413354981,33.004421108,31.591543316,31.76427943,32.249786349,33.392272495,31.934607426,30.708056486,32.061471006,31.45082343,48.719286872], Geom.line)
  # # plot(x = [31,9,23,10,17,16,35,5,40,27,18,4,34,32,12,7,2,11,8,1,13,3,33,26,21,28,14,38,22,36,39,25,24,15,30,37,19,20,6,29], y = [0.663763885,0.795551611,0.715931824,0.731583818,0.974614338,0.65141427,0.881060887,0.716975545,0.694817954,0.812678691,0.773986244,0.696794102,0.725649725,0.793091758,0.665638965,0.794086698,1.991063103,0.664803417,0.679153679,4.684491724,0.6393919,1.776576016,0.719184689,0.76449543,0.906088166,0.667785573,0.753942407,0.759590093,0.631475561,0.736668336,0.704332105,0.742905691,0.851447618,0.646210437,0.67289623,0.784782175,0.698064856,0.743893172,0.907585122,0.795386601],Geom.line)


text:
mutable struct WorkSpace
    sdg::Graphs.SequenceDistanceGraph{LongDNASeq}
    paired_reads_datastores::Vector{Union{PairedReads{DNAAlphabet{2}},PairedReads{DNAAlphabet{4}}}}
    long_reads_datastores::Vector{LongReads}
    linked_reads_datastores::Vector{LinkedReads}
    mer_count_stores::Vector{IndexedCounts}
end

"Create an empty workspace"
function WorkSpace()
    return WorkSpace(Graphs.SequenceDistanceGraph{LongDNASeq}(),
                     Vector{PairedReads}(),
                     Vector{LongReads}(),
                     Vector{LinkedReads}(),
                     Vector{IndexedCounts}())
end

graph(ws::WorkSpace) = ws.sdg

"""
    paired_reads(ws::WorkSpace, id::Symbol)

Get a reference to the `PairedReads` in the workspace that has the name
specified by `nm`.
"""
function paired_reads(ws::WorkSpace, nm::Symbol)
    for ds in ws.paired_reads_datastores
        if name(ds) == nm
            return ds
        end
    end
    error(string("WorkSpace has no paired read datastore called ", nm))
end

function add_paired_reads!(ws::WorkSpace, reads::PairedReads)
    push!(ws.paired_reads_datastores, reads)
    return ws
end

function add_paired_reads!(ws::WorkSpace, file::String, name::Union{String,Nothing} = nothing)
    reads = open(PairedReads, file, name)
    add_paired_reads!(ws, reads)
    return ws
end

function add_mer_counts!(ws::WorkSpace, ::Type{M}, mode::KmerCountMode, name::Symbol) where {M<:AbstractMer}
    push!(ws.mer_count_stores, MerCounts{M}(name, ws.sdg, mode))
    return ws
end

function mer_counts(ws::WorkSpace, name::Symbol)
    for mc in ws.mer_count_stores
        if mc.name == name
            return mc
        end
    end
    error(string("WorkSpace has no mer count datastore datastore called ", name))
end

function Base.show(io::IO, ws::WorkSpace)
    println("Graph Genome workspace")

    println(io, ' ', summary(ws.sdg))

    println(" Paired Read Datastores (", length(ws.paired_reads_datastores), "):")
    for ds in ws.paired_reads_datastores
        println(io, "  ", summary(ds)[23:end])
    end

    println(" Long Read Datastores (", length(ws.long_reads_datastores), "):")
    for ds in ws.long_reads_datastores
        println(io, "  ", summary(ds)[23:end])
    end

    println(" Linked Read Datastores (", length(ws.linked_reads_datastores), "):")
    for ds in ws.linked_reads_datastores
        println(io, "  ", summary(ds)[23:end])
    end

    println(" Mer Count Datastores (", length(ws.mer_count_stores), "):")
    for ds in ws.mer_count_stores
        msg = summary(ds)
        i = first(findfirst("'", msg))
        println(io, "  ", msg[i:end])
    end
end

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
using Spark
using ArgParse

include("op_tests.jl")
include("fault_tests.jl")

function parse_cli()
    s = ArgParseSettings("run_tests.jl")

    @add_arg_table s begin
        "tests"
            help = "Names of test, 'operation', 'fault', or 'all'"
            arg_type = ASCIIString
            default = "all"
    end

    return parse_args(s)
end

function start_worker(port::Int64)
    worker = Spark.Worker("127.0.0.1", port)
    Spark.start(worker)
    println("done")
end

function main()
    master = Spark.Master("127.0.0.1", 3333)
    for i in 0:2
        @async start_worker(6666+i)
    end
    # fill in master.workers
    Spark.load(master, "default_workers.json")
    # start master listener
    Spark.initserver(master)

    function run_test(test::ASCIIString)
        println(test, " test running")
        try
            eval(Expr(:call, symbol(test), master))
            print_with_color(:green, string(test, " test passed\n"))
        catch e
            println(e)
            print_with_color(:red, string(test, " test failed\n"))
        end
    end

    args = parse_cli()
    test = args["tests"]
    op_tests = {"input_test", "collect_test", "count_test", "lookup_test","partition_by_test", "filter_test", "map_test", "group_by_key_test", "join_test"}
    fault_tests = {"basic_disc_fault", "join_disc_fault"}

    tests = {}
    if test == "operation"
        tests = cat(1, tests, op_tests)
    elseif test == "fault"
        tests = cat(1, tests, fault_tests)
    elseif test == "all"
        tests = cat(1, tests, op_tests)
        tests = cat(1, tests, fault_tests)
    else
        tests = cat(1, tests, test)
    end

    for t in tests
        run_test(t)
    end
end

main()


text:
println("==============flag! MOI_wrapper>variables.jl======================================")
"""
MOI variables
"""

"""
Getters for variables
"""
MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variable_info)

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

"""
``MOI.add_variable()`` overloads and safety functions
"""
function MOI.add_variable(model::Optimizer)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end

function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function check_inbounds(model::Optimizer, vi::VI)
	num_variables = length(model.variable_info)
	if !(1 <= vi.value <= num_variables)
	    @error "Invalid variable index $vi. ($num_variables variables in the model.)"
	end
	return
end

check_inbounds(model::Optimizer, var::SVF) = check_inbounds(model, var.variable)

function check_inbounds(model::Optimizer, aff::SAF)
	for term in aff.terms
	    check_inbounds(model, term.variable_index)
	end
	return
end

function check_inbounds(model::Optimizer, quad::SQF)
	for term in quad.affine_terms
	    check_inbounds(model, term.variable_index)
	end
	for term in quad.quadratic_terms
	    check_inbounds(model, term.variable_index_1)
	    check_inbounds(model, term.variable_index_2)
	end
	return
end

has_upper_bound(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].has_upper_bound

has_lower_bound(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].has_lower_bound

is_fixed(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].is_fixed

"""
Primal-start support for Juniper
"""
MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, vi::MOI.VariableIndex, value::Union{Real, Nothing})
    check_inbounds(model, vi)
    if value === nothing
        value = 0.0
    end
    model.variable_info[vi.value].start = value
    return
end

"""
MOI.VariableName attribute support
"""
MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{VI}) = true

function MOI.set(model::Optimizer, ::MOI.VariableName, vi::VI, name::String)
	model.variable_info[vi.value].name = name
	return
end

function MOI.get(model::Optimizer, ::MOI.VariableName, vi::VI)
	return model.variable_info[vi.value].name
end

"""
Populating Variable bounds
"""
function MOI.add_constraint(model::Optimizer, v::SVF, lt::MOI.LessThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(lt.upper)
        @error "Invalid upper bound value $(lt.upper)."
    end
    if has_upper_bound(model, vi)
        @error "Upper bound on variable $vi already exists."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is fixed. Cannot also set upper bound."
    end
    model.variable_info[vi.value].upper_bound = lt.upper
    model.variable_info[vi.value].has_upper_bound = true
    return MOI.ConstraintIndex{SVF, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::SVF, gt::MOI.GreaterThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(gt.lower)
        @error "Invalid lower bound value $(gt.lower)."
    end
    if has_lower_bound(model, vi)
        @error "Lower bound on variable $vi already exists."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is fixed. Cannot also set lower bound."
    end
    model.variable_info[vi.value].lower_bound = gt.lower
    model.variable_info[vi.value].has_lower_bound = true
    return MOI.ConstraintIndex{SVF, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::SVF, eq::MOI.EqualTo{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(eq.value)
        @error "Invalid fixed value $(eq.value)."
    end
    if has_lower_bound(model, vi)
        @error "Variable $vi has a lower bound. Cannot be fixed."
    end
    if has_upper_bound(model, vi)
        @error "Variable $vi has an upper bound. Cannot be fixed."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is already fixed."
    end
    model.variable_info[vi.value].lower_bound = eq.value
    model.variable_info[vi.value].upper_bound = eq.value
    model.variable_info[vi.value].is_fixed = true
    return MOI.ConstraintIndex{SVF, MOI.EqualTo{Float64}}(vi.value)
end


text:
BigInt(1234)

BigInt(typemax(Int64) + 1)    # Incorrect - overflow

BigInt(typemax(Int64)) + 1    # Correct - evaluated using higher precision

parse(BigInt, "12345678901234567890")

factorial(BigInt(30))            # Would have given overflow with Int64

base_exp = [519432 525806; 632382 518061; 78864 613712; 466580 530130; 780495 510032; 525895 525320; 15991 714883; 960290 502358; 760018 511029; 166800 575487; 210884 564478; 555151 523163; 681146 515199; 563395 522587; 738250 512126; 923525 503780; 595148 520429; 177108 572629; 750923 511482; 440902 532446; 881418 505504; 422489 534197; 979858 501616; 685893 514935; 747477 511661; 167214 575367; 234140 559696; 940238 503122; 728969 512609; 232083 560102; 900971 504694; 688801 514772; 189664 569402; 891022 505104; 445689 531996; 119570 591871; 821453 508118; 371084 539600; 911745 504251; 623655 518600; 144361 582486; 352442 541775; 420726 534367; 295298 549387; 6530 787777; 468397 529976; 672336 515696; 431861 533289; 84228 610150; 805376 508857; 444409 532117; 33833 663511; 381850 538396; 402931 536157; 92901 604930; 304825 548004; 731917 512452; 753734 511344; 51894 637373; 151578 580103; 295075 549421; 303590 548183; 333594 544123; 683952 515042; 60090 628880; 951420 502692; 28335 674991; 714940 513349; 343858 542826; 549279 523586; 804571 508887; 260653 554881; 291399 549966; 402342 536213; 408889 535550; 40328 652524; 375856 539061; 768907 510590; 165993 575715; 976327 501755; 898500 504795; 360404 540830; 478714 529095; 694144 514472; 488726 528258; 841380 507226; 328012 544839; 22389 690868; 604053 519852; 329514 544641; 772965 510390; 492798 527927; 30125 670983; 895603 504906; 450785 531539; 840237 507276; 380711 538522; 63577 625673; 76801 615157; 502694 527123; 597706 520257; 310484 547206; 944468 502959; 121283 591152; 451131 531507; 566499 522367; 425373 533918; 40240 652665; 39130 654392; 714926 513355; 469219 529903; 806929 508783; 287970 550487; 92189 605332; 103841 599094; 671839 515725; 452048 531421; 987837 501323; 935192 503321; 88585 607450; 613883 519216; 144551 582413; 647359 517155; 213902 563816; 184120 570789; 258126 555322; 502546 527130; 407655 535678; 401528 536306; 477490 529193; 841085 507237; 732831 512408; 833000 507595; 904694 504542; 581435 521348; 455545 531110; 873558 505829; 94916 603796; 720176 513068; 545034 523891; 246348 557409; 556452 523079; 832015 507634; 173663 573564; 502634 527125; 250732 556611; 569786 522139; 216919 563178; 521815 525623; 92304 605270; 164446 576167; 753413 511364; 11410 740712; 448845 531712; 925072 503725; 564888 522477; 7062 780812; 641155 517535; 738878 512100; 636204 517828; 372540 539436; 443162 532237; 571192 522042; 655350 516680; 299741 548735; 581914 521307; 965471 502156; 513441 526277; 808682 508700; 237589 559034; 543300 524025; 804712 508889; 247511 557192; 543486 524008; 504383 526992; 326529 545039; 792493 509458; 86033 609017; 126554 589005; 579379 521481; 948026 502823; 404777 535969; 265767 554022; 266876 553840; 46631 643714; 492397 527958; 856106 506581; 795757 509305; 748946 511584; 294694 549480; 409781 535463; 775887 510253; 543747 523991; 210592 564536; 517119 525990; 520253 525751; 247926 557124; 592141 520626; 346580 542492; 544969 523902; 506501 526817; 244520 557738; 144745 582349; 69274 620858; 292620 549784; 926027 503687; 736320 512225; 515528 526113; 407549 535688; 848089 506927; 24141 685711; 9224 757964; 980684 501586; 175259 573121; 489160 528216; 878970 505604; 969546 502002; 525207 525365; 690461 514675; 156510 578551; 659778 516426; 468739 529945; 765252 510770; 76703 615230; 165151 575959; 29779 671736; 928865 503569; 577538 521605; 927555 503618; 185377 570477; 974756 501809; 800130 509093; 217016 563153; 365709 540216; 774508 510320; 588716 520851; 631673 518104; 954076 502590; 777828 510161; 990659 501222; 597799 520254; 786905 509727; 512547 526348; 756449 511212; 869787 505988; 653747 516779; 84623 609900; 839698 507295; 30159 670909; 797275 509234; 678136 515373; 897144 504851; 989554 501263; 413292 535106; 55297 633667; 788650 509637; 486748 528417; 150724 580377; 56434 632490; 77207 614869; 588631 520859; 611619 519367; 100006 601055; 528924 525093; 190225 569257; 851155 506789; 682593 515114; 613043 519275; 514673 526183; 877634 505655; 878905 505602; 1926 914951; 613245 519259; 152481 579816; 841774 507203; 71060 619442; 865335 506175; 90244 606469; 302156 548388; 399059 536557; 478465 529113; 558601 522925; 69132 620966; 267663 553700; 988276 501310; 378354 538787; 529909 525014; 161733 576968; 758541 511109; 823425 508024; 149821 580667; 269258 553438; 481152 528891; 120871 591322; 972322 501901; 981350 501567; 676129 515483; 950860 502717; 119000 592114; 392252 537272; 191618 568919; 946699 502874; 289555 550247; 799322 509139; 703886 513942; 194812 568143; 261823 554685; 203052 566221; 217330 563093; 734748 512313; 391759 537328; 807052 508777; 564467 522510; 59186 629748; 113447 594545; 518063 525916; 905944 504492; 613922 519213; 439093 532607; 445946 531981; 230530 560399; 297887 549007; 459029 530797; 403692 536075; 855118 506616; 963127 502245; 841711 507208; 407411 535699; 924729 503735; 914823 504132; 333725 544101; 176345 572832; 912507 504225; 411273 535308; 259774 555036; 632853 518038; 119723 591801; 163902 576321; 22691 689944; 402427 536212; 175769 572988; 837260 507402; 603432 519893; 313679 546767; 538165 524394; 549026 523608; 61083 627945; 898345 504798; 992556 501153; 369999 539727; 32847 665404; 891292 505088; 152715 579732; 824104 507997; 234057 559711; 730507 512532; 960529 502340; 388395 537687; 958170 502437; 57105 631806; 186025 570311; 993043 501133; 576770 521664; 215319 563513; 927342 503628; 521353 525666; 39563 653705; 752516 511408; 110755 595770; 309749 547305; 374379 539224; 919184 503952; 990652 501226; 647780 517135; 187177 570017; 168938 574877; 649558 517023; 278126 552016; 162039 576868; 658512 516499; 498115 527486; 896583 504868; 561170 522740; 747772 511647; 775093 510294; 652081 516882; 724905 512824; 499707 527365; 47388 642755; 646668 517204; 571700 522007; 180430 571747; 710015 513617; 435522 532941; 98137 602041; 759176 511070; 486124 528467; 526942 525236; 878921 505604; 408313 535602; 926980 503640; 882353 505459; 566887 522345; 3326 853312; 911981 504248; 416309 534800; 392991 537199; 622829 518651; 148647 581055; 496483 527624; 666314 516044; 48562 641293; 672618 515684; 443676 532187; 274065 552661; 265386 554079; 347668 542358; 31816 667448; 181575 571446; 961289 502320; 365689 540214; 987950 501317; 932299 503440; 27388 677243; 746701 511701; 492258 527969; 147823 581323; 57918 630985; 838849 507333; 678038 515375; 27852 676130; 850241 506828; 818403 508253; 131717 587014; 850216 506834; 904848 504529; 189758 569380; 392845 537217; 470876 529761; 925353 503711; 285431 550877; 454098 531234; 823910 508003; 318493 546112; 766067 510730; 261277 554775; 421530 534289; 694130 514478; 120439 591498; 213308 563949; 854063 506662; 365255 540263; 165437 575872; 662240 516281; 289970 550181; 847977 506933; 546083 523816; 413252 535113; 975829 501767; 361540 540701; 235522 559435; 224643 561577; 736350 512229; 328303 544808; 35022 661330; 307838 547578; 474366 529458; 873755 505819; 73978 617220; 827387 507845; 670830 515791; 326511 545034; 309909 547285; 400970 536363; 884827 505352; 718307 513175; 28462 674699; 599384 520150; 253565 556111; 284009 551093; 343403 542876; 446557 531921; 992372 501160; 961601 502308; 696629 514342; 919537 503945; 894709 504944; 892201 505051; 358160 541097; 448503 531745; 832156 507636; 920045 503924; 926137 503675; 416754 534757; 254422 555966; 92498 605151; 826833 507873; 660716 516371; 689335 514746; 160045 577467; 814642 508425; 969939 501993; 242856 558047; 76302 615517; 472083 529653; 587101 520964; 99066 601543; 498005 527503; 709800 513624; 708000 513716; 20171 698134; 285020 550936; 266564 553891; 981563 501557; 846502 506991; 334 1190800; 209268 564829; 9844 752610; 996519 501007; 410059 535426; 432931 533188; 848012 506929; 966803 502110; 983434 501486; 160700 577267; 504374 526989; 832061 507640; 392825 537214; 443842 532165; 440352 532492; 745125 511776; 13718 726392; 661753 516312; 70500 619875; 436952 532814; 424724 533973; 21954 692224; 262490 554567; 716622 513264; 907584 504425; 60086 628882; 837123 507412; 971345 501940; 947162 502855; 139920 584021; 68330 621624; 666452 516038; 731446 512481; 953350 502619; 183157 571042; 845400 507045; 651548 516910; 20399 697344; 861779 506331; 629771 518229; 801706 509026; 189207 569512; 737501 512168; 719272 513115; 479285 529045; 136046 585401; 896746 504860; 891735 505067; 684771 514999; 865309 506184; 379066 538702; 503117 527090; 621780 518717; 209518 564775; 677135 515423; 987500 501340; 197049 567613; 329315 544673; 236756 559196; 357092 541226; 520440 525733; 213471 563911; 956852 502490; 702223 514032; 404943 535955; 178880 572152; 689477 514734; 691351 514630; 866669 506128; 370561 539656; 739805 512051; 71060 619441; 624861 518534; 261660 554714; 366137 540160; 166054 575698; 601878 519990; 153445 579501; 279899 551729; 379166 538691; 423209 534125; 675310 515526; 145641 582050; 691353 514627; 917468 504026; 284778 550976; 81040 612235; 161699 576978; 616394 519057; 767490 510661; 156896 578431; 427408 533714; 254849 555884; 737217 512182; 897133 504851; 203815 566051; 270822 553189; 135854 585475; 778805 510111; 784373 509847; 305426 547921; 733418 512375; 732087 512448; 540668 524215; 702898 513996; 628057 518328; 640280 517587; 422405 534204; 10604 746569; 746038 511733; 839808 507293; 457417 530938; 479030 529064; 341758 543090; 620223 518824; 251661 556451; 561790 522696; 497733 527521; 724201 512863; 489217 528217; 415623 534867; 624610 518548; 847541 506953; 432295 533249; 400391 536421; 961158 502319; 139173 584284; 421225 534315; 579083 521501; 74274 617000; 701142 514087; 374465 539219; 217814 562985; 358972 540995; 88629 607424; 288597 550389; 285819 550812; 538400 524385; 809930 508645; 738326 512126; 955461 502535; 163829 576343; 826475 507891; 376488 538987; 102234 599905; 114650 594002; 52815 636341; 434037 533082; 804744 508880; 98385 601905; 856620 506559; 220057 562517; 844734 507078; 150677 580387; 558697 522917; 621751 518719; 207067 565321; 135297 585677; 932968 503404; 604456 519822; 579728 521462; 244138 557813; 706487 513800; 711627 513523; 853833 506674; 497220 527562; 59428 629511; 564845 522486; 623621 518603; 242689 558077; 125091 589591; 363819 540432; 686453 514901; 656813 516594; 489901 528155; 386380 537905; 542819 524052; 243987 557841; 693412 514514; 488484 528271; 896331 504881; 336730 543721; 728298 512647; 604215 519840; 153729 579413; 595687 520398; 540360 524240; 245779 557511; 924873 503730; 509628 526577; 528523 525122; 3509 847707; 522756 525555; 895447 504922; 44840 646067; 45860 644715; 463487 530404; 398164 536654; 894483 504959; 619415 518874; 966306 502129; 990922 501212; 835756 507474; 548881 523618; 453578 531282; 474993 529410; 80085 612879; 737091 512193; 50789 638638; 979768 501620; 792018 509483; 665001 516122; 86552 608694; 462772 530469; 589233 520821; 891694 505072; 592605 520594; 209645 564741; 42531 649269; 554376 523226; 803814 508929; 334157 544042; 175836 572970; 868379 506051; 658166 516520; 278203 551995; 966198 502126; 627162 518387; 296774 549165; 311803 547027; 843797 507118; 702304 514032; 563875 522553; 33103 664910; 191932 568841; 543514 524006; 506835 526794; 868368 506052; 847025 506971; 678623 515342; 876139 505726; 571997 521984; 598632 520198; 213590 563892; 625404 518497; 726508 512738; 689426 514738; 332495 544264; 411366 535302; 242546 558110; 315209 546555; 797544 509219; 93889 604371; 858879 506454; 124906 589666; 449072 531693; 235960 559345; 642403 517454; 720567 513047; 705534 513858; 603692 519870; 488137 528302; 157370 578285; 63515 625730; 666326 516041; 619226 518883; 443613 532186; 597717 520257; 96225 603069; 86940 608450; 40725 651929; 460976 530625; 268875 553508; 270671 553214; 363254 540500; 384248 538137; 762889 510892; 377941 538833; 278878 551890; 176615 572755; 860008 506412; 944392 502967; 608395 519571; 225283 561450; 45095 645728; 333798 544090; 625733 518476; 995584 501037; 506135 526853; 238050 558952; 557943 522972; 530978 524938; 634244 517949; 177168 572616; 85200 609541; 953043 502630; 523661 525484; 999295 500902; 840803 507246; 961490 502312; 471747 529685; 380705 538523; 911180 504275; 334149 544046; 478992 529065; 325789 545133; 335884 543826; 426976 533760; 749007 511582; 667067 516000; 607586 519623; 674054 515599; 188534 569675; 565185 522464; 172090 573988; 87592 608052; 907432 504424; 8912 760841; 928318 503590; 757917 511138; 718693 513153; 315141 546566; 728326 512645; 353492 541647; 638429 517695; 628892 518280; 877286 505672; 620895 518778; 385878 537959; 423311 534113; 633501 517997; 884833 505360; 883402 505416; 999665 500894; 708395 513697; 548142 523667; 756491 511205; 987352 501340; 766520 510705; 591775 520647; 833758 507563; 843890 507108; 925551 503698; 74816 616598; 646942 517187; 354923 541481; 256291 555638; 634470 517942; 930904 503494; 134221 586071; 282663 551304; 986070 501394; 123636 590176; 123678 590164; 481717 528841; 423076 534137; 866246 506145; 93313 604697; 783632 509880; 317066 546304; 502977 527103; 141272 583545; 71708 618938; 617748 518975; 581190 521362; 193824 568382; 682368 515131; 352956 541712; 351375 541905; 505362 526909; 905165 504518; 128645 588188; 267143 553787; 158409 577965; 482776 528754; 628896 518282; 485233 528547; 563606 522574; 111001 595655; 115920 593445; 365510 540237; 959724 502374; 938763 503184; 930044 503520; 970959 501956; 913658 504176; 68117 621790; 989729 501253; 567697 522288; 820427 508163; 54236 634794; 291557 549938; 124961 589646; 403177 536130; 405421 535899; 410233 535417; 815111 508403; 213176 563974; 83099 610879; 998588 500934; 513640 526263; 129817 587733; 1820 921851; 287584 550539; 299160 548820; 860621 506386; 529258 525059; 586297 521017; 953406 502616; 441234 532410; 986217 501386; 781938 509957; 461247 530595; 735424 512277; 146623 581722; 839838 507288; 510667 526494; 935085 503327; 737523 512167; 303455 548204; 992779 501145; 60240 628739; 939095 503174; 794368 509370; 501825 527189; 459028 530798; 884641 505363; 512287 526364; 835165 507499; 307723 547590; 160587 577304; 735043 512300; 493289 527887; 110717 595785; 306480 547772; 318593 546089; 179810 571911; 200531 566799; 314999 546580; 197020 567622; 301465 548487; 237808 559000; 131944 586923; 882527 505449; 468117 530003; 711319 513541; 156240 578628; 965452 502162; 992756 501148; 437959 532715; 739938 512046; 614249 519196; 391496 537356; 62746 626418; 688215 514806; 75501 616091; 883573 505412; 558824 522910; 759371 511061; 173913 573489; 891351 505089; 727464 512693; 164833 576051; 812317 508529; 540320 524243; 698061 514257; 69149 620952; 471673 529694; 159092 577753; 428134 533653; 89997 606608; 711061 513557; 779403 510081; 203327 566155; 798176 509187; 667688 515963; 636120 517833; 137410 584913; 217615 563034; 556887 523038; 667229 515991; 672276 515708; 325361 545187; 172115 573985; 13846 725685]

bases = base_exp[:,1]
exponents = base_exp[:,2]
values = BigInt.(bases) .^ BigInt.(exponents)
max_value, max_idx = findmax(values)
println("Largest value on row ", max_idx, ": ", bases[max_idx], "^", exponents[max_idx])