text:
# include("MeanFieldTheory.jl")
# include("Configuration.jl")
# include("VMCMatrix.jl")
# include("Random.jl")
# include("Move.jl")
# include("delegate.jl")

# having the type parameter here is actually pretty stupid, since some wavefunctions can have multiple types (e.g., CFL)
abstract Wavefunction

# spinless fermions or hard-core bosons:
abstract HardcoreWavefunction <: Wavefunction

################################################################
#################### Free spinless fermions ####################
################################################################

type FreeFermionWavefunction{T} <: HardcoreWavefunction
    orbitals::FreeFermionOrbitals{T}
    config::LatticeParticlesConfiguration
    cmat::CeperleyMatrix{T}
    neighbor_table::Vector{Vector{Int}}
    proposed_move::ParticlesMove{LatticeParticleMove}

    function FreeFermionWavefunction(orbs::FreeFermionOrbitals{T}, nbs::Vector{Vector{Int}})
        Nsites = get_Nsites(orbs)
        Nparticles = get_Nfilled(orbs)
        if ~isempty(nbs)
            @assert length(nbs) == Nsites
        end

        positions = random_combination(Nsites, Nparticles)
        config = LatticeParticlesConfiguration(Vector{Int}[positions], Nsites)

        cmat = CeperleyMatrix(fill_detmat(orbs, config))

        return new(orbs, config, cmat, nbs)
    end
end

function FreeFermionWavefunction{T}(orbs::FreeFermionOrbitals{T}, nbs::Vector{Vector{Int}})
    return FreeFermionWavefunction{eltype(orbs.filled)}(orbs, nbs)
end

function FreeFermionWavefunction{T}(orbs::FreeFermionOrbitals{T})
    return FreeFermionWavefunction{eltype(orbs.filled)}(orbs, Vector{Int}[])
end

get_Nspecies(wf::HardcoreWavefunction) = 1 # I think this is what we want
get_Nparticles(wf::HardcoreWavefunction) = get_Nparticles(wf.config)[1]
get_Nsites(wf::HardcoreWavefunction) = get_Nsites(wf.config)

function check_for_numerical_error!(wf::FreeFermionWavefunction)
    check_and_reset_inverse_and_det!(wf.cmat, "cmat")
end

function random_empty_site_move(wf::HardcoreWavefunction)
    particle_to_move = Particle(rand(1:get_Nparticles(wf)), 1)

    Nempty = get_Nsites(wf) - get_Nparticles(wf)
    proposed_position = setdiff(collect(1:get_Nsites(wf)), wf.config.particle_positions[1])[rand(1:Nempty)]

    return ParticlesMove([LatticeParticleMove(particle_to_move, proposed_position)])
end

function random_neighbor_move(wf::HardcoreWavefunction)
    particle_to_move = Particle(rand(1:get_Nparticles(wf)), 1)

    current_site = wf.config.particle_positions[1][particle_to_move.index]
    proposed_position = wf.neighbor_table[current_site][rand(1:end)]

    return ParticlesMove([LatticeParticleMove(particle_to_move, proposed_position)])
end

function assess_move!(wf::FreeFermionWavefunction, move::ParticlesMove{LatticeParticleMove})
    # @assert is_valid_move(move, wf)

    @assert length(move) == 1 # for now
    if is_occupied(wf.config, move[1].destination, 1)
        # finish_move! could technically still be called after this (i.e., we could try to accept the move),
        # but in practice it never will since probrat = 0.
        # This logic could afford to be a little tighter however..
        wf.cmat.current_state = :COLUMNS_UPDATE_PENDING
        probrat = 0.
    else
        detrat_from_columns_update!(wf.cmat, move[1].particle.index,
            columns_for_detmat(wf.orbitals, move[1].destination))
        probrat = abs2(get_detrat(wf.cmat))
    end

    # nAtest = 4
    # subsystemA = collect(1:4)
    # if get_Nparticles_in_subsystem(update_configuration(wf.config, move), subsystemA)[1] == nAtest
    #     printlnf("Proposing a move with $nAtest particles in subsystem A:  probrat = $probrat")
    # end

    wf.proposed_move = move

    return probrat
end

# FIXME: this actually returns a CeperleyMatrix object..
# Does this then have to be garbage collected etc., leading to slow-down?
# Issue probably crops up everywhere..
function finish_move!(wf::FreeFermionWavefunction)
    update_configuration!(wf.config, wf.proposed_move)
    finish_columns_update!(wf.cmat)
end

function cancel_move!(wf::FreeFermionWavefunction)
    cancel_columns_update!(wf.cmat)
end

##########################################################
#################### SBM wavefunction ####################
##########################################################

type SBMWavefunction{T} <: Wavefunction
    orbitals_up::FreeFermionOrbitals{T}
    orbitals_down::FreeFermionOrbitals{T}

    config::LatticeParticlesConfiguration

    cmat_up::CeperleyMatrix{T}
    cmat_down::CeperleyMatrix{T}

    neighbor_table::Vector{Vector{Int}} # this could be a Vector of neighbor tables for different neighbors

    proposed_move::ParticlesMove{LatticeParticleMove} # empty upon initialization

    function SBMWavefunction(orbs_up::FreeFermionOrbitals{T}, orbs_down::FreeFermionOrbitals{T}, nbs::Vector{Vector{Int}})
        @assert get_Nsites(orbs_up) == get_Nsites(orbs_down)
        Nsites = get_Nsites(orbs_up)
        @assert get_Nfilled(orbs_up) + get_Nfilled(orbs_down) == Nsites # half filling
        if !isempty(nbs)
            @assert length(nbs) == Nsites
        end

        # get random spin configuration
        positions_up = random_combination(get_Nsites(orbs_up), get_Nfilled(orbs_up))
        config = GutzSpinhalfConfiguration(positions_up, get_Nsites(orbs_up))
        @assert is_valid_Gutz_spinhalf_configuration(config) # being anal

        # compute Ceperley matrices
        detmat_up = fill_detmat(orbs_up, config, 1)
        detmat_down = fill_detmat(orbs_down, config, 2)
        cmat_up = CeperleyMatrix(detmat_up)
        cmat_down = CeperleyMatrix(detmat_down)

        return new(orbs_up, orbs_down, config, cmat_up, cmat_down, nbs)
    end
end

# This is annoying, but now I think I understand the issue:
# After defining a parameterized type, no constructors are available without type specification (in the call),
# e.g., mytype(...) is not an available method, so we need to define it as such with the type specified on the RHS.
# Also, presence of the type parameter T in this "outer constructor" ensures commonality amongst
# the orbital types of the different partons.
function SBMWavefunction{T}(orbs_up::FreeFermionOrbitals{T}, orbs_down::FreeFermionOrbitals{T}, nbs::Vector{Vector{Int}})
    return SBMWavefunction{eltype(orbs_up.filled)}(orbs_up, orbs_down, nbs)
end

# if not included in argument list, set neighbor_table to an empty array
function SBMWavefunction{T}(orbs_up::FreeFermionOrbitals{T}, orbs_down::FreeFermionOrbitals{T})
    return SBMWavefunction(orbs_up, orbs_down, Vector{Int}[])
end

@delegate SBMWavefunction.config [ is_valid_Gutz_spinhalf_configuration ]
export is_valid_Gutz_spinhalf_configuration

get_Nspecies(wf::SBMWavefunction) = 2
get_Nparticles(wf::SBMWavefunction) = get_Nparticles(wf.config)
get_Nup(wf::SBMWavefunction) = get_Nparticles(wf)[1]
get_Ndown(wf::SBMWavefunction) = get_Nparticles(wf)[2]
get_Nsites(wf::SBMWavefunction) = get_Nsites(wf.config)

function check_for_numerical_error!(wf::SBMWavefunction)
    check_and_reset_inverse_and_det!(wf.cmat_up, "cmat_up")
    check_and_reset_inverse_and_det!(wf.cmat_down, "cmat_down")
end

function random_exchange_move(wf::SBMWavefunction)
    @assert is_valid_Gutz_spinhalf_configuration(wf)

    particle_to_move_up = Particle(rand(1:get_Nup(wf)), 1)
    particle_to_move_down = Particle(rand(1:get_Ndown(wf)), 2)

    proposed_position_up = wf.config.particle_positions[2][particle_to_move_down.index] # current position of down particle
    proposed_position_down = wf.config.particle_positions[1][particle_to_move_up.index] # current position of up particle

    mv_up = LatticeParticleMove(particle_to_move_up, proposed_position_up)
    mv_down = LatticeParticleMove(particle_to_move_down, proposed_position_down)

    return ParticlesMove([mv_up, mv_down])
end

function neighbor_exchange_move(wf::SBMWavefunction)
    @assert is_valid_Gutz_spinhalf_configuration(wf)

    random_site = rand(1:get_Nsites(wf))
    random_neighbor = wf.neighbor_table[random_site][rand(1:end)]

    # printlnf("random_site = " * string(random_site))
    # printlnf("random_neighbor = " * string(random_neighbor))

    if is_occupied(wf.config, random_site, 1) && is_occupied(wf.config, random_neighbor, 1) ||
        is_occupied(wf.config, random_site, 2) && is_occupied(wf.config, random_neighbor, 2)
        return ParticlesMove{LatticeParticleMove}() # spins are parallel, so move is trivial
    else
        if is_occupied(wf.config, random_site, 1) # up particle at random_site
            @assert is_occupied(wf.config, random_neighbor, 2) # down particle at random_neighbor
            particle_to_move_up = Particle(wf.config.lattice_status[1][random_site], 1)
            proposed_position_up = random_neighbor
            particle_to_move_down = Particle(wf.config.lattice_status[2][random_neighbor], 2)
            proposed_position_down = random_site
        elseif is_occupied(wf.config, random_site, 2) # down particle at random_site
            @assert is_occupied(wf.config, random_neighbor, 1) # up particle at random_neighbor
            particle_to_move_up = Particle(wf.config.lattice_status[1][random_neighbor], 1)
            proposed_position_up = random_site
            particle_to_move_down = Particle(wf.config.lattice_status[2][random_site], 2)
            proposed_position_down = random_neighbor
        else
            @assert false
        end
        mv_up = LatticeParticleMove(particle_to_move_up, proposed_position_up)
        mv_down = LatticeParticleMove(particle_to_move_down, proposed_position_down)
        return ParticlesMove([mv_up, mv_down])
    end
end

# similar to Jim's perform_move functions
function assess_move!(wf::SBMWavefunction, move::ParticlesMove{LatticeParticleMove})
    # @assert is_valid_move(move, wf)

    detrat_from_columns_update!(wf.cmat_up, move[1].particle.index,
        columns_for_detmat(wf.orbitals_up, move[1].destination))
    detrat_from_columns_update!(wf.cmat_down, move[2].particle.index,
        columns_for_detmat(wf.orbitals_down, move[2].destination))

    wf.proposed_move = move

    probrat_up = abs2(get_detrat(wf.cmat_up))
    probrat_down = abs2(get_detrat(wf.cmat_down))
    return probrat_up * probrat_down
end

function finish_move!(wf::SBMWavefunction)
    update_configuration!(wf.config, wf.proposed_move)
    finish_columns_update!(wf.cmat_up)
    finish_columns_update!(wf.cmat_down)
end

function cancel_move!(wf::SBMWavefunction)
    cancel_columns_update!(wf.cmat_up)
    cancel_columns_update!(wf.cmat_down)
end

#########################################################################
#################### Boson CFL = det(nu=1) * det(FS) ####################
#########################################################################

type BosonCFLWavefunction{Td,Tf} <: HardcoreWavefunction
    orbitals_d::FreeFermionOrbitals{Td}
    orbitals_f::FreeFermionOrbitals{Tf}

    config::LatticeParticlesConfiguration

    cmat_d::CeperleyMatrix{Td}
    cmat_f::CeperleyMatrix{Tf}

    neighbor_table::Vector{Vector{Int}} # this could be a Vector of neighbor tables for different neighbors

    proposed_move::ParticlesMove{LatticeParticleMove} # empty upon initialization

    function BosonCFLWavefunction(orbs_d::FreeFermionOrbitals{Td}, orbs_f::FreeFermionOrbitals{Tf},
                                  nbs::Vector{Vector{Int}}, Ly_for_init::Int)
        @assert Ly_for_init >= 0
        @assert get_Nsites(orbs_d) == get_Nsites(orbs_f)
        @assert get_Nfilled(orbs_d) == get_Nfilled(orbs_f)
        Nsites = get_Nsites(orbs_d)
        Nparticles = get_Nfilled(orbs_d)
        @assert Nparticles <= Nsites
        if !isempty(nbs)
            @assert length(nbs) == Nsites
        end

        # get initial configuration of bosons
        if Ly_for_init == 0
            positions = random_combination(Nsites, Nparticles)
            config = LatticeParticlesConfiguration(Vector{Int}[positions], Nsites)
        else
            # put one particle per rung sequentially until we run out of particles
        end

        # compute Ceperley matrices
        cmat_d = CeperleyMatrix(fill_detmat(orbs_d, config))
        cmat_f = CeperleyMatrix(fill_detmat(orbs_f, config))

        return new(orbs_d, orbs_f, config, cmat_d, cmat_f, nbs)
    end
end

function BosonCFLWavefunction{Td,Tf}(orbs_d::FreeFermionOrbitals{Td}, orbs_f::FreeFermionOrbitals{Tf},
                                     nbs::Vector{Vector{Int}}, Ly_for_init::Int)
    return BosonCFLWavefunction{eltype(orbs_d.filled),eltype(orbs_f.filled)}(orbs_d, orbs_f, nbs, Ly_for_init)
end

function BosonCFLWavefunction{Td,Tf}(orbs_d::FreeFermionOrbitals{Td}, orbs_f::FreeFermionOrbitals{Tf}, nbs::Vector{Vector{Int}})
    return BosonCFLWavefunction{eltype(orbs_d.filled),eltype(orbs_f.filled)}(orbs_d, orbs_f, nbs, 0)
end

function BosonCFLWavefunction{Td,Tf}(orbs_d::FreeFermionOrbitals{Td}, orbs_f::FreeFermionOrbitals{Tf}, Ly_for_init::Int)
    return BosonCFLWavefunction{eltype(orbs_d.filled),eltype(orbs_f.filled)}(orbs_d, orbs_f, Vector{Int}[], Ly_for_init)
end

function BosonCFLWavefunction{Td,Tf}(orbs_d::FreeFermionOrbitals{Td}, orbs_f::FreeFermionOrbitals{Tf})
    return BosonCFLWavefunction{eltype(orbs_d.filled),eltype(orbs_f.filled)}(orbs_d, orbs_f, Vector{Int}[], 0)
end

function check_for_numerical_error!(wf::BosonCFLWavefunction)
    check_and_reset_inverse_and_det!(wf.cmat_d, "cmat_d")
    check_and_reset_inverse_and_det!(wf.cmat_f, "cmat_f")
end

function assess_move!(wf::BosonCFLWavefunction, move::ParticlesMove{LatticeParticleMove})
    # @assert is_valid_move(move, wf)

    @assert length(move) == 1 # for now
    if is_occupied(wf.config, move[1].destination, 1)
        wf.cmat_d.current_state = :COLUMNS_UPDATE_PENDING
        wf.cmat_f.current_state = :COLUMNS_UPDATE_PENDING
        probrat = 0.
    else
        detrat_from_columns_update!(wf.cmat_d, move[1].particle.index,
            columns_for_detmat(wf.orbitals_d, move[1].destination))
        detrat_from_columns_update!(wf.cmat_f, move[1].particle.index,
            columns_for_detmat(wf.orbitals_f, move[1].destination))
        probrat = abs2(get_detrat(wf.cmat_d)) * abs2(get_detrat(wf.cmat_f))
    end

    wf.proposed_move = move

    return probrat
end

function finish_move!(wf::BosonCFLWavefunction)
    update_configuration!(wf.config, wf.proposed_move)
    finish_columns_update!(wf.cmat_d)
    finish_columns_update!(wf.cmat_f)
end

function cancel_move!(wf::BosonCFLWavefunction)
    cancel_columns_update!(wf.cmat_d)
    cancel_columns_update!(wf.cmat_f)
end

###############################################################################
#################### CFL = det(nu=1) * det(nu=1) * det(FS) ####################
###############################################################################

type CFLWavefunction{Td1,Td2,Tf} <: HardcoreWavefunction
    orbitals_d1::FreeFermionOrbitals{Td1}
    orbitals_d2::FreeFermionOrbitals{Td2}
    orbitals_f::FreeFermionOrbitals{Tf}

    config::LatticeParticlesConfiguration

    cmat_d1::CeperleyMatrix{Td1}
    cmat_d2::CeperleyMatrix{Td2}
    cmat_f::CeperleyMatrix{Tf}

    neighbor_table::Vector{Vector{Int}} # this could be a Vector of neighbor tables for different neighbors

    proposed_move::ParticlesMove{LatticeParticleMove} # empty upon initialization

    function CFLWavefunction(orbs_d1::FreeFermionOrbitals{Td1}, orbs_d2::FreeFermionOrbitals{Td2}, orbs_f::FreeFermionOrbitals{Tf},
                             nbs::Vector{Vector{Int}}, Ly_for_init::Int)
        @assert Ly_for_init >= 0
        @assert get_Nsites(orbs_d1) == get_Nsites(orbs_d2) == get_Nsites(orbs_f)
        @assert get_Nfilled(orbs_d1) == get_Nfilled(orbs_d2) == get_Nfilled(orbs_f)
        Nsites = get_Nsites(orbs_d1)
        Nparticles = get_Nfilled(orbs_d1)
        @assert Nparticles <= Nsites
        if !isempty(nbs)
            @assert length(nbs) == Nsites
        end

        # get initial configuration of electrons
        if Ly_for_init == 0
            positions = random_combination(Nsites, Nparticles)
            config = LatticeParticlesConfiguration(Vector{Int}[positions], Nsites)
        else
            # put one particle per rung sequentially until we run out of particles
        end

        # compute Ceperley matrices
        cmat_d1 = CeperleyMatrix(fill_detmat(orbs_d1, config))
        cmat_d2 = CeperleyMatrix(fill_detmat(orbs_d2, config))
        cmat_f = CeperleyMatrix(fill_detmat(orbs_f, config))

        return new(orbs_d1, orbs_d2, orbs_f, config, cmat_d1, cmat_d2, cmat_f, nbs)
    end
end

function CFLWavefunction{Td1,Td2,Tf}(orbs_d1::FreeFermionOrbitals{Td1}, orbs_d2::FreeFermionOrbitals{Td2}, orbs_f::FreeFermionOrbitals{Tf},
                         nbs::Vector{Vector{Int}}, Ly_for_init::Int)
    return CFLWavefunction{eltype(orbs_d1.filled),eltype(orbs_d2.filled),eltype(orbs_f.filled)}(orbs_d1, orbs_d2, orbs_f, nbs, Ly_for_init)
end

function CFLWavefunction{Td1,Td2,Tf}(orbs_d1::FreeFermionOrbitals{Td1}, orbs_d2::FreeFermionOrbitals{Td2}, orbs_f::FreeFermionOrbitals{Tf},
                         nbs::Vector{Vector{Int}})
    return CFLWavefunction{eltype(orbs_d1.filled),eltype(orbs_d2.filled),eltype(orbs_f.filled)}(orbs_d1, orbs_d2, orbs_f, nbs, 0)
end

function CFLWavefunction{Td1,Td2,Tf}(orbs_d1::FreeFermionOrbitals{Td1}, orbs_d2::FreeFermionOrbitals{Td2}, orbs_f::FreeFermionOrbitals{Tf},
                         Ly_for_init::Int)
    return CFLWavefunction{eltype(orbs_d1.filled),eltype(orbs_d2.filled),eltype(orbs_f.filled)}(orbs_d1, orbs_d2, orbs_f, Vector{Int}[], Ly_for_init)
end

function CFLWavefunction{Td1,Td2,Tf}(orbs_d1::FreeFermionOrbitals{Td1}, orbs_d2::FreeFermionOrbitals{Td2}, orbs_f::FreeFermionOrbitals{Tf})
    return CFLWavefunction{eltype(orbs_d1.filled),eltype(orbs_d2.filled),eltype(orbs_f.filled)}(orbs_d1, orbs_d2, orbs_f, Vector{Int}[], 0)
end

function check_for_numerical_error!(wf::CFLWavefunction)
    check_and_reset_inverse_and_det!(wf.cmat_d1, "cmat_d1")
    check_and_reset_inverse_and_det!(wf.cmat_d2, "cmat_d2")
    check_and_reset_inverse_and_det!(wf.cmat_f, "cmat_f")
end

function assess_move!(wf::CFLWavefunction, move::ParticlesMove{LatticeParticleMove})
    # @assert is_valid_move(move, wf)

    @assert length(move) == 1 # for now
    if is_occupied(wf.config, move[1].destination, 1)
        wf.cmat_d1.current_state = :COLUMNS_UPDATE_PENDING
        wf.cmat_d2.current_state = :COLUMNS_UPDATE_PENDING
        wf.cmat_f.current_state = :COLUMNS_UPDATE_PENDING
        probrat = 0.
    else
        detrat_from_columns_update!(wf.cmat_d1, move[1].particle.index,
            columns_for_detmat(wf.orbitals_d1, move[1].destination))
        detrat_from_columns_update!(wf.cmat_d2, move[1].particle.index,
            columns_for_detmat(wf.orbitals_d2, move[1].destination))
        detrat_from_columns_update!(wf.cmat_f, move[1].particle.index,
            columns_for_detmat(wf.orbitals_f, move[1].destination))
        probrat = abs2(get_detrat(wf.cmat_d1)) * abs2(get_detrat(wf.cmat_d2)) * abs2(get_detrat(wf.cmat_f))
    end

    wf.proposed_move = move

    return probrat
end

function finish_move!(wf::CFLWavefunction)
    update_configuration!(wf.config, wf.proposed_move)
    finish_columns_update!(wf.cmat_d1)
    finish_columns_update!(wf.cmat_d2)
    finish_columns_update!(wf.cmat_f)
end

function cancel_move!(wf::CFLWavefunction)
    cancel_columns_update!(wf.cmat_d1)
    cancel_columns_update!(wf.cmat_d2)
    cancel_columns_update!(wf.cmat_f)
end

#########################################################################################
#################### RVB [single (1) species via p-h transformation] ####################
#########################################################################################

type RVBWavefunction1{T} <: HardcoreWavefunction
    bcs::SpinfulBCS{T}
    config::LatticeParticlesConfiguration
    cmat::CeperleyMatrix{T}
    neighbor_table::Vector{Vector{Int}}
    proposed_move::ParticlesMove{LatticeParticleMove}

    function RVBWavefunction1(bcs::SpinfulBCS{T}, nbs::Vector{Vector{Int}})
        Nsites = get_Nsites(bcs)
        @assert iseven(Nsites) # otherwise spin system will necessarily be polarized, which I'm assuming not to be the case below
        Nparticles = convert(Int, Nsites/2) # = number of effective hard-core bosons
        if !isempty(nbs)
            @assert length(nbs) == Nsites
        end

        positions = random_combination(Nsites, Nparticles) # positions of up spinons
        config = LatticeParticlesConfiguration(Vector{Int}[positions], Nsites)

        cmat = CeperleyMatrix(fill_detmat_for_SL(bcs, config))

        return new(bcs, config, cmat, nbs)
    end
end

function RVBWavefunction1(bcs::SpinfulBCS, nbs::Vector{Vector{Int}})
    return RVBWavefunction1{eltype(bcs.evecs)}(bcs, nbs)
end

function RVBWavefunction1(bcs::SpinfulBCS)
    return RVBWavefunction1{eltype(bcs.evecs)}(bcs, Vector{Int}[])
end

function check_for_numerical_error!(wf::RVBWavefunction1)
    check_and_reset_inverse_and_det!(wf.cmat, "cmat")
end

function assess_move!(wf::RVBWavefunction1, move::ParticlesMove{LatticeParticleMove})
    # @assert is_valid_move(move, wf)

    @assert length(move) == 1 # for now
    if is_occupied(wf.config, move[1].destination, 1)
        wf.cmat.current_state = :COLUMNS_UPDATE_PENDING
        probrat = 0.
    else
        detrat_from_columns_update!(wf.cmat, [2move[1].particle.index - 1, 2move[1].particle.index],
            columns_for_detmat_for_SL(wf.bcs, move[1].destination))
        probrat = abs2(get_detrat(wf.cmat))
    end

    wf.proposed_move = move

    return probrat
end

function finish_move!(wf::RVBWavefunction1)
    update_configuration!(wf.config, wf.proposed_move)
    finish_columns_update!(wf.cmat)
end

function cancel_move!(wf::RVBWavefunction1)
    cancel_columns_update!(wf.cmat)
end


text:
export outvtk_mesh,outvtk_scalar,outvtk_vector

# output mesh to VTK file for visualization

function outvtk_mesh(fname::String,g::Geometry,t::AbstractTopology)

  topo_ids    = Dict(VertexTopology=>1,IntervalTopology=>2,TriangleTopology=>3,TetrahedronTopology=>4)

  topo_id       = topo_ids[typeof(t)]
  paraview_type = [1;3;5;10][topo_id]
  paraview_len  = [1;2;3;4][topo_id]

  fid = open(fname,"w");

  # output header
  println(fid,"# vtk DataFile Version 2.0")
  println(fid,"Simplex mesh created from Julia")
  println(fid,"ASCII")

  # create unstructured grid with points
  println(fid,"DATASET UNSTRUCTURED_GRID")
  println(fid,"POINTS $(nvertices(g)) float")

  xx=zeros(Float64,3)

  for k=1:nvertices(g)
    for d=1:dim(g)
      xx[d] = xyz(g,k,d)
    end
    println(fid,"$(xx[1]) $(xx[2]) $(xx[3])")
  end

  # output cell connectivity
  println(fid,"CELLS $(nelement(t)) $((paraview_len+1)*nelement(t))")
  for e in t
    println(fid,"$(paraview_len) $(join(e[2:(1+paraview_len)].-1," "))")
  end

  # output cell type 10 = tetraeder
  println(fid,"CELL_TYPES $(nelement(t))")
  for i=1:nelement(t)
    println(fid,"$(paraview_type)")
  end

  # choose a) scalar or b) vector output
  # a) println(fid,"SCALARS $(nameofthatscalar) float 1") # scalar
  # b) println(fid,"VECTORS vectors float") # vector
  # println(fid,"LOOKUP_TABLE default")
  # for i=1:nvertices(g)
  #   a) println(fid,"$(u[i])") # scalar
  #   b) println(fid,"$(ux[i]) $(uy[i]) $(uz[i])") # vector
  # end

  close(fid)

end

function outvtk_scalar(fname::String,g::Geometry,t::AbstractTopology,sname::String,u::Vector{Float64})

  topo_ids    = Dict(VertexTopology=>1,IntervalTopology=>2,TriangleTopology=>3,TetrahedronTopology=>4)
  topo_id       = topo_ids[typeof(t)]
  paraview_type = [1;3;5;10][topo_id]
  paraview_len  = [1;2;3;4][topo_id]

  fid = open(fname,"w");

  # output header
  println(fid,"# vtk DataFile Version 2.0")
  println(fid,"Simplex mesh created from Julia")
  println(fid,"ASCII")

  # create unstructured grid with points
  println(fid,"DATASET UNSTRUCTURED_GRID")
  println(fid,"POINTS $(nvertices(g)) float")

  xx=zeros(Float64,3)

  for k=1:nvertices(g)
    for d=1:dim(g)
      xx[d] = xyz(g,k,d)
    end
    println(fid,"$(xx[1]) $(xx[2]) $(xx[3])")
  end

  # output cell connectivity
  println(fid,"CELLS $(nelement(t)) $((paraview_len+1)*nelement(t))")
  for e in t
    println(fid,"$(paraview_len) $(join(e[2:(1+paraview_len)].-1," "))")
  end

  # output cell type 10 = tetraeder
  println(fid,"CELL_TYPES $(nelement(t))")
  for i=1:nelement(t)
    println(fid,"$(paraview_type)")
  end
  println(fid,"POINT_DATA $(length(u))")

  println(fid,"SCALARS $(sname) float 1") # scalar
  println(fid,"LOOKUP_TABLE default")
  for i=1:length(u)
    println(fid,"$(u[i])")
  end

  close(fid)

end

function outvtk_vector(fname::String,g::Geometry,t::AbstractTopology,sname::String,ux::Vector{Float64},uy::Vector{Float64},uz::Vector{Float64})

  topo_ids    = Dict(VertexTopology=>1,IntervalTopology=>2,TriangleTopology=>3,TetrahedronTopology=>4)
  topo_id       = topo_ids[typeof(t)]
  paraview_type = [1;3;5;10][topo_id]
  paraview_len  = [1;2;3;4][topo_id]

  fid = open(fname,"w");

  # output header
  println(fid,"# vtk DataFile Version 2.0")
  println(fid,"Simplex mesh created from Julia")
  println(fid,"ASCII")

  # create unstructured grid with points
  println(fid,"DATASET UNSTRUCTURED_GRID")
  println(fid,"POINTS $(nvertices(g)) float")

  xx=zeros(Float64,3)

  for k=1:nvertices(g)
    for d=1:dim(g)
      xx[d] = xyz(g,k,d)
    end
    println(fid,"$(xx[1]) $(xx[2]) $(xx[3])")
  end

  # output cell connectivity
  println(fid,"CELLS $(nelement(t)) $((paraview_len+1)*nelement(t))")
  for e in t
    println(fid,"$(paraview_len) $(join(e[2:(1+paraview_len)].-1," "))")
  end

  # output cell type 10 = tetraeder
  println(fid,"CELL_TYPES $(nelement(t))")
  for i=1:nelement(t)
    println(fid,"$(paraview_type)")
  end
  println(fid,"POINT_DATA $(length(ux))")

  println(fid,"VECTORS $(sname) float") # scalar
  # println(fid,"LOOKUP_TABLE default")
  for i=1:length(ux)
    println(fid,"$(ux[i]) $(uy[i]) $(uz[i])")
  end

  close(fid)

end

ext: jl

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
# Luca Brugnolini
# Ph.D. in Economics, Law and Institutions
# I am starting programming in Julia
# Base my code on Kilian (2009) Local Projection
# I modify the code to implement it in Julia and to improve performance
# First Draft 30/11/2
# This Version 02/12/2014
# ** Denotes some performance improvements I have implemented. If you'll find some errors try to check these parts.

#cd("C:/Users/asus/Desktop/Luke/UniversitÃ /PhD/My Paper/About Local Projection/Local Projection/MATLAB/Killian Original Code")
z = readcsv("data.csv");
tic()
DY    = z[2:end,1]; # CFNAI
DP    = (log(z[2:end,2])-log(z[1:end-1,2]))*1200;  # D(log(CPI))
DPCOM = ((log(z[2:end,3])-log(z[2:end,2])) - (log(z[1:end-1,3])-log(z[1:end-1,2])))*1200; # D(log(real PCOM)): CRB
FF    = z[2:end,4];  #FedFound

# DEMINING FOR BC BOOTSTRAP
# FULL SAMPLE (1970-2007)
DY = DY - mean(DY);
DP = DP - mean(DP);
DPCOM = DPCOM - mean(DPCOM);
FF = FF - mean(FF);
z = [DY DP DPCOM FF];
p = 12;
H = 24;

# DGP: 4-Variable VAR(12)- Output, Inflation, Real Commodity-P inf, FFrate
(irftrue, std, P, Btrue) = irfvar(z, p, H);             # IRF for DGP

# SIMULATION
(T,K) = size(z);
nmc   = 1;                                 # number of Monte Carlo trials
throw = 300;                                 # initial obs to be discarded
pbar  = 12;                                  # max lag order for AIC

# Initialization
BIASv = zeros(K^2,H+1);   MSEv  = zeros(K^2,H+1);   # point estimates
BIASl = zeros(K^2,H+1);   MSEl  = zeros(K^2,H+1);

CVGvd = zeros(K^2,H+1);   LGTvd = zeros(K^2,H+1);   # inference - delta
CVGld = zeros(K^2,H+1);   LGTld = zeros(K^2,H+1);

CVGvb = zeros(K^2,H+1);   LGTvb = zeros(K^2,H+1);   # inference - bootstrap
CVGlb = zeros(K^2,H+1);   LGTlb = zeros(K^2,H+1);

# Monte-Carlo Trials
for n = 1:nmc                                       # loop for MC traials

#---GENERATE DATA
    y = zeros(T+throw,K);
    u = randn(T+throw,K);

    for i = p+1:T+throw
    y[i,:] = (Btrue[:,1]+P*u[i,:]')';

        for jj = 1:p
        y[i,:] = y[i,:] + (Btrue[:,(jj-1)*K+2:jj*K+1]*y[i-jj,:]')';

        end

     end

    y = y[throw+1:end,:]                          # discarding initial obs

    for i=1:K                                   # demean y
        y[:,i] = y[:,i] - mean(y[:,i]);
    end

#---IRF POINT ESTIMATE AND STDEV ESTIMATED BY VAR AND LP
    pv =  varlagorder(y,pbar);            # lag-order for VAR
    (irfv, stdv, Phat, B, SIGMA, U, rfirfv, COVsig) = irfvar(y, pv, H); # IRF by VAR

    pl = lplagorder(y,pbar,H);                                       # lag-order for LP
    (irfl, stdl, rfirfl) = irflp(y, pl, H, Phat, COVsig);            # IRF by LP

#---BIAS, MSE of IRF ESTIMATES
    BIASv = BIASv + (irfv-irftrue);                    # sum of BIAS
    BIASl = BIASl + (irfl-irftrue);

    MSEv  = MSEv  + (irfv-irftrue).^2;                 # sum of squared error
    MSEl  = MSEl  + (irfl-irftrue).^2;

#---BIAS-CORRECTED BOOTSTRAP INTERVAL FOR VAR: KILLIAN -- for VAR(12)
    # translating estimates into the companion form
    A = [B[:,2:end];[eye(K*(p-1)) zeros(K*(p-1),K)]];  # slope estimate
    V = [B[:,1];zeros(K*(p-1),1)];                     # intercept
    SIGMAc = [SIGMA zeros(K,K*p-K);zeros(K*p-K,K*p)];  # sigma

# Bias correction: if the largest root of the companion matrix
# is less than 1, do BIAS correction
   eigv = abs(eigvals(A));

   if ~ any(eigv.>=1)
        A = asybc(A,SIGMAc,T,K,p);

   end

   A = real(A);

# Bias-corrected estimate is the bootstrap DGP
   CIv = boot([V[1:K,1] A[1:K,1:K*p]], U, y, p, H);
   CILv  = reshape(CIv[1,:]',H+1,K^2)';               # lower bound
   CIHv  = reshape(CIv[2,:]',H+1,K^2)';               # upper bound

#---BLOCK BOOTSTRAP INTERVAL FOR LOCAL PROJECTIONS
   block = 8;                                         # block size
   CIl = bcbootlp(y, rfirfl, SIGMA, pl, block, H);    # bootstrap interval
   CILl = reshape(CIl[1,:]',H+1,K^2)';                # lower bound
   CIHl = reshape(CIl[2,:]',H+1,K^2)';                # upper bound

#---COVERAGE RATES
    for k=1:K^2

        for h=1:H+1

            if (irftrue[k,h]>=irfv[k,h]-1.96*stdv[k,h]) && (irftrue[k,h]<=irfv[k,h]+1.96*stdv[k,h])
               CVGvd[k,h] = CVGvd[k,h]+1;

            end

            if (irftrue[k,h]>=irfl[k,h]-1.96*stdl[k,h]) && (irftrue[k,h]<=irfl[k,h]+1.96*stdl[k,h])
               CVGld[k,h] = CVGld[k,h]+1;

            end

            if (irftrue[k,h]>=CILv[k,h]) && (irftrue[k,h]<=CIHv[k,h])
               CVGvb[k,h] = CVGvb[k,h]+1;

            end

            if (irftrue[k,h]>=CILl[k,h]) && (irftrue[k,h]<=CIHl[k,h])
               CVGlb[k,h] = CVGlb[k,h]+1;

            end

        end

    end

#---AVERAGE LENGTHS
    LGTvd = LGTvd + 1.96*2*stdv;
    LGTld = LGTld + 1.96*2*stdl;
    LGTvb = LGTvb + abs(CIHv-CILv);
    LGTlb = LGTlb + abs(CIHl-CILl);
end                                              # end loop for MC trials

BIASv = BIASv/nmc;                               # BIAS
BIASl = BIASl/nmc;
MSEv  = MSEv /nmc;                               # MSE
MSEl  = MSEl /nmc;
STDv  = sqrt(MSEv-BIASv.^2);                     # stdev
STDl  = sqrt(MSEl-BIASl.^2);

CVGvd  = CVGvd/nmc;                              # coverage rates
CVGld  = CVGld/nmc;
CVGvb  = CVGvb/nmc;
CVGlb  = CVGlb/nmc;

LGTvd  = LGTvd/nmc;                              # average lengths
LGTld  = LGTld/nmc;
LGTvb  = LGTvb/nmc;
LGTlb  = LGTlb/nmc;
toc()


text:
module Utility

import GreyDecision.GreyNumbers: GreyNumber

mean(x) = sum(x) / length(x)

# euclidean(g1, g2) -> GreyNumber ?
# euclidean(g1, g2) -> Real ?
function euclidean(g1::Array{GreyNumber{T}, 1}, g2::Array{GreyNumber{T}, 1})::GreyNumber where {T <: Real}
    return (g1 .- g2).^2.0 |> sum |> sqrt
end

function euclidean(v1::Array{GreyNumber{T}, 1})::GreyNumber where {T <: Real}
    n = length(v1)
    v2 = map(x -> GreyNumber(zero(T), zero(T)), 1:n)
    return euclidean(v1, v2)
end

function apply_columns(fs::Array{Function,1}, data::Array{GreyNumber{T},2}) where {T <: Real}
    _, m = size(data)
    return [fs[i](data[:,i]) for i in 1:m]
end

function apply_columns(f::Function, data::Array{GreyNumber{T},2}) where T <: Real
    _, p = size(data)
    return [f(data[:,i]) for i in 1:p]
end

function apply_rows(f::Function, data::Array{GreyNumber{T},2}) where T <: Real
    n, _ = size(data)
    return [f(data[i,:]) for i in 1:n]
end

function colmins(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_columns(minimum, data)
end

function colmaxs(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_columns(maximum, data)
end

function colsums(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_columns(sum, data)
end

function colmeans(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_columns(mean, data)
end

function rowmins(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_rows(minimum, data)
end

function rowmaxs(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_rows(maximum, data)
end

function rowsums(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_rows(sum, data)
end

function rowmeans(data::Array{GreyNumber{T}, 2}) where T <: Real
    return apply_rows(mean, data)
end

function Base.zero(::Type{GreyNumber{T}}) where T
    return GreyNumber(zero(T), zero(T))
end

function Base.zero(::Type{GreyNumber})
    return GreyNumber(0.0, 0.0)
end

function Base.zero(g::GreyNumber{T})::GreyNumber{T} where T
    return GreyNumber(zero(T), zero(T))
end

function Base.zeros(::Type{GreyNumber{T}}, n::Int64)::Array{GreyNumber{T}, 1} where T
    gs = Array{GreyNumber{T}, 1}(undef, n)
    for i in 1:n
        gs[i] = GreyNumber(zero(T), zero(T))
    end
    return gs
end

function Base.zeros(::Type{GreyNumber}, n::Int64, m::Int64)::Array{GreyNumber, 2}
    zeros(GreyNumber{Float64}, n, m)
end

function Base.zeros(::Type{GreyNumber}, t::Tuple{Int64, Int64})::Array{GreyNumber, 2}
    zeros(GreyNumber, first(t), last(t))
end

function Base.one(::Type{GreyNumber{T}}) where T
    return GreyNumber(one(T), one(T))
end

function Base.one(::Type{GreyNumber})
    return GreyNumber(1.0, 1.0)
end

function Base.one(g::GreyNumber{T})::GreyNumber{T} where T
    return GreyNumber(one(T), one(T))
end

function Base.ones(::Type{GreyNumber{T}}, n::Int64)::Array{GreyNumber{T}, 1} where T
    gs = Array{GreyNumber{T}, 1}(undef, n)
    for i in 1:n
        gs[i] = GreyNumber(one(T), one(T))
    end
    return gs
end

function Base.ones(::Type{GreyNumber}, n::Int64, m::Int64)::Array{GreyNumber{Float64}, 2}
    ones(GreyNumber{Float64}, n, m)
end

function Base.ones(::Type{GreyNumber}, t::Tuple{Int64, Int64})::Array{GreyNumber{Float64}, 2}
    ones(GreyNumber, first(t), last(t))
end

function Base.isvalid(g::GreyNumber)::Bool
    return g.a <= g.b
end

function Base.convert(::Type{Array{T, 1}}, g::GreyNumber)::Array{T, 1} where T
    arr = Array{T, 1}(undef, 2)
    arr[1] = min(g.a, g.b)
    arr[2] = max(g.a, g.b)
    return arr
end

function reverseminmax(fns::Array{Function,1})::Array{Function,1}
    newfs = map(x -> if x == minimum maximum else minimum end, fns)
    return newfs
end

function makeminmax(fns::Array{K,1} where K)::Array{Function,1}
    return convert(Array{Function,1}, fns)
end

function normalize(v1::Array{GreyNumber{T},1})::Array{GreyNumber{T},1} where {T <: Real}
    disttoorigin = euclidean(v1)
    return map(x -> x / disttoorigin, v1)
end

function normalize(datamat::Array{GreyNumber{T},2})::Array{GreyNumber{T},2} where {T <: Real}
    newdata = similar(datamat)
    _, p = size(datamat)
    for i in 1:p
        newdata[:, i] = normalize(datamat[:, i])
    end
    return newdata
end

function unitize(v::Array{T,1})::Array{Float64,1} where {T <: Number}
    return v ./ sum(v)
end

function weightize(data::Array{GreyNumber{T}, 2}, w::Array{Float64, 1})::Array{GreyNumber{T}, 2} where {T <: Real}
    n, p = size(data)
    newdata = Array{GreyNumber{T}, 2}(undef, n, p)
    for i in 1:p
        newdata[:, i] = map(x -> w[i] * x, data[:, i])
    end
    return newdata
end

function Base.isapprox(x::GreyNumber{T}, y::GreyNumber{T}; atol::Real=0)::Bool where {T <: Real}
    return isapprox(first(x), first(y), atol = atol) && isapprox(last(x), last(y), atol = atol)
end

end # end of module Utility


text:
# nonparametric regression, using pre-generated weights
# can be local constant, local linear, or local quadratic
#
# y is n X p
# x is n X dimx
# xeval is neval X dimx
# weights is n X neval (different weights for each eval. point)
#
# weights should sum to one by columns, they are typical
# nonparametric weights from a kernel.

function npreg()
    println("npreg(), with no arguments, runs Pkg.dir/Econometrics/examples/npreg_example.jl")
    println("examine that file for information on how to call npreg.jl")
    include(Pkg.dir()*"/Econometrics/examples/npreg_example.jl")
end

function npreg(y::Array{Float64,2}, x::Array{Float64,2}, xeval::Array{Float64,2}, weights::Array{Float64,2}, order::Int64=1, do_median::Bool=false, do_ci::Bool=false)
    weights = sqrt.(weights)
    neval, dimx = size(xeval)
    n, dimy = size(y)
    # local constant
    if order==0
        X = ones(n,1)
        Xeval = ones(neval,1)
    elseif order==1
    # local linear
        X = [ones(n,1) x]
        Xeval = [ones(neval,1) xeval]
    else
    # local quadratic
        stacked = [x; xeval]
        # cross products
        CP = Array{Float64}(n+neval, Int((dimx-1)*dimx/2))
        cpind = 0
        for i = 1:dimx-1
            for j = (i+1):dimx
                cpind += 1
                CP[:,cpind] = stacked[:,[i]].*stacked[:,[j]]
            end
        end
        ZZ = [ones(n+neval,1) stacked CP]
        X = view(ZZ,1:n,:)
        Xeval = view(ZZ,(n+1):n+neval,:)
    end
    # do the fit
    yhat = Array{Float64}(neval, dimy)
    for i = 1:neval
        WX = weights[:,i] .* X
        Wy = weights[:,i] .* y
        b = WX\Wy
        yhat[i,:] = Xeval[[i],:]*b
    end
    return yhat
end


text:
# This is not the usual style of trampoline (thunk-returning functions) and
# maybe there's a better name (pop-call-append?). We do not transform the
# recursion fully to a loop for work-stealing scheduler. Until the first
# base case (initial decent), the recursion is used as usual (so still consumes
# O(log2 n) stack space). The required amount of stack space is comparable to
# other divide-and-conquer implementations. The idea/assumption is that it
# _might_ be a good idea to keep the recursion so that the compiler can infer
# types in many cases (TODO: check this). Once the call stack is constructed by
# the recursion, the continuations in the chain (cactus stack) are evaluated in
# a loop ("trampoline"). That said, it'd be interesting to see if the standard
# trampoline has some performance/implementation advantages.

and_finally(@nospecialize f) = listof(Function, f)

before(@nospecialize(f), chain::List{Function}) = Cons{Function}(f, chain)

function trampoline_fallback(chain, x)
    chain isa Nil && return x
    while true
        f = car(chain)
        y = f(x)
        if y === nothing
            return
        elseif y isa Some
            x = something(y)
            chain2 = cdr(chain)
            chain2 isa Nil && return
            chain = chain2
        elseif y isa Tuple{Cons{Function},Any}
            chain, x = y
        else
            unexpected(y)
        end
    end
end

@noinline unexpected(x) = error("unexpected type: $x")

"""
Semantically equivalent to `trampoline_fallback` but tries to type-stabilize.
"""
function trampoline(chain, x)
    chain isa Nil && return x
    # return trampoline_fallback(chain, x)
    trampoline_stabilizing(chain, x, typeof(car(chain)), typeof(x))
end

@generated function trampoline_stabilizing(chain, x, ::Type{TF}, ::Type{TX}) where {TF,TX}
    function unwrap_types(@nospecialize T)
        types = Any[]
        while true
            if T isa Union
                push!(types, T.a)
                T = T.b
            else
                push!(types, T)
                break
            end
        end
        return types
    end
    ftypes = unwrap_types(TF)
    xtypes = unwrap_types(TX)

    if any(!Base.isconcretetype, ftypes) || any(!Base.isconcretetype, xtypes)
        # tuple of union? falling back...
        return :(trampoline_fallback(chain, x))
    end

    if length(ftypes) * length(xtypes) < 16
        # Recurse to re-stabilize:
        otherwise = quote
            return trampoline_stabilizing(chain, x, Union{typeof(f),TF}, Union{typeof(x),TX})
        end
    else
        # Too many recursions. Just keep going:
        otherwise = :(y = f(x))
    end

    fbranches = foldr(ftypes, init = otherwise) do ft, ex
        quote
            if f isa $ft
                y = f(x)
            else
                $ex
            end
        end
    end
    xbranches = foldr(xtypes, init = otherwise) do xt, ex
        quote
            if x isa $xt
                $fbranches
            else
                $ex
            end
        end
    end

    quote
        chain isa Nil && return x
        while true
            f = car(chain)
            $xbranches  # y = f(x)
            if y === nothing
                return
            elseif y isa Some
                x = something(y)
                chain2 = cdr(chain)
                chain2 isa Nil && return
                chain = chain2
            elseif y isa Tuple{Cons{Function},Any}
                chain, x = y
            else
                unexpected(y)
            end
        end
    end
end


text:
using Test, Yao.EasyBuild

@testset "block_extension" begin
    include("block_extension/blocks.jl")
end

@testset "variational_circuit" begin
    include("variational_circuit.jl")
end

@testset "hamiltonians" begin
    include("hamiltonians.jl")
end

@testset "supremacy_circuit" begin
    include("supremacy_circuit.jl")
end

@testset "google53" begin
    include("google53.jl")
end

@testset "general_U4" begin
    include("general_U4.jl")
end

@testset "phaseestimation" begin
    include("phaseestimation.jl")
end

@testset "hadamardtest" begin
    include("hadamardtest.jl")
end


text:
abstract type Proposal{P} end

struct StaticProposal{P} <: Proposal{P}
    proposal::P
end

struct RandomWalkProposal{P} <: Proposal{P}
    proposal::P
end

# Random draws
Base.rand(p::Proposal, args...) = rand(Random.GLOBAL_RNG, p, args...)
Base.rand(rng::Random.AbstractRNG, p::Proposal{<:Distribution}) = rand(rng, p.proposal)
function Base.rand(rng::Random.AbstractRNG, p::Proposal{<:AbstractArray})
    return map(x -> rand(rng, x), p.proposal)
end

# Densities
Distributions.logpdf(p::Proposal{<:Distribution}, v) = logpdf(p.proposal, v)
function Distributions.logpdf(p::Proposal{<:AbstractArray}, v)
    # `mapreduce` with multiple iterators requires Julia 1.2 or later
    return mapreduce(((pi, vi),) -> logpdf(pi, vi), +, zip(p.proposal, v))
end

###############
# Random Walk #
###############

function propose(rng::Random.AbstractRNG, p::RandomWalkProposal, m::DensityModel)
    return propose(rng, StaticProposal(p.proposal), m)
end

function propose(
    rng::Random.AbstractRNG,
    proposal::RandomWalkProposal{<:Union{Distribution,AbstractArray}},
    model::DensityModel,
    t
)
    return t + rand(rng, proposal)
end

function q(
    proposal::RandomWalkProposal{<:Union{Distribution,AbstractArray}},
    t,
    t_cond
)
    return logpdf(proposal, t - t_cond)
end

##########
# Static #
##########

function propose(
    rng::Random.AbstractRNG,
    proposal::StaticProposal{<:Union{Distribution,AbstractArray}},
    model::DensityModel,
    t=nothing
)
    return rand(rng, proposal)
end

function q(
    proposal::StaticProposal{<:Union{Distribution,AbstractArray}},
    t,
    t_cond
)
    return logpdf(proposal, t)
end

############
# Function #
############

# function definition with abstract types requires Julia 1.3 or later
for T in (StaticProposal, RandomWalkProposal)
    @eval begin
        (p::$T{<:Function})() = $T(p.proposal())
        (p::$T{<:Function})(t) = $T(p.proposal(t))
    end
end

function propose(
    rng::Random.AbstractRNG,
    proposal::Proposal{<:Function},
    model::DensityModel
)
    return propose(rng, proposal(), model)
end

function propose(
    rng::Random.AbstractRNG,
    proposal::Proposal{<:Function},
    model::DensityModel,
    t
)
    return propose(rng, proposal(t), model)
end

function q(
    proposal::Proposal{<:Function},
    t,
    t_cond
)
    return q(proposal(t_cond), t, t_cond)
end


text:
using LinearAlgebra, StaticArrays, Attitude, SatelliteDynamics
using Infiltrator, ProgressMeter

function initial_conditions()

    # epc0 = Epoch(2019, 1, 1, 12, 0, 0, 0.0)

    # Declare initial state in terms of osculating orbital elements
    rp = R_EARTH + 200e3;ra = R_EARTH + 35000e3
    a = (ra+rp)/2
    e = (ra-rp)/(ra+rp)
    i = deg2rad(27)
    # i = 0
    ω = deg2rad(34)
    Ω = deg2rad(133)
    M = deg2rad(0)
    oe0  = [a,e,i,Ω,ω,M]

    # Convert osculating elements to Cartesean state
    eci0 = sOSCtoCART(oe0, use_degrees=false)

    # modified equinoctial elements
    p = a*(1-e^2)
    f = e*cos(ω + Ω)
    g = e*sin(ω + Ω)
    h = tan(i/2)*cos(Ω)
    k = tan(i/2)*sin(Ω)
    L = wrap_to_2pi(ω + Ω + M)

    eoe0 = [p;f;g;h;k;L]

    r_eci = eci0[1:3]
    v_eci = eci0[4:6]

    p_0 = u_from_x(r_eci)                          # u
    p_prime_0 = uprime_from_xdot(v_eci,p_0)     # du/ds
    h_0 = GM_EARTH/norm(r_eci) - .5*norm(v_eci)^2         # m^2/s^2
    t_0 = 0

    ks0 = [p_0;p_prime_0;h_0;t_0]

    return eci0, eoe0, ks0
end
function rk4(f::Function, x_n, dt)
    k1 = dt*f(x_n)
    k2 = dt*f(x_n+k1/2)
    k3 = dt*f(x_n+k2/2)
    k4 = dt*f(x_n+k3)
    return (x_n + (1/6)*(k1+ 2*k2 + 2*k3 + k4))
end
function fode_j2(r_eci)
    r = norm(r_eci)
    x,y,z = r_eci
    J2 = J2_EARTH
    μ = GM_EARTH

    # precompute repeated stuff
    Re_r_sqr = 1.5*J2*(R_EARTH/r)^2
    five_z_sqr = 5*z^2/r^2

    return  SA[(μ/r^3)*x*(Re_r_sqr*(five_z_sqr - 1)),
               (μ/r^3)*y*(Re_r_sqr*(five_z_sqr - 1)),
               (μ/r^3)*z*(Re_r_sqr*(five_z_sqr - 3))]
end
function cartesian_dynamics(x)
    r_eci = SA[x[1],x[2],x[3]]
    v_eci = SA[x[4],x[5],x[6]]

    # relavent parameters
    J2 = J2_EARTH
    μ = GM_EARTH

    # eci position stuff
    r = norm(r_eci)
    x,y,z = r_eci

    # precompute repeated stuff
    Re_r_sqr = 1.5*J2*(R_EARTH/r)^2
    five_z_sqr = 5*z^2/r^2

    return  SA[v_eci[1],v_eci[2],v_eci[3],
               (-μ/r^3)*x*(1 - Re_r_sqr*(five_z_sqr - 1)),
               (-μ/r^3)*y*(1 - Re_r_sqr*(five_z_sqr - 1)),
               (-μ/r^3)*z*(1 - Re_r_sqr*(five_z_sqr - 3))]
end
function J2_from_EOE(x)
    """J2 acceleration in RTN from equinoctial orbital elements.
    Args:
        p: a*(1-e^2) semi-latus rectum       (m)
        f: e*cos(ω + Ω)                      (non-dimensional)
        g: e*sin(ω + Ω)                      (non-dimensional)
        h: tan(i/2)*cos(Ω)                   (non-dimensional)
        k: tan(i/2)*sin(Ω)                   (non-dimensional)
        L: ω + Ω + θ argument of latitude    (rad)
    Returns:
        a_rtn: acceleration in the RTN frame  (m/s)
    """

    # unpack state
    p,f,g,h,k,L = x
    J2 =  J2_EARTH

    # precompute trig functions
    sinL = sin(L)
    cosL = cos(L)

    # constants for EOE
    mu = GM_EARTH
    w = 1+f*cosL + g*sinL
    r = p/w
    Re = R_EARTH

    # accelerations
    a_r = -((3*mu*J2*Re^2)/(2*r^4))*(1 -   (  (12*(h*sinL-k*cosL)^2 )
                                                    /   ((1+h^2 + k^2)^2))  )
    a_t = -((12*mu*J2*Re^2)/(r^4))*((  ((h*sinL - k*cosL)*(h*cosL+k*sinL) )
                                                    /  ((1+h^2 + k^2)^2))  )
    a_n = -((6*mu*J2*Re^2)/(r^4))*((  ((h*sinL - k*cosL)*(1-h^2-k^2) )
                                                    /  ((1+h^2 + k^2)^2))  )
    return a_r, a_t, a_n
end
function eoe_dynamics(x)
    """Equations of motion for equinoctial elements in 2-body problem.
    Args:
        epc: time(::epoch)
        x: state
            - p: a*(1-e^2) semi-latus rectum       (m)
            - f: e*cos(ω + Ω)                      (non-dimensional)
            - g: e*sin(ω + Ω)                      (non-dimensional)
            - h: tan(i/2)*cos(Ω)                   (non-dimensional)
            - k: tan(i/2)*sin(Ω)                   (non-dimensional)
            - L: ω + Ω + θ argument of latitude    (rad)
            - mass: spacecraft mass                (kg)
        u: thruster acceleration                   (m/s^2)
        thruster_on: bool for thruster             ()
    Returns:
        ̇x: state derivative
    """

    # unpack state
    p,f,g,h,k,L = x

    # equinoctial parameters for GVE's
    sinL = sin(L)
    cosL = cos(L)
    alpha_2 = h^2 - k^2
    s2 = 1+h^2 + k^2
    w = 1+f*cosL + g*sinL
    r = p/w

    # precompute to save calculations
    sqrtpu = sqrt(p/GM_EARTH)

    # J2
    u_r, u_t, u_n = J2_from_EOE(x)

    # GVE's for equinoctial orbital elements
    p_dot = (2*p/w)*sqrtpu*u_t
    f_dot = sqrtpu*( u_r*sinL    +    ((w+1)*cosL + f)*u_t/w    -
                                            (h*sinL - k*cosL)*g*u_n/w  )
    g_dot = sqrtpu*(-u_r*cosL   +    ((w+1)*sinL + g)*u_t/w    +
                                            (h*sinL - k*cosL)*f*u_n/w  )
    h_dot = sqrtpu*(s2*u_n/(2*w))*cosL
    k_dot = sqrtpu*(s2*u_n/(2*w))*sinL
    L_dot = sqrt(GM_EARTH*p)*(w/p)^2 + (1/w)*sqrtpu*(h*sinL - k*cosL)*u_n

    # return state derivative
    return SA[p_dot,f_dot,g_dot,h_dot,k_dot,L_dot]
end
function rv_from_eoe(eoe_vec)
    """Position and Velocity resolved in ECI
    Args:
        p: a*(1-e^2) semi-latus rectum     (m)
        f: e*cos(ω + Ω)                    (non-dimensional)
        g: e*sin(ω + Ω)                    (non-dimensional)
        h: tan(i/2)*cos(Ω)                 (non-dimensional)
        k: tan(i/2)*sin(Ω)                 (non-dimensional)
        L: ω + Ω + θ argument of latitude  (rad)
    Returns:
        r: position vector from eci origin to spacecraft, resolved in eci. (m)
        v: velocity vector of spacecraft in eci, resolved in eci.          (m/s)
    """

    # unpack eoe
    p,f,g,h,k,L = eoe_vec

    # precompute sin and cosine
    sinL = sin(L)
    cosL = cos(L)

    # various eoe specific parameters
    alf2 = h^2 - k^2
    s2 = 1 + h^2 + k^2
    w = 1 + f*cosL + g*sinL
    r = p/w

    # eci position vector
    r = SA[(r/s2) * (  cosL + alf2*cosL + 2*h*k*sinL),
         (r/s2) * (  sinL - alf2*sinL + 2*h*k*cosL),
         (2*r/s2)*(h*sinL -   k *cosL)]

    # ECI velocity vector
    v = sqrt(GM_EARTH/p)*SA[-(1/(s2)) * (sinL + alf2*sinL - 2*h*k*cosL + g - 2*f*h*k + alf2*g),
         -(1/(s2))* (-cosL + alf2*cosL + 2*h*k*sinL - f + 2*g*h*k + alf2*f),
          (2/s2) * (h*cosL +    k*sinL + f*h + g*k)]

    # @infiltrate
    return SA[r[1],r[2],r[3],v[1],v[2],v[3]]
end

function ks_dynamics(x)

    # unpack state
    p = SVector(x[1],x[2],x[3],x[4])
    p_prime = SVector(x[5],x[6],x[7],x[8])
    h = x[9]
    t = x[10]

    # position in eci
    r_eci = x_from_u(p)

    # j2 accel
    a_j2 = fode_j2(r_eci)

    L = L_fx(p)

    # append 0 to end of P_j2 acceleration vector
    P_j2 = SA[a_j2[1],a_j2[2],a_j2[3],0]

    p_dp = -(h/2)*p + ((dot(p,p))/(2))*L'*P_j2

    return SVector(p_prime[1],p_prime[2],p_prime[3],p_prime[4],
                   p_dp[1],p_dp[2],p_dp[3],p_dp[4],
                   -2*dot(p_prime,L'*P_j2),
                   dot(p,p))
end
function rv_from_ks(p,p_prime)
    x = x_from_u(p)
    xdot = xdot_from_u(p,p_prime)
    return SA[x[1],x[2],x[3],xdot[1],xdot[2],xdot[3]]
end
function ϵ(rv)
    r = SA[rv[1],rv[2],rv[3]]
    v = SA[rv[4],rv[5],rv[6]]
    return dot(v,v)/2 - GM_EARTH/norm(r)
end
# function driver()
#     eci0, eoe0, ks0 = initial_conditions()
#
#     tf = 1*24*3600
#     N = 100000
#     dt = tf/N
#
#     X_eci = [@SVector zeros(6) for i = 1:N]
#     X_eci[1] = SVector{6}(eci0)
#     X_eoe = [@SVector zeros(6) for i = 1:N]
#     X_eoe[1] = eoe0
#     X_ks = [@SVector zeros(10) for i = 1:N]
#     X_ks[1] = ks0
#
#     # one day in fictitious time (this was guess and check)
#     sf = 0.003709551289
#     ds = sf/N
#
#     @inbounds for i = 1:(N-1)
#         X_eci[i+1] = rk4(cartesian_dynamics,X_eci[i],dt)
#         X_eoe[i+1] = rk4(eoe_dynamics,X_eoe[i],dt)
#         X_ks[i+1] = rk4(ks_dynamics,X_ks[i],ds)
#     end
#
#     ksm = mat_from_vec(X_ks)
#
#     ps = ksm[1:4,:]
#
#     # mat"
#     # figure
#     # hold on
#     # plot($ps')
#     # hold off
#     # "
#
#     # @show tf - X_ks[end][10]
#     X_eoe_eci = [rv_from_eoe(X_eoe[i]) for i = 1:length(X_eoe)]
#     X_ks_eci = [rv_from_ks(X_ks[i][1:4],X_ks[i][5:8]) for i = 1:length(X_ks)]
#
#     eci_m = mat_from_vec(X_eci)
#     ks_m = mat_from_vec(X_ks_eci)
#
#     # energies
#     ks_ϵ = ϵ(X_ks_eci[end])
#     eci_ϵ = ϵ(X_eci[end])
#     eoe_ϵ = ϵ(X_eoe_eci[end])
#
#     @show ks_ϵ
#     @show eci_ϵ
#     @show eoe_ϵ
#
#
#     # mat"
#     # figure
#     # hold on
#     # plot3($eci_m(1,:),$eci_m(2,:),$eci_m(3,:))
#     # plot3($ks_m(1,:),$ks_m(2,:),$ks_m(3,:))
#     # hold off
#     # "
#
#
#     # # @show "yes"
#     # eci_m = mat_from_vec(X_eci)
#     # eoe_m = mat_from_vec(X_eoe_eci)
#     # mat"
#     # figure
#     # hold on
#     # plot3($eci_m(1,:),$eci_m(2,:),$eci_m(3,:))
#     # plot3($eoe_m(1,:),$eoe_m(2,:),$eoe_m(3,:))
#     # hold off
#     # "
#     #
#     #
#     # @show norm(X_eci[end] - X_eoe_eci[end])
#     # @show norm(X_eci[end] - X_ks_eci[end])
# end
function driver(N,eci0, eoe0, ks0 )
    # eci0, eoe0, ks0 = initial_conditions()

    # N = 100

    X_eci = [@SVector zeros(6) for i = 1:N]
    X_eci[1] = SVector{6}(eci0)
    X_eoe = [@SVector zeros(6) for i = 1:N]
    X_eoe[1] = eoe0
    X_ks = [@SVector zeros(10) for i = 1:N]
    X_ks[1] = ks0

    # one day in fictitious time (this was guess and check)
    sf = 10*0.0003709551289
    ds = sf/(N)

    @inbounds for i = 1:(N-1)
        X_ks[i+1] = rk4(ks_dynamics,X_ks[i],ds)
    end
    tf = X_ks[end][10]
    # @show tf
    tf += tf/N
    dt = tf/N
    @inbounds for i = 1:(N-1)
        X_eci[i+1] = rk4(cartesian_dynamics,X_eci[i],dt)
        X_eoe[i+1] = rk4(eoe_dynamics,X_eoe[i],dt)
    end

    X_eoe_eci = [rv_from_eoe(X_eoe[i]) for i = 1:length(X_eoe)]
    X_ks_eci = [rv_from_ks(X_ks[i][1:4],X_ks[i][5:8]) for i = 1:length(X_ks)]

    # eci_m = mat_from_vec(X_eci)
    # ks_m = mat_from_vec(X_ks_eci)

    # true energy
    # true_ϵ = -8.336482081923751e6
    true_ϵ = -8.3364820836748555e6
    # energies
    ks_ϵ = ϵ(X_ks_eci[end])
    eci_ϵ = ϵ(X_eci[end])
    eoe_ϵ = ϵ(X_eoe_eci[end])

    # @show ks_ϵ
    # @show eci_ϵ
    # @show eoe_ϵ
    # error()

    ks_error = -100*abs(ks_ϵ - true_ϵ)/true_ϵ
    eci_error = -100*abs(eci_ϵ- true_ϵ)/true_ϵ
    eoe_error = -100*abs(eoe_ϵ- true_ϵ)/true_ϵ
    # ks_error = abs(ks_ϵ - true_ϵ)
    # eci_error = abs(eci_ϵ- true_ϵ)
    # eoe_error = abs(eoe_ϵ- true_ϵ)
    # @show ks_ϵ
    # @show eci_ϵ
    # @show eoe_ϵ
    #
    # @show ks_error
    # @show eci_error
    # @show eoe_error
    #
    #
    # mat"
    # figure
    # hold on
    # %plot($eci_m(1,:),$eci_m(2,:),'o')
    # plot($eci_m(1,:),$eci_m(2,:),'o')
    # plot($ks_m(1,:),$ks_m(2,:))
    # hold off
    # "
    #
    # @show norm(X_eci[end] - X_ks_eci[end])
    # @show norm(X_eci[end] - X_eoe_eci[end])
    #
    # @show X_eci[end]
    return SA[ks_error, eci_error, eoe_error]
end

function grid_search()
    eci0, eoe0, ks0 = initial_conditions()

    N_vec = 15:1:300
    errors = [@SVector zeros(3) for i = 1:length(N_vec)]

    @showprogress "simulating" for i = 1:length(N_vec)
        errors[i] = driver(N_vec[i],eci0, eoe0, ks0 )
        # errors[i] = driver(10000000,eci0, eoe0, ks0 )
    end

    errors = mat_from_vec(errors)
    mat"
    figure
    hold on
    plot($errors')
    set(gca, 'YScale', 'log')
    set(gca, 'XScale', 'log')
    ylim([0 10000])
    xlim([15 54])
    legend('KS','ECI','EOE')
    hold off
    "
end
grid_search()

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
using Test, Unitful, NamedDims, AxisIndices, MappedArrays

@testset "no units, no time" begin
    A = NIMArray(reshape(1:12, 3, 4), x=1:3, y =1:4)
    @test @inferred(timeaxis(A)) === nothing
    @test !has_time_axis(A)
    @test timedim(A) == 0
    @test nimages(A) == 1
    @test @inferred(pixelspacing(A)) === (1,1)
    @test @inferred(spacedirections(A)) === ((1,0),(0,1))
    @test @inferred(coords_spatial(A)) === (1,2)
    @test @inferred(spatialorder(A)) === (:x, :y)
    @test @inferred(size_spatial(A)) === (3,4)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3), Base.OneTo(4))
    @test_throws ErrorException assert_timedim_last(A)
    #@test map(istimeaxis, AxisArrays.axes(A)) == (false,false)

end

@testset "units, no time" begin
    mm = u"mm"     # in real use these should be global consts
    m = u"m"
    A = NIMArray(reshape(1:12, 3, 4), x = 1mm:1mm:3mm, y = 1m:2m:7m)
    @test @inferred(timeaxis(A)) === nothing
    @test !has_time_axis(A)
    @test timedim(A) == 0
    @test nimages(A) == 1
    @test @inferred(pixelspacing(A)) === (1mm,2m)
    @test spacedirections(A) === ((1mm,0m),(0mm,2m))   # TODO: make this inferrable
    @test @inferred(coords_spatial(A)) === (1,2)
    @test spatialorder(A) === (:x,:y)
    @test @inferred(size_spatial(A)) === (3,4)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),Base.OneTo(4))
    @test_throws ErrorException assert_timedim_last(A)
    #@test map(istimeaxis, AxisArrays.axes(A)) == (false,false)
end

@testset "units, time" begin
    s = u"s" # again, global const
    A = NIMArray(reshape(1:12, 3, 4), x = 1:3, time = 1s:1s:4s)
    @test @inferred(timeaxis(A)) === 1s:1s:4s
    @test has_time_axis(A)
    @test timedim(A) == 2
    @test nimages(A) == 4
    @test @inferred(pixelspacing(A)) === (1,)
    @test @inferred(spacedirections(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (1,)
    @test spatialorder(A) === (:x,)
    @test @inferred(size_spatial(A)) === (3,)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),)
    @test isnothing(assert_timedim_last(A))
    #@test map(istimeaxis, AxisArrays.axes(A)) == (false,true)
end

@testset "units, time first" begin
    s = u"s" # global const
    A = NIMArray(reshape(1:12, 4, 3), time = 1s:1s:4s, x = 1:3)
    @test @inferred(timeaxis(A)) === 1s:1s:4s
    @test has_time_axis(A)
    @test timedim(A) == 1
    @test nimages(A) == 4
    @test @inferred(pixelspacing(A)) === (1,)
    @test @inferred(spacedirections(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (2,)
    @test spatialorder(A) === (:x,)
    @test @inferred(size_spatial(A)) === (3,)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),)
    @test_throws ErrorException assert_timedim_last(A)
    #@test map(istimeaxis, AxisArrays.axes(A)) == (true,false)
end

# TODO
@testset "grayscale" begin
    A = NIMArray(rand(Gray{N0f8}, 4, 5), :y, :x)
    #@test summary(A) == "2-dimensional AxisArray{Gray{N0f8},2,...} with axes:\n    :y, Base.OneTo(4)\n    :x, Base.OneTo(5)\nAnd data, a 4×5 Array{Gray{N0f8},2} with eltype Gray{Normed{UInt8,8}}"
    cv = channel_view(A)
    @test axes(cv) == (1:4, 1:5)
    @test spatialorder(cv) == (:y, :x)
    @test colordim(cv) == 0
end

@testset "color" begin
    A = NIMArray(rand(RGB{N0f8}, 4, 5), :y, :x)
    cv = channel_iew(A)
    @test axes(cv) == (1:3, 1:4, 1:5)
    @test spatialorder(cv) == (:y, :x)
    @test colordim(cv) == 1
    p = permuteddimsview(cv, (2,3,1))
    @test dimnames(p) == (:y, :x, :color)
    @test axes(p) == (1:4, 1:5, 1:3)
    @test colordim(p) == 3
    @test has_color_axis(p)
    @test indices_color(p) == 1:3
end

@testset "frequency" begin
    a = rand(2,2)
    nia = NIArray(a, x = 2:3, frequency = 3:4)
    @test has_frequency_axis(nia)
    @test frequencyaxis(nia) == 3:4
    @test indices_frequency(nia) == 1:2
    @test frequencydim(nia) == 2
end

@testset "channel" begin
    a = rand(2,2)
    nia = NIArray(a, x = 2:3, channel = 3:4)
    @test has_channel_axis(nia)
    @test channelaxis(nia) == 3:4
    @test indices_channel(nia) == 1:2
    @test channeldim(nia) == 2
end

@testset "observation" begin
    a = rand(2,2)
    nia = NIArray(a, x = 2:3, observation = 3:4)
    @test has_observation_axis(nia)
    @test observationaxis(nia) == 3:4
    @test indices_observation(nia) == 1:2
    @test observationdim(nia) == 2
end

@testset "nested" begin
    A = NIMArray(rand(RGB{N0f8}, 4, 5), y = range(1, step=2, length=4), x = 1:5)
    P = permuteddimsview(A, (2, 1))
    @test @inferred(pixelspacing(P)) == (1, 2)
    #M = mappedarray(identity, A)
    #@test @inferred(pixelspacing(M)) == (2, 1)
    s = u"s" # global const
    μm = u"μm" # global const
    A = NIMArray(rand(N0f16, 4, 5, 11),
                 y = range(1μm, step=2μm, length=4),
                 x = range(1μm, step = 1μm, length=5),
                 time = range(0.0s, step=0.1s, length=11))
    P = permuteddimsview(A, (3, 1, 2))
    M = mappedarray(identity, A)
    @test @inferred(pixelspacing(P)) == @inferred(pixelspacing(M)) == (2μm, 1μm)
    @test @inferred(timeaxis(P)) == @inferred(timeaxis(M)) == range(0.0s, step=0.1s, length=11)
    @test has_time_axis(P)
    @test coords_spatial(P) == (2, 3)
    @test coords_spatial(M) == (1, 2)
    @test spatialorder(P) == spatialorder(M) == (:y, :x)
    @test sampling_rate(P) == 1 / step(timeaxis(P))
    @test time_end(P) == last(timeaxis(P))
    @test onset(P) == first(timeaxis(P))
    @test @inferred(size_spatial(P)) == @inferred(size_spatial(M)) == (4, 5)
    @test_throws ErrorException assert_timedim_last(P)
    assert_timedim_last(M)

    A = NIMArray(rand(N0f16, 11, 5, 4),
                  time = range(0.0s, step=0.1s, length=11),
                  x = range(1μm, step = 1μm, length=5),
                  y = range(1μm, step=2μm, length=4))
    P = permuteddimsview(A, (3, 2, 1))
    M = mappedarray(identity, A)
    @test @inferred(pixelspacing(P)) == (2μm, 1μm)
    @test @inferred(pixelspacing(M)) == (1μm, 2μm)
    @test @inferred(timeaxis(P)) == @inferred(timeaxis(M)) == range(0.0s, step=0.1s, length=11)
    @test has_time_axis(P)
    @test coords_spatial(P) == (1, 2)
    @test coords_spatial(M) == (2, 3)
    @test spatialorder(P) == (:y, :x)
    @test spatialorder(M) == (:x, :y)
    @test @inferred(size_spatial(P)) == (4, 5)
    @test @inferred(size_spatial(M)) == (5, 4)
    @test assert_timedim_last(P) == nothing
    @test_throws ErrorException assert_timedim_last(M)
end

# Possibly-ambiguous functions
@testset "ambig" begin
    A = NIMArray(rand(RGB{N0f8},3,5), :x, :y)
    @test isa(convert(Array{RGB{N0f8},2}, A), Array{RGB{N0f8},2})
    @test isa(convert(Array{Gray{N0f8},2}, A), Array{Gray{N0f8},2})
end

@testset "internal" begin
    A = NIMArray(rand(RGB{N0f8},3,5), :x, :y)
    @test axes(A) isa Tuple{Vararg{<:AbstractAxis}}
end

@testset "ImageMetadata" begin
    a = rand(2,2)
    nia = NIArray(a, x = 2:3, y = 3:4)
    ima = IMArray(a, 2:3, 3:4)
    nima = NIMArray(a, x = 2:3, y = 3:4)
    ma = MetaArray(a)

    @test isempty(properties(ma))
    @test isempty(properties(ima))
    @test isempty(properties(nima))

    ma.foo = true
    ima.foo = true
    nima.foo = true

    @test ma.foo
    @test ima.foo
    @test nima.foo

    @test first(properties(ma)) == Pair(:foo, true)
    @test first(properties(ima)) == Pair(:foo, true)
    @test first(properties(nima)) == Pair(:foo, true)
end


text:
@safetestset "canonical_one_form" begin
    using PoincareInvariants: canonical_one_form
    using LinearAlgebra: dot
    using Random: rand

    @test_throws ArgumentError canonical_one_form([1, 2, 3], 0.5, nothing)
    @test_throws ArgumentError canonical_one_form([0.3], 3.4, nothing)

    for mid in [3, 9, 22]
        n = mid * 2

        q = rand(mid)
        p = rand(mid)
        z = [q; p]

        @inferred canonical_one_form(z, 0.1, nothing)

        θ = canonical_one_form(z, 1.3, nothing)

        @test θ[1:mid] == p
        @test θ[mid+1:end] == zeros(mid)

        @test dot(θ, z) ≈ dot(z, θ) ≈ dot(p, q)

        v = rand(n)
        @test dot(θ, v) ≈ dot(v, θ) ≈ dot(p, v[1:mid])
    end
end

@safetestset "canonical_two_form" begin
    using PoincareInvariants: CanonicalSymplecticMatrix, canonical_two_form
    using LinearAlgebra: dot
    using Random: rand

    @test CanonicalSymplecticMatrix{Int}(2) == [0 -1; 1 0]
    @test CanonicalSymplecticMatrix{Int}(4) == [0 0 -1 0; 0 0 0 -1; 1 0 0 0; 0 1 0 0]
    @test CanonicalSymplecticMatrix{Int}(6) == [0  0  0 -1  0  0;
                                                0  0  0  0 -1  0;
                                                0  0  0  0  0 -1;
                                                1  0  0  0  0  0;
                                                0  1  0  0  0  0;
                                                0  0  1  0  0  0]

    @test_throws ArgumentError CanonicalSymplecticMatrix{Float64}(-2)
    @test_throws ArgumentError CanonicalSymplecticMatrix{Int}(0)

    function test_canonicalmatrix(mat, n, ::Type{T}) where T
        mid = n ÷ 2

        for i in 1:n, j in 1:n
            if i == j - mid
                mat[i, j] === T(-1) || return false
            elseif i - mid == j
                mat[i, j] === T(1) || return false
            else
                mat[i, j] === T(0) || return false
            end
        end

        return true
    end

    @testset "CanonicalSymplecticMatrix{$T}($n)" for T in [Int, Float64], n in [10, 64, 1002]
        @test CanonicalSymplecticMatrix(n) === CanonicalSymplecticMatrix{Int}(n)

        C = CanonicalSymplecticMatrix{T}(n)
        @test size(C) == (n, n)
        @test eltype(C) == T

        mid = n ÷ 2

        @test test_canonicalmatrix(C, n, T)

        @test_throws BoundsError C[0, 4]
        @test_throws BoundsError C[8, -2]
        @test_throws BoundsError C[2    , n + 1]
        @test_throws BoundsError C[n + 1,     3]

        @test C * collect(1:n) == [(-mid-1:-1:-n)..., 1:mid...]

        @test dot(ones(n), C, ones(n)) == 0

        v = rand(n); w = rand(n)
        vCw = dot(v, C, w)

        @test vCw ≈ dot(v[mid+1:end], w[1:mid]) - dot(v[1:mid], w[mid+1:end])
        @test vCw ≈ dot(v, C * w)

        @test_throws ArgumentError CanonicalSymplecticMatrix{T}(n + 1)

        @inferred CanonicalSymplecticMatrix{Float64} canonical_two_form(v, 0.1, nothing)

        @test canonical_two_form(v, 0.1, nothing) == CanonicalSymplecticMatrix(n)
    end
end


text:
# atom
a0 = Atom(() -> true; N=0)
a1 = Atom(x -> x > 1)
a2 = Atom((x, y) -> x > y; N=2)

@test a0 == a0 && a0 != a1
@test evaluate(a0) && a0()
@test evaluate(a1, 4) && a1(4)
@test !evaluate(a2, 4, 5) && !a2(4, 5)
@test evaluate(a2, 6, 5) && a2(6, 5)

# negation
n0 = Negation(a0)
n1 = Negation(a1)
n2 = Negation(a2)

@test n0 == n0 && n0 != n1
@test !evaluate(n0) && !n0()
@test !evaluate(n1, 4) && !n1(4)
@test evaluate(n2, 4, 5) && n2(4, 5)
@test !evaluate(n2, 6, 5) && !n2(6, 5)

# conjunction
c0 = Conjunction([(n0, Int[]), (a0, Int[])]; N=0)
c0_3 = Conjunction([(a0, Int[]), (a0, Int[]), (a0, Int[])]; N=0)
c1 = Conjunction([(a0, Int[]), (a1, Int[1])]; N=1)
c1_3 = Conjunction([(a0, Int[]), (a1, Int[1]), (a1, Int[1])]; N=1)
c2 = Conjunction([(a2, Int[2, 1]), (a2, Int[1, 2])]; N=2)
c2_3 = Conjunction([(a0, Int[]), (a1, Int[2]), (a2, Int[2, 1])]; N=2)

@test c0 == c0 && c0 != c1
@test !evaluate(c0) && !c0()
@test evaluate(c0_3) && c0_3()
@test evaluate(c1, 4) && c1(4)
@test evaluate(c1_3, 4) && c1_3(4)
@test !evaluate(c2, 4, 5) && !c2(4, 5)
@test evaluate(c2_3, 4, 5) && c2_3(4, 5)

# default constructor
@test Conjunction([(a0, Int[]), (a0, Int[]), (a0, Int[])]) isa Conjunction{Val{0}}
@test Conjunction([(a0, Int[]), (a1, Int[2]), (a2, Int[2, 1])]) isa Conjunction{Val{2}}

# disjunction
d0 = Disjunction([(n0, Int[]), (a0, Int[])]; N=0)
d0_3 = Disjunction([(n0, Int[]), (n0, Int[]), (n0, Int[])]; N=0)
d1 = Disjunction([(c0, Int[]), (a1, Int[1])]; N=1)
d1_3 = Disjunction([(a0, Int[]), (a1, Int[1]), (a1, Int[1])]; N=1)
d2 = Disjunction([(a2, Int[2, 1]), (a2, Int[1, 2])]; N=2)
d2_3 = Disjunction([(a0, Int[]), (a1, Int[2]), (a2, Int[2, 1])]; N=2)

@test d0 == d0 && d0 != d1
@test evaluate(d0) && d0()
@test !evaluate(d0_3) && !d0_3()
@test evaluate(d1, 4) && d1(4)
@test evaluate(d1_3, 4) && d1_3(4)
@test evaluate(d2, 4, 5) && d2(4, 5)
@test evaluate(d2_3, 4, 5) && d2_3(4, 5)

# default constructor
@test Disjunction([(a0, Int[]), (a0, Int[]), (a0, Int[])]) isa Disjunction{Val{0}}
@test Disjunction([(a0, Int[]), (a1, Int[2]), (a2, Int[2, 1])]) isa Disjunction{Val{2}}

# unary conjunction
cu = Conjunction([n1, a1])
cu_3 = Conjunction([c1, a1, d1])

@test !evaluate(cu, 5) && !cu(5)
@test evaluate(cu_3, 5) && cu_3(5)

# unary disjunction
du = Disjunction([n1, a1])
du_1 = Disjunction([n1])

@test evaluate(du, 5) && du(5)
@test !evaluate(du_1, 5) && !du_1(5)


text:
"""
bloom_best variant

    * use "bloom first" (bloom_v4) for small patterns

    * use struct for initialitation results directlye

benchmark result is: seems like optimizer cannot remove struct access,
so runtime is worse.
"""
function bloom_best3 end

struct SearchStructure
    t::SearchSequence
    bloom_mask::UInt64
    bloom_skip::Int
    bloom_bits::Int
    skip::Int
    tlast::UInt8
end

function bloom_best3(t::SearchSequence)
    n = sizeof(t)
    skip = n
    tlast = _nthbyte(t,n)
    bloom_mask = UInt64(_search_bloom_mask(tlast))
    bloom_bits = 1 # no. of bits set in bloom filter if <= 32
    bloom_skip = 1 # no. of bytes to skip if byte not in filter
    j = n
    while (j-=1)>=1  # j=n is already processed
        tj = _nthbyte(t,j)
        if tj == tlast && skip == n
            skip = n - j
        end
        if bloom_bits <= 32 # argument: is near max(p(bloom_skip)*bloom_skip)
            hash = _search_bloom_mask(tj)
            if hash&bloom_mask == 0
                if bloom_bits < 32
                    # put in bloom filter up to 32 bits
                    bloom_mask |= hash
                end
                bloom_bits += 1 # gets >32 to terminate bloom filter filling
            end
            if bloom_bits <= 32
                bloom_skip += 1
            end
        else
            if skip<n
                # we have finished adding hashes to bloom filter
                # and we have determined the skip distance if matching last byte
                # nothimg remains to be done in preprocessing - stop work.
                break
            end
        end
    end
    skip -=1
    return SearchStructure(t,bloom_mask,bloom_skip,bloom_bits,skip,tlast)
end

function bloom_best3(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_best3(s,bloom_best3(t),i)
end

function bloom_best3(s::SearchSequence,p::SearchStructure,i::Integer,sv::MaybeVector=nothing)
    DOSTATS = !(sv isa Nothing)
    n = sizeof(p.t)
    m = sizeof(s)
    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(p.t,1)), s, i), 0)
    end
    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i +=n-1
    if p.bloom_bits <= 12 # best guess from benchmarks
        # do bloom test first in loop
        while i < m
            if DOSTATS loops += 1 end
            if DOSTATS bloomtests += 1 end
            if p.bloom_mask & _search_bloom_mask(_nthbyte(s,i)) == 0
                if DOSTATS bloomskips += 1 end
                i += p.bloom_skip
                continue
            end
            if _nthbyte(s,i) == p.tlast
                # check candidate
                j = 1
                while j < n
                    if _nthbyte(s,i-n+j) != _nthbyte(p.t,j)
                        break
                    end
                    j += 1
                    # match found?
                    if j == n
                        if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(p.bloom_mask), p.skip) end
                        return i-n+1
                    end
                end
                # no match: skip and test bloom
                i += p.skip
            end
            i +=1
        end
    else
        # do byte test first
        while i < m
            if DOSTATS loops += 1 end
            if _nthbyte(s,i) == p.tlast
                # check candidate
                j = 1
                while j < n
                    if _nthbyte(s,i-n+j) != _nthbyte(p.t,j)
                        break
                    end
                    j += 1
                    # match found?
                    if j == n
                        if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(p.bloom_mask), p.skip) end
                        return i-n+1
                    end
                end
                # no match: skip and test bloom
                i += p.skip
                if i>=w
                    break
                end
            end
            i += 1
            if DOSTATS bloomtests += 1 end
            if p.bloom_mask & _search_bloom_mask(_nthbyte(s,i)) == 0
                if DOSTATS bloomskips += 1 end
                i += p.bloom_skip
            end
        end
    end
    if i==m
        # test end match
        j = 1
        while j <= n
            if _nthbyte(s,i-n+j) != _nthbyte(p.t,j)
                break # not found
            end
            if j == n
                if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(p.bloom_mask), p.skip) end
                return i-n+1
            end # match at the very end
            j += 1
        end
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(p.bloom_mask), p.skip) end
    0
end


text:
module naumer_Dimensionality_2022

# Include the other files
# Include systems.jl first because of abstractSystem
include("systems.jl")
include("core.jl")

end


text:
using SortingAlgorithms
using RCall
using StatPlots, DataFrames
using SortingLab, BenchmarkTools
using CSV
using Plots

function string_sort_vs_r(N, K)
    samplespace = "id".*dec.(1:N÷K,10)
    srand(1);
    svec = rand(samplespace, N);
    fsortperm_timing = @belapsed ss = fsortperm($svec);
    # @assert issorted(svec[ss])

    radixsort_timing = @belapsed radixsort($svec);

    r_timing = R"""
    memory.limit(40000)
    ss = sprintf("id%010d",1:$(N/K))
    x = sample(ss, $N, replace=T)
    st = system.time(order(x, method="radix"))
    rm(x); rm(ss); gc()
    st
    """[3]

    r_timing_sort = R"""
    memory.limit(40000)
    ss = sprintf("id%010d",1:$(N/K))
    x = sample(ss, $N, replace=T)
    st = system.time(sort(x, method="radix"))
    rm(x); rm(ss); gc()
    st
    """[3]

    DataFrame(N = N, K = K, fsortperm_timing = fsortperm_timing,
    radixsort_timing=radixsort_timing, r_timing=r_timing,
    r_timing_sort=r_timing_sort)
end

#warm up
string_sort_vs_r(1000, 10)

# @time a = string_sort_vs_r(10_000_000, 1)

# @time  for k in [2, 4, 10, 20, 50, 100, 200, 500, 1000]
#     a = vcat(a, string_sort_vs_r(10_000_000, k))
#     print(a)
# end
# CSV.write("string_sort_vs_r.csv", a)

@time b = string_sort_vs_r(100_000_000, 1)

@time  for k in [2, 4, 10, 20, 50, 100, 200, 500, 1000]
    b = vcat(b, string_sort_vs_r(100_000_000, k))
    print(b)
end

CSV.write("string_sort_vs_r_100m.csv", b)

# load the data and plot them
a = CSV.read("string_sort_vs_r.csv")
plot(a[:N]./a[:K], a[:radixsort_timing],
title="Id Strings Vector (10m) sorting perf by # of uniques",
ylabel="seconds",
label="Julia",
xlabel="# of uniques")
plot!(a[:N]./a[:K], a[:r_timing_sort], label="R")
savefig("string_sort_vs_r_10m.png")

a = CSV.read("string_sort_vs_r_100m.csv")
plot(a[:N]./a[:K], a[:radixsort_timing],
title="Id Strings Vector (100m) sorting perf by # of uniques",
ylabel="seconds",
label="Julia",
xlabel="# of uniques")
plot!(a[:N]./a[:K], a[:r_timing_sort], label="R")
savefig("string_sort_vs_r_100m.png")

max_stars_repo_path: benchmarks/at what proportion of uniques does R win out.jl


text:
# This file contains helper functions for the MCC DAQ HAT Python examples.

#includet("/home/pi/daqhatsJulia/daqhats.jl")    # general daqhats commands
#includet("/home/pi/daqhatsJulia/daqhats172.jl") # MCC172 commands

    # type: (HatIDs) -> int
    """
    function select_hat_device(filter_by_id)
    This function performs a query of available DAQ HAT devices and determines
    the address of a single DAQ HAT device to be used in an example.  If a
    single HAT device is present, the address for that device is automatically
    selected, otherwise the user is prompted to select an address from a list
    of displayed devices.

    Args:
        filter_by_id (int): If this is :py:const:`HatIDs.ANY` return all DAQ
            HATs found.  Otherwise, return only DAQ HATs with ID matching this
            value.

    Returns:
        int: The address of the selected device.

    Raises:
        Exception: No HAT devices are found or an invalid address was selected.

    """
function select_hat_device(filter_by_id::Integer)
    selected_hat_address = Integer[]

    # Get descriptors for all of the available HAT devices.
    hats = hat_list(filter_by_id)
    number_of_hats = length(hats)

    # Verify at least one HAT device is detected.
    if number_of_hats < 1
        error("Error: No HAT devices found")
    elseif number_of_hats == 1
        selected_hat_address = hats[1].address
    else
        # Display available HAT devices for selection.
        for h = 1:number_of_hats
            println("Address $(hats[h].address): $(hats[h].product_name)")
        end
        println()

        println("Select the address of the HAT device to use: ")
        address = parse(Integer, (readline()))

        # Verify the selected address is valid.
        for h in 1:number_of_hats
            if address == hats[h].address
                selected_hat_address = address
                break
            end
        end
    end

    if isempty(selected_hat_address)
        error("Error: Invalid HAT selection")
    end

    return selected_hat_address
end

    # This one not required? but taken from Python
    """
    function enum_mask_to_string(enum_type, bit_mask):
    This function converts a mask of values defined by an IntEnum class to a
    comma separated string of names corresponding to the IntEnum names of the
    values included in a bit mask.

    Args:
        enum_type (Enum): The IntEnum class from which the mask was created.
        bit_mask (int): A bit mask of values defined by the enum_type class.

    Returns:
        str: A comma separated string of names corresponding to the IntEnum
        names of the values included in the mask

    """
function enum_mask_to_string(enum_type::Integer, bit_mask::Integer)
    item_names = String[]
    if bit_mask == 0
        push!(item_names, "DEFAULT")
    end
    for item in enum_type
        if item & bit_mask
            push!(item_names, item.name)
        end
    end

    return join(item_names, ", ")
end

"""
    function chan_list_to_mask(chan_list::Vector{T})
The vector contains the number for each channel that is activated

channel_list = [0, 1]  means both channel 0 and 1 are activated.
Sorry about the throwback to c style numbering.
"""
function chan_list_to_mask(chan_list::Vector{T}) where T <: Integer
    chan_mask = zero(UInt8)
    for chan in chan_list
        chan_mask |= 0x01 << chan
    end
    return chan_mask
end

"""
    function validate_channels(channel_set, number_of_channels)

    Raises a ValueError exception if a channel number in the set of
    channels is not in the range of available channels.

    Args:
        channel_set (set): A set of channel numbers.
        number_of_channels (int): The number of available channels.

    Returns:
        None

    Raises:
        ValueError: If there is an invalid channel specified.

"""
function validate_channels(channel_set::Vector{<:Integer}, number_of_channels::Integer)
    # not sure if I have this right
    valid_chans = 0:number_of_channels - 1
    if !issubset(channel_set, valid_chans)
        error("Error: Invalid channel selected - must be $(min(valid_chans)) - $(max(valid_chans))")
    end
end

"""
    function deinterleave(data::AbstractVector, nc::Integer)
Convert interleaved vector to Matrix of nc columns
"""
deinterleave(data::AbstractVector, nc::Integer) = transpose(reshape(data, Int(nc), :))

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
function get_buy_value(config, trade_grid)
    cQty  = typeof(config["init_state"]["cummulativeQuoteQty"]) != Float64 ? parse(Float64, config["init_state"]["cummulativeQuoteQty"]) : config["init_state"]["cummulativeQuoteQty"]
    tQty  = typeof(config["total_value"]) != Float64 ? parse(Float64, config["total_value"]) : config["total_value"]

    total_range = typeof(config["grid"]["total_range"]) != Float64 ? parse(Float64, config["grid"]["total_range"]) : config["grid"]["total_range"]
    grid_step   = typeof(config["grid"]["step"]) != Float64 ? parse(Float64, config["grid"]["step"]) : config["grid"]["step"]

    total_steps =  total_range / grid_step
    sell_step_num = floor(total_steps * (cQty / tQty))
    buy_step_num  = floor(total_steps * (1 - cQty / tQty))

    current_volume_diff = Pinguu.get_current_volume(config) - tQty

    buy_coins  = tQty - cQty
    buy_step_val  = round((buy_coins  / buy_step_num) + (current_volume_diff/length(trade_grid)), digits=config["precision"])
    return buy_step_val
end
export get_buy_value


text:
using Mamba, Plots
pyplot(size=(800,800))

if !isdefined(Main, :traceplot)
  function traceplot(c::AbstractChains; legend::Bool=false, na...)
    nrows, nvars, nchains = size(c.value)
    plts = Array{Plots.Plot{Plots.PyPlotBackend}}(nvars)
    for i in 1:nvars
      plts[i] = plot(1:nrows, c.value[:, i, :])
    end
    plts
  end
end

if !isdefined(Main, :meanplot)
  function meanplot(c::AbstractChains; legend::Bool=false, na...)
    nrows, nvars, nchains = size(c.value)
    plts = Array{Plots.Plot{Plots.PyPlotBackend}}(nvars)
    val = Mamba.cummean(c.value)
    for i in 1:nvars
      plts[i] = plot(1:nrows, val[:, i, :], lab="Chain $i")
    end
    plts
  end
end

tp = traceplot(sim)

mp = meanplot(sim)


text:
"Radial Basis Function Kernel also called RBF or SE(Squared Exponential)"
struct RBFKernel{T<:Real,KT<:KernelType} <: Kernel{T,KT}
    fields::KernelFields{T,KT}
    function RBFKernel{T,KT}(θ::Vector{T};variance=one(T),dim::Integer=0) where {T<:Real,KT<:KernelType}
        if KT == ARDKernel
            if length(θ)==1 && dim ==0
                error("You defined an ARD RBF kernel without precising the number of dimensions or giving a vector for the lengthscale                   Please set dim in your kernel initialization")
            elseif dim!=0 && (length(θ)!=dim && length(θ)!=1)
                @warn "You did not use the same dimension for your params and dim, using the first value of params for all dimensions"
                θ = ones(dim,T=T)*θ[1]
            elseif length(θ)==1 && dim!=0
                θ = ones(dim)*θ[1]
            end
            dim = length(θ)
            return new{T,ARDKernel}(KernelFields{T,ARDKernel}(
                                        "Radial Basis",
                                        HyperParameter{T}(variance,interval(OpenBound(zero(T)),nothing),fixed=false),
                                        HyperParameters{T}(θ,[interval(OpenBound(zero(T)),NullBound{T}()) for _ in 1:dim]),dim,
                                        WeightedSqEuclidean(one(T)./(θ.^2))))
        else
            return new{T,IsoKernel}(KernelFields{T,IsoKernel}(
                                        "Radial Basis",
                                        HyperParameter{T}(variance,interval(OpenBound(zero(T)),nothing),fixed=false),
                                        HyperParameters{T}(θ,[interval(OpenBound(zero(T)),NullBound{T}())]),1,
                                        SqEuclidean()))
        end
    end
end

const SEKernel = RBFKernel

"""`RBFKernel(θ::T1=1.0;variance::T2=one(T1),dim::Integer=0,ARD::Bool=false)`

Create a Squared Exponential/Radial Basis Function kernel

``variance\\exp\\left(-\\frac{||x-x'||}{2\\theta^2}\\right)``

The kernel can be made with ARD (Automatic Relevance Determination) by giving `θ` as a vector or setting `ARD` as true, and in both cases setting the `dim` variable appropriately
"""
function RBFKernel(θ::T1=1.0;variance::T2=one(T1),dim::Integer=0,ARD::Bool=false) where {T1<:Real,T2<:Real}
    if ARD
        RBFKernel{T1,ARDKernel}([θ],variance=variance,dim=dim)
    else
        RBFKernel{T1,IsoKernel}([θ],variance=variance)
    end
 end

function RBFKernel(θ::Array{T1,1};variance::T2=one(T1),dim::Integer=0) where {T1<:Real,T2<:Real}
    RBFKernel{T1,ARDKernel}(θ,variance=variance,dim=dim)
end

@inline rbfkernel(z::Real, l::Real) = exp(-0.5*z/(l^2))

@inline rbfkernel(z::Real) = exp(-0.5*z)

function kappa(k::RBFKernel{T,IsoKernel}) where {T<:Real}
    return z->rbfkernel(z,getlengthscales(k))
end

function kappa(k::RBFKernel{T,ARDKernel}) where {T<:Real}
    return z->rbfkernel(z)
end

function updateweights!(k::RBFKernel{T,KT},w::Vector{T}) where {T,KT}
    k.fields.metric.weights .= 1.0./(w.^2)
end

function computeIndPointsJmm(k::RBFKernel{T,KT},X::Matrix{T},iPoint::Integer,K::Symmetric{T,Matrix{T}}) where {T,KT}
    l2 = (getlengthscales(k)).^2; v = getvariance(k)
    return -v.*((X[iPoint,:]'.-X)./l2').*K[:,iPoint]
end

function computeIndPointsJnm(k::RBFKernel{T,KT},X::Matrix{T},x::Vector{T},iPoint::Integer,K::Matrix{T}) where {T,KT}
    l2 = (getlengthscales(k)).^2; v = getvariance(k)
    return -v.*((x'.-X)./l2').*K[:,iPoint]
end

################# Matrix derivatives for the RBF kernel###################################
"Return the kernel matrix derivative for the Iso RBFKernel"
function kernelderivativematrix(X::AbstractArray{T},kernel::RBFKernel{T,IsoKernel}) where {T<:Real}
    v = getvariance(kernel); l = getlengthscales(kernel)
    P = pairwise(SqEuclidean(),X,dims=1); K = zero(P)
    map!(kappa(kernel),K,P)
    return Symmetric(lmul!(v./(l^3),P.*=K))
end

"""Return the matrix derivative for the Iso RBFKernel with the covariance matrix precomputed"""
function kernelderivativematrix_K(X::AbstractArray{T},K::Symmetric{T,AbstractArray{T,2}},kernel::RBFKernel{T,IsoKernel}) where {T<:Real}
    v = getvariance(kernel); l = getlengthscales(kernel)
    P = pairwise(SqEuclidean(),X,dims=1)
    return Symmetric(lmul!(v./(l^3),P.*=K))
end

"Return the matrix derivative for the ARD RBFKernel"
function kernelderivativematrix(X::AbstractArray{T},kernel::RBFKernel{T,ARDKernel}) where {T<:Real}
    v = getvariance(kernel); ls = getlengthscales(kernel)
    K = pairwise(getmetric(kernel),X,dims=1)
    Pi = [pairwise(SqEuclidean(), X[:,i]',dims=2) for i in 1:length(ls)]
    map!(kappa(kernel),K,K)
    return Symmetric.(map!((pi,l)->lmul!(v./(l^3),pi.*=K),Pi,Pi,ls))
end

"""Return the matrix derivative for the Iso RBFKernel with the covariance matrix precomputed"""
function kernelderivativematrix_K(X::AbstractArray{T},K::Symmetric{T,AbstractArray{T,2}},kernel::RBFKernel{T,ARDKernel}) where {T<:Real}
    v = getvariance(kernel); ls = getlengthscales(kernel)
    Pi = [pairwise(SqEuclidean(),X[:,i],dims=1) for i in 1:length(ls)]
    return Symmetric.(map((pi,l)->lmul!(v./(l^3),pi.*=K),Pi,ls))
end

########## DERIVATIVE MATRICES FOR TWO MATRICES #######

function kernelderivativematrix(X::AbstractArray{T},Y::AbstractArray{T},kernel::RBFKernel{T,IsoKernel}) where {T<:Real}
    v = getvariance(kernel); l = getlengthscales(kernel)
    P = pairwise(SqEuclidean(),X,Y,dims=1); K = zero(P);
    map!(kappa(kernel),K,P)
    return lmul!(v./(l^3),P.*=K)
end

"When K has already been computed"
function kernelderivativematrix_K(X::AbstractArray{T},Y::AbstractArray{T},K::AbstractArray{T},kernel::RBFKernel{T,IsoKernel}) where {T}
    v = getvariance(kernel); l = getlengthscales(kernel)
    P = pairwise(SqEuclidean(),X,Y,dims=1)
    return lmul!(v./(l^3),P.*=K)
end

function kernelderivativematrix(X::AbstractArray{T},Y::AbstractArray{T},kernel::RBFKernel{T,ARDKernel}) where {T}
    v = getvariance(kernel); ls = getlengthscales(kernel)
    K = pairwise(getmetric(kernel),X,Y,dims=1)
    Pi = [pairwise(SqEuclidean(),X[:,i]',Y[:,i]',dims=2) for i in 1:length(ls)]
    map!(kappa(kernel),K,K)
    return map((pi,l)->lmul!(v./(l^3),pi.*=K),Pi,ls)
end

"When K has already been computed"
function kernelderivativematrix_K(X::AbstractArray{T},Y::AbstractArray{T},K::AbstractArray{T},kernel::RBFKernel{T,ARDKernel}) where {T}
    v = getvariance(kernel); ls = getlengthscales(kernel)
    Pi = [pairwise(SqEuclidean(),X[:,i],Y[:,i],dims=1) for i in 1:length(ls)]
    return map((pi,l)->lmul!(v./(l^3),pi.*=K),Pi,ls)
end

############ DIAGONAL DERIVATIVES ###################

function kernelderivativediagmatrix(X::AbstractArray{T},kernel::RBFKernel{T,IsoKernel}) where {T<:Real}
    return zeros(T,size(X,1))
end

function kernelderivativediagmatrix(X::AbstractArray{T},kernel::RBFKernel{T,ARDKernel}) where {T}
    return [zeros(T,size(X,1)) for _ in 1:kernel.fields.Ndim]
end


text:
struct Enc{T<:Group}
    pk::T
    g::T
end

(enc::Enc{T})(m::T, r::Integer) where T <: Group = (enc.g^r, m*enc.pk^r) ### Message first?
(enc::Enc)(r::Integer) = (enc.g^r, enc.pk^r)

a(x::Tuple{T, T}) where T <: Group = x[1]
b(x::Tuple{T, T}) where T <: Group = x[2]

*(x::Tuple{G, G}, y::Tuple{G, G}) where G <: Group = (a(x)*a(y), b(x)*b(y))

(enc::Enc)(e::Tuple{G, G}, r::Integer) where G <: Group = e * enc(r)

struct ElGamal{G <: Group} <: AbstractVector{G}
    a::Vector{G}
    b::Vector{G}

    function ElGamal{G}(a::Vector{G}, b::Vector{G}) where {G <: Group}
        @assert length(a) == length(b)
        return new(a, b)
    end
end

ElGamal(a::Vector{G}, b::Vector{G}) where G <: Group = ElGamal{G}(a, b)

ElGamal(e::Vector{Tuple{G, G}}) where G <: Group = ElGamal([a(i) for i in e], [b(i) for i in e])

function ElGamal{G}(a::Vector{T}, b::Vector{T}) where {T, G<:Group}
    a′ = convert(Vector{G}, a)
    b′ = convert(Vector{G}, b)

    return ElGamal{G}(a′, b′)
end

a(e::ElGamal) = e.a
b(e::ElGamal) = e.b

Base.getindex(e::ElGamal, i::Int) = (a(e)[i], b(e)[i])
Base.getindex(e::ElGamal{G}, ivec::Vector) where G <: Group = ElGamal{G}(a(e)[ivec], b(e)[ivec])
Base.length(e::ElGamal) = length(a(e))
Base.size(e::ElGamal) = size(a(e))

Base.copy(e::ElGamal{G}) where G <: Group = ElGamal{G}(copy(a(e)), copy(b(e)))
Base.sort(e::ElGamal) = sort!(copy(e))

### NEW METHOD
function Base.setindex!(e::ElGamal{G}, val::Tuple{G, G}, i) where G <: Group
    a(e)[i] = a(val)
    b(e)[i] = b(val)
end

function *(x::ElGamal{G}, y::ElGamal{G}) where G <: Group

    @assert length(x) == length(y)

    a′ = a(x) .* a(y)
    b′ = b(x) .* b(y)

    return ElGamal(a′, b′)
end

function *(x::ElGamal{G}, y::Tuple{G, G}) where G <: Group

    a′ = a(x) .* a(y)
    b′ = b(x) .* b(y)

    return ElGamal(a′, b′)
end

*(x::Tuple{G, G}, y::ElGamal{G}) where G <: Group = y * x

import Base: ^

^(x::Tuple{G, G}, k::Integer) where G <: Group = (x[1]^k, x[2]^k)
^(x::ElGamal, k::AbstractVector{<:Integer}) = ElGamal(a(x) .^ k, b(x) .^ k)
^(x::ElGamal, k::Integer) = ElGamal(a(x) .^ k, b(x) .^ k)

Base.broadcasted(::typeof(^), x::ElGamal, y::AbstractVector{<:Integer}) = x^y

(enc::Enc)(e::ElGamal, r::Integer) = enc(r) * e

function (enc::Enc{G})(m::Vector{G}, r::AbstractVector{<:Integer}) where G <: Group

    a′ = enc.g .^ r
    b′ = m .* (enc.pk .^ r)

    return ElGamal(a′, b′)
end

function (enc::Enc{G})(𝐞::ElGamal{G}, 𝐫::AbstractVector{<:Integer}) where G <: Group

    ### Need to improve this
    r = ElGamal(enc.(𝐫))
    𝐞′ = 𝐞 * r

    return 𝐞′
end

struct Dec
    sk::Integer
end

(dec::Dec)(e::Tuple{G, G}) where G = b(e) * a(e)^(-dec.sk)

(dec::Dec)(e::ElGamal) = [dec(ei) for ei in e]

Base.isless(x::Tuple{G, G}, y::Tuple{G, G}) where G <: Group = x[1] == y[1] ? x[2] < y[2] : x[1] < y[1]

Base.prod(e::ElGamal{G}) where G <: Group = (prod(e.a), prod(e.b))


text:
@testset "Float32" begin
    for n in [1,3,10]

A = Array{Float32}(undef, n, n)
copyto!(A, rand(n,n))

r1 = polarfact(A, alg =:newton);

r2 = polarfact(A, alg =:qdwh);

r3 = polarfact(A, alg =:halley);

r4 = polarfact(A, alg =:svd);

r5 = polarfact(A, alg =:hybrid);

@test r1.U ≈ r2.U
@test r2.U ≈ r3.U
@test r3.U ≈ r4.U
@test r4.U ≈ r5.U

    end
end


text:
include("testPlotGrid.jl")
include("testViewOrthoSlices2D.jl")
include("testViewImage2D.jl")
include("testViewSlice2D.jl")


text:
@doc raw"""
    zernike(x, y, maxord)

Calculate  unit-normalized Zernike polynomials and their x,y derivatives

Translation of the code from
[1] T. B. Andersen, “Efficient and robust recurrence relations for the Zernike circle polynomials
and their derivatives in Cartesian coordinates,” _Opt. Express_, vol. 26, no. 15, p. 18878, Jul. 2018.

     Numbering scheme:
     Within a radial order, sine terms come first
             ...
           sin((n-2m)*theta)   for m = 0,..., [(n+1)/2]-1
             ...
              1                for n even, m = n/2
             ...
           cos((n-2m)*theta)   for m = [n/2]+1,...,n
             ...

     INPUT:
     x, y normalized (x,y) coordinates in unit circle
     MaxOrd: Maximum Zernike radial order

     OUTPUT:
     dUdx[...]   array to receive each derivative dU/dx at (x,y)
     dUdy[...]   array to receive each derivative dU/dy at (x,y)

"""
function zernike(x, y, maxord::Int)

    ztotnumber = Int((maxord + 2) * (maxord + 1) / 2)
    # print(ztotnumber)
    zern = zeros(ztotnumber)
    dudx = zeros(ztotnumber)
    dudy = zeros(ztotnumber)

    zern[1] = 1
    dudx[1] = 0
    dudy[1] = 0

    zern[2] = y
    zern[3] = x
    dudx[2] = 0
    dudx[3] = 1
    dudy[2] = 1
    dudy[3] = 0

    kndx = 1                # index for term from 2 orders down
    jbeg = 2                # start index for current radial order
    jend = 3                # end index for current radial order
    jndx = 3                # running index for current Zern
    even = -1

    #  Outer loop in radial order index
    for nn = 2:maxord

        even   = -even          # parity of radial index
        jndx1  = jbeg           # index for 1st ascending series in x
        jndx2  = jend           # index for 1st descending series in y
        jndx11 = jndx1 - 1      # index for 2nd ascending series in x
        jndx21 = jndx2 + 1      # index for 2nd descending series in y
        jbeg   = jend + 1       # end of previous radial order +1
        nn2    = nn / 2
        nn1    = (nn - 1) / 2

      # Inner loop in azimuthal index
        for mm = 0:nn
            jndx += 1                  # increment running index for current Zern

            if mm == 0
                zern[jndx] = x * zern[jndx1] + y * zern[jndx2]
                dudx[jndx] = zern[jndx1] * nn
                dudy[jndx] = zern[jndx2] * nn

            elseif mm == nn
                zern[jndx] = x * zern[jndx11] - y * zern[jndx21]
                dudx[jndx] = zern[jndx11] * nn
                dudy[jndx] = -zern[jndx21] * nn

            elseif even > 0 && mm == nn2
                zern[jndx] = 2 * (x * zern[jndx1] + y * zern[jndx2]) - zern[kndx]
                dudx[jndx] = 2 * nn * zern[jndx1] + dudx[kndx]
                dudy[jndx] = 2 * nn * zern[jndx2] + dudy[kndx]
                kndx += 1

            elseif even < 0 && mm == nn1
                qval = zern[jndx2] - zern[jndx21]
                zern[jndx] = x * zern[jndx11] + y * qval - zern[kndx]
                dudx[jndx] = zern[jndx11] * nn + dudx[kndx]
                dudy[jndx] = qval * nn + dudy[kndx]
                kndx += 1

            elseif even < 0 && mm == nn1 + 1
                pval = zern[jndx1] + zern[jndx11]
                zern[jndx] = x * pval + y * zern[jndx2] - zern[kndx]
                dudx[jndx] = pval * nn + dudx[kndx]
                dudy[jndx] = zern[jndx2] * nn + dudy[kndx]
                kndx += 1

            else
                pval = zern[jndx1] + zern[jndx11]
                qval = zern[jndx2] - zern[jndx21]
                zern[jndx] = x * pval + y * qval - zern[kndx]
                dudx[jndx] = pval * nn + dudx[kndx]
                dudy[jndx] = qval * nn + dudy[kndx]
                kndx += 1

            end

            jndx11 = jndx1                   # update indices
            jndx1 += 1
            jndx21 = jndx2
            jndx2 += -1
        end

        jend = jndx
    end

    (z = zern, zx = dudx, zy = dudy)

end

"""
    ZernikeBW(elements, ap, mask, norms)

Contains Zernike basis (in Born&Wolf norming) with the aperture, plotting mask, and element norms. Can be constructed as
    `ZernikeBW(dom::CartesianDomain2D, d::Real, maxorder::Integer) `
and
    `ZernikeBW(gridsize::Integer, maxorder::Integer)`

"""
struct ZernikeBW <: OrthogonalBasis
    elements::VectorOfArray
    ap::Array
    mask::Array
    norms::Vector
    ZernikeBW(elements, ap, mask) = new(elements, ap, mask, [sqrt.(inner(f, ap .* f)) for f in elements])
end

ZernikeBW(gridsize::Integer, maxorder::Integer)  = ZernikeBW(makezerniketable(gridsize, maxorder),makeaperture(gridsize)...)
ZernikeBW(dom::CartesianDomain2D, d::Real, maxorder::Integer)  = ZernikeBW(makezerniketable(dom, maxorder,d/2.),aperture(dom, d)...)

# function makezerniketable(gridsize::Integer, maxorder::Integer)
#     x = range(-1, 1, length=gridsize)
#     y = range(-1, 1, length=gridsize)
#     totalznum = Int((maxorder + 2) * (maxorder + 1) / 2)
#     ztable = [zernike(xc, yc, maxorder)[:z] for xc ∈ x,  yc ∈ y ]
#     zvec = VectorOfArray([zeros(gridsize, gridsize) for i = 1:totalznum])
#     [zvec[i,:] = ztable[i] for i = eachindex(ztable) ]
#     return zvec
# end

function makezerniketable(dom::CartesianDomain2D,  maxorder::Integer, scale = 1)
    x= dom.xrange/scale
    y = dom.yrange/scale
    totalznum = Int((maxorder + 2) * (maxorder + 1) / 2)
    ztable = [zernike(xc, yc, maxorder)[:z] for yc ∈ y, xc ∈ x ]
    zvec = VectorOfArray([zeros(length(y), length(x)) for i = 1:totalznum])
    [zvec[i,:] = ztable[i] for i = eachindex(ztable) ]
    return zvec
end

makezerniketable(gridsize::Integer, maxorder::Integer) =
    makezerniketable(CartesianDomain2D(range(-1, 1, length=gridsize), range(-1, 1, length=gridsize)), maxorder)

function makeaperture(gridsize::Integer, δ = 0.)
    x = range(-1, 1, length=gridsize)
    y = range(-1, 1, length=gridsize)
    # δ = 0. # tuning of the aperture size
    r =1 + δ  /gridsize
    ap = [ (xc^2 + yc^2) <= r^2 ? 1 : 0 for xc ∈ x, yc ∈ y]
    # area = +(ap[:]...)
    phmask = [ (xc^2 + yc^2) <= r^2 ? 1 : NaN for xc ∈ x, yc ∈ y]
    return(ap, phmask)
end