text:
# This file is a part of UnsafeArrays.jl, licensed under the MIT License (MIT).

using Base: @propagate_inbounds
using Base.MultiplicativeInverses: SignedMultiplicativeInverse

"""
    UnsafeArray{T,N} <: DenseArray{T,N}

An `UnsafeArray` is an bitstype wrapper around a memory pointer and a
size tuple. It's an intended to be used as a short-lived, allocation-free
alternative to an an `Array` returned by `unsafe_wrap` or a `SubArray`
returned by `view` or . Use with caution!

Constructors:

    UnsafeArray{T,N}(pointer::Ptr{T}, size::NTuple{N,Int}) where {T,N}
    UnsafeArray(pointer::Ptr{T}, size::NTuple{N,Int}) where {T,N}

UnsafeArray requires `isbitstype(T) == true`.

Note: It's safe to construct an empty multidimensional `UnsafeArray`:

    UnsafeArray(Ptr{Int}(0), (0,...))

Note: You *must* ensure that `A` is not garbage collected or reallocated
via (e.g.) `resize!`, `sizehint!` etc. while `U` is in use! Use only in
situations where you have full control over the life cycle of `A` and `U`.

`deepcopy(::UnsafeArray)` returns a standard `Array`.
"""
struct UnsafeArray{T,N} <: DenseArray{T,N}
    pointer::Ptr{T}
    size::NTuple{N,Int}

    UnsafeArray{T,N}(isbits_T::Val{true}, pointer::Ptr{T}, size::NTuple{N,Int}) where {T,N} =
        new{T,N}(pointer, size)
end

export UnsafeArray

UnsafeArray{T,N}(pointer::Ptr{T}, size::NTuple{N,Int}) where {T,N} =
    UnsafeArray{T,N}(Val{isbitstype(T)}(), pointer, size)

UnsafeArray(pointer::Ptr{T}, size::NTuple{N,Int}) where {T,N} =
    UnsafeArray{T,N}(pointer, size)

Base.size(A::UnsafeArray) = A.size

Base.elsize(::Type{UnsafeArray{T, N}}) where {T, N} = sizeof(T)

@inline function Base.getindex(A::UnsafeArray, i::Integer)
    @boundscheck checkbounds(A, i)
    unsafe_load(A.pointer, i)
end

@inline function Base.setindex!(A::UnsafeArray, x, i::Integer)
    @boundscheck checkbounds(A, i)
    unsafe_store!(A.pointer, x, i)
end

@inline Base.IndexStyle(A::UnsafeArray) = IndexLinear()

Base.unsafe_convert(::Type{Ptr{T}}, A::UnsafeArray{T}) where T = A.pointer

Base.iscontiguous(::UnsafeArray) = true
Base.iscontiguous(::Type{<:UnsafeArray}) = true

Base.@propagate_inbounds Base.view(A::UnsafeArray) = Base.unsafe_view(A)

Base.@propagate_inbounds function Base.view(A::UnsafeArray, idx, I...)
    J = Base.to_indices(A, (idx, I...))
    @boundscheck checkbounds(A, J...)
    Base.unsafe_view(A, J...)
end

Base.@propagate_inbounds Base.unsafe_view(A::UnsafeArray{T,N}, I::Vararg{Base.ViewIndex,N}) where {T,N} =
    _unsafe_view_impl((), A, I...)

Base.@propagate_inbounds Base.unsafe_view(A::UnsafeArray{T,N}, i::Base.ViewIndex) where {T,N} =
    _unsafe_view_impl((), A, i)

@inline function _unsafe_view_impl(IFwd::NTuple{N,Base.Slice}, A::UnsafeArray{T,N}) where {T,N}
    @assert IFwd == axes(A)
    A
end

Base.@propagate_inbounds _unsafe_view_impl(IFwd::NTuple{N,Base.ViewIndex}, A::UnsafeArray{T,N}) where {T,N} =
    SubArray(A, IFwd)

Base.@propagate_inbounds _unsafe_view_impl(IFwd::NTuple{M,Base.ViewIndex}, A::UnsafeArray{T,N}, i::Base.ViewIndex, I::Base.ViewIndex...) where {T,M,N} =
    _unsafe_view_impl((IFwd..., i), A, I...)

@inline function _unsafe_view_impl(IFwd::NTuple{M,Base.Slice}, A::UnsafeArray{T,N}, i::DenseIdx, I::Integer...) where {T,M,N}
    @assert IFwd == ntuple(i -> axes(A)[i], Val{M}())
    I_all = (IFwd..., i, I...)
    @boundscheck checkbounds(A, I_all...)
    startidxs = map(first, (IFwd..., i, I...))
    sub_s = _sub_size(I_all...)
    p = pointer(A, LinearIndices(size(A))[startidxs...])
    UnsafeArray(p, sub_s)
end

@inline function _unsafe_view_impl(IFwd::Tuple{}, A::UnsafeArray{T,N}, i::DenseIdx) where {T,N}
    @boundscheck checkbounds(A, i)
    p = pointer(A, first(i))
    sub_s = _sub_size(i)
    UnsafeArray(p, sub_s)
end

@inline Base.__reshape(p::Tuple{UnsafeArray,IndexLinear}, dims::Dims) =
    UnsafeArray(p[1].pointer, dims)

# From Julia Base (same implementation, with slight variations):

Base.copyto!(dest::Array{T}, src::UnsafeArray{T}) where {T} = copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::UnsafeArray{T}, soffs::Integer, n::Integer) where {T}
    n == 0 && return dest
    n > 0 || throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest)
        throw(BoundsError())
    end
    unsafe_copyto!(dest, doffs, src, soffs, n)
end

function Base.unsafe_copyto!(dest::Array{T}, doffs::Integer, src::UnsafeArray{T}, soffs::Integer, n::Integer) where {T}
    @assert isbitstype(T)
    unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
    return dest
end

Base.copyto!(dest::UnsafeArray{T}, src::Array{T}) where {T} = copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::UnsafeArray{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) where {T}
    n == 0 && return dest
    n > 0 || throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest)
        throw(BoundsError())
    end
    unsafe_copyto!(dest, doffs, src, soffs, n)
end

function Base.unsafe_copyto!(dest::UnsafeArray{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) where {T}
    @assert isbitstype(T)
    unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
    return dest
end

Base.deepcopy(A::UnsafeArray) = copyto!(similar(A), A)

# # Defining Base.unaliascopy results in very bad broadcast performance for
# # some reason, even when it shouldn't be called. By default, unaliascopy
# # results in an error for UnsafeArray.
#
# Base.unaliascopy(A::UnsafeArray) = begin
#     copy(A)
# end


text:
module MCHammer

#MC Hammer is a Monte-Carlo Tool Kit for quickly building/porting risk and decisision analysis in Julia for breakneck performance.

using LinearAlgebra
using Statistics, Random, Distributions, StatsBase, DataFrames, CSV, IterableTables
using TimeSeries
using Gadfly

import StatsBase: tiedrank
import Random: rand, rand!, seed!
import LinearAlgebra: BlasReal, BlasFloat
import Statistics: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
                   cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
median, middle

export
    cormat,
    covmat,
    corvar,
    GetCertainty,
    fractiles,
    cmd,
    truncate_digit,
    density_chrt,
    histogram_chrt,
    sensitivity_chrt,
    trend_chrt,
    GBMMfit,
    GBMM,
    GBMA_d,
    sip2csv,
    importxlsip,
    importsip,
    genmeta,
    marty,
    markov_a,
    markov_ts,
    RiskEvent

#Source files for MC_Hammer functions
include("mch_simtools.jl")
include("mch_charts.jl")
include("mch_timeseries.jl")
include("mch_SIPTools.jl")
include("mch_distributions.jl")
include("mch_stochastic.jl")

function test()
    println("Hello")

end
end # module

ext: jl


text:
using Yao
using QuDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using BitBasis
using Random
using Test
using YaoBlocks

#Linear Diff Equation Unitary M
function diffeqProblem(nbit::Int)
    siz = 1<<nbit
    Au = rand_unitary(siz)
    An = rand(ComplexF64,siz,siz)
    b = normalize!(rand(ComplexF64, siz))
    x = normalize!(rand(ComplexF64, siz))
    Au,An,b, x
end

@testset "QuLDE_Test" begin
    Random.seed!(2)
    N = 1
    k = 3
    tspan = (0.0,0.4)
    Au,An,b,x = diffeqProblem(N)
    qprob = QuLDEProblem(Au, b, x, tspan)
    f(u,p,t) = Au*u + b;
    prob = ODEProblem(f, x, tspan)

    sol = solve(prob, Tsit5(), dt = 0.1, adaptive = false)
    s = sol.u[end]

    out = solve(qprob, QuLDE(k))

    @test isapprox.(s, out, atol = 0.01) |> all

    qprob = QuLDEProblem(An, b, x, tspan)
    f(u,p,t) = An*u + b;
    prob = ODEProblem(f, x, tspan)

    sol = solve(prob, Tsit5(), dt = 0.1, adaptive = false)
    s = sol.u[end]

    out = solve(qprob, QuLDE(k))
    @test isapprox.(s, out, atol = 0.02) |> all

    # u0 equal to zero
    tspan = (0.0,0.1)

    qprob = QuLDEProblem(Au,b,tspan)
    out = solve(qprob,QuLDE(k))
    t = tspan[2] - tspan[1]
    r_out = (exp(Au*t) - Diagonal(ones(length(b))))*Au^(-1)*b
    @test isapprox.(r_out, out, atol = 1e-3) |> all

    qprob = QuLDEProblem(An,b,tspan)
    out = solve(qprob,QuLDE(k))
    t = tspan[2] - tspan[1]
    r_out = (exp(An*t) - Diagonal(ones(length(b))))*An^(-1)*b
    @test isapprox.(r_out, out, atol = 1e-3) |> all
end


text:
include("DataDictionary.jl")
include("LoadDictionaries.jl")
include("LoadSynthData.jl")
include("FluxDriver.jl")
include("Calculate_constraints.jl")
include("Bounds.jl")
include("TXTLDictionary.jl")
include("Utility.jl")
include("CalcError.jl")
include("Plot.jl")
include("PlotSeparate.jl")
include("PlotErrorbar.jl")
include("dFBA.jl")
include("RunFBA.jl")
using PyPlot
using Interpolations
using LaTeXStrings
using GLPK
term_out(GLPK.OFF)


text:
module TestMatrix
    using NullableArrays
    using Base.Test

    #----- test Base.diag -----#
    A = reshape([1:25...], 5, 5)
    m = fill(false, 5, 5)
    m[1] = true
    m[25] = true
    X = NullableArray(A)
    Y = NullableArray(A, m)

    @test isequal(diag(X), NullableArray([1, 7, 13, 19, 25]))
    @test isequal(diag(Y), NullableArray([1, 7, 13, 19, 25],
                                         [true, false, false, false, true]))

end


text:
#######################################################################################################################################################################################################
#
# Changes to the function
# General:
#     2021-Dec-24: move the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Jan-28: remove the complicated funtion to create var and dim at the same time
#
#######################################################################################################################################################################################################
"""
NCDatasets.jl does not have a convenient function (1 line command) to append or grow dataset into a file. Thus, we provide a few methods as supplements:

$(METHODLIST)

"""
function append_nc! end

#######################################################################################################################################################################################################
#
# Changes to the method
# General:
#     2022-Jan-28: add this method to add data to Dataset
#     2022-Jan-28: add documentation
#
#######################################################################################################################################################################################################
"""

    append_nc!(ds::Dataset, var_name::String, var_data::Array{T,N}, var_attributes::Dict{String,String}, dim_names::Vector{String}; compress::Int = 4) where {T<:Union{AbstractFloat,Int,String},N}

Append data to existing netcdf dataset, given
- `ds` A `NCDatasets.Dataset` type dataset
- `var_name` New variable name to write to
- `var_data` New variable data to write, can be integer, float, and string with N dimens
- `var_attributes` New variable attributes
- `dim_names` Dimension names in the netcdf file
- `compress` Compression level fro NetCDF, default is 4

---
# Examples
```julia
create_nc!("test.nc", String["lon", "lat", "ind"], [36, 18, 5]);
dset = Dataset("test.nc", "a");
append_nc!(dset, "str", ["A" for i in 1:18], Dict("longname" => "test strings"), ["lat"]);
append_nc!(dset, "lat", collect(1:18), Dict("longname" => "latitude"), ["lat"]);
append_nc!(dset, "lon", collect(1:36), Dict("longname" => "longitude"), ["lon"]; compress=4);
append_nc!(dset, "ind", collect(1:5), Dict("longname" => "index"), ["ind"]);
append_nc!(dset, "d2d", rand(36,18), Dict("longname" => "a 2d dataset"), ["lon", "lat"]);
append_nc!(dset, "d3d", rand(36,18,5), Dict("longname" => "a 3d dataset"), ["lon", "lat", "ind"]);
close(dset);
```
"""
append_nc!(ds::Dataset, var_name::String, var_data::Array{T,N}, var_attributes::Dict{String,String}, dim_names::Vector{String}; compress::Int = 4) where {T<:Union{AbstractFloat,Int,String},N} = (
    # only if variable does not exist create the variable
    @assert !(var_name in keys(ds)) "You can only add new variable to the dataset!";
    @assert length(dim_names) ==  N "Dimension must be match!";
    @assert 0 <= compress <= 9 "Compression rate must be within 0 to 9";

    _ds_var = defVar(ds, var_name, T, dim_names; attrib = var_attributes, deflatelevel = compress);
    _ds_var[:,:] = var_data;

    return nothing
);

#######################################################################################################################################################################################################
#
# Changes to the method
# General:
#     2021-Dec-24: add the case to write variables to the file directly
#     2022-Jan-28: add documentation
#
#######################################################################################################################################################################################################
"""

    append_nc!(file::String, var_name::String, var_data::Array{T,N}, var_attributes::Dict{String,String}, dim_names::Vector{String}; compress::Int = 4) where {T<:Union{AbstractFloat,Int,String},N}

Append data to existing netcdf dataset, given
- `file` Path of the netcdf dataset
- `var_name` New variable name to write to
- `var_data` New variable data to write, can be integer, float, and string with N dimens
- `var_attributes` New variable attributes
- `dim_names` Dimension names in the netcdf file
- `compress` Compression level fro NetCDF, default is 4

---
# Examples
```julia
create_nc!("test.nc", String["lon", "lat", "ind"], [36, 18, 5]);
append_nc!("test.nc", "str", ["A" for i in 1:18], Dict("longname" => "test strings"), ["lat"]);
append_nc!("test.nc", "lat", collect(1:18), Dict("longname" => "latitude"), ["lat"]);
append_nc!("test.nc", "lon", collect(1:36), Dict("longname" => "longitude"), ["lon"]; compress=4);
append_nc!("test.nc", "ind", collect(1:5), Dict("longname" => "index"), ["ind"]);
append_nc!("test.nc", "d2d", rand(36,18), Dict("longname" => "a 2d dataset"), ["lon", "lat"]);
append_nc!("test.nc", "d3d", rand(36,18,5), Dict("longname" => "a 3d dataset"), ["lon", "lat", "ind"]);
```
"""
append_nc!(file::String, var_name::String, var_data::Array{T,N}, var_attributes::Dict{String,String}, dim_names::Vector{String}; compress::Int = 4) where {T<:Union{AbstractFloat,Int,String},N} = (
    _dset = Dataset(file, "a");
    append_nc!(_dset, var_name, var_data, var_attributes, dim_names; compress = compress);
    close(_dset);

    return nothing
);


text:
############################################################
# L2 regularization
############################################################
function forward(backend::GPUBackend, regu :: L2Regu, global_regu::AbstractFloat, param :: Blob)
  return regu.coefficient * global_regu * CuBLAS.dot(backend.cublas_ctx, eltype(param), length(param),
      param.ptr, 1, param.ptr, 1)
end
function backward(backend::GPUBackend, regu :: L2Regu, global_regu::AbstractFloat, param :: Blob, gradient :: Blob)
    CuBLAS.axpy(backend.cublas_ctx, length(param),
        convert(eltype(param), 2 * regu.coefficient * global_regu), param.ptr, 1, gradient.ptr, 1)
end

############################################################
# L1 regularization
############################################################
function forward(backend::GPUBackend, regu :: L1Regu, global_regu::AbstractFloat, param :: Blob)
  loss_blob = make_zero_blob(backend, Float32, 1, 1, 1, 1)
  len = length(param)
  coef = convert(eltype(param), regu.coefficient * global_regu)
  x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
  if eltype(param) == Float32
    kernel = backend.mocha.l1_forward_float
  else
    kernel = backend.mocha.l1_forward_double
  end
  CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X, (param.ptr.p, len, coef, loss_blob.ptr.p))
  loss = Float32[0]
  copy!(loss, loss_blob)
  return loss[1]
end
function backward(backend::GPUBackend, regu :: L1Regu, global_regu::AbstractFloat, param :: Blob, gradient :: Blob)
  len = length(param)
  x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
  coef = convert(eltype(param), regu.coefficient * global_regu)
  if eltype(param) == Float32
    kernel = backend.mocha.l1_backward_float
  else
    kernel = backend.mocha.l1_backward_double
  end
  CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X, (param.ptr.p, gradient.ptr.p, len, coef))
end


text:
function schsolve_Emin_cg( pw::PWGrid, Vpot, psi::Array{ComplexF64,2};
                           NiterMax=1000 )
    #
    Npoints = size(psi)[1]
    Nstates = size(psi)[2]

    g = zeros(ComplexF64, Npoints, Nstates)
    Kg = zeros(ComplexF64, Npoints, Nstates)
    Kg_old = zeros(ComplexF64, Npoints, Nstates)
    psic = zeros(ComplexF64, Npoints, Nstates )

    d = zeros(ComplexF64, Npoints, Nstates)
    g_old = zeros(ComplexF64, Npoints, Nstates)
    d_old = zeros(ComplexF64, Npoints, Nstates)
    gt = zeros(ComplexF64, Npoints, Nstates)

    yk1 = zeros(ComplexF64, Npoints, Nstates) # for testing Hager-Zhang beta
    v1 = zeros(ComplexF64, Npoints, Nstates)

    alphat = 1.e-5
    beta = 0.0
    Etot_old = 0.0
    Etot = 0.0
    #
    for iter = 1:NiterMax

        g = calc_grad( pw, Vpot,  psi)
        Kg = Kprec(pw,g)

        if iter != 1

            #beta = real(tr(g'*Kg))/real(tr(g_old'*Kg_old))
            #beta = real(tr((g-g_old)'*Kprec(pw,g)))/real(tr(g_old'*Kprec(pw,g_old)))
            #beta = real(tr((g-g_old)'*Kg))/real(tr((g-g_old)'*d))

            yk1 = g - g_old
            tt = real(tr(d_old'*yk1))
            v1 = yk1 - 2*d_old*norm(yk1)/tt
            beta = real(tr(v1'*g)/tt)
        end
        #println("beta HZ = ", beta)

        if beta < 0.0
            println("beta is negative, resetting beta to zero")
            beta = 0.0
        end

        d = -Kprec(pw, g) + beta * d_old
        psic = ortho_gram_schmidt(psi + alphat*d)
        gt = calc_grad( pw, Vpot, psic )

        denum = real(tr((g-gt)'*d))
        if denum != 0.0
          alpha = abs(alphat*real(tr(g'*d))/denum )
        else
          alpha = 0.0
        end

        psi = psi[:,:] + alpha*d[:,:]

        psi = ortho_gram_schmidt(psi)
        Etot = calc_Etot( pw, Vpot, psi )

        diff = abs(Etot-Etot_old)
        @printf("E step %8d = %18.10f %18.10e\n", iter, Etot, diff)
        if diff < 1e-6
            @printf("CONVERGENCE ACHIEVED\n")
            break
        end
        g_old = copy(g)
        d_old = copy(d)
        Kg_old = copy(Kg)
        Etot_old = Etot
    end
    return psi, Etot
end


text:
using DEBTrait, Test

################################################################################
# test BatchC
###################
n_polymers      = 0
n_monomers      = 1
n_microbes      = 1
n_enzymes       = 1
n_minerals      = 0
#
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(n_monomers,n_microbes), 1.0*ones(n_microbes))
Setup           = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
Monomer         = BaseMonomer(1e-10*ones(n_monomers), 6*ones(n_monomers), 3*ones(n_monomers), 0.3*ones(n_monomers))
Soil            = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC     = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
AssimilationC   = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
p = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
p.PMetabolismC.y_EM = ones(n_microbes)
p.PMetabolismC.y_EX = ones(n_microbes)
p.PMetabolismC.α_n  = ones(n_microbes)*200.0
p.I_D = 1e-3*ones(n_monomers)
#
du = zeros(p.dim)
u  = 1e-8*ones(p.dim)
t = 0.0
#
out = rhs_BatchC!(du,u,p,t)
mass_balance = sum(out) .- p.I_D
@test mass_balance[1] <= 1e-10

###################
n_polymers      = 0
n_monomers      = 2
n_microbes      = 1
n_enzymes       = 1
n_minerals      = 0
#
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(n_monomers,n_microbes))
Setup           = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
Monomer         = BaseMonomer(1e-10*ones(n_monomers), 6*ones(n_monomers), 3*ones(n_monomers), 0.3*ones(n_monomers))
Soil            = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC     = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
AssimilationC   = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
p = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
p.PMetabolismC.y_EM = ones(n_microbes)
p.PMetabolismC.y_EX = ones(n_microbes)
p.PMetabolismC.α   = zeros(n_microbes)
p.I_D = 1e-3*ones(n_monomers)
#
du = zeros(p.dim)
u  = 1e-8*ones(p.dim)
t = 0.0
#
out = rhs_BatchC!(du,u,p,t)
mass_balance = sum(out) .- sum(p.I_D)
@test mass_balance[1] <= 1e-10

###################
n_polymers      = 0
n_monomers      = 1
n_microbes      = 2
n_enzymes       = 1
n_minerals      = 0
#
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(n_monomers,n_microbes))
Setup           = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
Monomer         = BaseMonomer(1e-10*ones(n_monomers), 6*ones(n_monomers), 3*ones(n_monomers), 0.3*ones(n_monomers))
Soil            = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC     = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
AssimilationC   = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
p = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
p.PMetabolismC.y_EM = ones(n_microbes)
p.PMetabolismC.y_EX = ones(n_microbes)
p.PMetabolismC.α   = zeros(n_microbes)
p.I_D = 1e-3*ones(n_monomers)
#
du = zeros(p.dim)
u  = 1e-8*ones(p.dim)
t = 0.0
#
out = rhs_BatchC!(du,u,p,t)
mass_balance = sum(out) .- sum(p.I_D)
@test mass_balance[1] <= 1e-10

#################
n_polymers      = 0
n_monomers      = 2
n_microbes      = 2
n_enzymes       = 1
n_minerals      = 0
#
Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(n_monomers,n_microbes))
Setup           = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
Monomer         = BaseMonomer(1e-10*ones(n_monomers), 6*ones(n_monomers), 3*ones(n_monomers), 0.3*ones(n_monomers))
Soil            = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC     = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
AssimilationC   = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
p = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
p.PMetabolismC.y_EM = ones(n_microbes)
p.PMetabolismC.y_EX = ones(n_microbes)
p.PMetabolismC.α   = zeros(n_microbes)
p.I_D = 1e-3*ones(n_monomers)
#
du = zeros(p.dim)
u  = 1e-8*ones(p.dim)
t = 0.0
#
out = rhs_BatchC!(du,u,p,t)
mass_balance = sum(out) .- sum(p.I_D)
@test mass_balance[1] <= 1e-10

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


text:
using Queryverse
using Test

@testset "Queryverse" begin

    df = load(joinpath(@__DIR__, "testdata.csv")) |>
@query(i, begin
        @select {i.Count, i.Year}
    end) |>
@tee(save("testoutput.csv")) |>
# @tee(begin
#     data_values() |>
#     mark_point() |>
#     encoding_x_ord(:Year) |>
#     encoding_y_quant(:Count) |>
#     save("testoutput.pdf")
# end) |>
DataFrame

    @test size(df) == (2, 2)

    io = IOBuffer()
    f(x) = print(io, sum(x + x))
    1:10 |> @tee(f) |> x->print(io, " ", sum(x))
    @test String(take!(io)) == "110 55"

end


text:
# equations.jl

# This file contains all implemented equations solved in the model

# Currently implemented:
# - Pressure equation
# - Darcy equation (velocity)

# This is the present data storage:
# grid       = (nx, ny, nz, Δx, Δy, Δz, x, y, z)
# state      = (h, u, v, w, hⁿ⁺¹, uⁿ⁺¹, vⁿ⁺¹, wⁿ⁺¹)
# model      = externals, recharge
# parameters = (K, Δt, tend)
# externals  = source
# source     = (i_src, j_src, k_src, n_src, externals)

function pressureEquation!(grid, model, state, parameters, timeData)

    source = model.source.externalSource
    rchFlux = model.recharge.rechargeFlux

    # Solve for the pressure/head
    deps = Event(device(grid.arch))

    toplayerNr = grid.nz
    event1 = launch!(
                    grid,               # grid
                    :xy,                # coverage
                    addToLayer!,        # kernel
                    source,             # kernel args
                    rchFlux,
                    toplayerNr,
                    dependencies=deps   # handle for sync.
                    )

    event2 = launch!(
                    grid,                # grid
                    :xyz,                # coverage
                    pressureKernel!,  # kernel
                    source,              # kernel args...
                    state.h,
                    state.hⁿ⁺¹,
                    grid.Δx,
                    grid.Δy,
                    grid.Δz,
                    timeData.Δt,
                    parameters.K,
                    parameters.specificStorage,
                    dependencies=event1  # handle for sync.
                    )

    wait(device(grid.arch), event2)

end

function darcyEquation!(grid, model, state, parameters, timeData)
    # Compute the velocities from the darcy equation using the pressure
    deps = Event(device(grid.arch))
    event = launch!(
                    grid,                # grid
                    :xyz,                # coverage
                    darcyKernel!,     # kernel
                    state.uⁿ⁺¹,          # kernel args...
                    state.vⁿ⁺¹,
                    state.wⁿ⁺¹,
                    state.h,
                    grid.Δx,
                    grid.Δy,
                    grid.Δz,
                    parameters.K,
                    dependencies=deps # event handle for sync.
                    )
    wait(device(grid.arch), event)
end


text:
# If there is no time dimension we return the same data for every timestep
_auxval(data::AbstractSimData, key::Union{Aux,Symbol}, I...) =
    _auxval(aux(data, key), data, key, I...)
_auxval(A::AbstractMatrix, data::SimData, key::Union{Aux,Symbol}, y, x) = A[y, x]
# function _auxval(A::AbstractDimArray{<:Any,2}, data::SimData, key::Union{Aux,Symbol}, y, x)
    # X = DD.basetypeof(dims(A, XDim))
    # Y = DD.basetypeof(dims(A, YDim))
    # A[y, x]
# end
function _auxval(A::AbstractDimArray{<:Any,3}, data::SimData, key::Union{Aux,Symbol}, y, x)
    # X = DD.basetypeof(dims(A, XDim))
    # Y = DD.basetypeof(dims(A, YDim))
    A[y, x, auxframe(data, key)]
end

function boundscheck_aux(data::SimData, auxkey::Aux{Key}) where Key
    a = aux(data)
    (a isa NamedTuple && haskey(a, Key)) || _auxmissingerror(Key, a)
    gsize = gridsize(data)
    asize = size(a[Key], 1), size(a[Key], 2)
    if asize != gsize
        _auxsizeerror(Key, asize, gsize)
    end
    return true
end

@noinline _auxmissingerror(key, aux) = error("Aux data $key is not present in aux")
@noinline _auxsizeerror(key, asize, gsize) =
    error("Aux data $key is size $asize does not match grid size $gsize")

_calc_auxframe(data) = _calc_auxframe(aux(data), data)
_calc_auxframe(aux::NamedTuple, data) = map(A -> _calc_auxframe(A, data), aux)
function _calc_auxframe(A::AbstractDimArray, data)
    hasdim(A, Ti) || return nothing
    curtime = currenttime(data)
    firstauxtime = first(dims(A, TimeDim))
    auxstep = step(dims(A, TimeDim))
    # Use julias range objects to calculate the distance between the
    # current time and the start of the aux
    i = if curtime >= firstauxtime
        length(firstauxtime:auxstep:curtime)
    else
        1 - length(firstauxtime-timestep(data):-auxstep:curtime)
    end
    _cyclic_index(i, size(A, 3))
end
_calc_auxframe(aux, data) = nothing

function _cyclic_index(i::Integer, len::Integer)
    return if i > len
        rem(i + len - 1, len) + 1
    elseif i <= 0
        i + (i ÷ len -1) * -len
    else
        i
    end
end


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

const KKT = LDLFactorizations.LDLFactorization{Float64,Int64}

"""
    AdmmParams

Parameters to the ADMM algorithm.

The fields represent:
- `P`: A sparse `n`x`n` positive semidefinite matrix
- `q`: an `n`-vector
- `A`: an `m`x`n` matrix
- `b`: an `m`-vector
- `g`: a vector of PiecewiseQuadratic functions
"""
mutable struct AdmmParams
    P::SparseMatrixCSC{Float64,Int64}
    q::Vector{Float64}
    A::SparseMatrixCSC{Float64,Int64}
    b::Vector{Float64}
    g::Vector{PiecewiseQuadratic}

    function AdmmParams(P, q, A, b, g)
        m, n = size(A)
        @assert size(P) == (n, n)
        @assert size(q) == (n,)
        @assert size(A) == (m, n)
        @assert size(b) == (m,)
        @assert size(g) == (n,)
        return new(P, q, A, b, g)
    end
end

function Base.show(io::IO, params::AdmmParams)
    print(io, "AdmmParams")
    print(io, "  P: ")
    print(io, params.P)
    print(io, "  q: ")
    print(io, params.q)
    print(io, "  A: ")
    print(io, params.A)
    print(io, "  b: ")
    print(io, params.b)
    print(io, "  g: ")
    print(io, params.g)
    return
end

"""
    get_num_vars(params::AdmmParams)

Return the number of variables in the problem that `params` was initialized with.
"""
get_num_vars(params::AdmmParams) = size(params.A, 2)

"""
    get_num_constrs(params::AdmmParams)

Return the number of constraints in the problem that `params` was initialized with.
"""
get_num_constrs(params::AdmmParams) = size(params.A, 1)

"""
    Settings

Settings for the ADMM algorithm.

The fields represent:
- `ρ`: Augmented Lagrangian parameter for constraints (must be positive)
- `σ`: Augmented Lagrangian parameter for variables (must be positive)
- `α`:  Over-relaxation parameter (must be ∈ [0, 2], typically ∈ [1.0, 1.8])
- `max_iters`: Maximum allowed ADMM iterations (must be positive)
- `ϵ`: Threshold for determining convergence (must be nonnegative)
- `term_cond_freq`: Frequency in iterations with which to check for termination
- `compute_stats`: Whether to compute objectives and residuals on every iteration
- `term_cond_type`:
    - `1` uses the `ConvergeTermCache`
    - `2` uses the `FirstVarsTermCache`
- `obj_tol`: Minimum amount by which the objective needs to improve to be considered and improvement
- `res_tol`: Maximum allowable residual for an iterate to be considered an improvement
- `non_improvement_iters`: Maximum allowed non-improvement iterations before termination
"""
@with_kw mutable struct Settings
    ρ::Vector{Float64}
    σ::Vector{Float64}
    α::Float64 = 1.0
    max_iters::Int64 = 1000
    ϵ::Float64 = 1e-4
    term_cond_freq::Int = 10
    compute_stats::Bool = false
    term_cond_type::Int64 = 1
    obj_tol::Float64 = 1e-5
    res_tol::Float64 = 3e-4
    non_improvement_iters::Int64 = 50

    @assert (α > 0) & (α < 2)
    @assert all(ρ .> 0)
    @assert all(σ .> 0)
    @assert max_iters > 0
    @assert ϵ >= 0
    @assert obj_tol > 0
    @assert res_tol > 0
    @assert non_improvement_iters > 0
end

"""
    Settings(n:: Int, m:: Int)

Construct ADMM Settings with default values given the number of variables `n` and constraints `m`
"""
Settings(n::Int, m::Int) = Settings(; ρ=ones(m), σ=ones(n))

"""
    Settings(n:: Int, m:: Int)

Construct ADMM Settings with default values in accordance with set of ADMM Params `params`
"""
function Settings(params::AdmmParams)
    n = get_num_vars(params::AdmmParams)
    m = get_num_constrs(params::AdmmParams)
    return Settings(n::Int, m::Int)
end

"""
    get_num_vars(settings::Settings)

Return the number of variables in the problem that `settings` was initialized with.
"""
get_num_vars(settings::Settings) = length(settings.σ)

"""
    get_num_constrs(settings::Settings)

Return the number of constraints in the problem that `settings` was initialized with.
"""
get_num_constrs(settings::Settings) = length(settings.ρ)

"""
    Stats

Stats from the ADMM algorithm per iteration

The fields represent:
- `obj`: Objective values
- `res`: Residuals
- `iters`: Number of iterations completed
"""
@with_kw mutable struct Stats
    obj::Vector{Float64}
    res::Vector{Float64}
    iters::Int64

    function Stats(max_iters::Int64)
        obj = fill(NaN, max_iters)
        res = fill(NaN, max_iters)
        return new(obj, res, 0)
    end
end

Stats(settings::Settings) = Stats(settings.max_iters)

"""
    Vars

Store variables that need to be updated on each iteration of ADMM.

The fields represent:
- `x`: Previous value of x
- `z`: Previous value of z
- `w`: Previous value of w
- `y`: Previous value of y
- `xt`: Auxiliary variable
- `zt`: Auxiliary variable

# Notes:
* We keep the previous values of `x`, `w`, `z`, `y` in order to check for convergence.
* The actual problem we are solving here is detailed in [`SeparableOptimization.factorize_kkt`](@ref)
* `w` and `y` are used when minimizing the augmented Lagrangian of the original problem (and consequently
  in ADMM's update rules)
"""
@with_kw mutable struct Vars
    x::Vector{Float64}
    z::Vector{Float64}
    w::Vector{Float64}
    y::Vector{Float64}
    xt::Vector{Float64}
    zt::Vector{Float64}

    @assert length(x) == length(xt)
    @assert length(x) == length(w)
    @assert length(z) == length(zt)
    @assert length(z) == length(y)
end

"""
    Vars(x::Vector{Float64}, z::Vector{Float64})

Construct ADMM Vars with the appropriate shapes given `x` and `z`
"""
function Vars(x::Vector{Float64}, z::Vector{Float64})
    n = length(x)
    m = length(z)
    return Vars(x, z, zeros(n), zeros(m), x, z)
end
"""
    Vars(n::Int, m::Int)

Construct ADMM Vars given the number of variables `n` and constraints `m`
"""
Vars(n::Int, m::Int) = Vars(zeros(n), zeros(m))
"""
    Vars(params::AdmmParams)

Construct ADMM Vars with the appropriate shapes given AdmmParams `params`
"""
function Vars(params::AdmmParams)
    n = get_num_vars(params)
    m = get_num_constrs(params)
    return Vars(n, m)
end

"""
    get_num_vars(vars::Vars)

Return the number of variables in the problem that `vars` was initialized with.
"""
get_num_vars(vars::Vars) = length(vars.x)

"""
    get_num_constrs(vars::Vars)

Return the number of constraints in the problem that `vars` was initialized with.
"""
get_num_constrs(vars::Vars) = length(vars.z)

"""
    assert_valid(params::AdmmParams, settings::Settings, vars::Vars)

Assert that the number of variables and constraints are consistent across a `params`-`settings`-`vars` combination
"""
function assert_valid(params::AdmmParams, settings::Settings, vars::Vars)
    n = get_num_vars(params)
    m = get_num_constrs(params)
    @assert n == get_num_vars(settings)
    @assert n == get_num_vars(vars)
    @assert m == get_num_constrs(settings)
    @assert m == get_num_constrs(vars)
end


text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset verbose=true "test_flows" begin

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 20    S1: 15  |----------------------|
      S2: 30    S2: 30  |         35           |
                        |                      |
                        |                      |
    (pilotable) prod_1_1|                      |(pilotable) prod_2_1
    Pmin=10, Pmax=100   |                      | Pmin=10, Pmax=100
    Csta=45k, Cprop=10  |                      | Csta=80k, Cprop=15
    ON                  |                      | ON
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10     S1: 15 |                      | S1: 40  S1: 50
      S2: 10     S2: 15 |                      | S2: 45  S2: 50

    =#

    TS = [DateTime("2015-01-01T11:00:00"), DateTime("2015-01-01T11:15:00")]
    ech = DateTime("2015-01-01T07:00:00")
    network = PSCOPF.Networks.Network()
    # Buses
    PSCOPF.Networks.add_new_bus!(network, "bus_1")
    PSCOPF.Networks.add_new_bus!(network, "bus_2")
    # Branches
    PSCOPF.Networks.add_new_branch!(network, "branch_1_2", 35.);
    # PTDF
    PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_1", 0.5)
    PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_2", -0.5)
    # Limitables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                            0., 100.,
                                            0., 1.,
                                            Dates.Second(0), Dates.Second(0))
    # Pilotables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                            0., 100.,
                                            0., 2.,
                                            Dates.Second(0), Dates.Second(0))
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                            0., 100.,
                                            0., 3.,
                                            Dates.Second(0), Dates.Second(0))
    # Uncertainties
    uncertainties = PSCOPF.Uncertainties()
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 20.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 30.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:15:00"), "S1", 15.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:15:00"), "S2", 30.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 10.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 10.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:15:00"), "S1", 15.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:15:00"), "S2", 15.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:00:00"), "S1", 40.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:00:00"), "S2", 45.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:15:00"), "S1", 50.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:15:00"), "S2", 50.)

    function create_schedule()
        schedule = PSCOPF.Schedule(PSCOPF.Utilitary(), DateTime("2015-01-01T06:00:00"), SortedDict(
            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                SortedDict(),
                SortedDict(TS[1]=> PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>20.,"S2"=>30.)),
                            TS[2] => PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>15.,"S2"=>30.)),
                            ),
                ),
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(),
                SortedDict(TS[1]=> PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>30.,"S2"=>25.)),
                            TS[2] => PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>50.,"S2"=>35.)),
                            ),
                ),
            "prod_2_1" => PSCOPF.GeneratorSchedule("prod_2_1",
                SortedDict(),
                SortedDict(TS[1]=> PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>0.,"S2"=>0.)),
                            TS[2] => PSCOPF.UncertainValue{Float64}(missing, SortedDict("S1"=>0.,"S2"=>0.)),
                            ),
                ),
            )
        )

        return  schedule
    end

    @testset "test_flows_simple" begin
        schedule = create_schedule()
        @test 40. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S1")
        @test 45. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S2")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S1")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S2")
    end

    @testset "test_flows_simple_change" begin
        schedule = create_schedule()
        PSCOPF.set_prod_value!(schedule, "prod_1_1", TS[1], "S1", 10.) # 30 -> 10 : -20
        PSCOPF.set_prod_value!(schedule, "prod_2_1", TS[1], "S1", 20.) # 0 -> 20 : +20

        # changed
        @test 20. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S1")

        # These are unchanged
        @test 45. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S2")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S1")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S2")
    end

    @testset "test_flows_slack" begin
        schedule = create_schedule()
        PSCOPF.set_prod_value!(schedule, "prod_1_1", TS[1], "S1", 10.) # 30 -> 10 : -20
        PSCOPF.set_prod_value!(schedule, "prod_2_1", TS[1], "S1", 10.) # 0 -> 10 : +10

        # slack effect
        @test ( PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S1")
                ≈ ( 0.5 * ( (20. + 10.) - (10.) ) #bus1
                   -0.5 * ( (10.) - (40.) ) #bus2
                )#25
            )

        # These are unchanged
        @test 45. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S2")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S1")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S2")
    end

    @testset "test_flows_capping" begin
        schedule = create_schedule()
        PSCOPF.set_capping_value!(schedule, "wind_1_1", TS[1], "S1", 5.)

        # capping 5MW from wind_1_1 => 20-5
        @test ( PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S1")
                ≈ ( 0.5 * ( ((20. - 5.) + 30.) - (10.) ) #bus1
                   -0.5 * ( (0.) - (40.) ) #bus2
                )#37.5
            )

        # These are unchanged
        @test 45. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S2")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S1")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S2")
    end

    @testset "test_flows_loss_of_load" begin
        schedule = create_schedule()
        PSCOPF.set_loss_of_load_value!(schedule, "bus_2", TS[1], "S1", 15.)

        # cut a 15MW conso on bus2 => 40-15
        @test ( PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S1")
                ≈ ( 0.5 * ( (20. + 30.) - (10.) ) #bus1
                   -0.5 * ( (0.) - (40. - 15.) ) #bus2
                )#32.5
            )

        # These are unchanged
        @test 45. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[1], "S2")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S1")
        @test 50. ≈ PSCOPF.compute_flow("branch_1_2", uncertainties, schedule, network, ech, TS[2], "S2")
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using LinearAlgebra
using Statistics
using PyCall
using Optim
using ASE
using JuLIP
using FiniteDifferences
using JuLIPMaterials.CLE
using JuLIPMaterials.Fracture

ase_units = pyimport("ase.units")

# Parameters

R = 50.0 # radius of domain
crack_surface = [0, 0, 1] # y-diretion, the surface which is opened
crack_front = [1, -1, 0] # z-direction, crack front line
relax_elasticity = false # if true, C_ij matrix computed with internal relaxation
relax_surface = true # if true, surface energy computed with internal relaxation

eam = JuLIP.Potentials.EAM(test_pots * "w_eam4.fs")

unitcell = bulk(:W, cubic=true)
@show energy(eam, unitcell)

variablecell!(unitcell)
set_calculator!(unitcell, eam)
minimise!(unitcell)
alat = mean(diag(cell(unitcell)))
@show alat
@test alat ≈ 3.143390027695279

# make a new bulk so it's perfectly cubic
unitcell = bulk(:W, cubic=true, a=alat)
set_calculator!(unitcell, eam)

# 6x6 elastic constant matrix
C = voigt_moduli(unitcell)
C11, C12, C44 = cubic_moduli(C, tol=1e-3)
@show C11, C12, C44
@test all(isapprox.((C11, C12, C44), (3.2644386308256323, 1.2618234185741113, 1.004134918612743), atol=1e-5))

bulk_at = bulk(:W, cubic=true, a=alat, y=crack_front, z=crack_surface)
slab = bulk_at * (1, 1, 10)
set_calculator!(slab, eam)
surface = deepcopy(slab)
X = positions(surface) |> mat
X[3, :] .+= 2.0
set_positions!(surface, X)
wrap_pbc!(surface)
c = Matrix(surface.cell)
c[3, :] += [0.0, 0.0, 10.0]
set_cell!(surface, c)
fixedcell!(surface)
set_calculator!(surface, eam)
area = norm(cross(slab.cell[:, 1], slab.cell[:, 2]))
γ = (energy(surface) - energy(slab)) / (2 * area)

@show γ / (ase_units.J / ase_units.m^2)
@test γ / (ase_units.J / ase_units.m^2) ≈ 2.950564005097423

# build the final crystal in correct orientation
r_cut = cutoff(eam) # cutoff distance for the potential

atu = bulk(:W, a=alat, y=crack_front, z=crack_surface)
cryst = cluster(atu, R, dims=[1,2], parity=[nothing,1,nothing])
X = positions(cryst) |> mat

rac = RectilinearAnisotropicCrack(PlaneStrain(), C11, C12, C44,
                                  crack_surface, crack_front)

k_G = k1g(rac, γ)
@show k_G

x0, y0, _ = diag(cell(cryst)) ./ 2

# check gradient wrt finite differnces
u, ∇u = u_CLE(rac, cryst, x0, y0)

X0 = copy(X)

U = [u(1.0, 0.0); zeros(length(cryst))']

clust = deepcopy(cryst)
X = X0 + U
set_positions!(clust, X)

@test maximum((central_fdm(2, 1))(α -> u(1.0, α), 0.0) - ∇u(1.0, 0.0)) < 1e-5

r = sqrt.((X0[1,:] .- x0).^2 + (X0[2,:] .- y0).^2)
set_calculator!(clust, eam)
F = norm.(forces(clust))

# exclude anything within cutoff of a surface atom
surface_atoms = findall([sum(norm.(R) .< 3.0) for (i,j,R) in sites(clust, 3.0)] .!= 8)
mask = [all([jj ∉ surface_atoms for jj in j]) for (i,j,R) in sites(clust, cutoff(eam))]

rm, Fm = r[mask], F[mask]

maxF = Float64[]
maxr = Float64[]
bins = range(minimum(rm), maximum(rm), length=10)
for (lower, upper) in zip(bins[1:end-1], bins[2:end])
   push!(maxr, (lower + upper)/2)
   push!(maxF, quantile(Fm[lower .< rm .< upper], 0.95))
end

# p = plot(; yscale=:log10, xscale=:log10, legend=:bottomleft, xlabel="r",
#             ylabel="|f|")
# scatter!(rm, Fm)
# plot!(r, 8*r.^(-2), label="r^-2", lw=2)
# plot!(maxr, maxF, marker=:o)

# fit a straight line to log-log plot to estimate decay
logR = log10.(maxr)
logF = log10.(maxF)
c = [ones(length(logR)) logR] \ logF
@show c
# plot!(10 .^ logR, 10 .^ (c[1] .+ c[2] * logR))

# decay should be at least r^{-2}
@test c[2] < -2


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# ----------------
# IMPLEMENTATIONS
# ----------------
include("statistics/data.jl")
include("statistics/ensembles.jl")


text:
function fibsq24()
#
#  test polynomial suggested by Goedecker
#  square of Fibocacci polynomial
#
    n = 24;
    p = fib(n)
    p = conv(p, p)
    z = [-.9558467190743565;
         -.9244083656460181-.2442575811931259*im;
         -.9244083656460181+.2442575811931259*im;
         -.8320677914606713-.4727784060949767*im;
         -.8320677914606713+.4727784060949767*im;
         -.6846250033032301-.6707962641470166*im;
         -.6846250033032301+.6707962641470166*im;
         -.4913433576989534-.8254218906657225*im;
         -.4913433576989534+.8254218906657225*im;
         -.2643730603203030-.9264240880183555*im;
         -.2643730603203030+.9264240880183555*im;
         -.1800084671534713e-1-.9668298465915332*im;
         -.1800084671534713e-1+.9668298465915332*im;
         .2322158170552186-.9432943625623267*im;
         .2322158170552186+.9432943625623267*im;
         .4703449187206489-.8561989307649443*im;
         .4703449187206489+.8561989307649443*im;
         .6808619003199168-.7094507040001145*im;
         .6808619003199168+.7094507040001145*im;
         .8490194865023282-.5100588510854889*im;
         .8490194865023282+.5100588510854889*im;
         .9602996918859324-.2681678908797330*im;
         .9602996918859324+.2681678908797330*im;
         1.999999940395313];
    z = [z 2*ones(n)]
    p, PolyZeros(z)
end


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using DifferentialEquations
using Plots
using GraphPlot

include("utils\\randutils.jl")
include("SI\\contact_process.jl")
include("SI\\metapopulation.jl")

# Test that contact process on a complete graph is equivalent to single node metapopulation
N = 100
g = complete_graph(N)
β = 0.3

frac_infected = 0.05
n_infected =  Int(N*frac_infected)

X = falses(N)
X[1:n_infected] .= true # complete graph. choice of nodes makes no difference

t, Xs = contact_process_gillespie(g, X, β)

plot(t, Xs/N,
    xlabel="time",
    ylabel="Fraction of infected individuals",
    label="",
    line=:steppre,
    title="Contact Process"
    )

t_mc, X_mc, sd_mc = contact_process_montecarlo(g, X, β, tmax = 0.3, nbins=400, nsims = 1000)

plot(t_mc, X_mc/N,
    label="Average (Contact Process)",
    legend=:bottomright,
    xlabel = "time",
    ylabel = "Fraction of infected individuals",
    )

mpp = Metapopulation(SimpleGraph(1), [1], β, 0.0)
s0 = [N-n_infected]
i0 = [n_infected]

t_mpp, s_mpp, i_mpp = metapopulation_gillespie(mpp,s0,i0, nmax = 100, tmax=0.3)

plot(t_mpp, sum(i_mpp, dims=2)/N,
    label="",
    line=:steppre,
    xlabel="Time",
    ylabel="Fraction of infected individuals",
    title="Metapopulation"
    )

t_mpp_mc, s_mpp_mc, i_mpp_mc = metapopulation_montecarlo(mpp,s0,i0, nmax=100, tmax=0.3, nbins=400)

plot(t_mc, X_mc/N,
    label="Average (Contact Process)",
    legend=:bottomright,
    xlabel = "time",
    ylabel = "Fraction of infected individuals",
    )
plot!(t_mpp_mc, i_mpp_mc/N,
    label = "Average (Metapopulation)",
    xlabel="Time",
    ylabel="Fraction of infected individuals"
    )

plot(t_mc, log10.(abs.((X_mc-i_mpp_mc)/N)),
    xlabel = "time",
    ylabel="log10 Deviation")

# Test that reusing the same seed for the rng yields the same result

N = 100
g = erdos_renyi(N, 0.5)
β = 0.3

frac_infected = 0.05
n_infected =  Int(N*frac_infected)

X = falses(N)
X[1:n_infected] .= true # complete graph. choice of nodes makes no difference

Random.seed!(42) # magic number
t1, Xs1, Es1 = contact_process_gillespie(g, X, β, record=true)

Random.seed!(42) # magic number
ts2, Xs2, Es2 = contact_process_gillespie(g, X, β, record=true)

ts1 == ts2
Xs1 == Xs2
Es1 == Es2

# Test the influence of network structure on epidemic

# other tests
N = 1000
g = complete_graph(N)
g = random_configuration_model(N,fill(3,N))
β = 0.3
X = falses(N)
X[rand_combination(N,50)] .= true
tmax=0.05
@time ts, Xs = contact_process_gillespie(g, X, β)
@time ts, Xs, Es = contact_process_gillespie(g, X, β, record=true)

plot(ts, Xs/N, xlabel="time", ylabel="Fraction of infected individuals", label="")

@time t_mc, X_mc, sd_mc = contact_process_montecarlo(g, X, β, tmax = tmax, nbins=400, nsims = 1000)

t_mf, X_mf, sol = contact_process_ode(g, X, β, tmax=0.05, saveat=t_mc)

plot(t_mf, X_mf, label="")

plot(t_mf, sum(X_mf, dims=2)/N, label="Mean Field")

plot(t_mc, X_mc/N,
    label="Average",
    legend=:bottomright,
    xlabel = "time",
    ylabel = "Fraction of infected individuals"
    )
plot!(t_mc, (X_mc+sd_mc)/N, linestyle=:dash, color=:lightblue, label="Standard Deviation")
plot!(t_mc, (X_mc-sd_mc)/N, linestyle=:dash, color=:lightblue, label="")
plot!(t_mf, sum(X_mf,dims=2)/N, label="Individual Based Mean Field")

using Statistics
using ProgressMeter

function compare_methods_time(n,m,step)
    N = 10:step:n
    m_a = zeros(length(N))
    m_s = zeros(length(N))
    m_t = zeros(length(N))
    t_a = zeros(length(N))
    t_s = zeros(length(N))
    t_t = zeros(length(N))
    qt_a = zeros(length(N),4)
    qt_s = zeros(length(N),4)
    qt_t = zeros(length(N),4)
    ta = zeros(m)
    ts = zeros(m)
    tt = zeros(m)
    β = 0.3
    f = 0.05 # at least 5 percent infected
    p = Progress(length(N)*m, dt=1.0)
    for i in 1:length(N)
        g = random_configuration_model(N[i], fill(4,N[i]))
        X = falses(N[i])
        n_inf = N[i] < 20 ? 1 : Int(floor(f*N[i]))
        X[1:n_inf] .= true
        for k in 1:m
            ta[k] = @elapsed contact_process_gillespie(g, X, β, nmax=N[i], method=:array)
            ts[k] = @elapsed contact_process_gillespie(g, X, β, nmax=N[i], method=:sparse)
            tt[k] = @elapsed contact_process_gillespie(g, X, β, nmax=N[i], method=:tree)
            ProgressMeter.next!(p; showvalues= [(:N,N[i]),(:iter,k)])
        end
        m_a[i] = mean(ta)
        m_s[i] = mean(ts)
        m_t[i] = mean(tt)
        t_a[i] = median(ta)
        t_s[i] = median(ts)
        t_t[i] = median(tt)
        qt_a[i,:] = quantile(ta, [0.05,0.25,0.75,0.95])
        qt_s[i,:] = quantile(ts, [0.05,0.25,0.75,0.95])
        qt_t[i,:] = quantile(tt, [0.05,0.25,0.75,0.95])
    end
    return N, m_a, m_s, m_t, t_a, t_s, t_t, qt_a, qt_s, qt_t
end

n = 3000
m = 1000
step = 500

N, ma, ms, mt, ta, ts, tt, qta, qts, qtt = compare_methods_time(n,m,step)

plot(N, ta, label="Median Array", legend=:topleft, linecolor=:blue)
plot!(N, ts, label="Median Sparse", linecolor=:red)
plot!(N, tt, label="Median Tree", linecolor=:green)
plot!(N, ma, label="Mean Array", linecolor=:blue, linestyle=:dash)
plot!(N, ms, label="Mean Sparse", linecolor=:red, linestyle=:dash)
plot!(N, mt, label="Mean Tree", linecolor=:green, linestyle=:dash)
plot!(N, qta, label="",
    linestyle=:dot, linecolor=:blue)
plot!(N, qts, label="",
    linestyle=:dot, linecolor=:red)
plot!(N, qtt, label="",
    linestyle=:dot, linecolor=:green)

scatter(N.^2, ta, label="", color=:blue)
scatter!(N.^2, ts, label="", color=:red)
scatter(N.*log2.(N), tt,
    label="",
    color=:green,
    xlabel="N log(N)",
    ylabel="t(N)")

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
function Damper(; name, b = 1.0)
    # 2 pin component
    @named oneport = OnePort()
    @unpack u, f = oneport
    ps = @parameters b=b
    eqs = [
        u ~ f / b
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function Mass(; name, m = 1.0)
    @named oneport = OnePort()
    @unpack u, f = oneport
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
        D(u) ~ f / m
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function Spring(; name, k = 1.0)
    @named oneport = OnePort()
    @unpack u, f = oneport
    ps = @parameters k=k
    D = Differential(t)
    eqs = [
        D(f) ~ k * u
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any objective functions supported by the
bridge `BT`. This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref)
which only bridges the objective functions that are unsupported by the internal model,
even if they are supported by one of its bridges.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: MOIB.AbstractBridgeOptimizer
    model::OT
    map::Map # `MOI.ObjectiveFunction` -> objective bridge
end
function SingleBridgeOptimizer{BT}(model::OT) where {BT, OT <: MOI.ModelLike}
    SingleBridgeOptimizer{BT, OT}(model, Map())
end

function bridges(bridge::MOI.Bridges.AbstractBridgeOptimizer)
    return EmptyMap()
end
bridges(bridge::SingleBridgeOptimizer) = bridge.map

MOIB.supports_constraint_bridges(::SingleBridgeOptimizer) = false
function MOIB.is_bridged(::SingleBridgeOptimizer, ::Type{<:MOI.AbstractSet})
    return false
end
function MOIB.is_bridged(::SingleBridgeOptimizer,
                         ::Type{<:MOI.AbstractFunction},
                         ::Type{<:MOI.AbstractSet})
    return false
end
function MOIB.supports_bridging_objective_function(
    ::SingleBridgeOptimizer{BT}, F::Type{<:MOI.AbstractScalarFunction}) where BT
    return supports_objective_function(BT, F)
end
function MOIB.is_bridged(
    b::SingleBridgeOptimizer, F::Type{<:MOI.AbstractScalarFunction})
    return MOIB.supports_bridging_objective_function(b, F)
end
function MOIB.bridge_type(::SingleBridgeOptimizer{BT},
                          ::Type{<:MOI.AbstractScalarFunction}) where BT
    return BT
end


text:
export result_nodal_reaction,result_nodal_displacement,result_nodal_time_history

"""
   result_nodal_displacement(assembly,lc_id,node_id)
获取已求解的assembly实例的单个节点位移
# Arguments
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `node_id`: 节点id
"""
function result_nodal_displacement(assembly,lc_id,node_id;substep=0)
    path=assembly.working_path
    node_id=string(node_id)
    if !(node_id in keys(assembly.structure.nodes))
        throw("node with id "*string(node_id)*" doesn't exist!")
    end
    lc_id=string(lc_id)
    if !(lc_id in get_all_id(assembly.lcset))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    hid=assembly.structure.nodes[node_id].hid
    DOFs=6*hid-5:6*hid
    if lc_id in keys(assembly.lcset.statics)
        u=read_vector(path*"/.analysis",lc_id*"_u.v")
        return u[DOFs]
    else
        u=read_matrix(path*"/.analysis",lc_id*"_u.m")
        if substep==0
            return u[DOFs,:]
        elseif substep<=size(u,2)
            return u[DOFs,substep]
        else
            throw("substep not exists!")
        end
    end
end

"""
   result_nodal_displacement(assembly,lc_id)
获取已求解的assembly实例的所有节点位移
# 参数
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `node_id`: 节点id
# 返回
- id=>位移字典
"""
function result_nodal_displacement(assembly,lc_id)
    path=assembly.working_path
    lc_id=string(lc_id)
    if !(lc_id in get_all_id(assembly.lcset))
        throw("loadcase with id "*string(lc_id)*" doesn't exist!")
    end
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    res=Dict{String,Array{Float64}}
    for node in values(assembly.structure.nodes)
        hid=node.hid
        DOFs=6*hid-5:6*hid
        res[node.id]=u[DOFs]
    end
    return res
end

"""
   result_nodal_reaction(assembly,lc_id,node_id)
获取已求解的assembly实例的单个节点反力
# Arguments
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `node_id`: 节点id
"""
function result_nodal_reaction(assembly,lc_id,node_id)
    path=assembly.working_path
    node_id=string(node_id)
    if !(node_id in keys(assembly.structure.nodes))
        throw("node with id "*string(node_id)*" doesn't exist!")
    end
    lc_id=string(lc_id)
    if !(lc_id in get_all_id(assembly.lcset))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    F=read_vector(path*"/.analysis",lc_id*"_F.v")
    F₀=zero(F)
    hid=assembly.structure.nodes[node_id].hid
    DOFs=6*hid-5:6*hid
    R=F-recursive_P(assembly,lc_id)
    return R[DOFs]
end

function recursive_P(assembly,lc_id)
    lc_id=string(lc_id)
    lc=assembly.lcset.statics[lc_id]
    if lc.plc==""
        return lc.P
    else
        return lc.P+recursive_P(assembly,lc.plc)
    end
end

"""
   result_nodal_time_history(assembly,lc_id,node_id)
获取已求解的assembly实例的单个节点响应时程
# Arguments
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `node_id`: 节点id
- `value`:值，0-位移；1-速度；2-加速度
- `dir`:方向，1-1方向，2-2方向，3-3方向，4-绕1轴转动方向，5-绕2轴转动方向，6-绕3轴转动方向
"""
function result_nodal_time_history(assembly,lc_id,node_id,value=0,dir=1)
    path=assembly.working_path
    node_id=string(node_id)
    if !(node_id in keys(assembly.structure.nodes))
        throw("node with id "*string(node_id)*" doesn't exist!")
    end
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.time_histories))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    hid=assembly.structure.nodes[node_id].hid
    DOFs=6*hid-5:6*hid
    if value==0
        u=read_matrix(path*"/.analysis",lc_id*"_u.m")
    elseif value==1
        u=read_matrix(path*"/.analysis",lc_id*"_v.m")
    elseif value==2
        u=read_matrix(path*"/.analysis",lc_id*"_a.m")
    end
    u[6*hid-6+dir,:]
end


text:
@testset "test_complex_systems" begin

C_1 = zpk([], [-1+im], 1.0+1im)
C_2 = zpk([-1+im], [], 1.0+1im)

# Basic arittmetic
@test C_1 + C_1 == zpk([], [-1+im], 2+2im)
@test (1.0+im)*C_1 == zpk([], [-1+im], 2.0im)
@test C_2 + C_2 == zpk([-1+im], [], 2+2im)

@test C_1 * C_1 == zpk([], [-1+im,-1+im], 2.0im)
@test C_2 * C_2 ≈ zpk([-1+im,-1+im], [], 2.0im)

@test im*ss(1) == ss(im)

@test poles(zpk([], [-1+im,-1+im,0], 2.0im)) == [-1+im,-1+im,0]
@test tzeros(zpk([-1+im,-1+im,0], [-2], 2.0im)) == [-1+im,-1+im,0]

@test zpk( tf([1.0, 1+im], [1.0, 2+im]) ) == zpk( [-1-im], [-2-im], 1.0+0im)

@test minreal(zpk([-1+im], [-1+im,-1+im],1+0im)) == zpk([], [-1+im],1+0im)
@test minreal(zpk([-1+im, -1+im], [-1+im],1+1im)) == zpk([-1+im], [], 1+1im)

@test_throws ArgumentError zpk([-1+im], [-1+im,-1+im],1) #  Given the type of k this should be a real-coefficient system, but poles and zeros don't come in conjugate pairs

@test zpk([-2+im], [-1+im],1+0im)*zpk([], [-1+im],1+0im) == zpk([-2+im], [-1+im, -1+im], 1+0im)
@test zpk([], [-2], 2) + zpk([], [-1], 1) == zpk([-4/3], [-2,-1], 3)

@test tf(zpk([-2+im], [-1+im],1+0im)) == tf([1, 2-im], [1, 1-im])

@test 1 / ( tf("s") + 1 + im ) == tf([1], [1, 1+im])

s = tf("s");
@test tzeros(ss(-1, 1, 1, 1.0im)) ≈ [-1.0 + im] rtol=1e-15
@test tzeros(ss([-1.0-im 1-im; 2 0], [2; 0], [-1+1im -0.5-1.25im], 1)) ≈ [-1-2im, 2-im]

@test tzeros(ss((s-2.0-1.5im)^3/(s+1+im)/(s+2)^3)) ≈ fill(2.0 + 1.5im, 3) rtol=1e-4
@test tzeros(ss((s-2.0-1.5im)*(s-3.0)/(s+1+im)/(s+2)^2)) ≈ [3.0, 2.0 + 1.5im] rtol=1e-14

end


text:
# To build and compile Theia
#
# Create a cmake/FindCeres.cmake with the following lines
#
# SET(CERES_FOUND TRUE)
# SET(CERES_LIBRARY "/usr/local/lib/libceres.so")
# SET(CERES_LIBRARIES "/usr/local/lib/libceres.so")
# SET(CERES_INCLUDE_DIRS "/usr/local/include/ceres/")
#
# Create cmake/FindRapidJSON.cmake
#
# SET(RapidJSON_FOUND TRUE)
# SET(RAPIDJSON_INCLUDE_DIRS "/usr/include/rapidjson")
#
# Run
#
# cmake -DSTATX_WITH_CERES=ON -DBUILD_SHARED_LIBS=ON -DBUILD_TESTING=OFF -DCMAKE_CXX_FLAGS="-shared -fPIC" ..
#
# Edit
#
# applications/CMakeFiles/view_reconstruction.dir/link.txt
#
# and remove the -lFALSE flags. This is the shared library version,
# therefore the Ceres library should also be compiled as a shared
# library, but it is easier to compile.
#
# There are other minor issues that can be easily solved.

using Cxx
using Libdl

push!(DL_LOAD_PATH, "/usr/local/lib")

function ransac_example()

    for header in ["/usr/include/eigen3/",
                   "/usr/local/include/theia/libraries/vlfeat",
                   "/usr/include/suitesparse",
                   "/usr/local/include/theia/libraries/",
                   "/usr/local/include/theia/libraries/statx",
                   "/usr/local/include/theia/libraries/optimo",
                   "/usr/local/include",
                   "/usr/local/include/theia",
                   "/usr/include"]

        addHeaderDir(header, kind=C_User)

    end

    Libdl.dlopen("libvlfeat.so", Libdl.RTLD_GLOBAL)
    Libdl.dlopen("libflann_cpp.so", Libdl.RTLD_GLOBAL)
    Libdl.dlopen("libakaze.so", Libdl.RTLD_GLOBAL)
    Libdl.dlopen("libceres.so", Libdl.RTLD_GLOBAL)
    Libdl.dlopen("libtheia.so", Libdl.RTLD_GLOBAL)

    cxx"""
#include<theia.h>

// Our "data".
struct Point {
  double x; double y;
};

// Our "model".
struct Line {
  double m; double b;
};

// Estimator class.
class LineEstimator: public theia::Estimator<Point, Line> {

public:

   LineEstimator(){};

  ~LineEstimator(){};

  // Number of points needed to estimate a line.
  double SampleSize() const override { return 2; }

  // Estimate a line from two points.
  bool EstimateModel(const std::vector<Point>& data,
                     std::vector<Line>* models) const override {
    Line model;
    model.m = (data[1].y - data[0].y)/(data[1].x - data[0].x);
    model.b = data[1].y - model.m*data[1].x;
    models->push_back(model);
    return true;
  }

  // Calculate the error as the y distance of the point to the line.
  double Error(const Point& point, const Line& line) const override {
    return point.y - (line.m*point.x + line.b);
  }

};

void run_ransac(void) {
  // Generate your input data using your desired method.
  // We put pseudo-code here for simplicity.
  std::vector<Point> input_data;

  // Add 700 inliers.
  for (int i = 0; i < 700; i++) {
    double x = i * 20.0 / 700.0;
    Point inlier_point = { x , 3.0 * x + 4.0 };
    input_data.push_back(inlier_point);
  }
  // Add 300 outliers.
  for (int i = 0; i < 300; i++) {
    double x = i * 20.0 / 300.0;
    Point outlier_point = { x , 4.0 };
    input_data.push_back(outlier_point);
  }

  // Specify RANSAC parameters.
  double error_threshold = 0.3;
  int min_num_inliers = 600;
  int max_iters = 1000;

  // Estimate the line with RANSAC.
  LineEstimator line_estimator;
  Line best_line;
  // Set the ransac parameters.
  theia::RansacParameters params;
  params.error_thresh = 0.1;

  // Create Ransac object, specifying the number of points to sample to
  // generate a model estimation.
  theia::Ransac<LineEstimator> ransac_estimator(params, line_estimator);
  // Initialize must always be called!
  ransac_estimator.Initialize();

  theia::RansacSummary summary;
  ransac_estimator.Estimate(input_data, &best_line, &summary);
  LOG(INFO) << "Line m = " << best_line.m << "*x + " << best_line.b;

}
    """

end

"""

Test the creating and call of a C++ function in Julia. Should print
10.

"""
function hello_world()

    cxx""" #include<iostream> """

    cxx"""
           void printme(int x) {
              std::cout << x << std::endl;
           }
    """

    @cxx printme(10)

end


text:
@testset "Memory: Static" begin

function memory_static_kernel(a,b)
    # Local
    ptr_local = alloc_local(:local, Float32, 1)
    unsafe_store!(ptr_local, a[1])
    b[1] = unsafe_load(ptr_local)

    # Region
    #= TODO: AMDGPU target cannot select
    ptr_region = alloc_special(Val(:region), Float32, Val(AS.Region), Val(1))
    unsafe_store!(ptr_region, a[2])
    b[2] = unsafe_load(ptr_region)
    =#

    # Private
    #= TODO
    ptr_private = alloc_special(Val(:private), Float32, Val(AS.Private), Val(1))
    unsafe_store!(ptr_private, a[3])
    b[3] = unsafe_load(ptr_private)
    =#

    nothing
end

A = ones(Float32, 1)
B = zeros(Float32, 1)

RA = ROCArray(A)
RB = ROCArray(B)

wait(@roc memory_static_kernel(RA, RB))

@test Array(RA) ≈ Array(RB)

end

@testset "Memory: Dynamic" begin

function malloc_kernel(X)
    ptr = AMDGPU.malloc(Csize_t(4))
    X[1] = reinterpret(UInt64, ptr)
    AMDGPU.free(ptr)
    nothing
end

RA = ROCArray(zeros(UInt64, 1))

wait(@roc malloc_kernel(RA))

@test Array(RA)[1] != 0

end

@testset "Memcpy/Memset" begin

function memcpy_kernel(X,Y)
    AMDGPU.memcpy!(Y.ptr, X.ptr, sizeof(Float32)*length(X))
    nothing
end

A = rand(Float32, 4)
B = zeros(Float32, 4)
RA, RB = ROCArray.((A,B))

wait(@roc memcpy_kernel(RA,RB))

@test A == collect(RA) == collect(RB)

function memset_kernel(X,y)
    AMDGPU.memset!(X.ptr, y, div(length(X),2))
    nothing
end

A = zeros(UInt8, 4)
RA = ROCArray(A)

wait(@roc memset_kernel(RA,0x3))

@test all(collect(RA)[1:2] .== 0x3)
@test all(collect(RA)[3:4] .== 0x0)

end


text:
applyblock(f, x::Union{Block,Loop}) = typeof(x)(f(x.body))
applyblock(f, x::If) = If(f(x.t), f(x.f))
applyblock(f, x) = x

walk(x, inner, outer) = outer(applyblock(xs -> map(inner, xs), x))

postwalk(f, x) = walk(x, x -> postwalk(f, x), f)
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)

# remove nops

nops(x) = prewalk(x -> applyblock(is -> filter(i -> i ≠ nop, is), x), x)

# dead code

function deadcode(x)
  prewalk(x) do x
    applyblock(x) do is
      i = findfirst(b -> b isa Branch && !b.cond, is)
      i == 0 ? is : is[1:i]
    end
  end
end

mapbranches(f, b, level = -1) = applyblock(is -> mapbranches.(f, is, level+1), b)
mapbranches(f, b::Branch, level = -1) = f(b, level)

function branches_to(b)
  result = false
  mapbranches(b) do b, l
    b.level==l && (result = true)
    b
  end
  return result
end

# After removing a block that has no branches to it, decrement all of the
# branches that point to any level above it.
function decbranches(i)
  mapbranches(i) do b, l
    b.level > l ? Branch(b.cond, b.level-1) : b
  end
end

# After adding a block increment all branches pointing to a level above it.
function incbranches(i)
  mapbranches(i) do b, l
    b.level >= l ? Branch(b.cond, b.level+1) : b
  end
end

# recover ifs

function makeifs(code)
  prewalk(code) do x
    x isa Block || return x
    i = findfirst(x -> x isa Branch, x.body)
    i != 0 && x.body[i].cond || return x
    cond = x.body[1:i-1]
    cond[end] == Op(i32, :eqz) ? pop!(cond) : push!(cond, Op(i32, :eqz))
    fixbranches = incbranches(Block(x.body[i+1:end])).body
    return Block([cond..., If(fixbranches, [])])
  end
end

# remove unused blocks
# TODO: collapse blocks with the same head

function rmblocks(code)
  prewalk(code) do x
    applyblock(x) do is
      is′ = []
      for i in is
        i isa Block && !branches_to(i) ?
          append!(is′, decbranches(i).body) :
          push!(is′, i)
      end
      is′
    end
  end
end

optimise(b) = b |> deadcode |> makeifs |> rmblocks


text:
#!/usr/local/bin/julia -q --color=yes

dataFile = Pkg.dir(homedir(), ".jtodo.jls")
pending = '✘'
isdone = '✔'

save(data) = begin
  open(f->serialize(f, data), dataFile, "w")
end

load() = begin
  open(deserialize, dataFile)
end

if !isfile(dataFile)
  save([])
end

todos = load()

ls(args) = begin
  for (i, todo) in enumerate(todos)
    item = "$i $todo \n"
    len = length(args)

    if search(todo, pending) > 0
      if len == 0 || args[1] != "-d"
        print_with_color(:red, item)
      end
    else
      if len > 0 && args[1] != "-p"
        print_with_color(:green, item)
      end
    end
  end
end

final() = begin
  ls([])
  save(todos)
end

intify(id) = begin
  parse(Int, id)
end

ad(args) = begin
  for i = 1:length(args)
    newTodo = "| $pending $(args[i])"
    push!(todos, newTodo)
  end

  final()
end

rm(args) = begin
  sort!(args, rev=true)

  for arg in args
    splice!(todos, parse(Int, arg))
  end

  final()
end

ch(args) = begin
  (id, cont) = args
  intid = intify(id)

  if search(todos[intid], pending) > 0
    todos[intid] = "| $pending $cont"
  else
    todos[intid] = "| $isdone $cont"
  end

  final()
end

sw(args) = begin
  (id1, id2) = map(intify, args)
  (todos[id1], todos[id2]) = [todos[id2], todos[id1]]

  final()
end

tg(args) = begin
  id = intify(args[1])

  if search(todos[id], pending) > 0
    todos[id] = replace(todos[id], pending, isdone)
  else
    todos[id] = replace(todos[id], isdone, pending)
  end

  final()
end

cl(args) = begin
  len = length(todos)

  while len > 0
    if(search(todos[len], isdone) > 0)
      splice!(todos, len)
    end
    len = len - 1
  end

  final()
end

usage(args) = begin
  str = """
    ls => list
    ad => add
    rm => remove
    ch => change
    sw => swap
    tg => toggle
    cl => clean
  """

  println(str)
end

handler(args) = begin
  cmd = args[1]
  newArgs = args[2:end]
  dict = Dict(
    "ls" => ls,
    "ad" => ad,
    "rm" => rm,
    "ch" => ch,
    "sw" => sw,
    "tg" => tg,
    "cl" => cl,
  )

  get(dict, cmd, usage)(newArgs)
end

main() = begin
  if length(ARGS) == 0
    usage()
  else
    handler(ARGS)
  end
end

main()

max_stars_repo_name: moldray/juliaTodo

max_issues_repo_name: moldray/juliaTodo

max_forks_repo_name: moldray/juliaTodo


text:
using Test
using ScottishTaxBenefitModel
using .Utils:coarse_match

using DataFrames, CSV
using .ExampleHelpers

#
# Donor and Recipient each have 2 fields `a` and `b` filled with random
# integers. 2 levels: _1 with exact numbers and 2 with 0/1 coarsened versions
# so matches can be, since both have to match & `a` is coarsened before `b`
# quality  a_1 a_2 b_1 b_2
#    1      X   -   X   -
#    2      -   X   X   -
#    3      -   -   X   X

printrows= false

n = 10_000
donor = DataFrame( sernum=collect(1:n), a_1=rand(1:50,n), b_1=rand(100:1500,n))
# coarsend
donor.a_2 = (donor.a_1 .<= 25)
donor.b_2 = (donor.b_1 .<= 600)

m = 5_000
recip = DataFrame( sernum=collect(1:m), a_1=rand(2:50,m), b_1=rand(110:1700,m))
# coarsend
recip.a_2 = (recip.a_1 .<= 25)
recip.b_2 = (recip.b_1 .<= 600)
max_matches = 2500

@testset "simple matching case" begin
    i = 0
    for r1 in eachrow(recip)
        # r1 = recip[1,:]
        i += 1
        if i > 250
            break
        end
        if printrows
            println("r1=$r1")
        end
        matches = coarse_match(
            r1,
            donor,
            [:a, :b],
            2 )
        # println(matches)
        nmatches = sum( matches.matches )
        println( "$nmatches ")
        donor.quality = matches.quality
        matchedrows = donor[matches.matches,:]
        # @test sum( matches.matches ) >= max_matches
        n = 0
        for match in eachrow( matchedrows )
            if printrows
                println( "on match=$(match)" )
            end
            # @test match.quality in 1:8
            if match.quality == 2 # 1+1
                @test (match.a_1 == r1.a_1)&&(match.b_1 == r1.b_1) # `a` and `b` should match fine.
            elseif match.quality == 5 # 1 + 2^2
                @test ! ((match.a_1 == r1.a_1)&&(match.b_1 == r1.b_1)) # can't be a q=1 match
                @test ((match.a_2 == r1.a_2)&&(match.b_1 == r1.b_1)||
                       (match.a_1 == r1.a_1)&&(match.b_2 == r1.b_2))# should be a_2 matches coarse, b_1 matches fine
            elseif match.quality == 8 # 2^2 + 2^2
                if printrows
                    println( "q=3; r1=$r1" )
                    println( "q=3; match=$match" )
                end
                @test ! ((match.a_1 == r1.a_1)&&(match.b_1 == r1.b_1)) # can't be a q=1 match
                @test ! ((match.a_2 == r1.a_2)&&(match.b_1 == r1.b_1)) # can't be a q=2 match
                @test (match.a_2 == r1.a_2)&&(match.b_2 == r1.b_2) # should be `a` and `b` match coarse

            end

        end
    end # each recipient
end # testset

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
"""
    Planar quadrotor
"""

mutable struct Quadrotor2D{T}
    n::Int
    m::Int
    d::Int

    L::T    # length
    mass::T # mass
    J::T    # inertia
    g::T    # gravity
end

function f(model::Quadrotor2D, x, u, w)
    qdd1 = -1.0 * sin(x[3]) / model.mass * (u[1] + u[2])
    qdd2 = -1.0 * model.g + cos(x[3]) / model.mass * (u[1] + u[2])
    qdd3 = model.L / model.J * (-1.0 * u[1] + u[2])

    @SVector [x[4],
              x[5],
              x[6],
              qdd1,
              qdd2,
              qdd3]
end

n, m, d = 6, 2, 0
model = Quadrotor2D(n, m, d, 1.0, 1.0, 1.0, 9.81)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Test
using FiniteDiff1D
using Symbolics
using Latexify
using LaTeXStrings

function post_process(sol, fʲ_i, k_range, n_unknowns)
    sol = fʲ_i ~ sol
    s = latexify(sol)

    s = LaTeXStrings.LaTeXString(s)
    for (i, k) in enumerate(k_range)
        str_to = if k==0
            "_{i}"
        elseif k > 0
            "_{i+$k}"
        elseif k < 0
            "_{i$k}"
        end

        str_from = if k==0 || k > 0
            "_{k{_$k}}"
        elseif k < 0
            "_{k\\_-{_$(abs(k))}}"
        end

        # Fixup k-subscripts
        s = replace(s, str_from => str_to)
    end

    for j in 1:n_unknowns
        # Fixup j-subscripts / unicode superscripts
        s = replace(s, "f\\^j_{i{_$j}}" => "f$(repeat("'", j))_{i}")
    end

    s = LaTeXStrings.LaTeXString(s)
    return s
end

@testset "Test" begin
    n_left, n_right = -1, 1
    solns, f_k, fʲ_i, k_range, n_unknowns = FiniteDiff1D.generate_stencil(n_left, n_right)
    s = post_process.(solns, fʲ_i, Ref(k_range), Ref(n_unknowns))
    @show s
end
nothing


text:
"""
    function pcdos(file; threshold = 0.02, gnz = 0.)
---
`pcdos` means PCD method for one-step.

This function:
- construct a `G` matrix with the genotypes in `file`.
- does a pivoted Cholesky decomposition, say, `C`.
- take individuals whose corresponding `C.U` is greater than `threshold`.
- calculate `G` again and return ID used and `G` inverse.
"""
function pcdos(file; threshold = 0.02, gnz = 0.)
    @info "Reading genotypes"
    dat = read_g012(file)
    @info "Calculating GRM with VanRaden method I"
    G = grm(dat.genotype)
    @info "Doing pivoted Cholesky decomposition"
    C = cholesky(G, Val(true), check=false)
    @info "Find individuals whose correspoding diagonals in `C.U > $threshold`"
    i = 0
    for t in diag(C.U)
        t < threshold && break
        i += 1
    end
    id = sort(C.p[1:i])
    @info "Calculating reduced G and its inverse"
    G = grm(dat.genotype[:, id])
    return (id = dat.ID[id], gi = inv(G))
end


text:
a=[]
# [try, catch](http://julia.readthedocs.org/en/latest/manual/control-flow/#the-try-catch-statement) can be used to deal with errors as with many other languages
try
    push!(a,1)
catch err
    showerror(STDOUT, err, backtrace());println()
end
println("Continuing after error")

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
using Distributed
addprocs(4)
@everywhere using Distributions, Turing, Random, Parameters, MCMCChains
ProjDir = @__DIR__
cd(ProjDir)
path = pwd()
@everywhere include($path*"/LNR.jl")

@everywhere @model model(data, Nr) = begin
    μ = Array{Real,1}(undef,Nr)
    μ ~ [Normal(0,3)]
    s ~ Uniform(0, pi/2)
    σ = tan(s)
    for i in 1:length(data)
        data[i] ~ LNR(μ, σ, 0.0)
    end
end

Random.seed!(343)
Nreps = 50
Nr = 3
Nobs = 10
μ = rand(Normal(0, 3), Nr)
σ = rand(Uniform(.2, 1))
dist = LNR(μ=μ, σ=σ, ϕ=0.0)
data = rand(dist, Nobs)

function estimate1(model,specs,Nsamples;Nchains=4)
    return reduce(chainscat, pmap(x->sample(model, specs, Nsamples), 1:Nchains))
end

pmap(x->sample(model(data, Nr), NUTS(1000, .8), 2000), 1:4)

reduce(chainscat, pmap(x->sample(model(data, Nr),NUTS(1000, .8), 2000), 1:4))

estimate1(model, NUTS(1000, .8), 2000)

max_stars_repo_name: itsdfish/Miscellaneous

max_issues_repo_name: itsdfish/Miscellaneous

max_forks_repo_name: itsdfish/Miscellaneous


text:
using QuantumControlBase: getcontrols, discretize
using Printf
using Dates

"""Result object returned by [`optimize_krotov`](@ref)."""
mutable struct KrotovResult{STST}
    tlist::Vector{Float64}
    iter_start::Int64  # the starting iteration number
    iter_stop::Int64 # the maximum iteration number
    iter::Int64  # the current iteration number
    secs::Float64  # seconds that the last iteration took
    tau_vals::Vector{ComplexF64}
    J_T::Float64  # the current value of the final-time functional J_T
    J_T_prev::Float64  # previous value of J_T
    guess_controls::Vector{Vector{Float64}}
    optimized_controls::Vector{Vector{Float64}}
    states::Vector{STST}
    start_local_time::DateTime
    end_local_time::DateTime
    records::Vector{Tuple}  # storage for info_hook to write data into at each iteration
    converged::Bool
    message::String

    function KrotovResult(problem)
        tlist = problem.tlist
        controls = getcontrols(problem.objectives)
        iter_start = get(problem.kwargs, :iter_start, 0)
        iter_stop = get(problem.kwargs, :iter_stop, 5000)
        iter = iter_start
        secs = 0
        tau_vals = Vector{ComplexF64}()
        guess_controls = [discretize(control, tlist) for control in controls]
        J_T = 0.0
        J_T_prev = 0.0
        optimized_controls = [copy(guess) for guess in guess_controls]
        states = [similar(obj.initial_state) for obj in problem.objectives]
        start_local_time = now()
        end_local_time = now()
        records = Vector{Tuple}()
        converged = false
        message = "in progress"
        new{eltype(states)}(
            tlist,
            iter_start,
            iter_stop,
            iter,
            secs,
            tau_vals,
            J_T,
            J_T_prev,
            guess_controls,
            optimized_controls,
            states,
            start_local_time,
            end_local_time,
            records,
            converged,
            message
        )
    end
end

Base.show(io::IO, r::KrotovResult) = print(io, "KrotovResult<$(r.message)>")
Base.show(io::IO, ::MIME"text/plain", r::KrotovResult) = print(
    io,
    """
Krotov Optimization Result
--------------------------
- Started at $(r.start_local_time)
- Number of objectives: $(length(r.states))
- Number of iterations: $(max(r.iter - r.iter_start, 0))
- Value of functional: $(@sprintf("%.5e", r.J_T))
- Reason for termination: $(r.message)
- Ended at $(r.end_local_time) ($(Dates.canonicalize(Dates.CompoundPeriod(r.end_local_time - r.start_local_time))))
"""
)


text:
using RunTests
using Base.Test

@testmodule FirstSetOfTests begin

  x = 100

  function test_apple()
      @test true
  end

  function test_aardvark()
      @test true
  end

  function test_andrew()
    @test true
  end

  function foo()
    @test x==100
  end

  test_hello() = foo()

  function test_world()
    println("seen")
    @test true
    println("also seen")
    @test false
    println("never seen")
    @test 100==100
  end

  function test_baz()
    @test false
  end

  @skipif false function test_not_skipped()
    @test true
  end

  @skipif true function test_skipped()
    @test true
  end

  const SKIP_FALSE = false

  @skipif SKIP_FALSE function test_not_skipped_constant()
    @test true
  end

  const SKIP_TRUE = true

  @skipif SKIP_TRUE function test_skipped_constant()
    @test true
  end

  @skipif (x->x)(100+25/5)==102+25/5 function test_not_skipped_complex()
    @test true
  end

  @skipif (x->x)(100+25/5)==100+25/5 function test_skipped_complex()
    @test true
  end

  @xfail function test_xfails()
    @test false
  end

  @xfail function test_xpasses()
    @test true
  end

  @parameterize 1:4 function test_parameterised(x)
    @test x<3
  end

  @parameterize [1,2,3,4] function test_parameterised_with_types(x::Integer)
    @test x<3
  end

  @parameterize [(1,4),(2,3),(3,2),(4,1)] function test_parameterised_2_args(x::Integer, y::Any)
    @test x<y
  end

  @parameterize zip(1:4,4:-1:1) function test_parameterised_zip(x::Integer, y::Any)
    @test x<y
  end

  const CONST_PARAMS = [1,2,3,4]

  @parameterize CONST_PARAMS function test_parameterised_constant(x)
    @test x<3
  end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libyaml"
version = v"0.2.1"

# Collection of sources required to build libyaml
sources = [
    "https://github.com/yaml/libyaml/archive/0.2.1.tar.gz" =>
    "1d2aeb87f7d317f1496e4c39410d913840714874a354970300f375eec9303dc4",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libyaml-0.2.1
./bootstrap
./configure --prefix=$prefix --host=$target
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms_base = [
    FreeBSD(:x86_64),
    Linux(:aarch64, libc=:glibc),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    Linux(:i686, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:x86_64, libc=:musl),
]
platforms_osx = [
    MacOS(:x86_64),
]

os = get(ENV, "TRAVIS_OS_NAME", nothing)
if os == "linux"
    platforms = platforms_base
elseif os == "osx"
    platforms = platforms_osx
else
    platforms = [platforms_base; platforms_osx]
end

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libyaml", :libyaml)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: JuliaTagBot/libyaml-binaries

max_issues_repo_name: JuliaTagBot/libyaml-binaries

max_forks_repo_name: JuliaTagBot/libyaml-binaries


text:
using ArgParse

include("functions_align.jl")
include("functions_collapse.jl")
include("functions_collapse_additive.jl")
include("functions_deletions.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--infile", "-i"
            help = "the alignment to collapse, in fasta format"
            required = true
        "--reference"
            help = "the reference sequence, in fasta format"
            required = true
        "--outfile", "-o"
            help = "the name of the reduced alignment to write"
            default = "out.fasta"
        "--rel1"
            help = "output csv with relationships between input sequences, one line per retained sequence"
            default = "tip_to_redundants.csv"
        "--rel2"
            help = "output csv with relationships between input sequences, one line per omitted sequence"
            default = "redundant_to_tips.csv"
        "--retain"
            help = "force the retention of this sequence in the output alignment"
            action = "append_arg"
        "--dels-file"
            help = "file of deletions to type and include as variants when collapsing"
        "--append-as-snp"
            action = "store_true"
            help = "if invoked, and a dels file is provided, append the deletion genotypes as SNPs in the output alignment"
        "--check"
            action = "store_true"
            help = "if invoked, check the sets for not containing any differences"
        "--additive"
            help = "provide a file of previous relationships along with the alignment to collapse
                    - new sequences will be assigned to existing groups if possible, or will form
                    new groups of their own. The full procedure will not be run."
    end

    return parse_args(s)
end

function deletions(parsed_args, ref_array, nucleotide_array, fasta_IDs)
    dels = read_del_file(parsed_args["dels-file"])
    nucleotide_array_dels, A_genotypes_table = type_deletions_and_append_as_SNP(nucleotide_array, dels)

    if Threads.nthreads() == 1
        pairwise_identity = reduced_is_same(nucleotide_array, ref_array)
    else
        pairwise_identity = parallel_reduced_is_same(nucleotide_array, ref_array)
    end

    println("sequence comparisons complete")

    final_sets = get_sets_in_one_go(pairwise_identity)

    levels = TEST_get_levels(final_sets)
    println("number of sequences in input dataset: ", length(fasta_IDs))
    println("number of sequences in output dataset: ", length(levels))

    if parsed_args["check"]
        safe_sets_test = TEST_compare_within_sets(collect(final_sets), nucleotide_array_dels)
    end

    if parsed_args["append-as-snp"]
        tip_to_seq_relationships = get_output_sets(pairwise_identity, final_sets, parsed_args["retain"], nucleotide_array_dels, fasta_IDs)
    else
        tip_to_seq_relationships = get_output_sets(pairwise_identity, final_sets, parsed_args["retain"], nucleotide_array, fasta_IDs)
    end

    write_fasta(tip_to_seq_relationships, nucleotide_array, fasta_IDs, parsed_args["outfile"], false)

    if parsed_args["check"]
        println("safe sets test (true/false): ", safe_sets_test)
    end

    println("number of sequences in output alignment: ", length(keys(tip_to_seq_relationships)))

    write_tip_to_redundant_seq_relationships(tip_to_seq_relationships, parsed_args["rel1"], parsed_args["retain"])

    seq_to_tip_relationships = get_redundant_seq_to_tip_relationships(tip_to_seq_relationships)

    write_redundant_seq_to_tip_relationships(seq_to_tip_relationships, parsed_args["rel2"])
end

function additive(parsed_args, ref_array, nucleotide_array, fasta_IDs)

    println("number of sequences in input dataset: ", length(fasta_IDs))

    old_sets_temp1, old_ID_list_temp1 = read_relationships(parsed_args["additive"])

    # should be a function here to rejig the previous relationships to deal with
    # the case where tips that were in the old dataset are missing in the new one
    old_sets_temp2, old_ID_list_temp2 = update_relationships_based_on_appearance(old_sets_temp1, old_ID_list_temp1, fasta_IDs)

    # Should also probably deal with the potentiality of the sequences (with the same ID)
    # changing from one iteration of the dataset to the next. The test could be the within-sets
    # test of all the old sets, and any sequences that break it could be treated like new sequences.
    # UPDATE: broken sets are thrown out and their sequences will end up in the new ID list:
    old_sets, old_ID_list = update_relationships_based_on_sequences(old_sets_temp2, old_ID_list_temp2, nucleotide_array, fasta_IDs, ref_array)

    old_array, new_array, old_IDs, new_IDs = split_data(nucleotide_array, fasta_IDs, old_ID_list)

    # println("additive mode, number of new sequences to add: ", length(new_IDs))

    # TODO: don't throw an error here - just return the original alignment (minus)
    # any missing seqs + the relationships?
    if length(new_IDs) == 0
        e = error("no new sequences in alignment - did you mean to use --additive?")
        throw(e)
    end

    m_n_pairwise_identity = m_n_reduced_is_same(new_array, old_array, ref_array)

    new_pairwise_identity = reduced_is_same(new_array, ref_array)

    updated_tip_to_seq_relationships, unplaced_new_ID_indices = assign_new_seqs(m_n_pairwise_identity, new_pairwise_identity, old_IDs, new_IDs, old_sets)

    unplaced_IDs = new_IDs[unplaced_new_ID_indices]

    if length(unplaced_new_ID_indices) > 1

        unplaced_pairwise_identity = new_pairwise_identity[unplaced_new_ID_indices, unplaced_new_ID_indices]

        unplaced_sets = get_sets_in_one_go(unplaced_pairwise_identity)

        unplaced_nuc_bit_array = new_array[:,unplaced_new_ID_indices]

        unplaced_tip_to_seq_relationships = get_output_sets(unplaced_pairwise_identity, unplaced_sets, [], unplaced_nuc_bit_array, unplaced_IDs)

        combined_tip_to_seq_relationships = merge(updated_tip_to_seq_relationships, unplaced_tip_to_seq_relationships)

    elseif length(unplaced_IDs) == 1

        combined_tip_to_seq_relationships = merge(updated_tip_to_seq_relationships, Dict{String, Array{String,1}}(unplaced_IDs[1] => []))

    else

        combined_tip_to_seq_relationships = updated_tip_to_seq_relationships
    end

    println("number of sequences in output dataset: ", length(TEST_get_dict_levels(combined_tip_to_seq_relationships)))
    println("number of sequences in output alignment: ", length(keys(combined_tip_to_seq_relationships)))

    if parsed_args["check"]
        final_sets = TEST_get_sets_from_relationships(fasta_IDs, combined_tip_to_seq_relationships)
        safe_sets_test = TEST_compare_within_sets(collect(final_sets), nucleotide_array)
        println("safe sets test (true/false): ", safe_sets_test)
    end

    write_fasta(combined_tip_to_seq_relationships, nucleotide_array, fasta_IDs, parsed_args["outfile"], false)

    write_tip_to_redundant_seq_relationships(combined_tip_to_seq_relationships, parsed_args["rel1"], [])

    combined_seq_to_tip_relationships = get_redundant_seq_to_tip_relationships(combined_tip_to_seq_relationships)

    write_redundant_seq_to_tip_relationships(combined_seq_to_tip_relationships, parsed_args["rel2"])

end

function vanilla(parsed_args, ref_array, nucleotide_array, fasta_IDs)
    if Threads.nthreads() == 1
        pairwise_identity = reduced_is_same(nucleotide_array, ref_array)
    else
        pairwise_identity = parallel_reduced_is_same(nucleotide_array, ref_array)
    end

    println("sequence comparisons complete")

    final_sets = get_sets_in_one_go(pairwise_identity)

    levels = TEST_get_levels(final_sets)
    println("number of sequences in input dataset: ", length(fasta_IDs))
    println("number of sequences in output dataset: ", length(levels))

    if parsed_args["check"]
        safe_sets_test = TEST_compare_within_sets(collect(final_sets), nucleotide_array)
    end

    tip_to_seq_relationships = get_output_sets(pairwise_identity, final_sets, parsed_args["retain"], nucleotide_array, fasta_IDs)

    write_fasta(tip_to_seq_relationships, nucleotide_array, fasta_IDs, parsed_args["outfile"], false)

    if parsed_args["check"]
        println("safe sets test (true/false): ", safe_sets_test)
    end

    println("number of sequences in output alignment: ", length(keys(tip_to_seq_relationships)))

    write_tip_to_redundant_seq_relationships(tip_to_seq_relationships, parsed_args["rel1"], parsed_args["retain"])

    seq_to_tip_relationships = get_redundant_seq_to_tip_relationships(tip_to_seq_relationships)

    write_redundant_seq_to_tip_relationships(seq_to_tip_relationships, parsed_args["rel2"])
end

function main()
    parsed_args = parse_commandline()

    ref_array, ref_id = populate_byte_array_get_names(parsed_args["reference"])
    nucleotide_array, fasta_IDs = populate_byte_array_get_names(parsed_args["infile"])

    if length(fasta_IDs) != length(Set(fasta_IDs))
        e = error("duplicate IDs in input fasta")
        throw(e)
    end

    if length(ref_array) != size(nucleotide_array, 1)
        e = error("reference and alignment do not have the same number of sites")
        throw(e)
    end

    if parsed_args["additive"] != nothing
        if parsed_args["dels-file"] != nothing
            e = error("can't use a dels file with --additive at the moment")
            throw(e)
        end

        additive(parsed_args, ref_array, nucleotide_array, fasta_IDs)

    elseif parsed_args["dels-file"] != nothing

        deletions(parsed_args, ref_array, nucleotide_array, fasta_IDs)

    else

        println("number of threads for sequence comparisons: ", Threads.nthreads())
        vanilla(parsed_args, ref_array, nucleotide_array, fasta_IDs)

    end
end

main()

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
using Luxor, Colors
using Printf
using JSON
using ArgParse

include("strange.jl")

function parse_commandline()
    settings = ArgParseSettings()

    @add_arg_table settings begin
        "--infile"
            help = "parameter input filename"
            arg_type = String
            default = "params.json"
        "--outfile"
            help = "strange image output filename"
            default = "regenerate.png"
	"--n"
	    help = "generate n variations of random colors on input file"
	    default = 1
	    arg_type = Int
	"--red"
	    help = "red midpoint"
	    default = 1.0
	    arg_type = Float64
	"--blue"
	    help = "blue midpoint"
	    default = 1.0
	    arg_type = Float64
	"--green"
	    help = "green midpoint"
	    default = 1.0
	    arg_type = Float64
	"--redvar"
	    help = "red midpoint"
	    default = .50
	    arg_type = Float64
	"--bluevar"
	    help = "blue midpoint"
	    default = .50
	    arg_type = Float64
	"--greenvar"
	    help = "green midpoint"
	    default = .50
	    arg_type = Float64
    end

    return parse_args(settings)
end

# generate n png images

function main()
	parsed_args = parse_commandline()

	infile = parsed_args["infile"]
	outfile = parsed_args["outfile"]

	n = parsed_args["n"]

	red_range = parsed_args["red"]
	green_range = parsed_args["green"]
	blue_range = parsed_args["blue"]

	red_var = parsed_args["redvar"]
	green_var = parsed_args["greenvar"]
	blue_var = parsed_args["bluevar"]

	println("input ctrl file: $infile")

	sp = JSON.parsefile(infile)

	ctrl_params = Dict([
			    :a => sp["ctrl_params"]["a"],
			    :b => sp["ctrl_params"]["b"],
			    :c => sp["ctrl_params"]["c"],
			    :d => sp["ctrl_params"]["d"],
			    :e => sp["ctrl_params"]["e"]])

	dot_size = get(sp, "dot_size", .50)
	res = get(sp, "res", 1000)

	color_profile = Dict([
			      :r => sp["color_profile"]["r"],
			      :g => sp["color_profile"]["g"],
			      :b => sp["color_profile"]["b"]])
	flip = get(sp, "flip", false)
	cmode = get(sp, "cmode", "rgb")

	params = @sprintf("Params: a=%-2.16f, b=%-2.16f, c=%-2.16f, d=%-2.16f, e=%-2.16f, res=%-2.2d, dotsize = %-2.16f flip = %d, cmode = %s\n",
			  ctrl_params[:a],
			  ctrl_params[:b],
			  ctrl_params[:c],
			  ctrl_params[:d],
			  ctrl_params[:e],
			  res,
			  dot_size,
			  flip,
			  cmode)
	println("$params")

	if n > 1

		# override cmode if generating more than one pic
		cmode = "rgb"
		for i in 1:n

			imagefile = @sprintf("%s_%-2.2d.png", outfile, i)
			jsonfile = @sprintf("%s_%-2.2d.json", outfile, i)

			println("outputfile: $imagefile")
			color_profile = Dict([
					      :r => red_range + rand()*rs()*red_range*red_var,
					      :g => green_range + rand()*rs()*green_range*green_var,
					      :b => blue_range + rand()*rs()*blue_range*blue_var
					     ])

			cf = open(jsonfile, "w")
                	sset = StrangeSet(ctrl_params, dot_size, res, color_profile, imagefile, flip, cmode)
                	JSON.print(cf, JSON.parse(JSON.json(sset)), 4)
                	close(cf)

			strange(imagefile, flip, cmode, dot_size, res, ctrl_params, color_profile)
		end
	else
		println("outputfile: $outfile")
		strange(outfile, flip, cmode, dot_size, res, ctrl_params, color_profile)
	end
	println("finished generating $outfile")
end

main()


text:
function geometry(params)

    Ny = params.Ny
    Ly = params.Ly

    Dy,y = cheb(Ny);
       y = y*(0.5*Ly);
      Dy = (2/Ly)*Dy;
     Dy2 = Dy*Dy;

    return Dy, Dy2, y
end