text:
# Loosely based on MiniZinc's fzn_lex_lesseq_float for MIP:
# https://github.com/MiniZinc/libminizinc/blob/master/share/minizinc/linear/fzn_lex_lesseq_float.mzn
# Major difference: that reformulation is only for two vectors,
# this implementation is for an arbitrary number of vectors (there is no
# dedicated set for just two vectors).

"""
Bridges `CP.Strictly{CP.LexicographicallyLessThan}` to indicators.
"""
struct StrictlyLexicographicallyLessThan2IndicatorBridge{T} <: MOIBC.AbstractBridge
    vars_eq::Matrix{MOI.VariableIndex}
    vars_lt::Matrix{MOI.VariableIndex}
    vars_eq_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    vars_lt_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_one_lt::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_move::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_indic_eq::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}}
    cons_indic_lt::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}}}
end

function MOIBC.bridge_constraint(
    ::Type{StrictlyLexicographicallyLessThan2IndicatorBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Strictly{CP.LexicographicallyLessThan, T},
) where {T}
    return MOIBC.bridge_constraint(
        StrictlyLexicographicallyLessThan2IndicatorBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{StrictlyLexicographicallyLessThan2IndicatorBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Strictly{CP.LexicographicallyLessThan, T},
) where {T}
    @assert s.set.row_dim >= 1
    @assert s.set.column_dim >= 2

    vars_eq = Matrix{MOI.VariableIndex}(undef, s.set.column_dim - 1, s.set.row_dim)
    vars_lt = Matrix{MOI.VariableIndex}(undef, s.set.column_dim - 1, s.set.row_dim)
    vars_eq_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.set.column_dim - 1, s.set.row_dim)
    vars_lt_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.set.column_dim - 1, s.set.row_dim)
    cons_one_lt = Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}(undef, s.set.column_dim - 1)
    cons_move = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}(undef, s.set.column_dim - 1, s.set.row_dim - 1)
    cons_indic_eq = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}}(undef, s.set.column_dim - 1, s.set.row_dim)
    cons_indic_lt = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}}}(undef, s.set.column_dim - 1, s.set.row_dim)

    f_scalars = MOIU.scalarize(f)
    f_matrix = reshape(f_scalars, s.set.column_dim, s.set.row_dim)
    # First index: column; second index: row. Rationale: columns are sorted
    # lexicographically, not rows.

    # Add the constraint between the columns i and i+1.
    for i in 1:(s.set.column_dim - 1)
        vars_eq[i, :], vars_eq_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.set.row_dim])
        vars_lt[i, :], vars_lt_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.set.row_dim])

        cons_one_lt[i] = MOI.add_constraint(
            model,
            sum(one(T) .* vars_lt[i, :]),
            MOI.EqualTo(one(T)),
        )

        cons_move[i, :] = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
            MOI.add_constraint(
                model,
                one(T) * vars_eq[i, j - 1] - one(T) * vars_eq[i, j] - one(T) * vars_lt[i, j],
                MOI.EqualTo(zero(T)),
            )
            for j in 2:s.set.row_dim
        ]

        cons_indic_eq[i, :] = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}[
            MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_eq[i, j],
                        f_matrix[i, j] - f_matrix[i + 1, j],
                    ]
                ),
                MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T))),
            )
            for j in 1:s.set.row_dim
        ]

        cons_indic_lt[i, :] = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}}[
            MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_lt[i, j],
                        f_matrix[i, j] - f_matrix[i + 1, j],
                    ]
                ),
                MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(zero(T))),
            )
            for j in 1:s.set.row_dim
        ]
    end

    return StrictlyLexicographicallyLessThan2IndicatorBridge(vars_eq, vars_lt, vars_eq_bin, vars_lt_bin, cons_one_lt, cons_move, cons_indic_eq, cons_indic_lt)
end

function MOI.supports_constraint(
    ::Type{StrictlyLexicographicallyLessThan2IndicatorBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Strictly{CP.LexicographicallyLessThan, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{StrictlyLexicographicallyLessThan2IndicatorBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{StrictlyLexicographicallyLessThan2IndicatorBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}),
    ]
end

function MOI.get(b::StrictlyLexicographicallyLessThan2IndicatorBridge{T}, ::MOI.NumberOfVariables) where {T}
    return length(b.vars_eq) + length(b.vars_lt)
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_eq_bin) + length(b.vars_lt_bin)
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons_move) + length(b.cons_one_lt)
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_indic_eq)
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}},
    },
) where {T}
    return length(b.cons_indic_lt)
end

function MOI.get(b::StrictlyLexicographicallyLessThan2IndicatorBridge{T}, ::MOI.ListOfVariableIndices) where {T}
    return vcat(vec(b.vars_eq), vec(b.vars_lt))
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return vcat(vec(b.vars_eq_bin), vec(b.vars_lt_bin))
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return vcat(b.cons_one_lt, vec(b.cons_move))
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T}
    return vec(b.cons_indic_eq)
end

function MOI.get(
    b::StrictlyLexicographicallyLessThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}},
    },
) where {T}
    return vec(b.cons_indic_lt)
end


text:
"""
    make_time_series(data, number_of_periods; loadprofile, genprofile)

Make a time series dict from profile matrices, to be used in `make_multinetwork`.

`number_of_periods` is by default the number of networks specified in `data`.
Profile matrices must have `number_of_periods` rows and one column for each component (load
or generator).

# Arguments
- `data`: a multinetwork data dictionary;
- `number_of_periods = dim_length(data)`;
- `loadprofile = ones(number_of_periods,length(data["load"]))`;
- `genprofile = ones(number_of_periods,length(data["gen"])))`.
"""
function make_time_series(data::Dict{String,Any}, number_of_periods::Int = dim_length(data); loadprofile = ones(number_of_periods,length(data["load"])), genprofile = ones(number_of_periods,length(data["gen"])))
    if size(loadprofile) ≠ (number_of_periods, length(data["load"]))
        right_size = (number_of_periods, length(data["load"]))
        Memento.error(_LOGGER, "Size of loadprofile matrix must be $right_size, found $(size(loadprofile)) instead.")
    end
    if size(genprofile) ≠ (number_of_periods, length(data["gen"]))
        right_size = (number_of_periods, length(data["gen"]))
        Memento.error(_LOGGER, "Size of genprofile matrix must be $right_size, found $(size(genprofile)) instead.")
    end
    return Dict{String,Any}(
        "load" => Dict{String,Any}(l => Dict("pd" => load["pd"] .* loadprofile[:, parse(Int, l)]) for (l,load) in data["load"]),
        "gen" => Dict{String,Any}(g => Dict("pmax" => gen["pmax"] .* genprofile[:, parse(Int, g)]) for (g,gen) in data["gen"]),
    )
end

ext: jl


text:
include("tests_ebayes_samples.jl")
include("tests_example_priors.jl")
include("tests_sure.jl")
include("tests_ebcf.jl")


text:
"""
    A = JopConvolve(dom, rng, h [, optional parameters])

`A` is an n-dimension convolution (using the filter `h::Array`) operator with domain
`dom::JetSpace` and range `rng::JetSpace`, and with the following optional named arguments:

* `x0` is a tuple defining the origin of the upper-left corner of `h`

* `dx` is a tuple defining the spacing along each dimension of `h`

# Examples

## 1D, causal
```julia
A = JopConvolve(JetSpace(Float64,128), JetSpace(Float64,128), rand(32))
m = zeros(domain(A))
m[64] = 1.0
d = A*m
```

## 1D, zero-phase
```julia
A = JopConvolve(JetSpace(Float64,128), JetSpace(Float64,128), rand(32), dx=(1.0,), x0=(-16.0,))
m = zeros(domain(A))
m[64] = 1.0
d = A*m
```

## 2D, zero-phase
```julia
A = JopConvolve(JetSpace(Float64,128,128), JetSpace(Float64,128,128), rand(32,32), dx=(1.0,1.0), x0=(-16.0,-16.0))
m = zeros(domain(A))
m[64,64] = 1.0
d = A*m
```

# Notes
* It is often the case that the domain and range of the convolution operator are the same.  For this use-case, we provide
a convenience method for construction the operator:

    A = JopConvolve(spc, h [, optional parameters])

where `spc::JetSpace` and is used for both `dom` and `rng`.

* Since smoothing is a common use-case for JopConvolve, we provide a convenience method for creating `A` specific
to n-dimensional smoothing:

    A = JopConvolve(spc [, optional arguments])

where the optional arguments and their default values are:

* `smoother=:gaussian` choose between `:gaussian`, `:triang` and `:rect`

* `n=(128,)` choose the size of the smoothing window in each dimension.  If `length(n)=1`, then we assume a square window.

* `sigma=(0.5,)` for a gaussian window choose the shape of the window.  If `length(sigma)=1`, then we assume the same shape in each dimension.

# 2D Smoothing Example
```julia
P = JopPad(JetSpace(Float64,256,256), -10:256+11, -10:256+11, extend=true)
S = JopConvolve(range(P), smoother=:rect, n=(1,1))
R = JopPad(JetSpace(Float64,256,256), -10:256+11, -10:256+11, extend=false)
m = rand(domain(P))
d = R'∘S∘P*m
```
"""
function JopConvolve(dom::JetSpace{T,N}, rng::JetSpace{T,N}, h::AbstractArray{T,N}; x0=(0.0,), dx=(1.0,)) where {T<:Real,N}
    x0 = (length(x0) == 1 ? ntuple(i->x0[1], N) : x0)::NTuple{N,T}
    dx = (length(dx) == 1 ? ntuple(i->dx[1], N) : dx)::NTuple{N,T}

    length(x0) == N || error("Expected length(x0)=1 or length(x0)=ndims(h), got length(x0)=$(length(x0))")
    length(dx) == N || error("Expected length(dx)=1 or length(dx)=ndims(h), got length(dx)=$(length(dx))")

    nker = size(h)
    ndom = size(dom)
    nrng = size(rng)
    ntot = ntuple(idim->nextprod((2,3,5,7), nker[idim] + ndom[idim] - 1), N)::NTuple{N,Int}

    ntot < nrng && error("Range must be greater than or equal size(model)+size(h)-1")

    nfft = ntuple(idim->idim == 1 ? div(ntot[idim],2) + 1 : ntot[idim], N)::NTuple{N,Int}
    dfft = ntuple(idim->2*T(pi)/dx[idim]/ntot[idim], N)::NTuple{N,T}

    k = Vector{Vector{T}}(undef, N)
    k[1] = dfft[1]*collect(0:nfft[1]-1)
    for idim = 2:N
        nfft_p = div(nfft[idim]+1, 2)
        nfft_n = div(nfft[idim], 2)
        k[idim] = dfft[idim]*[collect(0:nfft_p-1) ; collect(-nfft_n:1:-1)]
    end

    hpad = zeros(T, ntot)
    for i in CartesianIndices(nker)
        hpad[i] = h[i]
    end

    H = rfft(hpad)

    phaseshift!(H, k, x0)

    JopLn(dom = dom, rng = rng, df! = JopConvolve_df!, df′! = JopConvolve_df′!,
        s = (mpad=Array{T,N}(undef,ntot), dpad=Array{T,N}(undef,ntot), H=H))
end
JopConvolve(spc::JetSpace, h::Array; x0=(0.0,), dx=(1.0,)) = JopConvolve(spc, spc, h, x0=x0, dx=dx)

# Convenience constructor for n-dimensional smoothing
function JopConvolve(spc::JetSpace{T,N}; smoother=:gaussian, n=(128,), sigma=(0.5,)) where {T,N}
    if smoother != :gaussian && smoother != :triang && smoother != :rect
        error("expected smoother=:gaussian, smoother=:triang or smoother=:rect, got smoother=$(smoother)")
    end
    n = length(n) == 1 ? ntuple(idim->2*n[1]+1, N) : ntuple(idim->2*n[idim]+1, N)
    sigma = length(sigma) == 1 ? ntuple(idim->sigma[1], N) : sigma
    dx = ntuple(idim->one(T), N)
    x0 = ntuple(idim->-one(T)*div(n[idim]-1,2), N)

    w = Vector{Vector{T}}(undef, N)
    for idim = 1:N
        if smoother == :gaussian
            w[idim] = gaussian(n[idim],sigma[idim])
        elseif smoother == :triang
            w[idim] = triang(n[idim])
        elseif smoother == :rect
            w[idim] = rect(n[idim])
        end
    end

    h = filter!(ones(T, n), w)

    JopConvolve(spc, spc, h, x0=x0, dx=dx)
end

export JopConvolve

function phaseshift!(H::Array{T,N}, k, x0) where {T,N}
    for i in CartesianIndices(H)
        f = zero(T)
        for idim = 1:N
            f += k[idim][i[idim]]*x0[idim]
        end
        H[i] *= @fastmath exp(-im*f)
    end
end

function filter!(h::Array{T,N}, w) where {T,N}
    for i in CartesianIndices(h)
        h[i] = 1
        for idim = 1:N
            h[i] *= w[idim][i[idim]]
        end
    end
    h ./= sum(h)
    h
end

function JopConvolve_df!(d::AbstractArray{T,N}, m::AbstractArray{T,N}; mpad, dpad, H, kwargs...) where {T<:Real,N}
    ndom = size(m)
    nrng = size(d)
    nfft = size(H)
    ntot = size(mpad)
    mpad .= 0
    for i in CartesianIndices(ndom)
        mpad[i] = m[i]
    end
    M = rfft(mpad)
    D = Array{Complex{T}}(undef, nfft)
    for i = 1:length(D)
        D[i] = H[i]*M[i]
    end
    dd = irfft(D, ntot[1])
    for i in CartesianIndices(nrng)
        d[i] = dd[i]
    end
    d
end

function JopConvolve_df′!(m::AbstractArray{T,N}, d::AbstractArray{T,N}; mpad, dpad, H, kwargs...) where {T<:Real,N}
    ndom = size(m)
    nrng = size(d)
    nfft = size(H)
    ntot = size(mpad)
    dpad .= 0
    for i in CartesianIndices(nrng)
        dpad[i] = d[i]
    end
    D = rfft(dpad)
    M = Array{Complex{T}}(undef, nfft)
    for i = 1:length(M)
        M[i] = conj(H[i]) * D[i]
    end
    mm = irfft(M, ntot[1])
    for i in CartesianIndices(ndom)
        m[i] = mm[i]
    end
    m
end


text:
using ExoCalc
using Test

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

max_stars_repo_name: icweaver/exocalc

max_issues_repo_name: icweaver/exocalc

max_forks_repo_name: icweaver/exocalc


text:
using Revise
using Pkg
Pkg.activate(".")
using MCPhylo
using Random
using LinearAlgebra

Random.seed!(42)

s = rand(InverseWishart(10, diagm(ones(3))))
#s = rand(Exponential())
mu = [0.1,0.4,0.3,0.15,0.05]
μ = randn(3)
mu1 = mu + rand(5)
mu1 /= sum(mu1)
mu1 = rand(Dirichlet(5,1))
my_data = Dict{Symbol, Any}(
  :df => collect(transpose(rand(Multinomial(1, mu1), 500)))
  #collect(transpose(rand(MultivariateNormal(μ, s), 500)))
);

# model setup
model =  Model(
    df = Stochastic(2, (mu1) ->  MultivariateDistribution[Multinomial(1, mu1) for i in 1:500], false),
    #MultivariateDistribution[MultivariateNormal(mu1, cm) for i in 1:500], false),
    mu1 = Stochastic(1, () -> Dirichlet(5,1.0), true),
    cm = Stochastic(2, () -> InverseWishart(10, diagm(ones(3))), false)
     )
# intial model values
inits = [ Dict{Symbol, Union{Any, Real}}(
    :mu1 => rand(Dirichlet(5, 1.0)),
    :cm => rand(InverseWishart(3, diagm(ones(3)))),
    :df => my_data[:df],

    ) for i in 1:2
    ]

scheme = [#Slice([:mu1, :cm], 1.0),
          #RWM(:mu1, 1.0)
          SliceSimplex(:mu1, scale = 1)
          #NUTS([:mu1, :cm], tree_depth=10, target=0.7)
          #HMC([:mu1, :cm], 0.001, 10)
]

setsamplers!(model, scheme);

# do the mcmc simmulation. if trees=true the trees are stored and can later be
# flushed ot a file output.

sim = mcmc(model, my_data, inits, 25000, burnin=12500,thin=10, chains=2, trees=true)

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
module Xtal

using StaticArrays
using LinearAlgebra
using Printf
using ComputedFieldTypes
using FFTW

const ELEMENTS =
(
    "H",  "He", "Li", "Be", "B",  "C",  "N",  "O",  "F",  "Ne",
    "Na", "Mg", "Al", "Si", "P",  "S",  "Cl", "Ar", "K",  "Ca",
    "Sc", "Ti", "V",  "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
    "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y",  "Zr",
    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
    "Sb", "Te", "I",  "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
    "Lu", "Hf", "Ta", "W",  "Re", "Os", "Ir", "Pt", "Au", "Hg",
    "Tl", "Pb", "Bi", "Po", "As", "Rn", "Fr", "Ra", "Ac", "Th",
    "Pa", "U",  "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
    "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
    "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"
)

# TODO: can we modify this to make it unnecessary?
const ELEMENT_LOOKUP = Dict{String, Int}([(ELEMENTS[n] => n) for n in 1:length(ELEMENTS)])

# Matrices used to reduce conventional cells to corresponding primitive cells
# These can be accessed by calling the centering letter as a symbol
# So if you want the face-centered reduction matrix, use REDUCTION_MATRIX_3D[:F]
const REDUCTION_MATRIX_3D =
(
    # Transpose these matrices in your head, please
    A = SMatrix{3,3,Float64}( 1,    0,    0,
                              0,  1/2,  1/2,
                              0, -1/2,  1/2),
    B = SMatrix{3,3,Float64}( 1/2,  0,  1/2,
                                0,  1,    0,
                             -1/2,  0,  1/2),
    C = SMatrix{3,3,Float64}( 1/2,  1/2,  0,
                             -1/2,  1/2,  0,
                                0,    0,  1),
    F = SMatrix{3,3,Float64}( 1/2,  1/2,    0,
                              1/2,    0,  1/2,
                                0,  1/2,  1/2),
    I = SMatrix{3,3,Float64}(-1/2,  1/2,  1/2,
                              1/2, -1/2,  1/2,
                              1/2,  1/2, -1/2),
    P = SMatrix{3,3,Float64}(LinearAlgebra.I)
)

"""
    Xtal.TOL_DEF

Default tolerance for discrepancies in floating point values.
"""
const TOL_DEF = 1e-8

"""
    Xtal.BOHR2ANG

Converts lengths in bohr to angstrom
"""
const BOHR2ANG = 0.52917720859

#= Units of these two constants should be energy^-1 * length^-2
"""
    Xtal.CABINIT

The value of 2m_e/ħ^2 as used in abinit (hartree^-1 * angstrom^-2).
"""
const CABINIT = 7.142129652186264
=#

"""
    Xtal.CVASP

The value of 2m_e/ħ^2 as used in VASP (eV^-1 * angstrom^-2, but not exactly).
"""
const CVASP = 0.262465831
# the correct value is actually 0.26246842360751754

"""
    Xtal._allsame(itr)

Returns `true` if all the elements of an iterator are identical.
"""
_allsame(itr) = all(x -> x == first(itr), itr)

"""
    Xtal._selfdot(v)

Computes the dot product of a vector with itself.
"""
_selfdot(v) = dot(v,v)

# Methods used in vector operations that go beyond what's available in LinearAlgebra
include("vectors.jl")
# Abstract types used in type tree
include("types.jl")
export AbstractLattice, AbstractCrystal, AbstractCrystalData, AbstractRealSpaceData,
       AbstractReciprocalSpaceData, AbstractKPoints, AbstractDensityOfStates, AbstractPotential,
       AbstractPseudopotential
# Methods and structs for working with crystal lattices
include("lattices.jl")
export BasisVectors, RealLattice, ReciprocalLattice
export dual, prim, conv, cell_lengths, cell_volume, lengths, volume, lattice2D, lattice3D,
       maxHKLindex
# Methods and structs for working with atomic positions
include("atoms.jl")
export AtomPosition, AtomList
export atomname, atomicno, coord, natom, basis, cartesian, reduce_coords
# Methods and structs for working with crystal data
include("crystals.jl")
export Crystal, CrystalWithDatasets
export natom_template, data, prim, conv, volume
# Methods and structs for working with different types of data associated with crystals
include("data.jl")
export RealSpaceDataGrid, KPointGrid, KPointList, BandAtKPoint, BandStructure, HKLData,
       ReciprocalWavefunction, DensityOfStates, ProjectedDensityOfStates, FatBands, AtomicData,
       SphericalComponents
export grid, gridsize, volume, voxelsize, integrate, fft, nkpt, nband, bounds, fermi, smear,
       energies, nelectrons
include("potentials.jl")
export XCFunctional, HGHPseudopotential, zatom, zion
# Methods and structs for working with different file formats
include("filetypes.jl")
export readXYZ, readXSF3D, writeXSF, readCPcoeff, readCPgeo, readCPcell
export read_abinit_density, read_abinit_potential, read_abinit_wavefunction, readHGH
export readWAVECAR, readDOSCAR, readPROCAR
# Show methods for pretty printing this module's structs
include("show.jl")
# Precompilation directives
include("precompile.jl")

end # end of module


text:
#
# Description
#
#   Functions to configure the MGEO.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                              Public Functions
################################################################################

"""
    function conf_design_vars(T, Nv::Int64, bits::Int64, min::Number, max::Number)

This function configures `Nv` design variables in which all of them will have
the same number of bits `bits`, and maximum and minimum values `min` and `max`.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    function conf_design_vars(Val{:MGEO_Real}, Nv::Int64, min::Number, max::Number, σ::Number)

This function configures `Nv` design variables in which all of them will have
the same maximum and minimum values `min` and `max`, and the same standard
deviation `σ` used when modifying the variables. This is used for MGEO Real
algorithm.

# Args

* `T`: (OPTIONAL) Type of algorithm (**Default**: `MGEO_Canonical()`).
* `Nv`: The number of design variables.
* `bits`: The number of bits for each variable.
* `min`: The minimum value for all the design variables.
* `max`: The maximum value for all the design variables.

* `σ`: Standard deviation used when modifying a variable in MGEO Real.

# Returns

An array with the configured design variables.

"""
conf_design_vars(Nv::Int64, bits::Int64, min::Number, max::Number) =
    conf_design_vars(Val{:MGEO_Canonical}, Nv, bits, min, max)

function conf_design_vars(T::Type{Val{:MGEO_Canonical}},
                          Nv::Int64,
                          bits::Int64,
                          min::Number,
                          max::Number)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && @error("The number of design variables must be higher than 0.")

    # Check if the number of bits are greater than 0.
    ( bits <= 0 ) && @error("The number of bits must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && @error("The minimum value must be less than the maximum value for a variable.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Canonical}(undef, Nv)

    # Full scale for each variable.
    fullscale = (1 << bits) - 1

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Canonical(bits,
                                                        min,
                                                        max,
                                                        fullscale,
                                                        (i-1)*bits+1,
                                                        "Var. $i")
    end

    SVector{Nv, Design_Variable_MGEO_Canonical}(design_vars)
end

function conf_design_vars(T::Type{Val{:MGEO_Var}},
                          Nv::Int64,
                          bits::Int64,
                          min::Number,
                          max::Number)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && @error("The number of design variables must be higher than 0.")

    # Check if the number of bits are greater than 0.
    ( bits <= 0 ) && @error("The number of bits must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && @error("The minimum value must be less than the maximum value for a variable.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Var}(undef, Nv)

    # Full scale for each variable.
    fullscale = (1 << bits) - 1

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Var(bits,
                                                  min,
                                                  max,
                                                  fullscale,
                                                  (i-1)*bits+1,
                                                  "Var. $i")
    end

    SVector{Nv, Design_Variable_MGEO_Var}(design_vars)
end

function conf_design_vars(T::Type{Val{:MGEO_Real}},
                          Nv::Int64,
                          min::Number,
                          max::Number,
                          σ::Number)

    # Check if the number of design variables is greater than 1.
    ( Nv <= 0 ) && @error("The number of design variables must be higher than 0.")

    # Check if the minimum value is smaller than maximum value.
    ( min >= max ) && @error("The minimum value must be less than the maximum value for a variable.")

    # Check if standard deviation σ is higher than 0.
    ( σ <= 0 ) && @error("The standard deviation must be higher than 0.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Real}(undef, Nv)

    @inbounds for i=1:Nv
        design_vars[i] = Design_Variable_MGEO_Real(min, max, σ, "Var. $i")
    end

    SVector{Nv, Design_Variable_MGEO_Real}(design_vars)
end

"""
    function conf_design_vars(T, bits::Vector{Int64}, min::Vector{T1}, max::Vector{T2}, var_names::Vector{String}) where {T1<:Number, T2<:Number}

This function configures the design variables specifying for each one the number
of bits `bits`, the minimum `min` and maximum `max` values, and the variable
names `var_names`.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    function conf_design_vars(T::Type{Val{:MGEO_Real}}, min::Vector{T1}, max::Vector{T2}, σ::Vector{T3}, var_names::Vector{String}) where

This function configures the design variables specifying for each one the
minimum `min` and maximum `max` values, the standard deviations `σ` used when
modifying the variable, and the variable names `var_names`.

# Args

* `T`: (OPTIONAL) Type of algorithm (**Default**: `MGEO_Canonical()`).
* `bits`: List containing the number of bits for each design variable.
* `min`: List containing the minimum of each design variable.
* `max`: List containing the maximum of each design variable.
* `var_names`: List containing the name of each design variable.

* `σ`: List containing the standard deviations used when modifying a variable in
       MGEO Real.

# Returns

An array with the configured design variables.

"""
conf_design_vars(bits::Vector{Int64},
                 min::Vector{T1},
                 max::Vector{T2},
                 var_names::Vector{String}) where
    {T1<:Number, T2<:Number} =
    conf_design_vars(Val{:MGEO_Canonical}, bits, min, max, var_names)

function conf_design_vars(T::Type{Val{:MGEO_Canonical}},
                          bits::Vector{Int64},
                          min::Vector{T1},
                          max::Vector{T2},
                          var_names::Vector{String}) where
 	{T1<:Number, T2<:Number}

    # Check if the size of arrays is correct.
    Nv = length(bits)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(var_names) != Nv ) ) &&
    @error("The size of the vectors does not match.")

    # Number of bits.
    num_bits = 0

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Canonical}(undef, Nv)

    for i = 1:Nv
        # Check if the number of bits are larger than 0.
        ( bits[i] <= 0 ) && @error("The number of bits must be higher than 0.")

        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && @error("The minimum value must be less than the maximum value for a variable.")

        design_vars[i] = Design_Variable_MGEO_Canonical(bits[i],
                                                        min[i],
                                                        max[i],
                                                        (UInt64(1) << bits[i]) - 1,
                                                        num_bits+1,
                                                        var_names[i])
        num_bits += bits[i]
    end

    SVector{Nv, Design_Variable_MGEO_Canonical}(design_vars)
end

function conf_design_vars(T::Type{Val{:MGEO_Var}},
                          bits::Vector{Int64},
                          min::Vector{T1},
                          max::Vector{T2},
                          var_names::Vector{String}) where
 	{T1<:Number, T2<:Number}

    # Check if the size of arrays is correct.
    Nv = length(bits)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(var_names) != Nv ) ) &&
    @error("The size of the vectors does not match.")

    # Number of bits.
    num_bits = 0

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Var}(undef, Nv)

    for i = 1:Nv
        # Check if the number of bits are larger than 0.
        ( bits[i] <= 0 ) && @error("The number of bits must be higher than 0.")

        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && @error("The minimum value must be less than the maximum value for a variable.")

        design_vars[i] = Design_Variable_MGEO_Var(bits[i],
                                                  min[i],
                                                  max[i],
                                                  (UInt64(1) << bits[i]) - 1,
                                                  num_bits+1,
                                                  var_names[i])
        num_bits += bits[i]
    end

    SVector{Nv, Design_Variable_MGEO_Var}(design_vars)
end

function conf_design_vars(T::Type{Val{:MGEO_Real}},
                          min::Vector{T1},
                          max::Vector{T2},
                          σ::Vector{T3},
                          var_names::Vector{String}) where
 	{T1<:Number, T2<:Number, T3<:Number}

    # Check if the size of arrays is correct.
    Nv = length(min)

    ( ( length(min)       != Nv ) ||
      ( length(max)       != Nv ) ||
      ( length(σ)         != Nv ) ||
      ( length(var_names) != Nv ) ) &&
    @error("The size of the vectors does not match.")

    # Create the array of design variables.
    design_vars = Vector{Design_Variable_MGEO_Real}(undef, Nv)

    for i = 1:Nv
        # Check if the min value is smaller than max.
        ( min[i] >= max[i] ) && @error("The minimum value must be less than the maximum value for a variable.")

        # Check if standard deviation σ is higher than 0.
        ( σ[i] <= 0 ) && @error("The standard deviation must be higher than 0.")

        design_vars[i] = Design_Variable_MGEO_Real(min[i],
                                                   max[i],
                                                   σ[i],
                                                   var_names[i])
    end

    SVector{Nv, Design_Variable_MGEO_Real}(design_vars)
end

"""
    function conf_design_vars(bits::Vector{Int64}, min::Vector{T1}, max::Vector{T2}) where {T1<:Number, T2<:Number}

This function configures the design variables specifying for each one the number
of bits `bits`, and the minimum `min` and maximum `max` values. The variable
names will be selected as `Var. 1`, `Var. 2`, etc.

The parameter `T` selects which MGEO algorithm will be used and it can be
`MGEO_Canonical()` (**Default**) or `MGEO_Var()`.

    function conf_design_vars(T::Type{Val{:MGEO_Real}}, min::Vector{T1}, max::Vector{T2}, σ::Vector{T3}) where {T1<:Number, T2<:Number, T3<:Number}

This function configures the design variables specifying for each one the
minimum `min` and maximum `max` values, and the standard deviations `σ` used
when modifying the variable. The variable names will be selected as `Var. 1`,
`Var. 2`, etc.

# Args

* `T`: (OPTIONAL) Type of algorithm (**Default**: `MGEO_Canonical()`).
* `bits`: List containing the number of bits for each design variable.
* `min`: List containing the minimum of each design variable.
* `max`: List containing the maximum of each design variable.

* `σ`: List containing the standard deviations used when modifying a variable in
       MGEO Real.

# Returns

An array with the configured design variables.

"""
conf_design_vars(bits::Vector{Int64},
                 min::Vector{T1},
                 max::Vector{T2}) where {T1<:Number, T2<:Number} =
    conf_design_vars(Val{:MGEO_Canonical}, bits, min, max)

function conf_design_vars(T::Union{Type{Val{:MGEO_Canonical}},
								   Type{Val{:MGEO_Var}}},
                          bits::Vector{Int64},
                          min::Vector{T1},
                          max::Vector{T2}) where {T1<:Number, T2<:Number}

    # Create an array with variable names.
    var_names = Array(String, length(bits))

    for i=1:length(varNames)
        varNames[i] = "Var. $i"
    end

    conf_design_vars(T, bits, min, max, var_names)
end

function conf_design_vars(T::Type{Val{:MGEO_Real}},
                          min::Vector{T1},
                          max::Vector{T2},
                          σ::Vector{T3}) where
    {T1<:Number, T2<:Number, T3<:Number}

    # Create an array with variable names.
    var_names = Array(String, length(bits))

    for i=1:length(varNames)
        varNames[i] = "Var. $i"
    end

    conf_design_vars(T, min, max, σ, var_names)
end

"""
    function conf_mgeo(Nf::Int64, τ::Float64, ngen_max::Int64, run_max::Int64, design_vars::SVector{Nv, T}, mgeo_eps::Float64 = 1e-10) where {Nv,T}

Configure the MGEO.

The MGEO algorithm that will be used will be the same configured in
`design_vars`.

# Args

* `Nf`: Number of objective functions.
* `τ`: Parameter τ of MGEO.
* `ngen_max`: Maximum number of generations.
* `run_max`: Maximum number of independent runs.
* `design_vars`: Configuration of the design variables (see `conf_design_vars`).
* `mgeo_eps`: (OPTIONAL) Epsilon to compare objective functions (**Default**:
              `1e-10`)

# Returns

An instance of the structure `MGEO_Structure` with the MGEO configuration.

"""
function conf_mgeo(Nf::Int64,
                   τ::Float64,
                   ngen_max::Int64,
                   run_max::Int64,
                   design_vars::SVector{Nv, T},
                   mgeo_eps::Float64 = 1e-10) where {Nv,T}

    # Check if the number of objective functions is higher than 1.
    ( Nf <= 0 ) && @error("The number of objective functions must be higher than 0.")

    MGEO_Structure{Nv, Nf, T}(τ, ngen_max, run_max, design_vars, mgeo_eps)
end


text:
function (::Type{MT})(primitives::AbstractPrimitive{3, T}...;
                                         samples=(128,128,128),
                                         algorithm=MarchingCubes(),
                                         adf=false) where {T, MT <: AbstractMesh}

    f(x) = FRep(primitives[1], x)
    meshes = Vector{MT}(undef, length(primitives))
    for i = 1:length(primitives)
        b = HyperRectangle(primitives[i])
        meshes[i] = MT(x -> FRep(primitives[i], x), b, samples, algorithm)
    end
    return merge(meshes)
end

function piped_mesh(m::AbstractMesh,r)
    c = nothing
    for f in m.faces
        v1 = m.vertices[f[1]-O]
        v2 = m.vertices[f[2]-O]
        v3 = m.vertices[f[3]-O]
        c = CSGUnion(c, Piping(r, Point{3,Float64}[v1, v2, v3]))
    end
    c
end


text:
#=

  "PANDA: Phylogenetic ANalyses of DiversificAtion" module package

=#

module PANDA

__precompile__(true) #put back to true when checking done

#=
 Submodules
=#

# Utilities
include("Utils.jl")

# ESSE: Environmental and State Dependent Diversification
include("ESSE.jl")

# TRIBE: Trait and Range Interspecific Biogeographic Evolution
include("TRIBE.jl")

# ClaDS: Cladogenetic XX
include("ClaDS.jl")

#=
 Exported functions
=#

using .ESSE: esse, simulate_sse
export esse, simulate_sse

using .TRIBE: tribe, simulate_tribe
export tribe, simulate_tribe

using .ClaDS: Tree, tip_labels, n_tips, sample_tips,
    sim_ClaDS2_ntips, plot_ClaDS, infer_ClaDS,
    CladsOutput, print_CladsOutput, plot_CladsOutput, tip_rate,
    load_tree, save_ClaDS_in_R
export Tree, tip_labels, n_tips, sample_tips,
    sim_ClaDS2_ntips, plot_ClaDS, infer_ClaDS,
    CladsOutput, print_CladsOutput, plot_CladsOutput, tip_rate,
    load_tree, save_ClaDS_in_R
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
module Jeff

include("reflect.jl")
include("data.jl")
include("model.jl")
include("objective.jl")

end # module


text:
using Pkg
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")

const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

const is_TRAVIS = haskey(ENV,"TRAVIS")

is_CI = haskey(ENV,"CI")

@time begin
  if GROUP == "All" || GROUP == "Logging"
      @time @safetestset "AdaptiveLossLogNoImport" begin
        using Pkg
        neuralpde_dir = dirname(abspath(joinpath(@__DIR__, "..", "..", "..")))
        @info "loading neuralpde package at : $(neuralpde_dir)"
        neuralpde = Pkg.PackageSpec(path = neuralpde_dir)
        Pkg.develop(neuralpde)
        @info "making sure that there are no logs without having imported NeuralPDELogging"
        ENV["LOG_SETTING"] = "NoImport"
        include("adaptive_loss_log_tests.jl")
      end
      @time @safetestset "AdaptiveLossLogImportNoUse" begin
        using Pkg
        neuralpde_dir = dirname(abspath(joinpath(@__DIR__, "..", "..", "..")))
        @info "loading neuralpde package at : $(neuralpde_dir)"
        neuralpde = Pkg.PackageSpec(path = neuralpde_dir)
        Pkg.develop(neuralpde)
        @info "making sure that there are still no logs now that we have imported NeuralPDELogging"
        ENV["LOG_SETTING"] = "ImportNoUse"
        include("adaptive_loss_log_tests.jl")
      end
      @time @safetestset "AdaptiveLossLogImportUse" begin
        using Pkg
        neuralpde_dir = dirname(abspath(joinpath(@__DIR__, "..", "..", "..")))
        @info "loading neuralpde package at : $(neuralpde_dir)"
        neuralpde = Pkg.PackageSpec(path = neuralpde_dir)
        Pkg.develop(neuralpde)
        ENV["LOG_SETTING"] = "ImportUse"
        @info "making sure that logs are generated now if we use a logger"
        include("adaptive_loss_log_tests.jl")
      end
  end
end


text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setup.jl");

# periodic boundary conditions
pbc = [1; 1; 1]

# path to the database
datapath = "/Users/ngoccuongnguyen/Dropbox (MIT)/Research/Software/MDP/fitting/snap/Ta/JSON/";
dataformat = "json"
fileextension = "json"
atomspecies = ["Ta"];

# folders in the datapath
folders = ["Displaced_A15", "Displaced_BCC", "Displaced_FCC", "Elastic_BCC",
            "Elastic_FCC", "GSF_110", "GSF_112", "Liquid", "Surface",
            "Volume_A15", "Volume_BCC", "Volume_FCC"];

# weights for energies, forces, and stresses in the linear fit
weightinner=[100            1               1.00E-08
            100             1               1.00E-08
            100             1               1.00E-08
            1.00E-08        1.00E-08        0.0001
            1.00E-09        1.00E-09        1.00E-09
            100             1               1.00E-08
            100             1               1.00E-08
            4.67E+02        1               1.00E-08
            100             1               1.00E-08
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09]

# randomly selecting the configurations in the database
randomize = false;

# use all the data
percentage = 100.0;

# translate atom positions
translationvector = nothing

# rotate atom positions
rotationmatrix = nothing

# transpose lattice vectors because the lattice ordering is python-based
transposelattice = true

# training data
for i = 1:length(folders)
    traindata[i] = adddata(datapath * folders[i], dataformat, fileextension,
                percentage, randomize, atomspecies, weightinner[i,:], translationvector,
                rotationmatrix, transposelattice)
end

rcut = 5.0   # cut-off radiuss

# single descriptors to catpure energy of isolated pure elements
descriptors[1] = PotentialStruct("single", "nonbonded", [1], "single", [0.0], rcut);
include(descriptors[1].potentialfunction * ".jl");  # include the descriptors file

# Bessel descriptors
besselscale = 2.5;
eta = [rcut; besselscale];
kappa = [6];      # number of Bessel descriptors
descriptors[2] = Potential.PotentialStruct("pair", "nonbonded", [1, 1], "scaledbessel", [0.0], rcut);
include(descriptors[2].potentialfunction * ".jl");  # include the descriptors file

# SNAP parameters
twojmax = 6
ntypes = 1
rcutfac = 4.67637;
rfac0 = 0.99363;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)

# SNAP descriptors
descriptors[3] = initsna(snaparam, elemradius, elemweight)

# loss function style must be: "energy", "force", "energyforce", "energyforcestress"
lossfunc = "energyforcestress"
# use least-square method
method = "lsq"
# divide energies by the number of atoms
normalizeenergy = true
# divide virials by the box volume
normalizestress = true

# weights for energies, forces, and stresses in the outer nonlinear optimization
weightouter = [0.25, 0.75, 2e-7]

# bounds for cut-off radius to be optimized
rcutrange = [4.6, 5.4]

# bounds for scaling parameleter to be optimized
scalerange = [2.0, 4.0]

# define range for nonlinear parameters to be optimized
etarange =  vcat(reshape(rcutrange,(1,2)), reshape(scalerange,(1,2)))

# number of interpolation points for each nonlinear parameters
N = [8, 8]
etaspace = [etarange N]

# optimization parameters
optim = setoptim(lossfunc, method, normalizeenergy, normalizestress, eta, kappa, weightouter, etaspace)

# optimize the Bessel/SNAP potential
eta, coeff, fmin, iter, polycoeff = optimize2(traindata, traindata, descriptors, potentials, optim, pbc)

print("SNAP Coeffficients: "), show(stdout, "text/plain", coeff)

# compute unweighted MAE, RMSE, RSQ errors
energyerrors, forceerrors, stresserrors = validate2(traindata, descriptors, potentials, optim, coeff, pbc)

printerrors(folders, energyerrors, "Energy Errors")
printerrors(folders, forceerrors, "Force Errors")
printerrors(folders, stresserrors, "Stress Errors")

using DelimitedFiles
etaspace = [etarange N [32, 32]]
etaout, fout = mpolyeval(etaspace, polycoeff)
out = [etaout fout]
writedlm("optbesselsnaplossfunction.txt", out)

# ----------------------------------------------------------------------------------------
# Energy Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.04407201624074    |    0.13895234006755    |    0.99848863618353   |
# Displaced_A15 |   0.00064295683681    |    0.00084732580228    |    0.90326784300319   |
# Displaced_BCC |   0.00266815814537    |    0.00279961784024    |    0.99672203251388   |
# Displaced_FCC |   0.00159030140295    |    0.00172387396161    |    0.81493658950080   |
# Elastic_BCC   |   0.01091363883685    |    0.01100632694440    |    -588.43048075524   |
# Elastic_FCC   |   0.00302729308650    |    0.00306579881163    |    -1.3064502487690   |
# GSF_110       |   0.00366406390478    |    0.00400799472712    |    0.88510147584172   |
# GSF_112       |   0.00325432404242    |    0.00354119150458    |    0.96659981570708   |
# Liquid        |   0.00084617992731    |    0.00103532130901    |    0.99987144804575   |
# Surface       |   0.01079599316896    |    0.01500070756280    |    0.94230931245908   |
# Volume_A15    |   0.07848650524447    |    0.10962449687720    |    0.99885130263424   |
# Volume_BCC    |   0.24761485607900    |    0.34764026600377    |    0.99887102660169   |
# Volume_FCC    |   0.21855205786430    |    0.36344650749071    |    0.99444598091703   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Force Errors  |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.07363846374494    |    0.12575317788972    |    0.96777978926350   |
# Displaced_A15 |   0.09446578972838    |    0.11852562750111    |    0.96530855378974   |
# Displaced_BCC |   0.10660973760926    |    0.13590313156901    |    0.98968714163891   |
# Displaced_FCC |   0.09629424008693    |    0.12054105144872    |    0.94274401239914   |
# Elastic_BCC   |   0.05808574102001    |    0.07122021832773    |    -46968.652761763   |
# Elastic_FCC   |   0.03751218866409    |    0.04863178583668    |    -14439.203632284   |
# GSF_110       |   0.03187074058294    |    0.05363994853004    |    0.83204957019260   |
# GSF_112       |   0.04774908683784    |    0.07436887168531    |    0.84499362823215   |
# Liquid        |   0.26273384173361    |    0.33833743022054    |    0.95786271517765   |
# Surface       |   0.04000342054590    |    0.09161969982319    |    0.83889748233340   |
# Volume_A15    |   2.89993087276951    |    4.85293370653430    |    -Inf000000000000   |
# Volume_BCC    |   2.76832107868402    |    5.25683612049748    |    -Inf000000000000   |
# Volume_FCC    |   2.09738287182340    |    4.06096724650454    |    -Inf000000000000   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Stress Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   16575.0798602557    |    87940.0330737213    |    0.99740626554616   |
# Displaced_A15 |   15490.4558881970    |    21516.7061763517    |    -1.7547031528251   |
# Displaced_BCC |   11778.4200009532    |    16358.5489271636    |    0.63123812456294   |
# Displaced_FCC |   915.502519376297    |    1232.52113746799    |    0.99685768641770   |
# Elastic_BCC   |   680.783042315670    |    906.291087830642    |    0.98772277016923   |
# Elastic_FCC   |   8359.52276223574    |    11754.3177725659    |    0.44593828295390   |
# GSF_110       |   1828.42673919399    |    2901.92378103316    |    0.58993452570032   |
# GSF_112       |   1754.88638599766    |    2747.82369773695    |    0.73828261524334   |
# Liquid        |   43409.1051199719    |    60496.6513110316    |    -6.4577661295236   |
# Surface       |   2344.03298702141    |    4224.36580573058    |    0.80717237641878   |
# Volume_A15    |   25323.0104544580    |    56682.3988563785    |    0.99882817761392   |
# Volume_BCC    |   93126.9641553286    |    232208.518467598    |    0.99841748447201   |
# Volume_FCC    |   61878.6776523236    |    223359.641907389    |    0.99194235547032   |
# ----------------------------------------------------------------------------------------

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Printf

function Base.show(io::IO, organism::TierrianOrganism, model::TierraModel)
    println(io, "-"^10)

    println(io, organism.hash)
    println(io, organism.parent_hash)

    println(io, "*Last Instruction*")

    ip = _address_within(organism, organism.ip - 1)
    instruction_int = read_memory(model, ip)
    instruction = convert_instruction(instruction_int)
    @printf(io, "%16i \t %5i %s \n", ip, instruction_int, instruction)

    println(io, "*GENOME*")

    for line::UInt16 in organism.start_address:(organism.start_address - 1 + organism.length)
        instruction_int = read_memory(model, line)
        instruction = convert_instruction(instruction_int)
        if line == ip
            @printf(io, "=> %16i \t %5i %s \n", line, instruction_int, instruction)
        else
            @printf(io, "   %16i \t %5i %s \n", line, instruction_int, instruction)
        end
    end

    println(io, "*STACK*")

    for element in organism.stack
        @printf(io, "%16i", element)
    end

    println(io, "*Registers*")

    @printf(io, "a: %16i \t", organism.a)
    @printf(io, "b: %16i \t", organism.b)
    @printf(io, "c: %16i \t", organism.c)
    @printf(io, "d: %16i \n", organism.d)
end

function Base.show(io::IO, model::TierraModel)
    println(io, "-"^50)

    hash_count = Dict()

    for key in model.organism_keys
        hash = model.organisms[key].hash

        if haskey(hash_count, hash)
            hash_count[hash] += 1
        else
            hash_count[hash] = 1
        end
    end

    for (hash, count) in hash_count
        println(io, hash, ":\t ", count)
    end

    return

    for key in model.organism_keys[1:min(length(model.organism_keys), 5)]
        print(key)
        organism = model.organisms[key]
        show(io, organism, model)
    end
end

function Base.show(io::IO, free_block::FreeMemoryBlock)
    @printf(io, "%16i \t %16i \n", free_block.start_address, free_block.length)
end

function print_program(program::Vector{UInt8})
    for instruction_int in program
        instruction = convert_instruction(instruction_int)
        println(instruction)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function get_space_variables( fname::String )

    #Get Variables
    fid         = h5open( fname , "r")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Discretize space
    dx  = (xmax - xmin)/(Nx-1)
    x   = [xmin + (i-1)*dx for i in 1:Nx]

    dy  = (ymax - ymin)/(Ny-1)
    y   = [ymin + (i-1)*dy for i in 1:Ny]

    close(fid)

    return x , y

end


text:
#!jl # # EFTfitter.jl - Advanced Tutorial
#jl # EFTfitter.jl - Advanced Tutorial

#jl #~This tutorial introduces some of the advanced functionalities of EFTfitter.jl
#jl #~using a generic example. Please see the [tutorial] for basic usage of EFTfitter.
#md # Table of contents:
#md # ```@contents
#md # Pages = ["advanced_tutorial.md"]
#md # Depth = 3
#md # ```
                                                                                #src
#nb using EFTfitter
#nb using BAT           # for sampling
#nb using IntervalSets  # for specifying the prior
#nb using Distributions # for specifying the prior
#nb using Plots         # for plotting
                                                                                #src
#!jl include("lit_advanced_tutorial_inputs.jl")
#nb # ## File "runAdvancedTutorial.jl"
                                                                                #src
#jl using EFTfitter
#jl using BAT            # for sampling
#jl using IntervalSets   # for specifying the prior
#jl using Distributions  # for specifying the prior
#jl using Plots # for plotting
#jl
#jl include("advanced_inputs.jl")

#~We need to modify the definition of the `EFTfitterModel` by also passing the `nuisance_correlations`:
model = EFTfitterModel(parameters, measurements, correlations, nuisance_correlations)

#!md posterior = PosteriorDensity(model)
#!md
#!md algorithm = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)
#!md samples = bat_sample(posterior, algorithm).result;

#!md p = plot(samples)
#!nb savefig(p, "plot.pdf")
#md # ![plot with nuisances](plots/plot_nuisance.png)
#md  # !!! warning
#md  #     When using nuisance correlations, it may happen that the fit does not converge anymore.
#md  #     This happens mainly when the correlation values are close to -1 or +1.
#md  #     In such a case, you can try to reduce the allowed ranges in the priors for the Nuisance
#md  #     correlations to investigate at which values this happens.

#!md #~Warning:
#!md #~When using nuisance correlations, it may happen that the fit does not converge anymore.
#!md #~This happens mainly when the correlation values are close to -1 or +1.
#!md #~In such a case, you can try to reduce the allowed ranges in the priors for the Nuisance
#!md #~correlations to investigate at which values this happens.

#nb # ### Ranking of measurements and uncertainties
#md # ## Ranking of measurements and uncertainties
#jl #~============= Ranking of measurements and uncertainties ========================#
#~With the `rank_measurements` and `rank_uncertainties` functions, the influence of
#~the individual measurements or uncertainty types on the result of a fit can be estimated.
#~For the ranking, each active measurement (respectively uncertainty type) is deactivated
#~at a time and the fit is repeated. The results of all fits with one deactivated
#~measurement (or uncertainty type) are then compared to the fit result with all measurements (uncertainty types) activated.
#~A ranking is calculated based on a ranking criterion calculated from the posterior distributions of these fits.

#~The default ranking criterion is the relative increase of the total width of the
#~smallest interval containing 90% of the posterior probability when deactivating a measurement.
#~For models with more than one parameter, the sum of relative increases of all
#~one-dimensional smallest intervals is used. The default ranking criterion is therefore
#`SumOfSmallestIntervals(p=0.9, bins=200)`.
measurement_ranking = EFTfitter.rank_measurements(model)

#~For ranking the uncertainty types, the relative decrease is used.
uncertainty_ranking = EFTfitter.rank_uncertainties(model, criterion = SumOfSmallestIntervals(p=0.9, bins=200))

#~Please see the [ranking documentation] for further ranking criteria and keyword arguments.


text:
#!/usr/bin/env julia

const data = Dict{Tuple{String,String},Float64}()

for arg in ARGS, line in eachline(arg)
    lang, bench, time_str = split(line, ',')
    old_time = get(data, (lang, bench), Inf)
    new_time = parse(Float64, time_str)
    0 < new_time < old_time || continue
    data[lang, bench] = new_time
end

for ((lang, bench), min_time) in sort!(collect(data))
    println("$lang,$bench,$min_time")
end

max_stars_repo_name: MilesCranmer/Microbenchmarks

max_issues_repo_name: MilesCranmer/Microbenchmarks

max_forks_repo_name: MilesCranmer/Microbenchmarks


text:
# include(pwd()*"/src/SFFM.jl")
using Plots, SFFM

## define a model
T = [0.0]
C = [1]

rfun(x) = x.*0
Rfun(x) = r(x)

r = (
    r = function (x)
        rfun(x)
    end,
    R = function (x)
        Rfun(x)
    end
)

bounds = [0 12; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = bounds)

orders = [1;3;5;7;11;13;15;21;25;27;33;;37;41]
errors_1 = []

Δtemp = 1/2 # the grid size; must have kΔ = 1 for some k due to discontinuity in r at 1
nodes = collect(0:Δtemp:bounds[1,2])

for order in orders
    # order = 3
    println("order = "*string(order))
    dgmesh = SFFM.DGMesh(
        model,
        nodes,
        order,
        Basis = "lagrange",
    )
    frapmesh = SFFM.FRAPMesh(
        model,
        nodes,
        order,
    )
    fvmesh = SFFM.FVMesh(
        model,
        order,
        collect(0:Δtemp/order:bounds[1,2]),
    )
    simmesh = SFFM.FVMesh(
        model,
        1,
        nodes,
    )

    trueprobs = zeros(Float64,1,SFFM.NIntervals(frapmesh),SFFM.NPhases(model))
    truepos = convert(Int,ceil((1.2+eps())/Δtemp))
    trueprobs[truepos] = 1
    groundtruth = SFFM.SFFMProbability(
        [0],
        trueprobs,
        SFFM.CellNodes(frapmesh),
    )

    # DG
    B_DG = SFFM.MakeFullGenerator(model, dgmesh)
    #ME
    # me = SFFM.MakeME(SFFM.CMEParams[order], mean = Δtemp)
    B_ME = SFFM.MakeFullGenerator(model, frapmesh)
    # Erlang (this is the erlang which is equivalent to DG)
    erlang = SFFM.MakeErlang(order, mean = 1)
    B_Erlang = SFFM.MakeFullGenerator(model, frapmesh, erlang)
    # meph (this is the erlang treated as an ME)
    meph = SFFM.ME(erlang.a, erlang.S, erlang.s; D = SFFM.erlangDParams[string(order)])
    B_MEPH = SFFM.MakeFullGenerator(model, frapmesh, meph)
    # FVM
    B_FV = SFFM.MakeFullGenerator(model, fvmesh)

    point = 0+eps()
    pointIdx = convert(Int,ceil(point/Δtemp))
    begin
        V = SFFM.vandermonde(order)
        theNodes = SFFM.CellNodes(dgmesh)[:,pointIdx]
        basisValues = zeros(length(theNodes))
        for n in 1:length(theNodes)
            basisValues[n] = prod(point.-theNodes[[1:n-1;n+1:end]])./prod(theNodes[n].-theNodes[[1:n-1;n+1:end]])
        end
        initpm = [
            zeros(sum(model.C.<=0)) # LHS point mass
            zeros(sum(model.C.>=0)) # RHS point mass
        ]
        initprobs = zeros(Float64,SFFM.NBases(dgmesh),SFFM.NIntervals(dgmesh),SFFM.NPhases(model))
        initprobs[:,pointIdx,1] = basisValues'*V.V*V.V'.*2/Δtemp
        initdist = SFFM.SFFMDensity(
            initpm,
            initprobs,
            SFFM.CellNodes(dgmesh),
        ) # convert to a distribution object so we can apply Dist2Coeffs
        # convert to Coeffs α in the DG context
        x0_DG = SFFM.Dist2Coeffs(model, dgmesh, initdist)
    end
    begin
        initpm = [
            zeros(sum(model.C.<=0)) # LHS point mass
            zeros(sum(model.C.>=0)) # RHS point mass
        ]
        initprobs = zeros(Float64,SFFM.NBases(frapmesh),SFFM.NIntervals(frapmesh),SFFM.NPhases(model))
        initprobs[:,pointIdx,1] = me.a
        initdist = SFFM.SFFMDensity(
            initpm,
            initprobs,
            SFFM.CellNodes(frapmesh),
        ) # convert to a distribution object so we can apply Dist2Coeffs
        # convert to Coeffs α in the DG context
        x0_ME = SFFM.Dist2Coeffs( model, frapmesh, initdist)
    end
    begin
        initpm = [
            zeros(sum(model.C.<=0)) # LHS point mass
            zeros(sum(model.C.>=0)) # RHS point mass
        ]
        initprobs = zeros(Float64,SFFM.NBases(frapmesh),SFFM.NIntervals(frapmesh),SFFM.NPhases(model))
        initprobs[1,pointIdx,1] = 1
        initdist = SFFM.SFFMDensity(
            initpm,
            initprobs,
            SFFM.CellNodes(frapmesh),
        ) # convert to a distribution object so we can apply Dist2Coeffs
        # convert to Coeffs α in the DG context
        x0_Erlang = SFFM.Dist2Coeffs( model, frapmesh, initdist)
    end
    begin
        initpm = [
            zeros(sum(model.C.<=0)) # LHS point mass
            zeros(sum(model.C.>=0)) # RHS point mass
        ]
        initprobs = zeros(Float64,SFFM.NBases(fvmesh),SFFM.NIntervals(fvmesh),SFFM.NPhases(model))
        initprobs[1,pointIdx,1] = 1
        initdist = SFFM.SFFMDensity(
            initpm,
            initprobs,
            SFFM.CellNodes(fvmesh),
        ) # convert to a distribution object so we can apply Dist2Coeffs
        # convert to Coeffs α in the DG context
        x0_FV = SFFM.Dist2Coeffs( model, fvmesh, initdist)
    end

    euler(B,x0) = SFFM.EulerDG( B, 1.2, x0, h = 0.0001)
    x1_DG = euler(B_DG.B, x0_DG)
    x1_ME = euler(B_ME.B, x0_ME)
    x1_Erlang = euler(B_Erlang.B, x0_Erlang)
    x1_MEPH = euler(B_MEPH.B, x0_Erlang)
    x1_FV = euler(B_FV.B, x0_FV)

    x1_DG = SFFM.Coeffs2Dist(
        model,
        dgmesh,
        x1_DG,
        SFFM.SFFMProbability,
    )
    x1_ME = SFFM.Coeffs2Dist(
        model,
        frapmesh,
        x1_ME,
        SFFM.SFFMProbability,
    )
    x1_Erlang = SFFM.Coeffs2Dist(
        model,
        frapmesh,
        x1_Erlang,
        SFFM.SFFMProbability,
    )
    x1_MEPH = SFFM.Coeffs2Dist(
        model,
        frapmesh,
        x1_MEPH,
        SFFM.SFFMProbability,
    )

    x1_FV = SFFM.Coeffs2Dist(
        model,
        frapmesh,
        x1_FV,
        SFFM.SFFMProbability,
    )

    errVec_1 = (
        SFFM.starSeminorm(x1_DG, groundtruth),
        SFFM.starSeminorm(x1_ME, groundtruth),
        SFFM.starSeminorm(x1_Erlang, groundtruth),
        SFFM.starSeminorm(x1_MEPH, groundtruth),
        SFFM.starSeminorm(x1_FV, groundtruth),
    )
    p = SFFM.plot(model, frapmesh, x1_ME,
        color = 1, label = "DG")
    display(p)

    p = SFFM.plot(model, dgmesh, x1_DG,
        color = 1, label = "DG", alpha = 0.25)
    p = SFFM.plot!(p, model, frapmesh, x1_ME,
        color = 2, label = "ME", alpha = 0.25)
    # SFFM.plot!(p, model, frapmesh, x1_Erlang,
    #     color = 3, label = "Erlang", alpha = 0.25)
    # SFFM.plot!(p, model, frapmesh, x1_MEPH,
    #     color = 5, label = "ME-PH", alpha = 0.25)
    SFFM.plot!(p, model, frapmesh, x1_FV,
        color = 7, label = "FV", alpha = 0.25)
    p = plot!(title = "approx dist at t=1.2; order = "*string(order), subplot = 1)
    display(p)

    push!(errors_1, errVec_1)
end

q = plot(
    xlabel = "order",
    ylabel = "log10 error",
    title = "error for t=1.2",
    legend = :bottomleft
)
methodNames = ["DG", "ME", "Erlang", "MEPH", "FV"]
shapes = [:none, :none, :circle, :none, :none]

for whichOrder in 2:length(orders)
# order = 3
    for whichMethod in 1:length(methodNames)
        plot!(
            q,
            [orders[whichOrder-1];orders[whichOrder]],
            [
                log10(errors_1[whichOrder-1][whichMethod]);
                log10(errors_1[whichOrder][whichMethod])
            ],
            label = whichOrder == 2 && methodNames[whichMethod],
            colour = whichMethod,
            markershape = shapes[whichMethod],
            seriestype = :line,
        )
    end
end
display(q)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using ProblemParser
using Utils

using Chain
using StaticArrays

const file = find_input(@__FILE__)
const input = @chain file slurp parse(Lines(), _)

const opening = "([{<"
const closing = ")]}>"
const scores = SA[3,57,1197,25137]

function corrupted(line::AbstractString)
    stack = Char[]
    for c in line
        pos = findfirst(c, opening)
        if pos === nothing
            pos = findfirst(c, closing)
            (isempty(stack) || pop!(stack) != c) && return scores[pos]
        else
            push!(stack, closing[pos])
        end
    end
    0
end

function part1(input)
    @chain input begin
        map(corrupted, _)
        sum
    end
end

assertequal(part1(input), 358737)

function complete(line::AbstractString)
    stack = Char[]
    for c in line
        pos = findfirst(c, opening)
        if pos === nothing
            pos = findfirst(c, closing)
            @assert pop!(stack) == c
        else
            push!(stack, closing[pos])
        end
    end
    score = 0
    while !isempty(stack)
        score = score * 5 + findfirst(==(pop!(stack)), closing)
    end
    score
end
@assert complete("<{([{{}}[<[[[<>{}]]]>[]]") == 294

function part2(input)
    completed = @chain input begin
        filter(==(0) ∘ corrupted, _)
        map(complete, _)
    end
    middle = (length(completed) + 1) ÷ 2
    partialsort!(completed, middle)
end

assertequal(part2(input), 4329504793)

max_stars_repo_name: CmdQ/AoC


text:
using DiffResults

f = Dynamics.dubins_dynamics!
n,m = 3,2
f_aug! = f_augmented!(f,n,m)

fd_aug! = rk3(f_aug)
nm1 = n + m + 1

In = 1.0*Matrix(I,n,n)

# Initialize discrete and continuous dynamics Jacobians

# Get continuous dynamics jacobian
function gen_discrete_jacobians(fd_aug!,n,m,dt)
    nm1 = n+m+1
    Jd = zeros(nm1, nm1)
    s = zeros(nm1)
    sdot = zero(s)
    Fd!(Jd,Sdotd,Sd) = ForwardDiff.jacobian!(Jd,fd_aug!,Sdotd,Sd)
    s[end] = sqrt(dt)

    xi = 1:n
    ui = n .+ (1:m)

    Jdi = CartesianIndices(Jd)
    xx = Jdi[xi,xi]
    xu = Jdi[xi,ui]
    xu2 = Jdi[xi,1:m]

    function fd_jacobians!(fdx,fdu,x,u)
        # Assign state, control (and dt) to augmented vector
        s[1:n] = x
        s[n+1:n+m] = u[1:m]

        # Calculate Jacobian
        Fd!(Jd,sdot,s)

        fdx[1:n,1:n] = Jd[1:n,1:n]
        fdu[1:n,1:m] = Jd[1:n,n.+(1:m)]

    end

    function fd_jacobians2!(fdx,fdu,x,u)
        # Assign state, control (and dt) to augmented vector
        s[xi] = x
        s[ui] = u

        # Calculate Jacobian
        Fd!(Jd,sdot,s)

        fdx[xx] = Jd[xx]
        fdu[xu2] = Jd[xu]

    end
    return fd_jacobians!, fd_jacobians2!
end

function gen_continuous_jacobians(f_aug!,n,m,dt)
    xi = 1:n
    ui = n .+ (1:m)

    Jc = zeros(n+m,n+m)
    z = zeros(n+m)
    zdot = zeros(n+m)
    A = view(Jc,xi,xi)
    B = view(Jc,xi,ui)
    function fc_jacobians2!(x,u)
        z[xi] = x
        z[ui] = u
        ForwardDiff.jacobian!(Jc,f_aug!,zdot,z)
        return A,B
    end

    function fc_jacobians!(x,u)
        # infeasible = size(u,1) != m̄
        z[1:n] = x
        z[n+1:n+m] = u[1:m]
        Fc!(Jc,zdot,z)
        return Jc[1:n,1:n], Jc[1:n,n+1:n+m] # fx, fu
    end

    function rk3_jacobian!(x,u)
        z[xi] = x
        z[ui] = u
        A1,B1 = fc_jacobians2!(x,u)
        k1 = zdot*dt
        A2,B2 = fc_jacobians2!(x+k1[xi]/2,u)
        k2 = zdot*dt
        A3,B3 = fc_jacobians2!(x-k1[xi]+2*k2[xi],u)
        k3 = zdot*dt

        dk1_x = A1*dt
        dk2_x = (I + dk1_x/2)'A2*dt
        dk3_x = (I - dk1_x + 2*dk2_x)'A3*dt
        Ad = I + (dk1_x + 4*dk2_x + dk3_x)/6

        dk1_u = B1*dt
        dk2_u = A2'*(dk1_u/2)*dt + B2*dt
        dk3_u = A3'*(-dk1_u + 2*dk2_u)*dt + B3*dt
        Bd = (dk1_u + 4*dk2_u + dk3_u)/6
        return Ad,Bd
    end

    return fc_jacobians!, fc_jacobians2!, rk3_jacobian!
end

f_aug! = f_augmented!(f,n,m)
dt = 0.1
fc_jacob, fc_jacob2, fd_jacob = gen_continuous_jacobians(f_aug!,n,m,dt)

x,u = rand(n), rand(m)
xdot = zero(x)
z = [x;u]
zdot = zero(z)
s = [z;u;sqrt(dt)]
sdot = zero(s)
fc_jacob(x,u) == fc_jacob2(x,u)
fd_jacob(x,u)

# Discrete jacobians
fd! = rk3(f,dt)
fd_aug! = rk3(f_aug)
fd!(xdot,x,u)
fd_aug!(sdot,s)
xdot == sdot[1:n]

fd_jacob!, fd_jacob2! = gen_discrete_jacobians(fd_aug!,n,m,dt)
fdx = zeros(n,n)
fdu = zeros(n,m)
fd_jacob!(fdx,fdu,x,u) == fd_jacob2!(fdx,fdu,x,u)

@btime fc_jacob($x,$u)
@btime fc_jacob2($x,$u)
@btime fc_jacob2($x,$u)

@btime fd_jacob(x,u)
@btime fd_jacob!(fdx,fdu,x,u)

nm1 = n+m+1
Jd = zeros(nm1, nm1)

xi = 1:n
ui = n .+ (1:m)

Jdi = CartesianIndices(Jd)
xx = Jdi[xi,xi]
xu = Jdi[xi,ui]


text:
module TensorValues

using StaticArrays
using Base: @propagate_inbounds, @pure

export MultiValue
export TensorValue
export VectorValue

export inner, outer, meas
export det, inv, tr, dot, norm
export mutable
export trace
export symmetic_part

import Base: show
import Base: zero, one
import Base: +, -, *, /, \, ==, ≈
import Base: conj
import Base: sum, maximum, minimum
import Base: getindex, iterate, eachindex
import Base: size, length, eltype
import Base: reinterpret
import Base: convert
import Base: CartesianIndices
import Base: LinearIndices
import Base: adjoint

import LinearAlgebra: det, inv, tr, dot, norm

include("Types.jl")

include("Indexing.jl")

include("Operations.jl")

include("Reinterpret.jl")

end # module


text:
# Autogenerated wrapper script for BLST_jll for x86_64-apple-darwin
export libblst

JLLWrappers.@generate_wrapper_header("BLST")
JLLWrappers.@declare_library_product(libblst, "@rpath/libblst.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libblst,
        "lib/libblst.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
ENV["MPLBACKEND"]="agg" # no GUI

using PyPlot, PyCall, Test

@info("PyPlot is using Matplotlib $(PyPlot.version) with Python $(PyCall.pyversion)")

plot(1:5, 2:6, "ro-")

line = gca().lines[1]
@test line.get_xdata() == [1:5;]
@test line.get_ydata() == [2:6;]

fig = gcf()
@test isa(fig, PyPlot.Figure)
if PyPlot.version >= v"2"
    @test fig.get_size_inches() ≈ [6.4, 4.8]
else # matplotlib 1.3
    @test fig.get_size_inches() ≈ [8, 6]
end

# with Matplotlib 1.3, I get "UserWarning: bbox_inches option for ps backend is not implemented yet"
if PyPlot.version >= v"2"
    s = sprint(show, "application/postscript", fig);
    m = match(r"%%BoundingBox: *([0-9]+) +([0-9]+) +([0-9]+) +([0-9]+)", s)
    @test m !== nothing
    boundingbox = map(s -> parse(Int, s), m.captures)
    @info("got plot bounding box ", boundingbox)
    @test all([300, 200] .< boundingbox[3:4] - boundingbox[1:2] .< [450,350])
end

c = get_cmap("RdBu")
a = 0.0:0.25:1.0
rgba = pycall(pycall(PyPlot.ScalarMappable, PyObject, cmap=c,
                     norm=PyPlot.Normalize01)."to_rgba", PyArray, a)
@test rgba ≈ [  0.403921568627451   0.0                  0.12156862745098039  1.0
                0.8991926182237601  0.5144175317185697   0.4079200307574009   1.0
                0.9657054978854287  0.9672433679354094   0.9680891964628989   1.0
                0.4085351787773935  0.6687427912341408   0.8145328719723184   1.0
                0.0196078431372549  0.18823529411764706  0.3803921568627451   1.0 ]


text:
# Interface methods

function rastername(::Type{EarthEnv{LandCover}}, layer::Integer, discover::Bool)
    filetype = discover ? "complete" : "partial"
    "landcover_$(filetype)_$(layer).tif"
end

rasterpath(T::Type{EarthEnv{LandCover}}, layer::Integer, discover::Bool) =
    joinpath(rasterpath(T), rastername(T, layer, discover))

function rasterurl(T::Type{EarthEnv{LandCover}}, layer::Integer, discover::Bool)
    stem = discover ? "with_DISCover/consensus_full_class_$(layer).tif" :
                      "without_DISCover/Consensus_reduced_class_$(layer).tif"
    joinpath(rasterurl(T), stem)
end

function download_raster(T::Type{EarthEnv{LandCover}}; layer::Integer=1, discover::Bool=false)
    1 ≤ layer ≤ 12 || throw(ArgumentError("The layer must be between 1 and 12"))
    url = rasterurl(T, layer, discover)
    path = rasterpath(T, layer, discover)
    return _maybe_download(url, path)
end

_pathsegment(::Type{LandCover}) = "consensus_landcover"


text:
function n_extant_tips(tree::Tree)
    function aux(subtree)
        if subtree.n_nodes < 2
            if subtree.extant
                return 1
            else
                return  0
            end
        else
            return aux(subtree.offsprings[1]) + aux(subtree.offsprings[2])
        end
    end

    aux(tree)
end

function n_extant_tips(tree::Tree ,edge_trees::Array{EdgeTreeRates2,1})

    n::Int64 = 0

    for et in edge_trees
        n += et.tip_number
    end

    return n
end

function n_left(tree::Tree)
    function aux(subtree, x)
        if length(subtree.offsprings) == 0
            pushfirst!(x,0)
        else
            aux(subtree.offsprings[2], x)
            aux(subtree.offsprings[1], x)
            pushfirst!(x,subtree.offsprings[1].n_nodes)
        end
    end

    x = Array{Int64,1}(undef,0)
    aux(tree, x)
    return x
end


text:
module Interface

include("SpectralSolver.jl")
include("NewtonSolver.jl")
include("SchwarszchildModes.jl")

using .NewtonSolver
using .Schwarzschild
using .SpectralSolver

#P = ModeParameters((l=2,s=2,m=0,n=0,a=0.01,ω = 0.373672 - 0.0889623*im,Alm = 0, Nmax = 300, lmax = 10))

function GetModes(l,m,n,s; amax = 0.99, ϵ = 0.01, Nmax = 300, lmax = 10)
    ##Initialize
    as = [a for a in 0.0:ϵ:amax]
    ωs = zero(im*(as))
    Almss = zero(im*(as))

    #Get Schwarzschild Modes
    ωsch, Almsch = SchwarzschildModes(l,n,s)

    ϵs = 0.0000001 #Finite difference step size

    P = ModeParameters((l=l,s=s,m=m,n=n,a=0.0,ω = ωsch,Alm = Almsch, Nmax = Nmax, lmax = lmax))

    NmaxTrials = [100,300,500,1000,5_000,10_000,50_000]
    lmaxTrials = [10,12,15,15,20,25]

    for i =1:length(as)
        #print(i," ")
        P.a = as[i]
        ω₀ = P.ω; Alm₀ = P.Alm
        ωnew = P.ω + 1.0; Almnew = Alm₀ + 1.0;
        iterations = 1;
        while ((abs(ωnew-ω₀) > 1e-10) | (abs(Almnew-Alm₀) > 1e-10)) & (iterations < 6)
            P.Nmax = NmaxTrials[iterations]
            P.lmax = lmaxTrials[iterations]
            ω₀ = P.ω; Alm₀ = P.Alm;
            Newton!(P,ϵs)
            ωnew = P.ω; Almnew = P.Alm
            iterations += 1
            #println("ω diff = ",abs(ωnew-ω₀)," Alm diff =  ",abs(Almnew-Alm₀))
        end
        #println("Number of Iterations = ", iterations)
        ωs[i] = P.ω
        Almss[i] = P.Alm
    end

    final = [i for i in zip(as,ωs,Almss)]
end

Ms = GetModes(2,2,2,2)

#for i in Ms
#   println(i)
#end

export GetModes
end

max_stars_repo_name: Potatoasad/AttemptAtQNM

max_issues_repo_name: Potatoasad/AttemptAtQNM

max_forks_repo_name: Potatoasad/AttemptAtQNM


text:
export Product, LibraryProduct, FileProduct, ExecutableProduct, satisfied, locate, @write_deps_file

"""
A `Product` is an expected result after building or installation of a package.
"""
abstract type Product end

"""
A `LibraryProduct` is a special kind of `Product` that not only needs to exist,
but needs to be `dlopen()`'able.  You must know which directory the library
will be installed to, and its name, e.g. to build a `LibraryProduct` that
refers to `"/lib/libnettle.so"`, the "directory" would be "/lib", and the
"libname" would be "libnettle".
"""
struct LibraryProduct <: Product
    dir_path::String
    libname::String

    """
    `LibraryProduct(prefix::Prefix, libname::AbstractString)`

    Declares a `LibraryProduct` that points to a library located within the
    `libdir` of the given `Prefix`, with a name containing `libname`.  As an
    example, given that `libdir(prefix)` is equal to `usr/lib`, and `libname`
    is equal to `libnettle`, this would be satisfied by the following paths:

        usr/lib/libnettle.so
        usr/lib/libnettle.so.6
        usr/lib/libnettle.6.dylib
        usr/lib/libnettle-6.dll

    Libraries matching the search pattern are rejected if they are not
    `dlopen()`'able.
    """
    function LibraryProduct(prefix::Prefix, libname::AbstractString)
        return LibraryProduct(libdir(prefix), libname)
    end

    """
    `LibraryProduct(dir_path::AbstractString, libname::AbstractString)`

    For finer-grained control over `LibraryProduct` locations, you may directly
    pass in the `dir_path` instead of auto-inferring it from `libdir(prefix)`.
    """
    function LibraryProduct(dir_path::AbstractString, libname::AbstractString)
        return new(dir_path, libname)
    end
end

"""
locate(lp::LibraryProduct; verbose::Bool = false,
        platform::Platform = platform_key())

If the given library exists (under any reasonable name) and is `dlopen()`able,
(assuming it was built for the current platform) return its location.  Note
that the `dlopen()` test is only run if the current platform matches the given
`platform` keyword argument, as cross-compiled libraries cannot be `dlopen()`ed
on foreign platforms.
"""
function locate(lp::LibraryProduct; verbose::Bool = false,
                platform::Platform = platform_key())
    if !isdir(lp.dir_path)
        if verbose
            info("Directory $(lp.dir_path) does not exist!")
        end
        return nothing
    end
    for f in readdir(lp.dir_path)
        # Skip any names that aren't a valid dynamic library for the given
        # platform (note this will cause problems if something compiles a `.so`
        # on OSX, for instance)
        if !valid_dl_path(f, platform)
            continue
        end

        if verbose
            info("Found a valid dl path $(f) while looking for $(lp.libname)")
        end

        # If we found something that is a dynamic library, let's check to see
        # if it matches our libname:
        if startswith(basename(f), lp.libname)
            dl_path = abspath(joinpath(lp.dir_path), f)
            if verbose
                info("$(dl_path) matches our search criteria of $(lp.libname)")
            end

            # If it does, try to `dlopen()` it if the current platform is good
            if platform == platform_key()
                hdl = Libdl.dlopen_e(dl_path)
                if hdl == C_NULL
                    if verbose
                        info("$(dl_path) cannot be dlopen'ed")
                    end
                else
                    # Hey!  It worked!  Yay!
                    Libdl.dlclose(hdl)
                    return dl_path
                end
            else
                # If the current platform doesn't match, then just trust in our
                # cross-compilers and go with the flow
                return dl_path
            end
        end
    end

    if verbose
        info("Could not locate $(lp.libname) inside $(lp.dir_path)")
    end
    return nothing
end

"""
An `ExecutableProduct` is a `Product` that represents an executable file.

On all platforms, an ExecutableProduct checks for existence of the file.  On
non-Windows platforms, it will check for the executable bit being set.  On
Windows platforms, it will check that the file ends with ".exe", (adding it on
automatically, if it is not already present).
"""
struct ExecutableProduct <: Product
    path::AbstractString

    """
    `ExecutableProduct(prefix::Prefix, binname::AbstractString)`

    Declares an `ExecutableProduct` that points to an executable located within
    the `bindir` of the given `Prefix`, named `binname`.
    """
    function ExecutableProduct(prefix::Prefix, binname::AbstractString)
        return ExecutableProduct(joinpath(bindir(prefix), binname))
    end

    """
    `ExecutableProduct(binpath::AbstractString)`

    For finer-grained control over `ExecutableProduct` locations, you may directly
    pass in the full `binpath` instead of auto-inferring it from `bindir(prefix)`.
    """
    function ExecutableProduct(binpath::AbstractString)
        return new(binpath)
    end
end

"""
`locate(fp::FileProduct; platform::Platform = platform_key(),
                         verbose::Bool = false)`

If the given executable file exists and is executable, return its path.

On all platforms, an ExecutableProduct checks for existence of the file.  On
non-Windows platforms, it will check for the executable bit being set.  On
Windows platforms, it will check that the file ends with ".exe", (adding it on
automatically, if it is not already present).
"""
function locate(ep::ExecutableProduct; platform::Platform = platform_key(),
                verbose::Bool = false)
    # On windows, we always slap an .exe onto the end if it doesn't already
    # exist, as Windows won't execute files that don't have a .exe at the end.
    path = if platform isa Windows && !endswith(ep.path, ".exe")
        "$(ep.path).exe"
    else
        ep.path
    end

    if !isfile(path)
        if verbose
            info("$(ep.path) does not exist, reporting unsatisfied")
        end
        return nothing
    end

    # If the file is not executable, fail out (unless we're on windows since
    # windows doesn't honor these permissions on its filesystems)
    @static if !iswindows()
        if uperm(path) & 0x1 == 0
            if verbose
                info("$(path) is not executable, reporting unsatisfied")
            end
            return nothing
        end
    end

    return path
end

"""
A `FileProduct` represents a file that simply must exist to be satisfied.
"""
struct FileProduct <: Product
    path::AbstractString
end

"""
locate(fp::FileProduct; platform::Platform = platform_key(),
                        verbose::Bool = false)

If the given file exists, return its path.  The platform argument is ignored
here, but included for uniformity.
"""
function locate(fp::FileProduct; platform::Platform = platform_key(),
                                 verbose::Bool = false)
    if isfile(fp.path)
        if verbose
            info("FileProduct $(fp.path) does not exist")
        end
        return fp.path
    end
    return nothing
end

"""
`satisfied(p::Product; platform::Platform = platform_key(),
                       verbose::Bool = false)`

Given a `Product`, return `true` if that `Product` is satisfied, e.g. whether
a file exists that matches all criteria setup for that `Product`.
"""
function satisfied(p::Product; platform::Platform = platform_key(),
                               verbose::Bool = false)
    return locate(p; platform=platform, verbose=verbose) != nothing
end

"""
`@write_deps_file(products...)`

Helper macro to generate a `deps.jl` file out of a mapping of variable name
to  `Product` objects. Call using something like:

    fooifier = ExecutableProduct(...)
    libbar = LibraryProduct(...)
    @write_deps_file fooifier libbar

If any `Product` object cannot be satisfied (e.g. `LibraryProduct` objects must
be `dlopen()`-able, `FileProduct` objects must exist on the filesystem, etc...)
this macro will error out.  Ensure that you have used `install()` to install
the binaries you wish to write a `deps.jl` file for, and, optionally that you
have used `activate()` on the `Prefix` in which the binaries were installed so
as to make sure that the binaries are locatable.

The result of this macro call is a `deps.jl` file containing variables named
the same as the keys of the passed-in dictionary, holding the full path to the
installed binaries.  Given the example above, it would contain code similar to:

    global const fooifier = "<pkg path>/deps/usr/bin/fooifier"
    global const libbar = "<pkg path>/deps/usr/lib/libbar.so"

This file is intended to be `include()`'ed from within the `__init__()` method
of your package.  Note that all files are checked for consistency on package
load time, and if an error is discovered, package loading will fail, asking
the user to re-run `Pkg.build("package_name")`.
"""
macro write_deps_file(capture...)
    # props to @tshort for his macro wizardry
    names = :($(capture))
    products = esc(Expr(:tuple, capture...))

    # We have to create this dummy_source, because we cannot, in a single line,
    # have both `@__FILE__` and `__source__` interpreted by the same julia.
    dummy_source = VERSION >= v"0.7.0-" ? __source__.file : ""

    return quote
        # First pick up important pieces of information from the call-site
        const source = VERSION >= v"0.7.0-" ? $("$(dummy_source)") : @__FILE__
        const depsjl_path = joinpath(dirname(source), "deps.jl")
        const package_name = basename(dirname(dirname(source)))
        const platform = platform_key()
        escape_path = path -> replace(path, "\\", "\\\\")

        const rebuild = strip("""
        Please re-run Pkg.build(\\\"$(package_name)\\\"), and restart Julia.
        """)

        # Begin by ensuring that we can satisfy every product RIGHT NOW
        for product in $(products)
            # Check to make sure that we've passed in the right kind of
            # objects, e.g. subclasses of `Product`
            if !(typeof(product) <: Product)
                msg = "Cannot @write_deps_file for $product, which is " *
                        "of type $(typeof(product)), which is not a " *
                        "subtype of `Product`!"
                error(msg)
            end

            if !satisfied(product; verbose=true)
                error("$product is not satisfied, cannot generate deps.jl!")
            end
        end

        # If things look good, let's generate the `deps.jl` file
        open(depsjl_path, "w") do depsjl_file
            # First, dump the preamble
            println(depsjl_file, strip("""
            ## This file autogenerated by BinaryProvider.@write_deps_file.
            ## Do not edit.
            """))

            # Next, spit out the paths of all our products
            for idx in 1:$(length(capture))
                product = $(products)[idx]
                name = $(names)[idx]

                # Escape the location so that e.g. Windows platforms are happy
                # with the backslashes in a string literal
                escaped_path = escape_path(locate(product, platform=platform))
                println(depsjl_file, strip("""
                const $(name) = \"$(escaped_path)\"
                """))
            end

            # Next, generate a function to check they're all on the up-and-up
            println(depsjl_file, "function check_deps()")

            for idx in 1:$(length(capture))
                product = $(products)[idx]
                name = $(names)[idx]

                # Add a `global $(name)`
                println(depsjl_file, "    global $(name)");

                # Check that any file exists
                println(depsjl_file, """
                    if !isfile($(name))
                        error("\$($(name)) does not exist, $(rebuild)")
                    end
                """)

                # For Library products, check that we can dlopen it:
                if typeof(product) <: LibraryProduct
                    println(depsjl_file, """
                        if Libdl.dlopen_e($(name)) == C_NULL
                            error("\$($(name)) cannot be opened, $(rebuild)")
                        end
                    """)
                end
            end

            # Close the `check_deps()` function
            println(depsjl_file, "end")
        end
    end
end