text:
# Integers

# Literals generate a value of the "default" integer type (Int64 on my workstation)
1
1234

# Types of Literals
typeof(1)

# Internal variable, `Sys.WORD_SIZE`, indicates word architecture of target system
Sys.WORD_SIZE

# Synonyms for signed and unsigned default integers
Int
UInt

# Values that **do not** fit in 32-bits but **can** be represented in 64-bits **always** create a 64-bit integer
typeof(3_000_000_000)

# Similar promotion occurs from 64-bits to 128-bits on my workstation
typeof(3_000_000_000_000)
typeof(3_000_000_000_000_000)
typeof(3_000_000_000_000_000_000)
typeof(3_000_000_000_000_000_000_000)

# Unsigned integers are input and output using the `0x` prefix and hexadecimal digits '0-9a-f'. The size of the
# value is determined by the number of hexadecimal digits used.
x = 0xff
typeof(x)

x = 0xffff
typeof(x)

x = 0xfffffff
typeof(x)

x = 0xffffffff
typeof(x)

x = 0xffffffffff
typeof(x)

x = 0x1
typeof(x)

x = 0x123
typeof(x)

x = 0x1234567
typeof(x)

x = 0x123456789abcdef
typeof(x)

x = 0x11112222333344445555666677778888
typeof(x)

# Binary and octal literals are also supported
x = 0b10
typeof(x)

x = 0o010
typeof(x)

# 16 leading zeros
x = 0x0000000000000000111122223333444
typeof(x)

# 18 leading zeros
x = 0x000000000000000000111122223333444
typeof(x)

# If one preceeds a binary, octal or hexadecimal literal, the value is the two's complement of the literal without
# the sign
x = -0x2
typeof(x)

x = -0x0002
typeof(x)

# Use `typemin` and `typemax` to obtain the minimum and maximum representable values of the primitive numeric
# types such as integers.
(typemin(Int32), typemax(Int32))

for T in [Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]
    println("$(lpad(T, 7)): [$(typemin(T)), $(typemax(T))]")
end

# Overflow

# In Julia, exceeding the maximum representable value of a given type results in wrap around behavior.
x = typemax(Int64)
x + 1
x + 1 == typemin(Int64)

x = typemin(Int64)
x - 1
x - 1 == typemax(Int64)

# Beware. In applications where overflow is possible, explicit checking for wrap around behavior produced by
# overflow (and underflow) is essential; otherwise the `BigInt` type in _Arbitrary Precision Arithmetic_ is
# recommended instead.

# An example of overflow behavior and how to potentially resolve it is as follows:
10^19  # Does not produce expected value because of overflow

big(10)^19  # Use `BigInt` instead

# Division Errors

# Integer division, the `div` function, has two exceptional cases:
# - Dividing by zero
# - Dividing the lowest negative number (`typemin`) by -1
# Both of these cases produce a `DivideError`
div(1, 0)

div(typemin(Int64), -1)

# The remainder and modulus functions, `rem` and `mod`, respectively, throw a `DivideError` when their second
# argument is zero.
rem(1, 0)
mod(1, 0)

# Floating-point numbers

# Literal floating-point numbers are represented using the "standard" formats using e-notation when necessary:
1.0
typeof(1.0)

1.
typeof(1.)

0.5
typeof(0.5)

.5
typeof(.5)

-1.23
typeof(-1.23)

1e10
typeof(1e10)

2.5e-4
typeof(2.5e-4)

# The previous values are all of type `Float64`. Literal `Float32` values can be entered by writing an `f` in place
# of the `e`
x = 0.5f0
typeof(x)

x = 2.5f-4
typeof(x)

# One **cannot** simply use the `f` suffix
x = 0.5f

# One can easily convert values to `Float32` values
x = Float32(-1.5)
typeof(x)

# Hexadecimal floating-point literals are also valid, but only as `Float64` values, with `p` preceding the base-2
# exponent. One must understand how `Float64` values are stored using the IEEE standard. See
# [the Wikipedia article](https://en.wikipedia.org/wiki/Floating-point_arithmetic)
0x1p0

0x1.8p3

x = 0x.4p-1
typeof(x)

# Julia also supports half-precision floating-point numbers (`Float16`); however, they are implemented in software
# and use `Float32` for calculations
sizeof(Float16(4.))

2 * Float16(4.)

# Floating-point zero

# Floating-point numbers have **two** zeros: positive zero and negative zero. They are equal to each other, but
# actually have different binary representations. The different representations can be seen using the `bitstring`
# function.

0.0 == -0.0

bitstring(0.0)

bitstring(-0.0)

# Special floating-point values
Inf16  ## `Float16` infinity
Inf32  ## `Float32` infinity
Inf  ## `Float64` infinity

-Inf16  ## `Float16` negative infinity
-Inf32  ## `Float32` negative infinity
-Inf  ## `Float64` negative infinity

NaN16  ## `Float16` not a number
NaN32  ## `Float32` not a number
NaN  ## `Float64` not a number

# Comparing these values may not be obvious. See the section, "Numeric Comparisons."

# These special floating-point values, per the _IEEE-754 standard_, are the result of certain arithmetic
# operations.
1 / Inf

1 / 0

-5 / 0

0.000001 / 0

0 / 0

500 + Inf

500 - Inf

Inf + Inf

Inf - Inf

Inf * Inf

Inf / Inf

0 * Inf

# The `typemin` and `typemax` functions also apply to floating-point types.
(typemin(Float16), typemax(Float16))

(typemin(Float32), typemax(Float32))

(typemin(Float64), typemax(Float64))

# Machine epsilon

# Most real numbers **cannot** be represented with floating-point numbers. For many purposes, it is important to
# know the distance between two adjacent floating-point numbers. This distance is often known as the
# "machine epsilon."

# Julia provides `eps`, which gives the distance between 1.0 and the next larger representable floating-point
# value.
eps(Float32)  # Same as 2.0^-23 as a `Float32` value
Float32(2.0)^-23
eps(Float64)  # Same as 2.0^-52 as a `Float64` value
Float64(2.0^-52)
eps()  # Same as `eps(Float64)`

# The `eps` function also accepts a floating point value. In this situation, the result of `eps()` is the distance
# between the argument to `eps()` and the next (higher) representable floating-point number.
eps(1.0)
eps(1000.)
eps(1e-27)
eps(0.0)

# Note that the distance between two adjacent floating-point values is **not** uniform throughout the
# floating-point range; instead, floating-point numbers are densest close to 0.0 and grow exponentially sparser as
# one moves away from zero.

# Julia also provides `nextfloat` and `prevfloat` functions that return the next largest or smallest
x = 1.25f0
nextfloat(x)
prevfloat(x)
bitstring(x)
bitstring(nextfloat(x))
bitstring(prevfloat(x))

x = 1.25
nextfloat(x)
prevfloat(x)
bitstring(x)
bitstring(nextfloat(x))
bitstring(prevfloat(x))

# Arbitrary Precision Arithmetic

# Julia wraps the _GNU Multiple Precision Arithmetic Library (GMP)_ and the _GNU MPFR Library_, respectively. The
# `BigInt` and `BigFloat` types are available in Julia for arbitrary precision integer and floating point numbers,
# respectively.

# Use constructors to create these types from either primitive numerical types. Additionally, one can use the
# _string literal_ `@big_str` or `parse` to construct instances from `AbstractString` values. Values of type
# `BigInt` can be constructed from integer literals when the exceed the bounds of `Int128`.

# Note that **no** unsigned arbitrary-precision integer types exist in _Base_; instead, `BigInt` is sufficient in
# most cases. Additionally, binary, octal and hexadecimal integer literals can be used (in addition to decimal
# literals).

# Once created, arbitrary-precision integers and floating-point values participate in arithmetic with all other
# numeric types because of Julia's _type promotion and conversion mechansims_.

BigInt(typemax(Int64)) + 1

big"1234567890123456789012345678901234567890" + 1

parse(BigInt, "123456789012345678901234567890") + 1

string(big"2"^200, base=16)

typeof(0x1000_0000_0000_0000_0000_0000_0000_0000)
typeof(0x1000_0000_0000_0000_0000_0000_0000_0000_0)
0x1000_0000_0000_0000_0000_0000_0000_0000_0 - 1 == typemax(UInt128)

0x0000_0000_0000_0000_0000_0000_0000_0000_0
typeof(ans)

big"1.23456789012345678901"

parse(BigFloat, "1.23456789012345678901")

BigFloat(2.0^66) / 3

factorial(BigInt(40))

# However, type promotion between the primitive types above and `BigInt` / `BigFloat` is **not** automatic and must
# be explicity stated.

x = typemin(Int64)
x = x - 1
typeof(x)

y = BigInt(typemin(Int64))
y = y - 1
typeof(y)

# The default precision (in number of bits of the significand) and rounding mode of `BigFloat` operations can be
# changed globally by calling `setprecision` and `setrounding`, and all further calculations will take these
# changes into account. Alternatively, the precision or the rounding can be changed only within the execution of a
# particular block of code using the same functions with a `do` block.

setrounding(BigFloat, RoundUp) do
    BigFloat(1) + parse(BigFloat, "0.1")
end

setrounding(BigFloat, RoundDown) do
   BigFloat(1) + parse(BigFloat, "0.1")
end

setprecision(40) do
   BigFloat(1) + parse(BigFloat, "0.1")
end

# Numeric Literal Coefficents

# To make common numeric forulae and expressions clearer, Julia allows variables to be immediately preceded by a
# numeric literal, implying multiplication. This makes writing polynomial expressions much clearer.

x = 3
2x^2 - 3x + 1
1.5x^2 - .5x + 1

# It also makes writing exponential functions more elegant

2^2x

# The precedence of numeric literal coefficents is slightly lower than that of unary operators such as negation.
# So, -2x is parsed as (-2) * x and √2x is parsed as (√2) * x. However, numeric literal coefficients parse
# similarly to unary operators when combined with exponentiation. For example, 2^3x is parsed as 2 * (3x), and
# 2x^3 is parsed as 2 * (x ^ 3).

# Numeric literals also work as coefficients to parenthesized expressions.

2(x - 1)^2 - 3(x - 1) + 1

# **Note**
#
# The precedence of numeric literal coefficients used for implicit multiplication is higher than other binary
# operators such as multiplication (*) and division(/, \, and //). This means, for example, that 1 / 2im equals
# -0.5im and 6 // 2(2 + 1) equals 1 // 1.

# Additionally, parenthesized expressions can be used as coefficients to variables, implying multiplication of the
# expression by the variable:

(x - 1)x

# Neither juxtaposition of two parenthesized, nor placing a variable before a parenthesized expression, however,
# can be used to imply multiplication.

(x + 1)(x + 1)

x(x + 1)

# Instead, both expressions are interpreted as function application; that is, any expression that is **not** a
# numeric literal, when immediately followed by a parenthetical, is interpreted as a function application to the
# values in parentheses.

# The above syntatic enhancements significantly reduce the visual noise incurred when writing common mathematical
# formulae. Finally, note that **no** whitespace may come between a numeric literal coefficient and the identifier
# or parenthisized expression which it multiplies.


text:
## Genetic Algorithm
##

## fitness function for a subarray
## cfg: GASS parameters
## subarrid: subarray id (int)
## subind: subarray indices in the main Arr.
##
function fitness_subarray(cfg, subarrid, subind)
    subarr = cfg.arr[subind,:]

    bl= calc_baselines(subarr)
    uv= calc_uv(bl, cfg.obs.Source_Hour_Angle ,  cfg.obs.Source_Declination)
    h , dr=  calc_dirtybeam(uv , 512, 256, robust=0.5)
    b= fit_beam(h , dr)
    mrs= calc_mrs(uv)

    res= 0
    res += cfg.wei.Weight_Spatial_Resolution[subarrid]*abs(b.ar-cfg.sub.Spatial_Resolution[subarrid])
    res += cfg.wei.Weight_Elongation[subarrid]*abs(b.e-cfg.sub.Elongation[subarrid])*sign(b.e-
        cfg.sub.Elongation[subarrid])
    res += cfg.wei.Weight_Sidelobe_Levels[subarrid]*abs(b.sidelobe-
        cfg.sub.Sidelobe_Level[subarrid])*sign(b.sidelobe-cfg.sub.Sidelobe_Level[subarrid])
    res += cfg.wei.Weight_Maximum_Recoverable_Scale[subarrid]*abs(mrs-
        cfg.sub.Maximum_Recoverable_Scale[subarrid])*sign(cfg.sub.Maximum_Recoverable_Scale[subarrid]-mrs)

    if cfg.obs.Display_Verbose
      println("### Subarray fitness")
      println(subarr)
      @printf("## beam:")
      println(b)
      @printf("## MRS: %3.3f \n", mrs)
      @printf("## Fitness: %3.3f (smaller is better) \n", res)
    end

    return(res , [b , mrs])
end

### Creation of a population
###
function create_population(cfg)

    pop= Array{Array{Int,1},2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)
    fitness= Array{Float64,2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)
    score=  Array{Float64,1}(undef,cfg.ga.Population_Size)
    paramsub= Array{Dict{String,Float64},2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)

    subind= collect(1:cfg.obs.Antenna_Number)

    for i in 1:cfg.ga.Population_Size
        Random.shuffle!(subind)
        for j in 1:cfg.obs.Subarray_Number
            pop[i,j]= subind[cfg.sub.Subrange[j]]
            fitness[i,j] , res= fitness_subarray(cfg, j, pop[i,j])

            if cfg.obs.Display_Verbose
              @printf("### subarr %d, %d ; fitness: %f (smaller is better) \n" , i,j,fitness[i,j])
            end

            paramsub[i,j]= Dict("ar"=>res[1].ar,"e"=>res[1].e, "sidelobe"=>res[1].sidelobe, "mrs"=>res[2])
        end
        score[i]= -sum(cfg.wei.Weight_Subarray[:] .* fitness[i,:])
    end

    pop0= population(0, pop, fitness, score , paramsub)

    return(pop0)
end

### get best parents...
###
function get_elitism(cfg, pop::population)
    isort= reverse(sortperm(pop.score))

    elit= []
    for i in 1:cfg.ga.Number_Elitism
        push!(elit, [pop.subarr[isort[i],:], pop.fitness[isort[i],:], pop.score[isort[i]]])
    end
    return(elit)
end

## ... using tournament selection
function get_parents(cfg, pop::population)
    tour= cfg.ga.Tournament_Size

    popran= Random.shuffle!(collect(1:cfg.ga.Population_Size))
    isort= reverse(sortperm(pop.score[popran[1:tour]]))
    iparent1= popran[isort[1]]

    iparent2= iparent1
    while iparent1==iparent2
        popran= Random.shuffle!(collect(1:cfg.ga.Population_Size))
        isort= reverse(sortperm(pop.score[popran[1:tour]]))
        iparent2= popran[isort[1]]
    end
    return(iparent1 , iparent2)
end

#### auxiliary functions
function flatten_subarray(cfg, subarr)
    subflat= zeros(Int, cfg.obs.Antenna_Number)

    counter= 1
    for i in 1:cfg.obs.Subarray_Number
        for j in 1:cfg.sub.Pads_Per_Subarray[i]
            subflat[counter]= subarr[i][j]
            counter += 1
        end
    end
    return(subflat)
end

### recreate the subarray array..
function wrap_subarray(cfg, sub)
    subarray= Array{Array{Int,1}}(undef,0)

    for i in 1:cfg.obs.Subarray_Number
        push!(subarray, sub[cfg.sub.Subrange[i]])
    end

    return(subarray)
end

###
### crossover (#1) using the two parents
###
function get_crossover1(cfg, pop::population , parents)
    nmax= cfg.obs.Antenna_Number
    p1= flatten_subarray(cfg, pop.subarr[parents[1],:])
    p2= flatten_subarray(cfg, pop.subarr[parents[2],:])
    child= zeros(Int, nmax)

    pivot= rand(1:nmax,2)
    pivot1= minimum(pivot) ; pivot2= maximum(pivot)
    child[pivot1:pivot2] .= p1[pivot1:pivot2]

    irange1= 1:pivot1-1 ; irange2= pivot2+1:nmax

    cnterp2= 1
    for i in irange1
        while (p2[cnterp2] in child)
            cnterp2 += 1
        end
        child[i]= p2[cnterp2]
    end
    for i in irange2
        while (p2[cnterp2] in child)
            cnterp2 += 1
        end
        child[i]= p2[cnterp2]
    end

    c= wrap_subarray(cfg, child)
    return(child)
end

### mutation
### loop over all alleles and swap two if p is true
###
function get_mutation(cfg, child)
    nmax= cfg.obs.Antenna_Number
    si= copy(child)

    pmutation= 1-cfg.ga.Mutation_Rate
    rng = Random.MersenneTwister()
    rd= Random.rand(rng, nmax)
    # println(rd)
    for i in 1:nmax
        if rd[i] > pmutation
            iswap= Random.rand(rng, big.(1:nmax))
            value= si[iswap]
            si[iswap]= si[i]
            si[i]= value
        end
    end
    return(si)
end

###
### Evolve a population
###
function get_evolution(cfg, pi::population)
    npop= cfg.ga.Population_Size
    age= pi.age + 1
    nelit= cfg.ga.Number_Elitism

    ## Elitism
    pelit=  get_elitism(cfg, pi)

    ## parent selection, crossover and mutation
    ncross= npop-nelit

    crosspop= []
    for i in 1:ncross
        parent= get_parents(cfg, pi)
        child=  get_crossover1(cfg, pi , parent)
        mutated= get_mutation(cfg, child)
        push!(crosspop, mutated)
        println(parent)
        println(child)
    end

    ## create an evolved population
    pinew= Array{Array{Int,1},2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)
    fitness= Array{Float64,2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)
    score=  Array{Float64,1}(undef,cfg.ga.Population_Size)
    paramsub= Array{Dict{String,Float64},2}(undef,cfg.ga.Population_Size,  cfg.obs.Subarray_Number)

    for i in 1:nelit
        for j in 1:cfg.obs.Subarray_Number
            pinew[i,j]= pelit[i][1][j]
            fitness[i,j], res= fitness_subarray(cfg, j, pinew[i,j])
            paramsub[i,j]= Dict("ar"=>res[1].ar,"e"=>res[1].e, "sidelobe"=>res[1].sidelobe, "mrs"=>res[2])
            #println(i," ",j," ",fitness[i,j])
        end
        score[i]= -sum(cfg.wei.Weight_Subarray[:] .* fitness[i,:])
    end

    println("fitness of crossover..")
    for i in nelit+1:npop
        println(crosspop[i-nelit])
        pwrap= wrap_subarray(cfg, crosspop[i-nelit])
        println(pwrap)
        for j in 1:cfg.obs.Subarray_Number
            pinew[i,j]= pwrap[j]
            fitness[i,j] , res= fitness_subarray(cfg, j, pinew[i,j])
            paramsub[i,j]= Dict("ar"=>res[1].ar,"e"=>res[1].e, "sidelobe"=>res[1].sidelobe, "mrs"=>res[2])
        end
        score[i]= -sum(cfg.wei.Weight_Subarray[:] .* fitness[i,:])
    end

    popnew= population(age, pinew, fitness, score, paramsub)

    return(popnew)
end

### Run the full optimization process for a cfg set of params.
###
function gass_optimization(cfg)
    println("\n##")
    println("## Creating the population...")
    p0= create_population(cfg)
    pelit=  get_elitism(cfg, p0)
    parent= get_parents(cfg, p0)
    child=  get_crossover1(cfg, p0 , parent)
    mutated= get_mutation(cfg, child)

    ## evolution
    println("\n## Evolution...")

    species= []
    for i in 1:cfg.ga.Number_Iterations
        p1= get_evolution(cfg, p0)
        push!(species,p1)
        p0= p1
        #if (i % 10) == 0
            @printf("\n## Iteration: %d \n",i)
            @printf("### Best score: %3.3f \n", p1.score[1])
        #end
    end

   return(species)
end

ext: jl

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
using ComponentArrays
using Optim

rosen(u) = (1.0 - u.x)^2 + 100.0 * (u.y - u.x^2)^2

u₀ = ComponentArray(x=0.0, y=0.0)
lb = ComponentArray(x=-0.5, y=-0.5)
ub = ComponentArray(x=0.5, y=0.5)

df = TwiceDifferentiable(rosen, u₀; autodiff=:forward)
dfc = TwiceDifferentiableConstraints(lb, ub)

val = optimize(df, dfc, u₀, IPNewton())


text:
@doc """
	convert_to_reversible!(model)

Converts the model by combining duplicate reactions that previously flowed in
 opposite directions

""" convert_to_reversible!

@doc """
	convert_to_reversible(model)

Converts a copy of the model by combining duplicate reactions that previously flowed in
 opposite directions

""" convert_to_reversible

@doc """
    change_objective_coef(lp, objective, value)

Change the objective value for the entire lp model. ``objective`` may be
a single index or an array of indices

### Example
    change_objective_coef(lp, [2,5], [1, 0.5])
""" change_objective_coef

@doc """
    lb, ub = get_variable_bounds(lp, [index])

Return the column bounds for the lp model. Optionally, choose an index to return the bounds
of a single column

### Example
    lb, ub = get_variable_bounds(lp, 13)

""" get_variable_bounds

@doc """
    b_lb, b_ub = get_constraint_bounds(lp)

Return the row bounds for the lp model
""" get_constraint_bounds

@doc """
    S = get_constraint_matrix(lp)

Return a sparse matrix filled with the coefficients of the reactions
""" get_constraint_matrix

@doc """
    direction = get_direction(lp)

Return either "maximize" or "minimize"
""" get_direction

@doc """
    obj_val = get_objective_value(lp)

Return the objective value of the lp model. Note that the model must have been solved
""" get_objective_value

@doc """
    red_cost = get_reduced_costs(lp)

Return the reduced cost vector of the lp model. Note that the model must have been solved
""" get_reduced_costs

@doc """
    sol_vector = get_solution(lp)

Return the solution vector of the lp model. Note that the model must have been solved
""" get_solution

@doc """
    sol_stat = get_solution_status(lp)

Return the solution status of the lp model. Returns one of the following

* Optimal
* No feasible solution
* Infeasible
* Feasible
* Unbound
* Undefined

""" get_solution_status

@doc """
    stat_code = get_solution_status_code(lp)

Return the solution status code for a solved lp model. The solution status codes
differ between solvers.
""" get_solution_status_code

@doc """
    slack = get_slack(lp)

Return the slack vector of the lp model. Note that the model must have been solved

""" get_slack

@doc """
    objective = get_objective(lp)

Return the objective vector of the lp model. This corresponds to model.c in the Model
""" get_objective

@doc """
    set_direction(lp, direction)

Set the direction of optimization. Direction may be "max" or "min"
""" set_direction

@doc """
    set_col_bounds(lp, index, lb, ub)
    set_col_bounds(lp, index, fixed_value)

Set the bounds for a column specified by **index** to either a lower- and upper- bound or to
a **fixed value**

### Example
    set_col_bounds(lp, 11, 0, 10)

To set the bounds for column (reaction) 11 to 0 and 10
""" set_col_bounds

@doc """
    set_objective(lp, objective, value)

Change the objective value for a single column (reaction). This does not affect other reactions
""" set_objective

@doc """
    solve(lp)

Solve the lp model
""" solve

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
using Flux3D, BenchmarkTools

# mesh = load_trimesh(joinpath(@__DIR__, "../assets/teapot.obj"))
mesh = load_trimesh("./assets/teapot.obj")

trail = @benchmark $get_edges_packed($mesh; refresh = true)
minimum(trail.times) * 1.0e-6

trail = @benchmark $get_faces_to_edges_packed($mesh; refresh = true)
minimum(trail.times) * 1.0e-6

trail = @benchmark $get_laplacian_packed($mesh; refresh = true)
minimum(trail.times) * 1.0e-6

trail = @benchmark $compute_verts_normals_packed($mesh)
minimum(trail.times) * 1.0e-6

trail = @benchmark $compute_faces_areas_packed($mesh)
minimum(trail.times) * 1.0e-6

trail = @benchmark $compute_faces_normals_packed($mesh)
minimum(trail.times) * 1.0e-6

trail = @benchmark $laplacian_loss($mesh)
minimum(trail.times) * 1.0e-6

trail = @benchmark $edge_loss($mesh)
minimum(trail.times) * 1.0e-6

trail = @benchmark $chamfer_distance($mesh, $mesh, 10000)
minimum(trail.times) * 1.0e-6

trail = @benchmark $sample_points($mesh, 10000)
minimum(trail.times) * 1.0e-6


text:
@testset "SLIC" begin
  Z = [ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)]
  𝒮 = georef((Z=Z,))
  p = partition(𝒮, SLIC(4, 1.0))
  @test length(p) == 4
  @test all(nelements.(p) .== 100)
  p1, p2, p3, p4 = p
  @test mean(coordinates(centroid(p1, ind)) for ind in 1:nelements(p1)) == [5.0,5.0]
  @test mean(coordinates(centroid(p2, ind)) for ind in 1:nelements(p2)) == [15.0,5.0]
  @test mean(coordinates(centroid(p3, ind)) for ind in 1:nelements(p3)) == [5.0,15.0]
  @test mean(coordinates(centroid(p4, ind)) for ind in 1:nelements(p4)) == [15.0,15.0]

  C = cluster(𝒮, SLIC(4, 1.0))
  @test C.cluster == vec(Z')

  𝒮 = georef((z=[√(i^2+j^2) for i in 1:100, j in 1:100],))
  p = partition(𝒮, SLIC(50, 0.001))
  @test length(p) == 49

  if visualtests
    @test_reference "data/slic.png" plot(p)
  end

  # test SLIC with heterogeneous data
  Z = (a=rand(10), b=1:10, x=rand(10), y=rand(10))
  𝒮 = georef(Z, (:x, :y))
  C = cluster(𝒮, SLIC(2, 1.0))
  @test domain(C) == domain(𝒮)
  @test Set(C.cluster) ⊆ Set(1:2)

  # test SLIC for orphaned points
  a = [0,0,0,0,0,0,0,0,0,0]
  x = [0.4993029939801461, 0.14954882636793432, 0.23118957975519616, 0.6816610871344635, 0.6665309965318731, 0.691522274292691, 0.012495903053589608, 0.9831177095525963, 0.4445263730141056, 0.2175871587746574]
  y = [0.32721108209880256, 0.11427387079564899, 0.826401075107011, 0.6164294766961782, 0.6562529361193601, 0.43388375115444644, 0.7624847842129086, 0.1516623758764959, 0.07641616063237144, 0.8669098569279463]
  Z = (a=a, x=x, y=y)
  𝒮 = georef(Z, (:x, :y))
  C = cluster(𝒮, SLIC(2, 1.0))
  @test Set(C.cluster) ⊆ Set(1:2)

  # test SLIC with weights in attribute columns
  z1 = [√((i-0)^2+(j-0)^2) for i in 1:100, j in 1:100]
  z2 = [√((i-100)^2+(j-100)^2) for i in 1:100, j in 1:100]
  𝒮 = georef((z1=z1, z2=z2))
  w1 = Dict(:z1 => 10, :z2 => 0.1)
  w2 = Dict(:z1 => 0.1, :z2 => 10)
  p1 = partition(𝒮, SLIC(50, 0.001, weights=w1))
  p2 = partition(𝒮, SLIC(50, 0.001, weights=w2))
  @test length(p1) == 49
  @test length(p2) == 49

  if visualtests
    @test_reference "data/slic-w1.png" plot(p1)
    @test_reference "data/slic-w2.png" plot(p2)
  end
end


text:
module Day10

function scores(line)
    stack = Char[]
    for ch ∈ line
        if ch ∈ "([{<"
            push!(stack, ch)
        else
            match = pop!(stack)
            if abs(Int(ch) - Int(match)) > 2
                corruption_score = ch == ')' ? 3 :
                    ch == ']' ? 57 :
                    ch == '}' ? 1197 :
                    25137
                return (corruption_score, 0)
            end
        end
    end
    # We came to the end without finding any corrupted input.
    # Now calculate the autocompletion score.
    completion_score = 0
    while !isempty(stack)
        ch = pop!(stack)
        increment = ch == '(' ? 1 :
            ch == '[' ? 2 :
            ch == '{' ? 3 :
            4
        completion_score *= 5
        completion_score += increment
    end
    return (0, completion_score)
end

function sum_scores(input)
    corruption_score = 0
    completion_scores = []
    for line in filter(!isempty, split(input, "\n"))
        corruption, completion = scores(line)
        corruption_score += corruption
        if completion != 0
            push!(completion_scores, completion)
        end
    end
    sort!(completion_scores)
    (corruption_score, completion_scores[(length(completion_scores) + 1) ÷ 2])
end

function run()
    test = """[({(<(())[]>[[{[]{<()<>>
[(()[<>])]({[<{<<[]>>(
{([(<{}[<>[]}>{[]{[(<()>
(((({<>}<{<{<>}{[]{[]{}
[[<[([]))<([[{}[[()]]]
[{[{({}]{}}([{[{{{}}([]
{<[[]]>}<{[{[{[]{()[[[]
[<(<(<(<{}))><([]([]()
<{([([[(<>()){}]>(<<{{
<{([{{}}[<[[[<>{}]]]>[]]"""

    @assert sum_scores(test) == (26397, 288957)

    input = read("resources/day10_input.txt", String)
    display(sum_scores(input))
end

end


text:
module ReferenceTests

using Test
using ImageCore
using Distances
using FileIO
using ImageInTerminal
using SHA
using DeepDiffs
using Random

export
    @withcolor,
    @io2str,
    @test_reference,
    psnr_equality

include("utils.jl")
include("test_reference.jl")
include("fileio.jl")
include("equality_metrics.jl")
include("render.jl")

if Base.VERSION >= v"1.4.2" && ccall(:jl_generating_output, Cint, ()) == 1
    @assert precompile(Tuple{Core.kwftype(typeof(test_reference)),NamedTuple{(:by,), Tuple{typeof(isequal)}},typeof(test_reference),String,Matrix{RGB{Float64}}})   # time: 0.4697069
    @assert precompile(Tuple{Core.kwftype(typeof(test_reference)),NamedTuple{(:by,), Tuple{typeof(isequal)}},typeof(test_reference),String,Matrix{RGBA{Float64}}})   # time: 0.24671553
    @assert precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG}},Matrix{RGB{Float64}},Function,Nothing})   # time: 0.029915236
    @assert precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG}},Matrix{RGBA{Float64}},Function,Nothing})   # time: 0.013356352
end

end # module


text:
# include("../../src/SFFM.jl")
using LinearAlgebra, Plots, SFFM
# #
# # ## define the model(s)
# include("exampleModelDef.jl")

## section 4.4: the sensitivity of the stationary distribution of X to rates r
Tfun(γ₂) = [
    -(γ₁ + γ₂) γ₂ γ₁ 0;
    β₂ -(γ₁ + β₂) 0 γ₁;
    β₁ 0 -(γ₂ + β₁) γ₂;
    0 β₁ β₂ -(β₂ + β₂);
    ]
let
    c = 0
    colours = [:green;:blue;:red]
    shapes = [:diamond,:x,:+]
    styles = [:solid,:dash,:dot]
    q = plot(layout = (1,2))
    for sp in 1:2
        q = plot!(windowsize = (600,250), subplot = sp)
    end
    for γ₂ in [11;16;22]
        c = c+1
        Ttemp = Tfun(γ₂)
        tempModel = SFFM.Model( Ttemp, C, r, Bounds = approxModel.Bounds)
        println("created tempModel with upper bound x=", tempModel.Bounds[1,end])
        ## mesh
        Δtemp = 0.4
        Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])

        Basis = "lagrange"
        nBases = 2
        mesh = SFFM.DGMesh(
            tempModel,
            Nodes,
            nBases,
            Basis=Basis,
        )

        # compute the marginal via DG
        All = SFFM.MakeAll( tempModel, mesh, approxType = "projection")
        Ψ = SFFM.PsiFun( All.D)

        # the distribution of X when Y first returns to 0
        ξ = SFFM.MakeXi( All.B, Ψ)

        marginalX, p, K = SFFM.MakeLimitDistMatrices(
            All.B,
            All.D,
            All.R.RDict,
            Ψ,
            ξ,
            mesh,
            tempModel,
        )
        println("For γ₂ = ",
            γ₂, ", χ⁰ = ",
            sum(p),
            ", χ¹ = ",
            sum(marginalX)-sum(p),
            " and total prob is ",
            sum(marginalX),
            ".",
        )
        tempDist = SFFM.Coeffs2Dist(
            tempModel,
            mesh,
            marginalX,
            SFFM.SFFMDensity,
        )
        temp = zeros(SFFM.NBases(mesh),SFFM.NIntervals(mesh),2)
        temp[:,:,1] = tempDist.distribution[:,:,1]+tempDist.distribution[:,:,2]
        temp[:,:,2] = tempDist.distribution[:,:,3]+tempDist.distribution[:,:,4]
        q = plot!(
            tempDist.x,
            temp[:,:,1],
            subplot = 1,
            legend = false,
            color = colours[c],
            seriestype = :line,
            linestyle = styles[c],
            markershape = shapes[c],
            xlabel = "x",
        )
        q = scatter!(
            [-0.2 + (c-1)*0.2],
            [sum(tempDist.pm)],
            subplot = 2,
            color = colours[c],
            label = :none,
        )
        q = plot!(
            tempDist.x[:,1],
            temp[:,1,2],
            subplot = 2,
            color = colours[c],
            label = "γ₂: "*string(γ₂),
            linestyle = styles[c],
            markershape = shapes[c],
        )
        q = plot!(
            tempDist.x[:,2:end],
            temp[:,2:end,2],
            subplot = 2,
            color = colours[c],
            label = :none,
            markershape = shapes[c],
            linestyle = styles[c],
            xlabel = "x",
        )

        println("")
    end
    titles = ["Phases 11 + 10" "Phases 01 + 00"]
    for sp in 1:2
        q = plot!(
            subplot = sp,
            xlims = (-0.5,8),
            title = titles[sp],
            ylabel = "Density / Probability",
            grid = false,
        )
    end
    display(q)
    # savefig(pwd()*"/examples/paperNumerics/dump/sensitivityMarginalStationaryDistX.png")
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using DifferentialEquations
using Plots
using Distributions

n_parasites = 100

c = 4.0;
ux = 0.2;
ux1 = fill(0.2, n_parasites); #le ux et uy 1000 à cause de la β # une autre façon :[0.2 for x in 1:1000]
Random.seed!(1234);
uy = rand(200:1000, n_parasites)/1000;
Random.seed!(1235);
r1 = rand(Float64, n_parasites);
Random.seed!(1236);
r2 = rand(Float64, n_parasites);
Random.seed!(1237);
r3 = rand()*r1
bx = 1.0;
by = bx .* r1 .* (1 .- r1 .* r2);
V0 = (by .* ux) ./ (bx .* uy);
α = 1 .- V0;
βy = r1 .-(α.* by)/bx;
ey = bx .* (1 .- r3) .* (1 .- (r1 .* r2));

Y = zeros(Float64, length(ey));
Y[1] = 1.0;
X0 = 10.0;

# differential equations formula for finding densities of each pop at next time step
function fonction(u, p, t)
    x = u[1]
    y = u[2:end]
    regul = (1-(sum(u))/(p.K))
    dx = (p.bx*x + sum(p.ey.*y))*regul - p.ux*x-p.c*sum(p.βy.*y)*x
    dy = p.by .* y .* regul .- p.uy .* y .+ p.c .* p.βy .* x .* y
    return vcat(dx, dy)
end

debut = 0.0;
duree = 1000.0;
fin = debut + duree
N = zeros(Float64, (n_parasites+1, (n_parasites-1)*Int(duree)+1));
new_U = vcat(X0, Y);
parameters = (bx = bx, βy = βy, ey = ey, c = c, K = 80.0, ux = ux, by = by, uy = uy);

# each strain introduction (1000x)
@progress "Simulation" for i in 2:length(Y)
    # initial conditions
    prob = ODEProblem(fonction, new_U, (debut,fin), parameters)
    solution = solve(prob, saveat=debut:1.0:fin)
    for t in eachindex(solution.t)
        pop = solution.u[t]
        for i in 1:100
            if (pop.<0)[i]
                pop[i] = 0
            end
        end
        N[:,Int(solution.t[t]+1)] = pop
    end

    # set conditions & new parasite for next loop
    global new_U = solution[end]
    new_y = findfirst(x -> x == 0.0, new_U)
    new_U[new_y] = 1.0

    # set limits for next loop
    global debut = fin
    #global fin = Float64((i+1).* 1000.0)
    global fin = debut + duree
end

Np = N';

survival = (Np.>0.0)[:,2:end];
survived_ui = survival.*uy';
avg_survived = sum(survived_ui; dims=2)./sum(survival; dims=2);
avg_w_survived = sum(Np[:,2:end].*uy'; dims=2)./sum(Np[:,2:end]; dims=2);
uy_avg = avg_w_survived;

survived_by = survival.*by';
avg_by_survived = sum(survived_by; dims=2)./sum(survival; dims=2);
avg_w_by = sum(Np[:,2:end].*by'; dims=2)./sum(Np[:,2:end]; dims=2);
by_avg = avg_w_by;

#calculating R0

V0= by_avg  .*ux./(bx*uy_avg);
plot(V0,c=:black, title = "Vertical cases", xlabel = "Time", ylabel = "Mean V0", leg = false, ylims = (0.0,1.0))

# png("Figure 4/graph_4f.png")


text:
# Functions for the Vector Autoregressive model
import Pkg;
using ProgressMeter;
using LinearAlgebra, Statistics;
using Gadfly, Colors;

eye(n) = Matrix{Float64}(I, n, n)
chol(mat) = convert(Array{Float64}, cholesky(Hermitian(mat)).U');

# (1) Function to calculate impulse response
function func_VAR(data, p)
    # x: partial equilibrium data
    # p: lag order
    T, k = size(data)
    # Set up L.H.S.: set t-by-k array
    y = data';
    Y = y[:, p:T];

    for i = 1:(p-1)
        Y = [Y; y[:, (p-i):(T-i)]];
    end
    # Set up R.H.S.
    X = [ones(1, T-p); Y[:, 1:(T-p)]];
    Y = y[:, ((p+1):T)];

    # Least square estimation
    𝚩 = (Y*X')/(X*X');
    𝞄 = (Y - 𝚩 * X)';
    𝝨 = 𝞄' * 𝞄 ./ (T-p-p*k-1);
    return 𝚩, 𝞄, 𝝨;
end

function func_IRFvar(data, p, h, ℏ)
    T, k = size(data);

    𝚩_ols, 𝞄_ols, 𝝨_ols = func_VAR(data, p);
    V_ols = 𝚩_ols[:, 1];
    𝚩_ols = 𝚩_ols[:, 2:size(𝚩_ols,2)];
    # For standardization of the impulse response,
    ρ = convert(Array{Float64}, cholesky(𝝨_ols).U');

    𝝖 = [𝚩_ols; [kron(eye(k), eye(p - 1)) zeros(k*(p-1), k)]];
    J = [eye(k) zeros(k, k*(p-1))];

    # Compute impulse response
    ψ = zeros(k, k, (h + 1));

    for i_irf = 0:1:h
        ψ[:,:, (i_irf + 1)] = J * 𝝖^(i_irf) * J';
        # If needed to be standardized, multiply ρ
    end

    # Include recursive-design wild bootstrap calculation of confidence band
    #   Follow Kilian (2009, AER)
    #srand(8301993);
    n_repli = 2000; # Number of iterations

    ψ_boot = zeros(k, k, (h + 1), n_repli);

    y = data';
    Y = y[:, p:T];
    for i = 1:(p-1)
        Y = [Y; y[:, (p-i):(T-i)]];
    end

    Ur = zeros(k*p, T-p+1);
    Yr = zeros(k*p, T-p+1);

    @showprogress 1 "Wild-bootstrap confidence band " for i_iter = 1:n_repli
        pos = Int64(round(rand(1)[]*(T-p+1)));
        if pos == 0 #randomization hits zero, default to 1
            pos = 1
        end
        Yr[:, 1] = Y[:, pos];

        # Recrusive design bootstrap
        η = randn(1, size(𝞄_ols, 1));
        η = kron(η, ones(k, 1));

        Ur[1:k, 2:(T-p+1)] = 𝞄_ols' .* η;

        for ii_iter = 2:(T-p+1)
            Yr[:, ii_iter] = [V_ols; zeros(k*(p-1),1)] +
                                    𝝖 * Yr[:, (ii_iter - 1)] +
                                    Ur[:, ii_iter];
        end

        yr = Yr[1:k, :];
        for i=2:p
    		yr = [Yr[(((i-1)*k)+1):(i*k), 1] yr];
        end;
        yr = yr';

        # Iterated partial equilibrium estimation.;
        𝚩, 𝞄, 𝝨 = func_VAR(yr, p);

        𝚩 = 𝚩[:, 2:size(𝚩,2)];
        # For standardization of the impulse response,
        ρ = cholesky(𝝨).L[:,:]; # = cholesky()';

        𝝖 = [𝚩; [kron(eye(k), eye(p - 1)) zeros(k*(p-1), k)]];
        J = [eye(k) zeros(k, k*(p-1))];
        # Computaiton of impulse response
        temp_ψ_prev = [];
        for i_irf = 0:1:h
            temp_ψ = J * 𝝖^(i_irf) * J';
            ψ_boot[:,:, (i_irf + 1), i_iter] = temp_ψ ;
            # If needed to be standardized, multiply ρ
        end
    end
    # Need to calculate percentile
    α = 1; # One standard deviaiton confidence band of bootstrap iterations
    # Thus, this measure of confidence band does not carry literal definition
    #   as in classical model.

    # Calculate the standard deviation, across differen iteration (4th dim)
    σ_ψ_boot = std(ψ_boot, dims = 4);

    ψ_ub_1sd = ψ + σ_ψ_boot
    ψ_lb_1sd = ψ - σ_ψ_boot
    ψ_ub_2sd = ψ + σ_ψ_boot * 2;
    ψ_lb_2sd = ψ - σ_ψ_boot * 2;

    ############################################################################
    # FORECAST ERROR VARIANCE DECOMPOSITION
    T, k = size(data);

    𝚩_ols, 𝞄_ols, 𝝨_ols = func_VAR(data, p);
    V_ols = 𝚩_ols[:, 1];
    𝚩_ols = 𝚩_ols[:, 2:size(𝚩_ols,2)];
    # For standardization of the impulse response,
    ρ = convert(Array{Float64}, cholesky(𝝨_ols).U');

    𝝖 = [𝚩_ols; [kron(eye(k), eye(p - 1)) zeros(k*(p-1), k)]];
    J = [eye(k) zeros(k, k*(p-1))];

    # forecast error variance decomposition
    TH1 = J * 𝝖^(0) *J';
    TH = TH1 * ρ;
    TH = TH';
    TH2 = (TH .* TH);
    TH3 = TH2;

    for i=2:ℏ
        TH = J * 𝝖^(i-1) * J' * ρ;
        TH = TH'; TH2 = (TH.*TH);
        TH3 = TH3 + TH2;
    end;
    TH4 = sum(TH3, dims = 1); VC = zeros(k, k);
    for j=1:k
        VC[j,:] = TH3[j, :] ./ TH4';
    end;

    FEVDC = VC'*100;

    return ψ,
           ψ_lb_1sd, ψ_ub_1sd, ψ_lb_2sd, ψ_ub_2sd,
           FEVDC;
end

function func_contemp(bmat, 𝝨_ols, k, p)
    Ik = eye(k);
    i_block = 1:k; # initial, block of coefficients for same lag.
    for i = 1:p
        Ik = Ik .- bmat[:, i_block];
        # update the block index.
        i_block = i_block .+ k;
    end
    Ik = inv(Ik);

    Bk = eye(k);
    i_block = 1:k; # initial, block of coefficients for same lag.
    for i = 1:p
        Bk = Bk .- bmat[:, i_block]';
        # update the block index.
        i_block = i_block .+ k;
    end
    Bk = inv(Bk);

    # Long run impact matrix
    ψ0 = Ik * 𝝨_ols * Bk;
    #return Ik, Bk, ψ0
    ψ0 = chol(ψ0);

    # Contemporaneous shock identificaiton.
    ρ = Ik\ψ0';
    return ρ
end

function func_IRFvar_LR(data, p, h)
    # Implementation of Blanchard and Quah (1989)
    T, k = size(data);

    𝚩_ols, 𝞄_ols, 𝝨_ols = func_VAR(data, p);
    V_ols = 𝚩_ols[:, 1];
    𝚩_ols = 𝚩_ols[:, 2:size(𝚩_ols,2)];

    𝝖 = [𝚩_ols; [kron(eye(k), eye(p - 1)) zeros(k*(p-1), k)]];
    J = [eye(k) zeros(k, k*(p-1))];

    # For contemporaneous impact matrix - long run identificaiton
    # Contemporaneous shock identificaiton.
    ρ = func_contemp(𝚩_ols, 𝝨_ols, k, p);

    # Compute impulse response
    ψ = zeros(k, k, (h + 1));

    for i_irf = 0:1:h
        ψ[:,:, (i_irf + 1)] = J * 𝝖^(i_irf) * J' * ρ;
        # If needed to be standardized, multiply ρ
    end

    # Include recursive-design wild bootstrap calculation of confidence band
    #   Follow Kilian (2009, AER)
    #srand(8301993);
    n_repli = 2000; # Number of iterations

    ψ_boot = zeros(k, k, (h + 1), n_repli);

    y = data';
    Y = y[:, p:T];
    for i = 1:(p-1)
        Y = [Y; y[:, (p-i):(T-i)]];
    end

    Ur = zeros(k*p, T-p+1);
    Yr = zeros(k*p, T-p+1);

    @showprogress 1 "Wild-bootstrap confidence band " for i_iter = 1:n_repli
        pos = Int64(round(rand(1)[]*(T-p+1)));
        if pos == 0 #randomization hits zero, default to 1
            pos = 1
        end
        Yr[:, 1] = Y[:, pos];

        # Recrusive design bootstrap
        η = randn(1, size(𝞄_ols, 1));
        η = kron(η, ones(k, 1));

        Ur[1:k, 2:(T-p+1)] = 𝞄_ols' .* η;

        for ii_iter = 2:(T-p+1)
            Yr[:, ii_iter] = [V_ols; zeros(k*(p-1),1)] +
                                    𝝖 * Yr[:, (ii_iter - 1)] +
                                    Ur[:, ii_iter];
        end

        yr = Yr[1:k, :];
        for i=2:p
    		yr = [Yr[(((i-1)*k)+1):(i*k), 1] yr];
        end;
        yr = yr';

        # Iterated partial equilibrium estimation.;
        𝚩, 𝞄, 𝝨 = func_VAR(yr, p);

        𝚩 = 𝚩[:, 2:size(𝚩,2)];
        # For standardization of the impulse response,
        # For contemporaneous impact matrix - long run identificaiton
        # Contemporaneous shock identificaiton.
        ρ = func_contemp(𝚩, 𝝨, k, p);

        𝝖 = [𝚩; [kron(eye(k), eye(p - 1)) zeros(k*(p-1), k)]];
        J = [eye(k) zeros(k, k*(p-1))];
        # Computaiton of impulse response
        temp_ψ_prev = [];
        for i_irf = 0:1:h
            temp_ψ = J * 𝝖^(i_irf) * J' * ρ;
            ψ_boot[:,:, (i_irf + 1), i_iter] = temp_ψ ;
            # If needed to be standardized, multiply ρ
        end
    end
    # Need to calculate percentile
    α = 1; # One standard deviaiton confidence band of bootstrap iterations
    # Thus, this measure of confidence band does not carry literal definition
    #   as in classical model.

    # Calculate the standard deviation, across differen iteration (4th dim)
    σ_ψ_boot = std(ψ_boot, dims = 4);

    ψ_ub_1sd = ψ + σ_ψ_boot
    ψ_lb_1sd = ψ - σ_ψ_boot
    ψ_ub_2sd = ψ + σ_ψ_boot * 2;
    ψ_lb_2sd = ψ - σ_ψ_boot * 2;

    return ψ,
           ψ_lb_1sd, ψ_ub_1sd, ψ_lb_2sd, ψ_ub_2sd
end


text:
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../test/useful_routines.jl")

function cell_interface_error(quad, cellmap, facescale, distance_to_interface)
    err = 0.0
    for (idx, (p, w)) in enumerate(quad)
        x = cellmap(p)
        err += (distance_to_interface(x)[1])^2 * facescale[idx] * w
    end
    return err
end

function interface_error(levelset, levelsetcoeffs, cutmesh, numqp, distance_to_interface)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    cellids = findall(CutCell.cell_sign(cutmesh) .== 0)

    err = 0.0
    for cellid in cellids
        quad = interfacequads[+1, cellid]
        if length(quad) > 0
            normals = CutCell.interface_normals(interfacequads, cellid)
            cellmap = CutCell.cell_map(cutmesh, cellid)
            facescale = CutCell.scale_area(cellmap, normals)

            err += cell_interface_error(quad, cellmap, facescale, distance_to_interface)
        end
    end
    return sqrt(err)
end

function shrinking_circle_distance(x, xc, initialradius, speed, time)
    r = initialradius - speed * time
    return circle_distance_function(x, xc, r)
end

function interface_error_over_timesteps(
    levelset,
    mesh,
    speed,
    dt,
    nsteps,
    numqp,
    xc,
    initialradius,
)
    err = zeros(nsteps + 1)

    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> shrinking_circle_distance(x, xc, initialradius, speed, 0.0),
        mesh,
    )
    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)

    err[1] = interface_error(
        levelset,
        levelsetcoeffs,
        cutmesh,
        numqp,
        x -> shrinking_circle_distance(x, xc, initialradius, speed, 0.0),
    )/(2pi*initialradius)

    levelsetspeed = speed * ones(length(levelsetcoeffs))

    for t = 1:nsteps
        paddedmesh = CutCell.BoundaryPaddedMesh(cutmesh, 1)
        refseedpoints, spatialseedpoints, seedcellids =
            CutCell.seed_zero_levelset(2, levelset, levelsetcoeffs, cutmesh)
        paddedlevelset = CutCell.BoundaryPaddedLevelset(
            paddedmesh,
            refseedpoints,
            spatialseedpoints,
            seedcellids,
            levelset,
            levelsetcoeffs,
            cutmesh,
            1e-10,
        )
        levelsetcoeffs =
            CutCell.step_first_order_levelset(paddedlevelset, levelsetspeed, dt)

        cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
        err[t+1] = interface_error(
            levelset,
            levelsetcoeffs,
            cutmesh,
            numqp,
            x -> shrinking_circle_distance(x,xc,initialradius,speed,t*dt),
        )/(2pi*(initialradius - speed*t*dt))
    end
    return err
end

function final_interface_error(
    levelset,
    mesh,
    speed,
    dt,
    nsteps,
    numqp,
    xc,
    initialradius,
)

    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> shrinking_circle_distance(x, xc, initialradius, speed, 0.0),
        mesh,
    )
    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    levelsetspeed = speed * ones(length(levelsetcoeffs))

    for t = 1:nsteps
        paddedmesh = CutCell.BoundaryPaddedMesh(cutmesh, 1)
        refseedpoints, spatialseedpoints, seedcellids =
            CutCell.seed_zero_levelset(2, levelset, levelsetcoeffs, cutmesh)
        paddedlevelset = CutCell.BoundaryPaddedLevelset(
            paddedmesh,
            refseedpoints,
            spatialseedpoints,
            seedcellids,
            levelset,
            levelsetcoeffs,
            cutmesh,
            1e-10,
        )
        levelsetcoeffs =
            CutCell.step_first_order_levelset(paddedlevelset, levelsetspeed, dt)

        cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    end
    err = interface_error(
        levelset,
        levelsetcoeffs,
        cutmesh,
        numqp,
        x -> shrinking_circle_distance(x,xc,initialradius,speed,nsteps*dt),
    )/(2pi*(initialradius - speed*nsteps*dt))
    return err
end

function grid_size(mesh)
    w = CutCell.widths(mesh)
    nn = CutCell.nodes_per_mesh_side(mesh)
    return w ./ (nn .- 1)
end

function time_step_size(speed, mesh; CFL = 0.5)
    dx = minimum(grid_size(mesh))
    return CFL * dx / abs(speed)
end

function numsteps(stoptime,dt)
    r = stoptime/dt
    @assert isinteger(r)
    return round(Int,r)
end

function error_for_numelmts(nelmts)
    x0 = [0.0, 0.0]
    L, W = 1.0, 1.0
    numghostlayers = 1
    polyorder = 2
    numqp = 3

    xc = [0.5, 0.5]
    initialradius = 0.50
    speed = 1.0
    stoptime = 0.4

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)

    mesh = CutCell.Mesh(x0, [L, W], [nelmts, nelmts], basis)
    dt = time_step_size(speed, mesh)
    nsteps = numsteps(stoptime,dt)

    return final_interface_error(
        levelset,
        mesh,
        speed,
        dt,
        nsteps,
        numqp,
        xc,
        initialradius,
    )
end

nelmts = [5,10,20,40,80]

# err80 = error_for_numelmts(80)
err = [error_for_numelmts(ne) for ne in nelmts]

dx = 1.0 ./ nelmts
rate = diff(log.(err)) ./ diff(log.(dx))

using CSV,DataFrames
df = DataFrame([dx,err],["Element Size", "Error"])
CSV.write("examples/levelset-propagate/circle-convergence.csv",df)

#
# using Plots
# fig = plot(legend=:topleft)
# plot!(range(0,stop=0.4,length=length(err[1])),err[1],label="5x5")
# plot!(range(0,stop=0.4,length=length(err[2])),err[2],label="10x10")
# plot!(range(0,stop=0.4,length=length(err[3])),err[3],label="20x20")
# plot!(range(0,stop=0.4,length=length(err[4])),err[4],label="40x40")


text:
module Triangles

export Triangle, transform, hittriangle

using CthulhuVision.Math
using CthulhuVision.Light
using CthulhuVision.Materials
using CUDAnative

struct Triangle
    a::Vector3
    b::Vector3
    c::Vector3
    material::Material
end

transform(tri::Triangle, form::Transform) :: Triangle = Triangle(form * tri.a, form * tri.b, form * tri.c, tri.material)

##########
# Device #
##########

@inline function maxdim(v::Vector3)
    if v.x > v.y
        if v.x > v.z
            0
        else
            2
        end
    else
        if v.y > v.z
            1
        else
            2
        end
    end
end

@inline function hittriangle(tri::Triangle, tmin::Float32, tmax::Float32, ray::Ray) :: HitRecord
    v0v1 = tri.b - tri.a
    v0v2 = tri.c - tri.a

    direct = unit(direction(ray))

    normal = unit(cross(v0v1, v0v2))
    nr = dot(normal, direct)

    # Very nearly parallel, so don't report a hit.
    if abs(nr) < 0.0001f0
        return HitRecord()
    end

    no = dot(normal, origin(ray))
    d = -dot(normal, tri.a)
    t = -(no + d) / nr

    if t < tmin || t > tmax
        return HitRecord()
    end

    p = origin(ray) + t * direct

    # Determine if p is inside the triangle.
    edge0 = tri.b - tri.a
    v0p = p - tri.a
    e0 = cross(edge0, v0p)

    edge1 = tri.c - tri.b
    v1p = p - tri.b
    e1 = cross(edge1, v1p)

    edge2 = tri.a - tri.c
    v2p = p - tri.c
    e2 = cross(edge2, v2p)

    # The hit point p is on the wrong side of one of the edges.
    if dot(e0, normal) < 0f0 || dot(e1, normal) < 0f0 || dot(e2, normal) < 0f0
        return HitRecord()
    end

    if nr > 0f0
        normal = -normal
    end

    HitRecord(t, p, normal, tri.material)
end

end


text:
@testset "TopologicalRelation" begin
  @testset "HalfEdgeTopology" begin
    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    struc = HalfEdgeTopology(elems)
    ∂ = Boundary{2,0}(struc)
    @test ∂(1) == [2,3,1]
    @test ∂(2) == [3,2,4]
    ∂ = Boundary{2,1}(struc)
    @test ∂(1) == [1,3,2]
    @test ∂(2) == [1,4,5]
    ∂ = Boundary{1,0}(struc)
    @test ∂(1) == [3,2]
    @test ∂(2) == [1,2]
    @test ∂(3) == [3,1]
    @test ∂(4) == [2,4]
    @test ∂(5) == [4,3]
    𝒞 = Coboundary{0,1}(struc)
    @test 𝒞(1) == [2,3]
    @test 𝒞(2) == [4,1,2]
    @test 𝒞(3) == [3,1,5]
    @test 𝒞(4) == [5,4]
    𝒞 = Coboundary{0,2}(struc)
    @test 𝒞(1) == [1]
    @test 𝒞(2) == [2,1]
    @test 𝒞(3) == [1,2]
    @test 𝒞(4) == [2]
    𝒞 = Coboundary{1,2}(struc)
    @test 𝒞(1) == [2,1]
    @test 𝒞(2) == [1]
    @test 𝒞(3) == [1]
    @test 𝒞(4) == [2]
    @test 𝒞(5) == [2]
    𝒜 = Adjacency{0}(struc)
    @test 𝒜(1) == [2,3]
    @test 𝒜(2) == [4,3,1]
    @test 𝒜(3) == [1,2,4]
    @test 𝒜(4) == [3,2]

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    struc = HalfEdgeTopology(elems)
    ∂ = Boundary{2,0}(struc)
    @test ∂(1) == [1,2,6,5]
    @test ∂(2) == [6,2,4]
    @test ∂(3) == [6,4,3,5]
    @test ∂(4) == [3,1,5]
    ∂ = Boundary{2,1}(struc)
    @test ∂(1) == [1,3,5,6]
    @test ∂(2) == [3,9,4]
    @test ∂(3) == [4,7,8,5]
    @test ∂(4) == [2,6,8]
    ∂ = Boundary{1,0}(struc)
    @test ∂(1) == [1,2]
    @test ∂(2) == [3,1]
    @test ∂(3) == [6,2]
    @test ∂(4) == [4,6]
    @test ∂(5) == [5,6]
    @test ∂(6) == [1,5]
    @test ∂(7) == [4,3]
    @test ∂(8) == [3,5]
    @test ∂(9) == [2,4]
    𝒞 = Coboundary{0,1}(struc)
    @test 𝒞(1) == [1,6,2]
    @test 𝒞(2) == [9,3,1]
    @test 𝒞(3) == [2,8,7]
    @test 𝒞(4) == [7,4,9]
    @test 𝒞(5) == [5,8,6]
    @test 𝒞(6) == [3,4,5]
    𝒞 = Coboundary{0,2}(struc)
    @test 𝒞(1) == [1,4]
    @test 𝒞(2) == [2,1]
    @test 𝒞(3) == [4,3]
    @test 𝒞(4) == [3,2]
    @test 𝒞(5) == [3,4,1]
    @test 𝒞(6) == [2,3,1]
    𝒞 = Coboundary{1,2}(struc)
    @test 𝒞(1) == [1]
    @test 𝒞(2) == [4]
    @test 𝒞(3) == [2,1]
    @test 𝒞(4) == [2,3]
    @test 𝒞(5) == [3,1]
    @test 𝒞(6) == [4,1]
    @test 𝒞(7) == [3]
    @test 𝒞(8) == [3,4]
    @test 𝒞(9) == [2]
    𝒜 = Adjacency{0}(struc)
    @test 𝒜(1) == [2,5,3]
    @test 𝒜(2) == [4,6,1]
    @test 𝒜(3) == [1,5,4]
    @test 𝒜(4) == [3,6,2]
    @test 𝒜(5) == [6,3,1]
    @test 𝒜(6) == [2,4,5]
  end
end


text:
module MapperCore

using ..Mapper2.Helper
Helper.@SetupDocStringTemplates

using ..Mapper2.MapperGraphs

using IterTools
using DataStructures

using Logging
using Serialization

"""
    RuleSet

Abstract supertype for controlling dispatch to specialized functions for
architecture interpretation. Create a custom concrete subtype of this if you
want to use custom methods during placement or routing.
"""
abstract type RuleSet end

struct DefaultRuleSet <: RuleSet end

include("Taskgraphs.jl")

# Architecture includes
include("Architecture/Paths.jl")
include("Architecture/Architecture.jl")
include("Architecture/Methods.jl")
include("Architecture/Constructors.jl")

# Map Includes
include("Map/Map.jl")
include("Map/Verification.jl")
include("Map/Inspection.jl")

#-------------------------------------------------------------------------------
# Mapping methods
#-------------------------------------------------------------------------------

const TN = TaskgraphNode
const TE = TaskgraphEdge
const PLC = Union{Port,Link,Component}

# Placement Queries
"""
    isequivalent(ruleset::RuleSet, a::TaskgraphNode, b::TaskgraphNode) :: Bool

Return `true` if `TaskgraphNodes` `a` and `b` are semantically equivalent for
placement.

Default: `true`
"""
isequivalent(::RuleSet, a::TN, b::TN) = true

"""
    ismappable(ruleset::RuleSet, component::Component) :: Bool

Return `true` if some task can be mapped to `component` under `ruleset`.

Default: `true`
"""
ismappable(::RuleSet, c::Component) = true

"""
    canmap(ruleset::RuleSet, t::TaskgraphNode, c::Component) :: Bool

Return `true` if `t` can be mapped to `c` under `ruleset`.

Default: `true`
"""
canmap(::RuleSet, t::TN, c::Component) = true

# Routing Queries

"""
    canuse(ruleset::RuleSet, item::Union{Port,Link,Component}, edge::TaskgraphEdge)::Bool

Return `true` if `edge` can use `item` as a routing resource under `ruleset`.

Default: `true`
"""
canuse(::RuleSet, item::PLC, edge::TE) = true

"""
    getcapacity(ruleset::RuleSet, item::Union{Port,Link,Component})

Return the capacity of routing resource `item` under `ruleset`.

Default: `1`
"""
getcapacity(::RuleSet, item) = 1

"""
    is_source_port(ruleset::RuleSet, port::Port, edge::TaskgraphEdge)::Bool

Return `true` if `port` is a valid source port for `edge` under `ruleset`.

Default: `true`
"""
is_source_port(::RuleSet, port::Port, edge::TE) = true

"""
    is_sink_port(ruleset::RuleSet, port::Port, edge::TaskgraphEdge)::Bool

Return `true` if `port` is a vlid sink port for `edge` under `ruleset`.

Default: `true`
"""
is_sink_port(::RuleSet, port::Port, edge::TE) = true

"""
    needsrouting(ruleset::RuleSet, edge::TaskgraphEdge)::Bool

Return `true` if `edge` needs to be routed under `ruleset`.

Default: `true`
"""
needsrouting(::RuleSet, edge::TE) = true

#-------------------------------------------------------------------------------
# Exports
#-------------------------------------------------------------------------------

export isequivalent,
    ismappable, canmap, canuse, is_sink_port, is_source_port, getcapacity, needsrouting

### Taskgraph Exports###

export TaskgraphNode,
    TaskgraphEdge,
    Taskgraph,
    # Methods
    getsources,
    getsinks,
    getnodes,
    getedges,
    getnode,
    getedge,
    nodenames,
    num_nodes,
    num_edges,
    add_node,
    add_edge,
    out_edges,
    out_edge_indices,
    in_edges,
    in_edge_indices,
    hasnode,
    outnode_names,
    innode_names,
    out_nodes,
    in_nodes,

    ### Architecture Exports ###
    # Path Types
    AbstractPath,
    Path,
    AddressPath,
    catpath,
    striplast,
    splitpath,
    stripfirst,

    # Architecture stuff
    Port,
    Link,
    # Port Methods
    Input,
    Output,
    invert,
    # Link Methods
    isaddresslink,
    sources,
    dests,
    # Components
    AbstractComponent,
    TopLevel,
    Component,
    isaddress,
    @port_str,
    @link_str,
    @component_str,
    # Verification
    check_routing,
    # Methods
    getoffset,
    checkclass,
    ports,
    portpaths,
    addresses,
    pathtype,
    children,
    walk_children,
    connected_components,
    search_metadata,
    search_metadata!,
    get_metadata!,
    check_connectivity,
    get_connected_port,
    isfree,
    isgloballink,
    isglobalport,
    getaddress,
    hasaddress,
    getchild,
    getname,
    mappables,
    # Asserts
    assert_no_children,
    assert_no_intrarouting,
    # Constructor Types
    OffsetRule,
    PortRule,
    # Constructor Functions
    add_port,
    add_child,
    add_link,
    connection_rule,
    build_mux,
    check,
    ConnectionRule,
    Offset,
    # Analysis methods
    build_distance,
    build_component_table,

    ### Map ###
    RuleSet,
    Mapping,
    Map,
    NodeMap,
    EdgeMap,
    rules,
    save,
    load,
    getpath
end


text:
abstract type _AbstractFunction end

"""
    _Operation

An enum representing different node types in a symbolic expression graph.

 * `_OP_COEFFICIENT`: a numerical `Float64`
 * `_OP_VARIABLE`: a decision variable
 * `_OP_OPERATION`: a function call
 * `_OP_INTEGER_EXPONENT`: a special case for handling integer polynomials such
   as `x^2`, since the `2.0` is usually structural rather than a coefficient
   that differs between expressions.
"""
@enum(
    _Operation,
    _OP_COEFFICIENT,
    _OP_VARIABLE,
    _OP_OPERATION,
    _OP_INTEGER_EXPONENT,
)

"""
    _Node

A type to represent a node in the symbolic expression graph.

!!! warning
    `_Node`s can only be used in conjunction with `_Function`.

`_Node` has the following fields:

 * `.head::_Operation`: this describes the node type
 * `.hash::UInt64`: a hash of the expression tree rooted at this node. See the
   Fast hashing section for more details.
 * `.index::Int` : interpretation depends on `.head`. If head is:
   * `_OP_COEFFICIENT`, the index into `f.data`
   * `_OP_VARIABLE`, the index into `f.ordered_variables`
   * `_OP_INTEGER_EXPONENT`, the integer exponent `children[1]^index`
 * `.operation::Symbol` : the symbolic name of the function call, if `.head` is
   `_OP_OPERATION`, otherwise this can be ignored.
 * `.children::Union{Nothing,Vector{_Node}}` : a vector of children of the node,
   if `.head` is `_OP_OPERATION`, otherwise `nothing`.

## Fast hashing

The expression graph is structured as a Merkle tree (a.k.a. a hash tree). That
is, each node contains a field `.hash`, which is the hash of the node's
operation with the hash of its children. Since children are also nodes, hashing
a child is an O(1) lookup of `child.hash`. Moreover, two expressions are
structurally identical if the have the same hash when converted into `_Node`
form.

## Other considerations

Through experience, the representation of the symbolic expression graph is not
the bottleneck in a nonlinear optimization model. Therefore, we choose the
somewhat inefficient storage rpresentation of
`children::Union{Nothing,Vector{_Node}}` for the list of children. This results
in a GC tracked object for each node that has children. If memory issues are
identified in future, changing how we represent the expression graph is one
promising avenue of investigation.
"""
struct _Node
    head::_Operation
    hash::UInt
    index::Int
    operation::Symbol
    children::Union{Nothing,Vector{_Node}}

    function _Node(f::_AbstractFunction, coefficient::Float64)
        # Store the coefficient
        push!(f.data, coefficient)
        # Compute the hash of the node. Because all coefficients are alike, this
        # is just the hash of `_OP_COEFFICIENT`.
        h = hash(_OP_COEFFICIENT)
        # Index is the length of `f.data` when we need to look it up again.
        index = length(f.data)
        return new(_OP_COEFFICIENT, h, index, :NONE, nothing)
    end
    function _Node(f::_AbstractFunction, variable::MOI.VariableIndex)
        # We need to convert the variable to the order in which it will appear
        # in .ordered_variables. Currently, the order is stored in
        # `f.variables`...
        lookup = f.variables::Dict{MOI.VariableIndex,Int}
        # And if the variable hasn't been seen yet, add a new index
        index = get!(lookup, variable, length(lookup) + 1)
        # The hash of this node depends on the variable and it's position in
        # ordered_variables.
        h = hash(_OP_VARIABLE, hash(index))
        return new(_OP_VARIABLE, h, index, :NONE, nothing)
    end
    function _Node(::_AbstractFunction, operation::Symbol, children::_Node...)
        # Okay, a function call. The hash depends on the operation and the
        # children.
        h = hash(operation, _hash(children...))
        return new(_OP_OPERATION, h, 0, operation, collect(children))
    end
    # Special cased performance optimization: x^N where N::Int is a common
    # operation that would otherwise get re-written to x[i]^p[j]. Doing so
    # disables a lot of potential optimizations, so it's worth special-casing
    # this. Other operations aren't as important.
    function _Node(::_AbstractFunction, ::typeof(^), x::_Node, N::Int)
        # The hash depends on ^, as well as the child node `x` and the
        # exponent `N`
        h = hash(_OP_INTEGER_EXPONENT, hash(_hash(x), hash(N)))
        return new(_OP_INTEGER_EXPONENT, h, N, :^, [x])
    end
end

# The hash of a child is an O(1) lookup
_hash(child::_Node) = child.hash

# Recursively has nodes together
_hash(a::_Node, args::_Node...) = hash(a.hash, _hash(args...))

"""
    _Function(expr::Expr)

This function parses the Julia expression `expr` into a `_Function`.

`_Function` has the following fields:

 * `ordered_variables`: a vector containing the `.value` indices of each
   variable in the order that they are discovered in the expression tree.
 * `coefficients::Vector{Float64}`: a vector to store the numerical input values
   of `x`.
 * `data::Vector{Float64}`: the list of numerical coefficients that appear in
   `expr`
 * `expr::_Node`: the symbolic representation of `expr`.

For more information on how to interpret the fields, see `_Node`.
"""
mutable struct _Function <: _AbstractFunction
    ordered_variables::Vector{Int}
    coefficients::Vector{Float64}
    data::Vector{Float64}
    expr::_Node
    # A field that is only used when constructing the _Function
    variables::Union{Nothing,Dict{MOI.VariableIndex,Int}}
    function _Function(expr::Expr)
        f = new()
        f.variables = Dict{MOI.VariableIndex,Int}()
        f.data = Float64[]
        f.expr = _Node(f, expr)
        f.coefficients = zeros(length(f.variables))
        f.ordered_variables = zeros(Int, length(f.variables))
        for (k, v) in f.variables
            f.ordered_variables[v] = k.value
        end
        f.variables = nothing
        return f
    end
end

_is_integer(::Any) = false
_is_integer(x::Real) = Base.isinteger(x)

"""
    _Node(f::_Function, expr::Expr)

Parse `expr` into `f` and return a `_Node`.

!!! warning
    This function gets called recursively.
"""
function _Node(f::_Function, expr::Expr)
    if isexpr(expr, :call)
        # Performance optimization: most calls will be unary or binary
        # operators. Therefore, we can specialize an if-statement to handle the
        # common cases without needing to splat.
        if length(expr.args) == 2
            return _Node(f, expr.args[1], _Node(f, expr.args[2]))
        elseif length(expr.args) == 3
            # Special cased performance optimization for _OP_INTEGER_EXPONENT
            if expr.args[1] == :^ && _is_integer(expr.args[3])
                return _Node(f, ^, _Node(f, expr.args[2]), Int(expr.args[3]))
            end
            return _Node(
                f,
                expr.args[1],
                _Node(f, expr.args[2]),
                _Node(f, expr.args[3]),
            )
        else
            nodes = (_Node(f, expr.args[i]) for i in 2:length(expr.args))
            return _Node(f, expr.args[1], nodes...)
        end
    else
        @assert isexpr(expr, :ref)
        @assert expr.args[1] == :x
        return _Node(f, expr.args[2])
    end
end

struct _SymbolicsFunction{F,G,H}
    f::F
    ∇f::G
    ∇²f::H
    g::Vector{Float64}
    H::Vector{Float64}
    H_structure::Vector{Tuple{Int,Int}}
end

"""
    _expr_to_symbolics(expr::_Node, p, x)

Convert a `_Node` into a `Symbolics.jl` expression via operator overloading.
"""
function _expr_to_symbolics(expr::_Node, p, x)
    if expr.head == _OP_COEFFICIENT
        return p[expr.index]
    elseif expr.head == _OP_VARIABLE
        return x[expr.index]
    elseif expr.head == _OP_OPERATION
        # TODO(odow): improve this lookup of expr.operation to function.
        f = getfield(Base, expr.operation)
        return f((_expr_to_symbolics(c, p, x) for c in expr.children)...)
    else
        @assert expr.head == _OP_INTEGER_EXPONENT
        return _expr_to_symbolics(expr.children[1], p, x)^expr.index
    end
end

function _SymbolicsFunction(f::_Function, features::Vector{Symbol})
    d, n = length(f.data), length(f.ordered_variables)
    Symbolics.@variables(p[1:d], x[1:n])
    p, x = collect(p), collect(x)
    f_expr = _expr_to_symbolics(f.expr, p, x)
    f = Symbolics.build_function(f_expr, p, x; expression = Val{false})
    if :Jac in features || :Grad in features
        ∇f_expr = Symbolics.gradient(f_expr, x)
        _, g! = Symbolics.build_function(∇f_expr, p, x; expression = Val{false})
        g_cache = zeros(length(x))
    else
        g!, g_cache = nothing, Float64[]
    end
    if :Hess in features
        ∇²f_expr_square = Symbolics.sparsejacobian(∇f_expr, x)
        # ∇²f_expr_square is sparse but it is also symmetrical. MathOptInterface
        # needs only the lower-triangular (technically, it doesn't matter which
        # triangle, but let's choose the lower-triangular for simplicity).
        I, J, V = SparseArrays.findnz(∇²f_expr_square)
        # Rather than using a SparseArray as storage, convert to MOI's
        # ((i, j), v) datastructure. This simplifies
        # `hessian_lagrangian_structure` calls later.
        ∇²f_expr = [V[i] for i in 1:length(I) if I[i] >= J[i]]
        ∇²f_structure = [(I[i], J[i]) for i in 1:length(I) if I[i] >= J[i]]
        _, h! =
            Symbolics.build_function(∇²f_expr, p, x; expression = Val{false})
        h_cache = zeros(length(∇²f_structure))
    else
        h!, h_cache, ∇²f_structure = nothing, Float64[], Tuple{Int,Int}[]
    end
    return _SymbolicsFunction(f, g!, h!, g_cache, h_cache, ∇²f_structure)
end

"""
    AbstractNonlinearOracleBackend

An abstract type to help dispatch on different implementations of the function
evaluations.
"""
abstract type AbstractNonlinearOracleBackend end

"""
    DefaultBackend() <: AbstractNonlinearOracleBackend

A simple implementation that loops over the different constraints.
"""
struct DefaultBackend <: AbstractNonlinearOracleBackend end

mutable struct _NonlinearOracle{B} <: MOI.AbstractNLPEvaluator
    backend::B
    objective::Union{Nothing,_Function}
    constraints::Vector{_Function}
    symbolic_functions::Dict{UInt,_SymbolicsFunction}
    hessian_sparsity_map::Vector{Int}
    eval_objective_timer::Float64
    eval_objective_gradient_timer::Float64
    eval_constraint_timer::Float64
    eval_constraint_jacobian_timer::Float64
    eval_hessian_lagrangian_timer::Float64
    function _NonlinearOracle(backend::B, objective, constraints) where {B}
        return new{B}(
            backend,
            objective,
            constraints,
            Dict{UInt,_SymbolicsFunction}(),
            Int[],
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
        )
    end
end

function MOI.initialize(
    ::_NonlinearOracle,
    ::AbstractNonlinearOracleBackend,
    ::Vector{Symbol},
)
    return
end

function MOI.initialize(oracle::_NonlinearOracle, features::Vector{Symbol})
    oracle.eval_constraint_timer = 0.0
    oracle.eval_constraint_jacobian_timer = 0.0
    oracle.eval_hessian_lagrangian_timer = 0.0
    hashes = map(f -> f.expr.hash, oracle.constraints)
    for h in unique(hashes)
        f = oracle.constraints[findfirst(isequal(h), hashes)]
        oracle.symbolic_functions[h] = _SymbolicsFunction(f, features)
    end
    if oracle.objective !== nothing
        f = oracle.objective
        h = f.expr.hash
        if !haskey(oracle.symbolic_functions, h)
            oracle.symbolic_functions[h] = _SymbolicsFunction(f, features)
        end
    end
    MOI.initialize(oracle, oracle.backend, features)
    return
end

MOI.features_available(::_NonlinearOracle) = [:Grad, :Jac, :Hess]

function _update_coefficients(f::_Function, x::AbstractVector{Float64})
    for (i, v) in enumerate(f.ordered_variables)
        @inbounds f.coefficients[i] = x[v]
    end
    return
end

function _eval_function(
    oracle::_NonlinearOracle,
    func::_Function,
    x::AbstractVector{Float64},
)
    _update_coefficients(func, x)
    @inbounds f = oracle.symbolic_functions[func.expr.hash]
    return f.f(func.data, func.coefficients)::Float64
end

function _eval_gradient(
    oracle::_NonlinearOracle,
    func::_Function,
    x::AbstractVector{Float64},
)
    _update_coefficients(func, x)
    @inbounds f = oracle.symbolic_functions[func.expr.hash]
    f.∇f(f.g, func.data, func.coefficients)
    return f.g
end

function _eval_hessian(
    oracle::_NonlinearOracle,
    func::_Function,
    x::AbstractVector{Float64},
)
    _update_coefficients(func, x)
    @inbounds f = oracle.symbolic_functions[func.expr.hash]
    f.∇²f(f.H, func.data, func.coefficients)
    return f.H
end

function MOI.eval_objective(
    oracle::_NonlinearOracle,
    x::AbstractVector{Float64},
)
    start = time()
    @assert oracle.objective !== nothing
    objective_value = _eval_function(oracle, oracle.objective, x)
    oracle.eval_objective_timer += time() - start
    return objective_value
end

function MOI.eval_objective_gradient(
    oracle::_NonlinearOracle,
    g::AbstractVector{Float64},
    x::AbstractVector{Float64},
)
    start = time()
    @assert oracle.objective !== nothing
    ∇f = _eval_gradient(oracle, oracle.objective, x)
    fill!(g, 0.0)
    for (i, v) in enumerate(oracle.objective.ordered_variables)
        @inbounds g[v] = ∇f[i]
    end
    oracle.eval_objective_gradient_timer += time() - start
    return
end

function MOI.eval_constraint(
    oracle::_NonlinearOracle,
    g::AbstractVector{Float64},
    x::AbstractVector{Float64},
)
    start = time()
    for (row, c) in enumerate(oracle.constraints)
        g[row] = _eval_function(oracle, c, x)
    end
    oracle.eval_constraint_timer += time() - start
    return
end

function MOI.jacobian_structure(oracle::_NonlinearOracle)
    structure = Tuple{Int,Int}[]
    for (row, c) in enumerate(oracle.constraints)
        for col in c.ordered_variables
            push!(structure, (row, col))
        end
    end
    return structure
end

function MOI.eval_constraint_jacobian(
    oracle::_NonlinearOracle,
    J::AbstractVector{Float64},
    x::AbstractVector{Float64},
)
    start = time()
    k = 1
    for c in oracle.constraints
        g = _eval_gradient(oracle, c, x)
        for gi in g
            J[k] = gi
            k += 1
        end
    end
    oracle.eval_constraint_jacobian_timer += time() - start
    return
end

function _hessian_lagrangian_structure(
    ::Vector{Tuple{Int,Int}},
    ::_NonlinearOracle,
    ::Any,
    ::Nothing,
)
    return
end

function _hessian_lagrangian_structure(
    structure::Vector{Tuple{Int,Int}},
    oracle::_NonlinearOracle,
    map::Dict{Tuple{Int,Int},Int},
    c::_Function,
)
    f = oracle.symbolic_functions[c.expr.hash]
    for (i, j) in f.H_structure
        row, col = c.ordered_variables[i], c.ordered_variables[j]
        row_col = row >= col ? (row, col) : (col, row)
        index = get(map, row_col, nothing)
        if index === nothing
            push!(structure, row_col)
            map[row_col] = length(structure)
            push!(oracle.hessian_sparsity_map, length(structure))
        else
            push!(oracle.hessian_sparsity_map, index::Int)
        end
    end
    return
end

function MOI.hessian_lagrangian_structure(oracle::_NonlinearOracle)
    structure = Tuple{Int,Int}[]
    map = Dict{Tuple{Int,Int},Int}()
    _hessian_lagrangian_structure(structure, oracle, map, oracle.objective)
    for c in oracle.constraints
        _hessian_lagrangian_structure(structure, oracle, map, c)
    end
    return structure
end

function MOI.eval_hessian_lagrangian(
    oracle::_NonlinearOracle,
    H::AbstractVector{Float64},
    x::AbstractVector{Float64},
    σ::Float64,
    μ::AbstractVector{Float64},
)
    fill!(H, 0.0)
    start = time()
    k = 1
    if oracle.objective !== nothing && !iszero(σ)
        h = _eval_hessian(oracle, oracle.objective, x)
        for (i, hi) in enumerate(h)
            @inbounds H[oracle.hessian_sparsity_map[i]] += σ * hi
        end
        k += length(h)
    end
    for (μi, constraint) in zip(μ, oracle.constraints)
        if iszero(μi)
            continue
        end
        h = _eval_hessian(oracle, constraint, x)
        for nzval in h
            @inbounds H[oracle.hessian_sparsity_map[k]] += μi * nzval
            k += 1
        end
    end
    oracle.eval_hessian_lagrangian_timer += time() - start
    return
end

function _nonlinear_constraint(expr::Expr)
    lower, upper, body = -Inf, Inf, :()
    if isexpr(expr, :call, 3)
        if expr.args[1] == :(>=)
            lower, upper, body = expr.args[1], expr.args[5], expr.args[3]
        elseif expr.args[1] == :(<=)
            lower, upper, body = -Inf, expr.args[3], expr.args[2]
        else
            @assert expr.args[1] == :(==)
            lower, upper, body = expr.args[3], expr.args[3], expr.args[2]
        end
    else
        @assert isexpr(expr, :comparison, 5)
        lower, upper, body = expr.args[1], expr.args[5], expr.args[3]
    end
    return _Function(body), MOI.NLPBoundsPair(lower, upper)
end

"""
    _nlp_block_data(
        d::MOI.AbstractNLPEvaluator;
        backend::AbstractNonlinearOracleBackend,
    )

Convert an AbstractNLPEvaluator into a SymbolicAD instance of MOI.NLPBlockData.
"""
function _nlp_block_data(
    d::MOI.AbstractNLPEvaluator;
    backend::AbstractNonlinearOracleBackend,
)
    MOI.initialize(d, [:ExprGraph])
    objective = d.has_nlobj ? _Function(MOI.objective_expr(d)) : nothing
    n = length(d.constraints)
    functions = Vector{_Function}(undef, n)
    bounds = Vector{MOI.NLPBoundsPair}(undef, n)
    for i in 1:n
        f, bound = _nonlinear_constraint(MOI.constraint_expr(d, i))
        functions[i], bounds[i] = f, bound
    end
    oracle = _NonlinearOracle(backend, objective, functions)
    return MOI.NLPBlockData(bounds, oracle, d.has_nlobj)
end

###
### ThreadedBackend
###

struct _ConstraintOffset
    index::Int
    ∇f_offset::Int
    ∇²f_offset::Int
end

struct ThreadedBackend <: AbstractNonlinearOracleBackend
    offsets::Vector{Vector{_ConstraintOffset}}
    ThreadedBackend() = new(Vector{_ConstraintOffset}[])
end

function MOI.initialize(
    oracle::_NonlinearOracle,
    backend::ThreadedBackend,
    ::Vector{Symbol},
)
    obj_hess_offset = 0
    if oracle.objective !== nothing
        h = oracle.objective.expr.hash
        obj_hess_offset = length(oracle.symbolic_functions[h].H)
    end
    offsets = Dict{UInt,Vector{_ConstraintOffset}}(
        h => _ConstraintOffset[] for h in keys(oracle.symbolic_functions)
    )
    index, ∇f_offset, ∇²f_offset = 1, 1, obj_hess_offset + 1
    for f in oracle.constraints
        push!(
            offsets[f.expr.hash],
            _ConstraintOffset(index, ∇f_offset, ∇²f_offset),
        )
        index += 1
        symbolic_function = oracle.symbolic_functions[f.expr.hash]
        ∇f_offset += length(symbolic_function.g)
        ∇²f_offset += length(symbolic_function.H)
    end
    for v in values(offsets)
        push!(backend.offsets, v)
    end
    return
end

function MOI.eval_constraint(
    oracle::_NonlinearOracle{ThreadedBackend},
    g::AbstractVector{Float64},
    x::AbstractVector{Float64},
)
    start = time()
    Threads.@threads for offsets in oracle.backend.offsets
        for c in offsets
            func = oracle.constraints[c.index]
            @inbounds g[c.index] = _eval_function(oracle, func, x)
        end
    end
    oracle.eval_constraint_timer += time() - start
    return
end

function MOI.eval_constraint_jacobian(
    oracle::_NonlinearOracle{ThreadedBackend},
    J::AbstractVector{Float64},
    x::AbstractVector{Float64},
)
    start = time()
    Threads.@threads for offset in oracle.backend.offsets
        for c in offset
            func = oracle.constraints[c.index]
            g = _eval_gradient(oracle, func, x)
            for i in 1:length(g)
                @inbounds J[c.∇f_offset+i-1] = g[i]
            end
        end
    end
    oracle.eval_constraint_jacobian_timer += time() - start
    return
end

function _hessian_lagrangian_structure(
    structure::Vector{Tuple{Int,Int}},
    oracle::_NonlinearOracle{ThreadedBackend},
    ::Dict{Tuple{Int,Int},Int},
    c::_Function,
)
    f = oracle.symbolic_functions[c.expr.hash]
    for (i, j) in f.H_structure
        row, col = c.ordered_variables[i], c.ordered_variables[j]
        push!(structure, row >= col ? (row, col) : (col, row))
    end
    return
end

function MOI.eval_hessian_lagrangian(
    oracle::_NonlinearOracle{ThreadedBackend},
    H::AbstractVector{Float64},
    x::AbstractVector{Float64},
    σ::Float64,
    μ::AbstractVector{Float64},
)
    fill!(H, 0.0)
    start = time()
    hessian_offset = 0
    if oracle.objective !== nothing && !iszero(σ)
        h = _eval_hessian(oracle, oracle.objective, x)
        for i in 1:length(h)
            H[i] = σ * h[i]
        end
        hessian_offset += length(h)
    end
    Threads.@threads for offset in oracle.backend.offsets
        for c in offset
            if iszero(μ[c.index])
                continue
            end
            func = oracle.constraints[c.index]
            h = _eval_hessian(oracle, func, x)
            for i in 1:length(h)
                @inbounds H[c.∇²f_offset+i-1] = μ[c.index] * h[i]
            end
        end
    end
    oracle.eval_hessian_lagrangian_timer += time() - start
    return
end

###
### JuMP integration
###

function _to_expr(
    nlp_data::JuMP._NLPData,
    data::JuMP._NonlinearExprData,
    subexpressions::Vector{Expr},
)
    tree = Any[]
    for node in data.nd
        expr = if node.nodetype == JuMP._Derivatives.CALL
            Expr(:call, JuMP._Derivatives.operators[node.index])
        elseif node.nodetype == JuMP._Derivatives.CALLUNIVAR
            Expr(:call, JuMP._Derivatives.univariate_operators[node.index])
        elseif node.nodetype == JuMP._Derivatives.SUBEXPRESSION
            subexpressions[node.index]
        elseif node.nodetype == JuMP._Derivatives.MOIVARIABLE
            MOI.VariableIndex(node.index)
        elseif node.nodetype == JuMP._Derivatives.PARAMETER
            nlp_data.nlparamvalues[node.index]
        else
            @assert node.nodetype == JuMP._Derivatives.VALUE
            data.const_values[node.index]
        end
        if 1 <= node.parent <= length(tree)
            push!(tree[node.parent].args, expr)
        end
        push!(tree, expr)
    end
    return tree[1]
end

_to_expr(::JuMP._NLPData, ::Nothing, ::Vector{Expr}) = nothing

function _nlp_block_data(
    model::JuMP.Model;
    backend::AbstractNonlinearOracleBackend,
)
    return _nlp_block_data(model.nlp_data; backend = backend)
end

function _nlp_block_data(
    nlp_data::JuMP._NLPData;
    backend::AbstractNonlinearOracleBackend,
)
    subexpressions = map(nlp_data.nlexpr) do expr
        return _to_expr(nlp_data, expr, Expr[])::Expr
    end
    nlobj = _to_expr(nlp_data, nlp_data.nlobj, subexpressions)
    objective = nlobj === nothing ? nothing : _Function(nlobj)
    functions = map(nlp_data.nlconstr) do c
        return _Function(_to_expr(nlp_data, c.terms, subexpressions))
    end
    return MOI.NLPBlockData(
        [MOI.NLPBoundsPair(c.lb, c.ub) for c in nlp_data.nlconstr],
        _NonlinearOracle(backend, objective, functions),
        objective !== nothing,
    )
end

function optimize_hook(
    model::JuMP.Model;
    backend::AbstractNonlinearOracleBackend = DefaultBackend(),
)
    nlp_block = _nlp_block_data(model; backend = backend)
    nlp_data = model.nlp_data
    model.nlp_data = nothing
    MOI.set(model, MOI.NLPBlock(), nlp_block)
    JuMP.optimize!(model; ignore_optimize_hook = true)
    model.nlp_data = nlp_data
    return
end


text:
export distance², distance, lerp, permute

# Point functions.

@inline distance²(p::Point, q::Point) = magnitude²(p - q)

@inline distance(p::Point, q::Point) = magnitude(p - q)

@inline Base.floor(p::Point) = typeof(p)(floor.(p.coordinates)...)

@inline Base.ceil(p::Point) = typeof(p)(ceil.(p.coordinates)...)

@inline Base.min(a::Point, b::Point) = promote_type(typeof(a), typeof(b))(min.(a.coordinates, b.coordinates)...)

@inline Base.max(a::Point, b::Point) = promote_type(typeof(a), typeof(b))(max.(a.coordinates, b.coordinates)...)

@inline lerp(t::T, p₀::Point, p₁::Point) where T<:Number = (one(T) - t) * p₀ + t * p₁

@inline permute(a::Point3{T}, (x, y, z)) where T = typeof(a)(a[x], a[y], a[z])

@inline Base.:(-)(a::Point{N,T}, b::Point{N,U}) where {N,T,U} = Vect{N, promote_type(T,U)}(a.coordinates - b.coordinates)


text:
using TreeTools

include("$(dirname(pathof(TreeTools)))/../test/methods/test.jl")
include("$(dirname(pathof(TreeTools)))/../test/prunegraft/test.jl")
include("$(dirname(pathof(TreeTools)))/../test/iterators/test.jl")
include("$(dirname(pathof(TreeTools)))/../test/splits/test.jl")

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools

max_forks_repo_name: PierreBarrat/TreeTools


text:
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
@parameters t
D = Differential(t)

# Modelica example
begin
    @named mass1 = HeatCapacitor(C=15, T_start=373.15)
    @named mass2 = HeatCapacitor(C=15, T_start=273.15)
    @named conduction = ThermalConductor(G=10)
    @named Tsensor1 = TemperatureSensor()
    @named Tsensor2 = TemperatureSensor()

    connections = [
        connect(mass1.port, conduction.port_a),
        connect(conduction.port_b, mass2.port),
        connect(mass1.port, Tsensor1.port),
        connect(mass2.port, Tsensor2.port),
    ]

    @named model = ODESystem(connections, t, systems=[mass1, mass2, conduction, Tsensor1, Tsensor2])
    sys = structural_simplify(model)
    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, [mass1.der_T => 1.0, mass2.der_T => 1.0], (0, 3.0))
    sol = solve(prob, DFBDF())
    @test sol.retcode == :Success
end


text:
using Gtk
using Test
using DrugInteractions

DrugInteractions._apps_should_persist[1] = false

drug_interactions_app()

@test DrugInteractions._apps[end] isa GtkWindow


text:
""" Write model in DynamicsTrainer to verilog-A netlist
    input/output config format: list of tuples, (header, type, port number)
    * header: either "" or "ddt".
    * type: either I or V
    * port number: must be between 1 and port number
    One tuple per model input/output. For example, ("ddt","V",2) interprets the input as the 2nd port's voltage derivative
    Similarly, ("", "I", 4) interprets the input as the 4th port's input current """
function generate_netlist_task(;
    model_name = "model",
    num_ports = 0,
    input_config = [],
    output_config = [],
    time_scale = 1.0,
    input_scale,
    input_bias,
    output_scale,
    output_bias,
)
    # create a function that can be directly run by a DynamicsTrainer as a callback task
    function gen_netlist(f)
        # extract parameters
        _, _, _, p_feed, p_ode, p_link = extract_params(f)

        write_netlist(
            f.ode,
            f.link,
            f.feedthrough,
            model_name,
            p_ode,
            p_link,
            p_feed,
            num_ports,
            input_config,
            output_config;
            time_scale = time_scale,
            input_scale = input_scale,
            input_bias = input_bias,
            output_scale = output_scale,
            output_bias = output_bias,
            dir = f.savepath,
        )
    end
    return gen_netlist
end

""" Produces model called <model_name>_model """
function write_netlist(
    dyn::O,
    link::L,
    feedthrough::F,
    model_name::String,
    p_ode,
    p_link,
    p_feed,
    num_ports,
    input_config,
    output_config;
    time_scale = 1.0,
    input_scale = ones(dyn.input_dim),
    input_bias = zeros(dyn.input_dim),
    output_scale = ones(link.output_dim),
    output_bias = zeros(link.output_dim),
    dir = ".",
) where {O<:Ode} where {L<:Link} where {F<:Feedthrough}
    # dimension data
    input_dim = dyn.input_dim
    state_dim = dyn.state_dim

    # TODO: currently assumes σ_u is either identity or tanh
    if dyn.σ_u == identity
        gate_u = ""
    elseif dyn.σ_u == tanh
        gate_u = "tanh"
    else
        gate_u = ""
    end

    # odelist: list of strings, one for each dxdt
    # hidden_layers: dict of strings, one for each hidden layer unit
    # param_dict: dict of param name, param value
    odelist, hidden_layers, param_dict_ode = convert_verilogA(dyn, p_ode)

    # param_dict_output: dict of param name, param value
    # output_list: list of strings, one for each output dim
    param_dict_output, output_list = convert_verilogA(
        link,
        feedthrough,
        p_link,
        p_feed;
        i_scale = output_scale,
        i_bias = -output_scale .* output_bias,
    )

    # write verilog A file
    # implements architecture-specific connections
    # instance specific parameter values written to the .inc netlist
    open(dir * "/" * model_name * ".va", "w") do io
        # write includes
        write(io, "`include \"constants.vams\"\n")
        write(io, "`include \"disciplines.vams\"\n")
        write(io, "(* compact_module *)\n")
        # write(io, "(* ignore_hidden_state=\"all\" *)\n")

        # write inputs & outputs: ports
        iolist = ""
        for i = 1:num_ports
            iolist *= "a$(i),"
        end

        # write states as inouts as well in order to set initial condition
        for i = 1:state_dim
            iolist *= "x$(i),"
        end

        write(io, "module " * model_name * "(" * iolist * "c);\n")

        # define inouts, with branches
        write(io, "inout " * iolist * "c;\n")
        inlist = ""
        for i = 1:input_dim
            inlist *= "_u$(i),"
            inlist *= "u$(i),"
        end
        write(io, "electrical " * iolist * inlist * "c;\n")

        for i = 1:num_ports
            write(io, "branch (a$(i),c) Ia$(i);\n")
        end

        # define state variables: electrical & branch
        write(io, "electrical ")
        for i = 1:state_dim-1
            write(io, "dx$(i),")
            # write(io, "x$(i),dx$(i),")
        end
        write(io, "dx$(state_dim);\n")
        # write(io, "x$(i),dx$(i);\n")

        for i = 1:state_dim
            write(io, "branch (x$(i),c) Ix$(i);\n")
        end

        # write parameter definitions
        for param in keys(param_dict_ode)
            if param != "tauinv"
                write(io, "parameter real " * param * " = 0 from (-inf:inf);\n")
            else
                write(io, "parameter real " * param * " = 0 from(0:inf);\n")
            end
        end

        for param in keys(param_dict_output)
            write(io, "parameter real " * param * " = 0 from (-inf:inf);\n")
        end

        # initialize hidden layer units
        if length(hidden_layers) > 0
            write(io, "real")
            key_list = [key for key in keys(hidden_layers)]
            for i = 1:length(key_list)
                write(io, " " * key_list[i])
                if i == length(hidden_layers)
                    write(io, ";\n")
                else
                    write(io, ",")
                end
            end
        end

        # begin writing connections
        write(io, "analog begin\n")

        # connect inputs
        for i = 1:length(input_config)
            header = input_config[i][1]
            type = input_config[i][2]
            port = input_config[i][3]

            # write(io, "V(u$(i),c)<+V(a$(i),c);\n")
            write(
                io,
                "V(_u$(i),c) <+ $(1.0/input_scale[i])*$(header)($(type)(a$(port),c)) + $(input_bias[i]);\n",
            )
            write(io, "V(u$(i),c) <+ $(gate_u)(V(_u$(i),c));\n")
        end

        # connect derivatives
        for i = 1:state_dim
            write(io, "V(dx$(i),c)<+(" * string(time_scale) * ")*ddt(V(x$(i),c));\n")
        end

        # define hidden layers
        for key in keys(hidden_layers)
            write(io, key * " = " * hidden_layers[key] * "\n")
        end

        # write ode
        for ode in odelist
            write(io, ode * "\n")
        end

        # write output layer
        for i = 1:length(output_config)
            type = output_config[i][2]
            port = output_config[i][3]

            if type == "V"
                header = "V(a$(port),c) <+ "
            elseif type == "I"
                header = "I(Ia$(port)) <+ "
            elseif type == "O"
                write(io, "I(Ia$(port)) <+ 0.0; \n")
                continue
            else
                header = "$(type)(Ia$(port)) <+ "
            end

            write(io, header * output_list[i] * "\n")
        end

        # end
        write(io, "end // end analog\n")
        write(io, "endmodule\n")

    end # io for .va file

    #
    # write netlist file
    # instance specific parameter values written to the .inc netlist
    open(dir * "/" * model_name * ".inc", "w") do io
        write(io, "// subcircuit for learned Verilog-A model\n")
        write(io, "simulator lang=spectre\n")
        write(io, "ahdl_include \"" * model_name * ".va\"\n")

        # terminal definition
        iolist = ""
        for i = 1:num_ports
            iolist *= "a$(i) "
        end
        write(io, "subckt " * model_name * "_model " * iolist * "c\n")

        # # set initial condition guess to zero
        # write(io, "nodeset ")
        # for i=1:state_dim
        #     write(io, "x$(i)=0 ")
        # end
        # write(io, "\n")

        # include states in inner verilog model
        for i = 1:state_dim
            iolist *= "x$(i) "
        end
        write(io, model_name * "0 (" * iolist * "c) " * model_name)

        for k in sort([string(val) for val in keys(param_dict_ode)])
            write(io, " " * k * "=" * string(Float64(param_dict_ode[k])))
        end
        for k in sort([string(val) for val in keys(param_dict_output)])
            write(io, " " * k * "=" * string(Float64(param_dict_output[k])))
        end
        write(io, "\n")
        write(io, "ends\n")
    end # io for .inc file
end

""" Convert CTRNN0 to VerilogA format """
function convert_verilogA(f::CTRNN0, p)
    # unpack params
    logτ, A, B, μ, ν = unpack_params(f, p)

    # write dictionary of parameters with values
    param_dict = Dict()
    for i = 1:f.state_dim
        for j = 1:f.state_dim
            param_dict["A_$(i)_$(j)"] = A[i, j]
        end
        for j = 1:f.input_dim
            param_dict["B_$(i)_$(j)"] = B[i, j]
        end
        param_dict["mu_$(i)"] = μ[i]
    end
    for i = 1:f.state_dim
        param_dict["nu_$(i)"] = ν[i]
    end
    param_dict["tauinv"] = exp.(-logτ)[1]

    # hidden layer
    hidden_layers = Dict()
    # hidden_layer = []
    for i = 1:f.state_dim
        layer_name = "h_$(i)"
        layer_string = "(mu_$(i))"
        for j = 1:f.state_dim
            layer_string *= "+(A_$(i)_$(j))*(V(x$(j),c))"
        end
        for j = 1:f.input_dim
            layer_string *= "+(B_$(i)_$(j))*(V(u$(j),c))"
        end

        # TODO: currently assumes σ is tanh or identity
        if f.σ == tanh
            hidden_layers[layer_name] = "tanh($(layer_string));"
        elseif f.σ == relu
            hidden_layers[layer_name] = "max(0,$(layer_string));"
        else
            hidden_layers[layer_name] = "$(layer_string);"
        end
        # push!(hidden_layer, "$(layer_name) = $(activation)($(layer_string));")
    end

    # write ode derivative assignments
    odelist = []
    for i = 1:f.state_dim
        string_i = "I(Ix$(i)) <+ (V(dx$(i),c)) + (tauinv)*(V(x$(i),c)) - (nu_$(i)) - h_$(i)"

        push!(odelist, string_i * ";")
    end

    return odelist, hidden_layers, param_dict
end

""" Convert CTRNN1 to VerilogA format """
function convert_verilogA(f::CTRNN1, p)
    # unpack parameters
    logτ, W, A, B, μ, ν = unpack_params(f, p)

    # write dictionary of parameters with values
    param_dict = Dict()
    for i = 1:f.hidden_dim
        for j = 1:f.state_dim
            param_dict["A_$(i)_$(j)"] = A[i, j]
            param_dict["W_$(j)_$(i)"] = W[j, i]
        end
        for j = 1:f.input_dim
            param_dict["B_$(i)_$(j)"] = B[i, j]#/v_scale[j]
        end
        param_dict["mu_$(i)"] = μ[i]
    end
    for i = 1:f.state_dim
        param_dict["nu_$(i)"] = ν[i]
    end
    param_dict["tauinv"] = exp.(-logτ)[1]

    # hidden layer
    hidden_layers = Dict()
    # hidden_layer = []
    for i = 1:f.hidden_dim
        layer_name = "h_$(i)"
        layer_string = "(mu_$(i))"
        for j = 1:f.state_dim
            layer_string *= "+(A_$(i)_$(j))*(V(x$(j),c))"
        end
        for j = 1:f.input_dim
            layer_string *= "+(B_$(i)_$(j))*(V(u$(j),c))"
        end

        # TODO: currently assumes σ is tanh or identity
        if f.σ == tanh
            hidden_layers[layer_name] = "tanh($(layer_string));"
        elseif f.σ == relu
            hidden_layers[layer_name] = "max(0,$(layer_string));"
        else
            hidden_layers[layer_name] = "$(layer_string);"
        end
        # push!(hidden_layer, "$(layer_name) = $(activation)($(layer_string));")
    end

    # write ode derivative assignments
    odelist = []
    for i = 1:f.state_dim
        string_i = "I(Ix$(i)) <+ (V(dx$(i),c)) + (tauinv)*(V(x$(i),c)) - (nu_$(i))"
        for j = 1:f.hidden_dim
            string_i *= " - (W_$(i)_$(j))*h_$(j)"
        end

        push!(odelist, string_i * ";")
    end

    return odelist, hidden_layers, param_dict
end

""" Convert AffineLayer to VerilogA format """
function convert_verilogA(
    f::AffineLayer,
    g::F,
    p_link,
    p_feed;
    v_scale = ones(f.output_dim),
    i_scale = ones(f.output_dim),
    i_bias = zeros(f.output_dim),
) where {F<:Feedthrough}
    # TODO: assumes either tanh or identity
    if f.σ == identity
        gate_link = ""
    elseif f.σ == tanh
        gate_link = "tanh"
    else
        gate_link = ""
    end

    # unpack parameters
    C, D, b = unpack_params(f, p_link)
    D = unpack_params(g, p_feed)

    # write dictionary of parameters with values for output layer
    param_dict_output = Dict()
    for i = 1:f.output_dim
        for j = 1:f.state_dim
            param_dict_output["C_$(i)_$(j)"] = C[i, j] * i_scale[i]
        end
        if D != false
            for j = 1:g.input_dim
                param_dict_output["D_$(i)_$(j)"] = D[i, j] * i_scale[i]
            end
        end
        param_dict_output["b_$i"] = b[i] * i_scale[i] + i_bias[i]
    end

    # write state,input -> output assignments
    output_list = []
    for i = 1:f.output_dim
        string_i = "(b_$(i))"
        for j = 1:f.state_dim
            string_i *= "+(C_$(i)_$(j))*$(gate_link)(V(x$(j),c))"
        end
        if D != false
            for j = 1:g.input_dim
                if D[i, j] != 0.0
                    string_i *= "+(D_$(i)_$(j))*(V(u$(j),c))"
                end
            end
        end
        push!(output_list, string_i * ";")
    end

    return param_dict_output, output_list
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
function MOI.set(
    model::MOI.ModelLike,
    attr::BackwardInVariablePrimal,
    bridge::MOI.Bridges.Variable.VectorizeBridge,
    value,
)
    MOI.set(model, attr, bridge.variable, value)
end


text:
using ValidatedNumerics
using FactCheck

setprecision(Interval, Float64)

facts("displaymode tests") do

    context("Interval") do

        a = 1..2
        b = -1.1..1.3
        c = Interval(pi)
        d = @interval(π)

        context("6 sig figs") do
            displaymode(format=:standard, sigfigs=6)

            @fact string(a) --> "[1, 2]"
            @fact string(b) --> "[-1.10001, 1.30001]"
            @fact string(c) --> "[3.14159, 3.1416]"
            @fact string(d) --> "[3.14159, 3.1416]"
        end

        context("10 sig figs") do
            displaymode(sigfigs=10)

            @fact string(a) --> "[1, 2]"
            @fact string(b) --> "[-1.100000001, 1.300000001]"
            @fact string(c) --> "[3.141592653, 3.141592654]"
            @fact string(d) --> "[3.141592653, 3.141592654]"
        end

        context("20 sig figs") do
            displaymode(sigfigs=20)

            @fact string(a) --> "[1, 2]"
            @fact string(b) --> "[-1.1000000000000000889, 1.3000000000000000445]"
            @fact string(c) --> "[3.1415926535897931159, 3.141592653589793116]"
            @fact string(d) --> "[3.1415926535897931159, 3.1415926535897935601]"
        end

        context("Full") do
            displaymode(format=:full)

            @fact string(a) --> "Interval(1.0, 2.0)"
            @fact string(b) --> "Interval(-1.1, 1.3)"
            @fact string(c) --> "Interval(3.141592653589793, 3.141592653589793)"
            @fact string(d) --> "Interval(3.141592653589793, 3.1415926535897936)"
        end

        context("Midpoint") do
            displaymode(format=:midpoint, sigfigs=6)

            @fact string(a) --> "1.5 ± 0.5"
            @fact string(b) --> "0.1 ± 1.20001"
            @fact string(c) --> "3.14159 ± 0"
            @fact string(d) --> "3.14159 ± 4.4409e-16"
        end

    end

    context("DecoratedInterval") do
        a = @decorated(1, 2)

        displaymode(format=:standard, decorations=false)

        @fact string(a) --> "[1, 2]"

        displaymode(format=:standard, decorations=true)

        @fact string(a) --> "[1, 2]_com"

    end

    setprecision(Interval, 128)

    context("BigFloat intervals") do
        displaymode(format=:standard, decorations=false)

        a = @interval big(1)
        @fact string(a) --> "[1, 1]₁₂₈"

        displaymode(format=:full)
        @fact string(a) --> "Interval(1.000000000000000000000000000000000000000, 1.000000000000000000000000000000000000000)"

        a = DecoratedInterval(big(2), big(3), com)

        displaymode(format=:standard, decorations=false)
        @fact string(a) --> "[2, 3]₁₂₈"

        displaymode(format=:standard, decorations=true)
        @fact string(a) --> "[2, 3]₁₂₈_com"

        displaymode(format=:full)
        @fact string(a) --> "DecoratedInterval(Interval(2.000000000000000000000000000000000000000, 3.000000000000000000000000000000000000000), com)"

    end

    setprecision(Interval, Float64)

    context("IntervalBox") do

        displaymode(format=:standard, sigfigs=6)

        X = IntervalBox(1..2, 3..4)
        @fact string(X) --> "[1, 2] × [3, 4]"

        X = IntervalBox(1.1..1.2, 2.1..2.2)
        @fact string(X) --> "[1.09999, 1.20001] × [2.09999, 2.20001]"

        X = IntervalBox(-Inf..Inf, -Inf..Inf)
        @fact string(X) --> "[-∞, ∞] × [-∞, ∞]"

        displaymode(format=:full)
        @fact string(X) --> "IntervalBox(Interval(-Inf, Inf), Interval(-Inf, Inf))"

    end
end