text:
@testset "Double Description" begin
    @test Polyhedra.polytypefor(Float32) == Float64
    @testset "H-representation -> V-representation" begin
        @testset "Intersection" begin
            @testset "Vector" begin
                @testset "Exact" begin
                    h = hrep([HalfSpace([ 1,  1], 1),
                              HalfSpace([ 1, -1], 0),
                              HalfSpace([-1,  0], 0)])
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Rational{BigInt},Vector{Rational{BigInt}}}
                    @test collect(points(v)) == [[1//2, 1//2], [0, 0], [0, 1]]
                    @test !hasallrays(v)
                end
                @testset "Numerical" begin
                    h = hrep([HalfSpace([ 1.,  1], 1),
                              HalfSpace([ 1., -1], 0),
                              HalfSpace([-1.,  0], 0)])
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Float64,Vector{Float64}}
                    @test collect(points(v)) == [[1/2, 1/2], [0.0, 0.0], [0.0, 1.0]]
                    @test !hasallrays(v)
                end
            end
            @testset "SVector" begin
                @testset "Exact" begin
                    h = hrep([HalfSpace((@SVector [ 1,  1]), 1),
                              HalfSpace((@SVector [ 1, -1]), 0),
                              HalfSpace((@SVector [-1,  0]), 0)])
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Rational{BigInt},SVector{2,Rational{BigInt}}}
                    @test collect(points(v)) == [(@SVector [1//2, 1//2]), (@SVector [0, 0]), (@SVector [0, 1])]
                    @test !hasallrays(v)
                end
                @testset "Numerical" begin
                    h = hrep([HalfSpace((@SVector [ 1.,  1]), 1),
                              HalfSpace((@SVector [ 1., -1]), 0),
                              HalfSpace((@SVector [-1.,  0]), 0)])
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Float64,SVector{2,Float64}}
                    @test collect(points(v)) == [(@SVector [1/2, 1/2]), (@SVector [0.0, 0.0]), (@SVector [0.0, 1.0])]
                    @test !hasallrays(v)
                end
            end
        end
        @testset "MixedMatHRep" begin
            @testset "Exact" begin
                h = hrep([ 1  1
                           1 -1
                          -1  0],
                         [1, 0, 0])
                #v = @inferred doubledescription(h)
                v = doubledescription(h)
                @test v isa Polyhedra.MixedMatVRep{Rational{BigInt}}
                @test v.V == [1//2 1//2; 0//1 0//1; 0//1 1//1]
                @test v.R == zeros(Rational{BigInt}, 0, 2)
                @test isempty(v.Rlinset)
            end
            @testset "Numerical" begin
                h = hrep([ 1  1
                           1 -1
                          -1  0],
                         [1., 0, 0])
                #v = @inferred doubledescription(h)
                v = doubledescription(h)
                @test v isa Polyhedra.MixedMatVRep{Float64}
                @test v.V == [1/2 1/2; 0 0; 0 1]
                @test v.R == zeros(0, 2)
                @test isempty(v.Rlinset)
            end
        end
        @testset "Empty Intersection" begin
            h = HyperPlane([0, 1], 1) ∩ HyperPlane([0, 1], -1)
            @testset "0x_1 + 0x_2 = 1" begin
                @testset "Exact" begin
                    h = intersect(HyperPlane([0, 0], 1))
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Rational{BigInt},Vector{Rational{BigInt}}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
                @testset "Numerical" begin
                    h = intersect(HyperPlane([0., 0], 1))
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Float64,Vector{Float64}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
            end
            @testset "-1 = 0x_1 + x_2 = 1" begin
                @testset "Exact" begin
                    h = HyperPlane([0, 1], 1) ∩ HyperPlane([0, 1], -1)
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Rational{BigInt},Vector{Rational{BigInt}}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
                @testset "Numerical" begin
                    h = HyperPlane([0, 1.], 1) ∩ HyperPlane([0, 1.], -1)
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Float64,Vector{Float64}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
            end
            @testset "1 ≤ x_1 + 0x_2 ≤ -1" begin
                @testset "Exact" begin
                    h = HalfSpace([1, 0], -1) ∩ HalfSpace([-1, 0], -1)
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Rational{BigInt},Vector{Rational{BigInt}}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
                @testset "Numerical" begin
                    h = HalfSpace([1, 0.], -1) ∩ HalfSpace([-1, 0.], -1)
                    #v = @inferred doubledescription(h)
                    v = doubledescription(h)
                    @test v isa Polyhedra.Hull{Float64,Vector{Float64}}
                    @test !haspoints(v)
                    @test !hasallrays(v)
                end
            end
        end
    end
    @testset "V-representation -> H-representation" begin
        @testset "Hull" begin
            @testset "Vector" begin
                @testset "Exact" begin
                    v = conichull([1, 0],
                                  [0, 1])
                    h = doubledescription(v)
                    #h = @inferred doubledescription(v)
                    @test !hashyperplanes(h)
                    @test collect(halfspaces(h)) == [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0), HalfSpace([0, 0], 1)] # FIXME get rid of (0, 0) 1
                end
                @testset "Numerical" begin
                    v = conichull([1., 0.],
                                  [0., 1.])
                    h = doubledescription(v)
                    #h = @inferred doubledescription(v)
                    @test !hashyperplanes(h)
                    @test collect(halfspaces(h)) == [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0), HalfSpace([0, 0], 1)] # FIXME get rid of (0, 0) 1
                end
            end
        end
    end
end


text:
"""
  `zminsv(A::Array, tol) -> sval, svec`

  zminsv calculates the smallest singular value
  of matrix A and the associated right singular vector
  via an implicit inverse iteration in the form of the
  Gauss-Newton iteration

    input  A --- the matrix
           tol --- the error tolerence

   output  s --- the smallest singular value
           x --- the associated right singular vector
"""
function zminsv(A::AbstractArray{S,2}, tol::U) where {S<:Number,U<:AbstractFloat}
    E = one(S)
    m, n = size(A)           	# get the dimensions of A
    if m < n-1 throw(ArgumentError("zminsv only if m >= n-1 but $m < $n-1")) end

    Q, R = qr(A)        # QR decomp. of A, maybe input
    if m == n - 1				# Corner case when s == 0 is guaranteed
        y = normalize([backsub(R[1:m,1:m], R[1:m,end]); -E])
        ss = zero(real(S))
	    # println("zminsv (m=$m, n=$n) returns ss= $ss y = $y")
        return ss, y
    end

    F = svd(A)
    return F.S[end], F.V[:,end]

    cr = U(0)
    scale = norm(A, Inf)     		    # get the magnitude of rows as scaler
	a = scale * normalize(randn(S,n))	# random initial vector (row)
    b = [scale; zeros(S, n)]

    T, trans = hessqr([a'; R])
    z = hqrt(trans, b)

    x = backsub(T[1:n,1:n], z[1:n])
    normalize!(x)

    r = [scale * x'; R] * x
    r = r - b

    for k = 1:50
		D = [ x' * scale * 2; R]		# Least squares problem min |D z - b|²
        T, trans = hessqr(D) 			# Hessenberg QR decomp. of stacked matrix
        z = hqrt(trans, r)        		# same QQ on b
        u = backsub(T[1:n,1:n], z[1:n])	# backward substitution to solve RR z = QQ' b
        y = normalize(x - u)			# new iteration
        r = [scale * y'; R] * y - b
        ss = norm(r[2:n])
        crk = norm(y-x)
		println("k = $k ss = $ss crk = $crk")

        if (k == 1 && crk < tol ) || (k > 1 &&  (crk >= cr || crk^2 / (cr - crk) < tol))
            break
        end
        x = y
        cr = crk
    end
	# println("zminsv (m=$m, n=$n) returns ss= $ss x = $x")
    ss, y
end

ext: jl


text:
# ------------------------------------------
# Julia caller
# ------------------------------------------

cd(@__DIR__)

w = [1.0, 0.0, 1.0]
γ = 1.4

###
# C
###

prim = zeros(3)

ccall(
    #(:hello, :libkit), # name of C function and library
    (:hello, "../c/libkit.so"), # name of C function and library
    Nothing,              # output type
    (),                        # tuple of input types
                                  # name of Julia variable to pass in
)
ccall(
    (:conserve_prim_1d, "../c/libkit.so"), # name of C function and library
    #(:conserve_prim_1d, "./c/kinetic.so"), # name of C function and library
    Nothing,              # output type
    (Ref{Cdouble}, Ref{Cdouble}, Cdouble,),                        # tuple of input types
    prim, w, γ                                  # name of Julia variable to pass in
)

println(prim)

###
# C++
###

prim = zeros(3)

ccall(
    #(:hello, :libkit), # name of C function and library
    (:hello, "../c++/libkit.so"), # name of C function and library
    Nothing,              # output type
    (),                        # tuple of input types
                                  # name of Julia variable to pass in
)
ccall(
    (:conserve_prim_1d, "../c++/libkit.so"), # name of C function and library
    #(:conserve_prim_1d, "./c/kinetic.so"), # name of C function and library
    Nothing,              # output type
    (Ref{Cdouble}, Ref{Cdouble}, Cdouble,),                        # tuple of input types
    prim, w, γ                                  # name of Julia variable to pass in
)

println(prim)

"""
C++ doesn't work yet
"""

max_stars_repo_name: vavrines/Kic

max_issues_repo_name: vavrines/Kic

max_forks_repo_name: vavrines/Kic


text:
using HDF5
using Healpix
using Random
using Random: MersenneTwister

"""
    read_halo_catalog_hdf5(filename)

Utility function to read an HDF5 table with x, y, z, M_h as the four rows.
The hdf5 record is "halos". This is a format we use to distribute Websky halos.

## Example
```julia-repl
julia> halo_pos, halo_mass = read_halo_catalog_hdf5(
    "/global/cfs/cdirs/sobs/www/users/Radio_WebSky/websky_halos-light.hdf5")
```
"""
function read_halo_catalog_hdf5(filename)
    hdata = h5open(filename, "r") do file
        read(file, "halos")
    end
    pos = hdata[1:3,:]
    halo_mass = hdata[4,:]
    return pos, halo_mass
end

"""
    chunk(arr_len, chunksize::Integer)

Generates a list of tuples which describe starting and ending chunk indices.
Useful for parallelizing an array operation.
"""
function chunk(arr_len, chunksize::Integer)
    return [(i, min(i + chunksize-1, arr_len))
        for i in range(1, arr_len, step=chunksize)]
end

function getrange(n)
    tid = Threads.threadid()
    nt = Threads.nthreads()
    d , r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid ≤ r ? 1 : 0)
    from:to
end

function threaded_rand!(arr::Array{T,1};
      chunksize=4096) where T

   num = size(arr,1)
   Threads.@threads for (i1, i2) in chunk(num, chunksize)
      @views rand!(arr[i1:i2])
   end
end

"""
Generate an array where the value at index i corresponds to the index of the
first source of halo i. Takes an array where the value at index i corresponds
to the number of subhalos that halo i has.
"""
function generate_subhalo_offsets(num_subhalos)
    result = cumsum(num_subhalos)
    prepend!(result, 0)
    return result
end

"""
    get_basic_halo_properties(halo_pos::Array{T,2}, model::AbstractForegroundModel,
                              cosmo::Cosmology.FlatLCDM{T}, res::Resolution) where T

Compute distance, redshift, and healpix indices from halo positions and a cosmology.
"""
function get_basic_halo_properties(halo_pos::Array{T,2}, model::AbstractForegroundModel,
                                   cosmo::Cosmology.FlatLCDM{T}, res::Resolution) where T
    N_halos = size(halo_pos, 2)
    hp_ind = Array{Int64}(undef, N_halos)  # healpix index of halo
    redshift = Array{T}(undef, N_halos)
    dist = Array{T}(undef, N_halos)

    r2z = build_r2z_interpolator(
        model.min_redshift, model.max_redshift, cosmo)
    Threads.@threads for i in 1:N_halos
        dist[i] = sqrt(halo_pos[1,i]^2 + halo_pos[2,i]^2 + halo_pos[3,i]^2)
        redshift[i] = r2z(dist[i])
        hp_ind[i] = Healpix.vec2pixRing(res, halo_pos[1,i], halo_pos[2,i], halo_pos[3,i])
    end

    return dist, redshift, hp_ind
end

"""
Compute angles of halos
"""
function get_angles(halo_pos::Array{T,2}) where T
    N_halos = size(halo_pos, 2)
    θ = Array{T}(undef, N_halos)
    ϕ = Array{T}(undef, N_halos)

    Threads.@threads for i in 1:N_halos
        θ[i], ϕ[i] = Healpix.vec2ang(halo_pos[1,i], halo_pos[2,i], halo_pos[3,i])
    end

    return θ, ϕ
end

"""
Utility function which prepends some zeros to an array. It makes a copy instead
of modifying the input.
"""
function ellpad(arr::Array{T,N}; nzeros=1) where {T,N}
    result = arr[:]
    pushfirst!(result, zeros(T, nzeros)...)
    return result
end

function catalog2map!(m::Map{T,RingOrder}, flux, theta, phi) where T
    res = m.resolution
    pixel_array = m.pixels
    N_halo = length(flux)

    # try to prevent thread issues by sorting by theta
    perm = sortperm(theta, rev=true, alg=ThreadsX.MergeSort)
    Threads.@threads for i_perm in 1:N_halo
        i_halo = perm[i_perm]
        hp_ind = Healpix.ang2pixRing(res, theta[i_halo], phi[i_halo])
        pixel_array[hp_ind] += flux[i_halo]
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(res.nside)
    pixel_array .*= per_pixel_steradian
end


text:
# CompactBitVector
# ----------------
#
# Classical indexable bit vector implementation with 1/4 bits per bit additional
# space.  Exactly speaking, this data structure is not succinct: n/4 (= Θ(n))
# additional space is required rather than o(n).  But in practice, this is
# almost the same as theoretically optimal additional space and efficient in
# actual computers.
#
# A compact bit vector has two types of blocks: large blocks and small blocks.
# Four small blocks are associated to a large block. Large blocks hold the rank
# value at the beginning of the block while small blocks hold the offset of the
# rank value from a large block the small block belongs to. If we naively store
# the offset values of small blocks, the first small blocks for each large block
# is always filled with zero. To expand the maximum length of a bit vector, we
# exploit these small blocks: the rank value at the beginning of a large block
# is smallblock[1] << 32 + largeblock. As a result, about 2^(8 + 32) = 2^40 bits
# = 128 GiB can be stored in a compact bit vector.
#
# Data Layout
# bits: ..|................|................|................|................|..
# sbs:    |  smallblock 1  |  smallblock 2  |  smallblock 3  |  smallblock 4  |
# lbs:    |                            large block                            |

type CompactBitVector <: AbstractIndexableBitVector
    # data
    bits::BitVector
    # large blocks
    lbs::Vector{Uint32}
    # small blocks
    sbs::Vector{Uint8}
end

maxlength(::Type{CompactBitVector}) = 2^40 - 1

function CompactBitVector()
    return CompactBitVector(convert(BitVector, Bool[]), Uint32[], Uint8[])
end

function CompactBitVector(v::Union(BitVector,Vector{Bool}))
    sv = CompactBitVector()
    for bit in v
        push!(sv, bit != 0)
    end
    return sv
end

function convert(::Type{CompactBitVector}, v::Union(BitVector,Vector{Bool}))
    return CompactBitVector(v)
end

function show(io::IO, v::CompactBitVector)
    write(io, "CompactBitVector(bits=[...],")
    write(io, "large_blocks=$(v.lbs),")
    write(io, "small_blocks=$(v.sbs))")
end

immutable LargeBlock end
immutable SmallBlock end

blocksizeof(::Type{SmallBlock}) =  64
blocksizeof(::Type{LargeBlock}) = 256
const n_smallblocks_per_largeblock = div(blocksizeof(LargeBlock), blocksizeof(SmallBlock))

function push!(v::CompactBitVector, bit::Bool)
    len = length(v) + 1
    if len > maxlength(CompactBitVector)
        error("overflow")
    end
    ensureroom!(v, len)
    if len % blocksizeof(LargeBlock) == 1 && length(v.lbs) > 1
        # the first bit of a large block
        rank  = convert(Int, v.lbs[end-1])
        rank += v.sbs[end-1]
        rank += count_ones(v.bits.chunks[end])
        v.lbs[end] = rank & typemax(Uint32)
        v.sbs[end] = rank >> bitsof(Uint32)
    elseif len % blocksizeof(SmallBlock) == 1 && length(v.sbs) > 1
        # the first bit of a small block
        v.sbs[end] = v.sbs[end-1] + count_ones(v.bits.chunks[end])
    end
    push!(v.bits, bit)
    return v
end
push!(v::CompactBitVector, b::Integer) = push!(v::CompactBitVector, b != 0)

length(v::CompactBitVector) = length(v.bits)
endof(v::CompactBitVector)  = length(v.bits)

getindex(v::CompactBitVector, i::Integer) = v.bits[i]

function rank1(v::CompactBitVector, i::Integer)
    if !(0 ≤ i ≤ endof(v))
        throw(BoundsError())
    end
    return unsafe_rank1(v, i)
end

function unsafe_rank1(v::CompactBitVector, i::Integer)
    if i == 0
        return 0
    end
    lbi = div(i - 1, blocksizeof(LargeBlock)) + 1
    sbi = div(i - 1, blocksizeof(SmallBlock)) + 1
    @inbounds byte = v.bits.chunks[sbi]
    r = rem(i, 64)
    if r != 0
        byte &= rmask(Uint64, r)
    end
    @inbounds ret = (
          convert(Int, v.sbs[sbi-rem(sbi, n_smallblocks_per_largeblock)+1]) << bitsof(Uint32)
        + v.lbs[lbi]
        + v.sbs[sbi]
        + count_ones(byte)
    )
    return ret
end

# ensure room to store `len` bits
function ensureroom!(v::CompactBitVector, len::Int)
    @assert len ≥ length(v)
    n_required_sbs = div(len - 1, blocksizeof(SmallBlock)) + 1
    if n_required_sbs ≤ length(v.sbs)
        # enough space to hold `len` bits
        return v
    end
    len_sbs = length(v.sbs)
    len_lbs = length(v.lbs)
    resize!(v.sbs, n_required_sbs)
    resize!(v.lbs, div(n_required_sbs - 1, 4) + 1)
    for i in len_sbs+1:endof(v.sbs) v.sbs[i] = 0 end
    for i in len_lbs+1:endof(v.lbs) v.lbs[i] = 0 end
    @assert length(v.sbs) ≤ length(v.lbs) * n_smallblocks_per_largeblock
    return v
end


text:
#EasyPlot: show functions
#-------------------------------------------------------------------------------

const SHOW_INDENTSTR = "   "

Base.show(io::IO, ::Axis{T}) where T = print(io, "Axis{$T}")

showcompact_lengthinfo(io::IO, d) = print(io, typeof(d))
showcompact_lengthinfo(io::IO, d::DataF1) = print(io, "DataF1: ", length(d), "pts")
showcompact_lengthinfo(io::IO, d::DataRS) = print(io, "DataRS:", paramlist(d))
showcompact_lengthinfo(io::IO, d::DataHR) = print(io, "DataHR:", paramlist(d))

function showshorthand(io::IO, wfrm::Waveform)
	print(io, "Waveform(")
	show(io, wfrm.label) #With "quotes"
	print(io, ", ")
	showcompact_lengthinfo(io, wfrm.data)
	print(io, ")")
end

function showindented(io::IO, p::Plot, indent::String)
	xaxis = p.xaxis
	print(io, indent, "Plot(")
	show(io, p.title) #With "quotes"
	print(io, ", $xaxis, [")
	rmg = length(p.ystriplist)
	for strip in p.ystriplist
		scale = strip.scale
		print(io, scale)
		rmg -= 1
		rmg > 0 ? print(io, ", ") : nothing
	end
	print(io, "])[\n")
	wfrmindent = indent * SHOW_INDENTSTR
	for wfrm in p.wfrmlist
		print(io, wfrmindent); showshorthand(io, wfrm); println(io)
	end
	print(io, indent, "]")
end

Base.show(io::IO, p::Plot) = showindented(io, p, "")

function Base.show(io::IO, pc::PlotCollection)
	print(io, "PlotCollection(")
	show(io, pc.title) #With "quotes"
	print(io, ")[\n")
	for p in pc.plotlist
		showindented(io, p, SHOW_INDENTSTR); println(io)
	end
	print(io, "]")
end

#Show dispatchable symbols more succinctly (Errors easier to debug)
Base.show(io::IO, t::Type{DS}) = print(io, "DS")
function Base.show(io::IO, t::Type{DS{DT}}) where DT
	sym = @isdefined(DT) ? DT : :T #Used in declarations like: cons(::DT{T}, ...) where T
	print(io, "DS{$sym}")
end

#Last line


text:
function test_diffusion_simple(fld)
    Nx, Ny, Nz = 1, 1, 16
    Lx, Ly, Lz = 1, 1, 1
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if fld ∈ (:u, :v, :w)
        field = getfield(model.velocities, fld)
    else
        field = getfield(model.tracers, fld)
    end

    value = π
    field.data .= value

    Δt = 0.01 # time-step much less than diffusion time-scale
    Nt = 10
    time_step!(model, Nt, Δt)

    !any(@. !isapprox(value, field.data))
end

function test_diffusion_budget(field_name)
    Nx, Ny, Nz = 1, 1, 16
    Lx, Ly, Lz = 1, 1, 1
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if field_name ∈ (:u, :v, :w)
        field = getfield(model.velocities, field_name)
    else
        field = getfield(model.tracers, field_name)
    end

    half_Nz = floor(Int, Nz/2)
    @views @. field.data[:, :, 1:half_Nz]  = -1
    @views @. field.data[:, :, half_Nz:end] =  1

    mean_init = mean(field.data)
    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 0.0001 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)
    isapprox(mean_init, mean(field.data))
end

function test_diffusion_cosine(fld)
    Nx, Ny, Nz = 1, 1, 128
    Lx, Ly, Lz = 1, 1, π/2
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if fld == :w
        throw("There are no boundary condition tests for w yet.")
    elseif fld ∈ (:u, :v)
        field = getfield(model.velocities, fld)
    else
        field = getfield(model.tracers, fld)
    end

    zC = model.grid.zC
    m = 2
    @views @. field.data[1, 1, :] = cos.(m*zC)

    diffusing_cosine(κ, m, z, t) = exp(-κ*m^2*t) * cos(m*z)

    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 1e-6 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)

    field_numerical = dropdims(field.data, dims=(1, 2))

    !any(@. !isapprox(field_numerical, diffusing_cosine(κ, m, zC, model.clock.time), atol=1e-6, rtol=1e-6))
end

function test_z_boundary_condition_simple(field_name, bctype, bc, Nx, Ny, Nz)
    Lx, Ly, Lz = 0.1, 0.2, 0.3
    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz))

    bc = BoundaryCondition(bctype, bc)
    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.right = bc # just set a boundary condition somewhere

    time_step!(model, 1, 1e-16)
    typeof(model) <: Model
end

function test_flux_budget(field_name)
    Nx, Ny, Nz = 1, 1, 16
    Lx, Ly, Lz = 1, 1, 0.7
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if field_name ∈ (:u, :v, :w)
        field = getfield(model.velocities, field_name)
    else
        field = getfield(model.tracers, field_name)
    end

    @. field.data = 0

    bottom_flux = 0.3
    flux_bc = BoundaryCondition(Flux, bottom_flux)
    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.right = flux_bc # "right" = "bottom" in the convention where k=Nz is the bottom.

    mean_init = mean(field.data)

    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 1e-6 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)

    # budget is Lz * ∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right);
    # therefore <ϕ> = bottom_flux * t / Lz
    isapprox(mean(field.data) - mean_init, bottom_flux*model.clock.time/Lz)
end


text:
const Returns = TimeArray

const Signal = TimeArray

const Weights = TimeArray

struct Brokerage
	transaction_costs::Float64
	# maxSlippage::Float64
	# borrowing_costs::TimeSeries
end

struct PortfolioReturns
	composition::Weights # the weights that were actually used (potentially after some normalization)
	returns::Returns # for the whole of the portfolio
	individual_returns::Returns # (for every asset, separately, this is useful for calculating risk contribution, etc)
end


text:
function draw_sphere(r, k, ambient, light)
    shades = ('.', ':', '!', '*', 'o', 'e', '&', '#', '%', '@')
    for i in floor(Int, -r):ceil(Int, r)
        x = i + 0.5
        line = IOBuffer()
        for j in floor(Int, -2r):ceil(2r)
            y = j / 2 + 0.5
            if x ^ 2 + y ^ 2 ≤ r ^ 2
                v = normalize([x, y, sqrt(r ^ 2 - x ^ 2 - y ^ 2)])
                b = dot(light, v) ^ k + ambient
                intensity = ceil(Int, (1 - b) * (length(shades) - 1))
                if intensity < 1
                    intensity = 1 end
                if intensity > length(shades)
                    intensity = length(shades) end
                print(shades[intensity])
            else
                print(' ')
            end
        end
        println()
    end
end

light = normalize([30, 30, -50])
draw_sphere(20, 4, 0.1, light)
draw_sphere(10, 2, 0.4, light)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export find_intersections

"""
    find_intersections(segments::Vector{Segment{T}}; tol=1e-9) where {T<:AbstractFloat}

Compute all the intersections between the segments using the Bentley-Ottmann algorithm.

Seting a tolerance value to work for all situations is tricky, so try to play with it if you
are finding that not all intersections are found.

# Examples
```jldoctest
julia> segments = [Segment(0,0,5,5), Segment(0,5,5,0)];
julia> intersections = find_intersections(segments)
julia> [Point(2.5, 2.5)]
```
"""
function find_intersections(segments::Vector{Segment{T}}; tol=1e-9) where {T<:AbstractFloat}
    Q = EventQueue(segments)
    y0 = Q.tree[1].y
    ys = [y0]
    status = Status(y0, tol)
    intersections = []
    while length(Q) != 0
        event = fetch_event!(Q)
        handle_event_point(Q, event, status, intersections)
    end
    return intersections
end

function handle_event_point(
    Q::EventQueue,
    event::Event,
    status::Status,
    intersections,
)
    p = event.point
    status.y_sweep = p.y
    U = event.segments
    L = Set{Segment}()
    C = Set{Segment}()
    for (key, segment) in status.dict
        if is_lower_end(segment, p)
            push!(L, segment)
        elseif contains(segment, p, status.tol)
            push!(C, segment)
        end
    end
    CL = union(C, L)
    total = union(CL, U)
    if length(total) > 1
        push!(intersections, p)
    end
    UC = union(U, C)
    update!(status, insert=UC, delete=CL)
    if length(UC) == 0
        sl = find_left(status, p.x)
        sr = find_right(status, p.x)
        find_new_event(Q, sl, sr, p)
    else
        sp, xp = find_leftmost(UC, p.y, status.tol)
        sl = find_left(status, xp)
        find_new_event(Q, sl, sp, p)
        spp, xpp = find_rightmost(UC, p.y, status.tol)
        sr = find_right(status, xpp)
        find_new_event(Q, spp, sr, p)
    end
    return
end

function find_new_event(Q::EventQueue, s1::Segment, s2::Segment, p::Point, tol=1e-9)
    if s1 == s2
        return
    end
    do_intersect, intersection = intersect!(s1, s2)
    if do_intersect
        if intersection.y < p.y
            insert!(Q, intersection)
        elseif (intersection.y ≈ p.y) && (intersection.x > p.x)
            insert!(Q, intersection)
        else
            return
        end
    end
end

# for when s1 or s2 are nothing
find_new_event(Q::EventQueue, s1, s2, p::Point) = nothing


text:
using DataStructures
using Distributed
using CuboidGrids; Lar = CuboidGrids

using Random
using BenchmarkTools
using Base.Threads

#=
using CuboidGrids; Lar = CuboidGrids
mod1D = Lar.grid(repeat([.1,-.1],outer=5)...)

using ViewerGL; GL = ViewerGL

GL.VIEW([ GL.GLFrame2, GL.GLGrid(mod1D..., GL.COLORS[1],1) ])

mod3D = Lar.INSR(Lar.larModelProduct)([mod1D,mod1D,mod1D])

GL.VIEW([ GL.GLFrame2, GL.GLPol(mod3D..., GL.COLORS[1],1) ])
=#

map(x->x*2, [1,2,3])
reduce(*, [2;3;4])

#=il problema svolto con mapreduce
    mapreduce(x-> rem(33333331, x), max, 1:16666665)    => 11111109
  diventa:
@parallel max for x in 1:16666665
    rem(33333331, x)                                    => 11111109
end

io devo fare:
@parallel <risultato ottenuto?> for x in <collezizone>
    <funzione da applicare>  (si fa su funzioni che hanno bisogno di un solo parametro?)
end
=#

function serial_add()
	s=0.0
	for i = 1:100000
		s=s+Random.randn()
	end
	return s
end

function parallel_add()
	return @distributed (+) for i=1:100000
	Random.randn()
	end
end

@btime serial_add()
@btime parallel_add()

@threads for i in 1:nthreads()
	a[threadid()] = threadid()
	end


text:
module D3Magic

# ported from https://github.com/ResidentMario/py_d3

export @d3_str, display_d3, d3update_display, cur_cell

const d3id = Ref(0)
d3version = "4.9.1"

setD3version(version::String) = global d3version = version

function d3code_with_id(userd3code, id)
    userd3code = replace(userd3code, r"d3.select\((?!this)", "d3.select$id(")
    return replace(userd3code, r"d3.selectAll\((?!this)", "d3.selectAll$id(")
end

function d3string(userd3code, id)
    userd3code = d3code_with_id(userd3code, id)
    """
    <g id="d3-cell-$id">
        <script>
            d3.select$id = function(selection) {
                return d3.select("#d3-cell-$id").select(selection);
            }
            d3.selectAll$id = function(selection) {
                return d3.select("#d3-cell-$id").selectAll(selection);
            }
            d3._select = d3.select
            d3._selectAll = d3.selectAll
        </script>
        $userd3code
    </g>
    """
end

d3update_string(userd3code, id) = d3code_with_id(userd3code, id)

macro d3_str(sraw)
    s = parse(string("\"\"\"", sraw, "\"\"\""))
    quote
        display_d3($(esc(s)))
    end
end

function display_d3(s, id = (d3id[] += 1))
    display("text/html", d3string(s, id))
end

function d3update_display(s, id = d3id[])
    display("text/html", d3update_string(s, id))
end

cur_cell() = d3id[]

function __init__()
    display("text/html","""<script>
        requirejs.config({
            paths: {
                d3: "//cdnjs.cloudflare.com/ajax/libs/d3/$d3version/d3"
            }
        });
        require(['d3'], function(d3) {
            window.d3 = d3;
        })
    </script>""")
end

end #module


text:
using Revise
using POMDPTesting
include("mdp_definition.jl")

function testindexing(mdp::MergingMDP)
    ss = states(mdp)
    @showprogress for (i, s) in enumerate(ss)
        si = stateindex(mdp, s)
        @assert si == i "si = $si , i = $i, s = $s"
        for a in actions(mdp)
            d = transition(mdp, s, a)
            for (sp, w) in weighted_iterator(d)
                try
                    spi = stateindex(mdp, sp)
                catch
                    throw("Error s = $s, a = $a, sp=$sp")
                end
            end
        end
    end
end

testindexing(mdp)

trans_prob_consistency_check(mdp)

function POMDPTesting.trans_prob_consistency_check(pomdp::Union{MDP, POMDP})
    # initalize spaces
    sspace = states(pomdp)
    aspace = actions(pomdp)
    # iterate through all s-a pairs
    @showprogress for s in sspace
        for a in aspace
            tran = transition(pomdp, s, a)
            p = 0.0
            for sp in sspace
                p += pdf(tran, sp)
            end
            @assert isapprox(p, 1.0) "Transition probability does not sum to unity for state: ", s, " action: ", a, " sums to", p
        end
    end
end


text:
import matplotlib as mplt
mplt.use('Agg')

import matplotlib.pyplot as plt
import sys, argparse
import numpy as np
from scipy import signal
from pprint import pprint

def mavg(y, span):

    N = len(y)
    yy = np.zeros((N,))
    if span == 1:
        yy[:] = y

    else:
        for i in range(N):
            if i < span:
                rng = slice(0, i+1)
                yy[i] = np.nan
            else:
                rng = slice(i-span,i)
                yy[i] = np.mean(y[rng])

    return yy

parser = argparse.ArgumentParser()
parser.add_argument('--input-dir')
parser.add_argument('--output-dir')
parser.add_argument('--res')
parser.add_argument('--casenames')
parser.add_argument('--data-file')
parser.add_argument('--varname')
parser.add_argument('--mavg', type=int, default=1)
parser.add_argument('--yscale', type=float, default=1.0)
parser.add_argument('--ylabel', default="")
parser.add_argument('--extra-title', default="")

args = parser.parse_args()

pprint(args)

casenames = args.casenames.split(",")

print("Going to compare these models:")
pprint(casenames)

tss = []

new_casenames = []

for i in range(len(casenames)):

    try:

        f = Nio.open_file("%s/%s_%s/%s" % (args.input_dir, args.res, casenames[i], args.data_file), "r")

    except Exception as e:

        print("Error happens when doing casename %s. Going to ignore this one..." % casenames[i])

        continue

    new_casenames.append(casenames[i])
    ts = mavg(f.variables[args.varname][:] / args.yscale, args.mavg)

    tss.append(ts)

    f.close()

casenames = new_casenames

N = len(tss[0])
time = np.arange(N) / 12
nyears = N / 12.0

fig, ax = plt.subplots(1, 1, figsize=(12, 8))

ax.set_title("%s (%d years) %s" % (args.varname, nyears, args.extra_title))
ax.set_xlabel("Time [years]")
ax.set_ylabel(args.ylabel)
ax.grid()
for i in range(len(casenames)):
    ax.plot(time, tss[i], linewidth=2, label=casenames[i])

ax.legend()
ax.grid()

fig.savefig("%s/%s_%s_mc_timeseries_%s.png" % (args.output_dir, args.label, args.res, args.varname), dpi=200)

#plt.show()

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
macro portfoliomodels()
  println("")
  reload(Pkg.dir("PortfolioModels/test/runtests.jl"))
end


text:
@reexport module LShaped

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using MathOptInterface
using StochasticPrograms
using StochasticPrograms: AcceptableTermination
using StochasticPrograms: UnspecifiedInstantiation, VerticalStructure, AbstractScenarioProblems, ScenarioProblems, DistributedScenarioProblems, DecisionChannel
using StochasticPrograms: AbstractExecution, Serial, Synchronous, Asynchronous
using StochasticPrograms: AbstractStructuredOptimizer, RelativeTolerance, MasterOptimizer, SubproblemOptimizer
using StochasticPrograms: SingleDecisionSet, update_decision_constraint!
using StochasticPrograms: set_known_decision!, update_known_decisions!, SingleKnownSet, KnownModification, KnownValuesChange
using StochasticPrograms: add_subscript
using StochasticPrograms: AbstractPenaltyterm, Quadratic, InfNorm, ManhattanNorm, initialize_penaltyterm!, update_penaltyterm!, remove_penalty!, remove_penalty_variables!, remove_penalty_constraints!
using ProgressMeter
using Clustering

import Base: show, put!, wait, isready, take!, fetch, zero, +, length, size
import StochasticPrograms: supports_structure, default_structure, check_loadable, load_structure!, restore_structure!, optimize!, optimizer_name, master_optimizer, subproblem_optimizer, num_subproblems, remove_penalty_variables!, remove_penalty_constraints!

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const CI = MOI.ConstraintIndex
const CutConstraint = CI{AffineDecisionFunction{Float64}, MOI.GreaterThan{Float64}}

export
    AbstractLShapedAttribute,
    LShapedAlgorithm,
    num_cuts,
    num_iterations,
    FeasibilityCuts,
    Regularizer,
    RegularizationParameter,
    RawRegularizationParameter,
    set_regularization_attribute,
    set_regularization_attributes,
    DontRegularize,
    NoRegularization,
    RegularizedDecomposition,
    WithRegularizedDecomposition,
    RD,
    WithRD,
    TrustRegion,
    WithTrustRegion,
    TR,
    WithTR,
    LevelSet,
    WithLevelSets,
    LV,
    WithLV,
    Consolidator,
    DontConsolidate,
    Consolidate,
    Consolidation,
    ConsolidationParameter,
    RawConsolidationParameter,
    set_consolidation_attribute,
    set_consolidation_attributes,
    at_tolerance,
    num_iters,
    tolerance_reached,
    Aggregator,
    AggregationParameter,
    RawAggregationParameter,
    set_aggregation_attribute,
    set_aggregation_attributes,
    DontAggregate,
    NoAggregation,
    PartialAggregate,
    PartialAggregation,
    Aggregate,
    FullAggregation,
    DynamicAggregate,
    DynamicAggregation,
    ClusterAggregate,
    ClusterAggregation,
    HybridAggregate,
    HybridAggregation,
    SelectUniform,
    SelectDecaying,
    SelectRandom,
    SelectClosest,
    SelectClosestToReference,
    StaticCluster,
    ClusterByReference,
    Kmedoids,
    Hierarchical,
    absolute_distance,
    angular_distance,
    spatioangular_distance

# Include files
include("types/types.jl")
include("consolidators/consolidation.jl")
include("aggregators/aggregation.jl")
include("regularizers/regularization.jl")
include("execution/execution.jl")
include("solver.jl")
include("MOI_wrapper.jl")

end # module


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

mutable struct Value_Size <: ProtoType
    width::Int64
    height::Int64
    Value_Size(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Value_Size

mutable struct ArrayValue <: ProtoType
    values::Base.Any
    ArrayValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ArrayValue (has cyclic type dependency)
const __ftype_ArrayValue = Dict(:values => "Base.Vector{Value}")
meta(t::Type{ArrayValue}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, __ftype_ArrayValue)

mutable struct Value <: ProtoType
    boolean_value::Bool
    integer_value::Int64
    double_value::Float64
    entity_value::Base.Any
    array_value::ArrayValue
    timestamp_value::ProtoBuf.google.protobuf.Timestamp
    string_value::AbstractString
    blob_value::Array{UInt8,1}
    size_value::Value_Size
    Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Value (has cyclic type dependency)
const __fnum_Value = Int[1,2,3,6,9,10,17,18,19]
const __ftype_Value = Dict(:entity_value => "Entity")
const __oneofs_Value = Int[1,1,1,1,1,1,1,1,1]
const __oneof_names_Value = [Symbol("value_type")]
meta(t::Type{Value}) = meta(t, ProtoBuf.DEF_REQ, __fnum_Value, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_Value, __oneof_names_Value, __ftype_Value)

mutable struct Entity_PropertiesEntry <: ProtoType
    key::AbstractString
    value::Value
    Entity_PropertiesEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Entity_PropertiesEntry (mapentry) (has cyclic type dependency)

mutable struct Entity <: ProtoType
    properties::Base.Dict # map entry
    Entity(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Entity (has cyclic type dependency)
const __ftype_Entity = Dict(:properties => "Base.Dict{AbstractString,Value}")
meta(t::Type{Entity}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, __ftype_Entity)

export ArrayValue, Value_Size, Value, Entity_PropertiesEntry, Entity, ArrayValue, Value, Entity_PropertiesEntry, Entity
# mapentries: "Entity_PropertiesEntry" => ("AbstractString", "Value")


text:
## Dyad-based network measures ##

"""
    density(g::Graph)

Graph density. It is a proportion of number of existing edges to the maximum
possible number of edges. Density of an empty graph is 0, and density of
a full graph is 1. Graphs having self-loops and/or multiple edges can have
density larger than 1.
"""
density(g::Graph) = edgecount(g) / maxedgecount(g)

"""
    mutuality(g::Graph)

Number of mutual dyads in the graph. A dyad is mutual if for two nodes
`A` and `B` there are both edge `A → B` and `B → A`.
Edges in undirected graphs always constitute mutual dyads.

# References

Moreno, J.L., and Jennings, H.H. (1938). "Statistics of Social Configurations." Sociometry, 1, 342-374.
"""
function mutuality(g::DirectedGraph)
    count = 0
    for n = g.nodes
        i, j = 1, 1
        while i < length(n.forward) && j < length(n.reverse)
            if n.forward[i] < n.reverse[j]
                i += 1
            elseif n.forward[i] > n.reverse[j]
                j += 1
            else
                count += 1
                i += 1
                j += 1
            end
        end
    end
    return count
end

mutuality(g::UndirectedGraph) = edgecount(g)

"""
    reciprocity(g::Graph)

Graph reciprocity. It is a proportion of mutual dyads to a total
number of non-null dyads.
[Wiki](https://en.wikipedia.org/wiki/Reciprocity_(network_science)#Traditional_definition)
"""
reciprocity(g::Graph) = mutuality(g) / edgecount(g)

"""
    dyadcensus(g::Graph)

Dyad census of the graph. It is a 3-tuple of form (M, A, N), standing for
the numbers of mutual, asymmetric, and null dyads.
"""
function dyadcensus(g::Graph)
    m = mutuality(g)
    return (m, edgecount(g) - m, maxedgecount(g) - edgecount(g))
end

"""
    selfloopcount(g::Graph)

Number of self-loops (edges that start and end at the same node) in the graph.
"""
function selfloopcount(g::Graph)
    cnt = 0
    @inbounds for i = nodes(g)
        @inbounds for j = neighbors(g, i)
            if i == j; cnt += 1; end
        end
    end
    return cnt
end

function selfloopnodes(g::Graph)
    ns = Int[]
    @inbounds for i = nodes(g)
        @inbounds for j = neighbors(g, i)
            if i == j
                push!(ns, i)
                break
            end
        end
    end
    return ns
end


text:
type WeaveDoc
    source::AbstractString
    basename::AbstractString
    path::AbstractString
    chunks::Array
    cwd::AbstractString
    format
    doctype::AbstractString
    function WeaveDoc(source, chunks)
        path, fname = splitdir(abspath(source))
        basename = splitext(fname)[1]
        new(source, basename, path, chunks, "", nothing, "")
    end
end

immutable ChunkOutput
    code::AbstractString
    stdout::AbstractString
    displayed::AbstractString
    rich_output::AbstractString
    figures::Array{AbstractString}
end

type CodeChunk
    content::AbstractString
    number::Int
    result_no::Int
    start_line::Int
    optionstring::AbstractString
    options::Dict{Symbol, Any}
    output::AbstractString
    rich_output::AbstractString
    figures::Array{AbstractString}
    result::Array{ChunkOutput}
    function CodeChunk(content, number, start_line, optionstring, options)
        new(content, number, 0, start_line, optionstring, options, "","", AbstractString[], ChunkOutput[])
    end
end

type DocChunk
    content::AbstractString
    number::Int
    start_line::Int
end

type TermResult
end

type ScriptResult
end

type CollectResult
end


text:
# # Serializing PowerSystem Data

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSystems.jl supports serializing/deserializing data with JSON. This notebook
# provides an example of how to write and read a `System` to/from disk.

# ### Dependencies
# Let's use a dataset from the [tabular data parsing example](../../notebook/2_PowerSystems_examples/parse_matpower.ipynb)
using SIIPExamples

pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "02_parse_matpower.jl"))

# ### Write data to a temporary directory

folder = mktempdir()
path = joinpath(folder, "system.json")
println("Serializing to $path")
to_json(sys, path)

filesize(path) / (1024 * 1024) #MiB

# ### Read the JSON file and create a new `System`
sys2 = System(path)


text:
using QuantumStatistics, AbstractTrees
# using NewickTree
include("parameter.jl")

chan = [Parquet.T, Parquet.U, Parquet.S]

para = Parquet.Para(chan, [1, 2])

ver4 = Parquet.Ver4{Weight}(2, 1, para)

########## use AbstractTrees interface to print/manipulate the tree
print_tree(ver4)

# [println(node) for node in Leaves(ver4)]  #print all loopNum=0 ver4

# [println(node) for node in PostOrderDFS(ver4)]  # iterator ver4 in depth-first search (children before parents)
# [println(node) for node in PreOrderDFS(ver4)]  # iterator ver4 (parents before children)

########## print tree to a newick format file  ##############
# io = open("./test.newick", "w")
# write(io, Parquet.newick(ver4))
# close(io)

########## use ete3 package to visualize tree
# Parquet.showTree(ver4, para, verbose=1, depth=3)  # visualize tree using python3 package ete3


text:
"""
# Paths in graphs

- Distances
- Breadth-first search
    - Correctness and efficiency
- Lengths on edges
- Dijkstra's algorithm
    - An adaptation of breadth-first search
    - A more convenient graph
    - Alarm clocks
    - Dijkstra's algorithm
    - An alternative derivation
    - Running time
    - Which heap is best?
- Priority queue implementations
    - Array
    - Binary heap
    - d-ary heap
- Shortest paths in the presence of negative edges
    - Negative edges
    - Negative cycles
- Shortest paths in dags
"""
module Chapter4

using DataStructures

export
    dijkstra

import
    ..Graph

"""
    dijkstra(G, s)

Run Dijkstra's shortest path algorithm from node `s` for a graph
using an adjacency list representation. Return the distances to
each node and their previous nodes.

# Examples

```jldoctest
julia> G = Dict(
           :A=>[(:B,6), (:D,1)],
           :B=>[(:A,6), (:C,5), (:D,2), (:E,2)],
           :C=>[(:B,5), (:E,5)],
           :D=>[(:A,1), (:B,2), (:E,1)],
           :E=>[(:B,2), (:C,5), (:D,1)]
           )
Dict{Symbol, Vector{Tuple{Symbol, Int64}}} with 5 entries:
  :A => [(:B, 6), (:D, 1)]
  :D => [(:A, 1), (:B, 2), (:E, 1)]
  :B => [(:A, 6), (:C, 5), (:D, 2), (:E, 2)]
  :E => [(:B, 2), (:C, 5), (:D, 1)]
  :C => [(:B, 5), (:E, 5)]

julia> dijkstra(G, :A)[:dist]
Dict{Symbol, Float64} with 5 entries:
  :A => 0.0
  :D => 1.0
  :B => 3.0
  :E => 2.0
  :C => 7.0

julia> dijkstra(G, :A)[:prev]
Dict{Symbol, Symbol} with 5 entries:
  :A => :nil
  :D => :A
  :B => :D
  :E => :D
  :C => :E
```
"""
dijkstra(G::Graph{T,U}, s::T) where {T<:Symbol, U<:Real} = dijkstra(G, s, :nil::Union{T,Nothing})
dijkstra(G::Graph{T,U}, s::T) where {T<:Real, U<:Real} = dijkstra(G, s, 0::Union{T,Nothing})
dijkstra(G::Graph{T,U}, s::T) where {T<:String, U<:Real} = dijkstra(G, s, ""::Union{T,Nothing})

function dijkstra(G::Graph{T,U}, s::T, nil::Union{T,Nothing}) where{T, U<:Real}
    dist = Dict(collect(keys(G)) .=> Inf)
    prev = Dict(collect(keys(G)) .=> nil)
    dist[s] = 0

    H = PriorityQueue([(s, dist[s])])
    while !isempty(H)
        (u, du) = dequeue_pair!(H)
        for (v, length) ∈ G[u]
            alt = dist[u] + length
            if dist[v] > alt
                dist[v] = alt
                prev[v] = u
                if !haskey(H, v)
                    enqueue!(H, v => alt)
                end
            end
        end
    end

    return (dist=dist, prev=prev)
end

end # module


text:
# Sum of the positive components

"""
  SumPositive()

Returns the function `g(x) = sum(max(0, x))`.
"""

immutable SumPositive <: ProximableFunction end

function (f::SumPositive){T <: Real}(x::AbstractArray{T})
  return sum(max.(0.0, x))
end

function prox!{T <: Real}(f::SumPositive, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  fsum = 0.0
  for i in eachindex(x)
    y[i] = x[i] < gamma ? (x[i] > 0.0 ? 0.0 : x[i]) : x[i]-gamma
    fsum += y[i] > 0.0 ? y[i] : 0.0
  end
  return fsum
end

fun_name(f::SumPositive) = "Sum of the positive coefficients"
fun_dom(f::SumPositive) = "AbstractArray{Real}"
fun_expr(f::SumPositive) = "x ↦ sum(max(0, x))"

function prox_naive{T <: Real}(f::SumPositive, x::AbstractArray{T}, gamma::Real=1.0)
  y = copy(x)
  indpos = x .> 0.0
  y[indpos] = max.(0.0, x[indpos]-gamma)
  return y, sum(max.(0.0, y))
end

# ######################### #
# Prox with multiple gammas #
# ######################### #

function prox!{T <: Real}(f::SumPositive, x::AbstractArray{T}, y::AbstractArray{T}, gamma::AbstractArray{T})
  fsum = 0.0
  for i in eachindex(x)
    y[i] = x[i] < gamma[i] ? (x[i] > 0.0 ? 0.0 : x[i]) : x[i]-gamma[i]
    fsum += y[i] > 0.0 ? y[i] : 0.0
  end
  return fsum
end

function prox_naive{T <: Real}(f::SumPositive, x::AbstractArray{T}, gamma::AbstractArray{T})
  y = copy(x)
  indpos = x .> 0.0
  y[indpos] = max.(0.0, x[indpos]-gamma[indpos])
  return y, sum(max.(0.0, y))
end


text:
using Preprocessing, Potential

function linearsolve(Ae, Af, As, be, bf, bs)

    cefs = zeros(0)
    if (length(Ae)>0) & (length(Af)>0) & (length(As)>0) & (length(be)>0) & (length(bf)>0) & (length(bs)>0)
        cefs = cat(cat(Ae, Af, dims=1), As, dims=1)\[be[:]; bf[:]; bs[:]]
    end

    cef = zeros(0)
    if (length(Ae)>0) & (length(Af)>0) & (length(be)>0) & (length(bf)>0)
        cef = cat(Ae, Af, dims=1)\[be[:]; bf[:]]
    end

    ce = zeros(0)
    if (length(Ae)>0) & (length(be)>0)
        ce = Ae\be
    end

    cf = zeros(0)
    if (length(Af)>0) & (length(bf)>0)
        cf = Af\bf
    end

    return ce, cf, cef, cefs
end

function lsqsolve(Ae, Af, As, be, bf, bs)

    cefs = zeros(0)
    cef = zeros(0)
    ce = zeros(0)
    cf = zeros(0)

    sumAe = sum(abs.(Ae[:])) > 1e-15
    sumAf = sum(abs.(Af[:])) > 1e-15
    sumAs = sum(abs.(As[:])) > 1e-15
    sumbe = sum(abs.(be[:])) > 1e-15
    sumbf = sum(abs.(bf[:])) > 1e-15
    sumbs = sum(abs.(bs[:])) > 1e-15

    if (sumAe) & (sumAf) & (sumAs) & (sumbe) & (sumbf) & (sumbs)
        A = Ae+Af+As
        for i = 1:size(A,1)
            A[i,i] = A[i,i]*(1 + 1e-15)
        end
        cefs = (A)\(be+bf+bs)
    end
    if (sumAe) & (sumAf) & (sumbe) & (sumbf)
        A = Ae+Af
        for i = 1:size(A,1)
            A[i,i] = A[i,i]*(1 + 1e-15)
        end
        cef = A\(be+bf)
    end
    if (sumAe) & (sumbe)
        A = 1.0*Ae
        for i = 1:size(A,1)
            A[i,i] = A[i,i]*(1 + 1e-15)
        end
        ce = A\be
    end
    if (sumAf) & (sumbf)
        A = 1.0*Af
        for i = 1:size(A,1)
            A[i,i] = A[i,i]*(1 + 1e-15)
        end
        if abs(A[1,1])<=1e-10
            A[1,1] = 1e-10
        end
        cf = A\bf
    end

    return ce, cf, cef, cefs
end

function insertarray(ce, cei, i)
    if (length(cei)>0)
        ce[:,i] = cei
    end
    return ce
end

function catarray(Ae, be, Aei, bei)
    if (length(Aei) > 0) & (length(bei) > 0)
        Ae = cat(Ae, Aei, dims=1)
        be = [be; bei]
    end
    return Ae, be
end

function calmae(Aei, bei, cei)
    mae = -1.0
    if (length(bei) > 0) & (length(cei)>0) & (length(Aei)>0)
        mae = sum(abs.(bei - Aei*cei))/length(bei)
    end
    return mae
end

function calerrors(Aei, bei, cei)
    mae = -1.0
    rmse = - 1.0
    rsq = -1.0
    ssr = -1.0
    avg = -1.0
    if (length(bei) > 0) & (length(cei)>0) & (length(Aei)>0)
        res = bei - Aei*cei
        mae = sum(abs.(res))/length(bei)
        ssr = sum(res.^2)
        mse = ssr /length(bei);
        rmse = sqrt(mse)
        avg = sum(bei)/length(bei);
        rsq = 1.0 - ssr/sum((bei .- avg).^2);
    end
    return mae, rmse, rsq, ssr, avg
end

function applyweight(Aei, bei, we)
    if (length(bei) > 0) & (length(Aei) > 0) & (length(we) > 0)
        Aei = we[1]*Aei
        bei = we[1]*bei
    end
    return Aei, bei
end

function applytranspose(Aei, bei, m, n)
    Cei = zeros(m,n)
    if (length(bei) > 0) & (length(Aei) > 0)
        Cei = Aei'*bei
    end
    return Cei
end

function linearfit(data, emdescriptors, mldescriptors, eta, kappa, pbc=nothing, a=nothing, b=nothing, c=nothing, method=nothing, normalizeenergy=nothing)

    if normalizeenergy === nothing
        normalizeenergy = 0 # divide energy by number of atoms
    end
    if method === nothing
        method = "lsq"
    end

    n = length(data)
    emae = -ones(Float64,n,4)
    fmae = -ones(Float64,n,4)
    smae = -ones(Float64,n,4)
    ce = 1.0
    cf = 1.0
    cef = 1.0
    cefs = 1.0
    Ae = 1.0
    Af = 1.0
    As = 1.0
    be = 1.0
    bf = 1.0
    bs = 1.0
    De = 1.0
    Df = 1.0
    Ds = 1.0
    de = 1.0
    df = 1.0
    ds = 1.0
    Le = 1.0
    Lf = 1.0
    Ls = 1.0
    he = 1.0
    hf = 1.0
    hs = 1.0
    for i = 1:n
        # training data set
        if typeof(data[i]) == Preprocessing.DataStruct
            config, indices = Preprocessing.readconfigdata(data[i], pbc, a, b, c)
        else
            config = data[i]
        end
        training = Array(1:config.nconfigs)

        Aei, Afi, Asi, bei, bfi, bsi, nfi = Potential.emldescriptors(config, training, eta, kappa, normalizeenergy, emdescriptors, mldescriptors)

        m = size(Aei,2) # total number of descriptors
        if i == 1
            ce = ce*zeros(Float64,m,n)
            cf = cf*zeros(Float64,m,n)
            cef = cef*zeros(Float64,m,n)
            cefs = cefs*zeros(Float64,m,n)
            De = De*zeros(m,m,n)
            Df = Df*zeros(m,m,n)
            Ds = Ds*zeros(m,m,n)
            de = de*zeros(m,n)
            df = df*zeros(m,n)
            ds = ds*zeros(m,n)
            Le = Le*zeros(m,m)
            Lf = Lf*zeros(m,m)
            Ls = Ls*zeros(m,m)
            he = he*zeros(m)
            hf = hf*zeros(m)
            hs = hs*zeros(m)
        end

        De[:,:,i] = applytranspose(Aei, Aei, m, m)
        Df[:,:,i] = applytranspose(Afi, Afi, m, m)
        Ds[:,:,i] = applytranspose(Asi, Asi, m, m)
        de[:,i] = applytranspose(Aei, bei, m, 1)
        df[:,i] = applytranspose(Afi, bfi, m, 1)
        ds[:,i] = applytranspose(Asi, bsi, m, 1)

        Aei, bei = applyweight(Aei, bei, config.we)
        Afi, bfi = applyweight(Afi, bfi, config.wf)
        Asi, bsi = applyweight(Asi, bsi, config.ws)

        Le += applytranspose(Aei, Aei, m, m)
        Lf += applytranspose(Afi, Afi, m, m)
        Ls += applytranspose(Asi, Asi, m, m)
        he += applytranspose(Aei, bei, m, 1)
        hf += applytranspose(Afi, bfi, m, 1)
        hs += applytranspose(Asi, bsi, m, 1)

        if length(data)==1
            if method == "lsq"
                cei, cfi, cefi, cefsi = lsqsolve(Le, Lf, Ls, he, hf, hs)
            else
                cei, cfi, cefi, cefsi = linearsolve(Aei, Afi, Asi, bei, bfi, bsi)
            end
            ce = insertarray(ce, cei, i)
            cf = insertarray(cf, cfi, i)
            cef = insertarray(cef, cefi, i)
            cefs = insertarray(cefs, cefsi, i)

            emae[i,1] = calmae(Aei, bei, cei)
            emae[i,2] = calmae(Aei, bei, cfi)
            emae[i,3] = calmae(Aei, bei, cefi)
            emae[i,4] = calmae(Aei, bei, cefsi)
            fmae[i,1] = calmae(Afi, bfi, cei)
            fmae[i,2] = calmae(Afi, bfi, cfi)
            fmae[i,3] = calmae(Afi, bfi, cefi)
            fmae[i,4] = calmae(Afi, bfi, cefsi)
            smae[i,1] = calmae(Asi, bsi, cei)
            smae[i,2] = calmae(Asi, bsi, cfi)
            smae[i,3] = calmae(Asi, bsi, cefi)
            smae[i,4] = calmae(Asi, bsi, cefsi)

            return ce, cf, cef, cefs, emae, fmae, smae, De, Df, Ds, de, df, ds
        end

        if (i == 1) | (method == "lsq")
            Ae = Aei
            Af = Afi
            As = Asi
            be = bei
            bf = bfi
            bs = bsi
        else
            Ae, be = catarray(Ae, be, Aei, bei)
            Af, bf = catarray(Af, bf, Afi, bfi)
            As, bs = catarray(As, bs, Asi, bsi)
        end

        if (method == "lsq")
            cei, cfi, cefi, cefsi = lsqsolve(Le, Lf, Ls, he, hf, hs)
        else
            cei, cfi, cefi, cefsi = linearsolve(Ae, Af, As, be, bf, bs)
        end

        ce = insertarray(ce, cei, i)
        cf = insertarray(cf, cfi, i)
        cef = insertarray(cef, cefi, i)
        cefs = insertarray(cefs, cefsi, i)

        emae[i,1] = calmae(Ae, be, cei)
        emae[i,2] = calmae(Ae, be, cfi)
        emae[i,3] = calmae(Ae, be, cefi)
        emae[i,4] = calmae(Ae, be, cefsi)
        fmae[i,1] = calmae(Af, bf, cei)
        fmae[i,2] = calmae(Af, bf, cfi)
        fmae[i,3] = calmae(Af, bf, cefi)
        fmae[i,4] = calmae(Af, bf, cefsi)
        smae[i,1] = calmae(As, bs, cei)
        smae[i,2] = calmae(As, bs, cfi)
        smae[i,3] = calmae(As, bs, cefi)
        smae[i,4] = calmae(As, bs, cefsi)
    end

    return ce, cf, cef, cefs, emae, fmae, smae, De, Df, Ds, de, df, ds
end

function calmean(emae, szbe)
    n,p,m = size(emae)
    em = -ones(p,m)
    if (sum(szbe) > 0) & (sum(emae[:]) > 0)
        for k = 1:p
            for j = 1:m
                em[k,j] = 0.0
                for i = 1:n
                    em[k,j] = em[k,j] + emae[i,k,j]*szbe[i]
                end
                em[k,j] =  em[k,j]/sum(szbe)
            end
        end
    end
    return em
end

function validation(data, emdescriptors, mldescriptors, ce, cf, cef, cefs, eta, kappa, pbc=nothing, a=nothing, b=nothing, c=nothing)

    normalizeenergy = 0
    m = size(ce,2)
    n = length(data)
    emae = -ones(Float64,n,4,m)
    fmae = -ones(Float64,n,4,m)
    smae = -ones(Float64,n,4,m)
    szbe = zeros(Int64, n)
    szbf = zeros(Int64, n)
    szbs = zeros(Int64, n)
    for i = 1:n
        if typeof(data[i]) == Preprocessing.DataStruct
            config, indices = Preprocessing.readconfigdata(data[i], pbc, a, b, c)
        else
            config = data[i]
        end
        valid = Array(1:config.nconfigs)

        Aei, Afi, Asi, bei, bfi, bsi, nfi = Potential.emldescriptors(config, valid, eta, kappa, normalizeenergy, emdescriptors, mldescriptors)
        szbe[i] = length(bei)
        szbf[i] = length(bfi)
        szbs[i] = length(bsi)

        for j = 1:m
            cei = ce[:,j]
            cfi = cf[:,j]
            cefi = cef[:,j]
            cefsi = cefs[:,j]
            emae[i,1,j] = calmae(Aei, bei, cei)
            emae[i,2,j] = calmae(Aei, bei, cfi)
            emae[i,3,j] = calmae(Aei, bei, cefi)
            emae[i,4,j] = calmae(Aei, bei, cefsi)
            fmae[i,1,j] = calmae(Afi, bfi, cei)
            fmae[i,2,j] = calmae(Afi, bfi, cfi)
            fmae[i,3,j] = calmae(Afi, bfi, cefi)
            fmae[i,4,j] = calmae(Afi, bfi, cefsi)
            smae[i,1,j] = calmae(Asi, bsi, cei)
            smae[i,2,j] = calmae(Asi, bsi, cfi)
            smae[i,3,j] = calmae(Asi, bsi, cefi)
            smae[i,4,j] = calmae(Asi, bsi, cefsi)
        end
    end

    eme = calmean(emae, szbe)
    fme = calmean(fmae, szbf)
    sme = calmean(smae, szbs)

    return eme, fme, sme, emae, fmae, smae, szbe, szbf, szbs
end

function linearfit2(data, descriptors, potential, optim, pbc=nothing, a=nothing, b=nothing, c=nothing)

    data = deletenothing(data)
    descriptors = deletenothing(descriptors)
    potential = deletenothing(potential)

    if pbc === nothing
        pbc = [1; 1; 1] # periodic boundary conditions
    end

    lossfunc, method, normalizeenergy, normalizestress, eta, kappa,
        weightouter, etaspace, kappalist = getoptim(optim)

    n = length(data)
    emae = -ones(Float64,n,4)
    fmae = -ones(Float64,n,4)
    smae = -ones(Float64,n,4)
    ce = 1.0
    cf = 1.0
    cef = 1.0
    cefs = 1.0
    Ae = 1.0
    Af = 1.0
    As = 1.0
    be = 1.0
    bf = 1.0
    bs = 1.0
    De = 1.0
    Df = 1.0
    Ds = 1.0
    de = 1.0
    df = 1.0
    ds = 1.0
    Le = 1.0
    Lf = 1.0
    Ls = 1.0
    he = 1.0
    hf = 1.0
    hs = 1.0
    for i = 1:n
        # training data set
        if typeof(data[i]) == Preprocessing.DataStruct
            config, indices = Preprocessing.readconfigdata(data[i], pbc, a, b, c)
        else
            config = data[i]
        end
        training = Array(1:config.nconfigs)

        Aei, Afi, Asi, bei, bfi, bsi, nfi = Potential.emldescriptors2(config, training, eta, kappa,
                normalizeenergy, normalizestress, descriptors, potential)

        m = size(Aei,2) # total number of descriptors
        if i == 1
            ce = ce*zeros(Float64,m,n)
            cf = cf*zeros(Float64,m,n)
            cef = cef*zeros(Float64,m,n)
            cefs = cefs*zeros(Float64,m,n)
            De = De*zeros(m,m,n)
            Df = Df*zeros(m,m,n)
            Ds = Ds*zeros(m,m,n)
            de = de*zeros(m,n)
            df = df*zeros(m,n)
            ds = ds*zeros(m,n)
            Le = Le*zeros(m,m)
            Lf = Lf*zeros(m,m)
            Ls = Ls*zeros(m,m)
            he = he*zeros(m)
            hf = hf*zeros(m)
            hs = hs*zeros(m)
        end

        De[:,:,i] = applytranspose(Aei, Aei, m, m)
        Df[:,:,i] = applytranspose(Afi, Afi, m, m)
        Ds[:,:,i] = applytranspose(Asi, Asi, m, m)
        de[:,i] = applytranspose(Aei, bei, m, 1)
        df[:,i] = applytranspose(Afi, bfi, m, 1)
        ds[:,i] = applytranspose(Asi, bsi, m, 1)

        Aei, bei = applyweight(Aei, bei, config.we)
        Afi, bfi = applyweight(Afi, bfi, config.wf)
        Asi, bsi = applyweight(Asi, bsi, config.ws)

        Le += applytranspose(Aei, Aei, m, m)
        Lf += applytranspose(Afi, Afi, m, m)
        Ls += applytranspose(Asi, Asi, m, m)
        he += applytranspose(Aei, bei, m, 1)
        hf += applytranspose(Afi, bfi, m, 1)
        hs += applytranspose(Asi, bsi, m, 1)

        if length(data)==1
            if optim.method == "lsq"
                cei, cfi, cefi, cefsi = lsqsolve(Le, Lf, Ls, he, hf, hs)
            else
                cei, cfi, cefi, cefsi = linearsolve(Aei, Afi, Asi, bei, bfi, bsi)
            end
            ce = insertarray(ce, cei, i)
            cf = insertarray(cf, cfi, i)
            cef = insertarray(cef, cefi, i)
            cefs = insertarray(cefs, cefsi, i)

            emae[i,1] = calmae(Aei, bei, cei)
            emae[i,2] = calmae(Aei, bei, cfi)
            emae[i,3] = calmae(Aei, bei, cefi)
            emae[i,4] = calmae(Aei, bei, cefsi)
            fmae[i,1] = calmae(Afi, bfi, cei)
            fmae[i,2] = calmae(Afi, bfi, cfi)
            fmae[i,3] = calmae(Afi, bfi, cefi)
            fmae[i,4] = calmae(Afi, bfi, cefsi)
            smae[i,1] = calmae(Asi, bsi, cei)
            smae[i,2] = calmae(Asi, bsi, cfi)
            smae[i,3] = calmae(Asi, bsi, cefi)
            smae[i,4] = calmae(Asi, bsi, cefsi)

            return ce, cf, cef, cefs, emae, fmae, smae, De, Df, Ds, de, df, ds
        end

        if (i == 1) | (optim.method == "lsq")
            Ae = 1.0*Aei
            Af = 1.0*Afi
            As = 1.0*Asi
            be = 1.0*bei
            bf = 1.0*bfi
            bs = 1.0*bsi
        else
            Ae, be = catarray(Ae, be, Aei, bei)
            Af, bf = catarray(Af, bf, Afi, bfi)
            As, bs = catarray(As, bs, Asi, bsi)
        end

        if (optim.method == "lsq")
            #display(Le)
            cei, cfi, cefi, cefsi = lsqsolve(Le, Lf, Ls, he, hf, hs)
        else
            cei, cfi, cefi, cefsi = linearsolve(Ae, Af, As, be, bf, bs)
        end

        # if i==n
        #     display(i)
        #     A = Le+Lf+Ls
        #     b = he+hf+hs
        #     fileID = open("Ajl","w"); write(fileID,A[:]); close(fileID);
        #     fileID = open("bjl","w"); write(fileID,b[:]); close(fileID);
        #     display(A)
        #     display(b)
        #     coeff = (0.5*(A+A'))\b
        #     display(coeff)
        #     # display(size(Ae))
        #     # display(size(Af))
        #     # display(size(As))
        #     coeff = cat(cat(Ae, Af, dims=1), As, dims=1)\[be[:]; bf[:]; bs[:]]
        #     # A = Ae+Af+As
        #     # b = be+bf+bs
        #     # coeff = A\b
        #     display(coeff)
        #     error("here")
        # end

        ce = insertarray(ce, cei, i)
        cf = insertarray(cf, cfi, i)
        cef = insertarray(cef, cefi, i)
        cefs = insertarray(cefs, cefsi, i)

        emae[i,1] = calmae(Ae, be, cei)
        emae[i,2] = calmae(Ae, be, cfi)
        emae[i,3] = calmae(Ae, be, cefi)
        emae[i,4] = calmae(Ae, be, cefsi)
        fmae[i,1] = calmae(Af, bf, cei)
        fmae[i,2] = calmae(Af, bf, cfi)
        fmae[i,3] = calmae(Af, bf, cefi)
        fmae[i,4] = calmae(Af, bf, cefsi)
        smae[i,1] = calmae(As, bs, cei)
        smae[i,2] = calmae(As, bs, cfi)
        smae[i,3] = calmae(As, bs, cefi)
        smae[i,4] = calmae(As, bs, cefsi)
    end

    if lossfunc == "energy"
        coeff = ce[:,end]
    elseif lossfunc == "force"
        coeff = cf[:,end]
    elseif lossfunc == "energyforce"
        coeff = cef[:,end]
    elseif lossfunc == "energyforcestress"
        coeff = cefs[:,end]
    end

    return coeff, ce, cf, cef, cefs, emae, fmae, smae, De, Df, Ds, de, df, ds
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# Quadratic decision function #
# ========================== #
mutable struct QuadraticDecisionConstraintBridge{T, S, F, LQ} <: MOIB.Constraint.AbstractBridge
    constraint::CI{F, S}
    decision_function::QuadraticDecisionFunction{T, LQ}
    set::S
end

function MOIB.Constraint.bridge_constraint(::Type{QuadraticDecisionConstraintBridge{T,S,F,LQ}},
                                           model,
                                           f::QuadraticDecisionFunction{T, LQ},
                                           set::S) where {T, S, F <: MOI.ScalarAffineFunction{T}, LQ <: LinearPart{T}}
    # All decisions have been mapped to either the decision part constant
    # or the variable part terms at this point
    lq = f.linear_quadratic_terms
    # Calculate total constant
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.ScalarAffineFunction(
                                        lq.variable_part.terms,
                                        zero(T)),
                                    MOIU.shift_constant(set, -constant))
    # Save the constraint index, the decision function, and the set, to allow modifications
    return QuadraticDecisionConstraintBridge{T, S, F, LQ}(constraint, f, set)
end

function MOIB.Constraint.bridge_constraint(::Type{QuadraticDecisionConstraintBridge{T,S,F,LQ}},
                                           model,
                                           f::QuadraticDecisionFunction{T, LQ},
                                           set::S) where {T, S, F <: MOI.ScalarQuadraticFunction{T}, LQ <: QuadraticPart{T}}
    # All decisions have been mapped to either the decision part constant
    # or the variable part terms at this point
    lq = f.linear_quadratic_terms
    # Calculate total constant
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.ScalarQuadraticFunction(
                                        lq.variable_part.affine_terms,
                                        lq.variable_part.quadratic_terms,
                                        zero(T)),
                                    MOIU.shift_constant(set, -constant))
    # Save the constraint index, the decision function, and the set, to allow modifications
    return QuadraticDecisionConstraintBridge{T, S, F, LQ}(constraint, f, set)
end

function MOI.supports_constraint(::Type{<:QuadraticDecisionConstraintBridge{T}},
                                 ::Type{<:QuadraticDecisionFunction{T}},
                                 ::Type{<:MOI.AbstractScalarSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:QuadraticDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{QuadraticDecisionConstraintBridge{T, S, F, LQ}}) where {T, S, F, LQ}
    return [(F, S)]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:QuadraticDecisionConstraintBridge{T}},
                                              ::Type{QuadraticDecisionFunction{T, LinearPart{T}}},
                                              S::Type{<:MOI.AbstractScalarSet}) where T
    return QuadraticDecisionConstraintBridge{T, S, MOI.ScalarAffineFunction{T}, LinearPart{T}}
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:QuadraticDecisionConstraintBridge{T}},
                                              ::Type{QuadraticDecisionFunction{T, QuadraticPart{T}}},
                                              S::Type{<:MOI.AbstractScalarSet}) where T
    return QuadraticDecisionConstraintBridge{T, S, MOI.ScalarQuadraticFunction{T}, QuadraticPart{T}}
end

MOI.get(b::QuadraticDecisionConstraintBridge{T, S, F}, ::MOI.NumberOfConstraints{F, S}) where {T, S, F} = 1
MOI.get(b::QuadraticDecisionConstraintBridge{T, S, F}, ::MOI.ListOfConstraintIndices{F, S}) where {T, S, F} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarAffineFunction{T}}) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     lq.decision_part.terms)
        result |= any(t -> t.variable_index == v,
                      f.known_decision_terms.affine_terms)
    end
    g = QuadraticDecisionFunction(
        LinearPart(
            MOIU.filter_variables(v -> !from_decision(v), lq.variable_part),
            copy(lq.decision_part)),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
    return g
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarQuadraticFunction{T}}) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     lq.decision_part.affine_terms)
        result |= any(t -> t.variable_index_1 == v || t.variable_index_2 == v,
                      lq.decision_part.quadratic_terms)
        result |= any(t -> t.variable_index_1 == v,
                      lq.cross_terms.affine_terms)
        result |= any(t -> t.variable_index == v,
                      f.known_decision_terms.affine_terms)
    end
    g = QuadraticDecisionFunction(
        QuadraticPart(
            MOIU.filter_variables(v -> !from_decision(v), lq.variable_part),
            copy(lq.decision_part),
            copy(lq.cross_terms)),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
    return g
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::QuadraticDecisionConstraintBridge{T,S}) where {T,S}
    return bridge.set
end

function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual},
                 bridge::QuadraticDecisionConstraintBridge{T,S}) where {T,S}
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.delete(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::QuadraticDecisionConstraintBridge{T,S,F,LinearPart{T}},
                 f::QuadraticDecisionFunction{T,LinearPart{T}}) where {T,S,F}
    # Update bridge functions and function constant
    bridge.decision_function = f
    lq = f.linear_quadratic_terms
    # Change the function of the bridged constraints
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.ScalarAffineFunction(
                lq.variable_part.terms,
                zero(T)))
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::QuadraticDecisionConstraintBridge{T,S,F,QuadraticPart{T}},
                 f::QuadraticDecisionFunction{T,QuadraticPart{T}}) where {T,S,F}
    # Update bridge functions and function constant
    bridge.decision_function = f
    lq = f.linear_quadratic_terms
    # Change the function of the bridged constraints
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.ScalarQuadraticFunction(
                lq.variable_part.affine_terms,
                lq.variable_part.quadratic_terms,
                zero(T)))
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::QuadraticDecisionConstraintBridge{T,S}, change::S) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    bridge.set = change
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S}, change::MOI.ScalarConstantChange) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Modify constant of variable part
    lq.variable_part.constant = change.new_constant
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarAffineFunction{T}}, change::MOI.ScalarCoefficientChange) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Update coefficient in variable part
    modify_coefficient!(lq.variable_part.terms, change.variable, change.new_coefficient)
    # Modify decision part if variable originated from a decision
    i = something(findfirst(t -> t.variable_index == change.variable,
                            lq.decision_part.terms), 0)
    if !iszero(i)
        modify_coefficient!(lq.decision_part.terms, change.variable, change.new_coefficient)
    end
    # Modify the variable part of the constraint as usual
    MOI.modify(model, bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarQuadraticFunction{T}}, change::MOI.ScalarCoefficientChange) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Update variable part
    modify_coefficient!(lq.variable_part.affine_terms, change.variable, change.new_coefficient)
    # Modify decision part if variable originated from a decision
    i = something(findfirst(t -> t.variable_index == change.variable,
                            lq.decision_part.affine_terms), 0)
    if !iszero(i)
        modify_coefficient!(lq.decision_part.affine_terms, change.variable, change.new_coefficient)
    end
    # Complete rebridge required
    MOI.modify(model, bridge.constraint, DecisionsStateChange())
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S}, change::DecisionCoefficientChange) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    variable_terms = lq isa LinearPart ? lq.variable_part.terms : lq.variable_part.affine_terms
    decision_terms = lq isa LinearPart ? lq.decision_part.terms : lq.decision_part.affine_terms
    # Check if decision has been mapped
    i = something(findfirst(t -> t.variable_index == change.decision,
                            variable_terms), 0)
    if !iszero(i)
        # Update mapped variable through ScalarCoefficientChange
        MOI.modify(model, bridge, MOI.ScalarCoefficientChange(change.decision, change.new_coefficient))
    else
        i = something(findfirst(t -> t.variable_index == change.decision,
                                decision_terms), 0)
        # Query the fixed decision value
        unbridged = MOIB.unbridged_variable_function(model, change.decision)::MOI.SingleVariable
        decision_value = MOI.get(model, MOI.VariablePrimal(), unbridged.variable)
        # Update decision part of constraint constant
        coefficient = iszero(i) ? zero(T) : decision_terms[i].coefficient
        lq.decision_part.constant +=
            (change.new_coefficient - coefficient) * decision_value
        # Recalculate total constant and shift constraint set
        constant = lq.variable_part.constant +
            lq.decision_part.constant +
            f.known_part.constant +
            f.known_decision_terms.constant
        MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
                MOIU.shift_constant(bridge.set, -constant))
    end
    # Update coefficient in decision part
    modify_coefficient!(decision_terms, change.decision, change.new_coefficient)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S}, change::KnownCoefficientChange) where {T,S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    i = something(findfirst(t -> t.variable_index == change.known,
                            f.known_part.affine_terms), 0)
    # Update known part of objective constant
    coefficient = iszero(i) ? zero(T) : f.known_part.affine_terms[i].coefficient
    known_value = MOI.get(model, MOI.VariablePrimal(), change.known)
    f.known_part.constant +=
        (change.new_coefficient - coefficient) * known_value
    # Update coefficient in known part
    modify_coefficient!(f.known_part.affine_terms, change.known, change.new_coefficient)
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarAffineFunction{T}}, ::DecisionStateChange) where {T, S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Switch on state transition
    if change.new_state == NotTaken
        # Check if decision is included in the decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.decision_part.terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = lq.decision_part.terms[i].coefficient
            # Update fixed decision value
            change.value_difference < 0 ||
                error("Decision value update should be negative when transitioning to NotTaken state.")
            lq.decision_part.constant +=
                coefficient * change.value_difference
            # Modify variable part
            push!(lq.variable_part.terms,
                  MOI.ScalarAffineTerm(coefficient, change.decision))
            # Modify the coefficient of the mapped variable
            MOI.modify(model, bridge.constraint,
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
        # Check if decision is included in the known decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                f.known_decision_terms.affine_terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = f.known_decision_terms.affine_terms[i].coefficient
            # Update fixed decision value
            change.value_difference < 0 ||
                error("Decision value update should be negative when transitioning to NotTaken state.")
            f.known_decision_terms.constant +=
                coefficient * change.value_difference
            # Modify variable part
            add_term!(lq.variable_part.terms,
                      MOI.ScalarAffineTerm(coefficient, change.decision))
            # Get current coefficient in variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                    lq.variable_part.terms), 0)
            coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
            # Modify the coefficient of the mapped variable
            MOI.modify(model, bridge.constraint,
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
    end
    if change.new_state == Taken
        # Check if decision is included in the decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.decision_part.terms), 0)
        if !iszero(i)
            # Update fixed decision value
            coefficient = lq.variable_part.terms[i].coefficient
            lq.decision_part.constant +=
                coefficient * change.value_difference
            # Modify variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.variable_part.terms), 0)
            !iszero(i) && deleteat!(lq.variable_part.terms, i)
            # Modify the coefficient of the mapped variable
            MOI.modify(model, bridge.constraint,
                       MOI.ScalarCoefficientChange(change.decision, zero(T)))
        end
        # Check if decision is included in the known decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                f.known_decision_terms.affine_terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = f.known_decision_terms.affine_terms[i].coefficient
            # Update fixed decision value
            f.known_decision_terms.constant +=
                coefficient * change.value_difference
            # Modify variable part
            remove_term!(lq.variable_part.terms,
                         MOI.ScalarAffineTerm(coefficient, change.decision))
            # Get current coefficient in variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                    lq.variable_part.terms), 0)
            coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
            # Modify the coefficient of the mapped variable
            MOI.modify(model, bridge.constraint,
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
    end
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarQuadraticFunction{T}}, ::DecisionStateChange) where {T, S}
    # Complete rebridge only option because quadratic coefficients cannot be modified
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarAffineFunction{T}}, ::DecisionsStateChange) where {T, S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Unbridge the constraint function and cache
    # all occuring decisions.
    unbridged_func = zero(QuadraticDecisionFunction{T,LinearPart{T}})
    variables_to_remove = Vector{Int}()
    # Known/variable cross terms
    for term in f.known_variable_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(lq.variable_part.terms, term)
    end
    empty!(f.known_variable_terms.affine_terms)
    for (i, term) in enumerate(f.known_variable_terms.quadratic_terms)
        push!(unbridged_func.known_variable_terms.quadratic_terms, copy(term))
    end
    empty!(f.known_variable_terms.quadratic_terms)
    # Known/decision cross terms
    for term in f.known_decision_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(lq.variable_part.terms, term)
        remove_term!(lq.decision_part.terms, term)
    end
    empty!(f.known_decision_terms.affine_terms)
    for (i, term) in enumerate(f.known_decision_terms.quadratic_terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_2)
        push!(unbridged_func.known_decision_terms.quadratic_terms,
              MOI.ScalarQuadraticTerm(term.coefficient,
                                      term.variable_index_1,
                                      unbridged.variable))
    end
    empty!(f.known_decision_terms.quadratic_terms)
    # Affine terms
    for (i,term) in enumerate(lq.decision_part.terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index)
        push!(unbridged_func.linear_quadratic_terms.decision_part.terms,
              MOI.ScalarAffineTerm(term.coefficient, unbridged.variable))
        # Check if a mapped variable exists and remove it as well if so
        j = something(findfirst(t -> t.variable_index == term.variable_index,
                                lq.variable_part.terms), 0)
        j != 0 && push!(variables_to_remove, j)
    end
    empty!(lq.decision_part.affine_terms)
    # Remove terms that come from bridged decisions
    deleteat!(f.variable_part.affine_terms, variables_to_remove)
    # Reset constants
    f.decision_part.constant = zero(T)
    f.known_decision_terms.constant = zero(T)
    # Rebridge
    MOIU.operate!(+, T, f, MOIB.bridged_function(model, unbridged_func))
    # All decisions have been mapped to either the decision part constant
    # or the variable part terms at this point.
    # Delete previous constraint
    MOI.delete(model, bridge.constraint)
    # Recalculate total constant and readd constraint
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    bridge.constraint = MOI.add_constraint(model,
                                           MOI.ScalarQuadraticFunction(
                                               lq.variable_part.affine_terms,
                                               lq.variable_part.quadratic_terms,
                                               zero(T)),
                                           MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarQuadraticFunction{T}}, ::DecisionsStateChange) where {T, S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Unbridge the constraint function and cache
    # all occuring decisions.
    unbridged_func = zero(QuadraticDecisionFunction{T,QuadraticPart{T}})
    aff_variables_to_remove = Vector{Int}()
    quad_variables_to_remove = Vector{Int}()
    # Cross terms
    for (i,term) in enumerate(lq.cross_terms.affine_terms)
        # Subtract any added variable/decision cross terms
        remove_term!(lq.variable_part.affine_terms, term)
    end
    empty!(lq.cross_terms.affine_terms)
    for (i, term) in enumerate(lq.cross_terms.quadratic_terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_1)
        push!(unbridged_func.linear_quadratic_terms.cross_terms.quadratic_terms,
              MOI.ScalarQuadraticTerm(term.coefficient,
                                      unbridged.variable,
                                      term.variable_index_2))
        # Check if a mapped variable exists and remove it as well if so
        j = something(findfirst(t -> t.variable_index_2 == term.variable_index_2,
                                lq.variable_part.quadratic_terms), 0)
        j != 0 && push!(quad_variables_to_remove, j)
    end
    empty!(lq.cross_terms.quadratic_terms)
    # Known/variable cross terms
    for term in f.known_variable_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(lq.variable_part.affine_terms, term)
    end
    empty!(f.known_variable_terms.affine_terms)
    for (i, term) in enumerate(f.known_variable_terms.quadratic_terms)
        push!(unbridged_func.known_variable_terms.quadratic_terms, copy(term))
    end
    empty!(f.known_variable_terms.quadratic_terms)
    # Known/decision cross terms
    for term in f.known_decision_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(lq.variable_part.affine_terms, term)
        remove_term!(lq.decision_part.affine_terms, term)
    end
    empty!(f.known_decision_terms.affine_terms)
    for (i, term) in enumerate(f.known_decision_terms.quadratic_terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_2)
        push!(unbridged_func.known_decision_terms.quadratic_terms,
              MOI.ScalarQuadraticTerm(term.coefficient,
                                      term.variable_index_1,
                                      unbridged.variable))
    end
    empty!(f.known_decision_terms.quadratic_terms)
    # Affine decision terms
    for (i,term) in enumerate(lq.decision_part.affine_terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index)
        push!(unbridged_func.linear_quadratic_terms.decision_part.affine_terms,
              MOI.ScalarAffineTerm(term.coefficient, unbridged.variable))
        # Check if a mapped variable exists and remove it as well if so
        j = something(findfirst(t -> t.variable_index == term.variable_index,
                                lq.variable_part.affine_terms), 0)
        j != 0 && push!(aff_variables_to_remove, j)
    end
    empty!(lq.decision_part.affine_terms)
    # Quadratic decision terms
    for (i,term) in enumerate(lq.decision_part.quadratic_terms)
        unbridged_1 = MOIB.unbridged_variable_function(model, term.variable_index_1)
        unbridged_2 = MOIB.unbridged_variable_function(model, term.variable_index_2)
        push!(unbridged_func.linear_quadratic_terms.decision_part.quadratic_terms,
              MOI.ScalarQuadraticTerm(term.coefficient,
                                      unbridged_1.variable,
                                      unbridged_2.variable))
        # Check if a mapped variable exists and remove it as well if so
        j = something(findfirst(t -> t.variable_index_1 == term.variable_index_1 ||
                                t.variable_index_2 == term.variable_index_2,
                                lq.variable_part.quadratic_terms), 0)
        j != 0 && push!(quad_variables_to_remove, j)
    end
    empty!(lq.decision_part.quadratic_terms)
    # Remove terms that come from bridged decisions
    deleteat!(lq.variable_part.affine_terms, aff_variables_to_remove)
    deleteat!(lq.variable_part.quadratic_terms, quad_variables_to_remove)
    # Reset constants
    lq.decision_part.constant = zero(T)
    f.known_decision_terms.constant = zero(T)
    # Rebridge
    MOIU.operate!(+, T, f, MOIB.bridged_function(model, unbridged_func))
    # All decisions have been mapped to either the decision part constant
    # or the variable part terms at this point.
    # Delete previous constraint
    MOI.delete(model, bridge.constraint)
    # Recalculate total constant and readd constraint
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    bridge.constraint = MOI.add_constraint(model,
                                           MOI.ScalarQuadraticFunction(
                                               lq.variable_part.affine_terms,
                                               lq.variable_part.quadratic_terms,
                                               zero(T)),
                                           MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T}, ::KnownValueChange) where T
    # Complete rebridge only option because quadratic coefficients cannot be modified
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarAffineFunction{T}}, change::KnownValuesChange) where {T, S}
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    F = QuadraticDecisionFunction{T}
    # Recalculate known constant
    known_val = zero(T)
    for term in f.known_part.affine_terms
        known_val += term.coefficient * MOI.get(model, MOI.VariablePrimal(), term.variable_index)
    end
    for term in f.known_part.quadratic_terms
        val_1 = MOI.get(model, MOI.VariablePrimal(), term.variable_index_1)
        val_2 = MOI.get(model, MOI.VariablePrimal(), term.variable_index_2)
        coeff = val_1 == val_2 ? term.coefficient / 2 : term.coefficient
        known_val += coeff * val_1 * val_2
    end
    # Update known part of objective
    f.known_part.constant = known_val
    # Handle cross terms
    for term in f.known_variable_terms.quadratic_terms
        coefficient = term.coefficient
        known_index = term.variable_index_1
        var_index = term.variable_index_2
        i = something(findfirst(t -> t.variable_index == var_index,
                                f.known_variable_terms.affine_terms), 0)
        old_coefficient = iszero(i) ? zero(T) : f.known_variable_terms.affine_terms[i].coefficient
        new_coefficient = coefficient * MOI.get(model, MOI.VariablePrimal(), known_index)
        # Update cross term
        add_term!(f.known_variable_terms.affine_terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, var_index))
        # Update variable part
        add_term!(lq.variable_part.terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, var_index))
        # Get current coefficient in variable part
        i = something(findfirst(t -> t.variable_index == var_index,
                                lq.variable_part.terms), 0)
        coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
        # Modify the coefficient of the mapped variable
        MOI.modify(model, bridge.constraint,
                   MOI.ScalarCoefficientChange(var_index, coefficient))
    end
    for term in f.known_decision_terms.quadratic_terms
        coefficient = term.coefficient
        known_index = term.variable_index_1
        dvar_index = term.variable_index_2
        i = something(findfirst(t -> t.variable_index == dvar_index,
                                f.known_decision_terms.affine_terms), 0)
        old_coefficient = iszero(i) ? zero(T) : f.known_decision_terms.affine_terms[i].coefficient
        new_coefficient = coefficient * MOI.get(model, MOI.VariablePrimal(), known_index)
        # Update cross term
        add_term!(f.known_decision_terms.affine_terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, dvar_index))
        # Update decision part
        add_term!(lq.decision_part.terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, dvar_index))
        # Get current coefficient in decision part
        i = something(findfirst(t -> t.variable_index == dvar_index,
                                lq.decision_part.terms), 0)
        coefficient = iszero(i) ? zero(T) : lq.decision_part.terms[i].coefficient
        # Proper decision modification
        unbridged = MOIB.unbridged_variable_function(model, dvar_index)::MOI.SingleVariable
        MOIB.modify_bridged_change(model, bridge,
                                   DecisionCoefficientChange(unbridged.variable, coefficient))
    end
    # Recalculate total constant and shift constraint set
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionConstraintBridge{T,S,MOI.ScalarQuadraticFunction{T}}, change::KnownValuesChange) where {T, S}
    # Complete rebridge only option because quadratic coefficients cannot be modified
    MOI.modify(model, bridge, DecisionsStateChange())
    return nothing
end