text:
module TestVariable

using Test

import MutableArithmetics
const MA = MutableArithmetics

using MultivariatePolynomials
const MP = MultivariatePolynomials
import MultivariateBases
const MB = MultivariateBases
using DynamicPolynomials
#using TypedPolynomials

using JuMP
using PolyJuMP

function test_variable_macro_with_Poly(var)
    x = MP.similarvariable(var, Val{:x})
    y = MP.similarvariable(var, Val{:y})
    X = [1, y, x^2, x, y^2]

    m = Model()
    @variable(m, a)
    var_poly   = polynomial(a * x * y)
    var_poly_x = polynomial(a * x)
    aff_poly   = polynomial((a + 1) * x * y)
    aff_poly_x = polynomial((a + 1) * x)
    @test_throws ErrorException @variable m p Poly(X) unknown_kw=1
    @test_throws ErrorException @variable m p == 1 Poly(X)
    @test_throws ErrorException @variable m p Poly(X) start=1
    @test_throws ErrorException @variable m 0 <= p <= 1 Poly(X)
    @test_throws ErrorException @variable m 0 >= p >= 1 Poly(X)
    @test_throws ErrorException @variable m p <= 0 Poly(X)
    @test_throws ErrorException @variable m p >= 1 Poly(X)
    @test_throws ErrorException @variable m p == 1 Poly(X)
    @test_throws ErrorException @variable m p >= 0 Poly(X)
    @test_throws ErrorException @variable(m, p3[2:3] >= 0, Poly(X))
end

function _test_variable(m, p, monos, binary = false, integer = false,
                 vars = true, use_y = true)
    PT = polynomialtype(use_y ? first(monos) : first(variables(first(monos))), vars ? JuMP.VariableRef : JuMP.AffExpr)
    @test isa(p, PT)
    @test monomials(p) == monovec(monos)
    if vars
        @test all(α -> JuMP.is_binary(α) == binary, coefficients(p))
        @test all(α -> JuMP.is_integer(α) == integer, coefficients(p))
    else
        @test all(α -> JuMP.is_binary(first(keys(α.terms))) == binary, coefficients(p))
        @test all(α -> JuMP.is_integer(first(keys(α.terms))) == integer, coefficients(p))
    end
end

function test_MonomialBasis(var)
    x = MP.similarvariable(var, Val{:x})
    y = MP.similarvariable(var, Val{:y})
    X = [1, y, x^2, x, y^2]
    m = Model()
    @variable m p1[1:3] Poly(X)
    PT = polynomialtype(x * y, JuMP.VariableRef)
    @test isa(p1, Vector{PT})
    _test_variable(m, p1[1], X)
    @variable(m, p2, Poly(X), integer=true)
    _test_variable(m, p2, X, false, true)
    @variable(m, p3[2:3], Poly(X))
    @test isa(p3, JuMP.Containers.DenseAxisArray{PT,1,Tuple{UnitRange{Int}}})
    _test_variable(m, p3[2], X)
    @variable(m, p4[i=2:3,j=i:4], Poly(X), binary=true)
    _test_variable(m, p4[2,3], X, true)

    X = [x^2, y^2]
    @variable m p5[1:3] Poly(X)
    @test isa(p5, Vector{PT})
    _test_variable(m, p5[1], X)
    @variable(m, p6, Poly(MB.MonomialBasis(X)), integer=true)
    _test_variable(m, p6, X, false, true)
end

function test_ScaledMonomialBasis(var)
    x = MP.similarvariable(var, Val{:x})
    m = Model()
    @variable(m, p1, Poly(MB.ScaledMonomialBasis([1, x, x^2])), Int)
    _test_variable(m, p1, monovec([1, x, x^2]), false, true, false, false)
end

function test_FixedPolynomialBasis(var)
    x = MP.similarvariable(var, Val{:x})
    y = MP.similarvariable(var, Val{:y})
    m = Model()
    @variable(m, p1, Poly(MB.FixedPolynomialBasis([1 - x^2, x^2 + 2])), Bin)
    _test_variable(m, p1, monovec([x^2, 1]), true, false, false, false)
    @variable(m, p2[1:2], Poly(MB.FixedPolynomialBasis([1 - x^2, x^2 + 2])))
    _test_variable(m, p2[1], monovec([x^2, 1]), false, false, false, false)
    # Elements of the basis have type monomial
    @variable(m, p3[2:3], Poly(MB.FixedPolynomialBasis([x, x^2])))
    _test_variable(m, p3[2], monovec([x^2, x]), false, false, false, false)
    # Elements of the basis have type term
    @variable(m, p4[1:2], Poly(MB.FixedPolynomialBasis([1, x, x^2])), Int)
    _test_variable(m, p4[1], monovec([x^2, x, 1]), false, true, false, false)
    # Elements of the basis have type variable
    @variable(m, p5[-1:1], Poly(MB.FixedPolynomialBasis([x, y])), integer=true)
    _test_variable(m, p5[0], monovec([x, y]), false, true, false)
    @variable(m, p6[-1:1], Poly(MB.FixedPolynomialBasis([x])), integer=true)
    _test_variable(m, p6[0], monovec([x]), false, true, true, false)
end

function test_value_function(var)
    x = MP.similarvariable(var, Val{:x})
    y = MP.similarvariable(var, Val{:y})
    m = Model()
    @variable m α
    @variable m β
    p = α * x*y + β * x^2
    JuMP.fix(α, 2)
    JuMP.fix(β, 3)
    @test_broken JuMP.value(p) == 2x*y + 3x^2
    @test JuMP.value(fix_value, p) == 2x*y + 3x^2
end

function runtests(var)
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name) $(typeof(var))" begin
                getfield(@__MODULE__, name)(var)
            end
        end
    end
end

end

import DynamicPolynomials
DynamicPolynomials.@polyvar(x)
TestVariable.runtests(x)
import TypedPolynomials
TypedPolynomials.@polyvar(y)
TestVariable.runtests(y)


text:
# FASTA Reader
# ============

type Reader <: Bio.IO.AbstractReader
    state::Bio.Ragel.State
    index::Nullable{Index}

    function Reader(input::BufferedStreams.BufferedInputStream, index)
        return new(Bio.Ragel.State(file_machine.start_state, input), index)
    end
end

"""
    FASTA.Reader(input::IO; index=nothing)

Create a data reader of the FASTA file format.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *fai* is supported)
"""
function Reader(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = Index(index)
    else
        if index != nothing
            throw(ArgumentError("index must be a filepath or nothing"))
        end
    end
    return Reader(BufferedStreams.BufferedInputStream(input), index)
end

function Base.eltype(::Type{Reader})
    return Record
end

function Bio.IO.stream(reader::Reader)
    return reader.state.stream
end

function Base.getindex(reader::Reader, name::AbstractString)
    if isnull(reader.index)
        throw(ArgumentError("no index attached"))
    end
    seekrecord(reader.state.stream, get(reader.index), name)
    reader.state.cs = file_machine.start_state
    reader.state.finished = false
    return read(reader)
end

info("compiling FASTA")
const record_machine, file_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    rep1 = Automa.RegExp.rep1
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt
    any = Automa.RegExp.any
    space = Automa.RegExp.space

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(opt('\r'), lf)
    end

    hspace = re"[ \t\v]"

    identifier = rep(any() \ space())
    identifier.actions[:enter] = [:mark]
    identifier.actions[:exit]  = [:identifier]

    description = cat(any() \ hspace, re"[^\r\n]*")
    description.actions[:enter] = [:mark]
    description.actions[:exit]  = [:description]

    header = cat('>', identifier, opt(cat(rep1(hspace), description)))
    header.actions[:enter] = [:anchor]
    header.actions[:exit]  = [:header]

    # '*': terminal, `-': gap
    letters = re"[A-Za-z*\-]+"
    letters.actions[:enter] = [:movable_anchor]
    letters.actions[:exit]  = [:letters]

    sequence = opt(cat(letters, rep(cat(rep1(alt(hspace, newline)), letters))))
    sequence.actions[:enter] = [:sequence_start]

    record = cat(header, rep1(newline), sequence, rep1(newline))
    record.actions[:exit] = [:record]

    record_trailing = cat(header, rep1(newline), sequence)
    record_trailing.actions[:exit] = [:record]

    file = cat(rep(newline), rep(record), opt(record_trailing))

    return map(Automa.compile, (record, file))
end)()

#=
write("fasta.dot", Automa.machine2dot(file_machine))
run(`dot -Tsvg -o fasta.svg fasta.dot`)
=#

const record_actions = Dict(
    :identifier  => :(record.identifier  = (mark:p-1)),
    :description => :(record.description = (mark:p-1)),
    :header => quote
        @assert record.data === data
        copy!(record.data, 1, record.data, 1, p - 1)
        filled += p - 1
        if p ≤ endof(data)
            data[p] = UInt8('\n')
            filled += 1
        end
    end,
    :letters => quote
        let len = p - mark
            copy!(record.data, filled + 1, record.data, mark, p - mark)
            filled += len
            record.sequence = first(record.sequence):last(record.sequence)+len
        end
    end,
    :sequence_start => :(record.sequence = filled+1:filled),
    :record => :(record.filled = 1:filled),
    :anchor => :(),
    :movable_anchor => :(mark = p),
    :mark => :(mark = p),
    :countline => :(#= linenum += 1 =#))
eval(
    Bio.ReaderHelper.generate_index_function(
        Record,
        record_machine,
        :(filled = mark = 0),
        record_actions))
eval(
    Bio.ReaderHelper.generate_read_function(
        Reader,
        file_machine,
        :(filled = mark = 0),
        merge(record_actions, Dict(
            :identifier  => :(record.identifier  = (mark:p-1) - stream.anchor + 1),
            :description => :(record.description = (mark:p-1) - stream.anchor + 1),
            :header => quote
                range = Bio.ReaderHelper.upanchor!(stream):p-1
                Bio.ReaderHelper.resize_and_copy!(record.data, filled + 1, data, range)
                filled += length(range)
                if filled + 1 ≤ endof(record.data)
                    record.data[filled+1] = UInt8('\n')
                else
                    push!(record.data, UInt8('\n'))
                end
                filled += 1
            end,
            :sequence_start => :(record.sequence = filled+1:filled),
            :letters => quote
                let len = p - stream.anchor
                    Bio.ReaderHelper.append_from_anchor!(record.data, filled + 1, stream, p - 1)
                    filled += len
                    record.sequence = first(record.sequence):last(record.sequence)+len
                end
            end,
            :record => quote
                record.filled = 1:filled
                found_record = true
                @escape
            end,
            :countline => :(linenum += 1),
            :anchor => :(Bio.ReaderHelper.anchor!(stream, p)),
            :movable_anchor => :(Bio.ReaderHelper.anchor!(stream, p, false))))))

ext: jl


text:
using DelimitedFiles
using LOWDER

include("quad_generator.jl")

function runtest_quad(
                        list::Vector{Int64},
                        filename::String;
                        n::Int64=2
                        )

    nprob = length(list)

    file = open( filename, "w" )

    for i = 1:nprob

        try

            # Generates the problem
            ( fmin, l, u ) = genprob( list[i], 10.0 * ones( Float64, n ) )

            x = ones( Float64, n)
            p = length(fmin)
            n_points = n + 1

            # Solves the problem using 'lowder'
            sol = LOWDER.lowder( fmin, x, l, u; m = n_points )

            # Saves info
            text = "$( n ) $( p ) $( sol.status ) $( sol.true_val ) $( sol.iter ) $( sol.nf ) $( sol.nf / p ) $( sol.stationarity ) $( sol.sample_radius ) $( sol.tr_radius ) $( sol.index ) $( sol.f )"
            print( file, text )
            for j = 1 : ( n - 1 )

                print( file, " $(sol.solution[ j ])" )

            end
            println( file, " $( sol.solution[ n ] )" )

            println("succes!")

        catch

            println( file, "execution_fail" )

            println("fail!")

        end

    end

    close( file )

end

# -----------------------------------------------
# Path to file
# -----------------------------------------------

filename = "../data_files/quad_cons.dat"

# -----------------------------------------------
# Funtion call
# -----------------------------------------------

list = [2, 5, 10, 15, 20, 50, 100]

runtest_quad(list, filename)


text:
include("header.jl")

@testset "distributions" begin
    @test isa(gaussian(10),Array)
    @test isa(xavier(10,10),Array)
    @test isa(bilinear(Float32,2,2,128,128),Array)
end

nothing


text:
#=

=#

########################################
# Model Cost / Error
########################################

# With regularisation
function loss(net, X, y)

    # Data type
    F = net.dtype

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # Cost
    J = F(0)

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Cost without regularisation
        J += L.cost(L.A, L.y)
    end

    # Divide by number of observations
    J /= (net.batch_size * num_batches)

    # Add regularisation cost
    J += net.reg_cost(net)

    # Add sparsity cost
    J += sparse_cost(net)

    return J
end

# Without regularisation
function val_loss(net, X, y)

    # Data type
    F = net.dtype

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # Cost
    J = F(0)

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Cost without regularisation
        J += L.cost(L.A, L.y)
    end

    # Divide by number of observations
    J /= (net.batch_size * num_batches)

    return J
end

#loss(net) = loss(net, net.data.X_train, net.data.y_train)

########################################
# Model Cost / Error and Accuracy
########################################

# With regularisation
function loss_and_accuracy(net, X, y)

    # Data type
    F = net.dtype

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # Cost and Count of Correct Predictions
    J = F(0)
    num_correct = 0

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Cost without regularisation
        J += L.cost(L.A, L.y)

        # Count correct predictions
        num_correct += count_correct(L)
    end

    # Divide by number of observations
    num_used = net.batch_size * num_batches
    J /= num_used

    # Add regularisation cost
    J += net.reg_cost(net)

    # Add sparsity cost
    J += sparse_cost(net)

    # Accuracy (%)
    accuracy = num_correct / num_used * 100

    return J, accuracy
end

# Without regularisation
function val_loss_and_accuracy(net, X, y)

    # Data type
    F = net.dtype

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # Cost and Count of Correct Predictions
    J = F(0)
    num_correct = 0

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Cost without regularisation
        J += L.cost(L.A, L.y)

        # Count correct predictions
        num_correct += count_correct(L)
    end

    # Divide by number of observations
    num_used = net.batch_size * num_batches
    J /= num_used

    # Accuracy (%)
    accuracy = num_correct / num_used * 100

    return J, accuracy
end

########################################
# Model Accuracy
########################################

function accuracy(net, X, y)

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # Count of Correct Predictions
    num_correct = 0

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Count correct predictions
        num_correct += count_correct(L)
    end

    # Accuracy (%)
    return num_correct / (net.batch_size * num_batches) * 100
end

########################################
# Count Correct Predictions
########################################

# Counts the correct predictions in Softmax Output Layer
function count_correct(ℓ::SoftmaxOutputLayer)

    num_correct = 0
    for i in indices(ℓ.A, 2)
        indmax(ℓ.A[:, i]) == ℓ.y[i] && (num_correct += 1)
    end

    return num_correct
end

# Counts the correct predictions in Logistic Output Layer
function count_correct(ℓ::LogisticOutputLayer)

    num_correct = 0
    for i in indices(ℓ.A, 2)
        round(ℓ.A[i]) == ℓ.y[i] && (num_correct += 1)
    end

    return num_correct
end

########################################
# Model Predictions
########################################

# Outputs the predicted label from Softmax Output Layer
function predictions(ℓ::SoftmaxOutputLayer)

    # Number of observations / size of vector
    num_obs = size(ℓ.A, 2)

    # Prediction vector
    pred = Array(Int, num_obs)
    for i in 1:num_obs
        pred[i] = indmax(ℓ.A[:, i])
    end

    return pred
end

# Predictions for Logistic Output Layer - hypothesis closest to 0 or 1
predictions(ℓ::LogisticOutputLayer) = round(Int, ℓ.A)

# Predicted value for Linear Output Layer = Hypothesis
predictions(ℓ::LinearOutputLayer) = ℓ.A

# Predictions for data set
function predictions(net::NeuralNet, X, y)

    # Number of batches
    num_obs = size(X)[end]
    num_batches = div(num_obs, net.batch_size)

    # ℓ1 = first layer, L = final layer
    ℓ1 = net.layers[1]
    L = net.layers[end]

    # Create empty array for predictions
    pred = Array(eltype(L.y), 0)

    # Loop over mini batches
    for batch = 1:num_batches

        # Update network with next batch
        stop = batch * net.batch_size
        start = stop - net.batch_size + 1
        ℓ1.A[:] = viewbatch(X, start, stop)
        L.y[:] = viewbatch(y, start, stop)

        # Feed forward
        fwd_prop!(net)

        # Push predictions to array
        append!(pred, predictions(L))
    end

    # Accuracy (%)
    return pred
end


text:
#----------------------------------------------------------------------------------------------------------------------------------------

# Implements solution to the diffusion equation for the fluence in the steady-state and time-domain in a slab [1][2].
#
# [1] Contini et. al. "Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory."
#     Applied optics 36.19 (1997): 4587-4599.
# [2] Martelli et al. Light propagation through biological tissue and other diffusive media: theory, solutions and software.
#     Vol. 10. No. 3.824746. Bellingham: SPIE press, 2010.
#----------------------------------------------------------------------------------------------------------------------------------------
""" Structure containing inputs for simulating the fluence under the diffusion approximation in the slab space."""
struct DAslab{T <: AbstractFloat} <: DiffusionParameters
    ρ::T                                 # distance away from isotropic point source (cm)
    μa::T                                # absorption coefficient (cm⁻¹)
    μsp::T                               # reduced scattering coefficient (cm⁻¹)
    n_med::T                             # medium's index of refraction
    n_ext::T                             # external medium's index of refraction
    ω::T                                 # modulation frequency (1/ns)
    z::T                                 # z depth of detector within medium (cm)
    t::Union{T, AbstractVector{T}}       # time vector (ns)
    s::T                                 # slab thickness (cm) in z direction
    xs::Int                              # the number of sources to compute in the series

    # do not need to provide these arguments (calculated from previous inputs)
    D::T                                 # Diffusion coefficient
    ν::T                                 # Speed of light in medium (cm/ns)
    z0::T                                # isotroptic source depth
    zb::T                                # Extrapolation length
    function DAslab{T}(ρ::T, μa::T, μsp::T, n_med::T, n_ext::T, ω::T, z::T, t::Union{T, AbstractVector{T}}, s::T, xs::Int) where {T <: AbstractFloat}
        @assert ρ >= zero(T) "ρ must greater than or equal to 0"
        @assert μa >= zero(T) "μa must greater than or equal to 0"
        @assert μsp > zero(T) "μsp must greater than 0"
        @assert s > zero(T) "s must greater than 0"
        @assert xs > 0 "xs must greater than 0"
        @assert all(t .> zero(T)) "t must be positive"
        D = D_coeff(μsp)
        return new{T}(ρ, μa, μsp, n_med, n_ext, ω, z, t, s, xs, D, ν_coeff(n_med), z0_coeff(μsp), zb_coeff(A_coeff(n_med / n_ext), D))
    end
end

"""
    Generator function for DAslab structure.

Provides default parameters to use in the semi-infinite space fluence calculation in either the CW, FD, or TD.

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_med`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `ω`: modulation frequency (1/ns)
- `z`: the z-depth orthogonal from the boundary (cm)
- `t`: the time vector (ns).
- `s`: the thickness (z) of the slab (cm)
- `xs`: the number of sources to compute in the series

# Examples
```
julia> data = DAslab() # return default parameters
julia> data = DAslab(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults
julia> fluence_DA_slab_CW(data) # can then be used with corresponding functions
```
"""
function DAslab(;
    ρ::T = 1.0,
    μa::T = 0.1,
    μsp::T = 10.0,
    n_med::T = 1.0,
    n_ext::T = 1.0,
    ω::T = 0.0,
    z::T = 0.0,
    t::Union{T, AbstractVector{T}} = 1.0,
    s::T = 2.0,
    xs::Int = 10
) where {T<:AbstractFloat}
    return DAslab{T}(ρ, μa, μsp, n_med, n_ext, ω, z, t, s, xs)
end

#--------------------------------------
# Steady-State Fluence
#--------------------------------------
"""
    fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10)

Compute the steady-state fluence from a slab geometry (x, y -> inf, z -> finite).
The fluence is calculated from an infinite summation that converges rather rapidly.
Generally only 5 terms are needed (xs = 5), but more terms will be required if ρ >> s.

# Arguments
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `s`: the thickness (z-depth) of the slab (cm)
- `z`: the z-depth within slab (cm)
- `xs`: the number of sources to compute in the series

# Examples
```
julia> fluence_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0)
```
"""
function fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0, xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    μeff = sqrt(3 * μa * μsp)
    return _sum_fluence_DA_slab_CW(xs, s, params.zb, params.z0, ρ, z, μeff, params.D)
end

"""
    fluence_DA_slab_CW(data::DiffusionParameters)

Wrapper to `fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0, xs = 10)`

# Examples
```
julia> data = DAslab(ρ = 1.0, s = 3.0) # use structure to generate inputs
julia> fluence_DA_slab_CW(data) # then call the function
```
"""
function fluence_DA_slab_CW(data::DiffusionParameters)
    μeff = sqrt(3 * data.μa * data.μsp)
    return _sum_fluence_DA_slab_CW(data.xs, data.s, data.zb, data.z0, data.ρ, data.z, μeff, data.D)
end

function _sum_fluence_DA_slab_CW(xs, s, zb, z0, ρ, z, μeff, D)
    ϕ = zero(eltype(μeff))

    for m in xs:-1:1
        ϕ += _kernel_fluence_DA_slab_CW(-m, s, zb, z0, ρ, z, μeff)
        ϕ += _kernel_fluence_DA_slab_CW(m, s, zb, z0, ρ, z, μeff)
    end
    ϕ += _kernel_fluence_DA_slab_CW(0, s, zb, z0, ρ, z, μeff)

    return ϕ / (4 * π * D)
end

@inline function _kernel_fluence_DA_slab_CW(m, s, zb, z0, ρ, z, μeff)
    tmp1 = 2 * m * (s + 2 * zb)
    zmp = tmp1 + z0
    zmm = tmp1 - 2 * zb - z0

    tmp2 = sqrt(ρ^2 + (z - zmp)^2)
    tmp3 = sqrt(ρ^2 + (z - zmm)^2)

    ϕ = exp(-μeff * tmp2) / tmp2
    ϕ -= exp(-μeff * tmp3) / tmp3
    return ϕ
end

#--------------------------------------
# Time-Domain Fluence
#--------------------------------------
"""
    fluence_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 10)

Compute the time-domain fluence from a slab geometry (x, y -> inf, z -> finite).
The fluence is calculated from an infinite summation that converges rather rapidly.
Generally only 5 terms are needed (xs = 5), but more terms will be required if ρ >> s.
At long time values (low fluence values), the solution struggles with limited precision.
Using double floats or BigFloats is advised at these times.

# Arguments
- `t`: the time vector (ns).
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `s`: the thickness (z-depth) of the slab (cm)
- `z`: the z-depth coordinate (cm)
- `xs`: the number of sources to compute in the series

# Examples
julia> `fluence_DA_slab_TD(0.1, 1.0, 0.1, 10.0, s = 40.0)`
"""
function fluence_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t-> _kernel_fluence_DA_slab_TD(params.D, params.ν, t, ρ, μa, s, params.zb, params.z0, z, xs), t)
end

"""
    fluence_DA_slab_TD(t, data::DiffusionParameters)

Wrapper to `fluence_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 10)`

# Examples
```
julia> data = DAslab() # use structure to generate inputs
julia> fluence_DA_slab_TD(data) # then call the function
```
"""
function fluence_DA_slab_TD(data::DiffusionParameters)
    return map(t -> _kernel_fluence_DA_slab_TD(data.D, data.ν, t, data.ρ, data.μa, data.s, data.zb, data.z0, data.z, data.xs), data.t)
end

function _kernel_fluence_DA_slab_TD(D, ν, t, ρ, μa, s, zb, z0, z, xs)
    tmp1 = 4 * D * ν * t
    ϕ = ν * exp(-(ρ^2 / tmp1) - μa * ν * t)
    ϕ /= (π * tmp1)^(3/2)
    ϕ_tmp = zero(eltype(ϕ))

    for m in xs:-1:1
        ϕ_tmp += _images_fluence_DA_slab_TD(-m, s, zb, z0, z, tmp1)
        ϕ_tmp += _images_fluence_DA_slab_TD(m, s, zb, z0, z, tmp1)
    end
    ϕ_tmp += _images_fluence_DA_slab_TD(0, s, zb, z0, z, tmp1)
    ϕ *= ϕ_tmp
end

@inline function _images_fluence_DA_slab_TD(m, s, zb, z0, z, tmp1)
    tmp2 = 2 * m * (s + 2 * zb)
    zmp = tmp2 + z0
    zmm = tmp2 - 2 * zb - z0
    return exp(-((z - zmp)^2 / tmp1)) - exp(-((z - zmm)^2 / tmp1))
end

#--------------------------------------
# Time-Domain Flux
#--------------------------------------
"""
    flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 15)

    Compute the time-domain flux, D*∂ϕ(t)/∂z for z = 0 and -D*∂ϕ(t)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite).
    If z != 0 or s will default to compute D*∂ϕ(t)/∂z for the z given.

# Arguments
- `t`: the time vector (ns).
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `s`: the thickness (z-depth) of the slab (cm)
- `z`: the z-depth coordinate (cm) must be equal
- `xs`: the number of sources to compute in the series

# Examples
julia> `flux_DA_slab_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0)`
"""
function flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0, xs = 15)
    if z == zero(eltype(z))
        return _refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s, xs)
    elseif z == s
        return _trans_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s, xs)
    else
        return D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, z = dz, xs = xs), z)
    end
end

"""
    flux_DA_slab_TD(t, data::DiffusionParameters)

Wrapper to `flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0, xs = 15)`

# Examples
```
julia> data = DAslab(ρ = 1.0) # use structure to generate inputs
julia> flux_DA_slab_TD(data) # then call the function
```
"""
function flux_DA_slab_TD(data)
    return flux_DA_slab_TD(data.t, data.ρ, data.μa, data.μsp, n_ext = data.n_ext, n_med = data.n_med, z = data.z, s = data.s)
end

"""
    _refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15)

Compute the time-domain flux or reflectance at the top surface from a slab geometry (x,y->inf, z-> finite).

"""
function _refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s, xs)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t -> _kernel_refl_DA_slab_TD(ρ, params.D, params.ν, t, μa, xs, s, params.zb, params.z0), t)
end

@inline function _kernel_refl_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0)
    Rt1 = -exp(-(ρ^2 / (4 * D * ν * t)) - μa * ν * t) / (2 * (4 * π * D * ν)^(3/2) * t^(5/2))
    Rt2 = zero(eltype(Rt1))
	for m in -xs:xs
	    z3m = -2 * m * s - 4 * m * zb - z0
	    z4m = -2 * m * s - (4m - 2) * zb + z0
	    Rt2 += z3m * exp(-(z3m^2 / (4 * D * ν * t))) - z4m * exp(-(z4m^2 / (4 * D * ν * t)))
    end

	return Rt1 * Rt2
end

"""
    trans_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15)

Compute the time-domain transmittance (flux) from a slab geometry (x,y -> inf, z -> finite) with Eqn. 39 from Contini 97.
"""
function _trans_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s, xs)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t -> _kernel_trans_DA_slab_TD(ρ, params.D, params.ν, t, μa, xs, s, params.zb, params.z0), t)
end

@inline function _kernel_trans_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0)
    Rt1 = exp(-(ρ^2 / (4 * D * ν * t)) - μa * ν * t) / (2 * (4 * π * D * ν)^(3/2) * t^(5/2))
    Rt2 = zero(eltype(Rt1))
	for m in -xs:xs
	    z1m = (1 - 2 * m) * s - 4 * m * zb - z0
		z2m = (1 - 2 * m) * s - (4 * m - 2) * zb + z0
	    Rt2 += z1m * exp(-(z1m^2 / (4 * D * ν * t))) - z2m * exp(-(z2m^2 / (4 * D * ν * t)))
    end

	return Rt1 * Rt2
end

#--------------------------------------
# Steady-State Flux
#--------------------------------------
"""
    fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10)

Compute the steady-state flux, D*∂ϕ(ρ)/∂z for z = 0 and -D*∂ϕ(ρ)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite).
If z != 0 or s will default to compute D*∂ϕ(ρ)/∂z for the z given.

# Arguments
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `s`: the thickness (z-depth) of the slab (cm)
- `z`: the z-depth within slab (cm)
- `xs`: the number of sources to compute in the series

# Examples
```
julia> flux_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0)
```
"""
function flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    if z == zero(eltype(z))
        return params.D * ForwardDiff.derivative(dz -> fluence_DA_slab_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, s = s, z = dz, xs = xs), z)
    elseif z == s
        return -params.D * ForwardDiff.derivative(dz -> fluence_DA_slab_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, s = s, z = dz, xs = xs), z)
    else
        return params.D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, z = dz, xs = xs), z)
    end
end

"""
    flux_DA_slab_CW(data::DiffusionParameters)

Wrapper to `flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0)`

# Examples
```
julia> data = DAslab(ρ = 1.0) # use structure to generate inputs
julia> flux_DA_slab_CW(data) # then call the function
```
"""
function flux_DA_slab_CW(data)
    return flux_DA_slab_CW(data.ρ, data.μa, data.μsp, n_ext = data.n_ext, n_med = data.n_med, s = data.s, z = data.z, xs = data.xs)
end

max_issues_repo_name: heltonmc/LightPropagation


text:
using Base: depwarn

# issue #23: swap argument order
# Deprecated in ImageBinarization v0.3
function binarize(f::AbstractImageBinarizationAlgorithm, img::AbstractArray, args...; kwargs...)
    depwarn("binarize(alg, img) is deprecated, use binarize(img, alg) instead", :binarize)
    binarize(img, f, args...; kwargs...)
end

# unexport recommend_size
# Deprecated in ImageBinarization v0.3
function recommend_size(img)
    depwarn("recommend_size is deprecated, use AdaptiveThreshold(img) to directly generate the algorithm", :recommend_size)
    return round(Int, mean(size(img))/8)
end


text:
function ssd_write(dest::AbstractString, s::Simulation)
    if isfile(dest) @warn "Destination `$dest` already exists. Overwriting..." end
    HDF5.h5open(dest, "w") do h5f
        LegendHDF5IO.writedata( h5f, "SSD_Simulation", NamedTuple(s)  )
    end
end

function ssd_read(src::AbstractString, ::Type{Simulation})
    HDF5.h5open(src, "r") do h5f
        Simulation(LegendHDF5IO.readdata(h5f, "SSD_Simulation"));
    end
end


text:
"""
$(TYPEDEF)

A structure used to store information related to the bounds assigned to each
variable.

$(TYPEDFIELDS)
"""
Base.@kwdef struct VariableInfo{T<:AbstractFloat}
    "Is the variable integer valued?"
    is_integer::Bool                                            = false
    "Boolean indicating whether finite lower bound exists."
    has_lower_bound::Bool                                       = false
    "Boolean indicating whether finite upper bound exists."
    has_upper_bound::Bool                                       = false
    "Boolean indicating variable is fixed to a finite value."
    is_fixed::Bool                                              = false
    "Indicates that no constraints have been set"
    has_constraints::Bool                                       = false
    "Lower bounds. May be -Inf."
    lower_bound::T                                              = -Inf
    "Upper bounds. May be Inf."
    upper_bound::T                                              = Inf
    ""
end
_is_integer(x::VariableInfo) = x.is_integer
_has_lower_bound(x::VariableInfo) = x.has_lower_bound
_has_upper_bound(x::VariableInfo) = x.has_upper_bound
_lower_bound(x::VariableInfo{T}) where {T <: AbstractFloat} = x.lower_bound
_upper_bound(x::VariableInfo{T}) where {T <: AbstractFloat} = x.upper_bound

_is_fixed(x::VariableInfo) = x.is_fixed
function _is_less_than(x::VariableInfo)
    flag = x.has_upper_bound
    flag &= !x.has_lower_bound
    flag &= !x.is_integer
    return flag
end
function _is_greater_than(x::VariableInfo)
    flag = x.has_lower_bound
    flag &= !x.has_upper_bound
    flag &= !x.is_integer
    return flag
end
function _is_int_interval(x::VariableInfo)
    flag = x.has_lower_bound
    flag &= x.has_upper_bound
    flag &= x.is_integer
    return flag
end
function _is_real_interval(x::VariableInfo)
    flag = x.has_lower_bound
    flag &= x.has_upper_bound
    flag &= !x.is_integer
    return flag
end
function _is_zero_one(x::VariableInfo{T}) where {T <: AbstractFloat}
    flag = iszero(x.lower_bound)
    flag &= isone(x.upper_bound)
    flag &= x.is_integer
    return flag
end

mid(x::VariableInfo{T}) where {T <: AbstractFloat} = 0.5*(_upper_bound(x) - _lower_bound(x))
empty_variable_info(::Type{T}) where T = VariableInfo{T}(lower_bound = Inf,
                                                         upper_bound = -Inf)

Base.isempty(v::VariableInfo{T}) where {T <: AbstractFloat} = _lower_bound(v) > _upper_bound(v)
function check_isempty(l, u, b)
    flag = l < u
    if b
        flag &= (l <= 0.0) | (u >= 1.0)
    end
    return !flag
end

function VariableInfo(::Type{T}, ::ZO) where {T <: AbstractFloat}
    return VariableInfo{T}(is_integer = true,
                           has_lower_bound = true,
                           has_upper_bound = true,
                           has_constraints = true,
                           lower_bound = zero(T),
                           upper_bound = one(T))
end

function VariableInfo(it::MOI.Interval{T}) where {T <: AbstractFloat}
    VariableInfo{T}(has_lower_bound = !isinf(it.lower),
                    has_upper_bound = !isinf(it.upper),
                    has_constraints = !isinf(it.lower) | !isinf(it.upper),
                    is_fixed = it.lower == it.upper,
                    lower_bound = it.lower,
                    upper_bound = it.upper)
end

function VariableInfo(v::VariableInfo{T}, ::ZO) where {T <: AbstractFloat}
    isempty(v) && (return v)
    l = max(0.0, _lower_bound(v))
    u = min(1.0, _upper_bound(v))
    if check_isempty(l, u, _is_integer(v))
        return empty_variable_info()
    end
    return VariableInfo{T}(is_integer = true,
                        has_lower_bound = true,
                        has_upper_bound = true,
                        has_constraints = true,
                        is_fixed = l == u,
                        lower_bound = l,
                        upper_bound = u)
end

function VariableInfo(v::VariableInfo{T}, it::MOI.Interval{T}) where {T <: AbstractFloat}
    isempty(v) && (return v)
    l = max(it.lower, _lower_bound(v))
    u = max(it.upper, _upper_bound(v))
    if check_isempty(l, u, _is_integer(v))
        return empty_variable_info(T)
    end
    return VariableInfo(is_integer = _is_integer(v),
                        has_lower_bound = !isinf(l),
                        has_upper_bound = !isinf(u),
                        has_constraints = !isinf(l) | !isinf(u),
                        is_fixed = l == u,
                        lower_bound = l,
                        upper_bound = u)
end

_ZO(v::VariableInfo)  = MOI.ZeroOne()
_ET(v::VariableInfo{T}) where {T <: AbstractFloat}  = MOI.EqualTo{T}(v.lower_bound)
_IT(v::VariableInfo{T}) where {T <: AbstractFloat}  = MOI.Interval{T}(v.lower_bound)
_GT(v::VariableInfo{T}) where {T <: AbstractFloat}  = MOI.GreaterThan{T}(v.lower_bound)
_LT(v::VariableInfo{T}) where {T <: AbstractFloat}  = MOI.LessThan{T}(v.upper_bound)
_INT(v::VariableInfo{T}) where {T <: AbstractFloat} = MOI.Semiinteger{T}(v.lower_bound, v.upper_bound)


text:
"""
Compiled templating language for Genie.
"""
module Flax

using Revise, Gumbo, SHA, Reexport, JSON, OrderedCollections, Markdown, YAML
using Genie, Genie.Loggers, Genie.Configuration
@reexport using HttpCommon

export HTMLString, JSONString, JSString
export doctype, var_dump, @vars, @yield, el
export foreachvar, @foreach, foreachstr
export partial, template

import Base.string
import Base.show
import Base.==
import Base.hash

const NORMAL_ELEMENTS = [ :html, :head, :body, :title, :style, :address, :article, :aside, :footer,
                          :header, :h1, :h2, :h3, :h4, :h5, :h6, :hgroup, :nav, :section,
                          :dd, :div, :d, :dl, :dt, :figcaption, :figure, :li, :main, :ol, :p, :pre, :ul, :span,
                          :a, :abbr, :b, :bdi, :bdo, :cite, :code, :data, :dfn, :em, :i, :kbd, :mark,
                          :q, :rp, :rt, :rtc, :ruby, :s, :samp, :small, :spam, :strong, :sub, :sup, :time,
                          :u, :var, :wrb, :audio, :map, :void, :embed, :object, :canvas, :noscript, :script,
                          :del, :ins, :caption, :col, :colgroup, :table, :tbody, :td, :tfoot, :th, :thead, :tr,
                          :button, :datalist, :fieldset, :form, :label, :legend, :meter, :optgroup, :option,
                          :output, :progress, :select, :textarea, :details, :dialog, :menu, :menuitem, :summary,
                          :slot, :template, :blockquote, :center]
const VOID_ELEMENTS   = [:base, :link, :meta, :hr, :br, :area, :img, :track, :param, :source, :input]
const BOOL_ATTRIBUTES = [:checked, :disabled, :selected]

const FILE_EXT      = [".flax.jl", "html.jl"]
const TEMPLATE_EXT  = [".flax.html", ".jl.html"]
const JSON_FILE_EXT = ".json.jl"
const MARKDOWN_FILE_EXT = [".md", ".jl.md"]

const SUPPORTED_HTML_OUTPUT_FILE_FORMATS = TEMPLATE_EXT

const HTMLString = String
const JSONString = String

const BUILD_NAME    = "FlaxViews"

const MD_SEPARATOR_START = "---\n"
const MD_SEPARATOR_END   = "---\n"

init_task_local_storage() = (haskey(task_local_storage(), :__vars) || task_local_storage(:__vars, Dict{Symbol,Any}()))
init_task_local_storage()

task_local_storage(:__yield, "")

"""
    prepare_template(s::String)
    prepare_template{T}(v::Vector{T})

Cleans up the template before rendering (ex by removing empty nodes).
"""
@inline function prepare_template(s::String) :: String
  s
end
@inline function prepare_template(v::Vector{T})::String where {T}
  filter!(v) do (x)
    ! isa(x, Nothing)
  end

  join(v)
end

"""
    attributes(attrs::Vector{Pair{Symbol,String}} = Vector{Pair{Symbol,String}}()) :: Vector{String}

Parses HTML attributes.
"""
function attributes(attrs::Vector{Pair{Symbol,Any}} = Vector{Pair{Symbol,Any}}()) :: String
  a = IOBuffer()

  for (k,v) in attrs
    sk = string(k)
    sk == "typ" && (k = "type")
    startswith(sk, "_") && (k = sk = sk[2:end])
    k = replace(sk, "_"=>"-")

    print(a, "$(k)=\"$(v)\" ")
  end

  String(take!(a))
end

"""
    normalize_element(elem::String)

Cleans up problematic characters or DOM elements.
"""
@inline function normalize_element(elem::String)
  elem == "d" && (elem = "div")
  replace(string(lowercase(elem)), "_"=>"-")
end

"""
    normal_element(f::Function, elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString

Generates a HTML element in the form <...></...>
"""
@inline function normal_element(f::Function, elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  normal_element(f(), elem, attrs...)
end
@inline function normal_element(children::Union{String,Vector{String}}, elem::String, attrs::Pair{Symbol,Any}) :: HTMLString
  normal_element(children, elem, Pair{Symbol,Any}[attrs])
end
@inline function normal_element(children::Union{String,Vector{String}}, elem::String, attrs...) :: HTMLString
  normal_element(children, elem, Pair{Symbol,Any}[attrs...])
end
@inline function normal_element(children::Union{String,Vector{String}}, elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  children = join(children)
  elem = normalize_element(elem)

  string("<", elem, " ", attributes(attrs), ">", prepare_template(children), "</", elem, ">")
end
@inline function normal_element(elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  normal_element("", elem, attrs...)
end
@inline function normal_element(elems::Vector, elem::String, attrs...) :: HTMLString
  io = IOBuffer()

  for e in elems
    e === nothing && continue

    if isa(e, Function)
      print(io, e(), "\n")
    else
      print(io, e, "\n")
    end
  end

  normal_element(String(take!(io)), elem, attrs...)
end
@inline function normal_element(_::Nothing, __::Any) :: HTMLString
  ""
end

"""
    void_element(elem::String, attrs::Vector{Pair{Symbol,String}} = Vector{Pair{Symbol,String}}()) :: HTMLString

Generates a void HTML element in the form <...>
"""
@inline function void_element(elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  string("<", normalize_element(elem), " ", attributes(attrs), ">")
end

"""
    skip_element(f::Function) :: HTMLString
    skip_element() :: HTMLString

Cleans up empty elements.
"""
@inline function skip_element(f::Function) :: HTMLString
  "$(prepare_template(f()))"
end
@inline function skip_element() :: HTMLString
  ""
end

"""
    partial(path::String; mod::Module = @__MODULE__, vars...) :: String

Renders (includes) a view partial within a larger view or layout file.
"""
function partial(path::String; mod::Module = @__MODULE__, vars...) :: String
  for (k,v) in vars
    if k == :context && isa(v, Module)
      (mod = v)
      continue
    end

    try
      task_local_storage(:__vars)[k] = v
    catch
      init_task_local_storage()
      task_local_storage(:__vars)[k] = v
    end
  end

  template(path, partial = true, mod = mod)
end

"""
    view_file_info(path::String) :: Tuple{String,String}

Extracts path and extension info about a file
"""
function view_file_info(path::String) :: Tuple{String,String}
  _path, _extension = "", ""

  if isfile(path)
    _path, _extension = relpath(path), "." * split(path, ".", limit = 2)[end]
  else
    for file_extension in SUPPORTED_HTML_OUTPUT_FILE_FORMATS
      if isfile(path * file_extension)
        _path, _extension = path * file_extension, file_extension
        break
      end
    end
  end

  _path, _extension
end

"""
    include_markdown(path::String; mod::Module = @__MODULE__)

Includes and renders a markdown view file
"""
function include_markdown(path::String; mod::Module = @__MODULE__)
  md = read(path, String)

  vars_injection = ""

  if startswith(md, MD_SEPARATOR_START)
    close_sep_pos = findfirst(MD_SEPARATOR_END, md[length(MD_SEPARATOR_START)+1:end])
    metadata = md[length(MD_SEPARATOR_START)+1:close_sep_pos[end]] |> YAML.load

    vars_injection = ""
    for (k,v) in metadata
      task_local_storage(:__vars)[Symbol(k)] = v
      vars_injection *= """@vars($(repr(Symbol(k)))) = $(v)\n"""
    end

    md = replace(md[close_sep_pos[end]+length(MD_SEPARATOR_END)+1:end], "\"\"\""=>"\\\"\\\"\\\"")
  end

  content = string( "\"\"\"", md, "\"\"\"")

  vars_injection, (include_string(mod, content) |> Markdown.parse |> Markdown.html)
end

"""
    get_template(path::String; partial::Bool = true, mod::Module = @__MODULE__) :: Function

Resolves the inclusion and rendering of a template file
"""
function get_template(path::String; partial::Bool = true, mod::Module = @__MODULE__) :: Function
  orig_path = path

  path, extension = view_file_info(path)

  isfile(path) || error("Template file $orig_path does not exist")

  extension in FILE_EXT && return (() -> Base.include(mod, path))

  f_name = Symbol(function_name(path))
  f_path = joinpath(Genie.BUILD_PATH, BUILD_NAME, m_name(path) * ".jl")
  f_stale = build_is_stale(path, f_path)

  if f_stale || ! isdefined(mod, f_name)
    content = if extension in MARKDOWN_FILE_EXT
      vars_injection, md = include_markdown(path, mod = mod)
      string_to_flax(md, partial = partial, f_name = f_name, prepend = vars_injection)
    else
      html_to_flax(path, partial = partial)
    end

    f_stale && build_module(content, path)

    return Base.include(mod, joinpath(Genie.BUILD_PATH, BUILD_NAME, m_name(path) * ".jl"))
  end

  getfield(mod, f_name)
end

"""
    parse_view(data::String; partial = false, mod::Module = @__MODULE__) :: Function

Parses a view file, returning a rendering function. If necessary, the function is JIT-compiled, persisted and loaded into memory.
"""
function parse_view(data::String; partial = false, mod::Module = @__MODULE__) :: Function
  path = "Flax_" * string(hash(data))

  func_name = function_name(data) |> Symbol
  f_path = joinpath(Genie.BUILD_PATH, BUILD_NAME, m_name(path) * ".jl")
  f_stale = build_is_stale(f_path, f_path)

  if f_stale || ! isdefined(mod, func_name)
    f_stale && build_module(string_to_flax(data, partial = partial), path)

    return Base.include(mod, joinpath(Genie.BUILD_PATH, BUILD_NAME, m_name(path) * ".jl"))
  end

  getfield(mod, func_name)
end

"""
    template(path::String; partial::Bool = true, mod::Module = @__MODULE__) :: String

Renders a template file.
"""
@inline function template(path::String; partial::Bool = true, mod::Module = @__MODULE__) :: String
  get_template(path, partial = partial, mod = mod) |> Base.invokelatest
end

"""
    build_is_stale(file_path::String, build_path::String) :: Bool

Checks if the view template has been changed since the last time the template was compiled.
"""
function build_is_stale(file_path::String, build_path::String) :: Bool
  isfile(file_path) || return true

  file_mtime = stat(file_path).mtime
  build_mtime = stat(build_path).mtime
  status = file_mtime > build_mtime

  Genie.config.log_views && status && log("🚨  Flax view $file_path build $build_path is stale")

  status
end

"""
    register_vars(vars...) :: Nothing

Loads the rendering vars into the task's scope
"""
@inline function register_vars(vars...) :: Nothing
  init_task_local_storage()
  task_local_storage(:__vars, merge(Dict{Symbol,Any}(vars), task_local_storage(:__vars)))

  nothing
end

"""
    html_renderer(resource::Union{Symbol,String}, action::Union{Symbol,String}; layout::Union{Symbol,String} = Genie.config.renderer_default_layout_file, mod::Module = @__MODULE__, vars...) :: Function

Renders data as HTML
"""
function html_renderer(resource::Union{Symbol,String}, action::Union{Symbol,String}; layout::Union{Symbol,String} = Genie.config.renderer_default_layout_file, mod::Module = @__MODULE__, vars...) :: Function
  register_vars(vars...)
  task_local_storage(:__yield, get_template(joinpath(Genie.RESOURCES_PATH, string(resource), Genie.VIEWS_FOLDER, string(action)), partial = true, mod = mod) |> Base.invokelatest)

  layout = Base.get(task_local_storage(:__vars), :layout, layout)

  get_template(joinpath(Genie.APP_PATH, Genie.LAYOUTS_FOLDER, string(layout)), partial = false, mod = mod)
end
function html_renderer(data::String; mod::Module = @__MODULE__, layout::Union{Symbol,String,Nothing} = nothing, vars...) :: Function
  register_vars(vars...)

  if layout != nothing
    task_local_storage(:__yield, parse_view(data, partial = true, mod = mod))
    get_template(joinpath(Genie.APP_PATH, Genie.LAYOUTS_FOLDER, string(layout)), partial = false, mod = mod)
  else
    parse_view(data, partial = false, mod = mod)
  end
end

"""
    json_renderer(resource::Union{Symbol,String}, action::Union{Symbol,String}; mod::Module = @__MODULE__, vars...) :: Function

Renders data as JSON
"""
@inline function json_renderer(resource::Union{Symbol,String}, action::Union{Symbol,String}; mod::Module = @__MODULE__, vars...) :: Function
  register_vars(vars...)

    () -> (Base.include(mod, joinpath(Genie.RESOURCES_PATH, string(resource), Genie.VIEWS_FOLDER, string(action) * JSON_FILE_EXT)) |> JSON.json)
end

"""
    function_name(file_path::String)

Generates function name for generated Flax views.
"""
@inline function function_name(file_path::String) :: String
  file_path = relpath(file_path)
  "func_$(sha1(file_path) |> bytes2hex)"
end

"""
    m_name(file_path::String)

Generates module name for generated Flax views.
"""
@inline function m_name(file_path::String) :: String
  file_path = relpath(file_path)
  "$(sha1(file_path) |> bytes2hex)"
end

"""
    html_to_flax(file_path::String; partial = true) :: String

Converts a HTML document to a Flax document.
"""
@inline function html_to_flax(file_path::String; partial = true) :: String
  to_flax(file_path, parse_template, partial = partial)
end

"""
    string_to_flax(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts string view data to Flax code
"""
@inline function string_to_flax(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String
  to_flax(content, parse_string, partial = partial, f_name = f_name, prepend = prepend)
end

"""
    to_flax(input::String, f::Function; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts an input file to Flax code
"""
@inline function to_flax(input::String, f::Function; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String
  f_name = f_name === nothing ? function_name(input) : f_name

  string("function $(f_name)() \n",
          prepend,
          f(input, partial = partial),
          "\nend \n")
end

"""
    build_module(content::String, path::String) :: Bool

Persists compiled Flax view data to file
"""
function build_module(content::String, path::String) :: Bool
  module_path = joinpath(Genie.BUILD_PATH, BUILD_NAME, m_name(path) * ".jl")
  isdir(joinpath(Genie.BUILD_PATH, BUILD_NAME)) || mkpath(joinpath(Genie.BUILD_PATH, BUILD_NAME))
  open(module_path, "w") do io
    write(io, "# $path \n\n")
    write(io, content)
  end

  true
end

"""
    read_template_file(file_path::String) :: String

Reads `file_path` template from disk.
"""
function read_template_file(file_path::String) :: String
  io = IOBuffer()
  open(file_path) do f
    for line in enumerate(eachline(f))
      print(io, parse_tags(line), "\n")
    end
  end

  String(take!(io))
end

"""
    parse_template(file_path::String; partial = true) :: String

Parses a HTML file into Flax code.
"""
@inline function parse_template(file_path::String; partial = true) :: String
  parse(read_template_file(file_path), partial = partial)
end

"""
    parse_string(data::String; partial = true) :: String

Parses a HTML string into Flax code.
"""
@inline function parse_string(data::String; partial = true) :: String
  parse(parse_tags(data), partial = partial)
end

@inline function parse(input::String; partial = true) :: String
  parse_tree(Gumbo.parsehtml(input).root, "", 0, partial = partial)
end

"""
    parse_tree(elem, output, depth; partial = true) :: String

Parses a Gumbo tree structure into a `string` of Flax code.
"""
function parse_tree(elem::Union{HTMLElement,HTMLText}, output::String = "", depth::Int = 0; partial = true) :: String
  io = IOBuffer()

  if isa(elem, HTMLElement)
    tag_name = replace(lowercase(string(tag(elem))), "-"=>"_")

    if Genie.config.flax_autoregister_webcomponents && ! isdefined(@__MODULE__, Symbol(tag_name))
      log("Autoregistering HTML element $tag_name", :debug)

      register_element(Symbol(tag_name))
      print(io, "Genie.Flax.register_element(Symbol(\"$tag_name\")) \n")
    end

    invalid_tag = partial && (tag_name == "html" || tag_name == "head" || tag_name == "body")

    if tag_name == "script" && in("type", collect(keys(attrs(elem))))
      if attrs(elem)["type"] == "julia/eval"
        if ! isempty(children(elem))
          print(io, repeat("\t", depth), string(children(elem)[1].text), "\n")
        end
      end

    else
      print(io, repeat("\t", depth), ( ! invalid_tag ? "Flax.$(tag_name)(" : "Flax.skip_element(" ))

      attributes = IOBuffer()
      for (k,v) in attrs(elem)
        x = v

        if startswith(v, "<\$") && endswith(v, "\$>")
          v = (replace(replace(replace(v, "<\$"=>""), "\$>"=>""), "'"=>"\"") |> strip)
          x = v
          v = "\$($v)"
        end

        if in(Symbol(lowercase(k)), BOOL_ATTRIBUTES)
          if x == true || x == "true" || x == :true || x == ":true" || x == ""
            # push!(attributes, "$k = \"$k\"") # boolean attributes can have the same value as the attribute -- or be empty
            print(attributes, "$k = \"$k\"", ", ") # boolean attributes can have the same value as the attribute -- or be empty
          end
        else
          print(attributes, """$(replace(lowercase(string(k)), "-"=>"_")) = "$v" """, ", ")
        end
      end

      attributes_string = String(take!(attributes))
      endswith(attributes_string, ", ") && (attributes_string = attributes_string[1:end-2])
      print(io, attributes_string, ") ")

      inner = ""
      if ! isempty(children(elem))
        children_count = size(children(elem))[1]

        print(io, "do;[\n")

        idx = 0
        for child in children(elem)
          idx += 1
          inner *= parse_tree(child, "", depth + 1, partial = partial)
          if idx < children_count
            if isa(child, HTMLText) ||
                ( isa(child, HTMLElement) && ( ! in("type", collect(keys(attrs(child)))) || ( in("type", collect(keys(attrs(child)))) && (attrs(child)["type"] != "julia/eval") ) ) )
                ! isempty(inner) && (inner = repeat("\t", depth) * inner * "\n")
            end
          end
        end
        isempty(inner) || (print(io, inner, "\n", repeat("\t", depth)))

        print(io, "]end\n")
      end
    end

  elseif isa(elem, HTMLText)
    content = replace(elem.text |> strip |> string, "\""=>"\\\"")
    print(io, repeat("\t", depth), "\"$(content)\"")
  end

  String(take!(io))
end

"""
    parse_tags(line::Tuple{Int64,String}, strip_close_tag = false) :: String

Parses special Flax tags.
"""
@inline function parse_tags(line::Tuple{Int64,String}) :: String
  parse_tags(line[2])
end
@inline function parse_tags(code::String) :: String
  code = replace(code, "<%"=>"""<script type="julia/eval">""")
  replace(code, "%>"=>"""</script>""")
end

"""
Outputs document's doctype.
"""
@inline function doctype(doctype::Symbol = :html) :: HTMLString
  "<!DOCTYPE $doctype>"
end

"""
Outputs document's doctype.
"""
@inline function doc(html::String) :: HTMLString
  doctype() * "\n" * html
end
@inline function doc(doctype::Symbol, html::String) :: HTMLString
  doctype(doctype) * "\n" * html
end

"""
    register_elements() :: Nothing

Generated functions that represent Flax functions definitions corresponding to HTML elements.
"""
@inline function register_elements() :: Nothing
  for elem in NORMAL_ELEMENTS
    register_normal_element(elem)
  end

  for elem in VOID_ELEMENTS
    register_void_element(elem)
  end

  nothing
end

@inline function register_element(elem::Symbol, elem_type::Symbol = :normal) :: Nothing
  elem_type == :normal ? register_normal_element(elem) : register_void_element(elem)
end

function register_normal_element(elem::Symbol) :: Nothing
  Core.eval(@__MODULE__, """
    function $elem(f::Function; attrs...) :: HTMLString
      \"\"\"\$(normal_element(f, "$(string(elem))", Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)
  Core.eval(@__MODULE__, """
    function $elem(children::Union{String,Vector{String}} = ""; attrs...) :: HTMLString
      \"\"\"\$(normal_element(children, "$(string(elem))", Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  nothing
end

function register_void_element(elem::Symbol) :: Nothing
  Core.eval(@__MODULE__, """
    function $elem(; attrs...) :: HTMLString
      \"\"\"\$(void_element("$(string(elem))", Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  nothing
end

push!(LOAD_PATH,  abspath(Genie.HELPERS_PATH))

"""
    @foreach(f, arr)

Iterates over the `arr` Array and applies function `f` for each element.
The results of each iteration are concatenated and the final string is returned.
"""
macro foreach(f, arr)
  quote
    isempty($(esc(arr))) && return ""

    mapreduce(*, $(esc(arr))) do _s
      $f(_s) * "\n"
    end
  end
end

function foreachstr(f, arr)
  isempty(arr) && return ""

  mapreduce(*, arr) do _s
    f(_s)
  end
end

"""
    foreachvar(f::Function, key::Symbol, v::Vector) :: String

Utility function for looping over a vector `v` in the view layer.
"""
function foreachvar(f::Function, key::Symbol, v::Vector) :: String
  isempty(v) && return ""

  output = mapreduce(*, v) do (value)
    vars = task_local_storage(:__vars)
    vars[key] = value
    task_local_storage(:__vars, vars)

    f(value)
  end

  vars = task_local_storage(:__vars)
  delete!(vars, key)
  task_local_storage(:__vars, vars)

  output
end

register_elements()

"""
    var_dump(var, html = true) :: String

Utility function for dumping a variable into the view.
"""
function var_dump(var, html = true) :: String
  iobuffer = IOBuffer()
  show(iobuffer, var)
  content = String(take!(iobuffer))

  html ? replace(replace("<code>$content</code>", "\n"=>"<br>"), " "=>" ") : content
end

"""
    @vars

Utility macro for accessing view vars
"""
macro vars()
  :(task_local_storage(:__vars))
end

"""
    @vars(key)

Utility macro for accessing view vars stored under `key`
"""
macro vars(key)
  :(task_local_storage(:__vars)[$key])
end

"""
    @vars(key, value)

Utility macro for setting a new view var, as `key` => `value`
"""
macro vars(key, value)
  quote
    try
      task_local_storage(:__vars)[$key] = $(esc(value))
    catch
      init_task_local_storage()
      task_local_storage(:__vars)[$key] = $(esc(value))
    end
  end
end

"""
    @yield

Outputs the rendering of the view within the template.
"""
macro yield()
  quote
    try
      task_local_storage(:__yield)
    catch
      task_local_storage(:__yield, "")
    end
  end
end
macro yield(value)
  :(task_local_storage(:__yield, $value))
end

function el(; vars...)
  OrderedDict(vars)
end

"""
    prepare_build() :: Bool

Sets up the build folder and the build module file for generating the compiled views.
"""
function prepare_build(subfolder) :: Bool
  build_path = joinpath(Genie.BUILD_PATH, subfolder)
  @ifdev rm(build_path, force = true, recursive = true)
  if ! isdir(build_path)
    log("Creating build folder at $(build_path)", :info)
    mkpath(build_path)
  end

  true
end

"""
    create_build_folders()

Sets up build folders.
"""
@inline function create_build_folders()
  prepare_build(BUILD_NAME)
end

end


text:
using FiniteDifferences
import NbodyGradient: set_state!, zero_out!, amatrix

@testset "Transit Series" begin
    N = 3
    t0 = 7257.93115525 - 7300.0 - 0.5 # Initialize IC before first transit
    h = 0.04
    itime = 10.0
    tmax = itime

    # Setup initial conditions:
    elements = readdlm("elements.txt", ',')[1:N,:]
    elements[2:end,3] .-= 7300.0
    elements[:,end] .= 0.0 # Set Ωs to 0
    elements[2,1] *= 10.0
    elements[3,1] *= 10.0
    ic = ElementsIC(t0, N, elements)

    # Calculate transit times and b,vsky
    function calc_times(h;grad=false)
        intr = Integrator(ahl21!, h, 0.0, tmax)
        s = State(ic)
        ttbv = TransitParameters(itime, ic)
        intr(s, ttbv;grad=grad)
        return ttbv
    end

    # Calculate b,vsky for specified times
    function calc_pd(h, times;grad=false)
        intr = Integrator(ahl21!, h, 0.0, tmax)
        s = State(ic)
        ts = TransitSeries(times, ic)
        intr(s, ts; grad=grad)
        return ts
    end

    ttbv = calc_times(h; grad=true)
    mask = ttbv.ttbv[1,2,:] .!= 0.0
    ts = calc_pd(h, ttbv.ttbv[1,2,mask]; grad=true)

    tol = 1e-10
    # Test whether TransitSeries can reproduce b,vsky at transit time from TransitParameters
    @test all(abs.(ttbv.ttbv[2,2,mask] .- ts.vsky[2,:]) .< tol)
    @test all(abs.(ttbv.ttbv[3,2,mask] .- ts.bsky2[2,:]) .< tol)

    # Now make set of times to calc b and vsky and compare derivatives
    times = [t0, t0 + 1.0]
    ts = calc_pd(h, times; grad=true)

    # method for finite diff
    function calc_b(θ,i=1;b=true,times=ttbv.ttbv[1,2,mask])
        elements = reshape(θ,3,7)
        intr = Integrator(ahl21!, h, 0.0, 10.0)
        ic = ElementsIC(t0,N,elements)
        s = State(ic)
        pd = Photodynamics(times, ic)
        intr(s, pd; grad=true)
        if b
            return pd.bsky2[i,:]
        else
            return pd.vsky[i,:]
        end
    end

    # Now test derivatives
    function calc_finite_diff(h, t0, times, tmax, elements)
        # Now do finite difference with big float
        dq0 = big(1e-10)
        ic_big = ElementsIC(big(t0), N, big.(elements))
        elements_big = copy(ic_big.elements)
        s_big = State(ic_big)
        pdp = TransitSeries(big.(times), ic_big);
        pdm = TransitSeries(big.(times), ic_big);
        dbvde_num = zeros(BigFloat, size(pdp.dbvdelements));
        intr_big = Integrator(big(h), zero(BigFloat), big(tmax))

        for jq in 1:N
            for iq in 1:7
                zero_out!(pdp); zero_out!(pdm);
                ic_big.elements .= elements_big
                ic_big.m .= elements_big[:,1]
                if iq == 7; ivary = 1; else; ivary = iq + 1; end
                ic_big.elements[jq,ivary] += dq0
                if ivary == 1; ic_big.m[jq] += dq0; amatrix(ic_big); end # Masses don't update with elements array
                sp = State(ic_big);
                intr_big(sp, pdp; grad=false)
                ic_big.elements[jq,ivary] -= 2dq0
                if ivary == 1; ic_big.m[jq] -= 2dq0; amatrix(ic_big); end
                sm = State(ic_big)
                intr_big(sm, pdm; grad=false)
                for i in 2:N
                    for k in 1:ts.nt
                        # Compute double-sided derivative for more accuracy:
                        dbvde_num[1,i,k,iq,jq] = (pdp.vsky[i,k] - pdm.vsky[i,k]) / (2dq0)
                        dbvde_num[2,i,k,iq,jq] = (pdp.bsky2[i,k] - pdm.bsky2[i,k]) / (2dq0)
                    end
                end
            end
        end
        return dbvde_num
    end

    mask = zeros(Bool, size(ts.dbvdelements));
    for jq = 1:N
        for iq = 1:7
            if iq == 7; ivary = 1; else; ivary = iq + 1; end  # Shift mass variation to end
            for i = 2:N
                for k = 1:ts.nt
                    for itbv = 1:2
                        # Ignore inclination & longitude of nodes variations:
                        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
                            mask[2,i,k,iq,jq] = true
                        end
                    end
                end
            end
        end
    end

    dbvde_num = calc_finite_diff(h, t0, ts.times, tmax, elements)
    @test isapprox(asinh.(ts.dbvdelements[mask]), asinh.(dbvde_num[mask]); norm=maxabs)
    #mask_times = ttbv.ttbv[1,2,:] .!= 0.0
    #@test isapprox(asinh.(ts.dbvdelements[mask]), asinh.(ttbv.dtbvdelements[2:3,:,mask_times,:,:][mask]); norm=maxabs)
end


text:
function equals(obj::FieldUnivariateDerivative2Field, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_one(obj::FieldUnivariateDerivative2Field)
    return jcall(obj, "getOne", FieldUnivariateDerivative2, ())
end

function get_runtime_class(obj::FieldUnivariateDerivative2Field)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_univariate_derivative2_field(::Type{FieldUnivariateDerivative2Field}, arg0::Field)
    return jcall(FieldUnivariateDerivative2Field, "getUnivariateDerivative2Field", FieldUnivariateDerivative2Field, (Field,), arg0)
end

function get_zero(obj::FieldUnivariateDerivative2Field)
    return jcall(obj, "getZero", FieldUnivariateDerivative2, ())
end

function hash_code(obj::FieldUnivariateDerivative2Field)
    return jcall(obj, "hashCode", jint, ())
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using Lasso
using GLM
using LinearAlgebra

##

Ψ = rand(100_000, 200)
Y = rand(100_000)
qrΨ = qr(Ψ)
Q = Matrix(qrΨ.Q)
η = norm(Y - Q * (Q' * Y)) / norm(Y)

R = Matrix(qrΨ.R)
y = Q' * Y
η1 = norm(y) / norm(Y)

target = sqrt(1.2^2 - 1) * η

c = nothing

for n = 1:length(p.λ)
   ηn = norm(R * p.coefs[:,n] - y) / norm(Y)
   @show ηn
   if ηn < target
      global c = p.coefs[:,n]
      break
   end
end

@show length(c.nzind)
@show norm(Ψ * c - Y)/norm(Y)
@show 1.2 * η

## ------------------------------------------------------------------------


text:
"""
Содержит функции:
searchinrange(begs::Vector, ends::Vector, from, to) - поиск диапазона сегментов,пересекающихся с from-to

seg_in_range(segs::StructArray, from,to) - аналог searchinrange для StructArray

delete_seg(segs::StructArray, from, to) - удаление сегментов, попавших под диапазон rom-to

add_seg(segs::StructArray, newSeg::StructArray, mode) - добавление сегментов по 3-м режимам:"rewrite"(удаление всего, что пересекалось),
"simpleAdd"(тупо добавить и отсортировать), "raw"(слить,расширяя границы при пересечении)"

addRawSeg(seg1::StructArray,seg2::StructArray) - слияние seg1 с seg2 ,расширяя границы при пересечении

mergeIntersections(seg::StructArray) - удаление внутренних пересечений в seg

segCompare(seg1::StructArray, seg2::StructArray) - сравнение двух наборов на пересечение сегментов мд ними

segMinus(seg1::StructArray, seg2::StructArray) - удаление из seg1 сегментов, наложившихся на seg2. (в 2х версиях)
"""

using StructArrays
#сегмент - StructArray с обязательными ibeg, iend и необязательным набором признаков
#сегменты упорялочены по возрастанию
function searchinrange(begs::Vector{T}, ends::Vector{T}, from, to) where T
    i1 = searchsortedfirst(ends, from)
    i2 = searchsortedlast(begs, to)
    i1:i2 # соответствие с searchsorted
end

function seg_in_range(segs::StructArray, from,to)
    #определение сегментов, пересекающихся с диапазоном [from-to]
    s_range = searchinrange(segs.ibeg, segs.iend, from, to)
    segsInRange = segs[s_range]
    # segsInRange = filter(x->sign(x.ibeg-to) != sign(x.iend-from), LazyRows(segs))
    return segsInRange
end

"""
удаление сегментов, которые наложились на диапазон [from-to]
"""
function delete_seg(segs::StructArray, from, to)
    s_range = searchinrange(segs.ibeg, segs.iend, from, to)
    Ns = length(segs)
    i=1:Ns
    segs_reduced = StructArray(ibeg=Int[], iend=Int[], type = Int[])
    #segs_reduced =  StructArray{eltype(Int64)}(undef, 0) #так не пашет
    if !isempty(s_range)
    #if range - > ошибку дает, что небулевое значение в булевом контексте
    #результат выходил как вектор Struct-Arrays-ов, так что сделала так
        append!(segs_reduced,segs[1:s_range.start-1])
        append!(segs_reduced,segs[s_range.stop+1:end])
        #segs_reduced = segs[idTake]
    else
        segs_reduced = copy(segs) # !!! все массивы типа mutable - присваиваются и передаются в функции по ссылке
    end
    return segs_reduced
end

"""
добавление сегментов
mode = "rewrite"(удаление всего, что пересекалось),
"simpleAdd"(тупо добавить и отсортировать), "raw"(слить,расширяя границы при пересечении)
"""
function add_seg(segs::StructArray, newSeg::StructArray, mode)
    #seg_out = StructArray{eltype(Int)}(undef, 0) # вот так можно создать пустой массив
    #seg_out = StructArray(ibeg=typeof(newSeg.ibeg)(), iend=typeof(newSeg.iend)(), type = typeof(newSeg.type)())#type = typeof(newSeg[1].type)())
    seg_out = StructArray(ibeg=Int[], iend=Int[], type=Int[])
    if mode=="simpleAdd"

        NsAdd = length(newSeg)
        NsAdd > 0 || return copy(segs)

        i0 = i1 = 1
        for i=1:NsAdd
            #индекс первого бОльшего сегмента
            i1 = searchsortedfirst(segs.ibeg, newSeg[i].ibeg)
            temp_t_right = @view segs[i0:i1-1]
            append!(seg_out, temp_t_right) # ничего страшного, если пустой добавится - ничего не произойдет
            push!(seg_out, newSeg[i]) # push - для элементов, append - для векторов
            i0 = i1
        end
        temp_t_right = @view segs[i1:end]
        append!(seg_out, temp_t_right)

    elseif mode=="rewrite"
        NsAdd = length(newSeg)
        seg_out = copy(segs)
        for i = 1:NsAdd
            seg_out = delete_seg(seg_out, newSeg[i].ibeg, newSeg[i].iend)
        end
        seg_out = add_seg(seg_out, newSeg, "simpleAdd")
    elseif mode=="raw"
        seg_out = addRawSeg(segs,newSeg)
    end

    return seg_out
end

"""
Вспомогательные функции для addRawSeg()
"""
function addTempSeg(tempSeg,newSeg)
    if tempSeg.ibeg != 0
        if tempSeg.iend-tempSeg.ibeg > 0
            if !isa(newSeg.ibeg,UndefInitializer)
                # @info newSeg, tempSeg
                #так плохо, но я не знаю, как иначе
                ib = [newSeg.ibeg tempSeg.ibeg]
                ie = [newSeg.iend tempSeg.iend]
                tp = [newSeg.type tempSeg.type]
                newSeg = (ibeg = ib, iend = ie, type = tp)
            else
                newSeg = tempSeg
            end
        end
        tempSeg = (ibeg = 0, iend = 0, type = 0)
    end
    return newSeg, tempSeg
end

function sortSeg(seg::StructArray)
    return seg[sortperm(seg.ibeg)]
end

"""
Добавление сегментов из набора 2 в набор 1 с расширением сегментов при пересечении
"""
function addRawSeg(seg1::StructArray,seg2::StructArray)
    #проверка на внутренние пересечения. ИЛИ ЕЕ НАРУЖУ ВЫНЕСТИ??
    seg1 = mergeIntersections(seg1)
    seg2 = mergeIntersections(seg2)

    N1 = length(seg1)
    N2 = length(seg2)
    stop = false
    toDel1 = falses(N1)
    toDel2 = falses(N2)

    newSeg =  (ibeg = undef, iend = undef, type=undef) # StructArray(ibeg=Int[], iend=Int[], type=Int[])
    tempSeg = (ibeg = 0, iend = 0, type=0)
    i=j=1
    while !stop
        if seg1[i].iend < seg2[j].ibeg # Seg1 кончается раньше, чем начинается Seg2
            newSeg, tempSeg = addTempSeg(tempSeg, newSeg)
            if i<N1; i+=1;continue
            else stop = true end
        elseif seg2[j].iend < seg1[i].ibeg #Seg2 кончается раньше, чем начинается Seg1
            newSeg, tempSeg = addTempSeg(tempSeg, newSeg)
            if j<N2; j+=1;continue
            else stop = true end
        else #сегменты пересекаются
            #и имеют одинаковые типы
            equalTypes = seg1[i].type == seg2[j].type
            if equalTypes
                toDel1[i] = toDel2[j] = true
                if tempSeg.ibeg == 0
                    tempSeg = (ibeg = minimum([seg1[i].ibeg,seg2[j].ibeg]),
                    iend = maximum([seg1[i].iend,seg2[j].iend]),
                    type = seg1[i].type)
                else
                    tempSeg = (ibeg = tempSeg.ibeg,
                    iend = maximum([seg1[i].iend,seg2[j].iend]),
                    type = seg1[i].type)
                end
            end
                if seg1[i].iend >= seg2[j].iend
                    if j<N2; j+=1; continue
                    else stop = true end
                else
                    if i<N1; i+=1;continue
                    else stop = true end
            end
        end
        newSeg, tempSeg = addTempSeg(tempSeg, newSeg)
    end
    clearSegs1 = seg1[.!toDel1]
    clearSegs2 = seg2[.!toDel2]
    addSeg = []
    if !isempty(clearSegs1)
        addSeg = clearSegs1
    end
    if !isempty(clearSegs2)
        if !isempty(addSeg)
            append!(addSeg,clearSegs2)
        else
            addSeg = clearSegs2
        end
    end
    if !isempty(newSeg) && !isa(newSeg.ibeg,UndefInitializer)
        #newSeg был туплом, надо сделать страктэрреем
        if !isa(newSeg.ibeg,Array)
            newSegStr = StructArray(ibeg = [newSeg.ibeg], iend = [newSeg.iend], type = [newSeg.type])
        else
            newSegStr = StructArray(ibeg = newSeg.ibeg, iend = newSeg.iend, type = newSeg.type)
        end
        if !isempty(addSeg)
            append!(addSeg, newSegStr)
        else
            addSeg = newSegStr
        end
    end
    # @info addSeg
    outSeg = sortSeg(addSeg)
    return outSeg
end

"""
Схлопывание пересекающихся сегментов
ВНИМАНИЕ: если есть параметры, то они похерятся
"""
function mergeIntersections(seg::StructArray)
    N = length(seg)
    if N==1
        return seg
    end
    toDel = falses(N)
    idMerged = 0 #флаг, было ли на прошлом шаге слияние
    for i=1:N-1
        if seg[i].iend>=seg[i+1].ibeg #если следующий начинается раньше конца текущего
                if seg[i].type == seg[i+1].type #если одинаковые типы пересекаются
                    toDel[i] = true # удалим текущий
                    seg[i+1] = (ibeg= seg[i].ibeg, iend = max(seg[i].iend[1], seg[i+1].iend[1]), type = seg[i].type) #у следующего удлиним начало
                end
        end
    end
    seg[.!toDel]
end

"""
сравнение двух наборов сегментов на пересечения, аналог ф-ции в MATLAB
onlySeg1, вектор сегментов Seg1, которые не пересекаются с onlySeg2 и наоборот
сегменты не должны накладываться др на др
"""
function segCompare(seg1::StructArray, seg2::StructArray)
    N1 = length(seg1)
    N2 = length(seg2)
    i=j=1
    stop = false
    commonSeg1 = falses(N1)
    commonSeg2 = falses(N2);

    if N1 == 0 || N2 == 0
        only1ind = .!commonSeg1
        only2ind = .!commonSeg2
    else
        while !stop
            #Seg1 раньше Seg2
            if seg1[i].iend < seg2[j].ibeg # Seg1 кончается раньше, чем начинается Seg2
                if i<N1
                    i+=1
                    continue
                else
                    stop = true
                end
            end
            #Seg2 раньше Seg1
            if seg2[j].iend < seg1[i].ibeg  #Seg2 кончается раньше, чем начинается Seg1
                if j<N2
                    j+=1
                    continue
                else
                    stop = true
                end
            end
            # Seg2 пересекается с Seg1
            if sign(seg1[i].ibeg -seg2[j].iend) != sign(seg1[i].iend-seg2[j].ibeg)
                commonSeg1[i] = commonSeg2[j] = true
                if seg1[i].iend >= seg2[j].iend
                    if j<N2
                        j+=1
                        continue
                    else
                        stop = true
                    end
                else
                 	if i<N1
                        i=i+1
                        continue
                    else
                        stop = true
                    end
                end
            end
        end
        only1ind =  .!commonSeg1
        only2ind =  .!commonSeg2

    end
    return only1ind, only2ind
end
"""
удаляет любые сегменты из Seg1, которые пересекаются с Seg2 (Seg1 = Seg1 - Seg2)
"""
function segMinus_(seg1::StructArray, seg2::StructArray)
    only1, only2 = segCompare(seg1, seg2)
    return seg1[only1]
end
"""
версия 2 - удаляет любые сегменты из Seg1, которые пересекаются с Seg2 (Seg1 = Seg1 - Seg2)
"""
function segMinus(seg1::StructArray, seg2::StructArray)
    N1 = length(seg1)
    N2 = length(seg2)
    i=j=1
    stop = false
    takeID = trues(N1)

    if !(N1 == 0 || N2 == 0)
        while !stop
            #Seg1 раньше Seg2
            if seg1[i].ibeg<= seg2[j].iend # начало 1 меньше конца 2
                if seg1[i].iend >= seg2[j].ibeg
                    takeID[i] = false
                end
                if i<N1
                    i+=1
                    continue
                else
                    stop = true
                end
            else
                if j<N2
                    j+=1
                    continue
                else
                    stop = true
                end
            end
        end
    end
    return seg1[takeID]
end


text:
using Jecco.AdS5_3_1

grid = SpecCartGrid3D(
    x_min            = -10.0,
    x_max            =  10.0,
    x_nodes          =  128,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  6,
    u_outer_min      =  0.1,
    u_outer_max      =  1.01,
    u_outer_domains  =  1,
    u_outer_nodes    =  48,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

id = BlackBranePert(
    energy_dens = 1.0,
    a4_ampx  = 5.e-2,
    a4_kx    = 1,
    AH_pos   = 1.0,
    xmax     = grid.x_max,
    xmin     = grid.x_min,
    ymin     = grid.y_min,
    ymax     = grid.y_max,
)

evoleq = AffineNull(
    phi0           = 0.0,
    potential      = ZeroPotential(),
    gaugecondition = ConstantAH(u_AH = 1.0),
)

io = InOut(
    out_boundary_every  = 200,
    out_bulk_every      = 1000,
    out_gauge_every     = 200,
)

integration = Integration(
    dt              = 0.001,
    tmax            = 2.0,
    ODE_method      = AdS5_3_1.AB4(),
    filter_poststep = true,
)

run_model(grid, id, evoleq, integration, io)


text:
module Mongoc

using Dates

# load libmongoc
const libmongocpath = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(libmongocpath)
    error("Mongoc.jl is not installed properly, run Pkg.build(\"Mongoc\") and restart Julia.")
end
include(libmongocpath)

# utility functions for date conversion
const ISODATE_OFFSET = Dates.UNIXEPOCH + 31536000000 # offsets an additional year from UNIXEPOCH. Number of millis in a year: 24 * 60 * 60 * 1000 * 365
isodate2datetime(millis::Int64) = Dates.epochms2datetime(millis + ISODATE_OFFSET)
datetime2isodate(dt::DateTime) = Dates.datetime2epochms(dt) - ISODATE_OFFSET

include("bson.jl")
include("types.jl")
include("c_api.jl")
include("api.jl")
include("session.jl")

function __init__()
    check_deps()
    mongoc_init()
end

end # module Mongoc


text:
using ITensors, ITensorNetworkAD, AutoHOOT, Zygote, OptimKit
using ITensorNetworkAD.ITensorNetworks:
  PEPS,
  inner_network,
  inner_networks,
  Models,
  flatten,
  insert_projectors,
  split_network,
  GeneralMPSTensor,
  TreeTensor,
  broadcast_inner,
  tree
using ITensorNetworkAD.Optimizations:
  gradient_descent, backtracking_linesearch, loss_grad_wrap
using ITensorNetworkAD.ITensorAutoHOOT: batch_tensor_contraction

@testset "test MPSTensor" begin
  Nx, Ny = 3, 3
  num_sweeps = 20
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  H = Models.mpo(Models.Model("tfim"), sites; h=1.0)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  params = Dict(:maxdim => 1000, :cutoff => 1e-15)
  loss_w_grad_mps = loss_grad_wrap(peps, H_line, GeneralMPSTensor; params...)
  loss_w_grad = loss_grad_wrap(peps, H_line)
  loss_mps, grad_mps = loss_w_grad_mps(peps)
  loss, grad = loss_w_grad(peps)
  g_mps_nrm = broadcast_inner(grad_mps, grad_mps)
  g_nrm = broadcast_inner(grad, grad)
  @test isapprox(loss, loss_mps)
  @test isapprox(g_nrm, g_mps_nrm)
end

@testset "test TreeTensor" begin
  ITensors.set_warn_order(100)
  Nx, Ny = 3, 3
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  params = Dict(:maxdim => 1000, :cutoff => 1e-15)
  loss_w_grad_tree = loss_grad_wrap(peps, H_line, TreeTensor, tree; params...)
  loss_w_grad = loss_grad_wrap(peps, H_line)
  loss_tree, grad_tree = loss_w_grad_tree(peps)
  loss, grad = loss_w_grad(peps)
  g_tree_nrm = broadcast_inner(grad_tree, grad_tree)
  g_nrm = broadcast_inner(grad, grad)
  @test isapprox(loss, loss_tree)
  @test isapprox(g_nrm, g_tree_nrm)
end

@testset "test the loss of optimization" begin
  Nx, Ny = 2, 3
  num_sweeps = 20
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=10)
  randn!(peps)
  psi0 = randomMPS(vec(sites); linkdims=10)
  sweeps = Sweeps(10)
  maxdim!(sweeps, 10)

  H = Models.mpo(Models.Model("tfim"), sites; h=1.0)
  energy, _ = dmrg(H, psi0, sweeps)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  losses_gd = gradient_descent(peps, H_line; stepsize=0.005, num_sweeps=num_sweeps)
  losses_ls = optimize(peps, H_line; num_sweeps=num_sweeps, method="GD")
  losses_lbfgs = optimize(peps, H_line; num_sweeps=num_sweeps, method="LBFGS")
  losses_cg = optimize(peps, H_line; num_sweeps=num_sweeps, method="CG")
  for i in 3:(length(losses_gd) - 1)
    @test losses_gd[i] >= losses_gd[i + 1]
    @test losses_ls[i] >= losses_ls[i + 1]
    @test losses_lbfgs[i] >= losses_lbfgs[i + 1]
    @test losses_cg[i] >= losses_cg[i + 1]
  end
  @test abs(energy - losses_gd[end]) < 0.5
  @test abs(energy - losses_ls[end]) < 0.5
  @test abs(energy - losses_lbfgs[end]) < 0.5
  @test abs(energy - losses_cg[end]) < 0.5
end

@testset "test the loss of optimization" begin
  Nx, Ny = 3, 3
  num_sweeps = 20
  cutoff = 1e-15
  maxdim = 100
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)

  psi0 = randomMPS(vec(sites); linkdims=3)
  sweeps = Sweeps(10)
  maxdim!(sweeps, 10)

  H = Models.mpo(Models.Model("tfim"), sites; h=1.0)
  energy, _ = dmrg(H, psi0, sweeps)
  print(energy)

  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  losses_gdls = gradient_descent(
    peps,
    H_line,
    insert_projectors,
    backtracking_linesearch;
    beta=0.5,
    num_sweeps=num_sweeps,
    cutoff=cutoff,
    maxdim=maxdim,
  )
  for i in 3:(length(losses_gdls) - 1)
    @test losses_gdls[i] >= losses_gdls[i + 1]
  end
  @test abs(energy - losses_gdls[end]) < 0.5
end

@testset "test the equivalence of local and line hamiltonian" begin
  Nx, Ny = 2, 3
  num_sweeps = 3
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=10)
  randn!(peps)
  H_local = Models.localham(Models.Model("tfim"), sites; h=1.0)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  losses_line = gradient_descent(peps, H_line; stepsize=0.005, num_sweeps=num_sweeps)
  losses_local = gradient_descent(peps, H_local; stepsize=0.005, num_sweeps=num_sweeps)
  @test isapprox(losses_local, losses_line)
end

@testset "test monotonic loss decrease of optimization with inserting projectors" begin
  Nx, Ny = 3, 3
  num_sweeps = 20
  cutoff = 1e-15
  maxdim = 100
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  losses_gd = gradient_descent(
    peps,
    H_line,
    insert_projectors;
    stepsize=0.005,
    num_sweeps=num_sweeps,
    cutoff=cutoff,
    maxdim=maxdim,
  )
  losses_ls = optimize(
    peps,
    H_line,
    insert_projectors;
    num_sweeps=num_sweeps,
    method="GD",
    cutoff=cutoff,
    maxdim=maxdim,
  )
  losses_lbfgs = optimize(
    peps,
    H_line,
    insert_projectors;
    num_sweeps=num_sweeps,
    method="LBFGS",
    cutoff=cutoff,
    maxdim=maxdim,
  )
  losses_cg = optimize(
    peps,
    H_line,
    insert_projectors;
    num_sweeps=num_sweeps,
    method="CG",
    cutoff=cutoff,
    maxdim=maxdim,
  )
  for i in 3:(length(losses_gd) - 1)
    @test losses_gd[i] >= losses_gd[i + 1]
    @test losses_ls[i] >= losses_ls[i + 1]
    @test losses_lbfgs[i] >= losses_lbfgs[i + 1]
    @test losses_cg[i] >= losses_cg[i + 1]
  end
end

@testset "test inner product gradient" begin
  Nx = 2
  Ny = 2
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  function loss(peps::PEPS)
    peps_prime = prime(linkinds, peps)
    peps_prime_ham = prime(peps)
    network_list = [inner_network(peps, peps_prime)]
    variables = flatten([peps, peps_prime])
    inners = batch_tensor_contraction(network_list, variables...)
    return sum(inners)[]
  end
  g = gradient(loss, peps)
  inner = inner_network(peps, prime(linkinds, peps))
  g_true_first_site = contract(inner[2:length(inner)])
  g_true_first_site = 2 * g_true_first_site
  @test isapprox(g[1].data[1, 1], g_true_first_site)
end

@testset "test split network" begin
  ITensors.set_warn_order(40)
  Ny, Nx = 3, 4
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  H_line = [H for H in H_line if H.coord[2] isa Colon]
  tn_split_row, tn_split_column, projectors_row, projectors_column = insert_projectors(
    peps, 1e-15, 1000
  )
  for (i, H) in enumerate(H_line)
    function loss(peps::PEPS)
      peps_bra = addtags(linkinds, peps, "bra")
      peps_ket = addtags(linkinds, peps, "ket")
      sites = commoninds(peps_bra, peps_ket)
      peps_bra_split = split_network(peps_bra)
      peps_ket_split = split_network(peps_ket)
      peps_ket_split_ham = prime(sites, peps_ket_split)
      network_list = inner_networks(
        peps_bra_split, peps_ket_split, peps_ket_split_ham, projectors_row[i], [H]
      )
      variables = flatten([peps_bra_split])
      inners = batch_tensor_contraction(network_list, variables...)
      return sum(inners)[]
    end
    g = gradient(loss, peps)
    inner = inner_networks(peps, prime(linkinds, peps), prime(peps), [H])[1]
    g_true_first_site = contract(vcat(inner[2:length(inner)]))
    @test isapprox(g[1].data[1, 1], g_true_first_site)
  end
end

@testset "test inner product gradient with tagging" begin
  Nx, Ny = 3, 3
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  function loss(peps::PEPS)
    peps_bra = addtags(linkinds, peps, "bra")
    peps_ket = addtags(linkinds, peps, "ket")
    sites = commoninds(peps_bra, peps_ket)
    projectors = [ITensor(1.0)]
    network_list = [inner_network(peps_bra, peps_ket, projectors)]
    variables = flatten([peps_bra, peps_ket])
    inners = batch_tensor_contraction(network_list, variables...)
    return sum(inners)[]
  end
  g = gradient(loss, peps)
  inner = inner_network(peps, prime(linkinds, peps))
  g_true_first_site = contract(inner[2:length(inner)])
  g_true_first_site = 2 * g_true_first_site
  @test isapprox(g[1].data[1, 1], g_true_first_site)
  @test isapprox(loss(peps), contract(inner)[])
end


text:
module MLMC
export mlmc_test_short,mlmc,MLMCResult,MLMCResult_test,statistics, plot_test
using DocOpt
    using Gadfly # package for plotting
    using DataFrames

type MLMCResult
    eps::Float64
    estimate::Float64
    tictoc::Float64
    nums::Array{Float64,1} # would like this to be int
    suml::Array{Float64,2}
    L::Int64
    alpha::Float64
    beta::Float64
end

type MLMCResult_test
    alpha::Float64
    beta::Float64
    gamma::Float64
    kurtosis::Array{Float64,1}
    variance_l::Array{Float64,1}
    mean_l::Array{Float64,1}
    cost_l::Array{Float64,1}
    L::Int64
end

function mlmc(id,eps,mlmc_l, alpha_0=0,beta_0=0,gamma=0,Lmin=4,Lmax=32,N0=1000)

    doc= """ Follows Giles 2008 """
    srand(id);
    if (Lmin<2)
        error("error: needs Lmin >= 2");
    end

    if (Lmax<Lmin)
        error("error: needs Lmax >= Lmin");
    end

    if (N0<=0 || eps<=0 || gamma <= 0)
        error("error: needs N>0, eps>0, gamma>0 \n");
    end

    alpha = alpha_0;
    beta  = beta_0;
    M    = 2.0;
    dNl  = zeros(Lmax+1);
    Nl   = zeros(Lmax+1); # should be int
    dNl[1:(Lmin+1)] = N0;
    suml = zeros(3, Lmax+1)
    L = Lmin;
    # assigns the proportion of the error between Variance and Bias
    theta = 0.25;
    tic()
    while sum(dNl) > 0
        #update sample sums
        for l=0:L
            if dNl[l+1] > 0
#                 @show (ceil(Int64,dNl[l+1]),l);
                start_time=time()
                samples     =  mlmc_l(ceil(Int64,dNl[l+1]),l);
                suml[3,l+1]+= time()-start_time
                Nl[l+1]     = Nl[l+1] +length(samples);
                suml[1,l+1] = suml[1,l+1] + sum(samples);
                suml[2,l+1] = suml[2,l+1] + sum(samples.^2);
            end
        end

        #compute absolute average and variance

        mean_l = abs(vec(suml[1,1:(L+1)])./Nl[1:(L+1)]);
        Vars_l = max(0, vec(suml[2,1:(L+1)])./Nl[1:(L+1)] - mean_l.^2);

        #fix to cope with possible zero values for ml and Vl
        #(can happen in some applications when there are few samples)

        for l = 3:(L+1)
            mean_l[l] = max(mean_l[l], 0.5*mean_l[l-1]/2.0^alpha);
            Vars_l[l] = max(Vars_l[l], 0.5*Vars_l[l-1]/2.0^beta);
        end

        if alpha_0==0
            pa = linreg(vec(float([2:(L+1);])) ,vec(log2(mean_l[2:(L+1)])));
            alpha = -pa[2];
        end
        if beta_0==0
            pb = linreg(vec(float([2:(L+1);])) ,vec(log2(Vars_l[2:(L+1)])));
            beta = -pb[2];
        end

        # Cost of the level. Should be estimated with regression
        Cl  = 2.^(gamma*(0:L));

        #set optimal number of additional samples
        Ns   = zeros(L+1);

        Ns  = ceil(Integer, sqrt(Vars_l./Cl) * sum(sqrt(Vars_l.*Cl))/ ((1-theta)*eps^2) );
        dNl = max(zeros(L+1), Ns-Nl[1:(L+1)]);

        #if (almost) converged, estimate remaining error and decide
        #whether a new level is required

        if maximum(dNl./Nl[1:(L+1)]) < 0.01
            rem = mean_l[L+1] / (2^alpha - 1);
            if rem > sqrt(theta)*eps
                if (L==Lmax)
                    println("*** failed to achieve weak convergence *** ")
                else
                    L       = L+1;
                    Vars_l = push!(Vars_l,Vars_l[L] / 2^beta);
                    Nl[L+1] = 0;
                    Cl  = 2.^(gamma*(0:L));
                    Ns   = zeros(L+1);
                    Ns  = ceil( sqrt(Vars_l./Cl) * sum(sqrt(Vars_l.*Cl)) / ((1-theta)*eps^2) );
                    dNl = max(0, Ns-Nl[1:(L+1)]);
                end
            end
        end
  end
# finally, evaluate multilevel estimator
        tictoc=toq()
    P = sum(vec(suml[1,1:(L+1)])./Nl[1:(L+1)]);

    return MLMCResult(eps,P,tictoc, Nl,suml, L, alpha, beta)
end

function statistics(res::MLMCResult)
    levels=res.L
    nums=res.nums[1:levels]
    means=vec(res.suml[1,1:levels])./nums
    vars= vec(res.suml[2,1:levels])./nums -means.^2
    (means, vars)
end

function mlmc_test_short(mlmc_l,N,L)

    del1 = zeros(L+1);
    var1 = zeros(L+1);
    kur1 = zeros(L+1);
    cost = zeros(L+1);

    for l = 0:L
        tic();
        samples     = mlmc_l(N,l);
        s4   = sum(samples.^4)/N;
        s3   = sum(samples.^3)/N;
        s2   = sum(samples.^2)/N;
        s1   = sum(samples   )/N;
        cost[l+1] = toq();

        if (l==0)
            kurt = 0.0;
        else
            kurt = ( s4 - 4*s3*s1 + 6*s2*s1^2 - 3*s1^4 ) / (s2-s1^2)^2;
        end

        del1[l+1] = s1;
        var1[l+1] = s2-s1^2;
        kur1[l+1] = kurt;
    end

    if ( kur1[L+1] > 100.0 )
        println("\n WARNING: kurtosis on finest level = $(kur1[L+1])\n",);
    end

    L1 = ceil(Integer, 0.4*L)+1;
    L2 = L+1;
    X = float(collect(L1:L2));
    pa = linreg(X,log2(abs(del1[L1:L2])));  alpha = -pa[2];
    pb = linreg(X,log2(abs(var1[L1:L2])));  beta  = -pb[2];
    gamma = log2(cost[L+1]/cost[L]);

#    println("\n******************************************************\n");
#    println("*** Linear regression estimates of MLMC parameters ***\n");
#    println("******************************************************\n");
#    println("\n alpha = $(alpha)  (exponent for MLMC weak convergence)\n");
#    println(" beta  = $(beta)  (exponent for MLMC variance) \n");
#    println(" gamma = $(gamma)  (exponent for MLMC cost) \n");
    return MLMCResult_test(alpha, beta, gamma, kur1, var1, del1, cost, L)

end

function plot_test(res::MLMCResult_test)
    df1 = DataFrame(x=0:res.L, y = vec(res.variance_l) , label="Variance")
    df2 = DataFrame(x=0:res.L, y = vec(10.^(-(0:res.L))), label="First order decay")
    df3 = DataFrame(x=0:res.L, y = vec(abs(res.mean_l)), label="Mean")
    dfp1 = vcat(df1, df2, df3);
    #dfp2 = vcat(df3, df2);

    p1 = plot(dfp1, x="x", y="y", color="label", Scale.y_log,  Guide.xlabel("level"), Guide.ylabel("\log(Variance)"), Guide.title("Variance decay"), Geom.point, Geom.line)
#    p2 = plot(dfp2, x="x", y="y", color="label", Scale.y_log,  Guide.xlabel("level"), Guide.ylabel("\log(Mean)"),     Guide.title("Mean decay"), Geom.point, Geom.line)
#Scale.discrete_color_manual(colorant"blue",colorant"red", colorant"green", colorant"orange",colorant"black"))
    draw(PDF("mlmc_test.pdf", 12inch, 6inch),p1)

end
end

max_stars_repo_name: vollmersj/MLMC

max_issues_repo_name: vollmersj/MLMC

max_forks_repo_name: vollmersj/MLMC