text:
RuntimeGeneratedFunctions.init(@__MODULE__)
"""
Algorithm for solving Physics-Informed Neural Networks problems.

Arguments:
* `chain` is a Flux.jl chain with a d-dimensional input and a 1-dimensional output,
* `init_params` is the initial parameter of the neural network,
* `phi` is a trial solution,
* `autodiff` is a boolean variable that determines whether to use automatic, differentiation (not supported while) or numerical,
* `derivative` is method that calculates the derivative,
* `strategy` determines which training strategy will be used.
"""

struct PhysicsInformedNN{C,P,PH,DER,T,K}
  chain::C
  initθ::P
  phi::PH
  autodiff::Bool
  derivative::DER
  strategy::T
  kwargs::K
end

function PhysicsInformedNN(chain,
                           init_params = nothing;
                           _phi = nothing,
                           autodiff=false,
                           _derivative = nothing,
                           strategy = GridTraining(),
                           kwargs...)
    if init_params === nothing
        initθ = DiffEqFlux.initial_params(chain)
    else
        initθ = init_params
    end

    if _phi == nothing
        phi = get_phi(chain)
    else
        phi = _phi
    end

    if _derivative == nothing
        derivative = get_derivative(autodiff)
    else
        derivative = _derivative
    end

    PhysicsInformedNN(chain,initθ,phi,autodiff,derivative, strategy, kwargs)
end

abstract type TrainingStrategies  end

"""
* `dx` is the discretization of the grid
"""
struct GridTraining <: TrainingStrategies
    dx
end
function GridTraining(;dx= 0.1)
    GridTraining(dx)
end

"""
* `number_of_points` is number of points in random select training set
"""
struct StochasticTraining <:TrainingStrategies
    number_of_points:: Int64
end
function StochasticTraining(;number_of_points=100)
    StochasticTraining(number_of_points)
end

    # sampling_method::
struct QuasiRandomTraining <:TrainingStrategies
    sampling_method::QuasiMonteCarlo.SamplingAlgorithm
    number_of_points:: Int64
    number_of_minibatch:: Int64
end
function QuasiRandomTraining(;sampling_method = UniformSample(), number_of_points=100,number_of_minibatch=10)
    QuasiRandomTraining(sampling_method,number_of_points,number_of_minibatch)
end
"""
* `algorithm`: quadrature algorithm,
* `reltol`: relative tolerance,
* `abstol` absolute tolerance,
* `maxiters`: the maximum number of iterations in quadrature algorithm,
* `batch`: the preferred number of points to batch.

For more information look: Quadrature.jl https://github.com/SciML/Quadrature.jl
"""
struct QuadratureTraining <: TrainingStrategies
    algorithm::DiffEqBase.AbstractQuadratureAlgorithm
    reltol::Float64
    abstol::Float64
    maxiters::Int64
    batch::Int64
end

function QuadratureTraining(;algorithm=HCubatureJL(),reltol= 1e-8,abstol= 1e-8,maxiters=1e3,batch=0)
    QuadratureTraining(algorithm,reltol,abstol,maxiters,batch)
end

"""
Create dictionary: variable => unique number for variable

# Example 1

Dict{Symbol,Int64} with 3 entries:
  :y => 2
  :t => 3
  :x => 1

# Example 2

 Dict{Symbol,Int64} with 2 entries:
  :u1 => 1
  :u2 => 2
"""
get_dict_vars(vars) = Dict( [Symbol(v) .=> i for (i,v) in enumerate(vars)])

# Wrapper for _transform_derivative
function transform_derivative(ex,dict_indvars,dict_depvars)
    if ex isa Expr
        ex.args = _transform_derivative(ex.args,dict_indvars,dict_depvars)
    end
    return ex
end

function get_ε(dim, der_num)
    epsilon = cbrt(eps(Float32))
    ε = zeros(Float32, dim)
    ε[der_num] = epsilon
    ε
end

θ = gensym("θ")

"""
Transform the derivative expression to inner representation

# Examples

1. First compute the derivative of function 'u(x,y)' with respect to x.

Take expressions in the form: `derivative(u(x,y), x)` to `derivative(phi, u, [x, y], εs, order, θ)`,
where
 phi - trial solution
 u_d - derived function
 x,y - coordinates of point
 εs - epsilon mask
 order - order of derivative
 θ - weight in neural network
"""
function _transform_derivative(_args,dict_indvars,dict_depvars)
    dim = length(dict_indvars)
    εs = [get_ε(dim,d) for d in 1:dim]
    for (i,e) in enumerate(_args)
        if !(e isa Expr)
            if e in keys(dict_depvars)
                push!(_args,:($θ), :phi)
            elseif e == :derivative
                derivative_variables = Symbol[]
                order = 0
                while (_args[1] == :derivative)
                    order += 1
                    push!(derivative_variables, _args[end])
                    _args = _args[2].args
                end
                depvar = _args[1]
                indvars = _args[2:end]
                undv = [dict_indvars[d_p] for d_p  in derivative_variables]
                εs_dnv = [εs[d] for d in undv]
                _args = [:derivative, :phi, Symbol(:($depvar),:_d), :([$(indvars...)]), εs_dnv, order, :($θ)]
                break
            end
        else
            _args[i].args = _transform_derivative(_args[i].args,dict_indvars,dict_depvars)
        end
    end
    return _args
end

"""
Parse ModelingToolkit equation form to the inner representation.

Example:

1)  1-D ODE: Dt(u(t)) ~ t +1

    Take expressions in the form: 'Equation(derivative(u(t), t), t + 1)' to 'derivative(phi, u_d, [t], [[ε]], 1, θ) - (t + 1)'

2)  2-D PDE: Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

    Take expressions in the form:
     Equation(derivative(derivative(u(x, y), x), x) + derivative(derivative(u(x, y), y), y), -(sin(πx)) * sin(πy))
    to
     (derivative(phi,u_d, [x, y], [[ε,0],[ε,0]], 2, θ) + derivative(phi, u_d, [x, y], [[0,ε],[0,ε]], 2, θ)) - -(sin(πx)) * sin(πy)

3)  System of PDEs: [Dx(u1(x,y)) + 4*Dy(u2(x,y)) ~ 0,
                    Dx(u2(x,y)) + 9*Dy(u1(x,y)) ~ 0]

    Take expressions in the form:
    2-element Array{Equation,1}:
        Equation(derivative(u1(x, y), x) + 4 * derivative(u2(x, y), y), ModelingToolkit.Constant(0))
        Equation(derivative(u2(x, y), x) + 9 * derivative(u1(x, y), y), ModelingToolkit.Constant(0))
    to
      [(derivative(phi, u1_d, [x, y], [[ε,0]], 1, θ) + 4 * derivative(phi, u2_d, [x, y], [[0,ε]], 1, θ)) - 0,
       (derivative(phi, u2_d, [x, y], [[ε,0]], 1, θ) + 9 * derivative(phi, u1_d, [x, y], [[0,ε]], 1, θ)) - 0]
"""
function parse_equation(eq,dict_indvars,dict_depvars)
    left_expr = transform_derivative(toexpr(eq.lhs),
                                     dict_indvars,dict_depvars)
    right_expr = transform_derivative(toexpr(eq.rhs),
                                     dict_indvars,dict_depvars)
    loss_func = :($left_expr - $right_expr)
end

"""
Build a loss function for a PDE or a boundary condition

# Examples: System of PDEs:

Take expressions in the form:

[Dx(u1(x,y)) + 4*Dy(u2(x,y)) ~ 0,
 Dx(u2(x,y)) + 9*Dy(u1(x,y)) ~ 0]

to

:((phi, cord, θ, derivative)->begin
          #= none:35 =#
          #= none:35 =#
          begin
              begin
                  u1 = ((x, y, θ)->begin
                              #= none:16 =#
                              (phi([x, y], θ))[1]
                          end)
                  u2 = ((x, y, θ)->begin
                              #= none:16 =#
                              (phi([x, y], θ))[2]
                          end)
                  u1_d = ((cord, θ)->begin
                              #= none:20 =#
                              (phi(cord, θ))[1]
                          end)
                  u2_d = ((cord, θ)->begin
                              #= none:20 =#
                              (phi(cord, θ))[2]
                          end)
              end
              let (x, y) = (cord[1], cord[2])
                  [(derivative(phi, u1_d, [x, y], [[ε, 0]], 1, θ) + 4 * derivative(phi, u2_d, [x, y], [[0, ε]], 1, θ)) - 0,
                   (derivative(phi, u2_d, [x, y], [[ε, 0]], 1, θ) + 9 * derivative(phi, u1_d, [x, y], [[0, ε]], 1, θ)) - 0]
              end
          end
      end)
"""
function build_symbolic_loss_function(eqs,_indvars,_depvars, phi, derivative;bc_indvars=nothing)
    # dictionaries: variable -> unique number
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_symbolic_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars, phi, derivative,bc_indvars = bc_indvars)
end
function build_symbolic_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars, phi, derivative;
                             bc_indvars = indvars)
    if !(eqs isa Array)
        eqs = [eqs]
    end
    loss_functions= Expr[]
    for eq in eqs
        push!(loss_functions,parse_equation(eq,dict_indvars,dict_depvars))
    end

    vars = :(cord, $θ, phi, derivative,u_arr,u_d_arr)
    ex = Expr(:block)

    depvars_d = Symbol[]
    for v in depvars
        push!(depvars_d,:($(Symbol(:($v),:_d))))
    end

    expr_u_arr = Expr[]
    expr_u_d_arr = Expr[]
    for i in eachindex(depvars)
        push!(expr_u_arr, :(u_arr[$i]) )
        push!(expr_u_d_arr, :(u_d_arr[$i]) )
    end

    vars_u = Expr(:(=), build_expr(:tuple, depvars), build_expr(:tuple, expr_u_arr))
    push!(ex.args,  vars_u)

    vars_u_d = Expr(:(=), build_expr(:tuple, depvars_d), build_expr(:tuple, expr_u_d_arr))
    push!(ex.args,  vars_u_d)

    indvars_ex = [:($:cord[$i]) for (i, u) ∈ enumerate(bc_indvars)]

    left_arg_pairs, right_arg_pairs = bc_indvars,indvars_ex

    vars_eq = Expr(:(=), build_expr(:tuple, left_arg_pairs), build_expr(:tuple, right_arg_pairs))
    let_ex = Expr(:let, vars_eq, build_expr(:vect, loss_functions))
    push!(ex.args,  let_ex)

    expr_loss_function = :(($vars) -> begin $ex end)
end

function build_loss_function(eqs,_indvars,_depvars, phi, derivative;bc_indvars=nothing)

    # dictionaries: variable -> unique number
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars, phi, derivative,bc_indvars = bc_indvars)
end
function build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars, phi, derivative;
                             bc_indvars = indvars)
    depvars_d = Symbol[]
    for v in depvars
        push!(depvars_d,:($(Symbol(:($v),:_d))))
    end

    us = Expr[]
    for v in depvars
        var_num = dict_depvars[v]
        push!(us,:(($(indvars...), $θ, phi) -> phi([$(indvars...)],$θ)[$var_num]))
    end
    u_ds = Expr[]
    for (v,v_d) in zip(depvars, depvars_d)
        var_num = dict_depvars[v]
        push!(u_ds,:((cord, $θ, phi) -> phi(cord,$θ)[$var_num]))
    end

    expr_loss_function = build_symbolic_loss_function(eqs,indvars,depvars,
                                                      dict_indvars,dict_depvars,
                                                      phi, derivative;
                                                      bc_indvars = bc_indvars)
    u_arr = [@RuntimeGeneratedFunction(u) for u in us]
    u_d_arr = [@RuntimeGeneratedFunction(u_d) for u_d in u_ds]
    _loss_function = @RuntimeGeneratedFunction(expr_loss_function)
    loss_function = (cord, θ) -> _loss_function(cord, θ, phi, derivative, u_arr, u_d_arr)
    return loss_function
end

# Get arguments from boundary condition functions
function get_bc_argument(bcs,dict_indvars,dict_depvars)
    bc_args = []
    for _bc in bcs
        _bc isa Array && error("boundary conditions must be represented as a one-dimensional array")
        left_expr = transform_derivative(toexpr(_bc.lhs),
                                         dict_indvars,dict_depvars)
        bc_arg = if (left_expr.args[1] == :derivative)
            left_expr.args[4].args
        else
            left_expr.args[2:end-2]
        end
        push!(bc_args, bc_arg)
    end

    return bc_args
end

function get_bc_varibles(bcs,_indvars::Array,_depvars::Array)
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return get_bc_varibles(bcs,dict_indvars,dict_depvars)
end

function get_bc_varibles(bcs,dict_indvars,dict_depvars)
    bc_args = get_bc_argument(bcs,dict_indvars,dict_depvars)
    return map(barg -> filter(x -> x isa Symbol, barg), bc_args)
end

function generate_training_sets(domains,dx,bcs,_indvars::Array,_depvars::Array)
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return generate_training_sets(domains,dx,bcs,dict_indvars,dict_depvars)
end
# Generate training set in the domain and on the boundary
function generate_training_sets(domains,dx,bcs,dict_indvars::Dict,dict_depvars::Dict)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end

    bound_args = get_bc_argument(bcs,dict_indvars,dict_depvars)
    bound_vars = get_bc_varibles(bcs,dict_indvars,dict_depvars)

    spans = [d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)]
    dict_var_span = Dict([Symbol(d.variables) => d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)])

    dif = [Float64[] for i=1:size(domains)[1]]
    for _args in bound_args
        for (i,x) in enumerate(_args)
            if x isa Number
                push!(dif[i],x)
            end
        end
    end
    collect_spans = collect.(spans)
    bc_data = map(zip(dif,collect_spans)) do (d,c)
        setdiff(c, d)
    end

    train_set = vec(map(points -> collect(points), Iterators.product(spans...)))
    pde_train_set = vec(map(points -> collect(points), Iterators.product(bc_data...)))

    bcs_train_set = map(bound_vars) do bt
        span = map(b -> dict_var_span[b], bt)
        _set = vec(map(points -> collect(points), Iterators.product(span...)))
    end

    [pde_train_set,bcs_train_set,train_set]
end

function get_phi(chain)
    # The phi trial solution
    if chain isa FastChain
        phi = (x,θ) -> chain(adapt(typeof(θ),x),θ)
    else
        _,re  = Flux.destructure(chain)
        phi = (x,θ) -> re(θ)(adapt(typeof(θ),x))
    end
    phi
end

# the method to calculate the derivative
function get_derivative(autodiff)
    epsilon = 2*cbrt(eps(Float32))
    if autodiff # automatic differentiation (not implemented yet)
        error("automatic differentiation is not implemented yet)")
    else # numerical differentiation
        derivative = (phi,u,x,εs,order,θ) ->
        begin
            ε = εs[order]
            if order > 1
                return (derivative(phi,u,x+ε,εs,order-1,θ)
                      - derivative(phi,u,x-ε,εs,order-1,θ))/epsilon
            else
                return (u(x+ε,θ,phi) - u(x-ε,θ,phi))/epsilon
            end
        end
    end
    derivative
end

function get_loss_function(loss_functions, train_sets, strategy::GridTraining)
    # norm coefficient for loss function
    τ_ = loss_functions isa Array ? sum(length(train_set) for train_set in train_sets) : length(train_sets)
    τ = 1.0f0 / τ_

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        train_sets = [train_sets]
    end
    f = (loss,train_set,θ) -> sum(abs2,[inner_loss(loss,x,θ) for x in train_set])
    loss = (θ) ->  τ * sum(f(loss_function,train_set,θ) for (loss_function,train_set) in zip(loss_functions,train_sets))
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::StochasticTraining)
    number_of_points = strategy.number_of_points
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end
    τ = 1.0f0 / number_of_points

    loss = (θ) -> begin
        total = 0.
        for (lb, ub,l) in zip(lbs, ubs, loss_functions)
            len = length(lb)
            for i in 1:number_of_points
                r_point = lb .+ ub .* rand(len)
                total += inner_loss(l,r_point,θ)^2
            end
        end
        return τ * total
    end
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::QuasiRandomTraining)
    sampling_method = strategy.sampling_method
    number_of_points = strategy.number_of_points
    number_of_minibatch = strategy.number_of_minibatch
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end
    τ = number_of_points
    ss =[]
    for (lb, ub) in zip(lbs, ubs)
        s = QuasiMonteCarlo.generate_design_matrices(number_of_points,lb,ub,sampling_method,number_of_minibatch)
        # s =QuasiMonteCarlo.sample(number_of_points,lb,ub,UniformSample())
        push!(ss,s)
    end
    loss = (θ) -> begin
        total = 0.
        for (lb, ub,s_,l) in zip(lbs,ubs,ss,loss_functions)
            s =  s_[rand(1:number_of_minibatch)]
            ste_ = size(lb)[1]
            k = size(lb)[1]-1
            for i in 1:ste_:ste_*number_of_points
                r_point = s[i:i+k]
                total += inner_loss(l,r_point,θ)^2
            end
        end
        return (1.0f0/τ) * total
    end
    return loss
end

function get_bounds(domains,bcs,_indvars::Array,_depvars::Array)
    depvars = [nameof(value(d)) for d in _depvars]
    indvars = [nameof(value(i)) for i in _indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return get_bounds(domains,bcs,dict_indvars,dict_depvars)
end

function get_bounds(domains,bcs,dict_indvars,dict_depvars)
    bound_vars = get_bc_varibles(bcs,dict_indvars,dict_depvars)

    pde_lower_bounds = [d.domain.lower for d in domains]
    pde_upper_bounds = [d.domain.upper for d in domains]
    pde_bounds= [pde_lower_bounds,pde_upper_bounds]

    dict_lower_bound = Dict([Symbol(d.variables) => d.domain.lower for d in domains])
    dict_upper_bound = Dict([Symbol(d.variables) => d.domain.upper for d in domains])

    bcs_lower_bounds = map(bound_vars) do bt
        map(b -> dict_lower_bound[b], bt)
    end
    bcs_upper_bounds = map(bound_vars) do bt
        map(b -> dict_upper_bound[b], bt)
    end
    bcs_bounds= [bcs_lower_bounds,bcs_upper_bounds]

    [pde_bounds, bcs_bounds]
end

function get_loss_function(loss_functions, bounds, strategy::QuadratureTraining)
    lbs,ubs = bounds

    function inner_loss(loss_function,x,θ)
        sum(loss_function(x, θ))
    end

    if !(loss_functions isa Array)
        loss_functions = [loss_functions]
        lbs = [lbs]
        ubs = [ubs]
    end

    τ = 1.0f0 / ((10)^length(ubs[1])*length(ubs))

    f = (lb,ub,loss_,θ) -> begin
        _loss = (x,θ) -> sum(abs2,inner_loss(loss_, x, θ))
        prob = QuadratureProblem(_loss,lb,ub,θ;batch = strategy.batch)
        solve(prob,
              strategy.algorithm,
              reltol = strategy.reltol,
              abstol = strategy.abstol,
              maxiters = strategy.maxiters)[1]
    end
    loss = (θ) -> τ*sum(f(lb,ub,loss_,θ) for (lb,ub,loss_) in zip(lbs,ubs,loss_functions))
    return loss
end
function symbolic_discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)

    depvars = [nameof(value(d)) for d in pde_system.depvars]
    indvars = [nameof(value(i)) for i in pde_system.indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)

    chain = discretization.chain
    initθ = discretization.initθ
    phi = discretization.phi
    autodiff = discretization.autodiff
    derivative = discretization.derivative
    autodiff == true && error("Automatic differentiation is not support yet")
    strategy = discretization.strategy

    symbolic_pde_loss_function = build_symbolic_loss_function(eqs,indvars,depvars,
                                                 dict_indvars,dict_depvars,phi, derivative)

    bc_indvars = get_bc_varibles(bcs,dict_indvars,dict_depvars)
    symbolic_bc_loss_functions = [build_symbolic_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars,
                                                  phi, derivative;
                                                  bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
    symbolic_pde_loss_function,symbolic_bc_loss_functions
end
# Convert a PDE problem into an OptimizationProblem
function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)

    depvars = [nameof(value(d)) for d in pde_system.depvars]
    indvars = [nameof(value(i)) for i in pde_system.indvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)

    chain = discretization.chain
    initθ = discretization.initθ
    phi = discretization.phi
    autodiff = discretization.autodiff
    derivative = discretization.derivative
    autodiff == true && error("Automatic differentiation is not support yet")
    strategy = discretization.strategy

    _pde_loss_function = build_loss_function(eqs,indvars,depvars,
                                                 dict_indvars,dict_depvars,phi, derivative)

    bc_indvars = get_bc_varibles(bcs,dict_indvars,dict_depvars)
    _bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars,
                                                  phi, derivative;
                                                  bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    pde_loss_function, bc_loss_function =
    if strategy isa GridTraining
        dx = strategy.dx

        train_sets = generate_training_sets(domains,dx,bcs,
                                            dict_indvars,dict_depvars)

        # the points in the domain and on the boundary
        pde_train_set,bcs_train_set,train_set = train_sets

        pde_loss_function = get_loss_function(_pde_loss_function,
                                                        pde_train_set,
                                                        strategy)

        bc_loss_function = get_loss_function(_bc_loss_functions,
                                                       bcs_train_set,
                                                       strategy)

        (pde_loss_function, bc_loss_function)
    elseif strategy isa StochasticTraining
          bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
          pde_bounds, bcs_bounds = bounds
          pde_loss_function = get_loss_function(_pde_loss_function,
                                                          pde_bounds,
                                                          strategy)
          lbs,ubs = bcs_bounds
          number_of_points = length(lbs[1]) == 0 ? 1 : strategy.number_of_points^(1/length(lbs[1]))
          strategy = StochasticTraining(number_of_points = number_of_points)

          bc_loss_function = get_loss_function(_bc_loss_functions,
                                                         bcs_bounds,
                                                         strategy)
          (pde_loss_function, bc_loss_function)
    elseif strategy isa QuasiRandomTraining
         bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
         pde_bounds, bcs_bounds = bounds
         pde_loss_function = get_loss_function(_pde_loss_function,
                                                        pde_bounds,
                                                        strategy)
         plbs,pubs = pde_bounds
         blbs,bubs = bcs_bounds
         pl = length(plbs)
         bl = length(blbs[1])
         number_of_points = Int(round(strategy.number_of_points^(bl/pl)))
         strategy = QuasiRandomTraining(number_of_points = number_of_points)

         bc_loss_function = get_loss_function(_bc_loss_functions,
                                                       bcs_bounds,
                                                       strategy)
         (pde_loss_function, bc_loss_function)
    elseif strategy isa QuadratureTraining
        dim<=1 && error("QuadratureTraining works only with dimensionality more than 1")

        (dim<=2 && strategy.algorithm in [CubaCuhre(),CubaDivonne()]
        && error("$(strategy.algorithm) works only with dimensionality more than 2"))

        bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
        pde_bounds, bcs_bounds = bounds

        pde_loss_function = get_loss_function(_pde_loss_function,
                                              pde_bounds,
                                              strategy)

        bc_loss_function = get_loss_function(_bc_loss_functions,
                                             bcs_bounds,
                                             strategy)
        (pde_loss_function, bc_loss_function)
    end

    function loss_function_(θ,p)
        return pde_loss_function(θ) + bc_loss_function(θ)
    end

    f = OptimizationFunction(loss_function_, GalacticOptim.AutoZygote())
    prob = GalacticOptim.OptimizationProblem(f, initθ)
end


text:
#==============================================================================#
# IoTJobsDataPlane.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/iot-jobs-data-2017-09-29.normal.json
#==============================================================================#

__precompile__()

module IoTJobsDataPlane

using AWSCore

"""
    using AWSSDK.IoTJobsDataPlane.describe_job_execution
    describe_job_execution([::AWSConfig], arguments::Dict)
    describe_job_execution([::AWSConfig]; jobId=, thingName=, <keyword arguments>)

    using AWSCore.Services.iot_jobs_data
    iot_jobs_data([::AWSConfig], "GET", "/things/{thingName}/jobs/{jobId}", arguments::Dict)
    iot_jobs_data([::AWSConfig], "GET", "/things/{thingName}/jobs/{jobId}", jobId=, thingName=, <keyword arguments>)

# DescribeJobExecution Operation

Gets details of a job execution.

# Arguments

## `jobId = ::String` -- *Required*
The unique identifier assigned to this job when it was created.

## `thingName = ::String` -- *Required*
The thing name associated with the device the job execution is running on.

## `includeJobDocument = ::Bool`
Optional. When set to true, the response contains the job document. The default is false.

## `executionNumber = ::Int`
Optional. A number that identifies a particular job execution on a particular device. If not specified, the latest job execution is returned.

# Returns

`DescribeJobExecutionResponse`

# Exceptions

`InvalidRequestException`, `ResourceNotFoundException`, `ThrottlingException`, `ServiceUnavailableException`, `CertificateValidationException` or `TerminalStateException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-jobs-data-2017-09-29/DescribeJobExecution)
"""
@inline describe_job_execution(aws::AWSConfig=default_aws_config(); args...) = describe_job_execution(aws, args)

@inline describe_job_execution(aws::AWSConfig, args) = AWSCore.Services.iot_jobs_data(aws, "GET", "/things/{thingName}/jobs/{jobId}", args)

@inline describe_job_execution(args) = describe_job_execution(default_aws_config(), args)

"""
    using AWSSDK.IoTJobsDataPlane.get_pending_job_executions
    get_pending_job_executions([::AWSConfig], arguments::Dict)
    get_pending_job_executions([::AWSConfig]; thingName=)

    using AWSCore.Services.iot_jobs_data
    iot_jobs_data([::AWSConfig], "GET", "/things/{thingName}/jobs", arguments::Dict)
    iot_jobs_data([::AWSConfig], "GET", "/things/{thingName}/jobs", thingName=)

# GetPendingJobExecutions Operation

Gets the list of all jobs for a thing that are not in a terminal status.

# Arguments

## `thingName = ::String` -- *Required*
The name of the thing that is executing the job.

# Returns

`GetPendingJobExecutionsResponse`

# Exceptions

`InvalidRequestException`, `ResourceNotFoundException`, `ThrottlingException`, `ServiceUnavailableException` or `CertificateValidationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-jobs-data-2017-09-29/GetPendingJobExecutions)
"""
@inline get_pending_job_executions(aws::AWSConfig=default_aws_config(); args...) = get_pending_job_executions(aws, args)

@inline get_pending_job_executions(aws::AWSConfig, args) = AWSCore.Services.iot_jobs_data(aws, "GET", "/things/{thingName}/jobs", args)

@inline get_pending_job_executions(args) = get_pending_job_executions(default_aws_config(), args)

"""
    using AWSSDK.IoTJobsDataPlane.start_next_pending_job_execution
    start_next_pending_job_execution([::AWSConfig], arguments::Dict)
    start_next_pending_job_execution([::AWSConfig]; thingName=, <keyword arguments>)

    using AWSCore.Services.iot_jobs_data
    iot_jobs_data([::AWSConfig], "PUT", "/things/{thingName}/jobs/\$next", arguments::Dict)
    iot_jobs_data([::AWSConfig], "PUT", "/things/{thingName}/jobs/\$next", thingName=, <keyword arguments>)

# StartNextPendingJobExecution Operation

Gets and starts the next pending (status IN_PROGRESS or QUEUED) job execution for a thing.

# Arguments

## `thingName = ::String` -- *Required*
The name of the thing associated with the device.

## `statusDetails = ::Dict{String,String}`
A collection of name/value pairs that describe the status of the job execution. If not specified, the statusDetails are unchanged.

# Returns

`StartNextPendingJobExecutionResponse`

# Exceptions

`InvalidRequestException`, `ResourceNotFoundException`, `ThrottlingException`, `ServiceUnavailableException` or `CertificateValidationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-jobs-data-2017-09-29/StartNextPendingJobExecution)
"""
@inline start_next_pending_job_execution(aws::AWSConfig=default_aws_config(); args...) = start_next_pending_job_execution(aws, args)

@inline start_next_pending_job_execution(aws::AWSConfig, args) = AWSCore.Services.iot_jobs_data(aws, "PUT", "/things/{thingName}/jobs/\$next", args)

@inline start_next_pending_job_execution(args) = start_next_pending_job_execution(default_aws_config(), args)

"""
    using AWSSDK.IoTJobsDataPlane.update_job_execution
    update_job_execution([::AWSConfig], arguments::Dict)
    update_job_execution([::AWSConfig]; jobId=, thingName=, status=, <keyword arguments>)

    using AWSCore.Services.iot_jobs_data
    iot_jobs_data([::AWSConfig], "POST", "/things/{thingName}/jobs/{jobId}", arguments::Dict)
    iot_jobs_data([::AWSConfig], "POST", "/things/{thingName}/jobs/{jobId}", jobId=, thingName=, status=, <keyword arguments>)

# UpdateJobExecution Operation

Updates the status of a job execution.

# Arguments

## `jobId = ::String` -- *Required*
The unique identifier assigned to this job when it was created.

## `thingName = ::String` -- *Required*
The name of the thing associated with the device.

## `status = "QUEUED", "IN_PROGRESS", "SUCCEEDED", "FAILED", "REJECTED", "REMOVED" or "CANCELED"` -- *Required*
The new status for the job execution (IN_PROGRESS, FAILED, SUCCESS, or REJECTED). This must be specified on every update.

## `statusDetails = ::Dict{String,String}`
Optional. A collection of name/value pairs that describe the status of the job execution. If not specified, the statusDetails are unchanged.

## `expectedVersion = ::Int`
Optional. The expected current version of the job execution. Each time you update the job execution, its version is incremented. If the version of the job execution stored in Jobs does not match, the update is rejected with a VersionMismatch error, and an ErrorResponse that contains the current job execution status data is returned. (This makes it unnecessary to perform a separate DescribeJobExecution request in order to obtain the job execution status data.)

## `includeJobExecutionState = ::Bool`
Optional. When included and set to true, the response contains the JobExecutionState data. The default is false.

## `includeJobDocument = ::Bool`
Optional. When set to true, the response contains the job document. The default is false.

## `executionNumber = ::Int`
Optional. A number that identifies a particular job execution on a particular device.

# Returns

`UpdateJobExecutionResponse`

# Exceptions

`InvalidRequestException`, `ResourceNotFoundException`, `ThrottlingException`, `ServiceUnavailableException`, `CertificateValidationException` or `InvalidStateTransitionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-jobs-data-2017-09-29/UpdateJobExecution)
"""
@inline update_job_execution(aws::AWSConfig=default_aws_config(); args...) = update_job_execution(aws, args)

@inline update_job_execution(aws::AWSConfig, args) = AWSCore.Services.iot_jobs_data(aws, "POST", "/things/{thingName}/jobs/{jobId}", args)

@inline update_job_execution(args) = update_job_execution(default_aws_config(), args)

end # module IoTJobsDataPlane

#==============================================================================#
# End of file
#==============================================================================#

ext: jl


text:
"""
    PointCloud{T<:AbstractFloat, I<:Integer}(dataset::DataSet{T, I}, hyperrect::HyperRectVolume{T}, searchpts::Bool = false)::PointCloud

*AHMI-internal, not part of stable public API.*

creates a point cloud by searching the data tree for points which are inside the hyper-rectangle
The parameter searchpts determines if an array of the point IDs is created as well
"""
function PointCloud(
    dataset::DataSet{T, I},
    hyperrect::HyperRectVolume{T},
    searchpts::Bool)::PointCloud{T, I} where {T<:AbstractFloat, I<:Integer}

    result = PointCloud(T, I)

    PointCloud!(result, dataset, hyperrect, searchpts)

    result
end

function PointCloud!(
    cloud::PointCloud{T, I},
    dataset::DataSet{T, I},
    hyperrect::HyperRectVolume{T},
    searchpts::Bool) where {T<:AbstractFloat, I<:Integer}

    search!(cloud.searchres, dataset, hyperrect, searchpts)

    cloud.points = cloud.searchres.points

    resize!(cloud.pointIDs, cloud.points)
    copyto!(cloud.pointIDs, cloud.searchres.pointIDs)

    cloud.maxWeightProb = cloud.searchres.maxWeightProb
    cloud.minWeightProb = cloud.searchres.minWeightProb
    cloud.maxLogProb = cloud.searchres.maxLogProb
    cloud.minLogProb = cloud.searchres.minLogProb

    cloud.probfactor = exp(cloud.maxLogProb - cloud.minLogProb)
    cloud.probweightfactor = exp(cloud.maxWeightProb - cloud.minWeightProb)

    nothing
end

function create_pointweights(
    dataset::DataSet{T, I},
    volumes::Vector{IntegrationVolume{T, I, V}},
    ids::Array{I})::Vector{T} where {T<:AbstractFloat, I<:Integer, V<:SpatialVolume}

    pweights = zeros(T, dataset.N)

    for i in ids
        for p in eachindex(volumes[i].pointcloud.pointIDs)
            pweights[volumes[i].pointcloud.pointIDs[p]] += 1
        end
    end

    pweights
end

function Base.copy!(
    target::PointCloud{T, I},
    src::PointCloud{T, I}) where {T<:AbstractFloat, I<:Integer}

    target.points = src.points

    resize!(target.pointIDs, length(src.pointIDs))
    copyto!(target.pointIDs, src.pointIDs)

    target.maxLogProb = src.maxLogProb
    target.minLogProb = src.minLogProb

    target.maxWeightProb = src.maxWeightProb
    target.minWeightProb = src.minWeightProb

    target.probfactor = src.probfactor
    target.probweightfactor = src.probweightfactor

    nothing
end


text:
# mutable struct Point2DH <: FieldVector{3,Float64}
#     x::Float64
#     y::Float64
#     h::Float64
# end
#
# mutable struct Point3DH <: FieldVector{4,Float64}
#     x::Float64
#     y::Float64
#     z::Float64
#     h::Float64
# end

#const Point2DH = MMatrix{3,1,Float64,3}
#const Point3DH = MMatrix{4,1,Float64,4}

const Point2DH = SVector{3,Float64}
const Point3DH = SVector{4,Float64}

const Point2D = SVector{2,Float64}
const Point3D = SVector{3,Float64}

# mutable struct Point2DH <: MMatrix{3,1,Float64,3}
#     x::Float64
#     y::Float64
#     h::Float64
# end
#
# mutable struct Point3DH <: MMatrix{4,1,Float64,4}
#     x::Float64
#     y::Float64
#     z::Float64
#     h::Float64
# end

struct HomogeneousPoint{T <: AbstractFloat,N}
    coords::NTuple{N, T}
end

abstract type ProjectiveEntity end

abstract type CameraModel end

abstract type EstimationAlgorithm end

abstract type FactorisationAlgorithm end

abstract type CostFunction end

abstract type CoordinateSystemTransformation end

abstract type CovarianceEstimationScheme end

abstract type NoiseModel end

abstract type TotalViews end

abstract type PlanarCoordinateSystem end

mutable struct FundamentalMatrix <: ProjectiveEntity
end

mutable struct HomographyMatrix <: ProjectiveEntity
end

mutable struct HomographyMatrices <: ProjectiveEntity
end

mutable struct EssentialMatrix <: ProjectiveEntity
end

mutable struct  ProjectionMatrix <: ProjectiveEntity
end

mutable struct  ProjectionMatrices <: ProjectiveEntity
    src::ProjectiveEntity
    method::FactorisationAlgorithm
    views::TotalViews
end

mutable struct  LatentVariables <: ProjectiveEntity
    src::ProjectiveEntity
end

mutable struct  Chojnacki <: FactorisationAlgorithm
end

mutable struct  HomogeneousCoordinates <: ProjectiveEntity
end

struct CartesianSystem <: PlanarCoordinateSystem
    𝐞₁::Vec{2,Float64}
    𝐞₂::Vec{2,Float64}
end
# Convention I:
# CartesianSystem() = CartesianSystem(Vec(1.0, 0.0), Vec(0.0, -1.0))
CartesianSystem() = CartesianSystem(Vec(-1.0, 0.0), Vec(0.0, 1.0))

struct RasterSystem <: PlanarCoordinateSystem
    𝐞₁::Vec{2,Float64}
    𝐞₂::Vec{2,Float64}
end
#TODO Will have to change this to accomodate different conventions.
# Convention I:  RasterSystem() = RasterSystem(Vec(1.0, 0.0), Vec(0.0, 1.0))
RasterSystem() = RasterSystem(Vec(-1.0, 0.0), Vec(0.0, -1.0))

struct OpticalSystem <: PlanarCoordinateSystem
    𝐞₁::Vec{2,Float64}
    𝐞₂::Vec{2,Float64}
end
# Convention I: = OpticalSystem(Vec(1.0, 0.0), Vec(0.0, 1.0))
OpticalSystem() = OpticalSystem(Vec(-1.0, 0.0), Vec(0.0, -1.0))

mutable struct  Pinhole <: CameraModel
    image_width::Int
    image_height::Int
    focal_length::Float64
    # Center of projection.
    𝐜::Point{3,Float64}
    # Basis vectors that characterise the pose of the camera
    𝐞₁::Vec{3,Float64}
    𝐞₂::Vec{3,Float64}
    𝐞₃::Vec{3,Float64}
    # Origin of the picture plane (the image).
    𝐨::Point{2,Float64}
    # Basis vectors that characterise the coordinate system of the
    # picture plane (the image).
    𝐞₁′::Vec{2,Float64}
    𝐞₂′::Vec{2,Float64}
end

mutable struct  CanonicalLens <: CameraModel
end

mutable struct DirectLinearTransform <: EstimationAlgorithm
end

mutable struct BundleAdjustment{A<:AbstractVector} <: EstimationAlgorithm
    𝛉₀::A
    max_iter::Int64
    toleranceθ::Float64
end

mutable struct FundamentalNumericalScheme{A<:AbstractVector} <: EstimationAlgorithm
    𝛉₀::A
    max_iter::Int64
    toleranceθ::Float64
end

mutable struct CanonicalApproximation <: CovarianceEstimationScheme
end

mutable struct HessianApproximation <: CovarianceEstimationScheme
end

mutable struct Taubin <: EstimationAlgorithm
end

mutable struct ApproximateMaximumLikelihood <: CostFunction
end
const AML = ApproximateMaximumLikelihood

mutable struct CanonicalToHartley <: CoordinateSystemTransformation
end

mutable struct HartleyToCanonical <: CoordinateSystemTransformation
end

mutable struct CovarianceMatrices
end

mutable struct GaussianNoise <: NoiseModel

end

mutable struct TwoViews <: TotalViews
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    boundgrid(spatialdata, dims)

Returns a `RegularGrid` of given `dims` covering all the
locations in `spatialdata`.
"""
function boundgrid(spatialdata::AbstractSpatialData, dims::Dims)
  # retrieve data coordinates
  datacoords = coordinates(spatialdata)
  N = length(datacoords)

  @assert length(dims) == N "dimensions must match number of coordinates in data"
  @assert all(dims .> 0) "dimensions must be positive"

  # determine coordinate type
  T = promote_type([T for (var,T) in datacoords]...)

  start  = fill(typemax(T), N)
  finish = fill(typemin(T), N)

  for ind in 1:npoints(spatialdata)
    x = coordinates(spatialdata, ind)
    for d in 1:N
      x[d] < start[d] && (start[d] = x[d])
      x[d] > finish[d] && (finish[d] = x[d])
    end
  end

  RegularGrid(tuple(start...), tuple(finish...), dims=dims)
end

"""
    readgeotable(args; coordnames=[:x,:y,:z], kwargs)

Read data from disk using `CSV.read`, optionally specifying
the columns `coordnames` with spatial coordinates.

The arguments `args` and keyword arguments `kwargs` are
forwarded to the `CSV.read` function, please check their
documentation for more details.

This function returns a [`GeoDataFrame`](@ref) object.
"""
function readgeotable(args...; coordnames=[:x,:y,:z], kwargs...)
  data = CSV.read(args...; kwargs...)
  GeoDataFrame(data, coordnames)
end


text:
@doc raw"""
    get_reason(c)

return the current reason stored within a [`StoppingCriterion`](@ref) `c`.
This reason is empty if the criterion has never been met.
"""
get_reason(c::sC) where {sC<:StoppingCriterion} = c.reason

@doc raw"""
    get_stopping_criteria(c)
return the array of internally stored [`StoppingCriterion`](@ref)s for a
[`StoppingCriterionSet`](@ref) `c`.
"""
function get_stopping_criteria(c::S) where {S<:StoppingCriterionSet}
    return error("get_stopping_criteria() not defined for a $(typeof(c)).")
end

@doc raw"""
    StopAfterIteration <: StoppingCriterion

A functor for an easy stopping criterion, i.e. to stop after a maximal number
of iterations.

# Fields
* `maxIter` – stores the maximal iteration number where to stop at
* `reason` – stores a reason of stopping if the stopping criterion has one be
  reached, see [`get_reason`](@ref).

# Constructor

    StopAfterIteration(maxIter)

initialize the stopafterIteration functor to indicate to stop after `maxIter`
iterations.
"""
mutable struct StopAfterIteration <: StoppingCriterion
    maxIter::Int
    reason::String
    StopAfterIteration(mIter::Int) = new(mIter, "")
end
function (c::StopAfterIteration)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if i > c.maxIter
        c.reason = "The algorithm reached its maximal number of iterations ($(c.maxIter)).\n"
        return true
    end
    return false
end
"""
    StopWhenGradientNormLess <: StoppingCriterion

stores a threshold when to stop looking at the norm of the gradient from within
a [`GradientProblem`](@ref).
"""
mutable struct StopWhenGradientNormLess <: StoppingCriterion
    threshold::Float64
    reason::String
    StopWhenGradientNormLess(ε::Float64) = new(ε, "")
end
function (c::StopWhenGradientNormLess)(p::P, o::O, ::Int) where {P<:Problem,O<:Options}
    if norm(p.M, o.x, get_gradient(p, o.x)) < c.threshold
        c.reason = "The algorithm reached approximately critical point; the gradient norm ($(norm(p.M,o.x,get_gradient(p,o.x)))) is less than $(c.threshold).\n"
        return true
    end
    return false
end
"""
    StopWhenChangeLess <: StoppingCriterion

stores a threshold when to stop looking at the norm of the change of the
optimization variable from within a [`Options`](@ref), i.e `o.x`.
For the storage a [`StoreOptionsAction`](@ref) is used

# Constructor

    StopWhenChangeLess(ε[, a])

initialize the stopping criterion to a threshold `ε` using the
[`StoreOptionsAction`](@ref) `a`, which is initialized to just store `:x` by
default.
"""
mutable struct StopWhenChangeLess <: StoppingCriterion
    threshold::Float64
    reason::String
    storage::StoreOptionsAction
    function StopWhenChangeLess(ε::Float64, a::StoreOptionsAction=StoreOptionsAction((:x,)))
        return new(ε, "", a)
    end
end
function (c::StopWhenChangeLess)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if has_storage(c.storage, :x)
        xOld = get_storage(c.storage, :x)
        if distance(p.M, o.x, xOld) < c.threshold && i > 0
            c.reason = "The algorithm performed a step with a change ($(distance(p.M, o.x, xOld))) less than $(c.threshold).\n"
            c.storage(p, o, i)
            return true
        end
    end
    c.storage(p, o, i)
    return false
end
@doc raw"""
    StopWhenAll <: StoppingCriterion

store an array of [`StoppingCriterion`](@ref) elements and indicates to stop,
when _all_ indicate to stop. The `reseason` is given by the concatenation of all
reasons.

# Constructor
    StopWhenAll(c::NTuple{N,StoppingCriterion} where N)
    StopWhenAll(c::StoppingCriterion,...)
"""
mutable struct StopWhenAll{TCriteria<:Tuple} <: StoppingCriterionSet
    criteria::TCriteria
    reason::String
    StopWhenAll(c::Vector{StoppingCriterion}) = new{typeof(tuple(c...))}(tuple(c...), "")
    StopWhenAll(c...) = new{typeof(c)}(c, "")
end
function (c::StopWhenAll)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if all(subC -> subC(p, o, i), c.criteria)
        c.reason = string([get_reason(subC) for subC in c.criteria]...)
        return true
    end
    return false
end
get_stopping_criteria(c::StopWhenAll) = c.criteria

@doc raw"""
    StopWhenAny <: StoppingCriterion

store an array of [`StoppingCriterion`](@ref) elements and indicates to stop,
when _any_ single one indicates to stop. The `reseason` is given by the
concatenation of all reasons (assuming that all non-indicating return `""`).

# Constructor
    StopWhenAny(c::NTuple{N,StoppingCriterion} where N)
    StopWhenAny(c::StoppingCriterion...)
"""
mutable struct StopWhenAny{TCriteria<:Tuple} <: StoppingCriterionSet
    criteria::TCriteria
    reason::String
    StopWhenAny(c::Vector{StoppingCriterion}) = new{typeof(tuple(c...))}(tuple(c...), "")
    StopWhenAny(c::StoppingCriterion...) = new{typeof(c)}(c)
end
function (c::StopWhenAny)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if any(subC -> subC(p, o, i), c.criteria)
        c.reason = string([get_reason(subC) for subC in c.criteria]...)
        return true
    end
    return false
end
get_stopping_criteria(c::StopWhenAny) = c.criteria
"""
    StopWhenCostLess <: StoppingCriterion

store a threshold when to stop looking at the cost function of the
optimization problem from within a [`Problem`](@ref), i.e `get_cost(p,o.x)`.

# Constructor

    StopWhenCostLess(ε)

initialize the stopping criterion to a threshold `ε`.
"""
mutable struct StopWhenCostLess <: StoppingCriterion
    threshold::Float64
    reason::String
    StopWhenCostLess(ε::Float64) = new(ε, "")
end
function (c::StopWhenCostLess)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if i > 0 && get_cost(p, o.x) < c.threshold
        c.reason = "The algorithm reached a cost function value ($(get_cost(p,o.x))) less then the threshold ($(c.threshold)).\n"
        return true
    end
    return false
end

"""
    StopAfter <: StoppingCriterion

store a threshold when to stop looking at the complete runtime. It uses
`time_ns()` to measure the time and you provide a `Period` as a time limit,
i.e. `Minute(15)`

# Constructor

    StopAfter(t)

initialize the stopping criterion to a `Period t` to stop after.
"""
mutable struct StopAfter <: StoppingCriterion
    threshold::Period
    reason::String
    start::Nanosecond
    function StopAfter(t::Period)
        return if value(t) < 0
            error("You must provide a positive time period")
        else
            new(t, "", Nanosecond(0))
        end
    end
end
function (c::StopAfter)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if value(c.start) == 0 || i <= 0 # (re)start timer
        c.start = Nanosecond(time_ns())
    else
        cTime = Nanosecond(time_ns()) - c.start
        if i > 0 && (cTime > Nanosecond(c.threshold))
            c.reason = "The algorithm ran for about $(floor(cTime, typeof(c.threshold))) and has hence reached the threshold of $(c.threshold).\n"
            return true
        end
    end
    return false
end
@doc raw"""
    get_active_stopping_criteria(c)

returns all active stopping criteria, if any, that are within a
[`StoppingCriterion`](@ref) `c`, and indicated a stop, i.e. their reason is
nonempty.
To be precise for a simple stopping criterion, this returns either an empty
array if no stop is incated or the stopping criterion as the only element of
an array. For a [`StoppingCriterionSet`](@ref) all internal (even nested)
criteria that indicate to stop are returned.
"""
function get_active_stopping_criteria(c::sCS) where {sCS<:StoppingCriterionSet}
    c = get_active_stopping_criteria.(get_stopping_criteria(c))
    return vcat(c...)
end
# for non-array containing stopping criteria, the recursion ends in either
# returning nothing or an 1-element array contianing itself
function get_active_stopping_criteria(c::sC) where {sC<:StoppingCriterion}
    if c.reason != ""
        return [c] # recursion top
    else
        return []
    end
end


text:
# -------------------------------------------------------------------------------------
_tryparse(T, v) = tryparse(T, v)
_tryparse(::Type{String}, v) = string(v)

_parse_key(k) = string(k)
function _parse_val(v::AbstractString)
    for T in _OUTPUT_VAL_TYPES
        r = _tryparse(T, v)
        !isnothing(r) && return r
    end
    return string(v)
end

function _parse_pair(p::AbstractString)
    k, v = split(p, _SEPS[:PAIR_SEP])
    return (_parse_key(k) => _parse_val(v))
end

# -------------------------------------------------------------------------------------
function isvalid_dfname(dfn::String)
    dfn = basename(dfn)
    m = _parse_regex(dfn)
    isempty(m.digest)
end

# -------------------------------------------------------------------------------------
function _parse_dfname(dfn::String, ondigest::Function)
    dfn = basename(dfn)
    head_str, params_str, ext_str, digest = _parse_regex(dfn)
    !isempty(digest) && return ondigest(dfn, digest)

    keepempty = false
    head_split = split(head_str, _SEPS[:ELEMT_SEP]; keepempty) .|> string
    params_split = split(params_str, _SEPS[:ELEMT_SEP]; keepempty) .|> string

    return (;
        head = _parse_val.(head_split),
        params = Dict{String, Any}(_parse_pair.(params_split)...),
        ext = string(ext_str)
    )
end

_error_ondigest(dfn, digest) = error("Invalid name '", dfn, "'. Digest: '", digest, "'")
parse_dfname(dfn::String) = _parse_dfname(dfn::String, _error_ondigest)

_nothing_ondigest(dfn, digest) = nothing
tryparse_dfname(dfn::String) = _parse_dfname(dfn::String, _nothing_ondigest)

## ------------------------------------------------------
dfheads(dfname::String) = parse_dfname(dfname)[1]
dfparams(dfname::String) = parse_dfname(dfname)[2]
dfparam(dfname::String, k::String) = dfparams(dfname)[k]
dfparam(dfname::String, k::String, dft) = get(dfparams(dfname), k, dft)
dfext(dfname::String) = parse_dfname(dfname)[3]


text: 1 || false

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# Online MM algorithm via quadratic upper bound
# Majorizing function: f(u) + ∇f(u)'(Θ - u) + .5 * (Θ - u)'H(Θ - u)
# Q_t(Θ) = (Θ' * A * Θ) + (b' * Θ) + c
# Where:
#   A = (1 - γ) * A + γ * H
#   b = (1 - γ) * b + γ * (∇f(Θ) - H * Θ)

immutable MMQuadUpBound
    A::MatF
    b::VecF
end
function fit!(o::MMQuadUpBound, A, b, γ)
    smooth!(o.A, A, γ)
    smooth!(o.b, b, γ)
end
value(o::MMQuadUpBound) = o.A \ o.b


text:
push!(ARGS, "../../input_files/dynamic/basin_refinement/8GPa/12n.dat")
include("../../Basin.jl")
pop!(ARGS)


text:
""" Mesh Tools
This file includes tools for importing delta sets from mesh
files supported by MeshIO and for converting delta sets to
meshes (for the purposes of plotting.

Meshes here are stored in the GeometryBasics.Mesh object.
"""
module MeshInterop
export EmbeddedDeltaSet2D, Mesh

using FileIO, MeshIO, GeometryBasics
import GeometryBasics: Mesh

using Catlab.CategoricalAlgebra: copy_parts!
using ..SimplicialSets, ..DualSimplicialSets
import ..SimplicialSets: EmbeddedDeltaSet2D

# Helper Functions (should these be exposed?)
##################

""" Construct a GeometryBasics.Mesh object from an embedded delta set
"""
function Mesh(ds::EmbeddedDeltaSet2D)
  points = Point{3, Float64}[point(ds)...]
  tris = TriangleFace{Int}[zip(triangle_vertices(ds)...)...]
  GeometryBasics.Mesh(points, tris)
end

""" Construct a GeometryBasics.Mesh object from a dual embedded delta set
"""
function Mesh(ds::EmbeddedDeltaDualComplex2D{O,R,P}; primal=false) where {O,R,P}
  if(primal)
    ds′ = EmbeddedDeltaSet2D{O,P}()
    copy_parts!(ds′, ds)
    Mesh(ds′)
  else
    points = Point{3, Float64}[dual_point(ds)...]
    tris = TriangleFace{Int}[zip(dual_triangle_vertices(ds)...)...]
    GeometryBasics.Mesh(points, tris)
  end
end

# Import Tooling
################

""" Constructor for EmbeddedDeltaSet2D from GeometryBasics.Mesh object

This operator should work for any triangular mesh object. Note that it will
not preserve any normal, texture, or other attributes from the Mesh object.
"""
function EmbeddedDeltaSet2D(m::GeometryBasics.Mesh)
  coords = coordinates(m)
  tris = faces(m)
  s = EmbeddedDeltaSet2D{Bool, eltype(coords)}()
  add_vertices!(s, length(coords), point=coords)
  for tri in tris
    tri = convert.(Int64, tri)
    glue_sorted_triangle!(s, tri...)
  end
  # Properly orient the delta set
  orient!(s)
  s
end

""" Constructor for EmbeddedDeltaSet2D from mesh file

This operator should work for any file support for import from MeshIO. Note
that it will not preserve any normal, texture, or other data beyond points and
triangles from the mesh file.
"""
function EmbeddedDeltaSet2D(fn::String)
	EmbeddedDeltaSet2D(FileIO.load(fn))
end
end


text:
# This file includes two kinds of codes
#   - Codes for backward compatibility
#   - Glue codes that might no longer be necessary in the future

@static if VERSION < v"1.1"
    @inline isnothing(x) = x === nothing
end

if hasmethod(nan, Tuple{Type{HSV{Float32},}})
    # requires ColorTypes v0.11 and ColorVectorSpace v0.9.4
    # https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/75
    @inline _nan(::Type{T}) where T = nan(T)
else
    @inline _nan(::Type{HSV{Float16}}) = HSV{Float16}(NaN16,NaN16,NaN16)
    @inline _nan(::Type{HSV{Float32}}) = HSV{Float32}(NaN32,NaN32,NaN32)
    @inline _nan(::Type{HSV{Float64}}) = HSV{Float64}(NaN,NaN,NaN)
    @inline _nan(::Type{T}) where {T} = nan(T)
end