text:
import TranscodingStreams
import Dates

function compressdat(vtuheader,dat,headertype,level=5)
	compr_dat = Vector{UInt8}[]
	bytedata = reinterpret(UInt8,dat)
	for i = 1:length(vtuheader.compressed_blocksizes)-1
		_start = (i-1)*vtuheader.blocksize+1
		_stop = i*vtuheader.blocksize
		tmp = bytedata[_start:_stop]
		buf = IOBuffer(append=true)
		zWriter = ZlibCompressorStream(buf,level=level)
		write(zWriter, tmp)
		write(zWriter, TranscodingStreams.TOKEN_END)
		flush(zWriter)
		push!(compr_dat,read(buf))
	end
	_start = (vtuheader.num_blocks-1)*vtuheader.blocksize+1
	tmp = bytedata[_start:end]
	buf = IOBuffer(append=true)
	zWriter = ZlibCompressorStream(buf,level=level)
	write(zWriter, tmp)
	write(zWriter, TranscodingStreams.TOKEN_END)
	flush(zWriter)
	push!(compr_dat,read(buf))
	last_blocksize = headertype(length(tmp))
	compressed_blocksizes = map(x->headertype(length(x)),compr_dat)
	nvtuhead = VTUHeader(vtuheader.num_blocks,vtuheader.blocksize,last_blocksize,compressed_blocksizes)
	return nvtuhead, vcat(compr_dat...)
end

function updateappendeddata!(i,j,el,appendeddata,data,offsets,type,headertype,compressed_dat)
	vtuheader = data.header[j]
	dat = data.data[j].dat
	rangemin = minimum(filter(!isnan,dat))
	rangemax = maximum(filter(!isnan,dat))
	rawdat = appendeddata[1].content[1]
	if compressed_dat
		nvtuhead, compr_dat = compressdat(vtuheader,dat,headertype)
		head = bytes(nvtuhead)
		base_head = Base64.base64encode(head)
		base_dat = Base64.base64encode(compr_dat)
		if i < length(offsets)
			newoff = length(base_head)+length(base_dat)
			diffoff = offsets[i]+newoff-offsets[i+1]
			appendeddata[1].content[1] = rawdat[1:offsets[i]+1]*base_head*base_dat*rawdat[2+offsets[i+1]:end]
		else
			newoff = 0
			diffoff = 0
			appendeddata[1].content[1] = rawdat[1:offsets[i]+1]*base_head*base_dat
		end
	else
		a = headertype(length(dat)*sizeof(type))
		b = dat
		aa = reinterpret(UInt8,[a])
		bb = reinterpret(UInt8,b)
		base = Base64.base64encode(vcat(aa,bb))
		if i < length(offsets)
			newoff = length(base)
			diffoff = offsets[i]+newoff-offsets[i+1]
			appendeddata[1].content[1] = rawdat[1:offsets[i]+1]*base*rawdat[2+offsets[i+1]:end]
		else
			newoff = 0
			diffoff = 0
			appendeddata[1].content[1] = rawdat[1:offsets[i]+1]*base
		end
	end
	if i < length(offsets)
		offsets[i+1:end] .+= diffoff
	end
	setAttribute(el,"RangeMin",rangemin)
	setAttribute(el,"RangeMax",rangemax)
	setAttribute(el,"offset",offsets[i])
	return nothing
end

function updatedataarray!(i,j,el,data,offsets,type,headertype,compressed_dat)
	vtuheader = data.header[j]
	dat = data.data[j].dat
	bytedata = reinterpret(UInt8,dat)
	rangemin = minimum(filter(!isnan,dat))
	rangemax = maximum(filter(!isnan,dat))
	if compressed_dat
		nvtuhead, compr_dat = compressdat(vtuheader,dat,headertype,-1)
		head = bytes(nvtuhead)
		base_head = Base64.base64encode(head)
		base_dat = Base64.base64encode(compr_dat)
		el.content[1] = base_head*base_dat
	else
		a = headertype(length(dat)*sizeof(type))
		b = dat
		aa = reinterpret(UInt8,[a])
		bb = reinterpret(UInt8,b)
		el.content[1] = Base64.base64encode(vcat(aa,bb))
	end
	setAttribute(el,"RangeMin",rangemin)
	setAttribute(el,"RangeMax",rangemax)
	return nothing
end

function update_xml!(vtufile::VTUFile)
	dataarrays,appendeddata,headertype,offsets,data,compressed_dat = vtufile.dataarrays,vtufile.appendeddata,vtufile.headertype,vtufile.offsets,vtufile.data,vtufile.compressed_dat
	println("compressed ",compressed_dat)
	j = 0
	update_data!(data)
	if !isempty(appendeddata) && !isempty(appendeddata[1].content)
		println("length appended = ",length(appendeddata[1].content[1]))
	end
	for (i,el) in enumerate(dataarrays)
		_type = replace(getAttribute(el,"type"),"\""=>"")
		type = eval(Meta.parse(_type))
		_format = replace(getAttribute(el,"format"),"\""=>"")
		_name = replace(getAttribute(el,"Name"),"\""=>"")
		if findfirst(x->x==_name,vtukeywords.interpolation_keywords) != nothing || findfirst(x->x==_name,vtukeywords.uncompress_keywords) != nothing
			j+=1
			if _format == "appended"
				updateappendeddata!(i,j,el,appendeddata,data,offsets,type,headertype,compressed_dat)
			else
				updatedataarray!(i,j,el,data,offsets,type,headertype,compressed_dat)
			end
		else
			if _format == "appended"
				setAttribute(el,"offset",offsets[i])
			end
		end
	end
end

function timestamp()
	str = string(Dates.unix2datetime(time()))
	str = replace(str,"-"=>"_")
	str = replace(str,":"=>"_")
	str = replace(str,"."=>"_")
end

"""
	write(vtufile::VTUFile,add_timestamp::Bool=true)

Writes a [`VTUFile`](@ref) to destination `vtufile.name`

# Arguments
- `vtufile::VTUFile`: VTU file
- `add_timestamp::Bool`: adds a timestamp to `vtufile.name` if `add_timestamp==true`
"""
function Base.write(vtufile::VTUFile, add_timestamp::Bool=true)
	update_xml!(vtufile)
	name = vtufile.name
	if add_timestamp
		splitstr = split(name,".vtu")
		@assert length(splitstr) == 2 && isempty(splitstr[end])
		name = splitstr[1] * "_" * timestamp() * ".vtu"
	end
	f = open(name,"w")
	writeXMLElement(f,vtufile.xmlroot)
	close(f)
end

"""
	rename!(vtuf::VTUFile,name::String)

Renames a [`VTUFile`](@ref) to `name`. This determines the write-location.

# Arguments
- `vtuf::VTUFile`: VTU file
- `name::String`: Name, or write location, of vtu file
"""
function rename!(vtuf::VTUFile,name::String)
	vtuf.name = name
	return nothing
end

"""
	deepcopy(vtuf::VTUFile)

Deepcopy of a [`VTUFile`](@ref).
Only data that will be altered by algebraic operations is deepcopied at the moment due to efficiency.

# Arguments
- `vtuf::VTUFile`: VTU file
"""
function Base.deepcopy(vtuf::VTUFile)
	return VTUFile(vtuf.name,vtuf.xmlroot,vtuf.dataarrays,vtuf.appendeddata,vtuf.headertype,vtuf.offsets,similar(vtuf.data),vtuf.compressed_dat)
end

"""
	similar(vtuf::VTUFile)

Congruent to [`deepcopy`](@ref).

# Arguments
- `vtuf::VTUFile`: VTU file
"""
function Base.similar(vtuf::VTUFile)
	ret = deepcopy(vtuf)
	#foreach(fill_zeros,ret.data.data)
	return ret
end

"""
	fill(ret::VTUFile, c::Float64)

Fills all entries of all interpolation data fields with a constant value `c::Float64`.

# Arguments
- `ret::VTUFile`: VTU file
- `c::Float64`: Float all data is filled with
"""
function Base.fill!(ret::VTUFile, c::Float64)
	fill!(ret.data,c)
	return nothing
end

"""
	zero(vtu::VTUFile)

Shortcut for [`fill!`](@ref)(vtu,0.0).

# Arguments
- `vtu::VTUFile`: VTU file
"""
function Base.zero(vtu::VTUFile)
	ret = similar(vtu)
	fill!(ret,0.0)
	return ret
end

"""
	one(vtu::VTUFile)

Shortcut for [`fill!`](@ref)(vtu,1.0).

# Arguments
- `vtu::VTUFile`: VTU file
"""
function Base.one(vtu::VTUFile)
	ret = similar(vtu)
	fill!(ret,1.0)
	return ret
end

"""
	empty(vtu::VTUFile)

Empties each data field in [`VTUData`](@ref).

# Arguments
- `vtu::VTUFile`: VTU file
"""
function Base.empty!(vtu::VTUFile)
	for i = 1:length(vtu.data.data)
		empty!(vtu.data.data[i].dat)
	end
	return nothing
end

"""
	empty(vtu::VTUFile)

Returns a empty copy of `vtu`.

# Arguments
- `vtu::VTUFile`: VTU file
"""
function Base.empty(vtu::VTUFile)
	ret = similar(vtu)
	empty!(ret)
	return ret
end

function addPointData!(vtu::VTUFile,name::String,dat,interp_dat) #does only work for scalar fields and last field in result has to be scalar
	pointdata = getElements(vtu.xmlroot,"PointData")
	n = length(vtu.data.idat)+1
	n_pd = length(pointdata[1].content)+1
	apd = deepcopy(pointdata[1].content[end])
	push!(pointdata[1].content,apd)
	setAttribute(apd,"Name",name)
	if hasAttributekey(apd,"NumberOfComponents")
		@assert parse(Int,replace(getAttribute(apd,"NumberOfComponents"),"\""=>"")) == 1
	end
	insert!(vtu.data.data,n,dat)
	off = vtu.offsets[n_pd-1]
	diffoff = vtu.offsets[n_pd]-off
	vtu.appendeddata[1].content[1] = vtu.appendeddata[1].content[1][1:off+1]*vtu.appendeddata[1].content[1][off+2:off+1+diffoff]*vtu.appendeddata[1].content[1][off+2:end];
	insert!(vtu.offsets,n_pd,off)
	vtu.offsets[n_pd:end] .+= diffoff
	insert!(vtu.data.header,n,deepcopy(vtu.data.header[n-1]))
	push!(vtu.data.idat,n)
	push!(vtu.data.interp_data,interp_dat)
	insert!(vtu.data.names,n,name)
	push!(vtukeywords.interpolation_keywords,name)
	push!(vtukeywords.interpolation_keywords,"\""*name*"\"")
	vtu.dataarrays = getElements(vtu.xmlroot, "DataArray")
	return nothing
end

function addIntegrityTemp!(vtu::VTUFile,init=nothing,name::String="temperature_interpolated")
	ind = findall(x->replace(x,"\""=>"")==name,vtu.data.names)[1]
	dat = deepcopy(vtu.data.data[ind])
	interp_dat = deepcopy(vtu.data.interp_data[vtu.data.idat[ind]])
	if init != nothing
		interp_dat.dat .+= init.data.interp_data[ind].dat
	end
	interp_dat.dat ./= 373.15
	addPointData!(vtu,"integrity_temp",dat,interp_dat)
	return nothing
end

function addIntegrityDilat!(vtu::VTUFile,init=nothing,name::String="sigma")
	ind = findall(x->replace(x,"\""=>"")==name,vtu.data.names)[1]
	sigar = deepcopy(vtu.data.interp_data[ind].dat)
	if init != nothing
		sigar .+= init.data.interp_data[ind].dat
	end
	dat = deepcopy(vtu.data.data[end])
	interp_dat = deepcopy(vtu.data.interp_data[vtu.data.idat[end]])
	for (ig,i) = enumerate(1:4:length(sigar))
		_sig = SMatrix{3,3}(sigar[i], sigar[i+3], 0.0, sigar[i+3], sigar[i+1], 0.0, 0.0, 0.0, sigar[i+2])
		evals = sort(eigvals(_sig))
		r = (evals[end]-evals[1])/2
		phi = 1
		c = 1
		R = c*cos(phi)-(evals[end]+evals[1])/2*sin(phi)
		interp_dat.dat[ig] = r/R
	end
	addPointData!(vtu,"integrity_dilat",dat,interp_dat)
	return nothing
end

function addIntegrityFluid!(vtu::VTUFile,init=nothing,name1::String="sigma",name2::String="pressure_interpolated")
	ind1 = findall(x->replace(x,"\""=>"")==name1,vtu.data.names)[1]
	ind2 = findall(x->replace(x,"\""=>"")==name2,vtu.data.names)[1]
	sigar = deepcopy(vtu.data.interp_data[ind1].dat)
	par = deepcopy(vtu.data.interp_data[ind2].dat)
	if init != nothing
		sigar .+= init.data.interp_data[ind1].dat
		par .+= init.data.interp_data[ind2].dat
	end
	dat = deepcopy(vtu.data.data[end])
	interp_dat = deepcopy(vtu.data.interp_data[vtu.data.idat[end]])
	for (ig,i) = enumerate(1:4:length(sigar))
		_sig = SMatrix{3,3}(sigar[i], sigar[i+3], 0.0, sigar[i+3], sigar[i+1], 0.0, 0.0, 0.0, sigar[i+2])
		evals = sort(eigvals(_sig))
		p = par[ig]
		interp_dat.dat[ig] = evals[end]+p
	end
	addPointData!(vtu,"integrity_fluid",dat,interp_dat)
	return nothing
end

function addIntegrityChecks!(vtu::VTUFile,init=nothing)
	addIntegrityTemp!(vtu,init)
	addIntegrityDilat!(vtu,init)
	addIntegrityFluid!(vtu,init)
	return nothing
end

"""
	getindex(vtu::VTUFile, str::String)

Return the `str`-named data field of vtu.

# Arguments
- `vtu::VTUFile`: VTU file
- `str::String`: Name of data field

# Example
```julia
displacementfield = vtu["displacement"]
```
"""
function Base.getindex(vtu::VTUFile, str::String)
	ind = findfirst(x->replace(x,"\""=>"")==str,vtu.data.names)
	if ind ∈ vtu.data.idat
		ind = findfirst(x->x==ind,vtu.data.idat)
		return vtu.data.interp_data[ind].dat
	else
		return vtu.data.data[ind].dat
	end
end

max_stars_repo_name: baxmittens/VTUFileHandler

max_issues_repo_name: baxmittens/VTUFileHandler

max_forks_repo_name: baxmittens/VTUFileHandler


text:
"""
    create_sum_invested_variable(model::JuMP.Model, w)

Creates decision variable representing the invested money in variable-return assets: ``||w||``.

Arguments:
 - `model::JuMP.Model`: JuMP upper level portfolio optimization model.
 - `w`: portfolio optimization investment variable ("weights").
"""
function create_sum_invested_variable(model::JuMP.Model, w)
    if !haskey(object_dictionary(model), :sum_invested)
        sum_invested = @variable(model, sum_invested)
        @constraint(model, [sum_invested; w] in MOI.NormOneCone(length(w) + 1))
        return sum_invested
    else
        return model[:sum_invested]
    end
end

"""
    base_model(numA) -> model

Creates a JuMP model with generic PO variable and constraints:
    - Investment vector of variables `w` (portfolio weights if `current_wealth = 1`).
    - Invested monney should be lower current wealth.

Returns the model (which has a reference to the decision variable named `w`).

Arguments:
 - `numA::Integer`: number of assets in portfolio.

Optional Keywork Arguments:
 - `current_wealth::Real = 1.0`: Current available wealth to be invested.
 - `allow_borrow::Bool=false`: Short variables work as borrowing money.
"""
function base_model(numA::Integer; current_wealth::Real=1.0, allow_borrow::Bool=false)
    model = Model()
    w = @variable(model, w[i=1:numA])
    @variable(model, sum_invested)
    if allow_borrow
        @constraint(model, sum_invested == sum(w))
    else
        @constraint(model, [sum_invested; w] in MOI.NormOneCone(length(w) + 1))
    end
    @constraint(model, sum_invested <= current_wealth)
    return model
end

"""
    readjust_volumes(decision_values::Array{Float64,1}; current_wealth::Real = 1.0)

Adjust volumes to be feasible under current wealth.

Arguments:
 - `decision_values::Array{Float64,1}`: Decided values to be invested.

Optional Keywork Arguments:
 - `current_wealth::Real = 1.0`: Current available wealth to be invested.
"""
function readjust_volumes!(decision_values::Array{Float64,1}; current_wealth::Real = 1.0)
    decision_values_adjusted = deepcopy(decision_values)
    if norm(decision_values, 1) > current_wealth
        decision_values_adjusted = decision_values / norm(decision_values, 1)
    end
    return decision_values_adjusted
end

"""
    compute_solution(model, solver) -> Array{Float64,1}

Solves optimization model using provides solver and returns solutions adjusted to current wealth.

Arguments:
 - `model::JuMP.Model`: JuMP portfolio optimization model.
 - `solver`: Appropriate solver for the provided model.

Optional Keywork Arguments:
 - `current_wealth::Real = 1.0`: Current available wealth to be invested.
 - `decision_variable = model[:w]`: Portfolio optimization investment variable reference.
"""
function compute_solution(
    model::JuMP.Model, solver;
    current_wealth::Real = 1.0,
    decision_variable = model[:w]
)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"
    decision_values = value.(decision_variable)
    decision_values = readjust_volumes!(decision_values; current_wealth = current_wealth)
    return decision_values
end

ext: jl


text:
# Generate random unitaries and density matrices for testing

export random_unitary, random_densitymatrix

""" `U = random_unitary(dim, real=true/false, complex=true/false)`

Generate a random unitary of dimension *dim*. Depending on *real* and *complex* the matrix has real and complex parts.

"""

function random_unitary(dim::Integer; real::Bool = true, complex::Bool = false)

  	@assert dim > 1 "Dimension must be at least 2."

  	generated = randn(dim, dim)

  	if complex || !real
    		generated += im * randn(dim, dim)
  	end

	# Compute the QR decomposition
  	(Q, R) = qr(generated)

	# Get the upper diagonal part of R
  	F = diagm(diag(R))

  	R = sign(diagm(diag(R)))
  	R[find((x)->  x == 0, R)] = 1
  	F = diagm(diag(R))

  	return Q*F
end

""" `rho = random_densitymatrix(dim)`

Generate a random density matrix in dimension *dim*.

"""

function random_densitymatrix(dim::Integer)

	@assert dim > 1 "Dimension must be at least 2."

        # Draw random numbers
	rList = rand(dim);

	# Normalize to form a distribution
	s = sum(rList);
	nList = rList/s;

	# Construct a matrix with the numbers on the diagonal
	M = diagm(nList);

	# Draw a random unitary
   	U = random_unitary(dim, real = true, complex = true)

	# Roate the diagonal matrix
	return(U*M*U');
end


text:
function KS_solve_SCF_smearing( pw::PWGrid,
                       V_ionic, Focc_in, Nstates::Int64, Nelectrons::Int64;
                       β = 0.5, E_NN = nothing,
                       rho0 = nothing,
                       Potentials0 = nothing,
                       NiterMax=100,
                       verbose=false )

    Ngwx = pw.gvecw.Ngwx
    Ns = pw.Ns
    Npoints = prod(Ns)
    ΔV = pw.Ω/Npoints

    # starting rho
    if rho0 == nothing
        srand(1234)
        v = rand( ComplexF64, Ngwx, Nstates )
        v = ortho_gram_schmidt( v )
        rho = calc_rho( pw, Focc_in, v )
    else
        srand(1234)
        v = rand( ComplexF64, Ngwx, Nstates ) # What if v0 is provided ?
        rho = copy( rho )
    end

    # starting Potentials
    if Potentials0 == nothing
        # initialize potential: Ionic, Hartree, XC
        V_Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        V_xc = excVWN( rho ) + rho .* excpVWN( rho )
        Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )
    else
        Potentials = PotentialsT( Potentials0.Ionic,
                                  Potentials0.Hartree,
                                  Potentials0.XC )
    end

    Etot_old = 0.0
    rho_old = copy(rho)

    if( E_NN != nothing )
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, E_NN )
    else
        Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
    end

    λ = zeros(Nstates)

    const kT = 0.01

    MIXDIM = 4
    df = zeros(Float64,Npoints,MIXDIM)
    dv = zeros(Float64,Npoints,MIXDIM)

    for iter = 1:150

        λ, v = diag_lobpcg( pw, Potentials, v, verbose_last=false )

        Focc, E_fermi = calc_Focc(λ, Nelectrons, kT)
        Entropy = calc_entropy(Focc, kT)

        for ist = 1:Nstates
            @printf("%5d %18.10f %18.10f\n", ist, λ[ist], Focc[ist] )
        end

        println("sum(Focc) = ", sum(Focc))
        println("integFocc = ", sum_upto_E_fermi(Focc, λ, E_fermi))
        println("E_fermi = ", E_fermi)

        rho_new = calc_rho( pw, Focc, v )
        diffRho = norm(rho_new - rho)

        #@printf("Using simple mix with β = %10.5f\n", β)
        #rho = β*rho_new[:] + (1-β)*rho[:]

        @printf("Using andersonmix with β = %f, MIXDIM = %d\n", β, MIXDIM)
        rho = andersonmix!( rho, rho_new, β, df, dv, iter, MIXDIM )

        for ip = 1:Npoints
            if rho[ip] < 1e-12
                rho[ip] = 1e-12
            end
        end
        #
        integRho = sum(rho)*ΔV
        #@printf("integRho = %18.10f\n", integRho)
        diffNelec = abs(integRho - Nelectrons)
        if diffNelec > 1e-12
            println("WARNING: Rescaling electron density, diffNelec = ", diffNelec)
            rho[:] = rho[:]*Nelectrons/integRho
        end

        V_Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        V_xc = excVWN( rho ) + rho .* excpVWN( rho )
        Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )

        # Calculate energies and update potetials
        Energies = calc_Energies( pw, Potentials, Focc, v, Energies.NN )

        Etot = Energies.Total + Entropy
        diffE = abs( Etot - Etot_old )

        #
        @printf("SCF: %8d %18.10f %18.10e %18.10e\n", iter, Etot, diffE, diffRho )

        if diffE < 1e-7
            @printf("SCF is is converged: iter: %d , diffE = %10.7e\n", iter, diffE)
            @printf("Etot           = %18.10f\n", Etot)
            @printf("Entropy        = %18.10f\n", Entropy)
            @printf("Energies.Total = %18.10f\n", Etot)
            break
        end
        #
        Etot_old = Etot
    end

    exit()

    return Energies, Potentials, v, λ

end


text:
export ncores, nphysicalcores, systemload, normload, setnumworkers

@static if Sys.isapple() ncores() = asint(read(`sysctl -n hw.ncpu`, String)) end
@static if Sys.isapple() nphysicalcores() = asint(read(`sysctl -n hw.physicalcpu`, String)) end
@static if Sys.islinux() ncores() = asint(read(`nproc --all`, String)) end
@static if Sys.islinux() nphysicalcores() = div(ncores(),2) end

rmcomma(a) = a[end]==',' ? a[1:end-1] : a
systemload() = parsefloat(rmcomma(split(read(`uptime`, String))[end-2]))
normload() = systemload()/ncores()

function setnumworkers(n)
    oldn = nprocs()-1
    if n>oldn
        return addprocs(n-oldn)
    elseif n<oldn
        rmthese = workers()[end-(oldn-n)+1:end]
    end
end


text:
__precompile__()

module DataArrays
    using Base.Cartesian, Compat, Reexport
    import Compat.String
    @reexport using StatsBase

    const DEFAULT_POOLED_REF_TYPE = UInt32

    import Base: ==, !=, >, <, >=, <=, +, -, *, !, &, |, $, ^, /,
      .==, .!=, .>, .<, .>=, .<=, .+, .-, .*, .%, ./, .\, .^

    export @data,
           @pdata,
           AbstractDataArray,
           AbstractDataMatrix,
           AbstractDataVector,
           allna,
           anyna,
           array,
           autocor,
           compact,
           cut,
           data,
           DataArray,
           DataMatrix,
           DataVector,
           dropna,
           each_failna,
           each_dropna,
           each_replacena,
           EachFailNA,
           EachDropNA,
           EachReplaceNA,
           FastPerm,
           getpoolidx,
           gl,
           head,
           isna,
           levels,
           NA,
           NAException,
           NAtype,
           padNA,
           pdata,
           percent_change,
           PooledDataArray,
           PooledDataMatrix,
           PooledDataVecs,
           PooledDataVector,
           reldiff,
           reorder,
           rep,
           replace!,
           setlevels!,
           setlevels,
           tail,
           xtab,
           xtabs

    include("utils.jl")
    include("natype.jl")
    include("abstractdataarray.jl")
    include("dataarray.jl")
    include("pooleddataarray.jl")
    include("datavector.jl")
    include("indexing.jl")
    include("datamatrix.jl")
    include("linalg.jl")
    include("operators.jl")
    include("reduce.jl")
    include("reducedim.jl")
    include("broadcast.jl")
    include("sort.jl")
    include("extras.jl")
    include("grouping.jl")
    include("statistics.jl")
    include("predicates.jl")
    include("literals.jl")
    include("deprecated.jl")

    Base.@deprecate removeNA dropna
    Base.@deprecate each_failNA each_failna
    Base.@deprecate each_replaceNA each_replacena
    Base.@deprecate set_levels setlevels
    Base.@deprecate set_levels! setlevels!
end


text:
"""
    Cylinder(x, y, zmin, zmax, radius)

Create a `Cylinder` that can be used for spatial lookup.

See also `AbstractRegion`
"""
struct Cylinder{T <: Real} <: AbstractRegion{T}
    x::T
    y::T
    zmin::T
    zmax::T
    radius::T
end

"""
    volume(cylinder::Cylinder)

Returns the volume of the `cylinder`.
"""
@inline function volume(cylinder::Cylinder)
    π * cylinder.radius * cylinder.radius * (cylinder.zmax - cylinder.zmin)
end

function boundingbox(cylinder::Cylinder{T}) where T
    xmin = cylinder.x - cylinder.radius
    xmax = cylinder.x + cylinder.radius
    ymin = cylinder.y - cylinder.radius
    ymax = cylinder.y + cylinder.radius
    return BoundingBox{T}(xmin, ymin, cylinder.zmin, xmax, ymax, cylinder.zmax)
end

function in(p::StaticVector{3, <:Real}, cylinder::Cylinder)
    r² = cylinder.radius * cylinder.radius
    Δx = p[1] - cylinder.x
    Δy = p[2] - cylinder.y
    return Δx*Δx + Δy*Δy <= r² && p[3] >= cylinder.zmin && p[3] <= cylinder.zmax
end


text:
using AstroImages, FITSIO, Random, WCS, ImageBase, Test, WCS, Statistics

import AstroImages: _float, render

@testset "Conversion to float and fixed-point" begin
    @testset "Float" begin
        for T in (Float16, Float32, Float64)
            @test _float(T(-9.8)) === T(-9.8)
            @test _float(T(12.3)) === T(12.3)
        end
    end
    @testset "Integers" begin
        for (UIT, SIT) in ((UInt8,  Int8),
                           (UInt16, Int16),
                           (UInt32, Int32),
                           (UInt64, Int64))
            N = sizeof(UIT) * 8
            NT = Normed{UIT, N}
            maxint = UIT(big(2) ^ (N - 1))
            @test _float(typemin(UIT)) === _float(typemin(SIT)) === NT(0)
            @test _float(UIT(85)) === reinterpret(NT, UIT(85))
            @test _float(SIT(-85)) === _float(UIT(-85 + big(maxint)))
            @test _float(SIT(115)) === _float(UIT(115) + maxint)
            @test _float(typemax(UIT)) === _float(typemax(SIT)) === NT(1)
        end
    end
end

@testset "FITS and images" begin
    fname = tempname() * ".fits"
    for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
              Float32, Float64]
        data = reshape(T[1:100;], 5, 20)
        FITS(fname, "w") do f
            write(f, data)
        end
        @test load(fname, 1) == data
        @test load(fname, (1, 1)) == (data, data)
        img = AstroImage(fname)
        rendered_img = imview(img)
        @test eltype(rendered_img) <: RGBA
    end
    rm(fname, force = true)
end

@testset "default handler" begin
    fname = tempname() * ".fits"
    @testset "less dimensions than 2" begin
        data = rand(2)
        FITS(fname, "w") do f
            write(f, data)
        end
        @test ndims(AstroImage(fname)) == 1
    end

    @testset "no ImageHDU" begin
        ## Binary table
        indata = Dict{String, Array}()
        i = length(indata) + 1
        indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
        i += 1
        indata["col$i"] = ones(Bool, 20)  # Bool column
        i += 1
        indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
        i += 1
        indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
        indata["vcol"] = [randstring(j) for j=1:20]  # variable length column
        indata["VCOL"] = [collect(1.:j) for j=1.:20.] # variable length

        FITS(fname, "w") do f
            write(f, indata; varcols=["vcol", "VCOL"])
            @test_throws Exception AstroImage(f)
        end
    end

    @testset "Opening AstroImage in different ways" begin
        data = rand(2,2)
        FITS(fname, "w") do f
            write(f, data)
        end
        f = FITS(fname)
        header = read_header(f[1])
        @test AstroImage(fname, 1) isa AstroImage
        @test AstroImage(f, 1) isa AstroImage
        @test AstroImage(data, header) isa AstroImage
        close(f)
    end

    @testset "Image HDU is not at 1st position" begin
        ## Binary table
        indata = Dict{String, Array}()
        i = length(indata) + 1
        indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
        i += 1
        indata["col$i"] = ones(Bool, 20)  # Bool column
        i += 1
        indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
        i += 1
        indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
        indata["vcol"] = [randstring(j) for j=1:20]  # variable length column
        indata["VCOL"] = [collect(1.:j) for j=1.:20.] # variable length

        FITS(fname, "w") do f
            write(f, indata; varcols=["vcol", "VCOL"])
            write(f, rand(2, 2))
        end

        @test @test_logs (:info, "Image was loaded from HDU 3") AstroImage(fname) isa AstroImage
    end
    rm(fname, force = true)
end

@testset "Utility functions" begin
   @test size(AstroImage(rand(10,10))) == (10,10)
   @test length(AstroImage(rand(10,10))) == 100
end

@testset "multi wcs AstroImage" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")
    inhdr = FITSHeader([
            "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", "HISTORY",
            "CRVAL1a",
            "CRVAL2a",
            "CRPIX1a",
            "CRPIX2a",
            "CDELT1a",
            "CDELT2a",
            "CTYPE1a",
            "CTYPE2a",
            "CUNIT1a",
            "CUNIT2a",

            "CRVAL1b",
            "CRVAL2b",
            "CRPIX1b",
            "CRPIX2b",
            "CDELT1b",
            "CDELT2b",
            "CTYPE1b",
            "CTYPE2b",
            "CUNIT1b",
            "CUNIT2b",
        ],
        [
            1.0, 1, true, "string value", nothing, nothing,
            0.5,
            89.5,
            1,
            1,
            1,
            -1,
            "RA---TAN",
            "DEC--TAN",
            "deg     ",
            "deg     ",

            0.5,
            89.5,
            1,
            1,
            1,
            -1,
            "RA---TAN",
            "DEC--TAN",
            "deg     ",
            "deg     ",
        ],
        [
            "floating point keyword", "", "boolean keyword", "string value", "this is a comment", "this is a history",
            "",
            "",
            "",
            "",
            "",
            "",
            "Terrestrial East Longitude",
            "Terrestrial North Latitude",
            "",
            "",

            "",
            "",
            "",
            "",
            "",
            "",
            "Terrestrial East Longitude",
            "Terrestrial North Latitude",
            "",
            "",
        ])

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)
    close(f)

    img = AstroImage(fname)
    f = FITS(fname)
    @test length(wcs(img)) == 2
    @test WCS.to_header(wcs(img,1)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1])
    @test WCS.to_header(wcs(img,2)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[2])

    img = AstroImage(f)
    @test length(wcs(img)) == 2
    @test WCS.to_header(wcs(img,1)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1])
    @test WCS.to_header(wcs(img,2)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[2])
    close(f)
end

##
@testset "imview" begin

    arr1 = permutedims(reshape(1:9,3,3))
    img = AstroImage(arr1)

    @test imview(arr1) == imview(img)
    @test imview(img) isa AstroImage
    @test !(imview(arr1) isa AstroImage)

    img_rendered_1 = imview(img, clims=(1,9), stretch=identity, contrast=1, bias=0.5, cmap=nothing)

    # Image Orientation
    @test CartesianIndex(3,1) == argmin(Gray.(img_rendered_1))
    @test CartesianIndex(1,3) == argmax(Gray.(img_rendered_1))

    # Rendering Basics
    @test allunique(img_rendered_1)
    # It is intended that the rendered image is flipped vs it's data
    @test img_rendered_1[3,1] == RGBA(0,0,0,1)
    @test img_rendered_1[1,3] == RGBA(1,1,1,1)
    @test all(p -> p.r==p.g==p.b && p.alpha==1, img_rendered_1)

    # Limits
    img_rendered_2 = imview(img, clims=(3,7), stretch=identity, contrast=1, bias=0.5, cmap=nothing)
    @test length(unique(img_rendered_2)) == 5
    @test count(==(RGBA(0,0,0,1)), img_rendered_2) == 3
    @test count(==(RGBA(1,1,1,1)), img_rendered_2) == 3

    # Calculated limits
    @test img_rendered_1 == imview(img, clims=extrema, stretch=identity, contrast=1, bias=0.5, cmap=nothing)
    img_rendered_3 = imview(img, clims=Zscale(), stretch=identity, contrast=1, bias=0.5, cmap=nothing)
    img_rendered_4 = imview(img, clims=Percent(100), stretch=identity, contrast=1, bias=0.5, cmap=nothing)
    @test img_rendered_1 == img_rendered_3
    @test img_rendered_1 == img_rendered_4

    # Stretching
    for stretchfunc in (sqrtstretch, asinhstretch, powerdiststretch, logstretch, powstretch, squarestretch, sinhstretch)
        img_rendered_5 = imview(arr1, clims=(1,9), stretch=stretchfunc, contrast=1, bias=0.5, cmap=nothing)
        @test extrema(Gray.(img_rendered_5)) == (0,1)
        manual_stretch = stretchfunc.(AstroImages.clampednormedview(arr1,(1,9)))
        @test Gray.(img_rendered_5) ≈
            N0f8.((manual_stretch.-minimum(manual_stretch)) ./
                (maximum(manual_stretch)-minimum(manual_stretch)))'[end:-1:begin,:]
    end

    # Contrast/Bias
    @test Gray.(imview(img, clims=extrema, stretch=identity, contrast=1, bias=0.6, cmap=nothing)) ==
            N0f8.(clamp.(N0f8.(Gray.(img_rendered_1)) .- 0.1,false,true))

    img_rendered_5 = imview(arr1, clims=(1,9), stretch=sqrtstretch, contrast=0.5, bias=0.5, cmap=nothing)

    # Missing/NaN
    for m in (NaN, missing)
        arr2 = [
            1 2 3
            4 m 6
            7 8 9
        ]
        @test imview(arr2) == imview(AstroImage(arr2))
        @test imview(arr2)[2,2].alpha == 0
        @test 8 == count(img_rendered_1 .== imview(arr2, clims=(1,9), stretch=identity, contrast=1, bias=0.5, cmap=nothing))
    end

    img_rendered_6 = imview([1, 2, NaN, missing, -Inf, Inf], clims=extrema)
    img_rendered_6b = imview([1, 2], clims=extrema)

    @test img_rendered_6[1] == img_rendered_6b[1]
    @test img_rendered_6[2] == img_rendered_6b[2]
    @test img_rendered_6[1].alpha == 1
    @test img_rendered_6[2].alpha == 1
    @test img_rendered_6[3].alpha == 0
    @test img_rendered_6[4].alpha == 0
    @test img_rendered_6[5] == RGBA(0,0,0,1)
    @test img_rendered_6[6] == RGBA(1,1,1,1)
end

##

1

##

# @testset "multi file AstroImage" begin
#     fname1 = tempname() * ".fits"
#     f = FITS(fname1, "w")
#     inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
#                         "HISTORY"],
#                     ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing],
#                     ["",
#                         "",
#                         "",
#                         "floating point keyword",
#                         "",
#                         "boolean keyword",
#                         "string value",
#                         "this is a comment",
#                         "this is a history"])

#     indata1 = reshape(Int[1:100;], 5, 20)
#     write(f, indata1; header=inhdr)
#     close(f)

#     fname2 = tempname() * ".fits"
#     f = FITS(fname2, "w")
#     indata2 = reshape(Int[1:100;], 5, 20)
#     write(f, indata2; header=inhdr)
#     close(f)

#     fname3 = tempname() * ".fits"
#     f = FITS(fname3, "w")
#     indata3 = reshape(Int[1:100;], 5, 20)
#     write(f, indata3; header=inhdr)
#     close(f)

#     img = AstroImage((fname1, fname2, fname3))
#     f1 = FITS(fname1)
#     f2 = FITS(fname2)
#     f3 = FITS(fname3)

#     @test length(img.data) == length(img.wcs) == 3
#     @test img.data[1] == indata1
#     @test img.data[2] == indata2
#     @test img.data[3] == indata3
#     @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) ==
#         WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1])
#     @test eltype(eltype(img.data)) == Int

#     img = AstroImage(Gray, (f1, f2, f3), (1,1,1))
#     @test length(img.data) == length(img.wcs) == 3
#     @test img.data[1] == indata1
#     @test img.data[2] == indata2
#     @test img.data[3] == indata3
#     @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) ==
#         WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1])
#     @test eltype(eltype(img.data)) == Int
#     close(f1)
#     close(f2)
#     close(f3)
#     rm(fname1, force = true)
#     rm(fname2, force = true)
#     rm(fname3, force = true)
# end

# include("plots.jl")
# include("ccd2rgb.jl")


text:
@testset "parallel_PARSDMM" begin
Random.seed!(123)

#test parallel PARSDMM (test is for 2 workers)
    options = PARSDMM_options()
    default_PARSDMM_options(options,Float32)
    options.parallel     = true
    options.evol_rel_tol = 1e-5
    options.feas_tol     = 1e-5
    options.obj_tol      = 1e-5
    options.maxit        = 10000

    x         = randn(Float32,201,100)
    comp_grid = compgrid((1.0, 1.0),(201, 100))
    x         = vec(x)
    m         = deepcopy(x)

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #total variation
    (TD_OP, AtA_diag, dense, TD_n) = get_TD_operator(comp_grid,"TV",options.FL)
    m_min     = 0.0
    m_max     = 0.5f0*norm(TD_OP*x,1)
    set_type  = "l1"
    TD_OP     = "TV"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
    (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

    (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

    result = deepcopy(x1)
    for i=1:length(TD_OP)-1
      copy!(result,x1)
      @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
    end

#test parallel PARSDMM without many explicit BLAS calls (test is for 2 workers)
    #default_PARSDMM_options(options,Float64)
    options.parallel    = true
    options.Blas_active = false

    (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
    (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

    (x2,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
    result = deepcopy(x2)
    for i=1:length(TD_OP)-1
      copy!(result,x2)
      @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
    end

#test if with and without expicit Blas calls results are the same
@test isapprox(x1,x2,rtol=1f-4)

#run serial PARSDMM and compare results
    options.parallel    = false

    (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
    (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

    (x3,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
    result = deepcopy(x3)
    for i=1:length(TD_OP)-1
      copy!(result,x3)
      @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
    end

@test isapprox(x1,x3,rtol=5*1f-4)

#test parallel PARDMM with a JOLI operator in a projector, but not in TD_OP
  options = PARSDMM_options()
  default_PARSDMM_options(options,Float32)
  options.parallel     = false
  options.evol_rel_tol = 1e-5
  options.feas_tol     = 1e-5
  options.obj_tol      = 1e-5
  options.maxit        = 10000

  constraint = Vector{SetIntersectionProjection.set_definitions}()

  #DFT l1 constraints
  (TD_OP, AtA_diag, dense, TD_n)=get_TD_operator(comp_grid,"DFT",options.FL)
  m_min     = 0.0
  m_max     = 0.5f0*norm(TD_OP*x,1)
  set_type  = "l1"
  TD_OP     = "DFT"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

  (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

  result = deepcopy(x1)
  for i=1:length(TD_OP)-1
    copy!(result,x1)
    @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
  end

  #test in parallel and see if results are the same
  options.parallel     = true

  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

  (x2,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

  result = deepcopy(x2)
  for i=1:length(TD_OP)-1
    copy!(result,x2)
    @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
  end

  @test isapprox(x1,x2,rtol=5f-4)

# #test parallel PARDMM with a JOLI operator in TD_OP
# default_PARSDMM_options(options,Float64)
# options.parallel    = false

# constraint=Dict()

# #DFT l1 constraints
# (C, AtA_diag, dense, TD_n)=get_TD_operator(comp_grid,"curvelet",options.FL)
# constraint["use_TD_l1_1"]      = true
# constraint["TD_l1_operator_1"] = "curvelet"
# constraint["TD_l1_sigma_1"]    = 0.5*norm(C*x,1)

# (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,options.FL);
# (TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,TD_Prop,comp_grid,options)

# (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
# result=Vector{typeof(x[1])}(length(x1))
# for i=1:length(TD_OP)-1
#   copy!(result,x1)
#   @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= options.feas_tol
# end

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using MixedModels

include("utilities.jl")
include("statschol.jl")
include("UniformBlockDiagonal.jl")
include("linalg.jl")
include("matrixterm.jl")
include("FactorReTerm.jl")
include("pls.jl")
include("pirls.jl")
include("gausshermite.jl")
include("fit.jl")
include("missing.jl")
include("likelihoodratiotest.jl")


text:
#loads the value weighted crsp index
#also does some light preprocessing
function loadidxcrsp(idxcrspname::String, retainedidxcrspcols::Vector{Symbol};
  idxcrsppath::String = IDXCRSP_PATH,
  idxcrspdateformat::DateFormat = IDXCRSP_DATE_FORMAT,
  idxdatecol = IDX_DATE_COL)

  idxcrsp::DataFrame = CSV.read("$idxcrsppath\\$idxcrspname.csv") |> DataFrame
  idxcrsp.date = (i::Int->Date("$i", idxcrspdateformat)).(idxcrsp[!,idxdatecol])

  #TEST_OUTPUT && (idxcrsp |> CSV.write("$OUT_PATH\\$idxcrspname.csv"))

  select!(idxcrsp, retainedidxcrspcols)
  rename(idxcrsp, Dict(s => Symbol(:P_,s) for s ∈ retainedidxcrspcols))

  return idxcrsp
end

#load the ff data
function loadff(ffname::String;
  ffpath::String = FF_PATH,
  ffdateformat::DateFormat = FF_DATE_FORMAT,
  fieldprefix::String = lowercase("$ffname"),
  rescale::Bool = true)

  ff::DataFrame = CSV.read("$ffpath\\$ffname.csv") |> DataFrame

  #rename the columns
  ffnames::Vector{Symbol} = names(ff)
  ffnames[1] = :date
  for (i,s) ∈ enumerate(ffnames)
    local str::String = string(s)
    if i==1 #date column is a special case
      ffnames[1] = :date
      continue
    end

    str = replace(str, "-"=>"")
    str = replace(str, " "=>"")
    str = replace(str, "_"=>"")
    str = lowercase(str)
    ffnames[i] = Symbol("P_$(fieldprefix)_", str)
  end
  names!(ff, ffnames)
  ff.date = (i::Int->Date("$i", ffdateformat)).(ff.date)

  #now make the colulmns for the index
  idxnames::Vector{Symbol} = setdiff(ffnames, [:date])
  idxnames .= (s->Symbol(s,"_idx")).(idxnames)
  for n ∈ idxnames
    ff[!,n] = Vector{MFloat64}(undef, size(ff,1))
  end

  @mpar for i ∈ 1:length(idxnames) #iterate the focal columns
    idx::Symbol = idxnames[i]
    ret::Symbol = Symbol(replace(string(idx), "_idx"=>""))

    base::Float64 = 1.0
    for r ∈ eachrow(ff) #now iterate down the rows
      if (!ismissing(r[ret]))
        if (r[ret] == -99.99) || (r[ret] < -100.) #missing codes for ff
          r[ret] = missing
        else
          rescale && (r[ret] /= 100)
          base *= 1.0 + r[ret]
          r[idx] = base
        end
      end
    end
  end

  TEST_OUTPUT && (ff |> CSV.write("$OUT_PATH\\$ffname.csv"))
  select!(ff, [:date; idxnames;])

  rename!(ff, Dict(n=>Symbol(replace(string(n), "_idx"=>"")) for n ∈ idxnames))

  return ff
end

#interpolates the index to every day in an interval using the appropriate compounding
function interpolateindex(idxcol::AbstractVector{T},
  datecol::AbstractVector{Date}; tol::Float64 = maximum(idxcol) *10^-13)::Dict{MDate,T} where T<:Any

  local idf::DataFrame

  #check a couple basic types of errors
  (sum((ismissing).(idxcol))>0) && (error("interpolateindex does not support missing values.
    Consider supplying a view w/out missing vals."))
  @assert length(idxcol) == length(datecol)

  #accumulate dates
  begindate::Date = minimum(datecol)
  enddate::Date = maximum(datecol)
  alldates::Vector{Date} = collect(begindate:Day(1):enddate)
  actuals::Dict{Date,T} = Dict(datecol[i]=>idxcol[i] for i ∈ 1:length(datecol))

  #set the values to be interpolated between
  #first allocate space to work
  allN = length(alldates)
  allvals::Vector{T} = Vector{T}(undef, allN)
  nextvalids = Vector{T}(undef, allN)
  prevvalids = Vector{T}(undef, allN)
  nextvaliddates = Vector{MDate}(undef, allN)
  prevvaliddates = Vector{MDate}(undef, allN)

  #identify the most recent previous values
  for i ∈ 1:allN
    (i==1) && continue

    dprev::Date = alldates[i-1]
    if haskey(actuals, dprev)
      prevvalids[i] = actuals[dprev]
      prevvaliddates[i] = dprev
    else
      prevvalids[i] = prevvalids[i-1]
      prevvaliddates[i] = prevvaliddates[i-1]
    end
  end

  #identify the nextvalues by counting backwards
  for i ∈ allN:-1:1
    dnext::Date = alldates[i]
    if haskey(actuals, dnext)
      nextvalids[i] = actuals[dnext]
      nextvaliddates[i] = dnext
    else
      nextvalids[i] = nextvalids[i+1]
      nextvaliddates[i] = nextvaliddates[i+1]
    end
  end

  intervalperiods::Vector{MInt} = ((prev,next)->
    (!ismissing(prev)) && (!ismissing(next)) ?
    (next - prev).value : missing).(prevvaliddates, nextvaliddates)
  ccret::Vector{MFloat64} = (log).(nextvalids ./ prevvalids) ./ intervalperiods

  #finally do the interpolation
  for i ∈ 1:allN
    d::Date = alldates[i]
    if (!haskey(actuals, d)) #interpolation case
      allvals[i] = exp(ccret[i])*allvals[i-1]
    else #no need to interpolate if we have an actual value....
      allvals[i] = actuals[d]
      if i > 1 #but still check the answer if possible
        Δ::Float64 = abs(exp(ccret[i])*allvals[i-1] - allvals[i])
        (Δ > tol) && @warn("Tolerance $tol exceeded with $Δ on row $i.
          Value(t): $(allvals[i]), date: $d, return: $(ccret[i])")
      end
    end
  end

  interpolated::Dict{Date,T} = Dict(alldates[i] => allvals[i] for i ∈ 1:allN)
  return interpolated
end

#derives the returns from an outside cumulative return index
function returnsfromindex!(df::DataFrame, Fidxcol::Symbol; Fretcol::Symbol=Fidxcol,
  sortdf::Bool = true,
  Fdatecol::Symbol = :date, uselogs::Bool = false)::Nothing

  #pre-allocate and sort if needed
  N = size(df,1)
  sortdf && (sort!(df, Fdatecol))
  Fretcol_temp_ = Symbol(Fretcol,:_temp_)
  df[!,Fretcol_temp_] = Vector{MFloat64}(undef, N)

  local base::MFloat64 = missing
  for (i,r) ∈ enumerate(eachrow(df))
    if ismissing(base)
      base = df[i, Fidxcol]
      continue #otherwise we are just adding a 0 return row
    end

    base = ismissing(df[i-1, Fidxcol]) ? base : df[i-1, Fidxcol]

    #compute the return
    if !ismissing(r[Fidxcol])
      r[Fretcol_temp_] = r[Fidxcol]/base - 1.
      (uselogs) && (r[Fretcol_temp_] = log(1. + r[Fretcol_temp_]))
    end
  end

  #most likely scenario for the below is the idx and ret col have the same name
  if Fretcol ∈ names(df)
    select!(df, Not(Fretcol))
  end

  rename!(df, Fretcol_temp_ => Fretcol)

  return nothing
end

#helper function for multiple values
function returnsfromindex!(df::DataFrame, Fidxcols::Vector{Symbol};
  Fretcols::Vector{Symbol} = deepcopy(Fidxcols),
  sortdf::Bool = true, Fdatecol::Symbol = :date, uselogs::Bool = false)::Nothing

  sortdf && (sort!(df, Fdatecol)) #no need to sort multiple times

  for i ∈ 1:length(Fidxcols)
    returnsfromindex!(df, Fidxcols[i], Fretcol=Fretcols[i],
      sortdf = false,  Fdatecol = Fdatecol, uselogs=uselogs)
  end

  return nothing

end

function mergeoutsidedf!(::Type{T}, port::DataFrame, source::AbstractDataFrame,
  Fidxcol::Symbol; Fdatecol::Symbol = :date, interpolate::Bool = true)::Nothing where T<:Any

  local indexvals::Dict{Date,T}

  #make the index column in the output dataframe
  port[!,Fidxcol] = Vector{Union{Missing, T}}(undef, size(port,1))

  #avoid missing data in the input
  ssource::SubDataFrame = view(source, (!ismissing).(source[!, Fidxcol]),:)

  if interpolate
    indexvals = interpolateindex(ssource[!,Fidxcol], ssource[!,Fdatecol])
    maxdate = maximum(ssource[!,Fdatecol])
    mindate = minimum(ssource[!,Fdatecol])
    for r ∈ eachrow(port) #now update port
      if (r[Fdatecol] ≥ mindate) && (r[Fdatecol] ≤ maxdate)
        r[Fidxcol] = indexvals[r[Fdatecol]]
      end
    end
  else
    indexvals = Dict(r[Fdatecol]=>r[Fidxcol] for r ∈ eachrow(ssource))
    for r ∈ eachrow(port) #now update port but check that each entry exists
      if haskey(indexvals, r[Fdatecol])
        r[Fidxcol] = indexvals[r[Fdatecol]]
      end
    end
  end

  return nothing
end

#helper function that works with an aaray of columns
function mergeoutsidedf!(port::DataFrame, source::AbstractDataFrame,
  indexsyms = setdiff!(names(source), [:date]); interpolate::Bool = true)

  for s ∈ indexsyms
    T::Type = eltype(source[!,s])
    mergeoutsidedf!(T, port, source, s, interpolate=interpolate)
  end

  return nothing
end

function mergeoutside(port::DataFrame)
  local vwcrsp::DataFrame
  local rfcrsp::DataFrame
  local ff5::DataFrame
  local ff3m::DataFrame

  #form an intermediate portfolio consisting of the unique panel dates
  outside::DataFrame = DataFrame(date=unique(port.date))

  #load the index portfolios
  vwcrsp = loadidxcrsp(VWCRSP_NAME, RETAINED_VWCRSP_COLS)
  rfcrsp = loadidxcrsp(RFCRSP_NAME, RETAINED_RFCRSP_COLS)
  ff3m = loadff(FF3M_NAME, rescale=false) #use the CRSP version, so no rescale
  ff5 = loadff(FF5_NAME)

  indicesadded::Vector{Symbol} = setdiff([names(vwcrsp); names(rfcrsp);
    names(ff3m); names(ff5)], [:date])

  mergectr::Int = 0
  for df::DataFrame ∈ [vwcrsp, rfcrsp, ff3m, ff5]
    mergectr+=1
    mergeoutsidedf!(outside, df)
  end

  #if TEST_OUTPUT
  #  CSV.write("output\\mergeoutsidetest.csv", port[:,
  #    filter(s::Symbol->string(s)[1:2] ≠ "S_", names(port))])
  #end

  #compute the index returns
  returnsfromindex!(outside, indicesadded, uselogs=false)

  #some housekeeping with the market index
  rename!(outside, Dict(:vwindd=>:mktgross, :t30ind=>:rfr))
  outside.mkt = outside.mktgross .- outside.rfr

  port = join(port, outside, on=:date)

  return port
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#Ost - "E6" - восток, Nord - "N8" - север, West - "W4" - запад, Sud - "S2" - юг
"""
main_function4!(r) - это главная функция которая закрашивает поле в ввиде пирамиды
"""
function main_function4!(r)
    arr = []
    side = Ost
    exit_to_the_lower_left_cell!(r, arr)
    diagonal!(r, West, Nord)
    checks_whether_to_put_a_marker_in_the_cell_and_go_back(r, side)
    return_to_starting_point!(r, arr)
end
"""
diagonal!(r, side1, side2) - функция ставит маркеры в диагональной линии в переданном напралениях
"""
function diagonal!(r, side1, side2)
    while !(isborder(r, side1)&&isborder(r, side2))
        putmarker!(r)
        if !isborder(r, side1)
            move!(r, side1)
        end
        if !isborder(r, side1)
            move!(r, side2)
        end
    end
end
"""
checks_whether_to_put_a_marker_in_the_cell_and_go_back(r, side) - это функция действует так пока переменная try_ch_side
не примет значение false она будет проверять поставить ли маркер в клетке, пока робот не будет находится
в клетке где стоит маркер она и возвращается назад и спускается вниз если пройдет проверку
"""
function checks_whether_to_put_a_marker_in_the_cell_and_go_back(r, side)
    try_ch_side = true
    while try_ch_side == true
        while !ismarker(r)
            putmarker!(r)
            move!(r, side)
        end
        while (ismarker(r) == true)&&(!isborder(r, inverse(side)))
            move!(r, inverse(side))
        end
        if (!isborder(r, Sud))
            move!(r, Sud)
        else
            try_ch_side = false
        end
    end
end
"""
inverse(side::HorizonSide) - получает переданую в функцию side преобразовывает её в Int значение то есть Nord = 0,
West = 1, Sud = 2, Ost = 3 потом преобразовывает в сторону противоположную : допустим, что в функцию передано напрвление West после преобразование в значения Int
West = 1 после к этого к этому числу прибавляется 2 и берется остаток от деления на число 4 то есть 3 и это сторона Ost, а эта сторона противоположна стороне West
"""
function inverse(side)
    HorizonSide(mod(Int(side)+2, 4))
end
"""
exit_to_the_lower_left_cell!(r, arr) - это функция заполняет в переданый массив  двумя сторонами West || Sud: принцип алгоритма такой
цикл while работет так: пока слева(West) && снизу(Sud) нет перегоротки каждая проверка будет выполнятся проверка двух if: первый проверяет
есть ли перегоротка в заданном направлении, если нет, то робот делает шаг в сторону West и заносит это сторону в переданный массив
принцип работы второго if такой же.Цикл while прекратит свою работу, когда окажется в левом нижнем углу
"""
function exit_to_the_lower_left_cell!(r, arr)
    while !(isborder(r, Ost)&&isborder(r, Sud))
        if !isborder(r, Ost)
            move!(r, Ost)
            push!(arr, Ost)
        end
        if !isborder(r, Sud)
            move!(r, Sud)
            push!(arr, Sud)
        end
    end
end
"""
return_to_starting_point!(r, arr) - это функция использует переданный массив и передвигается по элементам массива то есть пока длина массива
не будет равна нулю выполняется следующие переменная side принимает значение последнего элемента массива и двигается в противоположном направление за это отвечает фунrция inverse
дальше удаляется последний элемент массива и так пока длина массива не станет равна 0
"""
function return_to_starting_point!(r, arr)
    while length(arr)>0
        side=pop!(arr)
        move!(r, inverse(side))
    end
end


text:
const global_style = "position: absolute; top: 0px; left: 0px; visibility: visible; will-change: transform; opacity: 1; transform: translate(525px, 1220px) scale(1); transition-duration: 250ms; transition-timing-function: cubic-bezier(0.4, 0, 0.2, 1); transition-property: transform, opacity;"

struct MediaItem
  title::String
  media::String
  preview::String
  link::String
  tags::Vector{String}
end

function master_url(root, path)
    urlbase = "https://simondanisch.github.io/ReferenceImages/gallery/"
    if Sys.iswindows()
        path = replace(path, "\\" => "/")
        root = replace(root, "\\" => "/")
    end
    url = replace(
        path,
        root => urlbase
    )
end

function MediaItem(path, example)
    media_file = ""
    media_folder = joinpath(path, "media")
    files = readdir(media_folder)
    for ext in ("mp4", "gif", "jpg", "png")
        idx = findfirst(x-> endswith(x, ext), files)
        if idx !== nothing
            media_file = joinpath(media_folder, files[idx])
            break
        end
    end
    if isempty(media_file)
        error("No media for $(example.title)")
    end
    root = dirname(path)
    MediaItem(
        example.title,
        embed_media(master_url(root, media_file)),
        "",
        master_url(root, joinpath(path, "index.html")),
        string.(example.tags)
    )
end

function create_item(item::MediaItem; style = global_style)
    create_item(item.title, item.media, item.preview, item.link, item.tags; style = style)
end

function create_item(title, media, preview, link, tags; style = global_style)
  """

          $(media)

$(repr(title))

  """
end

function create_page(
    items, tags;
    shufflejs = "https://cdnjs.cloudflare.com/ajax/libs/Shuffle/5.2.0/shuffle.min.js"
  )
  tag_buttons = join(map(tags) do tag
    """<button class="btn btn--primary" data-group=$(repr(tag))>$tag</button>"""
  end, "\n")
  groups = join(create_item.(items), "\n")
  html_groups = """

  <link href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" rel="stylesheet" />
  <link rel="stylesheet" href="https://rawgit.com/Vestride/Shuffle/master/docs/css/prism.css">
  <link rel="stylesheet" href="https://rawgit.com/Vestride/Shuffle/master/docs/css/normalize.css">
  <link rel="stylesheet" href="https://rawgit.com/Vestride/Shuffle/master/docs/css/style.css">
  <link rel="stylesheet" href="https://rawgit.com/Vestride/Shuffle/master/docs/css/shuffle-styles.css">

Example

            Search
            <input class="textfield filter__search js-shuffle-search" type="search" id="filters-search-input">

Filter

              $(tag_buttons)

          <fieldset class="filters-group">
            <legend class="filter-label">Sort</legend>

                <input type="radio" name="sort-value" value="dom"> Default

                <input type="radio" name="sort-value" value="title"> Title

                <input type="radio" name="sort-value" value="date-created"> Date Created

          </fieldset>

        $groups

  """
  js_script = read(joinpath(@__DIR__, "gallery.js"), String)
  return """
  <!doctype html>

      <meta charset="UTF-8">

      $html_groups

  """
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcomment

"""
    'feffery'_antdcomment(;kwargs...)
    'feffery'_antdcomment(children::Any;kwargs...)
    'feffery'_antdcomment(children_maker::Function;kwargs...)

An AntdComment component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `action` (a value equal to: 'liked', 'disliked', 'default'; optional)
- `authorName` (String; optional)
- `authorNameHref` (String; optional)
- `avatarProps` (Dict; optional)
- `className` (String; optional)
- `commentContent` (String; optional)
- `commentId` (String; optional)
- `defaultAction` (a value equal to: 'liked', 'disliked', 'default'; optional)
- `deleteClicks` (Real; optional)
- `dislikesCount` (Real; optional)
- `fromNow` (Bool; optional)
- `likesCount` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-cn', 'en-us'; optional)
- `publishTime` (required): . publishTime has the following type: lists containing elements 'value', 'format'.
Those elements have the following types:
  - `value` (String; required)
  - `format` (String; optional)
- `replyClicks` (Real; optional)
- `showDelete` (Bool; optional)
- `showLikeDislike` (Bool; optional)
- `showReply` (Bool; optional)
- `style` (Dict; optional)
"""
function 'feffery'_antdcomment(; kwargs...)
        available_props = Symbol[:children, :id, :action, :authorName, :authorNameHref, :avatarProps, :className, :commentContent, :commentId, :defaultAction, :deleteClicks, :dislikesCount, :fromNow, :likesCount, :loading_state, :locale, :publishTime, :replyClicks, :showDelete, :showLikeDislike, :showReply, :style]
        wild_props = Symbol[]
        return Component("'feffery'_antdcomment", "AntdComment", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdcomment(children::Any; kwargs...) = 'feffery'_antdcomment(;kwargs..., children = children)
'feffery'_antdcomment(children_maker::Function; kwargs...) = 'feffery'_antdcomment(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
"""
    DeterministicEquivalent

Deterministic equivalent memory structure. Stochastic program is stored as one large optimization problem. Supported by any standard `AbstractOptimizer`.

"""
struct DeterministicEquivalent{N, M, S <: NTuple{M, Scenarios}} <: AbstractStochasticStructure{N}
    decisions::NTuple{N, Decisions}
    scenarios::S
    sub_objectives::NTuple{N, Vector{Tuple{MOI.OptimizationSense, MOI.AbstractScalarFunction}}}
    model::JuMP.Model
    variable_map::Dict{Tuple{VI, Int}, VI}
    constraint_map::Dict{Tuple{CI, Int}, CI}

    function DeterministicEquivalent(decisions::NTuple{N, Decisions}, scenarios::NTuple{M, Scenarios}) where {N, M}
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        sub_objectives = ntuple(Val(N)) do i
            Vector{Tuple{MOI.ObjectiveSense, MOI.AbstractScalarFunction}}()
        end
        S = typeof(scenarios)
        return new{N,M,S}(decisions, scenarios, sub_objectives, Model(), Dict{Tuple{VI,Int}, VI}(), Dict{Tuple{CI,Int}, CI}())
    end
end

function StochasticStructure(decisions::NTuple{N, Decisions}, scenario_types::ScenarioTypes{M}, ::Deterministic) where {N, M}
    scenarios = ntuple(Val(M)) do i
        Vector{scenario_types[i]}()
    end
    return DeterministicEquivalent(decisions, scenarios)
end

function StochasticStructure(decisions::NTuple{N, Decisions}, scenarios::NTuple{M, Vector{<:AbstractScenario}}, ::Deterministic) where {N, M}
    return DeterministicEquivalent(decisions, scenarios)
end

# Base overloads #
# ========================== #
function Base.print(io::IO, structure::DeterministicEquivalent)
    print(io, "Deterministic equivalent problem\n")
    print(io, structure.model)
end
# ========================== #

# MOI #
# ========================== #
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractModelAttribute)
    return MOI.get(backend(structure.model), attr)
end
function MOI.get(structure::DeterministicEquivalent, attr::ScenarioDependentModelAttribute)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    if attr.attr isa MOI.ObjectiveFunction
        return structure.sub_objectives[attr.stage][attr.scenario_index][2]
    elseif attr.attr isa MOI.ObjectiveFunctionType
        return typeof(structure.sub_objectives[attr.stage][attr.scenario_index][2])
    elseif attr.attr isa MOI.ObjectiveSense
        return structure.sub_objectives[attr.stage][attr.scenario_index][1]
    elseif attr.attr isa MOI.ObjectiveValue || attr.attr isa MOI.DualObjectiveValue
        return MOIU.eval_variables(structure.sub_objectives[attr.stage][attr.scenario_index][2]) do idx
            return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
        end
    else
        # Most attributes are shared with the deterministic equivalent
        return MOI.get(backend(structure.model), attr.attr)
    end
end
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(backend(structure.model), attr, index)
end
function MOI.get(structure::DeterministicEquivalent, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, attr.scenario_index)]
    return MOI.get(backend(structure.model), attr.attr, mapped_index)
end
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractConstraintAttribute, cindex::MOI.ConstraintIndex)
    return MOI.get(backend(structure.model), attr, cindex)
end
function MOI.get(structure::DeterministicEquivalent, attr::ScenarioDependentConstraintAttribute, ci::CI{F,S}) where {F,S}
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_ci = structure.constraint_map[(ci, attr.scenario_index)]
    return MOI.get(backend(structure.model), attr.attr, mapped_ci)
end
function MOI.get(structure::DeterministicEquivalent, attr::ScenarioDependentConstraintAttribute, ci::CI{F,S}) where {F <: Union{MOI.SingleVariable, SingleDecision}, S}
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_ci = CI{F,S}(structure.variable_map[(MOI.VariableIndex(ci.value), attr.scenario_index)].value)
    return MOI.get(backend(structure.model), attr.attr, mapped_ci)
end

function MOI.set(structure::DeterministicEquivalent, attr::MOI.Silent, flag)
    MOI.set(backend(structure.model), attr, flag)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractModelAttribute, value)
    if attr isa MOI.ObjectiveFunction
        # Get full objective+sense
        dep_obj = value
        obj_sense = objective_sense(structure.model)
        # Update first-stage objective
        structure.sub_objectives[1][1] = (obj_sense, value)
        # Update main objective
        for (i, sub_objective) in enumerate(structure.sub_objectives[2])
            (sense, func) = sub_objective
            if obj_sense == sense
                dep_obj += probability(structure, 2, i) * func
            else
                dep_obj -= probability(structure, 2, i) * func
            end
        end
        MOI.set(backend(structure.model), attr, dep_obj)
    elseif attr isa MOI.ObjectiveSense
        # Get full objective+sense
        dep_obj, prev_sense = structure.sub_objective[1][1]
        # Update first-stage objective
        structure.sub_objectives[1][1] = (value, dep_obj)
        # Update main objective (if necessary)
        if value != prev_sense
            for (i, sub_objective) in enumerate(structure.sub_objectives[2])
                (sense, func) = sub_objective
                if value == sense
                    dep_obj += probability(structure, 2, i) * func
                else
                    dep_obj -= probability(structure, 2, i) * func
                end
            end
            MOI.set(backend(structure.model), MOI.ObjectiveFunction{typeof(dep_obj)}(), dep_obj)
        end
        MOI.set(backend(structure.model), MOI.ObjectiveSense(), value)
    else
        # Most attributes are shared with the deterministic equivalent
        MOI.set(backend(structure.model), attr, value)
    end
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::ScenarioDependentModelAttribute, value)
    if attr.attr isa MOI.ObjectiveFunction
        # Get full objective+sense
        obj_sense = objective_sense(structure.model)
        dep_obj = objective_function(structure.model)
        # Update subobjective
        (sub_sense, prev_func) = structure.sub_objectives[attr.stage][attr.scenario_index]
        structure.sub_objectives[attr.stage][attr.scenario_index] = (sub_sense, value)
        prev_sub_obj = jump_function(structure.model, prev_func)
        sub_obj = jump_function(structure.model, value)
        # Update main objective
        if obj_sense == sub_sense
            dep_obj += probability(structure, attr.stage, attr.scenario_index) * (sub_obj - prev_sub_obj)
        else
            dep_obj -= probability(structure, attr.stage, attr.scenario_index) * (sub_obj - prev_sub_obj)
        end
        MOI.set(backend(structure.model), attr.attr, moi_function(dep_obj))
    elseif attr.attr isa MOI.ObjectiveSense
        # Get full objective+sense
        obj_sense = objective_sense(structure.model)
        dep_obj = moi_function(objective_function(structure.model))
        # Update subobjective sense
        (prev_sense, func) = structure.sub_objectives[attr.stage][attr.scenario_index]
        structure.sub_objectives[attr.stage][attr.scenario_index] = (value, func)
        # Update main objective (if necessary)
        if value != prev_sense
            if value == obj_sense
                dep_obj += 2 * probability(structure, attr.stage, attr.scenario_index) * func
            else
                dep_obj -= 2 * probability(structure, attr.stage, attr.scenario_index) * func
            end
            MOI.set(backend(structure.model), MOI.ObjectiveFunction{typeof(dep_obj)}(), dep_obj)
        end
    else
        # Most attributes are shared with the deterministic equivalent
        MOI.set(backend(structure.model), attr.attr, value)
    end
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(backend(structure.model), attr, index, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::ScenarioDependentVariableAttribute,
                 index::MOI.VariableIndex, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, attr.scenario_index)]
    MOI.set(backend(structure.model), attr.attr, mapped_index, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractConstraintAttribute,
                 cindex::MOI.ConstraintIndex, value)
    MOI.set(backend(structure.model), attr, cindex, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::ScenarioDependentConstraintAttribute,
                 cindex::CI{F,S}, value) where {F, S}
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_index = structure.constraint_map[(cindex, attr.scenario_index)]
    MOI.set(backend(structure.model), attr.attr, mapped_index, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::ScenarioDependentConstraintAttribute,
                 cindex::CI{F,S}, value) where {F <: SingleDecision, S}
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    mapped_index = CI{F,S}(structure.variable_map[(MOI.VariableIndex(cindex.value), attr.scenario_index)].value)
    MOI.set(backend(structure.model), attr.attr, mapped_index, value)
    return nothing
end

function MOI.is_valid(structure::DeterministicEquivalent, index::MOI.VariableIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.model), index)
end
function MOI.is_valid(structure::DeterministicEquivalent, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer)
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, scenario_index)]
    return MOI.is_valid(backend(structure.model), mapped_index)
end

function MOI.is_valid(structure::DeterministicEquivalent, ci::MOI.ConstraintIndex{F,S}, stage::Integer) where {F, S}
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.model), ci)
end
function MOI.is_valid(structure::DeterministicEquivalent, ci::MOI.ConstraintIndex{F,S}, stage::Integer, scenario_index::Integer) where {F, S}
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_ci = structure.constraint_map[(ci, scenario_index)]
    return MOI.is_valid(backend(structure.model), mapped_ci)
end
function MOI.is_valid(structure::DeterministicEquivalent, ci::MOI.ConstraintIndex{F,S}, stage::Integer, scenario_index::Integer) where {F <: SingleDecision, S}
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_ci = CI{F,S}(structure.variable_map[(MOI.VariableIndex(ci.value), scenario_index)].value)
    return MOI.is_valid(backend(structure.model), mapped_ci)
end

function MOI.add_constraint(structure::DeterministicEquivalent, f::SingleDecision, s::MOI.AbstractSet)
    return MOI.add_constraint(backend(structure.model), f, s)
end
function MOI.add_constraint(structure::DeterministicEquivalent, f::SingleDecision, s::MOI.AbstractSet, stage::Integer, scenario_index::Integer)
    stage == 1 && error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    g = MOIU.map_indices(f) do index
        return structure.variable_map[(index, scenario_index)]
    end
    return MOI.add_constraint(backend(structure.model), g, s)
end

function MOI.delete(structure::DeterministicEquivalent, index::MOI.VariableIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    JuMP.delete(structure.model, DecisionRef(structure.model, index))
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, scenario_index)]
    JuMP.delete(structure.model, DecisionRef(structure.model, mapped_index))
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent, indices::Vector{MOI.VariableIndex}, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    JuMP.delete(structure.model, DecisionRef.(structure.model, indices))
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, indices::Vector{MOI.VariableIndex}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_indices = map(indices) do index
        structure.variable_map[(index, scenario_index)]
    end
    JuMP.delete(structure.model, DecisionRef.(structure.model, mapped_indices))
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent, ci::MOI.ConstraintIndex{F,S}, stage::Integer) where {F, S}
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.model), ci)
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent, cis::Vector{MOI.ConstraintIndex{F,S}}, stage::Integer) where {F, S}
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.model), cis)
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, ci::MOI.ConstraintIndex{F,S}, stage::Integer, scenario_index::Integer) where {F, S} where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_ci = structure.constraint_map[(ci, scenario_index)]
    MOI.delete(backend(structure.model), mapped_ci)
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, ci::MOI.ConstraintIndex{F,S}, stage::Integer, scenario_index::Integer) where {F <: SingleDecision, S} where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_ci = CI{F,S}(structure.variable_map[(MOI.VariableIndex(ci.value), scenario_index)].value)
    MOI.delete(backend(structure.model), mapped_ci)
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, cis::Vector{MOI.ConstraintIndex{F,S}}, stage::Integer, scenario_index::Integer) where {N, F, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_cis = map(cis) do ci
        structure.constraint_map[(ci, scenario_index)]
    end
    MOI.delete(backend(structure.model), mapped_cis)
    return nothing
end
function MOI.delete(structure::DeterministicEquivalent{N}, cis::Vector{MOI.ConstraintIndex{F,S}}, stage::Integer, scenario_index::Integer) where {N, F <: SingleDecision, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_cis = map(cis) do ci
        CI{F,S}(structure.variable_map[(MOI.VariableIndex(ci.value), scenario_index)].value)
    end
    MOI.delete(backend(structure.model), mapped_cis)
    return nothing
end

# JuMP #
# ========================== #
function JuMP.fix(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer, val::Number) where N
    dref = DecisionRef(structure.model, index)
    fix(dref, val)
    return nothing
end
function JuMP.fix(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer, val::Number) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    # Fix mapped decision
    mapped_index = structure.variable_map[(index, scenario_index)]
    dref = DecisionRef(structure.model, mapped_index)
    fix(dref, val)
    return nothing
end
function JuMP.unfix(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer) where N
    dref = DecisionRef(structure.model, index)
    unfix(dref)
    return nothing
end
function JuMP.unfix(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    # Unfix mapped decision
    mapped_index = structure.variable_map[(index, scenario_index)]
    dref = DecisionRef(structure.model, mapped_index)
    unfix(dref)
    return nothing
end

function JuMP.set_objective_sense(structure::DeterministicEquivalent, stage::Integer, sense::MOI.OptimizationSense)
    if stage == 1
        # Changes the first-stage sense modifies the whole objective as usual
        MOI.set(structure, MOI.ObjectiveSense(), sense)
    else
        # Every sub-objective in the given stage should be changes
        for scenario_index in 1:num_scenarios(structure, stage)
            # Use temporary model to apply modification
            attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.ObjectiveSense())
            MOI.set(structure, attr, sense)
        end
    end
    return nothing
end

function JuMP.objective_function_type(structure::DeterministicEquivalent)
    return jump_function_type(structure.model,
                              MOI.get(structure, MOI.ObjectiveFunctionType()))
end
function JuMP.objective_function_type(structure::DeterministicEquivalent{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.ObjectiveFunctionType())
    return jump_function_type(structure.model,
                              MOI.get(structure, attr))
end

function JuMP.objective_function(structure::DeterministicEquivalent, FunType::Type{<:AbstractJuMPScalar})
    MOIFunType = moi_function_type(FunType)
    func = MOI.get(structure,
                   MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
    return jump_function(structure.model, func)
end
function JuMP.objective_function(structure::DeterministicEquivalent, proxy::JuMP.Model, stage::Integer, FunType::Type{<:AbstractJuMPScalar})
    if stage == 1
        obj::FunType = jump_function(structure.model, structure.sub_objectives[stage][1][2])
        return obj
    else
        return objective_function(proxy, FunType)
    end
end
function JuMP.objective_function(structure::DeterministicEquivalent{N},
                                 proxy::JuMP.Model,
                                 stage::Integer,
                                 scenario_index::Integer,
                                 FunType::Type{<:AbstractJuMPScalar}) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOIFunType = moi_function_type(FunType)
    attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.ObjectiveFunction{MOIFunType}())
    func = MOI.get(structure, attr)::MOIFunType
    return jump_function(structure.model, func)
end

function JuMP._moi_optimizer_index(structure::DeterministicEquivalent, index::VI)
    return decision_index(backend(structure.model), index)
end
function JuMP._moi_optimizer_index(structure::DeterministicEquivalent, index::VI, scenario_index::Integer)
    mapped_index = structure.variable_map[(index, scenario_index)]
    return decision_index(backend(structure.model), mapped_index)
end
function JuMP._moi_optimizer_index(structure::DeterministicEquivalent, ci::CI)
    return decision_index(backend(structure.model), ci)
end
function JuMP._moi_optimizer_index(structure::DeterministicEquivalent, ci::CI, scenario_index::Integer)
    mapped_ci = structure.constraint_map[(ci, scenario_index)]
    return decision_index(backend(structure.model), mapped_ci)
end
function JuMP._moi_optimizer_index(structure::DeterministicEquivalent, ci::CI{F,S}, scenario_index::Integer) where {F <: SingleDecision, S}
    mapped_ci = CI{F,S}(structure.variable_map[(MOI.VariableIndex(ci.value), scenario_index)].value)
    return decision_index(backend(structure.model), mapped_ci)
end

function JuMP.set_objective_coefficient(structure::DeterministicEquivalent{N}, index::VI, var_stage::Integer, stage::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    var_stage <= stage || error("Can only modify coefficient in current stage of decision or subsequent stages from where decision is taken.")
    if var_stage == 1 && stage == 1
        # Use temporary model to apply modification
        obj = structure.sub_objectives[1][1][2]
        m = Model()
        MOI.set(backend(m), MOI.ObjectiveFunction{typeof(obj)}(), obj)
        dref = DecisionRef(m, index)
        set_objective_coefficient(m, dref, coeff)
        obj = objective_function(m)
        F = moi_function_type(typeof(obj))
        # Modify full objective
        MOI.set(structure, MOI.ObjectiveFunction{F}(), moi_function(obj))
    elseif (var_stage == 1 && stage > 1) || var_stage > 1
        for scenario_index in 1:num_scenarios(structure, stage)
            # Use temporary model to apply modification
            obj = structure.sub_objectives[stage][scenario_index][2]
            m = Model()
            F = moi_function_type(typeof(obj))
            MOI.set(backend(m), MOI.ObjectiveFunction{F}(), obj)
            dref = DecisionRef(m, index)
            set_objective_coefficient(m, dref, coeff)
            obj = objective_function(m)
            attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.ObjectiveFunction{typeof(obj)}())
            MOI.set(structure, attr, obj)
        end
    end
    dref = DecisionRef(structure.model, index)
    set_objective_coefficient(structure.model, dref, coeff)
    return nothing
end
function JuMP.set_objective_coefficient(structure::DeterministicEquivalent{N}, index::VI, var_stage::Integer, stage::Integer, scenario_index::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    if var_stage == 1
        # Use temporary model to apply modification
        obj = structure.sub_objectives[stage][scenario_index][2]
        m = Model()
        MOI.set(backend(m), MOI.ObjectiveFunction{typeof(obj)}(), obj)
        dref = DecisionRef(m, index)
        set_objective_coefficient(m, dref, coeff)
        obj = objective_function(m)
        F = moi_function_type(typeof(obj))
        attr = ScenarioDependentModelAttribute(stage, scenario_index, MOI.ObjectiveFunction{F}())
        MOI.set(structure, attr, moi_function(obj))
    else
        # Use temporary model to apply modification
        sense, obj = structure.sub_objectives[stage][scenario_index]
        m = Model()
        F = typeof(obj)
        MOI.set(backend(m), MOI.ObjectiveFunction{F}(), obj)
        mapped_index = structure.variable_map[(index, scenario_index)]
        dref = DecisionRef(m, mapped_index)
        set_objective_coefficient(m, dref, coeff)
        obj = objective_function(m)
        structure.sub_objectives[stage][scenario_index] = (sense, moi_function(obj))
        # Set coefficient of mapped second-stage variable
        dref = DecisionRef(structure.model, mapped_index)
        set_objective_coefficient(structure.model, dref, probability(structure, stage, scenario_index) * coeff)
    end
    return nothing
end

function JuMP.set_normalized_coefficient(structure::DeterministicEquivalent,
                                         ci::CI{F,S},
                                         index::VI,
                                         value) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    MOI.modify(backend(structure.model), ci,
               DecisionCoefficientChange(index, convert(T, value)))
    return nothing
end
function JuMP.set_normalized_coefficient(structure::DeterministicEquivalent{N},
                                         ci::CI{F,S},
                                         index::VI,
                                         stage::Integer,
                                         scenario_index::Integer,
                                         value) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, scenario_index)]
    mapped_ci = structure.constraint_map[(ci, scenario_index)]
    MOI.modify(backend(structure.model), mapped_ci,
               DecisionCoefficientChange(mapped_index, convert(T, value)))
    return nothing
end

function JuMP.normalized_coefficient(structure::DeterministicEquivalent{N},
                                     ci::CI{F,S},
                                     index::VI,
                                     stage::Integer,
                                     scenario_index::Integer) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_ci = structure.constraint_map[(ci, scenario_index)]
    f = MOI.get(structure, MOI.ConstraintFunction(), mapped_ci)::F
    mapped_index = structure.variable_map[(index, scenario_index)]
    dref = DecisionRef(structure.model, mapped_index)
    return JuMP._affine_coefficient(jump_function(structure.model, f), dref)
end

function JuMP.set_normalized_rhs(structure::DeterministicEquivalent,
                                 ci::CI{F,S},
                                 value) where {T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    MOI.set(backend(structure.model), MOI.ConstraintSet(), ci,
            S(convert(T, value)))
    return nothing
end
function JuMP.set_normalized_rhs(structure::DeterministicEquivalent{N},
                                 ci::CI{F,S},
                                 stage::Integer,
                                 scenario_index::Integer,
                                 value) where {N,
                                               T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_index = structure.constraint_map[(ci, scenario_index)]
    MOI.set(backend(structure.model), MOI.ConstraintSet(), mapped_index,
            S(convert(T, value)))
    return nothing
end

function DecisionRef(proxy::JuMP.Model, structure::DeterministicEquivalent, index::VI)
    return DecisionRef(structure.model, index)
end
function DecisionRef(proxy::JuMP.Model, structure::DeterministicEquivalent, index::VI, scenario_index::Integer)
    mapped_index = structure.variable_map[(index, scenario_index)]
    return DecisionRef(structure.model, mapped_index)
end
function DecisionRef(proxy::JuMP.Model, structure::DeterministicEquivalent{N}, index::VI, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return DecisionRef(structure.model, index)
end

function JuMP.jump_function(structure::DeterministicEquivalent{N},
                            proxy::JuMP.Model,
                            stage::Integer,
                            f::MOI.AbstractFunction) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return JuMP.jump_function(structure.model, f)
end

# Getters #
# ========================== #
function structure_name(structure::DeterministicEquivalent)
    return "Deterministic equivalent"
end
function scenario_types(structure::DeterministicEquivalent{N}) where N
    return ntuple(Val{N-1}()) do i
        eltype(structure.scenarios[i])
    end
end
function decision(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    return decision(structure.decisions[stage], index)
end
function decision(structure::DeterministicEquivalent{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    mapped_index = structure.variable_map[(index, scenario_index)]
    return decision(structure.decisions[stage], mapped_index)
end
function scenario(structure::DeterministicEquivalent{N}, stage::Integer, scenario_index::Integer) where N
    return structure.scenarios[stage-1][scenario_index]
end
function scenarios(structure::DeterministicEquivalent{N}, stage::Integer) where N
    return structure.scenarios[stage-1]
end
function subproblem(structure::DeterministicEquivalent, stage::Integer, scenario_index::Integer)
    error("The determinstic equivalent is not decomposed into subproblems.")
end
function subproblems(structure::DeterministicEquivalent, stage::Integer)
    error("The determinstic equivalent is not decomposed into subproblems.")
end
function num_subproblems(structure::DeterministicEquivalent, stage::Integer)
    return 0
end
function deferred(structure::DeterministicEquivalent)
    return num_variables(structure.model) == 0
end
# ========================== #

# Setters
# ========================== #
function update_known_decisions!(structure::DeterministicEquivalent)
    update_known_decisions!(structure.model)
    return nothing
end
function update_known_decisions!(structure::DeterministicEquivalent, stage::Integer, scenario_index::Integer)
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    update_known_decisions!(structure.model)
    return nothing
end

function add_scenario!(structure::DeterministicEquivalent, stage::Integer, scenario::AbstractScenario)
    push!(scenarios(structure, stage), scenario)
    return nothing
end
function add_worker_scenario!(structure::DeterministicEquivalent, stage::Integer, scenario::AbstractScenario, w::Integer)
    add_scenario!(structure, scenario, stage)
    return nothing
end
function add_scenario!(scenariogenerator::Function, structure::DeterministicEquivalent, stage::Integer)
    add_scenario!(structure, stage, scenariogenerator())
    return nothing
end
function add_worker_scenario!(scenariogenerator::Function, structure::DeterministicEquivalent, stage::Integer, w::Integer)
    add_scenario!(scenariogenerator, structure, stage)
    return nothing
end
function add_scenarios!(structure::DeterministicEquivalent, stage::Integer, _scenarios::Vector{<:AbstractScenario})
    append!(scenarios(structure, stage), _scenarios)
    return nothing
end
function add_worker_scenarios!(structure::DeterministicEquivalent, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)
    add_scenarios!(structure, scenarios, stasge)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, structure::DeterministicEquivalent, stage::Integer, n::Integer)
    for i = 1:n
        add_scenario!(structure, stage) do
            return scenariogenerator()
        end
    end
    return nothing
end
function add_worker_scenarios!(scenariogenerator::Function, structure::DeterministicEquivalent, stage::Integer, n::Integer, w::Integer)
    add_scenarios!(scenariogenerator, structure, n, stage)
    return nothing
end
function sample!(structure::DeterministicEquivalent, stage::Integer, sampler::AbstractSampler, n::Integer)
    sample!(scenarios(structure, stage), sampler, n)
    return nothing
end
# ========================== #


text:
struct Disturbances
	W
	w
end

function disturbances(W)
	T = length(W)
	d = size(W[1], 1)
	w = []
	for t = 1:T
		_w = sqrt(W[t])
		tmp = []
		for i = 1:d
			push!(tmp, _w[:, i])
			push!(tmp, -1.0 * _w[:, i])
		end

		push!(w, tmp)
	end

	return Disturbances(W, w)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
module SurrogateModelOptim

@doc read(joinpath(dirname(@__DIR__), "README.md"), String) SurrogateModelOptim

export  smoptimize,
        surrogate_model,
        model_infill,
        scaled_LHC_sampling_plan,
        best_fitness,
        best_candidate,
        f_calls

using LatinHypercubeSampling
import LatinHypercubeSampling.Continuous, LatinHypercubeSampling.Categorical
using ScatteredInterpolation
using BlackBoxOptim
import BlackBoxOptim: best_candidate, best_fitness
using Parameters
using Distances
using Statistics
using StatsBase
using LinearAlgebra
using Distributed
using StaticArrays
import NearestNeighbors: KDTree, knn
using Printf
using TimerOutputs

include("types.jl")
include("LHC_sampling_plan.jl")
include("model_infill_utilities.jl")
include("model_infill.jl")
include("surrogate_model_utilities.jl")
include("surrogate_model.jl")
include("smoptimize_utilities.jl")
include("smoptimize.jl")

end # module


text:
@testset "Movement models" begin
    @test eqdist_ap isa PopState
    @test eqdist_ap0.P ≈ eqdist_eigvecs
    @test all(eqdist_ap.P .> 0)
    @test all(sum(move.M; dims = 1) .≈ 1)
    @test all(0 .≤ eqdist_ap0.P .≤ 1)
    @test sum(eqdist_ap0) ≈ 1.0
    @test sum(eqdist_ap) ≈ 100.0
end