text:
module CuArrays

using CUDAapi, CUDAdrv, CUDAnative

using GPUArrays

export CuArray, CuVector, CuMatrix, CuVecOrMat, cu

import LinearAlgebra

using Adapt

using Requires

## discovery

let
    toolkit = find_toolkit()

    # required libraries that are part of the CUDA toolkit
    for name in ("cublas", "cusparse", "cusolver", "cufft", "curand")
        lib = Symbol("lib$name")
        path = find_cuda_library(name, toolkit)
        if path === nothing
            error("Could not find library '$name' (it should be part of the CUDA toolkit)")
        end
        Base.include_dependency(path)
        @eval global const $lib = $path
    end

    # optional libraries
    for name in ("cudnn", "cutensor")
        lib = Symbol("lib$name")
        path = find_cuda_library(name, toolkit)
        if path !== nothing
            Base.include_dependency(path)
        end

        # provide a global constant that returns the path to the library,
        # or nothing if the library is not available (for use in conditional expressions)
        @eval global const $lib = $path

        # provide a macro that either returns the path to the library,
        # or a run-time error if the library is not available (for use in ccall expressions)
        exception = :(error($"Your installation does not provide $lib, CuArrays.$(uppercase(name)) is unavailable"))
        @eval macro $lib() $lib === nothing ? $(QuoteNode(exception)) : $lib end

        # provide a function for external use (a la CUDAapi.has_cuda)
        fn = Symbol("has_$name")
        @eval (export $fn; $fn() = $lib !== nothing)
    end
end

## source code includes

include("memory.jl")
include("array.jl")
include("subarray.jl")
include("permuteddimsarray.jl")
include("utils.jl")
include("indexing.jl")
include("broadcast.jl")
include("matmul.jl")
include("mapreduce.jl")
include("accumulate.jl")

include("gpuarray_interface.jl")

# many libraries need to be initialized per-device (per-context, really, but we assume users
# of CuArrays and/or CUDAnative only use a single context), so keep track of the active one.
const active_context = Ref{CuContext}()

include("blas/CUBLAS.jl")
include("sparse/CUSPARSE.jl")
include("solver/CUSOLVER.jl")
include("fft/CUFFT.jl")
include("rand/CURAND.jl")
include("dnn/CUDNN.jl")
include("tensor/CUTENSOR.jl")

include("nnlib.jl")

include("deprecated.jl")

## initialization
function __init__()
    # package integrations
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")

    # update the active context when we switch devices
    callback = (::CuDevice, ctx::CuContext) -> begin
        active_context[] = ctx

        # wipe the active handles
        CUBLAS._handle[] = C_NULL
        CUBLAS._xt_handle[] = C_NULL
        CUSOLVER._dense_handle[] = C_NULL
        CUSOLVER._sparse_handle[] = C_NULL
        CUSPARSE._handle[] = C_NULL
        CURAND._generator[] = nothing
        CUDNN._handle[] = C_NULL
        CUTENSOR._handle[] = C_NULL
    end
    push!(CUDAnative.device!_listeners, callback)

    # a device might be active already
    existing_ctx = CUDAdrv.CuCurrentContext()
    if existing_ctx !== nothing
        active_context[] = existing_ctx
    end

    __init_memory__()
end

end # module


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Interpreter utility functions

"""
    read_at_line(io::IO, line_num::Int)

Read the `line_num`-th line of `io`.
"""
function read_at_line(io::IO, line_num::Int)
    seekstart(io)
    for _ ∈ 1:line_num - 1
        readline(io)
    end
    readline(io)
end

"""
    read_at_line(file_name::String, line_num::Int)

Read the `line_num`-th line of a file named `file_name`.
"""
function read_at_line(file_name::String, line_num::Int)
    open(file_name, "r") do io
        read_at_line(io, line_num)
    end
end

"""
    isnewline(c::Char)

Check if `c` is a newline character.
"""
isnewline(c::Char) = c ∈ ('\n', '\r')

"""
    issymbol(c::Char)

Check if `c` is in a partucular set of characters used in a SceneLang script.
"""
issymbol(c::Char) = c ∈ "{}()<>[],*"

"""
    valid_operations

A dictionary storing the operation `Symbol` and a function `::Int -> ::Bool` that return true when the number of arguments is under a threshold.

For example the pair `:+ => ( ::Int) -> true` indicates that the operation `+` is valid with any number of arguments,
while the pair `:floor => (n::Int) -> n == 1` indicates that the operation `floor` is valid only if it has one argument.

See also: [`Raytracer.isvalid`](@ref)
"""
const valid_operations = Dict(:+     => (; f = ( ::Int) -> true,        signature = "+(x...)"),
                              :-     => (; f = (n::Int) -> 1 <= n <= 2, signature = "-(x, [y])"),
                              :*     => (; f = ( ::Int) -> true,        signature = "*(x...)"),
                              :/     => (; f = (n::Int) -> n == 2,      signature = "/(x, y)"),
                              :%     => (; f = (n::Int) -> n == 2,      signature = "%(x, y)"),
                              :^     => (; f = (n::Int) -> n == 2,      signature = "^(x, y)"),
                              :div   => (; f = (n::Int) -> n == 2,      signature = "div(x, y)"),
                              :floor => (; f = (n::Int) -> n == 1,      signature = "floor(x)"),
                              :ceil  => (; f = (n::Int) -> n == 1,      signature = "ceil(x)"),
                              :round => (; f = (n::Int) -> n == 1,      signature = "round(x)"),
                              :exp   => (; f = (n::Int) -> n == 1,      signature = "exp(x)"),
                              :exp2  => (; f = (n::Int) -> n == 1,      signature = "exp2(x)"),
                              :exp10 => (; f = (n::Int) -> n == 1,      signature = "exp10(x)"),
                              :log   => (; f = (n::Int) -> n == 1,      signature = "log(x)"),
                              :log2  => (; f = (n::Int) -> n == 1,      signature = "log2(x)"),
                              :log10 => (; f = (n::Int) -> n == 1,      signature = "log10(x)"),
                              :log1p => (; f = (n::Int) -> n == 1,      signature = "log1p(x)"),
                              :sin   => (; f = (n::Int) -> n == 1,      signature = "sin(x)"),
                              :cos   => (; f = (n::Int) -> n == 1,      signature = "cos(x)"),
                              :tan   => (; f = (n::Int) -> n == 1,      signature = "tan(x)"),
                              :asin  => (; f = (n::Int) -> n == 1,      signature = "asin(x)"),
                              :acos  => (; f = (n::Int) -> n == 1,      signature = "acos(x)"),
                              :atan  => (; f = (n::Int) -> 1 <= n <= 2, signature = "atan(x, [y])"),
                              :Point => (; f = (n::Int) -> n == 3,      signature = "Point(x, y, z)"),
                              :RGB   => (; f = (n::Int) -> n == 3,      signature = "RGB(r, g, b)"),
                             )

"""
    Raytracer.isvalid(expr::Expr, str_len::Int, token_location::SourceLocation)

Return `true` if the given expression is valid in a SceneLang script, else throw an appropriate exception using `str_len` and `token_location`.

See also: [`valid_operations`](@ref)
"""
function Raytracer.isvalid(expr::Expr, str_len::Int, token_location::SourceLocation)
    expr.head == :call ||
        throw(InvalidExpression(token_location, "Invalid mathematical expression: expression head is not a call", str_len + 1))
    op_name = expr.args[begin]
    isa(op_name, Symbol) ||
        throw(InvalidExpression(token_location, "Invalid mathematical expression: operation name '$op_name' is not a symbol" , str_len + 1))
    op_name ∈ keys(valid_operations) ||
        throw(InvalidExpression(token_location, "Invalid mathematical expression: contains invalid operation $(expr.args[begin])\nValid operations are: " * join(valid_operations, ", "), str_len + 1))
    (invalid = findfirst(arg -> !isa(arg, Union{Integer, AbstractFloat, Expr, Symbol}), expr.args[begin + 1:end])) |> isnothing ||
        throw(InvalidExpression(token_location, "Invalid mathematical expression: contains invalid operand $(expr.args[invalid + 1])\nValid operands are instances of `Integer`, `AbstractFloat`, `Symbol` or `Expr`", str_len + 1))
    valid_operations[op_name].f(length(expr.args) - 1) ||
        throw(InvalidExpression(token_location, "Invalid mathematical expression: operation signature is '$(valid_operations[op_name].signature)', got '$(op_name)' with $(length(expr.args[begin+1:end])) arguments", str_len + 1))

    return all(arg -> (isa(arg, Expr) ? Raytracer.isvalid(arg, str_len, token_location) : true), expr.args[begin + 1:end])
end

ext: jl


text:
module Losses

using Statistics
using Zygote
using Zygote: @adjoint
using ..Flux: ofeltype, epseltype
using CUDA
using NNlib: logsoftmax, logσ
import Base.Broadcast: broadcasted

export mse, mae,nll, msle,softplus,logcosh,
    crossentropy, logitcrossentropy,
    # binarycrossentropy, logitbinarycrossentropy # export only after end deprecation
    kldivergence,
    huber_loss,
    tversky_loss,
    dice_coeff_loss,
    poisson_loss,
    hinge_loss, squared_hinge_loss

include("utils.jl")
include("functions.jl")

end #module


text:
# =============================================================================
# Spiking neuron
#
# system type: hybrid system with polynomial dynamics
# state dimension: 2
# modes: 1
# transitions: 1
#
# This model is adapted from Izhikevich's model [1].
# In [2, 3] this model has been formalized as a hybrid automaton.
#
# [1] Eugene M. Izhikevich: Dynamical systems in neuroscience. MIT press, 2007.
# [2] https://flowstar.org/examples/
# [3] Xin Chen: Reachability analysis of non-linear hybrid systems using Taylor
# models. PhD thesis, 2015.
# =============================================================================
using HybridSystems, MathematicalSystems, LazySets, TaylorIntegration

@taylorize function flow!(dx, x, params, t)
    local a = 0.02
    local b = 0.2
    local I = 40
    dx[1] = (0.04*(x[1]*x[1]) + 5*x[1]) + ((140+I)-x[2])
    dx[2] = a*((b*x[1]) - x[2])
    return dx
end

function spikingNeuron_model()
    # hybrid automaton
    HA = LightAutomaton(1)

    # mode 1
    X = HPolyhedron([HalfSpace([1.0, 0.0], 30.0)])  # x1 ≤ 30
    m1 = @system(x' = flow!(x), dim: 2, x ∈ X)

    # transition mode 1 → mode 1 (self loop)
    add_transition!(HA, 1, 1, 1)
    G = HPolyhedron([HalfSpace([-1.0, 0.0], -30.0)])  # x1 ≥ 30
    A = [0.0 0.0; 0.0 1.0]
    b = [-65.0, 8.0]
    R11 = ConstrainedAffineMap(A, b, G)  # x1 := -65, x2 := x2 + 8

    # hybrid system
    S = HybridSystem(HA, [m1], [R11], [AutonomousSwitching()])

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# Translated from the Golang Tour's select intro: https://tour.golang.org/concurrency/6
module Default
using Select

doevery(secs) = Channel{Int}() do ch
    while true
        sleep(secs)
        put!(ch, 1)
    end
end
doafter(secs) = Channel{Int}() do ch
    sleep(secs)
    put!(ch, 1)
end

function main()
    tick = doevery(0.1)
    boom = doafter(0.5)
    while true
        @select begin
            tick  => println("tick.")
            boom  => begin
                println("BOOM!")
                return
            end
            _  => begin
                println("    .")
                sleep(0.05)
            end
        end
    end
end

main()

end


text:
using Plots
using LinearAlgebra
using Polynomials

using SparseArrays
using Arpack

using ForwardDiff

using BenchmarkTools

#define integration functions
function Lagrange(x,xi,xj)
           p = poly(xj[:])
           p = p/p(xi)
           y = p(x)
           return y
end

function H(x,xi,xj)
    Lval(y) = Lagrange(y,xi,xj)
    Ldot = ForwardDiff.derivative(Lval,xi)
    H = (1-2*(x-xi)*Ldot)*Lagrange(x,xi,xj)^2
    return H
end
function Hbar(x,xi,xj)
    Hbar = (x-xi)*Lagrange(x,xi,xj)^2
    return Hbar
end
#produces (w,p)
function gauss(n)
    w = zeros(n,1)
    global coeffs = zeros(n+1,1);
    for k = Int64.(0:floor(n/2))
        coeffs[2*k+1] = (-1)^k*binomial(n,k)*binomial(2*n-2*k,n)/(2^n);
    end
    coeffs = Poly(reverse(coeffs,dims=1)[:])
    points = roots(coeffs);

    for i=1:n
        L = poly(points[[1:i-1; i+1:end]]);
        L = L/polyval(L,points[i]);
        L = polyint(L);
        w[i] = polyval(L,1)-polyval(L,-1);

    end
    #I = sum(w(i).*f(points));
    return w,points
end
gaussw,gausspoints = gauss(6) #used to integrate by gaussintegral
#results I
function gaussintegral(x1,x2,w,points,f)
    g(y) = 0.5*(x2-x1)*f.((x2-x1)*y/2 +(x2+x1)/2)
    I = sum(w.*g.(points))
    return I
end

function quad(ksi,eta,xs,ys)
    #x and y must be line vectors!!!!
    #returns x,y on selected ksi, eta points
    n = length(xs)
    xs = reshape(xs,(n,n))'
    ys = reshape(ys,(n,n))'
    nξ = n; nη = n;
    ξ = LinRange(-1,1,nξ); η = LinRange(-1,1,nη)
    Nij(i,j,ξi,ηi) = Lagrange(ξi,ξ[i],ξ[deleteat!(collect(1:nξ),i)])*Lagrange(ηi,η[j],η[deleteat!(collect(1:nη),j)])
    x = 0
    y = 0
    for i=1:n
        for j=1:n
            x = x + xs[i,j]*Nij(i,j,ksi,eta)
            y = y + ys[i,j]*Nij(i,j,ksi,eta)
        end
    end
    return x,y
end
mutable struct Foo2
    x::Int64
end

function setx(a::Foo2, v)
 a.x=v
end
#nodesξη =reverse(collect(Iterators.product(ξ, η)),dims=1)
#=
ξη
(-1.0, -1.0)       (-1.0, -0.333333)       (-1.0, 0.333333)       (-1.0, 1.0)
(-0.333333, -1.0)  (-0.333333, -0.333333)  (-0.333333, 0.333333)  (-0.333333, 1.0)
(0.333333, -1.0)   (0.333333, -0.333333)   (0.333333, 0.333333)   (0.333333, 1.0)
(1.0, -1.0)        (1.0, -0.333333)        (1.0, 0.333333)        (1.0, 1.0)
=#
#=
ij
-------------------------------
|(1, 1)  (1, 2)  (1, 3)  (1, 4)|
|(2, 1)  (2, 2)  (2, 3)  (2, 4)|
|(3, 1)  (3, 2)  (3, 3)  (3, 4)|
|(4, 1)  (4, 2)  (4, 3)  (4, 4)|
-------------------------------
=#

max_stars_repo_path: axnai julia/2d_tria_mess.jl

max_stars_repo_name: billyon/first

max_issues_repo_path: axnai julia/2d_tria_mess.jl

max_issues_repo_name: billyon/first

max_forks_repo_path: axnai julia/2d_tria_mess.jl

max_forks_repo_name: billyon/first


text:
# Implementations for Base/stdlib functions

isfrozenindex(::Type{<:Union{Array,SubArray,Base.ReshapedArray,BitArray}}) = true
isfrozenindex(::Type{<:Vector}) = false
isfrozenindex(::Type{<:BitArray}) = false


text:
module test_loltools_spectatorv4

using Test
using LOLTools.SpectatorV4

include("helpers.jl")

struct SpectatorController <: ApplicationController
    conn::Conn
end

function featured_games(c::SpectatorController)
    render(JSON, (a=1,))
end

routes() do
    get("/lol/spectator/v4/featured-games", SpectatorController, featured_games)
end

api_key = ""
platform = "kr"
SpectatorV4.featured_games(api_key, platform; action=mock_action)

end # module test_loltools_spectatorv4


text:
module Named

using Reexport
@reexport using ..LambdaCalculus

import Base: convert, show
import ..LambdaCalculus: boundvartype, freevartype, reify

export Term,
    Var,
    App,
    Abs

"""Named lambda terms, built using only the following rule:

    <Term> := <Name>              (variable)
            | λ <Name> . <Term>   (abstraction)
            | (<Term> <Term>)     (application)

"""
abstract type Term <: AbstractTerm end

struct Var <: Term
    name::Symbol
end

struct Abs <: Term
    boundname::Symbol
    body::Term
end

struct App <: Term
    car::Term
    cdr::Term
end

show(io::IO, t::Abs) = print(io, "(λ", t.boundname, ".", t.body, ")")
show(io::IO, t::App) =  print(io, "(", t.car, " ", t.cdr, ")")
show(io::IO, t::Var) = print(io, t.name)

# Conversions from expressions to terms
convert(::Type{Var}, v::Symbol) = Var(v)
convert(::Type{Term}, v::Symbol) = convert(Var, v)

function convert(::Type{App}, expr::Expr)
    @assert(length(expr.args) >= 2, "call must contain arguments")
    mapfoldl(e -> convert(Term, e), (e, arg) -> App(e, arg), expr.args)
end

function convert(::Type{Abs}, expr::Expr)
    @assert(isa(expr.args[1], Symbol), "only single-argument lambdas are allowed")
    # TODO: handle multiple arguments
    Abs(expr.args[1], convert(Term, expr.args[2]))
end

function convert(::Type{Term}, expr::Expr)
    if expr.head == :call
        # TODO: handle :* case
        convert(App, expr)
    elseif expr.head == :->
        convert(Abs, expr)
    elseif expr.head == :block
        # Such trivial blocks are used by the parser in lambdas to add metadata.
        # They consist of a LineNumberNode followed by the actual expression.
        convert(Term, expr.args[end])
    else
        @error "unhandled syntax: $expr"
    end
end

# Conversions from terms to expressions
convert(::Type{Expr}, v::Var) = :($(v.name))
convert(::Type{Expr}, t::Abs) = :($(t.boundname) -> $(convert(Expr, t.body)))
convert(::Type{Expr}, t::App) = :(($(convert(Expr, t.car))($(convert(Expr, t.cdr)))))

reify(v::Var) = :(Var($(Meta.quot(v.name))))
reify(t::Abs) = :(Abs($(Meta.quot(t.boundname)), $(reify(t.body))))
reify(t::App) = :(App($(reify(t.car)), $(reify(t.cdr))))

freevartype(::Type{<:Term}) = Var
boundvartype(::Type{<:Term}) = Var

include("syntactic.jl")
include("meta.jl")

end # module Named


text:
function demcz_sample(logobj, Zmat, N=4, K=10, Ngeneration=5000, Nblocks=1, blockindex=[1:size(Zmat,2)], eps_scale=1e-4*ones(size(Zmat,2)), γ=2.38; prevrun=nothing, verbose = true, print_step=100)
    M, d = size(Zmat)
    if prevrun == nothing
        # start chain at last N elements of Z
        X = Zmat[end-N+1:end, :]
        # initial obj values
        log_objcurrent = map(logobj, [X[i,:] for i = 1:N])
    else
        # start chain at last generation of prevrun
        X = prevrun.chain[:,:,end]
        log_objcurrent = prevrun.log_objcurrent[:,end]
    end
    # preallocate the chain
    mc = MC(zeros(N,  d, Ngeneration),zeros(N, Ngeneration), X, log_objcurrent)

    # print inital values
    if verbose
        print_status(mc, 0)
    end
    # run through generations
    for ig = 1:Ngeneration
        for ic = 1:N
            Xcurrent, current_logobj = update_blocks(mc.Xcurrent[ic, :], mc.log_objcurrent[ic], Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
            # update in chain
            mc.chain[ic, :, ig] .= Xcurrent
            mc.log_obj[ic, ig] = current_logobj
            mc.Xcurrent[ic, :] .= Xcurrent
            mc.log_objcurrent[ic] = current_logobj
        end
        if mod(ig, K) == 0.
            Zmat = vcat(Zmat, mc.Xcurrent)
            M += N
        end
        if verbose
            if mod(ig, print_step) == 0.
                print_status(mc, ig)
            end
        end
    end

    if prevrun != nothing
        return MC(cat(prevrun.chain, mc.chain, dims=3),cat(prevrun.log_obj, mc.log_obj, dims=2), mc.Xcurrent, mc.log_objcurrent), Zmat
    else
        return mc, Zmat
    end
end

function print_status(mc, ig)
    #bestval = maximum(mc.log_objcurrent)
    #bestpar = mc.Xcurrent[findfirst(bestval.==mc.log_objcurrent), :]
    avgpar = mean(mc.chain, dims = [1,3])
    avgval = mean(mc.log_obj, dims = [1,2])

    println("-----------------------")
    println("iteration $ig")
    println("average par = $avgpar")
    println("average val = $avgval")
    #println("bestval = $bestval")
    #println("bestpar = $bestpar")
    println("-----------------------")
end
#
# function demcz_sample_par(logobj, Zmat, N, K, Ngeneration, Nblocks, blockindex, eps_scale, γ; verbose=true, print_to_file=false, file = "./demcstat.txt")
#     wp = CachingPool(workers())
#     Mval, d = size(Zmat)
#     X = Zmat[end-N+1:end, :]
#     log_objcurrent = pmap(wp, logobj, [X[i,:] for i = 1:N])
#     mc = MC(Array{Float64}(N,  d, Ngeneration),Array{Float64}(N, Ngeneration), X, log_objcurrent)
#     global Xcurrent = copy(mc.Xcurrent)
#     global log_objcurrent = copy(mc.log_objcurrent)
#     global Z = Zmat
#     global M = Mval
#     @everywhere global Z
#     @everywhere global M
#     @everywhere global Xcurrent
#     @everywhere global log_objcurrent
#
#     passobj(myid(), workers(), [:Xcurrent, :log_objcurrent], from_mod = DEMC, to_mod = DEMC)
#     passobj(myid(), workers(), [:Z, :M], from_mod = DEMC, to_mod = DEMC)
#
#     if verbose
#         bestval = maximum(mc.log_objcurrent)
#         bestpar = mc.Xcurrent[findfirst(bestval.==mc.log_objcurrent), :]
#         if print_to_file == true
#             f = open(file, "w")
#             println(f, "iteration 0")
#             println(f, "bestval = $bestval")
#             println(f, "bestpar = $bestpar")
#             close(f)
#         else
#             println("iteration 0")
#             println("bestval = $bestval")
#             println("bestpar = $bestpar")
#         end
#     end
#     for ig = 1:Ngeneration
#         passobj(myid(), workers(), [:Xcurrent, :log_objcurrent], from_mod = DEMC, to_mod = DEMC)
#         res = pmap(wp, ic -> update_blocks(Xcurrent[ic,:], log_objcurrent[ic], Z, M, logobj, blockindex, eps_scale, γ, Nblocks), 1:N)
#         for ic = 1:N
#             # update in chain
#             mc.chain[ic, :, ig] = res[ic][1]
#             mc.log_obj[ic, ig] = res[ic][2]
#             mc.Xcurrent[ic, :] = res[ic][1]
#             mc.log_objcurrent[ic] = res[ic][2]
#             Xcurrent[ic, :] = res[ic][1]
#             log_objcurrent[ic] = res[ic][2]
#         end
#
#         if mod(ig, K) == 0.
#             Z = vcat(Z, mc.Xcurrent)
#             M += N
#             passobj(myid(), workers(), [:Z, :M], from_mod = DEMC, to_mod = DEMC)
#             if verbose
#                 bestval = maximum(mc.log_objcurrent)
#                 bestpar[:] = mc.Xcurrent[findfirst(bestval.==mc.log_objcurrent), :]
#                 avglast100 = mean(mc.chain[:, :, max(1, ig-100):ig], [1,3])
#                 avglast250 = mean(mc.chain[:, :, max(1, ig-250):ig], [1,3])
#                 accept_ratio = mean(mc.log_obj[:,max(1, ig-250):ig-1].!=mc.log_obj[:,max(2, ig-249):ig],2)
#                 accept_ratio100 = mean(mc.log_obj[:,max(1, ig-100):ig-1].!=mc.log_obj[:,max(2, ig-99):ig],2)
#                 if print_to_file
#                     f = open(file, "w")
#                     println(f, "iteration $ig")
#                     println(f, "bestval = $bestval")
#                     println(f, "bestpar = $bestpar")
#                     println(f, "avg last 100 = $avglast100")
#                     println(f, "accept_ratio last 100 = $accept_ratio100")
#                     println(f, "avg last 250 = $avglast250")
#                     println(f, "accept_ratio last 250 = $accept_ratio")
#                     close(f)
#                 else
#                     println("iteration $ig")
#                     println("bestval = $bestval")
#                     println("bestpar = $bestpar")
#                     println("avg last 100 = $avglast100")
#                     println("accept_ratio last 100 = $accept_ratio100")
#                     println("avg last 250 = $avglast250")
#                     println("accept_ratio last 250 = $accept_ratio")
#                 end
#             end
#         end
#     end
#     return mc
# end

function update_blocks(Xcurrent, current_logobj, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    for ib in 1:Nblocks
        Xcurrent, current_logobj = update_demcz_chain_block(Xcurrent, current_logobj, ib, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    end
    return Xcurrent, current_logobj
end

function update_demcz_chain_block(Xcurrent, current_logobj, ib, Zmat, M, logobj, blockindex, eps_scale, γ, Nblocks)
    # generate proposal
    set = collect(1:M)
    i1 = rand(set)
    deleteat!(set, i1)
    i2 = rand(set)
    de_diffvec = zeros(eltype(Xcurrent), length(Xcurrent))
    block = blockindex[ib]
    blocklen = length(block)
    if blocklen == 1
        de_diffvec[block] = γ*(Zmat[i1,block]- Zmat[i2,block]) .+ eps_scale[block] .* randn()
    else
        de_diffvec[block] = γ/sqrt(2*blocklen)*(Zmat[i1,block]-Zmat[i2,block]) .+ eps_scale[block] .* randn(blocklen)
    end
    Xproposal = Xcurrent .+ de_diffvec
    log_objXprop = logobj(Xproposal)
    if accept(log_objXprop, current_logobj)
        return Xproposal, log_objXprop
    else
        return Xcurrent, current_logobj
    end
end

function accept(objprop, objprev)
    if log(rand()) < objprop - objprev
        return true
    else
        return false
    end
end


text:
module PkgCmd

using IonBase

include(joinpath(pkgdir(IonBase), "test", "utils.jl"))

with_test_ion() do
    cd(joinpath(test_dir, "Basic")) do
        IonBase.add("OhMyREPL")
        IonBase.rm("OhMyREPL")
        IonBase.instantiate()
        IonBase.status()
    end
end

end


text:
using InvariantMeasures
using ValidatedNumerics

m = 30
m_extend = 100

D = ApproxInducedLSV(0.5, 15)
B = Ulam(1024)
Q = DiscretizedOperator(B, D)

normQ = opnormbound(weak_norm(B), Q)

trivial_norms = norms_of_powers_trivial(weak_norm(B), Q, m)
computed_norms = norms_of_powers(weak_norm(B), m, Q, integral_covector(B))

(dfly_strongs, dfly_norms) = norms_of_powers_dfly(B, D, m)

norms = min.(trivial_norms, computed_norms, dfly_norms) # in the current version, dfly_norms are always larger and can be omitted
@info norms

better_norms = refine_norms_of_powers(norms, m_extend)

w = invariant_vector(B, Q)
@show distance_from_invariant(B, D, Q, w, better_norms)

B_fine = Ulam(8192)
Q_fine = DiscretizedOperator(B_fine, D)
norm_Q_fine = opnormbound(weak_norm(B_fine), Q_fine)

trivial_norms_fine = norms_of_powers_trivial(weak_norm(B_fine), Q_fine, m_extend)
twogrid_norms_fine = norms_of_powers_from_coarser_grid(B_fine, B, D, better_norms, norm_Q_fine)

(dfly_strongs_fine, dfly_norms_fine) = norms_of_powers_dfly(B_fine, D, m_extend)

norms_fine = min.(trivial_norms_fine, twogrid_norms_fine, dfly_norms_fine)

better_norms_fine = refine_norms_of_powers(norms_fine, m_extend)

w_fine = invariant_vector(B_fine, Q_fine)
@show distance_from_invariant(B_fine, D, Q_fine, w_fine, better_norms_fine)

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
(::Type{CA})(::UndefInitializer) where {CA<:Caggregate} = init(bitstype(CA), Val(sizeof(CA)))
@generated function init(::Type{CA}, ::Val{size}) where {CA<:Caggregate, size}
	init = ((zero(UInt8) for _ in 1:size)...,)
	return :(CA(undef, $(init)))
end

function (::Type{CA})(init::Union{NTuple{N, UInt8} where {N}, CA, Cbinding{CA}, Cptrs{CA}, UndefInitializer} = undef; kwargs...) where {CA<:Caggregate}
	result =
		init isa UndefInitializer ? bitstype(CA)(undef) :
		init isa Tuple ? bitstype(CA)(undef, init) :
		init isa Cbinding ? init[] :
		init isa Cptrs ? init[] :
		init

	for kwarg in kwargs
		result = initproperty(result, kwarg...)
	end

	return result
end

Base.read(io::IO, ::Type{CA}) where {CA<:Caggregate} = CA(((read(io, UInt8) for _ in 1:sizeof(CA))...,))

Base.sizeof(::Type{CA}) where {CA<:Caggregate} = invoke(sizeof, Tuple{Any}, bitstype(CA))
Base.zero(::Type{CA}) where {CA<:Caggregate} = CA()
Base.convert(::Type{CA}, nt::NamedTuple) where {CA<:Caggregate} = CA(; nt...)

Base.fieldcount(::Type{CA}) where {CA<:Caggregates} = length(fieldnames(CA))
Base.fieldtypes(::Type{CA}) where {CA<:Caggregates} = ((Base.tuple_type_head(t) for (n, t) in fields(CA))...,)
Base.fieldnames(::Type{CA}) where {CA<:Caggregates} = ((n for (n, t) in fields(CA))...,)
Base.propertynames(::CA) where {CA<:Caggregates} = fieldnames(CA)

Base.getproperty(ca::CA, sym::Symbol) where {CA<:Caggregates} = field(CA, sym) === nothing ? error("type $(CA) has no field $(sym)") : load(ca, field(CA, sym))
Base.setproperty!(ca::CA, sym::Symbol, val) where {CA<:Caggregate} = error("cannot set value of field $(sym) in type $(CA)")
initproperty(ca::CA, sym::Symbol, val) where {CA<:Caggregate} = field(CA, sym) === nothing ? error("type $(CA) has no field $(sym)") : store!(ca, field(CA, sym), val)


text:
#=
Created on Wednesday 29 January 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Ticket to ride dataset

Choo choo
🚂
=#

module TicketToRide

tickettoride_edges = [(1, "Vancouver", "Seattle"),
          (1, "Seattle", "Portland"),
          (3, "Vancouver", "Calgary"),
          (6, "Calgary", "Winnipeg"),
          (6, "Winnipeg", "Sault Ste. Marie"),
          (4, "Winnipeg", "Helena"),
          (4, "Calgary", "Helena"),
          (6, "Seattle", "Helena"),
          (4, "Seattle", "Calgary"),
          (6, "Portland", "Salt Lake City"),
          (5, "Portland", "San Francisco"),
          (5, "San Francisco", "Salt Lake City"),
          (3, "San Francisco", "Los Angeles"),
          (2, "Los Angeles", "Las Vegas"),
          (3, "Los Angeles", "Phoenix"),
          (6, "Los Angeles", "El Paso"),
          (3, "Phoenix", "El Paso"),
          (3, "Phoenix", "Santa Fe"),
          (3, "Las Vegas", "Salt Lake City"),
          (5, "Phoenix", "Denver"),
          (3, "Salt Lake City", "Denver"),
          (3, "Helena", "Salt Lake City"),
          (6, "Helena", "Duluth"),
          (4, "Winnipeg", "Duluth"),
          (4, "Helena", "Denver"),
          (5, "Helena", "Omaha"),
          (4, "Denver", "Omaha"),
          (4, "Denver", "Kansas City"),
          (2, "Denver", "Santa Fe"),
          (2, "Santa Fe", "El Paso"),
          (3, "Santa Fe", "Oklahoma City"),
          (4, "Denver", "Oklahoma City"),
          (6, "El Paso", "Houston"),
          (4, "El Paso", "Dallas"),
          (5, "El Paso", "Oklahoma City"),
          (1, "Dallas", "Houston"),
          (2, "Dallas", "Oklahoma City"),
          (2, "Kansas City", "Oklahoma City"),
          (1, "Omaha", "Kansas City"),
          (2, "Omaha", "Duluth"),
          (3, "Duluth", "Chicago"),
          (4, "Omaha", "Chicago"),
          (6, "Duluth", "Toronto"),
          (3, "Duluth", "Sault Ste. Marie"),
          (5, "Sault Ste. Marie", "Montreal"),
          (2, "Montreal", "Boston"),
          (2, "Boston", "New York"),
          (3, "Montreal", "New York"),
          (3, "Montreal", "Toronto"),
          (4, "Toronto", "Chicago"),
          (3, "Chicago", "Pittsburgh"),
          (2, "Chicago", "Saint Louis"),
          (2, "Pittsburgh", "Toronto"),
          (2, "Toronto", "Sault Ste. Marie"),
          (2, "Pittsburgh", "New York"),
          (2, "Pittsburgh", "Washington"),
          (2, "Washington", "New York"),
          (2, "Washington", "Raleigh"),
          (2, "Pittsburgh", "Raleigh"),
          (5, "Pittsburgh", "Saint Louis"),
          (2, "Kansas City", "Saint Louis"),
          (2, "Nashville", "Saint Louis"),
          (2, "Little Rock", "Saint Louis"),
          (2, "Oklahoma City", "Little Rock"),
          (2, "Little Rock", "Dallas"),
          (3, "Little Rock", "Nashville"),
          (2, "Houston", "New Orleans"),
          (3, "Little Rock", "New Orleans"),
          (4, "New Orleans", "Atlanta"),
          (1, "Atlanta", "Nashville"),
          (4, "Nashville", "Pittsburgh"),
          (2, "Atlanta", "Raleigh"),
          (3, "Nashville", "Raleigh"),
          (2, "Raleigh", "Charleston"),
          (2, "Charleston", "Atlanta"),
          (6, "New Orleans", "Miami"),
          (5, "Atlanta", "Miami"),
          (4, "Charleston", "Miami")
          ]

cities_coordinates = Dict("Atlanta" => (-84.3901849, 33.7490987),
                               "Boston" => (-71.0595678, 42.3604823),
                               "Calgary" => (-114.0625892, 51.0534234),
                               "Charleston" => (-79.9402728, 32.7876012),
                               "Chicago" => (-87.6244212, 41.8755546),
                               "Dallas" => (-96.7968559, 32.7762719),
                               "Denver" => (-104.9847034, 39.7391536),
                               "Duluth" => (-92.1251218, 46.7729322),
                               "El Paso" => (-106.501349395577, 31.8111305),
                               "Helena" => (-112.036109, 46.592712),
                               "Houston" => (-95.3676974, 29.7589382),
                               "Kansas City" => (-94.5630298, 39.0844687),
                               "Las Vegas" => (-115.149225, 36.1662859),
                               "Little Rock" => (-92.2895948, 34.7464809),
                               "Los Angeles" => (-118.244476, 34.054935),
                               "Miami" => (-80.1936589, 25.7742658),
                               "Montreal" => (-73.6103642, 45.4972159),
                               "Nashville" => (-86.7743531, 36.1622296),
                               "New Orleans" => (-89.9750054503052, 30.03280175),
                               "New York" => (-73.9866136, 40.7306458),
                               "Oklahoma City" => (-97.5170536, 35.4729886),
                               "Omaha" => (-95.9378732, 41.2587317),
                               "Phoenix" => (-112.0773456, 33.4485866),
                               "Pittsburgh" => (-79.99, 40.42),
                               "Portland" => (-122.6741949, 45.5202471),
                               "Raleigh" => (-78.6390989, 35.7803977),
                               "Saint Louis" => (-90.12954315, 38.60187637),
                               "Salt Lake City" => (-111.8904308, 40.7670126),
                               "San Francisco" => (-122.49, 37.75),
                               "Santa Fe" => (-105.9377997, 35.6869996),
                               "Sault Ste. Marie" => (-84.320068, 46.52391),
                               "Seattle" => (-122.3300624, 47.6038321),
                               "Toronto" => (-79.387207, 43.653963),
                               "Vancouver" => (-123.1139529, 49.2608724),
                               "Washington" => (-77.0366456, 38.8949549),
                               "Winnipeg" => (-97.168579, 49.884017))

cities = [k for k in keys(cities_coordinates)]

export tickettoride_edges, cities_coordinates, cities
end  # module TicketToRide

max_stars_repo_name: mferrodo/STMO

max_issues_repo_name: mferrodo/STMO

max_forks_repo_name: mferrodo/STMO


text:
# Load Julia packages (libraries)

using StanModels, CSV

df = convert(DataFrame,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "rugged.csv"), delim=';'));

dcc = filter(row -> !(ismissing(row[:rgdppc_2000])), df)
dcc[!, :log_gdp] = log.(dcc[!, :rgdppc_2000])
dcc[!, :cont_africa] = Array{Float64}(convert(Array{Int}, dcc[!, :cont_africa]))
dcc[!, :rugged] = convert(Array{Float64}, dcc[!, :rugged])

# Define the Stan language model

m8_1s = "
data{
    int N;
    vector[N] log_gdp;
    vector[N] cont_africa;
    vector[N] rugged;
    vector[N] rugged_cont_africa;
}
parameters{
    real a;
    real bR;
    real bA;
    real bAR;
    real sigma;
}
model{
    vector[N] mu = a + bR * rugged + bA * cont_africa + bAR * rugged_cont_africa;
    sigma ~ uniform( 0 , 10 );
    bAR ~ normal( 0 , 10 );
    bA ~ normal( 0 , 10 );
    bR ~ normal( 0 , 10 );
    a ~ normal( 0 , 100 );
    log_gdp ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m8.1s", m8_1s);

# Input data for cmdstan

m8_1_data = Dict(
  "N" => size(dcc, 1),
  "log_gdp" => dcc[!, :log_gdp],
  "cont_africa" => dcc[!, :cont_africa],
  "rugged" => dcc[!, :rugged],
  "rugged_cont_africa" => dcc[!, :rugged] .* dcc[!, :cont_africa] );

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m8_1_data);

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.22 0.14  9.00  9.46   282    1
bR    -0.21 0.08 -0.33 -0.08   275    1
bA    -1.94 0.24 -2.33 -1.59   268    1
bAR    0.40 0.14  0.18  0.62   271    1
sigma  0.96 0.05  0.87  1.04   339    1
"

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
# see related https://discourse.julialang.org/t/nested-forwarddiff-jacobian-calls-with-inplace-function/21232/7
function make_out_of_place_func(f!, x, make_buffer)
    y = make_buffer(f!, x) # make_buffer overloads will be defined below
    x -> f!(y, x)
end

struct OutOfPlace{F,D}
    f!::F
    out_of_place_funcs::Dict{Type, D}
    make_buffer::D

    function OutOfPlace(f!::T, dict::Dict{Type, D}, shape::Tuple) where {T,D}

        mk(::T, p) = similar(p, shape)
        new{T,D}(f!, dict, mk)
    end
end

OutOfPlace(f!, shape) = OutOfPlace(f!, Dict{Type, Function}(), shape) #no sure how to remove the `Function` here

eval_f(f::F, x) where {F} = f(x) # function barrier
function (oop::OutOfPlace{F})(x) where {F}
    T = eltype(x)
    f = get!(() -> make_out_of_place_func(oop.f!, x, oop.make_buffer), oop.out_of_place_funcs, T)
    eval_f(f, x)
end

function make_hessian(f!,x0,p0)

    f = OutOfPlace(f!, (length(x0),))
    g! = (outjac, p) -> (ForwardDiff.jacobian!(outjac, f, p); outjac = outjac')
    g = OutOfPlace(g!,(length(x0),length(p0)))
    h! = (outjac2, p) -> (ForwardDiff.jacobian!(outjac2, g, p); reshape(outjac2',length(p0),length(p0),length(x0)))
    h!
end

struct Avv
	h!::Function
    hessians::Array{Float64}

    function Avv(h!::Function, n::Int, m::Int)
    	hessians = Array{Float64}(undef, m*n, n)
        new(h!, hessians)
    end
end

function (avv::Avv)(dir_deriv::AbstractVector, p::AbstractVector, v::AbstractVector)
    hess = avv.h!(avv.hessians, p) #half of the runtime
    vHv!(dir_deriv, hess,v) #half of the runtime, almost all the memory
end

function vHv!(dir_deriv::AbstractVector, hessians::AbstractArray,v::AbstractVector)
    tmp = similar(v) #v shouldn't be too large in general, so I kept it here
    vt = v'

    for i=1:length(dir_deriv)
        @views mul!(tmp, hessians[:,:,i], v) #this line is particularly expensive memory-wise
        dir_deriv[i] = vt*tmp
    end

end


text:
module Formats

using TranscodingStreams

export
		infer, specify, formatted, resource,
		isspecified, isinferred, isunknown, isambiguous, getformat, getcoding,
		openf, readf, readf!, writef

include("registry.jl")
include("inference.jl")
include("formatted.jl")
include("io.jl")

"""
Read, write and detect formatted data, based on MIME types

```julia
using Formats       # Framework for data formats and codings (this package)
using FormatCodecs  # Codecs for common codings (separate package)
```

To infer the format of a file from its filename:

```julia
f = infer("insulin.pdb")
```

The coding (such as data compression) is also automatically detected:

```julia
f = infer("myoglobin.gro.gz")
```

To infer a format based on a stream’s content:

```julia
f = infer(open("insulin.pdb"))
```

To specify a known file format and, optionally, a coding:

```julia
f1 = specify("insulin.dat", "structure/x-pdb")
f2 = specify("myoglobin.dat", "structure/x-gro", "application/gzip")
```

To check if a format was specified or inferred:

```julia
f = infer("kitten.png")
isspecified(f)  # -> False
isinferred(f)    # -> True
isunknown(f)    # -> False
isambiguous(f)  # -> False
```

To get the detected or specified format/coding:

```julia
h = infer("myoglobin.gro.gz")
getformat(h)  # -> "structure/x-gro"
getcoding(h)  # -> "application/gzip"
```

Functions `infer` and `specify` return a `Formatted` object which can be passed
to `read` and `write`:

```julia
mol = read(specify("insulin.pdb", "structure/x-pdb"))
write(infer("insulin.gro.gz"), mol)
```

Function `read!` targets a pre-allocated output object:

```julia
read!(infer("insulin.gro.gz"), mol))
```

`Formatted` objects created from a filename can be `open`ed, returning a new
`Formatted` object wrapping the underlying IO stream, and those wrapping IO
streams can be `close`d:

```julia
f1 = infer("insulin.pdb")
f2 = open(f1)
close(f2)
```

Functions `infer` and `specify` can be called with an existing `Formatted`
object. This will infer format and coding again, or override previous guesses
with the specified information:

```julia
f1 = infer("insulin.dat")
f2 = specify(f1, "structure/x-pdb")
f3 = infer(open(f2))
```

Convenience functions `openf`, `readf`, `readf!` and `writef` automatically
infer format and coding if passed a filename or IO stream, but preserve the
existing format/coding information if passed a `Formatted` object:

```julia
mol1 = readf("insulin.pdb")
writef("insulin.gro", mol1)
mol2 = read(openf("myoglobin.gro"))
writef(specify("myoglobin.dat", "structure/x-pdb"), mol2)
readf!("lysozyme.pdb", mol3)
```

Complete documentation is provided in `README.md`.
"""
Formats

end # module


text:
import gmsh

gmsh.initialize()

gmsh.model.add("t3")

lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(0, [1, 2], 1)
gmsh.model.addPhysicalGroup(1, [1, 2], 2)
gmsh.model.addPhysicalGroup(2, [1], 6)
gmsh.model.setPhysicalName(2, 6, "My surface")

h = 0.1
angle = 90.
ov = gmsh.model.geo.extrude([(2,1)], 0, 0, h, [8,2], [0.5,1])

ov = gmsh.model.geo.revolve([(2,28)], -0.1,0,0.1, 0,1,0, -pi/2, [7])

ov = gmsh.model.geo.twist([(2,50)], 0,0.15,0.25, -2*h,0,0, 1,0,0,
                          angle*pi/180., [10], [], true)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(3, [1, 2, ov[2][2]], 101)

gmsh.model.mesh.generate(3)
gmsh.write("t3.msh")

gmsh.finalize()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# TODO make it yield enodes only? ping pavel and marisa
# TODO STAGE IT! FASTER!
module EMatchCompiler

using AutoHashEquals
using ..Patterns

abstract type Instruction end
export Instruction

const Register = Int32

mutable struct Program
    instructions::Vector{Instruction}
    first_nonground::Int
    memsize::Int
    regs::Vector{Register}
    ground_terms::Dict{Pattern, Register}
end
export Program

function Program()
    Program(Instruction[], 0, 0, [], Dict{Pattern, Register}())
end

hasregister(prog::Program, i) = (prog.regs[i] != -1)
getregister(prog::Program, i) = prog.regs[i]
setregister(prog::Program, i, v) = (prog.regs[i] = v)
increment(prog::Program, i) = (prog.memsize += i)
memsize(prog::Program) = prog.memsize

Base.getindex(p::Program, i) = p.instructions[i]
Base.length(p::Program) = length(p.instructions)

@auto_hash_equals struct ENodePat
    head::Any
    # args::Vector{Register}
    args::UnitRange{Register}
end
export ENodePat

@auto_hash_equals struct Bind <: Instruction
    reg::Register
    enodepat::ENodePat
end
export Bind

@auto_hash_equals struct CheckClassEq <: Instruction
    left::Register
    right::Register
end
export CheckClassEq

@auto_hash_equals struct CheckType <: Instruction
    reg::Register
    type::Any
end
export CheckType

@auto_hash_equals struct Yield <: Instruction
    yields::Vector{Register}
end
export Yield

@auto_hash_equals struct Filter <: Instruction
    reg::Register
    head::Any
    arity::Int
end
export Filter

@auto_hash_equals struct Lookup <: Instruction
    reg::Register
    p::Pattern
end
export Lookup

@auto_hash_equals struct Fail <: Instruction
    err::Exception
end
export Lookup

# =============================================
# ========= GROUND patterns ================
# =============================================

function compile_ground!(reg, p::PatTerm, prog)
    if haskey(prog.ground_terms, p)
        # push!(prog.instructions, CheckClassEq(reg, prog.ground_terms[p]))
        return nothing
    end

    if isground(p)
        prog.ground_terms[p] = reg
        push!(prog.instructions, Lookup(reg, p))
        increment(prog, 1)
    else
        for p2 in p.args
            compile_ground!(prog.memsize, p2, prog)
        end
    end
end

function compile_ground!(reg, p::PatVar, prog)
    nothing
end

function compile_ground!(reg, p::PatTypeAssertion, prog)
    nothing
end

function compile_ground!(reg, p::PatLiteral, prog)
    if haskey(prog.ground_terms, p)
        return nothing
    end
    prog.ground_terms[p] = reg
    push!(prog.instructions, Lookup(reg, p))
    increment(prog, 1)
end

function compile_ground!(reg, p::PatEquiv, prog)
    compile_ground!(reg, p.left, prog)
    compile_ground!(reg, p.right, prog)
end

function compile_ground!(reg, p::PatSplatVar, prog)
    Fail(ErrorException("Splatting not yet supported"))
end

# =============================================
# ========= NONGROUND patterns ================
# =============================================

function compile_pat!(reg, p::PatTerm, prog)
    if haskey(prog.ground_terms, p)
        push!(prog.instructions, CheckClassEq(reg, prog.ground_terms[p]))
        return nothing
    end
    # a = [gensym() for i in 1:length(p.args)]
    c = memsize(prog)
    nargs = length(p.args)
    a = c:(c + nargs - 1)

    # println(p)
    # println(a)
    # filters = []
    # for (i, child) in enumerate(p.args)
    #     if child isa PatTerm
    #         push!(filters, Filter(a[i], child.head, arity(child)))
    #     end
    # end
    # println(filters)

    increment(prog, nargs)
    push!(prog.instructions, Bind(reg, ENodePat(p.head, a)))
    for (reg, p2) in zip(a, p.args)
        compile_pat!(reg, p2, prog)
    end
end

function compile_pat!(reg, p::PatVar, prog)
    if hasregister(prog, p.idx)
        push!(prog.instructions, CheckClassEq(reg, getregister(prog, p.idx)))
    else
        setregister(prog, p.idx, reg)
    end
end

function compile_pat!(reg, p::PatTypeAssertion, prog)
    if hasregister(prog, p.var.idx)
        push!(prog.instructions, CheckClassEq(reg, getregister(prog, p.var.idx)))
    else
        setregister(prog, p.var.idx, reg)
        push!(prog.instructions, CheckType(reg, p.type))
    end
end

function compile_pat!(reg, p::PatLiteral, prog)
    if haskey(prog.ground_terms, p)
        push!(prog.instructions, CheckClassEq(reg, prog.ground_terms[p]))
        return nothing
    end
    push!(prog.instructions, Bind(reg, ENodePat(p.val, 0:-1)))
end

function compile_pat!(reg, p::PatEquiv, prog)
    compile_pat!(reg, p.left, prog)
    compile_pat!(reg, p.right, prog)
end

function compile_pat!(reg, p::PatSplatVar, prog)
    Fail(ErrorException("Splatting not yet supported"))
end

#========================================================================================#

# EXPECTS INDEXES OF PATTERN VARIABLES TO BE ALREADY POPULATED
function compile_pat(p::Pattern)
    pvars = patvars(p)
    nvars = length(pvars)

    # The program will try to match against ground terms first
    prog = Program(Instruction[], 1, 1, fill(-1, nvars), Dict{Pattern, Register}())
    # println("compiling pattern $p")
    compile_ground!(1, p, prog)
    # println("compiled ground pattern \n $prog")
    prog.first_nonground = prog.memsize
    prog.memsize+=1

    # And then try to match against other patterns
    compile_pat!(prog.first_nonground, p, prog)
    push!(prog.instructions, Yield(prog.regs))
    # println("compiled pattern $p to \n $prog")
    return prog
end

export compile_pat

end


text:
@testset "solve_analytical.jl" begin

end


text:
using Test, WAVI, LinearAlgebra
@testset  "WAVI tests" begin
    @testset "Iceberg" begin
        include("verification_tests/iceberg_test.jl")
        sim=iceberg_test(end_time = 1000.)
        #Steady state iceberg thickness and velocity
        #a = 0.3 m/yr, A=2.0e-17 Pa^-3 a^-1
        #ice density 918 kg/m3 ocean density 1028.0 kg/m3, Glen law n=3.
        h0=((36.0*0.3/(2.0e-17))*(1.0/(9.81*918.0*(1-918.0/1028.0)))^3)^(1.0/4.0)
        u0=sim.model.grid.xxu*0.3/(2.0*h0)
        v0=sim.model.grid.yyv*0.3/(2.0*h0)
        relerr_h=norm(sim.model.fields.gh.h[sim.model.fields.gh.mask].-h0)/
                    norm(h0*ones(length(sim.model.fields.gh.h[sim.model.fields.gh.mask])))
        relerr_u=norm(sim.model.fields.gu.u[sim.model.fields.gu.mask]-u0[sim.model.fields.gu.mask])/
                         norm(u0[sim.model.fields.gu.mask])
        relerr_v=norm(sim.model.fields.gv.v[sim.model.fields.gv.mask]-v0[sim.model.fields.gv.mask])/
                         norm(v0[sim.model.fields.gv.mask])
        @test relerr_h < 1.0e-4
        @test relerr_u < 3.0e-4
        @test relerr_v < 3.0e-4
    end

end
if true
@testset "MISMIP+ verification experiments" begin
    @testset "MISMIP+ Ice0 verification experiments" begin
        @info "Performing MISMIP+ Ice0 verification experiments"
        include(joinpath("verification_tests","MISMIP_PLUS_Ice0.jl"))
        simulation=MISMIP_PLUS_Ice0()
        glx=WAVI.get_glx(simulation.model)
        glxtest=glx[[1,div(simulation.model.grid.ny,2),div(simulation.model.grid.ny,2)+1,simulation.model.grid.ny]]
        @test length(glx) == simulation.model.grid.ny #check that the grounding line covers the whole domain in the y-direction
        @test (glxtest[4]-glxtest[1])/(glxtest[4]+glxtest[1]) < 1e-4
        @test (glxtest[2]-glxtest[3])/(glxtest[2]+glxtest[3]) < 1e-4
        @test 480000<glxtest[1]<540000
        @test 480000<glxtest[4]<540000
        @test 430000<glxtest[2]<460000
        @test 430000<glxtest[3]<460000

        #check the melt rate for ice_1r is doing something sensible
        function m1(h, b)
            draft = -(918.0 / 1028.0) * h
            cavity_thickness = draft .- b
            cavity_thickness = max.(cavity_thickness, 0)
            m =  0.2*tanh.(cavity_thickness./75).*max.((-100 .- draft), 0)
            return m
        end
        melt = m1.(simulation.model.fields.gh.h, simulation.model.fields.gh.h)
        @test all(melt .>= 0)
        @test maximum(melt) < 100

    end
end
end


text:
# We do not inherit from Grid on purpose.
# This object cannot implement the Grid interface in a strict sense
struct DistributedGrid{Dc,Dp,A} <: GridapType
  grids::A
  function DistributedGrid(grids::AbstractPData{<:Grid{Dc,Dp}}) where {Dc,Dp}
    A = typeof(grids)
    new{Dc,Dp,A}(grids)
  end
end

local_views(g::DistributedGrid) = g.grids

function Geometry.OrientationStyle(
  ::Type{<:DistributedGrid{Dc,Dp,A}}) where {Dc,Dp,A}
  OrientationStyle(eltype(A))
end

function Geometry.RegularityStyle(
  ::Type{<:DistributedGrid{Dc,Dp,A}}) where {Dc,Dp,A}
  RegularityStyle(eltype(A))
end

Geometry.num_cell_dims(::DistributedGrid{Dc,Dp}) where {Dc,Dp} = Dc
Geometry.num_cell_dims(::Type{<:DistributedGrid{Dc,Dp}}) where {Dc,Dp} = Dc
Geometry.num_point_dims(::DistributedGrid{Dc,Dp}) where {Dc,Dp} = Dp
Geometry.num_point_dims(::Type{<:DistributedGrid{Dc,Dp}}) where {Dc,Dp} = Dp

# We do not inherit from GridTopology on purpose.
# This object cannot implement the GridTopology interface in a strict sense
struct DistributedGridTopology{Dc,Dp,A} <: GridapType
  topos::A
  function DistributedGridTopology(topos::AbstractPData{<:GridTopology{Dc,Dp}}) where {Dc,Dp}
    A = typeof(topos)
    new{Dc,Dp,A}(topos)
  end
end

local_views(a::DistributedGridTopology) = a.topo

function Geometry.OrientationStyle(
  ::Type{<:DistributedGridTopology{Dc,Dp,A}}) where {Dc,Dp,A}
  OrientationStyle(eltype(A))
end

function Geometry.RegularityStyle(
  ::Type{<:DistributedGridTopology{Dc,Dp,A}}) where {Dc,Dp,A}
  RegularityStyle(eltype(A))
end

Geometry.num_cell_dims(::DistributedGridTopology{Dc,Dp}) where {Dc,Dp} = Dc
Geometry.num_cell_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Dc
Geometry.num_point_dims(::DistributedGridTopology{Dc,Dp}) where {Dc,Dp} = Dp
Geometry.num_point_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Dp

struct DistributedFaceLabeling{A<:AbstractPData{<:FaceLabeling}}
  labels::A
end

local_views(a::DistributedFaceLabeling) = a.labels

function Geometry.add_tag_from_tags!(labels::DistributedFaceLabeling, name, tags)
  map_parts(labels.labels) do labels
    add_tag_from_tags!(labels, name, tags)
  end
end

# We do not inherit from DiscreteModel on purpose.
# This object cannot implement the DiscreteModel interface in a strict sense
struct DistributedDiscreteModel{Dc,Dp,A,B} <: GridapType
  models::A
  gids::B
  function DistributedDiscreteModel(
    models::AbstractPData{<:DiscreteModel{Dc,Dp}}, gids::PRange) where {Dc,Dp}
    A = typeof(models)
    B = typeof(gids)
    new{Dc,Dp,A,B}(models,gids)
  end
end

local_views(a::DistributedDiscreteModel) = a.models

function Geometry.num_cells(model::DistributedDiscreteModel)
 num_gids(model.gids)
end

function Geometry.get_grid(model::DistributedDiscreteModel)
  DistributedGrid(map_parts(get_grid,model.models))
end

function Geometry.get_grid_topology(model::DistributedDiscreteModel)
  DistributedGridTopology(map_parts(get_grid_topology,model.models))
end

function Geometry.get_face_labeling(model::DistributedDiscreteModel)
  DistributedFaceLabeling(map_parts(get_face_labeling,model.models))
end

# CartesianDiscreteModel

function Geometry.CartesianDiscreteModel(
  parts::AbstractPData{<:Integer},args...;kwargs...)

  desc = CartesianDescriptor(args...; kwargs...)
  nc = desc.partition
  msg = """
  A CartesianDiscreteModel needs a Cartesian subdomain partition
  of the right dimensions.
  """
  @assert length(size(parts)) == length(nc) msg
  gcids = PCartesianIndices(parts,nc,PArrays.with_ghost)
  gids = PRange(parts,nc,PArrays.with_ghost)
  models = map_parts(parts,gcids) do part, gcids
    cmin = first(gcids)
    cmax = last(gcids)
    CartesianDiscreteModel(desc,cmin,cmax)
  end
  DistributedDiscreteModel(models,gids)
end

## Helpers to partition a serial model
# Not very scalable but useful in moderate
# cell and proc counts

function compute_cell_graph(model::DiscreteModel,d::Integer=0)
  D = num_cell_dims(model)
  topo = get_grid_topology(model)
  cell_to_dfaces = get_faces(topo,D,d)
  dface_to_cells = get_faces(topo,d,D)
  _cell_graph(cell_to_dfaces,dface_to_cells)
end

function _cell_graph(cell_to_dfaces,dface_to_cells)
  # This can be improved using CSRR format
  ncells = length(cell_to_dfaces)
  c1 = array_cache(cell_to_dfaces)
  c2 = array_cache(dface_to_cells)
  ndata = 0
  for icell in 1:ncells
    ndata += 1
    dfaces = getindex!(c1,cell_to_dfaces,icell)
    for dface in dfaces
      jcells = getindex!(c2,dface_to_cells,dface)
      for jcell in jcells
        if jcell != icell
          ndata += 1
        end
      end
    end
  end
  I = zeros(Int32,ndata)
  J = zeros(Int32,ndata)
  p = 0
  for icell in 1:ncells
    p += 1
    I[p] = icell
    J[p] = icell
    dfaces = getindex!(c1,cell_to_dfaces,icell)
    for dface in dfaces
      jcells = getindex!(c2,dface_to_cells,dface)
      for jcell in jcells
        if jcell != icell
          p += 1
          I[p] = icell
          J[p] = jcell
        end
      end
    end
  end
  m = ncells
  n = ncells
  V = ones(Int8,ndata)
  g = sparse(I,J,V,m,n)
  fill!(g.nzval,Int8(1))
  g
end

function Geometry.DiscreteModel(
  parts::AbstractPData,
  model::DiscreteModel,
  cell_to_part::AbstractArray,
  cell_graph::SparseMatrixCSC = compute_cell_graph(model))

  ncells = num_cells(model)
  @assert length(cell_to_part) == ncells
  @assert size(cell_graph,1) == ncells
  @assert size(cell_graph,2) == ncells

  lcell_to_cell, lcell_to_part, gid_to_part = map_parts(parts) do part
    cell_to_mask = fill(false,ncells)
    icell_to_jcells_ptrs = cell_graph.colptr
    icell_to_jcells_data = cell_graph.rowval
    for icell in 1:ncells
      if cell_to_part[icell] == part
        cell_to_mask[icell] = true
        pini = icell_to_jcells_ptrs[icell]
        pend = icell_to_jcells_ptrs[icell+1]-1
        for p in pini:pend
          jcell = icell_to_jcells_data[p]
          cell_to_mask[jcell] = true
        end
      end
    end
    lcell_to_cell = findall(cell_to_mask)
    lcell_to_part = zeros(Int32,length(lcell_to_cell))
    lcell_to_part .= cell_to_part[lcell_to_cell]
    lcell_to_cell, lcell_to_part, cell_to_part
  end

  partition = map_parts(IndexSet,parts,lcell_to_cell,lcell_to_part)
  exchanger = Exchanger(partition;reuse_parts_rcv=true)
  gids = PRange(ncells,partition,exchanger,gid_to_part)

  models = map_parts(lcell_to_cell) do lcell_to_cell
    DiscreteModelPortion(model,lcell_to_cell)
  end

  DistributedDiscreteModel(models,gids)
end

# Triangulation

# We do not inherit from Triangulation on purpose.
# This object cannot implement the Triangulation interface in a strict sense
struct DistributedTriangulation{Dc,Dp,A,B} <: GridapType
  trians::A
  model::B
  function DistributedTriangulation(
    trians::AbstractPData{<:Triangulation{Dc,Dp}},
    model::DistributedDiscreteModel) where {Dc,Dp}
    A = typeof(trians)
    B = typeof(model)
    new{Dc,Dp,A,B}(trians,model)
  end
end

local_views(a::DistributedTriangulation) = a.trians

Geometry.num_cell_dims(::DistributedTriangulation{Dc,Dp}) where {Dc,Dp} = Dc
Geometry.num_cell_dims(::Type{<:DistributedTriangulation{Dc,Dp}}) where {Dc,Dp} = Dc
Geometry.num_point_dims(::DistributedTriangulation{Dc,Dp}) where {Dc,Dp} = Dp
Geometry.num_point_dims(::Type{<:DistributedTriangulation{Dc,Dp}}) where {Dc,Dp} = Dp

function Geometry.get_background_model(a::DistributedTriangulation)
  a.model
end

# Triangulation constructors

function Geometry.Triangulation(
  model::DistributedDiscreteModel;kwargs...)
  Triangulation(no_ghost,model;kwargs...)
end

function Geometry.BoundaryTriangulation(
  model::DistributedDiscreteModel;kwargs...)
  BoundaryTriangulation(no_ghost,model;kwargs...)
end

function Geometry.SkeletonTriangulation(
  model::DistributedDiscreteModel;kwargs...)
  SkeletonTriangulation(no_ghost,model;kwargs...)
end

function Geometry.Triangulation(
  portion,model::DistributedDiscreteModel;kwargs...)
  trians = map_parts(model.models,model.gids.partition) do model,gids
    Triangulation(portion,gids,model;kwargs...)
  end
  DistributedTriangulation(trians,model)
end

function Geometry.BoundaryTriangulation(
  portion,model::DistributedDiscreteModel;kwargs...)
  trians = map_parts(model.models,model.gids.partition) do model,gids
    BoundaryTriangulation(portion,gids,model;kwargs...)
  end
  DistributedTriangulation(trians,model)
end

function Geometry.SkeletonTriangulation(
  portion,model::DistributedDiscreteModel;kwargs...)
  trians = map_parts(model.models,model.gids.partition) do model,gids
    SkeletonTriangulation(portion,gids,model;kwargs...)
  end
  DistributedTriangulation(trians,model)
end

function Geometry.Triangulation(
  portion,gids::AbstractIndexSet,args...;kwargs...)
  trian = Triangulation(args...;kwargs...)
  filter_cells_when_needed(portion,gids,trian)
end

function Geometry.BoundaryTriangulation(
  portion,gids::AbstractIndexSet,args...;kwargs...)
  trian = BoundaryTriangulation(args...;kwargs...)
  filter_cells_when_needed(portion,gids,trian)
end

function Geometry.SkeletonTriangulation(
  portion,gids::AbstractIndexSet,args...;kwargs...)
  trian = SkeletonTriangulation(args...;kwargs...)
  filter_cells_when_needed(portion,gids,trian)
end

function Geometry.InterfaceTriangulation(
  portion,gids::AbstractIndexSet,args...;kwargs...)
  trian = InterfaceTriangulation(args...;kwargs...)
  filter_cells_when_needed(portion,gids,trian)
end

function filter_cells_when_needed(
  portion::PArrays.WithGhost,
  cell_gids::AbstractIndexSet,
  trian::Triangulation)

  trian
end

function filter_cells_when_needed(
  portion::PArrays.NoGhost,
  cell_gids::AbstractIndexSet,
  trian::Triangulation)

  remove_ghost_cells(trian,cell_gids)
end

function filter_cells_when_needed(
  portion::FullyAssembledRows,
  cell_gids::AbstractIndexSet,
  trian::Triangulation)

  trian
end

function filter_cells_when_needed(
  portion::SubAssembledRows,
  cell_gids::AbstractIndexSet,
  trian::Triangulation)

  remove_ghost_cells(trian,cell_gids)
end

# For the moment remove_ghost_cells
# refers to the triangulation faces
# pointing into the ghost cells
# in the underlying background discrete
# model. This might change when solving
# multi-field PDEs with one of the fields
# defined on the boundary (e.g. a Lagrange multiplier)
function remove_ghost_cells(trian,gids)
  model = get_background_model(trian)
  D = num_cell_dims(model)
  glue = get_glue(trian,Val(D))
  remove_ghost_cells(glue,trian,gids)
end

function remove_ghost_cells(glue::FaceToFaceGlue,trian,gids)
  tcell_to_mcell = glue.tface_to_mface
  mcell_to_part = gids.lid_to_part
  tcell_to_part = view(mcell_to_part,tcell_to_mcell)
  tcell_to_mask = tcell_to_part .== gids.part
  view(trian, findall(tcell_to_mask))
end

function remove_ghost_cells(glue::SkeletonPair,trian::SkeletonTriangulation,gids)
  ofacets = _find_owned_skeleton_facets(glue,gids)
  plus = view(trian.plus,ofacets)
  minus = view(trian.minus,ofacets)
  SkeletonTriangulation(plus,minus)
end

function remove_ghost_cells(glue::SkeletonPair,trian,gids)
  ofacets = _find_owned_skeleton_facets(glue,gids)
  view(trian,ofacets)
end

function _find_owned_skeleton_facets(glue,gids)
  glue_p = glue.plus
  glue_m = glue.minus
  T = eltype(gids.lid_to_part)
  ntfaces = length(glue_p.tface_to_mface)
  tface_to_part = zeros(T,ntfaces)
  for tface in 1:ntfaces
    mface_p = glue_p.tface_to_mface[tface]
    mface_m = glue_m.tface_to_mface[tface]
    gcell_p = gids.lid_to_gid[mface_p]
    gcell_m = gids.lid_to_gid[mface_m]
    if gcell_p > gcell_m
      part = gids.lid_to_part[mface_p]
    else
      part = gids.lid_to_part[mface_m]
    end
    tface_to_part[tface] = part
  end
  findall(part->part==gids.part,tface_to_part)
end


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#

using HOODESolver
using LinearAlgebra
using Plots
using Plots.PlotMeasures

function henon_heiles(u, p, t)
    return [0, u[4], -2u[1] * u[2], -u[2] - u[1]^2 + u[2]^2]
end

function fctmain(n_tau, prec)
    setprecision(prec)
    u0=BigFloat.([90, -44, 83, 13]//100)
    tab_eps = zeros(BigFloat,8)
    epsilon=big"0.19"
    for i=1:size(tab_eps,1)
        tab_eps[i] = epsilon
        epsilon /= 10
    end
    nbmaxtest=12
    order=6
    t_max = big"1.0"
    y = ones(Float64, nbmaxtest, size(tab_eps,1) )
    x=zeros(Float64,nbmaxtest)
    ind=1
    A=[0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    for epsilon in tab_eps
        prob = HOODEProblem(henon_heiles, u0, (big"0.0",t_max), missing, A, epsilon)
        tabsol = Array{Array{BigFloat,1},1}(undef,0)
        eps_v = convert(Float32,epsilon)
        nb = 60*2^nbmaxtest
        res = solve(prob, nb_tau=n_tau, order=order, nb_t=nb, dense=false)
        par_u0=res.par_u0
        solRef = res[end]
        nb = 100
        indc =1
        labels=Array{String,2}(undef, 1, ind)
        while indc <= nbmaxtest
            res = solve(prob, nb_tau=n_tau, order=order, nb_t=nb,par_u0=par_u0, dense=false)
            sol = res[end]
            nm = norm(sol - solRef, Inf)
            y[indc,ind] = (nm < 1) ? nm : NaN
            x[indc] = 1.0/nb
            nb *= 2
            indc += 1
        end
        for i=1:ind
            labels[1,i] = "  ε = $(convert(Float32,tab_eps[i])) "
        end
        p=Plots.plot(
    x,
    view(y,:,1:ind),
    xlabel="Δt",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:bottomright,
    label=labels,
    marker=2
)
        pp=Plots.plot(
    x,
    view(y,:,1:ind),
    xlabel="Δt",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:bottomright,
    label=labels,
    marker=2,
    bottom_margin=30px,
    left_margin=60px
)
        prec_v = precision(BigFloat)
        Plots.savefig(p,"out/r2p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.pdf")
        Plots.savefig(p,"out/r2p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.png")
        Plots.savefig(pp,"out/r2pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.pdf")
        Plots.savefig(pp,"out/r2pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.png")
        ind+= 1
    end
end

fctmain(32,256)


text:
export CondPlanar

"""
Implementation of Planar (Conditional Version)
"""
struct CondPlanar{T <: AbstractFloat} <: AbstractCondICNF{T}
    re::Optimisers.Restructure
    p::AbstractVector{T}

    nvars::Integer
    basedist::Distribution
    tspan::Tuple{T, T}

    solver_test::SciMLBase.AbstractODEAlgorithm
    solver_train::SciMLBase.AbstractODEAlgorithm

    sensealg_test::SciMLBase.AbstractSensitivityAlgorithm
    sensealg_train::SciMLBase.AbstractSensitivityAlgorithm

    acceleration::AbstractResource
    array_mover::Function

    # trace_test
    # trace_train
end

function CondPlanar{T}(
        nn::PlanarNN,
        nvars::Integer,
        ;
        basedist::Distribution=MvNormal(zeros(T, nvars), Diagonal(ones(T, nvars))),
        tspan::Tuple{T, T}=convert(Tuple{T, T}, (0, 1)),

        solver_test::SciMLBase.AbstractODEAlgorithm=default_solver_test,
        solver_train::SciMLBase.AbstractODEAlgorithm=default_solver_train,

        sensealg_test::SciMLBase.AbstractSensitivityAlgorithm=default_sensealg,
        sensealg_train::SciMLBase.AbstractSensitivityAlgorithm=default_sensealg,

        acceleration::AbstractResource=default_acceleration,
        ) where {T <: AbstractFloat}
    array_mover = make_mover(acceleration, T)
    p, re = destructure(nn)
    CondPlanar{T}(
        re, p |> array_mover, nvars, basedist, tspan,
        solver_test, solver_train,
        sensealg_test, sensealg_train,
        acceleration, array_mover,
    )
end

function augmented_f(icnf::CondPlanar{T}, ys::Union{AbstractMatrix{T}, CuArray{T, 2}}, sz::Tuple{T2, T2})::Function where {T <: AbstractFloat, T2 <: Integer}
    o_ = ones(T, sz) |> icnf.array_mover

    function f_aug(u, p, t)
        m_ = icnf.re(p)
        m = Chain(
            x -> vcat(x, ys),
            m_,
        )
        z = u[1:end - 1, :]
        mz, back = Zygote.pullback(m, z)
        J = only(back(o_))
        trace_J = transpose(m_.u) * J
        vcat(mz, -trace_J)
    end
    f_aug
end

function inference(icnf::CondPlanar{T}, mode::TestMode, xs::AbstractMatrix{T}, ys::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat}
    xs = xs |> icnf.array_mover
    ys = ys |> icnf.array_mover
    zrs = zeros(T, 1, size(xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, ys, size(xs))
    prob = ODEProblem{false}(f_aug, vcat(xs, zrs), icnf.tspan, p)
    sol = solve(prob, icnf.solver_test; sensealg=icnf.sensealg_test)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    Δlogp = fsol[end, :]
    logp̂x = logpdf(icnf.basedist, z) - Δlogp
    logp̂x
end

function inference(icnf::CondPlanar{T}, mode::TrainMode, xs::AbstractMatrix{T}, ys::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat}
    xs = xs |> icnf.array_mover
    ys = ys |> icnf.array_mover
    zrs = zeros(T, 1, size(xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, ys, size(xs))
    prob = ODEProblem{false}(f_aug, vcat(xs, zrs), icnf.tspan, p)
    sol = solve(prob, icnf.solver_train; sensealg=icnf.sensealg_train)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    Δlogp = fsol[end, :]
    logp̂x = logpdf(icnf.basedist, z) - Δlogp
    logp̂x
end

function generate(icnf::CondPlanar{T}, mode::TestMode, ys::AbstractMatrix{T}, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat}
    ys = ys |> icnf.array_mover
    new_xs = isnothing(rng) ? rand(icnf.basedist, n) : rand(rng, icnf.basedist, n)
    new_xs = new_xs |> icnf.array_mover
    zrs = zeros(T, 1, size(new_xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, ys, size(new_xs))
    prob = ODEProblem{false}(f_aug, vcat(new_xs, zrs), reverse(icnf.tspan), p)
    sol = solve(prob, icnf.solver_test; sensealg=icnf.sensealg_test)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    z
end

function generate(icnf::CondPlanar{T}, mode::TrainMode, ys::AbstractMatrix{T}, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat}
    ys = ys |> icnf.array_mover
    new_xs = isnothing(rng) ? rand(icnf.basedist, n) : rand(rng, icnf.basedist, n)
    new_xs = new_xs |> icnf.array_mover
    zrs = zeros(T, 1, size(new_xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, ys, size(new_xs))
    prob = ODEProblem{false}(f_aug, vcat(new_xs, zrs), reverse(icnf.tspan), p)
    sol = solve(prob, icnf.solver_train; sensealg=icnf.sensealg_train)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    z
end

Flux.@functor CondPlanar (p,)

function loss(icnf::CondPlanar{T}, xs::AbstractMatrix{T}, ys::AbstractMatrix{T}, p::AbstractVector=icnf.p; agg::Function=mean) where {T <: AbstractFloat}
    logp̂x = inference(icnf, TrainMode(), xs, ys, p)
    agg(-logp̂x)
end


text:
using CFMG, Test
using Pkg

if (in("JLD", keys(Pkg.installed())) == false)
  Pkg.add("JLD")
end

if (in("Statistics", keys(Pkg.installed())) == false)
  Pkg.add("Statistics")
end

using JLD

function CFMG_test1()
  # try
    # generate model with CFMG
    generate_default_project("./default1")
    cp("config/Network.vff", "./default1/Network.vff", force=true)
    make_julia_model("./default1/Network.vff", "./default2")

    # run the model
    println("change pwd to make model run normal")
    println(pwd())
    cd("./default2/")
    println(pwd())
    println("run the model")
    print("__DIR__ at test1.jl ")
    println(@__DIR__)
    include("default2/Static.jl")  # this is interesting!

    # verify simulation results
    println("load std result")
    test1_result_path = normpath(joinpath(@__DIR__, "test1_result.jld"))
    test1_result = load(test1_result_path)["soln_std"]
    println(static_soln_object)
    @assert test1_result.exit_flag == static_soln_object.exit_flag "check exit_flag"
    @assert test1_result.status_flag == static_soln_object.status_flag "check status_flag"
    @assert norm(test1_result.objective_value - static_soln_object.objective_value) < 1e-8 "Check objective_value"
    @assert norm(test1_result.flux - static_soln_object.flux) < 1e-8 "Check flux"

    # remove generated files
    # sleep(60)
    cd("../")
    rm("./default1", recursive=true)
    rm("./default2", recursive=true)
    # @assert 0 == 1
    # println(pwd())
    # default1 = normpath(joinpath(@__DIR__, "default1/"))
    # default2 = normpath(joinpath(@__DIR__, "default2/"))
    # rm(default1, recursive=true)
    # rm(default2, recursive=true)

  # catch e
  #   return false
  # end

  return true
end

# CFMG_test1()

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
using Base.Test

using OPCSPs
using OPCSPs.MVNTools

positions = Vector{Float64}[[0,0], [1,0], [1,1], [0,1], [0.5,0.5]]
r = [0.0, 10.0, 100.0, 1.0, 20.0]

covariance = diagm(ones(length(r)))

rng = MersenneTwister(1)
d = randn(rng, length(r))

op = OPCSP(r, d, positions, covariance, 3.0, 1, 4)
path = solve_op(op)
@test path==[1,5,3,4]

path = solve_opcsp_feedback(op)
@test path==[1,5,3,4]

covariance[3,2] = covariance[2,3] = 0.8
rng = MersenneTwister(2)
d = rand!(rng, Array(Float64, length(r)), MVN(zeros(length(r)), covariance))
op = OPCSP(r, d, positions, covariance, 3.0, 1, 4)
path = solve_opcsp_feedback(op)
@test path==[1,5,3,4]


text:
export Normal

import Distributions

"""
    mutable struct Normal <: Dist{Float64}

`Normal` defines an *sfunc* representing unconditional Gaussian distributions.  Its
has no input, its output type is `Float64`, and its parameters are `(mean,standard deviation)`.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `bounded_probs`
- `compute_pi`
"""
mutable struct Normal <: Dist{Float64}
    params :: Tuple{Float64, Float64}
    """
        Normal(m::Float64, sd::Float64)

    `Normal`'s constructor

    # Arguments
    - `m::Float`: the mean of the `Normal`
    - `sd::Float64`: the standard deviation of the `Normal`
    """
    Normal(m::Float64, sd::Float64) = new((m,sd))
end

"returns the mean of the `Normal`"
mean(n::Normal) = n.params[1]
"returns the standard deviation of the `Normal`"
sd(n::Normal) = n.params[2]

"returns a `Distributions.Normal` from a `Scruff.SFuncs.Normal`"
dist(n) = Distributions.Normal(mean(n), sd(n))

@impl begin
    struct NormalSupport end
    function support(
            sf::Normal,
            parranges::NTuple{N,Vector},
            size::Integer,
            curr::Vector{Float64}) where N

        if isempty(curr)
            oldl = mean(sf)
            oldu = mean(sf)
        else
            oldl = minimum(filter(x -> x > -Inf, curr))
            oldu = maximum(filter(x -> x < Inf, curr))
        end

        if oldl == oldu
            if size <= 1
                result = [oldl]
                return result
            else
                l = oldl - sd(sf)
                u = oldl + sd(sf)
                result = Vector{Float64}()
                gap = (u-l) / (size-1)
                push!(result, l)
                for i = 1:(size-1)
                    push!(result, l + i*gap)
                end
                return result
            end
        else
            len = length(curr)
            if size <= len
                return curr
            else # we need to extend and thicken the current range while keeping existing points
                olddiff = oldu - oldl
                oldgap = olddiff / (len-1)
                newl = oldl - olddiff
                newu = oldu + olddiff
                result = copy(curr)
                i = newl
                while i < oldl
                    push!(result, i)
                    i += oldgap
                end
                i = newu
                while i > oldu
                    push!(result, i)
                    i -= oldgap
                end
                # Now we have 3 * len points covering twice the range
                # If size is greater, we will add points in between
                numperinterval = ceil(size / (3*len))
                if numperinterval > 1
                    newgap = oldgap / numperinterval
                    for x in copy(result)
                        for j = 1:(numperinterval-1)
                            push!(result, x + newgap * j)
                        end
                    end
                end
                sort!(result)
                return result
            end
        end
    end
end

@impl begin
    struct NormalSupportQuality end
    function support_quality(::Normal, parranges)
        :IncrementalSupport
    end
end

@impl begin
    struct NormalSample end
    function sample(sf::Normal, x::Tuple{})::Float64
        rand(dist(sf))
    end
end

@impl begin
    struct NormalLogcpdf end
    function logcpdf(sf::Normal, i::Tuple{}, o::Float64)::AbstractFloat
        Distributions.logpdf(dist(sf), o)
    end
end

@impl begin
    struct NormalBoundedProbsBoundedProbs
        numpartitions::Int64 = 10
    end

    function bounded_probs(
            sf::Normal,
            range::__OptVec{Float64},
            parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {N}

        intervals = make_intervals(range)
        lower = Array{Float64}(undef, length(range))
        upper = Array{Float64}(undef, length(range))
        d = dist(sf)
        ls = []
        us = []
        for interval in intervals
            (l, u) = interval
            if l == -Inf || u == Inf
                push!(ls, 0.0)
                push!(us, 1.0)
            else
                diff = (u - l) / numpartitions
                mn = 0.0
                mx = 0.0
                for i = 1:numpartitions
                    il = l + (i-1) * diff
                    iu = l + i * diff
                    p1 = Distributions.pdf(d, il) * diff
                    p2 = Distributions.pdf(d, iu) * diff
                    mn += min(p1,p2)
                    mx += max(p1,p2)
                end
                push!(ls, mn)
                push!(us, mx)
            end
        end
        # We get better bounds by considering the bounds on other intervals.
        # For any point, the lower bound can be normalized using the upper bounds of every other point, and vice versa.
        # This is especially important for intervals with -Inf or Inf as an
        # endpoint.
        for i = 1:length(intervals)
            otherls = 0.0
            otherus = 0.0
            for j = 1:i-1
                otherls += ls[j]
                otherus += us[j]
            end
            for j = i+1:length(intervals)
                otherls += ls[j]
                otherus += us[j]
            end
            l = max(ls[i], 1 - otherus)
            u = min(us[i], 1 - otherls)
            lower[i] = l
            upper[i] = u
        end
        return(lower, upper)
    end
end

# TODO: Replace this with a lazy implicit representation that doesn't require enumerating until the last minute
@impl begin
    struct NormalComputePi end

    function compute_pi(sf::Normal, range::__OptVec{Float64}, parranges::NTuple{N,Vector},
            incoming_pis::Tuple)::Dist{Float64} where {N}
        Cat(range, collect(map(x -> Distributions.pdf(dist(sf), x), range)))
    end

end
#=
@impl begin
    struct NormalExpectedStats end

    function expected_stats(sf::Normal, range::__OptVec{Float64}, parranges::NTuple{N,Vector},
            pis::NTuple{M,Dist},
            child_lambda::Score{Float64}) where {N,M}

        pis = [Distributions.pdf(dist(sf), x) for x in range]
        ls = [get_score(child_lambda, r) for r in range]
        prob = pis .* ls
        let totalX = 0.0, totalX2 = 0.0, count = 0.0
            for (i, x) in enumerate(range)
                totalX += x * prob[i]
                totalX2 += x^2 * prob[i]
                count += prob[i]
            end
            return (count,totalX,totalX2)
        end
    end
end

@impl begin
    struct NormalAccumulateStats end

    function accumulate_stats(sf::Normal, existing_stats, new_stats)
        existing_stats .+ new_stats
    end
end

@impl begin
    struct NormalInitialStats end

    function initial_stats(sf::Normal)
        return (0.0,0.0,0.0)
    end
end

@impl begin
    struct NormalMaximizeStats end

    function maximize_stats(sf::Normal, stats)
        (count, totalX, totalX2) = stats
        mean = totalX/count
        std = sqrt(totalX2/count - mean^2)
        return (mean, std)
    end
end
=#

Base.hash(n::Normal, h::UInt) = hash(n.params[2], hash(n.params[1], hash(:Normal, h)))
Base.isequal(a::Normal, b::Normal) = Base.isequal(hash(a), hash(b))
Base.:(==)(a::Normal, b::Normal) = Base.isequal(hash(a), hash(b))
Base.isless(a::Normal, b::Normal) = a.params[1] < b.params[1]
Base.:<(a::Normal, b::Normal) = Base.isless(a, b)


text:
@testset "detect_gradient_orientation" begin
    @info "Test: detect_gradient_orientation"

    @testset "Numerical" begin
        # Points on a compass in the sequence: E, NE, N, NW, W, SW, S, SE
        # expressed in raster coordinates (increasing rows, increasing columns)
        g₁ = [0, -1, -1, -1, 0, 1, 1, 1]
        g₂ = [1, 1, 0, -1, -1, -1, 0, 1]

        # Measuring from East in counter-clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = false,
                                                       compass_direction = 'E')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [0.0,  45.0,  90.0,  135.0,  180.0,  225.0,  270.0,  315.0])

        # Measuring from East in counter-clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = false,
                                                       compass_direction = 'E')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([0.0,  45.0,  90.0,  135.0,  180.0,  225.0,  270.0,  315.0]))

        # Measuring from South in counter-clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = false,
                                                       compass_direction = 'S')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [90.0,  135.0,  180.0,  225.0,  270.0,  315.0,  0.0,  45.0])

        # Measuring from South in counter-clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = false,
                                                       compass_direction = 'S')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([90.0,  135.0,  180.0,  225.0,  270.0,  315.0,  0.0,  45.0]))

        # Measuring from West in counter-clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = false,
                                                       compass_direction = 'W')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [180.0,  225.0,  270.0,  315.0,  360.0,  45.0,  90.0,  135.0])

        # Measuring from West in counter-clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = false,
                                                       compass_direction = 'W')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([180.0,  225.0,  270.0,  315.0,  360.0,  45.0,  90.0,  135.0]))

        # Measuring from N in counter-clockwise manner degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = false,
                                                       compass_direction = 'N')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [270.0,  315.0,  0.0,  45.0,  90.0,  135.0,  180.0,  225.0])

        # Measuring from N in counter-clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = false,
                                                       compass_direction = 'N')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([270.0,  315.0,  0.0,  45.0,  90.0,  135.0,  180.0,  225.0]))

        # Measuring from East in clockwise manner.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = true,
                                                       compass_direction = 'E')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [0.0,  315.0,  270.0,  225.0,  180.0,  135.0,  90.0,  45.0])

        # Measuring from East in clockwise manner.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = true,
                                                       compass_direction = 'E')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([0.0,  315.0,  270.0,  225.0,  180.0,  135.0,  90.0,  45.0]))

        # Measuring from South in clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = true,
                                                       compass_direction = 'S')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [270.0,  225.0,  180.0,  135.0,  90.0,  45.0,  0.0,  315.0])

        # Measuring from South in clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = true,
                                                       compass_direction = 'S')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([270.0,  225.0,  180.0,  135.0,  90.0,  45.0,  0.0,  315.0]))

        # Measuring from West in clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = true,
                                                       compass_direction = 'W')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [180.0,  135.0,  90.0,  45.0,  0.0,  315.0,  270.0,  225.0])

        # Measuring from West in clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = true,
                                                       compass_direction = 'W')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([180.0,  135.0,  90.0,  45.0,  0.0,  315.0,  270.0,  225.0]))

        # Measuring from North in clockwise manner in degrees.
        orientation_convention = OrientationConvention(in_radians = false,
                                                       is_clockwise = true,
                                                       compass_direction = 'N')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== [90.0,  45.0,  0.0,  315.0,  270.0,  225.0,  180.0,  135.0])

        # Measuring from North in clockwise manner in radians.
        orientation_convention = OrientationConvention(in_radians = true,
                                                       is_clockwise = true,
                                                       compass_direction = 'N')
        out = detect_gradient_orientation(g₁, g₂, orientation_convention)
        @test all(out .== deg2rad.([90.0,  45.0,  0.0,  315.0,  270.0,  225.0,  180.0,  135.0]))

        # Verify default setting.
        # Measuring from South in counter-clockwise manner in radians.
        out = detect_gradient_orientation(g₁, g₂)
        @test all(out .== deg2rad.([90.0,  135.0,  180.0,  225.0,  270.0,  315.0,  0.0,  45.0]))
    end

    @testset "Undefined Angles" begin
        # When the gradient magnitude is close to zero the gradient orientation
        # is undefined and we return a "sentinel" value of 360 degrees or
        # 2π.
        out₁ = detect_gradient_orientation([0.0], [0.0], OrientationConvention(in_radians = false))
        out₂ = detect_gradient_orientation([0.0], [0.0], OrientationConvention(in_radians = true))
        @test out₁ == [360]
        @test out₂ == [2π]

        detect_gradient_orientation!(out₂, [0.0], [0.0])
        @test out₂ == [2π]
    end

    @testset "Warning" begin
        convention = OrientationConvention(compass_direction = 'X')
        @test_logs (:warn, "Unrecognised compass_direction... using a default direction of south (S).")  detect_gradient_orientation([0.0], [0.0], convention)
    end

end


text:
@testset "Webscale" begin
    seed!(42)
    nInd = 20
    kernel = SqExponentialKernel()
    alg = Webscale(nInd)
    @test repr(alg) == "Webscale (m = $nInd)"
    test_Zalg(alg; kernel=kernel)
end


text:
using FileIO

include("tiffxml.jl")
include("Segmentation3D.jl")

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200316_result"
#cell_name = filter(x->occursin(r"^s",x), readdir(data_dir) )
cell_name = filter(x->occursin(r"^s[0-9][0-9]?-c[0-9][0-9]?.tiff$",x), readdir(data_dir) )

for i in cell_name
	print("processing $i")
	@time cell_raw = load(File(format"TIFF", "$data_dir/$i"))
	@time nucleus_all, threshold_all = extract3dnucleus(cell_raw)
    save(File(format"TIFF", "$data_dir/$i-clear.ome.tiff"), N0f16.(nucleus_all))
	save("$data_dir/$i.jld", "threshold", threshold_all)
	x, y, z_all = size(nucleus_all)
	@time embedxml(x, y, 20, z_all÷20, "$data_dir/$i-clear.ome.tiff")
end


text:
using Reexport

# reexport gives the Levenshtein function in scope of file with "include" call
@reexport using Levenshtein
using BioSequences
import BioSequences: reverse_complement
using Distributions
using StatsBase
using Distances
using MultivariateStats
#using NextGenSeqUtils
using PyPlot

include("phreds.jl")
include("utils.jl")
include("io.jl")
include("kmers.jl")
include("simulation.jl")
include("align.jl")
include("hmm.jl")
include("orient.jl")
include("paths.jl")
include("wrappers.jl")
include("demux.jl")
include("FAD.jl")
include("evodist.jl")