text:
export VolumeMesh
export points, name, groups, materials, material_names, ishomogeneous, points_in_vtk_type

struct VolumeMesh{Dim,T,U} <: AbstractMesh
    points::Vector{Point{Dim,T}}
    offsets::Vector{U}              # First index in connectivity
    connectivity::Vector{U}         # Point IDs that compose each element
    types::Vector{UInt8}            # VTK integer type of the element
    materials::Vector{UInt8}        # ID of the element's material
    material_names::Vector{String}
    name::String
    groups::Dict{String,BitSet}     # "Label"=>{IDs of elements with this label}
end

points(mesh::VolumeMesh) = mesh.points
name(mesh::VolumeMesh) = mesh.name
groups(mesh::VolumeMesh) = mesh.groups
materials(mesh::VolumeMesh) = mesh.materials
material_names(mesh::VolumeMesh) = mesh.material_names

const VTK_TRIANGLE::UInt8 = 5
const VTK_QUAD::UInt8 = 9
const VTK_QUADRATIC_TRIANGLE::UInt8 = 22
const VTK_QUADRATIC_QUAD::UInt8 = 23
const VTK_TETRA::UInt8 = 10
const VTK_HEXAHEDRON::UInt8 = 12
const VTK_QUADRATIC_TETRA::UInt8 = 24
const VTK_QUADRATIC_HEXAHEDRON::UInt8 = 25

ishomogeneous(mesh::VolumeMesh) = all(i->mesh.types[1] === mesh.types[i], 2:length(mesh.types))

function points_in_vtk_type(vtk_type::I) where {I<:Integer}
    if vtk_type == VTK_TRIANGLE
        return 3
    elseif vtk_type == VTK_QUAD
        return 4
    elseif vtk_type == VTK_QUADRATIC_TRIANGLE
        return 6
    elseif vtk_type == VTK_QUADRATIC_QUAD
        return 8
    elseif vtk_type == VTK_TETRA
        return 4
    elseif vtk_type == VTK_HEXAHEDRON
        return 8
    elseif vtk_type == VTK_QUADRATIC_TETRA
        return 10
    elseif vtk_type == VTK_QUADRATIC_HEXAHEDRON
        return 12
    else
        error("Unsupported type.")
        return nothing
    end
end

function vtk_alias_string(vtk_type::I) where {I<:Integer}
    if vtk_type == VTK_TRIANGLE
        return "Triangle"
    elseif vtk_type == VTK_QUAD
        return "Quadrilateral"
    elseif vtk_type == VTK_QUADRATIC_TRIANGLE
        return "QuadraticTriangle"
    elseif vtk_type == VTK_QUADRATIC_QUAD
        return "QuadraticQuadrilateral"
    elseif vtk_type == VTK_TETRA
        return "Tetrahedron"
    elseif vtk_type == VTK_HEXAHEDRON
        return "Hexahedron"
    elseif vtk_type == VTK_QUADRATIC_TETRA
        return "QuadraticTetrahedron"
    elseif vtk_type == VTK_QUADRATIC_HEXAHEDRON
        return "QuadraticHexahedron"
    else
        error("Unsupported type.")
        return nothing
    end
end

function Base.show(io::IO, mesh::VolumeMesh{Dim,T,U}) where {Dim,T,U}
    println(io, "VolumeMesh{",Dim, ", ",T,", ",U,"}")
    println(io, "  ├─ Name      : ", mesh.name)
    size_B = Base.summarysize(mesh)
    if size_B < 1e6
        println(io, "  ├─ Size (KB) : ", string(@sprintf("%.3f", size_B/1000)))
    else
        println(io, "  ├─ Size (MB) : ", string(@sprintf("%.3f",size_B/1e6)))
    end
    println(io, "  ├─ Points    : ", length(mesh.points))
    if Dim === 3
        println(io, "  ├─ Faces     : ", length(mesh.offsets) - 1)
    else
        println(io, "  ├─ Cells     : ", length(mesh.offsets) - 1)
    end
    unique_types = unique(mesh.types)
    nunique_types = length(unique_types)
    for i = 1:nunique_types
        type = unique_types[i]
        nelements = count(x->x === type,  mesh.types)
        if i === nunique_types
            println(io, "  │  └─ ", rpad(vtk_alias_string(type), 22), ": ", nelements)
        else
            println(io, "  │  ├─ ", rpad(vtk_alias_string(type), 22), ": ", nelements)
        end
    end
    println(io, "  ├─ Materials : ", length(mesh.material_names))
    ngroups = length(mesh.groups)
    println(io, "  └─ Groups    : ", ngroups)
    if 0 < ngroups ≤ 5
        group_keys = sort!(collect(keys(mesh.groups)))
        for i = 1:ngroups
            if i === ngroups
                println(io, "     └─ ", group_keys[i])
            else
                println(io, "     ├─ ", group_keys[i])
            end
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Combinatorics

function masked(mask, lis)
    combos = []
    idx = 1
    for step in mask
        if(step < 1)
            push!(combos, Array{Int,1}[])
        else
            push!(combos, sort(lis[idx:idx+step-1]))
            idx += step
        end
    end
    Array{Array{Int, 1}, 1}(combos)
end

function orderedpartitions(mask)
    tostring(masklis) = replace("$masklis", r"Array{Int\d?\d?,1}|Int\d?\d?", "")
    join([tostring(lis) for lis in unique([masked(mask, p)
                        for p in permutations(1:sum(mask))])], "\n")
end

println(orderedpartitions([2, 0, 2]))
println(orderedpartitions([1, 1, 1]))

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Julia wrapper for header: cheader/tox_old.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function tox_callback_group_invite(tox::Ptr{Cint},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_invite,libtoxcore),Void,(Ptr{Cint},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_message(tox::Ptr{Cint},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_message,libtoxcore),Void,(Ptr{Cint},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_action(tox::Ptr{Cint},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_action,libtoxcore),Void,(Ptr{Cint},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_title(tox::Ptr{Cint},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_title,libtoxcore),Void,(Ptr{Cint},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_namelist_change(tox::Ptr{Cint},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_namelist_change,libtoxcore),Void,(Ptr{Cint},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_add_groupchat(tox::Ptr{Cint})
    ccall((:tox_add_groupchat,libtoxcore),Cint,(Ptr{Cint},),tox)
end

function tox_del_groupchat(tox::Ptr{Cint},groupnumber::Cint)
    ccall((:tox_del_groupchat,libtoxcore),Cint,(Ptr{Cint},Cint),tox,groupnumber)
end

function tox_group_peername(tox::Ptr{Cint},groupnumber::Cint,peernumber::Cint,name::Ptr{Cint})
    ccall((:tox_group_peername,libtoxcore),Cint,(Ptr{Cint},Cint,Cint,Ptr{Cint}),tox,groupnumber,peernumber,name)
end

function tox_group_peer_pubkey(tox::Ptr{Cint},groupnumber::Cint,peernumber::Cint,public_key::Ptr{Cint})
    ccall((:tox_group_peer_pubkey,libtoxcore),Cint,(Ptr{Cint},Cint,Cint,Ptr{Cint}),tox,groupnumber,peernumber,public_key)
end

function tox_invite_friend(tox::Ptr{Cint},friendnumber::Cint,groupnumber::Cint)
    ccall((:tox_invite_friend,libtoxcore),Cint,(Ptr{Cint},Cint,Cint),tox,friendnumber,groupnumber)
end

function tox_join_groupchat(tox::Ptr{Cint},friendnumber::Cint,data::Ptr{Cint},length::Cint)
    ccall((:tox_join_groupchat,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Cint),tox,friendnumber,data,length)
end

function tox_group_message_send(tox::Ptr{Cint},groupnumber::Cint,message::Ptr{Cint},length::Cint)
    ccall((:tox_group_message_send,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Cint),tox,groupnumber,message,length)
end

function tox_group_action_send(tox::Ptr{Cint},groupnumber::Cint,action::Ptr{Cint},length::Cint)
    ccall((:tox_group_action_send,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Cint),tox,groupnumber,action,length)
end

function tox_group_set_title(tox::Ptr{Cint},groupnumber::Cint,title::Ptr{Cint},length::Cint)
    ccall((:tox_group_set_title,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Cint),tox,groupnumber,title,length)
end

function tox_group_get_title(tox::Ptr{Cint},groupnumber::Cint,title::Ptr{Cint},max_length::Cint)
    ccall((:tox_group_get_title,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Cint),tox,groupnumber,title,max_length)
end

function tox_group_peernumber_is_ours(tox::Ptr{Cint},groupnumber::Cint,peernumber::Cint)
    ccall((:tox_group_peernumber_is_ours,libtoxcore),Uint32,(Ptr{Cint},Cint,Cint),tox,groupnumber,peernumber)
end

function tox_group_number_peers(tox::Ptr{Cint},groupnumber::Cint)
    ccall((:tox_group_number_peers,libtoxcore),Cint,(Ptr{Cint},Cint),tox,groupnumber)
end

function tox_group_get_names(tox::Ptr{Cint},groupnumber::Cint,names::Ptr{Cint},lengths::Ptr{Cint},length::Cint)
    ccall((:tox_group_get_names,libtoxcore),Cint,(Ptr{Cint},Cint,Ptr{Cint},Ptr{Cint},Cint),tox,groupnumber,names,lengths,length)
end

function tox_count_chatlist()
    ccall((:tox_count_chatlist,libtoxcore),Cint,())
end

function tox_get_chatlist()
    ccall((:tox_get_chatlist,libtoxcore),Cint,())
end

function tox_group_get_type(tox::Ptr{Cint},groupnumber::Cint)
    ccall((:tox_group_get_type,libtoxcore),Cint,(Ptr{Cint},Cint),tox,groupnumber)
end


text:
module LinearOps
import FFTW
import Hankel
import Luna: Modes, Grid, PhysData, Maths
import Luna.PhysData: wlfreq

#=================================================#
#===============    FREE SPACE     ===============#
#=================================================#
"""
    make_const_linop(grid, xygrid, n, frame_vel)

Make constant linear operator for full 3D propagation. `n` is the refractive index (array)
and β1 is 1/velocity of the reference frame.
"""
function make_const_linop(grid::Grid.RealGrid, xygrid::Grid.FreeGrid,
                          n::AbstractArray, β1::Number)
    kperp2 = @. (xygrid.kx^2)' + xygrid.ky^2
    idcs = CartesianIndices((length(xygrid.ky), length(xygrid.kx)))
    k2 = zero(grid.ω)
    k2[grid.sidx] .= (n[grid.sidx] .* grid.ω[grid.sidx] ./ PhysData.c).^2
    out = zeros(ComplexF64, (length(grid.ω), length(xygrid.ky), length(xygrid.kx)))
    _fill_linop_xy!(out, grid, β1, k2, kperp2, idcs)
    return out
end

function make_const_linop(grid::Grid.RealGrid, xygrid::Grid.FreeGrid, nfun)
    n = zero(grid.ω)
    n[grid.sidx] = nfun.(2π*PhysData.c./grid.ω[grid.sidx])
    β1 = PhysData.dispersion_func(1, nfun)(grid.referenceλ)
    make_const_linop(grid, xygrid, n, β1)
end

function make_const_linop(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid,
                          n::AbstractArray, β1::Number, β0ref::Number; thg=false)
    kperp2 = @. (xygrid.kx^2)' + xygrid.ky^2
    idcs = CartesianIndices((length(xygrid.ky), length(xygrid.kx)))
    k2 = zero(grid.ω)
    k2[grid.sidx] .= (n[grid.sidx].*grid.ω[grid.sidx]./PhysData.c).^2
    out = zeros(ComplexF64, (length(grid.ω), length(xygrid.ky), length(xygrid.kx)))
    _fill_linop_xy!(out, grid, β1, k2, kperp2, idcs, β0ref; thg=thg)
    return out
end

function make_const_linop(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid, nfun,
                          thg=false)
    n = zero(grid.ω)
    n[grid.sidx] = nfun.(wlfreq.(grid.ω[grid.sidx]))
    β1 = PhysData.dispersion_func(1, nfun)(grid.referenceλ)
    if thg
        β0const = 0.0
    else
        β0const = grid.ω0/PhysData.c * nfun(wlfreq(grid.ω0))
    end
    make_const_linop(grid, xygrid, n, β1, β0const; thg=thg)
end

"""
    make_linop(grid, xygrid, nfun)

Make z-dependent linear operator for free-space propagation. `nfun(ω; z)` should return the
refractive index as a function of frequency `ω` and (kwarg) propagation distance `z`.
"""
function make_linop(grid::Grid.RealGrid, xygrid::Grid.FreeGrid, nfun)
    kperp2 = @. (xygrid.kx^2)' + xygrid.ky^2
    idcs = CartesianIndices((length(xygrid.ky), length(xygrid.kx)))
    k2 = zero(grid.ω)
    nfunλ(z) = λ -> nfun(wlfreq(λ), z=z)
    function linop!(out, z)
        β1 = PhysData.dispersion_func(1, nfunλ(z))(grid.referenceλ)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z) .* grid.ω[grid.sidx] ./ PhysData.c).^2
        _fill_linop_xy!(out, grid, β1, k2, kperp2, idcs)
    end
end

# Internal routine -- function barrier aids with JIT compilation
function _fill_linop_xy!(out, grid::Grid.RealGrid, β1::Float64, k2, kperp2, idcs)
    for ii in idcs
        for iω in eachindex(grid.ω)
            βsq = k2[iω] - kperp2[ii]
            if βsq < 0
                # negative βsq -> evanescent fields -> attenuation
                out[iω, ii] = -im*(-β1*grid.ω[iω]) - min(sqrt(abs(βsq)), 200)
            else
                out[iω, ii] = -im*(sqrt(βsq) - β1*grid.ω[iω])
            end
        end
    end
end

function make_linop(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid, nfun; thg=false)
    kperp2 = @. (xygrid.kx^2)' + xygrid.ky^2
    idcs = CartesianIndices((length(xygrid.ky), length(xygrid.kx)))
    k2 = zero(grid.ω)
    nfunλ(z) = λ -> nfun(wlfreq(λ), z=z)
    function linop!(out, z)
        β1 = PhysData.dispersion_func(1, nfunλ(z))(grid.referenceλ)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z).*grid.ω[grid.sidx]./PhysData.c).^2
        βref = thg ? 0.0 : grid.ω0/PhysData.c * nfun(grid.ω0; z=z)
        _fill_linop_xy!(out, grid, β1, k2, kperp2, idcs, βref; thg=thg)
    end
end

function _fill_linop_xy!(out, grid::Grid.EnvGrid, β1::Float64, k2, kperp2, idcs, βref; thg)
    for ii in idcs
        for iω in eachindex(grid.ω)
            βsq = k2[iω] - kperp2[ii]
            if βsq < 0
                # negative βsq -> evanescent fields -> attenuation
                out[iω, ii] = -im*(-β1*grid.ω[iω]) - min(sqrt(abs(βsq)), 200)
            else
                out[iω, ii] = -im*(sqrt(βsq) - β1*grid.ω[iω])
            end
            if !thg
                out[iω, ii] -= -im*βref
            end
        end
    end
end

#=================================================#
#==============   RADIAL SYMMETRY   ==============#
#=================================================#
"""
    make_const_linop(grid, q::QDHT, n, frame_vel)

Make constant linear operator for radial free-space. `n` is the refractive index (array)
and β1 is 1/velocity of the reference frame.
"""
function make_const_linop(grid::Grid.RealGrid, q::Hankel.QDHT,
                          n::AbstractArray, β1::Number)
    out = Array{ComplexF64}(undef, (length(grid.ω), q.N))
    k2 = @. (n*grid.ω/PhysData.c)^2
    kr2 = q.k.^2
    _fill_linop_r!(out, grid, β1, k2, kr2, q.N)
    return out
end

function make_const_linop(grid::Grid.RealGrid, q::Hankel.QDHT, nfun)
    n = zero(grid.ω)
    n[grid.sidx] = nfun.(2π*PhysData.c./grid.ω[grid.sidx])
    β1 = PhysData.dispersion_func(1, nfun)(grid.referenceλ)
    make_const_linop(grid, q, n, β1)
end

function make_const_linop(grid::Grid.EnvGrid, q::Hankel.QDHT, nfun; thg=false)
    n = zero(grid.ω)
    n[grid.sidx] = nfun.(2π*PhysData.c./grid.ω[grid.sidx])
    β1 = PhysData.dispersion_func(1, nfun)(grid.referenceλ)
    if thg
        β0const = 0.0
    else
        β0const = grid.ω0/PhysData.c * nfun(2π*PhysData.c./grid.ω0)
    end
    make_const_linop(grid, q, n, β1, β0const; thg=thg)
end

function make_const_linop(grid::Grid.EnvGrid, q::Hankel.QDHT,
                          n::AbstractArray, β1::Number, β0ref::Number; thg=false)
    out = Array{ComplexF64}(undef, (length(grid.ω), q.N))
    k2 = @. (n*grid.ω/PhysData.c)^2
    kr2 = q.k.^2
    _fill_linop_r!(out, grid, β1, k2, kr2, q.N, β0ref, thg)
    return out
end

"""
    make_linop(grid, q::QDHT, nfun)

Make z-dependent linear operator for radial free-space propagation. `nfun(ω; z)` should
return the refractive index as a function of frequency `ω` and (kwarg) propagation
distance `z`.
"""
function make_linop(grid::Grid.RealGrid, q::Hankel.QDHT, nfun)
    kr2 = q.k.^2
    k2 = zero(grid.ω)
    nfunλ(z) = λ -> nfun(wlfreq(λ), z=z)
    function linop!(out, z)
        β1 = PhysData.dispersion_func(1, nfunλ(z))(grid.referenceλ)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z) .* grid.ω[grid.sidx]./PhysData.c).^2
        _fill_linop_r!(out, grid, β1, k2, kr2, q.N)
    end
end

function _fill_linop_r!(out, grid::Grid.RealGrid, β1, k2, kr2, Nr)
    for ir = 1:Nr
        for iω = 1:length(grid.ω)
            βsq = k2[iω] - kr2[ir]
            if βsq < 0
                # negative βsq -> evanescent fields -> attenuation
                out[iω, ir] = -im*(-β1*grid.ω[iω]) - min(sqrt(abs(βsq)), 200)
            else
                out[iω, ir] = -im*(sqrt(βsq) - β1*grid.ω[iω])
            end
        end
    end
end

function make_linop(grid::Grid.EnvGrid, q::Hankel.QDHT, nfun; thg=false)
    kr2 = q.k.^2
    k2 = zero(grid.ω)
    nfunλ(z) = λ -> nfun(wlfreq(λ), z=z)
    function linop!(out, z)
        β1 = PhysData.dispersion_func(1, nfunλ(z))(grid.referenceλ)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z) .* grid.ω[grid.sidx]./PhysData.c).^2
        βref = thg ? 0.0 : grid.ω0/PhysData.c * nfun(grid.ω0; z=z)
        _fill_linop_r!(out, grid, β1, k2, kr2, q.N, βref, thg)
    end
end

function _fill_linop_r!(out, grid::Grid.EnvGrid, β1, k2, kr2, Nr, βref, thg)
    for ir = 1:Nr
        for iω = 1:length(grid.ω)
            βsq = k2[iω] - kr2[ir]
            if βsq < 0
                # negative βsq -> evanescent fields -> attenuation
                out[iω, ir] = -im*(-β1*grid.ω[iω]) - min(sqrt(abs(βsq)), 200)
            else
                out[iω, ir] = -im*(sqrt(βsq) - β1*grid.ω[iω])
            end
            if !thg
                out[iω, ir] -= -im*βref
            end
        end
    end
end

#=================================================#
#===============   MODE AVERAGE   ================#
#=================================================#

"""
    αlim!(α)

Limit α so that we do not get overflow in exp(α*dz)
"""
function αlim!(α)
    # magic number: this is 130 dB/cm
    # a test script sensitive to this is test_main_rect_env.jl
    clamp!(α, 0.0, 3000.0)
end

"""
    conj_clamp(n, ω)

Simultaneously conjugate and clamp the effective index `n` to safe levels.

The real part is lower-bounded at 1e-3 and the imaginary part upper-bounded at an attenuation
coefficient `α` of 3000 (130 dB/cm). The limits are somewhat arbitrary and chosen empirically
from previous bugs. See https://github.com/LupoLab/Luna/pull/142.

See also [`αlim!`](@ref).
"""
conj_clamp(n, ω) = clamp(real(n), 1e-3, Inf) - im*clamp(imag(n), 0, 3000*PhysData.c/ω)

function make_const_linop(grid::Grid.RealGrid, βfun!, αfun!, β1)
    β = similar(grid.ω)
    βfun!(β, 0)
    α = similar(grid.ω)
    αfun!(α, 0)
    αlim!(α)
    linop = @. -im*(β-β1*grid.ω) - α/2
    linop[.!grid.sidx] .= 0
    return linop
end

function make_const_linop(grid::Grid.EnvGrid, βfun!, αfun!, β1, β0ref)
    β = similar(grid.ω)
    βfun!(β, 0)
    α = similar(grid.ω)
    αfun!(α, 0)
    αlim!(α)
    linop = -im.*(β .- β1.*(grid.ω .- grid.ω0) .- β0ref) .- α./2
    linop[.!grid.sidx] .= 0
    return linop
end

"""
    make_const_linop(grid, mode, λ0)

Make constant linear operator for mode-averaged propagation in mode `mode` with a reference
wavelength `λ0`.
"""
function make_const_linop(grid::Grid.EnvGrid, mode::Modes.AbstractMode, λ0; thg=false)
    β1const = Modes.dispersion(mode, 1, wlfreq(λ0))
    if thg
        β0const = 0.0
    else
        β0const = Modes.β(mode, wlfreq(λ0))
    end
    βconst = zero(grid.ω)
    βconst[grid.sidx] = Modes.β.(mode, grid.ω[grid.sidx])
    βconst[.!grid.sidx] .= 1
    function βfun!(out, z)
        out .= βconst
    end
    αconst = zero(grid.ω)
    αconst[grid.sidx] = Modes.α.(mode, grid.ω[grid.sidx])
    function αfun!(out, z)
        out .= αconst
    end
    make_const_linop(grid, βfun!, αfun!, β1const, β0const), βfun!, β1const, αfun!
end

function make_const_linop(grid::Grid.RealGrid, mode::Modes.AbstractMode, λ0)
    β1const = Modes.dispersion(mode, 1, wlfreq(λ0))
    βconst = zero(grid.ω)
    βconst[grid.sidx] = Modes.β.(mode, grid.ω[grid.sidx])
    βconst[.!grid.sidx] .= 1
    function βfun!(out, z)
        out .= βconst
    end
    αconst = zero(grid.ω)
    αconst[grid.sidx] = Modes.α.(mode, grid.ω[grid.sidx])
    function αfun!(out, z)
        out .= αconst
    end
    make_const_linop(grid, βfun!, αfun!, β1const), βfun!, β1const, αfun!
end

"""
    neff_β_grid(grid, mode, λ0; ref_mode=1)

Create closures which return the effective index and propagation constant
as a function of the frequency grid **index**, rather than the frequency itself.
Any [`Modes.AbstractMode`](@ref) may define its own method for `neff_β_grid` to
accelerate repeated calculation on the same frequency grid.
"""
function neff_β_grid(grid, mode, λ0)
    let grid=grid, mode=mode
        _neff(iω; z) = Modes.neff(mode, grid.ω[iω]; z=z)
        _β(iω; z) = Modes.β(mode, grid.ω[iω]; z=z)
        _neff, _β
    end
end

function make_linop(grid::Grid.RealGrid, mode::Modes.AbstractMode, λ0)
    sidcs = (1:length(grid.ω))[grid.sidx]
    neff, β = neff_β_grid(grid, mode, λ0)
    linop! = let neff=neff, ω=grid.ω, mode=mode, ω0=wlfreq(λ0)
        function linop!(out, z)
            fill!(out, 0.0)
            β1 = Modes.dispersion(mode, 1, ω0, z=z)::Float64
            for iω in sidcs
                nc = conj_clamp(neff(iω; z=z), ω[iω])
                out[iω] = -im*(ω[iω]/PhysData.c*nc - ω[iω]*β1)
            end
        end
    end
    βfun! = let β=β, ω=grid.ω
        function βfun!(out, z)
            fill!(out, 1.0)
            for iω in sidcs
                out[iω] = β(iω; z=z)
            end
        end
    end
    return linop!, βfun!
end

function make_linop(grid::Grid.EnvGrid, mode::Modes.AbstractMode, λ0; thg=false)
    sidcs = (1:length(grid.ω))[grid.sidx]
    neff, β = neff_β_grid(grid, mode, λ0)
    linop! = let neff=neff, ω=grid.ω, mode=mode, ω0=wlfreq(λ0), sidcs=sidcs
        function linop!(out, z)
            fill!(out, 0.0)
            β1 = Modes.dispersion(mode, 1, ω0, z=z)::Float64
            if !thg
                βref = Modes.β(mode, ω0, z=z)
            end
            for iω in sidcs
                nc = conj_clamp(neff(iω; z=z), ω[iω])
                out[iω] = -im*(ω[iω]/PhysData.c*nc - (ω[iω] - grid.ω0)*β1)
                if !thg
                    out[iω] -= -im*βref
                end
            end
        end
    end
    βfun! = let β=β, sidcs=sidcs
        function βfun!(out, z)
            fill!(out, 1.0)
            for iω in sidcs
                out[iω] = β(iω, z=z)
            end
        end
    end
    return linop!, βfun!
end

#=================================================#
#=================   MULTIMODE   =================#
#=================================================#
"""
    make_const_linop(grid, modes, λ0; ref_mode=1)

Make constant (z-invariant) linear operator for multimode propagation. The frame velocity is
taken as the group velocity at wavelength `λ0` in the mode given by `ref_mode` (which
indexes into `modes`)
"""
function make_const_linop(grid::Grid.RealGrid, modes, λ0; ref_mode=1)
    β1 = Modes.dispersion(modes[ref_mode], 1, wlfreq(λ0))
    nmodes = length(modes)
    linops = zeros(ComplexF64, length(grid.ω), nmodes)
    for i = 1:nmodes
        βconst = zero(grid.ω)
        βconst[grid.sidx] = Modes.β.(modes[i], grid.ω[grid.sidx])
        βconst[.!grid.sidx] .= 1
        α = zeros(length(grid.ω))
        α[grid.sidx] .= Modes.α.(modes[i], grid.ω[grid.sidx])
        αlim!(α)
        linops[:,i] = im.*(-βconst .+ grid.ω.*β1) .- α./2
    end
    linops
end

function make_const_linop(grid::Grid.EnvGrid, modes, λ0; ref_mode=1, thg=false)
    β1 = Modes.dispersion(modes[ref_mode], 1, wlfreq(λ0))
    if thg
        βref = 0.0
    else
        βref = Modes.β(modes[ref_mode], wlfreq(λ0))
    end
    nmodes = length(modes)
    linops = zeros(ComplexF64, length(grid.ω), nmodes)
    for i = 1:nmodes
        βconst = zero(grid.ω)
        βconst[grid.sidx] = Modes.β.(modes[i], grid.ω[grid.sidx])
        βconst[.!grid.sidx] .= 1
        α = Modes.α.(modes[i], grid.ω)
        αlim!(α)
        linops[:,i] = -im.*(βconst .- (grid.ω .- grid.ω0).*β1 .- βref) .- α./2
    end
    linops
end

"""
    neff_grid(grid, modes, λ0; ref_mode=1)

Create a closure that returns the effective index as a function of the frequency grid and mode
**index**, rather than the mode and frequency themselves. Any [`Modes.AbstractMode`](@ref)
may define its one method for `neff_grid` to accelerate repeated calculation on the same
frequency grid.
"""
function neff_grid(grid, modes, λ0; ref_mode=1)
    _neff = let grid=grid, modes=modes
        _neff(iω, iim; z) = Modes.neff(modes[iim], grid.ω[iω]; z=z)
    end
    _neff
end

function make_linop(grid::Grid.RealGrid, modes, λ0; ref_mode=1)
    sidcs = (1:length(grid.ω))[grid.sidx]
    neff = neff_grid(grid, modes, λ0; ref_mode=ref_mode)
    linop! = let neff=neff, ω=grid.ω, modes=modes, ω0=wlfreq(λ0), ref_mode=ref_mode
        function linop!(out, z)
            β1 = Modes.dispersion(modes[ref_mode], 1, ω0, z=z)::Float64
            fill!(out, 0.0)
            for i in eachindex(modes)
                for iω in sidcs
                    nc = conj_clamp(neff(iω, i; z=z), ω[iω])
                    out[iω, i] = -im*(ω[iω]/PhysData.c*nc - ω[iω]*β1)
                end
            end
        end
    end
end

function make_linop(grid::Grid.EnvGrid, modes, λ0; ref_mode=1, thg=false)
    sidcs = (1:length(grid.ω))[grid.sidx]
    neff = neff_grid(grid, modes, λ0; ref_mode=ref_mode)
    linop! = let neff=neff, ω=grid.ω, modes=modes, ω0=wlfreq(λ0), ref_mode=ref_mode
        function linop!(out, z)
            β1 = Modes.dispersion(modes[ref_mode], 1, ω0, z=z)::Float64
            fill!(out, 0.0)
            if !thg
                βref = Modes.β(modes[ref_mode], ω0, z=z)
            end
            for i in eachindex(modes)
                for iω in sidcs
                    nc = conj_clamp(neff(iω, i; z=z), ω[iω])
                    out[iω, i] = -im*(ω[iω]/PhysData.c*nc - (ω[iω] - grid.ω0)*β1)
                    if !thg
                        out[iω, i] -= -im*βref
                    end
                end
            end
        end
    end
end

end


text:
@generated function construct_complex_seeds(::Type{Partials{N,V}}) where {N,V}
    return Expr(:tuple, Expr(:tuple,[:(single_seed(Partials{N,V}, Val{2*$i-1}())) for i in 1:N÷2]...),
                        Expr(:tuple,[:(single_seed(Partials{N,V}, Val{2*$i-0}())) for i in 1:N÷2]...))
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x,
               rseed::Partials{M,V} = zero(Partials{M,V}),iseed::Partials{M,V} = zero(Partials{M,V})) where {T,V,N,M}
    for i in eachindex(duals)
        duals[i] = ComplexDual{T,V,M}(x[i],rseed,iseed)
    end
    return duals
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x,
               rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}}) where {T,V,N,M}
    for i in 1:N
        duals[i] = ComplexDual{T,V,M}(x[i],rseeds[i],iseeds[i])
    end
    return duals
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x, index,
               rseed::Partials{M,V} = zero(Partials{M,V}),iseed::Partials{M,V} = zero(Partials{M,V}),chunksize=M÷2) where {T,V,M}
    offset = index - 1
    seed_inds = 1:chunksize
    dual_inds = seed_inds .+ offset
    duals[dual_inds] .= ComplexDual{T,V,M}.(view(x, dual_inds), Ref(rseed),Ref(iseed))
    return duals
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x, index,
               rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}}, chunksize = N) where {T,V,N,M}
    offset = index - 1
    seed_inds = 1:chunksize
    dual_inds = seed_inds .+ offset
    duals[dual_inds] .= ComplexDual{T,V,M}.(view(x, dual_inds), getindex.(Ref(rseeds), seed_inds), getindex.(Ref(iseeds), seed_inds))
    return duals
end

function dualize(::Type{T},x::AbstractArray{<:Complex{V}}) where {T,V}
    xdual = similar(x,ComplexDual{T,V,2*length(v)})
    seed!(xdual,x)
end

function ForwardDiff.vector_mode_dual_eval(f::F, z, cfg::ComplexGradientConfig) where {F}
    zduals = cfg.duals
    seed!(zduals,z,cfg.rseeds,cfg.iseeds)
    return f(zduals)
end

# =====  other seed functions ==== ##

"""
    seed(v::Number,cfg::ComplexGradientConfig)

Given a number `v`, create a dual-valued copy of `v` with
a unit seed. The use of complex duals (even if `v` is real-valued)
ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed(v::Complex,cfg::ComplexGradientConfig)
    vduals = copy(cfg.duals)
    seed!(vduals,[v],cfg.rseeds,cfg.iseeds)
    return vduals[1]
end

function seed(v::Real,cfg::ComplexGradientConfig)
    vduals = copy(cfg.duals)
    seed!(vduals,complex([v]),cfg.rseeds,cfg.iseeds)
    return real(vduals[1])
end


text:
module SparsityOperators

import Base: length, iterate, \
using LinearAlgebra
import LinearAlgebra.BLAS: gemv, gemv!
import LinearAlgebra: BlasFloat, normalize!, norm, rmul!, lmul!
using LinearOperators
using SparseArrays
using Random

using FFTW
using Wavelets

const Trafo = Union{AbstractMatrix, AbstractLinearOperator, Nothing}
const FuncOrNothing = Union{Function, Nothing}

include("FFTOp.jl")
include("DCTOp.jl")
include("DSTOp.jl")
include("WaveletOp.jl")
include("SamplingOp.jl")
include("WeightingOp.jl")

export linearOperator, linearOperatorList

linearOperator(op::Nothing,shape) = nothing

"""
  returns a list of currently implemented `LinearOperator`s
"""
function linearOperatorList()
  return ["DCT-II", "DCT-IV", "FFT", "DST", "Wavelet"]
end

"""
    linearOperator(op::AbstractString, shape)

returns the `LinearOperator` with name `op`.

# valid names
* `"FFT"`
* `"DCT-II"`
* `"DCT-IV"`
* `"DST"`
* `"Wavelet"`
"""
function linearOperator(op::AbstractString, shape)
  shape_ = tuple(shape...)
  if op == "FFT"
    trafo = FFTOp(ComplexF32, shape_, false) #FFTOperator(shape)
  elseif op == "DCT-II"
    shape_ = tuple(shape[shape .!= 1]...)
    trafo = DCTOp(ComplexF32, shape_, 2)
  elseif op == "DCT-IV"
    shape_ = tuple(shape[shape .!= 1]...)
    trafo = DCTOp(ComplexF32, shape_, 4)
  elseif op == "DST"
    trafo = DSTOp(ComplexF32, shape_)
  elseif op == "Wavelet"
    trafo = WaveletOp(shape_)
  else
    error("Unknown transformation")
  end
  trafo
end

end # module


text:
set_problem!(u::Updater, ::Union{POMDP,MDP}) = u

abstract type BehaviorBelief end

mutable struct DiscreteBehaviorBelief <: BehaviorBelief
    physical::MLPhysicalState
    models::AbstractVector
    weights::Vector{Vector{Float64}}
end
DiscreteBehaviorBelief(ps::MLPhysicalState, models::AbstractVector) = DiscreteBehaviorBelief(ps, models, Vector{Vector{Float64}}(undef, length(ps.cars)))

function rand(rng::AbstractRNG,
              b::DiscreteBehaviorBelief,
              s::MLState=MLState(b.physical, Vector{CarState}(undef, length(b.physical.cars))))
    s.crashed = b.physical.crashed
    resize!(s.cars, length(b.physical.cars))
    for i in 1:length(s.cars)
        s.cars[i] = CarState(b.physical.cars[i], sample(rng, b.models, Weights(b.weights[i])))
    end
    return s
end

@with_kw mutable struct WeightUpdateParams
    smoothing::Float64 = 0.02 # value between 0 and 1, adds this fraction of the max to each entry in the vector
    wrong_lane_factor::Float64 = 0.1
end

function weights_from_particles!(b::DiscreteBehaviorBelief, problem::NoCrashProblem, o::MLPhysicalState, particles, p::WeightUpdateParams)
    b.physical = o
    resize!(b.weights, length(o.cars))
    for i in 1:length(o.cars)
        b.weights[i] = zeros(length(b.models)) #XXX lots of allocation
    end
    for sp in particles
        isp = 1
        io = 1
        while io <= length(o.cars) && isp <= length(sp.cars)
            co = o.cars[io]
            csp = sp.cars[isp]
            if co.id == csp.id
                # sigma_acc = dmodel.vel_sigma/dt
                # dv = acc*dt
                # sigma_v = sigma_dv = dmodel.vel_sigma
                if abs(co.x-csp.x) < 0.3*problem.dmodel.phys_param.lane_length
                    proportional_likelihood = proportional_normal_pdf(csp.vel,
                                                                      co.vel,
                                                                      problem.dmodel.vel_sigma)
                    if co.y == csp.y
                        b.weights[io][csp.behavior.idx] += proportional_likelihood
                    elseif abs(co.y - csp.y) < 1.0
                        b.weights[io][csp.behavior.idx] += p.wrong_lane_factor*proportional_likelihood
                    end # if greater than one lane apart, do nothing
                end
                io += 1
                isp += 1
            elseif co.id < csp.id
                io += 1
            else
                @assert co.id > csp.id
                isp += 1
            end
        end
    end

    # smoothing
    for i in 1:length(b.weights)
        if sum(b.weights[i]) > 0.0
            b.weights[i] .+= p.smoothing * sum(b.weights[i])
        else
            b.weights[i] .+= 1.0
        end
    end

    return b
end

mutable struct BehaviorBeliefUpdater <: Updater
    problem::NoCrashProblem
    smoothing::Float64
    nb_particles::Float64
end
function set_problem!(u::BehaviorBeliefUpdater, problem::Union{POMDP,MDP})
    u.problem = problem
end

function update(up::BehaviorBeliefUpdater,
                b_old::DiscreteBehaviorBelief,
                a::MLAction,
                o::MLPhysicalState,
                b_new::DiscreteBehaviorBelief=DiscreteBehaviorBelief(o,
                                                       up.problem.dmodel.behaviors,
                                                       Vector{Vector{Float64}}(undef, 0)))
    # particles = SharedArray(MLState, )
    # @parallel
    # # run simulations

    # weights_from_particles!(b_new, )
    error("not implemented")
end

# TODO clean up
#=
import DESPOT: fringe_upper_bound, lower_bound, init_lower_bound

function fringe_upper_bound(pomdp::NoCrashPOMDP, s::MLState)
  """
  Assume there are no cars in the environment: just go straight to target lane, and stay there
    for the rest of the time horizon
  """
  # NOTE: This function makes assumptions about the length of the time horizon!!!!
  T = 100 #time horizon
  pp = pomdp.dmodel.phys_param
  nb_lanes = pp.nb_lanes
  current_lane = s.cars[1].y
  t_single_lane_change = 1.0 / (pp.dt * pomdp.dmodel.lane_change_rate) #nb timesteps
  t_target_lane = (nb_lanes - current_lane) * t_single_lane_change
  return sum([pomdp.rmodel.reward_in_target_lane * discount(pomdp)^t for t = 0:(T-t_target_lane-1)])
end

type NoCrashLB <: DESPOTLowerBound end

init_lower_bound(lb::NoCrashLB, pomdp::NoCrashPOMDP, config::DESPOTConfig) = begin end

function lower_bound(lb::NoCrashLB, pomdp::NoCrashPOMDP, particles::Vector{DESPOTParticle{MLState}}, config:DESPOTConfig)
  """
  Just rollout the heuristic policy--it only really depends on observations, so
    bmodel (eg true behavior states) don't matter
  """

  s = particles[1].state

  sim = RolloutSimulator(initial_state=s, max_steps=100)

  policy = Simple(pomdp) # TODO check if this signature works with NoCrashPOMDP and stuff

  return simulate(sim, pomdp, policy)
end
=#


text:
abstract type AbstractBlockDiagonal{M,N,T} <: AbstractMatrix{T} end

struct BlockDiagonalColumnView{M,N,T,P,L,A<:AbstractFixedSizeMatrix{M,N,T,P,L}} <: AbstractBlockDiagonal{M,N,T}
    data::A#ConstantFixedSizeMatrix{M,N,T,P,L}
end
# function Base.getindex(bd::BlockDiagonalColumnView{M,N,T,P,L}) where {M,N,T,P,L}
#
# end
Base.size(::BlockDiagonalColumnView{M,N}) where {M,N} = (M*N,N)
@inline VectorizationBase.vectorizable(A::BlockDiagonalColumnView) = VectorizationBase.vectorizable(A.data)
function Base.getindex(A::BlockDiagonalColumnView{M,N,T}, i::Int, j::Int) where {M,N,T} # primary purpse is supporting show methods
    @boundscheck begin
        (i < 1 || i > (M*N) || j < 1 || j > N) && throw(BoundsError("Attempted to access index ($i,$j) of size $(size(A)) BlockDiagonalColumnView."))
    end
    i -= (j-1) * M
    if i < 1 || i > M
        return zero(T)
    end
    @inbounds A.data[i, j]
end

function block_diagonal_column_view_quote(M,N,T,PA,PB,increment::Bool = false)
    unroll = N > 7 ? 4 : N
    if N == unroll
        loop_body = quote end
        for n ∈ 0:N-1
            push!(loop_body.args, :($(Symbol(:s_,n+1)) += A[m + $(PA*n) ] * BD[m + $(PB*n)]))
        end
        #        return quote
        q = quote
            Base.Cartesian.@nexprs $unroll i -> s_i = zero($T)
            @vectorize $T for m ∈ 1:$M #$P # Can we assume the padding is uncontaminated?
                $loop_body
            end
        end
        if increment
            push!(q.args, quote
                  @inbounds begin
                  Base.Cartesian.@nexprs $unroll i -> c[i] = d[i] + s_i
                  end
                  end)
        else
            push!(q.args, quote
                  @inbounds begin
                  Base.Cartesian.@nexprs $unroll i -> c[i] = s_i
                  end
                  end)
        end
        push!(q.args, :c)
#        println(q)
        return q
    end
    reps, rem = divrem(N, unroll)
    if rem > 0
        # catch cases such as N = 9,
        # where we can evenly split the iterations 3 x 3
        # or N = 10, where iterations are 5 x 2
        reps_candidate, rem_candidate = divrem(N, unroll - 1)
        if rem_candidate == 0
            unroll -= 1
            reps, rem = reps_candidate, rem_candidate
        else
            reps_candidate, rem_candidate = divrem(N, unroll - 2)
            if rem_candidate == 0
                unroll -= 2
                reps, rem = reps_candidate, rem_candidate
            end
        end
    end
    loop_body = quote end
    for n ∈ 0:unroll-1
        push!(loop_body.args, :($(Symbol(:s_,n+1)) += A[m + $(unroll*PA)*n + $(PA*n) ] * BD[m + $(unroll*PB)*n + $(PB*n)]))
    end
    q = quote
        for n ∈ 0:$((N÷unroll)-1)
            Base.Cartesian.@nexprs $unroll i -> s_i = zero($T)
            @vectorize $T for m ∈ 1:$M #$P # Can we assume the padding is uncontaminated?
                $loop_body
            end
        end
        if increment
            push!(q.args, quote
                  @inbounds begin
                  Base.Cartesian.@nexprs $unroll i -> c[i + $unroll*n] = d[i + $unroll*n] + s_i
                  end
                  end)
        else
            push!(q.args, quote
                  @inbounds begin
                  Base.Cartesian.@nexprs $unroll i -> c[i + $unroll*n] = s_i
                  end
                  end)
        end
    end
    if rem > 0
        loop_body = quote end
        assignments = quote end
        for n ∈ N-rem+1:N
            ssym = Symbol(:s_,n)
            push!(q.args, :($ssym = zero($T)))
            push!(loop_body.args, :($ssym += A[m + $((n-1)*PA) ] * BD[m + $((n-1)*PB)]))
            if increment
                push!(assignments.args, :( c[$n] = d[$n] + $ssym ) )
            else
                push!(assignments.args, :( c[$n] = $ssym ) )
            end
        end
        push!(q.args, quote
            @vectorize $T for m ∈ 1:$M
                $loop_body
            end
            @inbounds begin
                $assignments
            end
        end)
    end
 #   push!(q.args, :(c))
    q

end

@generated function LinearAlgebra.mul!(
    c::PaddedMatrices.AbstractMutableFixedSizeVector{N,T},
    A::AbstractFixedSizeMatrix{M,N,T,PA},
    BD::BlockDiagonalColumnView{M,N,T,PB}
) where {M,N,T,PA,PB}
    quote
        $(block_diagonal_column_view_quote(M,N,T,PA,PB,false))
        c
    end
end
function Base.:*(
    A::AbstractFixedSizeMatrix{M,N,T,PA},
    BD::BlockDiagonalColumnView{M,N,T,PB}
) where {M,N,T,PA,PB}
    c  = FixedSizeVector{N,T}(undef)
    mul!(c, A, BD)'
end
@generated function Base.:*(
    sp::PaddedMatrices.StackPointer,
    A::AbstractFixedSizeMatrix{M,N,T,PA},
    BD::BlockDiagonalColumnView{M,N,T,PB}
) where {M,N,T,PA,PB}
    P = min(PA,PB)
    quote
        c  = PtrVector{$N,$T,$N}(pointer(sp,$T))
        sp + $(sizeof(T)*N), mul!(c, A, BD)'
    end
end
@generated function ReverseDiffExpressionsBase.RESERVED_INCREMENT_SEED_RESERVED!(
    d::PaddedMatrices.AbstractMutableFixedSizeVector{N,T,PC},
    BD::BlockDiagonalColumnView{M,N,T,PB},
    A::AbstractFixedSizeMatrix{M,N,T,PA}
) where {M,N,T,PA,PB,PC}
    quote
        c = PtrVector{$N,$T,$PC}(pointer(d))
        $(block_diagonal_column_view_quote(M,N,T,PA,PB,true))
        nothing
    end
end
# @generated function PaddedMatrices.RESERVED_INCREMENT_SEED_RESERVED(
    # sp::PaddedMatrices.StackPointer,
    # A::AbstractFixedSizeMatrix{M,N,T,PA},
    # BD::BlockDiagonalColumnView{M,N,T,PB},
    # d′::LinearAlgebra.Adjoint{T,<:PaddedMatrices.AbstractMutableFixedSizeVector{N,T,PC}}
# ) where {M,N,T,PA,PB,PC}
    # quote
        # d = d′'
        # c = PtrVector{$N,$T,$N}(pointer(sp,$T))
        # $(block_diagonal_column_view_quote(M,N,T,PA,PB,true))
        # sp + $(N*sizeof(T)), c'
    # end
# end


text:
using MotifScanner
using Test

@testset "MotifScanner.jl" begin
    # Write your tests here.
end


text:
abstract type AbstractPosteriors end

#holds the posterior distributions
mutable struct Posteriors <: AbstractPosteriors
  σ²G!::Function
  Z!::Function
  SG!::Function
  SGP!::Function
end

#NOTE: Each posterior distribution as a corresponding test version

function Posteriors(;
  σ²G!::Function = posteriorσ²G!,
  Z!::Function = posteriorZ!,
  SG!::Function = posteriorSG!,
  SGP!::Function = posteriorSGP!)::Posteriors

  return Posteriors(σ²G!, Z!, SG!, SGP!)
end

function PosteriorsTest(;
  σ²G!::Function = posteriorσ²Gtest!,
  Z!::Function = posteriorZtest!,
  SG!::Function = posteriorSGtest!,
  SGP!::Function = posteriorSGPtest!)::Posteriors

  return Posteriors(σ²G!, Z!, SG!, SGP!)
end

#############
function posteriorσ²G!(Θ::Parameters, Π::Priors)::Nothing

  μG::Float64 = ((-Θ.SGP)*(-Θ.ζ²G + Θ.ζGP) - Θ.SG*(Θ.ζGP - Θ.ζ²P))/(Θ.ζ²G - 2*Θ.ζGP + Θ.ζ²P)
  γG2::Float64 = (-Θ.ζGP^2 + Θ.ζ²G*Θ.ζ²P)/(Θ.ζ²G - 2*Θ.ζGP + Θ.ζ²P)
  ζ²Gstar::Float64 = 1/(1/γG2+1/Π.δ²G)
  μ::Float64 = ζ²Gstar * (μG/γG2+Π.θG/Π.δ²G)

  ζ²Gstar = max(ζ²Gstar, SMALL)

  try
    Θ.σ²G = rand(TruncatedNormal(μ,ζ²Gstar, 0.0, ∞))
  catch err
    println("ζ²Gstar: $ζ²Gstar, μ=$μ, Θ: $Θ")
    error(err)
  end

  return nothing
end

###########PosteriorsZ
function posteriorZ!(Θ::Parameters, Π::Priors)::Nothing

  s11::Float64 = (Θ.SG - Θ.σ²G)^2
  s12::Float64 = (Θ.SGP - Θ.σ²G)*(Θ.SG - Θ.σ²G)
  s22::Float64 = (Θ.SGP - Θ.σ²G)^2
  Σ::Matrix{Float64} = [s11 s12; s12 s22]

  Z::Matrix{Float64} = rand(InverseWishart(Π.ν + 1, Π.Ψ .+ Σ))
  Θ.ζ²G = Z[1,1]
  Θ.ζGP = Z[1,2]
  Θ.ζ²P = Z[2,2]
  return nothing
end

#################
function posteriorSG!(Θ::Parameters, Π::Priors)::Nothing
  μSG = Θ.σ²G+(Θ.ζGP/Θ.ζ²P)*(Θ.SGP-Θ.σ²G)
  γ²SG = Θ.ζ²G-(Θ.ζGP^2/Θ.ζ²P)
  ζ²SGstar::Float64 = 1.0 / (1.0/γ²SG+ 1.0/Π.δ²SG)

  μ::Float64 = (μSG/γ²SG + Π.θSG/Π.δ²SG) * ζ²SGstar

  ζ²SGstar = max(ζ²SGstar, SMALL)

  Θ.SG = rand(TruncatedNormal(μ,ζ²SGstar, 0.0, ∞))

  return nothing
end

##################
function posteriorSGP!(Θ::Parameters, Π::Priors)::Nothing
  μSGP = Θ.σ²G+(Θ.ζGP/Θ.ζ²G)*(Θ.SG-Θ.σ²G)
  γ²SGP = Θ.ζ²P-(Θ.ζGP^2/Θ.ζ²G)
  ζ²SGPstar::Float64 = 1.0 / (1.0/γ²SGP+ 1.0/Π.δ²SGP)

  μ::Float64 = (μSGP/γ²SGP + Π.θSGP/Π.δ²SGP) * ζ²SGPstar

  #println("SGTest## μSG: $μSGP  γ²SG: $γ²SGP ζ²SGstar: $ζ²SGPstar")
  ζ²SGPstar = max(ζ²SGstar, SMALL)
  Θ.SGP = rand(TruncatedNormal(μ, ζ²SGPstar, 0.0, ∞))

  return nothing
end

##############****************************#################################
##############****************************#################################
##############****************************#################################
##############****************************#################################

function posteriorσ²Gtest!(Θ::Parameters, Π::Priors)::Nothing

  μG::Float64 = ((-Θ.SGP)*(-Θ.ζ²G + Θ.ζGP) - Θ.SG*(Θ.ζGP - Θ.ζ²P))/(Θ.ζ²G - 2*Θ.ζGP + Θ.ζ²P)
  γG2::Float64 = (-Θ.ζGP^2 + Θ.ζ²G*Θ.ζ²P)/(Θ.ζ²G - 2*Θ.ζGP + Θ.ζ²P)

  ζ²Gstar::Float64 = 1/(1/γG2+1/Π.δ²G)
  μ::Float64 = ζ²Gstar * (μG/γG2+Π.θG/Π.δ²G)

  Θ.σ²G = rand(TruncatedNormal(μ, ζ²Gstar, 0.0, ∞))
  return nothing
end

function posteriorZtest!(Θ::Parameters, Π::Priors)::Nothing

  Σ::Matrix{Float64} =
    [(Θ.SG-Θ.σ²G)^2 (Θ.SG-Θ.σ²G)*(Θ.SGP-Θ.σ²G); (Θ.SG-Θ.σ²G)*(Θ.SGP-Θ.σ²G) (Θ.SGP-Θ.σ²G)^2]

  Z::Matrix{Float64} = rand(InverseWishart(Π.ν + 1, Π.Ψ + Σ))
  Θ.ζ²G = Z[1,1]
  Θ.ζGP = Z[1,2]
  Θ.ζ²P = Z[2,2]
  return nothing
end

function posteriorSGtest!(Θ::Parameters, Π::Priors)::Nothing
  ρSG::Float64 = Θ.ζGP/(Θ.ζ²G*Θ.ζ²P)^0.5
  μSG::Float64 = Θ.σ²G+(Θ.ζ²G/Θ.ζ²P)^0.5*ρSG*(Θ.SGP-Θ.σ²G)
  γ²SG::Float64 = (1-ρSG^2)*Θ.ζ²G

  ζ²SGstar::Float64 = 1/(1/γ²SG+1/Π.δ²SG)
  μ::Float64 = ζ²SGstar*(μSG/γ²SG+Π.θSG/Π.δ²SG)

  Θ.SG = rand(TruncatedNormal(μ, ζ²SGstar, 0.0, ∞))

  return nothing
end

function posteriorSGPtest!(Θ::Parameters, Π::Priors)::Nothing
  ρSGP::Float64 = Θ.ζGP/(Θ.ζ²G*Θ.ζ²P)^0.5
  μSGP::Float64 = Θ.σ²G+(Θ.ζ²P/Θ.ζ²G)^0.5*ρSGP*(Θ.SG-Θ.σ²G)
  γ²SGP::Float64 = (1-ρSGP^2)*Θ.ζ²P

  ζ²SGPstar::Float64 = 1/(1/γ²SGP+1/Π.δ²SGP)
  μ::Float64 = ζ²SGPstar*(μSGP/γ²SGP+Π.θSGP/Π.δ²SGP)

  #println("SGTest## μSG: $μSGP  γ²SG: $γ²SGP ζ²SGstar: $ζ²SGPstar")

  Θ.SGP = rand(TruncatedNormal(μ, ζ²SGPstar, 0.0, ∞))

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Changepoints
using Test

# write your own tests here
include("test_misc.jl")
include("test_ChangepointSampler.jl")
include("test_PELT.jl")
include("test_BS.jl")
include("test_R_comparison.jl")

# Testing code need to be verified
# include("test_CROPS.jl")


text:
function marAIC(m::LinearMixedModel)

    rho = length(MixedModels.θ(m))+length(fixef(m))+1
    llk = -objective(m)/2.0

    cr=[1.0, 0.0]

   val=-2.0*llk+2.0*(cr[1]*rho+cr[2])

    IC(llk, rho, "MarAIC", cr, val  )
end

function marBIC(m::LinearMixedModel)

    rho = length(MixedModels.θ(m))+length(fixef(m))+1
    llk = -objective(m)/2.0
    n,p,q,k = size(m)
    #(nn,pp,qq,tt)=size(m)

    cr=[1.0, 0.0]

    val=-2.0*llk+log(n)*(cr[1]*rho+cr[2])

    IC(llk, rho, "MarBIC", cr, val )
end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
@withkw @mixin immutable NoisyIFParameter <: IFParameter
    σ::Float = 0
end

@withkw @mixin type NoisyIF <: IF
    param::NoisyIFParameter = NoisyIFParameter()
    randncache::Vector{Float} = randn(N)
end

@replace function integrate!(p::NoisyIF, param::NoisyIFParameter, dt::Float)
    randn!(randncache)
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i] + σ / √dt * randncache[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
using NonNegLeastSquares

## compare two algorithms
function compare(alg1,alg2;var1=:none,var2=:none,m=200,n=200,k=3)
	println("\nComparing $alg1:$var1 to $alg2:$var2 with A = randn($m,$k) and B = randn($m,$n)")
	println("-------------------------------------------------------------------------------------")

	A,B = randn(m,k),randn(m,n)

	print(uppercase(string(alg1)),":",var1," → ")
	@time X1 = nonneg_lsq(A,B,alg=alg1,variant=var1)

	print(uppercase(string(alg2)),":",var2," → ")
	@time X2 = nonneg_lsq(A,B,alg=alg2,variant=var2)

	if vecnorm(X1-X2)/sqrt(length(X1)) > 1e-10
		warn("algorithms did not converge on the same answer.")
	end

end

compare(:nnls, :fnnls; m=100, n=1, k=100);

compare(:fnnls, :pivot; m=1000, n=1, k=1000);

compare(:fnnls, :admm; m=1000, n=1, k=1000);

compare(:pivot, :admm; m=1000, n=1, k=1000);

compare(:fnnls, :pivot; m=1000, n=1000, k=3);

compare(:pivot, :pivot; var2=:cache, m=100, n=100, k=100);

compare(:pivot, :pivot; var2=:comb, m=100, n=100, k=100);

compare(:pivot, :pivot; var1=:cache, var2=:comb, m=100, n=100, k=100);

compare(:pivot, :pivot; var1=:none, var2=:comb, m=5000, n=5000, k=3);

compare(:pivot, :pivot; var1=:cache, var2=:comb, m=10000, n=10000, k=10);


text:
"""
Atomic solution for an SDE.

### Parameters

* `DT`: data type
* `TT`: time step type
* `AT`: array type
* `IT`: internal variable types

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `ΔW`: Wiener process driving the stochastic process q
* `ΔZ`: Wiener process driving the stochastic process q
* `K`:  integer parameter defining the truncation of the increments of the Wiener process (for strong solutions)
"""
mutable struct AtomicSolutionSDE{DT <: Number, TT <: Real, AT <: AbstractArray{DT}, IT <: NamedTuple} <: AtomicSolution{DT,TT,AT}
    t::TT
    t̅::TT

    q::AT
    q̅::AT
    q̃::AT

    ΔW::AT
    ΔZ::AT

    K::Int

    internal::IT

    function AtomicSolutionSDE{DT,TT,AT,IT}(nd, nm, internal::IT) where {DT,TT,AT,IT}
        new(zero(TT), zero(TT),
            AT(zeros(DT, nd)), AT(zeros(DT, nd)), AT(zeros(DT, nd)),
            AT(zeros(DT, nm)), AT(zeros(DT, nm)), 0,
            internal)
    end

    function AtomicSolutionSDE{DT,TT,AT,IT}(t::TT, q::AT, ΔW::AT, ΔZ::AT, internal::IT) where {DT,TT,AT,IT}
        new(zero(t), zero(t),
            zero(q), zero(q), zero(q),
            zero(ΔW), zero(ΔZ), 0,
            internal)
    end
end

AtomicSolutionSDE(::Type{DT}, ::Type{TT}, ::Type{AT}, nd, nm, internal::IT=NamedTuple()) where {DT,TT,AT,IT} = AtomicSolutionSDE{DT,TT,AT,IT}(nd, nm, internal)
AtomicSolutionSDE(t::TT, q::AT, ΔW::AT, ΔZ::AT, internal::IT=NamedTuple()) where {DT, TT, AT <: AbstractArray{DT}, IT} = AtomicSolutionSDE{DT,TT,AT,IT}(t, q, ΔW, ΔZ, internal)

function Solutions.set_solution!(asol::AtomicSolutionSDE, sol)
    t, q = sol
    asol.t  = t
    asol.q .= q
end

function set_increments!(asol::AtomicSolutionSDE, incs)
    ΔW, ΔZ = incs
    asol.ΔW .= ΔW
    asol.ΔZ .= ΔZ
end

function Solutions.get_solution(asol::AtomicSolutionSDE)
    (asol.t, asol.q)
end

function get_increments(asol::AtomicSolutionSDE)
    (asol.ΔW, asol.ΔZ)
end

function get_increments!(asol::AtomicSolutionSDE, ΔW, ΔZ)
    ΔW .= asol.ΔW
    ΔZ .= asol.ΔZ
end

function Common.reset!(asol::AtomicSolutionSDE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.t += Δt
end

function Solutions.update!(asol::AtomicSolutionSDE{DT}, y::Vector{DT}) where {DT}
    for k in eachindex(y)
        update!(asol, y[k], k)
    end
end

function Solutions.update!(asol::AtomicSolutionSDE{DT}, y::DT, k::Union{Int,CartesianIndex}) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(y, asol.q[k], asol.q̃[k])
end


text:
import Parameters.@with_kw
import Statistics.mean
import Statistics.std
import Statistics.norm

"""
"""
@with_kw mutable struct ScaleMinMax
    data_min::Float64 = 0.0
    data_max::Float64 = 1.0
    is_fit::Bool = false
end

"""
"""
function fit_transform!(scale::ScaleMinMax, data::AbstractArray{<:Real})
    scale.data_min = minimum(Float64,data)
    scale.data_max = maximum(Float64,data)
    scale.is_fit = true
    scale
end

"""
"""
function apply_transform(scale::ScaleMinMax, data::AbstractArray{<:Real})
    if !scale.is_fit
        error("`scale` hasn't been fit")
    end
    span = scale.data_max - scale.data_min
    (data .- scale.data_min) ./ span
end

"""
"""
@with_kw mutable struct ScaleVariance
    data_mean::Float64 = 0.0
    data_std::Float64 = 1.0
    is_fit::Bool = false
end

"""
"""
function fit_transform!(scale::ScaleVariance, data::AbstractArray{<:Real})
    scale.data_mean = mean(data)
    scale.data_std = std(data)
    scale.is_fit = true
    scale
end

"""
"""
function apply_transform(scale::ScaleVariance, data::AbstractArray{<:Real})
    if !scale.is_fit
        error("`scale` hasn't been fit")
    end
    (data .- scale.data_mean) ./ scale.data_std
end

"""
"""
@with_kw mutable struct ScaleL2
    data_norm::Float64 = 1.0
    is_fit::Bool = false
end

"""
"""
function fit_transform!(scale::ScaleL2, data::AbstractArray{<:Real})
    scale.data_norm = norm(data)
    scale.is_fit = true
    scale
end

"""
"""
function apply_transform(scale::ScaleL2, data::AbstractArray{<:Real})
    if !scale.is_fit
        error("`scale` hasn't been fit")
    end
    data ./ scale.data_norm
end


text:
@testset "broadcast" begin
  @test testf((x)       -> fill!(x, 1),  rand(3,3))
  @test testf((x, y)    -> map(+, x, y), rand(2, 3), rand(2, 3))
  @test testf((x)       -> sin.(x),      rand(2, 3))
  @test testf((x)       -> log.(x) .+ 1, rand(2, 3))
  @test testf((x)       -> 2x,           rand(2, 3))
  @test testf((x)       -> x .^ 0,      rand(2, 3))
  @test testf((x)       -> x .^ 1,      rand(2, 3))
  @test testf((x)       -> x .^ 2,      rand(2, 3))
  @test testf((x)       -> x .^ 3,      rand(2, 3))
  @test testf((x)       -> x .^ 5,      rand(2, 3))
  @test testf((x)       -> (z = Int32(5); x .^ z),      rand(2, 3))
  @test testf((x)       -> (z = Float64(π); x .^ z),      rand(2, 3))
  @test testf((x)       -> (z = Float32(π); x .^ z),      rand(Float32, 2, 3))
  @test testf((x, y)    -> x .+ y,       rand(2, 3), rand(1, 3))
  @test testf((z, x, y) -> z .= x .+ y,  rand(2, 3), rand(2, 3), rand(2))
  @test (CuArray{Ptr{Cvoid}}(undef, 1) .= C_NULL) == CuArray([C_NULL])
  @test CuArray([1,2,3]) .+ CuArray([1.0,2.0,3.0]) == CuArray([2,4,6])

  @eval struct Whatever{T}
      x::Int
  end
  @test Array(Whatever{Int}.(CuArray([1]))) == Whatever{Int}.([1])
end
# https://github.com/JuliaGPU/CUDA.jl/issues/223
@testset "Ref Broadcast" begin
  foobar(idx, A) = A[idx]
  @test CuArray([42]) == foobar.(CuArray([1]), Base.RefValue(CuArray([42])))
end

@testset "Broadcast Fix" begin
  @test testf(x -> log.(x), rand(3,3))
  @test testf((x,xs) -> log.(x.+xs), Ref(1), rand(3,3))
end


text:
using RecipesBase
using GraphRecipes

export plotabm, plotabm!

mutable struct PlotABM{AbstractSpace}
    args::Any
end

"""
    plotabm(model::ABM{<: ContinuousSpace}; ac, as, am, kwargs...)
    plotabm(model::ABM{<: DiscreteSpace}; ac, as, am, kwargs...)

Plot the `model` as a `scatter`-plot, by configuring the agent shape, color and size
via the keywords `ac, as, am`.
These keywords can be constants, or they can be functions, each accepting an agent
and outputting a valid value for color/shape/size.

The keyword `scheduler = model.scheduler` decides the plotting order of agents
(which matters only if there is overlap).

The keyword `offset` is a function with argument `offest(a::Agent)`.
It targets scenarios where multiple agents existin within
a grid cell as it adds an offset (same type as `agent.pos`) to the plotted agent position.

All other keywords are propagated into `Plots.scatter` and the plot is returned.

    plotabm(model::ABM{<: GraphSpace}; ac, as, am, kwargs...)
This function is the same as `plotabm` for `ContinuousSpace`, but here the three key
functions `ac, as, am` do not get an agent as an input but a vector of agents at
each node of the graph. Their output is the same.

Here `as` defaults to `length`. Internally, the `graphplot` recipe is used, and
all other `kwargs...` are propagated there.
"""
plotabm(args...; kw...) = RecipesBase.plot(PlotABM{typeof(args[1].space)}(args); kw...)

"""
    plotabm!(model)
    plotabm!(plt, model)

Functionally the same as [`plotabm`](@ref), however this method appends to the active
plot, or one identified as `plt`.
"""
plotabm!(args...; kw...) = RecipesBase.plot!(PlotABM{typeof(args[1].space)}(args); kw...)
plotabm!(plt::RecipesBase.AbstractPlot, args...; kw...) =
    RecipesBase.plot!(plt, PlotABM{typeof(args[1].space)}(args); kw...)

@recipe function f(
    h::PlotABM{<:Union{GridSpace,ContinuousSpace}};
    scheduler = nothing,
    offset = nothing,
    ac = "#765db4",
    as = 10,
    am = :circle,
)
    if length(h.args) != 1 || !(typeof(h.args[1]) <: ABM)
        error("plotabm should be given a model::ABM.  Got: $(typeof(h.args))")
    end

    model = h.args[1]
    if scheduler == nothing
        scheduler = model.scheduler
    end

    ids = scheduler(model)
    colors = typeof(ac) <: Function ? [ac(model[i]) for i in ids] : ac
    sizes = typeof(as) <: Function ? [as(model[i]) for i in ids] : as
    markers = typeof(am) <: Function ? [am(model[i]) for i in ids] : am

    if offset == nothing
        pos = [model[i].pos for i in ids]
    else
        pos = [model[i].pos .+ offset(model[i]) for i in ids]
    end
    x := first.(pos)
    y := last.(pos)

    seriestype := :scatter
    markercolor := colors
    markersize := sizes
    markershape := markers
    legend --> false
    markerstrokewidth --> 0.5
    markerstrokecolor --> :black
    ()
end

@recipe function f(
    h::PlotABM{<:GraphSpace};
    ac = x -> "#765db4",
    as = length,
    am = x -> :circle,
)
    if length(h.args) != 1 || !(typeof(h.args[1]) <: ABM)
        error("plotabm should be given a model::ABM.  Got: $(typeof(h.args))")
    end
    model = h.args[1]
    N = nodes(model)
    ncolor = Vector(undef, length(N))
    weights = zeros(length(N))
    markers = Vector(undef, length(N))
    for (i, n) in enumerate(N)
        a = get_node_agents(n, model)
        ncolor[i] = ac(a)
        weights[i] = as(a)
        markers[i] = am(a)
    end

    node_weights := weights
    nodeshape := markers
    nodecolor := ncolor
    seriescolor --> :black
    markerstrokecolor --> :black
    markerstrokewidth --> 1.5
    model.space.graph
end


text:
module BioMotifInference
    using BioBackgroundModels, BioSequences, Distributed, Distributions, Serialization, UnicodePlots
    import DataFrames:DataFrame
    import ProgressMeter: AbstractProgress, Progress, @showprogress, next!, move_cursor_up_while_clearing_lines, printover, durationstring
    import Printf: @sprintf
    import StatsFuns: logaddexp, logsumexp, logsubexp
    import Random: rand, seed!, shuffle!
    import Distances: euclidean
    import Measurements: measurement

    #CONSTANTS AND PERMUTE FUNCTION ARGUMENT DEFAULTS GIVING RISE TO IMPLEMENTATION-SPECIFIC SAMPLING EFFECTS
    global MOTIF_EXPECT=1. #motif expectation- this value to be divided by obs lengths to obtain penalty factor for scoring

    global REVCOMP=true #calculate scores on both strands?

    global TUNING_MEMORY=20 #coefficient multiplied by Permute_Instruct function call limit to give total number of calls remembered by tuner per function
    global CONVERGENCE_MEMORY=500 #number of iterates to display for convergence interval history

    global SRC_PERM_FREQ=.5 #frequency with which random_decorrelate will permute a source

    global PWM_SHIFT_DIST=Weibull(.5,.1) #distribution of weight matrix permutation magnitudes
    global PWM_SHIFT_FREQ=.2 #proportion of positions in source to permute weight matrix
    global PWM_LENGTHPERM_FREQ=.2 #proportion of sources to permute length
    global LENGTHPERM_RANGE=1:3

    global PRIOR_WT=3. #estimate prior dirichlets from product of this constant and sample "mle" wm
    global PRIOR_LENGTH_MASS=.8

    global EROSION_INFO_THRESH=1.

    global CONSOLIDATE_THRESH=.035

    include("IPM/ICA_PWM_Model.jl")
    export Model_Record
    export ICA_PWM_Model
    include("IPM/IPM_likelihood.jl")
    include("IPM/IPM_prior_utilities.jl")
    assemble_source_priors
    include("ensemble/IPM_Ensemble.jl")
    export IPM_Ensemble, assemble_IPMs
    include("permutation/permute_utilities.jl")
    include("permutation/orthogonality_helper.jl")
    include("permutation/permute_functions.jl")
    export full_perm_funcvec
    include("permutation/permute_control.jl")
    export Permute_Instruct
    include("permutation/Permute_Tuner.jl")
    include("ensemble/ensemble_utilities.jl")
    export ensemble_history, reset_ensemble!, move_ensemble!, copy_ensemble!, rewind_ensemble, complete_evidence, get_model, show_models, reestimate_ensemble!
    include("utilities/model_display.jl")
    include("utilities/worker_diagnostics.jl")
    include("utilities/ns_progressmeter.jl")
    include("utilities/synthetic_genome.jl")
    include("utilities/worker_sequencer.jl")
    export synthetic_sample
    include("nested_sampler/nested_step.jl")
    include("nested_sampler/converge_ensemble.jl")
    export converge_ensemble!

end # module


text:
"""
    function elephant_in_snake()

Do we need to explain this?

# Example

```jldoctests
julia> cowsay("Help!", cow=Cowsay.elephant_in_snake)
 _______
< Help! >
 -------
   \\
    \\              ....
           ........    .
          .            .
         .             .
.........              .......
..............................

Elephant inside ASCII snake

```
"""
function elephant_in_snake(;eyes="oo", tongue="  ", thoughts="\\")
the_cow = """
   $thoughts
    $thoughts              ....
           ........    .
          .            .
         .             .
.........              .......
..............................

Elephant inside ASCII snake
"""
return the_cow
end


text:
using Base.Broadcast: broadcastable
using ChainRules
using ChainRulesCore
using ChainRulesTestUtils
using ChainRulesTestUtils: rand_tangent, _fdm
using Compat: hasproperty, only
using FiniteDifferences
using FiniteDifferences: rand_tangent
using SpecialFunctions
using LinearAlgebra
using LinearAlgebra.BLAS
using LinearAlgebra: dot
using Random
using Statistics
using Test

Random.seed!(1) # Set seed that all testsets should reset to.

function include_test(path)
    print("Testing $path:\t")  # print so TravisCI doesn't timeout due to no output
    @time include(path)  # show basic timing, (this will print a newline at end)
end

println("Testing ChainRules.jl")
@testset "ChainRules" begin
    @testset "rulesets" begin
        @testset "Base" begin
            include_test("rulesets/Base/base.jl")
            include_test("rulesets/Base/fastmath_able.jl")
            include_test("rulesets/Base/evalpoly.jl")
            include_test("rulesets/Base/array.jl")
            include_test("rulesets/Base/arraymath.jl")
            include_test("rulesets/Base/indexing.jl")
            include_test("rulesets/Base/mapreduce.jl")
        end
        println()

        @testset "Statistics" begin
            include_test("rulesets/Statistics/statistics.jl")
        end
        println()

        @testset "LinearAlgebra" begin
            include_test("rulesets/LinearAlgebra/dense.jl")
            include_test("rulesets/LinearAlgebra/norm.jl")
            include_test("rulesets/LinearAlgebra/matfun.jl")
            include_test("rulesets/LinearAlgebra/structured.jl")
            include_test("rulesets/LinearAlgebra/symmetric.jl")
            include_test("rulesets/LinearAlgebra/factorization.jl")
            include_test("rulesets/LinearAlgebra/blas.jl")
        end
        println()

        @testset "Random" begin
            include_test("rulesets/Random/random.jl")
        end
        println()

        @testset "packages" begin
            include_test("rulesets/packages/NaNMath.jl")
            # Note: drop SpecialFunctions dependency in next breaking release
            # https://github.com/JuliaDiff/ChainRules.jl/issues/319
            if !isdefined(SpecialFunctions, :ChainRulesCore)
                include_test("rulesets/packages/SpecialFunctions.jl")
            end
        end
        println()
    end
end


text:
function somefunc()
  42
end

max_stars_repo_name: org-arl/jajub

max_issues_repo_name: org-arl/jajub

max_forks_repo_name: org-arl/jajub


text:
function diffusion(r,x,eta)

#   constructs sparse matrix Msp for "diffusion" or damping term
#
#   input
#   r        - radial grid (includes boundaries)
#   theta    - theta grid (excludes poles and equator)
#   eta      - dimensionless diffusivity

   nr = length(r);
   nx = length(x);
   dr = r[2]-r[1];
   dr2 = dr * dr;

#  initialize arrays
   ndim = nx * nr;      # number of grid points per variable
   nelements =  3 * ndim;
   row = Array{Int64}(undef,nelements);
   col = Array{Int64}(undef,nelements);
   val = Array{ComplexF64}(undef,nelements);

   #
   # fill elements of matrix (excluding boundaries)
   count = 0;
   for l = 1 : nx
       for k = 1 : nr

           ipos = ndim + (l-1) * nr +  k;

           # second derivative (central difference)
           if (k>1 && k < nr)
             count = count + 1;
             row[count] = ipos;
             col[count] = ipos;
             val[count] = -2 * eta /dr2;

             count = count + 1;
             row[count] = ipos;
             col[count] = ipos - 1;
             val[count] = 1 * eta /dr2;

             count = count + 1;
             row[count] = ipos;
             col[count] = ipos + 1;
             val[count] = 1 * eta /dr2;
           end

           if (k == 1)    # impose bc (uneven dr)

               count = count+1;
               row[count] = ipos;
               col[count] = ipos;
               val[count] = -5 * eta  / dr2;

               count = count + 1;
               row[count] = ipos;
               col[count] = ipos + 1;
               val[count] = 2 * eta /dr2;

               count = count + 1;
               row[count] = ipos;
               col[count] = ipos + 2;
               val[count] = (-1/5) * eta/ dr2;

           end

           if (k == nr)    # impose bc (uneven dr)

               count = count + 1;
               row[count] = ipos;
               col[count] = ipos;
               val[count] = -5 * eta / dr2;

               count = count + 1;
               row[count] = ipos;
               col[count] = ipos-1;
               val[count] = 2 * eta /dr2;

               count = count + 1;
               row[count] = ipos;
               col[count] = ipos - 2;
               val[count] = (-1/5) * eta / dr2;

           end

       end
   end

   dsp = sparse(row,col,val,2*ndim,2*ndim);

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
function rkqs(y::Vector{Float64},dy::Vector{Float64},t::Float64,htry::Float64,
    hdid::Float64,hnext::Float64,ys::Vector{Float64},sparams::WK3.SolverParams,
    mparams::WK3.ModelParams,k::Float64)
    # quality control parameters
    safety = 0.9;
    pgrow = -0.2;
    pshr = -0.2;
    errcon = 1.89e-4;
    errmax = 0.0;

    # allocators for step size, new final time, error estimate, and solution
    htemp = 0.0;
    tnew = 0.0;
    yerr = zeros(sparams.nvar);
    ytemp = zeros(sparams.nvar);

    # first guess for step size
    h = htry;

    while true
        # attempt an integration step
        # println("RKCK:")
        rkck(y,dy,sparams.nvar,t,h,ytemp,yerr,k,mparams);
        # check for error < tolerance (and possibly exit)
        dP = abs.(elastancefn(t*ts,mparams,k)[1].*(y[2].*Vs .- mparams.V0)./Ps .- y[1]);
        if (dP >= sparams.eps && sparams.h0 == sparams.h0nom)
            errmax = 0.0;
            for i = 1:sparams.nvar
                errmax = max(errmax,abs.(yerr[i]/ys[i]));
            end
            errmax /= sparams.eps;
            if errmax <= 1.0
                break
            end
            # resize step if error > tolerance
            htemp = safety*h*errmax^pshr;
            h = h > 0 ? max(htemp,0.1*h) : min(htemp,0.1*h);
            tnew = t+h;
            if tnew == t
                error("Step size underflow in rkqs.")
            end
        elseif dP < sparams.eps
            warn("Transvalvular pressure difference under tolerance.
                Plv - Pa = $dP at t = $t. Taking reduced-order time step.")
            break
        elseif sparams.h0 != sparams.h0nom # assimilating data, don't change stepsize
            break
        end
    end

    # guess next step size
    if errmax > errcon
        hnext = safety*h*errmax^pgrow;
    else
        hnext = 5*h;
    end

    # output solution, new time, and step size taken
    y .= ytemp;
    t += h;
    hdid = h;
    return t,hdid,hnext
end


text:
module RandomFerns

using Distributions
using StatsBase

export Ensemble, kfoldCV_fern, kfoldCV_fernEnsemble, build_fern, apply_fern, build_ensemble, apply_ensemble

if VERSION.minor >= 4
  typealias Range1{Int} Range{Int}
  _int(x) = round(Int, x)
else
  _int(x) = int(x)
end

include("utils.jl")

immutable Part
  dist::MvNormal
  mean::Any
end

immutable Fern
  parts::Array{Part}
  proj_matrix::Matrix
  cuts::Vector
end

immutable Ensemble
  ferns::Array{Fern}
end

function random_projection_matrix(method::Symbol, n_features::Int, n_components::Int)

  P = zeros(n_components, n_features)

  if method == :gaussian
    dist = Normal(0, (1/sqrt(n_components)))
    rand!(dist, P)

  elseif method == :sparse

    # todo

  end

  return P

end

function _get_index{T<:FloatingPoint}(x::Matrix{T}, cuts::Vector{T}, n::Int)
  idx = @parallel (+) for i=1:n
    c = x[:,i] .< cuts[i]
    b = c .* 2^i
  end

  idx = idx + 1
  return idx
end

### regression model ###

function apply_fern{T<:FloatingPoint}(model::Fern, features::Matrix{T})

  yhat = zeros(size(features)[1])
  w = zeros(size(features)[1])

  # project unseen data
  x = features * model.proj_matrix'

  # get indices
  idx = _get_index(x, model.cuts, size(model.proj_matrix)[1])

  # we only need to consider those parts
  uidx = unique(idx)
  for i = 1:size(uidx)[1]
    j = uidx[i]
    s = idx .== j

    if model.parts[j].mean == 0.0
      yhat[s] = 0.0
      w[s] = 0.0
    else
      yhat[s] = model.parts[j].mean
      w[s] = pdf(model.parts[j].dist, x[s,:]') ./ pdf(model.parts[j].dist, mean(model.parts[j].dist))
    end
  end

  return yhat, w
end

function apply_ensemble{T<:FloatingPoint}(model::Ensemble, features::Matrix{T})

  n = size(model.ferns)[2]

  yhat = zeros(size(features)[1])
  wsum = zeros(size(features)[1])

  for i=1:n
    yy, w = apply_fern(model.ferns[i], features)
    yhat = yhat + (yy .* w)
    wsum = wsum + w
  end

  r = yhat ./ wsum

  b = isnan(r)
  r[b] = 0

  return r

end

function build_fern{T<:FloatingPoint, U<:Real}(labels::Vector{T}, features::Matrix{U}, n_components=5)

  # sanity check
  if size(labels)[1] != size(features)[1]
    println("number of instances in features and labels matrix not equal!")
    return None
  end

  # get random projection and project features
  p = random_projection_matrix(:gaussian, size(features)[2], n_components)
  x = features * p'

  # get cutting values
  min_x = minimum(x, 1)
  max_x = maximum(x, 1)

  cuts = zeros(n_components)

  for i = 1:n_components
    d = Uniform(min_x[i], max_x[i])
    cuts[i] = rand(d)
  end

  # estimate index
  idx = _get_index(x, cuts, n_components)

  max_idx = 2^(n_components+1)
  # fill all parts with Distributions
  parts = Array(Part, 1, max_idx)
  for i = 1:max_idx
    s = idx .== i

    if (sum(s) == 0)
      x_cov = eye(n_components)
      x_mean = zeros(n_components)
      y_mean = 0
    else
      c = mean_and_cov(x[s,:])
      x_cov = c[2] .* eye(n_components)
      x_mean = vec(c[1])
      y_mean = mean(labels[s])

      # do something if covariance matrix is not positive definite
      if (det(x_cov) <= 0.0)
        x_cov = eye(n_components) # todo: this is probably a bad idea
        y_mean = 0.0
      end

    end

    dist = MvNormal(x_mean, x_cov)
    parts[i] = Part(dist, y_mean)
  end

  return Fern(parts, p, cuts)
end

function build_ensemble{T<:FloatingPoint, U<:Real}(labels::Vector{T}, features::Matrix{U}, n_components=5, n_ferns=10)

  # sanity check
  if size(labels)[1] != size(features)[1]
    println("number of instances in features and labels matrix not equal!")
    return None
  end

  ferns = Array(Fern, 1, n_ferns)

  for i=1:n_ferns
    ferns[i] = build_fern(labels, features, n_components)
  end

  return Ensemble(ferns)
end

end # module


text:
abstract type AbstractUtility end

# Euler equation deviation. Generic!
function euler_dev(
    u :: AbstractUtility,
    cV :: AbstractVector{Float64},
    betaR :: Float64
    )
    muV = marg_utility(u, cV);
    T = length(cV);
    devV = muV[1 : (T-1)] ./ muV[2 : T] .- betaR;
    return devV
end

## ----------  Log

struct UtilityLog <: AbstractUtility end

# Note the broadcasting dot. This now works for scalars and arrays.
utility(u :: UtilityLog, c) = log.(c);
marg_utility(u :: UtilityLog, c) = 1.0 ./ c;
inv_utility(u :: UtilityLog, util) = exp.(util);
inv_marg_utility(u :: UtilityLog, mu) = 1.0 ./ mu;
c_growth(u :: UtilityLog, betaR) = betaR;

## -----------  CRRA

struct UtilityCRRA <: AbstractUtility
    sigma :: Float64
end

utility(u :: UtilityCRRA, c) = (c .^ (1.0 - u.sigma)) ./ (1.0 - u.sigma) .- 1.0;
marg_utility(u :: UtilityCRRA, c) = c .^ (-u.sigma);
inv_utility(u :: UtilityCRRA, util) = ((util .+ 1.0) .* (1.0 - u.sigma)) .^ (1.0 / (1.0 - u.sigma));
inv_marg_utility(u :: UtilityCRRA, mu) = mu .^ (-1.0 / u.sigma);
c_growth(u :: UtilityCRRA, betaR) = betaR .^ (1.0 / u.sigma);

# -------------


text:
module PaddedBlocks

using DocStringExtensions
using JSON3

include("blocks.jl")

export Blocks, slices

end


text:
# Parser
function parse_array!(tokens::Vector{Any})
    jsonArray = Vector{Any}()

    if tokens[1] == JSON_RIGHTBRACKET
        popfirst!(tokens)
        return jsonArray
    end # if

    @inbounds while true
        rest = @views parser!(tokens)
        push!(jsonArray, rest)

        if tokens[1] == JSON_RIGHTBRACKET
            popfirst!(tokens)
            return jsonArray
        elseif tokens[1] != JSON_COMMA
            error("No comma found in array")
        else
            popfirst!(tokens)
        end # if
    end # while

end # parse_array

function parse_object!(tokens::Vector{Any})
    jsonObject = Dict{String, Any}()

    if tokens[1] == JSON_RIGHTBRACE
        popfirst!(tokens)
        return jsonObject
    end # if

    @inbounds while true
        jsonKey = tokens[1]
        if typeof(jsonKey) == String
            popfirst!(tokens)
        else
            error("Expected key string, but got ", jsonKey)
        end # if

        if tokens[1] != JSON_COLON
            error("Expected colon but got ", tokens[1])
        end # if

        popfirst!(tokens)

        jsonObject[jsonKey] = @views parser!(tokens)

        if tokens[1] == JSON_RIGHTBRACE
            popfirst!(tokens)
            return jsonObject
        elseif tokens[1] != JSON_COMMA
            error("Expected comma in object but got", tokens[1])
        end # if

        try
            popfirst!(tokens)
        catch
        end
    end # while
end # parse_object

function parser!(tokens::Vector{Any})
    t = popfirst!(tokens)
    if t == JSON_LEFTBRACKET
        return parse_array!(tokens)
    elseif t == JSON_LEFTBRACE
        return parse_object!(tokens)
    else
        return t
    end # if

end # parser

max_stars_repo_name: Iaggelis/jjson

max_issues_repo_name: Iaggelis/jjson

max_forks_repo_name: Iaggelis/jjson


text:
using Test

function check(n)
    d = reverse(digits(n))
    double = false

    for i in 1:length(d)-1
        if d[i] > d[i+1]
            return false
        end
        if  d[i] == d[i+1]  # we have a double
            if i == 1
                if d[i+1] ≠ d[i+2]
                    double = true
                end
            elseif i == length(d) - 1
                if d[i-1] ≠ d[i]
                    double = true
                end
            else
                if d[i-1] ≠ d[i] && d[i+1] ≠ d[i+2]
                    double = true
                end
            end
        end
    end
    return double
end

function run(s)
    # println(s)
    n1, n2 = parse.(Int, (split(s, "-")))
    total = 0
    return sum(check.(n1:n2))
end

@test check(112233) == true
@test check(223334) == true
@test check(123444) == false
@test check(111122) == true

##############################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
# simple scan into a list of free buffers

## tunables

# how much larger a buf can be to fullfil an allocation request.
# lower values improve efficiency, but increase pressure on the underlying GC.
function max_oversize(sz)
    if sz <= 2^20       # 1 MiB
        # small buffers are fine no matter
        return typemax(Int)
    elseif sz <= 2^20   # 32 MiB
        return 2^20
    else
        return 2^22
    end
end

## pooling

Base.@kwdef struct SimplePool <: AbstractPool
    stream_ordered::Bool

    lock::ReentrantLock = ReentrantLock()
    cache::Set{Block} = Set{Block}()

    freed_lock::NonReentrantLock = NonReentrantLock()
    freed::Vector{Block} = Vector{Block}()
end

function scan(pool::SimplePool, sz)
    @lock pool.lock for block in pool.cache
        if sz <= sizeof(block) <= max_oversize(sz)
            delete!(pool.cache, block)
            return block
        end
    end
    return
end

function repopulate(pool::SimplePool)
    blocks = @lock pool.freed_lock begin
        isempty(pool.freed) && return
        blocks = Set(pool.freed)
        empty!(pool.freed)
        blocks
    end

    @lock pool.lock begin
        for block in blocks
            @assert !in(block, pool.cache)
            push!(pool.cache, block)
        end
    end

    return
end

function reclaim(pool::SimplePool, sz::Int=typemax(Int))
    repopulate(pool)

    @lock pool.lock begin
        freed_bytes = 0
        while freed_bytes < sz && !isempty(pool.cache)
            block = pop!(pool.cache)
            freed_bytes += sizeof(block)
            actual_free(block; pool.stream_ordered)
        end
        return freed_bytes
    end
end

function alloc(pool::SimplePool, sz; stream::CuStream)
    block = nothing
    for phase in 1:3
        if phase == 2
            GC.gc(false)
        elseif phase == 3
            GC.gc(true)
        end

        repopulate(pool)

        block = scan(pool, sz)
        block === nothing || break

        block = actual_alloc(sz; pool.stream_ordered)
        block === nothing || break

        reclaim(pool, sz)
        block = actual_alloc(sz, phase==3; pool.stream_ordered)

        block === nothing || break
    end

    return block
end

function free(pool::SimplePool, block; stream::CuStream)
    # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
    # and to simplify locking (preventing concurrent access during GC interventions)
    @spinlock pool.freed_lock begin
        push!(pool.freed, block)
    end
end

function cached_memory(pool::SimplePool)
    sz = @lock pool.freed_lock mapreduce(sizeof, +, pool.freed; init=0)
    sz += @lock pool.lock mapreduce(sizeof, +, pool.cache; init=0)
    return sz
end


text:
#Code written by Skipper af Rosenborg and Esben Bager
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap

include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

## LOAD DATA
#Skipper's path
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"

#HPC path
#inputArg = parse(Int64, ARGS[1]) #Should range from 0 to 106
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"
industry = "Utils"
mainData = loadIndexDataLOGReturn(industry, path)

## PIPELINE
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results"
#path = "/zhome/9f/d/88706/SpecialeCode/Results"
#### MUST CHANGE ####
fileName = path*"/IndexData/AALRTest/"

#Reset HPC path
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML"
#cd(path)

dateAndReseccion = Array(mainData[:,end-1:end])
mainDataArr = Array(mainData[:,1:end-2])

colNames = names(mainData)

expandedColNamesToString(expandColNamesTimeToString(colNames),1,1)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]
mainYarr = mainDataArr[:, nCols:nCols]

# Transform with time elements
mainXarr = expandWithTime3612(mainXarr)

# Transform #
mainXarr = expandWithTransformations(mainXarr)
mainXarr = expandWithMAandMomentum(mainXarr, mainYarr, (nCols-1))

# Standardize
standX = zScoreByColumn(mainXarr)
standY = mainDataArr[:, nCols:nCols]
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]

## OUTPUT OF DATA TO OTHER PROGRAMS
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/PythonCode/Data"
cd(path)
writedlm("allData"*industry*".csv", allData,",")
allData[1,end]
standX[1,end]

trainingSize = 12
predictions = 1
testRuns = nRows-trainingSize-predictions

Xtrain = allData[:, 1:bCols]
trainingData = Xtrain[:,1:nCols-1]

## STAGE1 OF AALR
function solveStage1(trainingData)
	### STAGE 1 ###
	println("STAGE 1 INITIATED")
	#Define solve
	m = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads = 1))

	#Add binary variables variables
	@variable(m, 0 <= z[1:nCols-1] <= 1, Bin )

	#Calculate highly correlated matix
	HC = cor(trainingData)

	#Define objective function
	@objective(m, Max, sum(z[i] for i=1:length(z)))

	#Define max correlation and constraints
	rho = 0.8
	for i=1:length(z)
		for j=1:length(z)
			if i != j
				if HC[i,j] >= rho
					@constraint(m,z[i]+z[j] <= 1)
				end
			end
		end
	end

	#Get solution status
	status = solve(m);

	#Get objective value
	println("Solved stage 1 with kMax = ", getobjectivevalue(m))
	return kmax = getobjectivevalue(m)
end

kmax = solveStage1(trainingData)
#Get solution value
#zSolved = getvalue(z)
println("STAGE 1 DONE")

## STAGE 2
println("STAGE 2 INITIATED")
bigM = 5
tau = 2
amountOfGammas = 5

bSolved = []
bestBeta = []
prevSolutions = zeros(Int64(kmax*amountOfGammas), bCols)

## DEFINE STAGE2 AND THE SOLVING OF IT.
function buildAndSolveStage2(standX, standY, curKmax, gamma, warmstartBool, warmStartSol, cuts, fixedZ, HC)
	HCPairCounter = 0

	#Define parameters and model
	stage2Model = JuMP.Model(solver = GurobiSolver(TimeLimit = 40, OutputFlag = 0, Threads = 1, PreMIQCPForm=0, MIPFocus=1, ImproveStartTime=20));

	#Define variables
	@variable(stage2Model, b[1:bCols]) #Beta values
	@variable(stage2Model, 0 <= z[1:bCols] <= 1, Bin ) #binary z variable
	@variable(stage2Model, v[1:bCols]) # auxiliary variables for abs
	@variable(stage2Model, T) #First objective term
	@variable(stage2Model, G) #Second objective term

	#Initialise warmstart
	if warmstartBool == true
		for j in find(warmStartSol)
			setvalue(b[j], warmStartSol[j])
			setvalue(v[j], abs(warmStartSol[j]))
			setvalue(z[j], 1)
		end
	end

	#Define objective function (5a)
	@objective(stage2Model, Min, T+G)
	@constraint(stage2Model, soc, norm( [1-T;sqrt(2)*(standX*b-standY)] ) <= 1+T)
	@constraint(stage2Model, gammaConstr, gamma*ones(bCols)'*v <= G) #4bCols

	#Define constraints (5c)
	@constraint(stage2Model, conBigMN, -1*b .<= bigM*z) #from 0 to bCols -1
	@constraint(stage2Model, conBigMP,  1*b .<= bigM*z) #from bCols to 2bCols-1

	#Second objective term
	@constraint(stage2Model, 1*b .<= v) #from 2bCols to 3bCols-1
	@constraint(stage2Model, -1*b .<= v) #from 3bCols to 4bCols-1

	#Define kmax constraint (5d)
	@constraint(stage2Model, kMaxConstr, sum(z[j] for j=1:bCols) <= curKmax) #4bCols+1

	#Constraint 5f - can only select one of a pair of highly correlated features
	rho = 0.8
	for k=1:bCols
		for j=1:bCols
			if k != j
				if HC[k,j] >= rho
					HCPairCounter += 1
					@constraint(stage2Model,z[k]+z[j] <= 1) #from 4bCols+1 to (4bCols+1+HCPairCounter)
				end
			end
		end
	end

	#Constraint (5g) - only one transformation allowed (x, x^2, log(x) or sqrt(x))
	for j=1:(nCols-1)
		@constraint(stage2Model, z[j]+z[j+(nCols-1)]+z[j+2*(nCols-1)]+z[j+3*(nCols-1)] <= 1) #from (4bCols+1+HCPairCounter) to (4bCols+1+HCPairCounter+nCols)
	end

	#Implement cuts
	if !isempty(cuts)
		#nEmptyCuts = 19
		#zOnes = ones(bCols)
		for c = 1:size(cuts)[1]
			@constraint(stage2Model, sum(z[i] for i=find(cuts[c,1:end])) <= countnz(cuts[c,1:end])-1)
			#@constraint(stage2Model, sum(zOnes[i]*z[i] for i=1:bCols) <= bCols) #(4bCols+1+HCPairCounter+nCols) to (4bCols+1+HCPairCounter+nCols+18)
		end
	end

	# Fixes "z" in case of bootstrapping
	if !isempty(fixedZ)
		@constraint(stage2Model, sum(z[i] for i=find(fixedZ)) == countnz(fixedZ))
	end

	# Solve model
	status = solve(stage2Model)

	if status == :InfeasibleOrUnbounded
		bSolved = zeros(bCols)
        obj = 9999
        objBound = 9999
	else
		bSolved = getvalue(b)
		if any(isnan,bSolved)
			bSolved = zeros(bCols)
		end
		obj = getobjectivevalue(stage2Model)
		objBound = getobjectivebound(stage2Model)
	end
	return bSolved, obj, objBound
end

## DEFINING STAGE 3 AND THE BOOTSTRAPPING
function stageThree(best3Beta, X, Y, allCuts)
	bCols = size(X)[2]
	nRows = size(X)[1]
	cuts = Matrix(0, bCols+1)
	rowsPerSample = nRows #All of rows in training data to generate beta estimates, but selected with replacement
	totalSamples = 50 #50 different times we will get a beta estimate
	nBoot = 1000

	#For loop start
	for i = 1:size(best3Beta)[1]
		if signifBoolean[i] == 1 #if the previous solution was already completely significant, skip the bootstrapping
			continue
		end
		bestBeta = best3Beta[i,4:bCols+3]
		bestK = best3Beta[i,2]
		bestGamma = best3Beta[i,3]

		xColumns = []
		bSample = Matrix(totalSamples, bCols)

		bZeros = zeros(bCols)
		for j = 1:bCols
			if !isequal(bestBeta[j],0)
				push!(xColumns, j)
			end
		end

		selectedX = X[:,xColumns]
		condNumber = cond(selectedX)
		if condNumber >= 15
			bZeros[xColumns] = 1
			subsetSize = size(xColumns)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on Condition number = $condNumber has been created from Beta$i")
		end

		# test significance
		bestZ = zeros(bestBeta)
		for l=1:size(bestBeta)[1]
			if bestBeta[l] != 0
				bestZ[l] = 1
			end
		end

		bZeros = zeros(bCols)
		bSample = createBetaDistribution(bSample, X, Y, bestK, totalSamples, rowsPerSample,  bestGamma, allCuts, bestZ) #standX, standY, k, sampleSize, rowsPerSample

		bSample = zeros(50,bCols)
		confArray99 = createConfidenceIntervalArray(bSample, nBoot, 0.99)
		confArray95 = createConfidenceIntervalArray(bSample, nBoot, 0.95)
		confArray90 = createConfidenceIntervalArray(bSample, nBoot, 0.90)

		bestBeta = zeros(bCols)

		significanceResult = testSignificance(confArray99, confArray95, confArray90, bestBeta)
		significanceResultNONSI = [] # = significanceResult[xColumns]
		subsetSize = size(xColumns)[1]
		for n = 1:size(significanceResult)[1]
			for s = 1:subsetSize
				if significanceResult[n] == 0 && xColumns[s] == n
					push!(significanceResultNONSI,xColumns[s])
					println("Parameter $n is selected, but NOT significant")
				elseif significanceResult[n] > 0 && xColumns[s] == n
					println("Parameter $n is significant with ", significanceResult[n])
				end
			end
		end

		if !isempty(significanceResultNONSI)
			bZeros[significanceResultNONSI] = 1
			subsetSize = size(significanceResultNONSI)[1]
			newCut = [bZeros' subsetSize]
			cuts = [cuts; newCut]
			println("A cut based on parameters being non-significant in Beta$i has been created")
		end
		if isempty(significanceResultNONSI)
			signifBoolean[i] = 1
		end
	end
	return cuts
end

## FUNCTION OT GET RMSE VALUE
function checkRMSE(standX, standY, bSolved, best3Beta, k, gamma)
	RMSE = getRMSE(standX,standY,bSolved)

	if RMSE < best3Beta[1,1] #largest than 3rd largest RMSE
		if RMSE < best3Beta[2,1] #largest than 2nd largest RMSE
			if RMSE < best3Beta[3,1] #largest than largest RMSE
				best3Beta[1,:] = best3Beta[2,:]
				best3Beta[2,:] = best3Beta[3,:]
				best3Beta[3,1] = RMSE
				best3Beta[3,2] = k
				best3Beta[3,3] = gamma #potentially store actual gamma value
				best3Beta[3,4:bCols+3] = bSolved
			else
				best3Beta[1,:] = best3Beta[2,:]
				best3Beta[2,1] = RMSE
				best3Beta[2,2] = k
				best3Beta[2,3] = gamma #potentially store actual gamma value
				best3Beta[2,4:bCols+3] = bSolved
			end
		else
			best3Beta[1,1] = RMSE
			best3Beta[1,2] = k
			best3Beta[1,3] = gamma #potentially store actual gamma value
			best3Beta[1,4:bCols+3] = bSolved
		end
	end

	return best3Beta
end

stage2Model = 0
best3Beta = zeros(3,bCols+3)
best3Beta[:,1] = 1e3
signifBoolean = zeros(3)

## MAIN LOOP TO RUN OVER DIFFERENT TIME PERIODS
@time(
for r = 1+inputArg*10:10+inputArg*10#(nRows-trainingSize-predictions)

	#SPLIT DATA
	trainingData = Xtrain[:,1:nCols-1]
	standX = allData[r:(trainingSize+r-1), 1:bCols]
	standY = allData[r:(trainingSize+r-1), bCols+1]
	standXTest  = allData[(trainingSize+r):(trainingSize+r+predictions-1), 1:bCols]
	standYTest  = allData[(trainingSize+r):(trainingSize+r+predictions-1), bCols+1]
	dateAndReseccionTest = dateAndReseccion[(trainingSize+r):(trainingSize+r+predictions-1), :]
	nRows = size(standX)[1]

	SSTO = sum((standY[i]-mean(standY))^2 for i=1:length(standY))
	amountOfGammas = 5

	gammaArray = zeros(amountOfGammas)
	maxVal = findmax(abs.(standX'*standY))[1]
	gammaArray = Array(linspace(maxVal,0,amountOfGammas))

	#Spaced between 0 and half the SSTO since this would then get SSTO*absSumOfBeta which would force everything to 0
	#gammaArray = log10.(logspace(SSTO/2, 0, amountOfGammas))

	#INITIALISE STORING ARRAYS
	bSample = []
	allCuts = []
	signifBoolean = zeros(3)
	best3Beta = zeros(3,bCols+3)
	best3Beta[:,1] = 1e3
	runCount = 1
	curRMSE = 1e3
	HC = cor(standX)

	#BUILD AND SOLVE MODEL
	bSolved = []
	warmstart = false
	for i = 1:kmax
		for g = 1:amountOfGammas
			println("\nSolving for k = $i and gamma index = $g\n")
			bSolved, obj, objBound = buildAndSolveStage2(standX, standY, i, gammaArray[g], warmstart, bSolved, [], [], HC)

			#println("Cur index is: ",Int64(g+(i-1)*amountOfGammas))
			prevSolutions[Int64(g+(i-1)*amountOfGammas),:] = bSolved

			best3Beta = checkRMSE(standX, standY, bSolved, best3Beta, i, g)

			warmstart = true
		 end
	end

	#CHECK FOR SIGNIFICANCE (THIS INCLUDES SOLVING THE MODEL WITH FIXED z[] VARIABLES)
	cuts = stageThree(best3Beta, standX, standY, allCuts)

	#Check for cuts
	while !isempty(cuts)
		if runCount > 4
			break
		end
		best3Beta[:,1] = 1e3
		allCuts = vcat(allCuts,cuts[:,1:end-1])

		# Based on the cuts, ONLY resolve the problems that the cut affects.
		# Find the minimum number of variables in over all cuts
		minNumberOfVar = kmax
		for j=1:size(allCuts)[1]
			if countnz(allCuts[1,:]) < kmax
				minNumberOfVar = countnz(allCuts[1,:])
			end
		end

		#Rebuild and solve model with cuts
		for i = 1:kmax
			for g = 1:amountOfGammas
				println("\nSolving for k = $i and gamma index = $g\n")
				bSolved = prevSolutions[Int64(g+(i-1)*amountOfGammas),:]
				if i >= minNumberOfVar
					#println("\nSolving for k = $i and gamma index = $g\n")
					bSolved, obj, objBound = buildAndSolveStage2(standX, standY, i, gammaArray[g], warmstart, prevSolutions[Int64(g+(i-1)*amountOfGammas),:], allCuts, [], HC)
				end

				#prevSolutions[Int64(i+(i-1)*amountOfGammas),:] = bSolved

				best3Beta = checkRMSE(standX, standY, bSolved, best3Beta, i, g)

				warmstart = true
			 end
		end

		#Check significance
		cuts = stageThree(best3Beta, standX, standY, allCuts)

		runCount += 1
	end

	best3Beta = cat(2,signifBoolean,best3Beta)

	prediction = Array{Float64}(3,1)
	prediction[1,1] = (standXTest*best3Beta[1,5:end])[1]
	prediction[2,1] = (standXTest*best3Beta[2,5:end])[1]
	prediction[3,1] = (standXTest*best3Beta[3,5:end])[1]

	best3Beta = cat(2,prediction,best3Beta)

	real = Array{Float64}(3,1)
	real[1,1] = standYTest[1]
	real[2,1] = standYTest[1]
	real[3,1] = standYTest[1]

	best3Beta = cat(2,real,best3Beta)

	dateAndReseccionOutput = Array{Int64}(3,2)
	dateAndReseccionOutput[1,:] = dateAndReseccionTest
	dateAndReseccionOutput[2,:] = dateAndReseccionTest
	dateAndReseccionOutput[3,:] = dateAndReseccionTest

	best3Beta = cat(2,dateAndReseccionOutput,best3Beta)

	### EXTRACT PREDICTION AND ISRS AND DATE###
	f = open(fileName*string(r)*"AALRBestK.csv", "w")
	write(f, "Date,Reseccion,Real,Prediciton,Significant,RMSE,k,gamma,"*expandedColNamesToString(colNames)*"\n")
	writecsv(f,best3Beta)
	close(f)
end
)
println("Complete")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/AALR_Time.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/AALR_Time.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
function include_all(folder::String)
  for file in readdir(folder)
    if endswith(file, ".jl")
        include(joinpath(folder, file))
    end
  end
end

function Base.split(condition::Function, associative::Associative)
  A = similar(associative)
  B = similar(associative)
  for (key, value) in associative
    if condition(key, value)
      A[key] = value
    else
      B[key] = value
    end
  end
  A, B
end


text:
struct ExceptionRethrow end

default_action(::ExceptionRethrow, mdp, s, ex) = rethrow(ex)
default_action(f::Function, mdp, s, ex) = f(s, ex)
default_action(p::POMDPs.Policy, mdp, s, ex) = action(p, s)
default_action(sol::POMDPs.Solver, mdp, s, ex) = action(solve(sol, mdp), s)
default_action(a, mdp, s, ex) = a

"""
    ReportWhenUsed(a)

When the planner fails, returns action `a`, but also prints the exception.
"""
struct ReportWhenUsed{T}
    a::T
end

function default_action(r::ReportWhenUsed, mdp, s, ex)
    showerror(stderr, ex)
    a = default_action(r.a, mdp, s, ex)
    @warn("Using default action $a")
    return a
end


text:
# Use version 1.5.3 and run in command line
#
#   julia --project=@. run_Dash.jl
#
# this will enable hot reloading, meaning the server refresh the page when
# learn_Dash.jl is updated.
#
# Run in interactive mode is fine by
#
#   include("learn_Dash.jl")
#
# but needs to shutdown and restart the server every time the file updated.

#
# Use this to run functions in learn_Dash.jl
# The separation is necessary due to hot reloading from Dash.  It automatically
# reload and this may not be desired sometimes.  The use of run_Dash.jl creates
# a buffer to get more control of when to reload.
#
# Apparently, hot redloading check included file so the idea of buffering is out.

include("learn_Dash.jl")

# https://dash-julia.plotly.com/getting-started
#hello_dash()
#dash_table()
dash_plotly() # RDatasets does not work
#dash_markdown()
#dash_core()
#dash_basic_callback()
#dash_layout_figure_slider()
#Ans = dash_app_multiple_inputs()
#Ans = dash_app_multiple_outputs()
#Ans = dash_app_chained_callbacks()
#Ans = dash_app_state()
#Ans = dash_intro_visual()
#Ans = dash_update_graphs_on_hover()
#Ans = dash_generic_crossfilter_recipe()