text:
"""
    $(SIGNATURES)

Converts units like kilojoule / mole into kJ/mol.
This function is not exported.
"""
function _parse_units(x)
    unit_conversions = [
        "kilo" => "k",
        "joule" => "J",
        "mole" => "mol",
        "milli" => "m",
        "kelvin" => "K",
        "molar" => "M",
        "micro" => "μ",
        "nano" => "n",
        "V" => "V",
    ]
    return foldl(replace, unit_conversions; init = x)
end

"""
    $(SIGNATURES)

Parse a reaction string and insert the prefix in front of each metabolite.
This function is not exported.
"""
function _parse_reaction_string(str::String, prefix::String)
    rxn_string_reagents_products = split(str, "=")
    if length(rxn_string_reagents_products) == 1
        return prefix * ":" * str # only metabolite
    else
        reagents = string(rxn_string_reagents_products[1])
        products = string(rxn_string_reagents_products[2])
        parsed_reagents = join(_parse_reaction_part(reagents, prefix * ":"), " + ")
        parsed_products = join(_parse_reaction_part(products, prefix * ":"), " + ")
        return parsed_reagents * " = " * parsed_products
    end
end

"""
    $(SIGNATURES)

Helper function for `_parse_reaction_string`.
This function is not exported.
"""
function _parse_reaction_part(str::String, prefix::String)
    isnothing(findfirst(isletter, str)) && return String[]
    nospace_str = replace(str, "*" => " ") # get rid of * if present
    mets = split(nospace_str, r"\+")
    parsed_str = String[]
    for met in mets
        met = replace(met, r"\s+$" => "")
        ind = findfirst(r"(\d+\s)|(\d+\.\s)|(\d+\.\d+\s)|(\.\d+\s)", met)
        if isnothing(ind)
            push!(parsed_str, prefix * replace(met, " " => ""))
        else
            stoich = met[ind]
            push!(parsed_str, stoich * prefix * replace(met[ind[end]+1:end], " " => ""))
        end
    end
    return parsed_str
end

"""
    $(SIGNATURES)

Return a reaction string that has `"bigg.metabolite:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.
A string literal version also exists: `bigg"..."`

See also [`kegg`](@ref), [`metanetx`](@ref), [`chebi`](@ref).

# Example
```
bigg("atp + h2o = adp + pi")
bigg("atp + h2o = adp + 2*pi")
bigg("atp + h2o = adp + 2 pi")
```
"""
function bigg(str)
    _parse_reaction_string(str, "bigg.metabolite")
end

"""
    $(SIGNATURES)

Return a reaction string that has `"bigg.metabolite:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.

See also [`kegg`](@ref), [`metanetx`](@ref), [`chebi`](@ref).

# Example
```
bigg"atp + h2o = adp + pi"
bigg"atp + h2o = adp + 2*pi"
bigg"atp + h2o = adp + 2 pi"
```
"""
macro bigg_str(str)
    return :(bigg($str))
end

"""
    $(SIGNATURES)

Return a reaction string that has `"kegg:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So 13pgm is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.
A string literal version also exists: `kegg"..."`

See also [`bigg`](@ref), [`metanetx`](@ref), [`chebi`](@ref).

# Example
```
kegg("C00002 + C00001 = C00008 + C00009")
kegg("C00002 + C00001 = C00008 + 2*C00009")
kegg("C00002 + C00001 = C00008 + 2 C00009")
```
"""
function kegg(str::String)
    _parse_reaction_string(str, "kegg")
end

"""
    $(SIGNATURES)

Return a reaction string that has `"kegg:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.

See also [`bigg`](@ref), [`metanetx`](@ref), [`chebi`](@ref).

# Example
```
kegg"C00002 + C00001 = C00008 + C00009"
kegg"C00002 + C00001 = C00008 + 2*C00009"
kegg"C00002 + C00001 = C00008 + 2 C00009"
```
"""
macro kegg_str(str::String)
    return :(kegg($str))
end

"""
    $(SIGNATURES)

Return a reaction string that has `"metanetx.chemical:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.
A string literal version also exists: `metanetx"..."`

See also [`bigg`](@ref), [`kegg`](@ref), [`chebi`](@ref)

# Example
```
metanetx("MNXM3 + MNXM2 = MNXM7 + MNXM9")
metanetx("MNXM3 + MNXM2 = MNXM7 + 2 MNXM9")
metanetx("MNXM3 + MNXM2 = MNXM7 + 2*MNXM9")
```
"""
function metanetx(str::String)
    _parse_reaction_string(str, "metanetx.chemical")
end

"""
    $(SIGNATURES)

Return a reaction string that has `"metanetx.chemical:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.

See also [`bigg`](@ref), [`kegg`](@ref), [`chebi`](@ref).

# Example
```
metanetx"MNXM3 + MNXM2 = MNXM7 + MNXM9"
metanetx"MNXM3 + MNXM2 = MNXM7 + 2 MNXM9"
metanetx"MNXM3 + MNXM2 = MNXM7 + 2*MNXM9"
```
"""
macro metanetx_str(str::String)
    return :(metanetx($str))
end

"""
    $(SIGNATURES)

Return a reaction string that has `"CHEBI:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.
A string literal version also exists: `chebi"..."`

See also [`bigg`](@ref), [`metanetx`](@ref), [`kegg`](@ref).

# Example
```
chebi(30616 + 33813 = 456216 + 43474")
chebi(30616 + 33813 = 456216 + 2*43474")
chebi(30616 + 33813 = 456216 + 2 43474")
```
"""
function chebi(str::String)
    _parse_reaction_string(str, "CHEBI")
end

"""
    $(SIGNATURES)

Return a reaction string that has `"CHEBI:"` appended in front of each metabolite.
Respects stoichiometric coefficients. `"*"` in front of metabolites do not matter, but spaces do matter.
So `13pgm` is different from `13 pgm` but `13*pgm` is the same as `13 pgm`.

See also [`bigg`](@ref), [`metanetx`](@ref), [`kegg`](@ref).

# Example
```
chebi"30616 + 33813 = 456216 + 43474"
chebi"30616 + 33813 = 456216 + 2*43474"
chebi"30616 + 33813 = 456216 + 2 43474"
```
"""
macro chebi_str(str::String)
    return :(chebi($str))
end

"""
    $(SIGNATURES)

Caps the lower bound of a concentration to ϵ.
Return ϵ if x is less than ϵ, otherwise x.
This function is not exported.
"""
function _lower_bound(x, ϵ = 0.001u"mM")
    x < ϵ ? ϵ : x
end


text:
using ChainRulesCore

#@non_differentiable leg_analysis(::Any...)
#@non_differentiable get_size_dict(::Any...)
#@non_differentiable optimaltree(::Any...)

ext: jl


text:
include("iPredict.jl/src/MixtureModelling.jl")
using DelimitedFiles
using Distributions
using StatsPlots
using DataFrames

test_data_1 = rand(LogNormal(1, 0.5), 20)
test_data_2 = rand(LogNormal(3, 0.5), 60)
append!(test_data_1, test_data_2)

relationships = Array{Float64}(undef, length(test_data_1), 2)
mu_arr = [1, 3]
sigma_arr = [0.5, 0.5]
k_arr = [0.5, 0.5]

expectation!(relationships, length(test_data_1), 2, mu_arr, sigma_arr, k_arr, test_data_1)
maximisation_k!(relationships, k_arr)

append!(test_data_1, rand(LogNormal(2, 0.7), 80))

relationships = Array{Float64}(undef, length(test_data_1), 3)
mu_arr = [1, 3, 2]
sigma_arr = [0.5, 0.5, 0.7]
k_arr = [1/3, 1/3, 1/3]

expectation!(relationships, length(test_data_1), 3, mu_arr, sigma_arr, k_arr, test_data_1)
maximisation_k!(relationships, k_arr)

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
module PdfFraud

using Flux
using ProtoBuf
using GenerativeModels
import GenerativeModels: elbo
using GMExtensions
using ValueHistories
using ProgressMeter
using ConditionalDists
using IPMeasures
using JLD2, FileIO, JSON, CSV, BSON
using MLDataPattern
using DrWatson

export PdfSampleRaw

include("data.jl")
include("experiments.jl")

end


text:
using JuliaTutor
using Test

@testset "JuliaTutor.jl" begin
    # Write your tests here.
end


text:
using Gen

@gen function my_model(xs)
    slope = @trace(normal(0, 2), :slope)
    intercept = @trace(normal(0, 10), :intercept)
    for (i, x) in enumerate(xs)
        @trace(normal(slope * x + intercept, 1), "y-$i")
    end
end

function my_inference_program(xs, ys, num_iters)
    constraints = choicemap()
    for (i, y) in enumerate(ys)
        constraints["y-$i"] = y
    end
    (trace, _) = generate(my_model, (xs,), constraints)
    slope_selection = select(:slope)
    intercept_selection = select(:intercept)
    for iter=1:num_iters
        (trace, _) = metropolis_hastings(trace, slope_selection)
        (trace, _) = metropolis_hastings(trace, intercept_selection)
    end
    return (trace[:slope], trace[:intercept])
end

xs = [1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]
ys = [8.23, 5.87, 3.99, 2.59, 0.23, -0.66, -3.53, -6.91, -7.24, -9.90]
(slope, intercept) = my_inference_program(xs, ys, 100)
println("slope: $slope, intercept: $intercept")


text:
# Packages. Note: make sure that all the packages that you use here are
# (1) available in either your global environment or your current page environment
# (2) are installed in the `.github/workflows/deploy.yml` file. See both the
# Project.toml and the `deploy.yml` file here as examples.
#
using DelimitedFiles
using TikzCDs
using Dates

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

###########
### 001 ###
###########

# approach 1; uses DelimitedFiles
function hfun_members_table(params::Vector{String})::String
    path_to_csv = params[1]
    members = readdlm(path_to_csv, ',', skipstart=1)
    # write a simple table
    io = IOBuffer()
    write(io, "
")
    write(io, "
Name	GitHub alias
")
    for (name, alias) in eachrow(members)
        write(io, "
")
        write(io, "$name	")
        write(io, """$alias	""")
        write(io, "
")
    end
    write(io, "
")
    return String(take!(io))
end

###########
### 007 ###
###########

# case 1
hfun_case_1() =
    """
var read from foo is $(pagevar("foo", "var"))
"""

# case 2, note the `@delay`
@delay function hfun_case_2()
    all_tags = globvar("fd_page_tags")
    (all_tags === nothing) && return ""
    all_tags = union(values(all_tags)...)
    tagstr = strip(prod("$t " for t in all_tags))
    return """
tags: { $tagstr }
"""
end

###########
### 008 ###
###########

function lx_capa(com, _)
    # this first line extracts the content of the brace
    content = Franklin.content(com.braces[1])
    output = replace(content, "a" => "A")
    return "**$output**"
end

function env_cap(com, _)
    content = Franklin.content(com)
    option = Franklin.content(com.braces[1])
    output = replace(content, option => uppercase(option))
    return "~~~~~~$output~~~~~~"
end

###########
### 009 ###
###########

function env_tikzcd(e, _)
  content = strip(Franklin.content(e))
  name = strip(Franklin.content(e.braces[1]))
  # save SVG at __site/assets/[path/to/file]/$name.svg
  rpath = joinpath("assets", splitext(Franklin.locvar(:fd_rpath))[1], "$name.svg")
  outpath = joinpath(Franklin.path(:site), rpath)
  # if the directory doesn't exist, create it
  outdir = dirname(outpath)
  isdir(outdir) || mkpath(outdir)
  # save the file and show it
  save(SVG(outpath), TikzCD(content))
  return "\\fig{/$(Franklin.unixify(rpath))}"
end


text:
"""
CLUSTERCOEFFS
-------------
    compute undirected clustering coefficients for a graph. clustercoeffs(A) computes a
    normalized, weighted clustering coefficients from a graph represented by a symmetric
    adjacency matrix A. clustercoeffs(A,weighted,normalized), with weighted and normalized
    boolean values indicating whether the computation has to be weighted and/or normalized.

Functions
---------
- cc = clustercoeffs(A::MatrixNetwork,weighted::Bool,normalized::Bool)
- cc = clustercoeffs{T}(A::SparseMatrixCSC{T,Int64},weighted::Bool,normalized::Bool)\n
If weighted and normalized are not specified, they are understood as true

Example
-------
~~~
A = load_matrix_network("clique-10")
cc = clustercoeffs(MatrixNetwork(A))
~~~
"""

function clustercoeffs(A::MatrixNetwork)
    return clustercoeffs(A, true, true);
end

function clustercoeffs(A::MatrixNetwork,weighted::Bool,normalized::Bool)
    donorm = true
    usew = true
    if !normalized
        donorm = false
    end
    if !weighted
        usew = false
    end

    if is_undirected(A) == false
        error("Only undirected (symmetric) inputs are allowed")
    end

    (rp,ci,ai) = (A.rp,A.ci,A.vals)
    if length(find(ai.<0)) != 0
        error("only positive edge weights allowed")
    end
    return clustercoeffs_phase2(donorm,rp,ci,ai,usew)
end

function clustercoeffs_phase2{T}(donorm::Bool,rp::Vector{Int64},ci::Vector{Int64},
                                          ai::Vector{T},usew::Bool)
    n = length(rp) - 1
    cc = zeros(Float64,n)
    # ind = falses(n,1)
    ind = zeros(Bool,n)
    cache=zeros(Float64,n)
    ew = 1
    ew2 = 1
    for v = 1:n
        for rpi = rp[v]:rp[v+1]-1
            w = ci[rpi]
            if usew
                ew = ai[rpi]
            end
            if v != w
                ind[w] = 1
                cache[w] = ew^(1/3)
            end
        end
        curcc = 0
        d = rp[v+1]-rp[v]
        # run two steps of bfs to try and find triangles.
        for rpi = rp[v]:rp[v+1]-1
            w = ci[rpi]
            if v == w
                d = d-1
                continue
            end #discount self-loop
            for rpi2 = rp[w]:rp[w+1]-1
                x = ci[rpi2]
                if x == w
                    continue
                end
                if ind[x]
                    if usew
                        ew = ai[rpi]
                        ew2 = ai[rpi2]
                    end
                    curcc = curcc+ew^(1/3)*ew2^(1/3)*cache[x]
                end
            end
        end
        if donorm && d>1
            cc[v] = curcc/(d*(d-1))
        elseif d>1
            cc[v] = curcc
        end

        for rpi = rp[v]:rp[v+1]-1
            w = ci[rpi]
            ind[w] = 0
        end # reset indicator
    end
    return cc
end

############################
### Additional functions ###
############################
# sparse matrices:
function clustercoeffs{T}(A::SparseMatrixCSC{T,Int64},weighted::Bool,normalized::Bool)
    donorm = true
    usew = true
    if !normalized
        donorm = false
    end
    if !weighted
        usew = false
    end

    At = A'
    if !(At == A)
        error("Only undirected (symmetric) inputs are allowed")
    end

    (rp,ci,ai) = (At.colptr,At.rowval,At.nzval);
    if length(find(ai.<0)) != 0
        error("only positive edge weights allowed")
    end
    return clustercoeffs_phase2(donorm,rp,ci,ai,usew)
end

function clustercoeffs{T}(A::SparseMatrixCSC{T,Int64})
    return clustercoeffs(A, true, true);
end

## Triplet Format:
function clustercoeffs(ei::Vector{Int64},ej::Vector{Int64})
    return clustercoeffs(MatrixNetwork(ei,ej))
end

function clustercoeffs(ei::Vector{Int64},ej::Vector{Int64},weighted::Bool,normalized::Bool)
    return clustercoeffs(MatrixNetwork(ei,ej),weighted,normalized)
end

## CSR sparse matrices - basically just like type MatrixNetwork
function clustercoeffs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64)
    return clustercoeffs(MatrixNetwork(n,rp,ci,vals))
end

function clustercoeffs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64,weighted::Bool,normalized::Bool)
    return clustercoeffs(MatrixNetwork(n,rp,ci,vals),weighted,normalized)
end


text:
using Documenter
using LoweredCodeUtils

makedocs(
    sitename = "LoweredCodeUtils",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [LoweredCodeUtils],
    pages = ["Home" => "index.md", "signatures.md", "edges.md", "api.md"]
)

deploydocs(
    repo = "github.com/JuliaDebug/LoweredCodeUtils.jl.git",
    push_preview = true
)


text:
# Set of allowed Julia operators (as given by MOI)
const OPS = [
    :+,       # n-ary
    :-,     # unary, binary
    :*,   # n-ary
    :/,       # unary
    :^, # binary (or INTPOWER)
    :sqrt,   # unary
    :exp,     # unary
    :log,     # unary
    :abs,     # unary
]

const UNARY_OPS_LOOKUP = (
    exp = SCIPcreateExprExp,
    log = SCIPcreateExprLog,
    abs = SCIPcreateExprAbs,
    cos = SCIPcreateExprCos,
    sin = SCIPcreateExprSin,
)

"""
Extract operators from Julia expr recursively and convert to SCIP expressions.
Returns the SCIP expression pointer, whether the expression is a pure value (without variables).
"""
function push_expr!(nonlin::NonlinExpr, scip::Ptr{SCIP_}, vars::Dict{VarRef, Ref{Ptr{SCIP_VAR}}}, expr::Expr)
    # Storage for SCIP_EXPR*
    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)
    num_children = length(expr.args) - 1
    pure_value = false

    if Meta.isexpr(expr, :comparison) # range constraint
        # args: [lhs, <=, mid, <=, rhs], lhs and rhs constant
        @assert length(expr.args) == 5
        @assert expr.args[2][1] == expr.args[4][1] == :<=

        # just call on middle expression, bounds are handled outside
        return push_expr!(nonlin, scip, vars, expr.args[3])

    elseif Meta.isexpr(expr, :call) # operator
        op = expr.args[1]
        if op in (:(==), :<=, :>=)
            # args: [op, lhs, rhs], rhs constant
            @assert length(expr.args) == 3

            # just call on lhs expression, bounds are handled outside
            return push_expr!(nonlin, scip, vars, expr.args[2])

        elseif op == :^
            # Special case: power with constant exponent.
            # The Julia expression considers the base and exponent to be
            # subexpressions. SCIP does in principle support constant
            # expressions, but in the case of SCIP_EXPR_REALPOWER, the exponent
            # value is stored directly, as the second child.
            @assert num_children == 2

            # Base (first child) is proper sub-expression.
            base, pure_value = push_expr!(nonlin, scip, vars, expr.args[2])

            # Exponent (second child) is stored as value.
            @assert isa(expr.args[3], Number)
            if !pure_value
                exponent = Cdouble(expr.args[3])
                @SCIP_CALL SCIPcreateExprPow(scip, expr__, base[], exponent, C_NULL, C_NULL)
            end
        elseif op in (:-, :+)
            @assert num_children >= 1

            if op === :- && num_children == 1
                # Special case: unary version of minus.
                # Then, insert the actual subexpression:
                right, pure_value = push_expr!(nonlin, scip, vars, expr.args[2])
                if !pure_value
                    # Finally, add the (binary) minus:
                    @SCIP_CALL SCIPcreateExprSum(
                        scip, expr__, Cint(1), [right[]],
                        [-1.0], 0.0, C_NULL, C_NULL,
                    )
                end
            else
                coef_mul = if op === :-
                    -1.0
                else
                    1.0
                end
                subexprs_pairs = [push_expr!(nonlin, scip, vars, expr.args[i + 1]) for i in 1:num_children]
                if all(pair -> pair[2], subexprs_pairs)
                    pure_value = true
                else
                    # replace pure-value expression by its evaluated value, since no expression has been evaluated yet
                    subexprs = map(eachindex(subexprs_pairs)) do idx
                        (subexpr, is_pure_value) = subexprs_pairs[idx]
                        if is_pure_value
                            expr_val = Meta.eval(expr.args[idx + 1])
                            expr_ptr, _ = push_expr!(nonlin, scip, vars, expr_val)
                            expr_ptr[]
                        else
                            subexpr[]
                        end
                    end
                    coefs = fill(coef_mul, num_children)
                    coefs[1] = 1.0
                    @SCIP_CALL SCIPcreateExprSum(scip, expr__, Cint(num_children), subexprs, coefs, 0.0, C_NULL, C_NULL)
                    pure_value = false
                end
            end
        elseif op == :*
            @assert num_children >= 1
            subexprs_pairs = [push_expr!(nonlin, scip, vars, expr.args[i + 1]) for i in 1:num_children]
            if all(pair -> pair[2], subexprs_pairs)
                pure_value = true
            else
                subexprs = map(eachindex(subexprs_pairs)) do idx
                    (subexpr, is_pure_value) = subexprs_pairs[idx]
                    if is_pure_value
                        expr_val = Meta.eval(expr.args[idx + 1])
                        expr_ptr, _ = push_expr!(nonlin, scip, vars, expr_val)
                        expr_ptr[]
                    else
                        subexpr[]
                    end
                end
                @SCIP_CALL SCIPcreateExprProduct(scip, expr__, num_children, subexprs, 1.0, C_NULL, C_NULL)
                pure_value = false
            end
        elseif op in keys(UNARY_OPS_LOOKUP)
            # Unary operators
            @assert num_children == 1
            # Insert child expression:
            child, pure_value = push_expr!(nonlin, scip, vars, expr.args[2])
            if !pure_value
                @SCIP_CALL UNARY_OPS_LOOKUP[op](scip, expr__, child[], C_NULL, C_NULL)
            end
        elseif op == :sqrt
            child, pure_value = push_expr!(nonlin, scip, vars, expr.args[2])
            if !pure_value
                @SCIP_CALL SCIPcreateExprPow(scip, expr__, child[], 0.5, C_NULL, C_NULL)
            end
        elseif op == :/
            @assert num_children == 2
            # Transform L / R => L * R^-1.0
            # Create left and right subexpression.
            left, pure_left = push_expr!(nonlin, scip, vars, expr.args[2])
            inverse_right_expr = :($(expr.args[3])^-1.0)
            right, pure_right = push_expr!(nonlin, scip, vars, inverse_right_expr)
            if pure_left && pure_right
                pure_value = true
            else
                if pure_left && expr.args[2] isa Expr
                    val = Meta.eval(expr.args[2])
                    left, _ = push_expr!(nonlin, scip, vars, val)
                end
                if pure_right
                    val = Meta.eval(inverse_right_expr)
                    right, _ = push_expr!(nonlin, scip, vars, val)
                end
                # TODO evaluate left or right
                @SCIP_CALL SCIPcreateExprProduct(scip, expr__, num_children, [left[], right[]], 1.0, C_NULL, C_NULL)
            end
        else
            # attempt computation of pure value
            try
                Meta.eval(expr)
            catch e
                error("Operator $op (in $expr) not supported by SCIP.jl!")
            end
            # if reaching this point, the expression is a pure value
            pure_value = true
        end

    elseif Meta.isexpr(expr, :ref) # variable
        # It should look like this:
        # :(x[MathOptInterface.VariableIndex(1)])
        @assert expr.args[1] == :x
        @assert num_children == 1
        vi = expr.args[2] # MOI.VariableIndex
        vr = VarRef(vi.value)
        v = vars[vr][]
        @SCIP_CALL SCIPcreateExprVar(scip, expr__, v, C_NULL, C_NULL)
        pure_value = false
    else
        error("Expression $expr not supported by SCIP.jl!")
    end

    if !pure_value
        # Return this expression to be referenced by parent expressions
        push!(nonlin.exprs, expr__)
    end
    return expr__, pure_value
end

function push_expr!(nonlin::NonlinExpr, scip::Ptr{SCIP_}, vars::Dict{VarRef, Ref{Ptr{SCIP_VAR}}}, expr::Number)
    # Storage for SCIP_EXPR*
    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)

    value = Cdouble(expr)
    @SCIP_CALL SCIPcreateExprValue(scip, expr__, value, C_NULL, C_NULL)

    @assert SCIPisExprValue(scip, expr__[]) == TRUE
    push!(nonlin.exprs, expr__)
    pure_value = true
    return expr__, pure_value
end

"""
Add nonlinear constraint to problem, return cons ref.

    lhs ≤ expression ≤ rhs

# Arguments
- `expr::Expr`: Julia expression of nonlinear function, as given by MOI.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

"""
function add_nonlinear_constraint(scipd::SCIPData, expr::Expr, lhs::Float64, rhs::Float64)
    nonlin = NonlinExpr()

    # convert expression recursively, extract root and variable pointers
    root_expr, pure_value = push_expr!(nonlin, scipd.scip[], scipd.vars, expr)
    if pure_value
        # TODO
        error("Constraint $expr with pure value")
    end
    # create and add cons_nonlinear
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicNonlinear(scipd, cons__, "", root_expr[], lhs, rhs)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])

    # register and return cons ref
    push!(scipd.nonlinear_storage, nonlin)
    return store_cons!(scipd, cons__)
end


text:
function get_coefs(margin::UD, n::Int)
    aₖ = zeros(Float64, n + 1)
    m = 2n
    tₛ, wₛ = gausshermite(m)
    tₛ    .= tₛ * sqrt2
    Xₛ = normal_to_margin(margin, tₛ)

    aₖ = [sum(wₛ .* _h(tₛ, k) .* Xₛ) for k in 0:n]

    return invsqrtpi * aₖ ./ factorial.(0:n)
end
get_coefs(margins::UD, n::Real) = get_coefs(margins, Int(n))

function Gn0d(n::Int, A::UnitRange{Int}, B::UnitRange{Int}, α::Vector{Float64}, β::Vector{Float64}, σAσB_inv::Float64)
    if n == 0
        return 0.0
    end
    M = length(A)
    N = length(B)
    accu = 0.0
    for r=1:M, s=1:N
        r11 = Hp(α[r+1], n-1) * Hp(β[s+1], n-1)
        r00 = Hp(α[r],   n-1) * Hp(β[s],   n-1)
        r01 = Hp(α[r],   n-1) * Hp(β[s+1], n-1)
        r10 = Hp(α[r+1], n-1) * Hp(β[s],   n-1)
        accu += A[r]*B[s] * (r11 + r00 - r01 - r10)
    end
    accu * σAσB_inv
end

function Gn0m(n::Int, A::UnitRange{Int}, α::Vector{Float64}, dB::UD, σAσB_inv::Float64)
    if n == 0
        return 0.0
    end
    M = length(A)
    accu = 0.0
    for r=1:M
        accu += A[r] * (Hp(α[r+1], n-1) - Hp(α[r], n-1))
    end
    m = n + 4
    t, w = gausshermite(m)
    t .= t * sqrt2
    X = normal_to_margin(dB, t)
    S = invsqrtpi * sum(w .* _h(t, n) .* X)
    return -σAσB_inv * accu * S
end

function solve_poly_pm_one(coef::Vector{<:Real})
    P = Polynomial(coef)
	dP = derivative(P)
    r = roots(x->P(x), x->dP(x), interval(-1, 1), Krawczyk, 1e-3)

    length(r) == 1 && return mid(r[1].interval)
    length(r) == 0 && return NaN
    length(r) >  1 && return [mid(rs.interval) for rs in r]
end
nearest_root(target::Real, roots::Vector{<:Real}) = roots[argmin(abs.(roots .- target))]

function _h(x::T, n::Int) where {T<:Real}
    if n == 0
        return one(T)
    elseif n == 1
        return x
    end

    Hkp1, Hk, Hkm1 = zero(T), x, one(T)
    for k in 2:n
        Hkp1 = x*Hk - (k-1) * Hkm1
        Hkm1, Hk = Hk, Hkp1
    end
    Hkp1
end
_h(X::Real, n::Real) = _h(X, Int(n))
_h(A::Array{<:Real, N}, n::Real) where N = _h.(A, Ref(n))

# We need to account for when x is ±∞ otherwise Julia will return NaN for 0×∞
function Hp(x::Real, n::Int)
    isinf(x) ? zero(x) : _h(x, n) * _normpdf(x)
end
Hp(x::Real, n::Real) = Hp(x, Int(n))


text:
# FIXME:
# Remove explicit names from these callbacks and add connectors which extract from NT
# modularse.  Have one cb that accumuluates
#

"""
```julia
using Lens
x = normal(0, 1)
@leval Loop => plotloss() rand(x, x >ₛ 0.0, 100; alg = SSMH)
"""
function plotloss(; title = "Loss vs Iteration", display_ = display, LTY::Type{TY} = Float64) where {TY}
  @show cb = Callbacks.plotscalar(; title = title, display_ = display_, LTX = Int, LTY = LTY)
  data -> cb((x = data.i, y = data.p))
end

"Scatter plot ω values with UnicodePlots"
function plotω(x::RandVar, y::RandVar, T = Float64)
  xωs = T[]
  yωs = T[]

  function innerplotω(data)
    color = :blue
    if isempty(xωs)
      color = :red
    end
    xω = data.ω.vals[x.id]
    yω = data.ω.vals[y.id]
    push!(xωs, xω)
    push!(yωs, yω)
    println(scatterplot(xωs, yωs, title="ω movement", color=color, xlim=[0, 1], ylim=[0, 1]))
  end
end

"Line Plot histogram of loss with UnicodePlots"
function plotscalar(key::Symbol, name, display_ = display, ::Type{T} = Float64) where T
  xs = T[]
  ys = Int[]
  maxseen = typemin(T)
  minseen = typemax(T)
  println("HELLO!!!!!")

  function innerplotrv(data)
    x_ = getfield(data, key)
    println("$name is:")
    display_(x_)
    push!(xs, x_)
    push!(ys, data.i)
    if !isempty(xs)
      println(UnicodePlots.lineplot(ys, xs, title="Time vs $name"))
    end
    if x_ > maxseen
      maxseen = x_
      printstyled("\nNew max at id $(data.i): $(x_)\n"; color = :light_blue)
    end
    if x_ < minseen
      minseen = x_
      printstyled("\nNew min at id $(data.i): $(data.p)\n"; color = :light_blue)
    end
  end
end

"Print acceptance statistics"
function printstats(data)
  printstyled("\nacceptance ratio: $(data.accepted/float(data.i))\n",
              "Last p: $(data.p)\n";
              color = :light_blue)
end

default_cbs_tpl(n) = (throttle(plotloss(), 0.1),
                      showprogress(n),
                      throttle(printstats, 1.0))

"Defautlt callbacks"
default_cbs(n) = runall([default_cbs_tpl(n)...])


text:
using CollectiveSpins

using PyCall
@pyimport matplotlib.pyplot as plt

function plot_geometry_xy(positions)
    for (i,p)=enumerate(positions)
        plt.text(p[1], p[2], "($i)", color="black", va="center", ha="center")
    end
    xmin = minimum([p[1] for p=positions])
    xmax = maximum([p[1] for p=positions])
    ymin = minimum([p[2] for p=positions])
    ymax = maximum([p[2] for p=positions])
    dx = xmax - xmin
    dy = ymax - ymin
    marginx = max(dx/5, 0.1)
    marginy = max(dy/5, 0.1)
    plt.xlim([xmin-marginx, xmax+marginx])
    plt.ylim([ymin-marginy, ymax+marginy])
end

d = 0.1
geomN = 5
positions = geometry.hexagonal(d; Nr=geomN)
positions = geometry.chain(d, 2*geomN+1)
plot_geometry_xy(positions)

plt.show()


text:
function show(io::IO, object::Node)
  if length(object.in) == 0
    if length(object.out) > 0
      for i in 1:length(object.out)
        if i == 1
          print(io, "\r\e[0m[\e[1mroot node\e[0m]-->[branch $(object.out[i])]")
        elseif i < length(object.out)
          print(io, "\n\e[0m$(repeat(" ", length("[root node]")))-->[branch $(object.out[i])]")
        else
          print(io, "\n\e[0m$(repeat(" ", length("[root node]")))-->[branch $(object.out[i])]")
        end
      end
    else
      print(io, "\r\e[0m[\e[1munattached node\e[0m]")
    end
  elseif length(object.in) == 1
    if length(object.out) == 0
      print(io, "\r\e[0m[branch $(object.in[1])]-->[\e[1mleaf node\e[0m]")
    elseif length(object.in) == 1
      for i in 1:length(object.out)
        if length(object.out) == 1
          print(io, "\r\e[0m[branch $(object.in[1])]-->[\e[1minternal node\e[0m]-->[branch $(object.out[i])]")
        elseif i == 1
          print(io, "\r\e[0m[branch $(object.in[1])]-->[\e[1minternal node\e[0m]-->[branch $(object.out[i])]")
        elseif i < length(object.out)
          print(io, "\n\e[0m$(repeat(" ", length("[branch $(object.in[1])]-->[internal node]")))-->[branch $(object.out[i])]")
        else
          print(io, "\n\e[0m$(repeat(" ", length("[branch $(object.in[1])]-->[internal node]")))-->[branch $(object.out[i])]")
        end
      end
    end
  end
end

function show(io::IO, object::Branch)
  print(io, "\r\e[0m[node $(object.source)]-->[\e[1m$(object.length) branch\e[0m]-->[node $(object.target)]")
end

function show(io::IO, object::Tree)
  print(io, "\r\e[0mPhylogenetic tree with $(length(object.nodes)) nodes and $(length(object.branches)) branches")
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 82467e1e-97e0-11eb-35bc-b70fe47f7296
using Revise, FourierTools, ImageCore, TestImages, PlutoUI, ImageShow

# ╔═╡ 8f0854f4-97e0-11eb-14b8-7958b0ab243f
begin
	img = Float32.(testimage("fabio_512_gray"))
	z = zeros(Float32, (768, 768))
	FourierTools.center_set!(z, img)
end

# ╔═╡ e43962ba-97e0-11eb-274f-fd6ed85b03ff
@bind ϕ Slider(-180:1:180)

# ╔═╡ dc4917a8-97e0-11eb-341d-25ec7207b748
md"ϕ=$ϕ"

# ╔═╡ a9e968a8-97e0-11eb-0486-3f75f68be3c7
Gray.(FourierTools.rotate(z, ϕ))

# ╔═╡ f1fecc74-97e1-11eb-203e-2949bafae27d
Gray.(FourierTools.rotate(img, ϕ))

# ╔═╡ Cell order:
# ╠═82467e1e-97e0-11eb-35bc-b70fe47f7296
# ╠═8f0854f4-97e0-11eb-14b8-7958b0ab243f
# ╠═dc4917a8-97e0-11eb-341d-25ec7207b748
# ╠═e43962ba-97e0-11eb-274f-fd6ed85b03ff
# ╠═a9e968a8-97e0-11eb-0486-3f75f68be3c7
# ╠═f1fecc74-97e1-11eb-203e-2949bafae27d


text:
#include("../dataset_reader.jl")
const dataset5_labels = Float64[0, 1, 2, 3, 4]
const dataset5_nfeatures = 13
const dataset5_filename = "original/processed/processed.cleveland.data"

function remove_null_data(data::Vector{Vector{String}})
    regex_number = r"\S*((\+|-)?([0-9]+(\.[0-9]*)?|\.[0-9]+))\S*"
    new_data = Vector{String}[]
    for line in data
        if all(occursin.(regex_number, line))
            push!(new_data, line)
        end
    end
    new_data
end

function get_numerical_dataset5(path::String="")
    data = get_data(joinpath(path, dataset5_filename); dlm=",", delete_first_line=false)
    data = remove_null_data(data)

    conversion_pattern = [Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64]
    convert_values(data, conversion_pattern)
end

function get_dataset5(path::String=""; shuffle::Bool=false, test_perc::Float64=0.2,
                      normalise::Bool=false, separate_test_data::Bool=true,
                      balance_classes::Bool=false)
    data = get_numerical_dataset5(path)
    separate_dataset(data, test_percentage=test_perc, shuffle=shuffle, normalise=normalise,
                     separate_test_data=separate_test_data, balance_classes=balance_classes)
end

function get_dataset5_setup(path::String=""; normalise::Bool=false, separate::Bool=false,
                            balance::Bool=false)
    dataset = get_dataset5(path, normalise=normalise, separate_test_data=separate, balance_classes=balance)
    if separate
        SeparatedDataset(dataset..., dataset5_labels, dataset5_nfeatures, 6)
    else
        Dataset(dataset..., dataset5_labels, dataset5_nfeatures, 6)
    end
end

#=
Data attributes:
#1 - Numerical (integer)
#2 - Categorical (0, 1)
#3 - Categorical (1, 2, 3, 4)
#4 - Numerical (fractional)
#5 - Numerical (fractional)
#6 - Categorical (0, 1)
#7 - Categorical (0, 1, 2)
#8 - Numerical (fractional)
#9 - Categorical (0, 1)
#10 - Numerical (fractional)
#11 - Categorical (1, 2, 3)
#12 - Categorical (0, 1, 2, 3)
#13 - Categorical (3, 6, 7)
#14 - Categorical (0, 1, 2, 3, 4)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
export BATResult
export batfit!
export get_marginalized_pdf
export get_histogram_pdf
export get_samples_inds_by_chain_id

struct HistogramModelLikelihood{H<:Histogram{<:Real, 1}, F<:FitFunction, T <: Real} <: BAT.AbstractDensity
    h::H
    f::F
    weights::Vector{T}
    midpoints::Vector{T}
    bin_volumes::Vector{T}
    logabsgamma::Vector{T}

    function HistogramModelLikelihood(h::H, f::FitFunction) where {H <: Histogram{<:Real, 1}}
        T = get_pricision_type(f)
        new{H, typeof(f), T}( h, f, h.weights,
                        StatsBase.midpoints(h.edges[1]),
                        T[StatsBase.binvolume(h, i) for i in eachindex(h.weights)],
                        T[logabsgamma(w + 1)[1] for w in h.weights]  )
    end
end

function log_pdf_poisson(λ::T, k::U, logabsgamma::T) where {T<:Real,U<:Real}
    R = float(promote_type(T,U))
    if λ >= 0 && k >= 0 && isinteger(k)
        result = (iszero(k) ? R(k) : R(log(λ) * k)) - λ - logabsgamma
        R(result)
    else
        R(-Inf)
    end
end

function BAT.density_logval( l::HistogramModelLikelihood{H, F, T}, pars) where {H, F, T}
    log_likelihood::Float64 = 0
    @inbounds for i in eachindex(l.weights)
        expected_counts::Float64 = l.bin_volumes[i] * l.f.model(l.midpoints[i], pars)
        if isnan(expected_counts) || expected_counts < 0
            expected_counts = T(Inf)
        end
        log_likelihood += log_pdf_poisson(expected_counts, l.weights[i], l.logabsgamma[i])
        if isnan(log_likelihood)
            println(pars)
            println(l.f.model(l.midpoints[i], pars))
            error("..")
        end
    end
    return log_likelihood
end

mutable struct BATResult{B, P}
    bat_result::B
    prior::P

    BATResult(b::B) where B = new{B, Missing}(b, missing)
    BATResult(b::B, p::P) where {B, P} = new{B, P}(b, p)
end

function batfit!(f::FitFunction{T, ND, NP}, h::Histogram{<:Real, 1};
                rng_seed = Philox4x((123, 456)),
                prior = missing,
                nsamples::Int = 10^5,
                nchains::Int = 4,
                pretunesamples::Int = length(f.parameter_bounds) * 1000,
                max_ncycles::Int = 30,
                tuning_r = 0.5,
                BGConvergenceThreshold::Real = sqrt(length(f.parameter_bounds))) where {T, ND, NP}
    first_bin = !isinf(first(f.fitranges[1])) ? StatsBase.binindex(h, first(f.fitranges[1])) : 1
    last_bin  = !isinf(last( f.fitranges[1])) ? StatsBase.binindex(h, last( f.fitranges[1])) : length(h.weights)
    if first_bin < 1 first_bin = 1 end
    if (last_bin > length(h.weights)) last_bin = length(h.weights) end
    h_sub = Histogram(h.edges[1][first_bin:last_bin], h.weights[first_bin:last_bin-1])
    h_bat = HistogramModelLikelihood(h_sub, f)

    if ismissing(prior)
        prior_dists = if !ismissing(f.parameter_bounds)
            Tuple(map(b -> Uniform(b.left, b.right), f.parameter_bounds))
        else
            error("`batfit!` needs a `FitFunction` with parameter_bounds.")
        end

        prior = BAT.NamedTupleDist( (;zip( f.parameter_names, prior_dists)...))
        prior = BAT.DistributionDensity(prior, bounds_type = BAT.reflective_bounds)
    end

    posterior = BAT.PosteriorDensity(h_bat, prior)

    algorithm = BAT.MetropolisHastings()

    tuning = BAT.AdaptiveMetropolisTuning(
        r = tuning_r,
        λ = 0.5,
        α = 0.15..0.35,
        β = 1.5,
        c = 1e-4..1e2
    )

    convergence = BAT.BrooksGelmanConvergence(
        threshold = T(BGConvergenceThreshold),
        corrected = false
    )

    init = BAT.MCMCInitStrategy(
        init_tries_per_chain = 8..128,
        max_nsamples_init = pretunesamples,
        max_nsteps_init = pretunesamples * 10,
        max_time_init = Inf
    )

    burnin = BAT.MCMCBurninStrategy(
        max_nsamples_per_cycle = pretunesamples,
        max_nsteps_per_cycle = pretunesamples * 10,
        max_time_per_cycle = Inf,
        max_ncycles = max_ncycles
    )

    bat_result = BATResult(BAT.bat_sample( rng_seed,
        posterior, (nsamples, nchains), algorithm,
        max_nsteps = 10 * nsamples,
        max_time = Inf,
        tuning = tuning,
        init = init,
        burnin = burnin,
        convergence = convergence,
        strict = false,
        filter = true),
        prior
    )

    f.backend_result = (bat_result, )
    try
        _set_fitted_parameters!(f, mode(RadiationSpectra.get_bat_result(f))[])
    catch err
        @warn "Determination of mode failed"
    end
    f
end

### Utilties

function rs_write(fn::AbstractString, b::BATResult)
    BAT.bat_write(fn, unshape(b))
    ser_bfn = begin
        bn = basename(fn)
        if occursin(".", bn)
            bn[1:first(findlast(".", bn))-1] * ".ser"
        else
            bn * ".ser"
        end
    end
    ser_fn = joinpath(dirname(fn), ser_bfn)
    serialize(ser_fn, b.prior)
    nothing
end
function rs_read(fn::AbstractString)
    ser_bfn = begin
        bn = basename(fn)
        if occursin(".", bn)
            bn[1:first(findlast(".", bn))-1] * ".ser"
        else
            bn * ".ser"
        end
    end
    ser_fn = joinpath(dirname(fn), ser_bfn)
    if !isfile(ser_fn) error("Prior was not serialized to file.") end
    prior = deserialize(ser_fn)
    s = varshape(prior)
    bat_samples = s.(BAT.bat_read(fn))
    return BATResult((result = bat_samples,), prior)
end

get_fit_backend_result(fr::Tuple) = fr

get_bat_result(ff::FitFunction)::BAT.DensitySampleVector =
    get_fit_backend_result(ff.backend_result)[1].bat_result.result

function BAT.NamedTupleShape(f::FitFunction)
    ts = Tuple(map(b -> ScalarShape{Real}(), f.parameter_names))
    return BAT.NamedTupleShape(NamedTuple{Tuple(f.parameter_names)}(ts))
end

function _set_fitted_parameters!(ff::FitFunction{T}, fitted_parameters::ShapedAsNT)::Nothing where {T <: AbstractFloat}
    _set_fitted_parameters!(ff, vcat(values(fitted_parameters[])...))
end

get_samples_inds_by_chain_id(ff::FitFunction, ichain::Int = 1) =
    get_samples_inds_by_chain_id(get_bat_result(ff), ichain)

function get_samples_inds_by_chain_id(bat_result, ichain::Int)
    chainids = Int[]
    sampleids = bat_result.info
    for s in sampleids
        if !(s.chainid in chainids) push!(chainids, s.chainid) end
    end
    chainid = chainids[ichain]
    inds = findall( sid -> sid.chainid == chainid, sampleids )
    return inds
end

function _get_standard_deviations(f::FitFunction{T}, fr::Tuple) where {T}
    uncertainties = zeros(T, length(f.fitted_parameters))
    for i in eachindex(f.fitted_parameters)
        h = get_marginalized_pdf(BAT.unshaped.(get_bat_result(f)), i, nbins = 100)
        d = fit(Normal, h.edges[1], h.weights)
        uncertainties[i] = d.σ
    end
    return uncertainties
end

unshape(samples::BAT.DensitySampleVector) =
    samples.v isa ShapedAsNTArray ? ValueShapes.unshaped.(samples) : samples
unshape(b::BATResult) =
    unshape(b.bat_result[1])

function get_marginalized_pdf(samples::BAT.DensitySampleVector, pars::NTuple{N, Int}; nbins = 100, sample_range::AbstractVector{Int} = Int[]) where {N}
    if isempty(sample_range) sample_range = 1:length(samples) end
    unshaped_samples = unshape(samples)
    pdf = fit(Histogram, Tuple([flatview(unshaped_samples.v)[ipar, sample_range] for ipar in pars]),
                FrequencyWeights(unshaped_samples.weight[sample_range]), closed=:left, nbins=nbins)
    return normalize(pdf)
end

get_marginalized_pdf(br::BATResult, args...; kwargs...) =
    get_marginalized_pdf(br.bat_result.result, args...; kwargs...)

get_marginalized_pdf(f::FitFunction, args...; nbins = 100, sample_range::AbstractVector{Int} = Int[]) =
    get_marginalized_pdf(get_bat_result(f), args..., nbins = nbins, sample_range = sample_range)

get_marginalized_pdf(samples::BAT.DensitySampleVector, ipar::Int; nbins = 100, sample_range::AbstractVector{Int} = Int[]) =
    get_marginalized_pdf(samples, (ipar,); nbins = nbins, sample_range = sample_range)

get_marginalized_pdf(samples::BAT.DensitySampleVector, var::ValueShapes.ValueAccessor; nbins = 100, sample_range::AbstractVector{Int} = Int[]) =
    get_marginalized_pdf(samples, Tuple([var.offset + i for i in 1:var.len ]); nbins = nbins, sample_range = sample_range)

function get_histogram_pdf(samples::BAT.DensitySampleVector; nbins = 10, sample_range::AbstractVector{Int} = Int[]) where {N}
    unshaped_samples = unshape(samples)
    n_params::Int = size(flatview(unshaped_samples.v), 1)
    get_marginalized_pdf(unshaped_samples, Tuple(1:n_params), nbins = nbins, sample_range = sample_range)
end
get_histogram_pdf(f::FitFunction, args...; kwargs...) =
    get_histogram_pdf(get_bat_result(f), args...; kwargs...)

function calculate_localmode(h::StatsBase.Histogram{<:Real, N}) where {N}
    cartesian_inds = CartesianIndices(h.weights)[findmax(h.weights)[2]]
    return [StatsBase.midpoints(h.edges[i])[cartesian_inds[i]] for i in 1:N]
end

function central_intervals(marginalized_pdf::Histogram{<:Real, 1}, intervals = BAT.standard_confidence_vals)
    sort!(intervals)
    splitted_histograms = reverse(BAT.split_central(marginalized_pdf, intervals)[1])
    r = Interval[]
    for ipar in 1:length(intervals)
        ifirst = findfirst(w -> w > 0, splitted_histograms[ipar].weights)
        ilast  = findlast( w -> w > 0, splitted_histograms[ipar].weights)
        int = splitted_histograms[ipar].edges[1][ifirst]..splitted_histograms[ipar].edges[1][ilast+1]
        push!(r, int)
    end
    return r
end
function smallest_intervals(marginalized_pdf::Histogram{<:Real, 1}, intervals = BAT.standard_confidence_vals)
    sort!(intervals)
    splitted_histograms = reverse(BAT.split_smallest(marginalized_pdf, intervals)[1])
    r = Interval[]
    for ipar in 1:length(intervals)
        ifirst = findfirst(w -> w > 0, splitted_histograms[ipar].weights)
        ilast  = findlast( w -> w > 0, splitted_histograms[ipar].weights)
        int = splitted_histograms[ipar].edges[1][ifirst]..splitted_histograms[ipar].edges[1][ilast+1]
        push!(r, int)
    end
    return r
end

@userplot Plot_Marginalized_PDF
@recipe function f(h::Plot_Marginalized_PDF; interval_type = :smallest)
    if (length(h.args) < 1) || !(typeof(h.args[1]) <: Histogram{<:Real, 1})
        error("Wrong arguments. Got: $(typeof(h.args))")
    end
    hmpdf = h.args[1]
    intervals = BAT.standard_confidence_vals

    @series begin
        seriestype := :step
        label --> "Marginalized PDF"
        hmpdf
    end

    h_ints = if interval_type == :smallest
        reverse(BAT.split_smallest(hmpdf, BAT.standard_confidence_vals)[1])
    elseif interval_type == :central
        reverse(BAT.split_central(hmpdf, BAT.standard_confidence_vals)[1])
    else
        error("`interval_type` must be either `:smallest` or `:central`")
    end
    for i in length(intervals):-1:1
        @series begin
            label := """$((interval_type == :smallest ? "SI" : "CI")): $(BAT.standard_confidence_vals[i])"""
            fillcolor := BAT.standard_colors[i]
            h_ints[i]
        end
    end
end


text:
immutable Givens{T} <: AbstractMatrix{T}
    size::Int
    i1::Int
    i2::Int
    c::T
    s::T
    r::T
end
type Rotation{T}
    rotations::Vector{T}
end
typealias AbstractRotation Union(Givens, Rotation)

function givensAlgorithm{T<:FloatingPoint}(f::T, g::T)
    zeropar = zero(T)
    onepar = one(T)
    twopar = 2one(T)

    safmin = realmin(T)
    epspar = eps(T)
    safmn2 = twopar^itrunc(log(safmin/epspar)/log(twopar)/twopar)
    safmx2 = onepar/safmn2

    if g == 0
        cs = onepar
        sn = zeropar
        r = f
    elseif f == 0
        cs = zeropar
        sn = onepar
        r = g
    else
        f1 = f
        g1 = g
        scalepar = max(abs(f1), abs(g1))
        if scalepar >= safmx2
            count = 0
            while true
                count += 1
                f1 *= safmn2
                g1 *= safmn2
                scalepar = max(abs(f1), abs(g1))
                if scalepar < safmx2 break end
            end
            r = hypot(f1, g1)
            cs = f1/r
            sn = g1/r
            for i = 1:count
                r *= safmx2
            end
        elseif scalepar <= safmn2
            count = 0
            while true
                count += 1
                f1 *= safmx2
                g1 *= safmx2
                scalepar = max(abs(f1), abs(g1))
                if scalepar > safmn2 break end
            end
            r = hypot(f1, g1)
            cs = f1/r
            sn = g1/r
            for i = 1:count
                r *= safmn2
            end
        else
            r = hypot(f1, g1)
            cs = f1/r
            sn = g1/r
        end
        if abs(f) > abs(g) && cs < 0
            cs = -cs
            sn = -sn
            r = -r
        end
    end
    return cs, sn, r
end

function givensAlgorithm{T<:FloatingPoint}(f::Complex{T}, g::Complex{T})
    twopar, onepar, zeropar = 2one(T), one(T), zero(T)
    czero = zero(Complex{T})

    abs1(ff) = max(abs(real(ff)), abs(imag(ff)))
    safmin = realmin(T)
    epspar = eps(T)
    safmn2 = twopar^itrunc(log(safmin/epspar)/log(twopar)/twopar)
    safmx2 = onepar/safmn2
    scalepar = max(abs1(f), abs1(g))
    fs = f
    gs = g
    count = 0
    if scalepar >= safmx2
        while true
            count += 1
            fs *= safmn2
            gs *= safmn2
            scalepar *= safmn2
            if scalepar < safmx2 break end
        end
    elseif scalepar <= safmn2
        if g == 0
            cs = onepar
            sn = czero
            r = f
            return cs, sn, r
        end
        while true
            count -= 1
            fs *= safmx2
            gs *= safmx2
            scalepar *= safmx2
            if scalepar > safmn2 break end
        end
    end
    f2 = abs2(fs)
    g2 = abs2(gs)
    if f2 <= max(g2, onepar)*safmin

     # This is a rare case: F is very small.

        if f == 0
            cs = zero
            r = hypot(real(g), imag(g))
        # do complex/real division explicitly with two real divisions
            d = hypot(real(gs), imag(gs))
            sn = complex(real(gs)/d, -imag(gs)/d)
            return cs, sn, r
        end
        f2s = hypot(real(fs), imag(fs))
     # g2 and g2s are accurate
     # g2 is at least safmin, and g2s is at least safmn2
        g2s = sqrt(g2)
     # error in cs from underflow in f2s is at most
     # unfl / safmn2 .lt. sqrt(unfl*eps) .lt. eps
     # if max(g2,one)=g2, then f2 .lt. g2*safmin,
     # and so cs .lt. sqrt(safmin)
     # if max(g2,one)=one, then f2 .lt. safmin
     # and so cs .lt. sqrt(safmin)/safmn2 = sqrt(eps)
     # therefore, cs = f2s/g2s / sqrt( 1 + (f2s/g2s)**2 ) = f2s/g2s
        cs = f2s/g2s
     # make sure abs(ff) = 1
     # do complex/real division explicitly with 2 real divisions
        if abs1(f) > 1
            d = hypot(real(f), imag(f))
            ff = complex(real(f)/d, imag(f)/d)
        else
            dr = safmx2*real(f)
            di = safmx2*imag(f)
            d = hypot(dr, di)
            ff = complex(dr/d, di/d)
        end
        sn = ff*complex(real(gs)/g2s, -imag(gs)/g2s)
        r = cs*f + sn*g
    else

     # This is the most common case.
     # Neither F2 nor F2/G2 are less than SAFMIN
     # F2S cannot overflow, and it is accurate

        f2s = sqrt(onepar + g2/f2)
     # do the f2s(real)*fs(complex) multiply with two real multiplies
        r = complex(f2s*real(fs), f2s*imag(fs))
        cs = onepar/f2s
        d = f2 + g2
     # do complex/real division explicitly with two real divisions
        sn = complex(real(r)/d, imag(r)/d)
        sn *= conj(gs)
        if count != 0
            if count > 0
                for i = 1:count
                    r *= safmx2
                end
            else
                for i = 1:-count
                    r *= safmn2
                end
            end
        end
    end
    return cs, sn, r
end

function givens{T}(f::T, g::T, i1::Integer, i2::Integer, size::Integer)
    i2 <= size || error("indices cannot be larger than size Givens rotation matrix")
    i1 < i2 || error("second index must be larger than the first")
    h = givensAlgorithm(f, g)
    Givens(size, i1, i2, convert(T, h[1]), convert(T, h[2]), convert(T, h[3]))
end

function givens{T}(A::AbstractMatrix{T}, i1::Integer, i2::Integer, col::Integer)
    i1 < i2 || error("second index must be larger than the first")
    h = givensAlgorithm(A[i1,col], A[i2,col])
    Givens(size(A, 1), i1, i2, convert(T, h[1]), convert(T, h[2]), convert(T, h[3]))
end

*{T}(G1::Givens{T}, G2::Givens{T}) = Rotation(push!(push!(Givens{T}[], G2), G1))
*(G::Givens, B::BitArray{2}) = error("method not defined")
*{TBf,TBi}(G::Givens, B::SparseMatrixCSC{TBf,TBi}) = error("method not defined")
*(R::AbstractRotation, A::AbstractMatrix) = A_mul_B!(R, copy(A))

A_mul_Bc(A::AbstractMatrix, R::Rotation) = A_mul_Bc!(copy(A), R)

getindex(G::Givens, i::Integer, j::Integer) = i == j ? (i == G.i1 || i == G.i2 ? G.c : one(G.c)) : (i == G.i1 && j == G.i2 ? G.s : (i == G.i2 && j == G.i1 ? -G.s : zero(G.s)))
size(G::Givens) = (G.size, G.size)
size(G::Givens, i::Integer) = 1 <= i <= 2 ? G.size : 1
A_mul_B!(G1::Givens, G2::Givens) = error("Operation not supported. Consider *")
function A_mul_B!(G::Givens, A::AbstractMatrix)
    m, n = size(A)
    for i = 1:n
        tmp = G.c*A[G.i1,i] + G.s*A[G.i2,i]
        A[G.i2,i] = G.c*A[G.i2,i] - conj(G.s)*A[G.i1,i]
        A[G.i1,i] = tmp
    end
    return A
end
function A_mul_Bc!(A::AbstractMatrix, G::Givens)
    m, n = size(A)
    for i = 1:m
        tmp = G.c*A[i,G.i1] + conj(G.s)*A[i,G.i2]
        A[i,G.i2] = G.c*A[i,G.i2] - G.s*A[i,G.i1]
        A[i,G.i1] = tmp
    end
    return A
end
function A_mul_B!(G::Givens, R::Rotation)
    push!(R.rotations, G)
    return R
end
function A_mul_B!(R::Rotation, A::AbstractMatrix)
    for i = 1:length(R.rotations)
        A_mul_B!(R.rotations[i], A)
    end
    return A
end
function A_mul_Bc!(A::AbstractMatrix, R::Rotation)
    for i = 1:length(R.rotations)
        A_mul_Bc!(A, R.rotations[i])
    end
    return A
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
@testset "Automatic differentiation rules" begin
    @testset "QDHT" begin
        rng = MersenneTwister(68)
        N = 64
        R = 10.0
        q = Hankel.QDHT{1,2}(R, N)
        @testset "rrule" begin
            q_rev, back = rrule(Hankel.QDHT{1,2}, R, N)
            @test typeof(q_rev) === typeof(q)
            for p in propertynames(q)
                @test getproperty(q, p) == getproperty(q_rev, p)
            end
            ∂QDHT, ∂R, ∂N = back(true)
            @test ∂QDHT === NoTangent()
            @test ∂R isa AbstractZero
            @test ∂N isa AbstractZero
        end
        @testset "frule" begin
            ΔQDHT, ΔR, ΔN = NoTangent(), true, true
            q_fwd, ∂q = frule((ΔQDHT, ΔR, ΔN), Hankel.QDHT{1,2}, R, N)
            @test typeof(q_fwd) === typeof(q)
            for p in propertynames(q)
                @test getproperty(q, p) == getproperty(q_fwd, p)
            end
            @test ∂q isa AbstractZero
        end
    end

    @testset "$f(::QDHT, ::Array{<:Real})" for f in (*, \)
        rng = MersenneTwister(86)
        N, M, K = 64, 5, 10
        @testset "Vector" begin
            q = Hankel.QDHT{1,2}(10, N)
            fr = randn(rng, N)
            test_rrule(f, q ⊢ NoTangent(), fr)
            test_frule(f, q ⊢ NoTangent(), fr)
        end

        @testset "Matrix" begin
            @testset for dim in (1, 2)
                q = Hankel.QDHT{1,2}(10, N; dim = dim)
                s = dim == 1 ? (N, M) : (M, N)
                fr = randn(rng, s)
                test_rrule(f, q ⊢ NoTangent(), fr)
                test_frule(f, q ⊢ NoTangent(), fr)
            end
        end

        @testset "Array{<:Real,3}" begin
            q = Hankel.QDHT{1,2}(10, N; dim = 2)
            s = (M, N, K)
            fr = randn(rng, s)
            test_rrule(f, q ⊢ NoTangent(), fr)
            test_frule(f, q ⊢ NoTangent(), fr)
        end
    end

    @testset "$f(::Array{<:Real}, ::QDHT, ::Array{<:Real})" for f in (mul!, ldiv!)
        rng = MersenneTwister(86)
        N, M, K = 64, 5, 10
        @testset "Vector" begin
            q = Hankel.QDHT{1,2}(10, N)
            fr, fk = ntuple(_ -> randn(rng, N), 2)
            test_frule(f, fk, q ⊢ NoTangent(), fr)
        end

        @testset "Matrix" begin
            @testset for dim in (1, 2)
                q = Hankel.QDHT{1,2}(10, N; dim = dim)
                s = dim == 1 ? (N, M) : (M, N)
                fr, fk = ntuple(_ -> randn(rng, s), 2)
                test_frule(f, fk, q ⊢ NoTangent(), fr)
            end
        end

        @testset "Array{<:Real,3}" begin
            q = Hankel.QDHT{1,2}(10, N; dim = 2)
            s = (M, N, K)
            fr, fk = ntuple(_ -> randn(rng, s), 2)
            test_frule(f, fk, q ⊢ NoTangent(), fr)
        end
    end

    @testset "$f(::Array{<:Real}, ::QDHT)" for f in (integrateR, integrateK)
        rng = MersenneTwister(27)
        N, M, K = 64, 5, 10
        q = Hankel.QDHT{1,2}(10, N)
        @testset "Vector" begin
            A = randn(rng, N)
            test_rrule(f, A, q ⊢ NoTangent())
            test_frule(f, A, q ⊢ NoTangent())
        end

        @testset "Matrix" begin
            @testset for dim in (1, 2)
                s = dim == 1 ? (N, M) : (M, N)
                sy = dim == 1 ? (1, M) : (M, 1)
                A = randn(rng, s)
                test_rrule(f, A, q ⊢ NoTangent(); fkwargs = (dim = dim,))
                test_frule(f, A, q ⊢ NoTangent(); fkwargs = (dim = dim,))
            end
        end

        @testset "Array{<:Real,3}" begin
            s = (M, N, K)
            sy = (M, 1, K)
            dim = 2
            A = randn(rng, s)
            test_rrule(f, A, q ⊢ NoTangent(); fkwargs = (dim = dim,))
            test_frule(f, A, q ⊢ NoTangent(); fkwargs = (dim = dim,))
        end
    end
end

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
function add_node!(bg::BondGraph, nodes)
    for node in nodes
        add_node!(bg, node)
    end
end

function add_node!(bg::BondGraph, node::AbstractNode)
    lg.add_vertex!(bg, node) || @warn "$node already in model"
end

function remove_node!(bg::BondGraph, nodes)
    for node in nodes
        remove_node!(bg, node)
    end
end

function remove_node!(bg::BondGraph, node::AbstractNode)
    lg.rem_vertex!(bg, node) || @warn "$node not in model"
    for bond in filter(bond -> node in bond, bg.bonds)
        lg.rem_edge!(bg, srcnode(bond), dstnode(bond))
    end
end

function connect!(bg::BondGraph, srcnode::AbstractNode, dstnode::AbstractNode;
        srcportindex=nextsrcport(srcnode), dstportindex=nextdstport(dstnode))
    srcnode in bg.nodes || error("$srcnode not found in bond graph")
    dstnode in bg.nodes || error("$dstnode not found in bond graph")
    srcport = Port(srcnode, srcportindex)
    dstport = Port(dstnode, dstportindex)
    return lg.add_edge!(bg, srcport, dstport)
end

function disconnect!(bg::BondGraph, node1::AbstractNode, node2::AbstractNode)
    # rem_edge! removes the bond regardless of the direction of the bond
    deleted_bond = lg.rem_edge!(bg, node1, node2)
    if isnothing(deleted_bond) # if returned nothing, try flipping node1 and node2
        deleted_bond = lg.rem_edge!(bg, node2, node1)
    end
    return deleted_bond
end

# TODO
# Flip bond function

function swap!(bg::BondGraph, oldnode::AbstractNode, newnode::AbstractNode)
    _check_port_number(oldnode,newnode)

    # may be a redundant check
    if !lg.has_vertex(bg, newnode)
        add_node!(bg, newnode)
    end

    srcnodes = lg.inneighbors(bg, oldnode)
    dstnodes = lg.outneighbors(bg, oldnode)
    remove_node!(bg, oldnode)

    for src in srcnodes
        connect!(bg, src, newnode)
    end
    for dst in dstnodes
        connect!(bg, newnode, dst)
    end
end

_check_port_number(oldnode::AbstractNode, newnode::AbstractNode) =
    numports(newnode) >= numports(oldnode) || error("New node must have a greater or equal number of ports to the old node")
_check_port_number(oldnode::AbstractNode, newnode::Junction) = true

# TODO implement according to https://bondgraphtools.readthedocs.io/en/latest/api.html#BondGraphTools.expose
# function expose!()

# end

# Inserts an AbstractNode between two connected (bonded) nodes
# The direction of the original bond is preserved by this action
function insert_node!(bg::BondGraph, bond::Bond, newnode::AbstractNode)
    src = srcnode(bond)
    dst = dstnode(bond)

    disconnect!(bg, src, dst)

    try
        add_node!(bg, newnode)
        connect!(bg, src, newnode)
        connect!(bg, newnode, dst)
    catch e
        # if connection fails, reconnect original bond
        disconnect!(bg, src, newnode)
        disconnect!(bg, newnode, dst)
        connect!(bg, src, dst)
        error(e)
    end
end
function insert_node!(bg::BondGraph, tuple::Tuple, newnode::AbstractNode)
    bonds = getbonds(bg, tuple)
    isempty(bonds) && error("$(tuple[1]) and $(tuple[2]) are not connected")
    insert_node!(bg, bonds[1], newnode)
end

function merge_nodes!(bg::BondGraph, node1::AbstractNode, node2::AbstractNode; junction=EqualEffort())
    node1.type == node2.type || error("$(node1.name) must be the same type as $(node2.name)")

    # node1 taken as the node to keep
    for src in lg.inneighbors(bg, node1)
        junc_src = deepcopy(junction)
        bond = getbonds(bg, src, node1)[1]
        insert_node!(bg, bond, junc_src)
        swap!(bg, node2, junc_src)
    end
    for dst in lg.outneighbors(bg, node1)
        junc_dst = deepcopy(junction)
        bond = getbonds(bg, node1, dst)[1]
        insert_node!(bg, bond, junc_dst)
        swap!(bg, node2, junc_dst)
    end
end
function merge_nodes!(bg::BondGraph, node1::Junction, node2::Junction)
    # node1 taken as the node to keep
    # remove conflicting connections between junctions if they exist
    disconnect!(bg, node1, node2)
    shared_neighbors = intersect(lg.all_neighbors(bg, node1), lg.all_neighbors(bg, node2))
    for shared_neighbor in shared_neighbors
        disconnect!(bg, node2, shared_neighbor)
    end
    swap!(bg, node2, node1)
end

function simplify_junctions!(bg::BondGraph; remove_redundant=true, squash_identical=true)
    junctions = filter(n -> n isa Junction, bg.nodes)

    # Removes junctions with 2 or less connected ports
    if remove_redundant
        for j in junctions
            n_nbrs = length(lg.all_neighbors(bg, j))
            if n_nbrs == 2
                #srcnode = lg.inneighbors(bg, j)[1]
                #dstnode = lg.outneighbors(bg, j)[1]
                node1, node2 = lg.all_neighbors(bg, j)
                remove_node!(bg, j)
                # bond direction may not be preserved here
                connect!(bg, node1, node2)
            elseif n_nbrs < 2
                remove_node!(bg, j)
            end
        end
    end

    # Squashes identical copies of the same junction type into one junction
    if squash_identical
        for j in junctions, nbr in lg.all_neighbors(bg, j)
            lg.has_vertex(bg, j) || continue # in case j was removed
            if type(j) == type(nbr)
                merge_nodes!(bg, j, nbr)
            end
        end
    end
    bg
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
        # isleft
        l = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(1, 0))
        @test isleft(Point2D{T}(0, 1) , l)
        @test !isleft(Point2D{T}(0, -1) , l)
        @test !isleft(Point2D{T}(0, -1e-6) , l)
        @test isleft(Point2D{T}(0, 1e-6) , l)
        @test isleft(Point2D{T}(0.5, 0) , l)

        # nearest_point
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 0)
        𝘅₃ = Point2D{T}(1, 1)
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        p = Point2D{T}(1, 1.1)
        r, p_c = nearest_point(p, q)
        @test r ≈ 0.5
        @test 𝘅₃ ≈ p_c

        # isleft
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 0)
        𝘅₃ = Point2D{T}(1, 1)
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        @test !isleft(Point2D{T}(1, 0), q)
        @test isleft(Point2D{T}(1, 2), q)
        @test !isleft(Point2D{T}(1, 0.9), q)

        # in
        @test Point2D{F}(0,0) ∈ plane
        @test Point2D{F}(4,4) ∈ plane
        @test Point2D{F}(1,2) ∉ plane

        # in_halfspace
        @test in_halfspace(Point2D{F}(0,0), plane)
        @test in_halfspace(Point2D{F}(0,1), plane)
        @test !in_halfspace(Point2D{F}(0,-1), plane)

        # in
        @test Point3D{F}(1,0,2) ∈ plane
        @test Point3D{F}(2,2,2) ∈ plane
        @test Point3D{F}(1,0,0) ∉ plane

        # in_halfspace
        @test in_halfspace(Point3D{F}(0,0,2), plane)
        @test in_halfspace(Point3D{F}(0,0,3), plane)
        @test !in_halfspace(Point3D{F}(0,0,-1), plane)

        # in
        @test Point2D{T}(2, 1) ∈ aab
        @test Point2D{T}(3, 1) ∈ aab
        @test !(Point2D{T}(4, 1) ∈ aab)
        @test !(Point2D{T}(2, 5) ∈ aab)

        # in
        @test Point3D{T}(2, 1, 0) ∈ aab
        @test Point3D{T}(3, 1, 0) ∈ aab
        @test !(Point3D{T}(4, 1, 0) ∈ aab)
        @test !(Point3D{T}(2, 5, 0) ∈ aab)
        @test !(Point3D{T}(2, 1, 2) ∈ aab)

        # in
        p = Point2D{T}(1//2, 1//10)
        @test p ∈ tri
        p = Point2D{T}(1//2, -1//10)
        @test p ∉ tri
        # in
        p = Point3D{T}(0, 1//2, 1//10)
        @test p ∈ tri
        p = Point3D{T}(0, 1//2, -1//10)
        @test p ∉ tri
        p = Point3D{T}(1//100, 1//2, 1//10)
        @test p ∉ tri

        # in
        p = Point2D{T}(1//2, 1//10)
        @test p ∈  quad
        p = Point2D{T}(1//2, -1//10)
        @test p ∉ quad

                    # in
            p₁ = Point2D{T}(0, 0)
            p₂ = Point2D{T}(2, 0)
            p₃ = Point2D{T}(2, 2)
            p₄ = Point2D{T}(3//2, 1//4)
            p₅ = Point2D{T}(3, 1)
            p₆ = Point2D{T}(1, 1)
            tri6 = QuadraticTriangle(p₁, p₂, p₃, p₄, p₅, p₆)
            @test Point2D{T}(1, 1//2) ∈  tri6
            @test Point2D{T}(1, 0) ∉  tri6

        # in
        p₁ = Point2D{T}(0, 0)
        p₂ = Point2D{T}(2, 0)
        p₃ = Point2D{T}(2, 3)
        p₄ = Point2D{T}(0, 3)
        p₅ = Point2D{T}(3//2, 1//2)
        p₆ = Point2D{T}(5//2, 3//2)
        p₇ = Point2D{T}(3//2, 5//2)
        p₈ = Point2D{T}(0, 1)
        quad8 = QuadraticQuadrilateral(p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈)
        @test Point2D{T}(1, 1) ∈  quad8
        @test Point2D{T}(1, 0) ∉  quad8

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using DelimitedFiles, Plots
# Training takes a while, so save results of simulations and
# have this separate plot script so can tweak plots without
#  rerunning simulations

# or whatever path to data
er_outfile = "er_data.txt"
plt_outfile = "plt_data.txt"

# er_dat = [Er(1,err_n) Er(2,err_n) Er(1,err_l) Er(2,err_l)]
# plt_dat = [[t for t in tsave]  σ_approx_n σ_approx_l σ_exact]

er_dat = readdlm(er_outfile, '\t', Float32, '\n')
plt_dat = readdlm(plt_outfile, '\t', Float32, '\n')

p1 = plot(er_dat[:,1],xscale=:log10,yscale=:log10,ylabel="Error",xlabel="Training steps",
        label="Training error, Neural net",legend = :bottomleft,linestyle = :dot,
        lw=4,titlefont = "Helvetica", legendfont = "Helvetica",color=2)
plot!(p1,er_dat[:,2],label="Testing error, Neural net",lw=4,color=2)
plot!(p1,er_dat[:,3], label="Training error, linear model",lw=4,color=1,linestyle = :dot)
plot!(p1,er_dat[:,4],label="Testing error, linear model",lw=4,color=1)

p2 = plot(plt_dat[:,1],plt_dat[:,3],m=:hexagon, lw=3, ms=3,
          titlefont = "Helvetica", legendfont = "Helvetica",
          label="Linear model",ylabel="stress",xlabel="time")
plot!(p2,plt_dat[:,1],plt_dat[:,2],m=:circle,lw=3,ms=3,
        label="NN solution")
plot!(p2,plt_dat[:,1],plt_dat[:,4],label="True solution",lw=1,lc=:black,
        xlims=(0,7),leg=:bottom)

pt = plot(p1,p2,layout=(2,1),size=(600,500))

# or path to wherever
savefig(pt, "FENEP.pdf")


text:
@testset "filelist.jl" begin
    dir = pwd()
    opts = (join = true, )
    fs = readdir(dir; opts...)
    @test filelist(dir; opts...) == fs[isfile.(fs)]
    @test isequal(filelist(dir; join=false), basename.(filelist(dir; join=true)))
    @test isequal(filelist(dir; sort=true), sort(filelist(dir; sort=false)))

    @test isequal(filelist(dir), FileTools.datalist(r".*",dir;find="file").fullpath)
end


text:
immutable Chi <: ContinuousUnivariateDistribution
    chisqd::Chisq

    Chi(df::Real) = new(Chisq(df))
end

@distr_support Chi 0.0 Inf

#### Parameters

dof(d::Chi) = dof(d.chisqd)
params(d::Chi) = (dof(d),)

#### Statistics

mean(d::Chi) = (k = dof(d); sqrt2 * gamma((k + 1.0) / 2.0) / gamma(k / 2.0))

var(d::Chi) = dof(d) - mean(d)^2

function skewness(d::Chi)
    μ, σ = mean(d), std(d)
    (μ / σ^3) * (1.0 - 2.0 * σ^2)
end

function kurtosis(d::Chi)
    μ, σ, γ = mean(d), std(d), skewness(d)
    (2.0 / σ^2) * (1 - μ * σ * γ - σ^2)
end

function entropy(d::Chi)
    k = dof(d)
    lgamma(k / 2.0) - log(sqrt(2.0)) -
        ((k - 1.0) / 2.0) * digamma(k / 2.0) + k / 2.0
end

function mode(d::Chi)
    k = dof(d)
    k >= 1.0 || error("Chi distribution has no mode when df < 1")
    sqrt(k - 1.0)
end

#### Evaluation

pdf(d::Chi, x::Float64) = exp(logpdf(d, x))

function logpdf(d::Chi, x::Float64)
    k = dof(d)
    (1.0 - 0.5 * k) * logtwo + (k - 1.0) * log(x) - 0.5 * x^2 - lgamma(0.5 * k)
end

gradlogpdf(d::Chi, x::Float64) = x >= 0.0 ? (dof(d) - 1.0) / x - x : 0.0

cdf(d::Chi, x::Float64) = cdf(d.chisqd, x^2)
ccdf(d::Chi, x::Float64) = ccdf(d.chisqd, x^2)
logcdf(d::Chi, x::Float64) = logcdf(d.chisqd, x^2)
logccdf(d::Chi, x::Float64) = logccdf(d.chisqd, x^2)

quantile(d::Chi, p::Float64) = sqrt(quantile(d.chisqd, p))
cquantile(d::Chi, p::Float64) = sqrt(cquantile(d.chisqd, p))
invlogcdf(d::Chi, p::Float64) = sqrt(invlogcdf(d.chisqd, p))
invlogccdf(d::Chi, p::Float64) = sqrt(invlogccdf(d.chisqd, p))

#### Sampling

rand(d::Chi) = sqrt(rand(d.chisqd))


text:
#### Travis

export Travis

"""
    Travis

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Travis
    "Julia version"
    julia_version::String
    "Checks"
    checks::String="continuous-integration/travis-ci/push"
end

function yaml(ci::Travis)
    return "language: julia
julia:
  - $(ci.julia_version)

script:
  - julia --project --color=yes --check-bounds=yes -e 'using Pkg;
                                                       Pkg.instantiate();
                                                       Pkg.build();
                                                      '

  - julia --project --color=yes --check-bounds=yes -e 'using Pkg;
                                                       Pkg.test();
                                                      '
"
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
TABLE_MIN = """id,height,age,name,occupation
1,10.0,1,string1,2012-06-15 00:00:00
2,10.1,2,string2,2013-06-15 01:00:00
3,10.2,3,string3,2014-06-15 02:00:00
4,10.3,4,string4,2015-06-15 03:00:00
5,10.4,5,string5,2016-06-15 04:00:00
"""

TABLE_WEIRD = """
a_dict,an_array,a_geopoint,a_date,a_time
{"test":3},"[1,2,3]","45.2,26.1","2014-06-15","02:00:00"
"""

@testset "Inference of a Schema from a table" begin

    @testset "A minimal Schema" begin
        t = Table(IOBuffer(TABLE_MIN))
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "integer"
        @test s.fields[2].typed == "number"
        @test s.fields[3].typed == "integer"
        @test s.fields[4].typed == "string"
        @test s.fields[5].typed == "string" # TODO: date
    end

    @testset "Inline with minimal schema" begin
        t = Table(IOBuffer(TABLE_MIN))
        TableSchema.infer(t)
        @test t.schema.fields[1].typed == "integer"
    end

    @testset "From a weird Table" begin
        t = Table(IOBuffer(TABLE_WEIRD))
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "object"
        @test s.fields[2].typed == "array"
        @test s.fields[3].typed == "geopoint"
        @test s.fields[4].typed == "date"
        @test s.fields[5].typed == "time"
    end

    @testset "One that does not meet constraints" begin
        t = Table("../data/data_constraints.csv")
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "integer"
        @test s.fields[2].typed == "number"
        @test s.fields[3].typed == "integer"
        @test s.fields[4].typed == "string"
        @test s.fields[5].typed == "string"
    end

    @testset "From data in a basic file" begin
        t = Table("../data/data_infer.csv")
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "integer"
        @test s.fields[2].typed == "integer"
        @test s.fields[3].typed == "string"
    end

    @testset "From a UTF8 file" begin
        t = Table("../data/data_infer_utf8.csv")
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "integer"
        @test s.fields[2].typed == "integer"
        @test s.fields[3].typed == "string"
    end

    @testset "From a ISO-8859-7 file" begin
        t = Table("../data/data_infer_iso-8859-7.csv")
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "integer"
        @test s.fields[2].typed == "integer"
        @test s.fields[3].typed == "string"
    end

    @testset "From 'row limit' file" begin
        t = Table("../data/data_infer_row_limit.csv")
        tr = TableSchema.read(t, cast=false)
        s = Schema()
        TableSchema.infer(s, tr, t.headers)
        @test s.fields[1].typed == "string"
        @test s.fields[2].typed == "string"
        @test s.fields[3].typed == "string"
    end

end