text:
using SpecialFunctions

type FastMslowDual
    # type to encapsulate the fast application of M = I + omega^2G*spadiagm(nu)
    GFFT :: Array{Complex128,2}
    nu :: Array{Float64,1}
    x  :: Array{Float64,1}
    y  :: Array{Float64,1}
    # number of points in the extended domain
    ne :: Int64
    me :: Int64
    # number of points in the original domain
    n  :: Int64
    m  :: Int64
    # frequency
    omega :: Float64
    # direction of the plane wave
    e :: Array{Float64,1}
    quadRule :: String
    function FastMslowDual(GFFT::Array{Complex128,2},
                           nu::Array{Float64,1},
                           x::Array{Float64,1},
                           y::Array{Float64,1},
                           ne::Int64,me::Int64,n::Int64,m::Int64,
                           k::Float64,e::Array{Float64,1};
                           quadRule::String = "trapezoidal")
      return new(GFFT,nu,x,y,ne,me,n, m, k,e, quadRule)
    end
end

# we need to overload these functions for gmres
import Base.*
import Base.A_mul_B!
import Base.size

function size(M::FastMslowDual, dim::Int64)
  # function to returns the size of the underliying matrix (M.m*M.n)^2
  return M.m*M.n
end

function A_mul_B!(Y,
                  M::FastMslowDual,
                  V)
    # in place matrix matrix multiplication
    @assert(size(Y) == size(V))
    # print(size(V))
    for ii = 1:size(V,2)
        Y[:,ii] = M*V[:,ii]
    end
end

function *(M::FastMslowDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # M using a Toeplitz reduction via a FFT
    # dummy function to call fastconvolution
    return fastconvolution(M,b)
end

@inline function fastconvolution(M::FastMslowDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # M using a Toeplitz reduction via a FFT

    # computing omega^2 nu G*(b)
    B = M.omega^2*(M.nu.*FFTconvolution(M,
                   exp(1im*M.omega*(M.e[1]*M.x + M.e[2]*M.y)).*b))

    # returning b + G*(b nu)
    return (-b + exp(-1im*M.omega*(M.e[1]*M.x + M.e[2]*M.y)).*B[:])
end

@inline function FFTconvolution(M::FastMslowDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # convolution of b times G

    if M.quadRule == "trapezoidal"

      #obtaining the middle index
      indMiddle = round(Integer, M.n)

      # Allocate the space for the extended B
      BExt = zeros(Complex128,M.ne, M.me);
      # Apply spadiagm(nu) and ented by zeros
      BExt[1:M.n,1:M.m]= reshape(M.nu.*b,M.n,M.m) ;

      # Fourier Transform
      BFft = fft(BExt)
      # Component-wise multiplication
      BFft = M.GFFT.*BFft
      # Inverse Fourier Transform
      BExt = ifft(BFft)

      # multiplication by omega^2
      B = (BExt[M.n:M.n+M.n-1, M.m:M.m+M.m-1]);

    elseif M.quadRule == "Greengard_Vico"
      # for this we use the Greengard Vico method in the
      # frequency domain

      # Allocate the space for the extended B
      BExt = zeros(Complex128,M.ne, M.me);
      # Apply spadiagm(nu) and ented by zeros
      BExt[1:M.n,1:M.m]= reshape(M.nu.*b,M.n,M.m) ;

      # Fourier Transform
      BFft = fftshift(fft(BExt))
     # Component-wise multiplication
      BFft = M.GFFT.*BFft
      # Inverse Fourier Transform
      BExt = ifft(ifftshift(BFft))

      # multiplication by omega^2
      B = BExt[1:M.n, 1:M.m];
    end
    return  B[:]
end

function buildFastConvolutionSlowDual(x::Array{Float64,1},y::Array{Float64,1},
                                      h::Float64,k::Float64,
                                      nu::Function,
                                      e::Array{Float64,1};
                                      quadRule::String = "trapezoidal")

  if quadRule == "trapezoidal"

    (ppw,D) = referenceValsTrapRule();
    D0      = D[round(Int,k*h)];
    (n,m) = length(x), length(y)
    Ge    = buildGConv(x,y,h,n,m,D0,k);
    GFFT  = fft(Ge);
    X = repmat(x, 1, m)[:]
    Y = repmat(y', n,1)[:]

    return FastMslowDual(GFFT,nu(X,Y),X,Y,2*n-1,2*m-1,n, m, k, e);

  elseif quadRule == "Greengard_Vico"

      Lp = 4*(abs(x[end] - x[1]) + h)
      L  =   (abs(x[end] - x[1]) + h)*1.5
      (n,m) = length(x), length(y)
      X = repmat(x, 1, m)[:]
      Y = repmat(y', n,1)[:]

      # this is depending if n is odd or not
      if mod(n,2) == 0
        kx = (-(2*n):1:(2*n-1));
        ky = (-(2*m):1:(2*m-1));

        KX = (2*pi/Lp)*repmat(kx, 1, 4*m);
        KY = (2*pi/Lp)*repmat(ky', 4*n,1);

        S = sqrt.(KX.^2 + KY.^2);

        GFFT = Gtruncated2D(L, k, S)
        return FastMslowDual(GFFT, nu(X,Y),X,Y, 4*n, 4*m,
                     n, m, k , e,quadRule="Greengard_Vico");
      else
        kx = (-2*n:1:2*n-1);
        ky = (-2*m:1:2*m-1);

        KX = (2*pi/Lp)*repmat( kx, 1,4*m);
        KY = (2*pi/Lp)*repmat(ky',4*n,  1);

        S = sqrt.(KX.^2 + KY.^2);

        GFFT = Gtruncated2D(L, k, S)

        return FastMslowDual(GFFT,nu(X,Y),X,Y,4*n,4*m,
                     n,m, k,e,quadRule = "Greengard_Vico");

    end
  end
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
# Ring Learning with Errors (RLWE) implementation in Julia
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# Reference:       https://eprint.iacr.org/2014/725.pdf
#                  https://eprint.iacr.org/2015/138.pdf
#                  https://eprint.iacr.org/2016/049.pdf
#                  https://securewww.esat.kuleuven.be/cosic/publications/article-2444.pdf
#                  https://www.microsoft.com/en-us/research/wp-content/uploads/2016/05/RLWE-1.pdf
# x86 C reference: https://github.com/ruandc/Ring-LWE-Encryption/blob/master/x86

using Krypto

abstract RLWEKey

type RLWEPubKey <: RLWEKey
    A::Array{Int16, 1}
    P::Array{Int16, 1}
    N::Integer
    Q::LatticeModuli
end

type RLWEPrivKey <: RLWEKey
    R::Array{Int16, 1}
    N::Integer
    Q::LatticeModuli
end

global const B = Int16(12)                # Global sampling bound B
global const Q = LatticeModuli(3, 12, 1)  # Global moduli (== 12289)
global const QBY2 = div(Int16(Q), 2)
global const QBY4 = div(Int16(Q), 4)

# Generates RLWE keypair
function RLWEKeyGen(N::Integer = 1024)
    A = genrandpoly(N, Q)
    R1, R2 = UniformSample(B, N), genrandpoly(N, Q)
    P = R1 - *(A, R2, Q)
    PUB = RLWEPubKey(NTTK(A, Q, N), NTTK(P, Q, N), N, Q)
    PRIV = RLWEPrivKey(NTTK(R2, Q, N), N, Q)
    return PUB, PRIV
end

# RLWE Encryption primitive
function RLWEEncrypt(K::RLWEPubKey, M::Array{UInt8, 1})
    ERR = [UniformSample(B, K.N) for i in 1:3]
    E1, E2 = NTTK(ERR[1], K.Q, K.N), NTTK(ERR[2], K.Q, K.N)
    M16 = [Int16(x * QBY2) % Int16(Q) for x in M]
    C1 = +(*(K.A, E1, K.Q), E2, K.Q)
    C2 = +(*(K.P, E1, K.Q), NTTK(+(ERR[3], M16), K.Q, K.N), K.Q)
    return C1, C2
end

# RLWE Decryption primitive
RLWEDecrypt(K::RLWEPrivKey, C1::Array{Int32, 1}, C2::Array{UInt8, 1}) = INTTK(+(*(C1, K.R, K.Q), C2, K.Q), K.Q, K.N)
RLWEDecrypt(K::RLWEPrivKey, C::Tuple{Array{UInt8, 1}, Array{UInt8, 1}}) = RLWEDecrypt(K, C[1], C[2])

function encrypt(::Type{Krypto.RLWE}, K::RLWEPubKey, M::Array{UInt8, 1})
    C1, C2 = RLWEEncrypt(K, M)
    # 0xff || n° of C1 || 0x00 || C1 || 0x00 × 8 || n° of C2 bytes || 0x00 || C2 || 0xff
    return vcat([0o377], [UInt8(degree(C1))], poly2bytes(C1), [0o0 for i in 1:8], [UInt8(degree(C2))], poly2bytes(C2), [0o377])
end

function decrypt(::Type{Krypto.RLWE}, K::RLWEPrivKey, M::Array{UInt8, 1})
    if M[1] != 0o377 || M[end] != 0o377 error("Decryption error.") end
    i, S1 =  2, 0
    while M[i] != 0o0 S1 += Int64(M[i]) * 256^(i - 1) end
    C1 = M[3:S1]
    if M[S1:(S1 + 8)] != [0o0 for i in 1:8] error("Decryption error.") end
    i, S2 =  S1 + 8, 0
    while M[i] != 0o0 S2 += Int64(M[i]) * 256^(i - S1 - 7) end
    C2 = M[(S1 + 8):(S1 + S2 + 8)]
    if M[end] != M[S1 + S2 + 8] || M[end] != 0o377 error("Decryption error.") end
    return poly2bytes(RLWEDecrypt(K, C1, C2))
end

ext: jl


text:
"""
Generic Boolean
"""
struct BooleanGeneric{T, S} <: AbstactBoolean
    m :: T
    dom :: BCube{Int}
    rng :: Val{S}
end

domain(f :: BooleanGeneric) = f.dom
codomain(f :: BooleanGeneric{T, S}) where {T, S} = S

(f :: BooleanGeneric)(v :: BoolVec) = f.m[value(v) + 1]

function BooleanGeneric(bc :: BCube)
    m = zeros(Int, length(bc))
    rng = Val(Int)
    return BooleanGeneric(m, bc, rng)
end

setindex!(f :: BooleanGeneric, v, k :: BoolVec) = (f.m[value(k) + 1] = v)

# given bool func, find its Walsh-Hadamard spectrum
function spectrum(f :: BoolFun)
    g = BooleanGeneric(domain(f))
    for y in domain(g)
        res = 0
        for x in domain(f)
            res += (-1)^(dot(x, y) + value(f(x)))
        end
        g[y] = res
    end
    return g
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# implementation of the GPUCompiler interfaces for generating GCN code

## target

export GCNCompilerTarget

Base.@kwdef struct GCNCompilerTarget <: AbstractCompilerTarget
    dev_isa::String
    features::String=""
end

llvm_triple(::GCNCompilerTarget) = "amdgcn-amd-amdhsa"

function llvm_machine(target::GCNCompilerTarget)
    triple = llvm_triple(target)
    t = Target(triple=triple)

    cpu = target.dev_isa
    feat = target.features
    optlevel = LLVM.API.LLVMCodeGenLevelDefault
    reloc = LLVM.API.LLVMRelocPIC
    tm = TargetMachine(t, triple, cpu, feat, optlevel, reloc)
    asm_verbosity!(tm, true)

    return tm
end

## job

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
runtime_slug(job::CompilerJob{GCNCompilerTarget}) = "gcn-$(job.target.dev_isa)$(job.target.features)"

const gcn_intrinsics = () # TODO: ("vprintf", "__assertfail", "malloc", "free")
isintrinsic(::CompilerJob{GCNCompilerTarget}, fn::String) = in(fn, gcn_intrinsics)

# we have to fake our target early in the pipeline because Julia's optimization passes
# weren't designed for a non-0 stack addrspace, and the AMDGPU target is very strict
# about which addrspaces are permitted for various code patterns
function process_module!(job::CompilerJob{GCNCompilerTarget}, mod::LLVM.Module)
    triple!(mod, llvm_triple(NativeCompilerTarget()))
    datalayout!(mod, julia_datalayout(NativeCompilerTarget()))
end

function process_entry!(job::CompilerJob{GCNCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    if job.source.kernel
        entry = lower_byval(job, mod, entry)

        # calling convention
        callconv!(entry, LLVM.API.LLVMAMDGPUKERNELCallConv)
    end

    entry
end

function add_lowering_passes!(job::CompilerJob{GCNCompilerTarget}, pm::LLVM.PassManager)
    add!(pm, ModulePass("LowerThrowExtra", lower_throw_extra!))
end
# We need to do alloca rewriting (from 0 to 5) after Julia's optimization
# passes because of two reasons:
# 1. Debug builds call the target verifier first, which would trip if AMDGPU
#    was the target at that time
# 2. We don't want any chance of messing with Julia's optimizations, since they
#    eliminate target-unsafe IR patterns
function optimize_module!(job::CompilerJob{GCNCompilerTarget}, mod::LLVM.Module)
    # revert back to the AMDGPU target
    triple!(mod, llvm_triple(job.target))
    datalayout!(mod, julia_datalayout(job.target))

    tm = llvm_machine(job.target)
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)

        add!(pm, FunctionPass("FixAllocaAddrspace", fix_alloca_addrspace!))

        run!(pm, mod)
    end
end

function lower_throw_extra!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw (extra)" begin

    throw_functions = [
        r"julia_bounds_error.*",
        r"julia_throw_boundserror.*",
        r"julia_error_if_canonical_getindex.*",
        r"julia_error_if_canonical_setindex.*",
        r"julia___subarray_throw_boundserror.*",
    ]

    for f in functions(mod)
        f_name = LLVM.name(f)
        for fn in throw_functions
            if occursin(fn, f_name)
                for use in uses(f)
                    call = user(use)::LLVM.CallInst

                    # replace the throw with a trap
                    let builder = Builder(ctx)
                        position!(builder, call)
                        emit_exception!(builder, f_name, call)
                        dispose(builder)
                    end

                    # remove the call
                    call_args = collect(operands(call))[1:end-1] # last arg is function itself
                    unsafe_delete!(LLVM.parent(call), call)

                    # HACK: kill the exceptions' unused arguments
                    for arg in call_args
                        # peek through casts
                        if isa(arg, LLVM.AddrSpaceCastInst)
                            cast = arg
                            arg = first(operands(cast))
                            isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                        end

                        if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                            unsafe_delete!(LLVM.parent(arg), arg)
                        end
                    end

                    changed = true
                end

                @compiler_assert isempty(uses(f)) job
            end
        end
    end

    end
    return changed
end
function fix_alloca_addrspace!(fn::LLVM.Function)
    changed = false
    alloca_as = 5
    ctx = context(fn)

    for bb in blocks(fn)
        for insn in instructions(bb)
            if isa(insn, LLVM.AllocaInst)
                ty = llvmtype(insn)
                ety = eltype(ty)
                addrspace(ty) == alloca_as && continue

                new_insn = nothing
                Builder(ctx) do builder
                    position!(builder, insn)
                    _alloca = alloca!(builder, ety, name(insn))
                    new_insn = addrspacecast!(builder, _alloca, ty)
                end
                replace_uses!(insn, new_insn)
                unsafe_delete!(LLVM.parent(insn), insn)
            end
        end
    end

    return changed
end

function emit_trap!(job::CompilerJob{GCNCompilerTarget}, builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    if Base.libllvm_version < v"9"
        rl_ft = LLVM.FunctionType(LLVM.Int32Type(ctx),
                                  [LLVM.Int32Type(ctx)])
        rl = if haskey(functions(mod), "llvm.amdgcn.readfirstlane")
            functions(mod)["llvm.amdgcn.readfirstlane"]
        else
            LLVM.Function(mod, "llvm.amdgcn.readfirstlane", rl_ft)
        end
        # FIXME: Early versions of the AMDGPU target fail to skip machine
        # blocks with certain side effects when EXEC==0, except when certain
        # criteria are met within said block. We emit a v_readfirstlane_b32
        # instruction here, as that is sufficient to trigger a skip. Without
        # this, the target will only attempt to do a "masked branch", which
        # only works on vector instructions (trap is a scalar instruction, and
        # therefore it is executed even when EXEC==0).
        rl_val = call!(builder, rl, [ConstantInt(Int32(32), ctx)])
        rl_bc = inttoptr!(builder, rl_val, LLVM.PointerType(LLVM.Int32Type(ctx)))
        store!(builder, rl_val, rl_bc)
    end
    call!(builder, trap)
end


text:
# lattice Boltzmann equation related stuff in one and two dimension
abstract type moments end
abstract type distribution end

mutable struct moments1d_64 <: moments
    height::Vector{Float64}
    velocity::Vector{Float64}
    energy::Vector{Float64}
end

mutable struct moments1d_32 <: moments
    height::Float32
    velocity::Float32
    energy::Float32
end

mutable struct distribution64_1d <: distribution
    f0::Vector{Float64}
    f1::Vector{Float64}
    f2::Vector{Float64}
end

mutable struct distribution32_1d <: distribution
    f0::Vector{Float32}
    f1::Vector{Float32}
    f2::Vector{Float32}
end

mutable struct distribution64_2d <: distribution
    f0::Vector{Float64}
    f1::Vector{Float64}
    f2::Vector{Float64}
    f3::Vector{Float64}
    f4::Vector{Float64}
    f5::Vector{Float64}
    f6::Vector{Float64}
    f7::Vector{Float64}
    f8::Vector{Float64}
end


text:
#EasyPlot plot builder interface (implemented by backends)
#-------------------------------------------------------------------------------

#=
Concerns with chosen EasyPlotBuilder implementation:
 - Garbage collection when you re-evaluate builder files?
   (evalfile() creates auto-named modules to house returned functions)

Alternative solutions considered:
 1. Implement functions for subtypes of EasyPlotBuilder (make EasyPlotBuilder abstract).
Concerns:
    - Too many `struct`s (basically one for each plot we define).
    - Objects might linger around more given loaded with evalfile() in load_plotbuilders()??
 2. Return references to modules created by load_plotbuilders() instead of EasyPlotBuilder objects
Concerns:
    - Users don't have a EasyPlotBuilder object to tell them how to structure their modules.
    - Garbage collection when you redefine module?
=#

#==Types
===============================================================================#
"""`struct EasyPlotBuilder`

Build `EasyPlot.Plot/PlotCollection` using data provided as `NamedTuple`.
"""
struct EasyPlotBuilder
	fnbuild::Function
end

#Dimension of entre plot canvas (distinct from PlotDim, but same values)
struct CanvasDim <: AbstractPlotDimensions
	w::Int; h::Int
end

mutable struct ShowOptions <: AbstractAttributeReceiver
	dim::AbstractPlotDimensions
end
_ShowOptions() = ShowOptions(plotautosize)
ShowOptions(args...; kwargs...) = _apply(_ShowOptions(), args, kwargs)

#Target application:
struct Target{T}; end
Target(id::Symbol) = Target{id}()

"""`abstract type AbstractBuilder`

Builds native plots (`NP`) from `EasyPlot.Plot/PlotCollection` objects.

Backend modules must implement a concrete builder object (`CB`), which generates
native plots.

# Must implement:
 - Member: `.postproc::Optional{Function}` with function signature `.postproc(plot::NP)`
 - `AbstractBuilder(::DS{:CB_ID}) = CB` (register symbol=>object mapping)
 - `getbuilder(::Target{T}, ::Type{CB}, args...; postproc=fn, kwargs...)` where T ∈ {:gui, :image}
   (args & kwargs are optional).
 - `build(b::CB, pcoll::PlotCollection)`
 - `_showable(mime::MIME, ::CB)`
 - `_show(io::IO, mime::MIME, opt::ShowOptions, plot::NP)`
"""
abstract type AbstractBuilder; end

#==Helper functions
===============================================================================#
function postproc(b::AbstractBuilder, plot)
	if !isnothing(b.postproc)
		b.postproc(plot)
	end
	return plot
end

#==EasyPlotBuilder interface:
===============================================================================#
function build(pf::EasyPlotBuilder, data::NamedTuple)
	plot = pf.fnbuild(data)
	@assert(isa(plot, PlotCollection), "EasyPlotBuilder.fnbuild() must return EasyPlot.PlotCollection")
	return plot
end

"""
`load_plotbuilders(@__DIR__, id1="gen1.jl", id2="gen2.jl", ...)`

Loads EasyPlotBuilder definitions from given files.

# Outputs
Returns a Dict of EasyPlotBuilder objects.

# NOTE
Must run from global scope to avoid world-age issue (code loading).
"""
function load_plotbuilders(rootpath::String; kwargs...)
#	outfile=splitext(filepath)[1] * ".svg"
	fmtlist = Dict{Symbol, EasyPlotBuilder}()
	for (k,filepath) in kwargs
		push!(fmtlist, k=>evalfile(joinpath(rootpath, filepath)))
	end
	return fmtlist
end

#==AbstractBuilder interface:
===============================================================================#
#(Includes catch-all functions/placeholders)

#Maps a symbol => concrete builder type:
AbstractBuilder(ds::DS{T}) where T = throw(MethodError(AbstractBuilder, (ds,)))
AbstractBuilder(bid::Symbol) = AbstractBuilder(DS(bid))

getbuilder(::Target{TS}, ::Type{T}, args...; kwargs...) where {TS, T} =
	throw(MethodError(getbuilder, (TS, T)))

"""`getbuilder(target::Symbol, builderid::Symbol, args...; kwargs...)`

Return a <: AbstractBuilder object specific to a `builderid` backend , yet
appropriate for use in a `target` application.

#Arugments
 - `target`: should be one of: `{:gui, :image}`.
 - `builderid`: should be one of: {:InspectDR, :Grace, :Grace_headless,
   :PyPlot, :PlotsJl}
"""
getbuilder(target::Symbol, builderid::Symbol, args...; kwargs...) where T =
	getbuilder(Target(target), AbstractBuilder(builderid), args...; kwargs...)

build(b::AbstractBuilder, pcoll::PlotCollection) = throw(MethodError(build, (b, plot)))

_showable(mime::MIME, b::AbstractBuilder) = throw(MethodError(Base.showable, (mime,b)))

#To be implemented by backend:
function _show(io::IO, mime::MIME, opt::ShowOptions, nativeplot::T) where T
	error("EasyPlot._show: no support for plots of type $T.")
end

#==Implement _apply() methods for ShowOptions
===============================================================================#
function _apply(opt::ShowOptions, ::DS{:plotdim}; w::Int=nothing, h::Int=nothing)
	if isnothing(w) || isnothing(h)
		throw(ArgumentError("set(`::ShowOptions`, plotdim=...): must specify w::Int & h::Int."))
	end
	opt.dim = PlotDim(w, h)
end

function _apply(opt::ShowOptions, ::DS{:dim}; w::Int=nothing, h::Int=nothing)
	if isnothing(w) || isnothing(h)
		throw(ArgumentError("set(`::ShowOptions`, dim=...): must specify w::Int & h::Int."))
	end
	opt.dim = CanvasDim(w, h)
end

#Last line


text:
import Base.show
using Random

mutable struct HelpModel <: CmdStanModels
    @shared_fields_stanmodels
    method::Help
end

function HelpModel(
  name::AbstractString,
  model::AbstractString,
  n_chains=[4];
  seed = RandomSeed(),
  init = Init(),
  output = Output(),
  tmpdir = mktempdir(),
  method = Help(),
  kwargs...)

  !isdir(tmpdir) && mkdir(tmpdir)

  StanBase.update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))
  sm = StanModel(joinpath(tmpdir, "$(name).stan"))

  output_base = StanRun.default_output_base(sm)
  exec_path = StanRun.ensure_executable(sm)

  stan_compile(sm)

  HelpModel(name, model, n_chains, seed, init, output,
    tmpdir, output_base, exec_path, String[], String[],
    Cmd[], String[], String[], String[], false, false, sm, method)
end

function help_model_show(io::IO, m, compact::Bool)
  println(io, "  name =                    \"$(m.name)\"")
  println(io, "  method =                  $(m.method)")
  println(io, "  n_chains =                $(get_n_chains(m))")
  println(io, "  output =                  Output()")
  println(io, "    file =                    \"$(split(m.output.file, "/")[end])\"")
  println(io, "    diagnostics_file =        \"$(split(m.output.diagnostic_file, "/")[end])\"")
  println(io, "    refresh =                 $(m.output.refresh)")
  println(io, "  tmpdir =                  \"$(m.tmpdir)\"")
end

show(io::IO, m::HelpModel) = help_model_show(io, m, false)


text:
@inline function _sampleParticlesSerial(rpio::RPIO{Particle}, μ!::F) where
  {Particle, F<:Function}
  particles::Vector{Particle} = rpio.particles
  N::Int64 = rpio.N
  rng::RNG = getRNG()
  for i in 1:N
    @inbounds μ!(particles[i], rng)
  end
end

@inline function _sampleParticlesParallel(rpio::RPIO{Particle}, μ!::F) where {Particle, F<:Function}
  particles::Vector{Particle} = rpio.particles
  N::Int64 = rpio.N
  nt::Int64 = rpio.threads
  b::Int64 = div(N, nt)
  Threads.@threads for i in 1:nt
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : N
    rng::RNG = getRNG()
    for j in start:finish
      @inbounds μ!(particles[j], rng)
    end
  end
end

@inline function _sampleParticles(model::Model, rpio::RPIO{Particle}) where Particle
  if rpio.threads > 1
    _sampleParticlesParallel(rpio, model.μ!)
  else
    _sampleParticlesSerial(rpio, model.μ!)
  end
end

@inline function _evaluatePotentialsSerial(rpio::RPIO{Particle}, p::Int64,
   G::F) where {Particle, F<:Function}
  particles::Vector{Particle} = rpio.particles
  potentials::Vector{Float64} = rpio.potentials
  v::Float64 = 0.0
  for i in 1:length(potentials)
    @inbounds potentials[i] = G(p, particles[i])
    @inbounds v += potentials[i]
  end
  return v
end

@inline function _evaluatePotentialsParallel(rpio::RPIO{Particle}, p::Int64,
   G::F) where {Particle, F<:Function}
  particles::Vector{Particle} = rpio.particles
  potentials::Vector{Float64} = rpio.potentials
  scratch::Vector{Float64} = rpio.scratch
  N::Int64 = rpio.N
  nt::Int64 = rpio.threads
  b::Int64 = div(N, nt)
  Threads.@threads for i in 1:nt
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : N
    v::Float64 = 0.0
    for j in start:finish
      @inbounds potentials[j] = G(p, particles[j])
      @inbounds v += potentials[j]
    end
    @inbounds scratch[i] = v
  end
  return sum(scratch)
end

 @inline function _evaluatePotentials(model::Model, rpio::RPIO{Particle},
   p::Int64) where Particle
   if rpio.threads > 1
     return _evaluatePotentialsParallel(rpio, p, model.G)
   else
     return _evaluatePotentialsSerial(rpio, p, model.G)
   end
 end

@inline function sampleOne(xs::Vector{Float64}, sumXs::Float64,
  start::Int64 = 1)
  rng = getRNG()
  u::Float64 = rand(rng)*sumXs
  i::Int64 = start
  v::Float64 = xs[i]
  while u > v
    i += 1
    @inbounds v += xs[i]
  end
  return i
end


text:
export i_sqdistance, i_dirtymul, i_factorial

"""
    i_sqdistance(dist!, x1, x2)

Squared distance between two points `x1` and `x2`.
"""
@i function i_sqdistance(dist!, x1::AbstractVector{T}, x2::AbstractVector) where T
    @inbounds for i=1:length(x1)
        x1[i] -= x2[i]
        dist! += x1[i] ^ 2
        x1[i] += x2[i]
    end
end

"""
    i_dirtymul(out!, x, anc!)

"dirty" reversible multiplication that computes `out! *= x` approximately for floating point numbers,
the `anc!` is anticipated as a number ~0.
"""
@i @inline function i_dirtymul(out!, x, anc!)
    anc! += out! * x
    out! -= anc! / x
    SWAP(out!, anc!)
end

@i @inline function i_dirtymul(out!::Int, x::Int, anc!::Int)
    anc! += out! * x
    out! -= anc! ÷ x
    SWAP(out!, anc!)
end

"""
    i_factorial(out!, n)

Compute the factorial `out! = factorial(n)`.
"""
@i function i_factorial(out!::Int, n::Int)
    INC(out!)
    @invcheckoff for i=1:n
        i_dirtymul(out!, i, 0)
    end
end


text:
using Statistics

Semantic = Union{<:Real,Array{<:Real,1}}
SemanticSet = Union{<:Real,Array{Semantic,1}}

@inline valueSemantic(semantic::Real, index::Int64) = semantic;
@inline valueSemantic(semantic::Array{<:Real,1}, index::Int64) = @inbounds semantic[index];

# Use "mutable" here so this struct is allocated in the heap, and the garbage collector can free the memory stored on a,b,c,d,S
mutable struct NodeEquation
    a::Semantic
    b::Semantic
    c::Semantic
    d::Semantic
    S::SemanticSet
    NodeEquation(a,b,c,d) = new(a,b,c,d,Semantic[])
    NodeEquation(a,b,c,d,S) = new(a,b,c,d,S)
end

#############################################################################################################
#
# Useful functions for evaluating the semantics
#

@inline allequal(x::Real) = true
@inline function allequal(x::Array{<:Real,1})
    length(x) < 2 && return true
    v = x[1]
    @inbounds for i=2:length(x)
        x[i] == v || return false
    end
    return true
end

@inline any0(x::Real) = (x==0)
@inline function any0(x::Array{<:Real,1})
    @inbounds for i=1:length(x)
        x[i] == 0 && return true
    end
    return false
end

# If vectors x and y have any zero in the same position
@inline any0(x::Real,y::Real) = (x==0) && (y==0)
@inline any0(x::Real,y::Array{<:Real,1}) = x==0 ? any0(y) : false
@inline any0(x::Array{<:Real,1},y::Real) = any0(y,x)
@inline function any0(x::Array{<:Real,1},y::Array{<:Real,1})
    @inbounds for i=1:length(x)
        (x[i]==0) && (y[i]==0) && return true
    end
    return false
end

@inline anyInf(x::Real) = isinf(x)
@inline function anyInf(x::Array{<:Real,1})
    @inbounds for i=1:length(x)
        isinf(x[i]) && return true
    end
    return false
end

@inline anyNaN(x::Real) = isnan(x)
@inline function anyNaN(x::Array{<:Real,1})
    @inbounds for i=1:length(x)
        isnan(x[i]) && return true
    end
    return false
end

@inline all0(x::Real) = (x==0)
@inline function all0(x::Array{<:Real,1})
    @inbounds for i=1:length(x)
        x[i]!=0 && return false
    end
    return true
end

@inline all1(x::Real) = (x==1)
@inline function all1(x::Array{<:Real,1})
    @inbounds for i=1:length(x)
        x[i]!=1 && return false
    end
    return true
end

@inline anyequal(x::Real, y::Real) = (x==y)
@inline function anyequal(x::Real, y::Array{<:Real,1})
    @inbounds for i=1:length(y)
        y[i]==x && return true
    end
    return false
end
@inline anyequal(x::Int64, y::Int64) = (x==y)
@inline function anyequal(x::Int64, y::Array{Int64,1})
    @inbounds for i=1:length(y)
        y[i]==x && return true
    end
    return false
end
@inline anyequal(x::Array{<:Real,1}, y::Real) = anyequal(y,x)
@inline function anyequal(x::Array{<:Real,1}, y::Array{<:Real,1})
    @inbounds for i=1:length(x)
        x[i]==y[i] && return true
    end
    return false
end

# Only for checking errors
equal(x,y; tolerance=NaN) = (x_y = abs.(x.-y); indNot0 = isa(x_y,Number) ? ((x_y!=0) ? 1 : []) : findall(x_y.!=0); if isnan(tolerance) tolerance = eltype(x)==Float32 ? 1e-3 : 1e-5; end; isempty(indNot0) ? true : all( ((x_y.<tolerance) .| (x_y.<abs.(x.*tolerance)) .| (x_y.<abs.(y.*tolerance)))[indNot0]) );

@inline cleanVector(v::Real) = v;
@inline cleanVector(v::Array{<:Real,1}) = (allequal(v)) ? v[1] : v

#############################################################################################################
#
# Functions for the operators
#

# These equations make the same computations, but allow to specify a buffer to write the results
# This avoids memory allocations and thus are faster
@inline addSemantics(x, y, buffer::Array{<:Real,1}) = (isa(x,Number) && isa(y,Number)) ? x+y : ( buffer .= x; buffer .+= y; cleanVector(buffer); )
@inline subSemantics(x, y, buffer::Array{<:Real,1}) = (isa(x,Number) && isa(y,Number)) ? x-y : ( buffer .= x; buffer .-= y; cleanVector(buffer); )
@inline mulSemantics(x, y, buffer::Array{<:Real,1}) = (isa(x,Number) && isa(y,Number)) ? x*y : ( buffer .= x; buffer .*= y; cleanVector(buffer); )
@inline divSemantics(x, y, buffer::Array{<:Real,1}) = (isa(x,Number) && isa(y,Number)) ? x/y : ( buffer .= x; buffer ./= y; cleanVector(buffer); )
@inline addSemantics(x, y, buffer::Nothing) = x.+y;
@inline subSemantics(x, y, buffer::Nothing) = x.-y;
@inline mulSemantics(x, y, buffer::Nothing) = x.*y;
@inline divSemantics(x, y, buffer::Nothing) = x./y;
@inline addSemantics(x, y, buffer::Real) = x.+y;
@inline subSemantics(x, y, buffer::Real) = x.-y;
@inline mulSemantics(x, y, buffer::Real) = x.*y;
@inline divSemantics(x, y, buffer::Real) = x./y;
@inline addSemantics(x, y) = x.+y;
@inline subSemantics(x, y) = x.-y;
@inline mulSemantics(x, y) = x.*y;
@inline divSemantics(x, y) = x./y;

#############################################################################################################
#
# Functions for calculating the equation coefficients a,b,c,d, and S
#

isNaN(S::SemanticSet) = isa(S,Number) && isnan(S)
semanticSetOpSemantic(S::SemanticSet, f::Function, semantic::Semantic) = (isempty(S) || isNaN(S)) ? S : Semantic[cleanVector(f.(s,semantic)) for s in S]
semanticOpSemanticSet(semantic::Semantic, f::Function, S::SemanticSet) = (isempty(S) || isNaN(S)) ? S : Semantic[cleanVector(f.(semantic,s)) for s in S]
function semanticSetDivSemantic(S::SemanticSet, semantic::Semantic)
    (isempty(S) || isNaN(S)) && return S;
    for s in S
        any0(s,semantic) && return NaN;
    end;
    return Semantic[cleanVector(s./semantic) for s in S];
end;
function semanticDivSemanticSet(semantic::Semantic, S::SemanticSet)
    isempty(S) && return Semantic[0.];
    isNaN(S) && return S;
    for s in S
        any0(s,semantic) && return NaN;
    end;
    newSemanticSet = Semantic[cleanVector(semantic./s) for s in S];
    for s in newSemanticSet
        s==0 && return newSemanticSet
    end;
    pushfirst!(newSemanticSet,0.);
    return newSemanticSet
end;

@inline equationAddChild1(y, a, b, c, d, S) = NodeEquation(                               a , all1(a) ? cleanVector(b.-y) : cleanVector(b.-(a.*y)) ,                               c , all0(c) ?  d : cleanVector(d.-(c.*y)), semanticSetOpSemantic(S,-,y) );
@inline equationAddChild2(x, a, b, c, d, S) = NodeEquation(                               a , all1(a) ? cleanVector(b.-x) : cleanVector(b.-(a.*x)) ,                               c , all0(c) ?  d : cleanVector(d.-(c.*x)), semanticSetOpSemantic(S,-,x) );
@inline equationSubChild1(y, a, b, c, d, S) = NodeEquation(                               a , all1(a) ? cleanVector(b.+y) : cleanVector(b.+(a.*y)) ,                               c , all0(c) ?  d : cleanVector(d.+(c.*y)), semanticSetOpSemantic(S,+,y) );
@inline equationSubChild2(x, a, b, c, d, S) = NodeEquation(                               a , all1(a) ? cleanVector(x.-b) : cleanVector((a.*x).-b) ,                               c , all0(c) ? -d : cleanVector((c.*x).-d), semanticOpSemanticSet(x,-,S) );
@inline equationMulChild1(y, a, b, c, d, S) = NodeEquation( all1(a) ? y : cleanVector(a.*y) ,                                                    b , all0(c) ? c : cleanVector(c.*y) ,                                     d, semanticSetDivSemantic(S,y) );
@inline equationMulChild2(x, a, b, c, d, S) = NodeEquation( all1(a) ? x : cleanVector(a.*x) ,                                                    b , all0(c) ? c : cleanVector(c.*x) ,                                     d, semanticSetDivSemantic(S,x) );
@inline equationDivChild1(y, a, b, c, d, S) = NodeEquation(                               a ,                                    cleanVector(b.*y) ,                               c , all0(d) ?  d :      cleanVector(d.*y), semanticSetOpSemantic(S,*,y) );
@inline equationDivChild2(x, a, b, c, d, S) = NodeEquation(                               b , all1(a) ?                 x :      cleanVector(a.*x) ,                               d , all0(c) ?  c :      cleanVector(c.*x), semanticDivSemanticSet(x,S) );

@inline equationAddChild1(y, eq::NodeEquation) = equationAddChild1(y, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationAddChild2(x, eq::NodeEquation) = equationAddChild2(x, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationSubChild1(y, eq::NodeEquation) = equationSubChild1(y, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationSubChild2(x, eq::NodeEquation) = equationSubChild2(x, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationMulChild1(y, eq::NodeEquation) = equationMulChild1(y, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationMulChild2(x, eq::NodeEquation) = equationMulChild2(x, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationDivChild1(y, eq::NodeEquation) = equationDivChild1(y, eq.a, eq.b, eq.c, eq.d, eq.S)
@inline equationDivChild2(x, eq::NodeEquation) = equationDivChild2(x, eq.a, eq.b, eq.c, eq.d, eq.S)

function semanticNotInDomain(semantic::Semantic, S::SemanticSet)
    isempty(S) && return false;
    isNaN(S) && return true;
    for s in S
        # anynan(s) && return false;
        anyequal(semantic,s) && return true;
    end;
    return false;
end;
semanticNotInDomain(semantic::Semantic, eq::NodeEquation) = semanticNotInDomain(semantic, eq.S)

function showEquation(eq::NodeEquation, semantic=nothing)
    println("a= ", eltype(eq.a), ".(", eq.a, ");");
    println("b= ", eltype(eq.b), ".(", eq.b, ");");
    println("c= ", eltype(eq.c), ".(", eq.c, ");");
    println("d= ", eltype(eq.d), ".(", eq.d, ");");
    # println("S= ", isempty(eq.S) ? "∅" : eq.S , ";");
    if isempty(eq.S)
        println("S= ", isempty(eq.S) ? "∅" : eq.S , ";");
    elseif isa(eq.S, Real) && isnan(eq.S)
        println("S= NaN;");
    else
        text = string(eq.S[1]);
        for i = 2:length(eq.S)
            text = string(text, ", ", string(eq.S[i]));
        end;

        println("S= {", text , "}");
    end;
    (semantic!=nothing) && println("semantic= ", semantic, ";");
end;

function checkEquation(eq::NodeEquation)
    @assert(isa(eq.a,Vector) || isa(eq.b,Vector));
    @assert(eltype(eq.a)==eltype(eq.b)==eltype(eq.c)==eltype(eq.d));
    # checkVector(v) = (isa(v,Number) || !allequal(v)) && !anyNaN(v) && !anyInf(v);
    checkVector(v) = !anyNaN(v) && !anyInf(v);
    @assert(checkVector(eq.a));
    @assert(checkVector(eq.b));
    @assert(checkVector(eq.c));
    @assert(checkVector(eq.d));
    if (isempty(eq.S))
        @assert(all0(eq.c))
    else
        if (isa(eq.S,Number))
            @assert(isnan(eq.S));
        else
            @assert(isa(eq.S,Array));
            @assert(sum(eq.S .== 0)<=1);
            for s in eq.S
                @assert(!anyNaN(s));
            end;
            poles = eq.S[end];
            if (isa(poles,Number) && isinf(poles))
                @assert(all0(eq.c));
            elseif (isa(poles,Number) && (poles==0))
                @assert(all0(eq.d));
            else
                if (anyInf(poles))
                    poles = poles[.!isinf.(poles)]
                end;
                if (!isempty(poles))
                    poles = cleanVector(poles);
                end;
                @assert(!anyNaN(poles));

                @assert(!any0(eq.c,eq.d));
                poles2 = eq.d./eq.c;

                if (isa(poles2,Vector))
                    if (anyInf(poles2))
                        poles2 = poles2[.!isinf.(poles2)]
                    end;
                    if (!isempty(poles2))
                        poles2 = cleanVector(poles2);
                    end;
                end;
                @assert(!anyNaN(poles2));
                @assert(equal(poles,poles2));
            end;
        end;
    end;
end;

function horizontalAsymptote(a,b,c,d,S; checkForErrors=false)
    if (isa(c,Number))
        checkForErrors && @assert(c!=0)
        return mean(a.^2)/(c^2);
    elseif (isa(a,Number))
        if (a==0)
            asymptote = b./d;
            if (isa(asymptote,Number))
                return (asymptote^2)*mean(c.==0)
            else
                asymptote[c.!=0] .= 0;
                return mean(asymptote.^2);
            end;
        else
            any0(c) && return NaN;
            return mean((a./c).^2);
        end;
    else
        # a and c are vectors
        asymptote = a./c;
        for i in 1:length(c)
            if (c[i]==0)
                (a[i]!=0) && return NaN
                asymptote[i] = (isa(b,Vector) ? b[i] : b) / (isa(d,Vector) ? d[i] : d);
            else
                if (a[i]==0)
                    asymptote[i] = 0;
                end;
            end;
        end;
    end;
    if (checkForErrors)
        @assert(!anyNaN(asymptote))
        @assert(!anyInf(asymptote))
    end;
    return mean(asymptote.^2);
end;
horizontalAsymptote(eq::NodeEquation; checkForErrors=false) = horizontalAsymptote(eq.a,eq.b,eq.c,eq.d,eq.S; checkForErrors=checkForErrors);

#############################################################################################################
#
# Functions for calculating the MSE from the equation values and the semantics of this node
#

calculateMSEFromEquation(semantic::Semantic, eq::NodeEquation; checkForErrors=false) = mean(( ((semantic.*eq.a).-eq.b)./((semantic.*eq.c).-eq.d) ).^2);

function calculateMSEFromEquationEfficient(semantic::Semantic, eq::NodeEquation; checkForErrors=false)
    a=eq.a; b=eq.b; c=eq.c; d=eq.d; S=eq.S;
    # Simplified common form of this equation:
    if all0(c)
        numerator = (semantic.*a).-b;
        if any0(d)
            any0(numerator,d) && return NaN;
            return Inf;
        end;
        if (isa(d,Number))
             # d is a constant
            if (checkForErrors) @assert(d!=0); end;
            result = mean(numerator.^2)/(d^2);
            # result = mean((mulFunction(semantic,a) .- b).^2) ./ (d^2)
        else
            # d is a vector
            result = mean((numerator./d).^2);
            # result = mean(((mulFunction(semantic,a).-b)./d).^2);
        end;
    elseif all1(c) && all0(d)
        if any0(semantic)
            any0(semantic,b) && return NaN;
            return Inf;
        end;
        # result = mean(((a.*semantic .- b)./semantic).^2);
        result = mean((a .- (b./semantic)).^2);
    elseif all0(d)
        if any0(semantic)
            any0(semantic,b) && return NaN;
            return Inf;
        end;
        any0(c) && return Inf;
        result = mean(((a.*semantic .- b)./(c.*semantic)).^2);
    else
        numerator = (semantic.*a).-b;
        denominator = (semantic.*c).-d;
        any0(numerator,denominator) && return NaN;
        any0(denominator) && return Inf;
        result = mean( (numerator./denominator).^2 );
    end;
    return result;
end;

function calculateMSEFromEquationMoreEfficient(semantic::Semantic, a,b,c,d,S; checkForErrors=false)
    # This function is similar to the one above, except for the fact that NaNs and Infs are returned as NaNs
    # Also, vectorized operations have been replaced by loops, so no memory allocations are needed
    #  Therefore, since less checks and memory allocations are done, it is slightly faster
    if all0(c)
        (any0(d)) && return NaN;
        if allequal(d)
            # d is a constant
            d = d[1];
            # return mean(((semantic.*a).-b).^2)/(d^2);

            n = max(length(a),length(b)); result = eltype(a)(0.);
            @inbounds for i in 1:n
                result += (valueSemantic(semantic,i)*valueSemantic(a,i) - valueSemantic(b,i))^2;
            end;
            result /= n*(d^2);

            if (checkForErrors) @assert(equal(result, mean(((semantic.*a).-b).^2)/(d^2))); end;
            return result;

            # result = mean((mulFunction(semantic,a) .- b).^2) ./ (d^2)
        else
            # d is a vector
            # return mean((((semantic.*a).-b)./d).^2);

            n = length(b); result = eltype(a)(0.);
            @inbounds for i in 1:n
                result += ((valueSemantic(semantic,i)*valueSemantic(a,i) - valueSemantic(b,i)) / valueSemantic(d,i)) ^2;
            end;
            result /= n;

            if (checkForErrors) @assert(equal(result, mean((((semantic.*a).-b)./d).^2))); end;
            return result;

            # result = mean(((mulFunction(semantic,a).-b)./d).^2);
        end;
    elseif all1(c) && all0(d)
        # any0(semantic) && return NaN;
        # result = mean(((a.*semantic .- b)./semantic).^2);
        # return mean((a .- (b./semantic)).^2);

        n = max(length(a),length(b)); result = eltype(a)(0.);
        @inbounds for i in 1:n
            semantic_i = valueSemantic(semantic,i);
            (semantic_i==0) && return NaN;
            result += (valueSemantic(a,i) - valueSemantic(b,i)/semantic_i) ^2;
        end;
        result /= n;

        if (checkForErrors) @assert(equal(result, mean((a .- (b./semantic)).^2))); end;
        return result;

    elseif all0(d)
        # any0(semantic) && return NaN;
        # any0(c) && return NaN;
        # return mean((((a.*semantic) .- b)./(c.*semantic)).^2);

        n = max(length(a),length(b)); result = eltype(a)(0.);
        @inbounds for i in 1:n
            semantic_i = valueSemantic(semantic,i);
            (semantic_i==0) && return NaN;
            c_i = valueSemantic(c,i);
            (c_i==0) && return NaN;
            result += ((valueSemantic(a,i)*semantic_i - valueSemantic(b,i))/(c_i*semantic_i) ) ^2;
        end;
        result /= n;

        if (checkForErrors) @assert(equal(result, mean((((a.*semantic) .- b)./(c.*semantic)).^2))); end;
        return result;
    else
        # denominator = (semantic.*c).-d;
        # any0(denominator) && return NaN;
        # return mean( (((semantic.*a).-b)./denominator).^2 );

        n = max(length(a),length(b),length(c),length(d)); result = eltype(a)(0.);
        @inbounds for i in 1:n
            denominator = (valueSemantic(c,i)*valueSemantic(semantic,i) - valueSemantic(d,i));
            (denominator==0) && return NaN;
            result += ((valueSemantic(a,i)*valueSemantic(semantic,i) - valueSemantic(b,i))/denominator)^2;
        end;
        result /= n;

        denominator = (semantic.*c).-d;
        any0(denominator) && return NaN;
        if (checkForErrors)
            denominator = (semantic.*c).-d;
            @assert(!any0(denominator))
            @assert(equal(result, mean( (((semantic.*a).-b)./denominator).^2 )));
        end;
        return result;
    end;
    return nothing;
end;
calculateMSEFromEquationMoreEfficient(semantic::Semantic, eq::NodeEquation; checkForErrors=false) = calculateMSEFromEquationMoreEfficient(semantic, eq.a,eq.b,eq.c,eq.d,eq.S; checkForErrors=checkForErrors);

#############################################################################################################
#
# Function for calculating the derivative of the equation on a value x
#

function derivativeEquation(eq::NodeEquation, x::Real)
    # return 2 .* mean( (b.*c .- a.*d) .* (a.*x .- b ) ./ ((c.*x .- d).^3) );
    a=eq.a; b=eq.b; c=eq.c; d=eq.d; S=eq.S;

    if anyDenominatorHas0Coefficients(c,d)
        return Inf;
    end;
    return 2 .* mean( (b.*c .- a.*d) .* (a.*x .- b ) ./ ((c.*x .- d).^3) );

    (a,b,c,d) = equation;
    bc_ad = (b.*c .- a.*d);
    result = bc_ad .* (a.*x .- b ) ./ ((c.*x .- d).^3);
    result[bc_ad.==0] .= 0.;
    return 2*mean(result);
end;

#############################################################################################################
#
# Function for calculating the reduction in MSE
#

function calculateMSEReduction(semantic::Semantic, equation::NodeEquation, mse::Real; checkForErrors=false)
    if (checkForErrors)
        @assert(!all0(equation.c) || !all0(equation.d));
        @assert(!anyNaN(semantic));
        @assert(!anyInf(semantic));
        @assert(!isnan(mse) && !isinf(mse));
    end;
    newMSE = calculateMSEFromEquationMoreEfficient(semantic, equation; checkForErrors=checkForErrors);
    if (checkForErrors)

        @assert(!isinf(newMSE));
        newMSE2 = calculateMSEFromEquationEfficient(semantic, equation; checkForErrors=checkForErrors);
        @assert(isnan(newMSE) == (isnan(newMSE2) || isinf(newMSE2)))
        if (!isnan(newMSE))
            @assert(equal(newMSE,newMSE2));
        end;
        newMSE3 = calculateMSEFromEquation(semantic, equation; checkForErrors=checkForErrors);
        @assert(isnan(newMSE2) == isnan(newMSE3));
        @assert(isinf(newMSE2) == isinf(newMSE3));
        if isinf(newMSE2)
            @assert(sign(newMSE2)==sign(newMSE3))
        elseif (!isnan(newMSE))
            @assert(equal(newMSE,newMSE2));
        end;
    end;
    return mse - newMSE;
end;

#############################################################################################################
#
# Functions for calculating the best constant for a node
#  Return the constant and the reduction in MSE
#

function calculateRoots(c,d; checkForErrors=false)
    (c==0) && return Array{eltype(c),1}();
    roots = d./c;
    if any0(c)
        # c must be a vector
        (checkForErrors) && @assert(isa(c,Vector));
        roots=roots[c.!=0];
    end;
    roots = unique(roots);
    if (checkForErrors)
        @assert(!any(isnan.(roots)));
        @assert(!any(isinf.(roots)));
        @assert(isa(roots,Vector));
    end;
    return roots;
end;

function calculateConstantGeneralCase(a,b,c,d,S; checkForErrors=false)

    if (checkForErrors)
        @assert(!isempty(S))
        @assert(!isNaN(S))
        @assert(!anyNaN(S[end]))
    end;

    zeros = calculateRoots(a,b; checkForErrors=checkForErrors);
    # These are not zeros of the whole equation, but zeros in each term

    isempty(zeros) && return NaN;
    zeros = zeros[.!semanticNotInDomain.(zeros, [S])];

    if (checkForErrors)
        @assert(!any(isnan.(zeros)));
    end;

    isempty(zeros) && return NaN;

    # Avoid those zeros which are too close to a root of the denominator
    # The roots of the denominator should already be calculated in the last term of S
    poles = S[end];
    if (isa(poles,Number) && isinf(poles))
        # In this case, any value is valid (there are no roots in the denominator)
        if (checkForErrors)
            @assert(all0(c));
            # @assert(isempty(cleanEquation(calculateRoots(c,d; checkForErrors=checkForErrors))));
            @assert(isempty(calculateRoots(c,d; checkForErrors=checkForErrors)));
        end;
    else

        # However, in order to avoid possible rounding errors, calculate the roots of the denominator instead of S[end]
        poles = calculateRoots(c,d; checkForErrors=checkForErrors);

        if (checkForErrors);
            @assert(!isempty(poles));
            @assert(!anyInf(poles));
            @assert(!anyNaN(poles));
        end;

        # poles = cleanVector(poles);

        function isNotClose(zero, poles, tolerance)
            for pole in poles
                limit = abs(pole*tolerance);
                (zero>=pole-limit) && (zero<=pole+limit) && return false;
            end;
            return true;
        end;
        zeros = zeros[isNotClose.(zeros,[poles],[eltype(a)(1e-2)])];

        isempty(zeros) && return NaN;

    end;

    yZeros = calculateMSEFromEquationMoreEfficient.(zeros, [a],[b],[c],[d],[S]; checkForErrors=checkForErrors);

    if (checkForErrors)
        @assert(!anyNaN(yZeros));
        @assert(!anyInf(yZeros));
    end;
    (_,index) = findmin(yZeros);
    zero = zeros[index];

    return zero;

end;
calculateConstantGeneralCase(eq::NodeEquation; checkForErrors=false) = ( a=eq.a; b=eq.b; c=eq.c; d=eq.d; S=eq.S; calculateConstantGeneralCase(a,b,c,d,S; checkForErrors=checkForErrors); )

function calculateConstant(a,b,c,d,S; checkForErrors=false)
    any0(c,d) && return NaN;
    lengthArray = maximum([length(a), length(b), length(c), length(d)]);
    # anyDenominatorHas0Coefficients(c,d) && return NaN;
    createVector(v::Vector) = v;
    createVector(x::Number) = ( v = Array{eltype(a),1}(undef,lengthArray); v.=x; return v; )
    a = createVector(a);
    b = createVector(b);
    c = createVector(c);
    d = createVector(d);

    if all0(c)
        any0(d) && return NaN;
        if (isa(a,Vector))
            index0 = (a.==0);
            num = (b.*a)./(d.^2);
            num[index0] .= eltype(a)(0.);
            den = ((a./d).^2);
            den[index0] .= eltype(a)(0.);
            constant = sum(num)./sum(den);
        else
            constant = sum((b.*a)./(d.^2))/sum((a./d).^2);
        end;
        semanticNotInDomain(constant, S) && return NaN;
    elseif all0(d)
        any0(c) && return NaN;
        if (isa(b,Vector))
            index0 = (b.==0);
            num = ((b./c).^2);
            num[index0] .= eltype(a)(0.);
            den = ((a.*b)./(c.^2));
            den[index0] .= eltype(a)(0.);
            constant = sum(num)/sum(den);
        else
            constant = sum((b./c).^2) ./ sum(a.*b./(c.^2));
        end;
        semanticNotInDomain(constant, S) && return NaN;
    elseif allequal(c) && allequal(d)
        mult = (b.*c .- a.*d);
        constant = sum(mult.*b)/sum(mult.*a);
        semanticNotInDomain(constant, S) && return NaN;
    else
        constant = calculateConstantGeneralCase(a,b,c,d,S; checkForErrors=checkForErrors);
    end;
    return constant;
end;
calculateConstant(eq::NodeEquation; checkForErrors=false) = ( a=eq.a; b=eq.b; c=eq.c; d=eq.d; S=eq.S; calculateConstant(a,b,c,d,S; checkForErrors=checkForErrors); )

function calculateConstantEfficient(a,b,c,d,S; checkForErrors=false)
    # Vectorial operations have been replaced by loops, so no memory allocation is needed and all the operations are much faster
    if all0(c)
        all0(a) && return NaN;
        any0(d) && return NaN;
        if allequal(a)
            # a is constant
            a = a[1];
            if (checkForErrors) @assert(isa(b,Vector)); end;
            # if isa(d,Number)
            if allequal(d)
                # d is constant
                d = d[1];
                # constant = mean(b)/a;
                constant = eltype(a)(0.);
                @inbounds for i in 1:length(b)
                    constant += b[i];
                end;
                constant /= (length(b)*a);

                if (checkForErrors) @assert(equal(constant,sum((b.*a)./(d.^2))/(length(b)*((a./d).^2)))); end;
            else
                # d is vector
                # inv_d = (1. ./ d).^2; constant = sum(b.*inv_d)/ (a*sum(inv_d));

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(d)
                    v = (1/d[i])^2;
                    numerator += valueSemantic(b,i)*v;
                    denominator += v;
                end;
                constant = numerator / (a*denominator);

                if (checkForErrors) @assert(equal(constant,sum((b.*a)./(d.^2))/sum((a./d).^2))); end;
            end;
        else
            # a is vector
            if allequal(d)
                # d is constant
                d = d[1];
                # constant = sum(b.*a)/sum(a.^2);

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(a)
                    a_i = valueSemantic(a,i);
                    numerator += a_i*valueSemantic(b,i);
                    denominator += (a_i^2);
                end;
                constant = numerator / denominator;

                if (checkForErrors) @assert(equal(constant,sum((b.*a)./(d.^2))/sum((a./d).^2))); end;
            else
                # d is vector
                # constant = sum((b.*a)./(d.^2))/sum((a./d).^2);

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(a)
                    a_i = valueSemantic(a,i);
                    d_i = valueSemantic(d,i);
                    numerator += a_i*valueSemantic(b,i)/(d_i^2);
                    denominator += (a_i/d_i)^2;
                end;
                constant = numerator / denominator;

                if (checkForErrors) @assert(equal(constant,sum((b.*a)./(d.^2))/sum((a./d).^2))); end;

            end;
        end;
        semanticNotInDomain(constant, S) && return NaN;

    elseif all0(d)
        any0(c) && return NaN;
        if allequal(b)
            # b is constant
            b = b[1];
            if (checkForErrors) @assert(isa(a,Vector)); end;
            # if isa(c,Number)
            if allequal(c)
                # c is constant
                c = c[1];
                # constant = b./mean(a);

                constant = eltype(a)(0.);
                @inbounds for i in 1:length(a)
                    constant += a[i];
                end;
                constant = length(a)*b/constant;

                if (checkForErrors) @assert(equal(constant,length(a)*((b./c).^2) ./ sum(a.*b./(c.^2)))); end;
            else
                # c is vector
                # inv_c = (1. ./ c).^2;
                # constant = b.*sum(inv_c)/sum(a.*inv_c);

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(c)
                    inv_c_i = (1/c[i])^2;
                    numerator += inv_c_i;
                    denominator += (valueSemantic(a,i)*inv_c_i);
                end;
                constant = b * numerator / denominator;

                if (checkForErrors) @assert(equal(constant,sum((b./c).^2) ./ sum(a.*b./(c.^2)))); end;
            end;
            semanticNotInDomain(constant, S) && return NaN;
        else
            # b is vector
            if allequal(c)
                # c is constant
                c = c[1];
                # constant = sum(b.^2)/sum(a.*b);

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(b)
                    b_i = b[i];
                    numerator += b_i^2;
                    denominator += (valueSemantic(a,i)*b_i);
                end;
                constant = numerator / denominator;

                if (checkForErrors) @assert(equal(constant,sum((b./c).^2) ./ sum(a.*b./(c.^2)))); end;
            else
                # c is vector
                # constant = sum((b./c).^2) ./ sum(a.*b./(c.^2));

                numerator = eltype(a)(0.); denominator = eltype(a)(0.);
                @inbounds for i in 1:length(c)
                    b_i = valueSemantic(b,i);
                    numerator += (b_i/c[i])^2;
                    denominator += (valueSemantic(a,i)*b_i/(c[i]^2));
                end;
                constant = numerator / denominator;

                if (checkForErrors) @assert(equal(constant,sum((b./c).^2) ./ sum(a.*b./(c.^2)))); end;
            end;
            semanticNotInDomain(constant, S) && return NaN;
        end;

    elseif allequal(c) && allequal(d)
        # constant = sum(b.*b.*c .- a.*b.*d)/sum(a.*b.*c .- a.*a.*d);

        # mult = (b.*c .- a.*d);
        # constant = sum(mult.*b)/sum(mult.*a);

        c = c[1]; d = d[1];
        numerator = eltype(a)(0.); denominator = eltype(a)(0.);
        @inbounds for i in 1:max(length(a),length(b))
            a_i = valueSemantic(a,i);
            b_i = valueSemantic(b,i);
            mult = b_i*c - a_i*d;
            numerator += b_i*mult;
            denominator += a_i*mult;
        end;
        constant = numerator / denominator;

        if (checkForErrors)
            mult = (b.*c .- a.*d);
            @assert(!all0(a));
            @assert(!all0(mult));
            @assert(equal(constant,sum(mult.*b)/sum(mult.*a)));
        end;

        semanticNotInDomain(constant, S) && return NaN;

    else
        any0(c,d) && return NaN;
        return calculateConstantGeneralCase(a,b,c,d,S; checkForErrors=checkForErrors);
    end;

    if (checkForErrors)
        @assert(!isnan(constant));
        @assert(!isinf(constant));
    end;
    return constant;
end;
calculateConstantEfficient(eq::NodeEquation; checkForErrors=false) = calculateConstantEfficient(eq.a,eq.b,eq.c,eq.d,eq.S; checkForErrors=checkForErrors);

function calculateConstantMinimizeEquation(equation::NodeEquation, mse::Real; checkForErrors=false)
    constant = calculateConstantEfficient(equation; checkForErrors=checkForErrors);
    if (checkForErrors)
        @assert(!isnan(mse) && !isinf(mse));
        @assert(!isinf(constant));
        newConstant = calculateConstant(equation; checkForErrors=checkForErrors);
        @assert(isnan(constant)==isnan(newConstant));
        if (!isnan(constant))
            @assert(equal(constant,newConstant));
        end;
        @assert(!isinf(constant))
    end;
    if (isnan(constant))
        reduction = -Inf;
    else
        reduction = calculateMSEReduction(constant, equation, mse; checkForErrors=checkForErrors);
    end;
    reduction = isnan(reduction) ? -Inf : reduction;
    return (constant, reduction);
end;


text:
push!(LOAD_PATH,".")

using DelimitedFiles
using CUDA
using Zygote
using Flux
using BSON
using Test
using Plots

using DrugRepositioning

data = readdlm(raw"E:\Projects\DrugRepositioning\data\origin\drug_disease_mat.csv", ',')
drugs = data[2:end,1] |> Array{String}
diseases = data[1,2:end] |> Array{String}
RD = data[2:end,2:end] |> Array{Bool}
data = readdlm(raw"E:\Projects\DrugRepositioning\data/origin/drug_pubchem_mat.csv", ',')
dp1 = data[2:end,2:end] |> Array{Float32}
data = readdlm(raw"E:\Projects\DrugRepositioning\data/origin/drug_target_go_mat.csv", ',')
dp2 = data[2:end,2:end] |> Array{Float32}
chara = Dict{Symbol,Array{Matrix}}(:drug => [dp1',dp2'], :disease => [])
data = readdlm(raw"E:\Projects\DrugRepositioning\data/origin/disease_similarity_mat.csv", ',')
sim = data[2:end,2:end] |> Array{Float32}

@testset "DrugRepositioning" begin
    local dataset
    local s
    @testset "dataset" begin
        sim=DrugRepositioning.get_similarity(RD)
        heatmap(sim)|>display
        @test true
        dataset = Dataset(drugs, diseases, RD, chara;α=0.3,step=4)
        @test true
        n = length(dataset.similaritys[:disease])
        add_similarity!(dataset, :disease, sim)
        @test length(dataset.similaritys[:disease]) - n == 1
        @test length(keys(dataset.similaritys))==2
        @test length(keys(dataset.similaritys[:drug]))==2
        @test length(keys(dataset.similaritys[:disease]))==1
        s = Samples(dataset, RD, RD, findall(RD .== 1);mode=:both)
        @test true
        s = Samples(dataset, RD, RD, findall(RD .== 1);use_gpu=true,mode=:both)
        @test true
        change_indices!(s,findall(RD.==0))
        @test Flux.Data._nobs(s)==sum(RD.==0)
    end
    @testset "experiment" begin
        labels = findall(RD .== 1)
        maxn = rand(3:7)
        n_folds = CrossValidation(maxn, labels, shuffle=true)
        @test true
        n = rand(1:maxn)
        train, test = n_folds(n)
        @test length(train) > length(test)
        @test length(train) + length(test) == length(labels)
        vr=ValidateResult(raw"E:\Projects\DrugRepositioning\result")
        @test true
        plot(vr)|>display
        @test true
        # println(vr)
        r=analyze(vr)
        @test true
        # println(r)
        ms=compared_scores(raw"E:\Projects\DrugRepositioning\data/all_result")
        pr_roc_curve(ms)|>display
        @test true
    end
    @testset "model" begin
        att = Attention(10)
        @test size(att(rand(Float32,10,10)))==(10,10)
        att2=Attention(10,1)
        @test size(att2(rand(Float32,10,10)))==(1,10)
        att3=AttrAttention(2888,500,3)|>gpu
        @test size(att3(CUDA.rand(Float32,2,1444,3,7,9)))==(500,7,9)
        att4=ScaleAttention(500)
        @test size(att4(rand(Float32,500,7,9)))==(500,7,9)
        x1 = rand(Float32, 2, 1444, 2, 16)
        m1 = RDConvolution(dataset,out=2)
        @test true
        ŷ1 = m1(x1)
        @test true
        @test size(ŷ1) == (2, 16)
        x2 = CUDA.rand(Float32, 2, 1444, 2, 16)
        m2 = RDConvolution(dataset;out=1024) |> Flux.gpu
        @test true
        ŷ2 = m2(x2)
        @test true
        @test size(ŷ2) == (1024, 16)
        DrugRepositioning.save("./weights.bson", m2)
        weights = DrugRepositioning.load("./weights.bson")
        m2 = RDConvolution(dataset) |> Flux.gpu
        Flux.loadparams!(m2, weights)
        ŷ3 = m2(x2)
        @test all(ŷ3 .== ŷ2)
        rm("./weights.bson")
        reset_model!(m2, init=Flux.kaiming_uniform)
        ŷ = m2(x2)
        @test true
        loader = Flux.Data.DataLoader(s, batchsize=3)
        for (i, (x1,x2), y) in loader
            ŷ = m2(x2)
            break
        end
        @test true
        m3=MSRD(dataset)|>gpu
        @test size(m3.attratt.encoders,3)==2
        for(i,(x1,x2),y) in loader
            ŷ=m3(x1)
            break
        end
        @test true
        m4=CombineModel(dataset)|>gpu
        for (i,(x1,x2),y) in loader
            ŷ=m4(x1,x2)
            @test size(ŷ)==(2,3)
            break
        end
        @test true
        # trainer = Trainer(m2, Flux.Losses.crossentropy, dataset, epochs=10)
        # training(trainer, Flux.Optimise.Descent())
        # @test true
    end
end

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
using FractionalDiffEq
using Test
using MittagLeffler

@testset "Test Diethelm PECE algorithms" begin
    fun(x, y) = 1-y
    prob=FDEProblem(fun, 1.8, 0, 5, 0.01)
    result=solve(prob, PECE())
    tspan=collect(0:0.01:5)

    target = []

    for i in 0:0.01:5
        push!(target, i^1.8*mittleff(1.8, 2.8,-i^1.8))
    end

    @test isapprox(result, target; atol=1)
end


text:
function roundtrip(x)
    io = IOBuffer()
    serialize(io, x)
    deserialize(seekstart(io)) == x
end

@testset "serialisation" begin
    @testset "basic types" for x in [
            true, 'x', 1, [1, 2],
        ]
        @test roundtrip(K(x))
    end
    @testset "tables" begin
        @test roundtrip(K_Table(a=1:5))
        @test roundtrip(K_KeyTable(1;a=[:a,:b], b=[1,2]))
    end
end


text:
module StorageServices

using ..Azure
using Azure.StorageManagementClient
using Azure.REST
using HTTP
using MbedTLS
using Dates
using Base64

include("common.jl")
include("blob.jl")
include("file.jl")
include("sas.jl")

end


text:
struct ShortTermIndices <: AbstractModelIndices
    hours::Vector{Int}
    plants::Vector{Plant}
    segments::Vector{Int}
end

struct ShortTermData{T <: AbstractFloat} <: AbstractModelData
    hydrodata::HydroPlantCollection{T,2}
    pricecurve::PriceCurve{T}

    function (::Type{ShortTermData})(plantdata::HydroPlantCollection{T,2},pricecurve::PriceCurve{T}) where T <: AbstractFloat
        return new{T}(plantdata,pricecurve)
    end
end
function ShortTermData(plantfilename::String,pricecurve::PriceCurve)
    ShortTermData(HydroPlantCollection(plantfilename),pricecurve)
end

function modelindices(data::ShortTermData, horizon::Horizon; areas::Vector{Area} = [0], rivers::Vector{River} = [:All])
    hours = collect(1:nhours(horizon))
    plants = plants_in_areas_and_rivers(hydrodata(data),areas,rivers)
    if isempty(plants)
        error("No plants in given set of price areas and rivers")
    end
    segments = collect(1:2)
    return ShortTermIndices(hours, plants, segments)
end

@hydromodel Deterministic ShortTerm = begin
    @unpack hours, plants, segments = indices
    hdata = hydrodata(data)
    λ = data.pricecurve
    λ̄ = expected(λ)
    HydroModels.horizon(λ) >= horizon || error("Not enough price data for chosen horizon")

    # Variables
    # ========================================================
    @variable(model, Q[p = plants, s = segments, t = hours], lowerbound = 0, upperbound = hdata[p].Q̄[s])
    @variable(model, S[p = plants, t = hours] >= 0)
    @variable(model, M[p = plants, t = hours], lowerbound = 0, upperbound = hdata[p].M̄)
    @variable(model, H[t = hours] >= 0)
    @variable(model, Qf[p = plants, t = hours] >= 0)
    @variable(model, Sf[p = plants, t = hours] >= 0)

    # Objectives
    # ========================================================
    # Net profit
    @expression(model, net_profit, sum(λ[t]*H[t] for t = hours))
    # Value of stored water
    @expression(model,value_of_stored_water,
                λ̄*sum(M[p,nhours(horizon)]*sum(hdata[i].μ[1]
                                                    for i = hdata.Qd[p])
                               for p = plants))
    # Define objective
    @objective(model, Max, net_profit + value_of_stored_water)

    # Constraints
    # ========================================================
    # Hydrological balance
    @constraint(model, hydro_constraints[p = plants, t = hours],
                # Previous reservoir content
                M[p,t] == (t > 1 ? M[p,t-1] : hdata[p].M₀)
                # Inflow
                + sum(Qf[i,t]
                      for i = intersect(hdata.Qu[p],plants))
                + sum(Sf[i,t]
                      for i = intersect(hdata.Su[p],plants))
                # Local inflow
                + hdata[p].V
                + (t <= hdata[p].Rqh ? 1 : 0)*sum(hdata[i].Q̃
                                                  for i = hdata.Qu[p])
                + (t == (hdata[p].Rqh + 1) ? 1 : 0)*sum(hdata[i].Q̃*(1-hdata[p].Rqm/60)
                                                         for i = hdata.Qu[p])
                # Outflow
                - sum(Q[p,s,t]
                      for s = segments)
                - S[p,t]
                )

    # Water flow: Discharge + Spillage
    @constraintref Qflow[1:length(plants),1:nhours(horizon)]
    @constraintref Sflow[1:length(plants),1:nhours(horizon)]
    for (pidx,p) = enumerate(plants)
        for t = hours
            if t - hdata[p].Rqh > 1
                Qflow[pidx,t] = @constraint(model,
                                            Qf[p,t] == (hdata[p].Rqm/60)*sum(Q[p,s,t-(hdata[p].Rqh+1)]
                                                                             for s = segments)
                                            + (1-hdata[p].Rqm/60)*sum(Q[p,s,t-hdata[p].Rqh]
                                                                      for s = segments)
                                            )
            elseif t - hdata[p].Rqh > 0
                Qflow[pidx,t] = @constraint(model,
                                            Qf[p,t] == (1-hdata[p].Rqm/60)*sum(Q[p,s,t-hdata[p].Rqh]
                                                                               for s = segments)
                                            )
            else
                Qflow[pidx,t] = @constraint(model,
                                            Qf[p,t] == 0
                                            )
            end
            if t - hdata[p].Rsh > 1
                Sflow[pidx,t] = @constraint(model,
                                            Sf[p,t] == (hdata[p].Rsm/60)*S[p,t-(hdata[p].Rsh+1)]
                                            + (1-hdata[p].Rsm/60)*S[p,t-hdata[p].Rsh]
                                            )
            elseif t - hdata[p].Rsh > 0
                Sflow[pidx,t] = @constraint(model,
                                            Sf[p,t] == (1-hdata[p].Rsm/60)*S[p,t-hdata[p].Rsh]
                                            )
            else
                Sflow[pidx,t] = @constraint(model,
                                            Sf[p,t] == 0
                                            )
            end
        end
    end

    # Power production
    @constraint(model, production[t = hours],
                H[t] == sum(hdata[p].μ[s]*Q[p,s,t]
                            for p = plants, s = segments)
                )
end


text:
using JuMP
using StructJuMP
using ECOS
using Cbc
using Base.Test

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

misocp_solver = CbcSolver()
socp_solver = ECOS.ECOSSolver(verbose=false)

@testset "[Benders] Empty scenario test" begin

    m = StructuredModel(num_scenarios=0)
    @variable(m, x, Int)
    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -20.0)

end

@testset "[Benders] Infeasible problem test" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 1)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 2)
    @variable(bl, y2 <= 2)
    @constraint(bl, x >= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Infeasible

end

@testset "[Benders] Infeasibility cut execution test #1" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 0)
    @variable(bl, y2 <= 2)
    @constraint(bl, x <= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -4.0)

end

@testset "[Benders] Optimality cut execution test #1" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 2)
    @variable(bl, y2 <= 4)
    @constraint(bl, x <= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -8.0)

end


text:
using diffprog
using Test


text:
"""
    ManyPencilArray{T,N,M}

Container holding `M` different [`PencilArray`](@ref) views to the same
underlying data buffer. All views share the same element type `T` and
dimensionality `N`.

This can be used to perform in-place data transpositions with
[`transpose!`](@ref Transpositions.transpose!).

---

    ManyPencilArray{T}(undef, pencils...; extra_dims=())

Create a `ManyPencilArray` container that can hold data of type `T` associated
to all the given [`Pencil`](@ref)s.

The optional `extra_dims` argument is the same as for [`PencilArray`](@ref).
"""
struct ManyPencilArray{
        T,  # element type of each array
        N,  # number of dimensions of each array (including extra_dims)
        M,  # number of arrays
        Arrays <: Tuple{Vararg{PencilArray,M}},
    }
    data :: Vector{T}
    arrays :: Arrays

    function ManyPencilArray{T}(
            init, pencils::Vararg{Pencil{Np},M};
            extra_dims::Dims=(),
        ) where {Np,M,T}
        data_length = max(length.(pencils)...) * prod(extra_dims)
        data = Vector{T}(init, data_length)
        arrays = _make_arrays(data, extra_dims, pencils...)
        A = typeof(arrays)
        N = Np + length(extra_dims)
        new{T,N,M,A}(data, arrays)
    end
end

function _make_arrays(data::Vector{T}, extra_dims::Dims, p::Pencil,
                      pens::Vararg{Pencil}) where {T}
    dims = (size_local(p, MemoryOrder())..., extra_dims...)
    n = prod(dims)
    @assert n == length_local(p) * prod(extra_dims)
    vec = view(data, Base.OneTo(n))
    arr = reshape(vec, dims)
    A = PencilArray(p, arr)
    (A, _make_arrays(data, extra_dims, pens...)...)
end

_make_arrays(::Vector, ::Dims) = ()

Base.eltype(A::ManyPencilArray{T}) where {T} = T
Base.ndims(A::ManyPencilArray{T,N}) where {T,N} = N

"""
    length(A::ManyPencilArray)

Returns the number of [`PencilArray`](@ref)s wrapped by `A`.
"""
Base.length(A::ManyPencilArray{T,N,M}) where {T,N,M} = M

"""
    first(A::ManyPencilArray)

Returns the first [`PencilArray`](@ref) wrapped by `A`.
"""
Base.first(A::ManyPencilArray) = first(A.arrays)

"""
    last(A::ManyPencilArray)

Returns the last [`PencilArray`](@ref) wrapped by `A`.
"""
Base.last(A::ManyPencilArray) = last(A.arrays)

"""
    getindex(A::ManyPencilArray, ::Val{i})
    getindex(A::ManyPencilArray, i::Integer)

Returns the i-th [`PencilArray`](@ref) wrapped by `A`.

If possible, the `Val{i}` form should be preferred, as it is more efficient and
it allows the compiler to know the return type.

See also [`first(::ManyPencilArray)`](@ref), [`last(::ManyPencilArray)`](@ref).

# Example

```julia
A = ManyPencilArray(pencil1, pencil2, pencil3)

# Get the PencilArray associated to `pencil2`.
# u2 = A[2]
u2 = A[Val(2)]  # faster!
```
"""
Base.getindex(A::ManyPencilArray, ::Val{i}) where {i} =
    _getindex(Val(i), A.arrays...)
@inline Base.getindex(A::ManyPencilArray, i) = A[Val(i)]

@inline function _getindex(::Val{i}, a, t::Vararg) where {i}
    i :: Integer
    i <= 0 && throw(BoundsError("index must be >= 1"))
    i == 1 && return a
    _getindex(Val(i - 1), t...)
end

# This will happen if the index `i` intially passed is too large.
@inline _getindex(::Val) = throw(BoundsError("invalid index"))


text:
## Data structure for tracking sparsity

"""
The sparsity pattern.

- `I`: Input index
- `J`: Ouput index

`(i, j)` means the `j`th element of the output depends on
the `i`th element of the input. Therefore `length(I) == length(J)`
"""
struct Sparsity
    m::Int
    n::Int
    I::Vector{Int} # Input
    J::Vector{Int} # Output
end

SparseArrays.sparse(s::Sparsity) = sparse(s.I, s.J, true, s.m, s.n)

Sparsity(m, n) = Sparsity(m, n, Int[], Int[])

function Base.push!(S::Sparsity, i::Int, j::Int)
    push!(S.I, i)
    push!(S.J, j)
end

struct ProvinanceSet
    set::Set{Int}
    ProvinanceSet(s::Set) = new(s)
    ProvinanceSet(s) = new(Set(s))
end

# note: this is not strictly set union, just some efficient way of concating
Base.union(p::ProvinanceSet, ::Cassette.NoMetaData) = p
Base.union(::Cassette.NoMetaData, p::ProvinanceSet) = p

Base.union(p::ProvinanceSet,
           q::ProvinanceSet) = ProvinanceSet(union(p.set, q.set))
Base.union(p::ProvinanceSet,
           q::ProvinanceSet,
           rs::ProvinanceSet...) = union(union(p, q), rs...)
Base.union(p::ProvinanceSet) = p

function Base.push!(S::Sparsity, i::Int, js::ProvinanceSet)
    for j in js.set
        push!(S, i, j)
    end
end

# Cassette

@proptagcontext JacobianSparsityContext

struct JacInput end
struct JacOutput end

const TagType = Union{JacInput,
                      JacOutput,
                      ProvinanceSet}

istainted(ctx::JacobianSparsityContext, val) = metatype(val, ctx) <: ProvinanceSet

Cassette.metadatatype(::Type{<:JacobianSparsityContext},
                      ::DataType) = TagType
# optimization
Cassette.metadatatype(::Type{<:JacobianSparsityContext},
                      ::Type{<:Number}) = ProvinanceSet

function propagate_tags(ctx::JacobianSparsityContext,
                        f, result, args...)

    # e.g. X .+ X[1]
    # XXX: This wouldn't be required if we didn't have Input
    #
    any(x->metatype(x, ctx) <: JacInput,
        args) && return result

    idxs = findall(args) do x
        metatype(x, ctx) <: ProvinanceSet
    end

    if isempty(idxs)
        return result
    else
        tag(untag(result, ctx),
            ctx,
            union(map(x->metadata(x, ctx), args[idxs])...))
    end
end
const TaggedOf{T} = Cassette.Tagged{A, T} where A

function Cassette.overdub(ctx::JacobianSparsityContext,
                          f::typeof(getindex),
                          X::Tagged,
                          idx::Union{TaggedOf{Int},Int}...)
    if metatype(X, ctx) <: JacInput
        i = LinearIndices(untag(X, ctx))[idx...]
        val = Cassette.fallback(ctx, f, X, idx...)
        tag(val, ctx, ProvinanceSet(i))
    else
        Cassette.recurse(ctx, f, X, idx...)
    end
end

# setindex! on the output
function Cassette.overdub(ctx::JacobianSparsityContext,
                          f::typeof(setindex!),
                          Y::Tagged,
                          val::Tagged,
                          idx::Int...)
    S = ctx.metadata[1]
    if metatype(Y, ctx) <: JacOutput
        set = metadata(val, ctx)
        if set isa ProvinanceSet
            i = LinearIndices(untag(Y, ctx))[idx...]
            push!(S, i, set)
        end
        Cassette.fallback(ctx, f, Y, val, idx...)
    else
        Cassette.recurse(ctx, f, Y, val, idx...)
    end
end

function jacobian_sparsity(f!, Y, X, args...;
                           sparsity=Sparsity(length(Y), length(X)),
                           verbose = true,
                           raw = false)

    ctx = JacobianSparsityContext(metadata=sparsity)
    ctx = Cassette.enabletagging(ctx, f!)
    ctx = Cassette.disablehooks(ctx)

    res = nothing
    abstract_run((result)->(res=result),
                 ctx,
                 f!,
                 tag(Y, ctx, JacOutput()),
                 tag(X, ctx, JacInput()),
                 map(arg -> arg isa Fixed ?
                     arg.val : tag(arg, ctx, ProvinanceSet(())), args)...;
                 verbose=verbose)

    if raw
        return (ctx, res)
    else
        return sparse(sparsity)
    end
end

function Cassette.overdub(ctx::JacobianSparsityContext,
                          f::typeof(Base.unsafe_copyto!),
                          X::Tagged,
                          xstart,
                          Y::Tagged,
                          ystart,
                          len)
    S = ctx.metadata[1]
    if metatype(Y, ctx) <: JacInput && metatype(X, ctx) <: JacOutput
        # Write directly to the output sparsity
        val = Cassette.fallback(ctx, f, X, xstart, Y, ystart, len)
        for (i, j) in zip(xstart:xstart+len-1, ystart:ystart+len-1)
            push!(S, i, j)
        end
        val
    elseif metatype(Y, ctx) <: JacInput
        # Keep around a ProvinanceSet
        val = Cassette.fallback(ctx, f, X, xstart, Y, ystart, len)
        nometa = Cassette.NoMetaMeta()
        rhs = (i->Cassette.Meta(ProvinanceSet(i), nometa)).(ystart:ystart+len-1)
        X.meta.meta[xstart:xstart+len-1] .= rhs
        val
    elseif metatype(X, ctx) <: JacOutput
        val = Cassette.fallback(ctx, f, X, xstart, Y, ystart, len)
        for (i, j) in zip(xstart:xstart+len-1, ystart:ystart+len-1)
            y = Cassette.@overdub ctx Y[j]
            set = metadata(y, ctx)
            if set isa ProvinanceSet
                push!(S, i, set)
            end
        end
        val
    else
        val = Cassette.fallback(ctx, f, X, xstart, Y, ystart, len)
        for (i, j) in zip(xstart:xstart+len-1, ystart:ystart+len-1)
            y = Cassette.@overdub ctx Y[j]
            set = metadata(y, ctx)
            nometa = Cassette.NoMetaMeta()
            X.meta.meta[i] = Cassette.Meta(set, nometa)
        end
        val
    end
end


text:
using Documenter, Weave

makedocs(modules = Weave)
include("make_examples.jl")
run(`mkdocs build`)


text:
@testset "Consider Kalman filter" begin

    @testset "Measurement update" begin

        x = randn(6)
        PL = randn(6,6)
        P = PL'PL
        P_chol = cholesky(P)
        RL = randn(3,3)
        R = RL'RL
        y = randn(3)
        R_chol = cholesky(R)
        H = randn(3,6)
        h(x) = H * x

        # Comparison measurement update
        # Nulling Kalman Gain for considered states
        # http://www.anuncommonlab.com/articles/how-kalman-filters-work/part2.html
        ỹ = KalmanFilter.calc_innovation(H, x, y)
        PHᵀ = KalmanFilter.calc_P_xy(P, H)
        S = KalmanFilter.calc_innovation_covariance(H, P, R)
        K = KalmanFilter.calc_kalman_gain(PHᵀ, S, [])
        K[5:6,:] .= 0.0
        x_post = KalmanFilter.calc_posterior_state(x, K, ỹ, [])
        # Use Joseph form, because of nulling Kalman gain
        P_post = (I - K * H) * P * (I - K * H)' + K * R * K'

        cons_mu = @inferred measurement_update(x, P, y, H, R, consider = 5:6)
        @test @inferred(get_covariance(cons_mu)) ≈ P_post
        @test @inferred(get_state(cons_mu)) ≈ x_post
        @test @inferred(get_kalman_gain(cons_mu)) ≈ K[1:4,:]

        cons_mu_chol = measurement_update(x, P_chol, y, H, R_chol, consider = 5:6)
        @test @inferred(get_covariance(cons_mu_chol)) ≈ P_post
        @test @inferred(get_state(cons_mu_chol)) ≈ x_post
        @test @inferred(get_kalman_gain(cons_mu_chol)) ≈ K[1:4,:]

        cons_mu_ukf = measurement_update(x, P, y, h, R, consider = 5:6)
        @test @inferred(get_covariance(cons_mu_ukf)) ≈ P_post
        @test @inferred(get_state(cons_mu_ukf)) ≈ x_post
        @test @inferred(get_kalman_gain(cons_mu_ukf)) ≈ K[1:4,:]

        cons_mu_srukf = measurement_update(x, P_chol, y, h, R_chol, consider = 5:6)
        @test @inferred(get_covariance(cons_mu_srukf)) ≈ P_post
        @test @inferred(get_state(cons_mu_srukf)) ≈ x_post
        @test @inferred(get_kalman_gain(cons_mu_srukf)) ≈ K[1:4,:]
    end
end


text:
function select(::Type{Neg}, y, s)
    inds = findall(y .== 0)
    return view(s, inds), inds
end

function select(::Type{Pos}, y, s)
    inds = findall(y .== 1)
    return view(s, inds), inds
end

# Thresholds
function threshold(tp::Threshold, y::AbstractMatrix, s::AbstractMatrix; kwargs...)
    ts, ~ = find_threshold(tp, y, s; kwargs...)
    return ts
end

function ChainRulesCore.rrule(::typeof(threshold), tp::Threshold, y, s; kwargs...)
    ts, inds = find_threshold(tp, y, s; kwargs...)

    function threshold_pullback(Δ)
        Δt = zeros(Bool, size(s)...)
        Δt[CartesianIndex.(1:length(inds), inds)] .= 1
        return NoTangent(), NoTangent(), NoTangent(), Δ .* convert(typeof(s), Δt)
    end
    return ts, threshold_pullback
end

function find_threshold(
    tp::Threshold,
    y::AbstractMatrix,
    s::AbstractMatrix;
    update_buffer = true
)

    k = size(s, 1)
    ts, inds = zeros(eltype(s), k), zeros(Int, k)

    @inbounds for i in 1:k
        ts[i], inds[i] = find_threshold(tp, Vector(y[i, :]), Vector(s[i, :]))
    end
    update_buffer && update_buffer!(ts, inds)
    return convert(typeof(similar(s, k)), ts), inds
end

function find_threshold(
    I::Type{<:Indices},
    find_func,
    y_in::AbstractVector,
    s_in::AbstractVector,
    args...,
)
    s, inds = select(I, y_in, s_in)
    t, ind = find_func(s, args...)
    return t, inds[ind]
end

function find_threshold(
    ::Type{All},
    find_func,
    ::AbstractVector,
    s_in::AbstractVector,
    args...,
)

    return find_func(s_in, args...)
end

function find_kth(x, k::Int, rev::Bool)
    ind = partialsortperm(x, k; rev)
    return x[ind], ind
end

function find_kth(x, ks, rev::Bool)
    inds = partialsortperm(x, 1:maximum(ks); rev)
    return x[inds[ks]], inds[ks]
end

function find_quantile(x, τ, rev::Bool)
    if !all(0 .<= τ .<= 1)
        throw(ArgumentError("input probability out of [0,1] range"))
    end

    n = length(x)
    i = rev ? 1 .- τ : τ
    ks = min.(max.(1, round.(Int64, n .* i)), n)

    if maximum(ks) <= n/2
        return find_kth(x, ks, false)
    else
        return find_kth(x, n .- ks .+ 1, true)
    end
end

# Maximum and minimum scores
struct Maximum{I<:Indices} <: Threshold
    Maximum(indices = Neg) = new{indices}()
end

function Base.show(io::IO, ::Maximum{I}) where {I}
    from = I <: All ? "" : "$(I) "
    return print(io, "largest $(from)score")
end

function find_threshold(::Maximum{I}, y::AbstractVector, s::AbstractVector) where {I}
    return find_threshold(I, findmax, y, s)
end

struct Minimum{I<:Indices} <: Threshold
    Minimum(indices = Pos) = new{indices}()
end

function Base.show(io::IO, ::Minimum{I}) where {I}
    from = I <: All ? "" : "$(I) "
    return print(io, "smallest $(from)score")
end

function find_threshold(::Minimum{I}, y::AbstractVector, s::AbstractVector) where {I}
    return find_threshold(I, findmin, y, s)
end

# K-th laregst/smallest scores
struct Kth{I<:Indices} <: Threshold
    K::Int64
    rev::Bool

    Kth(K, indices = Neg; rev = true) = new{indices}(K, rev)
end

function ordinal(i::Int)
    j = i % 10
    k = i % 100
    if j == 1 && k != 11
        return "$(i)st"
    elseif j == 2 && k != 12
        return  "$(i)nd"
    elseif j == 3 && k != 13
        return  "$(i)rd"
    else
        return "$(i)th"
    end
end

function Base.show(io::IO, tp::Kth{I}) where {I}
    ord = tp.rev ? "largest" : "smallest"
    from = I <: All ? "" : "$(I) "
    return print(io, "$(ordinal(tp.K)) $(ord) $(from)score")
end

function find_threshold(tp::Kth{I}, y::AbstractVector, s::AbstractVector) where {I}
    return find_threshold(I, find_kth, y, s, tp.K, tp.rev)
end

# Quantiles
struct Quantile{I<:Indices} <: Threshold
    τ::Float64
    rev::Bool

    Quantile(τ, indices = Neg; rev = true) = new{indices}(τ, rev)
end

function Base.show(io::IO, tp::Quantile{I}) where {I}
    ord = tp.rev ? "top" : "bottom"
    print(io, "$(ord) $(tp.τ)-quantile from $(I) scores")
end

function find_threshold(tp::Quantile{I}, y::AbstractVector, s::AbstractVector) where {I}
    return find_threshold(I, find_quantile, y, s, tp.τ, tp.rev)
end

struct SampledQuantile{I<:Indices} <: Threshold
    sampler
    rev::Bool

    SampledQuantile(sampler, indices = Neg; rev = true) = new{indices}(sampler, rev)
end

function Base.show(io::IO, tp::SampledQuantile{I}) where {I}
    ord = tp.rev ? "top" : "bottom"
    print(io, "$(ord) sampled quantile from $(I) scores")
end

function find_threshold(tp::SampledQuantile{I}, y::AbstractVector, s::AbstractVector) where {I}
    return find_threshold(I, find_quantile, y, s, tp.sampler(), tp.rev)
end


text:
"""Abstract type for container modules -Modules that specify an
execution structure for a given group of modules"""
abstract type KnetContainer <: KnetModule end

"""
Sequential <: KnetModule

# Constructor
    Sequential(ls...) adds layers in ls.
# Fields
    layers::Array{Union{KnetModule, Function}, 1}

# Usage
    s = Sequential(Linear(...), ReLU(), Linear(...),...)
    @mc s(x)
    @run r(x)

 Layers stored in layers field is executed in order, where layer
 n+1 takes layer n's output as its input. Layer inputs and outputs
 should be consistent.

See `add!` function to add additional layers after construction

# Note on storing functions in layers
    Currently, JLD throws error if you add functions in layers,
    but it is still kept for convenience, considering the future
    support. Currently, using functional modules is a better practice.
    See FnModule types.
"""
type Sequential <: KnetContainer
    layers::Array{Union{KnetModule, Function}, 1}
    #layers::Array{KnetModule, 1}
end

function Sequential(ls...)
    new = Sequential([l for l in ls])
    if any(x->isa(x, Function), new.layers)
        warn("Serialization functions like save_module and load_module",
             " may not work when you use functions in Sequential")
    end
    return new
end

"""
add!(s::Sequential, m...) adds modules m... to the end of layers.
Identical to push!(s.layers, m...)
"""
add!(s::Sequential, m...) = push!(s.layers, m...)

function (s::Sequential)(ctx, x)
    for l in s.layers
        if isa(l, Function) || isa(l, FnModule)
            x = l(x)
        else
            x = l(ctx, x)
        end
    end
    return x
end


text:
module SubTriangulationsTests

using Test

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using STLCutters

using STLCutters: Polyhedron
using STLCutters: STL
using STLCutters: restrict
using STLCutters: clip
using STLCutters: split
using STLCutters: surface, volume
using STLCutters: get_facet_planes
using STLCutters: get_reflex_planes
using STLCutters: compute_distances!
using STLCutters: get_original_reflex_faces
using STLCutters: get_cell_planes
using STLCutters: filter_face_planes
using STLCutters: refine
using STLCutters: get_cell_nodes_to_inout
using STLCutters: read_stl
using STLCutters: merge_nodes
using STLCutters: compute_stl_model
using STLCutters: compute_grid
using STLCutters: FACE_IN, FACE_OUT, FACE_CUT

vertices = [
  Point(0.1,-0.2,0.5),
  Point(1.0,-0.3,0.5),
  Point(-0.2,1.2,0.5),
  Point(-0.3,0.3,0.5),
  Point(0.5,1.3,0.5),
  Point(0.5,0.5,0.7),
  Point(1.2,1.2,0.5)]

facet_to_vertices =
[[4,6,1],
 [6,5,7],
 [4,3,6],
 [1,6,2],
 [3,5,6],
 [6,7,2] ]

stlmodel = compute_stl_model( Table(facet_to_vertices), vertices )
stl = STL(stlmodel)

# Global setup
Πf = get_facet_planes(stl)
Πr = get_reflex_planes(stl,Πf)

# Setup cell
stl_facets_k = 1:num_cells(stl)
K = Polyhedron(HEX)

Γ0 = Polyhedron(stl)
Γk0 = restrict(Γ0,stl,stl_facets_k)

stl_reflex_faces_k = get_original_reflex_faces(Γk0,stl)
stl_facets_k = stl_facets_k.+get_offset(stl,num_dims(stl))

Πk,Πk_ids,Πk_io = get_cell_planes(HEX,Point(0,0,0),Point(1,1,1))
Πf,Πf_ids = filter_face_planes(stl,Πr,stl_reflex_faces_k,Πf,stl_facets_k)

compute_distances!(Γk0,lazy_append(Πk,Πf),lazy_append(Πk_ids,Πf_ids))
compute_distances!(K,Πf,Πf_ids)

Γk = clip(Γk0,Πk_ids,inout=Πk_io)

Kn_in = refine(K,Γk,stl,stl_reflex_faces_k,inside=true)
Kn_out = refine(K,Γk,stl,stl_reflex_faces_k,inside=false)

T_in,X_in = simplexify(Kn_in)
T_out,X_out = simplexify(Kn_out)
T_Γ,X_Γ = simplexify(Γk)

n_to_io = get_cell_nodes_to_inout(Kn_in,Kn_out,HEX)

# Test and Write

@test n_to_io == [fill(FACE_OUT,num_vertices(HEX)÷2);fill(FACE_IN,num_vertices(HEX)÷2)]
mesh_in = compute_grid(T_in,X_in,TET)
mesh_out = compute_grid(T_out,X_out,TET)
mesh_Γ = compute_grid(T_Γ,X_Γ,TRI)

#  writevtk(Γk,"Gk")
#  writevtk(Kn_in,"Kin")
#  writevtk(Kn_out,"Kin")
#  writevtk(mesh_in,"mesh_in")
#  writevtk(mesh_out,"mesh_out")
#  writevtk(mesh_Γ,"mesh_G")

@test volume(mesh_in) + volume(mesh_out) ≈ 1

# Sharp corner

vertices = [
  Point(0.1,-0.2,-1.0),
  Point(1.0,-0.3,-1.0),
  Point(-0.2,1.2,-1.0),
  Point(-0.3,0.3,-1.0),
  Point(0.5,1.3,-1.0),
  Point(0.5,0.5,0.5),
  Point(1.2,1.2,-1.0)]

facet_to_vertices =
[[4,6,1],
 [6,5,7],
 [4,3,6],
 [1,6,2],
 [3,5,6],
 [6,7,2] ]

stlmodel = compute_stl_model( Table(facet_to_vertices), vertices )
stl = STL(stlmodel)

# Global setup
Πf = get_facet_planes(stl)
Πr = get_reflex_planes(stl,Πf)

# Setup cell
stl_facets_k = 1:num_cells(stl)
K = Polyhedron(HEX)

Γ0 = Polyhedron(stl)
Γk0 = restrict(Γ0,stl,stl_facets_k)

stl_reflex_faces_k = get_original_reflex_faces(Γk0,stl)
stl_facets_k = stl_facets_k.+get_offset(stl,num_dims(stl))

Πk,Πk_ids,Πk_io = get_cell_planes(HEX,Point(0,0,0),Point(1,1,1))
Πf,Πf_ids = filter_face_planes(stl,Πr,stl_reflex_faces_k,Πf,stl_facets_k)

compute_distances!(Γk0,lazy_append(Πk,Πf),lazy_append(Πk_ids,Πf_ids))
compute_distances!(K,Πf,Πf_ids)

Γk = clip(Γk0,Πk_ids,inout=Πk_io)

Kn_in = refine(K,Γk,stl,stl_reflex_faces_k,inside=true)
Kn_out = refine(K,Γk,stl,stl_reflex_faces_k,inside=false)

T_in,X_in = simplexify(Kn_in)
T_out,X_out = simplexify(Kn_out)
T_Γ,X_Γ = simplexify(Γk)

n_to_io = get_cell_nodes_to_inout(Kn_in,Kn_out,HEX)

# Test and Write

@test n_to_io == fill(FACE_IN,num_vertices(HEX))

mesh_in = compute_grid(T_in,X_in,TET)
mesh_out = compute_grid(T_out,X_out,TET)

#  writevtk(Γk,"Gk")
#  writevtk(Kn_in,"Kin")
#  writevtk(Kn_out,"Kin")
#  writevtk(mesh_in,"mesh_in")
#  writevtk(mesh_out,"mesh_out")
#  writevtk(mesh_Γ,"mesh_G")

@test volume(mesh_in) + volume(mesh_out) ≈ 1

## Real STL

X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))
stl = compute_stl_model(T,X)
stl = merge_nodes(stl)

#writevtk(get_grid(stl),"stl")

p = HEX
δ = 0.2
n = 10
D = 3

pmin,pmax = get_bounding_box(stl)
Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)
grid = get_grid(model)

@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=false)

#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])

bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
submesh_out_vol = volume(subcells,labels.cell_to_io,:out)

in_volume = bgmesh_in_vol + submesh_in_vol
out_volume = bgmesh_out_vol + submesh_out_vol
cut_volume = bgmesh_cut_vol

#celldata = [ "inoutcut" => labels.bgcell_to_ioc ]
#writevtk( grid, "bgcells"; celldata )
#
#celldata = [ "inout" => labels.cell_to_io, "bgcell" => labels.cell_to_bgcell ]
#writevtk( subcells, "subcells"; celldata )
#
#celldata = [ "bgcell" => labels.face_to_bgcell ]
#writevtk( subfaces, "subfaces" )

println("Num subcells: $(num_cells(subcells))")
@test surface(get_grid(stl)) ≈ surface(subfaces)
@test submesh_in_vol + submesh_out_vol ≈ cut_volume
@test in_volume + out_volume ≈ volume(grid)
@test in_volume ≈ 273280.03374196636

X,T,N = read_stl(joinpath(@__DIR__,"data/wine_glass.stl"))
stl = compute_stl_model(T,X)
stl = merge_nodes(stl)

#writevtk(get_grid(stl),"stl")

p = HEX
δ = 0.2
n = 20
D = 3

pmin,pmax = get_bounding_box(stl)
Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)
grid = get_grid(model)

@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=false)

#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])

bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
submesh_out_vol = volume(subcells,labels.cell_to_io,:out)

in_volume = bgmesh_in_vol + submesh_in_vol
out_volume = bgmesh_out_vol + submesh_out_vol
cut_volume = bgmesh_cut_vol

println("Num subcells: $(num_cells(subcells))")
@test surface(get_grid(stl)) ≈ surface(subfaces)
@test submesh_in_vol + submesh_out_vol ≈ cut_volume
@test in_volume + out_volume ≈ volume(grid)
@test in_volume ≈ 74.12595970214333

## Kd-Tree

X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))
stl = compute_stl_model(T,X)
stl = merge_nodes(stl)

#writevtk(get_grid(stl),"stl")

p = HEX
δ = 0.2
n = 10
D = 3

pmin,pmax = get_bounding_box(stl)
Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)
grid = get_grid(model)

@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=true)

#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])

bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
submesh_out_vol = volume(subcells,labels.cell_to_io,:out)

in_volume = bgmesh_in_vol + submesh_in_vol
out_volume = bgmesh_out_vol + submesh_out_vol
cut_volume = bgmesh_cut_vol

println("Num subcells: $(num_cells(subcells))")
@test surface(get_grid(stl)) ≈ surface(subfaces)
@test submesh_in_vol + submesh_out_vol ≈ cut_volume
@test in_volume + out_volume ≈ volume(grid)
@test in_volume ≈ 273280.03374196636

X,T,N = read_stl(joinpath(@__DIR__,"data/wine_glass.stl"))
stl = compute_stl_model(T,X)
stl = merge_nodes(stl)

#writevtk(get_grid(stl),"stl")

p = HEX
δ = 0.2
n = 20
D = 3

pmin,pmax = get_bounding_box(stl)
Δ = (pmax-pmin)*δ
pmin = pmin - Δ
pmax = pmax + Δ
partition = (n,n,n)

model = CartesianDiscreteModel(pmin,pmax,partition)
grid = get_grid(model)

@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=true)

#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])

bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
submesh_out_vol = volume(subcells,labels.cell_to_io,:out)

in_volume = bgmesh_in_vol + submesh_in_vol
out_volume = bgmesh_out_vol + submesh_out_vol
cut_volume = bgmesh_cut_vol

println("Num subcells: $(num_cells(subcells))")
@test surface(get_grid(stl)) ≈ surface(subfaces)
@test submesh_in_vol + submesh_out_vol ≈ cut_volume
@test in_volume + out_volume ≈ volume(grid)
@test in_volume ≈ 74.12595970214333
end # module