text:
# ------------------------------------
# Tests adapted from Git2Go Library
# ------------------------------------
context("test lookup git objects") do
    test_path = joinpath(pwd(), "testrepo")
    repo = create_test_repo(test_path)
    try
        cid, tid = seed_test_repo(repo)

        c1 = lookup(GitCommit, repo, cid)
        c2 = lookup_commit(repo, cid)
        @test c1 == c2
        @test git_otype(c1) == 1 # api.OBJ_COMMIT
        @test git_otype(c2) == 1 # api.OBJ_COMMIT

        # test repo has one commit
        ctree = GitTree(c1)
        @test length(ctree) == 1

        t = lookup(GitTree, repo, tid)
        @test Oid(t) == tid
        @test isa(t, GitTree)
        entry = entry_byname(t, "README")
        @test entry != nothing

        @test filemode(entry) == 33188 #api.FILEMODE_BLOB
        obj = lookup(repo, Oid(t))
        @test isa(obj, GitTree)

        obj = revparse(repo, "HEAD")
        @test isa(obj, GitCommit)
        @test Oid(obj) == cid

        obj = revparse(repo, "HEAD^{tree}")
        @test isa(obj, GitTree)
        @test Oid(obj) == Oid(t)
    finally
        close(repo)
        LibGit2.free!(repo)
        Base.gc()
        rm(test_path, recursive=true)
    end
end

# -----------------------------------------
# Tests adapted from Ruby's Rugged Library
# -----------------------------------------
with_repo_access() do test_repo, path

   #@test LibGit2.path(test_repo) == path

   context("test lookup any object type") do
       blob = test_repo[Oid("fa49b077972391ad58037050f2a75f74e3671e92")]
       @test isa(blob, GitBlob)

       commit = test_repo[Oid("8496071c1b46c854b31185ea97743be6a8774479")]
       @test isa(commit, GitCommit)

       tag = test_repo[Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d")]
       @test isa(tag, GitTag)

       tree = test_repo[Oid("c4dc1555e4d4fa0e0c9c3fc46734c7c35b3ce90b")]
       @test isa(tree, GitTree)
   end

   context("test fail to lookup non existant objects") do
       @test_throws LibGitError{:Odb,:NotFound} test_repo[Oid("a496071c1b46c854b31185ea97743be6a8774479")]
   end

   context("test lookup object") do
        obj = test_repo[Oid("8496071c1b46c854b31185ea97743be6a8774479")]
        @test isa(obj, GitCommit)
        @test Oid("8496071c1b46c854b31185ea97743be6a8774479") == Oid(obj)
   end

   context("test objects are the same") do
        obj1 = test_repo[Oid("8496071c1b46c854b31185ea97743be6a8774479")]
        obj2 = test_repo[Oid("8496071c1b46c854b31185ea97743be6a8774479")]
        @test obj1 == obj2
   end

   context("test read raw data") do
        obj = test_repo[Oid("8496071c1b46c854b31185ea97743be6a8774479")]
        @test isa(raw(obj), OdbObject)
   end

   context("test lookup by rev") do
        obj = revparse(test_repo, "v1.0")
        @test Oid(obj) == Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d")
        obj = revparse(test_repo, "v1.0^1")
        @test Oid(obj) == Oid("8496071c1b46c854b31185ea97743be6a8774479")
   end

   context("test lookup id by rev") do
       o = revparse_oid(test_repo, "v1.0")
       @test o == Oid("0c37a5391bbff43c37f0d0371823a5509eed5b1d")
       o = revparse_oid(test_repo, "v1.0^1")
       @test o == Oid("8496071c1b46c854b31185ea97743be6a8774479")
   end
end


text:
# indicator of the Cartesian product of real binary sets

"""
  IndBinary(low, high)

Returns the function `g = ind{x : x_i == low || x_i == high}`.
"""

immutable IndBinary{T <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}} <: IndicatorNonconvex
  low::T
  high::S
end

IndBinary{T <: Real}(low::T=0.0, high::T=1.0) = IndBinary{T, T}(low, high)

IndBinary_low{T <: Real, S}(f::IndBinary{T, S}, i) = f.low
IndBinary_low{T <: AbstractArray, S}(f::IndBinary{T, S}, i) = f.low[i]
IndBinary_high{T, S <: Real}(f::IndBinary{T, S}, i) = f.high
IndBinary_high{T, S <: AbstractArray}(f::IndBinary{T, S}, i) = f.high[i]

function (f::IndBinary){T <: Real}(x::AbstractArray{T})
  for k in eachindex(x)
    if x[k] != IndBinary_low(f, k) && x[k] != IndBinary_high(f, k)
      return +Inf
    end
  end
  return 0.0
end

function prox!{T <: Real}(f::IndBinary, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  low = 0.0
  high = 0.0
  for k in eachindex(x)
    low = IndBinary_low(f, k)
    high = IndBinary_high(f, k)
    if abs(x[k] - low) < abs(x[k] - high)
      y[k] = low
    else
      y[k] = high
    end
  end
  return 0.0
end

fun_name(f::IndBinary) = "indicator of binary array"
fun_dom(f::IndBinary) = "AbstractArray{Real}"

function prox_naive{T <: Real}(f::IndBinary, x::AbstractArray{T}, gamma::Real=1.0)
  distlow = abs.(x-f.low)
  disthigh = abs.(x-f.high)
  indlow = distlow .< disthigh
  indhigh = distlow .>= disthigh
  y = f.low.*indlow + f.high.*indhigh
  return y, 0.0
end

ext: jl


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/implicit_routines/preconditioner/dense.jl
# Definitions of dense preconditioners used with implicit relaxations.
#############################################################################

"""
$(TYPEDEF)

A dense LU preconditioner for implicit McCormick relaxation.
"""
struct DenseMidInv{S <: VecOrMat{Float64}} <: AbstractPreconditionerMC
    "Storage for the midpoint matrix and inverse"
    Y::S
    "Vector of length(1) used for intermediate calculation"
    YInterval::Vector{Interval{Float64}}
    "Number of state space variables"
    nx::Int
    "Number of decision space variables"
    np::Int
end
DenseMidInv(nx::Int, np::Int) = DenseMidInv(zeros(Float64,nx,nx), zeros(Interval{Float64},1), nx, np)
function precondition!(d::DenseMidInv{S}, H::Vector{MC{N,T}}, J::Array{MC{N,T},2}) where {N, T <: RelaxTag, S <: VecOrMat{Float64}}
    for i in eachindex(J)
        @inbounds d.YInterval[1] = J[i].Intv
        @inbounds d.Y[i] = 0.5*(d.YInterval[1].lo + d.YInterval[1].hi)
    end
    F = lu!(d.Y)
    H .= F\H
    J .= F\J
    return
end
function (d::DenseMidInv)(h!::FH, hj!::FJ, nx::Int, np::Int) where {FH <: Function, FJ <: Function}
    return DenseMidInv{Array{Float64,2}}(zeros(nx,nx), zeros(Interval{Float64},1), nx, np)
end
function preconditioner_storage(d::DenseMidInv, tag::T) where T <: RelaxTag
    zeros(MC{d.np, T}, (d.nx, d.nx))
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using MonteCarlo, MonteCarloObservable, StableDQMC
using Test
using LinearAlgebra, Random, Dates
using MonteCarlo: @bm, TimerOutputs

using MonteCarlo: vmul!, lvmul!, rvmul!, rdivp!, udt_AVX_pivot!, rvadd!, vsub!
using MonteCarlo: vmin!, vmininv!, vmax!, vmaxinv!, vinv!
using MonteCarlo: BlockDiagonal, CMat64, CVec64, StructArray

struct DummyModel <: MonteCarlo.Model
    lattice
end
DummyModel() = DummyModel(SquareLattice(2))
struct DummyField <: MonteCarlo.AbstractField end

# check elementwise, not matrix norm
function check(A::Array, B::Array, atol, rtol=atol)
    for (x, y) in zip(A, B)
        if !isapprox(x, y, atol=atol, rtol=rtol)
            @info "$x ≉ $y "
            return false
        end
    end
    true
end
function check(x::Number, y::Number, atol, rtol=atol)
    result = isapprox(x, y, atol=atol, rtol=rtol)
    result || @info "$x ≉ $y"
    result
end

# In case some test failed and left behind a .jld file
for f in readdir()
    if endswith(f, ".jld") || endswith(f, "jld2")
        @warn "Removing $f"
        rm(f)
    end
end

macro benchmark_test(name, code)
    TimerOutputs.timer_expr(MonteCarlo, true, name, code)
end

@testset "All Tests" begin
    println("Utilities")
    @time @testset "Utilities" begin
        @bm function test1(x, y)
            sleep(x+y)
        end
        @bm test2(x, y) = sleep(x+y)
        # @eval MonteCarlo ... will put the body of test3, test4 in the wrong scope
        # eval(timer_expr(...)) doesn't work
        function test3(x, y)
            @benchmark_test "test1(::Any, ::Any)" begin sleep(x+y) end
        end
        test4(x, y) = @benchmark_test "test2(::Any, ::Any)" begin sleep(x+y) end

        x = code_lowered(test1, Tuple{Float64, Float64})[1]
        y = code_lowered(test3, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        x = code_lowered(test2, Tuple{Float64, Float64})[1]
        y = code_lowered(test4, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        @test !MonteCarlo.timeit_debug_enabled()
        enable_benchmarks()
        @test MonteCarlo.timeit_debug_enabled()
        disable_benchmarks()
        @test !MonteCarlo.timeit_debug_enabled()
    end
    # 0.473232 seconds (778.10 k allocations: 47.050 MiB, 2.72% gc time, 99.50% compilation time)

    println("Lattices")
    @time @testset "Lattices" begin
        include("lattices.jl")
    end
    # 22.329261 seconds (51.70 M allocations: 2.487 GiB, 3.09% gc time, 63.07% compilation time)

    println("Model")
    @time @testset "Models" begin
        include("modeltests_IsingModel.jl")
        include("modeltests_HubbardModel.jl")
    end
    # 11.418479 seconds (29.34 M allocations: 1.301 GiB, 4.78% gc time, 7.60% compilation time)

    println("Linear Algebra")
    @time @testset "Linear Algebra" begin
        include("linalg.jl")
    end
    # 72.325652 seconds (119.03 M allocations: 5.922 GiB, 3.30% gc time, 33.24% compilation time)

    println("Slice Matrices")
    @time @testset "Slice Matrices" begin
        include("slice_matrices.jl")
    end
    # 5.212298 seconds (7.85 M allocations: 424.303 MiB, 1.66% gc time, 96.86% compilation time)

    println("DQMC")
    @time @testset "Flavors" begin
        include("flavortests_DQMC.jl")
    end
    # 27.297888 seconds (55.24 M allocations: 3.149 GiB, 2.60% gc time, 2.67% compilation time)

    println("Fields")
    @time @testset "Fields" begin
        include("fields.jl")
    end

    println("Scheduler")
    @time @testset "Scheduler & (DQMC) Updates" begin
        include("updates.jl")
    end
    # 17.197345 seconds (33.06 M allocations: 1.846 GiB, 2.53% gc time, 4.69% compilation time)

    println("Measurement")
    @time @testset "Measurements" begin
        include("measurements.jl")
    end
    # 20.091780 seconds (37.56 M allocations: 1.877 GiB, 2.46% gc time, 1.43% compilation time)

    println("Integration")
    @time @testset "Integration tests" begin
        include("integration_tests.jl")
    end
    # 11.359702 seconds (19.67 M allocations: 1.132 GiB, 2.00% gc time, 28.81% compilation time)

    println("ED")
    @time @testset "Exact Diagonalization" begin
        include("ED/ED_tests.jl")
    end
    # 71.714791 seconds (52.37 M allocations: 14.521 GiB, 1.82% gc time, 4.57% compilation time)

    println("File IO")
    @time @testset "File IO" begin
        include("FileIO.jl")
    end
    # 72.591676 seconds (86.55 M allocations: 4.172 GiB, 1.27% gc time, 9.25% compilation time)
end


text:
(Sys.islinux() || Sys.iswindows()) && import ImageMagick # work around libz issues
using ImageView
using Images, OffsetArrays
using Gtk
using Test

function imshow_now(args...; kwargs...)
    guidict = imshow(args...; kwargs...)
    Gtk.showall(guidict["gui"]["window"])
    sleep(0.01)
    guidict
end

include("simple.jl")
include("contrast.jl")
include("scalesigned.jl")
include("tile.jl")
include("annotations.jl")
include("statusbar.jl")
include("test4d.jl")
imshow(img)

include("newtests.jl")

ImageView.closeall()


text:
prodir = realpath(joinpath(dirname(dirname(@__FILE__))))
srcdir = joinpath(prodir, "src")
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)

using Documenter
using IntegerTriangles

# UpdateDocs()

makedocs(
    modules = [IntegerTriangles],
    sitename = "IntegerTriangles.jl",
    clean = true,
    checkdocs = :none,
    doctest = false,
    pages = [
        "Introduction" => "introduction.md",
        "Library" => "index.md",
        "Modules" => "modules.md",
        "Aitken" => "Aitken.md",
        "Bessel1" => "Bessel1.md",
        "Binomial" => "Binomial.md",
        "Catalan" => "Catalan.md",
        "DArcais" => "DArcais.md",
        "Delannoy" => "Delannoy.md",
        "Euler" => "Euler.md",
        "EulerSec" => "EulerS.md",
        "EulerTan" => "EulerT.md",
        "Eulerian" => "Eulerian.md",
        "EulerianSO2" => "EulerianSO2.md",
        "FallingFact" => "FallingFact.md",
        "Fibonacci" => "Fibonacci.md",
        "Fine" => "Fine.md",
        "Fubini" => "Fubini.md",
        "Hermite" => "Hermite.md",
        "Laguerre" => "Laguerre.md",
        "Lah" => "Lah.md",
        "Motzkin" => "Motzkin.md",
        "Narayana" => "Narayana.md",
        "Rencontres" => "Rencontres.md",
        "RisingFact" => "RisingFact.md",
        "SchroederBig" => "SchroederB.md",
        "SchroederLittle" => "SchroederL.md",
        "StirlingCycle" => "StirlingCycle.md",
        "StirlingSet" => "StirlingSet.md",
        "Trinomial" => "Trinomial.md",
        "TTree" => "TTree.md",
        "Uni" => "Uni.md",
        "Worpitzky" => "Worpitzky.md"
    ]
)

deploydocs(
    repo = "github.com/OpenLibMathSeq/IntegerTriangles.jl.git",
    target = "build"
)


text:
using LogExpFunctions
using ChainRulesTestUtils
using ChangesOfVariables
using InverseFunctions
using OffsetArrays

using Random
using Test

Random.seed!(1234)

include("basicfuns.jl")
include("chainrules.jl")
include("inverse.jl")
include("with_logabsdet_jacobian.jl")


text:
ROS.init("testBase")
nh = ROS.NodeHandle()

@test ROS.ok()
ROS.spinOnce()
ROS.info("j")
ROS.info_stream(3.2)
ROS.warn("u")
ROS.warn_stream(2)
ROS.error("l")
ROS.error_stream(16.07)
ROS.fatal("i")
ROS.fatal_stream(true)
ROS.debug("a")
ROS.debug_stream("!")

println("All $(basename(@__FILE__)) tests passed.")


text:
function HDF5.write(parent::Union{HDF5File, HDF5Group}, data::Dict, plists...)
    for (k, v) in data
        write(parent, k, v, plists...)
    end
end

tryreadmmap(f::Union{HDF5File, HDF5Group}) = Dict(s => tryreadmmap(f[s]) for s in names(f))

HDF5.readmmap(f::Union{HDF5File, HDF5Group}) = Dict(s => readmmap(f[s]) for s in names(f))

function tryreadmmap(obj::HDF5Dataset)
    T = hdf5_to_julia(obj)
    try
        readmmap(obj, T)
    catch e
        try
            HDF5DiskArray(obj)
        catch e
            read(obj, T)
        end
    end
end

Base.lastindex(dset::HDF5Dataset, i) = size(dset, i)

function h5readmmap(filename, name::String, pv...; mode = "r")
    local dat
    fid = h5open(filename, mode, pv...)
    try
        obj = fid[name, pv...]
        dat = readmmap(obj)
    finally
        close(fid)
    end
    dat
end

function HDF5.write(parent::Union{HDF5File, HDF5Group}, name::String, data::Dict, plists...)
    g = g_create(parent, name)
    for (k, v) in data
        write(g, k, v, plists...)
    end
end


text:
using YaoBase, StaticArrays

const C = Vector{Matrix{ComplexF64}}(undef, 0)

push!(C,  Const.I2)
push!(C,  Const.X)
push!(C,  Const.Y)
push!(C,  Const.Z)

push!(C,  sqrt(im * Const.X))
push!(C,  sqrt(-im * Const.X))
push!(C,  Const.Z*sqrt(im * Const.X))
push!(C,  Const.Z*sqrt(-im * Const.X))

push!(C,  sqrt(im * Const.Y))
push!(C,  sqrt(-im * Const.Y))
push!(C,  Const.Z * C[9])
push!(C,  Const.Z * C[10])

push!(C,  sqrt(im * Const.Z))
push!(C,  sqrt(-im * Const.Z))
push!(C,  Const.X * sqrt(im * Const.Z))
push!(C,  Const.X * sqrt(-im * Const.Z))

push!(C,  sqrt(im * Const.Z) * sqrt(im * Const.X))
push!(C,  sqrt(im * Const.Z) * sqrt(-im * Const.X))
push!(C,  sqrt(-im * Const.Z) * sqrt(im * Const.X))
push!(C,  sqrt(-im * Const.Z) * sqrt(-im * Const.X))

push!(C,  sqrt(im * Const.Z) * sqrt(im * Const.Y))
push!(C,  sqrt(im * Const.Z) * sqrt(-im * Const.Y))
push!(C,  sqrt(-im * Const.Z) * sqrt(im * Const.Y))
push!(C,  sqrt(-im * Const.Z) * sqrt(-im * Const.Y))

function is_equal_clifford(x::AbstractMatrix{T1}, y::AbstractMatrix{T2}) where {T1, T2}
    k = findfirst(a->!isapprox(a, zero(a); atol=1e-5), x)
    phi = x[k] / y[k]
    for (a, b) in zip(x, y)
        if isapprox(abs(a), zero(T1); atol=1e-5) && isapprox(abs(b), zero(T1); atol=1e-5)
            continue
        else
            isapprox(phi, a / b; atol=1e-5) || return false
        end
    end
    return true
end

function find_output(x::Int, y::Int)
    out = findall(σ->is_equal_clifford(C[x] * C[y], σ), C)
    if length(out) < 1
        error("product is not in Clifford{1} group: $x * $y")
    elseif length(out) > 1
        error("product is not unique")
    else
        return out[]
    end
end

const MulTable = SMatrix{24, 24}([find_output(i, j) for i in 1:24, j in 1:24])

max_stars_repo_name: QuantumBFS/YaoGraphs

max_issues_repo_name: QuantumBFS/YaoGraphs

max_forks_repo_name: QuantumBFS/YaoGraphs


text:
function testFourierDomainProperties(daughters, isAve)
    dMags = abs.(daughters)
    lowAprxAnalyt = daughters[1,2] ./ maximum(dMags[:,2])
    if abs(lowAprxAnalyt) >= .01
        @warn "the lowest frequency wavelet has more than 1% its max at zero, so it may not be analytic. Think carefully" lowAprxAnalyt
    end
    highAprxAnalyt = abs.(daughters[end,end]) ./ maximum(dMags[:,2])
    if highAprxAnalyt >= .01
        @warn "the highest frequency wavelet has more than 1% its max at the end, so it may not be analytic. Think carefully" highAprxAnalyt
    end

    netWeight = sum(dMags, dims=2)
    peakArgs = argmax(dMags[:, (isAve + 1):end], dims=1)
    peakFreqs = dMags[:,isAve + 1:end][peakArgs]
    maxPeak = maximum(peakFreqs)
    minPeak = minimum(peakFreqs)
    ratioOfCoverage = maxPeak / minPeak
    if ratioOfCoverage > 5
        @warn "there are frequencies which are significantly more covered than others, with a ratio of" ratioOfCoverage
    end
    peakLocs = [x[1] for x in peakArgs] # the location of the peaks, irrespective of which wavelet has that peak location
    first = minimum(peakLocs) # the first peak
    last =  maximum(peakLocs)
    minimalRegionComparedToLastPeak = maximum(peakFreqs) /
        minimum(abs.(netWeight[first:last]))
    if minimalRegionComparedToLastPeak > 2
        @warn "there are wavelets whose peaks are far enough apart that the trough between them is less than half the height of the highest frequency wavelet" minimalRegionComparedToLastPeak
    end
end

function testDualCoverage(canonicalDualFrame)
    dualCover = sum(canonicalDualFrame, dims=2)
    err = norm(dualCover .- 1)
    if err > 50
        @warn "the canonical dual frame is off by $(err), consider using one of the delta dual frames"
    end
end

function testDualCoverage(β, Ŵ)
    dualCover = sum(conj.(β) .* Ŵ, dims=2)
    err = norm(dualCover .- 1)
    if err > 50
        @warn "the dual frame is off by $(err), consider using a different dual frame"
    end
end


text:
function convert{S<:Sculpt,C<:Clay}(::Type{Flex{S,EXACT,C}}, x::Flex{S,INEXACT,C})
   Flex{S,EXACT,C}(widened(x)...)
end

function convert{S<:Sculpt,C<:Clay}(::Type{Flex{S,INEXACT,C}}, x::Flex{S,EXACT,C})
   Flex{S,INEXACT,C}(widened(x)...)
end

promote_rule{S<:Sculpt, C<:Clay}(::Type{Flex{S,INEXACT,C}}, ::Type{Flex{S,EXACT,C}}) = Flex{S,INEXACT,C}

#=
    Should you want comparisons to be influenced when the comparands differ in qualia,
    define these accordingly. (These, as given below, do not use qualia in comparisons).
=#
(==){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = false
(!=){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = true
(<=){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = opened(a) < opened(b)
(>=){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = opened(a) > opened(b)
(< ){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = closed(a) < closed(b)
(> ){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = closed(a) > closed(b)
(isless){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = (a < b)
(isequal){S<:Sculpt, Q<:Qualia, R<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{S,R,C}) = (a == b)


text:
module TestMultiReturningFunction10

using Modia.StaticArrays
using Modia
@usingModiaPlot

mutable struct MbsData
    phi1::Float64
    w1::Float64
    derw1::Float64
    u1::Float64

    phi2::Vector{Float64}
    w2::Vector{Float64}
    derw2::Vector{Float64}
    residuals2::Vector{Float64}

    function MbsData(; path = "", phi1=0.0, w1=0.0, phi2, equationInfo = missing, evaluatedParameters=missing)
        println("      TestMultiReturningFunction10: init Mbs with phi2=$phi2")
        ndof = length(phi2)
        obj = new(phi1, w1, 0.0, 0.0, phi2, zeros(ndof), zeros(ndof), zeros(ndof))
        # get_sizes(evaluatedParameters, ...)
        # add_states(obj,equationInfo, copy_from_x, copy_to_der_x)
        return obj
    end
end

function myBuildFunction(model::AbstractDict, FloatType::Type, TimeType::Type, buildDict::AbstractDict,
                         modelPath::Union{Expr,Symbol,Nothing}; buildOption = "Default")
    modelPathAsString = if isnothing(modelPath); "" else string(modelPath) end
    println("  TestMultiReturningFunction10: Test output from function myBuildFunction at modelPath = \"$modelPathAsString\":\n  Code could be constructed here and merged to the model with buildOption=$buildOption")
    return nothing
end

MyModelWithBuild(; kwargs...) = Model(; _buildFunction = :myBuildFunction, kwargs...)

Mbs(; kwargs...) = Par(; _constructor = :(MbsData), _path = true, kwargs...)

function copy_from_x(mbs::MbsData, x, ibeg, iend)
end

function copy_to_der_x(mbs::MbsData, der_x, ibeg, iend)
end

function setStates(mbs::MbsData, phi1, w1)::MbsData
    mbs.phi1 = phi1
    mbs.w1   = w1
    return mbs
end

function setAccelerations1(mbs::MbsData,derw1)::MbsData
    mbs.derw1 = derw1
    return mbs
end

function setAccelerations2(mbs::MbsData,derw2)::MbsData
    mbs.derw2 = derw2
    return mbs
end

function computeForcesAndResiduals(mbs::MbsData,time)
    # m*L^2*derw + m*g*l*sin(phi) ) = u
    #   m=L=g=1
    # derw + sin(phi) = 0
    mbs.u1         = mbs.derw1 + sin(mbs.phi1)
    mbs.residuals2 = (1+time)*mbs.derw2 + sin.(mbs.phi2)
    return mbs
end

function getForces(mbs::MbsData)
    return mbs.u1
end

function getResiduals(mbs::MbsData)
    return mbs.residuals2
end

Pendulum = MyModelWithBuild(_buildOption = "MyBuildOption",
    phi1 = Var(init=pi/2),
    w1   = Var(init=0.0),
    qdd  = Var(start=zeros(2)),
    mbs  = Mbs(phi2=[1.0,2.0]),
    mbs1 = Var(hideResult=true),
    mbs2 = Var(hideResult=true),
    mbs3 = Var(hideResult=true),
    mbs4 = Var(hideResult=true),
    equations = :[
        w1   = der(phi1)
        mbs1 = setStates(mbs,phi1,w1)
        mbs2 = setAccelerations1(mbs1,der(w1))
        mbs3 = setAccelerations2(mbs2,qdd)
        mbs4 = computeForcesAndResiduals(mbs3,time)
        tau1 = implicitDependency(getForces(mbs4), der(w1))
        0    = implicitDependency(getResiduals(mbs4), qdd)
        tau1 = -0.1*w1
    ]
)

pendulum = @instantiateModel(Pendulum , unitless=true, log=false, logDetails=false, logCode=true, logStateSelection=false)
simulate!(pendulum, stopTime = 2.0, log=true)
printResultInfo(pendulum)

plot(pendulum, [("phi1", "w1"), "der(w1)", "qdd"])

simulate!(pendulum, stopTime = 2.0, log=true, merge=Map(mbs = Mbs(phi2=[10.0,20.0,30.0]), qdd = zeros(3)))
plot(pendulum, [("phi1", "w1"), "der(w1)", "qdd"], figure=2)

end


text:
using JuMP, JuMPeR, Gurobi, Random, Distributions, LinearAlgebra, DataFrames, Plots

n = 10 # Number of facilities
m = 50 # Number of customers

# Generating random data (please don't change the seeds.)
facilities = 0.6.*rand(MersenneTwister(5), n,2) .+ 0.2;
customers = rand(MersenneTwister(2), m, 2);
c = [LinearAlgebra.norm(customers[i, :] .- facilities[j, :])[1] for j=1:n, i=1:m];
f = rand(MersenneTwister(3), n)*1 .+ 5;
s = rand(MersenneTwister(4), n)*2 .+ 15;
d = rand(MersenneTwister(5), m)*0.5 .+ 0.75

# P matrix
R_D = 0.25
P = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P = (P .>= 0.2*exp(-1/R_D .* R_D)) .* P

""" Adaptive facility location model, with JuMP using robust counterpart,
    under budget uncertainty.
    Note: We restrict u and V to be elementwise positive, for tractability! """
function adaptive_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = Model(solver = GurobiSolver())

    # VARIABLES
    @variable(model, x[1:n], Bin)           # Facility locations
    @variable(model, V[1:n, 1:m, 1:m] >= 0) # Affine policy w.r.t each z
    @variable(model, u[1:n, 1:m] >= 0)

    # CONSTRAINTS
    for i = 1:n # Positivity constraints
        for j = 1:m
            infdummy = @variable(model)
            normdummy = @variable(model, [1:m])
            y = @variable(model, [1:m])
            @constraint(model, [k=1:m], normdummy[k] >= y[k])
            @constraint(model, [k=1:m], normdummy[k] >= -y[k])
            @constraint(model, [k=1:m], infdummy >= -y[k] + V[i,j,k])
            @constraint(model, [k=1:m], infdummy >= y[k] - V[i,j,k])
            @constraint(model, u[i,j] >= rho*sum(normdummy) + Gamm * infdummy)
        end
    end
    for j=1:m # Demand constraints
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [i=1:m], normdummy[i] >= y[i])
        @constraint(model, [i=1:m], normdummy[i] >= -y[i])
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @constraint(model, [i=1:m], infdummy >= -y[i] - P[i,j] + sum(V[:,:,i].*dt))
        @constraint(model, [i=1:m], infdummy >= y[i] + P[i,j] - sum(V[:,:,i].*dt))
        @constraint(model, sum(dt .* u) >= d[j] + rho*sum(normdummy) + Gamm * infdummy)
    end
    for i=1:n # Capacity constraints
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [j=1:m], normdummy[j] >= y[j])
        @constraint(model, [j=1:m], normdummy[j] >= -y[j])
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @constraint(model, [j=1:m], infdummy >= -y[j] + sum(V[:,:,j].*dt))
        @constraint(model, [j=1:m], infdummy >= y[j] - sum(V[:,:,j].*dt))
        @constraint(model, sum(dt .* u) + rho*sum(normdummy) + Gamm * infdummy <=
                    s[i] * x[i]) # capacity constraint
    end
    # OBJECTIVE
    @variable(model, F)
    infdummy = @variable(model)
    normdummy = @variable(model, [1:m])
    y = @variable(model, [1:m])
    @constraint(model, [j=1:m], normdummy[j] >= y[j])
    @constraint(model, [j=1:m], normdummy[j] >= -y[j])
    @constraint(model, [j=1:m], infdummy >= - y[j] + sum(c.* V[:,:,j]))
    @constraint(model, [j=1:m], infdummy >= y[j] - sum(c.* V[:,:,j]))
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m) +
                                rho*sum(normdummy) + Gamm * infdummy)
    @objective(model, Min, F)
    return model, x, u, V
end

""" Plots the solution of the facility location model.
Blue circles are active facilities with different capacities.
Orange plus signs are other potential facility locations.
Rays describe connections between facilities and demand nodes.
"""
function plot_solution(model, x, y, cost = nothing)
    plt = scatter(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s .* getvalue(x))
    scatter!(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s, markershape = :+)
    for i=1:n
        for j=1:m
            if getvalue(y[i,j]) >= 1e-10
                plot!([customers[j, 1], facilities[i,1]], [customers[j,2], facilities[i,2]], linewidth = getvalue(y[i,j]), legend=false)

            end
        end
    end
    if cost == nothing
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
                title = "Total cost: $(round(getobjectivevalue(model), sigdigits=5))")
    else
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
        title = "Total cost: $(round(cost,sigdigits=5))")
    end
    println("Facility cost: $(getvalue(sum(f[j] * x[j] for j = 1:n)))")
    println("Transportation cost: $(getvalue(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)))")
    return plt
end

# Adaptive model through robust counterpart
rho = 1
Gamm = 5
model, x, u, V = adaptive_facility_model(c, f, rho, Gamm)
@time solve(model)
plt = plot_solution(model, x, u)
@show plt


text:
# using Flux: conv
using LinearAlgebra
function FieldProd(l1, l2, l)
    # l1, l2, l = ranks
    if l1 == 0
        return prod0__
    elseif l2 == 0
        return prod_0_
    elseif l1 == l2 && l == 0
        return prod__0
    elseif l1 == l2 == l == 1
        return prod111
    end

    # (x, y) -> yield(f(x, y), x.grid)
    # (x,y)->yield
end

function getl(x)
    n = size(x)[end]
    if n == 1
        return 0
    end

    dims = ndims(x) - 1
    if dims == 2
        return 1
    elseif dims == 3
        return (n - 1) ÷ 2
    end

end
function tp(x, f, l, p)
    r = l .+ size(f).-1
    xl = max.(l, 1)
    xr = min.(r, size(x))
    fl = 1 .+ xl .- l
    fr = size(f) .- r .+ xr
    sum(p.(
        x[[a:b for (a, b) in zip(xl, xr)]...],
        f[[a:b for (a, b) in zip(fl, fr)]...],
    ))
end
function xcor(x, f; product = *, stride = 1, pad = 0)
    l = Iterators.product([
        a:stride:b for (a, b) in zip(ones(Int,ndims(x)) .- pad, size(x) .- size(f) .+ 1 .+ pad)
    ]...)
    [tp(x, f, l, product) for l in l]
end
function conv(x, f; product = *, boundary = :extend)
    if boundary == :smooth
        r = xcor(x, reverse(f); product)
        ix = [Int.([1, (1:a)..., a]) for a in size(r)]
        return r[ix...]
    end

    if boundary == :extend
        pad = size(f) .- 1
    elseif boundary == :same
        pad = (size(f) .- 1) .÷ 2
    end
    xcor(x, reverse(f); product, pad)

end
function Δ(x, y)
    sum(abs.(x .- y)) / sum(abs.(y))
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ c02664e7-2046-4103-8a59-dca4998638df
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
    Pkg.instantiate()

	# Pkg.status()
end

# ╔═╡ 8a90d8a0-eb33-417c-8ac3-440822ae99f3
LOAD_PATH

# ╔═╡ 3103370e-488a-4cac-9540-1ef4bec5503b
using PlutoPkgTestA

# ╔═╡ 82d46919-ccb0-4ad2-bb3d-77adfafebc4d
PlutoPkgTestA.MY_VERSION |> Text

# ╔═╡ c44e23b8-3101-11ec-2112-df6ef8652469
@__DIR__

# ╔═╡ a6cee179-f82c-4ef5-8091-cf0be115ec92
pwd()

# ╔═╡ 26283153-f597-44b9-8a17-8018ca7ca34c
Base.active_project()

# ╔═╡ 3b96bb61-08f0-4ba8-90d2-2a2be9902c2d
Base.current_project()

# ╔═╡ 8cc87597-f0c5-4902-9f3e-4ed8c6798ee9
Base.current_project(@__DIR__)

# ╔═╡ 8559b034-f7d2-4eea-a080-557f22ed98d9
Base.current_project(joinpath(@__DIR__, "..")) |> normpath

# ╔═╡ Cell order:
# ╠═c02664e7-2046-4103-8a59-dca4998638df
# ╠═8a90d8a0-eb33-417c-8ac3-440822ae99f3
# ╠═82d46919-ccb0-4ad2-bb3d-77adfafebc4d
# ╠═3103370e-488a-4cac-9540-1ef4bec5503b
# ╠═c44e23b8-3101-11ec-2112-df6ef8652469
# ╠═a6cee179-f82c-4ef5-8091-cf0be115ec92
# ╠═26283153-f597-44b9-8a17-8018ca7ca34c
# ╠═3b96bb61-08f0-4ba8-90d2-2a2be9902c2d
# ╠═8cc87597-f0c5-4902-9f3e-4ed8c6798ee9
# ╠═8559b034-f7d2-4eea-a080-557f22ed98d9


text:
using GPFlux
using LinearAlgebra
using Test, Random

@testset "square_metric" begin
  Random.seed!(3)
  T = Float64
	for N in [2, 4, 7]
    @show N
    X = rand(T, 4, N)
		v = rand(T, N)

    function tfunc(x)
			D = square_metric(x)
			(v'*Symmetric(D, :L)*v)[] |> real
    end
    @test gradient_check(tfunc, X)
  end
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
# ### Loading some packages

using GeoPhyInv
using SparseArrays
using StatsBase
using LinearAlgebra
using Random
using LinearAlgebra
using Test
using ForwardDiff
using Calculus
#src # include("core.jl")
#src # include("expt.jl")

# ### Solve for ``ψ`` in a `PoissonExpt`
# This module represents an explicit, direct sparse 2D finite-difference Poisson solver for heterogeneous media,
# i.e. media having spatially varying (space-dependent) medium parameters.
# Current implementation assumes Neumann boundary conditions at all the boundaries.

# Consider the following Poisson experiment:
# ```math
# ∇⋅(σ(x,z)∇) ψ(t) = ∇⋅(Q(x,z)∇) p(t),
# ```
# ```math
# Q = k * Q_v / η.
# ```
# We start with the dimensions and spatial grids are allocated as follows.

nx=21
nz=21
nt=4
nznx=nz*nx
mgrid=[range(-div(nz,2), step=1.0, length=nz), range(-div(nx,2), step=1.0, length=nx)]
tgrid=range(0.0,step=0.5, length=nt)
@info "Grids are all set."

# Now lets allocate the inputs for a toy experiment.
# These medium parameters are used to generate the *observed* field ``ψ``.

Qv=abs.(randn(nz,nx))
η=abs.(randn(nz,nx))
k=abs.(randn(nz,nx))
σ=abs.(randn(nz,nx))
p=randn(nz,nx,nt)
@info "Medium parameters allocated."

# ### Acquisition
# Now, we will generate an acquisition ageometry and allocate a projection matrix `ACQ`.

ageom=AGeom(mgrid, SSrcs(1), Srcs(1), Recs(30))
update!(ageom, SSrcs(), [0,0], 5, [0,2π])
update!(ageom, Recs(), [0,0], 5, [0,2π])
ACQ=SparseMatrixCSC(ageom[1],mgrid);
@info "ACQ will be used to project ψ onto receivers."

# ### Generate `PoissonExpt` and then applying `mod!`
# This will first
# * apply operator ``A=∇⋅(Q(x,z)∇)`` on a field ``p``;
# * then apply ``(∇⋅(σ(x,z)∇))^{-1}`` in order to solve for ``ψ``;
# * finally, records ``ψ`` at the receiver locations to generate data.

paE=PoissonExpt(p, tgrid, mgrid, Qv, k, η, σ, ACQ)
GeoPhyInv.mod!(paE)

# ### Extracting data from `Expt`
data=paE[:data]
@info string("The dimensions of data are (nt,nr)=",size(data))


text:
using Base.Test
using OffsetArrays

@test isempty(detect_ambiguities(OffsetArrays, Base, Core))

# Basics
for n = 0:5
    a = OffsetArray(ones(Int,ntuple(d->1,n)), ntuple(x->x-1,n))
    @test length(linearindices(a)) == 1
    @test indices(a) == ntuple(x->x:x,n)
    @test a[1] == 1
end
a0 = reshape([3])
a = OffsetArray(a0)
@test indices(a) == ()
@test ndims(a) == 0
@test a[] == 3

y = OffsetArray(Float64, -1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
@test indices(y) == (-1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
y[-1,-7,-128,-5,-1,-3,-2,-1] = 14
y[-1,-7,-128,-5,-1,-3,-2,-1] += 5
@test y[-1,-7,-128,-5,-1,-3,-2,-1] == 19

r = -2:5
y = OffsetArray(r, r)
@test indices(y) == (r,)
y = OffsetArray(r, (r,))
@test indices(y) == (r,)

A0 = [1 3; 2 4]
A = OffsetArray(A0, (-1,2))                   # LinearFast
S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))   # LinearSlow
@test indices(A) == indices(S) == (0:1, 3:4)
@test_throws ErrorException size(A)
@test_throws ErrorException size(A, 1)
@test A == OffsetArray(A0, 0:1, 3:4)
@test_throws DimensionMismatch OffsetArray(A0, 0:2, 3:4)
@test_throws DimensionMismatch OffsetArray(A0, 0:1, 2:4)

# Scalar indexing
@test A[0,3] == A[1] == S[0,3] == S[1] == 1
@test A[1,3] == A[2] == S[1,3] == S[2] == 2
@test A[0,4] == A[3] == S[0,4] == S[3] == 3
@test A[1,4] == A[4] == S[1,4] == S[4] == 4
@test_throws BoundsError A[1,1]
@test_throws BoundsError S[1,1]

# Vector indexing
@test A[:, 3] == S[:, 3] == OffsetArray([1,2], (A.offsets[1],))
@test A[:, 4] == S[:, 4] == OffsetArray([3,4], (A.offsets[1],))
@test_throws BoundsError A[:, 1]
@test_throws BoundsError S[:, 1]
@test A[0, :] == S[0, :] == OffsetArray([1,3], (A.offsets[2],))
@test A[1, :] == S[1, :] == OffsetArray([2,4], (A.offsets[2],))
@test_throws BoundsError A[2, :]
@test_throws BoundsError S[2, :]
@test A[0:1, 3] == S[0:1, 3] == [1,2]
@test A[[1,0], 3] == S[[1,0], 3] == [2,1]
@test A[0, 3:4] == S[0, 3:4] == [1,3]
@test A[1, [4,3]] == S[1, [4,3]] == [4,2]
@test A[:, :] == S[:, :] == A

# CartesianIndexing
@test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
@test_throws BoundsError A[CartesianIndex(1,1)]
@test_throws BoundsError S[CartesianIndex(1,1)]
@test eachindex(A) == 1:4
@test eachindex(S) == CartesianRange((0:1,3:4))

# view
S = view(A, :, 3)
@test S == OffsetArray([1,2], (A.offsets[1],))
@test S[0] == 1
@test S[1] == 2
@test_throws BoundsError S[2]
@test indices(S) === (0:1,)
S = view(A, 0, :)
@test S == OffsetArray([1,3], (A.offsets[2],))
@test S[3] == 1
@test S[4] == 3
@test_throws BoundsError S[1]
@test indices(S) === (3:4,)
S = view(A, 0:0, 4)
@test S == [3]
@test S[1] == 3
@test_throws BoundsError S[0]
@test indices(S) === (Base.OneTo(1),)
S = view(A, 1, 3:4)
@test S == [2,4]
@test S[1] == 2
@test S[2] == 4
@test_throws BoundsError S[3]
@test indices(S) === (Base.OneTo(2),)
S = view(A, :, :)
@test S == A
@test S[0,3] == S[1] == 1
@test S[1,3] == S[2] == 2
@test S[0,4] == S[3] == 3
@test S[1,4] == S[4] == 4
@test_throws BoundsError S[1,1]
@test indices(S) === (0:1, 3:4)

# iteration
for (a,d) in zip(A, A0)
    @test a == d
end

# show
io = IOBuffer()
show(io, A)
str = takebuf_string(io)
@test str == "[1 3; 2 4]"
show(io, S)
str = takebuf_string(io)
@test str == "[1 3; 2 4]"
show(io, MIME("text/plain"), A)
strs = split(strip(takebuf_string(io)), '\n')
@test strs[2] == " 1  3"
@test strs[3] == " 2  4"
v = OffsetArray(rand(3), (-2,))
show(io, v)
str = takebuf_string(io)
show(io, parent(v))
@test str == takebuf_string(io)
function cmp_showf(printfunc, io, A)
    ioc = IOContext(io, limit=true, compact=true)
    printfunc(ioc, A)
    str1 = takebuf_string(io)
    printfunc(ioc, parent(A))
    str2 = takebuf_string(io)
    @test str1 == str2
end
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,5), (10,-9)))       # rows&cols fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,5), (10,-9)))    # columns fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,10^3), (10,-9)))    # rows fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neither fits

# Similar
B = similar(A, Float32)
@test isa(B, OffsetArray{Float32,2})
@test indices(B) === indices(A)
B = similar(A, (3,4))
@test isa(B, Array{Int,2})
@test size(B) == (3,4)
@test indices(B) === (Base.OneTo(3), Base.OneTo(4))
B = similar(A, (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test indices(B) === (-3:3, 1:4)
B = similar(parent(A), (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test indices(B) === (-3:3, 1:4)

# Reshape
B = reshape(A0, -10:-9, 9:10)
@test isa(B, OffsetArray{Int,2})
@test parent(B) === A0
@test indices(B) == (-10:-9, 9:10)
B = reshape(A, -10:-9, 9:10)
@test isa(B, OffsetArray{Int,2})
@test parent(B) === A0
@test indices(B) == (-10:-9, 9:10)
b = reshape(A, -7:-4)
@test indices(b) == (-7:-4,)
@test isa(parent(b), Vector{Int})
@test parent(b) == A0[:]
a = OffsetArray(rand(3,3,3), -1:1, 0:2, 3:5)
@test_throws ArgumentError reshape(a, Val{2})
@test_throws ArgumentError reshape(a, Val{4})

# Indexing with OffsetArray indices
i1 = OffsetArray([2,1], (-5,))
i1 = OffsetArray([2,1], -5)
b = A0[i1, 1]
@test indices(b) === (-4:-3,)
@test b[-4] == 2
@test b[-3] == 1
b = A0[1,i1]
@test indices(b) === (-4:-3,)
@test b[-4] == 3
@test b[-3] == 1
v = view(A0, i1, 1)
@test indices(v) === (-4:-3,)
v = view(A0, 1:1, i1)
@test indices(v) === (Base.OneTo(1), -4:-3)

# logical indexing
@test A[A .> 2] == [3,4]

# copy!
a = OffsetArray{Int}((-3:-1,))
fill!(a, -1)
copy!(a, (1,2))   # non-array iterables
@test a[-3] == 1
@test a[-2] == 2
@test a[-1] == -1
fill!(a, -1)
copy!(a, -2, (1,2))
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copy!(a, 1, (1,2))
fill!(a, -1)
copy!(a, -2, (1,2,3), 2)
@test a[-3] == -1
@test a[-2] == 2
@test a[-1] == 3
@test_throws BoundsError copy!(a, -2, (1,2,3), 1)
fill!(a, -1)
copy!(a, -2, (1,2,3), 1, 2)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2

b = 1:2    # copy between AbstractArrays
bo = OffsetArray(1:2, (-3,))
@test_throws BoundsError copy!(a, b)
fill!(a, -1)
copy!(a, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
fill!(a, -1)
copy!(a, -2, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copy!(a, -4, bo)
@test_throws BoundsError copy!(a, -1, bo)
fill!(a, -1)
copy!(a, -3, b, 2)
@test a[-3] == 2
@test a[-2] == a[-1] == -1
@test_throws BoundsError copy!(a, -3, b, 1, 4)
am = OffsetArray{Int}((1:1, 7:9))  # for testing linear indexing
fill!(am, -1)
copy!(am, b)
@test am[1] == 1
@test am[2] == 2
@test am[3] == -1
@test am[1,7] == 1
@test am[1,8] == 2
@test am[1,9] == -1

dest = similar(am)
map!(+, dest, am, am)
@test dest[1,7] == 2
@test dest[1,8] == 4
@test dest[1,9] == -2

A = OffsetArray(rand(4,4), (-3,5))
@test maximum(A) == maximum(parent(A))
@test minimum(A) == minimum(parent(A))
@test extrema(A) == extrema(parent(A))
C = similar(A)
cumsum!(C, A, 1)
@test parent(C) == cumsum(parent(A), 1)
@test parent(cumsum(A, 1)) == cumsum(parent(A), 1)
cumsum!(C, A, 2)
@test parent(C) == cumsum(parent(A), 2)
R = similar(A, (1:1, 6:9))
maximum!(R, A)
@test parent(R) == maximum(parent(A), 1)
R = similar(A, (-2:1, 1:1))
maximum!(R, A)
@test parent(R) == maximum(parent(A), 2)
amin, iamin = findmin(A)
pmin, ipmin = findmin(parent(A))
@test amin == pmin
@test A[iamin] == amin
@test amin == parent(A)[ipmin]
amax, iamax = findmax(A)
pmax, ipmax = findmax(parent(A))
@test amax == pmax
@test A[iamax] == amax
@test amax == parent(A)[ipmax]

v  = OffsetArray([1,1e100,1,-1e100], (-3,))*1000
v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000
@test isa(v, OffsetArray)
cv  = OffsetArray([1,1e100,1e100,2], (-3,))*1000
cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000
@test isequal(cumsum_kbn(v), cv)
@test isequal(cumsum_kbn(v2), cv2)
@test isequal(sum_kbn(v), sum_kbn(parent(v)))

io = IOBuffer()
writedlm(io, A)
seek(io, 0)
@test readdlm(io, eltype(A)) == parent(A)

amin, amax = extrema(parent(A))
@test clamp(A, (amax+amin)/2, amax) == clamp(parent(A), (amax+amin)/2, amax)

@test unique(A, 1) == parent(A)
@test unique(A, 2) == parent(A)
v = OffsetArray(rand(8), (-2,))
@test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
@test sortrows(A) == OffsetArray(sortrows(parent(A)), A.offsets)
@test sortcols(A) == OffsetArray(sortcols(parent(A)), A.offsets)
@test sort(A, 1) == OffsetArray(sort(parent(A), 1), A.offsets)
@test sort(A, 2) == OffsetArray(sort(parent(A), 2), A.offsets)

@test mapslices(v->sort(v), A, 1) == OffsetArray(mapslices(v->sort(v), parent(A), 1), A.offsets)
@test mapslices(v->sort(v), A, 2) == OffsetArray(mapslices(v->sort(v), parent(A), 2), A.offsets)

@test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]])
@test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]])
@test flipdim(A, 1) == OffsetArray(flipdim(parent(A), 1), A.offsets)
@test flipdim(A, 2) == OffsetArray(flipdim(parent(A), 2), A.offsets)

@test A+1 == OffsetArray(parent(A)+1, A.offsets)
@test 2*A == OffsetArray(2*parent(A), A.offsets)
@test A+A == OffsetArray(parent(A)+parent(A), A.offsets)
@test A.*A == OffsetArray(parent(A).*parent(A), A.offsets)

B = fill(5, 1:3, -1:1)
@test indices(B) == (1:3,-1:1)
@test all(B.==5)

# @unsafe
a = OffsetArray(zeros(7), -3:3)
unsafe_fill!(x) = @unsafe(for i in indices(x,1); x[i] = i; end)
function unsafe_sum(x)
    s = zero(eltype(x))
    @unsafe for i in indices(x,1)
        s += x[i]
    end
    s
end
unsafe_fill!(a)
for i = -3:3
    @test a[i] == i
end
@test unsafe_sum(a) == 0


text:
using MLJBase
using MLJSerialization
using Test
import XGBoost
using MLJXGBoostInterface
using Distributions
import StableRNGs
rng = StableRNGs.StableRNG(123)

@test_logs (:warn, r"Changing ") XGBoostClassifier(objective="wrong")
@test_logs (:warn, r"Changing ") XGBoostCount(objective="wrong")
@test_logs (:warn, r"Your") XGBoostRegressor(objective="binary:logistic")

@testset "Binary Classification" begin
    plain_classifier = XGBoostClassifier(num_round=100, seed=0)

    N=2
    X = (x1=rand(1000), x2=rand(1000), x3=rand(1000))
    ycat = map(X.x1) do x
        string(mod(round(Int, 10*x), N))
    end |> categorical

    train, test = partition(eachindex(ycat), 0.6)

    #fitresult, cache, report = MLJBase.fit(plain_classifier, 1, X, ycat;)

    m = machine(plain_classifier,X,ycat)
    fit!(m,verbosity = 0)
end

## REGRESSOR
plain_regressor = XGBoostRegressor()
n,m = 10^3, 5 ;
features = rand(rng, n,m);
weights = rand(rng, -1:1,m);
labels = features * weights;
features = MLJBase.table(features)
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 1, features, labels);
rpred = predict(plain_regressor, fitresultR, features);

plain_regressor.objective = "gamma"
labels = abs.(labels)
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 1, features, labels);
rpred = predict(plain_regressor, fitresultR, features);

importances = reportR.feature_importances

# serialization:
serializable_fitresult =
    MLJBase.save("mymodel", plain_regressor, fitresultR)

restored_fitresult = MLJBase.restore("mymodel",
                                     plain_regressor,
                                     serializable_fitresult)
@test predict(plain_regressor, restored_fitresult, features) ≈ rpred

## COUNT

count_regressor = XGBoostCount(num_round=10)

X = randn(rng, 100, 3) .* randn(rng, 3)'
Xtable = table(X)

α = 0.1
β = [-0.3, 0.2, -0.1]
λ = exp.(α .+ X * β)
ycount = [rand(rng, Poisson(λᵢ)) for λᵢ ∈ λ]

fitresultC, cacheC, reportC = MLJBase.fit(count_regressor, 1, Xtable, ycount);
cpred = predict(count_regressor, fitresultC, Xtable);

importances = reportC.feature_importances

## CLASSIFIER

plain_classifier = XGBoostClassifier(num_round=100, seed=0)

# test binary case:
N=2
X = (x1=rand(rng, 1000), x2=rand(rng, 1000), x3=rand(rng, 1000))
ycat = map(X.x1) do x
    string(mod(round(Int, 10*x), N))
end |> categorical
y = identity.(ycat) # make plain Vector with categ. elements
train, test = partition(eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(plain_classifier, 1,
                                            selectrows(X, train), y[train];)
yhat = mode.(predict(plain_classifier, fitresult, selectrows(X, test)))
misclassification_rate = sum(yhat .!= y[test])/length(test)
@test misclassification_rate < 0.01

importances = report.feature_importances

# Multiclass{10} case:
N=10
X = (x1=rand(rng, 1000), x2=rand(rng, 1000), x3=rand(rng, 1000))
ycat = map(X.x1) do x
    string(mod(round(Int, 10*x), N))
end |> categorical
y = identity.(ycat) # make plain Vector with categ. elements

train, test = partition(eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(plain_classifier, 1,
                                            selectrows(X, train), y[train];)
yhat = mode.(predict(plain_classifier, fitresult, selectrows(X, test)))
misclassification_rate = sum(yhat .!= y[test])/length(test)
@test misclassification_rate < 0.01

# check target pool preserved:
X = (x1=rand(rng, 400), x2=rand(rng, 400), x3=rand(rng, 400))
ycat = vcat(fill('x', 100), fill('y', 100), fill('z', 200)) |>categorical
y = identity.(ycat)
train, test = partition(eachindex(y), 0.5)
@test length(unique(y[train])) == 2
@test length(unique(y[test])) == 1
fitresult, cache, report = MLJBase.fit(plain_classifier, 1,
                                            selectrows(X, train), y[train];)
yhat = predict_mode(plain_classifier, fitresult, selectrows(X, test))
@test Set(MLJBase.classes(yhat[1])) == Set(MLJBase.classes(y[train][1]))

# serialization:
serializable_fitresult =
    MLJBase.save("mymodel", plain_classifier, fitresult)

restored_fitresult = MLJBase.restore("mymodel",
                                     plain_classifier,
                                     serializable_fitresult)

@test predict_mode(plain_classifier, restored_fitresult, selectrows(X, test)) ==
    yhat

## MACHINE INTEGRATION

# count regressor (`count_regressor`, `Xtable` and `ycount`
# defined above):

mach = machine(count_regressor, Xtable, ycount) |> fit!
yhat = predict(mach, Xtable)

# serialize:
io = IOBuffer()
MLJBase.save(io, mach)

# deserialize:
seekstart(io)
mach2 = machine(io)
close(io)

# compare:
@test predict(mach2, Xtable) ≈ yhat

# classifier
mach = machine(plain_classifier, X, y) |> fit!
yhat = predict_mode(mach, X);

# serialize:
io = IOBuffer()
MLJBase.save(io, mach)

# deserialize:
seekstart(io)
mach2 = machine(io)

# compare:
@test predict_mode(mach2, X) == yhat


text:
typealias IntegerVector{I<:Integer} Vector{I}
typealias RealVector{N<:Real} Vector{N}
typealias RealMatrix{N<:Real} Matrix{N}

typealias RealLowerTriangular{T<:Real, S<:AbstractMatrix} LowerTriangular{T, S}

multivecs{T}(::Type{T}, n::Int) = [T[] for _ =1:n]


text:
mutable struct SystemCatalog
    data::Dict{DataType, Dict{String, SystemDescriptor}}
end

function get_system_descriptor(
    category::Type{<:SystemCategory},
    catalog::SystemCatalog,
    name::String,
)
    data = catalog.data
    if haskey(data, category) && haskey(data[category], name)
        return data[category][name]
    else
        error("System $(name) of Category $(category) not found in current SystemCatalog")
    end
end

function get_system_descriptors(category::Type{<:SystemCategory}, catalog::SystemCatalog)
    data = catalog.data
    if haskey(data, category)
        array = SystemDescriptor[descriptor for descriptor in values(data[category])]
        return array
    else
        error("Category $(category) not found in SystemCatalog")
    end
end

function list_categories()
    catalog = SystemCatalog()
    return list_categories(catalog)
end

list_categories(c::SystemCatalog) = sort!([x for x in (keys(c.data))], by = x -> string(x))

function SystemCatalog(system_catalogue::Array{SystemDescriptor} = SYSTEM_CATALOG)
    data = Dict{DataType, Dict{String, SystemDescriptor}}()
    for descriptor in system_catalogue
        category = get_category(descriptor)
        if haskey(data, category)
            push!(data[category], (descriptor.name => descriptor))
        else
            push!(data, (category => Dict{String, SystemDescriptor}()))
            push!(data[category], (descriptor.name => descriptor))
        end
    end
    return SystemCatalog(data)
end

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
using StaticArrays

const σz = SA[1.0+0.0im 0.0+0.0im; 0.0+0.0im -1.0+0.0im]
const σx = SA[0.0+0.0im 1.0+0.0im; 1.0+0.0im  0.0+0.0im]
const σy = SA[0.0+0.0im 0.0-1.0im; 0.0+1.0im  0.0+0.0im]
const zm = SA[0.0+0.0im 0.0+0.0im; 0.0+0.0im  0.0+0.0im]
const eye = SA[1.0+0.0im 0.0+0.0im; 0.0+0.0im  1.0+0.0im]

@inline function comm(A,B)
    A*B - B*A
end

@inline function acomm(A,B)
    A*B + B*A
end


text:
module OpticSimVis

greet() = println("Hello World!")

using CoordinateTransformations
using Rotations
import GeometryBasics
# using GeometryBasics: HyperRectangle, Vec, Point, Mesh
using Colors: RGBA, RGB
import Colors
import ColorSchemes
import FileIO
import MeshIO
import UUIDs

using OpticSim, OpticSim.Geometry, OpticSim.Emitters
using StaticArrays

using Glimmer, Glimmer.FlexUI

include("Misc.jl")

include("OpticSimAPI/General.jl")
include("OpticSimAPI/Emitters.jl")
include("OpticSimAPI/HeadEyeModel.jl")

# export UIVariables
# export FlexUI

#------------------------------------------------------------------------------
# EXports
#------------------------------------------------------------------------------
export  Glimmer,
        FlexUI,
        App,
        set_camera!,
        prop, prop!,

        Scene,
        set_Y_up!,
        set_Z_up!,
        grid!,
        axes!,
        root, root!,
        parent, parent!,
        material, material!,
        Material,
        transform, lookAt,
        rotation,
        local_tr, local_tr!,
        cameraTransform!, cameraPlanes!,
        Box,
        Mesh,
        Arrow,
        Axes,
        LineSegments,
        PointCloud,
        dataPath,
        cameraTransform!,
        url,
        clear,
        DummyExport

export  Container,
        Slider,
        Button,
        MeshCatViewer,
        Label,
        Image,
        PanZoom,
        Field,
        ButtonToggle,
        RadioGroup,
        CheckBox,
        ExpansionPanel,
        Accordion,
        Tabs,
        Tab,
        Divider,
        Card,

        VContainer,
        HContainer,
        HContainerSpace,
        H1Label,
        H2Label,
        H3Label,
        H4Label,

        DummyExport

export  BasicValidation, on, Variable
export  FlexUI, addVariable!, prop, prop!, controls, controls!, renderFunction, renderFunction!

end # module


text:
function Base.:+(G1::GMTgrid, G2::GMTgrid)
# Add two grids, element by element. Inherit header parameters from G1 grid
	(size(G1.z) != size(G2.z)) && error("Grids have different sizes, so they cannot be added.")
	G3 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .+ G2.z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	grd_min_max!(G3)		# Also take care of NaNs
	return G3
end

# ---------------------------------------------------------------------------------------------------
Base.:+(shift::Real, G1::GMTgrid) = Base.:+(G1::GMTgrid, shift::Real)
function Base.:+(G1::GMTgrid, shift::Real)
	_shift = convert(eltype(G1.z), shift)
	G2 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .+ _shift, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	G2.range[5:6] .+= shift
	return G2
end

# ---------------------------------------------------------------------------------------------------
function Base.:-(G1::GMTgrid, G2::GMTgrid)
# Subtract two grids, element by element. Inherit header parameters from G1 grid
	(size(G1.z) != size(G2.z)) && error("Grids have different sizes, so they cannot be subtracted.")
	G3 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .- G2.z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	grd_min_max!(G3)		# Also take care of NaNs
	return G3
end

# ---------------------------------------------------------------------------------------------------
function Base.:-(G1::GMTgrid, shift::Real)
	_shift = convert(eltype(G1.z), shift)
	G2 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .- _shift, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	G2.range[5:6] .-= shift
	return G2
end

# ---------------------------------------------------------------------------------------------------
function Base.:*(G1::GMTgrid, G2::GMTgrid)
# Multiply two grids, element by element. Inherit header parameters from G1 grid
	(size(G1.z) != size(G2.z)) && error("Grids have different sizes, so they cannot be multiplied.")
	G3 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .* G2.z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	grd_min_max!(G3)		# Also take care of NaNs
	return G3
end

# ---------------------------------------------------------------------------------------------------
Base.:-(G1::GMTgrid) = -1 * G1
Base.:*(scale::Real, G1::GMTgrid) = Base.:*(G1::GMTgrid, scale::Real)
function Base.:*(G1::GMTgrid, scale::Real)
	_scale = convert(eltype(G1.z), scale)
	G2 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z .* _scale, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	G2.range[5:6] .*= scale
	return G2
end

# ---------------------------------------------------------------------------------------------------
Base.:^(G1::GMTgrid, scale::Int) = Base.:^(G1::GMTgrid, Float64(scale))
function Base.:^(G1::GMTgrid, scale::Real)
	_scale = convert(eltype(G1.z), scale)
	G2 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z.^_scale, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	G2.range[5:6] .^= scale
	return G2
end

# ---------------------------------------------------------------------------------------------------
function Base.:/(G1::GMTgrid, G2::GMTgrid)
# Divide two grids, element by element. Inherit header parameters from G1 grid
	if (size(G1.z) != size(G2.z))  error("Grids have different sizes, so they cannot be divided.")  end
	G3 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z ./ G2.z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	grd_min_max!(G3)		# Also take care of NaNs
	return G3
end

# ---------------------------------------------------------------------------------------------------
function Base.:/(G1::GMTgrid, scale::Real)
	_scale = convert(eltype(G1.z), scale)
	G2 = GMTgrid(G1.proj4, G1.wkt, G1.epsg, G1.range, G1.inc, G1.registration, G1.nodata, G1.title, G1.remark,
				 G1.command, G1.x, G1.y, G1.v, G1.z ./ _scale, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, G1.pad)
	G2.range[5:6] ./= scale
	return G2
end

# ---------------------------------------------------------------------------------------------------
function grd_min_max!(G::GMTgrid)
	# The non-nan version is way faster so use it as a proxy of NaNs and recompute if needed.
	min = minimum(G.z);
	!isnan(min) ? (G.range[5:6] = [min, maximum(G.z)]) : (G.range[5:6] = [minimum_nan(G.z), maximum_nan(G.z)])
	#=
	if !isnan(min)
		G.range[5:6] = [min, maximum(G.z)]
	else
		tic()
		t = G.layout;	G.layout = "TRBa"
		info = grdinfo(G, C=:n, L=0, Vd=1)
		G.range[5:6] = info[1][5:6]
		G.layout = t
		toc()
	end
	=#
end


text:
# Define the new measure data type with <: AbstractMeasureData
struct NewMeasureData <: AbstractMeasureData
    attr1::String # REPLACE WITH ACTUAL ATTRIBUTE
    attr2::DiscreteMeasureData # REPLACE WITH ACTUAL ATTRIBUTE
    # ADD MORE ATTRIBUTES AS NEEDED
    # constructor
    function NewMeasureData(attr1::String, attr2::DiscreteMeasureData)
        # INSERT CHECKS AND OTHER CONSTRUCTOR METHODS HERE
        return new(attr1, attr2) # REPLACE WITH ACTUAL ATTRIBUTES
    end
end

# Extend parameter_refs to return the parameter(s) being measured by a measure using NewMeasureData
function InfiniteOpt.parameter_refs(data::NewMeasureData)::Union{ParameterRef, AbstractArray{<:ParameterRef}}
    return data.attr2.parameter_ref # REPLACE WITH ACTUAL PARAMETER LOCATION
end

# Extend expand_measure to return the finite reformulation of a measure using NewMeasureData
function InfiniteOpt.expand_measure(expr::JuMP.AbstractJuMPScalar,
                                    data::NewMeasureData,
                                    write_model::JuMP.AbstractModel,
                                    point_mapper::Function)::JuMP.AbstractJuMPScalar
    # INSERT APPROPRIATE METHODS HERE USING point_mapper AS APPROPRIATE
    # USING make_point_variable_ref AND make_reduced_variable_ref MAY BE USEFUL
    return expand_measure(expr, data.attr2, write_model, point_mapper) # REPLACE ACTUAL RESULT
end

# Extend supports to return any infinite parameter supports employed by NewMeasureData
# This is only optional if the new abstraction doesn't use supports at all
function InfiniteOpt.supports(data::NewMeasureData)::Vector
    return data.attr2.supports # REPLACE WITH ACTUAL LOCATION
end

# Extend measure_data_in_hold_bounds to determine if NewMeasureData is in the
# domain of hold variable bounds. (Enables hold variable error checking)
function InfiniteOpt.measure_data_in_hold_bounds(data::NewMeasureData,
                                                 bounds::ParameterBounds)::Bool
    # INSERT ACTUAL CHECK HERE
    in_bounds = measure_data_in_hold_bounds(data.attr2, bounds) # REPLACE WITH ACTUAL RESULT
    return in_bounds
end

# Extend measure_name to return the name of a measure using NewMeasureData
# This is optional (uses "measure" otherwise)
function InfiniteOpt.measure_name(data::NewMeasureData)::String
    return data.attr1 # REPLACE WITH ACTUAL NAME LOCATION
end

# Make a convenient measure constructor function for our new measure type
# This should employ measure(expr, data)
function new_measure(expr::JuMP.AbstractJuMPScalar, param::ParameterRef,
                     lb::Number, ub::Number; name::String = "NewMeas",
                     num_supports::Int = 10)::MeasureRef # REPLACE ARGS WITH ACTUAL DESIRED
    # INSERT RELAVENT CHECKS AND OPERATIONS HERE
    # REPLACE BELOW WITH ACTUAL CONSTRUCTION
    attr2 = generate_measure_data(param, num_supports, lb, ub) # just an example
    data = NewMeasureData(name, attr2) # REPLACE WITH ACTUAL
    # built the measure using the built-in constructor
    return measure(expr, data)
end


text:
try
    import Gurobi
    open(joinpath(@__DIR__, "USESOLVER"), "w") do io
        write(io, "Gurobi")
    end
catch
    import ECOS, Cbc
    open(joinpath(@__DIR__, "USESOLVER"), "w") do io
        write(io, "Free")
    end
end


text:
module TestConstructors
	using Base.Test
	using DataArrays
	using DataFrames

	#
	# DataFrame
	#

	df = DataFrame()
	@assert isequal(df.columns, {})
	# TODO: Get this to work
	#@assert isequal(df.colindex, Index())

	df = DataFrame({DataArray(zeros(3)), DataArray(ones(3))}, Index(["x1", "x2"]))
	@assert nrow(df) == 3
	@assert ncol(df) == 2

	# TODO: Make isequal fail if colnames don't match
	@assert isequal(df, DataFrame({DataArray(zeros(3)), DataArray(ones(3))}))
	@assert isequal(df, DataFrame(quote x1 = [0.0, 0.0, 0.0]; x2 = [1.0, 1.0, 1.0] end))

	@assert isequal(df, DataFrame([0.0 1.0; 0.0 1.0; 0.0 1.0], ["x1", "x2"]))
	@assert isequal(df, DataFrame([0.0 1.0; 0.0 1.0; 0.0 1.0]))
	@assert isequal(df, DataFrame(DataArray(zeros(3)), DataArray(ones(3))))

	# TODO: Fill these in
	# From (Associative): ???
	# From (Vector, Vector, Groupings): ???

	@assert isequal(df, DataFrame({"x1" => [0.0, 0.0, 0.0],
		                           "x2" => [1.0, 1.0, 1.0]}))
	@assert isequal(df, DataFrame({"x1" => [0.0, 0.0, 0.0],
		                           "x2" => [1.0, 1.0, 1.0],
		                           "x3" => [2.0, 2.0, 2.0]},
		                          ["x1", "x2"]))

	df = DataFrame(Int, 2, 2)
	@assert size(df) == (2, 2)
	@assert all(types(df) .== {Int, Int})
	@assert all(isna(df))

	df = DataFrame(2, 2)
	@assert size(df) == (2, 2)
	@assert all(types(df) .== {Float64, Float64})
	@assert all(isna(df))

	df = DataFrame({Int, Float64}, ["x1", "x2"], 2)
	@assert size(df) == (2, 2)
	@assert all(types(df) .== {Int, Float64})
	@assert all(isna(df))

	@assert isequal(df, DataFrame({Int, Float64}, 2))
end


text:
using Printf
using Statistics
using Test
using CUDAdrv, CuArrays, CUDAnative

using CuDiff

function derivative_num(f, x::T) where T<:Number
    u = f(x)
    eps = T(1e-3)
    du = (derivative(f, x+eps)[1] - derivative(f, x-eps)[1]) / (2*eps)
    return u, du
end

function kernel(f, d_x, d_a, d_da, d_b, d_db)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    d_a[i], d_da[i] = derivative(f, d_x[i])
    d_b[i], d_db[i] = derivative_num(f, d_x[i])
    #d_b[i], d_db[i] = d_a[i], d_da[i]
    return nothing
end

function test_proc(f)
    x = rand(Float32, 1000)

    a = similar(x)
    da = similar(x)
    b = similar(x)
    db = similar(x)

    d_x = CuArray(x)
    d_a = similar(d_x)
    d_da = similar(d_x)
    d_b = similar(d_x)
    d_db = similar(d_x)

    CUDAnative.@cuda threads=1000 kernel(f, d_x, d_a, d_da, d_b, d_db)

    copy!(a, d_a)
    copy!(da, d_da)
    copy!(b, d_b)
    copy!(db, d_db)

    return x, a, da, b, db
end

function test_all()
    x = sort(rand(Float32, 1000))

    a = similar(x)
    da = similar(x)
    b = similar(x)
    db = similar(x)

    d_x = CuArray(x)
    d_a = similar(d_x)
    d_da = similar(d_x)
    d_b = similar(d_x)
    d_db = similar(d_x)

    function run(f)
        CUDAnative.@cuda threads=1000 kernel(f, d_x, d_a, d_da, d_b, d_db)
        copy!(a, d_a)
        copy!(da, d_da)
        copy!(b, d_b)
        copy!(db, d_db)
        err1 = quantile!(abs.(a .- b), 0.99)
        err2 = quantile!(abs.(da .- db), 0.99)
        passed =  err1 < 1e-5 && err2 < 1e-2
        #Printf.@printf "%.6f\t%.6f\n" err1 err2
        if !passed
            Printf.@printf "%.6f\t%.6f\n" err1 err2
        end
        return passed
    end

    @printf "+:\t%s\n" @test run(x -> x + 2.0)
    @printf "rationals:\t%s\n" @test run(x -> (2*x-0.65f0)/(4.3-3*x))
    @printf "poly:\t%s\n" @test run(x -> CUDAnative.exp(0.04*(x+85)))
    @printf "sin:\t%s\n" @test run(x -> CUDAnative.sin(x))
    @printf "sinpi:\t%s\n" @test run(x -> CUDAnative.sinpi(x))
    # @printf "sin_fast:\t%s\n" @test run(x -> CUDAnative.sin_fast(x))
    @printf "cos:\t%s\n" @test run(x -> CUDAnative.cos(x))
    @printf "cospi:\t%s\n" @test run(x -> CUDAnative.cospi(x))
    #@printf "cos_fast:\t%s\n" @test run(x -> CUDAnative.cos_fast(x))
    @printf "tan:\t%s\n" @test run(x -> -CUDAnative.tan(x))
    #@printf "tan_fast:\t%s\n" @test run(x -> -CUDAnative.tan_fast(x))
    @printf "acos:\t%s\n" @test run(x -> CUDAnative.acos(0.9*x))
    @printf "asin:\t%s\n" @test run(x -> CUDAnative.asin(0.9*x))
    @printf "atan:\t%s\n" @test run(x -> CUDAnative.atan(x))
    @printf "cosh:\t%s\n" @test run(x -> CUDAnative.cosh(x))
    @printf "sinh:\t%s\n" @test run(x -> CUDAnative.sinh(x))
    @printf "tanh:\t%s\n" @test run(x -> CUDAnative.tanh(x))
    @printf "acosh:\t%s\n" @test run(x -> CUDAnative.acosh(1.1+x))
    @printf "asinh:\t%s\n" @test run(x -> CUDAnative.asinh(x))
    @printf "atanh:\t%s\n" @test run(x -> CUDAnative.atanh(0.9*x))
    #@printf "atanh2:\t%s\n" @test run(x -> CUDAnative.atanh2(x,x))
    @printf "exp:\t%s\n" @test run(x -> CUDAnative.exp(x))
    @printf "exp_fast:\t%s\n" @test run(x -> CUDAnative.exp_fast(x))
    @printf "exp2:\t%s\n" @test run(x -> CUDAnative.exp2(x))
    @printf "exp10:\t%s\n" @test run(x -> CUDAnative.exp10(x))
    @printf "expm1:\t%s\n" @test run(x -> CUDAnative.expm1(x))
    @printf "log:\t%s\n" @test run(x -> CUDAnative.log(0.1+x))
    #@printf "log_fast:\t%s\n" @test run(x -> CUDAnative.log_fast(0.1+x))
    @printf "log10:\t%s\n" @test run(x -> CUDAnative.log10(0.1+x))
    #@printf "log10_fast:\t%s\n" @test run(x -> CUDAnative.log10_fast(0.1+x))
    @printf "log2:\t%s\n" @test run(x -> CUDAnative.log2(0.1+x))
    #@printf "log2_fast:\t%s\n" @test run(x -> CUDAnative.log2_fast(0.1+x))
    @printf "log1p:\t%s\n" @test run(x -> CUDAnative.log1p(x+0.1))
    @printf "abs:\t%s\n" @test run(x -> CUDAnative.abs(x-0.5))
    #@printf "saturate:\t%s\n" @test run(x -> CUDAnative.saturate(2*(x-0.5)))
    @printf "sqrt:\t%s\n" @test run(x -> CUDAnative.sqrt(x+0.1))
    @printf "rsqrt:\t%s\n" @test run(x -> CUDAnative.rsqrt(x+0.1))
    @printf "cbrt:\t%s\n" @test run(x -> CUDAnative.cbrt(x+0.1))
    @printf "rcbrt:\t%s\n" @test run(x -> CUDAnative.rcbrt(x+0.1))
    @printf "pow:\t%s\n" @test run(x -> CUDAnative.pow(2+x,1.5f0))
    @printf "erf:\t%s\n" @test run(x -> CUDAnative.erf(10*x-1))
    @printf "erfc:\t%s\n" @test run(x -> CUDAnative.erfc(10*x-1))
    @printf "erfinv:\t%s\n" @test run(x -> CUDAnative.erfinv(0.01+x*0.5))
    @printf "erfcinv:\t%s\n" @test run(x -> CUDAnative.erfcinv(0.2+x*0.5))

    @printf "<:\t%s\n" @test run(x -> x < 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf ">:\t%s\n" @test run(x -> x > 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf "<=:\t%s\n" @test run(x -> x <= 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf ">=:\t%s\n" @test run(x -> x >= 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf "==:\t%s\n" @test run(x -> x == 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf "!=:\t%s\n" @test run(x -> x != 0.5 ? CUDAnative.sin(x) : CUDAnative.cos(x))
    @printf "isapprox:\t%s\n" @test run(x -> isapprox(x, 0.5) ? CUDAnative.sin(x) : CUDAnative.cos(x))
end

test_all()


text:
"""
`SolutionDAE`: Solution of a differential algebraic equation

Contains all fields necessary to store the solution of an DAE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nm`: dimension of the constraint submanifold
* `nt`: number of time steps to store
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ni]` with `q[:,0,:]` the initial conditions
* `λ`:  Lagrange multiplier `λ[nd, nt+1, ni]`
* `ntime`: number of time steps to compute
* `nsave`: store every nsave'th time step (default: 1)
* `nwrite`: save data to disk after every nwrite'th time step (default: ntime)
* `counter`: counter for copied solution entries
* `woffset`: counter for file offset
* `h5`: HDF5 file for storage
"""
mutable struct SolutionDAE{dType, tType, N} <: DeterministicSolution{dType, tType, N}
    nd::Int
    nm::Int
    nt::Int
    ni::Int
    t::TimeSeries{tType}
    q::SDataSeries{dType,N}
    λ::SDataSeries{dType,N}
    ntime::Int
    nsave::Int
    nwrite::Int
    counter::Vector{Int}
    woffset::Int
    h5::HDF5File

    function SolutionDAE{dType, tType, N}(nd, nm, nt, ni, t, q, λ, ntime, nsave, nwrite) where {dType <: Number, tType <: Real, N}
        new(nd, nm, nt, ni, t, q, λ, ntime, nsave, nwrite, zeros(Int, ni), 0)
    end
end

function SolutionDAE(equation::DAE{DT,TT}, Δt::TT, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE; filename=nothing) where {DT,TT}
    @assert nsave > 0
    @assert ntime == 0 || ntime ≥ nsave
    @assert nwrite == 0 || nwrite ≥ nsave
    @assert mod(ntime, nsave) == 0

    if nwrite > 0
        @assert mod(nwrite, nsave) == 0
        @assert mod(ntime, nwrite) == 0
    end

    N  = equation.n > 1 ? 3 : 2
    nd = equation.d
    nm = equation.m
    ni = equation.n
    nt = div(ntime, nsave)
    nt = (nwrite == 0 ? nt : div(nwrite, nsave))
    nw = (nwrite == 0 ? ntime : nwrite)

    @assert nd > 0
    @assert nm > 0
    @assert ni > 0
    @assert nw > 0

    t = TimeSeries{TT}(nt, Δt, nsave)
    q = SDataSeries(DT, nd, nt, ni)
    λ = SDataSeries(DT, nm, nt, ni)
    s = SolutionDAE{DT,TT,N}(nd, nm, nt, ni, t, q, λ, ntime, nsave, nw)
    set_initial_conditions!(s, equation)

    if !isnothing(filename)
        isfile(filename) ? @warn("Overwriting existing HDF5 file.") : nothing
        create_hdf5(s, filename)
    end

    return s
end

function SolutionDAE(t::TimeSeries{TT}, q::SDataSeries{DT,N}, λ::SDataSeries{DT,N}, ntime::Int) where {DT,TT,N}
    @assert q.nd >= λ.nd
    @assert q.nt == λ.nt
    @assert q.ni == λ.ni

    # extract parameters
    nd = q.nd
    nm = λ.nd
    ni = q.ni
    nt = t.n
    ns = div(ntime, nt)

    @assert mod(ntime, nt) == 0

    # create solution
    SolutionDAE{DT,TT,N}(nd, nm, nt, ni, t, q, λ, ntime, ns, 0)
end

function SolutionDAE(file::String)
    # open HDF5 file
    get_config(:verbosity) > 1 ? @info("Reading HDF5 file ", file) : nothing
    h5 = h5open(file, "r")

    # read attributes
    ntime = read(attrs(h5)["ntime"])
    nsave = read(attrs(h5)["nsave"])

    # reading data arrays
    t = TimeSeries(read(h5["t"]), nsave)
    q = SDataSeries(read(h5["q"]))
    λ = SDataSeries(read(h5["λ"]))

    # need to close the file
    close(h5)

    # create solution
    SolutionDAE(t, q, λ, ntime)
end

hdf5(sol::SolutionDAE)  = sol.h5
time(sol::SolutionDAE)  = sol.t.t
ntime(sol::SolutionDAE) = sol.ntime
nsave(sol::SolutionDAE) = sol.nsave
offset(sol::SolutionDAE) = sol.woffset

function set_initial_conditions!(sol::SolutionDAE{DT,TT}, equ::DAE{DT,TT}) where {DT,TT}
    set_initial_conditions!(sol, equ.t₀, equ.q₀, equ.λ₀)
end

function set_initial_conditions!(sol::SolutionDAE{DT,TT}, t₀::TT, q₀::Array{DT}, λ₀::Array{DT}) where {DT,TT}
    set_data!(sol.q, q₀, 0)
    set_data!(sol.λ, λ₀, 0)
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function get_initial_conditions!(sol::SolutionDAE{DT,TT}, asol::AtomicSolutionDAE{DT,TT}, k, n=1) where {DT,TT}
    get_solution!(sol, asol.q, asol.λ, n-1, k)
    asol.t  = sol.t[n-1]
    asol.q̃ .= 0
end

function get_initial_conditions!(sol::SolutionDAE{DT,TT}, q::Vector{DT}, λ::Vector{DT}, k) where {DT,TT}
    get_data!(sol.q, q, 0, k)
    get_data!(sol.λ, λ, 0, k)
end

function CommonFunctions.set_solution!(sol::SolutionDAE, t, q, λ, n, k)
    set_solution!(sol, q, λ, n, k)
end

function CommonFunctions.set_solution!(sol::SolutionDAE{DT,TT}, asol::AtomicSolutionDAE{DT,TT}, n, k) where {DT,TT}
    set_solution!(sol, asol.t, asol.q, asol.λ, n, k)
end

function CommonFunctions.set_solution!(sol::SolutionDAE{DT,TT}, q::Vector{DT}, λ::Vector{DT}, n, k) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.λ, λ, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function CommonFunctions.reset!(sol::SolutionDAE)
    reset!(sol.q)
    reset!(sol.λ)
    compute_timeseries!(sol.t, sol.t[end])
    sol.counter .= 1
    sol.woffset += sol.nt
end

function Base.close(solution::SolutionDAE)
    # close(solution.t)
    # close(solution.q)
    # close(solution.λ)
    close(solution.h5)
end

"Creates HDF5 file and initialises datasets for DAE solution object."
function create_hdf5(solution::SolutionDAE{DT,TT,2}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution)

    # create datasets
    nt = div(solution.ntime, solution.nsave)
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, nt+1), "chunk", (solution.nd,1))
    λ = d_create(solution.h5, "λ", datatype(DT), dataspace(solution.nm, nt+1), "chunk", (solution.nm,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1] = solution.q[:, 0]
    λ[:, 1] = solution.λ[:, 0]

    return solution.h5
end

"Creates HDF5 file and initialises datasets for DAE solution object."
function create_hdf5(solution::SolutionDAE{DT,TT,3}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution)

    # create datasets
    nt = div(solution.ntime, solution.nsave)
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ni), "chunk", (solution.nd,1,1))
    λ = d_create(solution.h5, "λ", datatype(DT), dataspace(solution.nm, nt+1, solution.ni), "chunk", (solution.nm,1,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1, :] = solution.q[:, 0, :]
    λ[:, 1, :] = solution.λ[:, 0, :]

    return solution.h5
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionDAE{DT,TT,2}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # copy data from solution to HDF5 dataset
    h5["q"][:, j1:j2] = solution.q[:, 1:n]
    h5["λ"][:, j1:j2] = solution.λ[:, 1:n]

    return nothing
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionDAE{DT,TT,3}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # copy data from solution to HDF5 dataset
    h5["q"][:, j1:j2, :] = solution.q[:, 1:n, :]
    h5["λ"][:, j1:j2, :] = solution.λ[:, 1:n, :]

    return nothing
end


text:
using DataFrames, DataFramesMeta, Mimi

include("../src/main.jl")

df_diff = DataFrame([UTF8String, Float64], [:name, :maxdiff],0)
for c in Mimi.components(m), v in Mimi.variables(m, c)
    filename = "csharp_output/$c.$v.csv"
    if isfile(filename)
        t = typeof(m[c,v])
        if t==Float64
            s = readall(filename)
            csharp_val = parse(Float64, s)
            diff = csharp_val - m[c,v]
            push!(df_diff, (string(v), diff))
        elseif t==Array{Float64,1}
            csharp_val = vec(readdlm(filename))
            jl_val = m[c,v]
            if length(jl_val)!=16
                jl_val = jl_val[1:end-1]
            end
            diff = csharp_val .- jl_val
            push!(df_diff, (string(v), maxabs(diff)))
        elseif t==Array{Float64,2}
            csharp_val = readcsv(filename)
            jl_val = m[c,v]
            if size(jl_val, 1)!=16
                jl_val = jl_val = m[c,v][1:end-1,:]
            end
            diff = csharp_val .- jl_val
            push!(df_diff, (string(v), maxabs(diff)))
        else
            println(c, ".", v)
        end
    else
        println(c, ".", v)
    end
end

@where(df_diff, :maxdiff .> 0.00000000001)


text:
struct InducingPoints{M,O<:Union{Optimizer,Nothing}}
    Z::M
    opt::O
end

Base.size(Z::InducingPoints) = size(Z.Z)
Base.size(Z::InducingPoints,i::Int) = size(Z.Z,i)
Base.convert(::Type{<:AbstractMatrix},Z::InducingPoints) = Z.Z


text:
__precompile__()
module Coinigy

  export  CoinigyHandler,
          getUserInfo,
          getActivityLog,
          getNotifications,
          getAccounts,
          getBalances,
          getBalanceHistory,
          getOrders,
          getOrderTypes,
          getAlerts,
          getUserWatchList,
          getNewsFeed,
          getExchanges,
          getMarkets,
          getTradeHistory,
          getAsks,
          getBids,
          getOrders,
          getMarketData,
          getTicker,
          activateTradingKey,
          addOrder

  include("client.jl")

end


text:
#############################################
###### Structs ##############################
#############################################
struct Reaction ## Reaction Object
	index::Int64 # List index of reaction
	reactants::Array{Int64} # List of reactant indices
	react_coef::Array{Int64} # List of reactant coefficients
	products::Array{Int64} # List of product indices
	prod_coef::Array{Int64} # List of product coefficients
	catalysts::Array{Int64} # Catalyst indicies
	cat_coef::Array{Float64} # Catalysts coefficients
	propensity::String # Propensity identifier
	rate_constant::Float64 # Rate constant
end

struct CRS # Chemical reaction system object
	molecule_list::Array{String} # List of molcule strings
	molecule_dict::Dict{String, Int64} # Dict maps molecule strings to index
	reaction_list::Array{Reaction} # List of reaction objects
	parameters::Dict{Any,Any} # Every other parameter
end

#############################################
###### CRS functions ########################
#############################################
function generate_all_binary_seqs(max_L)
    ### Generate all possible sequences
    all_seqs = []
    for j in 1:max_L
        comps = Array{String}(undef,j + 1)
        for l in 0:j
            comps[l+1] = "B"^l * "A"^(j-l)
        end

        for c in comps
            p = unique(collect(permutations(c)))
            all_seqs = vcat(p,all_seqs)
        end
    end
    seqs = [join(s) for s in all_seqs]
    sort!(seqs, by = length, rev= false)
    return seqs
end

### Count occurance of t in s Function ##
function countall(s, t)
    n = length(collect(eachmatch(Regex(t), s)))
    return n
end

function all_splits(seq)
    l = length(seq)
    left = Array{String}(undef,(l-1))
    right =Array{String}(undef,(l-1))

    for i in 1:(l-1)
        left[i] = seq[1:i]
        right[i] = seq[(i+1):end]
    end
    return left,right
end

function calculate_reduced_mass(m1_string, m2_string, mA,mB)
    m1_mass = countall(m1_string,"A")*mA + countall(m1_string,"B")*mB
    m2_mass = countall(m2_string,"A")*mA + countall(m2_string,"B")*mB
    reduce = reduced_mass(m1_mass,m2_mass)
    return reduce
end

function reduced_mass(m1,m2)
    return (m1*m2)/(m1+m2)
end

function bimolecular_coef(k, T, R, V,m)
    return (k/V)*sqrt(R*T/m)
end

function generate_binary_ligation_CRS(max_L, kf, kb, volume::Float64 = 1.0, T::Float64 = 300, R::Float64= 8.3144598, mA::Int64= 100, mB::Int64 = 100)
    seqs = generate_all_binary_seqs(max_L)
    ## Generate Molecule List and Dictionary
    molecule_dict = Dict{String,Int64}()
    for (n, m) in enumerate(seqs)
        molecule_dict[m] = n
    end

    ## Generate Reaction List
    Aindex = molecule_dict["A"]
    Bindex = molecule_dict["B"]

    ### Generate all forward and backward reactions
    n_rxn = sum([(i-1)*2^(i +1) for i in 1:max_L  ])
    rxn_list = Array{Reaction}(undef,n_rxn )
    rID = 1
    for s in seqs
        ## Get all splits of the sequence
        left_side, right_side = all_splits(s)
        for i in 1:length(left_side)
            ### Forward Reaction
            ## Check if both sides of sequence are identical
            if left_side[i] != right_side[i]
                reactants = [ molecule_dict[ left_side[i] ], molecule_dict[ right_side[i] ] ]
                reactant_coef = [1,1]
            else
                reactants = [molecule_dict[ left_side[i] ] ]
                reactant_coef = [2]
            end

            products = [ molecule_dict[s] ]
            product_coef = [1]

            m_red = calculate_reduced_mass(left_side[i], right_side[i],mA, mB)
            k = bimolecular_coef(kf, T, R, volume, m_red)

            forward_rxn = Reaction(rID, reactants, reactant_coef, products, product_coef, Array{Int64}(undef,0), Array{Float64}(undef,0), "STD", k)
            rxn_list[rID] = forward_rxn
            rID += 1

            ## Backward Reaction
            backward_rxn = Reaction(rID, products, product_coef,reactants, reactant_coef, Array{Int64}(undef,0), Array{Float64}(undef,0), "STD", kb)
            rxn_list[rID] = backward_rxn
            rID += 1

        end

    end

    parameters = Dict(:max_L => max_L, :kf => kf, :kb => kb, :volume => volume, :T => T, :R => R, :mA=> mA, :mB => mB)
    return(CRS(seqs, molecule_dict, rxn_list, parameters))
end

function copy_CRS(original_CRS, volume_multiplier)
    reaction_list = similar(original_CRS.reaction_list)
    for rID in 1:length(reaction_list)
        rxn = original_CRS.reaction_list[rID]
        if sum(rxn.react_coef) ==1
            reaction_list[rID] = deepcopy(original_CRS.reaction_list[rID])
        elseif sum(rxn.react_coef) == 2
            new_rxn = Reaction(rxn.index, rxn.reactants, rxn.react_coef, rxn.products, rxn.prod_coef, rxn.catalysts, rxn.cat_coef, rxn.propensity, (1.0/volume_multiplier)*rxn.rate_constant)
            reaction_list[rID] = new_rxn
        end

    end
    new_parameters = original_CRS.parameters
    new_parameters[:volume] *= volume_multiplier

    new_CRS = CRS(original_CRS.molecule_list, original_CRS.molecule_dict, reaction_list, new_parameters )
    return new_CRS
end


text:
using MPI

using ClimateMachine: Settings

using ClimateMachine.Mesh.Topologies: StackedBrickTopology

import ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid

#####
##### RectangularDomain
#####

struct RectangularDomain{FT} <: AbstractDomain
    Np::Int
    Ne::NamedTuple{(:x, :y, :z), NTuple{3, Int}}
    L::NamedTuple{(:x, :y, :z), NTuple{3, FT}}
    x::NTuple{2, FT}
    y::NTuple{2, FT}
    z::NTuple{2, FT}
    periodicity::NamedTuple{(:x, :y, :z), NTuple{3, Bool}}
end

Base.eltype(::RectangularDomain{FT}) where {FT} = FT

function Base.show(io::IO, domain::RectangularDomain{FT}) where {FT}
    Np = domain.Np
    Ne = domain.Ne
    L = domain.L

    first = "RectangularDomain{$FT}\n"
    second = "    Np = $Np, Ne = $Ne\n"
    third = "    L = $L\n"

    return print(io, first, second, third)
end

name_it(Ne::NamedTuple{(:x, :y, :z)}) = Ne
name_it(Ne) = (x = Ne[1], y = Ne[2], z = Ne[3])

"""
    RectangularDomain(FT=Float64;
                      Ne,
                      Np,
                      x = (-1, 1),
                      y = (-1, 1),
                      z = (-1, 1),
                      periodicity = (true, true, false))

Returns a `RectangularDomain` representing the product of `x, y, z` intervals,
specified by 2-tuples.

The `RectangularDomain` is meshed with a simple `DiscontinuousSpectralElementGrid`
with an isotropic polynomial order `Np` and a 3-tuple of `Ne`lements
giving the number of elements in `x, y, z`.

Additional arguments are:

- `periodicity`: a 3-tuple that indicates periodic dimensions with `true`,

- `boundary`: specifies the boundary condition on each boundary with a
              boundary condition `tag`

- `array_type`: either `Array` for CPU computations or `CuArray` for
                GPU computations

- `mpicomm`: communicator for sending data across nodes in a distributed memory
             configuration using the Message Passing Interface (MPI).
             See https://pages.tacc.utexas.edu/~eijkhout/pcse/html/mpi-comm.html

Example
=======

```julia
julia> using ClimateMachine; ClimateMachine.init()

julia> using ClimateMachine.Ocean.RectangularDomains: RectangularDomain

julia> domain = RectangularDomain(Ne=(7, 8, 9), Np=4, x=(0, 1), y=(0, 1), z=(0, 1))
RectangularDomain{Float64}:
    Np = 4, Ne = (x = 7, y = 8, z = 9)
    L = (x = 1.00e+00, y = 1.00e+00, z = 1.00e+00)
    x = (0.00e+00, 1.00e+00), y = (0.00e+00, 1.00e+00), z = (1.00e+00, 0.00e+00)
```
"""
function RectangularDomain(
    FT = Float64;
    Ne,
    Np,
    x::Tuple{<:Number, <:Number},
    y::Tuple{<:Number, <:Number},
    z::Tuple{<:Number, <:Number},
    periodicity = (true, true, false),
)

    Ne = name_it(Ne)
    periodicity = name_it(periodicity)

    west, east = FT.(x)
    south, north = FT.(y)
    bottom, top = FT.(z)

    east > west || error("Domain x-limits must be increasing!")
    north > south || error("Domain y-limits must be increasing!")
    top > bottom || error("Domain z-limits must be increasing!")

    L = (x = east - west, y = north - south, z = top - bottom)

    return RectangularDomain(
        Np,
        Ne,
        L,
        (west, east),
        (south, north),
        (bottom, top),
        periodicity,
    )
end

array_type(domain::RectangularDomain) = Settings.array_type
eltype(::RectangularDomain{FT}) where {FT} = FT

function DiscontinuousSpectralElementGrid(
    domain::RectangularDomain{FT};
    boundary_tags = ((0, 0), (0, 0), (1, 2)),
    array_type = Settings.array_type,
    mpicomm = MPI.COMM_WORLD,
) where {FT}

    west, east = domain.x
    south, north = domain.y
    bottom, top = domain.z

    element_coordinates = (
        range(west, east, length = domain.Ne.x + 1),
        range(south, north, length = domain.Ne.y + 1),
        range(bottom, top, length = domain.Ne.z + 1),
    )

    topology = StackedBrickTopology(
        mpicomm,
        element_coordinates;
        periodicity = tuple(domain.periodicity...),
        boundary = boundary_tags,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = domain.Np,
    )

    return grid
end

max_forks_repo_name: jm-c/CLIMA


text:
function constructportfolios(;
  refreshcomp::Bool = true,
  refreshtretcrsp::Bool = true,
  refreshtvolmcapcrsp::Bool = true,
  refreshdataseries::Bool = true,
  refreshportfolios::Bool = true,
  portname::String = PORT_NAME,
  panelname::String = PANEL_NAME,
  datapath::String=DATA_PATH,
  injlsstream::Function = IN_JLS_STREAM,
  outjlsstream::Function = OUT_JLS_STREAM,
  outpath::String = OUT_PATH,
  trialrunonly::Bool = false)::DataFrame

  local panel::DataFrame
  local port::DataFrame

  #load the leverage data set
  local dailyunlesstrial::Symbol = trialrunonly ? :monthly : :daily
  local cqports::CQuartileConstructions = leverageconstructions()
  panelname = "$panelname-$dailyunlesstrial"
  portname = "$portname-$dailyunlesstrial"

  if refreshportfolios
    local panel = formdataseries(
      refreshcomp = refreshcomp,
      refreshtvolmcapcrsp = refreshtvolmcapcrsp,
      refreshtretcrsp = refreshtretcrsp,
      refreshdataseries=refreshdataseries,
      trialrunonly = trialrunonly)

    panel = makepanelmonthly(panel, crsptype=dailyunlesstrial)
    port = formfactors(panel, crsptype=dailyunlesstrial)
    port = mergeoutside(port)
    println(describe(port))

    #merge and verify integrity
    println("Size of panel: $(size(panel))")
    #panel = join(panel, port, on=:date, kind=:inner)
    checkdfkey(panel, [:date, :lpermno])
    println("Size of panel: $(size(panel))")

    #WARNING These rows save space, but are not essential
    retainedpanelcols = unique([:gvkey; :lpermno; :date; :linkeffdate; :linkenddate;
      :adate; :begindate; :enddate; names(port); Fvals(cqport);
      :ret; :lret; T1_CONTROL10;])
    #select!(panel, retainedpanelcols)

    checksum::Float64 = computechecksum(panel) + computechecksum(port)
    @info "Panel constructed with checksum $(checksum)"

    if TEST_OUTPUT
      d = Date(2015,12,31)
      CSV.write("output\\$(panelname)_$d.csv", panel[(cd->cd==d).(panel.date), :])
      CSV.write("output\\$portname.csv", port[:,
        filter(s::Symbol->string(s)[1:2] ≠ "S_", names(port))])
    end

    outjlsstream("$datapath\\$panelname.jls.lz4") do s
      serialize(s, panel)
    end
    CSV.write("$datapath\\$portname.csv")
  #error("A-OK!")
  else
    injlsstream("$datapath\\$panelname.jls.lz4") do s
      panel = deserialize(s) #NOTE: DEBUG ONLY
    end
  end

  return panel
end

#computes a checksum of a dataframe. Used to check consistency.
function computechecksum(df::DataFrame; checksum::Float64 = 0.0)::Float64
  for c ∈ eachcol(df)
    if eltype(c) <: Union{Missing, Real}
      checksum += sum((x->ismissing(x) ? 0.0 : isfinite(x) ? x : 0.0).(c))
    end
  end
  checksum = log(checksum)

  return checksum
end

#this set of functions computes the returns from the panel
#accomplished via dispath
const NandRet = Tuple{Int, MFloat64}

function NandRet(panel::AbstractDataFrame, weighton::Symbol; Fret::Symbol=FRET)::NandRet
  spanel::SubDataFrame = view(
    panel, completecases(spanel[!, [Fret, Fweighton]]), [Fret, Fweighton])

  N::Int = size(spanel,1)
  (N==0) && return (0,missing)
  return NandRet(N, spanel[!,Fret], spanel[!,Fweighton])
end

function NandRet(panel::AbstractDataFrame, ::Nothing; Fret::Symbol=FRET)::NandRet

  rets::T where T<:SubArray = view(panel[!, Fret], (!ismissing).(panel[!, Fret]))
  N::Int = length(rets)
  (N==0) && (return (0,missing))
  return NandRet(N, rets)
end

NandRet(N::Int, rets::AbstractVector)::NandRet = (N, sum(rets)/N)
NandRet(N::Int, rets::AbstractVector, weights::AbstractVector
  )::NandRet = (N, sum(rets .* weights)/sum(weights))

#gets a flat return file
#assumes the existance of a return file
function formfactors(panel::DataFrame; Fret::Symbol=FRET, crsptype::Symbol)::DataFrame

  local port::DataFrame
  local cqports::CQuartileConstructions = leverageconstructions()
  local fgrps::Vector{Symbol} = Fgroups(cqports)
  local debugsuffix="$crsptype"

  #get the monthly returns
  #print("PROFILING: make panel monthly - ")

  #this will hold all of the portfolios
  #get the returns
  #print("PROFILING: unstack returns - ")
  #port = unstack(view(panel, :, [:permno, :date, Fret]), :date, :permno,
  #  Fret, renamecols=(s->Symbol("S_", s)))
  port = DataFrame(date = unique(panel.date))
  sort!(port, :date)

  Nport::Int = size(port, 1)

  #form the necessary columns
  for pc ∈ cqports
    port[!, pc.name] = missings(MFloat64, Nport)
    port[!, pc.Fn] = missings(MInt, Nport)
  end
  #index the return dates
  portrows::Dict = Dict(r.date=>r for r ∈ eachrow(port))

  #good point to check our answers
  if TEST_OUTPUT
    d = Date(2015,12,31)
    CSV.write("output\\panel$(debugsuffix)_$d-mid.csv", panel[(cd->cd==d).(panel.date), :])
  end

  #form the portfolio returns
  #print("PROFILING: portfolio returns loop - ")
  spanels::GroupedDataFrame = groupby(panel, :date)
  @mpar for i ∈ 1:length(spanels)
    spanel::SubDataFrame = spanels[i]
    r::DataFrameRow = portrows[spanel.date[1]]

    #iterate through our portfolio constructions
    for pc ∈ cqports
      r[pc.name] = 0.0
      r[pc.Fn] = 0
      for grp::Int ∈ 1:pc.Ngroups #add the contribution of each group to the factor return
        sspanel = view(spanel, (g::MInt->(!ismissing(g)) && (g == grp)).(spanel[!,pc.Fgroup]),:)

        local N::Int
        local val::MFloat64
        (N, val) = NandRet(sspanel, pc.Fweighton)
        r[pc.Fn] += N
        r[pc.name] += val * pc.groupweights[grp]
      end
    end
  end

  #flatten the quartiles too
  #print("PROFILING: quartile flattening and joining - ")
  #=for pc ∈ cqports
    local portQ::DataFrame = unstack(view(panel, :, [:permno, :date, pc.Fgroup]), :date, :permno,
      pc.Fgroup, renamecols=(s->Symbol("Q_", pc.Fgroup, "_", s)))
    port = join(port, portQ, on=:date)
  end=#

  return port
end

#assigns the group for a portfolio
function assigngroups!(vals::AbstractVector{T},
  groups::AbstractVector{MInt},
  breakpoints::NTuple) where T

  N::Int = length(vals)
  quantiles::Vector{MFloat64} = Vector{MFloat64}(undef, N)

  #compute the quantiles
  qdist = StatsBase.ecdf((val::T->val).(vals))
  @inbounds @simd for i ∈ 1:N
    f::Float64 = vals[i]
    quantiles[i] = qdist(f)
  end

  #compute the groups
  groups .= (q::Float64-> findfirst(b::Float64->
      q≤b, breakpoints)).(quantiles)

  return nothing
end

#assigns the group for a portfolio
function assigngroups!(vals::AbstractVector{T},
  control::AbstractVector{V},
  groups::AbstractVector{MInt},
  breakpoints::NTuple) where {T,V}

  N::Int = length(vals)
  cgroups::Vector{MInt} = deepcopy(groups)

  assigngroups!(control, cgroups, breakpoints)

  for i ∈ 1:length(breakpoints)
    local svals::SubArray{T} = view(vals, cgroups .== i)
    local sgroups::SubArray{MInt} = view(groups, cgroups .== i)
    assigngroups!(svals, sgroups, breakpoints)
  end

  return nothing
end

#use this to create the default portoflio constructions
function leverageconstructions()
  local cqports::CQuartileConstructions
  local cqportsalt::CQuartileConstructions

  local cqports = CQuartileConstructions(CQUART_FVALS, CQUART_FCONTROLS,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES)
  local cqportsalt = CQuartileConstructions(CQUART_FVALS, CQUART_FCONTROLS_ALT,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES_ALT)

  cqports = merge(cqports, cqportsalt)

  return cqports
end

function constructportfolios!(lev::DataFrame)
  #first construct the portfolios
  local cqports::CQuartileConstructions = leverageconstructions()

  N::Int = size(lev,1)

  for pc ∈ cqports
    lev[!, pc.Fgroup] = Vector{MInt}(undef, N)
    #lev[!, pc.Fweight] = Vector{MFloat64}(undef, N)
  end

  #assign the gorups and weights
  slevs = groupby(lev, :fyear)
  @mpar for i ∈ 1:length(slevs)
    local slev::SubDataFrame = slevs[i]
    for pc ∈ cqports #for each portfolio construction

      #don't want to worry about missings here
      local essentialfields::Vector{Symbol} = filter!(
        s->!isnothing(s), [pc.Fval, pc.Fcontrol, pc.Fweighton])
      local sslev::SubDataFrame = view(slev,completecases(slev[:,essentialfields]), :)
      if size(sslev,1)>length(pc.breakpoints) #check that we have enough complete cases

        #assign the groups
        assigngroups!(sslev[!, pc.Fval], sslev[!, pc.Fcontrol], sslev[!, pc.Fgroup], pc.breakpoints)

        #=now assign the weights
        for s3lev ∈ groupby(sslev, pc.Fgroup)
          assignweightswithin!(s3lev, pc.Fweight)
        end=#
      end
    end
  end

end

#=function assignportfolios()
  returncontrolled::CQuartileConstructions = CQuartileConstructions(
    CQUART_FVALS, CQUART_FCONTROLS,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES)
  assetcontrolled::CQuartileConstructions = CQuartileConstructions(
    CQUART_FVALS, CQUART_FCONTROLS_ALT,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES_ALT)
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("MicrobiomeTests.jl")
MicrobiomeTests.runtests()

using Microbiome
Microbiome.runtests()


text:
@testset "Covariance Matrix for GuaranteedEllipseFit" begin
    @info "Test: Guaranteed Ellipse Fit Covariance Matrix"

    @testset "Numerical" begin
        𝚺_ref = [9.26328248166718          -2.0798069234395        -0.583656167803226        0.0481226621171429        0.0018586684115733;
                -2.07980692343951            6.006823008749        0.0378537156164336       0.00295257680184317     -0.000656827363252417;
                -0.583656167803228        0.0378537156164377          7.37460418632127         0.158598345149172      0.000644886416916763;
                 0.0481226621171417       0.00295257680184355         0.158598345149172          3.72310516565684       -0.0007869916229663;
                 0.00185866841157333     -0.000656827363252431      0.000644886416916745     -0.000786991622966293      0.000179354045153511]

        ℳ = GuaranteedEllipseFitting.svectors(transpose(data), Val{2}())
        observations = Observations(tuple(ℳ))

        𝛉₀ = fit_ellipse(observations, DirectEllipseFit())
        𝛉₁ = fit_ellipse(observations, GuaranteedEllipseFit(𝛉₀))
        alg_to_geo = AlgebraicToGeometric()
        𝛈 = alg_to_geo(𝛉₁)
        𝚲 = determine_algebraic_covariance(observations, 𝛉₁, GuaranteedEllipseFit())
        𝚺 = determine_geometric_covariance(observations, 𝛈, GuaranteedEllipseFit())
        @test isapprox(norm(𝚺_ref - 𝚺), 0.0, atol= 1e-5)
    end

end


text:
@testset "GPCP coorShift" begin
	# convert from -180to180 to "0to360"
	lon = collect(-179.5:1:179.5);
	datain = repeat(-179.5:1:179.5,1,180);
	lonout,dataout = coorShift_lon(lon,datain,to="0to360")
	@test size(datain) == size(dataout)
	@test lonout == collect(0.5:1:359.5)
	@test dataout[1,1] == 0.5
	@test dataout[end,1] == -0.5
	@test dataout[181,2] == -179.5

	# convert from 0to360 to "-180to180" + 3D
	lon = collect(0.125:0.25:359.875);
	datain = Array{Float64}(undef,720,length(lon),2);
	datain[:,:,1] = transpose(repeat(0.125:0.25:359.875,1,720));
	datain[:,:,2] = transpose(repeat(0.125:0.25:359.875,1,720));
	lonout,dataout = coorShift_lon(lon,datain,to="-180to180")
	@test size(datain) == size(dataout)
	@test lonout == collect(-179.875:0.25:179.875)
	@test dataout[1,1,1] == 180.125
	@test dataout[end,1,2] == 180.125
	@test dataout[2,181,1] == dataout[2,181,2]
	@test dataout[2,end,1] == 179.875
end

@testset "GPCP read (head, lon/lat, time, data)" begin
	header_read = readgpcpd_head(joinpath(dirname(@__DIR__),"test","input","gpcpd_data"));
	@test length(header_read) == 1440
	lon,lat = readgpcpd_lonlat(joinpath(dirname(@__DIR__),"test","input","gpcpd_data"))
	@test lon == collect(0.5:1:359.5)
	@test lat == collect(89.5:-1:-89.5)
	timeout = readgpcpd_time(joinpath(dirname(@__DIR__),"test","input","gpcpd_data"));
	@test timeout == collect(DateTime(2011,6,1):Dates.Day(1):DateTime(2011,6,30))
	dataout = readgpcpd(joinpath(dirname(@__DIR__),"test","input","gpcpd_data"));
	@test size(dataout) == (360,180,30)
end