text:
#================================================================

Definition of SignALSH, an LSH function for hashing on inner products.

================================================================#

import LinearAlgebra: norm

#========================
Typedefs
========================#

mutable struct SignALSH{T <: Union{Float32,Float64}} <: AsymmetricLSHFunction
    coeff_A :: Matrix{T}
    coeff_B :: Matrix{T}
    P_shift :: Vector{T}
    m :: Int64

    # An upper bound on the norm of the data points this hash function will
    # process
    maxnorm :: T

    # Whether or not SignALSH should round up to the next power of 2 when
    # resizing its coefficient array.
    resize_pow2 :: Bool
end

### External SignALSH constructors

@doc """
    SignALSH(n_hashes::Integer = $(DEFAULT_N_HASHES),
             dtype::Type = $(DEFAULT_DTYPE),
             maxnorm::Union{Nothing,Real} = nothing,
             m::Integer = 3,
             resize_pow2::Bool = $(DEFAULT_RESIZE_POW2))

Create a `SignALSH` hash function for hashing on inner product similarity.

# Arguments
- $(N_HASHES_DOCSTR())

# Keyword parameters
- $(DTYPE_DOCSTR(SignALSH))
- `maxnorm::Union{Nothing,Real}` (default: `nothing`): an upper bound on the ``\\ell^2``-norm of the data points.

!!! warning "Warning: maxnorm must be set"
    The `maxnorm` keyword parameter must be explicitly specified. If it is left unspecified (or set to `nothing`), `SignALSH()` will raise an error.

- `m::Integer` (default: `3`): parameter `m` that affects the probability of a hash collision.
- $(RESIZE_POW2_DOCSTR(SignALSH))

# Examples
`SignALSH` is an [`AsymmetricLSHFunction`](@ref), and hence hashes must be computed using [`index_hash`](@ref) and [`query_hash`](@ref).

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SignALSH(12; maxnorm=10);

julia> x = rand(4);

julia> ih = index_hash(hashfn, x); qh = query_hash(hashfn, x);

julia> length(ih) == length(qh) == 12
true

julia> typeof(ih) == typeof(qh) == BitArray{1}
true
```

You need to explicitly specify the `maxnorm` keyword parameter when constructing `SignALSH`, otherwise you will get an error.

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SignALSH(12)
ERROR: maxnorm must be specified for SignALSH
```

You'll also get an error if you try to hash a vector that has norm greater than the `maxnorm` that you specified.

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SignALSH(; maxnorm=1);

julia> index_hash(hashfn, ones(4))
ERROR: norm 2.0 exceeds maxnorm (1.0)
```

# References
- Anshumali Shrivastava and Ping Li. *Improved Asymmetric Locality Sensitive Hashing (ALSH) for Maximum Inner Product Search (MIPS)*. In Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence, UAI'15, page 812–821, Arlington, Virginia, USA, 2015. AUAI Press. 10.5555/3020847.3020931. [arXiv:1405.5869](https://arxiv.org/abs/1410.5410v2)

See also: [`inner_prod`](@ref inner_prod(::AbstractVector, ::AbstractVector)), [`ℓ2_norm`](@ref ℓp_norm)
"""
@generated function SignALSH{T}(n_hashes::Integer = DEFAULT_N_HASHES;
                                maxnorm::Union{Nothing,Real} = nothing,
                                m::Integer = 3,
                                resize_pow2::Bool = DEFAULT_RESIZE_POW2) where T

    if maxnorm <: Nothing
        :("maxnorm must be specified for SignALSH" |> ErrorException |> throw)
    else
        quote
            if maxnorm < 0
                "maxnorm must be non-negative" |> ErrorException |> throw
            elseif n_hashes ≤ 0
                "n_hashes must be positive" |> ErrorException |> throw
            elseif m ≤ 0
                "m must be positive" |> ErrorException |> throw
            end

            coeff_A = Matrix{T}(undef, n_hashes, 0)
            coeff_B = randn(T, n_hashes, m)
            P_shift = coeff_B * fill(T(1/2), m)

            SignALSH(coeff_A, coeff_B, P_shift, Int64(m),
                     T(maxnorm), resize_pow2)
        end
    end
end

SignALSH(args...; dtype=DEFAULT_DTYPE, kws...) =
    SignALSH{dtype}(args...; kws...)

#============
SignALSH helper functions
=============#

function Base.resize!(hashfn::SignALSH{T}, n::Integer) where T
    n = (hashfn.resize_pow2) ? nextpow(2, n) : n

    # Note: the only field of SignALSH that is dependent on the input size is
    # coeff_A
    n_hashes, old_n = size(hashfn.coeff_A)
    old_coeff_A = hashfn.coeff_A
    new_coeff_A = similar(old_coeff_A, n_hashes, n)

    new_coeff_A[1:end,1:min(n,old_n)] .= old_coeff_A

    if n > old_n
        new_coeff_slice = @views new_coeff_A[1:end,old_n+1:end]
        @views map!(x -> randn(T), new_coeff_slice, new_coeff_slice)
    end

    hashfn.coeff_A = new_coeff_A
end

current_max_input_size(hashfn::SignALSH) =
    size(hashfn.coeff_A, 2)

#============
Hash computation implementation
=============#

#=
h(P(x)) definitions
=#
SignALSH_P(hashfn::SignALSH{T}, x::AbstractArray) where T =
    SignALSH_P(hashfn, T.(x))

function SignALSH_P(hashfn::SignALSH{T}, x::AbstractArray{T}) where {T}
    # SignALSH_P is essentially just SimHash on
    #
    #	P(x) = [x; 1/2-norms^2; 1/2-norms^4; ... 1/2-norms^(2^m)]
    #
    # after dividing through x by the largest norm of all of the columns
    # of x.
    norms = col_norms(x)

    for norm_ii in norms
        if norm_ii > hashfn.maxnorm
            "norm $(norm_ii) exceeds maxnorm ($(hashfn.maxnorm))" |>
            ErrorException |>
            throw
        end
    end

    norms .*= 1/hashfn.maxnorm

    n = size(x,1)
    if n > current_max_input_size(hashfn)
        resize!(hashfn, n)
    end

    Ax = @views hashfn.coeff_A[1:end,1:n] * x .* (1/hashfn.maxnorm)

    # Perform the transformation P(x) on x, except that instead of actually
    # allocating memory for it and computing it, pile it onto Ax
    @. Ax += hashfn.P_shift

    for ii = 1:hashfn.m
        # Iteratively compute norms^2, norms^4, ... norms^(2^m). Multiply them by
        # columns of coeff_B and add to aTx, so that we end up with
        #
        # 	Ax = [coeff_A, coeff_B] * P(x)
        # 		= [coeff_A, coeff_B] * [x; 1/2-norms^2; ...; 1/2-norms^(2^m)]
        #
        # Note that we don't need to account for the 1/2 terms, which are
        # accounted for by the earlier addition of P_shift.
        norms .^= 2
        SignALSH_P_update_Ax!(hashfn.coeff_B[:,ii], norms, Ax)
    end

    return Ax .≥ 0
end

SignALSH_P_update_Ax!(coeff::Vector{T}, norms::Vector{T}, Ax::Matrix{T}) where T =
	BLAS.ger!(T(-1), coeff, norms, Ax)

function SignALSH_P_update_Ax!(coeff::AbstractVector,
                               norms::AbstractVector,
                               Ax::AbstractVector)

    Ax .-= coeff * norms' |> vec
end

# When the coefficients or norms are AbstractVectors, cascade through to reshape
# them into AbstractMatrix before updating Ax.
SignALSH_P_update_Ax!(coeff::AbstractVector, norms::AbstractVector, Ax) =
    SignALSH_P_update_Ax!(reshape(coeff, length(coeff), 1), norms, Ax)

SignALSH_P_update_Ax!(coeff, norms::AbstractVector, Ax) =
    SignALSH_P_update_Ax!(coeff, reshape(norms, length(norms), 1), Ax)

SignALSH_P_update_Ax!(coeff, norms, Ax) =
    (Ax .-= coeff * norms')

#=
h(Q(x)) definitions
=#

SignALSH_Q(hashfn::SignALSH{T}, x::AbstractArray) where T =
    SignALSH_Q(hashfn, T.(x))

SignALSH_Q(hashfn::SignALSH{T}, x::AbstractVector{T}) where T =
    invoke(SignALSH_Q, Tuple{SignALSH{T},AbstractArray{T}}, hashfn, x) |> vec

function SignALSH_Q(hashfn::SignALSH{T}, x::AbstractArray{T}) where {T}
    n = size(x,1)

    if n > current_max_input_size(hashfn)
        resize!(hashfn, n)
    end

    Ax = @views hashfn.coeff_A[1:end,1:n] * x
    norms = col_norms(x)

    for norm_ii in norms
        if norm_ii > hashfn.maxnorm
            "norm $(norm_ii) exceeds maxnorm ($(hashfn.maxnorm))" |>
            ErrorException |>
            throw
        end
    end

    map!(inv, norms, norms)
    @. Ax * norms' ≥ T(0)
end

#========================
LSHFunction and AsymmetricLSHFunction API compliance
========================#
index_hash(h::SignALSH, x) = SignALSH_P(h, x)
query_hash(h::SignALSH, x) = SignALSH_Q(h, x)

n_hashes(h::SignALSH) = size(h.coeff_A, 1)
similarity(::SignALSH) = inner_prod
hashtype(::SignALSH) = Bool


text:
using my package here
using Test

# these are subject to change as they develop, let me know when it breaks
@testset "Types" begin
    @test AbstractLinearCode <: AbstractCode
    @test AbstractCyclicCode <: AbstractLinearCode
    @test AbstractBCHCode <: AbstractCyclicCode
    @test AbstractReedSolomonCode <: AbstractBCHCode
    @test AbstractAdditiveCode <: AbstractCode
    @test AbstractQuantumCode <: AbstractAdditiveCode
    @test AbstractCSSCode <: AbstractQuantumCode
end

@testset "utils.jl" begin

end

@testset "linearcode.jl" begin
    # Hammming code
    F, _ = FiniteField(2, 1, "α");
    G = matrix(F, [1 0 0 0 0 1 1;
           0 1 0 0 1 0 1;
           0 0 1 0 1 1 0;
           0 0 0 1 1 1 1]);
    C = LinearCode(G);
    @test length(C) == 7
    @test rank(G) == dimension(C)
    @test dimension(C) == 4
    @test ismissing(minimumdistance(C))
    setminimumdistance!(C, 3)
    @test minimumdistance(C) == 3
    @test !isMDS(C)
    @test G == generatormatrix(C)
    @test G == originalgeneratormatrix(C)
    H = paritycheckmatrix(C)
    @test iszero(G * H')
    @test iszero(H * G')
    @test C ⊆ C
    D = dual(C)
    @test !(C ⊆ D)
    @test !issubcode(C, D)
    @test !isequivalent(C, D)
    @test !isselfdual(C)
    @test !isselforthogonal(C)
    cw = matrix(F, [1 0 0 0 0 1 1]);
    @test encode(C.G[:, 1], C) == cw
    # these v's are C.G[:, 1], just testing different formats
    v = [1, 0, 0, 0];
    @test encode(v, C) == cw
    v2 = [1; 0; 0; 0];
    @test encode(v2, C) == cw
    # this vector is the first row of the generator matrix and should
    # therefore have zero syndrome
    v = [1, 0, 0, 0, 0, 1, 1];
    @test iszero(syndrome(v, C))
    v = [1; 0; 0; 0; 0; 1; 1];
    @test iszero(syndrome(v, C))

    # lower rank test
    GandG = vcat(G, G);
    CGandG = LinearCode(G);
    @test rank(GandG) == dimension(CGandG)
    @test G == generatormatrix(CGandG)
    @test GandG == originalgeneratormatrix(CGandG)

end

ext: jl

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
# Test variational information

using Test
using Clustering

@testset "varinfo() (variational information)" begin

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]

@test varinfo(3, a1, 3, a1) ≈ 0.0 atol=1.0e-12
@test varinfo(2, a2, 2, a2) ≈ 0.0 atol=1.0e-12

v = varinfo(3, a1, 2, a2)
v_ = varinfo(2, a2, 3, a1)
@test 0.0 < v < log(3)
@test v ≈ v_

a1 = [1, 2, 3, 4, 5]
a2 = [1, 1, 1, 1, 1]

@test varinfo(5, a1, 1, a2) ≈ log(5)
@test varinfo(1, a2, 5, a1) ≈ log(5)

a1 = [1, 1, 1, 2, 2, 2]
a2 = [2, 2, 2, 1, 1, 1]
@test varinfo(2, a1, 2, a2) ≈ 0.0 atol=1.0e-12
@test varinfo(2, a1, 3, a2) ≈ 0.0 atol=1.0e-12
@test varinfo(4, a1, 3, a2) ≈ 0.0 atol=1.0e-12

end


text:
"""
A file with structs and functions for sampling hops and updating hopping rates
"""

# TODO to simplify the design can do two things:
# 1. make an abstract type HopRatesGridOptim (grid-optimized) and dispatch on it. Or maybe dispatch on the type of hop_const_cumulative_sums instead (Array{F} vs Array{Vector{F}})
# 2. make an abstract type HopRatesGraphDsi, which has hop rates uniform wrt neighbors, and treat things not of that type as having hop_const_cumulative_sums field

### spatial hop rates ###
abstract type AbstractHopRates end

HopRates(hopping_constants::Matrix{F}, spatial_system) where F <: Number = HopRatesGraphDsi(hopping_constants)
HopRates(hopping_constants::Matrix{F}, grid::Union{CartesianGridRej, CartesianGridIter}) where F <: Number = HopRatesGraphDsi(hopping_constants)

HopRates(hopping_constants::Matrix{Vector{F}}, spatial_system) where F <: Number = HopRatesGraphDsij(hopping_constants)
HopRates(hopping_constants::Matrix{Vector{F}}, grid::Union{CartesianGridRej, CartesianGridIter}) where F <: Number = HopRatesGridDsij(hopping_constants, grid)

HopRates(p::Pair{SpecHop, SiteHop}, spatial_system) where {SpecHop <: Vector{F}, SiteHop <: Vector{Vector{F}}} where F <: Number = HopRatesGraphDsLij(p...)
HopRates(p::Pair{SpecHop, SiteHop}, grid::Union{CartesianGridRej, CartesianGridIter}) where {SpecHop <: Vector{F}, SiteHop <: Vector{Vector{F}}} where F <: Number = HopRatesGridDsLij(p..., grid)

HopRates(p::Pair{SpecHop, SiteHop}, spatial_system) where {SpecHop <: Matrix{F}, SiteHop <: Vector{Vector{F}}} where F <: Number = HopRatesGraphDsiLij(p...)
HopRates(p::Pair{SpecHop, SiteHop}, grid::Union{CartesianGridRej, CartesianGridIter}) where {SpecHop <: Matrix{F}, SiteHop <: Vector{Vector{F}}} where F <: Number = HopRatesGridDsiLij(p..., grid)

"""
    update_hop_rates!(hop_rates::AbstractHopRates, species::AbstractArray, u, site, spatial_system)

update rates of all specs in species at site
"""
function update_hop_rates!(hop_rates::AbstractHopRates, species::AbstractArray, u, site, spatial_system)
    @inbounds for spec in species
        update_hop_rate!(hop_rates, spec, u, site, spatial_system)
    end
end

"""
    update_hop_rate!(hop_rates::HopRatesGraphDsi, species, u, site, spatial_system)

update rates of single species at site
"""
function update_hop_rate!(hop_rates::AbstractHopRates, species, u, site, spatial_system)
    rates = hop_rates.rates
    @inbounds old_rate = rates[species, site]
    @inbounds rates[species, site] = evalhoprate(hop_rates, u, species, site, spatial_system)
    @inbounds hop_rates.sum_rates[site] += rates[species, site] - old_rate
    old_rate
end

"""
    total_site_hop_rate(hop_rates::AbstractHopRates, site)

return total hopping rate out of site
"""
total_site_hop_rate(hop_rates::AbstractHopRates, site) = @inbounds hop_rates.sum_rates[site]

"""
    reset!(hop_rates::AbstractHopRates)

make all rates zero
"""
function reset!(hop_rates::AbstractHopRates)
    hop_rates.rates .= zero(eltype(hop_rates.rates))
    hop_rates.sum_rates .= zero(eltype(hop_rates.sum_rates))
    nothing
end

"""
    sample_species(hop_rates::AbstractHopRates, site, rng)

sample species to hop from site
"""
sample_species(hop_rates::AbstractHopRates, site, rng) = @inbounds linear_search((@view hop_rates.rates[:,site]), rand(rng) * total_site_hop_rate(hop_rates, site))

"""
    sample_hop_at_site(hop_rates::AbstractHopRates, site, rng, spatial_system)

sample a reaction at site, return (species, target_site)
"""
function sample_hop_at_site(hop_rates::AbstractHopRates, site, rng, spatial_system)
    species = sample_species(hop_rates, site, rng)
    return species, sample_target_site(hop_rates, site, species, rng, spatial_system)
end

############## hopping rates of form D_{s,i} ################
struct HopRatesGraphDsi{F} <: AbstractHopRates
    "hopping_constants[i,j] is the hop constant of species i at site j"
    hopping_constants::Matrix{F}

    "rates[i,j] is total hopping rate of species i at site j"
    rates::Matrix{F}

    "sum_rates[j] is the sum of hopping rates at site j"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsi)
    num_specs, num_sites = size(hop_rates.rates)
    println(io, "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} where s is species, and i is source.")
end

"""
    HopRatesGraphDsi(hopping_constants::Matrix{F}) where F <: Number

initializes HopRatesGraphDsi with zero rates
"""
function HopRatesGraphDsi(hopping_constants::Matrix{F}) where F <: Number
    rates = zeros(F, size(hopping_constants))
    HopRatesGraphDsi{F}(hopping_constants, rates, zeros(F, size(rates, 2)))
end

sample_target_site(hop_rates::HopRatesGraphDsi, site, species, rng, spatial_system) = rand_nbr(rng, spatial_system, site)

"""
return hopping rate of species at site
"""
function evalhoprate(hop_rates::HopRatesGraphDsi, u, species, site, spatial_system)
    @inbounds u[species,site]*hop_rates.hopping_constants[species,site]*outdegree(spatial_system, site)
end

############## hopping rates of form D_{s,i,j} ################
struct HopRatesGraphDsij{F} <: AbstractHopRates
    "hop_const_cumulative_sums[s,i] is the vector of cumulative sums of hopping constants of species s at site i"
    hop_const_cumulative_sums::Matrix{Vector{F}}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsij)
    num_specs, num_sites = size(hop_rates.rates)
    println(io, "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i,j} where s is species, i is source and j is destination.")
end

"""
    HopRatesGraphDsij(hopping_constants::Matrix{Vector{F}}) where F <: Number

initializes HopRates with zero rates
"""
function HopRatesGraphDsij(hopping_constants::Matrix{Vector{F}}; do_cumsum = true) where F <: Number
    do_cumsum && (hopping_constants = map(cumsum, hopping_constants))
    rates = zeros(F, size(hopping_constants))
    sum_rates = zeros(F, size(rates, 2))
    HopRatesGraphDsij{F}(hopping_constants, rates, sum_rates)
end

function sample_target_site(hop_rates::HopRatesGraphDsij, site, species, rng, spatial_system)
    @inbounds cumulative_hop_constants = hop_rates.hop_const_cumulative_sums[species, site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_nbr(spatial_system, site, n)
end

function evalhoprate(hop_rates::HopRatesGraphDsij, u, species, site, spatial_system)
    @inbounds u[species, site] * hop_rates.hop_const_cumulative_sums[species, site][end]
end

#################  hopping rates of form L_{s,i,j} optimized for cartesian grid  ######################
"""
Analogue of HopRatesGraphDsij, optimized for CartesianGrid
"""
struct HopRatesGridDsij{F} <: AbstractHopRates
    "hop_const_cumulative_sums[:,s,i] is the vector of cumulative sums of hopping constants of species s at site i"
    hop_const_cumulative_sums::Array{F, 3}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsij)
    num_specs, num_sites = size(hop_rates.rates)
    println(io, "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form L_{s,i,j} where s is species, i is source and j is destination.")
end

"""
    HopRatesGridDsij(hopping_constants::Matrix{Vector{F}}) where F <: Number

initializes HopRates with zero rates
"""
function HopRatesGridDsij(hopping_constants::Array{F, 3}; do_cumsum = true) where F <: Number
    do_cumsum && (hopping_constants = mapslices(cumsum, hopping_constants, dims = 1))
    rates = zeros(F, size(hopping_constants)[2:3])
    sum_rates = zeros(F, size(rates, 2))
    HopRatesGridDsij{F}(hopping_constants, rates, sum_rates)
end

function HopRatesGridDsij(hopping_constants::Matrix{Vector{F}}, grid) where F <: Number
    new_hopping_constants = Array{F, 3}(undef, 2*dimension(grid), size(hopping_constants)...)
    for ci in CartesianIndices(hopping_constants)
        species, site = Tuple(ci)
        nb_constants = @view new_hopping_constants[:, species, site]
        pad_hop_vec!(nb_constants, grid, site, hopping_constants[ci])
        cumsum!(nb_constants, nb_constants)
    end
    HopRatesGridDsij(new_hopping_constants, do_cumsum = false)
end

function sample_target_site(hop_rates::HopRatesGridDsij, site, species, rng, grid)
    @inbounds cumulative_hop_constants = @view hop_rates.hop_const_cumulative_sums[:,species, site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_potential_nbr(grid, site, n)
end

function evalhoprate(hop_rates::HopRatesGridDsij, u, species, site, spatial_system)
    @inbounds u[species, site] * hop_rates.hop_const_cumulative_sums[end, species, site]
end

############## hopping rates of form D_s * L_{i,j} ################
struct HopRatesGraphDsLij{F} <: AbstractHopRates
    "hopping constants of species -- D_s"
    species_hop_constants::Vector{F}

    "nbs_cumulative[i] is the vector of cumulative sums of hopping constants from site i to its neighbors"
    hop_const_cumulative_sums::Vector{Vector{F}}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsLij)
    num_specs, num_sites = length(hop_rates.species_hop_constants), length(hop_rates.hop_const_cumulative_sums)
    println(io, "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.")
end

"""
    HopRatesGraphDsLij(species_hop_constants::Vector{F}, site_hop_constants::Vector{Vector{F}}) where F <: Number

initializes HopRates with zero rates
"""
function HopRatesGraphDsLij(species_hop_constants::Vector{F}, site_hop_constants::Vector{Vector{F}}; do_cumsum = true) where F <: Number
    do_cumsum && (site_hop_constants = map(cumsum, site_hop_constants))
    rates = zeros(F, length(species_hop_constants), length(site_hop_constants))
    sum_rates = zeros(size(rates, 2))
    HopRatesGraphDsLij{F}(species_hop_constants, site_hop_constants, rates, sum_rates)
end

function sample_target_site(hop_rates::HopRatesGraphDsLij, site, species, rng, spatial_system)
    @inbounds cumulative_hop_constants = hop_rates.hop_const_cumulative_sums[site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_nbr(spatial_system, site, n)
end

function evalhoprate(hop_rates::HopRatesGraphDsLij, u, species, site, spatial_system)
    @inbounds u[species, site] * hop_rates.species_hop_constants[species] * hop_rates.hop_const_cumulative_sums[site][end]
end

############## hopping rates of form D_s * L_{i,j} optimized for cartesian grid ################
struct HopRatesGridDsLij{F} <: AbstractHopRates
    "hopping constants of species -- D_s"
    species_hop_constants::Vector{F}

    "nbs_cumulative[i] is the vector of cumulative sums of hopping constants from site i to its neighbors"
    hop_const_cumulative_sums::Matrix{F}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsLij)
    num_specs, num_sites = length(hop_rates.species_hop_constants), size(hop_rates.hop_const_cumulative_sums, 2)
    println(io, "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.")
end

function HopRatesGridDsLij(species_hop_constants::Vector{F}, site_hop_constants::Matrix{F}; do_cumsum = true) where F <: Number
    do_cumsum && (site_hop_constants = mapslices(cumsum, site_hop_constants, dims = 1))
    rates = zeros(F, length(species_hop_constants), size(site_hop_constants, 2))
    sum_rates = zeros(size(rates, 2))
    HopRatesGridDsLij{F}(species_hop_constants, site_hop_constants, rates, sum_rates)
end

function HopRatesGridDsLij(species_hop_constants::Vector{F}, site_hop_constants::Vector{Vector{F}}, grid) where F <: Number
    new_hopping_constants = Matrix{F}(undef, 2*dimension(grid), length(site_hop_constants))
    for site in 1:length(site_hop_constants)
        nb_constants = @view new_hopping_constants[:, site]
        pad_hop_vec!(nb_constants, grid, site, site_hop_constants[site])
        cumsum!(nb_constants, nb_constants)
    end
    HopRatesGridDsLij(species_hop_constants, new_hopping_constants, do_cumsum = false)
end

function sample_target_site(hop_rates::HopRatesGridDsLij, site, species, rng, grid)
    @inbounds cumulative_hop_constants = @view hop_rates.hop_const_cumulative_sums[:, site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_potential_nbr(grid, site, n)
end

function evalhoprate(hop_rates::HopRatesGridDsLij, u, species, site, grid)
    @inbounds u[species, site] * hop_rates.species_hop_constants[species] * hop_rates.hop_const_cumulative_sums[end, site]
end

############## hopping rates of form D_{s,i} * L_{i,j} ################
struct HopRatesGraphDsiLij{F} <: AbstractHopRates
    "hopping constants of species -- D_s"
    species_hop_constants::Matrix{F}

    "nbs_cumulative[i] is the vector of cumulative sums of hopping constants from site i to its neighbors"
    hop_const_cumulative_sums::Vector{Vector{F}}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsiLij)
    num_specs, num_sites = size(hop_rates.species_hop_constants)
    println(io, "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.")
end

function HopRatesGraphDsiLij(species_hop_constants::Matrix{F}, site_hop_constants::Vector{Vector{F}}; do_cumsum = true) where F <: Number
    @assert size(species_hop_constants, 2) == length(site_hop_constants)
    do_cumsum && (site_hop_constants = map(cumsum, site_hop_constants))
    rates = zeros(F, length(species_hop_constants), length(site_hop_constants))
    sum_rates = zeros(size(rates, 2))
    HopRatesGraphDsiLij{F}(species_hop_constants, site_hop_constants, rates, sum_rates)
end

function sample_target_site(hop_rates::HopRatesGraphDsiLij, site, species, rng, spatial_system)
    @inbounds cumulative_hop_constants = hop_rates.hop_const_cumulative_sums[site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_nbr(spatial_system, site, n)
end

function evalhoprate(hop_rates::HopRatesGraphDsiLij, u, species, site, spatial_system)
    @inbounds u[species, site] * hop_rates.species_hop_constants[species, site] * hop_rates.hop_const_cumulative_sums[site][end]
end

############## hopping rates of form D_{s,i} * L_{i,j} optimized for cartesian grid ################
struct HopRatesGridDsiLij{F} <: AbstractHopRates
    "hopping constants of species -- D_s"
    species_hop_constants::Matrix{F}

    "nbs_cumulative[i] is the vector of cumulative sums of hopping constants from site i to its neighbors"
    hop_const_cumulative_sums::Matrix{F}

    "rates[s,i] is the total hopping rate of species s at site i"
    rates::Matrix{F}

    "sum_rates[i] is the total hopping rate out of site i"
    sum_rates::Vector{F}
end

function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsiLij)
    num_specs, num_sites = length(hop_rates.species_hop_constants), size(hop_rates.hop_const_cumulative_sums, 2)
    println(io, "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.")
end

function HopRatesGridDsiLij(species_hop_constants::Matrix{F}, site_hop_constants::Matrix{F}; do_cumsum = true) where F <: Number
    @assert size(species_hop_constants, 2) == size(site_hop_constants, 2)
    do_cumsum && (site_hop_constants = mapslices(cumsum, site_hop_constants, dims = 1))
    rates = zeros(F, size(species_hop_constants))
    sum_rates = zeros(size(rates, 2))
    HopRatesGridDsiLij{F}(species_hop_constants, site_hop_constants, rates, sum_rates)
end

function HopRatesGridDsiLij(species_hop_constants::Matrix{F}, site_hop_constants::Vector{Vector{F}}, grid) where F <: Number
    new_hopping_constants = Matrix{F}(undef, 2*dimension(grid), length(site_hop_constants))
    for site in 1:length(site_hop_constants)
        nb_constants = @view new_hopping_constants[:, site]
        pad_hop_vec!(nb_constants, grid, site, site_hop_constants[site])
        cumsum!(nb_constants, nb_constants)
    end
    HopRatesGridDsiLij(species_hop_constants, new_hopping_constants, do_cumsum = false)
end

function sample_target_site(hop_rates::HopRatesGridDsiLij, site, species, rng, grid)
    @inbounds cumulative_hop_constants = @view hop_rates.hop_const_cumulative_sums[:, site]
    @inbounds n = searchsortedfirst(cumulative_hop_constants, rand(rng) * cumulative_hop_constants[end])
    return nth_potential_nbr(grid, site, n)
end

function evalhoprate(hop_rates::HopRatesGridDsiLij, u, species, site, grid)
    @inbounds u[species, site] * hop_rates.species_hop_constants[species, site] * hop_rates.hop_const_cumulative_sums[end, site]
end


text:
module VariationalInequalitySolver

using FastClosures, LinearAlgebra, Logging, NLPModels, SparseArrays, Stopping

import NLPModels:
  residual!,
  jac_structure_residual!,
  jac_coord_residual!,
  jprod_residual!,
  jtprod_residual!,
  jac_op_residual!,
  hess_structure_residual!,
  hess_coord_residual!,
  hprod_residual!,
  hess_op_residual!

abstract type AbstractVIModel{T, S} end

abstract type AbstractVIMeta{T, S} end

mutable struct VIMeta{T, S} <: AbstractVIMeta{T, S}
  x0::S
  nvar
  nnzj
  nnzh
end

#largely inspired from https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/main/src/nls/api.jl
include("projector/proj_nls.jl")
include("model/api.jl")

export NLSVIModel

mutable struct NLSVIModel{S, T, NLP <: AbstractNLSModel{T, S}} <: AbstractVIModel{T, S}
  meta::VIMeta{T, S}
  nls::NLP
  function NLSVIModel(nls::AbstractNLSModel{T, S}) where {T, S}
    @lencheck nls.nls_meta.nequ nls.meta.x0 #test that nls.meta.nvar == nls.nls_meta.nequ
    return new{S, T, typeof(nls)}(
      VIMeta{T, S}(nls.meta.x0, nls.meta.nvar, nls.nls_meta.nnzj, nls.nls_meta.nnzh),
      nls,
    )
  end
end

for meth in [
  :residual!,
  :jac_structure_residual!,
  :jac_coord_residual!,
  :jprod_residual!,
  :jtprod_residual!,
  :jac_op_residual!,
  :hess_structure_residual!,
  :hess_coord_residual!,
  :hprod_residual!,
  :hess_op_residual!,
]
  @eval begin
    $meth(model::NLSVIModel, args...; kwargs...) = $meth(model.nls, args...; kwargs...)
  end
end

using CaNNOLeS

function project!(model::NLSVIModel, d::AbstractVector{T}, Px::AbstractVector{T}) where {T, S}
  # Here we need to solve the optimization problem
  proj = NLSProjector(model.nls, d)
  stats = with_logger(NullLogger()) do
    cannoles(proj)
  end
  if stats.status != :first_order
    @warn "There was an error in the projection computation"
  end
  Px .= stats.solution
  return Px
end

include("solvers/projectionVI.jl")
include("solvers/projected_reflected_gradient.jl")

end


text:
function GCD(a::Int64,b::Int64)

	while !(mod(a,b) == 0)

		c = mod(a,b);
		a = b;
		b = c;
	end
	return b
end

function LCM(a::Int64, b::Int64)
	return div(a*b,GCD(a,b))
end

function LCMMulti(N::Array{Int64,1})
	return reduce(LCM, N)
end

function Euler5()
	println(LCMMulti([1:20]))
end

@time Euler5()

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
# Model
include_model("hopper")

# Dimensions
nq = 4 # configuration dimension
nu = 2 # control dimension
nc = 1 # number of contact points
nf = 2 # number of faces for friction cone
nb = nc * nf
ns = nq

# Parameters
g = 9.81 # gravity
μ = 1.0  # coefficient of friction
mb = 10.0 # body mass
ml = 1.0  # leg mass
Jb = 2.5 # body inertia
Jl = 0.25 # leg inertia

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

model = Hopper(n, m, d,
			   mb, ml, Jb, Jl,
			   μ, g,
			   qL, qU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

# Free-time model
model_ft = free_time_model(model)

function fd(model::Hopper, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	s = view(u, model.idx_s)
	h = u[end]

	[q2⁺ - q2⁻;
	((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
	- M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
	+ transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{2}(b)
	- h * G_func(model, q2⁺)) + s]
end

function maximum_dissipation(model::Hopper, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(model.nb)
	η = u[model.idx_η]
	h = u[end]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

# Horizon
T = 21

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= 25.0
_uu[model_ft.idx_s] .= 0.0
_uu[end] = 1.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -25.0
_ul[model_ft.idx_s] .= 0.0
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.5 , 0.0, 0.5]
q_right = [-0.25, 0.5 + 0.5, pi / 2.0, 0.25]
q_top = [-0.5, 0.5 + 1.0, pi, 0.25]
q_left = [-0.75, 0.5 + 0.5, 3.0 * pi / 2.0, 0.25]
qT = [-1.0, 0.5,  2.0 * pi, 0.5]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; q1],
		xT = [Inf * ones(model_ft.nq); qT])

q_ref = [linear_interpolation(q1, q_right, 6)...,
         linear_interpolation(q_right, q_top, 6)[2:end]...,
         linear_interpolation(q_top, q_left, 6)[2:end]...,
         linear_interpolation(q_left, qT, 6)[2:end]...]

x_ref = configuration_to_state(q_ref)

# Objective
include_objective(["velocity"])
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(1.0 * [1.0; 1.0; 0.0; 0.0; 1.0; 1.0; 0.0; 0.0]) for t = 1:T],
    [Diagonal([1.0e-1, 1.0e-1, zeros(model_ft.m - model_ft.nu)...]) for t = 1:T-1],
    [x_ref[t] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T-1],
    1.0)

obj_velocity = velocity_objective([Diagonal(1.0 * ones(model_ft.nq)) for t = 1:T],
	model_ft.nq, idx_angle = (3:3))
obj = MultiObjective([obj_tracking, obj_velocity])

# Constraints
include_constraints(["contact_al", "free_time"])
con_contact = contact_al_constraints(model_ft, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_free_time, con_contact])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = deepcopy(x_ref) # linear interpolation on state
u0 = [[1.0e-5 * rand(model_ft.m - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

using LinearAlgebra, ForwardDiff, SparseArrays, Optim, LineSearches
include(joinpath(pwd(),"src/solvers/augmented_lagrangian.jl"))

function f(x)
    MOI.eval_objective(prob, x)
end

function g!(G, x)
    MOI.eval_objective_gradient(prob, G, x)
    nothing
end

function c!(c, x)
    MOI.eval_constraint(prob, c, x)
    c .*= -1.0
    nothing
end

spar = sparsity_jacobian(prob)
global jac = zeros(length(spar))
global ii = [s[1] for s in spar]
global jj = [s[2] for s in spar]
function d!(D, x)
    MOI.eval_constraint_jacobian(prob, jac, x)
    D .= sparse(ii, jj, -1.0 .* jac)
end

function f_al(x, al::AugmentedLagrangian)
    # evaluate constraints
    c!(al.c, x)
    bounds!(al, x)
    active_set_update!(al)

    # compute objective
    J = f(x)

    # add augmented Lagrangian terms
    J += al.λ' * al.c + 0.5 * sum(al.as .* al.ρ .* (al.c.^2.0))
    J += al.λl' * al.cl + 0.5 * sum(al.asl .* al.ρl .* (al.cl.^2.0))
    J += al.λu' * al.cu + 0.5 * sum(al.asu .* al.ρu .* (al.cu.^2.0))
end

function g_al!(G, x, al::AugmentedLagrangian)
    # compute objective gradient
    g!(G, x)

    # evaluate constraints
    # ForwardDiff.jacobian!(al.∇c, c!, al.c, x)
    c!(al.c, x)
    d!(al.∇c, x)

    bounds!(al, x)
    active_set_update!(al)

    # add augmented Lagrangian gradient terms
    G .+= al.∇c' * (al.λ + al.as .* al.ρ .* al.c)
    G[al.idx_l] -= (al.λl + al.asl .* al.ρl .* al.cl)
    G[al.idx_u] += (al.λu + al.asu .* al.ρu .* al.cu)
    return nothing
end

function solve(x, al; alg = :LBFGS, max_iter = 5, c_tol = 1.0e-3)
	println("*augmented Lagrangian solve*")
    # reset augmented Lagrangian
    reset!(al)
    println("	solving...")
    for i = 1:max_iter

        # update augmented Lagrangian methods
        _f(z) = f_al(z, al)
        _g!(G, z) = g_al!(G, z, al)

        # solve
        sol = Optim.optimize(_f, _g!, x, @eval $alg()) # linesearch = LineSearches.BackTracking()

        # evaluate constraints
        x = sol.minimizer
        c_max = constraint_violation(al, x)
        println("	iter: $i")
        println("	c_max: $c_max")

        # check for convergence -> update augmented Lagrangian
        if c_max < c_tol
			println("solve: success")
            return x, sol
        else
			if i >= max_iter
				println("solve: failed")
				return x, sol
			end

            c!(al.c, x)
            bounds!(al, x)
            update!(al)
        end
    end
end

n = prob.num_var
m = prob.num_con
xl, xu = prob.primal_bounds
cl, cu = prob.constraint_bounds
idx_ineq = (1:m)[cu .> cl]
sum(isfinite.(xl))
sum(isfinite.(xu))
sum(cu .> cl)

al = augmented_lagrangian(n, m,
    xl = xl, xu = xu, ρ0 = 1.0, s = 10.0,
    idx_ineq = idx_ineq)

f_al(z0, al)
_c0 = zero(z0)
g_al!(_c0, z0, al)
_c0
@time x_sol_al, sol , info = solve(copy(z0), al,
    alg = :BFGS, max_iter = 2, c_tol = 1.0e-2)

# Visualize
using Plots
x̄, ū = unpack(x_sol_al, prob)
plot(hcat(x̄...)', width = 2.0)
plot(hcat(ū...)[1:model.nu, :]', width = 2.0, linetype = :steppost)
s̄ = [ū[t][model.idx_s] for t = 1:T-1]
plot(hcat(s̄...)', width = 2.0, linetype = :steppost)

# using Plots
tf, t, h = get_time(ū)
plot(t[1:end-1], hcat(ū...)[1:2,:]', linetype=:steppost,
	xlabel="time (s)", ylabel = "control",
	label = ["angle" "length"],
	width = 2.0, legend = :top)
plot(t[1:end-1], h, linetype=:steppost)
#
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model_ft, state_to_configuration(x̄), Δt = h[1])


text:
import FLOWFarm; const ff = FLOWFarm

turbine_x = [0.0]
turbine_y = [0.0]
turbine_z = [0.0]
turbine_yaw = [0.0]

nturbines = length(turbine_x)
diam = 126.4

# set turbine design parameters
rotor_diameter = zeros(nturbines) .+ diam # m
hub_height = zeros(nturbines) .+ 90.0   # m
cut_in_speed = zeros(nturbines) .+3.  # m/s
cut_out_speed = zeros(nturbines) .+25.  # m/s
rated_speed = zeros(nturbines) .+11.4  # m/s
rated_power = zeros(nturbines) .+5.0E6  # W
generator_efficiency = zeros(nturbines) .+ 0.944

# load power curve
cpctdata = readdlm("inputfiles/NREL5MWCPCT.txt",  ' ', skipstart=1)
velpoints = cpctdata[:,1]
cppoints = cpctdata[:,2]

# initialize power model
power_model = ff.PowerModelCpPoints(velpoints, cppoints)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# load thrust curve
ctpoints = cpctdata[:,3]

# initialize thurst model
ct_model = ff.ThrustModelCtPoints(velpoints, ctpoints)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end

sorted_turbine_index = [1]

wind_speed = 8.1
air_density = 1.1716  # kg/m^3
winddirections = [270.0*pi/180.0]
windspeeds = [wind_speed]
windprobabilities = [1.0]
measurementheights = [hub_height[1]]
wtvelocities = [wind_speed]
ambient_tis = [0.1]
shearexponent = 0.0
wind_shear_model = ff.PowerLawWindShear(shearexponent)
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)

# set up wake and related models
wakedeficitmodel = ff.GaussYawVariableSpread()
wakedeflectionmodel = ff.GaussYawVariableSpreadDeflection()
wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()

model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)


text:
using LinearAlgebra

function nearest_neighbors(X, p, r)
    neighbors = Int[]
    for i in 1:size(X, 2)
        if norm(X[:,i] - p) < r
            push!(neighbors, i)
        end
    end
    return neighbors
end


text:
# Function to parse instruction line
# into inst, operator, and val
function parse_instruction(code)
    inst, num = split(code, " ")
    operator = string(num[1])
    # operator = Meta.parse(operator)
    val = parse(Int, num[2:length(num)])

    return inst, operator, val
end

# run the program as given
function run_program(file)
    global i = 1                # line counter
    global acc = 0              # accumulator
    global read_lines = []      # storage for already read lines
    global len = length(file)   # length of file
    global valid_exit = false   # false if infinite loop

    # Run this loop until a break
    while true
        if i in read_lines  # check for infinite loop
            valid_exit = false
            break
        elseif i == (len+1) # check for valid program halt
            valid_exit = true
            break
        end

        inst, operator, val = parse_instruction(file[i])

        append!(read_lines, i)  # add line to read_lines

        # Check each inst and change i/acc accordingly
        if inst == "acc"
            exp = Meta.parse(join([acc, operator, val]))
            global acc = eval(exp)
            global i += 1

        elseif inst == "nop"
            global i += 1

        elseif inst == "jmp"
            exp = Meta.parse(join([i, operator, val]))
            i = eval(exp)
        end

    end
    return acc, valid_exit
end

# Function to change the program and test
# all possible single changes of "jmp" <=> "nop"
function change_program(file)

    new_file = copy(file)   # temporary instruction set
    for i in  1:length(new_file)
        newline = new_file[i]

        if newline[1:3] == "jmp"
            newline = replace(newline, "jmp" => "nop")
        elseif newline[1:3] == "nop"
            newline = replace(newline, "nop" => "jmp")
        end

        new_file[i] = newline
        results = run_program(new_file)

        if results[2] == true  # if run_program was valid
            return results
        else
            new_file = copy(file)
        end
    end
end
code_file = readlines("./code.txt")
# code_file = readlines("./test_code.txt")

println("Part 1: acc = ", run_program(code_file), " before infinite loop")
println("Part 2: acc = ", change_program(code_file), " after successful completion")
# println(run_program(code_file, change_inst=true))


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0 1 0]
dims = [0 0 0]
MPI.Dims_create!(nnodes, dims)

cperiods = Cint.(periods[:])
cdims = Cint.(dims[:])
comm_cart = MPI.Cart_create(comm, ndims, cdims, cperiods, reorder)
@test MPI.Comm_size(comm_cart) == nnodes

comm_cart2 = MPI.Cart_create(comm, dims, periods, reorder)
@test MPI.Comm_size(comm_cart2) == nnodes

MPI.Finalize()


text:
# simple periodic system runner.
#
# This file defines the periodic system problem
# with background noise and sparse large
# deviations. Both (global) L2 and Robust fitting are
# run for several background noise levels
# and the quality of the resulting eigenvalues is
# recorded and written to file
#

using RobustDMD, Printf, Random, JLD

################################################################
################################################################
################################################################
# Define some helper routines

function periodic_system_data(;nt=128)
################################################
# Generate data for simple periodic linear system
#
# z(t) = exp(At)z0
#
# where A has eigenvalues equal to +/- i
#
################################################

    z0 = [1.0; 0.1]
    A = [1 -2; 1 -1]
    dt = 0.1

    evals = [im; -im]

    z = Array{Complex{Float64}}(undef,nt,2)
    t = Array{Complex{Float64}}(undef,nt)

    for m = 1:nt
        t[m] = dt*m
        ztemp = exp(A*t[m])*z0
        z[m,:] = transpose(ztemp)
    end

    return z, evals, t
end

function make_corrupt_data(xclean,sigma,mu,p)
    #------------------------------------------
    # Generate Corrupted Data - Spikes
    #------------------------------------------
    #
    # sigma - background noise
    # mu - spike size
    # p - spike frequency

    m, n = size(xclean)

    noise = sigma*randn(m,n) .+ mu*randn(m,n).*( rand(m,n) .< p )
    xdat = xclean .+ noise

    return xdat
end

function dofits(xdat,t,eigs,kappa,optionsl2,optionshuber)
    # this wraps up the calls to the DMD solver
    # Here we use the robust but generally slower
    # proximal gradient routine for both L2 and Huber
    # fits. On return, we provide the error for these
    # fits and the exact DMD error.

    k = length(eigs)

    #-------------------------------------------------------------
    # Initial Guess
    #-------------------------------------------------------------
    m,n = size(xdat);

    a0, B0 = dmdexactestimate(m, n, k, xdat, t);

    ainit, Binit = dmdexactestimate(m,n,k,xdat,t,dmdtype="trap")

    #----------------------------------------------------
    # l2 Trial
    #----------------------------------------------------
    #println("running l2 experiment...");
    # loss functions
    lossf = (z) -> l2_func(z);
    lossg = (z) -> l2_grad!(z);

    l2params = DMDParams(k, xdat, t, lossf, lossg,
                         inner_directl2=true);
    copyto!(l2params.a,ainit)

    # apply solver
    obj_his1, err_his1 = solveDMD_withPG(l2params, optionsl2);

    #----------------------------------------------------
    # huber Trial
    #----------------------------------------------------
    #println("running huber experiment...");
    # loss functions
    lossf = (z) -> huber_func(z,kappa);
    lossg = (z) -> huber_grad!(z,kappa);

    huberparams = DMDParams(k, xdat, t, lossf, lossg);
    copyto!(huberparams.a,ainit)

    # apply solver
    obj_his2, err_his2 = solveDMD_withPG(huberparams, optionshuber);

    #--------------------------------------------------------------------
    # Evaluate Errors and Get Best Eigenvalue Permutation
    #--------------------------------------------------------------------

    err1, p1 = besterrperm_wi(l2params.a, eigs);
    err2, p2 = besterrperm_wi(huberparams.a, eigs);
    err3, p3 = besterrperm_wi(a0, eigs);

    return err1, err2, err3, l2params.a[p1], huberparams.a[p2],
    a0[p3], l2params.B[p1,:], huberparams.B[p2,:], B0[p3,:], obj_his1,
    err_his1, obj_his2, err_his2

end

################################################################
################################################################
################################################################
# Actually run the tests...

# make data

xclean, eigs, t = periodic_system_data()

# solver parameters

function prox_mr(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    for i = 1:k
        ar[i<<1-1] = min(1.0, ar[i<<1-1]);
    end
end

itm = 1000;
tol = 1e-7;
ptf = itm+1;
optsl2 = DMD_PG_Options(itm=itm, tol=tol, ptf=ptf, prox=prox_mr);
optshub = DMD_PG_Options(itm=itm, tol=tol, ptf=ptf, prox=prox_mr);

# test parameters

seed = 8675309
Random.seed!(seed)

p = 0.05 # add spikes at this percentage of data points
k = 2 # rank of fit

nj = 5 # number of different levels of background noise
ntest = 200

errs = Array{Float64}(undef,3,ntest,nj)
eigs_comp = Array{Complex{Float64}}(undef,k,3,ntest,nj)
params = Array{Float64}(undef,4,nj)

@printf "==================================================\n\n"
@printf "Running the simple periodic example... \n\n"

for j = 1:nj
    @printf "starting noise level %d of %d\n" j nj
    sigma = 10.0^(-6+j) # noise floor
    mu = 1.0 # size of spikes
    kappa = 5*sigma
    nu = sqrt(mu)
    params[1,j] = sigma
    params[2,j] = mu
    params[3,j] = p
    params[4,j] = kappa

    Threads.@threads for i = 1:ntest
        xdat = make_corrupt_data(xclean,sigma,mu,p)

        (err1, err2, err3, eig1, eig2, eig3, b1, b2, b3,
         obj_his1, err_his1, obj_his2, err_his2) = (
             dofits(xdat,t,eigs,kappa,optsl2,optshub))

        #@printf "test %d %d %d\n"  j l i
        #@printf "errs %e %e %e\n" err1 err2 err3
        errs[1,i,j] = err1
        errs[2,i,j] = err2
        errs[3,i,j] = err3

        eigs_comp[:,1,i,j] = eig1[:]
        eigs_comp[:,2,i,j] = eig2[:]
        eigs_comp[:,3,i,j] = eig3[:]

    end
end

dirname = @__DIR__
fname = dirname * "/results/simple_periodic_pg_out.jld"

@printf "writing results to file %s\n" fname

# save the results to file

mkpath(dirname * "/results")
file = jldopen(fname,"w")
file["errs"] = errs
file["eigs_comp"] = eigs_comp
file["params"] = params
close(file)

@printf "\n ...done\n"

@printf "==================================================\n"


text:
import IntervalArithmetic: diam, isinterior, bisect, isnan

export branch_and_prune, Bisection, Newton
export BreadthFirstSearch, DepthFirstSearch

diam(r::Root) = diam(interval(r))
isnan(X::IntervalBox) = any(isnan.(X))
isnan(r::Root) = isnan(interval(r))

"""
    BreadthFirstSearch <: BreadthFirstBBSearch

Type implementing  the `BreadthFirstBBSearch` interface for interval roots finding.

# Fields:
  - `initial`: region (as a `Root` object) in which roots are searched.
  - `contractor`: contractor to use (`Bisection`, `Newton` or `Krawczyk`)
  - `tol`: tolerance of the search
"""
struct BreadthFirstSearch{R <: Region, C <: Contractor, T <: Real} <: BreadthFirstBBSearch{Root{R}}
    initial::Root{R}
    contractor::C
    tol::T
end

"""
    DepthFirstSearch <: DepthFirstBBSearch

Type implementing the `DepthFirstBBSearch` interface for interval roots finding.

# Fields:
  - `initial`: region (as a `Root` object) in which roots are searched.
  - `contractor`: contractor to use (`Bisection`, `Newton` or `Krawczyk`)
  - `tol`: tolerance of the search
"""
struct DepthFirstSearch{R <: Region, C <: Contractor, T <: Real} <: DepthFirstBBSearch{Root{R}}
    initial::Root{R}
    contractor::C
    tol::T
end

BreadthFirstSearch(X::Region, C::Contractor, tol::Real) = BreadthFirstSearch(Root(X, :unknown), C, tol)
DepthFirstSearch(X::Region, C::Contractor, tol::Real) = DepthFirstSearch(Root(X, :unknown), C, tol)

root_element(search::BBSearch{Root{R}}) where {R <: Region} = search.initial

function bisect(r::Root)
    Y1, Y2 = bisect(interval(r))
    return Root(Y1, :unknown), Root(Y2, :unknown)
end

bisect(::BBSearch, r::Root) = bisect(r::Root)

function process(search::Union{BreadthFirstSearch, DepthFirstSearch}, r::Root)
    contracted_root = search.contractor(r, search.tol)
    status = root_status(contracted_root)

    status == :unique && return :store, contracted_root
    status == :empty && return :discard, contracted_root

    if status == :unknown
        # Avoid infinite division of intervals with singularity
        isnan(contracted_root) && diam(r) < search.tol && return :store, r
        diam(contracted_root) < search.tol && return :store, contracted_root

        return :bisect, r
    else
        error("Unrecognized root status: $status")
    end
end

"""
    branch_and_prune(X, contractor, strategy, tol)

Generic branch and prune routine for finding isolated roots using the given
contractor to determine the status of a given box `X`.

See the documentation of the `roots` function for explanation of the other
arguments.
"""
function branch_and_prune(r::Root, contractor, search, tol)
    iter = search(r, contractor, tol)
    local endstate
    # complete iteration
    for state in iter
        endstate = state
    end
    return data(endstate)
end

const NewtonLike = Union{Type{Newton}, Type{Krawczyk}}
const default_strategy = DepthFirstSearch
const default_tolerance = 1e-7
const default_contractor = Newton

"""
    roots(f, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15)
    roots(f, deriv, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15)
    roots(f, X, contractor, tol)
    roots(f, deriv, X, contractor, tol)

Uses a generic branch and prune routine to find in principle all isolated roots
of a function `f:R^n → R^n` in a region `X`, if the number of roots is finite.

Inputs:
  - `f`: function whose roots will be found
  - `X`: `Interval` or `IntervalBox` in which roots are searched
  - `contractor`: function that, when applied to the function `f`, determines
    the status of a given box `X`. It returns the new box and a symbol indicating
    the status. Current possible values are `Bisection`, `Newton` and `Krawczyk`
  - `deriv`: explicit derivative of `f` for `Newton` and `Krawczyk`
  - `strategy`: `SearchStrategy` determining the order in which regions are
    processed.
  - `tol`: Absolute tolerance. If a region has a diameter smaller than `tol`, it
    is returned with status `:unknown`.

"""
function roots(f::Function, X, contractor::Type{C}=default_contractor,
               strategy::Type{S}=default_strategy,
               tol::Float64=default_tolerance) where {C <: Contractor, S <: BBSearch}

    _roots(f, X, contractor, strategy, tol)
end

function roots(f::Function, deriv::Function, X, contractor::Type{C}=default_contractor,
               strategy::Type{S}=default_strategy,
               tol::Float64=default_tolerance) where {C <: Contractor, S <: BBSearch}

    _roots(f, deriv, X, contractor, strategy, tol)
end

function roots(f::Function, X, contractor::Type{C},
               tol::Float64) where {C <: Contractor}

    _roots(f, X, contractor, default_strategy, tol)
end

function roots(f::Function, deriv::Function, X, contractor::Type{C},
               tol::Float64) where {C <: Contractor}

    _roots(f, deriv, X, contractor, default_strategy, tol)
end

#===
    More specific `roots` methods (all parameters are present)
    These functions are called `_roots` to avoid recursive calls.
===#

# For `Bisection` method
function _roots(f, r::Root{T}, ::Type{Bisection},
               strategy::Type{S}, tol::Float64) where {T, S <: BBSearch}

    branch_and_prune(r, Bisection(f), strategy, tol)
end

# For `NewtonLike` acting on `Interval`
function _roots(f, r::Root{Interval{T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {T, S <: BBSearch}

    deriv = x -> ForwardDiff.derivative(f, x)
    _roots(f, deriv, r, contractor, strategy, tol)
end

function _roots(f, deriv, r::Root{Interval{T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {T, S <: BBSearch}

    branch_and_prune(r, contractor(f, deriv), strategy, tol)
end

# For `NewtonLike` acting on `IntervalBox`
function _roots(f, r::Root{IntervalBox{N, T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {N, T, S <: BBSearch}

    deriv = x -> ForwardDiff.jacobian(f, x)
    _roots(f, deriv, r, contractor, strategy, tol)
end

function _roots(f, deriv, r::Root{IntervalBox{N, T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {N, T, S <: BBSearch}

    branch_and_prune(r, contractor(f, deriv), strategy, tol)
end

# Acting on `Interval`
function _roots(f, X::Region, contractor::Type{C},
               strategy::Type{S}, tol::Float64) where {C <: Contractor, S <: BBSearch}

    _roots(f, Root(X, :unknown), contractor, strategy, tol)
end

function _roots(f, deriv, X::Region, contractor::Type{C},
               strategy::Type{S}, tol::Float64) where {C <: Contractor, S <: BBSearch}

    _roots(f, deriv, Root(X, :unknown), contractor, strategy, tol)
end

# Acting on `Vector` of `Root`
function _roots(f, V::Vector{Root{T}}, contractor::Type{C},
               strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch}

    vcat(_roots.(f, V, contractor, strategy, tol)...)
end

function _roots(f, deriv, V::Vector{Root{T}}, contractor::Type{C},
               strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch}

    vcat(_roots.(f, deriv, V, contractor, strategy, tol)...)
end

# Acting on complex `Interval`
function _roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
               strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch}

    g = realify(f)
    Y = IntervalBox(reim(Xc)...)
    rts = _roots(g, Root(Y, :unknown), contractor, strategy, tol)

    return [Root(Complex(root.interval...), root.status) for root in rts]
end

function _roots(f, Xc::Complex{Interval{T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {T, S <: BBSearch}

    g = realify(f)
    g_prime = x -> ForwardDiff.jacobian(g, x)
    Y = IntervalBox(reim(Xc)...)
    rts = _roots(g, g_prime, Root(Y, :unknown), contractor, strategy, tol)

    return [Root(Complex(root.interval...), root.status) for root in rts]
end

function _roots(f, deriv, Xc::Complex{Interval{T}}, contractor::NewtonLike,
               strategy::Type{S}, tol::Float64) where {T, S <: BBSearch}

    g = realify(f)
    g_prime = realify_derivative(deriv)
    Y = IntervalBox(reim(Xc)...)
    rts = _roots(g, g_prime, Root(Y, :unknown), contractor, strategy, tol)

    return [Root(Complex(root.interval...), root.status) for root in rts]
end


text:
function get_draw(s::Symbol, custom_draw::Function)
    return (args...) -> begin
        custom_draw(args...)
        get_draw(s)(args...)
    end
end

function get_draw(s::Symbol)
    return (video, object, frame; kwargs...) -> begin
        l = s == :graph ? GRAPHS[object.opts[:_graph_idx]] :
                          s == :vertex ? GRAPH_VERTICES[object.opts[:_vertex_idx]] :
                                          GRAPH_EDGES[object.opts[:_edge_idx]]
        kw = merge(object.opts, Dict(collect(kwargs)...))
        for (_, style) in l.opts[:styles]
            style(video, object, frame; kw...)
            if frame == first(get_frames(object))
                kw = merge(object.opts[:_style_opts_cache], kw)
            end
        end
        Luxor.clipreset()
        object.opts = merge(object.opts, object.opts[:_style_opts_cache])
    end
end

"""
    @add_styles(component::AbstractJavisGraphElement, draw)

Aggregate all the drawing styles and store as part of the current object.
"""
macro add_styles(component, draw...)
    # Process drawing functions
    return quote
        # Inside a quote - end block
        # $(d) will turn Expr([a, b, c]) to [a, b, c]
        # $(d...) will turn [Expr(a()), Expr(b())] to [a(), b()]
        # $(d...)... further expanded the generated list above
        # since component is a variable which might be created inside a function
        # we need to escape (with esc()) it since it is not a global variable
        # if !haskey($(esc(component)).opts, :styles)
        #     $(esc(component)).opts[:styles] = Function[(args...; kw...) -> Luxor.clipreset()]
        # end
        # append!($(esc(component)).opts[:styles], [$(esc(draw...))...])
        merge!($(esc(component)).opts[:styles], OrderedDict(collect([$(esc(draw...))...])))
    end
end

"""
    @register_style_opts(component::AbstractObject, styles...)

Register the style for the current object
"""
macro register_style_opts(component, style_opts...)
    return quote
        # Get the name and the value for the arguments passed
        for (n, v) in zip([$(style_opts)...], [$(map(esc, style_opts)...)])
            $(esc(component)).opts[:_style_opts_cache][n] = v
        end
    end
end

function degreetoradians(angle::Float64)
    return angle/180 * pi
end

function create_graph(g::AbstractGraph, frames=:automatic; animate=false)
#    return JGraph()
end

"""
Use julia matrix
>>> matrix = [[1 1 1]; [2 2 2]]
"""
function graph_from_matrix(adjacency_matrix)
    sparse_mat = SparseMatrixCSC(adjacency_matrix)
    nv = size(sparse_mat)[1]
    g = transpose(sparse_mat) == sparse_mat ? SimpleGraph(nv) : SimpleDiGraph(nv)
    for ij in findall(!iszero, sparse_mat)
        add_edge!(g, ij[1], ij[2])
    end
    return g
end

"""
    check_frames_within(parent_frames, child_frames)

Checks if child frames are contained within parent frame to ensure parent-child ownership.
If not, throws a mild warning.
"""
function check_frames_within(parent_frames, child_frames, child_id="")
    child_frame_range = (child_frames isa Frames ? child_frames.frames : child_frames)
    parent_frame_range = (parent_frames isa Frames ? parent_frames.frames : parent_frames)
    if parent_frame_range isa UnitRange && child_frame_range isa UnitRange
        if !(child_frame_range ⊆ parent_frame_range)
            @warn "Child $(child_id) frame range $(child_frame_range) out of parent graph range $(parent_frame_range)"
        end
    end
end

#=
Debug utlities - frame amplify and frame marker with Javis liveviewer mode turned on.
=#


text:
# test null hypothesis cases

using Test
using IncrementalInference
# going to introduce two new constraint types
import IncrementalInference: getSample

mutable struct DevelopPriorNH <: IncrementalInference.FunctorSingletonNH
  x::Distribution
  nullhypothesis::Distributions.Categorical
end
function getSample(dpl::DevelopPriorNH, N::Int=1)
  return (reshape(rand(dpl.x, N),1,N), )
end

@testset "test null hypothesis singletons..." begin

global N  = 100
global fg = initfg()

global v1 = addVariable!(fg, :x1, ContinuousScalar, N=N)

global pr = DevelopPriorNH(Normal(10.0,1.0), Categorical([0.5;0.5]))
global f1 = addFactor!(fg,[:x1],pr, autoinit=true)

# ensureAllInitialized!(fg)
# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl",121)

global pts = evalFactor2(fg, f1, v1.label, N=N)

@test sum(abs.(pts .- 1.0) .< 5) > 30
@test sum(abs.(pts .- 10.0) .< 5) > 30

end

#


text:
using AmbientForcing
using OrdinaryDiffEq, Distributions

# Adopting the Robertson Example from the DifferentialEquation Docs
# https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Handling-Mass-Matrices
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  du[1] = -0.04 * y₁ + 1e4 * y₂ * y₃
  du[2] =  0.04 * y₁ - 1e4 * y₂ * y₃ - 3e7 * y₂^2
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end

# Creating the mass matrix M
# The last row depicts our constraint
M = [1. 0  0
     0  1. 0
     0  0  0];

# Setting up the DAE as an ODE in mass matrix form
ode_rober = ODEFunction(rober, mass_matrix = M)

# Choosing an inital condition which fulfills the constraint
u0 = [1.0,0.0,0.0];

# Acsessing the constraint equations g
g_rober = constraint_equations(ode_rober)
g_rober(u0)

# Randomly perturbing all variables
Frand = random_force(ode_rober, [0.0,1], Uniform)
z_new = ambient_forcing(ode_rober, u0, 2.0, Frand)

# g ≈ 0 means the constraint is fulfilled
g_rober(z_new)

# Only perturbing the second variable y₂
h = [0, 1 ,0]
z_new = ambient_forcing(ode_rober, u0, 2.0, h)

g_rober(z_new)


text:
#=
##############################
#
# mtl-Files
#
##############################

export WavefrontMtlMaterial,
    readMtlFile

type WavefrontMtlMaterial{T}
    name::String
    ambient::Vector3{T}
    specular::Vector3{T}
    diffuse::Vector3{T}
    transmission_filter::Vector3{T}
    illum::Int
    dissolve::T
    specular_exponent::T
    ambient_texture::String
    specular_texture::String
    diffuse_texture::String
    bump_map::String
end

function WavefrontMtlMaterial(colortype=Float64)
    return WavefrontMtlMaterial{colortype}(
                "",
                Vector3(zero(colortype)),
                Vector3(zero(colortype)),
                Vector3(zero(colortype)),
                Vector3(zero(colortype)),
                0,
                zero(colortype),
                zero(colortype),
                "",
                "",
                "",
                ""
            )
end

function readMtlFile(fn::String; colortype=Float64)
    str = open(fn,"r")
    mesh = readMtlFile(str, colortype=colortype)
    close(str)
    return mesh
end

function parseMtlColor(s::String, colortype=Float64)
    colorconvfnc = colortype == Float64 ? float64 : (colortype == Float32 ? float32 : error("colortype: ", colortype, " not supported"))

    line_parts = split(s)

    if length(line_parts) == 3 # r b g
        return Vector3{colortype}(colorconvfnc(line_parts[1]), colorconvfnc(line_parts[2]), colorconvfnc(line_parts[3]))
    elseif line_parts[1] == "spectral"
        println("WARNING Parsing Mtl-File: spectral color type not supported")
        return Vector3{colortype}(colorconvfnc(1.0), colorconvfnc(0.412), colorconvfnc(0.705))
    else
        println("WARNING Parsing Mtl-File: CIEXYZ color space or wrong color type. Not supported")
        return Vector3{colortype}(colorconvfnc(1.0), colorconvfnc(0.412), colorconvfnc(0.705))
    end
end

function parseMtlTextureMap(s::String)
    line_parts = split(s)

    if length(line_parts) == 1 # no options
        return line_parts[1]
    else
        println("WARNING Parsing Mtl-File: texture map options or invalid texutre map command. Not supported")
        return ""
    end
end

# supports the non-standard keywords: map_bump, map_refl
function readMtlFile(io::IO; colortype=Float64)
    colorconvfnc = colortype == Float64 ? float64 : (colortype == Float32 ? float32 : error("colortype: ", colortype, " not supported"))

    materials = WavefrontMtlMaterial{colortype}[]

    lineNumber = 1
    while !eof(io)
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(readline(io)))
        @assert is_valid_ascii(line)

        if !startswith(line, "#") && !isempty(line) && !iscntrl(line) #ignore comments
            line_parts = split(line)
            command = line_parts[1]
            remainder = length(line_parts) > 1 ? line[searchindex(line, line_parts[2]):end] : ""

            # new material
            if command == "newmtl"
                push!(materials, WavefrontMtlMaterial(colortype))
                materials[end].name = line_parts[2]
            # abmient
            elseif command == "Ka"
                materials[end].ambient = parseMtlColor(remainder, colortype)
            # diffuse
            elseif command == "Kd"
                materials[end].diffuse = parseMtlColor(remainder, colortype)
            # specular
            elseif command == "Ks"
                materials[end].specular = parseMtlColor(remainder, colortype)
            # transmission filter
            elseif command == "Tf"
                materials[end].transmission_filter = parseMtlColor(remainder, colortype)
            # illumination model
            elseif command == "illum"
                materials[end].illum = int(line_parts[2])
            # dissolve
            elseif command == "d"
                if line_parts[2] == "-halo"
                    println("WARNING Parsing Mtl-File: d -halo not supported")
                else
                    materials[end].dissolve = colorconvfnc(line_parts[2])
                end
            # specular exponent
            elseif command == "Ns"
                materials[end].specular_exponent = colorconvfnc(line_parts[2])
            # sharpness
            elseif command == "sharpness"
                println("WARNING Parsing Mtl-File: sharpness not supported")
            # optical density
            elseif command == "Ni"
                println("WARNING Parsing Mtl-File: optical density not supported")
            # ambient texture map
            elseif command == "map_Ka"
                materials[end].ambient_texture = parseMtlTextureMap(remainder)
            # diffuse texture map
            elseif command == "map_Kd"
                materials[end].diffuse_texture = parseMtlTextureMap(remainder)
            # specular texture map
            elseif command == "map_Ks"
                materials[end].specular_texture = parseMtlTextureMap(remainder)
            # specular exponent texture map
            elseif command == "map_Ns"
                println("WARNING Parsing Mtl-File: map_Ns not supported")
            # dissolve texture map
            elseif command == "map_d"
                println("WARNING Parsing Mtl-File: map_d not supported")
            # ???
            elseif command == "disp"
                println("WARNING Parsing Mtl-File: disp not supported")
            # ???
            elseif command == "decal"
                println("WARNING Parsing Mtl-File: decal not supported")
            # bump map
            elseif command == "bump" || command == "map_bump"
                materials[end].bump_map = parseMtlTextureMap(remainder)
            # ???
            elseif command == "refl" || command == "map_refl"
                println("WARNING Parsing Mtl-File: refl not supported")
            # unknown line
            else
                println("WARNING: Unknown line while parsing wavefront .mtl: '$line' (line $lineNumber)")
            end
        end

        # read next line
        lineNumber += 1
    end

    return materials
end

=#


text:
using PrescriptionUtils
using Test

n = 100
p = 5
X = rand(n, p)
y = rand(n)
T = rand(1:3, n)

@testset "Methods" for method in (:knn, :lasso, :randomforest)
  getcounterfactuals(X, y, T, method)
end


text:
# test on generic type

# We used to be able to handle small generic cases w/ GenericSVD alone (not GenericSchur)
# but the switch to GenericLinearAlgebra provides half-baked schur methods which trip us up.
using Pseudospectra, Test, GenericLinearAlgebra, GenericSchur

@testset "Generic(Big)" begin

    A = Matrix{BigFloat}(Pseudospectra.grcar(8))
    # ax is needed w/o eigvals(::Matrix{BigFloat})
    opts = Dict{Symbol,Any}(:ax => [-1,3,-3,3], :npts => 20)
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)

    # Just big enough to get to the inverse-Lanczos branch
    # this is a stunt, so don't waste time with usual npts.
    A = Matrix{BigFloat}(Pseudospectra.grcar(56))
    opts = Dict{Symbol,Any}(:ax => [-1,3,-3,3], :npts => 10)
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
end


text:
export seed!, categorical, choice
# random variables
for op in [:uniform, :normal]
    @eval begin
        export $op
        function $op(shape...;kwargs...)
            kwargs = jlargs(kwargs)
            if !haskey(kwargs, :dtype); kwargs[:dtype] = tf.float64; end
            tf.random.$op(shape;kwargs...)
        end
        $op() = squeeze($op(1, dtype=Float64))
        Docs.getdoc(::typeof($op)) = Docs.getdoc(tf.random.$op)
    end
end

"""
categorical(n::Union{PyObject, Integer}; kwargs...)

`kwargs` has a keyword argument `logits`, a 2-D Tensor with shape `[batch_size, num_classes]`.
Each slice `[i, :]` represents the unnormalized log-probabilities for all classes.
"""
function categorical(n::Union{PyObject, Integer}; kwargs...)
    flag = false
    kwargs = jlargs(kwargs)
    if !haskey(kwargs, :dtype); kwargs[:dtype] = tf.float64; end
    if !haskey(kwargs, :logits)
        kwargs[:logits] = tf.ones(n, dtype=kwargs[:dtype])
    end
    logits = kwargs[:logits]; delete!(kwargs, :logits)
    if length(size(logits))==1
        flag = true
        logits = convert_to_tensor(reshape(logits, 1, length(logits)))
    end
    out = tf.random.categorical(logits, n)
    if flag
        out = squeeze(out)
    end
    return out+1
end

"""
choice(inputs::Union{PyObject, Array}, n_samples::Union{PyObject, Integer};replace::Bool=false)

Choose `n_samples` samples from `inputs` with/without replacement.
"""
function choice(inputs::Union{PyObject, Array}, n_samples::Union{PyObject, Integer};
    replace::Bool=false)
    inputs = convert_to_tensor(inputs)
    if replace
        Idx = categorical(n_samples, logits=ones(size(inputs,1)))
        tf.gather(inputs, Idx-1)
    else
        dist = uniform(size(inputs,1))
        _, indices_to_keep = tf.nn.top_k(dist, n_samples)
        indices_to_keep_sorted = tf.sort(indices_to_keep)
        tf.gather(inputs,indices_to_keep)
    end
end

for op in [:Beta, :Bernoulli,:Gamma, :TruncatedNormal, :Binomial, :Cauchy,
        :Chi, :Chi2, :Exponential, :Gumbel,
        :HalfCauchy, :HalfNormal, :Horseshoe, :InverseGamma,
        :InverseGaussian, :Kumaraswamy, :Pareto, :SinhArcsinh,
        :StudentT, :VonMises,
        :Poisson]
    opname = Symbol(lowercase(string(op)))
    @eval begin
        export $opname
        function $opname(shape...;kwargs...)
            tfp = pyimport("tensorflow_probability")
            if !haskey(kwargs, :dtype); T = Float64;
            else T=kwargs[:dtype]; end
            kwargs = jlargs(kwargs)
            if haskey(kwargs, :dtype); delete!(kwargs, :dtype); end
            out = tfp.distributions.$op(;kwargs...).sample(shape...)
            cast(out, T)
        end
        $opname(;kwargs...) = squeeze($opname(1; dtype =Float64, kwargs...))
        # Docs.getdoc(::typeof($opname)) = Docs.getdoc(tfp.distributions.$op)
    end
end

function seed!(k::Int64)
    tf.random.set_random_seed(k)
end


text:
module SymPoly

using Primes

using MultivariatePolynomials
using DynamicPolynomials

using SymbolicUtils
using SymbolicUtils: istree, operation, arguments
using Symbolics: value, get_variables, expand_derivatives
using SymbolicUtils.Rewriters
# using SymbolicUtils.Code

using Symbolics

include("symbolics.jl")

export poly, sym

include("factoredpoly.jl")
include("wrapper.jl")

export wrap, unwrap, var, leading, cont, prim
export deg, derivative
export rationalize, unrationalize
export gcd_extended

include("roots.jl")

export solve_newton, find_roots

include("factorization.jl")

export FactoredPoly, factors, factor_schubert_kronecker
export decompose, factor, expand_frac, power, frac

include("modular.jl")
include("roundabout.jl")

include("roots_comb.jl")

include("integral.jl")

export integrate, generate_basis, test_integrals

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

function Primes.factor(p; method=:roots_SSP)
    if method == :schubert_kronecker
        return factor_schubert_kronecker(p)
    # elseif method == :roundabout
    #    return factor_roundabout(p)
    elseif method == :roots_comb
        return factor_roots_comb(p)
    elseif method == :roots_SSP
        return factor_roots_SSP(p)
    else
        error("undefined factorization method")
    end
end

Primes.factor(p, q) = factor_rational(p, q)
Primes.factor(r::RationalPoly) = factor_rational(r)

end # module


text:
"""
This function calculates the derating data for existing and new renewable generation
based on top 100 net-load hour methodology.
"""
function calculate_derating_data(simulation_dir::String,
                                active_projects::Vector{Project})

    cap_mkt_params = read_data(joinpath(simulation_dir, "markets_data", "capacity_mkt_param.csv"))

    renewable_existing = filter(p -> typeof(p) == RenewableGenEMIS{Existing}, active_projects)
    renewable_options = filter(p -> typeof(p) == RenewableGenEMIS{Option}, active_projects)

    zones = unique(get_zone.(get_tech.(renewable_existing)))
    types = unique(get_type.(get_tech.(renewable_existing)))

    load_n_vg_data =  read_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data", "load_n_vg_data.csv"))
    availability_data = read_data(joinpath(simulation_dir, "timeseries_data_files", "Availability", "DAY_AHEAD_availability.csv"))

    num_hours = DataFrames.nrow(load_n_vg_data)
    num_top_hours = cap_mkt_params.num_top_hours[1]

    existing_vg_power = zeros(num_hours)

    load = vec(sum(Matrix(load_n_vg_data[:, r"load"]), dims=2))

    for g in renewable_existing
        existing_vg_power += load_n_vg_data[!,get_name(g)]
    end

    net_load_df = load_n_vg_data[:, 1:4]
    net_load_df[:, "net_load"] = load - existing_vg_power

    net_load_sorted_df = DataFrames.sort(net_load_df, "net_load", rev = true)

    derating_factors = read_data(joinpath(simulation_dir, "markets_data", "derating_dict.csv"))

    type_zone_max_cap = Dict{String, Float64}()
    for zone in zones
        for type in types
            type_zone_id = "$(type)_$(zone)"
            type_zone_max_cap[type_zone_id] = 0.0
            net_load_df[:, "net_load_w/o_existing_$(type_zone_id)"] = net_load_df[:, "net_load"]
            for g in renewable_existing
                gen_name = get_name(g)
                tech = get_tech(g)
                if "$(get_type(tech))_$(get_zone(tech))" == type_zone_id
                    net_load_df[:, "net_load_w/o_existing_$(type_zone_id)"] += load_n_vg_data[:, gen_name]
                    type_zone_max_cap[type_zone_id] += get_maxcap(g)
                end
            end
        end
    end

    for zone in zones
        for type in types
            type_zone_id = "$(type)_$(zone)"
            gen_sorted_df = DataFrames.sort(net_load_df, "net_load_w/o_existing_$(type_zone_id)", rev = true)

            load_reduction = gen_sorted_df[1:num_top_hours, "net_load_w/o_existing_$(type_zone_id)"] - net_load_sorted_df[1:num_top_hours, "net_load"]
            derating_factors[:, "existing_$(type_zone_id)"] .= sum(load_reduction) / type_zone_max_cap[type_zone_id] / num_top_hours
        end
    end

    for g in renewable_options
        gen_name = get_name(g)
        tech = get_tech(g)
        type_zone_id = "$(get_type(tech))_$(get_zone(tech))"
        gen_cap = get_maxcap(g)

        net_load_df[:, "net_load_with_$(gen_name)"] =  net_load_df[:, "net_load"] - availability_data[:, "$(type_zone_id)"] * gen_cap
        gen_sorted_df = DataFrames.sort(net_load_df, "net_load_with_$(gen_name)", rev = true)

        load_reduction = net_load_sorted_df[1:num_top_hours, "net_load"] - gen_sorted_df[1:num_top_hours, "net_load_with_$(gen_name)"]
        derating_factors[:, "new_$(type_zone_id)"] .= sum(load_reduction) / gen_cap / num_top_hours
    end
    write_data(joinpath(simulation_dir, "markets_data"), "derating_dict.csv", derating_factors)
    return
end

"""
This function does nothing is project is not of ThermalGenEMIS, HydroGenEMIS, RenewableGenEMIS or BatteryEMIS type.
"""
function update_derating_factor!(project::P,
                               simulation_dir::String,
                               ) where P <: Project{<:BuildPhase}
    return
end

"""
This function updates the derating factors of ThermalGenEMIS and HydroGenEMIS projects.
"""
function update_derating_factor!(project::Union{ThermalGenEMIS{<:BuildPhase}, HydroGenEMIS{<:BuildPhase}},
                               simulation_dir::String,
                               )

    derating_data = read_data(joinpath(simulation_dir, "markets_data", "derating_dict.csv"))
    derating_factor = derating_data[1, get_type(get_tech(project))]
    for product in get_products(project)
        set_derating!(product, derating_factor)
    end
    return
end

"""
This function updates the derating factors of existing RenewableGenEMIS projects.
"""
function update_derating_factor!(project::RenewableGenEMIS{Existing},
                               simulation_dir::String,
                               )

    derating_data = read_data(joinpath(simulation_dir, "markets_data", "derating_dict.csv"))
    name = get_name(project)
    tech = get_tech(project)
    type_zone_id = "$(get_type(tech))_$(get_zone(tech))"

    if  in("existing_$(type_zone_id)", names(derating_data))
        derating_factor = derating_data[1, "existing_$(type_zone_id)"]
    else
        error("Derating data not found")
    end

    for product in get_products(project)
        set_derating!(product, derating_factor)
    end

    return
end

"""
This function updates the derating factors of new RenewableGenEMIS projects.
"""
function update_derating_factor!(project::RenewableGenEMIS{<:BuildPhase},
                               simulation_dir::String,
                               )

    derating_data = read_data(joinpath(simulation_dir, "markets_data", "derating_dict.csv"))
    name = get_name(project)
    tech = get_tech(project)
    type_zone_id = "$(get_type(tech))_$(get_zone(tech))"

    if in("new_$(type_zone_id)", names(derating_data))
        derating_factor = derating_data[1, "new_$(type_zone_id)"]
    else
        error("Derating data not found")
    end

    for product in get_products(project)
        set_derating!(product, derating_factor)
    end

    return
end

"""
This function updates the derating factors of BatteryEMIS projects.
"""
function update_derating_factor!(project::BatteryEMIS{<:BuildPhase},
                               simulation_dir::String,
                               )
    tech = get_tech(project)
    duration = Int(get_storage_capacity(tech)[:max] / get_maxcap(project))
    project_type = "$(get_type(tech))_$(duration)"

    derating_data = read_data(joinpath(simulation_dir, "markets_data", "derating_dict.csv"))
    derating_factor = derating_data[1, project_type]
    for product in get_products(project)
        set_derating!(product, derating_factor)
    end
    return
end

"""
This function updates the derating factors of all active projects in the simulation.
"""
function update_simulation_derating_data!(simulation::Union{AgentSimulation, AgentSimulationData})
    data_dir = get_data_dir(get_case(simulation))
    active_projects = get_activeprojects(simulation)

    calculate_derating_data(data_dir, active_projects)
    for project in active_projects
        update_derating_factor!(project, data_dir)
    end
    return
end