text:
using XPORTA, Polyhedra

"""
    rotate_facet(F::HalfSpace, G::HalfSpace xbar::Vector{Int})

Performs a rotation of facet `F` relative to non-included vertex `xbar`. Returns
the rotated `HalfSpace`. `F` is a polytope facet, `G` is a facet of `F` and `xbar`
is a vertex not contained on the boundary of `F`.
"""
function rotate_facet(F::HalfSpace, G::HalfSpace, xbar::Vector{Int}) :: HalfSpace
    a_rot = (F.β - F.a'*xbar)*G.a - (G.β - G.a'*xbar)*F.a
    β_rot = (F.β - F.a'*xbar)*G.β - (G.β - G.a'*xbar)*F.β

    HalfSpace(a_rot, β_rot)
end

function adjacent_facets(facet, vertices)

end

function adjacency_decomposition(facet, vertices)

end


text:
struct Posterior{Tp, Tl}
    name::Symbol                    # Model name
    params_names::Vector{Symbol}    # Names of hyperparameters
    latent_names::Vector{Symbol}    # Names of latent variables
    hyperparameters::Tp             # Posterior distribution of hyperparameters
    latent::Tl                      # Posterior distribution of latent variables
    loglikelihood_func::Function    # Loglikelihood function
    nsamples::Int                   # Number of posterior samples
    nparams::Int                    # Number of hyperparameters
    ndata::Int                      # Number of data points

    Posterior(name, pnames, lnames, hp::Tp, latent::Tl, l_fun, ns, np, nd) where {Tp, Tl} = new{Tp, Tl}(
        name,
        pnames,
        lnames,
        hp,
        latent,
        l_fun,
        ns,
        np,
        nd
    )
end

# Constructor from chains
function posterior(mdl, chn::Chains, data)
    npc, _, nc = size(chn)
    nsmp = npc * nc
    ndata = data_size(data)

    # Extract parameters
    params_posterior, latent_posterior = extract_posterior_samples(mdl, chn)
    params_names = collect(keys(params_posterior[1]))
    latent_names = collect(keys(latent_posterior[1]))
    # Compute loglikelihoods for each sample at each data point
    l_fun(x) = loglikelihood(mdl, x, data)

    return Posterior(mdl.name,
                    params_names,
                    latent_names,
                    params_posterior,
                    latent_posterior,
                    l_fun,
                    nsmp,
                    length(params_names),
                    ndata
    )
end

# Show
function Base.show(io::IO, ::MIME"text/plain", P::AnimalBehavior.Posterior)
    table_conf = set_pt_conf(tf = tf_markdown, alignment = :c)
    println(io, "Posterior probability for ", P.name, " with $(P.nsamples) samples, from $(P.ndata) data points.")

    println(io)
    pretty_table_with_conf(table_conf,
        collect(values(expectation(P).hyperparameters))';
        header = P.params_names,
        title = "Expected hyperparameters values")

    println(io)
    DIC, pD = dic(P)
    WAIC, pW = waic(P)
    fit_header = ["", "Goodness of fit", "Complexity"]
    fit_vals = ["DIC" DIC pD;
                "WAIC" WAIC pW;
                "BIC" bic(P) P.nparams * log(P.ndata)]
    pretty_table_with_conf(table_conf, fit_vals; header=fit_header)
end

ext: jl

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
using Test
using DataDeps
using DataDepsGenerators
using ReferenceTests

@testset "CKAN Demo test" begin
    registration_code = generate(CKAN(), "https://demo.ckan.org/dataset/gold-prices")

    @testset "Integration Test" begin
        eval(Meta.parse(registration_code))
        @test length(collect(readdir(datadep"Gold Prices in London 1950-2008 (Monthly)"))) > 0
    end

    @test_reference "references/CKAN Gold Prices.txt" registration_code
end

@testset "data.gov test" begin
    @test_reference "references/CKAN Consumer Complaint.txt" generate(CKAN(), "https://catalog.data.gov/api/3/action/package_show?id=consumer-complaint-database")
end

@testset "data.gov.au test" begin
    @test_reference "references/CKAN Hazelwood Fire.txt" generate(CKAN(), "https://data.gov.au/api/3/action/package_show?id=2016-soe-atmosphere-hourly-co-and-24h-pm2-5-concentrations-measured-during-the-hazelwood-mine-fire")
end


text:
module FastIOBuffersBenchmarks

using Random
using BenchmarkTools
using FastIOBuffers

const N = 1000

suite = BenchmarkGroup()

suite["write"] = BenchmarkGroup()
suite["write"]["Float64"] = @benchmarkable write(buf, x) evals = N setup = begin
    x = rand(Float64)
    buf = FastWriteBuffer(Vector{UInt8}(undef, N * Core.sizeof(x)))
end
suite["write"]["String"] = @benchmarkable write(buf, x) evals = N setup = begin
    x = randstring(8)
    buf = FastWriteBuffer(Vector{UInt8}(undef, N * sizeof(x)))
end

suite["read"] = BenchmarkGroup()
suite["read"]["Float64"] = @benchmarkable read(buf, Float64) evals = N setup = begin
    rng = MersenneTwister(1)
    writebuf = IOBuffer()
    map(1 : N) do _
        write(writebuf, rand(rng, Float64))
    end
    buf = FastReadBuffer(take!(writebuf))
end
suite["read"]["String"] = @benchmarkable (seekstart(buf); read(buf, String)) setup = begin
    rng = MersenneTwister(1)
    writebuf = IOBuffer()
    write(writebuf, randstring(rng, N))
    buf = FastReadBuffer(take!(writebuf))
end

overhead = BenchmarkTools.estimate_overhead()
results = run(suite, verbose=true, overhead=overhead, gctrial=false)

for result in results["write"]
    println("$(first(result)):")
    display(last(result))
    println()
end
for result in results["read"]
    println("$(first(result)):")
    display(last(result))
    println()
end

end


text:
@everywhere include("../RandomWalk.jl")

numtraj = 1000000
time = vcat(0, [2^(i/3) for i=-15:50])
Nvals = [2, 4, 6, 10, 12, 15, 20, 24, 30, 40, 50, 60, 80, 100]

for N in Nvals
    println("Begin N = $N")
    path = "/data/niels/sim/rw/dom/N_$N"
    mkpath(path)
    traj = RandomWalk.dom_random_walk(numtraj, 2*N, 70000, time)
    RandomWalk.savedata("$path/traj.h5", time, traj, 2*N)
end

max_stars_repo_name: nielsmde/RandomWalk

max_issues_repo_name: nielsmde/RandomWalk

max_forks_repo_name: nielsmde/RandomWalk


text:
module UMAP

using Distances: evaluate, pairwise!, Euclidean, SqEuclidean, SemiMetric
using NearestNeighborDescent: DescentGraph
using LsqFit: curve_fit
using SparseArrays: SparseMatrixCSC, sparse, dropzeros, nzrange, rowvals, nonzeros
using LinearAlgebra: Symmetric, Diagonal, issymmetric, I
using Arpack: eigs

include("utils.jl")
include("umap_.jl")

export umap, UMAP_

end # module


text:
"""
The **pcl_io** library contains classes and functions for reading and writing
point cloud data (PCD) files, as well as capturing point clouds from a variety
of sensing devices

## Exports

$(EXPORTS)
"""
module PCLIO

export loadPCDFile, loadOBJFile, loadPLYFile, load, savePCDFile,
    saveOBJFile, savePLYFile

using DocStringExtensions
using LibPCL
using PCLCommon
using Cxx

libpcl_io = LibPCL.find_library_e("libpcl_io")
try
    Libdl.dlopen(libpcl_io, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

cxx"""
#include <pcl/io/pcd_io.h>
#include <pcl/io/ply_io.h>
#include <pcl/io/obj_io.h>
#include <pcl/io/auto_io.h>
"""

for f in [
        :loadPCDFile,
        :loadOBJFile,
        :loadPLYFile,
        :load,
        ]
    ex = Expr(:macrocall, Symbol("@icxx_str"),
        "pcl::io::$f(\$(pointer(s)), *\$(cloud.handle));")
    @eval begin
        function $f{T}(s::AbstractString, cloud::PointCloud{T})
            ret = $ex
            if ret != 0
                error("failed to $f: code $ret")
            end
            ret
        end
    end
end

for f in [
        :savePCDFile,
        :saveOBJFile,
        :savePLYFile,
        ]
    ex = Expr(:macrocall, Symbol("@icxx_str"),
        "pcl::io::$f(\$(pointer(s)), *\$(cloud.handle), \$binary_mode);")
    @eval begin
        function $f{T}(s::AbstractString, cloud::PointCloud{T};
                binary_mode::Bool=true)
            ret = $ex
            if ret != 0
                error("failed to $f: code $ret")
            end
            ret
        end
    end
end

# Create PointCloud instance and then load PCD data
function (::Type{PointCloud{T}}){T}(path::AbstractString)
    handle = @boostsharedptr "pcl::PointCloud<\$T>"
    cloud = PointCloud(handle)
    @assert !isempty(path)
    load(path, cloud)
    return cloud
end

end # module


text:
using ImageCore, Colors, ColorVectorSpace
using Base.Test

# Different access patterns (getindex)
function mysum_elt_boundscheck(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end
function mysum_index_boundscheck(A)
    s = zero(eltype(A))
    for I in eachindex(A)
        s += A[I]
    end
    s
end
function mysum_elt_inbounds(A)
    s = zero(eltype(A))
    @inbounds for a in A
        s += a
    end
    s
end
function mysum_index_inbounds_simd(A)
    s = zero(eltype(A))
    @inbounds @simd for I in eachindex(A)
        s += A[I]
    end
    s
end
# setindex!
function myfill1!(A, val)
    f = convert(eltype(A), val)
    for I in eachindex(A)
        A[I] = f
    end
    A
end
function myfill2!(A, val)
    f = convert(eltype(A), val)
    @inbounds @simd for I in eachindex(A)
        A[I] = f
    end
    A
end

# Rather than using BenchmarkTools (and thus run one test repeatedly,
# accumulating timings), we run the same test interleaving the two
# array types. This is designed to reduce the risk of spurious
# failure, particularly on shared machines like Travis where they may
# get "distracted" by other tasks
function test_getindex(f, ar, cv, n)
    t_ar = Array{Float64}(n)
    t_cv = Array{Float64}(n)
    for i = 1:n
        t_ar[i] = @elapsed f(ar)
        t_cv[i] = @elapsed f(cv)
    end
    median(t_ar), median(t_cv)
end
function test_setindex(f, ar, cv, n)
    t_ar = Array{Float64}(n)
    t_cv = Array{Float64}(n)
    for i = 1:n
        t_ar[i] = @elapsed f(ar, zero(eltype(ar)))
        t_cv[i] = @elapsed f(cv, zero(eltype(cv)))
    end
    median(t_ar), median(t_cv)
end

ssz = (1000,1000)
a = rand(3,ssz...)
c = reinterpret(RGB{Float64}, a, ssz)
vchan = ChannelView(c)
vcol = ColorView{RGB}(a)
cc_getindex_funcs = (mysum_elt_boundscheck,
                     mysum_index_boundscheck,
                     mysum_elt_inbounds,
                     mysum_index_inbounds_simd)
cc_setindex_funcs = (myfill1!,
                     myfill2!)
chanvtol = Dict(mysum_index_inbounds_simd => 20,   # @simd doesn't work for ChannelView :(
                mysum_elt_boundscheck => 20,
                myfill1! => 20,                    # crappy setindex! performance
                myfill2! => 20)
chanvdefault = 10
colvtol = Dict(mysum_elt_boundscheck=>5, mysum_index_boundscheck=>5)
colvdefault = 3

for (suite, testf) in ((cc_getindex_funcs, test_getindex),
                       (cc_setindex_funcs, test_setindex))
    for f in suite
        # Channelview
        t_ar, t_cv = testf(f, a, vchan, 10^2)
        tol = haskey(chanvtol, f) ? chanvtol[f] : chanvdefault
        if t_cv >= tol*t_ar
            println("ChannelView: failed on $f, time ratio $(t_cv/t_ar), tol $tol")
        end
        @test t_cv < tol*t_ar
        # ColorView
        t_ar, t_cv = testf(f, c, vcol, 10^2)
        tol = haskey(colvtol, f) ? colvtol[f] : colvdefault
        if t_cv >= tol*t_ar
            println("ColorView: failed on $f, time ratio $(t_cv/t_ar), tol $tol")
        end
        @test t_cv < tol*t_ar
    end
end


text:
using Documenter, JLSO

makedocs(
    modules=[JLSO],
    format=Documenter.HTML(
        assets=["assets/invenia.css"],
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages=[
        "Home" => "index.md",
        "Metadata" => "metadata.md",
        "Upgrading" => "upgrading.md",
        "API" => "api.md",
    ],
    repo="https://github.com/invenia/JLSO.jl/blob/{commit}{path}#L{line}",
    sitename="JLSO.jl",
    authors="Invenia Technical Computing Corporation",
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/invenia/JLSO.jl",
    devbranch = "master",
    push_preview = true,
)

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestSymbolicNode

using MXNet
using Test

using ..Main: mlp2, mlpchain, exec

################################################################################
# Test Implementations
################################################################################
function test_basic()
  @info("SymbolicNode::basic")

  model = mlp2()
  @test mx.list_arguments(model) == [:data,:fc1_weight,:fc1_bias,:fc2_weight,:fc2_bias]
  @test mx.list_outputs(model) == [:fc2_output]
  @test mx.list_auxiliary_states(model) == Symbol[]
end

function test_chain()
  @info("SymbolicNode::chain")

  model = mlpchain()
  @test mx.list_arguments(model) == [:data,:fc1_weight,:fc1_bias,:fc2_weight,:fc2_bias]
  @test mx.list_outputs(model) == [:fc2_output]
  @test mx.list_auxiliary_states(model) == Symbol[]

  let layerconfig = [20, 10, 6]
    model = @mx.chain mx.Variable(:data) =>
      mx.MLP(layerconfig, prefix=:magic_) =>
      mx.LinearRegressionOutput(mx.Variable(:label))

    @test mx.list_arguments(model) == [
      :data,
      :magic_fc1_weight, :magic_fc1_bias,
      :magic_fc2_weight, :magic_fc2_bias,
      :magic_fc3_weight, :magic_fc3_bias,
      :label]
  end
end

function test_internal()
  @info("SymbolicNode::internal")

  data  = mx.Variable(:data)
  oldfc = mx.FullyConnected(data, name=:fc1, num_hidden=10)
  net1  = mx.FullyConnected(oldfc, name=:fc2, num_hidden=100)

  @test mx.list_arguments(net1) == [:data,:fc1_weight,:fc1_bias,:fc2_weight,:fc2_bias]

  internal = mx.get_internals(net1)
  fc1      = internal[:fc1_output]
  @test mx.list_arguments(fc1) == mx.list_arguments(oldfc)
end

function test_get_children()
  @info("SymbolicNode::get_children")

  let x = mx.Variable(:x), y = mx.Variable(:y)
    z = x + y
    @test length(mx.list_outputs(z)) == 1
    @test length(mx.list_outputs(mx.get_children(z))) == 2
    @test mx.list_outputs(mx.get_children(z)) == [:x, :y]
  end

  @info("SymbolicNode::get_children::on leaf")
  let x = mx.Variable(:x)
    @test mx.get_children(x) == nothing
  end
end  # test_get_children

function test_compose()
  @info("SymbolicNode::compose")

  data = mx.Variable(:data)
  net1 = mx.FullyConnected(data, name=:fc1, num_hidden=10)
  net1 = mx.FullyConnected(net1, name=:fc2, num_hidden=100)

  net2 = mx.FullyConnected(mx.SymbolicNode, name=:fc3, num_hidden=10)
  net2 = mx.Activation(net2, act_type=:relu)
  net2 = mx.FullyConnected(net2, name=:fc4, num_hidden=20)

  composed  = net2(fc3_data=net1, name=:composed)
  multi_out = mx.Group(composed, net1)
  @test mx.list_outputs(multi_out) == [:composed_output, :fc2_output]
end

function test_infer_shape()
  @info("SymbolicNode::infer_shape::mlp2")

  model = mlp2()
  data_shape = (100, 100)
  arg_shapes, out_shapes, aux_shapes = mx.infer_shape(model, data=data_shape)
  arg_shape_dict = Dict{Symbol,Tuple}(zip(mx.list_arguments(model), arg_shapes))
  @test arg_shape_dict == Dict{Symbol,Tuple}(:fc2_bias => (10,),:fc2_weight => (1000,10),
                                             :fc1_bias => (1000,), :fc1_weight => (100, 1000),
                                             :data => data_shape)
  @test length(out_shapes) == 1
  @test out_shapes[1] == (10, 100)
end

function test_infer_shape_error()
  @info("SymbolicNode::infer_shape::throws")

  model = mlp2()
  weight_shape = (100, 1)
  data_shape   = (100, 100)
  @test_throws mx.MXError mx.infer_shape(model, data=data_shape, fc1_weight=weight_shape)
end

function test_saveload()
  @info("SymbolicNode::saveload::mlp2")

  model = mlp2()
  fname = tempname()
  mx.save(fname, model)
  model2 = mx.load(fname, mx.SymbolicNode)
  @test mx.to_json(model) == mx.to_json(model2)

  rm(fname)
end

function test_attrs()
  @info("SymbolicNode::Attributes")

  data = mx.Variable(:data)

  @test mx.get_name(data) == :data
  result = mx.get_attr(data, :test)
  @test ismissing(result)
  mx.set_attr(data, :test, "1.0")
  result = mx.get_attr(data, :test)
  @test !ismissing(result)
  @test result == "1.0"

  data2 = mx.Variable(:data2, attrs = Dict(:test => "hallo!"))
  @test mx.get_attr(data2, :test) == "hallo!"

  conv = mx.Convolution(data2, kernel = (1,1), num_filter = 1)
  @test ismissing(mx.get_attr(conv, :b))
  @test mx.get_name(conv) isa Symbol

  @test_throws MethodError mx.Variable(:data3, attrs = Dict(:test => "1.0", :test2 => 1.0))
  @test_throws MethodError mx.Convolution(data2, kernel = (1,1), num_filter = 1, attrs = Dict(:test => "1.0", :test2 => 1.0))
end

function test_functions()
  @info("SymbolicNode::Functions")
  data = mx.Variable(:data)
  typeof(mx.sum(data)) == mx.SymbolicNode
end

function test_reshape()
  @info("SymbolicNode::reshape(sym, dim...)")

  A = mx.NDArray(collect(1:24))
  x = mx.Variable(:x)
  y = mx.reshape(x, 2, 3, 4)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (2, 3, 4)
  @test copy(out) == reshape(1:24, 2, 3, 4)

  @info("SymbolicNode::reshape(sym, dim)")

  A = mx.NDArray(collect(1:24))
  x = mx.Variable(:x)
  y = mx.reshape(x, (2, 3, 4))
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (2, 3, 4)
  @test copy(out) == reshape(1:24, 2, 3, 4)

  @info("SymbolicNode::reshape::reverse")

  A = mx.zeros(10, 5, 4)
  x = mx.Variable(:x)
  y = mx.reshape(x, -1, 0, reverse = true)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (50, 4)

  @info("SymbolicNode::reshape::0")

  A = mx.zeros(2, 3, 4)
  x = mx.Variable(:x)
  y = mx.reshape(x, 4, 0, 2)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (4, 3, 2)

  @info("SymbolicNode::reshape::-1")

  A = mx.zeros(2, 3, 4)
  x = mx.Variable(:x)
  y = mx.reshape(x, 6, 1, -1)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (6, 1, 4)

  @info("SymbolicNode::reshape::-2")

  A = mx.zeros(2, 3, 4, 2)
  x = mx.Variable(:x)
  y = mx.reshape(x, 3, 2, -2)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (3, 2, 4, 2)

  @info("SymbolicNode::reshape::-3")

  A = mx.zeros(2, 3, 4, 5)
  x = mx.Variable(:x)
  y = mx.reshape(x, -3, -3)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (6, 20)

  @info("SymbolicNode::reshape::-4")

  A = mx.zeros(2, 3, 4)
  x = mx.Variable(:x)
  y = mx.reshape(x, 0, 0, -4, 2, 2)
  e = mx.bind(y, mx.cpu(), Dict(:x => A))
  mx.forward(e)
  out = e.outputs[1]

  @test size(out) == (2, 3, 2, 2)
end

function test_dot()
  @info("SymbolicNode::dot")
  x = mx.Variable(:x)
  y = mx.Variable(:y)
  z = mx.dot(x, y)
  z_exec = mx.bind(z, context = mx.cpu(),
                   args = Dict(:x => mx.ones((100, 2)), :y => mx.ones((2, 200))))
  mx.forward(z_exec)

  ret = copy(z_exec.outputs[1])
  @test size(ret) == (100, 200)
  @test ret ≈ 2*ones(100, 200)
end

function test_print()
  @info("SymbolicNode::print")
  io = IOBuffer()
  print(io, mx.Variable(:x))
  @test !isempty(String(take!(io)))
end

function test_misc()
  @info("SymbolicNode::Miscellaneous")
  # Test for #189
  a = mx.Variable("a")
  b = mx.Variable("b")
  symb = mx.ElementWiseSum(a, b)
end

function test_add()
  @info("SymbolicNode::elementwise add")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = exec(x .+ 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == A .+ 42
    end

    let y = exec(42 .+ x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == 42 .+ A
    end

    let y = exec(-1 .+ x .+ 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == -1 .+ A .+ 42
    end
  end

  let A = Float32[1 2; 3 4], B = Float32[2 4; 6 8]
    x = mx.Variable(:x)
    y = mx.Variable(:y)

    let z = x .+ y
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == A .+ B
    end

    let z = y .+ x
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == B .+ A
    end
  end
end  # function test_add

function test_minus()
  @info("SymbolicNode::elementwise minus")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = exec(x .- 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == A .- 42
    end

    let y = exec(42 .- x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == 42 .- A
    end

    let y = exec(-1 .- x .- 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == -1 .- A .- 42
    end

    let y = exec(-x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == -A
    end
  end

  let A = Float32[1 2; 3 4], B = Float32[2 4; 6 8]
    x = mx.Variable(:x)
    y = mx.Variable(:y)

    let z = x .- y
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == A .- B
    end

    let z = y .- x
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == B .- A
    end
  end
end  # function test_minus

function test_mul()
  @info("SymbolicNode::elementwise mul")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = exec(x .* 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == A .* 42
    end

    let y = exec(42 .* x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == 42 .* A
    end

    let y = exec(-1 .* x .* 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) == -1 .* A .* 42
    end
  end

  let A = Float32[1 2; 3 4], B = Float32[2 4; 6 8]
    x = mx.Variable(:x)
    y = mx.Variable(:y)

    let z = x .* y
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == A .* B
    end

    let z = y .* x
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) == B .* A
    end
  end
end  # function test_mul

function test_div()
  @info("SymbolicNode::elementwise div")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = exec(x ./ 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) ≈ A ./ 42
    end

    let y = exec(42 ./ x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) ≈ 42 ./ A
    end

    let y = exec(-1 ./ x ./ 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) ≈ -1 ./ A ./ 42
    end
  end

  let A = Float32[1 2; 3 4], B = Float32[2 4; 6 8]
    x = mx.Variable(:x)
    y = mx.Variable(:y)

    let z = x ./ y
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) ≈ A ./ B
    end

    let z = y ./ x
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) ≈ B ./ A
    end
  end
end  # function test_div

function test_power()
  @info("SymbolicNode::elementwise power")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = exec(x .^ 42; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) ≈ A .^ 42
    end

    let y = exec(42 .^ x; :x => A)[]
      @test size(y) == size(A)
      @test copy(y) ≈ 42 .^ A
    end
  end

  let A = Float32[1 2; 3 4], B = Float32[2 4; 6 8]
    x = mx.Variable(:x)
    y = mx.Variable(:y)

    let z = x .^ y
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) ≈ A .^ B
    end

    let z = y .^ x
      z = exec(z; :x => A, :y => B)[]

      @test size(z) == size(A)
      @test copy(z) ≈ B .^ A
    end
  end

  @info("SymbolicNode::power::e .^ x::x .^ e")
  let x = mx.Variable(:x), A = [0 0 0; 0 0 0]
    y = exec(ℯ .^ x; :x => A)[]
    @test copy(y) ≈ fill(1, size(A))
  end

  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = ℯ .^ x
      z = exec(y; :x => A)[]
      @test copy(z) ≈ ℯ .^ A
    end

    let y = x .^ ℯ
      z = exec(y; :x => A)[]
      @test copy(z) ≈ A .^ ℯ
    end
  end

  @info("SymbolicNode::power::π .^ x::x .^ π")
  let x = mx.Variable(:x), A = Float32[1 2; 3 4]
    let y = π .^ x
      z = exec(y; :x => A)[]
      @test copy(z) ≈ π .^ A
    end

    let y = x .^ π
      z = exec(y; :x => A)[]
      @test copy(z) ≈ A .^ π
    end
  end
end  # function test_power

function test_get_name()
  @info("SymbolicNode::get_name::with get_internals")
  name = mx.get_name(mx.get_internals(mlp2()))  # no error
  @test occursin("Ptr", name)
end  # function test_get_name

function test_var()
  @info("SymbolicNode::var")
  x = @mx.var x
  @test x isa mx.SymbolicNode

  x′ = @mx.var x
  @test x.handle != x′.handle

  x, y, z = @mx.var x y z
  @test x isa mx.SymbolicNode
  @test y isa mx.SymbolicNode
  @test z isa mx.SymbolicNode
end  # test_var

################################################################################
# Run tests
################################################################################
@testset "SymbolicNode Test" begin
  test_basic()
  test_chain()
  test_internal()
  test_compose()
  test_infer_shape()
  test_infer_shape_error()
  test_saveload()
  test_attrs()
  test_functions()
  test_reshape()
  test_dot()
  test_print()
  test_misc()
  test_add()
  test_minus()
  test_mul()
  test_div()
  test_power()
  test_get_name()
  test_var()
end

end


text:
"""
    tfactors(x::Incremental)

Returns a vector containing the incremental development factors of
the accumulated triangle built out of the incremental triangle `x`.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> tfactors(ti)
3-element Vector{Float64}:
 2.0
 1.25
 1.0
```
"""
function tfactors(x::Incremental)
    tmp = taccum(x)
    tfactors = ones(Float64, size(x.claims)[2])
    for i in 1:(size(tmp.claims)[2] - 1)
        tfactors[i] = sum(tmp.claims[:,i + 1]) / sum(tmp.claims[1:(end - i), i])
    end
    return tfactors
end

"""
    tfactors(x::Accumulated)

Returns a vector containing the incremental development factors of the accumulated triangle `x`.

# Example
````
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> tfactors(ta)
3-element Vector{Float64}:
 2.0
 1.25
 1.0
```
"""
function tfactors(x::Accumulated)
    tfactors = ones(Float64, size(x.claims)[2])
    for i in 1:(size(x.claims)[2] - 1)
        tfactors[i] = sum(x.claims[:,i + 1]) / sum(x.claims[1:(end - i), i])
    end
    return tfactors
end

"""
    fillcl(x::Accumulated)

Returns a completed accumulated triangle, which contains both the
known accumulated claims of accumulated triangle `x` and the chain ladder
estimated accumulated claims.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> fillcl(ta)
Accumulated Completed Triangle {Float64}:
 3.0  8.0  10.0
 4.0  6.0   7.5
 1.0  2.0   2.5
```
"""
function fillcl(x::Accumulated)
    factores = tfactors(x)
    matriz_completa = convert(Matrix{Float64}, deepcopy(x.claims))
    for i in 1:size(matriz_completa, 2)
        for j in 1:size(matriz_completa, 1)
            if (i + j) > (size(matriz_completa, 1) + 1)
                matriz_completa[i, j] = matriz_completa[i, j - 1] * factores[j - 1]
            end
        end
    end
    return AccumulatedCompleted(matriz_completa)
end

"""
    fillcl(x::Incremental)

Returns a completed incremental triangle, which contains both the
known incremental claims of incremental triangle `x` and the chain ladder
estimated incremental claims.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> fillcl(ti)
Incremental Completed Triangle {Float64}:
 3.0  5.0  2.0
 4.0  2.0  1.5
 1.0  1.0  0.5
```
"""
function fillcl(x::Incremental)
    factores = tfactors(x)
    tri_acum = taccum(x)
    matriz_completa = convert(Matrix{Float64}, deepcopy(tri_acum.claims))
    for i in 1:size(matriz_completa, 2)
        for j in 1:size(matriz_completa, 1)
            if (i + j) > (size(matriz_completa, 1) + 1)
                matriz_completa[i, j] = matriz_completa[i, j - 1] * factores[j - 1]
            end
        end
    end
    des = zeros(Float64, size(matriz_completa, 1), size(matriz_completa, 1))
    des[:,1] = matriz_completa[:,1]
    for j in 2:size(matriz_completa, 1)
        for i in 1:size(matriz_completa, 1)
            des[i, j] = matriz_completa[i, j] - matriz_completa[i, j - 1]
        end
    end
    return IncrementalCompleted(des)
end

"""
    ibnrcl(x::Accumulated)

Returns the estimated incurred but not reported claim reserve for the
accumulated triangle `x` using the chain ladder methodology.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> ibnrcl(ta)
3.0
```
"""
function ibnrcl(x::Accumulated)
    tri_completado = fillcl(x)
    ib = 0
    for i in 1:size(tri_completado.claims, 2)
        j = (tri_completado.claims[i,end] - tri_completado.claims[i, end - i + 1])
        ib = ib + j
    end
    return ib
end

"""
    ibnrcl(x::Incremental)

Returns the estimated incurred but not reported claim reserve for the
incremental triangle `x` using the chain ladder methodology.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> ibnrcl(ti)
3.0
```
"""
function ibnrcl(x::Incremental)
    y = taccum(x)
    return ibnrcl(y)
end


text:
"""
    typeinf_lock(f)

Type inference lock. This prevents you from recursing into type inference when you don't want.
equivalent to the following code, which you may see in Julia compiler implementation.

```julia
ccall(:jl_typeinf_begin, Cvoid, ())
ret = f()
ccall(:jl_typeinf_end, Cvoid, ())
return ret
```
"""
function typeinf_lock(f)
    ccall(:jl_typeinf_begin, Cvoid, ())
    ret = f()
    ccall(:jl_typeinf_end, Cvoid, ())
    return ret
end


text:
module Optim

import Random
import ..Storage

function RandomLocalSearch(problem)

    function optimize(state, stateStorage, data, iterations, rng)
        solution = stateStorage[1]
        a = 0
        for i in 1:iterations
            direction = Tuple(Random.rand(rng, problem.neighborSpace))
            delta = problem.propose(solution, data, direction)
            if delta < 0
                problem.apply(solution, data, direction)
            end
        end
        problem.loss(stateStorage[1], data)
    end

    (
     stateType=nothing, # No need to allocate memory for this optimizer
     init=nothing, # No need to initialize this optimizer
     solutionStates=1, # How many solution states does this solver uses
     optimize=optimize, # What is doing the optimizer to solve the problem
    )

end

struct ExhaustiveLocalSearchState
    currentIteration::Int
    hasImproved::Bool
end

function ExhaustiveLocalSearch(problem)

    function init(state)
        state.currentIteration = 17
        state.hasImproved = false
    end

    function optimize(state, stateStorage, data, iterations, rng)
        solution = stateStorage[1]
        idx = state.currentIteration
        hasImproved = state.hasImproved

        for i in 1:iterations
            direction = problem.neighborSpace[idx]
            delta = problem.propose(solution, data, direction)
            if delta < 0
                state.hasImproved = true
                problem.apply(solution, data, direction)
            end

            idx += 1

            if idx > length(problem.neighborSpace)
                idx = 1
                if !hasImproved
                    # break
                end
            end
        end

        state.currentIteration = idx
        state.hasImproved = hasImproved

        r = problem.loss(stateStorage[1], data)
    end

    (
     stateType=ExhaustiveLocalSearchState, # The type of state required by this optimizer
     init=init, # How to init the state of this optimizer
     solutionStates=1, # How many solution states does this solver uses
     optimize=optimize, # What is doing the optimizer to solve the problem
    )

end

end


text:
function perturb_elements!(enum::EnumParameter, element::Int)
    perturb!(enum.values[element])
end

function perturb_elements!(enum::EnumParameter)
    for parameter in enum.values
        if !(typeof(parameter) <: StringParameter)
            perturb!(parameter)
        end
    end
end

function perturb_elements!(enum::EnumParameter, interval::Array)
    if length(interval) > length(enum.values)
        error("too many intervals.")
    end
    for i = 1:length(interval)
        parameter = enum.values[i]
        if typeof(parameter) <: NumberParameter
            @inbounds interv = interval[i]
            perturb!(parameter, interv)
        end
    end
end

function perturb_elements!(enum::EnumParameter, element::Int, interval::Number)
    perturb!(enum.values[element], interval)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9fb05e18-3480-11ec-372e-f36b56f00c89
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 99f842de-080e-43b8-9d44-63a4dff4a7f7
begin
	using AE
	using Plots
	plotly()
end

# ╔═╡ dc1724ff-096c-4d30-a811-bc0c1d92bf80
using TSne

# ╔═╡ 867597af-8bd0-4515-a114-c07202d4043b
md"
# Autoencoder Demo

JingYu Ning
"

# ╔═╡ 2874e9ab-e919-48b8-84a7-31353cd39be2
md"
## Autoencoder with 32-dim of embedding
"

# ╔═╡ 6f975bbf-be08-44ee-a5d6-4c7c32491945
md"
### Get data from MNIST dataset
"

# ╔═╡ 390a55c6-0846-4588-9686-05ebfb6d5e0b
data, label = AE.get_test_data_and_label();

# ╔═╡ ad014136-6970-46c4-9a78-88e1c10dada4
md"
### Trained model
"

# ╔═╡ 81184495-8fd8-491e-a761-ce371bb05adb
m = AE.get_model("model_h")

# ╔═╡ 243f552f-802f-4489-aa76-dce6cc8bab8c
md"
### Inference
"

# ╔═╡ 509bbbb6-d1e8-4593-b829-c8d670d885a2
function plot_data(i::Integer, size=(580, 280))
	p1 = heatmap(
		data[:, end:-1:1, i]',
		color=:coolwarm, clim=(-1, 1),
		axis=false, ticks=[], colorbar=false, xlabel="Ground Truth",
		size=(280, 280)
	)
	p2 = heatmap(
		m(data[:, :, i:i])[:, end:-1:1, 1]',
		color=:coolwarm, clim=(-1, 1),
		axis=false, ticks=[], colorbar=false, xlabel="Inferenced data",
		size=(280, 280)
	)

	return plot(p1, p2, size=size)
end;

# ╔═╡ 1f31a397-521c-419e-be91-3457e3d5e429
plot_data(4, (290, 140))

# ╔═╡ aa46c772-e624-46cf-af64-1821e2e1acea
plot_data(3, (290, 140))

# ╔═╡ b0ddc77d-13ca-493c-ab37-e9d9e90683ef
plot_data(2, (290, 140))

# ╔═╡ 8fd27e2b-e4c1-4708-bf3c-92560ba1da4f
plot_data(19, (290, 140))

# ╔═╡ 7b818fe7-0bcf-4199-87c8-3f1339764fd6
plot_data(5, (290, 140))

# ╔═╡ c6700ed5-8ca7-46bb-835f-ae1ef9495651
plot_data(9, (290, 140))

# ╔═╡ 38ab2928-1fd8-46a4-a82e-79c29857743e
plot_data(12, (290, 140))

# ╔═╡ 11bafa2b-7396-4b19-bf0e-c58dacada767
plot_data(1, (290, 140))

# ╔═╡ 5f5b90e6-7695-46ce-b598-cbb3e37be1b9
plot_data(62, (290, 140))

# ╔═╡ df959e23-199f-4410-9475-a77ea77c24f6
plot_data(8, (290, 140))

# ╔═╡ 81dcfe8e-1216-480d-8b0d-0550f7d17149
md"
### Dimention reduction for visualization
"

# ╔═╡ d17b1e5e-390d-4b6a-addd-4f2096646113
m_dim_reduced = AE.get_model()

# ╔═╡ 32538eb6-57af-4582-9854-a2089ae934fb
embed = m_dim_reduced[1:2]

# ╔═╡ 15d6fc1a-e489-41a4-adff-a276a39f8415
embeddiing = vcat(embed(data), reshape(label, 1, :))

# ╔═╡ e39a081f-27b5-4934-9259-1d87c9077d35
scatter(
	embeddiing[1, :], embeddiing[2, :], embeddiing[3, :],
	group=Int.(embeddiing[end, :]),
	markersize=2, size=(650, 850), legend=:right
)

# ╔═╡ 1a5e67b2-0074-40db-9c82-23f524a3c7d5
md"
## t-SNE
"

# ╔═╡ 63e1212d-f7a4-422c-bb1c-cf98160dab4a
begin
	rescale(A; dims=1) = (A .- mean(A, dims=dims)) ./ max.(std(A, dims=dims), eps())

	embedding32 = permutedims(m[1:3](data), (2, 1))
	embedding32_n = rescale(embedding32)

	tsne_dim_reduction = tsne(embedding32_n, 2)
end;

# ╔═╡ 722aa1f3-10b7-434b-9947-78b5d4fbff42
scatter(
	tsne_dim_reduction[:, 1], tsne_dim_reduction[:, 2],
	group=label, markersize=3, size=(650, 600),
	legend=:outertopright,
)

# ╔═╡ Cell order:
# ╟─867597af-8bd0-4515-a114-c07202d4043b
# ╟─9fb05e18-3480-11ec-372e-f36b56f00c89
# ╠═99f842de-080e-43b8-9d44-63a4dff4a7f7
# ╟─2874e9ab-e919-48b8-84a7-31353cd39be2
# ╟─6f975bbf-be08-44ee-a5d6-4c7c32491945
# ╠═390a55c6-0846-4588-9686-05ebfb6d5e0b
# ╟─ad014136-6970-46c4-9a78-88e1c10dada4
# ╠═81184495-8fd8-491e-a761-ce371bb05adb
# ╟─243f552f-802f-4489-aa76-dce6cc8bab8c
# ╟─509bbbb6-d1e8-4593-b829-c8d670d885a2
# ╟─1f31a397-521c-419e-be91-3457e3d5e429
# ╟─aa46c772-e624-46cf-af64-1821e2e1acea
# ╟─b0ddc77d-13ca-493c-ab37-e9d9e90683ef
# ╟─8fd27e2b-e4c1-4708-bf3c-92560ba1da4f
# ╟─7b818fe7-0bcf-4199-87c8-3f1339764fd6
# ╟─c6700ed5-8ca7-46bb-835f-ae1ef9495651
# ╟─38ab2928-1fd8-46a4-a82e-79c29857743e
# ╟─11bafa2b-7396-4b19-bf0e-c58dacada767
# ╟─5f5b90e6-7695-46ce-b598-cbb3e37be1b9
# ╟─df959e23-199f-4410-9475-a77ea77c24f6
# ╟─81dcfe8e-1216-480d-8b0d-0550f7d17149
# ╠═d17b1e5e-390d-4b6a-addd-4f2096646113
# ╠═32538eb6-57af-4582-9854-a2089ae934fb
# ╠═15d6fc1a-e489-41a4-adff-a276a39f8415
# ╟─e39a081f-27b5-4934-9259-1d87c9077d35
# ╟─1a5e67b2-0074-40db-9c82-23f524a3c7d5
# ╠═dc1724ff-096c-4d30-a811-bc0c1d92bf80
# ╠═63e1212d-f7a4-422c-bb1c-cf98160dab4a
# ╟─722aa1f3-10b7-434b-9947-78b5d4fbff42


text:
abstract type Instruction end

struct Mask <: Instruction
    mask::AbstractString
end

struct Write <: Instruction
    address::Int64
    value::Int64
end

function read_data(lines)
    program = Instruction[]
    for line in lines
        first, v = split(line, " = ")
        if first == "mask"
            push!(program, Mask(v))
        else
            address = parse(Int64, first[5:end-1])
            value = parse(Int64, v)
            push!(program, Write(address, value))
        end
    end
    program
end

function initialization(program)
    memory = Dict{Int64, Int64}()
    mask = nothing
    for instruction in program
        if isa(instruction, Mask)
            mask = instruction.mask
        else
            s = bitstring(instruction.value)
            masked = join(map(((bit, m),) -> m == 'X' ? bit : m, zip(s[length(s)-35:end], mask)), "")
            masked_value = parse(Int64, masked; base=2)
            memory[instruction.address] = masked_value
        end
    end
    memory
end

function main()
    lines = readlines("../resources/input_14.txt")
    program = read_data(lines)
    memory = initialization(program)
    println(sum(values(memory)))
end

main()


text:
module SimpleExponentialSmoothingTests

using Test

using TSeriesForecast: SES, loss, forecast, fit

ϵ = 1.0e-6

### Data

α = 0.83389014338874
l0 = 446.586793455216
years = collect(1996:2013)
observations = [445.364098092, 453.195010427, 454.409641012, 422.378905779, 456.037121728, 440.386604674, 425.194372519, 486.20517351, 500.429086073, 521.27590917, 508.947617045, 488.888857729, 509.870574978, 456.722912348, 473.816602915, 525.950870565, 549.833807597, 542.340469826]
y_pred = [446.586793455216, 445.567199843463, 451.927955904711, 453.997408654674, 427.631050757941, 451.318593352274, 442.202515745824, 428.019592751632, 476.539975033385, 496.460869263639, 517.153886449352, 510.31075927907, 492.447246723823, 506.976388420008, 465.070510052539, 472.363790683708, 517.049528409702, 544.388015682095]

y_forecast = [542.6805873746]

model = SES(α, l0)

### Tests

@testset "Simple Exponential Smoothing" begin
    @testset "SES struct" begin

        @testset "SES default constructor" begin
            model = SES()
            @test model.α ≈ 0
            @test model.l0 ≈ 0
        end

        @testset "SES parametric constructor" begin
            α = 3.14
            l0 = 42.0
            model = SES(α, l0)
            @test model.α ≈ α
            @test model.l0 ≈ l0
        end

        @testset "SES parametric constructor takes any Number type" begin
            α = 42
            l0 = Float32(3.1415)
            model = SES(α, l0)
            @test model.α ≈ α
            @test model.l0 ≈ l0
        end
    end

    @testset "loss" begin
        y = observations
        expected_loss = sum((y - y_pred) .^ 2)

        @test isapprox(loss(model, y), expected_loss)
    end

    @testset "fit" begin
        y = observations

        starting_point = SES()
        optimal_model = fit(starting_point, y)

        @test isapprox(optimal_model.α, α, atol=0.0001)
        @test isapprox(optimal_model.l0, l0, atol=0.1)
    end

    @testset "forecast" begin
        @testset "The forecast is correct" begin
            @test forecast(model, observations, 1) ≈ y_forecast
        end

        @testset "forecast returns an array of length `h`" begin
            @test length(forecast(model, observations, 1)) == 1
        end

        @testset "The forecast vector repeats the same value" begin
            @test forecast(model, observations, 42) ≈ repeat(y_forecast, 42)
        end
    end

end

end


text:
using Dates

"""
Runs a Monte Carlo simulation with the FUND model over it's distributional parameters.
`trials`: the number of trials to run.
`ntimesteps`: how many timesteps to run.
`output_dir`: an output directory; if none provided, will create and use "output/yyyy-mm-dd HH-MM-SS MCtrials".
`save_trials`: whether or not to generate and save the MC trial values up front to a file.
"""
function run_fund_mcs(trials = 10000; ntimesteps = MimiFUND.default_nsteps + 1, output_dir = nothing, save_trials = false)

    # Set up output directories
    output_dir = output_dir === nothing ? joinpath(@__DIR__, "../../output/", "SCC $(Dates.format(now(), "yyyy-mm-dd HH-MM-SS")) MC$trials") : output_dir
    mkpath("$output_dir/results")

    trials_output_filename = save_trials ?  joinpath("$output_dir/trials.csv") :  nothing

    # Get an instance of FUND's mcs
    mcs = getmcs()

    # run monte carlo trials
    res = run(mcs, get_model(), trials; trials_output_filename = trials_output_filename, ntimesteps = ntimesteps, results_output_dir = "$output_dir/results")

    return res
end


text:
using DifferentialEquations
using StochasticDiffEq
using DiffEqCallbacks
using Random
using SparseArrays
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_mbp!(du,u,p,t)
    (Ctilde,Rtilde) = u
    (β,c,γ,S₀,I₀,N) = p
    C = exp(Ctilde)-1.0
    R = exp(Rtilde)-1.0
    S = S₀-C
    I = I₀+C-R
    @inbounds begin
        du[1] = (exp(-Ctilde)-0.5*exp(-2.0*Ctilde))*(β*c*I/N*S)
        du[2] = (exp(-Rtilde)-0.5*exp(-2.0*Rtilde))*(γ*I)
    end
    nothing
end;

# Define a sparse matrix by making a dense matrix and setting some values as not zero
A = zeros(2,2)
A[1,1] = 1
A[2,2] = 1
A = SparseArrays.sparse(A);

# Make `g` write the sparse matrix values
function sir_noise!(du,u,p,t)
    (Ctilde,Rtilde) = u
    (β,c,γ,S₀,I₀,N) = p
    C = exp(Ctilde)-1.0
    R = exp(Rtilde)-1.0
    S = S₀-C
    I = I₀+C-R
    du[1,1] = exp(-Ctilde)*sqrt(β*c*I/N*S)
    du[2,2] = exp(-Rtilde)*sqrt(γ*I)
end;

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [0.0,0.0]; # C,R

p = [0.05,10.0,0.25,990.0,10.0,1000.0]; # β,c,γ,S₀,I₀,N

Random.seed!(1234);

function condition(u,t,integrator,p) # Event when event_f(u,t) == 0
    (Ctilde,Rtilde) = u
    (β,c,γ,S₀,I₀,N) = p
    C = exp(Ctilde)-1.0
    R = exp(Rtilde)-1.0
    S = S₀-C
    I = I₀+C-R
    I
end

function affect!(integrator)
    terminate!(integrator)
end

cb = ContinuousCallback(
        (u,t,integrator)->condition(u,t,integrator,p),
        affect!)

prob_mbp = SDEProblem(sir_mbp!,sir_noise!,u0,tspan,p,noise_rate_prototype=A)

sol_mbp = solve(prob_mbp,
            SRA1(),
            callback=cb,
            saveat=δt);

df_mbp = DataFrame(sol_mbp(sol_mbp.t)')
df_mbp[!,:C] = exp.(df_mbp[!,:x1]) .- 1.0
df_mbp[!,:R] = exp.(df_mbp[!,:x2]) .- 1.0
df_mbp[!,:S] = p[4] .- df_mbp[!,:C]
df_mbp[!,:I] = p[5] .+ df_mbp[!,:C] .- df_mbp[!,:R]
df_mbp[!,:t] = sol_mbp.t

@df df_mbp plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_mbp,SRA1(),callback=cb)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: snowdj/sir-julia

max_issues_repo_name: snowdj/sir-julia

max_forks_repo_name: snowdj/sir-julia


text:
const DevicesModelContainer = Dict{Symbol, DeviceModel}
const BranchModelContainer = Dict{Symbol, DeviceModelForBranches}
const ServicesModelContainer = Dict{Tuple{String, Symbol}, ServiceModel}

"""
    OperationsProblemTemplate(::Type{T}) where {T<:PM.AbstractPowerFormulation}
Creates a model reference of the PowerSimulations Optimization Problem.
# Arguments
- `model::Type{T<:PM.AbstractPowerFormulation}`:
# Example
```julia
template = OperationsProblemTemplate(CopperPlatePowerModel)
```
"""
mutable struct OperationsProblemTemplate
    transmission::Type{<:PM.AbstractPowerModel}
    devices::DevicesModelContainer
    branches::BranchModelContainer
    services::ServicesModelContainer
    function OperationsProblemTemplate(::Type{T}) where {T <: PM.AbstractPowerModel}
        new(T, DevicesModelContainer(), BranchModelContainer(), ServicesModelContainer())
    end
end

OperationsProblemTemplate() = OperationsProblemTemplate(CopperPlatePowerModel)
# TODO: make getter functions here
# Note: use the file test_operations_template to test the getter functions
get_transmission_model(template::OperationsProblemTemplate) = template.transmission

# Note to devs. PSY exports set_model! these names are choosen to avoid name clashes

"""Sets the transmission model in a template"""
function set_transmission_model!(
    template::OperationsProblemTemplate,
    model::Type{<:PM.AbstractPowerModel},
)
    template.transmission = model
    return
end

"""
    Sets the device model in a template using the component type and formulation.
    Builds a default DeviceModel
"""
function set_device_model!(
    template::OperationsProblemTemplate,
    component_type::Type{<:PSY.Device},
    formulation::Type{<:AbstractDeviceFormulation},
)
    set_device_model!(template, DeviceModel(component_type, formulation))
    return
end

"""
    Sets the device model in a template using a DeviceModel instance
"""
function set_device_model!(
    template::OperationsProblemTemplate,
    model::DeviceModel{<:PSY.Device, <:AbstractDeviceFormulation},
)
    _set_model!(template.devices, model)
    return
end

function set_device_model!(
    template::OperationsProblemTemplate,
    model::DeviceModel{<:PSY.Branch, <:AbstractDeviceFormulation},
)
    _set_model!(template.branches, model)
    return
end

"""
    Sets the service model in a template using a name and the service type and formulation. Builds a default ServiceModel with use_service_name set to true.
"""
function set_service_model!(
    template::OperationsProblemTemplate,
    service_name::String,
    service_type::Type{<:PSY.Service},
    formulation::Type{<:AbstractServiceFormulation},
)
    set_service_model!(
        template,
        service_name,
        ServiceModel(service_type, formulation, use_service_name = true),
    )
    return
end

"""
    Sets the service model in a template using a ServiceModel instance.
"""
function set_service_model!(
    template::OperationsProblemTemplate,
    service_type::Type{<:PSY.Service},
    formulation::Type{<:AbstractServiceFormulation},
)
    set_service_model!(template, ServiceModel(service_type, formulation))
    return
end

function set_service_model!(
    template::OperationsProblemTemplate,
    service_name::String,
    model::ServiceModel{<:PSY.Service, <:AbstractServiceFormulation},
)
    _set_model!(template.services, service_name, model)
    return
end

function set_service_model!(
    template::OperationsProblemTemplate,
    model::ServiceModel{<:PSY.Service, <:AbstractServiceFormulation},
)
    _set_model!(template.services, model)
    return
end


text:
"""
Tableau of one-stage, 1st order implicit (backward) Euler method

```julia
TableauImplicitEuler(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Augustin-Louis Cauchy.
    Équations différentielles ordinaires. Cours inédit (fragment). Douzième leçon.
    Ed. Christian Gilain, Etudes Vivantes, 1981.
    Page 102, Equation (5), Θ=1.

"""
function TableauImplicitEuler(::Type{T}=Float64) where {T}
    a = ones(BigFloat, 1, 1)
    b = ones(BigFloat, 1)
    c = ones(BigFloat, 1)
    o = 1

    Tableau{T}(:implicit_euler, o, a, b, c)
end

"Alias for [`TableauImplicitEuler`](@ref)"
TableauBackwardEuler = TableauImplicitEuler

"""
Tableau of two-stage, 2nd order implicit midpoint method

```julia
TableauImplicitMidpoint(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Augustin-Louis Cauchy.
    Équations différentielles ordinaires. Cours inédit (fragment). Douzième leçon.
    Ed. Christian Gilain, Etudes Vivantes, 1981.
    Page 102, Equation (5), Θ=1/2.

"""
function TableauImplicitMidpoint(::Type{T}=Float64) where {T}
    a = ones(BigFloat, 1, 1) ./ 2
    b = ones(BigFloat, 1)
    c = ones(BigFloat, 1) ./ 2
    o = 2

    Tableau{T}(:implicit_midpoint, o, a, b, c)
end

"""
Tableau of symmetric and symplectic three-stage, 4th order Runge-Kutta method

```julia
TableauSRK3(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Shan Zhao and Guo-Wei Wei.
    A unified discontinuous Galerkin framework for time integration.
    Mathematical Methods in the Applied Sciences, Volume 37, Issue 7, Pages 1042-1071, 2014.
    doi: 10.1002/mma.2863.

"""
function TableauSRK3(::Type{T}=Float64) where {T}
    a = @big [[ 5/36         2/9        5/36-√15/10 ]
              [ 5/36         2/9        5/36        ]
              [ 5/36+√15/10  2/9        5/36        ]]
    b = @big  [ 5/18,        4/9,       5/18        ]
    c = @big  [ 1/2-√15/10,  1/2,       1/2+√15/10  ]
    o = 4

    Tableau{T}(:SRK3, o, a, b, c)
end


text:
using Pythagorus
using Test

include("provided_example.jl")

@test length(Object.example_found)==0
@test length(Object.objects_to_try)==10
@test length(Object.unused_predicates)==3

explore_concepts(Singleton)

@test length(Object.example_found)==10
@test length(Object.objects_to_try)==0

## to lazy to write a proper test.


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(sigmoid(-0.8925450194703983 + 0.6110890427800215*sigmoid(-0.16626426792027038 + 0.44845618413236066*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.1740395626916258*sigmoid(0.8063063231762757 + 0.4905142959916642*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8017072131372296*sigmoid(0.9210510999614945 + -0.4935968211971775*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7966067754659663*sigmoid(-0.7328670881809516 + 0.6102429954796835*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + sigmoid(-0.7008872510227278 + -0.258173144464545*sigmoid(-0.16626426792027038 + 0.44845618413236066*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.49550787111793193*sigmoid(0.8063063231762757 + 0.4905142959916642*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8462943069272191*sigmoid(0.9210510999614945 + -0.4935968211971775*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.05869131996649157*sigmoid(-0.7328670881809516 + 0.6102429954796835*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + sigmoid(-0.2643911691309153 + -0.31660511055293883*sigmoid(-0.16626426792027038 + 0.44845618413236066*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.09136517532701127*sigmoid(0.8063063231762757 + 0.4905142959916642*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.36400896662330684*sigmoid(0.9210510999614945 + -0.4935968211971775*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.4876638787119827*sigmoid(-0.7328670881809516 + 0.6102429954796835*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + sigmoid(0.6679332747384503 + -0.6061375406574494*sigmoid(-0.16626426792027038 + 0.44845618413236066*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7532640510772808*sigmoid(0.8063063231762757 + 0.4905142959916642*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.9959397069462161*sigmoid(0.9210510999614945 + -0.4935968211971775*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8544693505755525*sigmoid(-0.7328670881809516 + 0.6102429954796835*sigmoid(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*sigmoid(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*sigmoid(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*sigmoid(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(element_type = Quad(),
             approximation_type = periodic_derivative_operator(
               derivative_order=1, accuracy_order=4, xmin=0.0, xmax=1.0, N=50),
             surface_flux = flux_hll,
             volume_integral = VolumeIntegralWeakForm())

equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

mesh = DGMultiMesh(dg, coordinates_min=(-1.0, -1.0),
                       coordinates_max=( 1.0,  1.0))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                    source_terms = source_terms)

tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
stepsize_callback = StepsizeCallback(cfl=1.0)
callbacks = CallbackSet(summary_callback, alive_callback, stepsize_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# NOTE:  backend should implement `html_body` and `html_head`

# CREDIT: parts of this implementation were inspired by @joshday's PlotlyLocal.jl
@require Revise begin
    Revise.track(Plots, joinpath(Pkg.dir("Plots"), "src", "backends", "web.jl"))
end

function standalone_html(plt::AbstractPlot; title::AbstractString = get(plt.attr, :window_title, "Plots.jl"))
    """
    <!DOCTYPE html>
    <html>
        <head>
            <title>$title</title>
            <meta http-equiv="content-type" content="text/html; charset=UTF-8">
            $(html_head(plt))
        </head>
        <body>
            $(html_body(plt))
        </body>
    </html>
    """
end

function open_browser_window(filename::AbstractString)
    @static if is_apple()
        return run(`open $(filename)`)
    end
    @static if is_linux() || is_bsd()    # is_bsd() addition is as yet untested, but based on suggestion in https://github.com/JuliaPlots/Plots.jl/issues/681
        return run(`xdg-open $(filename)`)
    end
    @static if is_windows()
        return run(`$(ENV["COMSPEC"]) /c start "" "$(filename)"`)
    end
    warn("Unknown OS... cannot open browser window.")
end

function write_temp_html(plt::AbstractPlot)
    html = standalone_html(plt; title = plt.attr[:window_title])
    filename = string(tempname(), ".html")
    output = open(filename, "w")
    write(output, html)
    close(output)
    filename
end

function standalone_html_window(plt::AbstractPlot)
    filename = write_temp_html(plt)
    open_browser_window(filename)
end

# uses wkhtmltopdf/wkhtmltoimage: http://wkhtmltopdf.org/downloads.html
function html_to_png(html_fn, png_fn, w, h)
    run(`wkhtmltoimage -f png -q --width $w --height $h --disable-smart-width $html_fn $png_fn`)
end

function show_png_from_html(io::IO, plt::AbstractPlot)
    # write html to a temporary file
    html_fn = write_temp_html(plt)

    # convert that html file to a temporary png file using wkhtmltoimage
    png_fn = tempname() * ".png"
    w, h = plt.attr[:size]
    html_to_png(html_fn, png_fn, w, h)

    # now read that file data into io
    pngdata = readall(png_fn)
    write(io, pngdata)
end


text:
# test partial pose3 constraints and evaluation

using Statistics
using RoME
using Test

##

N=50
fg = initfg()

v1 = addVariable!(fg,:x1, Pose3, N=N) # 0.001*randn(6,N)
f0 = addFactor!(fg, [:x1;], PriorPose3(MvNormal(zeros(6),1e-2*Matrix{Float64}(LinearAlgebra.I, 6,6))))

sigx = 0.01
sigy = 0.01
sigth = 0.0281
mu1 = [0.0;0.0; -10.0]
prpz = PartialPriorRollPitchZ(
  MvNormal( mu1[1:2], [sigx 0.0; 0.0 sigy]^2 ),
  Normal( mu1[3], sigth )
)

mu2 = [20.0,5.0,pi/2]
xyy = PartialPose3XYYaw(
  MvNormal(
    mu2[1:2],
    [sigx 0.0; 0.0 sigy]^2
  ),
  Normal( mu2[3], sigth )
)

v2 = addVariable!(fg,:x2, Pose3, N=N) # randn(6,N)

f1 = addFactor!(fg, [:x2], prpz, graphinit=false)
f2 = addFactor!(fg, [:x1;:x2], xyy, graphinit=false)

##

@testset "test PartialPriorRollPitchZ evaluations" begin

##

# ensure that at least the first pose is already initialized
doautoinit!(fg, :x1)
@test isInitialized(fg, :x1)

X1pts = getVal(fg, :x1)
@test sum(isnan.(X1pts)) == 0

ppts = approxConv(fg, :x1x2f1,:x2)
@test sum(isnan.(ppts)) == 0

ensureAllInitialized!(fg)
@test isInitialized(fg, :x2)

# get values and ensure that a re-evaluation produces consistent results
X2pts = getVal(fg, :x2)
@test sum(isnan.(X2pts)) == 0

# check that only partial states are updated
pts = IIF.approxConv(fg, :x2f1, :x2, N=N)

newdims = collect(DFG.getSolverData(f1).fnc.usrfnc!.partial)

olddims = setdiff(collect(1:6), newdims)

@test size(pts, 1) == 6
@test size(pts, 2) == N

# check that untouched dimensions (not in the partial list) truely remain untouched
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10

# check that the prior new dims are updated to new and correct values
# @show Statistics.mean(pts,dims=2)[newdims]
@test sum(abs.(Statistics.mean(pts,dims=2)[newdims]-mu1[[3;1;2]]) .< [0.5; 0.1; 0.1]) == 3

# ensure a forced re-evaluatoin
@test norm(X2pts[newdims,:] - pts[newdims,:]) < 1.0

# memcheck that the exact same values are used
@test norm(X2pts - getVal(v2)) < 1e-10

##

end

@testset "test residual function of PartialPose3XYYaw" begin

##

tfg = initfg()
X0 = addVariable!(tfg, :x0, Pose3)
X1 = addVariable!(tfg, :x1, Pose3)

xi = zeros(6)
xja = zeros(6)

res = testFactorResidualBinary(xyy, Pose3, Pose3, xi, xja)

@test abs(res[1]-mu2[1]) < 0.3
@test abs(res[2]-mu2[2]) < 0.3
@test abs(res[3]-mu2[3]) < 0.2

##

xjb = zeros(6,1)
xjb[collect(xyy.partial),1] = mu2
# res = zeros(3)
# xyy(res, fmd, idx, meas, xi, xjb)

res = testFactorResidualBinary(xyy, Pose3, Pose3, xi, xjb)

@test 0.0 < norm(res) < 0.3

##

addFactor!(tfg, [:x0;:x1], xyy, graphinit=false)

# meas = getSample(xyy,100)
ccw = IIF._getCCW(tfg, :x0x1f1)
meas = sampleFactor(ccw, 100)

@test norm(Statistics.std(meas[1],dims=2) - [0.01;0.01;0.002]) < 0.05

##

end

@testset "test PartialPose3XYYaw evaluations" begin

##

# get existing and predict new
X2pts = getBelief(fg, :x2) |> getPoints
pts = approxConv(fg, :x1x2f1, :x2, N=N)

# find which dimensions are and and are not updated by XYYaw partial
newdims = collect(DFG.getSolverData(f2).fnc.usrfnc!.partial)
olddims = setdiff(collect(1:6), newdims)

# check the number of points are correct
@test size(pts, 1) == 6
@test size(pts, 2) == N

# ensure the unchanged dimensions actually remain unchanged
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10

for i in 1:N
    pts[6,i] = wrapRad(pts[6,i])
end

# ensure the newly updated values match what is specified in mu2
@show mu2 # mu2 is used for XYYaw
# trying to compare world and body frame values -- not right!
@show Statistics.mean(pts[newdims,:],dims=2)
@test sum(abs.(Statistics.mean(pts[newdims,:],dims=2)-mu2) .< [1.5;1.5;0.3]) == 3

# ensure a re-evaluation of the partial factor updates the partial variable dimensions correclty
@test norm(X2pts[newdims,:] - pts[newdims,:]) < 1.0

# ensure that memory pointers are working correctly
memcheck = getVal(v2)
@test norm(X2pts - memcheck) < 1e-10

##

end

@testset "test predictbelief with two functions" begin

##

val, = predictbelief(fg, :x2, ls(fg, :x2), N=N)

for i in 1:N
  val[6,i] = wrapRad(val[6,i])
end

@test size(val, 1) == 6
@test size(val, 2) == N

estmu1mean = Statistics.mean(val[collect(DFG.getSolverData(f1).fnc.usrfnc!.partial),:],dims=2)
estmu2mean = Statistics.mean(val[collect(DFG.getSolverData(f2).fnc.usrfnc!.partial),:],dims=2)

@show estmu1mean
@show estmu2mean
@test sum(abs.(estmu1mean - mu1[[3;1;2]]) .< [0.7; 0.1; 0.1]) == 3
@test sum(abs.(estmu2mean - mu2) .< [1.0; 1.5; 0.2] ) == 3

memcheck = getVal(v2)
@test 1e-10 < norm(val - memcheck)

##

end

#


text:
function _bioclim_score(x::Number)
    if isnan(x)
        return NaN
    end
    @assert 0. ≤ x ≤ 1.
    if x ≥ 0.5
        return 2(1-x)
    else
        return 2x
    end
end

function bioclim(layer::T, records::GBIFRecords) where {T <: SimpleSDMLayer}
    bioclim_prediction = similar(layer)
    qf = StatsBase.ecdf(filter(!isnan, layer[records]))
    bioclim_prediction.grid = _bioclim_score.(qf.(layer.grid))
    return bioclim_prediction
end

function Base.min(l1::T, l2::T) where {T <: SimpleSDMLayer}
    SimpleSDMLayers._layers_are_compatible(l1, l2)
    min_layer = similar(l1)
    for i in eachindex(l1.grid)
        if !isnan(l1[i])
            min_layer[i] = min(l1[i], l2[i])
        end
    end
    return min_layer
end


text:
#directory info

SRCDIR = dirname(pathof(ThreeBodyTB))
EXAMPLESDIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "examples")
TESTDIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "test")
TEMPLATEDIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "template_inputs")
STRUCTDIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "reference_structures")

PSEUDODIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "pseudo", "gbrv_pbesol")

DATSDIR1 = joinpath(dirname(pathof(ThreeBodyTB)), "..", "dats", "pbesol", "v1.2")

#DATSDIR1 = joinpath(dirname(pathof(ThreeBodyTB)), "..", "dats", "pbesol", "v0.9")

#DATSDIR1 = joinpath(dirname(pathof(ThreeBodyTB)), "..", "dats", "pbesol", "v0.4")

DATSDIR2 = joinpath(dirname(pathof(ThreeBodyTB)), "..", "dats", "pbesol", "v1.0")

DOCSDIR = joinpath(dirname(pathof(ThreeBodyTB)), "..", "docs")

#DATSDIR1 = "/home/kfg/codes/TB_fit/fit_dimer/datab_2/"
#DATSDIR2 = ""

#DATSDIR2 = "/home/kfg/codes/TB_fit/binary_v10/datab/"

global MPI_STRING="mpirun -np"
global QE_BIN_DIR_STRING="/home/kfg/codes/q-e-qe-6.5/bin/"
global PSEUDOS=PSEUDODIR
global TEMPLATES=TEMPLATEDIR

global WANNIER_BIN_DIR_STRING="/users/kfg/codes/wannier90-2.1.0/" #not necessary


text:
export TFIM, hamiltonian
abstract type AbstractModel{D} end

nspin(model::AbstractModel) = prod(model.size)

struct TFIM{D} <:AbstractModel{D}
    size::NTuple{D, Int}
    Γ::Real
    periodic::Bool
    TFIM(size::Int...; Γ::Real, periodic::Bool) = new{length(size)}(size, Γ, periodic)
end

function get_bonds(model::AbstractModel{1})
    nbit, = model.size
    [(i, i%nbit+1) for i in 1:(model.periodic ? nbit : nbit-1)]
end

function get_bonds(model::AbstractModel{2})
    m, n = model.size
    cis = LinearIndices(model.size)
    bonds = Tuple{Int, Int}[]
    for i=1:m, j=1:n
        (i!=m || model.periodic) && push!(bonds, (cis[i,j], cis[i%m+1,j]))
        (j!=n || model.periodic) && push!(bonds, (cis[i,j], cis[i,j%n+1]))
    end
    bonds
end

function hamiltonian(model::TFIM)
    nbit = nspin(model)
    -sum([repeat(nbit, Z, (i,j)) for (i,j) in get_bonds(model)]) +
    -sum([put(nbit, i=>X) for i=1:nbit])*model.Γ
end

#https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.70.3339
struct SY <:AbstractModel{0}
    size::Int
    J::Real
end

function hamiltonian(model::SY)
    nbit = nspin(model)
    h = Add{nbit}()
    coupling = randn(nbit*(nbit-1)÷2)*model.J/sqrt(2*nbit)
    t = 1
    for i in 1:nbit
        for j in i+1:nbit
            for σ in (X, Y, Z)
                push!(h, repeat(nbit, σ, (i,j)) * coupling[t] )
            end
            t += 1
        end
    end
    h
end


text:
export ray_class_group

mutable struct MapRayClassGrp{T} <: Map{T, Hecke.NfOrdIdlSet}
  header::Hecke.MapHeader
  modulus_fin::NfOrdIdl
  modulus_inf::Array{InfPlc,1}

  function MapRayClassGrp{T}() where {T}
    return new{T}()
  end
end

#
# Modify the map of the class group so that the chosen representatives are coprime to m
#

function _coprime_ideal(C::GrpAbFinGen, mC::Map, m::NfOrdIdl)

  O=parent(m).order
  K=nf(O)
  L=NfOrdIdl[]
  for i=1:ngens(C)
    a=mC(C[i])
    if iscoprime(a,m)
      push!(L,a)
    else
      J=inv(a)
      s=K(rand(J.num,5))//J.den  # Is the bound acceptable?
      I=s*a
      simplify(I)
      I = numerator(I)
      while !iscoprime(I,m)
        s=K(rand(J.num,5))//J.den
        I=s*a
        simplify(I)
        I = numerator(I)
      end
      push!(L,I)
    end
  end

  function exp(a::GrpAbFinGenElem)
    I=ideal(O,1)
    for i=1:ngens(C)
      if Int(a.coeff[1,i])!= 0
        I=I*(L[i]^(Int(a.coeff[1,i])))
      end
    end
    return I
  end

  mp=Hecke.MapClassGrp{typeof(C)}()
  mp.header=Hecke.MapHeader(C,Hecke.NfOrdIdlSet(O),exp, mC.header.preimage)

  return mp

end

doc"""
***
    ray_class_group(m::NfOrdIdl, A::Array{InfPlc,1}=[]) -> FinGenGrpAb, Map

> Compute the ray class group of the maximal order $L$ with respect to the modulus given by $m$ (the finite part) and the infinite primes of $A$
> and return an abstract group isomorphic to the ray class group with a map
> from the group to the ideals of $L$

"""
function ray_class_group(m::NfOrdIdl, primes::Array{InfPlc,1}=InfPlc[])

#
# We compute the group using the sequence U -> (O/m)^* _> Cl^m -> Cl -> 1
# First of all, we compute all these groups with their own maps
#
  O=parent(m).order
  K=nf(O)

  C, mC = class_group(O)
  mC=_coprime_ideal(C,mC,m)
  U, mU = unit_group(O)
  M, pi= quo(O,m)
  G, mG=unit_group(M)

  p = [ x for x in primes if isreal(x) ]
  if !isempty(p)
    H,lH,eH=_infinite_primes(O,p,m)
    T=G
    G=direct_product(G,H)
  end

#
# We start to construct the relation matrix
#
  RG=rels(G)
  RC=rels(C)

  A=vcat(RC, MatrixSpace(FlintZZ, ngens(G)+ngens(U), cols(RC))())
  B=vcat(MatrixSpace(FlintZZ, ngens(C), cols(RG))(), RG)
  B=vcat(B, MatrixSpace(FlintZZ, ngens(U) , cols(RG))())

#
# We compute the relation matrix given by the image of the map U -> (O/m)^*
#
  for i=1:ngens(U)
    u=mU(U[i])
    a=(mG\(pi(u))).coeff
    if !isempty(primes)
      a=hcat(a, (lH(K(u))).coeff)
    end
    for j=1:ngens(G)
      B[i+rows(RC)+rows(RG),j]=a[1,j]
    end
  end

#
# We compute the relation between generators of Cl and (O/m)^* in Cl^m
#

  for i=1: ngens(C)
    if order(C[i])!=1
      y=Hecke.principal_gen((mC(C[i]))^(Int(order(C[i]))))
      b=(mG\(pi(y))).coeff
      if primes != []
        b=hcat(b, (lH(K(y))).coeff)
      end
      for j=1: ngens(G)
        B[i,j]=-b[1,j]
      end
    end
  end

  R=hcat(A,B)
  X=AbelianGroup(R)

#
# Discrete logarithm
#

  function disclog(J::NfOrdIdl)

    !iscoprime(J,m) && error("The ideal is not coprime to the modulus")
    if isone(J)
      return X([0 for i=1:ngens(X)])
    else
      L=mC\J
      s=mC(L)
      I= J // s
      simplify(I)
      gamma=Hecke.principal_gen(I.num)
      y=((mG\(pi(gamma)))-(mG\(pi(O(I.den))))).coeff
      if primes!=[]
        z=lH(K(gamma))
        y=hcat(y, z.coeff)
      end
      return X(hcat(L.coeff,y))
    end
  end

#
# Exp map
#

  function expo(a::GrpAbFinGenElem)
    b=C([a.coeff[1,i] for i=1:ngens(C)])
    if isempty(primes)
      c=G([a.coeff[1,i] for i=ngens(C)+1:ngens(X)])
      return mC(b)*(pi\(mG(c)))
    else
      c=T([a.coeff[1,i] for i=ngens(C)+1:ngens(T)+ngens(C)])
      d=H([a.coeff[1,i] for i=ngens(T)+ngens(C)+1: ngens(X)])
      el=pi\(mG(c))
      # I need to modify $el$ so that it has the correct sign at the embeddings contained in primes
      vect=(lH(K(el))).coeff
      if vect==d.coeff
        return el*mC(b)
      else
        correction=O(1)
        for i=1:ngens(H)
          if d.coeff[1,i]==1
            correction=correction*eH(H[i])
          end
        end
        while vect!=d.coeff
          el=el+correction
          vect=(lH(K(el))).coeff
        end
        return el*mC(b)
      end
    end
  end

  mp=MapRayClassGrp{typeof(X)}()
  mp.header = Hecke.MapHeader(X, Hecke.NfOrdIdlSet(O), expo, disclog)
  mp.modulus_fin=m
  mp.modulus_inf=p

  return X, mp

end


text:
module TypesParticles

export Particle, Atom
export SystemOfAtoms

abstract type Particle end

"""
ab-initio Information about an atom.

Use: Atom(x::Float64,y::Float64,z::Float64,Z::Integer)
""" struct Atom <: Particle
           x::Float64
           y::Float64
           z::Float64
           Z::Integer
           function Atom(x::Real=0.00e0,
                         y::Real=0.00e0,
                         z::Real=0.00e0,
                         Z::Real=1)
                   x0,y0,z0 = map(Float64,(x,y,z));
                   Z0 = Int(Z);
                   if Z0 < 1
                       error("Atomic number < 1 is not allowed.")
                   end
                   new(x0,y0,z0,Z0);
           end
end #Atom

"""
Data type for storing list of information about an atomic system.

Use: SystemOfAtoms(natoms::Integer,atoms::Array{Atom,1},cell::Array{Float64,2}(3x3))
""" struct SystemOfAtoms
    natoms::Integer;
    atoms::Array{Atom,1};
    cell::Array{Float64,2};
    function SystemOfAtoms(natoms=1,
                    atoms=[Atom()],
                    cell=zeros(3,3))
            if natoms < 0
                error("Simulation cell cannot have negative atoms.")
            end
            if size(cell) != (3,3)
                error("Simulation cell is not a 3x3 array.")
            end
            if typeof(atoms) != Array{Atom,1}
                error("Atomic entries are not of the typeof() == Atom.")
            end
            new(natoms,atoms,cell)
    end
end #SystemOfAtoms

end#module


text:
using MPI
using StaticArrays
using CLIMA
using CLIMA.VariableTemplates
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging
using GPUifyLoops

import CLIMA.DGmethods:
    BalanceLaw,
    vars_aux,
    vars_state,
    vars_gradient,
    vars_diffusive,
    vars_integrals,
    flux_nondiffusive!,
    flux_diffusive!,
    source!,
    wavespeed,
    update_aux!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    boundary_state!,
    init_aux!,
    init_state!,
    init_ode_state,
    LocalGeometry

struct VarsTestModel{dim} <: BalanceLaw end

vars_state(::VarsTestModel, T) = @vars(x::T, coord::SVector{3, T})
vars_aux(m::VarsTestModel, T) = @vars(coord::SVector{3, T}, polynomial::T)
vars_diffusive(m::VarsTestModel, T) = @vars()

flux_nondiffusive!(::VarsTestModel, _...) = nothing
flux_diffusive!(::VarsTestModel, _...) = nothing
source!(::VarsTestModel, _...) = nothing
boundary_state!(_, ::VarsTestModel, _...) = nothing
wavespeed(::VarsTestModel, _...) = 1

function init_state!(m::VarsTestModel, state::Vars, aux::Vars, coord, t::Real)
    @inbounds state.x = coord[1]
    state.coord = coord
end

function init_aux!(
    ::VarsTestModel{dim},
    aux::Vars,
    g::LocalGeometry,
) where {dim}
    x, y, z = aux.coord = g.coord
    aux.polynomial = x * y + x * z + y * z
end

using Test
function run(mpicomm, dim, Ne, N, FT, ArrayType)

    brickrange = ntuple(j -> range(FT(0); length = Ne[j] + 1, stop = 3), dim)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, dim),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    dg = DGModel(
        VarsTestModel{dim}(),
        grid,
        Rusanov(),
        CentralNumericalFluxDiffusive(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    @test Array(Q.x)[:, 1, :] == Array(Q.coord)[:, 1, :]
    @test Array(dg.auxstate.coord) == Array(Q.coord)
    x = Array(Q.coord)[:, 1, :]
    y = Array(Q.coord)[:, 2, :]
    z = Array(Q.coord)[:, 3, :]
    @test Array(dg.auxstate.polynomial)[:, 1, :] ≈ x .* y + x .* z + y .* z
end

let
    CLIMA.init()
    ArrayType = CLIMA.array_type()

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 5, 5)
    lvls = 1

    polynomialorder = 4

    for FT in (Float64,) #Float32)
        for dim in 2:3
            err = zeros(FT, lvls)
            for l in 1:lvls
                @info (ArrayType, FT, dim)
                run(
                    mpicomm,
                    dim,
                    ntuple(j -> 2^(l - 1) * numelem[j], dim),
                    polynomialorder,
                    FT,
                    ArrayType,
                )
            end
        end
    end
end

nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
@testset "4.5 Secant" begin
    include("4.5.0 (a sec)^m (b trg)^n.jl")
    include("4.5.1.2 (d sec)^n (a+b sec)^m.jl")
    include("4.5.1.3 (d sin)^n (a+b sec)^m.jl")
    include("4.5.1.4 (d tan)^n (a+b sec)^m.jl")
    include("4.5.10 (c+d x)^m (a+b sec)^n.jl")
    include("4.5.11 (e x)^m (a+b sec(c+d x^n))^p.jl")
    include("4.5.2.1 (a+b sec)^m (c+d sec)^n.jl")
    include("4.5.2.3 (g sec)^p (a+b sec)^m (c+d sec)^n.jl")
    include("4.5.3.1 (a+b sec)^m (d sec)^n (A+B sec).jl")
    include("4.5.4.1 (a+b sec)^m (A+B sec+C sec^2).jl")
    include("4.5.4.2 (a+b sec)^m (d sec)^n (A+B sec+C sec^2).jl")
    include("4.5.7 (d trig)^m (a+b (c sec)^n)^p.jl")
end

max_stars_repo_path: test/4 Trig functions/4.5 Secant/.jl

max_issues_repo_path: test/4 Trig functions/4.5 Secant/.jl

max_forks_repo_path: test/4 Trig functions/4.5 Secant/.jl


text:
module FlatGeobuf

import FlatBuffers

# Partly autogenerated
include("schema/header.jl")
include("schema/feature.jl")

include("flatgeobuffer.jl")
include("index.jl")
include("io.jl")
include("table.jl")

export findall, read_file

end # module


text:
## CSV

using .CSV

read_table(::TableIOInterface.CSVFormat, filename:: AbstractString; kwargs...) = CSV.File(filename; kwargs...)
read_table(::TableIOInterface.CSVFormat, io:: IO; kwargs...) = CSV.File(read(io); kwargs...)

function write_table!(::TableIOInterface.CSVFormat, output:: Union{AbstractString, IO}, table; kwargs...)
    _checktable(table)
    table |> CSV.write(output; kwargs...)
    nothing
end


text:
## print.jl : useful print statements
#
# A collection of useful print statements for Estimators.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

#
# convenience functions
#
l() = 81

l_short() = 29

n() = 15

spaces(n) = repeat(" ", n)

short(num) = @sprintf("%5.3f", num)

shorte(num) = @sprintf("%5.3e", num)

long(num) = @sprintf("%12.5e", num)

end_table_row(str) = println(string(str, spaces(l()-length(str)-1), "|"))

hline() = println(string("+", repeat("-", l()-2), "+"))

table_hline(m) = println(string("+", repeat(string(repeat("-", n()), "+"), m)))

#
# header and footer
#
function print_header(estimator::Estimator, ϵ::Real)
    hline()
    str = string("| *** MultilevelEstimators.jl @", now())
    end_table_row(str)
    str = string("| *** This is a ", estimator)
    end_table_row(str)
    str = string("| *** Simulating ", first(split(estimator[:name], ".")))
    end_table_row(str)
    str = string("| *** Tolerance on RMSE ϵ = ", shorte(ϵ))
    end_table_row(str)
    hline()
end

function print_footer()
    hline()
    str = string("| *** MultilevelEstimators.jl @", now())
    end_table_row(str)
    str = "| *** Successfull termination."
    end_table_row(str)
    hline()
end

#
# status
#
function print_status(estimator::Estimator)
    table_hline(7)
    header = "| "
    for name in [print_elname(estimator) "E" "dE" "V" "dV" "N" "W"]
        header = string(header, name, spaces(n()-length(name)-1), "| ")
    end
    println(header)
    table_hline(7)
    for index in keys(estimator)
		if nb_of_samples(estimator, index) > 0
			index_str = string(index)
			str = "| "
			str = string(str, index_str, spaces(n()-length(index_str)-2), " |")
			str = string(str, long(mean0(estimator, index)), spaces(n()-12-1), " |")
			str = string(str, long(mean(estimator, index)), spaces(n()-12-1), " |")
			str = string(str, long(var0(estimator, index)), spaces(n()-12-1), " |")
			str = string(str, long(var(estimator, index)), spaces(n()-12-1), " |")
			samples_str = print_nb_of_samples(estimator, index)
			str = string(str, " ", samples_str, spaces(n()-length(samples_str)-2), " |")
			str = string(str, long(cost(estimator, index)), spaces(n()-12-1), " |")
			println(str)
		end
    end
    table_hline(7)
end

#
# index or level
#
print_elname(::Estimator{<:SL}) = "level"

print_elname(::Estimator{<:AbstractML}) = "level"

print_elname(::Estimator{<:AbstractMI}) = "index"

#
# optimal number of samples
#
function print_optimal_nb_of_samples(estimator::Estimator, samples)
    println("Samples will be updated ", estimator isa U ? "with" : "to")
    table_hline(2)
    header = "| "
    for name in [print_elname(estimator) "N"]
        header = string(header, name, spaces(n()-length(name)-1), "| ")
    end
    println(header)
    table_hline(2)
    _print_optimal_nb_of_samples(estimator, samples)
    table_hline(2)
end

function _print_optimal_nb_of_samples(estimator::Estimator, samples::Dict)
    for index in sort(collect(keys(samples)))
		if samples[index] > 0
			index_str = string(index)
			str = "| "
			str = string(str, index_str, spaces(n()-length(index_str)-2), " |")
			samples_str = print_nb_of_samples(estimator, index, samples[index])
			str = string(str, " ", samples_str, spaces(n()-length(samples_str)-2), " |")
			println(str)
		end
    end
end

print_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index, n::Integer) = string(n)

print_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index, n::Integer) = string(estimator[:nb_of_shifts](index), " x ", n)

print_nb_of_samples(estimator::Estimator, index::Index) = print_nb_of_samples(estimator, index, nb_of_samples(estimator, index))

#
# profits
#
function print_largest_profit(estimator, max_index, max_profit, indices, profits)
    println("Profit indicators:")
    table_hline(2)
    header = "| "
    for name in [print_elname(estimator) "profit"]
        header = string(header, name, spaces(n()-length(name)-1), "| ")
    end
    println(header)
    table_hline(2)
    for i in 1:length(indices)
        index_str = string(indices[i])
        str = "| "
        str = string(str, index_str, spaces(n()-length(index_str)-2), " |")
        profit_str = shorte(profits[i])
        str = string(str, " ", profit_str, spaces(n()-length(profit_str)-2), " |")
        println(str)
    end
    table_hline(2)
    println("Index with max profit is ", max_index, " (value = ", long(max_profit), ").")
end

#
# probabilities
#
function print_pmf(estimator::Estimator)
	P = pmf(estimator)
    println("Using approximate probability mass function:")
    table_hline(2)
    header = "| "
	for name in [print_elname(estimator) "P"]
        header = string(header, name, spaces(n()-length(name)-1), "| ")
    end
    println(header)
    table_hline(2)
	for idx in sort(collect(keys(P)))
        index_str = string(idx)
        str = "| "
        str = string(str, index_str, spaces(n()-length(index_str)-2), " |")
        pmf_str = shorte(P[idx])
        str = string(str, " ", pmf_str, spaces(n()-length(pmf_str)-2), " |")
        println(str)
    end
    table_hline(2)
end

#
# convergence
#
function print_convergence(estimator::Estimator, converged::Bool)
    print_status(estimator)
    converged && println(string("Convergence reached. RMSE ≈", long(rmse(estimator)), "."))
    print_footer()
end

function print_mse_analysis(estimator::Estimator, ϵ::Real, θ::Real)
    println("Checking convergence...")
    print_rates(estimator)
    println(string("  ==> Variance of the estimator ≈", long(varest(estimator)), "."))
    println(string("  ==> Bias of the estimator ≈", long(bias(estimator)), "."))
    println(string("  ==> MSE splitting parameter ≈ ", short(θ), "."))
    if !converged(estimator, ϵ, θ) > ϵ
        println(string("No convergence yet. RMSE ≈", long(rmse(estimator)), " > ", shorte(ϵ), "."))
        println("Adding an extra level...")
    end
end

function print_qmc_convergence(estimator::Estimator{<:AbstractIndexSet, <:QMC}, ϵ::Real, θ::Real)
    println(string("Checking if variance of estimator is larger than θ*ϵ^2..."))
    converged = !(varest(estimator) > θ*ϵ^2)
    println(string("  ==> ", converged ? "no!" : "yes", " (", long(varest(estimator))[2:end], converged ? " < " : " > ", short(θ), " *", long(ϵ^2), ")"))
end

function print_unbiased_convergence(estimator::Estimator{<:U}, ϵ::Real)
    println(string("Checking if variance of estimator is larger than ϵ^2..."))
    converged = !(varest(estimator) > ϵ^2)
    println(string("  ==> ", converged ? "no!" : "yes", " (", long(varest(estimator))[2:end], converged ? " < " : " > ", long(ϵ^2)[2:end], ")"))
    print_rates(estimator)
end

#
# taking ...
#
print_sample!_header(estimator::Estimator, index::Index, n::Integer, warm_up::Bool) = print(string("Taking ", print_nb_of_samples(estimator, index, n), print_warm_up(Val(warm_up)), " sample", print_with_s(n), " at ", print_elname(estimator), " ", index, "..."))

print_warm_up(::Val{true}) = " warm-up"

print_warm_up(::Val{false}) = " additional"

print_with_s(n::Integer) = n > 1 ? "s" : ""

print_sample!_footer() = println(" done")

#
# rates
#
function print_rates(estimator::Estimator)
    str = string("  ==> Rates: α ≈ ", print_rate(estimator, α))
    str = string(str, ", β ≈ ", print_rate(estimator, β))
    str = string(str, ", γ ≈ ", print_rate(estimator, γ), ".")
    println(str)
end

print_rate(estimator::Estimator{<:AbstractML}, f::Function) = short(f(estimator)[1])

print_rate(estimator::Estimator{<:AbstractMI}, f::Function) = string("(", join(short.(f(estimator)), ", "), ")")

#
# max level/index
#
warn_max_level(estimator::Estimator) = @warn string("Maximum ", _warn_max_level_name(estimator), " L = ", estimator[:max_index_set_param], " reached, no convergence yet.")

_warn_max_level_name(estimator::Estimator{<:AbstractML}) = "level"

_warn_max_level_name(estimator::Estimator{<:AbstractMI}) = "index set parameter"

warn_max_index(estimator::Estimator{<:AD}, index) = @warn string("Trying to add index ", index, " outside search space to index set, ignoring...")

#
# level parameter
#
print_level(estimator::Estimator{<:SL}, L) = println(string("Currently running on finest level."))

print_level(estimator::Estimator{<:ML}, L) = println(string("Currently running on level ", L, "."))

print_level(estimator::Estimator, L) = println(string("Currently running with L = ", L, "."))

#
# index set
#
print_index_set(estimator::Estimator, index_set) = nothing

print_index_set(estimator::Estimator{<:AbstractMI}, index_set) = begin
    if ndims(estimator) == 2
        println("Shape of the index set:")
        print(union(keys(estimator), index_set))
    end
end

print_index_set(estimator::Estimator{<:U}, index_set) = begin
    if ndims(estimator) == 2
        println("Shape of the index set:")
		print(filter(i -> has_samples_at_index(estimator, i), keys(estimator)))
    end
end

print_index_set(estimator::Estimator{<:AD}, index_set) = begin
    if ndims(estimator) == 2
        println("Shape of the index set (\u25A3 = active set):")
        print(collect(old_set(estimator)), collect(active_set(estimator)))
    end
end