text:
using MAT # include only on main worker
using Test
using DivSigGrad
using jInv.InverseSolve
using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.Utils

using EikonalInv

plotting = false;
if plotting
	using jInvVis
	using PyPlot
end

@everywhere begin
	include("DriversDC/velToConductMod.jl");
end

include("DriversEikonalInv/readModelAndGenerateMeshMref.jl");
include("DriversEikonalInv/prepareTravelTimeDataFiles.jl");
include("DriversEikonalInv/setupTravelTimeTomography.jl");
include("DriversDC/getSourcesAndReceivers.jl")
include("DriversDC/setupDCAndGetData.jl");

invertDC  = true;
invertEik = true;
invertJoint = invertDC & invertEik

matfile   = matread("../3Dseg12812864.mat")
mtrue     = matfile["VELc"]
n         = [64;64;32]
nfwd      = [32;32;16]
domain    = [0;13.5;0;13.5;0;4.2]
mtrue*=1e-3; ## mtrue is in km per sec.

###########################################################################################################################
################ Parameters For Inversion #################################################################################
###########################################################################################################################
(mtrue,MinvNodal,mref,~,~)  = readModelAndGenerateMeshMref(mtrue,0,n,domain);
Minv    					= getRegularMesh(MinvNodal.domain,collect(size(mtrue)));
HesPrec     = getSSORRegularizationPreconditioner(1.0,1e-15,200);
alpha	   	= 1e-15;
cgit       	= 8;
pcgTol     	= 1e-1;
maxIter    	= 10;
minUpdate 	= 1e-2;
Iact 		= sparse(I,prod(n),prod(n));
boundsLow  	= (minimum(mtrue)-0.05)*ones(size(Iact,2)); # minimum(mtrue) is sea...
boundsHigh 	= (maximum(mtrue)+0.05)*ones(size(Iact,2));
maxStep		= 0.1*maximum(boundsHigh);

mtrue   	= Iact'*mtrue[:];
mref    	= Iact'*mref[:];
mback       = mref - Iact*(Iact'*mref)

regparams = [1.0,1.0,1.0,1e-6];
regfun(m,mref,Minv) = wdiffusionReg(m,mref,Minv,Iact=Iact,C = regparams);

if invertJoint
	filenamePrefix = "SEG_Joint_"
elseif invertEik
	filenamePrefix = "SEG_Eik_"
elseif invertDC
	filenamePrefix = "SEG_DC_"
end

function dumpResults(mc,Dc,iter,pInv,pMis)
	mDc = reshape(Iact*vec(mc)+mback,tuple(Minv.n...))
	if plotting
		figure(888,figsize = (22,10));
		clf();
	end
	file = string(filenamePrefix,Minv.n,iter,".dat");
	if plotting
		plotModel(mDc,false,[],0,[minimum(boundsLow),maximum(boundsHigh)],file);
	end
	writedlm(file,convert(Array{Float16},mDc[:]));
end
dumpResults(mref,[],"ref",[],[]);
dumpResults(mtrue,[],"true",[],[]);

###########################################################################################################################
###########################################################################################################################
################ DC Resistivity Alone Code ################################################################################
###########################################################################################################################

if invertDC
	Mfwd    = getRegularMesh(MinvNodal.domain,nfwd);
	sigtrue = velToConductMod(mtrue,3.0,0.1,1.0)[1]
	sigref  = velToConductMod(mref,3.0,0.1,1.0)[1]
	sigback      = sigref - Iact*(Iact'*sigref)

	@everywhere begin
		modfunDC(m) = velToConductMod(m,3.0,0.1,1.0);
	end

	println("------- setup pFor -----")
	Ainv = getMUMPSsolver()
	Mesh2Mesh  = getInterpolationMatrix(Minv,Mfwd)'
	redoSetup = true; # load setup if possible
	redoData =  true; # if possible, do not recompute data

	println("---------  sending out pMis ---------- ")

	if nworkers()>1
		# Julia sets BLAS #threads to be 1 if workers are used.
		set_num_threads(nworkers());
	end
	pFor,Sources,Receivers,Wd,dobsDC,dobs0 = setupDivSigGrad(redoSetup,redoData,Mfwd,Ainv,Mesh2Mesh,sigtrue,sigref);
	Wt   = 1./(mean(abs.(vec(dobsDC)))/2+abs.(dobsDC));
	dobsDC += 0.01*randn(size(dobsDC))*mean(abs.(vec(dobsDC)))
	println("adding noise of magnitude: ",0.01*randn()*mean(abs.(vec(dobsDC))), " Compared to ",mean(abs.(vec(dobsDC))));
	gloc = GlobalToLocal(Iact'*Mesh2Mesh,Mesh2Mesh'*sigback);

	if invertJoint
		worker = workers()[1];
		pMis = initRemoteChannel(getMisfitParam,worker, pFor,Wt,dobsDC,SSDFun,modfunDC,gloc);
		pMis = [pMis];
		modfun = identityMod;
	elseif invertDC
		pMis = getMisfitParam(pFor,Wt,dobsDC,SSDFun,identityMod,gloc);
		modfun = modfunDC;
	end
	pInv = getInverseParam(Minv,modfun,
                regfun,alpha,mref[:],
                boundsLow,boundsHigh,
                maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
                minUpdate=minUpdate, maxIter = maxIter,HesPrec=HesPrec)
	println("-------- solve DivSigGrad ---------")
end
###########################################################################################################################
###########################################################################################################################
################ Travel Time Tomography Code ##############################################################################
###########################################################################################################################

if invertEik
	jump 	= 10;
	offset 	= 128;
	pad 	= 2;

	slowref 	= velocityToSlowSquared(mref)[1];
	slowback	= slowref - Iact*(Iact'*slowref);

	dataFilenamePrefix = string("DATA_SEG",tuple((MinvNodal.n+1)...));
	resultsFilename    = string("travelTimeInvSEG",tuple((MinvNodal.n+1)...),".dat");
	prepareTravelTimeDataFiles(mtrue,MinvNodal,mref,boundsHigh,boundsLow,dataFilenamePrefix,pad,jump,offset);

	if invertJoint
		pMisEik = setupTravelTimeTomography(dataFilenamePrefix, MinvNodal,Iact,slowback,velocityToSlowSquared)[3];
		modfun = identityMod;
		pMis = [pMis;pMisEik];
	elseif invertEik
		pMis = setupTravelTimeTomography(dataFilenamePrefix, MinvNodal,Iact,slowback,identityMod)[3];
		modfun = velocityToSlowSquared;
	end
	pInv = getInverseParam(Minv,modfun,regfun,alpha,mref[:],boundsLow,boundsHigh,
                     maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate=minUpdate, maxIter = maxIter,HesPrec=HesPrec);
end
#### Projected Gauss Newton
mc,Dc,flag = projGNCG(copy(mref[:]),pInv,pMis,dumpResults = dumpResults);
#### Barrier Gauss Newton
# mc,Dc,flag = barrierGNCG(copy(mref[:]),pInv,pMis,dumpResults = dumpResults,epsilon = 0.3);
mv("jInv.out",string(filenamePrefix,".out"),remove_destination=true);
t = 1.0;


text:
@info "Create EOS table"

@info "Creating file \"eos.quantities\"..."

# regular_quantities = [1, 2, 12, 20, 21, 22, 23, 24]
# additional_quantities = []
# derivative_quantities = []
# pair_quantities = []
# quadruple_quantities = []
# microscopic_quantities = []
# error_quantities = []

# These are all the quantities (why not?)
regular_quantities = collect(1:24)
additional_quantities = collect(1:1)
derivative_quantities = collect(1:10)
pair_quantities = [10, 11, 4002, 0]
quadruple_quantities = [1]
microscopic_quantities = [10040, 11040, 10050, 11050]
error_quantities = collect(1:8)
output_format = 0               # 1=ASCII, 0=HDF5

open("eos.quantities", "w") do fh
    println(fh,
            "# number of regular, additional and derivative quantities (see table 7.1)")
    println(fh,
            join([length(regular_quantities), length(additional_quantities),
                  length(derivative_quantities)], " "))
    println(fh, "# indices of regular, additional and derivative quantities")
    println(fh,
            join([regular_quantities; additional_quantities;
                  derivative_quantities], " "))
    println(fh,
            "# number of pairs and quadruples for composition data (see table 3.2/3.3)")
    println(fh,
            join([length(pair_quantities), length(quadruple_quantities)], " "))
    println(fh, "# indices for pairs and quadruples for composition data")
    println(fh, join([pair_quantities; quadruple_quantities], " "))
    println(fh, "# number of microscopic quantities (see section 4.2.4)")
    println(fh, length(microscopic_quantities))
    println(fh, "# indices of microscopic quantities")
    println(fh, join(microscopic_quantities, " "))
    println(fh, "# number of error quantities")
    println(fh, length(error_quantities))
    println(fh, "# indices of error quantities")
    println(fh, join(error_quantities, " "))
    println(fh, "# format of output file, ASCII (1) or HDF5 (else)")
    println(fh, output_format)
    return nothing
end

@info "Creating file \"eos.parameters\"..."

# This selects all the points in the original table, presumably without interpolation (why not?)

file_t = split(read("eos.t", String))
log_t = parse(Int, file_t[1])
num_t = parse(Int, file_t[2])
min_t = parse(Float64, file_t[3])
max_t = parse(Float64, file_t[end])
@assert length(file_t) == num_t + 2

file_nb = split(read("eos.nb", String))
log_nb = parse(Int, file_nb[1])
num_nb = parse(Int, file_nb[2])
min_nb = parse(Float64, file_nb[3])
max_nb = parse(Float64, file_nb[end])
@assert length(file_nb) == num_nb + 2

file_yq = split(read("eos.yq", String))
log_yq = parse(Int, file_yq[1])
num_yq = parse(Int, file_yq[2])
min_yq = parse(Float64, file_yq[3])
max_yq = parse(Float64, file_yq[end])
@assert length(file_yq) == num_yq + 2

open("eos.parameters", "w") do fh
    println(fh, "# order of interpolation in first, second and third index")
    println(fh, "1 1 1")
    println(fh,
            "# calculation of beta-equilibrium (1: yes, else: no) and for given entropy (1: yes, else: no)")
    println(fh, "0 0")
    println(fh,
            "# tabulation scheme (0 = explicit listing, 1 = loops, see manual)")
    println(fh, "1")
    println(fh,
            "# parameter values (first, second and third index) depending on tabulation scheme")
    println(fh, join([min_t, min_nb, min_yq], " "))
    println(fh, join([max_t, max_nb, max_yq], " "))
    println(fh, join([num_t, num_nb, num_yq], " "))
    println(fh, join([log_t, log_nb, log_yq], " "))
    return nothing
end

@info "Generate EOS file..."
run(pipeline(`echo 3`, `/Users/eschnett/src/codehdf5/compose`))

# h5repack -v -f SHUF -f GZIP=9 eoscompose.h5 eoscompose_compressed.h5

@info "Done."

ext: jl


text:
import Base: rand!, randn!

const libcurand = Libdl.find_library(["libcurand"], [])
RNG=0

function rng(init=false)
    global RNG
    if RNG==0 || init
        ptr = Ptr{Void}[0]
        # CURAND_RNG_PSEUDO_DEFAULT = 100, ///< Default pseudorandom generator
        @assert 0==ccall((:curandCreateGenerator,libcurand),Cint,(Ptr{Void},Cint),ptr,100)
        RNG = ptr[1]
    end
    return RNG
end

# our version of srand sets both gpu and cpu
function setseed(n::Integer)
    # need to regenerate RNG for the seed to take effect for some reason
    @assert 0==ccall((:curandSetPseudoRandomGeneratorSeed,libcurand),Cint,(Ptr{Void},Culonglong),rng(true),n)
    srand(n)
end

rand!(a::CudaArray{Float32})=(@assert 0==ccall((:curandGenerateUniform,libcurand),Cint,(Ptr{Void},Ptr{Cfloat},Csize_t),rng(),a,length(a)); a)
rand!(a::CudaArray{Float64})=(@assert 0==ccall((:curandGenerateUniformDouble,libcurand),Cint,(Ptr{Void},Ptr{Cdouble},Csize_t),rng(),a,length(a)); a)
rand!(a::Union{CudaArray{UInt32},CudaArray{Int32}})=(@assert 0==ccall((:curandGenerate,libcurand),Cint,(Ptr{Void},Ptr{Cuint},Csize_t),rng(),a,length(a)); a)
rand!(a::Union{CudaArray{UInt64},CudaArray{Int64}})=(@assert 0==ccall((:curandGenerate,libcurand),Cint,(Ptr{Void},Ptr{Cuint},Csize_t),rng(),a,2*length(a)); a)

# These are a pain, curand insists that array length should be even!
function randn!(a::CudaArray{Float32},mean=0f0,stddev=1f0)
    if length(a) % 2 == 0
        @assert 0==ccall((:curandGenerateNormal,libcurand),Cint,(Ptr{Void},Ptr{Cfloat},Csize_t,Cfloat,Cfloat),rng(),a,length(a),mean,stddev)
    else
        @assert 0==ccall((:curandGenerateNormal,libcurand),Cint,(Ptr{Void},Ptr{Cfloat},Csize_t,Cfloat,Cfloat),rng(),a,length(a)-1,mean,stddev)
        copysync!(a,length(a),Float32[randn()*stddev+mean],1,1)
    end
    return a
end

function randn!(a::CudaArray{Float64},mean=0e0,stddev=1e0)
    if length(a) % 2 == 0
        @assert 0==ccall((:curandGenerateNormalDouble,libcurand),Cint,(Ptr{Void},Ptr{Cdouble},Csize_t,Cdouble,Cdouble),rng(),a,length(a),mean,stddev)
    else
        @assert 0==ccall((:curandGenerateNormalDouble,libcurand),Cint,(Ptr{Void},Ptr{Cdouble},Csize_t,Cdouble,Cdouble),rng(),a,length(a)-1,mean,stddev)
        copysync!(a,length(a),Float64[randn()*stddev+mean],1,1)
    end
    return a
end

# randn!{T}(a::Array{T}, mean=zero(T), std=one(T))=(for i=1:length(a); a[i] = mean + std * randn(); end; a)
# rand!(a::BaseArray, x0, x1)=(rand!(a); axpb!(x1-x0, x0, a); a)


text:
"""
Optimizes f using Newton's method and exact Hessians.  For now, it is
not clear whether this or BFGS is better, so it is kept as a separate function.

Args:
  - f: A function that takes elbo args and constrained coordinates
       (e.g. DeterministicVI.elbo)
  - ea: Constrained initial ElboArgs
  - lbs: An array of lower bounds (in the transformed space)
  - ubs: An array of upper bounds (in the transformed space)
  - omitted_ids: Omitted ids from the _unconstrained_ parameterization
      (i.e. elements of free_ids).
  - xtol_rel: X convergence
  - ftol_abs: F convergence
  - verbose: Print detailed output

Returns:
  - iter_count: The number of iterations taken
  - max_f: The maximum function value achieved
  - max_x: The optimal function input
  - ret: The return code of optimize()
"""
function maximize_f(f::Function,
                    ea::ElboArgs,
                    transform::DataTransform;
                    omitted_ids=Int[],
                    xtol_rel=1e-7,
                    ftol_abs=1e-6,
                    verbose=false,
                    max_iters=50,
                    fast_hessian=true)
    # Make sure the model parameters are within the transform bounds
    enforce_bounds!(ea.vp, ea.active_sources, transform)
    @assert ea.active_sources == transform.active_sources

    kept_ids = setdiff(1:length(UnconstrainedParams), omitted_ids)
    x_length = length(kept_ids) * transform.active_S
    x_size = (length(kept_ids), transform.active_S)

    Sa = length(ea.active_sources)

    f_evals = 0
    const print_every_n = 10

    function f_wrapped_nocache(x::Vector{Float64})
        f_evals += 1

        # Evaluate in the constrained space and then unconstrain again.
        transform.array_to_vp!(reshape(x, x_size), ea.vp, omitted_ids)
        f_res = f(ea)

        # TODO: Add an option to print either the transformed or
        # free parameterizations.

        if verbose || (f_evals % print_every_n == 0)
            Log.info("f_evals: $(f_evals) value: $(f_res.v[1])")
        end

        if verbose
            S = length(iter_vp)

            if length(iter_vp[1]) == length(ids_names)
                state_df = DataFrames.DataFrame(names=ids_names)
            elseif length(iter_vp[1]) == length(ids_free_names)
                state_df = DataFrames.DataFrame(names=ids_free_names)
            else
                state_df = DataFrames.DataFrame(
                    names=[ "x$i" for i=1:length(iter_vp[1, :])])
            end

            for s=1:S
                state_df[Symbol(string("val", s))] = iter_vp[s]
            end

            for s=1:S
                state_df[Symbol(string("grad", s))] = f_res.d[:, s]
            end

            Log.info(repr(state_df))
            Log.info("=======================================\n")
        end

        transform.transform_sensitive_float(f_res, ea.vp, ea.active_sources)
    end

    last_sf = zero_sensitive_float(UnconstrainedParams, Sa)
    last_x = [ NaN ]

    function f_wrapped_cached(x::Vector{Float64})
        if x != last_x
            last_x = deepcopy(x)
            last_sf = deepcopy(f_wrapped_nocache(x))
        end

        last_sf
    end

    function neg_f_value{T <: Number}(x::Vector{T})
        @assert length(x) == x_length
        -f_wrapped_cached(x).v[1]
    end

    function neg_f_grad!{T <: Number}(x::Vector{T}, grad::Vector{T})
        @assert length(x) == x_length
        res = f_wrapped_cached(x)
        if length(grad) > 0
            svs = [res.d[kept_ids, si] for si in 1:transform.active_S]
            grad[:] = reduce(vcat, svs)
        end
        grad .*= -1
    end

    function neg_f_hessian!{T <: Number}(x::Vector{T}, hess::Matrix{T})
        @assert length(x) == x_length
        res = f_wrapped_cached(x)
        all_kept_ids = Int[]
        for sa=1:transform.active_S
            append!(all_kept_ids, kept_ids + (sa - 1) * length(kept_ids))
        end
        sub_hess = res.h[all_kept_ids, all_kept_ids]
        hess[:, :] = -1 .* 0.5 * (sub_hess + sub_hess')
    end

    x0 = transform.vp_to_array(ea.vp, omitted_ids)
    tr_method = Optim.NewtonTrustRegion(
                    initial_delta=10.0,
                    delta_hat=1e9,
                    eta=0.1,
                    rho_lower=0.25,
                    rho_upper=0.75)

    options = Optim.OptimizationOptions(;
        x_tol = xtol_rel, f_tol = ftol_abs, g_tol = 1e-8,
        iterations = max_iters, store_trace = verbose,
        show_trace = false, extended_trace = verbose)

    nm_result = Optim.optimize(neg_f_value,
                               neg_f_grad!,
                               neg_f_hessian!,
                               x0[:],
                               tr_method,
                               options)

    transform.array_to_vp!(reshape(nm_result.minimum, size(x0)),
                           ea.vp, omitted_ids)
    max_f = -1.0 * nm_result.f_minimum
    max_x = nm_result.minimum

    Log.info(string("got $max_f at $max_x after $f_evals function evaluations ",
            "($(nm_result.iterations) Newton steps)\n"))
    f_evals, max_f, max_x, nm_result
end

function maximize_f(f::Function,
                    ea::ElboArgs;
                    loc_width=1.5e-3,
                    omitted_ids=Int[],
                    xtol_rel=1e-7,
                    ftol_abs=1e-6,
                    verbose=false,
                    max_iters=50,
                    fast_hessian=true)
    transform = get_mp_transform(ea.vp, ea.active_sources, loc_width=loc_width)

    maximize_f(f, ea, transform;
                omitted_ids=omitted_ids,
                xtol_rel=xtol_rel,
                ftol_abs=ftol_abs,
                verbose=verbose,
                max_iters=max_iters,
                fast_hessian=fast_hessian)
end


text:
function check_liquid(client, invest)
    account = client.get_account()
    for b in account["balances"]
        if b["asset"] == invest["symbol"][end-2:end]
            if parse(Float64, b["free"]) > parse(Float64, invest["amount"])
                return true
            else
                return false
            end
        end
    end
    return false
end
export check_liquid


text:
module MainMenuBars

Base.Experimental.@optlevel 2

using Redux
using CImGui
using ..Menus
import ..Menus: AbstractMenuAction

## actions
"""
    AbstractGenericMainMenuBarAction <: AbstractSyncAction
Abstract supertype for all kinds of main menu bar actions. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMainMenuBarAction <: AbstractSyncAction end

"""
    AbstractMainMenuBarAction <: AbstractGenericMainMenuBarAction
Abstract supertype for all main menu bar actions. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMainMenuBarAction <: AbstractGenericMainMenuBarAction end

# for determining whether we should apply the reducer function to certain action types.
struct Valid end
struct Invalid end

is_valid(a::AbstractSyncAction) = Invalid
is_valid(a::AbstractGenericMainMenuBarAction) = Invalid
is_valid(a::AbstractMainMenuBarAction) = Valid

_get_label(a::AbstractSyncAction, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_ACTION_LABEL"
_get_label(a::AbstractSyncAction, ::Type{Valid}) = a.label
get_label(a::AbstractSyncAction) = _get_label(a, is_valid(a))

"""
    Rename(label, new_label)
Change widget's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractMainMenuBarAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractMainMenuBarAction
    label::String
    is_triggered::Bool
end

"""
    EditMenus(label, action::AbstractMenuAction)
Edit manus.
"""
struct EditMenus{T<:AbstractMenuAction} <: AbstractMainMenuBarAction
    label::String
    action::T
end

"""
    Show(label)
Show the menubar.
"""
struct Show <: AbstractMainMenuBarAction
    label::String
end

"""
    Hide(label)
Hide the menubar.
"""
struct Hide <: AbstractMainMenuBarAction
    label::String
end

## state
"""
    AbstractGenericMainMenuBar <: AbstractImmutableState
Abstract supertype for all kinds of main menu bars. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMainMenuBar <: AbstractImmutableState end

"""
    AbstractMainMenuBar <: AbstractGenericMainMenuBar
Abstract supertype for all main menu bars. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMainMenuBar <: AbstractGenericMainMenuBar end

is_valid(s::AbstractImmutableState) = Invalid
is_valid(s::AbstractGenericMainMenuBar) = Invalid
is_valid(s::AbstractMainMenuBar) = Valid

_get_label(s::AbstractImmutableState, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_STATE_LABEL"
_get_label(s::AbstractImmutableState, ::Type{Valid}) = s.label
get_label(s::AbstractImmutableState) = _get_label(s, is_valid(s))

"""
    MainMenuBar(label::AbstractString, menus = [], is_hidden = false)
A main menu bar state which contains a label a.k.a the identifier, a list of menus,
a flag value `is_hidden` and a flag value `is_triggered` that records the state of
the latest poll events.
"""
struct MainMenuBar <: AbstractMainMenuBar
    label::String
    menus::Vector{Menus.Menu}
    is_hidden::Bool
    is_triggered::Bool
end
MainMenuBar(label::AbstractString, menus = [], is_hidden = false) =
    MainMenuBar(label, menus, is_hidden, false)

is_hidden(s::MainMenuBar) = s.is_hidden
is_triggered(s::MainMenuBar) = s.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::MainMenuBar, a::Rename) = MainMenuBar(a.new_label, s.menus, s.is_hidden, s.is_triggered)
reducer(s::MainMenuBar, a::SetTriggeredTo) = MainMenuBar(s.label, s.menus, s.is_hidden, a.is_triggered)
reducer(s::MainMenuBar, a::Show) = MainMenuBar(s.label, s.menus, false, s.is_triggered)
reducer(s::MainMenuBar, a::Hide) = MainMenuBar(s.label, s.menus, true, s.is_triggered)

reducer(s::MainMenuBar, a::EditMenus) =
    MainMenuBar(s.label, Menus.reducer(s.menus, a.action), s.is_hidden, s.is_triggered)

## UI
"""
    (::MainMenuBar)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity) -> Bool
Return `true` when triggered/activated.
`get_state` is a router function that tells how to find the target state from `store`.
`chain_action` is for chaining upstream actions.
"""
function (::MainMenuBar)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity)
    s = get_state(store)
    is_activated = CImGui.BeginMainMenuBar()
    dispatch!(
        store,
        SetTriggeredTo(s.label, is_activated) |> chain_action,
    )
    if is_activated && !s.is_hidden
        for (i, menu) in enumerate(s.menus)
            menu(
                store,
                x->get_state(x).menus[i],
                x->EditMenus(s.label, x) |> chain_action,
            )
        end
        CImGui.EndMainMenuBar()
    end
    return is_activated
end

end # module


text:
module ContinuumArrays
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, diff,
                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
                getproperty
import Base.Broadcast: materialize, BroadcastStyle, broadcasted
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
                        adjointlayout, LdivApplyStyle, arguments, _arguments, call, broadcastlayout, lazy_getindex,
                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles
import LinearAlgebra: pinv
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
import FillArrays: AbstractFill, getindex_value, SquareEye

import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, quasildivapplystyle, _factorize

export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine

####
# Interval indexing support
####
struct AlephInfinity{N} <: Integer end

==(::AlephInfinity, ::Int) = false
==(::Int, ::AlephInfinity) = false

*(::AlephInfinity{N}, ::AlephInfinity{N}) where N = AlephInfinity{N}()

const ℵ₁ = AlephInfinity{1}()

show(io::IO, F::AlephInfinity{1}) where N =
    print(io, "ℵ₁")

const QMul2{A,B} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B}}
const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}

cardinality(::AbstractInterval) = ℵ₁
*(ℵ::AlephInfinity) = ℵ

Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)

for find in (:findfirst, :findlast)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::Inclusion) = f.x in d.domain ? f.x : nothing
end

function findall(f::Base.Fix2{typeof(isequal)}, d::Inclusion)
    r = findfirst(f,d)
    r === nothing ? eltype(d)[] : [r]
end

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval})
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

BroadcastStyle(::Type{<:Inclusion{<:Any,<:AbstractInterval}}) = LazyQuasiArrayStyle{1}()
BroadcastStyle(::Type{<:QuasiAdjoint{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()
BroadcastStyle(::Type{<:QuasiTranspose{<:Any,<:Inclusion{<:Any,<:AbstractInterval}}}) = LazyQuasiArrayStyle{2}()

###
# Maps
###

# Affine map represents A*x .+ b
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: AbstractQuasiVector{T} end

struct AffineQuasiVector{T,AA,X,B} <: AbstractAffineQuasiVector{T,AA,X,B}
    A::AA
    x::X
    b::B
end

AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
    AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)

AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)

AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)

axes(A::AbstractAffineQuasiVector) = axes(A.x)
affine_getindex(A, k) = A.A*A.x[k] .+ A.b
getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
    @boundscheck A.x[k] # throws bounds error if k ≠ x
    A
end

getindex(A::AbstractAffineQuasiVector, ::Colon) = copy(A)

copy(A::AbstractAffineQuasiVector) = A

inbounds_getindex(A::AbstractAffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
isempty(A::AbstractAffineQuasiVector) = isempty(A.x)
==(a::AbstractAffineQuasiVector, b::AbstractAffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b

BroadcastStyle(::Type{<:AbstractAffineQuasiVector}) = LazyQuasiArrayStyle{1}()

for op in(:*, :\, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
end
for op in(:/, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
end

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(a, x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AbstractAffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)

function checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector)
    @_propagate_inbounds_meta
    isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r)))
end

affine_igetindex(d, x) = d.A \ (x .- d.b)
igetindex(d::AbstractAffineQuasiVector, x) = affine_igetindex(d, x)

for find in (:findfirst, :findlast, :findall)
    @eval $find(f::Base.Fix2{typeof(isequal)}, d::AbstractAffineQuasiVector) = $find(isequal(igetindex(d, f.x)), d.x)
end

minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)

union(d::AbstractAffineQuasiVector) = Inclusion(minimum(d)..maximum(d))

struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
    domain::D
    range::R
end

AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} =
    AffineMap{promote_type(T,V), typeof(domain),typeof(range)}(domain,range)

measure(x::Inclusion) = last(x)-first(x)

function getproperty(A::AffineMap, d::Symbol)
    d == :x && return A.domain
    d == :A && return measure(A.range)/measure(A.domain)
    d == :b && return (last(A.domain)*first(A.range) - first(A.domain)*last(A.range))/measure(A.domain)
    getfield(A, d)
end

function getindex(A::AffineMap, k::Number)
    # ensure we exactly hit range
    k == first(A.domain) && return first(A.range)
    k == last(A.domain) && return last(A.range)
    affine_getindex(A, k)
end

function igetindex(A::AffineMap, k::Number)
    # ensure we exactly hit range
    k == first(A.range) && return first(A.domain)
    k == last(A.range) && return last(A.domain)
    affine_igetindex(A, k)
end

first(A::AffineMap) = first(A.range)
last(A::AffineMap) = last(A.range)

affine(a::AbstractQuasiVector, b::AbstractQuasiVector) = AffineMap(a, b)
affine(a, b::AbstractQuasiVector) = affine(Inclusion(a), b)
affine(a::AbstractQuasiVector, b) = affine(a, Inclusion(b))
affine(a, b) = affine(Inclusion(a), Inclusion(b))

const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}

sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V

include("operators.jl")
include("bases/bases.jl")

end


text:
"Singly implicit Runge-Kutta integrator."
immutable IntegratorSIRK{T} <: Integrator{T}
    equation::ODE{T}
    tableau::TableauSIRK{T}
    Δt::T

    x::Array{T,1}
    X::Array{T,2}
    Y::Array{T,2}
    F::Array{T,2}

    function IntegratorSIRK(equation, tableau, Δt)
        D = equation.d
        S = tableau.s
        new(equation, tableau, Δt, zeros(T,D), zeros(T,D,S), zeros(T,D,S), zeros(T,D,S))
    end
end

function IntegratorSIRK(equation::Equation, tableau::TableauSIRK, Δt)
    T = eltype(equation.q₀)
    IntegratorSIRK{T}(equation, tableau, Δt)
end

"Integrate ODE with singly implicit Runge-Kutta integrator."
function integrate!(int::IntegratorSIRK, s::SolutionODE)
    # TODO
end

"Integrate partitioned ODE with singly implicit Runge-Kutta integrator."
function integrate!(int::IntegratorSIRK, s::SolutionPODE)
    # TODO
end


text:
using CSV
using DataFrames
using JLD
using Images
cd("Examples/Cactus")
input_file = "train.csv"
df = CSV.read(input_file)

first(df, 10)
describe(df)
df = sort(df, [:id])
images = df[!, :id]
train_Y = df[!, :has_cactus]
train_data = load("dataset.jld")
train_data

# d = train_data["test_X_orig"]
# size(d)[4]
# vcat(d[:, :, :, 1]...)
# for i = 0 : size(d)[4]
#     d[:, :, :, i] =
# end
# d2 = (reshape(d, :, size(d)[4]))
# d2[:, 1]

function reshape_image_array(image_array::Array{Float32,4})
    return reshape(image_array, :, size(image_array)[4])
end
function reshape_image_array_reverse(image_array::Array{Float32,2})
    return reshape(image_array, 32, 32, 3, :)
end

test_X = reshape_image_array(train_data["test_X_orig"])
train_X = reshape_image_array(train_data["train_X_orig"])
train_Y = reshape(train_Y, (1, size(train_Y)[1]))
# test_Y = reshape(test_Y, (1, size(test_Y)[1]))

include("../../NeuralNetwork.jl")
# neural_network_dense(train_X, train_Y, [3072, 1, 1], 0.05)
parameters, activations = neural_network_dense(train_X, train_Y, [3072, 10, 10, 1], 100, 0.1)
predicts, accuracy = predict(train_X, train_Y, parameters, activations)

a = reshape_image_array_reverse(train_X[:, [1,7]])
a2 = a[:,:,:,2]
a3 = permuteddimsview(a2, [3, 1, 2])
a3 = Array{Float32}(a3)
colorview(RGB, a3)


text:
@testset "Algebra and Properties" begin
    @test USD == FinancialDSL.Currencies.Currency(:USD)
    @test FinancialDSL.Currencies.currency_symbol(FinancialDSL.Currencies.Currency(:X)) === :X
    @test 1BRL == FinancialDSL.Currencies.Cash{FinancialDSL.Currencies.Currency{:BRL}, Int}(1)
    @test 10.2BRL > 10BRL
    @test 1BRL + 20.5BRL == 21.5BRL
    @test 100BRL * 0.5BRL == 50BRL
    @test 100BRL / 10BRL == 10
    @test 100BRL / 10 == 10BRL
    @test isapprox(10.334BRL, 10.334BRL)
    @test FinancialDSL.Currencies.cashcurrency(10.0BRL) == BRL
    @test FinancialDSL.Currencies.cashvalue(10.0BRL) == 10.0
end

@testset "exch" begin
    provider = FinancialDSL.Core.FixedScenario()
    provider[USD] = 3.5BRL
    provider[EUR] = 4.0BRL

    usd_price = FinancialDSL.Currencies.cashvalue(provider[USD])
    @test usd_price == 3.5

    eur_price = FinancialDSL.Currencies.cashvalue(provider[EUR])
    @test eur_price == 4.0

    @test FinancialDSL.Currencies.cashcurrency(provider[EUR]) == BRL

    @test FinancialDSL.Core.exch(provider, BRL, BRL) == 1.0
    @test FinancialDSL.Core.exch(provider, EUR, EUR) == 1.0
    @test FinancialDSL.Core.exch(provider, FinancialDSL.Currencies.Currency(:X), FinancialDSL.Currencies.Currency(:X)) == 1.0
    @test FinancialDSL.Core.exchcash(provider, BRL, BRL) == 1.0BRL
    @test FinancialDSL.Core.exchcash(provider, EUR, EUR) == 1.0EUR

    @test FinancialDSL.Core.exch(provider, BRL, USD) == 1.0 / usd_price
    @test FinancialDSL.Core.exch(provider, BRL, EUR) == 1.0 / eur_price
    @test FinancialDSL.Core.exchcash(provider, BRL, EUR) == (1.0 / eur_price)*EUR

    @test FinancialDSL.Core.exch(provider, USD, BRL) == usd_price
    @test FinancialDSL.Core.exch(provider, EUR, BRL) == eur_price
    @test FinancialDSL.Core.exchcash(provider, EUR, BRL) == eur_price*BRL

    @test FinancialDSL.Core.exch(provider, EUR, USD) == eur_price / usd_price
    @test FinancialDSL.Core.exchcash(provider, EUR, USD) == (eur_price / usd_price)*USD
end


text:
using Plots, Test

@testset "Series Attributes" begin
    pl = plot([[1, 2, 3], [2, 3, 4]], lw = 5)
    @test hline!(deepcopy(pl), [1.75])[1].series_list[3][:label] ==
          hline!(deepcopy(pl), [1.75], z_order = :front)[1].series_list[3][:label] ==
          "y3"
    @test hline!(deepcopy(pl), [1.75], z_order = :back)[1].series_list[1][:label] == "y3"
    @test hline!(deepcopy(pl), [1.75], z_order = 2)[1].series_list[2][:label] == "y3"
end

@testset "Axis Attributes" begin
    pl = @test_nowarn plot(; tickfont = font(10, "Times"))
    for axis in (:xaxis, :yaxis, :zaxis)
        @test pl[1][axis][:tickfontsize] == 10
        @test pl[1][axis][:tickfontfamily] == "Times"
    end
end


text:
using CxxWrap
using CxxWrap.CxxWrapCore: SafeCFunction
using BenchmarkTools

using Emporos
using Emporos.IBKR
using Emporos.IBKR: IBKRCore

const account = "U6403725"

function main()
    ib = IB(port=7496)
    ib() do
        IBKRCore.reqManagedAccounts(ib.client)
        IBKRCore.reqAccountSummary(ib.client, 1, "All", "TotalCashValue")
        sleep(1)
        IBKRCore.reqAccountSummary(ib.client, 1, "All", "TotalCashValue")
        IBKRCore.reqAccountSummary(ib.client, 1, "All", "TotalCashValue")
        sleep(5)
        IBKRCore.reqAccountSummary(ib.client, 1, "All", "TotalCashValue")
    end
    @info "DONE main"
end

main()


text:
module HighPrecisionNonsymmetricEigenproblem

using StructArrays, VectorizationBase, LinearAlgebra

"""
Passed into the Schur factorization function if you do not wish to have the Schur vectors.
"""
struct NotWanted end

abstract type SmallRotation end

"""
Given's rotation acting on rows i:i+1
"""
struct Rotation2{Tc,Ts} <: SmallRotation
    c::Tc
    s::Ts
    i::Int
end

"""
Two Given's rotations acting on rows i:i+2. This could also be implemented as one Householder
reflector!
"""
struct Rotation3{Tc,Ts} <: SmallRotation
    c₁::Tc
    s₁::Ts
    c₂::Tc
    s₂::Ts
    i::Int
end

include("tsqr.jl")
include("hessenberg.jl")
include("givens.jl")
include("soa_multifloat.jl")
include("schur.jl")
include("utils.jl")
include("multiprecision.jl")

end


text:
#================
1D seismic models
================#
"""
Abstract type for radially-symmetric models of a planet.  1D seismic models should
be subtypes of this.
"""
abstract type SeisModel1D <: SeisModel end

const model_variables_SeisModel1D = (:vp, :vs, :density, :vph, :vpv, :vsh, :vsv,
    :eta, :Qμ, :Qκ)
const model_names_SeisModel1D = ("isotropic average Vp", "isotropic average Vs",
    "density", "Vph", "Vpv", "Vsh", "Vsv", "η", "Qμ", "Qκ")
const model_units_SeisModel1D = (" (km/s)", " (km/s)", " (g/cm^3)", " (km/s)",
    " (km/s)", " (km/s)", " (km/s)", "", "", "")

"""
    findlayer(m::SeisModel1D, r) -> layer_index

Return the layer index in which radius `r` km is located in the SeisModel1D `m`.
"""
function findlayer(m::SeisModel1D, r::Real)
    r < 0. && throw(DomainError(r, "radius cannot be negative"))
    r > m.a && throw(DomainError(r, "radius is greater than Earth radius for model ($(m.a) km)"))
    for i in 1:length(m.r)
        @inbounds r < m.r[i] && return i
    end
    length(m.r)
end

# Implementations of vp, vs, etc., evaluate here
include("prempolymodel.jl")
include("steppedlayeredmodel.jl")
include("linearlayeredmodel.jl")

## Derived quantities
"""
    bulk_modulus(m::SeisModel1D, r; depth=false) -> K

Return the bulk modulus `K` in GPa at radius `r` km in the model `m`.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function bulk_modulus(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    density(m, r)*vp(m, r)^2 - 4/3*shear_modulus(m, r)
end

"Newton's gravitational constant G = 6.67428e-11"
const NewtonG = 6.67428e-11

"""
    gravity(m::SeisModel1D, r; depth=false) -> g

Return the acceleration due to gravity, `g`, in m/s^2 at radius `r` in km.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function gravity(m::SeisModel1D, r; depth=false)
    depth && (r = radius(m, r))
    (r == 0) ? 0. : NewtonG*mass(m, r)/(r*1.e3)^2
end

# Individual implementations of mass are in model types' respective files
"""
    mass(m::SeisModel1D, r; depth=false) -> mass

Return the mass in kg between the centre of the model and the radius `r` km.

If `depth` is `true`, `r` is treated as a depth in km instead.

    mass(m) -> mass

Return the mass for the whole body in kg (between the centre of the model
and the surface).
"""
mass(m::SeisModel1D) = mass(m, surface_radius(m))

"""
    pressure(m::SeisModel1D, r; depth=false) -> p

Return the pressure `p` in Pa at radius `r` km.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function pressure(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    f(r) = pressure_integration_func(m, r)
    1.e3*quadgk(f, r, m.a)[1]
end
pressure_integration_func(m::SeisModel1D, r) = 1.e3*density(m, r)*gravity(m, r)

"""
    poissons_ratio(m, r; depth=false) -> ν

Return the Poisson's ratio for the model `m` given a radius `r` in km.
The calculation uses the isotropic average velocities at `r`.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function poissons_ratio(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    G = shear_modulus(m, r)
    K = bulk_modulus(m, r)
    (3K - 2G)/(6K + 2G)
end

"""
    shear_modulus(m::SeisModel1D, r; depth=false) -> μ

Return the shear modulus `μ` (often also called G) in GPa at radius `r` in the
model `m`.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function shear_modulus(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    vs(m, r)^2*density(m, r)
end

"""
    surface_mass(m::SeisModel1D, r; depth=false) -> mass

Return the mass in kg betwen radius `r` km and the surface.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function surface_mass(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    mass(m, m.a) - mass(m, r)
end

"""
    youngs_modulus(m, r; depth=false) -> E

Return the Young's modulus `E` in GPa for the model `m` given a radius `r` in km.

If `depth` is `true`, `r` is treated as a depth in km instead.
"""
function youngs_modulus(m::SeisModel1D, r; depth::Bool=false)
    depth && (r = radius(m, r))
    2*shear_modulus(m, r)*(1 + poissons_ratio(m, r))
end

"""
    moment_of_inertia(m, r0=0, r1=surface_radius(m)) -> I

Return the moment of interia `I` in kg m² for the model `m` between radii `r0`
and `r1` in km.
"""
moment_of_inertia(m::SeisModel1D, r0=0, r1=surface_radius(m)) =
    8/3*π*quadgk(r->1e3*density(m, r/1e3)*r^4, r0*1e3, r1*1e3)[1]


text:
function get_char()
	term = REPL.Terminals.TTYTerminal("xterm",stdin,stdout,stderr)
	REPL.Terminals.raw!(term,true)
	Base.start_reading(stdin)

    bb = bytesavailable(stdin)
    if bb > 0
        data = read(stdin, bb)
        if data[1] == UInt(3)
            println("Ctrl+C - exiting")
            exit()
        end
        # println("Got $bb bytes: $(string(data))")
		empty_stream = false
		char = Char.(data[end])
    else
		empty_stream = true
		char = nothing
	end
	t = REPL.TerminalMenus.terminal
	REPL.TerminalMenus.disableRawMode(t)
	return char, empty_stream
end


text:
struct ICA_PWM_model
    name::String #designator for saving model to posterior
    sources::Vector{Tuple{Matrix{AbstractFloat},Integer}} #vector of PWM signal sources (LOG PROBABILITY!!!) tupled with an index denoting the position of the first PWM base on the prior matrix- allows us to permute length and redraw from the appropriate prior position
    informed_sources::Integer #number of sources with informative priors- these are not subject to frequency sorting in model mergers
    source_length_limits::UnitRange{Integer} #min/max source lengths for init and permutation
    mix_matrix::BitMatrix # obs x sources bool matrix
    log_Li::AbstractFloat
    flags::Vector{String} #optional flags for search patterns
end

#ICA_PWM_model FUNCTIONS
ICA_PWM_model(name::String, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, mix_prior::Tuple{BitMatrix,AbstractFloat}, bg_scores::AbstractArray{AbstractFloat}, observations::AbstractArray{Integer}, source_length_limits::UnitRange{Integer}) = init_IPM(name, source_priors,mix_prior,bg_scores,observations,source_length_limits)

#MODEL INIT
function init_IPM(name::String, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, mix_prior::Tuple{BitMatrix,AbstractFloat}, bg_scores::AbstractArray{AbstractFloat}, observations::AbstractArray{Integer}, source_length_limits::UnitRange{Integer})
    T,O = size(observations)
    S=length(source_priors)
    obs_lengths=[findfirst(iszero,observations[:,o])-1 for o in 1:size(observations)[2]]
    sources=init_logPWM_sources(source_priors, source_length_limits)
    mix=init_mix_matrix(mix_prior,O,S)
    log_lh = IPM_likelihood(sources, observations, obs_lengths, bg_scores, mix)

   return ICA_PWM_model(name, sources, size(mix_prior[1],2), source_length_limits, mix, log_lh, ["init"])
end

                #init_IPM SUBFUNCS
                function init_logPWM_sources(prior_vector::Vector{Vector{Dirichlet{AbstractFloat}}}, source_length_limits::UnitRange{Integer})
                    srcvec = Vector{Tuple{Matrix{AbstractFloat},Integer}}()
                    prior_coord = 1
                        for (p, prior) in enumerate(prior_vector)
                            min_PWM_length=source_length_limits[1]
                            PWM_length = rand(min_PWM_length:length(prior)) #generate a PWM from a random subset of the prior
                            PWM = zeros(PWM_length,4)
                            prior_coord = rand(1:length(prior)-PWM_length+1) #determine what position on the prior to begin sampling from based on the PWM length
                            for (position, dirichlet) in enumerate(prior)
                                if position >= prior_coord #skip prior positions that are before the selected prior_coord
                                    sample_coord = min(position-prior_coord+1,PWM_length) #sample_coord is the position on the sampled PWM
                                    PWM[sample_coord, :] = rand(dirichlet) #draw the position WM from the dirichlet
                                    !isprobvec(PWM[sample_coord, :]) && throw(DomainError("Bad weight vec produced by init_sources! $(PWM[sample_coord,:])"))#make sure it's a valid probvec
                                end
                            end
                            push!(srcvec, (log.(PWM), prior_coord)) #push the source PWM to the source vector with the prior coord idx to allow drawing from the appropriate prior dirichlets on permuting source length
                        end
                    return srcvec
                end

                function init_mix_matrix(mix_prior::Tuple{BitMatrix,AbstractFloat}, no_observations::Integer, no_sources::Integer)
                    inform,uninform=mix_prior
                    if size(inform,2) > 0
                        @assert size(inform,1)==no_observations && size(inform,2)<=no_sources "Bad informative mix prior dimensions!"
                    end
                    @assert 0.0 <= uninform <=1.0 "Uninformative mix prior not between 0.0 and 1.0!"
                    mix_matrix = falses(no_observations, no_sources)
                    if size(inform,2)>0
                        mix_matrix[:,1:size(inform,2)]=inform
                    end
                    for index in CartesianIndices((1:no_observations,size(inform,2)+1:no_sources))
                        rand() <= uninform && (mix_matrix[index] = true)
                    end
                    return mix_matrix
                end

#LIKELIHOOD SCORING FUNCS
function IPM_likelihood(sources::Vector{Tuple{Matrix{AbstractFloat},Integer}}, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::AbstractArray{AbstractFloat}, mix::BitMatrix, revcomp::Bool=true, returncache::Bool=false, cache::Vector{AbstractFloat}=zeros(0), clean::Vector{Bool}=Vector(falses(size(observations)[2])))
    source_wmls=[size(source[1])[1] for source in sources]
    O = size(bg_scores)[2]
    obs_lhs=Vector{Vector{AbstractFloat}}()
    nt=Threads.nthreads()
    for t in 1:nt-1
        push!(obs_lhs,zeros(Int(floor(O/nt))))
    end
    push!(obs_lhs, zeros(Int(floor(O/nt)+(O%nt))))

    Threads.@threads for t in 1:nt
        opt = floor(O/nt) #obs per thread
        for i in 1:Int(opt+(t==nt)*(O%nt))
            o=Int(i+(t-1)*opt)
            if clean[o]
                obs_lhs[t][i]=cache[o]
            else
                obsl = obs_lengths[o]
                revcomp ? log_motif_expectation = log(0.5/obsl) : log_motif_expectation = log(1/obsl)#log_motif_expectation-nMica has 0.5 per base for including the reverse complement, 1 otherwise
                mixview=view(mix,o,:)
                mixwmls=source_wmls[mixview]
                score_mat=score_obs_sources(sources[mixview], observations[1:obsl,o], obsl, mixwmls, revcomp=revcomp)
                obs_source_indices = findall(mixview)
                obs_cardinality = length(obs_source_indices) #the more sources, the greater the cardinality_penalty
                if obs_cardinality > 0
                    penalty_sum = sum(exp.(fill(log_motif_expectation,obs_cardinality)))
                    penalty_sum > 1. && (penalty_sum=1.)
                    cardinality_penalty=log(1.0-penalty_sum)
                else
                    cardinality_penalty=0.0
                end

                obs_lhs[t][i]=weave_scores(obsl, view(bg_scores,:,o), score_mat, obs_source_indices, mixwmls, log_motif_expectation, cardinality_penalty, revcomp)
            end
        end
    end

    returncache ? (return lps([lps(obs_lhs[t]) for t in 1:nt]), vcat(obs_lhs...)) : (return lps([lps(obs_lhs[t]) for t in 1:nt]))
end

                function score_obs_sources(sources::Vector{Tuple{Matrix{AbstractFloat},Integer}}, observation::Vector{Integer}, obsl::Integer, source_wmls::Vector{Integer}; revcomp=true)
                    scores=Vector{Matrix{AbstractFloat}}(undef, length(sources))

                    for (s,source) in enumerate(sources)
                        pwm = source[1] #get the PWM from the source tuple
                        wml = source_wmls[s] #weight matrix length
                        source_stop=obsl-wml+1 #stop scannng th source across the observation here

                        scores[s]=nnlearn.score_source(observation, pwm, source_stop, revcomp) #get the scores for this oxs
                    end

                    return scores
                end

                function score_source(observation::AbstractArray{Integer,1}, source::Matrix{AbstractFloat}, source_stop::Integer, revcomp::Bool=true)
                    revcomp ? (revsource = revcomp_pwm(source); score_matrix = zeros(source_stop,2)) : score_matrix = zeros(source_stop)
                    forward_score = 0.0
                    revcomp && (reverse_score = 0.0)

                    for t in 1:source_stop
                        forward_score = 0.0 #initialise scores as log(p=1)
                        revcomp && (reverse_score = 0.0)

                        for position in 1:size(source)[1]
                            score_loc = t+position-1 #score_loc is the position of the obs to be scored by PWM
                            forward_score += source[position,observation[score_loc]] #add the appropriate log PWM value from the source to the score
                            revcomp && (reverse_score += revsource[position,observation[score_loc]])
                        end

                        score_matrix[t,1] = forward_score #assign scores to the matrix
                        revcomp && (score_matrix[t,2] = reverse_score)
                    end

                    return score_matrix
                end
                                function revcomp_pwm(pwm::Matrix{AbstractFloat}) #in order to find a motif on the reverse strand, we scan the forward strand with the reverse complement of the pwm, reordered 3' to 5', so that eg. an PWM for an ATG motif would become one for a CAT motif
                                    return pwm[end:-1:1,end:-1:1]
                                end

                function weave_scores(obsl::Integer, bg_scores::SubArray, score_mat::Vector{Matrix{AbstractFloat}}, obs_source_indices::Vector{Integer}, source_wmls::Vector{Integer}, log_motif_expectation::AbstractFloat, cardinality_penalty::AbstractFloat,  revcomp::Bool=true)
                    L=obsl+1
                    lh_vec = zeros(L)#likelihood vector is one position (0 initialiser) longer than the observation
                    osi_emitting = Vector{Integer}()

                    @inbounds for i in 2:L #i=1 is ithe lh_vec initializing 0, i=2 is the score of the first background position (ie t=1)
                        t=i-1
                        score = lps(lh_vec[i-1], bg_scores[t], cardinality_penalty)

                        #logic: all observations are scored from t=wml to the end of the obs-therefore check at each position for new sources to add (indexed by vector position to retrieve source wml and score matrix)
                        if length(osi_emitting)<length(obs_source_indices)
                            for n in 1:length(obs_source_indices)
                                if !(n in osi_emitting)
                                    t>= source_wmls[n] && (push!(osi_emitting,n))
                                end
                            end
                        end

                        for n in osi_emitting
                            wml = source_wmls[n]
                            from_score = lh_vec[i-wml+1] #score at the first position of the PWM
                            score_array = score_mat[n] #get the source score matrix
                            score_idx = t - wml + 1 #translate t to score_array idx for emission score
                            emit_score = score_array[score_idx,1] #emission score at the last position of the PWM

                            score = logaddexp(score, lps(from_score, emit_score, log_motif_expectation))
                            if revcomp #repeat each source contribution with the revcomp score vector if required
                                emit_score = score_array[score_idx,2]
                                score = logaddexp(score, lps(from_score, emit_score, log_motif_expectation))
                            end
                        end
                        lh_vec[i] = score
                    end
                    return lh_vec[end]
                end


text:
using MarkImg

img = rand(100, 100)
imgc = convert.(RGB, img)
x, y, r = 10, 15, 5 #row, column, radius
markImg!(imgc, x, y, r)


text:
using Test
using SymDiff

@test simplify(:(a * b * 0)) == 0
@test simplify(:(a * 1 * b)) == :(a * b)


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # SIMD: The parallelism that can (sometimes) happen automatically
#
# SIMD: Single-instruction, multiple data
#
# (Also confusingly called vectorization)

#-

# ## The architecture
#
# Instead of computing four sums sequentially:
#
# \begin{align}
# x_1 + y_1 &\rightarrow z_1 \\
# x_2 + y_2 &\rightarrow z_2 \\
# x_3 + y_3 &\rightarrow z_3 \\
# x_4 + y_4 &\rightarrow z_4
# \end{align}
#
# Modern processors have vector processing units that can do it all at once:
#
# $$
# \left(\begin{array}{cc}
# x_1 \\
# x_2 \\
# x_3 \\
# x_4
# \end{array}\right)
# +
# \left(\begin{array}{cc}
# y_1 \\
# y_2 \\
# y_3 \\
# y_4
# \end{array}\right)
# \rightarrow
# \left(\begin{array}{cc}
# z_1 \\
# z_2 \\
# z_3 \\
# z_4
# \end{array}\right)
# $$

#-

# ## Making it happen

#-

# Simple task: compute the sum of a vector:

A = rand(100_000)
function simplesum(A)
    result = zero(eltype(A))
    for i in eachindex(A)
        @inbounds result += A[i]
    end
    return result
end

simplesum(A)

#-

using BenchmarkTools
@btime simplesum($A)

# So, is that good?

@btime sum($A)

# We're slower that the builtin `sum` — and we're getting a different answer, too! Let's look at what happens with a 32-bit float instead of a 64 bit one. Each element has half the number of bits, so lets also double the length (so the total number of bits processed remains constant).

A32 = rand(Float32, length(A)*2)
@btime simplesum($A32)
@btime sum($A32);

# That's even worse! What's going on here?  We're seeing an even multiple number
# difference in our performance — perhaps Julia's builtin sum is using some
# parallelism? Let's try using SIMD ourselves:

function simdsum(A)
    result = zero(eltype(A))
    @simd for i in eachindex(A)
        @inbounds result += A[i]
    end
    return result
end
@btime simdsum($A)
@btime simdsum($A32)

# What did that do and why don't we always use `@simd for` — or why doesn't Julia
# just always use `@simd` for every `for` loop automatically?  Look at the values:

simplesum(A), simdsum(A), sum(A)

#-

simplesum(A32), simdsum(A32), sum(A32)

# Why aren't they the same?
#
# Without `@simd`, Julia is doing _exactly_ what we told it to do: it's taking
# each element of our array and adding it to a big pile sequentially. Our answer
# is smaller than what Julia's builtin `sum` thinks it is: that's because as our
# pile gets bigger we begin losing the lower bits of each element that we're
# adding, and those small losses begin to add up!
#
# The `@simd` macro tells Julia that it can re-arrange floating point additions —
# even if it would change the answer. Depending on your CPU, this may lead to 2x or 4x
# or even 8x parallelism. Essentially, Julia is computing independent sums for
# the even indices and the odd indices simultaneously:
#
# \begin{align}
# odds &\leftarrow 0 \\
# evens &\leftarrow 0 \\
# \text{loop}&\ \text{odd}\ i: \\
#     &\left(\begin{array}{cc}
# odds \\
# evens
# \end{array}\right)
# \leftarrow
# \left(\begin{array}{cc}
# odds \\
# evens
# \end{array}\right)
# +
# \left(\begin{array}{cc}
# x_{i} \\
# x_{i+1}
# \end{array}\right) \\
# total &\leftarrow evens + odds
# \end{align}
#
# In many cases, Julia can and does know that a for-loop can be SIMD-ed and it
# will take advantage of this by default!

B = rand(1:10, 100_000)
@btime simplesum($B)
@btime sum($B)
B32 = rand(Int32(1):Int32(10), length(B)*2)
@btime simplesum($B32)
@btime simdsum($B32)

# How can we see if something is getting vectorized?

@code_llvm simdsum(A32)

# So what are the challenges?
#
# * Biggest hurdle is that you have to convince Julia and LLVM that it's able to
#   use SIMD instructions for your given algorithm. That's not always possible.
# * There are lots of limitations of what can and cannot be SIMD-ed:

@doc @simd

# * You do need to think through the consequences of re-ordering your algorithm.

#-

# ## A slightly trickier case

using BenchmarkTools

#-

function diff!(A, B)
    A[1] = B[1]
    for i in 2:length(A)
        @inbounds A[i] = B[i] - B[i-1]
    end
    return A
end
A = zeros(Float32, 100_000)
B = rand(Float32, 100_000)

diff!(A, B)
[B[1];diff(B)] == A

#-

@btime diff!($A, $B)
@btime diff($B);

# But what happens if we do it in-place?

Bcopy = copy(B)
@btime diff!($Bcopy, $Bcopy);

# What happened?

@code_llvm diff!(A, B)

# We can manually assert that arrays don't alias (or have any loop-dependencies),
# with the very special `@simd ivdep` flag, but this can be disastrous:

function unsafe_diff!(A, B)
    A[1] = B[1]
    @simd ivdep for i in 2:length(A)
        @inbounds A[i] = B[i] - B[i-1]
    end
    return A
end
@btime unsafe_diff!($A, $B)
[B[1];diff(B)] == A
Bcopy = copy(B)
unsafe_diff!(Bcopy, Bcopy)
[B[1];diff(B)] == Bcopy

# If you really want to get your hands dirty, you can use the [SIMD.jl](https://github.com/eschnett/SIMD.jl)
# package to manually specify those `<8 x float>` things that LLVM generates.
# BUT: this is tricky and a pain; often it's just to be aware of what makes
# Julia code automatically SIMD-able, some of the cases where it may fail, and
# how to check its work.

#-

# ## SIMD
#
# * Exploits built-in parallelism in a processor
# * Best for small, tight innermost loops
# * Often happens automatically if you're careful
#     * Follow the [perforance best practices](https://docs.julialang.org/en/v1/manual/performance-tips/)
#     * `@inbounds` any array acesses
#     * No branches or (non-inlined) function calls
# * Can use `@simd` to allow Julia to break some rules to make it happen
#     * But be careful, especially with `@simd ivdep`!
# * Depending on processor and types involved, can yield 2-16x gains with extraordinarily little overhead
#     * Smaller datatypes can improve this further; use `Float32` instead of `Float64`
#       if possible, `Int32` instead of `Int64`, etc.
#     * When buying a new processor, look for [AVX-512](https://en.wikichip.org/wiki/x86/avx-512) support

max_stars_repo_path: Courses/Parallel_Computing/030 SIMD.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_path: Courses/Parallel_Computing/030 SIMD.jl

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_path: Courses/Parallel_Computing/030 SIMD.jl

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "cargo_license"
version = v"0.4.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/onur/cargo-license.git", "f489e18b0fc45ab31d9766841e7f551e96ee6a8c"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cargo-license
cargo build --release --all
mkdir ${bindir}
cp target/${rust_target}/release/cargo-license${exeext} ${bindir}/
install_license $WORKSPACE/srcdir/cargo-license/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
# Rust toolchain for i686 Windows is unusable
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("cargo-license", :cargo_license),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
                compilers=[:c, :rust], julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Bounds

using Statistics
using LinearAlgebra

export SimpleBound

### Silverman Rule of Thumb for bandwidth selection ##
function silverman(x::Vector{T}) where T<:Real
    # Silverman rule of thumb
    return 1.06*length(x)^(-0.2)*std(x)
end

### Epanechnikov Kernel function ###
function epa(x::T,x0::T=0.0,h::T=1.0) where T<:Real
    temp = min(abs((x-x0)/h),1.0)
    return 0.75*(1.0-temp^2)
end

function epa(x::Vector{T},x0::T=0.0,h::T=1.0) where T<:Real
    temp = min.(norm.((x.-x0)./h),1.0)
    return 0.75*(1.0.-temp.^2)
end

### Kernel density estimation ###
function kdens(X::Vector{T},x0::T,h::T) where T<:Real
    f = sum(epa.((X .- x0) ./ h))
    return f/(length(X)*h)
end

### Kernel regression estimation ###
function kreg(Y::Vector{T},X::Vector{T},x0::T,h::T) where T<:Real
    f = 0.0
    m = 0.0
	nx=length(X)
	ny=length(Y)
	if nx != ny
		error("length of Y and X do not match")
	else
		df = epa.((X .- x0) ./ h)
		f = sum(df)
		m = sum(Y.*df)
	end
	if f >0
        return m/f
    else
        return 0.0
    end
end

function kreg(Y::Vector{T},X::Vector{T},x0::T,h::T,w::Vector{T}) where T<:AbstractFloat
    # This is a version of kreg where we can introduce weights.
    # It is done mostly so we can conditon on binary weights, i.e. exclude some observations

    f = 0.0
    m = 0.0
	nx=length(X)
	ny=length(Y)
    nw=length(w)

    if nx != ny || nx !=nw
		error("length of Y and X do not match")
	else
		df = w.*epa.((X .- x0) ./ h)
		f = sum(df)
		m = sum(w.*Y.*df)
	end
	if f >0
        return m/f
    else
        return 0.0
    end
end

### Worst Case Scenario Bounds
mutable struct Treatment{T}
    prob0 :: T
    prob1 :: T
    yhat0 :: T
    yhat1 :: T
    bound0L :: T
    bound0U :: T
    bound1L :: T
    bount1U :: T
    treatL :: T
    treatU :: T
end

function treatment(y::Vector{T},z::Vector{T},x::Vector{T},x0::T,cont::Bool,h::T=1.0) where T<:Real
    ny=length(y)
    nx=length(x)
    nz=length(z)
    if (ny !=nx) || (ny != nz)
        error("vector length not matching")
    end
    zer0 ::T = 0
    res::Vector{T} = Treatment(zer0,0,0,0,0,0,0,0,0,0)
    if cont
        res.prob1 = kreg(z,x,x0,h)
        res.prob0 = 1- res.prob1
        res.yhat1 = kreg(y,x,x0,h,z)
        res.yhat0 = kreg(y,x,x0,h,1 .-z)
    else
        nx0 = count(x .== x0)
        temp1 = z .* (x .== x0)
        nz1 = count(temp1)

        res.prob1 = nz1/nx0
        res.prob0 = 1 - res.prob1
        res.yhat1 = count(y .* z .* (x .== x0))/nz1
        res.yhat0 = count(y .* (1 .-z) .* (x .== x0))/nz1
    end

end

struct SimpleBound
	LB :: Float64
	UB :: Float64
	method :: String
end

end # module


text:
using Test

using Agtor, Dates, CSV

@testset "Loading YAML spec" begin
    data_dir::String="test/data/"
    pump_spec_dir::String = "$(data_dir)pumps/"
    pump_specs::Dict{String, Dict} = load_yaml(pump_spec_dir)

    @test pump_specs["groundwater"]["capital_cost"][1] == "RealParameter"
end

@testset "Updating an AgParameter" begin

    test = Agtor.RealParameter("Test", 0.0, 1.0, 0.5)

    @test test.value == 0.5
    set_params!(test, 0.75)
    @test test.value == 0.75

    data_dir = "test/data/"

    z1 = setup_zone(data_dir)
    tc_z1 = deepcopy(z1)  # test copy

    samples = DataFrame!(CSV.File("test/data/scenarios/sampled_params.csv"))

    @test z1.fields[1].irrigation.efficiency.value == tc_z1.fields[1].irrigation.efficiency.value
    update_model!(z1, samples[1, :])
    @test z1.fields[1].irrigation.efficiency.value != tc_z1.fields[1].irrigation.efficiency.value
end

@testset "No duplicate parameters" begin
    zone = setup_zone()

    all_params = collect_agparams!(zone)
    param_names = all_params[:name]
    uni_params = unique(param_names)

    @test length(param_names) == length(uni_params)
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
using Literate

# Copy tutorial sources
for t in readdir("$(@__DIR__)")
    if occursin(r".*_tut.jl", t)
        println("\nTutorial $t in $(@__DIR__)\n")
        Literate.markdown(joinpath(@__DIR__, t), "$(@__DIR__)"; documenter=false);
        # cp(t, "../../../src/" * t, force = true)
    end
end

# Copy ancillary files
# for a in ["fast_top_ref.txt"]
#     println("\nAncillary $a in $(pwd())\n")
#     cp(a, "../../../src/" * a, force = true)
# end


text:
#=

  Huey, Dewey and Louie problem in ConstraintSolver.jl

  From Marriott & Stuckey, Programming with Constraints, page 42
  """
  Huey, Dewey and Louie are being questioned by their uncle. These are the
  statements the make:
   Huey: Dewey and Louie has equal share in it; if one is quitly, so
         is the other.
   Dewey: If Huey is guilty, then so am I.
   Louie: Dewey and I are not both quilty.

  Their uncle, knowing that they are cub scouts, realises that they
  cannot tell a lie. Has he got sufficient information to decide who
  (if any) are quilty?
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function huey_dewey_louie(print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 3
    @variable(model, x[1:n], Bin)
    Huey,Dewey,Louie = x

    #  Huey: Dewey and Louie has equal share in it; if one is
    # quitly, so is the other.
    # Dewey == 1 <=> Louie == 1
    @variable(model, b1, Bin)
    @constraint(model, b1 := {Dewey ==1})
    @constraint(model, b1 := {Louie ==1})

    # Dewey: If Huey is guilty, then so am I.
    # Huey ==1 => Dewey ==1
    @variable(model, b2, Bin)
    @constraint(model, b2 := {Huey ==1})
    @constraint(model, b2 => {Dewey ==1})

    #  Louie: Dewey and I are not both quilty.
    # ~(Dewey ==1  /\ Louie ==1)
    @variable(model, b3, Bin)
    @variable(model, b4, Bin)
    @constraint(model, b3 := {Dewey == 1})
    @constraint(model, b4 := {Louie == 1})
    @constraint(model, b3+b4 < 2)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time huey_dewey_louie()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#
# test the linear program solver on a simple problem
#

include("../src/homogeneous_algorithm/homogeneous_algorithm.jl")
include_print("testing_tools.jl")
using Ipopt
#simple_tests();

# 54
# 59
#

function test_problem(name::String)
	file_name = "Problems/" * name * ".mat";
	A, b, c = get_netlib_problem(file_name);
	#Q = -speye(length(c));
	Q = tridiagonal(length(c),0.0,1.0);
	#Q = spzeros(length(c),length(c));
	println("Solving ", file_name, " with the homogeneous algorithm")
	println(size(A,2), " variables and ", size(A,1), " constraints")
	println("Non-zeros: ", length(nonzeros(A)))
	settings = class_settings();

	println("=================== Linear system solver is matlab ldl ==================")
	settings.linear_system_solver = linear_solver_MATLAB();
	settings.linear_system_solver.options.sym = 2;

	settings.newton_solver = class_newton_hsd();
	trivial_test(A, b, c, Q, 1, file_name, settings, true);

	#settings.newton_solver = class_newton_ip();
	#trivial_test(A, b, c, Q, 1, file_name, settings, true);

	println("=========================================================================")
	println("Calls IPOPT")
	solve_with_JuMP(A, b, c, Q, IpoptSolver(max_iter=1000));
end

#file_list = ["Problems/AFIRO.mat", "Problems/BANDM.mat", "Problems/BEACONFD.mat", "Problems/BLEND.mat", "Problems/BOEING2.mat", "Problems/BRANDY.mat", "Problems/CAPRI.mat", "Problems/BORED_3D.mat", "Problems/STANDGUB.mat"];

#main();
println("For example, type: test_problem(\"BORED_3D\")")

max_forks_repo_name: ohinder/Everest


text:
using GalacticMetaheuristics, GalacticOptim
using Test

@testset "GalacticMetaheuristics.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    optprob = OptimizationFunction(rosenbrock)
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[1.5, 1.5])
    sol = solve(prob, ECA())
    @test 10 * sol.minimum < l1

    sol = solve(prob, Metaheuristics.DE())
    @test 10 * sol.minimum < l1

    sol = solve(prob, PSO())
    @test 10 * sol.minimum < l1

    sol = solve(prob, ABC())
    @test 10 * sol.minimum < l1

    sol = solve(prob, CGSA(N=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, SA())
    @test 10 * sol.minimum < l1

    sol = solve(prob, WOA())
    @test 10 * sol.minimum < l1

    sol = solve(prob, ECA())
    @test 10 * sol.minimum < l1

    sol = solve(prob, Metaheuristics.DE(), use_initial=true)
    @test 10 * sol.minimum < l1

    sol = solve(prob, PSO(), use_initial=true)
    @test 10 * sol.minimum < l1

    sol = solve(prob, ABC(), use_initial=true)
    @test 10 * sol.minimum < l1

    sol = solve(prob, CGSA(N=100), use_initial=true)
    @test 10 * sol.minimum < l1

    sol = solve(prob, SA(), use_initial=true)
    @test 10 * sol.minimum < l1

    sol = solve(prob, WOA(), use_initial=true)
    @test 10 * sol.minimum < l1
end


text:
# Autogenerated wrapper script for licensecheck_jll for x86_64-linux-gnu
export licensecheck

JLLWrappers.@generate_wrapper_header("licensecheck")
JLLWrappers.@declare_library_product(licensecheck, "licensecheck.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        licensecheck,
        "lib/licensecheck.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
############ Tickers  ####################
"""
tickers(opts::PolyOpts, search::AbstractString;
        active=true, sort="ticker", order="asc", limit=10, kwargs...)

Query all ticker symbols which are supported by Polygon.io. This API currently includes Stocks/Equities, Crypto, and Forex.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * search::AbstractString - Search for terms within the ticker and/or company name.
 * active::Boolean - Specify if the tickers returned should be actively traded on the queried date. Default is true.
 * sort::String - The field to sort the results on. Default is ticker. If the search query parameter is present, sort is ignored and results are ordered by relevance.
 * order::String - The order to sort the results on. Default is asc (ascending).
 * limit::Integer - Limit the size of the response, default is 100 and max is 1000. If your query returns more than the max limit and you want to retrieve the next page of results, see the next_url response attribute.
 * kwargs::Dict - Additional query parameters.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> tickers(opts, "bitcoin")
```

# Returns
* A JSON3.Array or specified tabular representation of the JSON3.Array.
* See https://polygon.io/docs/get_v3_reference_tickers_anchor for documentation on response attributes and supported keyword arguments.
"""
function tickers(opts::PolyOpts, search::AbstractString;
                active=true, sort="ticker", order="asc", limit=10, kwargs...)

    params = Dict(
        "search" => search,
        "active" => active,
        "sort" => sort,
        "order" => order,
        "limit" => limit,
        "apiKey" => opts.api_key
    )
    # Extract kwargs and add to params
    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), tickers_base_url, params, opts.sink)

end

############ Ticker Types  ####################
"""
    ticker_types(opts::PolyOpts)

Get a mapping of ticker types to their descriptive names.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> ticker_types(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_types_anchor for documentation on response attributes and supported keyword arguments.
"""
function ticker_types(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(NoSinkYesResults(), ticker_types_base_url, params, opts.sink)
end

############ Tickers Details ####################
"""
    ticker_details(opts::PolyOpts, stocksTicker::AbstractString)

Get details for a ticker symbol's company/entity.
This provides a general overview of the entity with information such as name, sector, exchange, logo and similar companies.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString - The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> ticker_details(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_ticker_details_anchor for documentation on response attributes and supported keyword arguments.
"""
function ticker_details(opts::PolyOpts, stocksTicker::AbstractString)
    params = Dict("apiKey" => opts.api_key)
    ticker_details_url = "$ticker_details_base_url/$stocksTicker/company"

    return generate_output_from_url(YesSinkNoResults(), ticker_details_url, params, opts.sink)
end

############ Ticker Details vX ####################
"""
    ticker_details_vX(opts::PolyOpts, ticker::AbstractString, date::AbstractString)

Get a single ticker supported by Polygon.io.
This response will have detailed information about the ticker and the company behind it.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * ticker::AbstractString - The ticker symbol of the asset.
 * date::AbstractString - Specify a point in time to get information about the ticker available on that date. When retrieving information from SEC filings, we compare this date with the period of report date on the SEC filing.
For example, consider an SEC filing submitted by AAPL on 2019-07-31, with a period of report date ending on 2019-06-29. That means that the filing was submitted on 2019-07-31, but the filing was created based on information from 2019-06-29. If you were to query for AAPL details on 2019-06-29, the ticker details would include information from the SEC filing.
Defaults to the most recent available date.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> ticker_details_vX(opts, "AAPL", "2017-01-01")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_vX_reference_tickers__ticker__anchor for documentation on response attributes and supported keyword arguments.
"""
function ticker_details_vX(opts::PolyOpts, ticker::AbstractString, date::AbstractString)
    # TODO: Dispatch on proper Date type?
    ticker_details_vX_url = "$ticker_details_vX_base_url/$ticker"

    params = Dict(
        "apiKey" => opts.api_key,
        "date"   => date
    )

    return generate_output_from_url(YesSinkYesResults(), ticker_details_vX_url, params, opts.sink)
end

############ Ticker News  #######################
"""
    ticker_news(opts::PolyOpts, ticker::AbstractString;
                published_utc_gte="2021-04-26", limit=10,
                order="descending", sort="published_utc",kwargs...)

Get the most recent news articles relating to a stock ticker symbol, including a summary of the article and a link to the original source.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * ticker::AbstractString - Ticker symbol. Return results where this field equals the value.
 * published_utc_gte::AbstractString - Return results where this field is greater than or equal to the value.
 * limit::Integer - Limit the size of the response, default is 100 and max is 1000.
    If your query returns more than the max limit and you want to retrieve the next page of results, see the next_url response attribute.
 * order::String - Order the results in ascending or descending order.
 * sort::String - The field key to sort the results on.
 * kwargs::Dict - Additional query parameters

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> ticker_news(opts, "AAPL", limit=5)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_news_anchor for documentation on response attributes and supported keyword arguments.
"""
function ticker_news(opts::PolyOpts, ticker::AbstractString;
                    published_utc_gte="2021-04-26", limit=10,
                    order="descending", sort="published_utc",kwargs...)

    params = Dict(
        "apiKey" => opts.api_key,
        "ticker" => ticker,
        "published_utc.gte" => published_utc_gte,
        "limit" => limit,
        "order" => order,
        "sort" => sort
    )
    # Extract kwargs and add to params
    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), ticker_news_base_url, params, opts.sink)

end

############ Markets  ####################
"""
    markets(opts::PolyOpts)

Get a list of markets that are currently supported by Polygon.io.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> markets(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_markets_anchor for documentation on response attributes and supported keyword arguments.

"""
function markets(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkYesResults(), markets_base_url, params, opts.sink)
end

############ Locales  ####################
"""
    locales(opts::PolyOpts)

Get a list of locales currently supported by Polygon.io.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> locales(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_locales_anchor for documentation on response attributes and supported keyword arguments.
"""
function locales(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkYesResults(), locales_base_url, params, opts.sink)
end

############ Stock Splits  ####################
"""
    stock_splits(opts::PolyOpts, stocksTicker::AbstractString)

Get a list of historical stock splits for a ticker symbol, including the execution and payment dates of the stock split, and the split ratio.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString - The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stock_splits(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_splits__stocksTicker__anchor for documentation on response attributes and supported keyword arguments.
"""
function stock_splits(opts::PolyOpts, stocksTicker::AbstractString)
    stock_splits_url = "$stock_splits_base_url/$stocksTicker"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(YesSinkYesResults(), stock_splits_url, params, opts.sink)
end

############ Stock Dividends  ####################
"""
    stock_dividends(opts::PolyOpts, stocksTicker::AbstractString)

Get a list of historical dividends for a stock, including the relevant dates and the amount of the dividend.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString - The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stock_dividends(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_dividends__stocksTicker__anchor for documentation on response attributes and supported keyword arguments.
"""
function stock_dividends(opts::PolyOpts, stocksTicker::AbstractString)
    stock_dividends_url = "$stock_dividends_base_url/$stocksTicker"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(YesSinkYesResults(), stock_dividends_url, params, opts.sink)
end

############ Stock Financials  ####################
"""
    stock_financials(opts::PolyOpts, stocksTicker::AbstractString; limit=5, kwargs...)

Get historical financial data for a stock ticker.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString - The ticker symbol of the stock/equity.
 * limit::Integer - Limit the number of results.
 * kwargs::Any: A list of additional arguments to pass to the Polygon IO API.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stock_financials(opts, "AAPL", limit=5)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_reference_financials_anchor for documentation on response attributes and supported keyword arguments.
"""
function stock_financials(opts::PolyOpts, stocksTicker::AbstractString; limit=5, kwargs...)
    stock_financials_url = "$stock_financials_base_url/$stocksTicker"
    params = Dict(
        "apiKey" => opts.api_key,
        "limit" => limit
    )
    # Extract kwargs and add to params
    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), stock_financials_url, params, opts.sink)
end

############ Market Holidays  ####################
"""
    market_holidays(opts::PolyOpts)

Get upcoming market holidays and their open/close times.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> market_holidays(opts)
```
# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v1_marketstatus_upcoming_anchor for documentation on response attributes and supported keyword arguments.
"""
function market_holidays(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkNoResults(), market_holidays_base_url, params, opts.sink)
end

############ Market Status  ####################
"""
    market_status(opts::PolyOpts)

Get the current trading status of the exchanges and overall financial markets.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> market_status(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v1_marketstatus_now_anchor for documentation on response attributes and supported keyword arguments.
"""
function market_status(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkNoResults(), market_status_base_url, params, opts.sink)
end

############ Stock Exchanges  ####################
"""
    stock_exchanges(opts::PolyOpts)

Get a list of stock exchanges which are supported by Polygon.io.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stock_exchanges(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v1_meta_exchanges_anchor for documentation on response attributes and supported keyword arguments.
"""
function stock_exchanges(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkNoResults(), stock_exchanges_base_url, params, opts.sink)
end

############ Condition Mappings  ####################
"""
    condition_mappings(opts::PolyOpts, tickertype="trades")

Get a unified numerical mapping for conditions on trades and quotes.
Each feed/exchange uses its own set of codes to identify conditions, so the same condition may have a different code depending on the originator of the data.
Polygon.io defines its own mapping to allow for uniformly identifying a condition across feeds/exchanges.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.
 * tickertype::AbstractString - The type of ticks to return mappings for. Must be one of "trades" or "quotes".

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> condition_mappings(opts, "trades")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v1_condition_mappings_anchor for documentation on response attributes and supported keyword arguments.
"""
function condition_mappings(opts::PolyOpts, tickertype="trades")
    condition_mappings_url = "$condition_mappings_base_url/$tickertype"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(YesSinkNoResults(), condition_mappings_url, params, opts.sink)
end

############ Crypto Exchanges ####################
"""
    crypto_exchanges(opts::PolyOpts)

Get a list of cryptocurrency exchanges which are supported by Polygon.io.

# Arguments
 * opts::PolyOpts - The PolyOpts object used to configure the request.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> crypto_exchanges(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v1_meta_crypto-exchanges_anchor for documentation on response attributes and supported keyword arguments.
"""
function crypto_exchanges(opts::PolyOpts)
    params = Dict("apiKey" => opts.api_key)
    return generate_output_from_url(YesSinkNoResults(), crypto_exchanges_base_url, params, opts.sink)
end


text:
__precompile__()

module ospa

using Hungarian
using LinearAlgebra
using Statistics
using Random

export ospa_kmeans, ospa_barycenter,
       ospa_dist, Pointcloud,
       generateData
       #mean_initialization, sample_from_mean_initialization,
       #ospa_barycenter_sol2

include("types.jl")
include("distance.jl")
include("barycenter.jl")
include("mean_barycenter_initialization.jl")
include("recursive_barycenter_initialization.jl")
include("cluster_initialization.jl")
include("kmeans.jl")
include("testdata.jl")

end # module


text:
const USE_GPU = false  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using Plots, Printf, Statistics

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, P::Data.Array, dt::Data.Number, ρ::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @inn(Vx) = @inn(Vx) - dt/ρ*@d_xi(P)/dx
    @inn(Vy) = @inn(Vy) - dt/ρ*@d_yi(P)/dy
    @inn(Vz) = @inn(Vz) - dt/ρ*@d_zi(P)/dz
    return
end

@parallel function compute_P!(P::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, dt::Data.Number, k::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(P) = @all(P) - dt*k*(@d_xa(Vx)/dx + @d_ya(Vy)/dy + @d_za(Vz)/dz)
    return
end

##################################################
@views function acoustic3D()
    # Physics
    lx, ly, lz = 40.0, 40.0, 40.0  # domain extends
    k          = 1.0               # bulk modulus
    ρ          = 1.0               # density
    t          = 0.0               # physical time
    # Numerics
    nx, ny, nz = 255, 255, 255     # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    nt         = 1000              # number of timesteps
    nout       = 10                # plotting frequency
    # Derived numerics
    dx, dy, dz = lx/(nx-1), ly/(ny-1), lz/(nz-1) # cell sizes
    # Array allocations
    P          = @zeros(nx  ,ny  ,nz  )
    Vx         = @zeros(nx+1,ny  ,nz  )
    Vy         = @zeros(nx  ,ny+1,nz  )
    Vz         = @zeros(nx  ,ny  ,nz+1)
    # Initial conditions
    P         .= Data.Array([exp(-((ix-1)*dx-0.5*lx)^2 -((iy-1)*dy-0.5*ly)^2 -((iz-1)*dz-0.5*lz)^2) for ix=1:size(P,1), iy=1:size(P,2), iz=1:size(P,3)])
    dt         = min(dx,dy,dz)/sqrt(k/ρ)/6.1
    # Preparation of visualisation
    ENV["GKSwstype"]="nul"; if isdir("viz3D_out")==false mkdir("viz3D_out") end; loadpath = "./viz3D_out/"; anim = Animation(loadpath,String[])
    println("Animation directory: $(anim.dir)")
    y_sl       = Int(ceil(ny/2))
    X, Y, Z    = -lx/2:dx:lx/2, -ly/2:dy:ly/2, -lz/2:dz:lz/2
    # Time loop
    for it = 1:nt
        if (it==11)  global wtime0 = Base.time()  end
        @parallel compute_V!(Vx, Vy, Vz, P, dt, ρ, dx, dy, dz)
        @parallel compute_P!(P, Vx, Vy, Vz, dt, k, dx, dy, dz)
        t = t + dt
        # Visualisation
        if mod(it,nout)==0
            heatmap(X, Z, Array(P)[:,y_sl,:]', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Z[1],Z[end]), c=:viridis, title="Pressure"); frame(anim)
        end
    end
    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = (4*2)/1e9*nx*ny*nz*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    wtime_it = wtime/(nt-10)                           # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                          # Effective memory throughput [GB/s]
    @printf("Total steps=%d, time=%1.3e sec (@ T_eff = %1.2f GB/s) \n", nt, wtime, round(T_eff, sigdigits=2))
    gif(anim, "acoustic3D.gif", fps = 15)
    return
end

acoustic3D()


text:
using PETLION
using Test

"""
The tests work properly if they are run in a terminal, but running them in the @testset environment fails.
They are blanked out for now.
"""

@testset "PETLION.jl" begin
    #= SAVE_SYMBOLIC_FUNCTIONS = PETLION.options[:SAVE_SYMBOLIC_FUNCTIONS]
    PETLION.options[:SAVE_SYMBOLIC_FUNCTIONS] = false

    opts = (
        temperature = false,
        N_p=10,
        N_n=10,
        N_s=10,
    )

    p    = petlion(LCO; jacobian=:symbolic, opts...)
    p_AD = petlion(LCO; jacobian=:AD, opts...)

    # AD matches symbolic
    @test isapprox([simulate(model, 0:100:3600, I=-1,  SOC=1, outputs=(:t, :V)).V for model in (p, p_AD)]...)
    @test isapprox([simulate(model, 0:100, P=-10, SOC=1).V for model in (p, p_AD)]...)
    @test isapprox([simulate(model, 0:10,  V=3.5, SOC=1).V for model in (p, p_AD)]...)

    # all outputs work
    @test !isempty(simulate(p, I=-1, SOC=1, outputs=:all))

    # functions are working
    p.opts.SOC=0
    @test simulate(p, 1,I=1).V[end] ≠ simulate(p, 1,I=cos).V[end]
    P_fun = t -> 100cos(t)
    @test simulate(p, 1,P=100).V[end] ≠ simulate(p, 1,P=P_fun).V[end]

    # function matches CC
    p.opts.SOC=0
    @test all(simulate(p, 0:1000, I=1).V .=== simulate(p, 0:1000, I=(t)->1).V)

    # :hold and I_max stop conditions are working
    sol = simulate(p, 100, I=-0.1, SOC=1, outputs=(:t, :V, :I, :P))
    simulate!(sol, p, 100, I=-0.1)
    simulate!(sol, p, 100, V=:hold)
    simulate!(sol, p, 100, V=:hold, I_max=-0.05)
    @test sol.I[end] ≈ -0.05
    simulate!(sol, p, P=:hold)
    @test sol[end-1].P[end] ≈ sol[end].P[1]

    PETLION.options[:SAVE_SYMBOLIC_FUNCTIONS] = SAVE_SYMBOLIC_FUNCTIONS
    =#
end


text:
using HPAT
using MPI

@acc hpat function calcKmeans(k::Int64, file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points",file_name)
    clusters = Kmeans(points, k, 10)
    return clusters
end

K = 3

if (length(ARGS) > 0)
	K = parse(Int, ARGS[1])
end

c2 = calcKmeans(K, HPAT.getDefaultDataPath()*"/kmeans_large.hdf5")

if MPI.Comm_rank(MPI.COMM_WORLD)==0	println("Centroids: ", c2) end