text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Iterators.Enumerate

"""
    AsyncCollector(f, results, c...; ntasks=0) -> iterator

Apply `f` to each element of `c` using at most `ntasks` asynchronous
tasks.
If `ntasks` is unspecified, uses `max(100, nworkers())` tasks.
For multiple collection arguments, apply `f` elementwise.
Output is collected into `results`.

Note: `next(::AsyncCollector, state) -> (nothing, state)`

Note: `for task in AsyncCollector(f, results, c...) end` is equivalent to
`map!(f, results, c...)`.
"""
type AsyncCollector
    f
    results
    enumerator::Enumerate
    max_tasks::Function
    task_chnl::Channel{Tuple{Int, Any}}             # to communicate with the tasks

    AsyncCollector(f, r, en::Enumerate, mt::Function, c::Channel) = new(f, r, en, mt, c)
end

function AsyncCollector(f, results, c...; ntasks=0)
    if ntasks == 0
        ntasks = max(nworkers(), 100)
        max_tasks = ()->ntasks
    elseif isa(ntasks, Integer)
        max_tasks = ()->ntasks
    elseif isa(ntasks, Function)
        max_tasks = ntasks
    else
        throw(ArgumentError("ntasks must be an Integer or a zero-arg function returning the maximum number of tasks allowed."))
    end
    AsyncCollector(f, results, enumerate(zip(c...)), max_tasks, Channel{Tuple{Int, Any}}(typemax(Int)))
end

type AsyncCollectorState
    enum_state
    active_count::Int
    item_done::Condition
    done::Bool
    in_error::Bool
    nfree::Int                  # number of free tasks
end

isbusy(itr::AsyncCollector, state::AsyncCollectorState) = (state.nfree == 0)

# Wait for @async task to end.
wait(state::AsyncCollectorState) = wait(state.item_done)

function start_collector_task(itr::AsyncCollector, state::AsyncCollectorState)
    t = @async begin
        try
            for (i, args) in itr.task_chnl
                state.nfree -= 1

                itr.results[i] = itr.f(args...)
                notify(state.item_done, nothing)

                state.nfree += 1
            end
        catch e
            # The in_error flag causes done() to end the iteration early and call sync_end().
            # sync_end() then re-throws "e" in the main task.
            state.in_error = true
            clear_collector_channel(itr)
            notify(state.item_done, nothing)

            rethrow(e)
        end
    end

    state.active_count += 1
    t
end

function clear_collector_channel(itr::AsyncCollector)
    try
        # empty out the channel and close it, ignore any errors in doing this
        while isready(itr.task_chnl)
            take!(itr.task_chnl)
        end
        close(itr.task_chnl)
    catch
    end
    nothing
end

# Open a @sync block and initialise iterator state.
function start(itr::AsyncCollector)
    sync_begin()

    state = AsyncCollectorState(start(itr.enumerator),  0, Condition(), false, false, 0)

    for _ in 1:itr.max_tasks()
        start_collector_task(itr, state)
        state.nfree += 1
    end
    state
end

# Close @sync block when iterator is done.
function done(itr::AsyncCollector, state::AsyncCollectorState)
    if state.in_error
        @assert !isopen(itr.task_chnl)  # Channel should have been cleared and closed in the async task.

        sync_end()

        # state.in_error is only being set in the @async block (and an error thrown),
        # which in turn should have been caught and thrown by the sync_end() call above.
        # Control should not come here.
        @assert false "Error should have been captured and thrown previously."
    end

    if !state.done && done(itr.enumerator, state.enum_state)
        state.done = true
        close(itr.task_chnl)
        sync_end()
    end
    return state.done
end

function next(itr::AsyncCollector, state::AsyncCollectorState)
    # start a task if required.
    # Note: Shouldn't need to check on every iteration. Do this at a periodic interval?
    if state.active_count < itr.max_tasks()
        start_collector_task(itr, state)
    end

    while isbusy(itr, state)
        wait(state)
        if state.in_error
            # Stop processing immediately on error.
            return (nothing, state)
        end
    end

    # Get index and mapped function arguments from enumeration iterator.
    (i, args), state.enum_state = next(itr.enumerator, state.enum_state)
    put!(itr.task_chnl, (i, args))

    return (nothing, state)
end

"""
    AsyncGenerator(f, c...; ntasks=0) -> iterator

Apply `f` to each element of `c` using at most `ntasks` asynchronous tasks.
If `ntasks` is unspecified, uses `max(100, nworkers())` tasks.
For multiple collection arguments, apply f elementwise.
Results are returned by the iterator as they become available.
Note: `collect(AsyncGenerator(f, c...; ntasks=1))` is equivalent to
`map(f, c...)`.
"""
type AsyncGenerator
    collector::AsyncCollector
end

function AsyncGenerator(f, c...; ntasks=0)
    AsyncGenerator(AsyncCollector(f, Dict{Int,Any}(), c...; ntasks=ntasks))
end

type AsyncGeneratorState
    i::Int
    async_state::AsyncCollectorState
end

start(itr::AsyncGenerator) = AsyncGeneratorState(0, start(itr.collector))

# Done when source async collector is done and all results have been consumed.
function done(itr::AsyncGenerator, state::AsyncGeneratorState)
    done(itr.collector, state.async_state) && isempty(itr.collector.results)
end

# Pump the source async collector if it is not already busy.
function pump_source(itr::AsyncGenerator, state::AsyncGeneratorState)
    if !isbusy(itr.collector, state.async_state) &&
       !done(itr.collector, state.async_state)
        ignored, state.async_state = next(itr.collector, state.async_state)
        return true
    else
        return false
    end
end

function next(itr::AsyncGenerator, state::AsyncGeneratorState)
    state.i += 1

    results = itr.collector.results
    while !haskey(results, state.i)

        # Wait for results to become available.
        if !pump_source(itr,state) && !haskey(results, state.i)
            wait(state.async_state)
        end
    end
    r = results[state.i]
    delete!(results, state.i)

    return (r, state)
end

# pass-through iterator traits to the iterable
# on which the mapping function is being applied
iteratorsize(itr::AsyncGenerator) = iteratorsize(itr.collector.enumerator)
size(itr::AsyncGenerator) = size(itr.collector.enumerator)
length(itr::AsyncGenerator) = length(itr.collector.enumerator)

"""
    asyncmap(f, c...) -> collection

Transform collection `c` by applying `@async f` to each element.

For multiple collection arguments, apply f elementwise.
"""
asyncmap(f, c...) = collect(AsyncGenerator(f, c...))

"""
    asyncmap!(f, c)

In-place version of `asyncmap()`.
"""
asyncmap!(f, c) = (for x in AsyncCollector(f, c, c) end; c)

"""
    asyncmap!(f, results, c...)

Like `asyncmap()`, but stores output in `results` rather returning a collection.
"""
asyncmap!(f, r, c1, c...) = (for x in AsyncCollector(f, r, c1, c...) end; r)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
perf = ref(String)
timing = ref(Float64)
baseline = ref(Float64)
function test()
	tic()
	y,m,d,h,mi,s = 1982,6,30,23,59,59
	for i = 1:1000000
		b = datetime(y,m,d,h,mi,s) #UTC
	end
	return toq()
end
push!(perf,"Create UTC DateTime")
push!(timing,test())
push!(baseline,0.12)
function test()
	tic()
	y,m,d,h,mi,s = 1982,6,30,23,59,59
	for i = 1:1000000
		b = datetime(y,m,d,h,mi,s,VET)
	end
	return toq()
end
push!(perf,"Create VET DateTime")
push!(timing,test())
push!(baseline,2.0)
function test()
	tic()
	y,m,d,h,mi,s = 1982,6,30,23,59,59
	for i = 1:1000000
		b = datetime(y,m,d,h,mi,s,PST)
	end
	return toq()
end
push!(perf,"Create PST DateTime")
push!(timing,test())
push!(baseline,4.3)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = lastdayofmonth(dt)
	end
	return toq()
end
push!(perf,"lastdayofmonth DateTime")
push!(timing,test())
push!(baseline,0.12)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = dayofyear(dt)
	end
	return toq()
end
push!(perf,"dayofyear DateTime")
push!(timing,test())
push!(baseline,0.11)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = dayofweek(dt)
	end
	return toq()
end
push!(perf,"dayofweek DateTime")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = dayofweekinmonth(dt)
	end
	return toq()
end
push!(perf,"dayofweekinmonth DateTime")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = daysofweekinmonth(dt)
	end
	return toq()
end
push!(perf,"daysofweekinmonth DateTime")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = firstdayofweek(dt)
	end
	return toq()
end
push!(perf,"firstdayofweek DateTime")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = lastdayofweek(dt)
	end
	return toq()
end
push!(perf,"lastdayofweek DateTime")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	dt = datetime(2013,7,8,23,59,59)
	for i = 1:1000000
		t = week(dt)
	end
	return toq()
end
push!(perf,"week DateTime")
push!(timing,test())
push!(baseline,0.07)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = years(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ years(1) DateTime")
push!(timing,test())
push!(baseline,0.43)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = months(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ months(1) DateTime")
push!(timing,test())
push!(baseline,0.42)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = weeks(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ weeks(1) DateTime")
push!(timing,test())
push!(baseline,0.07)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = days(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ days(1) DateTime")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = hours(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ hours(1) DateTime")
push!(timing,test())
push!(baseline,0.07)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = minutes(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ minutes(1) DateTime")
push!(timing,test())
push!(baseline,0.08)
function test()
	tic()
	t = datetime(2013,7,8,23,59,59)
	tt = seconds(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ seconds(1) DateTime")
push!(timing,test())
push!(baseline,0.07)

#Date
function test()
	tic()
	y,m,d = 1982,6,30
	for i = 1:1000000
		b = date(y,m,d)
	end
	return toq()
end
push!(perf,"Create Date")
push!(timing,test())
push!(baseline,0.07)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = dayofyear(dt)
	end
	return toq()
end
push!(perf,"dayofyear Date")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = lastdayofmonth(dt)
	end
	return toq()
end
push!(perf,"lastdayofmonth Date")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = dayofweekinmonth(dt)
	end
	return toq()
end
push!(perf,"dayofweekinmonth Date")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = daysofweekinmonth(dt)
	end
	return toq()
end
push!(perf,"daysofweekinmonth Date")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = firstdayofweek(dt)
	end
	return toq()
end
push!(perf,"firstdayofweek Date")
push!(timing,test())
push!(baseline,0.03)
function test()
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = lastdayofweek(dt)
	end
	return toq()
end
push!(perf,"lastdayofweek Date")
push!(timing,test())
push!(baseline,0.03)
function test() ####
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = dayofweek(dt)
	end
	return toq()
end
push!(perf,"dayofweek Date")
push!(timing,test())
push!(baseline,0.0003)
function test() #####
	tic()
	dt = date(2013,7,8)
	for i = 1:1000000
		t = week(dt)
	end
	return toq()
end
push!(perf,"week Date")
push!(timing,test())
push!(baseline,0.02)
function test()
	tic()
	t = date(2013,7,8)
	tt = years(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ years(1) Date")
push!(timing,test())
push!(baseline,0.15)
function test()
	tic()
	t = date(2013,7,8)
	tt = months(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ months(1) Date")
push!(timing,test())
push!(baseline,0.23)
function test()
	tic()
	t = date(2013,7,8)
	tt = weeks(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ weeks(1) Date")
push!(timing,test())
push!(baseline,0.05)
function test()
	tic()
	t = date(2013,7,8)
	tt = days(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ days(1) Date")
push!(timing,test())
push!(baseline,0.06)
#results
results = [perf timing baseline]
println(results)

ext: jl


text:
"""
Bridges `CP.CountCompare` to `CP.Count`.
"""
struct CountCompare2CountBridge{T} <: MOIBC.AbstractBridge
    con::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{CountCompare2CountBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.CountCompare,
) where {T}
    return MOIBC.bridge_constraint(
        CountCompare2CountBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{CountCompare2CountBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.CountCompare,
) where {T}
    f_scalars = MOIU.scalarize(f)
    con = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                f_scalars[1],
                [
                    f_scalars[1 + i] - f_scalars[1 + s.dimension + i]
                    for i in 1:s.dimension
                ]...
            ]
        ),
        CP.Count(s.dimension, MOI.EqualTo(zero(T)))
    )

    return CountCompare2CountBridge(con)
end

function MOI.supports_constraint(
    ::Type{CountCompare2CountBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.CountCompare},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{CountCompare2CountBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{CountCompare2CountBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}),
    ]
end

function MOI.get(
    ::CountCompare2CountBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}},
    },
) where {T}
    return 1
end

function MOI.get(
    b::CountCompare2CountBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}},
    },
) where {T}
    return [b.con]
end


text:
#!/usr/bin/env julia

# Burma Shave, but evaluating multiple channels per core.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--p", "-p"
    help = "number of processes to add"
    arg_type = Int
    default = 0
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
    # default = 0
    "--chain"
    help = "Write the chain to ~/web/ directory?"
    action = :store_true
    "--opt"
    help = "Use NLOpt instead of LittleMC"
    action = :store_true
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

# Since we've made this a #!/usr/bin/env julia script, we can no longer specify the extra
# processors via the julia prompt, so we need to ghost this behavior back in.

if parsed_args["p"] > 0
    addprocs(parsed_args["p"])
end

import YAML
config = YAML.load(open(parsed_args["config"]))

outfmt(run_index::Int) = config["out_base"] * @sprintf("run%02d/", run_index)

# This code is necessary for multiple simultaneous runs on odyssey
# so that different runs do not write into the same output directory
if parsed_args["run_index"] == nothing
    run_index = 0
    outdir = outfmt(run_index)
    while ispath(outdir)
        println(outdir, " exists")
        run_index += 1
        outdir = outfmt(run_index)
    end
else
    run_index = parsed_args["run_index"]
    outdir = outfmt(run_index)
    println("Deleting old $outdir")
    run(`rm -rf $outdir`)
end

# make the output directory
println("Creating ", outdir)
mkdir(outdir)

using DiskJockey.visibilities
# load data and figure out how many channels
dvarr = DataVis(config["data_file"])
nchan = length(dvarr)

if haskey(config, "exclude")
    exclude = config["exclude"]
    const global keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan]) # which channels of the dset to fit
else
    const global keylist = Int[i for i=1:nchan]
end

# go through any previously created directories and remove them
function cleardirs!(keylist::Vector{Int})
    println("Removing old directories")
    for keys in keylist
        key = keys[1]
        keydir = basedir * "jud$key"
        run(`rm -rf $keydir`)
    end
    println("Removed directories")
end

# This program is meant to be started with the -p option.
nchild = nworkers()
println("Workers allocated ", nchild)

# Chunk up the keylist into an acceptable number of keys per process.
max_keys_per_process = iceil(length(keylist)/nchild)
partial_proc = max_keys_per_process * nchild - length(keylist) # num of unfull processors
full_proc = nchild - partial_proc # number of full processors

chunk_keylist = Array(Any, nchild)
start_key = 1
end_key = start_key + (max_keys_per_process - 1)
for i=1:full_proc
    chunk_keylist[i] = keylist[start_key:end_key]
    start_key += max_keys_per_process
    end_key += max_keys_per_process
end

# Inch back to prep for filling the partially filled processes
end_key -= 1
# Now fill the partial_proc
for i=(full_proc+1):nchild
    chunk_keylist[i] = keylist[start_key:end_key]
    start_key += (max_keys_per_process - 1)
    end_key += (max_keys_per_process - 1)
end

# make the values of run_index and config available on all processes
for process in procs()
    @spawnat process global run_id=run_index
    @spawnat process global cfg=config
end

@everywhere const global species = cfg["species"]

@everywhere basefmt(id::Int) = cfg["base_dir"] * @sprintf("run%02d/", id)

@everywhere const global basedir = basefmt(run_id)

# make the internal Judith directory, if it doesn't exist
if !ispath(basedir)
    println("Creating ", basedir)
    mkdir(basedir)
end

# Clear all directories
cleardirs!(keylist)

@everywhere using DiskJockey.constants
@everywhere using DiskJockey.parallel
@everywhere using DiskJockey.visibilities
@everywhere using DiskJockey.image
@everywhere using DiskJockey.gridding
@everywhere using DiskJockey.model

# Delete the old log file (if it exists)
const logfile = outdir * "log.log"
if isfile(logfile)
    rm(logfile)
end

using Logging
# change the default logger
Logging.configure(filename=logfile, level=DEBUG)

debug("Created logfile.")

@everywhere function initfunc(keys)

    # In contrast to `burma_shave.jl`, now each dataset receives several
    # channels, specified by keys

    # Load the relevant chunk of the dataset, return array
    dvarr = DataVis(cfg["data_file"], keys)
    for dset in dvarr
        # Conjugation is necessary for the SMA, methinks
        visibilities.conj!(dset) # Swap UV convention
    end

    # Create a directory where all RADMC files will reside and be driven from
    # using the first key as the directory index
    key = keys[1]
    keydir = basedir * "jud$key"
    println("Creating $keydir")
    mkdir(keydir)

    # Copy all relevant configuration scripts to this subdirectory
    # these are mainly setup files which will not change throughout the run
    run(`cp radmc3d.inp $keydir`)
    ag = basedir * "amr_grid.inp"
    run(`cp $ag $keydir`)
    run(`cp wavelength_micron.inp $keydir`)

    run(`cp lines.inp $keydir`)
    # run(`cp molecule_co.inp $keydir`)

    # ls = "lines_" * molnames[species] * ".inp"
    # lf = keydir * "/lines.inp"
    # run(`cp $ls $lf`)

    mf = "molecule_" * molnames[species] * ".inp"
    run(`cp $mf $keydir`)

    # change the subprocess to reside in this directory for the remainder of the run
    # where it will drive its own independent RADMC3D process for a subset of channels
    cd(keydir)

    # For each channel, also calculate the interpolation closures
    # npix = cfg["npix"]
    # pixsize = cfg["pixsize"] # [cm]
    # pix_AU = pixsize/AU # [AU]
    #
    # dl = sin(pix_AU/cfg["parameters"]["dpc"][1] * arcsec)
    #
    # uu = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]
    # vv = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]
    #
    # int_arr = Array(Function, length(keys))
    # for (i, dset) in enumerate(dvarr)
    #     int_arr[i] = plan_interpolate(dset, uu, vv)
    # end

    # return the array of datasets
    return dvarr

end

# This is the likelihood function called by each individual process
@everywhere function f(dvarr, keys::Vector{Int}, p::Parameters)

    # dvarr, int_arr = data

    nkeys = length(keys)

    # We are using the following conventions: inclination ranges from
    # 0 to 180 degrees. 0 means face on, angular momentum vector pointing
    # at observer; 90 means edge on; and 180 means face on, angular momentum
    # vector pointing away from observer.
    # These are the same as the RADMC-3D conventions.
    incl = p.incl # [deg]

    # We also adopt the RADMC-3D convention for position angle, which defines position angle
    # by the angular momentum vector.
    # A positive PA angle means the disk angular momentum vector will be
    # rotated counter clockwise (from North towards East).
    PA = p.PA # [deg]

    vel = p.vel # [km/s]
    npix = cfg["npix"] # number of pixels, can alternatively specify x and y separately

    # Doppler shift the dataset wavelengths to rest-frame wavelength
    beta = vel/c_kms # relativistic Doppler formula
    lams = Array(Float64, length(keys))
    for i=1:nkeys
        lams[i] =  dvarr[i].lam * sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    # key = keys[1]
    # keydir = basedir * "jud$key/"
    # println("Writing cameralambda into $keydir")
    write_lambda(lams, "") # write into current directory

    # Run RADMC-3D, redirect output to /dev/null
    run(`radmc3d image incl $incl posang $PA npix $npix loadlambda` |> DevNull)

    # Read the RADMC-3D images from disk (we should already be in sub-directory)
    im = imread()

    # Convert raw images to the appropriate distance
    skim = imToSky(im, p.dpc)

    # Apply the gridding correction function before doing the FFT
    # No shift needed, since we will shift the resampled visibilities
    corrfun!(skim)

    lnprobs = Array(Float64, nkeys)
    # Do the Fourier domain stuff per channel
    for i=1:nkeys
        dv = dvarr[i]
        # FFT the appropriate image channel
        vis_fft = transform(skim, i)

        # Interpolate the `vis_fft` to the same locations as the DataSet
        # mvis = int_arr[i](dv, vis_fft)
        mvis = ModelVis(dv, vis_fft)

        # Apply the phase shift here
        phase_shift!(mvis, p.mu_RA, p.mu_DEC)

        # Calculate chi^2 between these two
        lnprobs[i] = lnprob(dv, mvis)
    end

    # Sum them all together and feed back to the master process
    return sum(lnprobs)

end

# Create the model grid
grd = config["grid"]
global const grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# Regenerate all of the static files (e.g., amr_grid.inp)
# so that they may be later copied
debug("Writing grid")
write_grid(basedir, grid)
debug("Wrote grid")

debug("Initializing processes")
pipes = initialize(nchild, chunk_keylist, initfunc, f)
gather!(pipes)
debug("Initialized processes")

# Calculate the lnprior based upon the current parameter values
function lnprior(pars::Parameters)
    mu_d = 142. # [pc]
    sig_d = 6. # [pc]
    return -0.5 * (pars.dpc - mu_d)^2 / sig_d^2
end

# this function is called only on the main process, which proposes MCMC jumps
# to this function, and farms out the likelihood evaluation to all of the child
# processes
function fprob(p::Vector{Float64})

    # Here is where we make the distinction between a proposed vector of floats
    # (i.e., the parameters), and the object which defines all of the disk parameters
    # which every single subprocess will use

    # Parameters has the following definition (in model.jl)
    # M_star::Float64 # [g] stellar mass
    # r_c::Float64 # [cm] characteristic radius
    # T_10::Float64 # [K] temperature at 10 AU
    # q::Float64 # temperature gradient exponent
    # gamma::Float64 # surface temperature gradient exponent
    # M_gas::Float64 # [g] disk mass of gas
    # ksi::Float64 # [cm s^{-1}] microturbulence
    # dpc::Float64 # [pc] distance to system
    # incl::Float64 # [degrees] inclination 0 deg = face on, 90 = edge on.
    # PA::Float64 # [degrees] position angle (East of North)
    # vel::Float64 # [km/s] systemic velocity (positive is redshift/receeding)
    # mu_RA::Float64 # [arcsec] central offset in RA
    # mu_DEC::Float64 # [arcsec] central offset in DEC

    # Fix the following arguments: gamma, dpc
    gamma = 1.0 # surface temperature gradient exponent

    if config["fix_d"]
        dpc = cfg["parameters"]["dpc"][1] # [pc] distance
        M_star, r_c, T_10, q, logM_gas, ksi, incl, PA, vel, mu_RA, mu_DEC = p
    else
        M_star, r_c, T_10, q, logM_gas, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC = p
    end

    # Enforce hard priors on physical parameters
    # Short circuit evaluation if we know the RADMC won't be valid.
    if ksi <= 0. || T_10 <= 0. || r_c <= 0.0 || M_star <= 0.0
        return -Inf
    end

    if incl < 0. || incl > 180.
        return -Inf
    end

    M_gas = 10^logM_gas

    # If we are going to fit with some parameters dropped out, here's the place to do it
    pars = Parameters(M_star, r_c, T_10, q, gamma, M_gas, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC)

    # Compute parameter file using model.jl, write to disk in base directory
    write_model(pars, basedir, grid, species)

    # nd = basedir * "numberdens_co.inp"
    nd = basedir * "numberdens_" * molnames[species] * ".inp"
    gv = basedir * "gas_velocity.inp"
    gt = basedir * "gas_temperature.inp"
    mt = basedir * "microturbulence.inp"

    # Copy new parameter files to all subdirectories
    for keys in keylist
        key = keys[1]
        keydir = basedir * "jud$key"
        run(`cp $nd $keydir`)
        run(`cp $gv $keydir`)
        run(`cp $gt $keydir`)
        run(`cp $mt $keydir`)
    end

    distribute!(pipes, pars)
    return gather!(pipes) + lnprior(pars)# the summed lnprob
end

# wrapper for NLopt requires gradient as an argument (even if it's not used)
function fgrad(p::Vector, grad::Vector)
    val = fprob(p)
    debug(p, " : ", val)
    return val
end

function fp(p::Vector)
    val = fprob(p)
    debug(p, " : ", val)
    return val
end

debug("Initializing MCMC")
using Distributions
using PDMats

pp = config["parameters"]
# The parameters we'll be using
if config["fix_d"]
    params = ["M_star", "r_c", "T_10", "q", "logM_gas", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
else
    params = ["M_star", "r_c", "T_10", "q", "logM_gas", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
end
# params = ["M_star", "r_c", "T_10", "q", "logM_CO", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)
jumps = Array(Float64, nparam)

for i=1:nparam
    starting_param[i], jumps[i] = pp[params[i]]
end

if parsed_args["opt"]

    # Now try optimizing the function using NLopt
    using NLopt

    nparam = length(starting_param)
    opt = Opt(:LN_COBYLA, nparam)

    max_objective!(opt, fgrad)
    # ftol_abs!(opt, 0.05) # the precision we want lnprob to
    xtol_rel!(opt, 0.001)

    # params = ["M_star", "r_c", "T_10", "q", "logM_CO", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]

    lower = Float64[0.1, 3.0, 1.0, 0.01, -5.0, 0.01, 0.0, 0.0, 42., -0.5, -0.5]
    upper = Float64[10.0, 2000., 2000., 2.0, 1.0, 2.0, 180., 360., 43., 0.5, 0.5]

    lower_bounds!(opt, lower)
    upper_bounds!(opt, upper)

    (optf,optx,ret) = optimize(opt, starting_param)
    println(optf, " ", optx, " ", ret)

else

    # Use LittleMC to do the optimization

    jump_param = PDiagMat(jumps.^2)
    jump_param = full(jump_param)

    # Perturb the starting parameters
    proposal = MvNormal(jump_param)
    starting_param = starting_param .+ 3. * rand(proposal)

    # If we've provided an empirically measured covariance matrix for the MCMC
    # jump proposals, use that instead of our jumps without covariance
    if haskey(config, "opt_jump")
        using NPZ
        jump_param = config["jump_scale"]^2 * npzread(config["opt_jump"])
    end

    using DiskJockey.LittleMC

    # Code to hook into the chain.js plot generator
    if parsed_args["chain"]
        csv = open("/n/home07/iczekala/web/chain.js/mc.csv", "w")
    else
        csv = open(outdir * "mc.csv", "w")
    end

    #Write the parameter header
    writecsv(csv, params')

    mc = MC(fp, config["samples"], starting_param, jump_param, csv)
    debug("Initialized MCMC")

    DiskJockey.LittleMC.start(mc)

    println(mean(mc.samples, 2))
    println(std(mc.samples, 2))

    runstats(mc)

    debug("Acceptance: ", DiskJockey.LittleMC.acceptance(mc))

    DiskJockey.LittleMC.write(mc, outdir * "mc.hdf5")
    close(csv)

end

quit!(pipes)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const SAF{T} = MOI.ScalarAffineFunction{T}
const EQ{T} = MOI.EqualTo{T}
const LT{T} = MOI.LessThan{T}
const GT{T} = MOI.GreaterThan{T}

const SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

const SCALAR_TYPES = Union{
    MOI.ZeroOne,
    MOI.Integer,
}

@enum(
    VariableType,
    CONTINUOUS,
    BINARY,
    INTEGER,
    SEMIINTEGER,
    SEMICONTINUOUS,
)

@enum(
    BoundType,
    NONE,
    LESS_THAN,
    GREATER_THAN,
    LESS_AND_GREATER_THAN,
    INTERVAL,
    EQUAL_TO,
)

mutable struct VariableInfo
    upper::Float64
    lower::Float64
    bound::BoundType
    type::VariableType
    initial_precision::Float64
    # target_precision::Float64
    # start::Union{Float64, Nothing}
    function VariableInfo()
        return new(
            +Inf,
            -Inf,
            NONE,
            CONTINUOUS,
            NaN
            )
    end
end

# Supported Functions
const SF = Union{MOI.VariableIndex,
                 MOI.ScalarAffineFunction{Float64},
                 MOI.VectorOfVariables,
                 MOI.VectorAffineFunction{Float64}}

# Supported Sets
const SS = Union{MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64},
                 MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
                 MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
                 MOI.ExponentialCone, MOI.DualExponentialCone,
                 MOI.PowerCone, MOI.DualPowerCone,
                 MOI.PositiveSemidefiniteConeTriangle}

const QuadraticFunction{T} = Union{
    MOI.ScalarQuadraticFunction{T},
    MOI.VectorQuadraticFunction{T}
}

struct IndexDataCache{T}
    vi   ::Vector{VI}
    sa_eq::Vector{CI{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    sa_lt::Vector{CI{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    sa_gt::Vector{CI{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}
    sv_lt::Vector{CI{MOI.VariableIndex, MOI.LessThan{T}}}
    sv_gt::Vector{CI{MOI.VariableIndex, MOI.GreaterThan{T}}}
    sv_zo::Vector{CI{MOI.VariableIndex, MOI.ZeroOne}}
    function IndexDataCache{T}() where T
        new(
            VI[],
            CI{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[],
            CI{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}[],
            CI{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}[],
            CI{MOI.VariableIndex, MOI.LessThan{T}}[],
            CI{MOI.VariableIndex, MOI.GreaterThan{T}}[],
            CI{MOI.VariableIndex, MOI.ZeroOne}[],
        )
    end
end

mutable struct Optimizer{T, OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    # model to be solved
    optimizer::OT # integer programming solver

    global_initial_precision::Float64
    # global_target_precision::Float64

    # quadratic solver, usually NLP (has to accpet other conic constraints)

    # quadratic constraint cache
    # quadratic_cache::PureQuadraticModel{T}

    # map between quadratic cache and optimizer data
    # goes from ci{quad,set} to ci{aff,set}
    quad_obj::Union{Nothing, MOI.ScalarQuadraticFunction{T}}
    ci_to_quad_scalar::Dict{CI, MOI.ScalarQuadraticFunction{T}}
    ci_to_quad_vector::Dict{CI, MOI.VectorQuadraticFunction{T}}

    original_variables::Dict{VI, VariableInfo}
    ci_to_var::Dict{CI, VI}

    pair_to_var::Dict{Tuple{VI,VI}, VI} # wij variable

    index_cache::Union{Nothing, IndexDataCache{T}}

    has_quad_change::Bool

    fallback_lb::Float64
    fallback_ub::Float64

    allow_soc::Bool

    function Optimizer{T}(
        optimizer::OT;
        lb = -Inf,
        ub = +Inf,
        global_precision = 1e-4
        ) where {T, OT <: MOI.ModelLike}
        # TODO optimizer must support binary, and affine in less and greater
        return new{T, OT}(
            optimizer,
            global_precision,
            nothing,
            Dict{CI, MOI.ScalarQuadraticFunction{T}}(),
            Dict{CI, MOI.VectorQuadraticFunction{T}}(),
            Dict{VI, VariableInfo}(),
            Dict{CI, VI}(),
            Dict{Tuple{VI,VI}, VI}(),
            nothing,
            false,
            lb,
            ub,
            true,
            )
    end
end

function set_quad_change(model)
    model.has_quad_change = true
end
function unset_quad_change(model)
    model.has_quad_change = true
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
    model.quad_obj === nothing &&
    isempty(model.ci_to_quad_scalar) &&
    isempty(model.ci_to_quad_vector) &&
    isempty(model.original_variables) &&
    isempty(model.ci_to_var) &&
    isempty(model.pair_to_var) &&
    model.index_cache === nothing
end

function MOI.empty!(model::Optimizer{T}) where T
    MOI.empty!(model.optimizer)
    model.quad_obj = nothing
    model.ci_to_quad_scalar = Dict{CI, MOI.ScalarQuadraticFunction{T}}()
    model.ci_to_quad_vector = Dict{CI, MOI.VectorQuadraticFunction{T}}()
    model.original_variables = Dict{VI, VariableInfo}()
    model.ci_to_var = Dict{CI, VI}()
    model.pair_to_var = Dict{Tuple{VI,VI}, VI}()
    model.index_cache = nothing
    return
end

function store_quadratic!(model::Optimizer, ci::CI{F,S}, f::F
    ) where {S, F<:MOI.ScalarQuadraticFunction{T}} where T
    model.ci_to_quad_scalar[ci] = f
    return
end
function store_quadratic!(model::Optimizer, ci::CI{F,S}, f::F
    ) where {S, F<:MOI.VectorQuadraticFunction{T}} where T
    model.ci_to_quad_vector[ci] = f
    return
end
function get_quadratic(model::Optimizer, ci::CI{F,S}
    ) where {S, F<:MOI.ScalarQuadraticFunction{T}} where T
    return model.ci_to_quad_scalar[ci]
end
function get_quadratic(model::Optimizer, ci::CI{F,S}
    ) where {S, F<:MOI.VectorQuadraticFunction{T}} where T
    return model.ci_to_quad_vector[ci]
end
function delete_quadratic!(model::Optimizer, ci::CI{F,S}
    ) where {S, F<:MOI.ScalarQuadraticFunction{T}} where T
    delete!(model.ci_to_quad_scalar, ci)
end
function delete_quadratic!(model::Optimizer, ci::CI{F,S}
    ) where {S, F<:MOI.VectorQuadraticFunction{T}} where T
    delete!(model.ci_to_quad_vector, ci)
end
function get_indices(model::Optimizer, ::Type{F}
    ) where {F<:Union{
        MOI.ScalarQuadraticFunction{T},
        MOI.ScalarAffineFunction{T},
        }} where T
    keys(model.ci_to_quad_scalar)
end
function get_indices(model::Optimizer, ::Type{F}
    ) where {F<:Union{
        MOI.VectorQuadraticFunction{T},
        MOI.VectorAffineFunction{T},
        }} where T
    keys(model.ci_to_quad_vector)
end

# MOI.supports(model::Optimizer, args...) = MOI.supports(model.optimizer, args...)

# TODO, call this on inner solver
MOI.supports_incremental_interface(::Optimizer) = true

function MOI.supports(model::Optimizer, attr::MOI.VariableName, tp::Type{MOI.VariableIndex})
    MOI.supports(model.optimizer, attr, tp)
end
function MOI.supports(model::Optimizer, attr::MOI.ConstraintName,
    tb::Type{<:MOI.ConstraintIndex{F,S}}) where {
        F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    return MOI.supports(model.optimizer, attr, tb)
end
function MOI.supports(model::Optimizer, attr::MOI.ConstraintName,
    tb::Type{<:MOI.ConstraintIndex{F,S}}) where {
        F<:QuadraticFunction,S<:MOI.AbstractSet}
    return MOI.supports(model.optimizer, attr, MOI.ConstraintIndex{affine_type(F),S})
end
function MOI.supports(model::Optimizer, attr::Union{
    MOI.Name,
    MOI.Silent,
    MOI.NumberOfThreads,
    MOI.TimeLimitSec,
    MOI.RawOptimizerAttribute,
}
)
    return MOI.supports(model.optimizer, attr)
end
function MOI.get(model::Optimizer, attr::Union{
    MOI.Name,
    MOI.Silent,
    MOI.NumberOfThreads,
    MOI.TimeLimitSec,
    MOI.RawOptimizerAttribute,
}
)
    return MOI.get(model.optimizer, attr)
end
function MOI.set(model::Optimizer, attr::Union{
    MOI.Name,
    MOI.Silent,
    MOI.NumberOfThreads,
    MOI.TimeLimitSec,
}, val
)
    return MOI.set(model.optimizer, attr, val)
end

function MOI.supports(model::Optimizer,
    attr::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{MOI.VariableIndex},
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
            }) where T
    return MOI.supports(model.optimizer, attr)
end
function MOI.supports(model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}) where T
    # this is what reformulation will add to optimizer
    return MOI.supports(model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
end

function MOI.supports_constraint(model::Optimizer,
    F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    return MOI.supports_constraint(model.optimizer, F, S)
end
function MOI.supports_constraint(model::Optimizer,
    F::Type{<:QuadraticFunction{T}}, S::Type{<:MOI.AbstractSet}) where T
    # TODO additional constraints must be tested
    return MOI.supports_constraint(model.optimizer, affine_type(F), S)
end
affine_type(::MOI.ScalarQuadraticFunction{T}) where T = MOI.ScalarAffineFunction{T}
affine_type(::MOI.VectorQuadraticFunction{T}) where T = MOI.VectorAffineFunction{T}
affine_type(::Type{MOI.ScalarQuadraticFunction{T}}) where T = MOI.ScalarAffineFunction{T}
affine_type(::Type{MOI.VectorQuadraticFunction{T}}) where T = MOI.VectorAffineFunction{T}
quadratic_type(::Type{MOI.ScalarAffineFunction{T}}) where T = MOI.ScalarQuadraticFunction{T}
quadratic_type(::Type{MOI.VectorAffineFunction{T}}) where T = MOI.VectorQuadraticFunction{T}

function MOI.get(model::Optimizer, ::MOI.SolverName)
    return "Binary reformulation of quadratic model with solver " *
        MOI.get(model.optimizer, MOI.SolverName()) * " attached"
end

function MOI.supports_add_constrained_variables(
    model::Optimizer, ::Type{S}) where {S<:MOI.AbstractVectorSet}
    return MOI.supports_add_constrained_variables(model.optimizer, S)
end
function MOI.supports_add_constrained_variables(
    model::Optimizer, ::Type{S}) where {S<:MOI.AbstractScalarSet}
    return MOI.supports_add_constrained_variables(model.optimizer, S)
end
function MOI.supports_add_constrained_variables(
    model::Optimizer, ::Type{MOI.Reals})
    return MOI.supports_add_constrained_variables(model.optimizer, MOI.Reals)
end
function MOI.supports_add_constrained_variables(
    model::Optimizer, ::Type{S}) where {
        S<:Union{MOI.SecondOrderCone, MOI.RotatedSecondOrderCone}
    }
    if model.allow_soc
        return MOI.supports_add_constrained_variables(model.optimizer, S)
    else
        return false
    end
end
function MOI.supports_add_constrained_variables(
    model::Optimizer, ::Type{S}) where {
        S<:Union{MOI.SOS1{T}, MOI.SOS2{T}}
    } where T
    return MOI.supports_add_constrained_variables(model.optimizer, S)
end

# function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
#     # if in a subset of the q2b save it
#     # otherwise pass it forward
# end

# function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
# end

function MOI.Utilities.supports_default_copy_to(model::Optimizer, val::Bool)
    return true # val # MOI.Utilities.supports_default_copy_to(model.optimizer, val)
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kwargs...)
    return MOI.Utilities.automatic_copy_to(dest, src; kwargs...)
end

##
## variable
##

function MOI.add_variable(model::Optimizer)
    v = MOI.add_variable(model.optimizer)
    model.original_variables[v] = VariableInfo()
    return v
end

function MOI.add_variables(model::Optimizer, N::Int)
    vs = MOI.add_variables(model.optimizer, N)
    for v in vs
        model.original_variables[v] = VariableInfo()
    end
    return vs
end

function MOI.is_valid(model::Optimizer, v::MOI.VariableIndex)
    #@assert MOI.is_valid(model.quadratic_cache, v)
    return MOI.is_valid(model.optimizer, v)
end

function MOI.delete(model::Optimizer, v::MOI.VariableIndex)
    #MOI.delete(model.quadratic_cache, v)
    # set_quad_change(model)
    # TODO do not delete if in quad
    return MOI.delete(model.optimizer, v)
end

function MOI.get(model::Optimizer, attr::Type{MOI.VariableIndex}, name::String)
    MOI.get(model.optimizer, attr, name)
end
function MOI.get(model::Optimizer, attr::MOI.VariableName, v::MOI.VariableIndex)
    return MOI.get(model.optimizer, attr, v)
end
function MOI.set(
    model::Optimizer, attr::MOI.VariableName, v::MOI.VariableIndex, name::String
)
    return MOI.set(model.optimizer, attr, v, name)
end

##
## objective
##

function MOI.set(
    model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense
)
    return MOI.set(model.optimizer, attr, sense)
end
function MOI.get(model::Optimizer, attr::MOI.ObjectiveSense)
    return MOI.get(model.optimizer, attr)
end
function MOI.set(
    model::Optimizer, attr::MOI.ObjectiveFunction{F}, f::F
) where {F <: Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{T}
}} where T
    model.quad_obj = nothing
    MOI.set(model.optimizer, attr, f)
    return
end
function MOI.set(
    model::Optimizer, attr::MOI.ObjectiveFunction{F}, f::F
) where {F <: MOI.ScalarQuadraticFunction{T}} where T
    fc = MOIU.canonical(f)
    model.quad_obj = fc
    set_quad_change(model)
    f_new = convert_to_affine(model, fc)
    MOI.set(model.optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), f_new)
    return
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}
) where T
    if model.quad_obj === nothing
        f = MOI.get(model.optimizer,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}()
        )
        MOI.ScalarQuadraticFunction{T}(
            MOI.ScalarQuadraticTerm{T}[],
            f.terms,
            f.constant
        )
    else
        return copy(model.quad_obj)
    end
end
function MOI.get(
    model::Optimizer,
    attr::MOI.ObjectiveFunction{
        F
        }
) where F <: Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{T}
} where T
    if model.quad_obj !== nothing
        error()
    else
        return MOI.get(model.optimizer, attr)
    end
end

##
## constraints non-quadratic
##

function MOI.is_valid(
    model::Optimizer,
    ci::MOI.ConstraintIndex{F, S}
) where {F, S}
    return MOI.is_valid(model.optimizer, affine_index(ci))
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F, S}
) where {F, S}
    aci = affine_index(ci)
    if MOI.is_valid(model.optimizer, aci)
        return MOI.get(model.optimizer, attr, aci)
    else
        return throw(MOI.InvalidIndex(ci))
    end
end
function MOI.set(
    model::Optimizer, attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F, S}, s
) where {F, S}
    # TODO improve - only need for quads
    aci = affine_index(ci)
    if MOI.is_valid(model.optimizer, aci)
        return MOI.set(model.optimizer, attr, aci, s)
    else
        return throw(MOI.InvalidIndex(ci))
    end
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F, S}) where {F, S}
    return MOI.get(model.optimizer, attr, ci)
end
function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F, S}, f) where {F, S}
    return MOI.set(model.optimizer, attr, ci, f)
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F, S}) where {F<:QuadraticFunction{T}, S} where T
    # TODO improve
    aci = affine_index(ci)
    if MOI.is_valid(model.optimizer, aci)
        return get_quadratic(model, ci)
    else
        return throw(MOI.InvalidIndex(ci))
    end
end
function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F, S}, f) where {F<:QuadraticFunction{T}, S} where T
    # return MOI.set(model.optimizer, attr, ci, f)
    error("Operation not allowed. Quadratic functions cant be modified.")
end

function MOI.add_constrained_variable(model::Optimizer, s::S
) where {S<:MOI.AbstractScalarSet}
    v, c = MOI.add_constrained_variable(model.optimizer, s)
    model.original_variables[v] = VariableInfo()
    return v, c
end
function MOI.add_constrained_variables(model::Optimizer, s::Vector{S}
) where {S<:MOI.AbstractScalarSet}
    vs, cs = MOI.add_constrained_variables(model.optimizer, s)
    for v in vs
        model.original_variables[v] = VariableInfo()
    end
    return vs, cs
end
function MOI.add_constrained_variables(model::Optimizer, s::S
) where {S<:MOI.AbstractVectorSet}
    vs, cs = MOI.add_constrained_variables(model.optimizer, s)
    for v in vs
        model.original_variables[v] = VariableInfo()
    end
    return vs, cs
end

function MOI.add_constraint(model::Optimizer, f::F, s::S
) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    return MOI.add_constraint(model.optimizer, f, s)
end
function MOI.add_constraints(model::Optimizer, f::Vector{F}, s::Vector{S}
) where {F, S}
    return MOI.add_constraints(model.optimizer, f, s)
end

function MOI.add_constraint(model::Optimizer, f::F, s::S
) where {F<:QuadraticFunction{T}, S<:MOI.AbstractSet} where T
    fc = MOIU.canonical(f)
    ci = MOI.add_constraint(model.optimizer, convert_to_affine(model, fc), s)
    qci = quad_index(ci)
    store_quadratic!(model, qci, fc)
    set_quad_change(model)
    return qci
end
function MOI.add_constraints(model::Optimizer, f::Vector{F}, s::Vector{S}
) where {F<:QuadraticFunction{T}, S<:MOI.AbstractSet} where T
    fc = MOIU.canonical.(f)
    cis = MOI.add_constraints(model.optimizer, convert_to_affine.(model, fc), s)
    qcis = quad_index.(cis)
    for i in eachindex(cis)
        store_quadratic!(model, qcis[i], fc[i])
    end
    set_quad_change(model)
    return qcis
end

function MOI.delete(
    model::Optimizer,
    ci::MOI.ConstraintIndex{F, S}
) where {F<:QuadraticFunction, S<:MOI.AbstractSet}
    delete_quadratic!(model, ci)
    set_quad_change(model)
    aci = affine_index(ci)
    if MOI.is_valid(model.optimizer, aci)
        MOI.delete(model.optimizer, aci)
    else
        return throw(MOI.InvalidIndex(ci))
    end
end

function MOI.delete(
    model::Optimizer,
    ci::MOI.ConstraintIndex{F, S}
) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    MOI.delete(model.optimizer, ci)
end

# function MOI.get(model::Optimizer, attr::MOI.ConstraintName, ci::Type{MOI.ConstraintIndex})
#     return MOI.get(model.optimizer, attr, affine_index(ci))
# end
function MOI.get(model::Optimizer, attr::MOI.ConstraintName, ci::MOI.ConstraintIndex)
    return MOI.get(model.optimizer, attr, affine_index(ci))
end
# function MOI.set(model::Optimizer, attr::MOI.ConstraintName, ci::Type{MOI.ConstraintIndex}, name::String)
#     return MOI.set(model.optimizer, attr, affine_index(ci), name)
# end
function MOI.set(model::Optimizer, attr::MOI.ConstraintName, ci::MOI.ConstraintIndex, name::String)
    return MOI.set(model.optimizer, attr, affine_index(ci), name)
end
function MOI.get(model::Optimizer, attr::Type{<:MOI.ConstraintIndex}, name::String)
    return quad_index(MOI.get(model.optimizer, attr, name), attr)
end

function affine_index(ci::MOI.ConstraintIndex{F, S}) where {F,S}
    return ci
end
function affine_index(ci::MOI.ConstraintIndex{F, S}) where
    {F<:QuadraticFunction{T}, S} where T
    return MOI.ConstraintIndex{affine_type(F), S}(ci.value)
end
function quad_index(ci::MOI.ConstraintIndex{F, S}) where {F,S}
    return CI{quadratic_type(F),S}(ci.value)
end
function quad_index(ci::Nothing, G)
    return ci
end
function quad_index(ci::MOI.ConstraintIndex{F, S}, G) where {F,S}
    return ci
end
function quad_index(ci::MOI.ConstraintIndex{F, S},
    ci2::MOI.ConstraintIndex{F, S}) where {F, S}
    return MOI.ConstraintIndex{F, S}(ci.value)
end
function quad_index(ci::MOI.ConstraintIndex{F, S},
    ci2::MOI.ConstraintIndex{G, S}) where {F, G<:QuadraticFunction{T}, S} where T
    return MOI.ConstraintIndex{G, S}(ci.value)
end
function quad_index(ci::MOI.ConstraintIndex{F, S},
    attr::Type{MOI.ConstraintIndex{F, S}}) where {F, S}
    return MOI.ConstraintIndex{F, S}(ci.value)
end
function quad_index(ci::MOI.ConstraintIndex{F, S},
    attr::Type{MOI.ConstraintIndex{G, S}}) where {F, G<:QuadraticFunction{T}, S} where T
    return MOI.ConstraintIndex{G, S}(ci.value)
end

##
## VariableIndex-in-Set
##

_bounds(s::MOI.GreaterThan{Float64}) = (s.lower, Inf)
_bounds(s::MOI.LessThan{Float64}) = (-Inf, s.upper)
_bounds(s::MOI.EqualTo{Float64}) = (s.value, s.value)
_bounds(s::MOI.Interval{Float64}) = (s.lower, s.upper)

function _throw_if_existing_lower(
    bound::BoundType, var_type::VariableType,
    new_set::Type{<:MOI.AbstractSet},
    variable::MOI.VariableIndex
)
    existing_set = if bound == LESS_AND_GREATER_THAN || bound == GREATER_THAN
        MOI.GreaterThan{Float64}
    elseif bound == INTERVAL
        MOI.Interval{Float64}
    elseif bound == EQUAL_TO
        MOI.EqualTo{Float64}
    elseif var_type == SEMIINTEGER
        MOI.Semiinteger{Float64}
    elseif var_type == SEMICONTINUOUS
        MOI.Semicontinuous{Float64}
    else
        nothing  # Also covers `NONE` and `LESS_THAN`.
    end
    if existing_set !== nothing
        throw(MOI.LowerBoundAlreadySet{existing_set, new_set}(variable))
    end
end
function _throw_if_existing_upper(
    bound::BoundType,
    var_type::VariableType,
    new_set::Type{<:MOI.AbstractSet},
    variable::MOI.VariableIndex
)
    existing_set = if bound == LESS_AND_GREATER_THAN || bound == LESS_THAN
        MOI.LessThan{Float64}
    elseif bound == INTERVAL
        MOI.Interval{Float64}
    elseif bound == EQUAL_TO
        MOI.EqualTo{Float64}
    elseif var_type == SEMIINTEGER
        MOI.Semiinteger{Float64}
    elseif var_type == SEMICONTINUOUS
        MOI.Semicontinuous{Float64}
    else
        nothing  # Also covers `NONE` and `GREATER_THAN`.
    end
    if existing_set !== nothing
        throw(MOI.UpperBoundAlreadySet{existing_set, new_set}(variable))
    end
end
function cache_bounds(
    model::Optimizer, f::MOI.VariableIndex, s::MOI.LessThan{T}) where T
    info = model.original_variables[f]
    _throw_if_existing_upper(info.bound, info.type, MOI.LessThan{T}, f)
    info.bound = info.bound == GREATER_THAN ? LESS_AND_GREATER_THAN : LESS_THAN
    lb, ub = _bounds(s)
    info.upper = ub
    return
end
function cache_bounds(
    model::Optimizer, f::MOI.VariableIndex, s::MOI.GreaterThan{T}) where T
    info = model.original_variables[f]
    _throw_if_existing_lower(info.bound, info.type, MOI.GreaterThan{T}, f)
    info.bound = info.bound == LESS_THAN ? LESS_AND_GREATER_THAN : GREATER_THAN
    lb, ub = _bounds(s)
    info.lower = lb
    return
end
function cache_bounds(
    model::Optimizer, f::MOI.VariableIndex, s::MOI.EqualTo{T}) where T
    info = model.original_variables[f]
    _throw_if_existing_lower(info.bound, info.type, MOI.EqualTo{T}, f)
    _throw_if_existing_upper(info.bound, info.type, MOI.EqualTo{T}, f)
    info.bound = EQUAL_TO
    lb, ub = _bounds(s)
    info.lower = lb
    info.upper = ub
    return
end
function cache_bounds(
    model::Optimizer, f::MOI.VariableIndex, s::MOI.Interval{T}) where T
    info = model.original_variables[f]
    _throw_if_existing_lower(info.bound, info.type, MOI.Interval{T}, f)
    _throw_if_existing_upper(info.bound, info.type, MOI.Interval{T}, f)
    info.bound = INTERVAL
    lb, ub = _bounds(s)
    info.lower = lb
    info.upper = ub
    return
end

function MOI.add_constrained_variables(
    model::Optimizer, s::Vector{S}
) where {S <: SCALAR_SETS}
    vs, c = MOI.add_constrained_variables(model.optimizer, s)
    # from add var
    for v in vs
        model.original_variables[v] = VariableInfo()
    end
    # from add con
    for i in eachindex(vs)
        cache_bounds(model, vs[i], s[i])
        model.ci_to_var[c[i]] = vs[i]
    end
    return vs, c
end
function MOI.add_constrained_variable(
    model::Optimizer, s::S
) where {S <: SCALAR_SETS}
    v, ci = MOI.add_constrained_variable(model.optimizer, s)
    # from add var
    model.original_variables[v] = VariableInfo()
    # from add con
    cache_bounds(model, v, s)
    model.ci_to_var[ci] = v
    return v, ci
end

function MOI.add_constraints(
    model::Optimizer, f::Vector{MOI.VariableIndex}, s::Vector{S}
) where {S <: SCALAR_SETS}
    c = MOI.add_constraints(model.optimizer, f, s)
    for i in eachindex(f)
        cache_bounds(model, f[i], s[i])
        model.ci_to_var[c[i]] = f[i]
    end
    return c
end
function MOI.add_constraint(
    model::Optimizer, f::MOI.VariableIndex, s::S
) where {S <: SCALAR_SETS}
    ci = MOI.add_constraint(model.optimizer, f, s)
    cache_bounds(model, f, s)
    model.ci_to_var[ci] = f
    return ci
end

function MOI.set(
    model::Optimizer, ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.VariableIndex, <:Any}, ::MOI.VariableIndex
)
    return throw(MOI.SettingVariableIndexNotAllowed())
end
function MOI.set(
    model::Optimizer, ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S}, s::S
) where S<:Union{MOI.Interval{T}, MOI.EqualTo{T}} where T
    MOI.throw_if_not_valid(model.optimizer, c)
    MOI.set(model.optimizer, MOI.ConstraintSet(), c, s)
    lower, upper = _bounds(s)
    var = model.ci_to_var[c]
    info = model.original_variables[var]
    info.lower = lower
    info.upper = upper
    return
end
function MOI.set(
    model::Optimizer, ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S}, s::S
) where S<:Union{MOI.LessThan{T}} where T
    MOI.throw_if_not_valid(model.optimizer, c)
    MOI.set(model.optimizer, MOI.ConstraintSet(), c, s)
    lower, upper = _bounds(s)
    var = model.ci_to_var[c]
    info = model.original_variables[var]
    info.upper = upper
    return
end
function MOI.set(
    model::Optimizer, ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S}, s::S
) where S<:Union{MOI.GreaterThan{T}} where T
    MOI.throw_if_not_valid(model.optimizer, c)
    MOI.set(model.optimizer, MOI.ConstraintSet(), c, s)
    lower, upper = _bounds(s)
    var = model.ci_to_var[c]
    info = model.original_variables[var]
    info.lower = lower
    return
end

function MOI.delete(model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}}
) where T
    MOI.throw_if_not_valid(model.optimizer, c)
    var = model.ci_to_var[c]
    MOI.delete(model.optimizer, c)
    delete!(model.ci_to_var, c)
    info = model.original_variables[var]
    info.upper = Inf
    info.bound = info.bound == LESS_AND_GREATER_THAN ? GREATER_THAN : NONE
    return
end
function MOI.delete(model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}}
) where T
    MOI.throw_if_not_valid(model.optimizer, c)
    var = model.ci_to_var[c]
    MOI.delete(model.optimizer, c)
    delete!(model.ci_to_var, c)
    info = model.original_variables[var]
    info.lower = -Inf
    info.bound = info.bound == LESS_AND_GREATER_THAN ? LESS_THAN : NONE
    return
end
function MOI.delete(model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S}
) where S<:Union{MOI.Interval{T}, MOI.EqualTo{T}} where T
    MOI.throw_if_not_valid(model.optimizer, c)
    var = model.ci_to_var[c]
    MOI.delete(model.optimizer, c)
    delete!(model.ci_to_var, c)
    info = model.original_variables[var]
    info.lower = -Inf
    info.upper = +Inf
    info.bound = NONE
    return
end

scalar_type(S::Type{MOI.ZeroOne}) = BINARY
scalar_type(S::Type{MOI.Integer}) = INTEGER

function MOI.add_constrained_variables(
    model::Optimizer, s::Vector{S}
) where {S <: SCALAR_TYPES}
    vs, c = MOI.add_constrained_variables(model.optimizer, s)
    # from add var
    for v in vs
        model.original_variables[v] = VariableInfo()
    end
    # from add con
    for i in eachindex(vs)
        info = model.original_variables[vs[i]]
        info.type = scalar_type(S)
        model.ci_to_var[c[i]] = vs[i]
    end
    return vs, c
end
function MOI.add_constrained_variable(
    model::Optimizer, s::S
) where {S <: SCALAR_TYPES}
    v, ci = MOI.add_constrained_variable(model.optimizer, s)
    # from add var
    model.original_variables[v] = VariableInfo()
    # from add con
    info = model.original_variables[v]
    info.type = scalar_type(S)
    model.ci_to_var[ci] = v
    return v, ci
end

function MOI.add_constraints(
    model::Optimizer, f::Vector{MOI.VariableIndex}, s::Vector{S}
) where {S <: SCALAR_TYPES}
    c = MOI.add_constraints(model.optimizer, f, s)
    for i in eachindex(f)
        info = model.original_variables[f[i]]
        if info.type != BINARY && info.type != INTEGER
            info.type = scalar_type(S)
        else
            error("Variable $(f[i]) is already of type $(info.type)")
        end
    end
    for i in eachindex(c)
        model.ci_to_var[c[i]] = f[i]
    end
    return c
end
function MOI.add_constraint(
    model::Optimizer, f::MOI.VariableIndex, s::S
) where {S <: SCALAR_TYPES}
    ci = MOI.add_constraint(model.optimizer, f, s)
    info = model.original_variables[f]
    if info.type != BINARY && info.type != INTEGER
        info.type = scalar_type(S)
    else
        error("Variable $(f) is already of type $(info.type)")
    end
    model.ci_to_var[ci] = f
    return ci
end
function MOI.delete(model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S}
) where S<:SCALAR_TYPES
    MOI.throw_if_not_valid(model.optimizer, c)
    var = model.ci_to_var[c]
    MOI.delete(model.optimizer, c)
    delete!(model.ci_to_var, c)
    info = model.original_variables[var]
    info.type = CONTINUOUS
    return
end

function add_and_cache_variables(ch, optimizer, N)
    vis = MOI.add_variables(optimizer, N)
    append!(ch.vi, vis)
    return vis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.ScalarAffineFunction{T}}, s::Vector{MOI.EqualTo{T}}) where T
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sa_eq, cis)
    return cis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.ScalarAffineFunction{T}}, s::Vector{MOI.LessThan{T}}) where T
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sa_lt, cis)
    return cis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.ScalarAffineFunction{T}}, s::Vector{MOI.GreaterThan{T}}) where T
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sa_gt, cis)
    return cis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.VariableIndex}, s::Vector{MOI.LessThan{T}}) where T
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sv_lt, cis)
    return cis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.VariableIndex}, s::Vector{MOI.GreaterThan{T}}) where T
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sv_gt, cis)
    return cis
end
function add_and_cache_constraints(ch, optimizer,
    f::Vector{MOI.VariableIndex}, s::Vector{MOI.ZeroOne})
    cis = MOI.add_constraints(optimizer, f, s)
    append!(ch.sv_zo, cis)
    return cis
end

function delete_additional_constraints(model)
    ch = model.index_cache
    if true
        MOI.delete(model.optimizer, ch.sa_eq)
        MOI.delete(model.optimizer, ch.sa_lt)
        MOI.delete(model.optimizer, ch.sa_gt)
        MOI.delete(model.optimizer, ch.vi)
    else
        for ci in ch.sa_eq
            MOI.delete(model.optimizer, ci)
        end
        for ci in ch.sa_lt
            MOI.delete(model.optimizer, ci)
        end
        for ci in ch.sa_gt
            MOI.delete(model.optimizer, ci)
        end
        for vi in ch.vi
            MOI.delete(model.optimizer, vi)
        end
    end
end

function lower(model, info)
    if info.type == BINARY
        return zero(typeof(info.lower))
    elseif info.lower > -Inf
        return info.lower
    end
    return model.fallback_lb
end
function upper(model, info)
    if info.type == BINARY
        return one(typeof(info.lower))
    elseif info.upper < Inf
        return info.upper
    end
    return model.fallback_ub
end
function get_precision(model::Optimizer{T}, x) where T
    info = model.original_variables[x]
    a = model.global_initial_precision
    b = info.initial_precision
    pre = min(isnan(a) ? one(T) : a, isnan(b) ? one(T) : b)
    if info.type == BINARY
        pre = one(T)
    end
    @assert zero(T) < pre <= one(T)
    ceil(Int, -log2(pre))
end
function build_approximation!(model::Optimizer)

    T = Float64

    if model.index_cache !== nothing
        delete_additional_constraints(model)
    end
    ch = IndexDataCache{T}()
    model.index_cache = ch

    QT = keys(model.pair_to_var)

    # variables to expand
    DS = Set()
    # build graph
    neighbors = Dict()
    degree = Dict()
    for (i,j) in keys(model.pair_to_var)
        if i == j
            push!(DS, i)
        else
            for (a,b) in [(i,j), (j,i)]
                if haskey(neighbors, a)
                    push!(neighbors[a], b)
                    degree[a] += 1
                else
                    neighbors[a] = [b]
                    degree[a] = 1
                end
            end
        end
    end

    # Based in: https://github.com/JuliaGraphs/LightGraphs.jl/blob/v1.3.1/src/vertexcover/degree_vertex_cover.jl
    degree_queue = PriorityQueue(Base.Order.Reverse, degree)
    # Loops in the graph must be added to this cover
    for v in DS
        degree_queue[v] = 0
        if haskey(neighbors, v)
            @inbounds @simd for u in neighbors[v]
                if !(v in DS)#!in_cover[u]
                    degree_queue[u] -= 1
                end
            end
        end
    end
    # length_cover = length(DS)
    while !isempty(degree_queue) && peek(degree_queue)[2] > 0
        v = dequeue!(degree_queue)
        # in_cover[v] = true
        push!(DS, v)
        # length_cover += 1
        @inbounds @simd for u in neighbors[v]
            if !(v in DS)#!in_cover[u]
                degree_queue[u] -= 1
            end
        end
    end

    # Bounds detection
    for pair in QT
        for x in pair
            info = model.original_variables[x]
            if !(upper(model, info) < Inf)
                error("Variable $x has no upper bound.")
            end
            if !(lower(model, info) > -Inf)
                error("Variable $x has no lower bound.")
            end
        end
    end

    #
    preΔx = add_and_cache_variables(ch, model.optimizer, length(DS))
    Δx = Dict()
    for (i,v) in enumerate(DS)
        Δx[v] = preΔx[i]
    end
    #
    w = model.pair_to_var
    #
    preΔw = add_and_cache_variables(ch, model.optimizer, length(w))
    Δw = Dict()
    for (i,p) in enumerate(QT)
        Δw[p] = preΔw[i]
    end
    #
    z = Dict()
    for v in DS
        z[v] = add_and_cache_variables(ch, model.optimizer, get_precision(model, v))#VI[]
    end
    #
    xh = Dict()
    for (xa,xb) in QT
        if xb in DS
            xi = xa
            xj = xb # in DS
        else
            xi = xb
            xj = xa
        end
        xh[(xa,xb)] = add_and_cache_variables(ch, model.optimizer, get_precision(model, xj))#VI[]
    end

    # zeros are useless in the paper

    # eq 28 - binary expansion of xj
    f28 = SAF{T}[]
    s28 = EQ{T}[]
    for xj in DS
        info = model.original_variables[xj]
        Xu = upper(model, info)
        Xl = lower(model, info)
        Δxj = Δx[xj]
        zj = z[xj]
        terms = [
            MOI.ScalarAffineTerm(-T(1), xj),
            MOI.ScalarAffineTerm(Xu-Xl, Δxj)
        ]
        for l in 1:get_precision(model, xj)
            push!(terms, MOI.ScalarAffineTerm( (Xu-Xl)*T(2)^(-T(l)), zj[l]))
        end
        push!(f28, MOI.ScalarAffineFunction(terms, zero(T)))
        push!(s28, MOI.EqualTo{T}(-Xl))
    end
    c28 = add_and_cache_constraints(ch, model.optimizer, f28, s28)

    # eq 29 - binary expansion of the product
    f29 = SAF{T}[]
    for (xa,xb) in QT
        if xb in DS
            xi = xa
            xj = xb # in DS
        else
            xi = xb
            xj = xa
        end
        info = model.original_variables[xj]
        Xu = upper(model, info)
        Xl = lower(model, info)
        #order is important in the pairs
        Δwij = Δw[(xa,xb)]
        wij = w[(xa,xb)]
        xhij = xh[(xa,xb)] # vector
        terms = [
            MOI.ScalarAffineTerm(-one(T), wij),
            MOI.ScalarAffineTerm(Xu-Xl, Δwij),
            MOI.ScalarAffineTerm(Xl, xi) #missing in the paper
        ]
        for l in 1:get_precision(model, xj)
            push!(terms, MOI.ScalarAffineTerm( (Xu-Xl)*T(2)^(-T(l)), xhij[l]))
        end
        push!(f29, MOI.ScalarAffineFunction(terms, zero(T)))
    end
    s29 = MOI.EqualTo{T}[MOI.EqualTo{T}(zero(T)) for i in eachindex(f29)]
    c29 = add_and_cache_constraints(ch, model.optimizer, f29, s29)

    # eq 30 - bound on \Delta xj
    f30a = MOI.VariableIndex[]
    f30b = MOI.VariableIndex[]
    s30b = LT{T}[]
    for xj in DS
        push!(f30a, Δx[xj])
        push!(f30b, Δx[xj])
        push!(s30b, MOI.LessThan{T}(T(2)^(-T(get_precision(model, xj)))))
    end
    s30a = MOI.GreaterThan{T}[MOI.GreaterThan{T}(zero(T)) for i in eachindex(f30a)]
    c30a = add_and_cache_constraints(ch, model.optimizer, f30a, s30a)
    c30b = add_and_cache_constraints(ch, model.optimizer, f30b, s30b)

    # eq 31 - McCormick in \Delta wij
    f31a = SAF{T}[]
    f31b = SAF{T}[]
    s31a = GT{T}[]
    s31b = LT{T}[]
    # 32
    f32a = SAF{T}[]
    f32b = SAF{T}[]
    s32a = GT{T}[]
    s32b = LT{T}[]
    for (xa,xb) in QT
        if xb in DS
            xi = xa
            xj = xb # in DS
        else
            xi = xb
            xj = xa
        end
        info_i = model.original_variables[xi]
        Xu_i = upper(model, info_i)
        Xl_i = lower(model, info_i)
        Δwij = Δw[(xa,xb)]
        Δxj = Δx[xj]
        p = get_precision(model, xj)
        # 31
        terms_31a = [
            MOI.ScalarAffineTerm(one(T), Δwij),
            MOI.ScalarAffineTerm(-Xu_i, Δxj), # review here if both are indeed I
            MOI.ScalarAffineTerm(-T(2)^(-T(p)), xi), # review here if both are indeed I
        ]
        push!(f31a, MOI.ScalarAffineFunction(terms_31a, zero(T)))
        push!(s31a, MOI.GreaterThan{T}(-Xu_i*T(2)^(-T(p))))
        terms_31b = [
            MOI.ScalarAffineTerm(one(T), Δwij),
            MOI.ScalarAffineTerm(-Xl_i, Δxj), # review here if both are indeed I
            MOI.ScalarAffineTerm(-T(2)^(-T(p)), xi), # review here if both are indeed I
        ]
        push!(f31b, MOI.ScalarAffineFunction(terms_31b, zero(T)))
        push!(s31b, MOI.LessThan{T}(-Xl_i*T(2)^(-T(p))))
        # 32
        terms_32a = [
            MOI.ScalarAffineTerm(one(T), Δwij),
            MOI.ScalarAffineTerm(-Xl_i, Δxj), # review here if both are indeed I
        ]
        push!(f32a, MOI.ScalarAffineFunction(terms_32a, zero(T)))
        push!(s32a, MOI.GreaterThan{T}(zero(T)))
        terms_32b = [
            MOI.ScalarAffineTerm(one(T), Δwij),
            MOI.ScalarAffineTerm(-Xu_i, Δxj), # review here if both are indeed I
        ]
        push!(f32b, MOI.ScalarAffineFunction(terms_32b, zero(T)))
        push!(s32b, MOI.LessThan{T}(zero(T)))
    end
    c31a = add_and_cache_constraints(ch, model.optimizer, f31a, s31a)
    c31b = add_and_cache_constraints(ch, model.optimizer, f31b, s31b)
    c32a = add_and_cache_constraints(ch, model.optimizer, f32a, s32a)
    c32b = add_and_cache_constraints(ch, model.optimizer, f32b, s32b)

    # eq 33 - Disjunction
    f33a = SAF{T}[]
    f33b = SAF{T}[]
    s33a = GT{T}[]
    s33b = LT{T}[]
    # 34
    f34a = SAF{T}[]
    f34b = SAF{T}[]
    s34a = GT{T}[]
    s34b = LT{T}[]
    for (xa,xb) in QT
        if xb in DS
            xi = xa
            xj = xb # in DS
        else
            xi = xb
            xj = xa
        end
        info_i = model.original_variables[xi]
        Xu_i = upper(model, info_i)
        Xl_i = lower(model, info_i)
        zj = z[xj]
        xhij = xh[(xa,xb)]
        # 33
        for l in 1:get_precision(model, xj)
            terms_33a = [
                MOI.ScalarAffineTerm(one(T), xhij[l]),
                MOI.ScalarAffineTerm(-Xl_i, zj[l]), # review here if both are indeed I
            ]
            push!(f33a, MOI.ScalarAffineFunction(terms_33a, zero(T)))
            push!(s33a, MOI.GreaterThan{T}(zero(T)))
            terms_33b = [
                MOI.ScalarAffineTerm(one(T), xhij[l]),
                MOI.ScalarAffineTerm(-Xu_i, zj[l]), # review here if both are indeed I
            ]
            push!(f33b, MOI.ScalarAffineFunction(terms_33b, zero(T)))
            push!(s33b, MOI.LessThan{T}(zero(T)))
            # 34
            terms_34a = [
                MOI.ScalarAffineTerm( one(T), xi),
                MOI.ScalarAffineTerm(-one(T), xhij[l]),
                MOI.ScalarAffineTerm(Xl_i, zj[l]), # review here if both are indeed I
            ]
            push!(f34a, MOI.ScalarAffineFunction(terms_34a, zero(T)))
            push!(s34a, MOI.GreaterThan{T}(Xl_i))
            terms_34b = [
                MOI.ScalarAffineTerm( one(T), xi),
                MOI.ScalarAffineTerm(-one(T), xhij[l]),
                MOI.ScalarAffineTerm(Xu_i, zj[l]), # review here if both are indeed I
            ]
            push!(f34b, MOI.ScalarAffineFunction(terms_34b, zero(T)))
            push!(s34b, MOI.LessThan{T}(Xu_i))
        end
    end
    c33a = add_and_cache_constraints(ch, model.optimizer, f33a, s33a)
    c33b = add_and_cache_constraints(ch, model.optimizer, f33b, s33b)
    c34a = add_and_cache_constraints(ch, model.optimizer, f34a, s34a)
    c34b = add_and_cache_constraints(ch, model.optimizer, f34b, s34b)

    # 37 - z is binary
    f37 = MOI.VariableIndex[]
    s37 = MOI.ZeroOne[]
    for zj in values(z), zjl in zj
        push!(f37, zjl)
        push!(s37, MOI.ZeroOne())
    end
    c37 = add_and_cache_constraints(ch, model.optimizer, f37, s37)

    unset_quad_change(model)

    # TODO: deal with integers
    # TODO: use model to hold intercepted ctrs?
    return
end

function MOI.optimize!(model::Optimizer)
    if model.has_quad_change || model.index_cache === nothing
        build_approximation!(model)
    end
    return MOI.optimize!(model.optimizer)
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::CI{F,S}) where {F, S}
    error("No duals available")
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::CI{F,S}) where {F, S}
    return MOI.get(model.optimizer, attr, ci)
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
    ci::CI{F,S}) where {F <: QuadraticFunction{T}, S} where T
    return MOI.get(model.optimizer, attr, affine_index(ci))
end

function MOI.get(model::Optimizer, attr::T) where {
    T <: Union{
        MOI.TerminationStatus,
        MOI.ObjectiveValue,
        MOI.PrimalStatus,
    }
}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::Union{
    MOI.AbstractModelAttribute,
    MOI.AbstractOptimizerAttribute
    })
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.RawStatusString)
    return "Raw status from solver: " *
        MOI.get(model.optimizer, MOI.SolverName()) * ", is: " *
        MOI.get(model.optimizer, MOI.RawStatusString())
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ListOfConstraintIndices{F, S}
) where {S, F<:Union{
    MOI.VectorOfVariables,
    MOI.VariableIndex,
}}
    return MOI.get(model.optimizer, attr)
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F, S}
) where {S<:MOI.AbstractSet, F<:Union{
    MOI.ScalarAffineFunction{T},
    MOI.VectorAffineFunction{T},
}} where T
    in_opt = MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F, S}())
    quads = get_indices(model, F)
    return setdiff!(in_opt, affine_index.(quads))
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F, S}
) where {S<:MOI.AbstractSet, F<:Union{
    MOI.ScalarQuadraticFunction{Float64},
    MOI.VectorQuadraticFunction{Float64},
}}
    in_opt = MOI.get(model.optimizer, MOI.ListOfConstraintIndices{affine_type(F), S}())
    quads = get_indices(model, F)
    return intersect!(quad_index.(in_opt), quads)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveFunctionType)
    if model.quad_obj === nothing
        return MOI.get(model.optimizer, attr)
    else
        return typeof(model.quad_obj)
    end
end

function MOI.get(model::Optimizer, attr::MOI.NumberOfVariables)
    return MOI.get(model.optimizer, attr)
end
function MOI.get(model::Optimizer, attr::MOI.ListOfVariableIndices)
    return MOI.get(model.optimizer, attr)
end

MOI.get(model::Optimizer, attr::MOI.RawSolver) = MOI.get(model.optimizer, attr)

# we will need to propagate the primal start
function MOI.set(
    model::Optimizer,
    attr::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    value::Union{Nothing, Float64}
)
    MOI.set(model.optimizer, attr, x, value)
    return
end

function MOI.get(
    model::Optimizer, attr::MOI.VariablePrimalStart, x::MOI.VariableIndex
)
    return MOI.get(model.optimizer, attr, x)
end

struct VariablePrecision <: MOI.AbstractVariableAttribute end
struct GlobalVariablePrecision <: MOI.AbstractOptimizerAttribute end
struct FallbackUpperBound <: MOI.AbstractOptimizerAttribute  end
struct FallbackLowerBound <: MOI.AbstractOptimizerAttribute  end

function MOI.set(
    model::Optimizer,
    attr::VariablePrecision,
    x::MOI.VariableIndex,
    value::Union{Nothing, Float64}
)
    val = value === nothing ? NaN : value
    @assert 0 < val <= 1
    model.original_variables[x].initial_precision = val
    return
end
function MOI.get(
    model::Optimizer, attr::VariablePrecision, x::MOI.VariableIndex
)
    return model.original_variables[x].initial_precision
end

function MOI.set(
    model::Optimizer,
    ::GlobalVariablePrecision,
    value::Union{Nothing, Float64}
)
    val = value === nothing ? NaN : value
    @assert 0 < val <= 1
    model.global_initial_precision = val
    return
end
function MOI.get(
    model::Optimizer, ::GlobalVariablePrecision, x::MOI.VariableIndex
)
    return model.global_initial_precision
end

function MOI.set(model::Optimizer, ::FallbackLowerBound, val)
    model.fallback_lb = val
    return nothing
end
function MOI.set(model::Optimizer, ::FallbackLowerBound, ::Nothing)
    model.fallback_lb = -Inf
    return nothing
end
function MOI.get(model::Optimizer, ::FallbackLowerBound)
    return model.fallback_lb
end
function MOI.set(model::Optimizer, ::FallbackUpperBound, val)
    model.fallback_ub = val
    return nothing
end
function MOI.set(model::Optimizer, ::FallbackUpperBound, ::Nothing)
    model.fallback_ub = Inf
    return nothing
end
function MOI.get(model::Optimizer, ::FallbackUpperBound)
    return model.fallback_ub
end

function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{F, S}) where {F, S}
    # TODO: this could be more efficient.
    return length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}()))
end

function MOI.get(model::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    constraints = Set{Tuple{DataType, DataType}}()
    inner_ctrs = MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())
    for (F, S) in inner_ctrs
        if F <: Union{MOI.ScalarAffineFunction, MOI.VectorQuadraticFunction}
            if MOI.get(model, MOI.NumberOfConstraints{F, S}()) > 0
                push!(constraints, (F,S))
            end
            if MOI.get(model, MOI.NumberOfConstraints{quadratic_type(F), S}()) > 0
                push!(constraints, (quadratic_type(F),S))
            end
        else
            push!(constraints, (F,S))
        end
    end
    # error("TODO")
    # deal with deleted parameters
    collect(constraints)
end

function ordered_term_indices(t::MOI.ScalarQuadraticTerm)
    return VI.(minmax(t.variable_1.value, t.variable_2.value))
end
function ordered_term_indices(t::Union{MOI.VectorAffineTerm, MOI.VectorQuadraticTerm})
    return (t.output_index, ordered_term_indices(t.scalar_term)...)
end
function get_affine_part(f::MOI.ScalarQuadraticFunction{T}) where T
    return MOI.ScalarAffineFunction{T}(f.affine_terms, f.constant)
end
function get_affine_part(f::MOI.VectorQuadraticFunction{T}) where T
    return MOI.VectorAffineFunction{T}(f.affine_terms, f.constants)
end
function convert_to_affine(model::Optimizer, f::QuadraticFunction{T}) where T
    fc = copy(f)
    for t in fc.quadratic_terms
        aff_term = build_affine_term(model, t)
        push!(fc.affine_terms, aff_term)
    end
    return get_affine_part(fc)
end
function build_affine_term(model::Optimizer{T}, t::MOI.ScalarQuadraticTerm) where T
    pair = ordered_term_indices(t)
    var = add_or_get_variable(model, pair)
    m = pair[1] == pair[2] ? T(1/2) : T(1)
    return MOI.ScalarAffineTerm(m*MOI.coefficient(t), var)
end
function build_affine_term(model::Optimizer{T}, t::MOI.VectorQuadraticTerm) where T
    triple = ordered_term_indices(t)
    var = add_or_get_variable(model, (triple[2], triple[3]))
    m = triple[2] == triple[3] ? T(1/2) : T(1)
    return MOI.VectorAffineTerm(triple[1],
        MOI.ScalarAffineTerm(m*MOI.coefficient(t), var))
end
function add_or_get_variable(model::Optimizer, pair)
    # for each unique pair we need a single variable
    if haskey(model.pair_to_var, pair)
        return model.pair_to_var[pair]
    else
        var = MOI.add_variable(model.optimizer)
        model.pair_to_var[pair] = var
        return var
    end
end

##
## Modifications
##

function MOI.modify(
    model::Optimizer,
    mod::Union{
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any},
    },
    chg::Union{
        MOI.ScalarConstantChange{Float64},
        MOI.ScalarCoefficientChange{Float64},
    }
)
    MOI.modify(model.optimizer, mod, chg)
    return
end
function MOI.modify(
    model::Optimizer,
    mod::Union{
        MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, <:Any},
    },
    chg::Union{
        MOI.VectorConstantChange{Float64},
    }
)
    MOI.modify(model.optimizer, mod, chg)
    return
end


text:
using Test, Julianda

foo = Julianda.Config.loadConfig("../config")

bar = Julianda.Position.listPositions(foo)
@test length(bar) != 0

bar = Julianda.Position.listOpenPositions(foo)
@test bar != nothing

bar = Julianda.Position.getPosition(foo, "GBP_USD")
@test bar.instrument == "GBP_USD"

#@test Julianda.Position.closePosition(foo, "GBP_USD", 50)

#@test Julianda.Position.closePositionFull(foo, "GBP_USD")

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
const CKSP = Ptr{Cvoid}

struct PetscKSP
    ptr::Ref{CKSP}

    PetscKSP() = new(Ref{CKSP}())
end

# allows us to pass PetscKSP objects directly into CKSP ccall signatures
Base.cconvert(::Type{CKSP}, ksp::PetscKSP) = ksp.ptr[]

"""
    KSPCreate(comm::MPI.Comm, ksp::PetscKSP)

Wrapper for KSPCreate
"""
function KSPCreate(comm::MPI.Comm, ksp::PetscKSP)
    error = ccall((:KSPCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CKSP}), comm, ksp.ptr)
    @assert iszero(error)
end

function KSPCreate(comm::MPI.Comm = MPI.COMM_WORLD)
    ksp = PetscKSP()
    KSPCreate(comm, ksp)
    return ksp
end

"""
    KSPSetOperators(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat)

Wrapper for KSPSetOperators
"""
function KSPSetOperators(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat)
    error = ccall((:KSPSetOperators, libpetsc), PetscErrorCode, (CKSP, CMat, CMat), ksp, Amat, Pmat)
    @assert iszero(error)
end

"""
    KSPSolve(ksp::PetscKSP, b::PetscVec, x::PetscVec)

Wrapper for KSPSolve
"""
function KSPSolve(ksp::PetscKSP, b::PetscVec, x::PetscVec)
    error = ccall((:KSPSolve, libpetsc), PetscErrorCode, (CKSP, CVec, CVec), ksp, b, x)
    @assert iszero(error)
end

"""
    KSPSetUp(ksp::PetscKSP)

Wrapper to KSPSetUp
"""
function KSPSetUp(ksp::PetscKSP)
    error = ccall((:KSPSetUp, libpetsc), PetscErrorCode, (CKSP,), ksp)
    @assert iszero(error)
end

"""
    KSPSetFromOptions(ksp::PetscKSP)

Wrapper to KSPSetFromOptions
"""
function KSPSetFromOptions(ksp::PetscKSP)
    error = ccall((:KSPSetFromOptions, libpetsc), PetscErrorCode, (CKSP,), ksp)
    @assert iszero(error)
end

"""
    KSPDestroy(ksp::PetscKSP)

Wrapper to `KSPDestroy`
"""
function KSPDestroy(ksp::PetscKSP)
    error = ccall((:KSPDestroy, libpetsc), PetscErrorCode, (Ptr{CKSP},), ksp.ptr)
    @assert iszero(error)
end


text:
"""
    readhdf5(file::String)

Read Nastran HDF5 formated output results.
"""
function readhdf5(file::String)
    println("Did not read nastran h5 file. Still not implemented.")
end

max_stars_repo_name: mgcth/SDNastran

max_issues_repo_name: mgcth/SDNastran

max_forks_repo_name: mgcth/SDNastran