text:
# Canonical form of multivariate normal

import NumericExtensions
import NumericExtensions.ScalMat
import NumericExtensions.PDiagMat
import NumericExtensions.PDMat

using Distributions
using Base.Test

##### construction, basic properties, and evaluation

h = [1., 2., 3.]
dv = [1.2, 3.4, 2.6]
J = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]

x1 = [3.2, 1.8, 2.4]
x = rand(3, 100)

# SGauss

gs = IsoNormalCanon(h, 2.0)
@test isa(gs, IsoNormalCanon)
@test dim(gs) == 3
@test mean(gs) == mode(gs) == h / 2.0
@test !gs.zeromean
@test invcov(gs) == diagm(fill(2.0, 3))
@test cov(gs) == diagm(fill(0.5, 3))
@test var(gs) == diag(cov(gs))
@test_approx_eq entropy(gs) 0.5 * logdet(2π * e * cov(gs))

gsz = IsoNormalCanon(3, 2.0)
@test isa(gsz, IsoNormalCanon)
@test dim(gsz) == 3
@test mean(gsz) == zeros(3)
@test gsz.zeromean

# DGauss

gd = DiagNormalCanon(h, dv)
@test isa(gd, DiagNormalCanon)
@test dim(gd) == 3
@test_approx_eq mean(gd) h ./ dv
@test !gd.zeromean
@test invcov(gd) == diagm(dv)
@test_approx_eq cov(gd) diagm(1.0 / dv)
@test_approx_eq var(gd) diag(cov(gd))
@test_approx_eq entropy(gd) 0.5 * logdet(2π * e * cov(gd))

gdz = DiagNormalCanon(dv)
@test isa(gdz, DiagNormalCanon)
@test dim(gdz) == 3
@test mean(gdz) == zeros(3)
@test gdz.zeromean

# Gauss

gf = MvNormalCanon(h, J)
@test isa(gf, MvNormalCanon)
@test dim(gf) == 3
@test_approx_eq mean(gf) J \ h
@test !gf.zeromean
@test invcov(gf) == J
@test_approx_eq cov(gf) inv(J)
@test_approx_eq var(gf) diag(cov(gf))
@test_approx_eq entropy(gf) 0.5 * logdet(2π * e * cov(gf))

gfz = MvNormalCanon(J)
@test isa(gfz, MvNormalCanon)
@test dim(gfz) == 3
@test mean(gfz) == zeros(3)
@test gfz.zeromean

# conversion

us = convert(IsoNormal, gs)
gs2 = convert(IsoNormalCanon, us)
@test isa(us, IsoNormal)
@test isa(gs2, IsoNormalCanon)
@test dim(gs2) == dim(gs)
@test_approx_eq gs2.h gs.h
@test_approx_eq gs2.J.value gs.J.value

ud = convert(DiagNormal, gd)
gd2 = convert(DiagNormalCanon, ud)
@test isa(ud, DiagNormal)
@test isa(gd2, DiagNormalCanon)
@test dim(gd2) == dim(gd)
@test_approx_eq gd2.h gd.h
@test_approx_eq gd2.J.diag gd.J.diag

uf = convert(MvNormal, gf)
gf2 = convert(MvNormalCanon, uf)
@test isa(uf, MvNormal)
@test isa(gf2, MvNormalCanon)
@test dim(gf2) == dim(gf)
@test_approx_eq gf2.h gf.h
@test_approx_eq gf2.J.mat gf.J.mat

# logpdf evaluation

@test_approx_eq logpdf(gs, x1) logpdf(us, x1)
@test_approx_eq logpdf(gs, x)  logpdf(us, x)
@test_approx_eq logpdf(gsz, x1) logpdf(convert(IsoNormal, gsz), x1)
@test_approx_eq logpdf(gsz, x)  logpdf(convert(IsoNormal, gsz), x)

@test_approx_eq logpdf(gd, x1) logpdf(ud, x1)
@test_approx_eq logpdf(gd, x)  logpdf(ud, x)
@test_approx_eq logpdf(gdz, x1) logpdf(convert(DiagNormal, gdz), x1)
@test_approx_eq logpdf(gdz, x)  logpdf(convert(DiagNormal, gdz), x)

@test_approx_eq logpdf(gf, x1) logpdf(uf, x1)
@test_approx_eq logpdf(gf, x)  logpdf(uf, x)
@test_approx_eq logpdf(gfz, x1) logpdf(convert(MvNormal, gfz), x1)
@test_approx_eq logpdf(gfz, x)  logpdf(convert(MvNormal, gfz), x)

# sampling

x = rand(gs)
@test isa(x, Vector{Float64})
@test length(x) == dim(gs)

x = rand(gd)
@test isa(x, Vector{Float64})
@test length(x) == dim(gd)

x = rand(gf)
@test isa(x, Vector{Float64})
@test length(x) == dim(gf)

x = rand(gsz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gsz)

x = rand(gdz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gdz)

x = rand(gfz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gfz)

n = 10
x = rand(gs, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gs), n)

x = rand(gd, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gd), n)

x = rand(gf, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gf), n)

x = rand(gsz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gsz), n)

x = rand(gdz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gdz), n)

x = rand(gfz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gfz), n)


text:
__precompile__()

module MatFun

using Base.LinAlg, TaylorSeries

export schurparlett, aaa, ratkrylov, ratkrylovf

include("schurparlett.jl")
include("aaa.jl")
include("ratkrylov.jl")

end # module

ext: jl


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

"""
    struct MCMCSampling <: AbstractSamplingAlgorithm

Samples a probability density using Markov chain Monte Carlo.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MCMCSampling{
    AL<:MCMCAlgorithm,
    TR<:AbstractDensityTransformTarget,
    IN<:MCMCInitAlgorithm,
    BI<:MCMCBurninAlgorithm,
    CT<:MCMCConvergenceTest,
    CB<:Function
} <: AbstractSamplingAlgorithm
    mcalg::AL = MetropolisHastings()
    trafo::TR = bat_default(MCMCSampling, Val(:trafo), mcalg)
    nchains::Int = 4
    nsteps::Int = bat_default(MCMCSampling, Val(:nsteps), mcalg, trafo, nchains)
    init::IN = MCMCChainPoolInit(
        nsteps_init = max(div(nsteps, 100), 250)
    )
    burnin::BI = MCMCMultiCycleBurnin(
        nsteps_per_cycle = max(div(nsteps, 10), 2500)
    )
    convergence::CT = BrooksGelmanConvergence()
    strict::Bool = true
    store_burnin::Bool = false
    nonzero_weights::Bool = true
    callback::CB = nop_func
end

export MCMCSampling

function bat_sample_impl(
    rng::AbstractRNG,
    target::AnyDensityLike,
    algorithm::MCMCSampling
)
    density_notrafo = convert(AbstractDensity, target)
    density, trafo = bat_transform(algorithm.trafo, density_notrafo)

    mcmc_algorithm = algorithm.mcalg

    (chains, tuners, chain_outputs) = mcmc_init!(
        rng,
        mcmc_algorithm,
        density,
        algorithm.nchains,
        algorithm.init,
        get_mcmc_tuning(mcmc_algorithm),
        algorithm.nonzero_weights,
        algorithm.store_burnin ? algorithm.callback : nop_func
    )

    if !algorithm.store_burnin
        chain_outputs .= DensitySampleVector.(chains)
    end

    mcmc_burnin!(
        algorithm.store_burnin ? chain_outputs : nothing,
        tuners,
        chains,
        algorithm.burnin,
        algorithm.convergence,
        algorithm.strict,
        algorithm.nonzero_weights,
        algorithm.store_burnin ? algorithm.callback : nop_func
    )

    mcmc_iterate!(
        chain_outputs,
        chains;
        max_nsteps = algorithm.nsteps,
        nonzero_weights = algorithm.nonzero_weights,
        callback = algorithm.callback
    )

    output = DensitySampleVector(first(chains))
    isnothing(output) || append!.(Ref(output), chain_outputs)
    samples_trafo = varshape(density).(output)

    samples_notrafo = inv(trafo).(samples_trafo)

    (result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo, generator = MCMCSampleGenerator(chains))
end


text:
"""
    generate_id() :: Int64

Generate ID. It is unique in most instances.
"""
function generate_id()
    time_value = (now().instant.periods.value - 63749462400000) << 16
    rand_value = rand(UInt16)
    time_value + rand_value
end

format_stdxxx_file(::Nothing) = ""
format_stdxxx_file(x::String) = x
format_stdxxx_file(x::AbstractString) = convert(String, x)
format_stdxxx_file(x::IOStream) = x.name[7:end-1]
format_stdxxx_file(x) = ""

function check_ncpu_mem(ncpu::Int64, mem::Int64)
    if ncpu < 1
        error("ncpu < 1 is not supported for Job")
    end
    if mem < 0
        error("mem < 1 is not supported for Job")
    end
end

"""
    Job(command::Base.AbstractCmd; stdout=nothing, stderr=nothing, append::Bool=false, kwargs...)
    Job(f::Function; kwargs...)
    Job(task::Task; kwargs...)

# Arguments

- `command::Base.AbstractCmd`: the command to run.
- `f::Function`: the function to run without any arguments, like `f()`.
- `task::Task`: the task to run. Eg: `@task(1+1)`.

# Common Keyword Arguments (kwargs...)

 - `name::String = ""`: job name.
 - `user::String = ""`: user that job belongs to.
 - `ncpu::Int64 = 1`: number of CPU this job is about to use.
 - `mem::Int64 = 0`: number of memory this job is about to use (supports TB, GB, MB, KB, B=1).
 - `schedule_time::Union{DateTime,Period} = DateTime(0)`: The expected time to run.
 - `dependency::Vector{Pair{Symbol,Int64}}`: defer job until specified jobs reach specified state (QUEUING, RUNNING, DONE, FAILED, CANCELLED).
 - `wall_time::Period = Week(1)`: wall clock time limit.
 - `priority::Int = 20`: lower means higher priority.

# Experimental Keyword Arguments - Output Redirection:

 - `stdout=nothing`: redirect stdout to the file.
 - `stderr=nothing`: redirect stderr to the file.
 - `append::Bool=false`: append the stdout or stderr or not.

!!! note
    Redirecting in Julia are not thread safe, so unexpected redirection might be happen if you are running programs in different `Tasks` simultaneously (multi-threading).
"""
mutable struct Job
    id::Int64
    name::String
    user::String
    ncpu::Int64
    mem::Int64
    schedule_time::DateTime
    create_time::DateTime
    start_time::DateTime
    stop_time::DateTime
    wall_time::Period
    state::Symbol
    priority::Int
    dependency::Vector{Pair{Symbol,Int64}}
    task::Union{Task,Nothing}
    stdout_file::String
    stderr_file::String
    _thread_id::Int

    # original
    function Job(id::Int64, name::String, user::String, ncpu::Int64, mem::Int64, schedule_time::DateTime, create_time::DateTime, start_time::DateTime, stop_time::DateTime, wall_time::Period, state::Symbol, priority::Int, dependency::Vector{Pair{Symbol,Int64}}, task::Union{Task,Nothing}, stdout_file::String, stderr_file::String)
        if !(typeof(wall_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.wall_time is not one of Union{Millisecond,Second,Minute,Hour,Day,Week}")
        end
        check_ncpu_mem(ncpu, mem)
        new(id, name, user, ncpu, mem, schedule_time, create_time, start_time, stop_time, wall_time, state, priority, dependency, task, stdout_file, stderr_file, 0)
    end

    # schedule_time::Period
    function Job(id::Int64, name::String, user::String, ncpu::Int64, mem::Int64, schedule_time::Period, create_time::DateTime, start_time::DateTime, stop_time::DateTime, wall_time::Period, state::Symbol, priority::Int, dependency::Vector{Pair{Symbol,Int64}}, task::Union{Task,Nothing}, stdout_file::String, stderr_file::String)
        if !(typeof(wall_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.wall_time is not one of Union{Millisecond,Second,Minute,Hour,Day,Week}")
        end
        if !(typeof(schedule_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.schedule_time is not one of Union{DateTime,Millisecond,Second,Minute,Hour,Day,Week}")
        end
        check_ncpu_mem(ncpu, mem)
        new(id, name, user, ncpu, mem, now() + schedule_time, create_time, start_time, stop_time, wall_time, state, priority, dependency, task, stdout_file, stderr_file, 0)
    end

    # dependency::Pair{Symbol,Int64}
    function Job(id::Int64, name::String, user::String, ncpu::Int64, mem::Int64, schedule_time::DateTime, create_time::DateTime, start_time::DateTime, stop_time::DateTime, wall_time::Period, state::Symbol, priority::Int, dependency::Pair{Symbol,Int64}, task::Union{Task,Nothing}, stdout_file::String, stderr_file::String)
        if !(typeof(wall_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.wall_time is not one of Union{Millisecond,Second,Minute,Hour,Day,Week}")
        end
        check_ncpu_mem(ncpu, mem)
        new(id, name, user, ncpu, mem, schedule_time, create_time, start_time, stop_time, wall_time, state, priority, [dependency], task, stdout_file, stderr_file, 0)
    end

    # schedule_time::Period
    # dependency::Pair{Symbol,Int64}
    function Job(id::Int64, name::String, user::String, ncpu::Int64, mem::Int64, schedule_time::Period, create_time::DateTime, start_time::DateTime, stop_time::DateTime, wall_time::Period, state::Symbol, priority::Int, dependency::Pair{Symbol,Int64}, task::Union{Task,Nothing}, stdout_file::String, stderr_file::String)
        if !(typeof(wall_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.wall_time is not one of Union{Millisecond,Second,Minute,Hour,Day,Week}")
        end
        if !(typeof(schedule_time) <: Union{Millisecond,Second,Minute,Hour,Day,Week})
            error("Job.schedule_time is not one of Union{DateTime,Millisecond,Second,Minute,Hour,Day,Week}")
        end
        check_ncpu_mem(ncpu, mem)
        new(id, name, user, ncpu, mem, now() + schedule_time, create_time, start_time, stop_time, wall_time, state, priority, [dependency], task, stdout_file, stderr_file, 0)
    end
end

### Job

function Job(task::Task;
    name::AbstractString = "",
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing, stderr=nothing, append::Bool=false
)
    if ncpu > 1
        @warn "ncpu != 1 for Job(task::Task) is not fully supported by JobScheduler. If a task uses multi-threads, it is recommended to split it into different jobs. Job: $name." maxlog=1
    end
    task2 = @task Pipelines.redirect_to_files(stdout, stderr; mode = append ? "a+" : "w+") do
        task.code()
    end
    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task2, "", "")
end

function Job(f::Function;
    name::AbstractString = "",
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing, stderr=nothing, append::Bool=false
)
    if ncpu != 1
        @warn "ncpu != 1 for Job(task::Task) is not fully supported by JobScheduler. If a task uses multi-threads, it is recommended to split it into different jobs. Job: $name." maxlog=1
    end
    task2 = @task Pipelines.redirect_to_files(stdout, stderr; mode = append ? "a+" : "w+") do
        f()
    end
    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task2, "", "")
end

function Job(command::Base.AbstractCmd;
    name::AbstractString = "",
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing, stderr=nothing, append::Bool=false
)
    task = @task Pipelines.redirect_to_files(stdout, stderr; mode = append ? "a+" : "w+") do
        run(command)
    end
    stdout_file = format_stdxxx_file(stdout)
    stderr_file = format_stdxxx_file(stderr)
    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task, stdout_file, stderr_file)
end

function result(job::Job)
    if job.state !== DONE
        @warn "Getting result from a $(job.state) job: returned value might be unexpected."
    end
    job.task.result
end


text:
function binaryRepresentation(arrayTuple, maxSymbol)
    if (maximum(arrayTuple) > maxSymbol)
        # Element in the array exceeding the max symbol index, therefore not a valid card
        return 0;
    end
    if in(0, arrayTuple)
        # Array representation contains a zero, which is not a valid symbol
        return 0;
    end
	len = length(arrayTuple);
	ulen = length(unique(arrayTuple));
	if len != ulen
        # Duplicate symbols found on the card
		return 0;
    end
	return sum(2.^arrayTuple) >>> 1;
end

max_stars_repo_name: driesverhees/dobble

max_issues_repo_name: driesverhees/dobble

max_forks_repo_name: driesverhees/dobble


text:
using DiffEqBase

struct CTKAlg <: DiffEqBase.AbstractSteadyStateAlgorithm end

"""
    solve(prob::DiffEqBase.AbstractSteadyStateProblem,
          alg::CTKAlg;
          nrm=norm,
          τstop=1e12*365*24*60*60,
          preprint="",
          maxItNewton=50)

Solves `prob` using the modified C.T.Kelley Shamanskii algorithm.
"""
function DiffEqBase.solve(prob::DiffEqBase.AbstractSteadyStateProblem,
                          alg::CTKAlg;
                          nrm=norm,
                          τstop=ustrip(1e6u"Myr"|>u"s"),
                          preprint="",
                          maxItNewton=50)
    # Define the functions according to DiffEqBase.SteadyStateProblem type
    p = prob.p
    t = 0
    x0 = copy(prob.u0)
    dx, df = copy(x0), copy(x0)
    F(x) = prob.f(dx, x, p, t)
    ∇ₓF(x) = prob.f(df, dx, x, p, t)
    # Compute `u_steady` and `resid` as per DiffEqBase using my algorithm
    x_steady = NewtonChordShamanskii(F, ∇ₓF, nrm, x0, τstop; preprint=preprint, maxItNewton=maxItNewton)
    resid = F(x_steady)
    # Return the common DiffEqBase solution type
    DiffEqBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
end

"""
    SteadyStateProblem(F, ∇ₓF, x, p)

Returns the `SteadyStateProblem` defined by `F(x,p)=0`.
"""
function DiffEqBase.SteadyStateProblem(F, ∇ₓF, x, p)
    f(dx, x, p, t) = F(x, p)
    f(df, dx, x, p, t) = ∇ₓF(x, p)
    return DiffEqBase.SteadyStateProblem(f, x, p)
end

export solve, SteadyStateProblem, CTKAlg, CTKAlg2


text:
using NCDatasets

const INFILE = "SOCATv5.tsv"
const OUTFILE = "daily.nc"
const CELLSIZE = 2.5
const STARTYEAR = 1985
const ENDYEAR = 2016

# Zero-based day of year of each month
const MONTHSTARTS = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]
const LEAPMONTHSTARTS = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]

# Days of the year
DAYSTARTS = range(1, step=1, length=365)

# Calculate the conversion of leap year days to 365 'normal' year days
# Each leap year day lasts 1 1/365 normal days
LEAPDAYSTARTS = Array{Float64}(undef, 365)

for i in 1:365
    LEAPDAYSTARTS[i] = i + (((366 / 365) / 365) * (i - 1))
end

# Determine whether or not a given year is a leap year
function isleapyear(year::Int64)::Bool
    local leapyear::Bool = false
    if rem(year, 4) == 0
        leapyear = true
        if rem(year, 100) == 0
            if rem(year, 400) != 0
                leapyear = false
            end
        end
    end

    leapyear
end

# Calculate the Julian day from a year/month/day set
function calcjdate(year::Int64, month::Int64, day::Int64)::Int64
    local result::Int64 = 0

    if isleapyear(year)
        result = LEAPMONTHSTARTS[month]
    else
        result = MONTHSTARTS[month]
    end

    result = result + day
end

# Calculate the day of the year for a given date
function getdayindex(date, days)::Float64
    return findall(days .<= date)[end]
end

# Calcaulte the nth day of the complete data set
# Days are only calculated between the start and end years
function getdateindex(year::Int64, month::Int64, day::Int64)::Int64

    local jdate::Float64 = calcjdate(year, month, day)
    local index::Float64 = -1
    local dayindex::Float64 = 0

    if year >= STARTYEAR && year <= ENDYEAR
        if isleapyear(year)
            dayindex = getdayindex(jdate, DAYSTARTS)
        else
            dayindex = getdayindex(jdate, LEAPDAYSTARTS)
        end

        index = ((year - STARTYEAR) * 365) + dayindex
    end

    return convert(Int64, floor(index))
end

# Calculate the cell indices from a longitude and latitude
function getcellindex(longitude::Float64, latitude::Float64)::Tuple{Int64, Int64}
    local loncell::Int64 = floor(longitude / CELLSIZE) + 1
    local latcell::Int64 = floor((latitude + 90) / CELLSIZE) + 1

    if loncell == 145
        loncell = 1
    end

    return loncell, latcell
end

function applydataset(datasettotals::Array{Float64, 3}, datasetcounts::Array{Int64, 3}, overalltotals::Array{Float64, 3}, overallcounts::Array{Int64, 3})
    datasetmean::Array{Float64, 3} = datasettotals ./ datasetcounts
    datasetmean[isnan.(datasetmean)] .= 0

    overalltotals .= overalltotals .+ datasetmean
    datasetmean[datasetmean .> 0] .= 1
    overallcounts .= overallcounts .+ datasetmean
end

function run()

    totaldays::Int64 = (ENDYEAR - STARTYEAR + 1) * 365

    # Output data set
    overallcelltotals::Array{Float64, 3} = zeros(convert(Int64, 360 / CELLSIZE), convert(Int64, 180 / CELLSIZE), totaldays)
    overallcellcounts::Array{Int64, 3} = zeros(convert(Int64, 360 / CELLSIZE), convert(Int64, 180 / CELLSIZE), totaldays)

    # Open input file
    inchan::IOStream = open(INFILE)

    currentdataset::String = ""
    datasetcelltotals::Array{Float64, 3} = zeros(convert(Int64, 360 / CELLSIZE), convert(Int64, 180 / CELLSIZE), totaldays)
    datasetcellcounts::Array{Int64, 3} = zeros(convert(Int64, 360 / CELLSIZE), convert(Int64, 180 / CELLSIZE), totaldays)

    currentline::String = readline(inchan)
    linecount::Int64 = 1
    while length(currentline) > 0

        fields::Array{String, 1} = split(currentline, "\t")

        dataset::String = fields[1]

        if dataset != currentdataset
            if length(currentdataset) > 0
                applydataset(datasetcelltotals, datasetcellcounts, overallcelltotals, overallcellcounts)

                print("\033[1K\r$currentdataset ($linecount)")

                datasetcelltotals .= 0
                datasetcellcounts .= 0
            end

            currentdataset = dataset
        end

        year::Int64 = parse(Int64, fields[5])
        month::Int64 = parse(Int64, fields[6])
        day::Int64 = parse(Int64, fields[7])
        longitude::Float64 = parse(Float64, fields[11])
        latitude::Float64 = parse(Float64, fields[12])
        fco2::Float64 = parse(Float64, fields[24])

        cellindex::Tuple{Int64, Int64} = getcellindex(longitude, latitude)

        dateindex::Int64 = getdateindex(year, month, day)

        if dateindex != -1
            datasetcelltotals[cellindex[1], cellindex[2], dateindex] = datasetcelltotals[cellindex[1], cellindex[2], dateindex] + fco2
            datasetcellcounts[cellindex[1], cellindex[2], dateindex] = datasetcellcounts[cellindex[1], cellindex[2], dateindex] + 1
        end

        currentline = readline(inchan)
        linecount = linecount + 1
    end

    close(inchan)

    # The last dataset
    applydataset(datasetcelltotals, datasetcellcounts, overallcelltotals, overallcellcounts)

    # Overall cell means
    local meanfco2::Array{Float64, 3} = overallcelltotals ./ overallcellcounts
    meanfco2 = replace(meanfco2, NaN=>-1e35)

    # Write NetCDF
    nc = Dataset(OUTFILE, "c")
    defDim(nc, "longitude", trunc(Int, (360 / CELLSIZE)))
    defDim(nc, "latitude", trunc(Int, (180 / CELLSIZE)))
    defDim(nc, "time", totaldays)

    nclon = defVar(nc, "longitude", Float32, ("longitude",))
    nclon.attrib["units"] = "degrees_east"

    nclat = defVar(nc, "latitude", Float32, ("latitude",))
    nclat.attrib["units"] = "degrees_north"

    nctime = defVar(nc, "time", Float32, ("time",))
    nctime.attrib["calendar"] = "noleap"

    ncfco2 = defVar(nc, "fCO2", Float64, ("longitude", "latitude", "time"))
    ncfco2.attrib["_FillValue"] = -1e35

    nclon[:] = collect(range(CELLSIZE / 2, step=CELLSIZE, stop=360))
    nclat[:] = collect(range(-90 + CELLSIZE / 2, step=CELLSIZE, stop=90))
    nctime[:] = collect(range(STARTYEAR, step=(1/365), stop=(ENDYEAR + 1) - (1/365)))
    ncfco2[:,:,:] = meanfco2

    close(nc)
end

@time run()


text:
module LinAlgBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat
using Compat.Iterators

if VERSION >= v"0.7.0-DEV.3449"
    using LinearAlgebra
    using LinearAlgebra: UnitUpperTriangular
else
    using Base.LinAlg
    using Base.LinAlg: UnitUpperTriangular
    const LinearAlgebra = Base.LinAlg
end

const SUITE = BenchmarkGroup(["array"])

const SIZES = (2^8, 2^10)
const MATS = (Matrix, Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal, UpperTriangular, LowerTriangular)
const DIVMATS = Iterators.filter(x -> !(in(x, (Bidiagonal, Tridiagonal, SymTridiagonal))), MATS)

typename(::Type{T}) where {T} = string(isa(T,DataType) ? T.name : Base.unwrap_unionall(T).name)
typename(::Type{M}) where {M<:Matrix} = "Matrix"
typename(::Type{V}) where {V<:Vector} = "Vector"

const UPLO = VERSION >= v"0.7.0-DEV.884" ? :U : true

linalgmat(::Type{Matrix}, s) = randmat(s)
linalgmat(::Type{Diagonal}, s) = Diagonal(randvec(s))
linalgmat(::Type{Bidiagonal}, s) = Bidiagonal(randvec(s), randvec(s-1), UPLO)
linalgmat(::Type{Tridiagonal}, s) = Tridiagonal(randvec(s-1), randvec(s), randvec(s-1))
linalgmat(::Type{SymTridiagonal}, s) = SymTridiagonal(randvec(s), randvec(s-1))
linalgmat(::Type{UpperTriangular}, s) = UpperTriangular(randmat(s))
linalgmat(::Type{LowerTriangular}, s) = LowerTriangular(randmat(s))
linalgmat(::Type{UnitUpperTriangular}, s) = UnitUpperTriangular(randmat(s))

function linalgmat(::Type{Hermitian}, s)
    A = randmat(s)
    A = A + im*A
    return Hermitian(A + A')
end

# Non-positive-definite upper-triangular matrix
mutable struct NPDUpperTriangular
end
function linalgmat(::Type{NPDUpperTriangular}, s)
    A = linalgmat(UpperTriangular, s)
    rr = samerand(s)
    for i in 1:s
        A[i,i] = (2*rr[i]-1)*A[i,i]
    end
    return A
end
typename(::Type{NPDUpperTriangular}) = "NPDUpperTriangular"

# Hermitian Sparse Matrix With Nonzero Pivots
mutable struct HermitianSparseWithNonZeroPivots
end
function linalgmat(::Type{HermitianSparseWithNonZeroPivots}, s)
    A = samesprand(s, s, 1/s)
    A = A + A' + I
    A
end
typename(::Type{HermitianSparseWithNonZeroPivots}) = "HermitianSparseWithNonZeroPivots"

############################
# matrix/vector arithmetic #
############################

g = addgroup!(SUITE, "arithmetic")

for s in SIZES
    v = randvec(s)
    vstr = typename(Vector)
    g["-", vstr, vstr, s] = @benchmarkable -($v, $v)
    g["+", vstr, vstr, s] = @benchmarkable +($v, $v)
    for M in MATS
        mstr = typename(M)
        m = linalgmat(M, s)
        g["*", mstr, vstr, s]  = @benchmarkable *($m, $v)
        g["\\", mstr, vstr, s] = @benchmarkable \($m, $v)
        g["+", mstr, mstr, s]  = @benchmarkable +($m, $m)
        g["-", mstr, mstr, s]  = @benchmarkable -($m, $m)
        g["*", mstr, mstr, s]  = @benchmarkable *($m, $m)
    end
    for M in DIVMATS
        mstr = typename(M)
        m = linalgmat(M, s)
        g["/", mstr, mstr, s]  = @benchmarkable /($m, $m)
        g["\\", mstr, mstr, s] = @benchmarkable \($m, $m)
    end
    # Issue #14722
    C = zeros(Float32, s, s)
    A = randmat(s)
    B = randmat(s)
    if VERSION >= v"0.7.0-DEV.3204"
        g["mul!", "Matrix{Float32}", "Matrix{Float64}", "Matrix{Float64}", s] = @benchmarkable LinearAlgebra.mul!($C, $A, $B)
    else
        g["mul!", "Matrix{Float32}", "Matrix{Float64}", "Matrix{Float64}", s] = @benchmarkable A_mul_B!($C, $A, $B)
    end

    for T in [Int32, Int64, Float32, Float64]
        arr = samerand(T, s)
        g["cumsum!", T, s] = @benchmarkable cumsum!($arr, $(arr.÷length(arr)))
    end

    for M in (UpperTriangular, UnitUpperTriangular, NPDUpperTriangular, Hermitian)
        mstr = typename(M)
        m = linalgmat(M, s)
        if VERSION >= v"0.7.0-DEV.1599"
            g["sqrt", mstr, s] = @benchmarkable sqrt($m)
        else
            g["sqrt", mstr, s] = @benchmarkable sqrtm($m)
        end
        if M == Hermitian
            if VERSION >= v"0.7.0-DEV.1597"
                g["log", mstr, s] = @benchmarkable log($m)
            else
                g["log", mstr, s] = @benchmarkable logm($m)
            end
            if VERSION >= v"0.7.0-DEV.1486"
                g["exp", mstr, s] = @benchmarkable exp($m)
            else
                g["exp", mstr, s] = @benchmarkable expm($m)
            end
        end
    end

    # PR 21165
    for M in (HermitianSparseWithNonZeroPivots,)
        mstr = typename(M)
        m = linalgmat(M, s)
        g["\\", mstr, vstr, s] = @benchmarkable \($m, $v)
    end

end

for b in values(g)
    b.params.time_tolerance = 0.45
    b.params.samples = 100
    b.params.seconds = 20
end

##################
# factorizations #
##################

g = addgroup!(SUITE, "factorization", ["eig", "svd", "lu", "qr", "schur", "chol"])

for M in (Matrix, Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, LowerTriangular)
    mstr = typename(M)
    for s in SIZES
        m = linalgmat(M, s)
        if VERSION < v"0.7.0-DEV.5211"
            g["eigen", mstr, s] = @benchmarkable eigfact($m)
        else
            g["eigen", mstr, s] = @benchmarkable eigen($m)
        end
    end
end

for M in (Matrix, Diagonal, Bidiagonal, UpperTriangular, LowerTriangular)
    mstr = typename(M)
    for s in SIZES
        m = linalgmat(M, s)
        if VERSION < v"0.7.0-DEV.5211"
            g["svd", mstr, s] = @benchmarkable svdfact($m)
        else
            g["svd", mstr, s] = @benchmarkable svd($m)
        end
    end
end

for M in (Matrix, Tridiagonal)
    mstr = typename(M)
    for s in SIZES
        m = linalgmat(M, s)
        if VERSION < v"0.7.0-DEV.5211"
            g["lu", mstr, s] = @benchmarkable lufact($m)
        else
            g["lu", mstr, s] = @benchmarkable lu($m)
        end
    end
end

for s in SIZES
    mstr = typename(Matrix)
    m = randmat(s)
    arr = m' * m
    if VERSION < v"0.7.0-DEV.5211"
        g["cholesky", mstr, s] = @benchmarkable cholfact($arr)
        g["schur",    mstr, s] = @benchmarkable schurfact($m)
        g["qr",       mstr, s] = @benchmarkable qrfact($m)
    else
        g["cholesky", mstr, s] = @benchmarkable cholesky($arr)
        g["schur",    mstr, s] = @benchmarkable schur($m)
        g["qr",       mstr, s] = @benchmarkable qr($m)
    end
end

for b in values(g)
    b.params.time_tolerance = 0.45
end

########
# BLAS #
########

g = addgroup!(SUITE, "blas")

s = 1024
C = Complex{Float64}
n = samerand()
cn = samerand(C)
v = randvec(s)
cv = randvec(C, s)
m = randmat(s)
cm = randmat(C, s)

g["dot"]       = @benchmarkable BLAS.dot($s, $v, 1, $v, 1)
g["dotu"]      = @benchmarkable BLAS.dotu($s, $cv, 1, $cv, 1)
g["dotc"]      = @benchmarkable BLAS.dotc($s, $cv, 1, $cv, 1)
g["blascopy!"] = @benchmarkable BLAS.blascopy!($s, $v, 1, $v, 1)
g["nrm2"]      = @benchmarkable BLAS.nrm2($s, $v, 1)
g["asum"]      = @benchmarkable BLAS.asum($s, $v, 1)
g["axpy!"]     = @benchmarkable BLAS.axpy!($n, $v, fill!($v, $n))
g["scal!"]     = @benchmarkable BLAS.scal!($s, $n, fill!($v, $n), 1)
g["scal"]      = @benchmarkable BLAS.scal($s, $n, $v, 1)
g["ger!"]      = @benchmarkable BLAS.ger!($n, $v, $v, fill!($m, $n))
g["syr!"]      = @benchmarkable BLAS.syr!('U', 1.0, $v, fill!($m, $n))
g["syrk!"]     = @benchmarkable BLAS.syrk!('U', 'N', $n, $m, $n, fill!($m, $n))
g["syrk"]      = @benchmarkable BLAS.syrk('U', 'N', $n, $m)
g["her!"]      = @benchmarkable BLAS.her!('U', $n, $cv, fill!($cm, $cn))
g["herk!"]     = @benchmarkable BLAS.herk!('U', 'N', $n, $cm, $n, fill!($cm, $cn))
g["herk"]      = @benchmarkable BLAS.herk('U', 'N', $n, $cm)
g["gbmv!"]     = @benchmarkable BLAS.gbmv!('N', $s, 500, 500, $n, $m, $v, $n, fill!($v, $n))
g["gbmv"]      = @benchmarkable BLAS.gbmv('N', $s, 500, 500, $n, $m, $v)
g["sbmv!"]     = @benchmarkable BLAS.sbmv!('U', $s-1, $n, $m, $v, $n, fill!($v, $n))
g["sbmv"]      = @benchmarkable BLAS.sbmv('U', $s-1, $n, $m, $v)
g["gemm!"]     = @benchmarkable BLAS.gemm!('N', 'N', $n, $m, $m, $n, fill!($m, $n))
g["gemm"]      = @benchmarkable BLAS.gemm('N', 'N', $n, $m, $m)
g["gemv!"]     = @benchmarkable BLAS.gemv!('N', $n, $m, $v, $n, fill!($v, $n))
g["gemv"]      = @benchmarkable BLAS.gemv('N', $n, $m, $v)
g["symm!"]     = @benchmarkable BLAS.symm!('L', 'U', $n, $m, $m, $n, fill!($m, $n))
g["symm"]      = @benchmarkable BLAS.symm('L', 'U', $n, $m, $m)
g["symv!"]     = @benchmarkable BLAS.symv!('U', $n, $m, $v, $n, fill!($v, $n))
g["symv"]      = @benchmarkable BLAS.symv('U', $n, $m, $v)
g["trmm!"]     = @benchmarkable BLAS.trmm!('L', 'U', 'N', 'N', $n, $m, fill!($m, $n))
g["trmm"]      = @benchmarkable BLAS.trmm('L', 'U', 'N', 'N', $n, $m, $m)
g["trsm!"]     = @benchmarkable BLAS.trsm!('L', 'U', 'N', 'N', $n, $m, fill!($m, $n))
g["trsm"]      = @benchmarkable BLAS.trsm('L', 'U', 'N', 'N', $n, $m, $m)
g["trmv!"]     = @benchmarkable BLAS.trmv!('L', 'N', 'U', $m, fill!($v, $n))
g["trmv"]      = @benchmarkable BLAS.trmv('L', 'N', 'U', $m, $v)
g["trsv!"]     = @benchmarkable BLAS.trsv!('U', 'N', 'N', $m, fill!($v, $n))
g["trsv"]      = @benchmarkable BLAS.trsv('U', 'N', 'N', $m, $v)

for b in values(g)
    b.params.time_tolerance = 0.40
end

#############
# small exp #
#############

SUITE["small exp #29116"] = @benchmarkable exp([1. 0; 2 0])

end # module


text:
using Test, Interfaces, Tables

abstract type Iterable end

@interface Iterable begin

    iterate(::Iterable)
    iterate(::Iterable, st)

    Base.IteratorSize(::Type{Iterable})::Union{Base.HasLength, Base.HasShape, Base.IsInfinite, Base.SizeUnknown}

    if Base.IteratorSize(Iterable) == Base.HasLength()
        length(::Iterable)
    elseif Base.IteratorSize(Iterable) isa Base.HasShape
        length(::Iterable)
        size(::Iterable)
    end

    Base.IteratorEltype(::Type{Iterable})::Union{Base.HasEltype, Base.EltypeUnknown}

    if Base.IteratorEltype(Iterable) == Base.HasEltype()
        eltype(::Iterable) || eltype(::Type{Iterable})
    end

end

@test Interfaces.implements(Int, Iterable, [Base])

abstract type Indexable end

@interface Indexable begin
    getindex(::Indexable, i)
    setindex!(::Indexable, v, i)
    firstindex(::Indexable)
    lastindex(::Indexable)
end

@test Interfaces.implements(Array, Indexable, [Base])

@interface AbstractArray T begin

    T <: AbstractArray
    size(::T)
    Base.IndexStyle(::Type{T})::Union{Base.IndexLinear, Base.IndexCartesian}

    if Base.IndexStyle(T) == Base.IndexLinear()
        getindex(::T, i::Int)
        setindex!(::T, v, i::Int)
    elseif Base.IndexStyle(T) == Base.IndexCartesian()
        getindex(::T, I::Vararg{Int})
        setindex!(::T, v, I::Vararg{Int})
    end

end

@test Interfaces.implements(Array, AbstractArray, [Base])

struct Foo end

@test_throws InterfaceImplementationError Interfaces.implements(Foo, AbstractArray)

@test_throws InterfaceImplementationError Interfaces.implements(Foo, Iterable)

@interface Tables.AbstractColumns T begin
    Tables.getcolumn(::T, ::Int)
    Tables.getcolumn(::T, ::Symbol)
    Tables.columnnames(::T)
end

@interface Tables.AbstractRow T begin
    Tables.getcolumn(::T, ::Int)
    Tables.getcolumn(::T, ::Symbol)
    Tables.columnnames(::T)
end

abstract type AbstractTable end

@interface AbstractTable begin
    Tables.isrowtable(::Type{AbstractTable}) ||
        Tables.rowaccess(::Type{AbstractTable}) ||
            Tables.columnaccess(::Type{AbstractTable})
    if Tables.isrowtable(AbstractTable) || Tables.rowaccess(AbstractTable)
        Tables.rows(::AbstractTable)::RT where {RT <: Iterable}
        if Base.IteratorEltype(RT) == Base.HasEltype()
            eltype(::RT)::Type{Tables.AbstractRow} ||
                eltype(::Type{RT})::Type{Tables.AbstractRow}
        end
    elseif Tables.columnaccess(AbstractTable)
        Tables.istable(::Type{AbstractTable})
        Tables.columns(::AbstractTable)::Tables.AbstractColumns
    end

    @optional Tables.schema(::AbstractTable)::Union{Nothing, Tables.Schema}
end

@test Interfaces.implements(Tables.DictRowTable, AbstractTable)
@test Interfaces.implements(Tables.DictColumnTable, AbstractTable)

abstract type Example end

@interface Example T begin

    # required subtyping
    T <: Example

    # required method definition
    foo1(arg::Example, arg2::Int)

    # method def w/ return type
    foo2(arg::Example)::Float64

    # method def w/ interface return type
    foo3(arg::Example)::Iterable

    # one of many required
    foo4(arg::Example) || foo5(arg::Example) || foo6(arg::Example)

    # conditional
    if foo7(T)
        foo8(arg::Example)
        foo9(arg::Example)::RT where {RT <: Iterable}
        foo10(::RT)::Int
    elseif foo11(T)
        foo12(arg::Example)
    else
        foo13(arg::Example) || foo14(arg::Example)
    end

    @optional foo15(arg::Example)

end

@test repr("text/plain", Interfaces.interface(Example)) == """
Interface: Example
  * subtyping: `T <: Example`
  * method definition: `foo1(arg::Example, arg2::Int)`
  * method definition with required return type: `foo2(arg::Example)::Float64`
  * method definition with required return type: `foo3(arg::Example)::Iterable`
  * one of the following required:
    * method definition: `foo4(arg::Example)`
    * method definition: `foo5(arg::Example)`
    * method definition: `foo6(arg::Example)`
  * conditional requirements:
    * if foo7(Example)
      * method definition: `foo8(arg::Example)`
      * method definition with required return type: `foo9(arg::Example)::(RT where RT <: Iterable)`
      * method definition with required return type: `foo10(::RT)::Int`
    * elseif foo11(Example)
      * method definition: `foo12(arg::Example)`
    * else
      * one of the following required:
        * method definition: `foo13(arg::Example)`
        * method definition: `foo14(arg::Example)`
  * optional interface requirement:
    * method definition: `foo15(arg::Example)`"""

abstract type AbstractStridedArray end

@interface AbstractStridedArray A begin
    A <: AbstractArray
    strides(::A)
    Base.unsafe_convert(::Type{<:Ptr}, ::A)
    Base.elsize(::Type{<:A})
    @optional stride(::A, i::Int)
end

@test Interfaces.implements(Vector{Int}, AbstractStridedArray, [Base])


text:
using Test
using FormatSpecimens

@testset "FormatSpecimens" begin

@testset "Parse" begin

    path = joinpath(@__DIR__, "..")

    exclusions = [".git", "src", "test", "juliamnt", ".github"]

    formats = filter(str -> isdir(joinpath(path, str)) && !in(str, exclusions), readdir(path))

    for format in formats
        @info "Parsing:" format
        @test_nowarn FormatSpecimens.read_all_specimens(format)
    end

end #testest Parse

end #testest FormatSpecimens

max_stars_repo_name: BioJulia/FormatSpecimens

max_issues_repo_name: BioJulia/FormatSpecimens

max_forks_repo_name: BioJulia/FormatSpecimens


text:
abstract FieldVector{T} <: StaticVector{T}

# Is this a good idea?? Should people just define constructors that accept tuples?
@inline (::Type{FV}){FV<:FieldVector}(x::Tuple) = FV(x...)

@pure size{FV<:FieldVector}(::Union{FV,Type{FV}}) = (length(FV.types),)
@pure length{FV<:FieldVector}(::Union{FV,Type{FV}}) = length(FV.types)

@inline getindex(v::FieldVector, i::Integer) = getfield(v, i)
@inline setindex!(v::FieldVector, x, i::Integer) = setfield!(v, i, x)

@inline function Base.unsafe_convert{T}(::Type{Ptr{T}}, v::FieldVector{T})
    Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(v))
end

# For example:
#
# immutable/type Point3D <: FieldVector{Float64}
#     x::Float64
#     y::Float64
#     z::Float64
# end


text:
@testset "NumField/NonSimpleNumField" begin
  include("NonSimpleNumField/Conjugates.jl")
end


text:
const MOI = MathOptInterface
const SUCCESS = [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED]

mutable struct Moreau{SpecializedType}
    dynamics::Function
    gap::Function                   # Gap function -- provide the force matrices
    hcon::Function                  # Holonomic constraint function
    jac::Function                   # Jacobian of hcon: del hcon / del q
    jacdot::Function                # Time-derivative of jac
    M::Array{SpecializedType, 2}
    h::Array{SpecializedType, 1}
    ϕ::Array{SpecializedType, 1}
    J::Array{SpecializedType, 2}
    Jdot::Array{SpecializedType, 2}
    tA::SpecializedType
    tM::SpecializedType
    tE::SpecializedType
    qA::Array{SpecializedType, 1}
    uA::Array{SpecializedType, 1}
    qM::Array{SpecializedType, 1}
    qE::Array{SpecializedType, 1}
    uE::Array{SpecializedType, 1}
    extra_state::Array{SpecializedType, 1}
    g::Array{SpecializedType, 1}            # Normal distance of contact
    W::Array{SpecializedType, 2}            # Jacobian transpose in the normal direction (changes during integration)
    Wn::Array{SpecializedType, 2}           # Jacobian transpose in the normal direction (fixed, does not change)
    Λ::Array{SpecializedType, 1}            # Normal contact force
    μ::Array{SpecializedType, 1}            # Holonomic constraint forces
    H::SortedSet{Int64, Base.Order.ForwardOrdering} # Which contacts are active?
    Δt::SpecializedType
    ε::SpecializedType                      # Coefficient of restitution (normal dir.)
    ϵ::SpecializedType                      # Baumgarte stabilization constant
end

"""
n: degrees of freedom of the system
m: number of unilateral constraints
"""

function Moreau(gap::Function, dynamics::Function, q::AbstractArray, u::AbstractArray,
        Δt::Wildcard=1e-3, ε::Wildcard=0.5) where {Wildcard <: Real}
    qA = q
    uA = u
    extra_state = Wildcard[]
    n = length(qA)
    qM = _compute_mid_displacements(qA, uA, Δt)
    M, h = dynamics(qM, uA)
    ϕ = Wildcard[]
    J = Array{Wildcard, 2}(undef, 0, 0)
    Jdot = Array{Wildcard, 2}(undef, 0, 0)
    qE = zeros(Wildcard, n)
    uE = zeros(Wildcard, n)
    tA, tM, tE = zeros(Wildcard, 3)
    tM = _compute_mid_time(tA, Δt)
    tE = _compute_mid_time(tM, Δt)
    ϵ = 0.01
    g, Wn = gap(qA, uA)
    W = Wn
    H = SortedSet(Int[])
    Λ = zeros(Wildcard, length(H))
    μ = zeros(Wildcard, length(ϕ))

    Moreau{Wildcard}(dynamics, gap, x->x, x->x, x->x,
        M, h, ϕ, J, Jdot, tA, tM, tE, qA, uA, qM, qE, uE, extra_state,
        g, W, Wn, Λ, μ, H, Δt, ε, ϵ)
end

function Moreau(gap::Function, dynamics::Function, hcon::Function,
        jac::Function, jacdot::Function, q::AbstractArray, u::AbstractArray,
        Δt::Wildcard=1e-3, ε::Wildcard=0.5) where {Wildcard <: Real}
    qA = q
    uA = u
    extra_state = Wildcard[]
    n = length(qA)
    qM = _compute_mid_displacements(qA, uA, Δt)
    M, h = dynamics(qM, uA)
    ϕ = hcon(qM)
    J = jac(qM)
    Jdot = jacdot(qM, uA)
    qE = zeros(Wildcard, n)
    uE = zeros(Wildcard, n)
    tA, tM, tE = zeros(Wildcard, 3)
    tM = _compute_mid_time(tA, Δt)
    tE = _compute_mid_time(tM, Δt)
    ϵ = 0.01
    g, Wn = gap(qA, uA)
    W = Wn
    H = SortedSet(Int[])
    Λ = zeros(Wildcard, length(H))
    μ = zeros(Wildcard, length(ϕ))

    Moreau{Wildcard}(dynamics, gap, hcon, jac, jacdot,
        M, h, ϕ, J, Jdot, tA, tM, tE, qA, uA, qM, qE, uE, extra_state,
        g, W, Wn, Λ, μ, H, Δt, ε, ϵ)
end

function Moreau(gap::Function, dynamics::Function, q::AbstractArray, u::AbstractArray,
    extra_state::AbstractArray, Δt::Wildcard=1e-3, ε::Wildcard=0.1) where {Wildcard <: Real}
    qA = q
    uA = u
    n = length(qA)
    qM = _compute_mid_displacements(qA, uA, Δt)
    M, h = dynamics(qM, uA, extra_state)
    ϕ = Wildcard[]
    J = Array{Wildcard, 2}(undef, 0, 0)
    Jdot = Array{Wildcard, 2}(undef, 0, 0)
    qE = zeros(Wildcard, n)
    uE = zeros(Wildcard, n)
    tA, tM, tE = zeros(Wildcard, 3)
    tM = _compute_mid_time(tA, Δt)
    tE = _compute_mid_time(tM, Δt)
    ϵ = 0.01
    g, Wn = gap(qA, uA)
    W = Wn
    H = SortedSet(Int[])
    Λ = zeros(Wildcard, length(H))
    μ = zeros(Wildcard, length(ϕ))

    Moreau{Wildcard}(dynamics, gap, x->x, x->x, x->x,
        M, h, ϕ, J, Jdot, tA, tM, tE, qA, uA, qM, qE, uE, extra_state,
        g, W, Wn, Λ, μ, H, Δt, ε, ϵ)
end

function _compute_mid_time(t::Number, Δt::Real)
    return t + 1/2*Δt
end

function _compute_mid_displacements(q::AbstractArray, u::AbstractArray, Δt::Real)
    return q + 1/2*Δt*u
end

function _compute_index_set(m::Moreau)
    contact_threshold = 1e-4
    i=1
    for g in m.g
        if g <= contact_threshold
            push!(m.H, i)
        else
            try
                pop!(m.H, i)
            catch e
            end
        end
        i += 1
    end
    _update_force_matrix(m)
end

function _update_force_matrix(m::Moreau)
    n = length(m.qA)
    m.W = zeros(eltype(m.qA), n, 0)
    for i in m.H
        m.W = hcat(m.W, m.Wn[:,i])
    end
    m.Λ = zeros(eltype(m.qA), length(m.g))
end

function step(m::Moreau)
    if isempty(m.ϕ)
        step_unconstrained(m)
    else
        step_constrained(m)
    end
end

function step_unconstrained(m::Moreau)
    _compute_index_set(m)
    Minv = inv(m.M)

    model = Model(Mosek.Optimizer)
    set_optimizer_attribute(model, "QUIET", true)
    set_optimizer_attribute(model, "INTPNT_CO_TOL_DFEAS", 1e-7)
    @variable(model, q[1:length(m.qA)])
    @variable(model, u[1:length(m.uA)])
    @variable(model, λ[1:length(m.H)] >= 0)

    if length(m.W) != 0
        """ The following does not work because of positive semidefiniteness...
        ξ = m.W' * ( u + m.ε * m.uA )
        @constraint(model, ξ .>= 0)
        @objective(model, Min, dot(λ, ξ))
        """
        A = m.W' * Minv * m.W
        b = m.W' * Minv * m.h * m.Δt + (1+m.ε) * m.W' * m.uA
        @constraint(model, b .+ A*λ .>= 0)
        @constraint(model, m.M * (u - m.uA) .== m.W * λ + m.h * m.Δt)
        @objective(model, Min, dot(λ, b .+ A*λ))
    else
        @constraint(model, λ .== 0)
        @constraint(model, m.M * (u - m.uA) .== m.h * m.Δt)
        @objective(model, Min, 0)
    end
    @constraint(model, q .== _compute_mid_displacements(m.qM, u, m.Δt))
    optimize!(model)
    if !any( JuMP.termination_status(model) .== SUCCESS )
        @warn "termination status: $(JuMP.termination_status(model))."
    else
        m.qE = JuMP.value.(q)
        m.uE = JuMP.value.(u)
        if length(m.W) != 0
            i = 1
            for k in 1:length(m.g)
                if k in m.H
                    m.Λ[k] = JuMP.value.(λ)[i]
                    i += 1
                else
                    m.Λ[k] = 0.0
                end
            end
        else
            m.Λ = zeros(eltype(m.qA), length(m.g))
        end
    end
end

function step_constrained(m::Moreau)
    _compute_index_set(m)
    Minv = inv(m.M)
    hJ = -m.Jdot * m.uA - 2/m.ϵ*m.J*m.uA - 1/m.ϵ/m.ϵ*m.ϕ
    Mhat = inv( m.J * Minv * m.J' )

    model = Model(Mosek.Optimizer)
    set_optimizer_attribute(model, "QUIET", true)
    set_optimizer_attribute(model, "INTPNT_CO_TOL_DFEAS", 1e-7)
    @variable(model, q[1:length(m.qA)])
    @variable(model, u[1:length(m.uA)])
    @variable(model, λ[1:length(m.H)] >= 0)

    if length(m.W) != 0
        A = m.W' * Minv * ( I - m.J' * Mhat * m.J * Minv ) * m.W
        b = m.W' * Minv * ( m.h + m.J' * Mhat *(hJ - m.J * Minv * m.h) ) * m.Δt + (1+m.ε) * m.W' * m.uA
        μ = Mhat * ( hJ - m.J * Minv * m.h - 1/m.Δt * m.J * Minv * m.W * λ )
        @constraint(model, b .+ A*λ .>= 0)
        @constraint(model, m.M * (u - m.uA) .== m.W * λ + (m.h + m.J' * μ) * m.Δt)
        @objective(model, Min, dot(λ, b .+ A*λ))
    else
        μ = Mhat * ( hJ - m.J * Minv * m.h )
        @constraint(model, λ .== 0)
        @constraint(model, m.M * (u - m.uA) .== (m.h + m.J' * μ) * m.Δt)
        @objective(model, Min, 0)
    end
    @constraint(model, m.J*(u - m.uA) .== hJ * m.Δt)
    @constraint(model, q .== _compute_mid_displacements(m.qM, u, m.Δt))
    optimize!(model)
    if !any( JuMP.termination_status(model) .== SUCCESS )
        @warn "termination status: $(JuMP.termination_status(model))."
    else
        m.qE = JuMP.value.(q)
        m.uE = JuMP.value.(u)
        if length(m.W) != 0
            i = 1
            for k in 1:length(m.g)
                if k in m.H
                    m.Λ[k] = JuMP.value.(λ)[i]
                    i += 1
                else
                    m.Λ[k] = 0.0
                end
            end
            m.μ = Mhat * ( hJ - m.J * Minv * m.h - 1/m.Δt * m.J * Minv * m.W * m.Λ )
        else
            m.Λ = zeros(eltype(m.qA), length(m.g))
            m.μ = Mhat * ( hJ - m.J * Minv * m.h )
        end
    end
end

function set_state(m::Moreau, q::AbstractArray, u::AbstractArray)
    m.qA = q
    m.uA = u
    n = length(m.qA)
    m.qM = _compute_mid_displacements(m.qA, m.uA, m.Δt)
    m.M, m.h = m.dynamics(m.qM, m.uA)
    m.g, m.W = m.gap(m.qM, m.uA)
    if !isempty(m.ϕ)
        m.ϕ = m.hcon(m.qM)
        m.J = m.jac(m.qM)
        m.Jdot = m.jacdot(m.qM, m.uA)
    end
    _compute_index_set(m)
end

function set_state(m::Moreau, q::AbstractArray, u::AbstractArray, extra_state::AbstractArray)
    m.qA = q
    m.uA = u
    m.extra_state = extra_state
    n = length(m.qA)
    m.qM = _compute_mid_displacements(m.qA, m.uA, m.Δt)
    m.M, m.h = m.dynamics(m.qM, m.uA, m.extra_state)
    m.g, m.W = m.gap(m.qM, m.uA)
    if !isempty(m.ϕ)
        m.ϕ = m.hcon(m.qM)
        m.J = m.jac(m.qM)
        m.Jdot = m.jacdot(m.qM, m.uA)
    end
    _compute_index_set(m)
end


text:
using Cliffords
import LinearAlgebra

# helper to remove the first "+/-" char from Pauli string labels
snip(s::String) = s[nextind(s,1):end];

###############################################################################
## Useful tools for manipulation of matrix representations in a tomography ####
###############################################################################

"""
    unitary2choi(U::AbstractArray)

Returns the Choi super operator representation of a unitary.

# Arguments
- `U`: the unitary matrix representation of the process

# Returns
- `choi`: a choi matrix representation

# Examples
julia> choi = unitary2choi(str2unitary("1QX90p"))
4×4 Array{Complex{Float64},2}:
 0.25+0.0im    0.0+0.25im   0.0+0.25im  0.25+0.0im
  0.0-0.25im  0.25+0.0im   0.25+0.0im    0.0-0.25im
  0.0-0.25im  0.25+0.0im   0.25+0.0im    0.0-0.25im
 0.25+0.0im    0.0+0.25im   0.0+0.25im  0.25+0.0im
"""
function unitary2choi(U::AbstractArray)

    d=size(U,1);
    choi=zeros(ComplexF64,d^2,d^2);
    for ii in 1:d
        for jj in 1:d
            proj = zeros(ComplexF64,d,d);
            proj[ii,jj] = 1;
            choi = choi + kron(proj,U*proj*U')/d;
        end
    end
    return choi
end

"""
    choi2pauliMap(choi::AbstractArray)

Converts a Choi representation to a Pauli Map representation.

# Arguments
- `choi`: the matrix representation of the process

# Returns
- `pauli_map`: a Pauli representation of the choi matrix

# Examples
julia> choi = unitary2choi(str2unitary("1QX90p"))
julia> choi2pauliMap(choi)
4×4 Array{Float64,2}:
  1.0          0.0   0.0          0.0
  0.0          1.0   0.0          0.0
  0.0          0.0   3.33067e-16  1.0
 -5.55112e-17  0.0  -1.0          3.33067e-16
"""
function choi2pauliMap(choi::AbstractArray)

    #Dimension of superoperator
    d2 = size(choi,1);

    #Create the Pauli opearators for n qubits
    num_qubits = Int(log2(sqrt(d2)))
    pauliOps = map(x -> complex(x), vec(permutedims(allpaulis(num_qubits),
                                    reverse(collect(1:num_qubits)))));

    pauliMap = zeros(d2,d2);
    for ct1 in 1:d2
        p1 = pauliOps[ct1]
        for ct2 in 1:d2
            p2 = pauliOps[ct2]
            pauliMap[ct2,ct1] = real(LinearAlgebra.tr(choi*kron(transpose(p2), p1)));
        end
    end
    return pauliMap
end

"""
    unitary2pauli(U:AbstractArray)

Compute the Pauli representation of an arbitrary unitary

# Arguments
- `U`: the matrix representation of the unitary.

# Returns
- pauli_map: Matrix{ComplexF64}

# Examples
julia> unitary2pauli(str2unitary("1QY90p"))
4×4 Array{Float64,2}:
  1.0          0.0          0.0   0.0
  0.0          3.33067e-16  0.0  -1.0
  0.0          0.0          1.0   0.0
 -5.55112e-17  1.0          0.0   3.33067e-16
"""
function unitary2pauli(U::AbstractArray)
    return choi2pauliMap(unitary2choi(U));
end

"""
    str2unitary(strIn::String)

Converts a string description of a gate into a unitary matrix.  Supported
unitaries are: X_4I, IX_4, IX_2, InvCNOT12, CNOT12, CNOT21, 1QXp, 1QX22p,
IY_2, IY, ZXm90_2, 1QHad, XX, 1QId, IX_8, 1QT, YI, X90ZX90, CNOT12_y, X90X90,
1QX45p, ZXm90, 1QX90p, ZX90, 1QY90p, Y90Y90, 1QYp, X_8I, 1QY90m, X_2I, XZm90,
XI, XZ90, 1QX90m, 1QZ90, II, IX, Y_2I

# Arguments
- `strIn`: string id of the unitary

# Returns
- `U`: Matrix{ComplexF64} representing the unitary transformation

# Examples
julia> str2unitary("1QX90p")
2×2 Array{Complex{Float64},2}:
 0.707107+0.0im            0.0-0.707107im
      0.0-0.707107im  0.707107+0.0im
"""
function str2unitary(strIn::String)

    #Basic Pauli operators
    X = complex(Pauli(1));
    Y = complex(Pauli(3));
    Z = complex(Pauli(2));
    I = complex(Pauli(4));

    Uout = Dict("CNOT12" => Matrix{ComplexF64}([[1,0,0,0] [0,1,0,0] [0,0,0,1] [0,0,1,0]]),
                "InvCNOT12" => Matrix{ComplexF64}([[0,1,0,0] [1,0,0,0] [0,0,1,0] [0,0,0,1]]),
                "CNOT21" => Matrix{ComplexF64}([[1,0,0,0] [0,0,0,1] [0,0,1,0] [0,1,0,0]]),
                "1QId" => I,
                "1QX90p" => exp(-1im*pi*X/4),
                "1QX90m" => exp(1im*pi*X/4),
                "1QY90p" => exp(-1im*pi*Y/4),
                "1QY90m" => exp(1im*pi*Y/4),
                "1QXp" => exp(-1im*pi*X/2),
                "1QYp" => exp(-1im*pi*Y/2),
                "1QX45p" => exp(-1im*pi*X/8),
                "1QX22p" => exp(-1im*pi*X/16),
                "1QHad" => exp(-1im*(pi/2)*(1/sqrt(2))*(X+Z)),
                "1QZ90" => exp(-1im*(pi/4)*Z),
                "1QT" => exp(-1im*(pi/8)*Z),
                "II" => kron(I,I),
                "XI" => exp(-1im*kron(X,I)*pi/2),
                "IX" => exp(-1im*kron(I,X)*pi/2),
                "YI" => exp(-1im*kron(Y,I)*pi/2),
                "IY" => exp(-1im*kron(I,Y)*pi/2),
                "X_2I" => exp(-1im*kron(X,I)*pi/4),
                "IX_2" => exp(-1im*kron(I,X)*pi/4),
                "Y_2I" => exp(-1im*kron(Y,I)*pi/4),
                "IY_2" => exp(-1im*kron(I,Y)*pi/4),
                "XZ90" => exp(-1im*kron(X,Z)*pi/4),
                "ZX90" => exp(-1im*kron(Z,X)*pi/4),
                "XZm90" => exp(1im*kron(X,Z)*pi/4),
                "ZXm90" => exp(1im*kron(Z,X)*pi/4),
                "X_8I" => exp(-1im*kron(X,I)*pi/16),
                "X_4I" => exp(-1im*kron(X,I)*pi/8),
                "IX_8" => exp(-1im*kron(I,X)*pi/16),
                "IX_4" => exp(-1im*kron(I,X)*pi/8),
                "CNOT12_y" => Matrix{ComplexF64}([[1,0,0,0] [0,1,0,0] [0,0,-1im,0] [0,0,0,1im]]),
                "X90ZX90" => Matrix{ComplexF64}([[1,0,0,0] [0,0,0,-1im] [0,0,1,0] [0,-1im,0,0]]),
                "XX" => exp(-1im*kron(X,I)*pi/2)*exp(-1im*kron(I,X)*pi/2),
                "X90X90" => exp(-1im*kron(X,I)*pi/4)*exp(-1im*kron(I,X)*pi/4),
                "Y90Y90" => exp(-1im*kron(Y,I)*pi/4)*exp(-1im*kron(I,Y)*pi/4)
            );
    return Uout[strIn]
end

"""
    choi2chi(choi::AbstractArray)

Converts a Choi superoperator to a chi represtation in Pauli basis.

# Arguments
- `choi`: the matrix representation of the process

# Returns
- `chi`: chi matrix representation of the process

# Examples
julia> using RandomQuantum
julia> U = rand(RandomClosedEvolution(2,2))
julia> choi = unitary2choi(U)
julia> choi2chi(choi)
"""
function choi2chi(choi::AbstractArray)

    # Some dimensions
    d2 = size(choi,1);
    d = Int(sqrt(d2));

    # Get the Kraus operators from the eigen decomposition of the Choi
    vals, vecs = LinearAlgebra.eigen(choi);

    # note the ordering diff between Matlab
    vals = LinearAlgebra.diagm(reverse(vals));
    vecs = reverse(vecs, dims=2);

    chi = zeros(ComplexF64,d2,d2);

    #Create the Pauli opearators for n qubits
    num_qubits = Int(log2(sqrt(d2)));
    pauliOps = map(x -> complex(x), vec(permutedims(allpaulis(num_qubits),
                                    reverse(collect(1:num_qubits)))));

    # Transform from the Krauss basis to the Pauli basis
    for kraussct in 1:size(vals,1)
        tmpKrauss = reshape(vecs[:,kraussct],d,d)*sqrt(d); # Krauss operator should have norm of d
        for paulict1 in 1:d2
            pauliLeft = LinearAlgebra.tr(pauliOps[paulict1]*tmpKrauss)/d;
            for paulict2 in 1:d2
                pauliRight = LinearAlgebra.tr(pauliOps[paulict2]*transpose(tmpKrauss))/d;
                chi[paulict1, paulict2] = chi[paulict1, paulict2] + vals[kraussct]*pauliLeft*pauliRight;
            end
        end
    end
    return chi
end

###############################################################################
## Fidelity and plotting tools for tomography #################################
###############################################################################
"""
    getProcessFidelity(choi::AbstractArray, idealProcess::String)

Compute the gate and process fidelity for a give Choi matrix and a given ideal
process matrix.

# Arguments
- `choi`: the matrix representation of the process
- `idealProcess`: string representing the ideal process in str2unitary()

# Returns
- `processFidelity`: -> real(tr(choi*choiIdeal))
- `gateFidelity`: -> (d*processFidelity +1)/(d + 1)

# Examples
julia> getProcessFidelity(choi, "XI")
"""
function getProcessFidelity(choi::AbstractArray, idealProcess::String)
    nbrQubits = log2(sqrt(size(choi, 1)))
    # Calculate the overlaps with the ideal gate
    if typeof(idealProcess) == String
        unitaryIdeal = str2unitary(idealProcess);
    else
        # assume we are passed a matrix otherwise
        unitaryIdeal = idealProcess;
    end
    choiIdeal = unitary2choi(unitaryIdeal);

    processFidelity = real(LinearAlgebra.tr(choi*choiIdeal));
    gateFidelity = (2^nbrQubits*processFidelity+1)/(2^nbrQubits+1);

    return processFidelity, gateFidelity
end

"""
    getProcessFidelity(choi::AbstractArray, idealProcess::AbstractMatrix)

Compute the gate and process fidelity for a give Choi matrix and a given ideal
unitary.

# Arguments
- `choi`: the matrix representation of the process
- `idealProcess`: matrix representing the ideal process unitary

# Returns
- `processFidelity`: -> real(tr(choi*choiIdeal))
- `gateFidelity`: -> (d*processFidelity +1)/(d + 1)

# Examples
julia> getProcessFidelity(choi, U)
"""
function getProcessFidelity(choi::AbstractArray, idealProcess::AbstractMatrix)
    nbrQubits = log2(sqrt(size(choi, 1)))
    # Calculate the overlaps with the ideal gate
    choiIdeal = unitary2choi(idealProcess);

    processFidelity = real(LinearAlgebra.tr(choi*choiIdeal));
    gateFidelity = (2^nbrQubits*processFidelity+1)/(2^nbrQubits+1);

    return processFidelity, gateFidelity
end

"""
    PlotProcess(U::AbstractMatrix; p_map::String="pauli", save_fig::Bool=false)

Plot the Pauli transfer matrix of an operator U.  Only one and two qubit
operators are supported.

# Arguments
- `U`: the matrix representation of the unitary you want to plot.
- `p_map`: which process map to use {"pauli", "chi"}.  "pauli" will plot the
           Pauli transfer matrix and "chi" will plot the chi matrix.
- `save_fig`: set to true if you want a .png of the figure saved in the current
              directory.

# Returns
- void: The function outputs a plot

# Examples
julia> PauliPlot(str2unitary("1QY90p"), save_fig=true)
"""
function PlotProcess(U::AbstractMatrix;
                     p_map::String="pauli",
                     save_fig::Bool=false)

    num_qubits = Int(log2(size(U,1)))
    op_num = 0:(4^num_qubits - 1)
    # This code snippet works for n-qubit unitaries, it's just not that
    # verbose...
    labels = map(p -> snip(string(p)),
                 vec(permutedims(allpaulis(num_qubits),
                 reverse(collect(1:num_qubits)))))
    # a more verbose and brittle option would be something like:
    # if num_qubits == 1
    #     labels = ["I", "X", "Y", "Z"]
    # elseif num_qubits == 2
    #     labels = ["II","IX","IY","IZ",
    #               "XI","XX","XY","XZ",
    #               "YI","YX","YY","YZ",
    #               "ZI","ZX","ZY","ZZ"]
    # else
    #     labels = []
    # end
    fig = figure()
    ax = fig.add_subplot(111)
    if p_map == "pauli"
        cax = ax.matshow(unitary2pauli(U), cmap="bwr")
        ax.set_title("Pauli transfer matrix")
    elseif p_map == "chi"
        cax = ax.matshow(real(choi2chi(unitary2choi(U))), cmap="bwr")
        ax.set_title(L"$\chi$ process matrix")
    end

    fig.colorbar(cax)

    ax.set_xticklabels(labels)
    ax.set_yticklabels(labels)
    ax.set_xticks(op_num)
    ax.set_yticks(op_num)
    ax.xaxis.set_ticks_position("bottom")
    ax.grid(color="grey", linestyle="--", alpha=0.1)
    ax.spines["right"].set_visible(false)
    ax.spines["top"].set_visible(false)
    ax.spines["left"].set_visible(false)
    ax.spines["bottom"].set_visible(false)
    ax.set_xlabel("Input Pauli Operator")
    ax.set_ylabel("Output Pauli Operator")

    if save_fig
        fig.savefig("ProcessPlot.png")
        println("Please rename the saved figure!")
    end
    fig.show()
end

"""
    PauliToCounts(Z::AbstractArray,nbrQubits)

Convert a set of 2^nbrQubits-1 of SigmaZ_j and correlator SigmaZ_iSigmaZ_j... expectation values into state populations (expectation value of projection operators)

# Arguments
- `Z`: array of SigmaZ expectation values. Z[k] is the expectation value of a product of I and Z operators corresponding to the binary value of k.
        So for example if k=0110 , Z[k]=<IZZI> . Z needs to have size of 2^nbrQubits-1
- `nbrQubits`: number of Qubits.

# Returns
- Array of projector expectation values

# Examples
julia> PauliToCounts([1,1,1],2)
"""
function PauliToCounts(Z::AbstractArray,nbrQubits)
  P = ones(1,2^nbrQubits)

  for k in 0:2^nbrQubits-1
       kbin = digits!(zeros(Int64,nbrQubits),k,base=2)
       kbin = -2*kbin .+ 1

       for h in 1:2^nbrQubits-1
           hbin = digits!(zeros(Int64,nbrQubits),h,base=2)
           P[k+1] += prod(sign.(kbin.*hbin .+ 0.5))*Z[h]
       end

   end
   return P/2^nbrQubits
end


text:
using MPI
using CUDA
using TickTock

## lr=parse(Int64,ENV["OMPI_COMM_WORLD_LOCAL_RANK"])
## lr=parse(Int64,ENV["SLURM_LOCALID"])
## CUDA.device!(lr)
## dev=device()
## pool=CuMemoryPool(dev,alloc_type=CUDA.CU_MEM_ALLOCATION_TYPE_PINNED,handle_type=CUDA.CU_MEM_HANDLE_TYPE_POSIX_FILE_DESCRIPTOR)
## pool=CuMemoryPool(dev,alloc_type=CUDA.CU_MEM_ALLOCATION_TYPE_PINNED,handle_type=CUDA.CU_MEM_HANDLE_TYPE_NONE)
## memory_pool!(dev,pool)

MPI.Init()

      comm = MPI.COMM_WORLD
local_rank = MPI.Comm_rank(comm)
         R = MPI.Comm_size(comm)

if R != 2
 exit()
end

# See: https://github.com/CliMA/ClimateMachine.jl/blob/2def0d653fcef3d069e816844fa4ec5d745c84d1/src/Driver/Driver.jl#L40-L43
local_comm = MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED, MPI.Comm_rank(comm))
CUDA.device!(MPI.Comm_rank(local_comm) % length(devices()))
dev=device()
pool=CuMemoryPool(dev,alloc_type=CUDA.CU_MEM_ALLOCATION_TYPE_PINNED,handle_type=CUDA.CU_MEM_HANDLE_TYPE_POSIX_FILE_DESCRIPTOR)
memory_pool!(dev,pool)

nel=10000000
nm =50

sendArr = CUDA.zeros(Float64, nel)
recvArr = CUDA.zeros(Float64, nel)
fill!(sendArr, local_rank)

if local_rank == 0
 recv_from_rank = 1
 recv_from_tag  = 11
 send_to_rank   = 1
 sender_tag     = 10
end

if local_rank == 1
 recv_from_rank = 0
 recv_from_tag  = 10
 send_to_rank   = 0
 sender_tag     = 11
end

for nmesg = 1:4
 rreq  = MPI.Irecv!(recvArr, recv_from_rank, recv_from_tag, comm)
 sreq  = MPI.Isend( sendArr, send_to_rank  , sender_tag,    comm)
 stats = MPI.Waitall!([rreq, sreq])
 tick();tok()
end

function sr(nm)
 for nmesg = 1:nm
  rreq  = MPI.Irecv!(recvArr, recv_from_rank, recv_from_tag, comm)
  sreq  = MPI.Isend( sendArr, send_to_rank  , sender_tag,    comm)
  stats = MPI.Waitall!([rreq, sreq])
 end
end

sr(1)
sr(1)
sr(1)

tick()
sr(nm)
elap=tok()

println("elap = ",elap, " BW (GB/s) ", 2*nel*nm*8/elap/1e9)

MPI.Barrier(comm)


text:
using JLD2, DataFrames, Serialization
#, Distributions, Serialization, Dates, Statistics, DataFrames
import GEAS:breeding, sim_QTL, summarize

function vnqtl(nqtl)
    par = begin
        Parameters = Dict(
            :nSire=> 100,
            :nDam => 200,   # male:female = 1:2
            :nSib => 40,
            :nC7e => 30,    # number of sib to be challenged
            :nG8n => 10,    # number of generations
            :nQTL => [nqtl, nqtl],
            :h²   => [.5, .5],
            :p8e  => .5,    # percentage of death after challenge
            :e19e => 1.,    # edit_successsful_rate
            :w4t  => 1,     # weight on the binary trait EBV
            :nK3n => 10,    # number of known QTL on binary trait
            :dd   => 0.01,  # A small value to be added to diagonals
            #:log  => slog
        )
        (; Parameters...)       # named tuple.  contents as above
    end
end

function append_rst(rst, ped, nqtl, method, i)
    df = summarize(ped)
    ng = nrow(df)
    df.nqtl = ones(Int, ng) .* nqtl
    df.method = ones(Int, ng) .* method
    df.repeat = ones(Int, ng) .* i
    append!(rst, df)
end

function t_2021_08_16(file)
    @load "dat/run/base.jld2" base
    nrpt = 30
    rst = DataFrame()
    for nqtl in [100, 500]
        par = vnqtl(nqtl)
        for i in 1:nrpt
            @info "==== Repeat $i of $nrpt ===="
            qtl = sim_QTL(base, nqtl, nqtl)

            ped, snp = breeding(base, par, qtl)
            append_rst(rst, ped, nqtl, 1, i)

            ped, snp = breeding(base, par, qtl, fixed=true)
            append_rst(rst, ped, nqtl, 2, i)

            ped, snp = breeding(base, par, qtl, edit=true)
            append_rst(rst, ped, nqtl, 3, i)
        end
    end
    serialize(file, rst)
end


text:
"""Trial Balance?
book.jl
 WIP

"""

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

baremodule Base

using Core.Intrinsics
ccall(:jl_set_istopmod, Void, (Bool,), true)
function include(path::AbstractString)
    local result
    if INCLUDE_STATE === 1
        result = Core.include(path)
    elseif INCLUDE_STATE === 2
        result = _include(path)
    elseif INCLUDE_STATE === 3
        result = include_from_node1(path)
    end
    result
end
INCLUDE_STATE = 1 # include = Core.include

include("coreio.jl")

eval(x) = Core.eval(Base, x)
eval(m, x) = Core.eval(m, x)
(::Type{T}){T}(arg) = convert(T, arg)::T # Hidden from the REPL.
(::Type{VecElement{T}}){T}(arg) = VecElement{T}(convert(T, arg))
convert{T<:VecElement}(::Type{T}, arg) = T(arg)
convert{T<:VecElement}(::Type{T}, arg::T) = arg

# init core docsystem
import Core: @doc, @__doc__, @doc_str
if isdefined(Core, :Inference)
    import Core.Inference.CoreDocs
    Core.atdoc!(CoreDocs.docm)
end

include("exports.jl")

if false
    # simple print definitions for debugging. enable these if something
    # goes wrong during bootstrap before printing code is available.
    # otherwise, they just just eventually get (noisily) overwritten later
    global show, print, println
    show(io::IO, x::ANY) = Core.show(io, x)
    print(io::IO, a::ANY...) = Core.print(io, a...)
    println(io::IO, x::ANY...) = Core.println(io, x...)
end

## Load essential files and libraries
include("essentials.jl")
include("ctypes.jl")
include("base.jl")
include("generator.jl")
include("reflection.jl")
include("options.jl")

# core operations & types
include("promotion.jl")
include("tuple.jl")
include("pair.jl")
include("traits.jl")
include("range.jl")
include("twiceprecision.jl")
include("expr.jl")
include("error.jl")

# core numeric operations & types
include("bool.jl")
include("number.jl")
include("int.jl")
include("operators.jl")
include("pointer.jl")
include("refpointer.jl")
include("checked.jl")
importall .Checked

# vararg Symbol constructor
Symbol(x...) = Symbol(string(x...))

# Define the broadcast function, which is mostly implemented in
# broadcast.jl, so that we can overload broadcast methods for
# specific array types etc.
#  --Here, just define fallback routines for broadcasting with no arguments
broadcast(f) = f()
broadcast!(f, X::AbstractArray) = (@inbounds for I in eachindex(X); X[I] = f(); end; X)

# array structures
include("indices.jl")
include("array.jl")
include("abstractarray.jl")
include("subarray.jl")

# Array convenience converting constructors
(::Type{Array{T}}){T}(m::Integer) = Array{T,1}(Int(m))
(::Type{Array{T}}){T}(m::Integer, n::Integer) = Array{T,2}(Int(m), Int(n))
(::Type{Array{T}}){T}(m::Integer, n::Integer, o::Integer) = Array{T,3}(Int(m), Int(n), Int(o))
(::Type{Array{T}}){T}(d::Integer...) = Array{T}(convert(Tuple{Vararg{Int}}, d))

(::Type{Vector})() = Array{Any,1}(0)
(::Type{Vector{T}}){T}(m::Integer) = Array{T,1}(Int(m))
(::Type{Vector})(m::Integer) = Array{Any,1}(Int(m))
(::Type{Matrix{T}}){T}(m::Integer, n::Integer) = Matrix{T}(Int(m), Int(n))
(::Type{Matrix})(m::Integer, n::Integer) = Matrix{Any}(Int(m), Int(n))

# numeric operations
include("hashing.jl")
include("rounding.jl")
importall .Rounding
include("float.jl")
include("complex.jl")
include("rational.jl")
include("multinverses.jl")
using .MultiplicativeInverses
include("abstractarraymath.jl")
include("arraymath.jl")

# define MIME"foo/bar" early so that we can overload 3-arg show
struct MIME{mime} end
macro MIME_str(s)
    :(MIME{$(Expr(:quote, Symbol(s)))})
end

include("char.jl")
include("strings/string.jl")

# SIMD loops
include("simdloop.jl")
importall .SimdLoop

# map-reduce operators
include("reduce.jl")

## core structures
include("reshapedarray.jl")
include("bitarray.jl")
include("intset.jl")
include("associative.jl")
include("dict.jl")
include("set.jl")
include("iterators.jl")
using .Iterators: zip, enumerate
using .Iterators: Flatten, product  # for generators

# Definition of StridedArray
StridedReshapedArray{T,N,A<:DenseArray} = ReshapedArray{T,N,A}
StridedArray{T,N,A<:Union{DenseArray,StridedReshapedArray},
    I<:Tuple{Vararg{Union{RangeIndex, AbstractCartesianIndex}}}} =
    Union{DenseArray{T,N}, SubArray{T,N,A,I}, StridedReshapedArray{T,N}}
StridedVector{T,A<:Union{DenseArray,StridedReshapedArray},
    I<:Tuple{Vararg{Union{RangeIndex, AbstractCartesianIndex}}}} =
    Union{DenseArray{T,1}, SubArray{T,1,A,I}, StridedReshapedArray{T,1}}
StridedMatrix{T,A<:Union{DenseArray,StridedReshapedArray},
    I<:Tuple{Vararg{Union{RangeIndex, AbstractCartesianIndex}}}} =
    Union{DenseArray{T,2}, SubArray{T,2,A,I}, StridedReshapedArray{T,2}}
StridedVecOrMat{T} = Union{StridedVector{T}, StridedMatrix{T}}

# For OS specific stuff
include(string((length(Core.ARGS)>=2 ? Core.ARGS[2] : ""), "build_h.jl"))     # include($BUILDROOT/base/build_h.jl)
include(string((length(Core.ARGS)>=2 ? Core.ARGS[2] : ""), "version_git.jl")) # include($BUILDROOT/base/version_git.jl)

include("osutils.jl")
include("c.jl")
include("sysinfo.jl")

if !isdefined(Core, :Inference)
    include("docs/core.jl")
    Core.atdoc!(CoreDocs.docm)
end

# Core I/O
include("io.jl")
include("iostream.jl")
include("iobuffer.jl")

# strings & printing
include("intfuncs.jl")
include("strings/strings.jl")
include("parse.jl")
include("shell.jl")
include("regex.jl")
include("show.jl")

# multidimensional arrays
include("cartesian.jl")
using .Cartesian
include("multidimensional.jl")
include("permuteddimsarray.jl")
using .PermutedDimsArrays

# nullable types
include("nullable.jl")

include("broadcast.jl")
importall .Broadcast

# base64 conversions (need broadcast)
include("base64.jl")
importall .Base64

# version
include("version.jl")

# system & environment
include("libc.jl")
using .Libc: getpid, gethostname, time
include("libdl.jl")
using .Libdl: DL_LOAD_PATH
include("env.jl")

# Scheduling
include("libuv.jl")
include("event.jl")
include("task.jl")
include("lock.jl")
include("threads.jl")
include("weakkeydict.jl")

# I/O
include("stream.jl")
include("socket.jl")
include("filesystem.jl")
importall .Filesystem
include("process.jl")
include("multimedia.jl")
importall .Multimedia
include("grisu/grisu.jl")
import .Grisu.print_shortest
include("methodshow.jl")

# core math functions
include("floatfuncs.jl")
include("math.jl")
importall .Math
const (√)=sqrt
const (∛)=cbrt

let SOURCE_PATH = ""
    global function _include(path)
        prev = SOURCE_PATH
        path = joinpath(dirname(prev),path)
        SOURCE_PATH = path
        Core.include(path)
        SOURCE_PATH = prev
    end
end
INCLUDE_STATE = 2 # include = _include (from lines above)

# reduction along dims
include("reducedim.jl")  # macros in this file relies on string.jl

# basic data structures
include("ordering.jl")
importall .Order

# Combinatorics
include("sort.jl")
importall .Sort

function deepcopy_internal end

# BigInts and BigFloats
include("gmp.jl")
importall .GMP
include("mpfr.jl")
importall .MPFR
big(n::Integer) = convert(BigInt,n)
big(x::AbstractFloat) = convert(BigFloat,x)
big(q::Rational) = big(numerator(q))//big(denominator(q))

include("combinatorics.jl")

# more hashing definitions
include("hashing2.jl")

# random number generation
include("dSFMT.jl")
include("random.jl")
importall .Random

# (s)printf macros
include("printf.jl")
importall .Printf

# metaprogramming
include("meta.jl")

# enums
include("Enums.jl")
importall .Enums

# concurrency and parallelism
include("serialize.jl")
importall .Serializer
include("channels.jl")

# memory-mapped and shared arrays
include("mmap.jl")
import .Mmap

# utilities - timing, help, edit
include("datafmt.jl")
importall .DataFmt
include("deepcopy.jl")
include("interactiveutil.jl")
include("summarysize.jl")
include("replutil.jl")
include("test.jl")
include("i18n.jl")
using .I18n

# frontend
include("initdefs.jl")
include("Terminals.jl")
include("LineEdit.jl")
include("REPLCompletions.jl")
include("REPL.jl")
include("client.jl")

# Stack frames and traces
include("stacktraces.jl")
importall .StackTraces

# misc useful functions & macros
include("util.jl")

# dense linear algebra
include("linalg/linalg.jl")
importall .LinAlg
const ⋅ = dot
const × = cross

# statistics
include("statistics.jl")

# irrational mathematical constants
include("irrationals.jl")

# signal processing
include("dft.jl")
importall .DFT
include("dsp.jl")
importall .DSP

# Fast math
include("fastmath.jl")
importall .FastMath

# libgit2 support
include("libgit2/libgit2.jl")

# package manager
include("pkg/pkg.jl")

# profiler
include("profile.jl")
importall .Profile

# dates
include("dates/Dates.jl")
import .Dates: Date, DateTime, DateFormat, @dateformat_str, now

# sparse matrices, vectors, and sparse linear algebra
include("sparse/sparse.jl")
importall .SparseArrays

include("asyncmap.jl")

include("distributed/Distributed.jl")
importall .Distributed
include("sharedarray.jl")

# code loading
include("loading.jl")

# worker threads
include("threadcall.jl")

# deprecated functions
include("deprecated.jl")

# Some basic documentation
include("docs/helpdb.jl")
include("docs/basedocs.jl")

# Documentation -- should always be included last in sysimg.
include("markdown/Markdown.jl")
include("docs/Docs.jl")
using .Docs, .Markdown
isdefined(Core, :Inference) && Docs.loaddocs(Core.Inference.CoreDocs.DOCS)

function __init__()
    # Base library init
    reinit_stdio()
    Multimedia.reinit_displays() # since Multimedia.displays uses STDOUT as fallback
    early_init()
    init_load_path()
    Distributed.init_parallel()
    init_threadcall()
end

INCLUDE_STATE = 3 # include = include_from_node1
include("precompile.jl")

end # baremodule Base

using Base
importall Base.Operators

Base.isfile("userimg.jl") && Base.include("userimg.jl")

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
using DandelionWebSockets
using Logging

# As part of improving the tests in this package, I feel the need to differentiate between the
# different types of tests I have here.
#
# For want of better terms, I'll categorize the tests into one of the following:
#
# - Small
# - Medium
# - Large
#
# Below is a table of the different aspects of a test that are used to categorize them. The most
# important aspects are speed, scope, and determinism. It's important that most tests are
# entirely deterministic, that is, the test will either succeed 100% of the time, or fail 100% of
# the time. It should not succeed 99% of the time, and fail 1% of the time. Because of this,
# testing concurrent code is forbidden in small tests. Medium tests may be concurrent, but care
# should be taken to make the test deterministic. Large tests may be written against code that fails
# due to external factors, but care should be taken to minimize the influence of those external
# factors.
#
# ====================================================================
# Aspect      | Small          | Medium                 | Large      |
# ====================================================================
# Speed       | < 10s          | < 10s                  | < 1 minute |
# --------------------------------------------------------------------
# Concurrency | None           | Yes, but deterministic | Yes        |
# --------------------------------------------------------------------
# Scope       | Type or method | Type or method         | Package    |
# ====================================================================
#
# NOTE: Speed applies to _all_ tests in that category, not a single test case.

disable_logging(Logging.AboveMaxLevel)

# TODO Remove dependency on FactCheck
#include("mock_test.jl")

include("utilities.jl")

# include("network_test.jl")
# TODO Remove dependency on FactCheck
#include("ping_test.jl")
# TODO Remove dependency on FactCheck
#include("handshake_test.jl")
# TODO Remove dependency on FactCheck
#include("reconnect_test.jl")
# TODO Remove dependency on FactCheck
#include("integration_test.jl")

include("runsmalltests.jl")
include("runmediumtests.jl")


text:
using SEC, Test
import SEC.Emulate: mean, cov, GaussianProcess, condition

@testset "GaussianProcess unit test" begin
    μ(x) = 0.0
    k(x,y) = exp(-(x-y)^2/2)
    gp = GaussianProcess(mean = μ, covariance = k)
    @test μ === mean(gp)
    @test k === cov(gp)
end

@testset "1D GPR mean and cov test" begin
    μ(x) = 0.0
    k(x,y) = exp(-(x-y)^2/2)
    gp = GaussianProcess(mean = μ, covariance = k)
    N = 10
    X = [i for i in 1:N]
    μX = mean(gp, X)
    σX = cov(gp, X)
    @test all(μX .== zeros(Float64, N))
    @test σX[1,1] ≈ 1.0
    @test σX[1, N] ≈ k(1, N)
end

@testset "1D GPR conditional mean and cov test" begin
    μ(x) = 0.0
    k(x,y) = exp(-(x-y)^2/2)
    gp = GaussianProcess(mean = μ, covariance = k)
    N = 10
    X = [i for i in 1:N]
    Y = [2*i for i in 1:N]
    data = (; X, Y)
    cgp = condition(gp, data)
    cμX =  mean(cgp, X)
    tolerance = 1e-12
    @test all(abs.(cμX .- Y) .≤ tolerance)
    covXX =  [cgp.covariance(x,x) for x in X]
    @test all(abs.(covXX) .≤ tolerance)
end

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
using DataFrames
using CSV
using Statistics
using Base.Threads: @threads
using StatsBase: sample

using EvoTrees
using EvoTrees: get_gain, get_max_gain, update_grads!, grow_tree, grow_gbtree, SplitInfo, Tree, TrainNode, TreeNode, Params, predict, predict!, find_split!, SplitTrack, update_track!, sigmoid

# prepare a dataset
# prepare a dataset
features = rand(UInt8, 100000, 300)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X,1))
𝑗 = collect(1:size(X,2))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# idx
X_perm = zeros(Int, size(X))
@threads for feat in 1:size(X, 2)
    X_perm[:, feat] = sortperm(X[:, feat]) # returns gain value and idx split
    # idx[:, feat] = sortperm(view(X, :, feat)) # returns gain value and idx split
end

# set parameters
nrounds = 1
λ = 1.0
γ = 1e-15
η = 0.5
max_depth = 5
min_weight = 5.0
rowsample = 1.0
colsample = 1.0

# params1 = Params(nrounds, λ, γ, η, max_depth, min_weight, :linear)
params1 = Params(:linear, 1, λ, γ, 1.0, 2, min_weight, rowsample, colsample)

# initial info
δ, δ² = zeros(size(X, 1)), zeros(size(X, 1))
𝑤 = ones(size(X, 1))
pred = zeros(size(Y, 1))
# @time update_grads!(Val{params1.loss}(), pred, Y, δ, δ²)
update_grads!(Val{params1.loss}(), pred, Y, δ, δ², 𝑤)
∑δ, ∑δ², ∑𝑤 = sum(δ), sum(δ²), sum(𝑤)
gain = get_gain(∑δ, ∑δ², ∑𝑤, params1.λ)

# initialize train_nodes
train_nodes = Vector{TrainNode{Float64, Array{Int64,1}, Array{Int64, 1}, Int}}(undef, 2^params1.max_depth-1)
for feat in 1:2^params1.max_depth-1
    train_nodes[feat] = TrainNode(0, -Inf, -Inf, -Inf, -Inf, [0], [0])
end
# initializde node splits info and tracks - colsample size (𝑗)
splits = Vector{SplitInfo{Float64, Int}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    splits[feat] = SplitInfo{Float64, Int}(-Inf, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -Inf, -Inf, 0, 0, 0.0)
end
tracks = Vector{SplitTrack{Float64}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    tracks[feat] = SplitTrack{Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -Inf, -Inf, -Inf)
end

# placeholder for sort perm
perm_ini = zeros(Int, size(X))
for feat in 1:size(𝑗, 1)
    perm_ini[:, feat] = 𝑖
end

train_nodes[1] = TrainNode(1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)

function sortperm_test1(node, splits, perm_ini, params)
    node_size = size(node.𝑖, 1)
    @threads for feat in node.𝑗
        sortperm!(view(perm_ini, 1:node_size, feat), view(X, node.𝑖, feat), alg = QuickSort, initialized = true)
        find_split!(view(X, view(node.𝑖, view(perm_ini, 1:node_size, feat)), feat), view(δ, view(node.𝑖, view(perm_ini, 1:node_size, feat))) , view(δ², view(node.𝑖, view(perm_ini, 1:node_size, feat))), view(𝑤, view(node.𝑖, view(perm_ini, 1:node_size, feat))), node.∑δ, node.∑δ², node.∑𝑤, params.λ, splits[feat], tracks[feat])
    end
end

function sortperm_test2(node, splits, perm_ini, params)
    node_size = size(node.𝑖, 1)
    @threads for feat in node.𝑗
        perm_ini[1:node_size, feat] = sortperm(X[node.𝑖, feat], alg = QuickSort)
        find_split!(view(X, view(node.𝑖, view(perm_ini, 1:node_size, feat)), feat), view(δ, view(node.𝑖, view(perm_ini, 1:node_size, feat))) , view(δ², view(node.𝑖, view(perm_ini, 1:node_size, feat))), view(𝑤, view(node.𝑖, view(perm_ini, 1:node_size, feat))), node.∑δ, node.∑δ², node.∑𝑤, params.λ, splits[feat], tracks[feat])
    end
end

function sortperm_radix_test1(node, splits, perm_ini, params)
    node_size = size(node.𝑖, 1)
    @threads for feat in node.𝑗
        perm_ini[1:node_size, feat] = faster_sortperm_radix(X[node.𝑖, feat])
        find_split!(view(X, view(node.𝑖, view(perm_ini, 1:node_size, feat)), feat), view(δ, view(node.𝑖, view(perm_ini, 1:node_size, feat))) , view(δ², view(node.𝑖, view(perm_ini, 1:node_size, feat))), view(𝑤, view(node.𝑖, view(perm_ini, 1:node_size, feat))), node.∑δ, node.∑δ², node.∑𝑤, params.λ, splits[feat], tracks[feat])
    end
end

@time sortperm_test1(train_nodes[1], splits, perm_ini, params1)
@time sortperm_test2(train_nodes[1], splits, perm_ini, params1)
@time sortperm_radix_test1(train_nodes[1], splits, perm_ini, params1)

sizeof(X)/1024^2

using SortingAlgorithms
function faster_sortperm_radix(v)
  ai = [Pair(i, a) for (i,a) in enumerate(v)]
  sort!(ai, by=x->x.second, alg=RadixSort)
  [a.first for a in ai]
end

function faster_sortperm(v)
  ai = [Pair(i, a) for (i,a) in enumerate(v)]
  sort!(ai, by=x->x.second)
  [a.first for a in ai]
end

function test_faster_sortperm(node, splits, perm_ini, params)
    node_size = size(node.𝑖, 1)
    @threads for feat in node.𝑗
        perm_ini[1:node_size, feat] = faster_sortperm(X[node.𝑖, feat])
        find_split!(view(X, view(node.𝑖, view(perm_ini, 1:node_size, feat)), feat), view(δ, view(node.𝑖, view(perm_ini, 1:node_size, feat))) , view(δ², view(node.𝑖, view(perm_ini, 1:node_size, feat))), view(𝑤, view(node.𝑖, view(perm_ini, 1:node_size, feat))), node.∑δ, node.∑δ², node.∑𝑤, params.λ, splits[feat], tracks[feat])

    end
end

function test_faster_sortperm_radix(node, splits, perm_ini, params)
    node_size = size(node.𝑖, 1)
    @threads for feat in node.𝑗
        perm_ini[1:node_size, feat] = faster_sortperm_radix(X[node.𝑖, feat])
        find_split!(view(X, view(node.𝑖, view(perm_ini, 1:node_size, feat)), feat), view(δ, view(node.𝑖, view(perm_ini, 1:node_size, feat))) , view(δ², view(node.𝑖, view(perm_ini, 1:node_size, feat))), view(𝑤, view(node.𝑖, view(perm_ini, 1:node_size, feat))), node.∑δ, node.∑δ², node.∑𝑤, params.λ, splits[feat], tracks[feat])

    end
end

@time test_faster_sortperm(train_nodes[1], splits, perm_ini, params1)
@time test_faster_sortperm_radix(train_nodes[1], splits, perm_ini, params1)

x1 = rand(Int, 100000)

@time sortperm(x1)
@time faster_sortperm(x1)
@time faster_sortperm_radix(x1)


text:
using InfixFunctions
using Test

# Test short and long function definitions:
add(x, y) = x + y
function sub(x, y)
    return x - y
end

@infix add
@infix sub

x = 5

@test x |add| x == 10
@test x |sub| x == 0

# Test multiple dispatch and type parameters:
function foo(x::T, y::T) where {T<:Int}
    return Complex(x + y)
end

foo(x::T, y::S) where {T<:Int, S<:Float64} = x - y

@infix foo

@test 3 |foo| 5 == 8 + 0im
@test 3 |foo| 5.0 == -2.0

# Test infixifying normal functions:
@infix div
@infix div  # should work again

@test 10 |div| 5 == 2

# Test return types:
function bar(x::T, y::T)::T where {T<:Int}
    return x + y
end

(bar(x::T, y::S)::S) where {T<:Int, S<:Float64} = x - y

@infix bar

@test (5 |bar| 5)::Int == 10
@test (5 |bar| 5.0)::Float64 == 0.0
@test (10 |div| 5)::Int == 2

# Test normal function call syntax:
@test x |add| x == add(x, x)
@test x |sub| x == sub(x, x)
@test 3 |foo| 5 == foo(3, 5)
@test 3 |foo| 5.0 == foo(3, 5.0)
@test 10 |div| 5 == div(10, 5)
@test (5 |bar| 5)::Int == bar(5, 5)::Int
@test (5 |bar| 5.0)::Float64 == bar(5, 5.0)::Float64
@test (10 |div| 5)::Int == div(10, 5)::Int


text:
using DiffEqFlux, Flux

Zygote.@nograd diagm
include("eig.jl")

""" LDS condition for CTRNN1
        Enforce Ω(I - τAW) + (I - τW'A')Ω ≻ 0,
        parametrize sqrt(Ω) """
struct ISSCTRNN1{P,Q1,Q2,Q3,F} <: Modifier
    name::String
    param_size::Int
    state_dim::Int
    hidden_dim::Int
    initial_params::P
    get_Ω::Q1
    get_sqrtΩ::Q2
    get_sqrtΩinv::Q3
    ϵ::F
    lb::F
    logτ_idx::UnitRange{Int}
    W_idx::UnitRange{Int}
    A_idx::UnitRange{Int}
    enable::Bool
    scale_W::Bool
    function ISSCTRNN1(
        state_dim,
        hidden_dim;
        enable = true,
        is_fixed = false,
        scale_W = false,
        initb = Flux.zeros,
        ϵ = 1e-3,
        shift = -4.0,
    )
        # initial parameters
        initial_params() = vec(initb(hidden_dim))

        # function to fetch parameters: if fixed, ignore parameters
        if is_fixed
            Ω_init = ones(hidden_dim)
            get_Ω = (p) -> Ω_init
            get_sqrtΩ = (p) -> Ω_init
            get_sqrtΩinv = (p) -> Ω_init
        else
            # get_sqrtΩ = (p) -> 0.1 .+ log.(1.0 .+ p .^ 2) #1.0.+exp.(-p.+shift)
            get_sqrtΩ = (p) -> softmax(p)
            get_sqrtΩinv = (p) -> 1.0 ./ get_sqrtΩ(p)  #sigmoid.(p.-shift)

            get_Ω = (p) -> get_sqrtΩ(p) .^ 2
            # get_sqrtΩ = (p) -> 1.0.+exp.(-p.+shift)
            # get_sqrtΩinv = (p) -> sigmoid.(p.-shift)
        end

        # parameter indices: assume computed same way as by CTRNN1
        i = 1
        logτ_idx = i:i
        i += 1
        W_idx = i:i+state_dim*hidden_dim-1
        i += state_dim * hidden_dim
        A_idx = i:i+hidden_dim*state_dim-1

        new{
            typeof(initial_params),
            typeof(get_Ω),
            typeof(get_sqrtΩ),
            typeof(get_sqrtΩinv),
            typeof(ϵ),
        }(
            "ISSCTRNN1",
            hidden_dim,
            state_dim,
            hidden_dim,
            initial_params,
            get_Ω,
            get_sqrtΩ,
            get_sqrtΩinv,
            ϵ,
            shift,
            logτ_idx,
            W_idx,
            A_idx,
            enable,
            scale_W,
        )
    end
end

""" Modify p_sub to guarantee stability """
function (f::ISSCTRNN1)(p_sub, p_link, p; T = Float64, rhs = (one(T) - T(f.ϵ)))
    if f.enable
        sqrtΩ = f.get_sqrtΩ(p)
        invsqrtΩ = f.get_sqrtΩinv(p)

        τ = exp.(@view p_sub[f.logτ_idx])
        W = reshape(view(p_sub, f.W_idx), f.state_dim, f.hidden_dim)
        A = reshape(view(p_sub, f.A_idx), f.hidden_dim, f.state_dim)

        mat = T(0.5) * τ .* (sqrtΩ .* A) * (W' .* invsqrtΩ)'
        factor = one(T) / (one(T) + relu(maxeig(mat + mat') - one(T) + T(f.ϵ)))

        if f.scale_W
            _p_sub = vcat(
                p_sub[1],
                factor * view(p_sub, f.W_idx),
                view(p_sub, f.A_idx),
                p_sub[f.A_idx[end]+1:end],
            )
        else
            _p_sub = vcat(
                p_sub[1],
                view(p_sub, f.W_idx),
                factor * view(p_sub, f.A_idx),
                p_sub[f.A_idx[end]+1:end],
            )
        end

        return _p_sub, p_link
    else
        return p_sub, p_link
    end
end

"""  Debug: Computes the largest eigenvalue to ensure Ω(τAW - I) + (τW'A' - I)Ω ≺ 0 """
function debug(f::ISSCTRNN1, p_sub, p_link, p)
    Ω = f.get_Ω(p)
    τ = exp.(view(p_sub, f.logτ_idx))
    W = reshape(view(p_sub, f.W_idx), f.state_dim, f.hidden_dim)
    A = reshape(view(p_sub, f.A_idx), f.hidden_dim, f.state_dim)

    mat = Ω .* τ .* A * W - diagm(Ω)
    return maxeig(mat + mat')
end

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

""" LDS condition for CTRNN1
        Enforce Ω(I - τAW) + (I - τW'A')Ω ≻ 0,
        parametrize sqrt(Ω) """
struct ISSCTRNN1Alt{P,Q1,Q2,Q3,F} <: Modifier
    name::String
    param_size::Int
    state_dim::Int
    hidden_dim::Int
    initial_params::P
    get_Ω::Q1
    get_sqrtΩ::Q2
    get_sqrtΩinv::Q3
    ϵ::F
    lb::F
    logτ_idx::UnitRange{Int}
    W_idx::UnitRange{Int}
    A_idx::UnitRange{Int}
    enable::Bool
    scale_W::Bool
    function ISSCTRNN1Alt(
        state_dim,
        hidden_dim;
        enable = true,
        is_fixed = false,
        scale_W = false,
        initb = Flux.zeros,
        ϵ = 1e-3,
        lb = -3.0,
    )
        # initial parameters
        initial_params() = vec(initb(hidden_dim))

        # function to fetch parameters: if fixed, ignore parameters
        if is_fixed
            Ω_init = ones(hidden_dim)
            get_Ω = (p) -> Ω_init
            get_sqrtΩ = (p) -> Ω_init
            get_sqrtΩinv = (p) -> Ω_init
        else
            get_sqrtΩ = (p) -> 0.5 .+ exp.(p) #log.(1.0 .+ p.^2)
            get_Ω = (p) -> get_sqrtΩ(p) .^ 2
            get_sqrtΩinv = (p) -> 1.0 ./ get_sqrtΩ(p)
        end

        # parameter indices: assume computed same way as by CTRNN1
        i = 1
        logτ_idx = i:i
        i += 1
        W_idx = i:i+state_dim*hidden_dim-1
        i += state_dim * hidden_dim
        A_idx = i:i+hidden_dim*state_dim-1

        new{
            typeof(initial_params),
            typeof(get_Ω),
            typeof(get_sqrtΩ),
            typeof(get_sqrtΩinv),
            typeof(ϵ),
        }(
            "ISSCTRNN1Alt",
            hidden_dim,
            state_dim,
            hidden_dim,
            initial_params,
            get_Ω,
            get_sqrtΩ,
            get_sqrtΩinv,
            ϵ,
            lb,
            logτ_idx,
            W_idx,
            A_idx,
            enable,
            scale_W,
        )
    end
end

""" Modify p_sub to guarantee stability """
function (f::ISSCTRNN1Alt)(p_sub, p_link, p; T = Float64, rhs = (one(T) - T(f.ϵ)))
    if f.enable
        sqrtΩ = f.get_sqrtΩ(p)
        # Ω = get_Ω(p)

        τ = exp.(@view p_sub[f.logτ_idx])
        W = reshape(view(p_sub, f.W_idx), f.state_dim, f.hidden_dim)
        A = reshape(view(p_sub, f.A_idx), f.hidden_dim, f.state_dim)

        # mat = T(0.5)*τ.*(sqrtΩ.*A)*(W'.*sqrtΩ)'
        mat = T(0.5) * τ .* (sqrtΩ .* A) * W
        factor = one(T) / (one(T) + relu(maxeig(mat + mat') - one(T) + T(f.ϵ)))

        if f.scale_W
            _p_sub = vcat(
                p_sub[1],
                factor * vec((W' .* sqrtΩ)'),
                view(p_sub, f.A_idx),
                p_sub[f.A_idx[end]+1:end],
            )
        else
            _p_sub = vcat(
                p_sub[1],
                vec((W' .* sqrtΩ)'),
                factor * view(p_sub, f.A_idx),
                p_sub[f.A_idx[end]+1:end],
            )
        end

        return _p_sub, p_link
    else
        return p_sub, p_link
    end
end

"""  Debug: Computes the largest eigenvalue to ensure Ω(τAW - I) + (τW'A' - I)Ω ≺ 0 """
function debug(f::ISSCTRNN1Alt, p_sub, p_link, p)
    Ω = f.get_Ω(p)
    τ = exp.(view(p_sub, f.logτ_idx))
    W = reshape(view(p_sub, f.W_idx), f.state_dim, f.hidden_dim)
    A = reshape(view(p_sub, f.A_idx), f.hidden_dim, f.state_dim)

    mat = Ω .* τ .* A * W - diagm(Ω)
    return maxeig(mat + mat')
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
# @test_throw, with additional testing for the exception message
macro test_throws_message(f, typ, ex...)
    quote
        msg = ""
        @test_throws $(esc(typ)) try
            $(esc(ex...))
        catch err
            msg = sprint(showerror, err)
            rethrow()
        end

        if !$(esc(f))(msg)
            # @test should return its result, but doesn't
            @error "Failed to validate error message\n$msg"
        end
        @test $(esc(f))(msg)
    end
end

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# Run some code on-device, returning captured standard output
macro on_device(ex)
    quote
        let
            function kernel()
                $(esc(ex))
                return
            end

            @cuda kernel()
            synchronize()
        end
    end
end

function julia_cmd(cmd)
    return `
        $(Base.julia_cmd())
        --color=$(Base.have_color ? "yes" : "no")
        --compiled-modules=$(Base.JLOptions().use_compiled_modules != 0 ? "yes" : "no")
        --history-file=no
        --startup-file=$(Base.JLOptions().startupfile != 2 ? "yes" : "no")
        --code-coverage=$(["none", "user", "all"][1+Base.JLOptions().code_coverage])
        $cmd
    `
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::Int32) =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", Int32, Tuple{Int32}, i)
@inline sink(i::Int64) =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", Int64, Tuple{Int64}, i)

# a lightweight CUDA array type for testing purposes
## ctor & finalizer
mutable struct CuTestArray{T,N}
    buf::Mem.Buffer
    shape::NTuple{N,Int}
    function CuTestArray{T,N}(shape::NTuple{N,Int}) where {T,N}
        len = prod(shape)
        buf = Mem.alloc(len*sizeof(T))

        obj = new{T,N}(buf, shape)
        finalizer(unsafe_free!, obj)
        return obj
    end
end
function unsafe_free!(a::CuTestArray)
    CUDAdrv.isvalid(a.buf.ctx) && Mem.free(a.buf)
end
## memory copy operations
function CuTestArray(src::Array{T,N}) where {T,N}
    dst = CuTestArray{T,N}(size(src))
    Mem.upload!(dst.buf, pointer(src), length(src) * sizeof(T))
    return dst
end
function Base.Array(src::CuTestArray{T,N}) where {T,N}
    dst = Array{T,N}(undef, src.shape)
    Mem.download!(pointer(dst), src.buf, prod(src.shape) * sizeof(T))
    return dst
end
## conversions
function CUDAnative.cudaconvert(a::CuTestArray{T,N}) where {T,N}
    ptr = Base.unsafe_convert(Ptr{T}, a.buf)
    devptr = CUDAnative.DevicePtr{T,AS.Global}(ptr)
    CuDeviceArray{T,N,AS.Global}(a.shape, devptr)
end


text:
# Routines to compute surface temperatures and their jumps

@inline function surface_temperature_jump!(Δus::VectorData{N},sys::HeatConduction{NX,NY,N},t::Real) where {NX,NY,N}
    #assign_temperature!(sys.Sb,sys.bodies,sys.bodytemps,t)
    fill!(sys.Sb,0.0)
    surface_temperature_jump!(Δus,sys.Sb,sys)
    return Δus
  end

@inline surface_temperature_jump!(ΔTs::ScalarData{N},Ts::ScalarData{N},
                               ::HeatConduction{NX,NY,N,MT,InternalProblem}) where {NX,NY,N,MT} =
                               (ΔTs .= -Ts; ΔTs)

@inline surface_temperature_jump!(ΔTs::ScalarData{N},Ts::ScalarData{N},
                               ::HeatConduction{NX,NY,N,MT,ExternalProblem}) where {NX,NY,N,MT} =
                               (ΔTs .= Ts; ΔTs)

@inline surface_temperature_jump!(ΔTs::ScalarData{N},Ts::ScalarData{N},
                               ::HeatConduction{NX,NY,N,MT,ExternalInternalProblem}) where {NX,NY,N,MT} =
                               (ΔTs .= 0.0; ΔTs)

@inline function surface_temperature!(T̄s::ScalarData{N},sys::HeatConduction{NX,NY,N},t::Real) where {NX,NY,N}
    #assign_temperature!(sys.Sb,sys.bodies,sys.bodytemps,t)
    fill!(sys.Sb,0.0)
    surface_temperature!(T̄s,sys.Sb,sys)
    return T̄s
  end

@inline surface_temperature!(T̄s::ScalarData{N},Ts::ScalarData{N},
                               ::HeatConduction{NX,NY,N,MT,ExternalInternalProblem}) where {NX,NY,N,MT} =
                               (T̄s .= Ts; T̄s)

@inline surface_temperature!(T̄s::ScalarData{N},Ts::ScalarData{N},
                               ::HeatConduction{NX,NY,N,MT,SD}) where {NX,NY,N,MT,SD} =
                               (T̄s .= 0.5*Ts; T̄s)

@inline function relative_surface_temperature!(T̄sr::VectorData{N},sys::HeatConduction{NX,NY,N},t::Real) where {NX,NY,N}
  assign_temperature!(sys.Sb,sys.bodies,sys.bodytemps,t)
  surface_temperature!(T̄sr,sys.Sb,sys)
  T̄sr .-= sys.Sb
end