text:
import Base.Broadcast
Base.Broadcast.dottable(x::Function) = true
RuntimeGeneratedFunctions.init(@__MODULE__)
"""
Algorithm for solving Physics-Informed Neural Networks problems.

Arguments:
* `chain`: a Flux.jl chain with a d-dimensional input and a 1-dimensional output,
* `strategy`: determines which training strategy will be used,
* `init_params`: the initial parameter of the neural network,
* `phi`: a trial solution,
* `derivative`: method that calculates the derivative.

"""
abstract type AbstractPINN{isinplace} <: SciMLBase.SciMLProblem end

struct PhysicsInformedNN{isinplace,C,T,P,PH,DER,PE,AL,K} <: AbstractPINN{isinplace}
  chain::C
  strategy::T
  init_params::P
  phi::PH
  derivative::DER
  param_estim::PE
  additional_loss::AL
  kwargs::K

 @add_kwonly function PhysicsInformedNN{iip}(chain,
                                             strategy;
                                             init_params = nothing,
                                             phi = nothing,
                                             derivative = nothing,
                                             param_estim=false,
                                             additional_loss=nothing,
                                             kwargs...) where iip
        if init_params == nothing
            if chain isa AbstractArray
                initθ = DiffEqFlux.initial_params.(chain)
            else
                initθ = DiffEqFlux.initial_params(chain)
            end

        else
            initθ = init_params
        end

        if phi == nothing
            if chain isa AbstractArray
                _phi = get_phi.(chain)
            else
                _phi = get_phi(chain)
            end
        else
            _phi = phi
        end

        if derivative == nothing
            _derivative = get_numeric_derivative()
        else
            _derivative = derivative
        end
        new{iip,typeof(chain),typeof(strategy),typeof(initθ),typeof(_phi),typeof(_derivative),typeof(param_estim),typeof(additional_loss),typeof(kwargs)}(chain,strategy,initθ,_phi,_derivative,param_estim,additional_loss,kwargs)
    end
end
PhysicsInformedNN(chain,strategy,args...;kwargs...) = PhysicsInformedNN{true}(chain,strategy,args...;kwargs...)

SciMLBase.isinplace(prob::PhysicsInformedNN{iip}) where iip = iip

abstract type TrainingStrategies  end

"""
* `dx`: the discretization of the grid.
"""
struct GridTraining <: TrainingStrategies
    dx
end

"""
* `points`: number of points in random select training set.
"""
struct StochasticTraining <:TrainingStrategies
    points:: Int64
end

"""
* `points`:  the number of quasi-random points in minibatch,
* `sampling_alg`: the quasi-Monte Carlo sampling algorithm,
* `minibatch`: the number of subsets.

For more information look: QuasiMonteCarlo.jl https://github.com/SciML/QuasiMonteCarlo.jl
"""
struct QuasiRandomTraining <:TrainingStrategies
    points:: Int64
    sampling_alg::QuasiMonteCarlo.SamplingAlgorithm
    minibatch:: Int64
end
function QuasiRandomTraining(points;sampling_alg = UniformSample(),minibatch=500)
    QuasiRandomTraining(points,sampling_alg,minibatch)
end
"""
* `quadrature_alg`: 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
    quadrature_alg::DiffEqBase.AbstractQuadratureAlgorithm
    reltol::Float64
    abstol::Float64
    maxiters::Int64
    batch::Int64
end

function QuadratureTraining(;quadrature_alg=CubatureJLh(),reltol= 1e-6,abstol= 1e-3,maxiters=1e3,batch=100)
    QuadratureTraining(quadrature_alg,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_expression
function transform_expression(ex,dict_indvars,dict_depvars, initθ, strategy)
    if ex isa Expr
        ex = _transform_expression(ex,dict_indvars,dict_depvars,initθ,strategy)
    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 - function
 x,y - coordinates of point
 εs - epsilon mask
 order - order of derivative
 θ - weight in neural network
"""
function _transform_expression(ex,dict_indvars,dict_depvars, initθ, strategy)
    _args = ex.args
    for (i,e) in enumerate(_args)
        if !(e isa Expr)
            if e in keys(dict_depvars)
                depvar = _args[1]
                num_depvar = dict_depvars[depvar]
                indvars = _args[2:end]
                cord = :cord
                ex.args = if length(dict_depvars) == 1
                    [:u, cord, :($θ), :phi]
                else
                    [:u, cord, Symbol(:($θ),num_depvar), Symbol(:phi,num_depvar)]
                end
                break
            elseif e isa ModelingToolkit.Differential
                derivative_variables = Symbol[]
                order = 0
                while (_args[1] isa ModelingToolkit.Differential)
                    order += 1
                    push!(derivative_variables, toexpr(_args[1].x))
                    _args = _args[2].args
                end
                depvar = _args[1]
                num_depvar = dict_depvars[depvar]
                indvars = _args[2:end]
                cord = :cord
                dim_l = length(indvars)
                εs = [get_ε(dim_l,d) for d in 1:dim_l]
                undv = [dict_indvars[d_p] for d_p  in derivative_variables]
                εs_dnv = [εs[d] for d in undv]
                ex.args = if length(dict_depvars) == 1
                    [:derivative, :phi, :u, cord, εs_dnv, order, :($θ)]
                else
                    [:derivative, Symbol(:phi,num_depvar), :u, cord, εs_dnv, order, Symbol(:($θ),num_depvar)]
                end
                break
            end
        else
            ex.args[i] = _transform_expression(ex.args[i],dict_indvars,dict_depvars,initθ,strategy)
        end
    end
    return ex
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, [x, y], [[ε,0],[ε,0]], 2, θ) + derivative(phi, u, [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(phi1, u1, [x, y], [[ε,0]], 1, θ1) + 4 * derivative(phi2, u, [x, y], [[0,ε]], 1, θ2)) - 0,
       (derivative(phi2, u2, [x, y], [[ε,0]], 1, θ2) + 9 * derivative(phi1, u, [x, y], [[0,ε]], 1, θ1)) - 0]
"""

function build_symbolic_equation(eq,_indvars,_depvars,initθ,strategy)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    parse_equation(eq,dict_indvars,dict_depvars,initθ,strategy)
end

function parse_equation(eq,dict_indvars,dict_depvars,initθ, strategy)
    eq_lhs = isequal(expand_derivatives(eq.lhs), 0) ? eq.lhs : expand_derivatives(eq.lhs)
    eq_rhs = isequal(expand_derivatives(eq.rhs), 0) ? eq.rhs : expand_derivatives(eq.rhs)

    left_expr = transform_expression(toexpr(eq_lhs),dict_indvars,dict_depvars,initθ,strategy)
    right_expr = transform_expression(toexpr(eq_rhs),dict_indvars,dict_depvars,initθ,strategy)
    left_expr = Broadcast.__dot__(left_expr)
    right_expr = Broadcast.__dot__(right_expr)
    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

:((cord, θ, phi, derivative, u)->begin
          #= ... =#
          #= ... =#
          begin
              (θ1, θ2) = (θ[1:33], θ"[34:66])
              (phi1, phi2) = (phi[1], phi[2])
              let (x, y) = (cord[1], cord[2])
                  [(+)(derivative(phi1, u, [x, y], [[ε, 0.0]], 1, θ1), (*)(4, derivative(phi2, u, [x, y], [[0.0, ε]], 1, θ2))) - 0,
                   (+)(derivative(phi2, u, [x, y], [[ε, 0.0]], 1, θ2), (*)(9, derivative(phi1, u, [x, y], [[0.0, ε]], 1, θ1))) - 0]
              end
          end
      end)
"""
function build_symbolic_loss_function(eqs,_indvars,_depvars, phi, derivative,initθ,strategy; bc_indvars=nothing , eq_params = SciMLBase.NullParameters(), param_estim = false,default_p=nothing)
    # dictionaries: variable -> unique number
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_symbolic_loss_function(eqs,indvars,depvars,
                                        dict_indvars,dict_depvars,
                                        phi, derivative,initθ,strategy,
                                        bc_indvars = bc_indvars, eq_params = eq_params, param_estim = param_estim,default_p=default_p)
end

function get_indvars_ex(bc_indvars)
    i_=1
    indvars_ex = map(bc_indvars) do u
        if u isa Symbol
             ex = :($:cord[[$i_],:])
             i_+=1
             ex
        else
           :(fill($u,size($:cord[[1],:])))
        end
    end
    indvars_ex
end

function build_symbolic_loss_function(eqs,indvars,depvars,
                                      dict_indvars,dict_depvars,
                                      phi, derivative, initθ, strategy; eq_params = SciMLBase.NullParameters(), param_estim = param_estim,default_p=default_p,
                                      bc_indvars = indvars)

    loss_function = parse_equation(eqs,dict_indvars,dict_depvars,initθ,strategy)
    vars = :(cord, $θ, phi, derivative,u,p)
    ex = Expr(:block)
    if length(depvars) != 1
        θ_nums = Symbol[]
        phi_nums = Symbol[]
        for v in depvars
            num = dict_depvars[v]
            push!(θ_nums,:($(Symbol(:($θ),num))))
            push!(phi_nums,:($(Symbol(:phi,num))))
        end

        expr_θ = Expr[]
        expr_phi = Expr[]

        acum =  [0;accumulate(+, length.(initθ))]
        sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]

        for i in eachindex(depvars)
            push!(expr_θ, :($θ[$(sep[i])]))
            push!(expr_phi, :(phi[$i]))
        end

        vars_θ = Expr(:(=), build_expr(:tuple, θ_nums), build_expr(:tuple, expr_θ))
        push!(ex.args,  vars_θ)

        vars_phi = Expr(:(=), build_expr(:tuple, phi_nums), build_expr(:tuple, expr_phi))
        push!(ex.args,  vars_phi)
    end
    #Add an expression for parameter symbols
    if param_estim == true && eq_params != SciMLBase.NullParameters()
        param_len = length(eq_params)
        last_indx =  [0;accumulate(+, length.(initθ))][end]
        params_symbols = Symbol[]
        expr_params = Expr[]
        for (i , eq_param) in enumerate(eq_params)
            push!(expr_params, :($θ[$(i+last_indx:i+last_indx)]))
            push!(params_symbols, Symbol(:($eq_param)))
        end
        params_eq = Expr(:(=), build_expr(:tuple, params_symbols), build_expr(:tuple, expr_params))
        push!(ex.args,  params_eq)
    end

    if eq_params != SciMLBase.NullParameters() && param_estim == false
        params_symbols = Symbol[]
        expr_params = Expr[]
        for (i , eq_param) in enumerate(eq_params)
            push!(expr_params, :(ArrayInterface.allowed_getindex(p,$i:$i)))
            push!(params_symbols, Symbol(:($eq_param)))
        end
        params_eq = Expr(:(=), build_expr(:tuple, params_symbols), build_expr(:tuple, expr_params))
        push!(ex.args,  params_eq)
    end

    if strategy isa QuadratureTraining

        indvars_ex = get_indvars_ex(bc_indvars)

        left_arg_pairs, right_arg_pairs = indvars,indvars_ex
        vcat_expr =  :(cord = vcat($(indvars...)))
        vcat_expr_loss_functions = Expr(:block,vcat_expr,loss_function) #TODO rename
        vars_eq = Expr(:(=), build_expr(:tuple, left_arg_pairs), build_expr(:tuple, right_arg_pairs))

    else
        indvars_ex = [:($:cord[[$i],:]) for (i, u) ∈ enumerate(indvars)]
        left_arg_pairs, right_arg_pairs = indvars,indvars_ex
        vars_eq = Expr(:(=), build_expr(:tuple, left_arg_pairs), build_expr(:tuple, right_arg_pairs))
        vcat_expr_loss_functions = loss_function #TODO rename
    end

    let_ex = Expr(:let, vars_eq, vcat_expr_loss_functions)
    push!(ex.args,  let_ex)

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

function build_loss_function(eqs,_indvars,_depvars, phi, derivative,initθ,strategy;bc_indvars=nothing,eq_params=SciMLBase.NullParameters(),param_estim=false,default_p=nothing)
    # dictionaries: variable -> unique number
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    bc_indvars = bc_indvars==nothing ? indvars : bc_indvars
    return build_loss_function(eqs,indvars,depvars,
                               dict_indvars,dict_depvars,
                               phi, derivative,initθ,strategy,
                               bc_indvars = bc_indvars, eq_params=eq_params,param_estim=param_estim,default_p=default_p)
end

function build_loss_function(eqs,indvars,depvars,
                             dict_indvars,dict_depvars,
                             phi, derivative, initθ, strategy;
                             bc_indvars = indvars,eq_params=SciMLBase.NullParameters(),param_estim=false,default_p=nothing)
     expr_loss_function = build_symbolic_loss_function(eqs,indvars,depvars,
                                                       dict_indvars,dict_depvars,
                                                       phi, derivative, initθ,strategy;
                                                       bc_indvars = bc_indvars,eq_params = eq_params,param_estim=param_estim,default_p=default_p)
    u = get_u()
    _loss_function = @RuntimeGeneratedFunction(expr_loss_function)
    loss_function = (cord, θ) -> begin
        _loss_function(cord, θ, phi, derivative, u, default_p)
    end
    return loss_function
end

function get_vars(indvars_, depvars_)
    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 depvars,indvars,dict_indvars,dict_depvars
end

function get_variables(eqs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return get_variables(eqs,dict_indvars,dict_depvars)
end

function get_variables(eqs,dict_indvars,dict_depvars)
    bc_args = get_argument(eqs,dict_indvars,dict_depvars)
    return map(barg -> filter(x -> x isa Symbol, barg), bc_args)
end

function get_number(eqs,dict_indvars,dict_depvars)
    bc_args = get_argument(eqs,dict_indvars,dict_depvars)
    return map(barg -> filter(x -> x isa Number, barg), bc_args)
end

function find_thing_in_expr(ex::Expr, thing; ans = Expr[])
    for e in ex.args
        if e isa Expr
            if thing in e.args
                push!(ans,e)
            end
            find_thing_in_expr(e,thing; ans=ans)
        end
    end
    return collect(Set(ans))
end

# Get arguments from boundary condition functions
function get_argument(eqs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    get_argument(eqs,dict_indvars,dict_depvars)
end
function get_argument(eqs,dict_indvars,dict_depvars)
    exprs = toexpr.(eqs)
    vars = map(exprs) do expr
        _vars =  map(depvar -> find_thing_in_expr(expr,  depvar), collect(keys(dict_depvars)))
        f_vars = filter(x -> !isempty(x), _vars)
        map(x -> first(x), f_vars)
    end
    args_ = map(vars) do _vars
        map(var -> var.args[2:end] , _vars)
    end
    return first.(args_) #TODO for all arguments
end

function generate_training_sets(domains,dx,eqs,bcs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return generate_training_sets(domains,dx,eqs,bcs,dict_indvars,dict_depvars)
end
# Generate training set in the domain and on the boundary
function generate_training_sets(domains,dx,eqs,bcs,dict_indvars::Dict,dict_depvars::Dict)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end

    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)])

    bound_args = get_argument(bcs,dict_indvars,dict_depvars)
    bound_vars = get_variables(bcs,dict_indvars,dict_depvars)

    dif = [Float32[] 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
    cord_train_set = collect.(spans)
    bc_data = map(zip(dif,cord_train_set)) do (d,c)
        setdiff(c, d)
    end

    dict_var_span_ = Dict([Symbol(d.variables) => bc for (d,bc) in zip(domains,bc_data)])

    bcs_train_sets = map(bound_args) do bt
        span = map(b -> get(dict_var_span, b, b), bt)
        _set = Float32.(hcat(vec(map(points -> collect(points), Iterators.product(span...)))...))
    end

    pde_vars = get_variables(eqs,dict_indvars,dict_depvars)
    pde_args = get_argument(eqs,dict_indvars,dict_depvars)

    pde_train_set = Float32.(hcat(vec(map(points -> collect(points), Iterators.product(bc_data...)))...))

    pde_train_sets = map(pde_args) do bt
        span = map(b -> get(dict_var_span_, b, b), bt)
        _set = Float32.(hcat(vec(map(points -> collect(points), Iterators.product(span...)))...))
    end
    [pde_train_sets,bcs_train_sets]
end

function get_bounds(domains,eqs,bcs,_indvars::Array,_depvars::Array)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return get_bounds(domains,eqs,bcs,dict_indvars,dict_depvars)
end

function get_bounds(domains,bcs,_indvars::Array,_depvars::Array,strategy::QuadratureTraining)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
    return get_bounds(domains,bcs,dict_indvars,dict_depvars,strategy)
end

function get_bounds(domains,bcs,dict_indvars,dict_depvars,strategy::QuadratureTraining)
    bound_vars = get_variables(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_bounds(domains,eqs,bcs,dict_indvars,dict_depvars)
    dict_span = Dict([Symbol(d.variables) => [d.domain.lower, d.domain.upper] for d in domains])
    # pde_bounds = [[d.domain.lower,d.domain.upper] for d in domains]
    pde_args = get_argument(eqs,dict_indvars,dict_depvars)

    pde_bounds= map(pde_args) do pd
        span = map(p -> get(dict_span, p, p), pd)
    end

    bound_args = get_argument(bcs,dict_indvars,dict_depvars)
    dict_span = Dict([Symbol(d.variables) => [d.domain.lower, d.domain.upper] for d in domains])

    bcs_bounds= map(bound_args) do bt
        span = map(b -> get(dict_span, b, b), bt)
    end
    [pde_bounds,bcs_bounds]
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

function get_u()
    u = (cord, θ, phi)-> phi(cord, θ)
end

Base.Broadcast.broadcasted(::typeof(get_u()), cord, θ, phi) = get_u()(cord, θ, phi)

# the method to calculate the derivative
function get_numeric_derivative()
    _epsilon = 1 / (2*cbrt(eps(Float32)))
    derivative =
        (phi,u,x,εs,order,θ) ->
        begin
            ε = εs[order]
            ε = adapt(typeof(θ),ε)
            x = adapt(typeof(θ),x)
            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

Base.Broadcast.broadcasted(::typeof(get_numeric_derivative()), phi,u,x,εs,order,θ) = get_numeric_derivative()(phi,u,x,εs,order,θ)

function get_loss_function(loss_functions, train_sets, strategy::GridTraining;τ=nothing)
    if τ == nothing
        τs_ = [size(set)[2] for set in train_sets]
        τs = 1.0f0 ./ τs_
    end

    loss = (θ) ->  sum(τ * sum(abs2,loss_function(train_set, θ)) for (loss_function,train_set,τ) in zip(loss_functions,train_sets,τs))
    return loss
end

function generate_random_points(points, bound)
    function f(b)
      if b isa Number
           fill(Float32(b),(1,points))
       else
           lb, ub =  b[1], b[2]
           lb .+ ub .* Float32.(rand(1,points))
       end
    end
    vcat(f.(bound)...)
end

function get_loss_function(loss_functions, bounds, strategy::StochasticTraining;τ=nothing)
    points = strategy.points

    if τ == nothing
        τ = 1.0f0 / points
    end

    loss = (θ) -> begin
        total = 0.
        for (bound, loss_function) in zip(bounds, loss_functions)
            sets = generate_random_points(points, bound)
            sets_ = adapt(typeof(θ),sets)
            total += τ * sum(abs2,loss_function(sets_,θ))
        end
        return total
    end

    return loss
end

function get_loss_function(loss_functions, bounds, strategy::QuasiRandomTraining;τ=nothing)
    sampling_alg = strategy.sampling_alg
    points = strategy.points
    minibatch = strategy.minibatch

    if τ == nothing
        τ = 1.0f0 / points
    end

    sss = map(bounds) do bound
             map(bound) do b
                if !(b isa Number)
                    lb, ub =  [b[1]], [b[2]]
                    s = QuasiMonteCarlo.generate_design_matrices(points,lb,ub,sampling_alg,minibatch)
                else
                    s = fill(Float64(b),(1,points))
                end
             end
          end

    loss = (θ) -> begin
        total = 0.
        for (bound,ss_,loss_function) in zip(bounds,sss,loss_functions)
            ss__ =  [ss_[i] isa Array{Float64,2} ? ss_[i] : ss_[i][rand(1:minibatch)] for i in 1:length(ss_)]
            r_point = vcat(ss__...)
            r_point_ = adapt(typeof(θ),r_point)
            total += τ * sum(abs2,loss_function(r_point_,θ))
        end
        return total
    end
    return loss
end

function get_loss_function(loss_functions, bounds, strategy::QuadratureTraining;τ=nothing)
    lbs,ubs = bounds
    if τ == nothing
        τ = 1.0f0
    end

    f_ = (lb,ub,loss_,θ) -> begin
        function ___loss(x,θ)
            #x = adapt(typeof(θ),x)
            sum(abs2,loss_(x,θ), dims=2)
        end

        prob = QuadratureProblem(___loss,lb,ub,θ,batch = strategy.batch,nout=1)
        abs(solve(prob,
              strategy.quadrature_alg,
              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.eqs
    bcs = pde_system.bcs

    domains = pde_system.domain
    eq_params = pde_system.ps
    defaults = pde_system.defaults
    default_p = collect(values(defaults))

    param_estim = discretization.param_estim
    additional_loss = discretization.additional_loss

    # dimensionality of equation
    dim = length(domains)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)

    chain = discretization.chain
    initθ = discretization.init_params
    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy
    if !(eqs isa Array)
        eqs = [eqs]
    end
    symbolic_pde_loss_functions = [build_symbolic_loss_function(eq,indvars,depvars,
                                                              dict_indvars,dict_depvars,
                                                              phi, derivative,initθ,strategy;eq_params=eq_params,param_estim=param_estim,default_p=default_p) for eq in eqs]

    bc_indvars = if strategy isa QuadratureTraining
         get_argument(bcs,dict_indvars,dict_depvars)
    else
         get_variables(bcs,dict_indvars,dict_depvars)
    end
    symbolic_bc_loss_functions = [build_symbolic_loss_function(bc,indvars,depvars,
                                                               dict_indvars,dict_depvars,
                                                               phi, derivative,initθ,strategy,eq_params=eq_params,param_estim=param_estim,default_p=default_p;
                                                               bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
    symbolic_pde_loss_functions,symbolic_bc_loss_functions
end

# Convert a PDE problem into an OptimizationProblem
function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eqs
    bcs = pde_system.bcs

    domains = pde_system.domain
    eq_params = pde_system.ps
    defaults = pde_system.defaults
    default_p = collect(values(defaults))

    param_estim = discretization.param_estim
    additional_loss = discretization.additional_loss

    # dimensionality of equation
    dim = length(domains)
    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars,pde_system.depvars)

    chain = discretization.chain
    initθ = discretization.init_params
    flat_initθ = if length(depvars) != 1 vcat(initθ...) else  initθ end
    flat_initθ = if param_estim == false flat_initθ else vcat(flat_initθ, adapt(typeof(flat_initθ),default_p)) end

    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy
    if !(eqs isa Array)
        eqs = [eqs]
    end
    _pde_loss_functions = [build_loss_function(eq,indvars,depvars,
                                             dict_indvars,dict_depvars,
                                             phi, derivative, initθ,strategy,eq_params=eq_params,param_estim=param_estim,default_p=default_p) for eq in eqs]
    bc_indvars = if strategy isa QuadratureTraining
         get_argument(bcs,dict_indvars,dict_depvars)
    else
         get_variables(bcs,dict_indvars,dict_depvars)
    end

    _bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars,
                                                  phi, derivative, initθ, strategy,eq_params=eq_params,param_estim=param_estim,default_p=default_p;
                                                  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,eqs,bcs,
                                            dict_indvars,dict_depvars)

        # the points in the domain and on the boundary
        pde_train_sets, bcs_train_sets = train_sets

        pde_train_sets = adapt.(typeof(flat_initθ),pde_train_sets)
        bcs_train_sets =  adapt.(typeof(flat_initθ),bcs_train_sets)

        pde_loss_function = get_loss_function(_pde_loss_functions,
                                                        pde_train_sets,
                                                        strategy)

        bc_loss_function = get_loss_function(_bc_loss_functions,
                                                       bcs_train_sets,
                                                       strategy)
        (pde_loss_function, bc_loss_function)
    elseif strategy isa StochasticTraining
          bounds = get_bounds(domains,eqs,bcs,dict_indvars,dict_depvars)
          pde_bounds, bcs_bounds = bounds
          pde_loss_function = get_loss_function(_pde_loss_functions,
                                                          pde_bounds,
                                                          strategy)

          pde_dim = size(pde_bounds[1])[1] #TODO
          bcs_dim = isempty(maximum(size.(bcs_bounds[1]))) ? nothing : maximum(size.(bcs_bounds))[1]
          bcs_cond_size = size(bcs_bounds)[1]

          points = bcs_dim == nothing ? 1 : bcs_cond_size*Int(round(strategy.points^(bcs_dim/pde_dim)))
          strategy_ = StochasticTraining(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,eqs,bcs,dict_indvars,dict_depvars)
         pde_bounds, bcs_bounds = bounds
         pde_loss_function = get_loss_function(_pde_loss_functions,
                                                        pde_bounds,
                                                        strategy)

         pde_dim = size(pde_bounds[1])[1] #TODO
         bcs_dim = isempty(maximum(size.(bcs_bounds[1]))) ? nothing : maximum(size.(bcs_bounds))[1]
         bcs_cond_size = size(bcs_bounds)[1]

         points = bcs_dim == nothing ? 1 : bcs_cond_size*Int(round(strategy.points^(bcs_dim/pde_dim)))
         strategy_ = QuasiRandomTraining(points;
                                        sampling_alg = strategy.sampling_alg,
                                        minibatch = strategy.minibatch)

         bc_loss_function = get_loss_function(_bc_loss_functions,
                                                       bcs_bounds,
                                                       strategy_)
         (pde_loss_function, bc_loss_function)
    elseif strategy isa QuadratureTraining
        bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars,strategy)
        pde_bounds, bcs_bounds = bounds
        plbs,pubs = pde_bounds
        blbs,bubs = bcs_bounds
        pl = length(plbs[1])
        bl = length(blbs[1])
        bsl = length(blbs)

        τ_ = (10)^pl
        τp = 1.0f0 / τ_

        pde_loss_function = get_loss_function(_pde_loss_functions,
                                              pde_bounds,
                                              strategy;
                                              τ=τp)

        τb =  1.0f0 / (bsl * τ_^(bl/pl))

        if bl == 0
            _bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                          dict_indvars,dict_depvars,
                                                          phi, get_numeric_derivative(), initθ, GridTraining(0.1);
                                                          bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
            train_sets = generate_training_sets(domains,0.1,eqs,bcs,
                                                dict_indvars,dict_depvars)
            pde_train_set, bcs_train_sets = train_sets
            bc_loss_function = get_loss_function(_bc_loss_functions,
                                                 bcs_train_sets,
                                                 GridTraining(0.1))

        elseif bl == 1 && strategy.quadrature_alg in [CubaCuhre(),CubaDivonne()]
            @warn "$(strategy.quadrature_alg) does not work with one-dimensional
            problems, so for the boundary conditions loss function,
            the quadrature algorithm was replaced by HCubatureJL"

            strategy = QuadratureTraining(quadrature_alg = HCubatureJL(),
                                          reltol = bsl*(strategy.reltol)^(bl/pl),
                                          abstol = bsl*(strategy.abstol)^(bl/pl),
                                          maxiters = bsl*Int(round((strategy.maxiters)^(bl/pl))),
                                          batch = strategy.batch)
            bc_loss_function = get_loss_function(_bc_loss_functions,
                                                 bcs_bounds,
                                                 strategy;
                                                 τ=τb)
        else
            strategy = QuadratureTraining(quadrature_alg = strategy.quadrature_alg,
                                          reltol = bsl*(strategy.reltol)^(bl/pl),
                                          abstol = bsl*(strategy.abstol)^(bl/pl),
                                          maxiters = bsl*Int(round((strategy.maxiters)^(bl/pl))),
                                          batch = strategy.batch)
            bc_loss_function = get_loss_function(_bc_loss_functions,
                                             bcs_bounds,
                                             strategy;
                                             τ=τb)
        end
        (pde_loss_function, bc_loss_function)
    end

    function loss_function_(θ,p)
        if additional_loss isa Nothing
            return pde_loss_function(θ) + bc_loss_function(θ)
        else
            function _additional_loss(phi,θ)
                θ_ = θ[1:end - length(default_p)]
                p = θ[(end - length(default_p) + 1):end]
                return additional_loss(phi,θ_,p)
            end
            return pde_loss_function(θ) + bc_loss_function(θ) + _additional_loss(phi,θ)
        end
    end

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


text:
module mmmultest1
using Test
# using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using FinEtools.MatrixUtilityModule: mulCAB!

function gemmavx!(C, A, B)
	M, N = size(C); K = size(B,1)
	C .= 0
    @avx for n ∈ 1:N, k ∈ 1:K
    	Bkn = B[k,n]
        for m ∈ 1:M
            C[m,n] += A[m,k] * Bkn
        end
    end
    return C
end

function gemmblas!(C, A, B)
    mul!(C, A, B)
    return C
end

function test(n)
	# println("n = ", n)
	C, A, B = rand(n, n), rand(n, n), rand(n, n)
	@test norm(mulCAB!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	@test norm(gemmblas!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	@test norm(gemmavx!(C, A, B) - A * B) / norm(C) <= 1.0e-9

	BLAS.set_num_threads(1)
	# println("mulCAB!")

# @btime mulCAB!($C, $A, $B)
	# println("gemmblas!")

# @btime gemmblas!($C, $A, $B)
	# println("gemmavx!")

# @btime gemmavx!($C, $A, $B)
end
end
using .mmmultest1

mmmultest1.test(2^6)
mmmultest1.test(2^5)
mmmultest1.test(2^4)
mmmultest1.test(2^3)

mmmultest1.test(3)
mmmultest1.test(5)
mmmultest1.test(7)

module mmmultest2
using Test
# using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using FinEtools.MatrixUtilityModule: mulCAB!

function gemmavx!(C, A, B)
	M, N = size(C); K = size(B,1)
	C .= 0
    @avx  for n ∈ 1:N, k ∈ 1:K # for k ∈ 1:K, n ∈ 1:N #
    	Bkn = B[k,n]
        for m ∈ 1:M
            C[m,n] += A[m,k] * Bkn
        end
    end
    return C
end

function gemmblas!(C, A, B)
    mul!(C, A, B)
    return C
end

function test(n)
	# println("n = ", n)
	C, A, B = rand(n, n), rand(n, n), rand(n, n)
	@test norm(mulCAB!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	@test norm(gemmblas!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	@test norm(gemmavx!(C, A, B) - A * B) / norm(C) <= 1.0e-9

	BLAS.set_num_threads(1)
	# println("mulCAB!")

# @btime mulCAB!($C, $A, $B)
	# println("gemmblas!")

# @btime gemmblas!($C, $A, $B)
	# println("gemmavx!")

# @btime gemmavx!($C, $A, $B)
end
end
using .mmmultest2

mmmultest2.test(2^6)
mmmultest2.test(2^5)
mmmultest2.test(2^4)
mmmultest2.test(2^3)

module mmmultest3
using Test
# using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using FinEtools.MatrixUtilityModule: mulCAtB!

function gemmblas!(C, A, B)
    mul!(C, Transpose(A), B)
    return C
end

function test(n)
	# println("n = ", n)
	C, A, B = rand(n, n), rand(n, n), rand(n, n)
	mulCAtB!(C, A, B)
	@test norm(mulCAtB!(C, A, B) - A' * B) / norm(C) <= 1.0e-9
	@test norm(gemmblas!(C, A, B) - A' * B) / norm(C) <= 1.0e-9

	# println("mulCAtB!")

# @btime mulCAtB!($C, $A, $B)
	# println("gemmblas!")

# @btime gemmblas!($C, $A, $B)
end
end
using .mmmultest3
mmmultest3.test(2^6)
mmmultest3.test(2^5)
mmmultest3.test(2^4)
mmmultest3.test(2^3)

# mlatest1.test(3)
# mlatest1.test(5)
# mlatest1.test(7)

module mmmultest4
using Test
# using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using FinEtools.MatrixUtilityModule: mulCAtB!

function gemmblas!(C, A, B)
    mul!(C, Transpose(A), B)
    return C
end

function test(M, N, K)
	# println("M, N, K = ", (M, N, K))
	C, A, B = rand(M, N), rand(K, M), rand(K, N)
	mulCAtB!(C, A, B)
	@test norm(mulCAtB!(C, A, B) - A' * B) / norm(C) <= 1.0e-9
	@test norm(gemmblas!(C, A, B) - A' * B) / norm(C) <= 1.0e-9

	# println("mulCAtB!")

# @btime mulCAtB!($C, $A, $B)
	# println("gemmblas!")

# @btime gemmblas!($C, $A, $B)
end
end
using .mmmultest4
M, N, K = 2^6, 2^6, 2^6
mmmultest4.test(M, N, K)
M, N, K = 2^5, 2^5, 2^5
mmmultest4.test(M, N, K)
M, N, K = 2^4, 2^4, 2^4
mmmultest4.test(M, N, K)
M, N, K = 2^3, 2^3, 2^3
mmmultest4.test(M, N, K)

M, N, K = 3, 4, 5
mmmultest4.test(M, N, K)
M, N, K = 13, 4, 5
mmmultest4.test(M, N, K)
M, N, K = 6, 4, 9
mmmultest4.test(M, N, K)

module mmmultest6
using Test
# using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using FinEtools.MatrixUtilityModule: mulCABt!

function gemmblas!(C, A, B)
    mul!(C, A, Transpose(B))
    return C
end

function test(M, N, K)
	# println("M, N, K = ", (M, N, K))
	C, A, B = rand(M, N), rand(M, K), rand(N, K)
	@test norm(mulCABt!(C, A, B) - A * B') / norm(C) <= 1.0e-9
	@test norm(gemmblas!(C, A, B) - A * B') / norm(C) <= 1.0e-9

	# println("mulCABt!")

# @btime mulCABt!($C, $A, $B)
	# println("gemmblas!")

# @btime gemmblas!($C, $A, $B)
end
end
using .mmmultest6

M, N, K = 3, 4, 5
mmmultest6.test(M, N, K)
M, N, K = 13, 4, 5
mmmultest6.test(M, N, K)
M, N, K = 6, 4, 9
mmmultest6.test(M, N, K)
M, N, K = 7, 4, 9
mmmultest6.test(M, N, K)
M, N, K = 1, 4, 9
mmmultest6.test(M, N, K)
M, N, K = 7, 1, 9
mmmultest6.test(M, N, K)
M, N, K = 1, 1, 9
mmmultest6.test(M, N, K)
M, N, K = 1, 1, 1
mmmultest6.test(M, N, K)
M, N, K = 7, 4, 1
mmmultest6.test(M, N, K)

M, N, K = 2^6, 2^6, 2^6
mmmultest6.test(M, N, K)
M, N, K = 2^5, 2^5, 2^5
mmmultest6.test(M, N, K)
M, N, K = 2^4, 2^4, 2^4
mmmultest6.test(M, N, K)
M, N, K = 2^3, 2^3, 2^3
mmmultest6.test(M, N, K)

module mmvtest1
using FinEtools
using FinEtools.MatrixUtilityModule: mulCAB!
using LinearAlgebra
using Test
function test()
	M, N = 7, 9
	C, A, B = rand(M), rand(M, N), rand(N)
	@test norm(mulCAB!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	C, A, B = rand(M), rand(M, N), rand(N)
	@test norm(mulCAB!(C, A, B) - A * B) / norm(C) <= 1.0e-9
	C, A, B = rand(M), rand(M, N), rand(N)
	@test norm(mulCAB!(C, A, B) - A * B) / norm(C) <= 1.0e-9
end
end
using .mmvtest1
mmvtest1.test()

ext: jl


text:
using Statistics

function find_maxima(w, n)
    maxima = []
    for i in (1+n):(length(w)-n)
        if issorted(w[i-n:i], lt=isless) && issorted(w[i:i+n], lt=isless, rev = true)
            append!(maxima, i)
        end
        end
    return maxima
end

function get_regression_αβ(t, W, n)
    α = zeros(axes(W,2))
    β = zeros(axes(W,2))
    for i in axes(W,2)
        m = find_maxima(W[:,i], n)
        #!! bump on tail
        m = m[2:5]
        β[i] = cov(t[m], log.(W[m,i]), corrected=false) / var(t[m], corrected=false)
        α[i] = mean(log.(W[m,i])) - β[i]*mean(t[m])
    end
    return α, β
end


text:
"""
    Hopper
    	model inspired by "Dynamically Stable Legged Locomotion"
		s = (x, z, t, r)
			x - lateral position
			z - vertical position
			t - body orientation
			r - leg length
"""
struct Hopper{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    mb # mass of body
    ml # mass of leg
    Jb # inertia of body
    Jl # inertia of leg

    μ  # coefficient of friction
    g  # gravity

    qL::Vector
    qU::Vector

	uL::Vector
	uU::Vector

    nq
    nu
    nc
    nf
    nb
	ns

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 4 # configuration dimension
nu = 2 # control dimension
nc = 1 # number of contact points
nf = 2 # number of faces for friction cone
nb = nc * nf
ns = 1

# Parameters
gravity = 9.81 # gravity
μ = 1.0  # coefficient of friction
mb = 1.0 # body mass
ml = 0.1  # leg mass
Jb = 0.25 # body inertia
Jl = 0.025 # leg inertia

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = nq

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
kinematics(::Hopper, q) = [q[1] + q[4] * sin(q[3]), q[2] - q[4] * cos(q[3])]

# Methods
M_func(model::Hopper, q) = Diagonal(@SVector [
								 model.mb + model.ml,
								 model.mb + model.ml,
								 model.Jb + model.Jl,
								 model.ml])

C_func(model::Hopper, q, q̇) = @SVector [0.0,
									 (model.mb + model.ml) * model.g,
									 0.0,
									 0.0]

function ϕ_func(::Hopper, q)
    @SVector [q[2] - q[4] * cos(q[3])]
end

N_func(::Hopper, q) = @SMatrix [0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]))]

function _P_func(model, q)
	@SMatrix [1.0 0.0 (q[4] * cos(q[3])) sin(q[3])]
end

function P_func(::Hopper, q)
    @SMatrix [1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
    	-1.0 0.0 (-1.0 * q[4] * cos(q[3])) -1.0 * sin(q[3])]
end

B_func(::Hopper, q) = @SMatrix [0.0 0.0 1.0 0.0;
                                -sin(q[3]) cos(q[3]) 0.0 1.0]

function lagrangian_derivatives(model, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function fd(model::Hopper{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

	# evalutate at midpoint
	qm1 = 0.5 * (q1 + q2⁺)
    vm1 = (q2⁺ - q1) / h
    qm2 = 0.5 * (q2⁺ + q3)
    vm2 = (q3 - q2⁺) / h

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

    [q2⁺ - q2⁻;
	 (0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
    + (transpose(B_func(model, qm2)) * SVector{2}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{1}(λ)
    + transpose(P_func(model, q3)) * SVector{2}(b))
    + h * w)]
end

function maximum_dissipation(model::Hopper{Discrete, FixedTime}, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(model.nf)
	η = u[model.idx_η]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function fd(model::Hopper{Discrete, FreeTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	h = u[end]

	# evalutate at midpoint
	qm1 = 0.5 * (q1 + q2⁺)
    vm1 = (q2⁺ - q1) / h
    qm2 = 0.5 * (q2⁺ + q3)
    vm2 = (q3 - q2⁺) / h

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

    [q2⁺ - q2⁻;
	 (0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
	+ (transpose(B_func(model, qm2)) * SVector{2}(u_ctrl)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{2}(b))
	+ h * w)]
end

function maximum_dissipation(model::Hopper{Discrete, FreeTime}, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(model.nf)
	η = u[model.idx_η]
	h = u[end]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function no_slip(model::Hopper{Discrete, FixedTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	λ = view(u, model.idx_λ)
	s = view(u, model.idx_s)

	return s[1] - (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function no_slip(model::Hopper{Discrete, FreeTime}, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	λ = view(u, model.idx_λ)
	s = view(u, model.idx_s)
	h = u[end]
	return s[1] - (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function friction_cone(model::Hopper, u)
	λ = u[model.idx_λ]
	b = u[model.idx_b]
	return @SVector [model.μ * λ[1] - sum(b)]
end

r = 0.5
qL = -Inf * ones(nq)
qU = Inf * ones(nq)
qL[4] = 0.1
qU[4] = r

uL = -100.0 * ones(nu)
uU = 100.0 * ones(nu)

model = Hopper{Discrete, FixedTime}(n, m, d,
			   mb, ml, Jb, Jl,
			   μ, gravity,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

# Visualization
function visualize!(vis, model::Hopper, q;
		Δt = 0.1, scenario = :vertical)

    r_foot = 0.05
    r_leg = 0.5 * r_foot

	default_background!(vis)

    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))

    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics(model, q[t])[1], 0.0, kinematics(model, q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]

        z_shift = [0.0; 0.0; r_foot]

        MeshCat.atframe(anim, t) do
            settransform!(vis["body"], Translation(p_body + z_shift))
            settransform!(vis["foot"], Translation(p_foot + z_shift))

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

	if scenario == :vertical
		settransform!(vis["/Cameras/default"],
			compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))
	end

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using StatsPlots

# warm up (少し時間を取られる。その間に下のセルでコードを入力する)
histogram(randn(10^5); norm=true, alpha=0.3, size=(300, 200))
plot!(Normal())

# %%
# コインを20回投げたときの表が出る回数の分布は正規分布で近似される。

n = 20
p = 0.5
sample_of_means = rand(Binomial(n, p), 10^5)
histogram(sample_of_means; norm=true, alpha=0.3, bin=-0.5:n+0.5, label="rand(Binomial($n, $p))")
plot!(Normal(n*p, √(n*p*(1-p))), 0, 20; label="normal approx.", lw=2)

# %%
# 回数n、確率パラメータpの二項分布オブジェクトを次のように作れる。

bin = Binomial(n, p)

# %%
# 確率分布オブジェクトの乱数

rand(bin)

# %%
# 確率分布オブジェクトの乱数を複数個生成

@show rand(bin, 20);

# %%
# サンプルのヒストグラム
# norm = true で確率密度函数と比較できるようにする。
# alpha = 0.3 で半透明化pp
# bin = -0.5:n+0.5 で 1 刻みの適切なビンを設定

histogram(sample_of_means; norm=true, alpha=0.3, bin=-0.5:n+0.5, label="sample of means")

# %%
# -0.5:n+0.5 は [-0.5, 0.5, 1.5, 2.5, ..., 19.5, 20.5] とほぼ等価

-0.5:n+0.5

# %%
@show collect(-0.5:n+0.5);

# %%
# 平均μ、標準偏差σの正規分布オブジェクトは次にようにして作れる。

normal = Normal(n*p, √(n*p*(1-p)))

# %%
# 確率分布オブジェクト dist は plot(dist) でプロットできる。
# lw = 2 で2倍の太さでプロット

plot(normal, 0, 20; label="normal dist.", lw=2)

# %%
# 以上を重ねてプロット

histogram(sample_of_means; norm=true, alpha=0.3, bin=-0.5:n+0.5, label="rand(Binomial($n, $p))")
plot!(Normal(n*p, √(n*p*(1-p))), 0, 20; label="normal approx.", lw=2)

# %%
# サイコロを10回振って出た目の平均値の分布は正規分布でよく近似される。

n = 10
M = [mean(rand(1:6, n)) for _ in 1:10^5]
histogram(M; norm=true, alpha=0.3, bin=0.95:0.1:6.05, label="mean(rand(1:6, 10))")
plot!(Normal(mean(1:6), std(1:6; corrected=false)/√n), 1, 6; label="normal approx.", lw=2)

# %%
# 2つの確率パラメータ(比率パラメータ)が等しい二項分布で生成されたサンプルのZ統計量は
# 近似的に標準正規分布に従う。

""""""
function gensample(bin1, bin2)
    m, n = ntrials(bin1), ntrials(bin2)
    a, b = rand(bin1), rand(bin2)
    c, d = m - a, n - b
    a, b, c, d
end

"""Z statistic of 2×2 contingency table"""
function zstat(a, b, c, d)
    a*d - b*c == 0 && return zero(inv(a))
    m, n = a + c, b + d
    N = a + b + c + d
    p = (a + b)/N
    Z = (a/m - b/n) / √(p*(1-p)*(1/m + 1/n))
    Z
end

m, n, p = 30, 20, 0.3
bin1 = Binomial(m, p)
bin2 = Binomial(n, p)
Z = [zstat(gensample(bin1, bin2)...) for _ in 1:10^5]

histogram(Z; norm=true, alpha=0.3, bin=-4.25:0.5:4.25, label="Z statistics")
plot!(Normal(), -4, 4; label="std normal dist.", lw=2)

# %%
@code_warntype zstat(1, 2, 3, 4)

# %%
# 数式処理を使って、Z統計量の2乗がPearsonのχ²統計量に一致することを確認。
# 上の数値計算で使ったzstat函数をそのまま数式処理で利用する。

using SymPy

@syms a b c d
Z = zstat(a, b, c, d)

# %%
# Z^2 |> factor は factor(Z^2) に等価
# 以下の計算結果はPearsonのχ²統計量の有名な表示に一致

Z^2 |> factor

# %%
# (O - E)²/E の和で定義されたPearsonのχ²統計量

"""Pearson's χ² statistic of 2×2 contingency table"""
function chisqstat(a, b, c, d)
    a*d - b*c == 0 && return zero(inv(a))
    N = a + b + c + d
    s, f = a + b, c + d
    m, n = a + c, b + d
    Ea, Eb, Ec, Ed = m*s/N, n*s/N, m*f/N, n*f/N
    X² = (a - Ea)^2/Ea + (b - Eb)^2/Eb + (c - Ec)^2/Ec + (d - Ed)^2/Ed
end

X² = chisqstat(a, b, c, d)

# %%
# 以下の計算結果は2×2の分割表のPearson χ²統計量に関する有名な公式

X² |> factor

# %%
# 正確二項検定のP値の定義

# \lessapprox TAB → ⪅
# \approx TAB → ≈
x ⪅ y = x < y || x ≈ y

"""p-value function of exact binomial test"""
pval_exact(dist, k) = sum(pdf(dist, j) for j in support(dist) if pdf(dist, j) ⪅ pdf(dist, k))

# %%
# 正確二項検定のP値の例　(パラメータpを固定した場合(1))

n, p = 20, 0.3
bin = Binomial(n, p)

k = support(bin)
y = pval_exact.(bin, k)
plot(k, y; label="", marker=:o)
plot!(; xtick=0:20, ytick=0:0.05:1)
plot!(; xlabel="data k", ylabel="p-value for parameter n = $n, p = $p")

# %%
# 正確二項検定のP値の例　(データkを固定した場合)

n, k = 20, 6

p = 0:0.002:1
y = @. pval_exact(Binomial(n, p), k)
plot(p, y; label="")
plot!(; xtick=0:0.1:1, ytick=0:0.05:1)
plot!(; xlabel="parameter p", ylabel="p-value for data n = $n, k = $k")

# %%
# 信頼区間函数の定義

using Roots

"""generic confidence interval function"""
function confint(pvalfunc, mlefunc, n, k, pmin, pmax; α = 0.05)
    f(p) = pvalfunc(n, p, k) - α
    ci = find_zeros(f, pmin, pmax)
    length(ci) ≥ 2 && return ci
    pmle = mlefunc(n, k)
    first(ci) > pmle ? [pmin, first(ci)] : [first(ci), pmax]
end

"""confidence interval function for exact binomial test"""
function confint_bin(n, k; α = 0.05)
    pvalfunc(n, p, k) = pval_exact(Binomial(n, p), k)
    mlefunc(n, k) = k/n
    pmin, pmax = 0.0, 1.0
    confint(pvalfunc, mlefunc, n, k, pmin, pmax; α)
end

# %%
# 正確二項検定のデータkを固定した場合のP値函数のプロットに
# 95%信頼区間のプロットを追加

n, k = 20, 6
α = 0.05
ci = confint_bin(n, k; α)

p = 0:0.002:1
y = @. pval_exact(Binomial(n, p), k)
plot(p, y; label="")
plot!(ci, [α, α]; label="$(100(1 - α))% confidence interval", lw=3)
plot!(; xtick=0:0.1:1, ytick=0:0.05:1)
plot!(; xlabel="parameter p", ylabel="p-value for data n = $n, k = $k")

# %% tags=[]
# データkを動かして信頼区間をプロット (1)

n = 20
k = 0:n
α = 0.05
ci = confint_bin.(n, k; α)
ciL, ciR = first.(ci), last.(ci)

plot(; legend=:topleft)
plot!(k, ciL; label="min. of confidence interval")
plot!(k, ciR; label="max. of confidenceinterval")
plot!(k, ciL; label="$(100(1 - α))% confidence interval (n = $n)", lw=0, frange=ciR, fa=0.1)
plot!(; xtick=0:n, ytick=0:0.1:1)
plot!(; xlabel="data k", ylabel="parameter p")

# %% tags=[]
# データkを動かして信頼区間をプロット (2)

n = 20
k = 0:n
α = 0.05
ci = confint_bin.(n, k; α)

plot(; legend=:topleft)
for i in eachindex(ci)
    plot!([k[i], k[i]], ci[i]; label="", lw=5)
end
plot!(; xtick=0:n, ytick=0:0.1:1)
plot!(; xlabel="data k", ylabel="parameter p")
title!("$(100(1 - α))% confidence intervals (n = $n)"; titlefontsize=12)

# %% tags=[]
# データkを動かして信頼区間をプロット (3)

n = 20
k = 0:n
α = 0.05
ci = confint_bin.(n, k; α)

plot(; legend=:topleft)
for i in eachindex(ci)
    plot!([k[i], k[i]], ci[i]; label="", lw=5, c=1)
end
plot!(; xtick=0:n, ytick=0:0.1:1)
plot!(; xlabel="data k", ylabel="parameter p")
title!("$(100(1 - α))% confidence intervals (n = $n)"; titlefontsize=12)

# %%
# P値函数のヒートマップ

n = 20
p = 0:0.002:1
k = 0:n
α = 0.05
ci = confint_bin.(n, k; α)

z = @. pval_exact(Binomial(n, p), k')
heatmap(k, p, z; clim=(-0.2, 1.2), colorbar_title="p-value")
plot!(; xtick=0:n, ytick=0:0.1:1)
plot!(; xlabel="data k", ylabel="parameter p")
plot!(; size=(680, 400))

# %%
# P値函数のヒートマップと信頼区間

n = 20
p = 0:0.002:1
k = 0:n
α = 0.05
ci = confint_bin.(n, k; α)

z = @. pval_exact(Binomial(n, p), k')
heatmap(k, p, z; clim=(-0.2, 1.2), colorbar_title="p-value")
plot!(; xtick=0:n, ytick=0:0.1:1)
plot!(; xlabel="data k", ylabel="parameter p")
plot!(; size=(680, 400))
for i in eachindex(ci)
    plot!([k[i], k[i]], ci[i]; label="", lw=3, c=:cyan, alpha=0.8)
end
plot!(xlim=(-0.5, n+0.5), ylim=(0, 1))
title!("$(100(1 - α))% confidence intervals (n = $n)"; titlefontsize=12)

# %%

max_stars_repo_path: 0024/JuliaLang examples of probability and statistics.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@testset "check parameter order change" begin
    R, (x, y, z) = PolynomialRing(QQ, ["x","y","z"], ordering=:deglex)

    fs = [y^2, x]
    x, y, z = gens(parent(first(Groebner.groebner(fs, ordering=:lex))))
    @test Groebner.groebner(fs, ordering=:lex) == [y^2, x]
    @test ordering(parent(first(Groebner.groebner(fs, ordering=:lex)))) === :lex

    x, y, z = gens(parent(first(Groebner.groebner(fs, ordering=:degrevlex))))
    @test Groebner.groebner(fs, ordering=:degrevlex) == [x, y^2]
    @test ordering(parent(first(Groebner.groebner(fs, ordering=:degrevlex)))) === :degrevlex

    x, y, z = gens(parent(first(fs)))
    @test Groebner.groebner(fs, ordering=:deglex) == [x, y^2]
    @test ordering(parent(first(Groebner.groebner(fs, ordering=:deglex)))) === :deglex

    x, y, z = gens(parent(first(Groebner.groebner(fs, forsolve=true))))
    @test Groebner.groebner(fs, forsolve=true) == [y^2, x]
    @test ordering(parent(first(Groebner.groebner(fs, forsolve=true)))) === :lex

    noon = Groebner.change_ordering(Groebner.noonn(2), :lex)
    gb = Groebner.groebner(noon, ordering=:lex)
    x1, x2 = gens(parent(first(gb)))
    @test gb == [
            x2^5 - 10//11*x2^4 - 11//5*x2^3 + 2*x2^2 + 331//1100*x2 - 11//10,
            x1 + x2^4 - 10//11*x2^3 - 11//10*x2^2 + x2 - 10//11]

    gb = Groebner.groebner(noon, ordering=:degrevlex)
    x1, x2 = gens(parent(first(gb)))
    @test gb == [
            x1^2 - x2^2 - 10//11*x1 + 10//11*x2,
            x2^3 + 10//11*x1*x2 - 10//11*x2^2 - 11//10*x2 + 1,
            x1*x2^2 - 11//10*x1 + 1]

    gb = Groebner.groebner(noon, ordering=:deglex)
    x1, x2 = gens(parent(first(gb)))
    @test gb == [
            x1^2 - x2^2 - 10//11*x1 + 10//11*x2,
            x2^3 + 10//11*x1*x2 - 10//11*x2^2 - 11//10*x2 + 1,
            x1*x2^2 - 11//10*x1 + 1]

end


text:
using Paxos.Ballots
using Paxos.Configurations
using Paxos.Nodes

@testset CustomTestSet "Configurations" begin

@test isempty(configuration())

@test quorumSize(configuration()) == 1

@test length(addMember(configuration(), NodeID(), "test1").new.members) == 1

@test quorumSize(addMember(configuration(), NodeID(), "test1").new) == 1

@test quorumSize(
  addMember(
    addMember(
      configuration(),
      NodeID(),
      "test1"),
    NodeID(),
    "test2"
    ).new) == 2

@test quorumSize(3) == 2

end


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "MHD" begin

  @trixi_testset "elixir_mhd_alfven_wave.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.440611823425164e-15, 1.1373567770134494e-14, 3.024482376149653e-15, 2.0553143516814395e-15, 3.9938347410210535e-14, 3.984545392098788e-16, 2.4782402104201577e-15, 1.551737464879987e-15],
      linf = [1.9984014443252818e-15, 1.3405943022348765e-14, 3.3584246494910985e-15, 3.164135620181696e-15, 7.815970093361102e-14, 8.881784197001252e-16, 2.886579864025407e-15, 2.942091015256665e-15],
      initial_condition = initial_condition_constant,
      tspan = (0.0,1.0))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.0375628983659061e-5, 6.571144191446236e-7, 3.5833569836289104e-5, 3.583356983615859e-5, 5.084863194951084e-6, 1.1963224165731992e-16, 3.598916927583752e-5, 3.598916927594727e-5],
      linf = [2.614095879338585e-5, 9.577266731216823e-7, 0.00012406198007461344, 0.00012406198007509917, 1.5066209528846741e-5, 2.220446049250313e-16, 0.00012658678753942054, 0.00012658678753908748])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.4396053943470756e-5, 1.1211016739165248e-5, 3.577870687983967e-5, 3.577870687982181e-5, 1.967962220860377e-6, 1.1963224165731992e-16, 3.583562899483433e-5, 3.583562899486565e-5],
      linf = [5.830577969345718e-5, 3.280495696370357e-5, 0.00012279619948236953, 0.00012279619948227238, 6.978806516122482e-6, 2.220446049250313e-16, 0.00012564003648959932, 0.00012564003648994626],
      volume_flux = flux_derigs_etal)
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [5.86009540e-02, 8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02],
      linf = [1.10014999e-01, 1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02])
  end

  @trixi_testset "elixir_mhd_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"),
      l2   = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0],
      linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0])
  end

  @trixi_testset "elixir_mhd_torrilhon_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
      l2   = [0.45700904847931145, 0.4792535936512035, 0.340651203521865, 0.4478034694296928, 0.9204708961093411, 1.3216517820475193e-16, 0.28897419402047725, 0.25521206483145126],
      linf = [1.2185238171352286, 0.8913202384963431, 0.8488793580488431, 0.973083603686, 1.660723397705417, 2.220446049250313e-16, 0.6874726847741993, 0.65536978110274])
  end

  @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"),
      l2   = [2.34809441e-01, 3.92255943e-01, 8.23575546e-02, 1.75599624e-01, 9.61613519e-01, 6.60825891e-17, 2.15346454e-01, 1.07006529e-01],
      linf = [6.40732148e-01, 9.44889516e-01, 3.54932707e-01, 8.54060243e-01, 2.07757711e+00, 1.11022302e-16, 4.92584725e-01, 2.49526561e-01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_mhd_shu_osher_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.01126210e+00, 8.27157902e+00, 1.30882545e+00, 0.00000000e+00, 5.21930435e+01, 6.56538824e-16, 1.01022340e+00, 0.00000000e+00],
      linf = [2.87172004e+00, 2.26438057e+01, 4.16672442e+00, 0.00000000e+00, 1.35152372e+02, 3.44169138e-15, 2.83556069e+00, 0.00000000e+00],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_mhd_shu_osher_shock_tube.jl with flipped shock direction" begin
    # Include this elixir to make `initial_condition_shu_osher_shock_tube_flipped` available, which is used below
    trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"), tspan=(0.0, 0.0))
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.01539817e+00, 8.29625810e+00, 1.29548008e+00, 0.00000000e+00, 5.23565514e+01, 3.18641825e-16, 1.00485291e+00, 0.00000000e+00],
      linf = [2.92876280e+00, 2.28341581e+01, 4.11643561e+00, 0.00000000e+00, 1.36966213e+02, 1.55431223e-15, 2.80548864e+00, 0.00000000e+00],
      initial_condition = initial_condition_shu_osher_shock_tube_flipped,
      boundary_conditions=BoundaryConditionDirichlet(initial_condition_shu_osher_shock_tube_flipped),
      tspan = (0.0, 0.2))
  end
end

end # module


text:
"""
    GlmResp

The response vector and various derived vectors in a generalized linear model.
"""
struct GlmResp{V<:FPVector,D<:UnivariateDistribution,L<:Link} <: ModResp
    "`y`: response vector"
    y::V
    d::D
    "`devresid`: the squared deviance residuals"
    devresid::V
    "`eta`: the linear predictor"
    eta::V
    "`mu`: mean response"
    mu::V
    "`offset:` offset added to `Xβ` to form `eta`.  Can be of length 0"
    offset::V
    "`wts:` prior case weights.  Can be of length 0."
    wts::V
    "`wrkwt`: working case weights for the Iteratively Reweighted Least Squares (IRLS) algorithm"
    wrkwt::V
    "`wrkresid`: working residuals for IRLS"
    wrkresid::V
end

function GlmResp(y::V, d::D, l::L, η::V, μ::V, off::V, wts::V) where {V<:FPVector, D, L}
    if d == Binomial()
        for yy in y
            0 ≤ yy ≤ 1 || throw(ArgumentError("$yy in y is not in [0,1]"))
        end
    else
        all(x -> insupport(d, x), y) || throw(ArgumentError("y must be in the support of D"))
    end
    n = length(y)
    nη = length(η)
    nμ = length(μ)
    length(wts) == nη == nμ == n || throw(DimensionMismatch(
        "lengths of η, μ, y and wts ($nη, $nμ, $(length(wts)), $n) are not equal"))
    lo = length(off)
    lo == 0 || lo == n || error("offset must have length $n or length 0")
    res = GlmResp{V,D,L}(y, d, similar(y), η, μ, off, wts, similar(y), similar(y))
    updateμ!(res, η)
    res
end

deviance(r::GlmResp) = sum(r.devresid)

"""
    cancancel(r::GlmResp{V,D,L})

Returns `true` if dμ/dη for link `L` is the variance function for distribution `D`

When `L` is the canonical link for `D` the derivative of the inverse link is a multiple
of the variance function for `D`.  If they are the same a numerator and denominator term in
the expression for the working weights will cancel.
"""
cancancel(::GlmResp) = false
cancancel(::GlmResp{V,D,LogitLink}) where {V,D<:Union{Bernoulli,Binomial}} = true
cancancel(::GlmResp{V,D,IdentityLink}) where {V,D<:Normal} = true
cancancel(::GlmResp{V,D,LogLink}) where {V,D<:Poisson} = true

"""
    updateμ!{T<:FPVector}(r::GlmResp{T}, linPr::T)

Update the mean, working weights and working residuals, in `r` given a value of
the linear predictor, `linPr`.
"""
function updateμ!(r::GlmResp{T}, linPr::T) where T<:FPVector
    isempty(r.offset) ? copy!(r.eta, linPr) : broadcast!(+, r.eta, linPr, r.offset)
    updateμ!(r)
    if !isempty(r.wts)
        map!(*, r.devresid, r.devresid, r.wts)
        map!(*, r.wrkwt, r.wrkwt, r.wts)
    end
    r
end

function updateμ!(r::GlmResp{V,D,L}) where {V<:FPVector,D,L}
    y, η, μ, wrkres, wrkwt, dres = r.y, r.eta, r.mu, r.wrkresid, r.wrkwt, r.devresid

    @inbounds for i in eachindex(y, η, μ, wrkres, wrkwt, dres)
        μi, dμdη = inverselink(L(), η[i])
        μ[i] = μi
        yi = y[i]
        wrkres[i] = (yi - μi) / dμdη
        wrkwt[i] = cancancel(r) ? dμdη : abs2(dμdη) / glmvar(r.d, μi)
        dres[i] = devresid(r.d, yi, μi)
    end
end

function updateμ!(r::GlmResp{V,D,L}) where {V<:FPVector,D<:Union{Bernoulli,Binomial},L<:Link01}
    y, η, μ, wrkres, wrkwt, dres = r.y, r.eta, r.mu, r.wrkresid, r.wrkwt, r.devresid

    @inbounds for i in eachindex(y, η, μ, wrkres, wrkwt, dres)
        μi, dμdη, μomμ = inverselink(L(), η[i])
        μ[i] = μi
        yi = y[i]
        wrkres[i] = (yi - μi) / dμdη
        wrkwt[i] = cancancel(r) ? dμdη : abs2(dμdη) / μomμ
        dres[i] = devresid(r.d, yi, μi)
    end
end

"""
    wrkresp(r::GlmResp)

The working response, `r.eta + r.wrkresid - r.offset`.
"""
wrkresp(r::GlmResp) = wrkresp!(similar(r.eta), r)

"""
    wrkresp!{T<:FPVector}(v::T, r::GlmResp{T})

Overwrite `v` with the working response of `r`
"""
function wrkresp!(v::T, r::GlmResp{T}) where T<:FPVector
    broadcast!(+, v, r.eta, r.wrkresid)
    isempty(r.offset) ? v : broadcast!(-, v, v, r.offset)
end

abstract type AbstractGLM <: LinPredModel end

mutable struct GeneralizedLinearModel{G<:GlmResp,L<:LinPred} <: AbstractGLM
    rr::G
    pp::L
    fit::Bool
end

function coeftable(mm::AbstractGLM)
    cc = coef(mm)
    se = stderr(mm)
    zz = cc ./ se
    CoefTable(hcat(cc,se,zz,2.0 * ccdf.(Normal(), abs.(zz))),
              ["Estimate","Std.Error","z value", "Pr(>|z|)"],
              ["x$i" for i = 1:size(mm.pp.X, 2)], 4)
end

function confint(obj::AbstractGLM, level::Real)
    hcat(coef(obj),coef(obj)) + stderr(obj)*quantile(Normal(),(1. -level)/2.)*[1. -1.]
end
confint(obj::AbstractGLM) = confint(obj, 0.95)

deviance(m::AbstractGLM) = deviance(m.rr)

function loglikelihood(m::AbstractGLM)
    r = m.rr
    wts = r.wts
    y = r.y
    mu = r.mu
    ϕ = deviance(m)/sum(wts)
    d = r.d
    ll = zero(loglik_obs(d, y[1], mu[1], wts[1], ϕ))
    @inbounds for i in eachindex(y, mu, wts)
        ll += loglik_obs(d, y[i], mu[i], wts[i], ϕ)
    end
    ll
end

dof(x::GeneralizedLinearModel) = dispersion_parameter(x.rr.d) ? length(coef(x)) + 1 : length(coef(x))

function _fit!(m::AbstractGLM, verbose::Bool, maxIter::Integer, minStepFac::Real,
              convTol::Real, start)
    m.fit && return m
    maxIter >= 1 || throw(ArgumentError("maxIter must be positive"))
    0 < minStepFac < 1 || throw(ArgumentError("minStepFac must be in (0, 1)"))

    cvg, p, r = false, m.pp, m.rr
    lp = r.mu
    if start == nothing || isempty(start)
        delbeta!(p, wrkresp(r), r.wrkwt)
        linpred!(lp, p)
        updateμ!(r, lp)
        installbeta!(p)
    else
        copy!(p.beta0, start)
        fill!(p.delbeta, 0)
        linpred!(lp, p, 0)
        updateμ!(r, lp)
    end
    devold = deviance(m)
    for i = 1:maxIter
        f = 1.0
        local dev
        try
            delbeta!(p, r.wrkresid, r.wrkwt)
            linpred!(lp, p)
            updateμ!(r, lp)
            dev = deviance(m)
        catch e
            isa(e, DomainError) ? (dev = Inf) : rethrow(e)
        end
        while dev > devold
            f /= 2.
            f > minStepFac || error("step-halving failed at beta0 = $(p.beta0)")
            try
                updateμ!(r, linpred(p, f))
                dev = deviance(m)
            catch e
                isa(e, DomainError) ? (dev = Inf) : rethrow(e)
            end
        end
        installbeta!(p, f)
        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol || dev == 0
            cvg = true
            break
        end
        @assert isfinite(crit)
        devold = dev
    end
    cvg || throw(ConvergenceException(maxIter))
    m.fit = true
    m
end

StatsBase.fit!(m::AbstractGLM; verbose::Bool=false, maxIter::Integer=30,
              minStepFac::Real=0.001, convTol::Real=1.e-6, start=nothing) =
    _fit!(m, verbose, maxIter, minStepFac, convTol, start)

function initialeta!(dist::UnivariateDistribution, link::Link,
                     eta::AbstractVector, y::AbstractVector, wts::AbstractVector,
                     off::AbstractVector)
    length(eta) == length(y) == length(wts) || throw(DimensionMismatch("argument lengths do not match"))
    @inbounds @simd for i = eachindex(y, eta, wts)
        μ = mustart(dist, y[i], wts[i])
        eta[i] = linkfun(link, μ)
    end
    if !isempty(off)
        @inbounds @simd for i = eachindex(eta,off)
            eta[i] -= off[i]
        end
    end
    eta
end

function StatsBase.fit!(m::AbstractGLM, y; wts=nothing, offset=nothing, dofit::Bool=true,
                        verbose::Bool=false, maxIter::Integer=30, minStepFac::Real=0.001, convTol::Real=1.e-6,
                        start=nothing)
    r = m.rr
    V = typeof(r.y)
    r.y = copy!(r.y, y)
    isa(wts, Void) || copy!(r.wts, wts)
    isa(offset, Void) || copy!(r.offset, offset)
    initialeta!(r.d, r.l, r.eta, r.y, r.wts, r.offset)
    updateμ!(r, r.eta)
    fill!(m.pp.beta0, 0)
    m.fit = false
    if dofit
        _fit!(m, verbose, maxIter, minStepFac, convTol, start)
    else
        m
    end
end

function fit(::Type{M},
    X::Union{Matrix{T},SparseMatrixCSC{T}}, y::V,
    d::UnivariateDistribution,
    l::Link = canonicallink(d);
    dofit::Bool = true,
    wts::V = ones(y),
    offset::V = similar(y, 0), fitargs...) where {M<:AbstractGLM,T<:FP,V<:FPVector}

    size(X, 1) == size(y, 1) || throw(DimensionMismatch("number of rows in X and y must match"))
    n = length(y)
    length(wts) == n || throw(DimensionMismatch("length(wts) does not match length(y)"))
    if length(offset) != n && length(offset) != 0
        throw(DimensionMismatch("length(offset) does not match length(y)"))
    end

    wts = T <: Float64 ? copy(wts) : convert(typeof(y), wts)
    off = T <: Float64 ? copy(offset) : convert(Vector{T}, offset)
    eta = initialeta!(d, l, similar(y), y, wts, off)
    rr = GlmResp(y, d, l, eta, similar(y), offset, wts)
    res = M(rr, cholpred(X), false)
    dofit ? fit!(res; fitargs...) : res
end

fit(::Type{M},
X::Union{Matrix,SparseMatrixCSC},
y::AbstractVector,
d::UnivariateDistribution,
l::Link=canonicallink(d); kwargs...) where {M<:AbstractGLM} =
    fit(M, float(X), float(y), d, l; kwargs...)

glm(X, y, args...; kwargs...) = fit(GeneralizedLinearModel, X, y, args...; kwargs...)

GLM.Link(mm::AbstractGLM) = mm.l
GLM.Link(r::GlmResp{T,D,L}) where {T,D,L} = L()
GLM.Link(m::GeneralizedLinearModel) = Link(m.rr)

Distributions.Distribution(r::GlmResp{T,D,L}) where {T,D,L} = D
Distributions.Distribution(m::GeneralizedLinearModel) = Distribution(m.rr)

"""
    dispersion(m::AbstractGLM, sqr::Bool=false)

Return the estimated dispersion (or scale) parameter for a model's distribution,
generally written σ² for linear models and ϕ for generalized linear models.
It is, by definition, equal to 1 for the Bernoulli, Binomial, and Poisson families.

If `sqr` is `true`, the squared dispersion parameter is returned.
"""
function dispersion(m::AbstractGLM, sqr::Bool=false)
    r = m.rr
    if dispersion_parameter(r.d)
        wrkwt, wrkresid = r.wrkwt, r.wrkresid
        s = sum(i -> wrkwt[i] * abs2(wrkresid[i]), eachindex(wrkwt, wrkresid)) / dof_residual(m)
        sqr ? s : sqrt(s)
    else
        one(eltype(r.mu))
    end
end

"""
    predict(mm::AbstractGLM, newX::AbstractMatrix; offset::FPVector=Vector{eltype(newX)}(0))

Form the predicted response of model `mm` from covariate values `newX` and, optionally,
an offset.
"""
function predict(mm::AbstractGLM, newX::AbstractMatrix;
                 offset::FPVector=eltype(newX)[])
    eta = newX * coef(mm)
    if !isempty(mm.rr.offset)
        length(offset) == size(newX, 1) ||
            throw(ArgumentError("fit with offset, so `offset` kw arg must be an offset of length `size(newX, 1)`"))
        broadcast!(+, eta, eta, offset)
    else
        length(offset) > 0 && throw(ArgumentError("fit without offset, so value of `offset` kw arg does not make sense"))
    end
    mu = [linkinv(Link(mm), x) for x in eta]
end


text:
##########################################################################
#
#    Model Expression parsing
#      - transforms MCMC specific idioms (~) into regular Julia syntax
#      - calls the ReverseDiffSource module for gradient code generation
#      - creates function
#
##########################################################################

# naming conventions
const ACC_SYM = :__acc       # name of accumulator variable
const PARAM_SYM = :__beta    # name of parameter vector

#######################################################################
#   generates the log-likelihood function
#######################################################################
# - 'init' contains the dictionary of model params and their initial value
# - If 'debug' is set to true, the function returns only the function expression
#  that would have been created

function generateModelFunction(model::Expr; gradient=false, debug=false, init...)

	model.head != :block && (model = Expr(:block, model))  # enclose in block if needed
	length(model.args)==0 && error("model should have at least 1 statement")

	vsize, pmap, vinit = modelVars(;init...) # model param info

	model = translate(model) # rewrite ~ statements
	rv = symbol("$(ACC_SYM)v")  # final result in this variable
	model = Expr(:block, [ :($ACC_SYM = LLAcc(0.)), # add log-lik accumulator initialization
		                   model.args,
		                   # :( $ACC_SYM = $(Expr(:., ACC_SYM, Expr(:quote, :val)) ) )]... )
		                   :( $rv = $(Expr(:., ACC_SYM, Expr(:quote, :val)) ) )]... )

	## build function expression
	if gradient  # case with gradient
		head, body, outsym = ReverseDiffSource.reversediff(model,
			                                               rv, false, MCMC;
			                                               init...)

		body = [ vec2var(;init...),  # assigments beta vector -> model parameter vars
		         body.args,
		         :(($outsym, $(var2vec(;init...))))]

		# enclose in a try block
		body = Expr(:try, Expr(:block, body...),
				          :e,
				          quote
				          	if isa(e, OutOfSupportError)
				          		return(-Inf, zero($PARAM_SYM))
				          	else
				          		rethrow(e)
				          	end
				          end)

	else  # case without gradient
		head, body, outsym = ReverseDiffSource.reversediff(model,
			                                               rv, true, MCMC;
			                                               init...)

		body = [ vec2var(;init...),  # assigments beta vector -> model parameter vars
		         body.args,
		         outsym ]

		# enclose in a try block
		body = Expr(:try, Expr(:block, body...),
				          :e,
				          quote
				          	if isa(e, OutOfSupportError)
				          		return(-Inf)
				          	else
				          		rethrow(e)
				          	end
				          end)

	end

	# build and evaluate the let block containing the function and var declarations
	fn = gensym("ll")
	body = Expr(:function, Expr(:call, fn, :($PARAM_SYM::Vector{Float64})),	Expr(:block, body) )
	body = Expr(:let, Expr(:block, :(global $fn), head.args..., body))

	# println("#############\n$body\n############")

	debug ? body : (eval(body) ; (eval(fn), vsize, pmap, vinit) )
end


text:
## Documentation not attached to an Symata Head with apprules

# FIXME: All means other things as well
@sjdoc All """
   All

a symbol used in options.
"""

@sjdoc Regex """
    r\"expr\"

create a regular expression (PCRE). This is a Julia `DataType`.
"""

@sjdoc macros """
macros

If Symata encounters a macro call in input, it first Julia-evaluates all the
arguments then Julia-evaluates the macro and inserts it into the Symata expression
tree. For instance, big numbers and regular expressions are constructed this way.
"""

@sjdoc ans """
    ans

holds the result of the most recent evaluation.
"""

@sjdoc Sequence """
    Sequence

represents a sequence of arguments to be spliced into function arguments.
"""

@sjexamp(Sequence,
         ("f(a,b,Sequence(c,d),e,g)", "f(a,b,c,d,e,g)"))

@sjdoc bf """
    bf"str"

convert string `str` to a `BigFloat`.
"""

@sjdoc BF """
 BF(str)

converts `str` to type `BigFloat`.
"""

@sjdoc bi """
    bi"str"

convert string `str` to type `BigInt`.
"""

@sjdoc BI "
    bI(n)

converts `n` type `BigInt`.
"

@sjdoc Unicode_Input """
Unicode characters may be input at the Symata REPL just as they are in Julia. Eg. \\alpha [TAB] will enter the greek letter α.
"""

@sjseealso_group(ReturnSymPy,ToSymPy,ToSymata,ShowSymPyDocs)

@sjdoc Reals """
    Reals

represents the real numbers.
"""

@sjdoc Integers """
    Integers

represents the integer numbers.
"""

@sjdoc Complexes """
    Complexes

represents the complex numbers.
"""

@sjdoc Rationals """
    Rationals

represents the rational numbers.
"""

@sjdoc Primes """
    Primes

represents the prime numbers.
"""

@sjdoc Booleans """
    Booleans

represents the set of boolean values, `True` and `False`.
"""

@sjdoc Algebraics """
    Algebraics

represents the algebraic numbers.
"""

@sjdoc Repeated """
    Repeated
"""

@sjdoc RepeatedNull """
    RepeatedNull
"""

@sjdoc Condition """
    Condition(cond,expr)

a pattern that matches only if `cond` is true.
"""

@sjdoc Except """
    Except(expr)

a pattern than fails to match if `expr` matches.
"""


text:
using StochasticArithmetic

# proposed by Siegfried Rump (rump's royal pain)
# http://www.davidhbailey.com/dhbtalks/dhb-num-repro.pdf
#
# true result: -0.82739605994682136814116
# the result is astronomicaly wrong !
# but consistant across precisions
function rumpTest(T)
    x::T = 77617.0
    y::T = 33096.0
    x2 = x*x
    y2 = y*y
    y4 = y2*y2
    y6 = y4*y2
    y8 = y4*y4
    return 333.75*y6 + x2*(11.0*x2*y2 - y6 - 121.0*y4 - 2.0) + 5.5*y8 + x/(2.0*y)
end

# non-instrumented run
println("result float: ", rumpTest(Float64))

# instrumented (thus likely different from Float64 runs) but we always get the same result as we are using the same seed
println("result Sfloat: ", rumpTest(SFloat64).value)
println("result Sfloat: ", rumpTest(SFloat64).value)
println("result Sfloat: ", rumpTest(SFloat64).value)
println("result Sfloat: ", rumpTest(SFloat64).value)
println("result Sfloat: ", rumpTest(SFloat64).value)

# now varying the seed between runs (getting the usual behaviour of stochastic arithmetic)
StochasticArithmetic.Determinism.resetSeed()
println("result Sfloat (reset): ", rumpTest(SFloat64).value)
StochasticArithmetic.Determinism.resetSeed()
println("result Sfloat (reset): ", rumpTest(SFloat64).value)
StochasticArithmetic.Determinism.resetSeed()
println("result Sfloat (reset): ", rumpTest(SFloat64).value)
StochasticArithmetic.Determinism.resetSeed()
println("result Sfloat (reset): ", rumpTest(SFloat64).value)
StochasticArithmetic.Determinism.resetSeed()
println("result Sfloat (reset): ", rumpTest(SFloat64).value)

# number of digits, computed over 10 runs
println( @reliable_digits rumpTest(SFloat64) )

# this function is exact despite having imprecise components
# it is, wrongly, labelled as inacurate by traditional stochastic arithmetic
function detTest(T)
    x = rumpTest(T) # numerical noise
    y = rumpTest(T) # identical numerical noise
    zero = x - y # exact zero
    return 1. + zero # exact one
end

# this should have an infinity of digits
println( @reliable_digits detTest(SFloat64) )


text:
using Oceananigans.Grids: Face, Face, Center
using Oceananigans.Operators: ζ₃ᶠᶠᵃ
using Oceananigans.AbstractOperations: KernelFunctionOperation

"""
    VerticalVorticityField(model)

Returns a Field that `compute!`s vertical vorticity in a
manner consistent with the `VectorInvariant` momentum advection scheme
for curvilinear grids.

In particular, `VerticalVorticityField` uses `ζ₃ᶠᶠᵃ`, which in turn computes
the vertical vorticity by first integrating the velocity field around the borders
of the vorticity cell to find the vertical circulation, and then dividing by the area of
the vorticity cell to compute vertical vorticity.
"""
function VerticalVorticityField(model; kw...)
    grid = model.grid
    u, v, w = model.velocities
    vorticity_operation = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᵃ, grid, computed_dependencies=(u, v))
    return Field(vorticity_operation; kw...)
end


text:
using Test
using DeZero

function numerical_diff(f::Union{Function,Core.Function}, x::Variable, eps::Float64=1e-4)::Array
    x0 = Variable(x.data .- eps)
    x1 = Variable(x.data .+ eps)
    y0::Variable = f(x0)
    y1::Variable = f(x1)
    return (y1.data - y0.data) / (2 * eps)
end

x = Variable([2.0])
a = square(x)
y = add(square(a), square(a))
backward(y)

println(y.data)
println(x.grad)


text:
# parameters for the case study

# model and data generation files
const modelFile = "resource_20_30.jl"
const dataFile = "demand_model_iid.jl"

# parameters for estimating the true optimal objective value using MC sampling
const numMCScenarios = 1000
const numMCReplicates = 30

# number of covariates
const numCovariates = 3

# determine number of samples depending on the number of covariates
const numDataSamples = Int64[6,8,20,80,400]

# degree of nonlinearity in demand model
const degree = Float64[1,0.5,2]

# scaling factor for the additive demand errors
const demand_errors_scaling = 5.0

# number of folds for K-fold CV to determine optimal k in kNN
const numFolds = 5

# minimum and maximum exponents to determine optimal k in kNN
const minExponent = 0.1
const maxExponent = 0.9

# starting seed
const startingSeed = 1280

max_stars_repo_path: Parameter settings/params_case1_knn.jl

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_path: Parameter settings/params_case1_knn.jl

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_path: Parameter settings/params_case1_knn.jl

max_forks_repo_name: rohitkannan/DD-SAA


text:
import Base: kron, randn, vec, dot, vecdot, IntSet, exp, log, sqrt
import Base.LinAlg: trace, ctranspose, expm!, eigs, scale!
export sincos, eigen, adjoint

const adjoint = ctranspose
const qr! = qrfact!
const luf = lufact
const eigen = eig

exp(A::AbstractMatrix) = expm(A)
exp!(A::AbstractMatrix) = expm!(A)
log(A::AbstractMatrix) = logm(A)
sqrt(A::AbstractMatrix) = sqrtm(A)

mul!(C::AbstractVecOrMat,A::AbstractVecOrMat,B::AbstractVecOrMat) = A_mul_B!(C,A,B)
mul!(C::AbstractVecOrMat,A::AbstractSparseArray,B::AbstractVecOrMat,a::Number,b::Number) = A_mul_B!(a,A,B,b,C)
dot(A::AbstractMatrix,B::AbstractMatrix) = vecdot(A,B)
rmul!(A::AbstractArray,b::Number) = scale!(A,b)

vec(x::RowVector) = x.vec

sincos(x::Number) = (sin(x),cos(x))

parseb2(s::AbstractString) = Base.parse(Int,s,2)

print_array(io::IO, A::AbstractArray) = Base.showarray(io, A, false, header=false)

randn(rng::AbstractRNG,::Type{Complex{T}}) where {T<:AbstractFloat} = Complex{T}(SQRT_HALF*randn(rng,T), SQRT_HALF*randn(rng,T))

Base.@irrational SQRT_HALF 0.70710678118654752  sqrt(big(0.5))

function vecdot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T2,S1,S2}
    m, n = size(A)
    size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions"))
    r = zero(promote_type(T1,T2))
    @inbounds for j = 1:n
        ia = A.colptr[j]; ia_nxt = A.colptr[j+1]
        ib = B.colptr[j]; ib_nxt = B.colptr[j+1]
        if ia < ia_nxt && ib < ib_nxt
            ra = A.rowval[ia]; rb = B.rowval[ib]
            while true
                if ra < rb
                    ia += 1
                    ia < ia_nxt || break
                    ra = A.rowval[ia]
                elseif ra > rb
                    ib += 1
                    ib < ib_nxt || break
                    rb = B.rowval[ib]
                else # ra == rb
                    r += A.nzval[ia]' * B.nzval[ib]
                    ia += 1; ib += 1
                    ia < ia_nxt && ib < ib_nxt || break
                    ra = A.rowval[ia]; rb = B.rowval[ib]
                end
            end
        end
    end
    return r
end

function kron(x::SparseVector{T1,S1}, y::SparseVector{T2,S2}) where {T1,S1,T2,S2}
    nnzx = nnz(x); nnzy = nnz(y)
    nnzz = nnzx*nnzy # number of nonzeros in new vector
    nzind = Vector{promote_type(S1,S2)}(nnzz) # the indices of nonzeros
    nzval = Vector{typeof(one(T1)*one(T2))}(nnzz) # the values of nonzeros
    @inbounds for i = 1:nnzx, j = 1:nnzy
        this_ind = (i-1)*nnzy+j
        nzind[this_ind] = (x.nzind[i]-1)*y.n + y.nzind[j]
        nzval[this_ind] = x.nzval[i] * y.nzval[j]
    end
    return SparseVector(x.n*y.n, nzind, nzval)
end

kron(A::SparseVector, B::VecOrMat) = kron(A, sparse(B))
kron(A::VecOrMat, B::SparseVector) = kron(sparse(A), B)

function kron(A::Diagonal{T1}, B::Diagonal{T2}) where {T1<:Number, T2<:Number}
    valA = A.diag; nA = length(valA)
    valB = B.diag; nB = length(valB)
    valC = Vector{typeof(zero(T1)*zero(T2))}(nA*nB)
    @inbounds for i = 1:nA, j = 1:nB
        valC[(i-1)*nB+j] = valA[i] * valB[j]
    end
    return Diagonal(valC)
end

Base.reduce(::typeof(vcat), A::AbstractVector{<:AbstractVecOrMat}) =
    _typed_vcat(mapreduce(eltype, promote_type, A), A)

function _typed_vcat(::Type{T}, A::AbstractVecOrMat) where T
    nargs = length(A)
    nrows = sum(a->size(a, 1), A)::Int
    ncols = size(A[1], 2)
    for j = 2:nargs
        if size(A[j], 2) != ncols
            throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    B = similar(A[1], T, nrows, ncols)
    pos = 1
    for k=1:nargs
        Ak = A[k]
        p1 = pos+size(Ak,1)-1
        B[pos:p1, :] = Ak
        pos = p1+1
    end
    return B
end


text:
"""
Block
"""

"""
SystemBlock

The structure to represent a system block.

- inports: inports of system block
- outports: outports of system block
- blocks: inner blocks of system block
"""

# struct SystemBlock <: AbstractBlock
#     inports::Dict{Symbol,InPort}
#     outports::Dict{Symbol,OutPort}
#     blocks::Dict{Symbol,Any}
# end

"""
systemequation

Generate function
"""

# function systemequation(system::SystemBlock, label::Symbol = :model)
#     u = gensym()
#     du = Symbol("d", u)
#     states = sort(collect(keys(system.states)))
#     n = length(states)
#     dict = Dict(zip(states, [:($u[$i]) for i = 1:n]))
#     expr = [Expr(:(=), :($(du)[$i]), _replace(dict, _tosystemexpr(system.states[states[i]]))) for i = 1:n]
#     initvec = [_toexpr(system.states[x].initialcondition) for x = states]
#     init = Symbol(label, "init")
#     func = Symbol(label, "!")
#     quote
#         $init = $initvec
#         function $func($(du), $(u), p, t)
#             $(Expr(:block, expr...))
#         end
#     end
# end

# function _systemequation_func(system::SystemBlock)
#     u = gensym()
#     du = Symbol("d", u)
#     states = sort(collect(keys(system.states)))
#     n = length(states)
#     dict = Dict(zip(states, [:($u[$i]) for i = 1:n]))
#     expr = [Expr(:(=), :($(du)[$i]), _replace(dict, _tosystemexpr(system.states[states[i]]))) for i = 1:n]
#     func = quote
#         ($(du), $(u), p, t) -> begin
#             $(Expr(:block, expr...))
#         end
#     end
#     eval(func)
# end

# function _tosystemexpr(blk)
#     s = Set{AbstractBlock}()
#     _toexpr(s, blk.inport)
#     saturationlimits = [_toexpr(x) for x = blk.saturationlimits]
#     if length(saturationlimits) == 0
#         _toexpr(blk.inport[1])
#     else
#         lower, upper = saturationlimits
#         Expr(:if, Expr(:comparison, lower, :<=, blk.state, :<=, uupper), _toexpr(blk.inport[1]), 0)
#     end

# end

"""
AbstractBlock

An abstract object of block.
"""

abstract type AbstractBlock end

"""
Var

A struct of symbol and datatype
"""

struct Var
    label::Symbol
    datatype::DataType
end

"""
Line

An object of line with one symbol
"""

const Line = Var

"""
InPort

An object of an inport. InPort has only at most one line. If InPort has no line, line is Nothing

- line An object of Line. If it has no line, line is Nothing
- parent An object of parent block.
"""

mutable struct InPort
    line::Union{Line,Nothing}
    parent::AbstractBlock
end

InPort(parent::AbstractBlock) = InPort(nothing, parent)

"""
OutPort

An object of an outport. OutPort may have multiple lines.

- lines A vector of Line. If it has no line, line is length 0
- parent An object of parent block.
"""

mutable struct OutPort
    lines::Vector{Line}
    parent::AbstractBlock
end

OutPort(parent::AbstractBlock) = OutPort(Vector{Line}(), parent)

"""
connect

A function to connect OutPort and InPort with a line
"""

function connect(source::OutPort, target::InPort, line::Line)
    push!(source.lines, line)
    target.line = line
end

"""
toexpr
"""

_toexpr(v::Number) = v
_toexpr(v::Symbol) = v
_toexpr(v::Var) = v.label

function _toexpr(p::InPort)
    if isnothing(p.line)
        throw(ErrorException("inport is nothing"))
    else
        p.line.label
    end
end

function _toexpr(p::OutPort)
    Expr(:block, [Expr(:(=), Expr(:(::), x.label, Symbol(x.datatype)), _toexpr(p.parent)) for x = p.lines]...)
end

"""
Constant block
"""

mutable struct ConstantBlock <: AbstractBlock
    outport::OutPort
    value

    function ConstantBlock(value)
        blk = new()
        blk.outport = OutPort(blk)
        blk.value = value
        blk
    end
end

function _toexpr(blk::ConstantBlock)
    _toexpr(blk.value)
end

"""
Gain block
"""

mutable struct GainBlock <: AbstractBlock
    inport::InPort
    outport::OutPort
    K

    function GainBlock(K)
        blk = new()
        blk.inport = InPort(blk)
        blk.outport = OutPort(blk)
        blk.K = K
        blk
    end
end

function _toexpr(blk::GainBlock)
    Expr(:call, :*, _toexpr(blk.K), _toexpr(blk.inport))
end

max_stars_repo_name: JuliaReliab/JuliaMBD

max_issues_repo_name: JuliaReliab/JuliaMBD

max_forks_repo_name: JuliaReliab/JuliaMBD


text:
const order = 5
const a2 = Float64[1/5]
const a3 = Float64[3/40, 9/40]
const a4 = Float64[44/45, -56/15, 32/9]
const a5 = Float64[19372/6561, -25360/2187, 64448/6561, -212/729]
const a6 = Float64[9017/3168, -355/33, 46732/5247, 49/176, -5103/18656]
const a7 = Float64[35/384, 0., 500/1113, 125/192, -2187/6784, 11/84]
const a = {Float64[] a2 a3 a4 a5 a6 a7}
const bs = Float64[5179/57600, 0., 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]
const c = Float64[0., 1/5, 3/10, 4/5, 8/9, 1., 1.]

function substep{T}(x::Vector{T}, x0::Vector{T}, h::Float64, coeffs::Vector{Float64}, k::Vector{Vector{T}})
    @inbounds for m=1:length(x0)
        dx::T = 0.
        @inbounds for i=1:length(coeffs)
            dx += coeffs[i]::Float64*k[i][m]
        end
        x[m] = x0[m] + h*dx
    end
    return nothing
end

function step{T}(F::Function, t::Float64, h::Float64,
                x0::Vector{T}, xp::Vector{T}, xs::Vector{T}, k::Vector{Vector{T}})
    for i=2:length(c)
        substep(xp, x0, h, a[i], k)
        F(t + h*c[i], xp, k[i])
    end
    substep(xs, x0, h, bs, k)
    return nothing
end

function allocate_memory{T}(x::Vector{T})
    xp = zeros(T, length(x))
    xs = zeros(T, length(x))
    k = Vector{Vector{T}}[]
    for i=1:7
        push!(k, zeros(T, length(x))
    return xp, xs, k
end

function error_estimate(xp, xs, abstol, reltol)
    err::Float64 = 0.
    for i=1:length(xp)
        sc_i = abstol + reltol*max(abs(xp[i]), abs(xs[i]))
        err += abs2(xp[i]-xs[i])/sc_i^2
    return sqrt(err/length(xp))
end

function initial_stepsize(F, x, k, tmp1, tmp2)
    d0 = norm(x,2)
    d1 = norm(k[1],2)
    h0 = ((d0<1e-5 | d1<1e-5) ? 1e-6 : 0.01*d0/d1)
    substep(tmp1, x, h0, [1.], k)
    F(t+h0, tmp1, tmp2)
    for i=1:length(tmp2)
        tmp2[i] -= k[1][i]
    end
    d2 = norm(tmp2, 2)
    if max(d1,d2)<1e-15
        h1 = max(1e-6, h0*1e-3)
    else
        h1 = (0.01/max(d1,d2))^(1./order)
    return min(100*h0, h1)
end

function ode{T}(F, tspan::Vector{Float64}, x0::Vector{T};
                    reltol = 1.0e-5,
                    abstol = 1.0e-8,
                    h0 = 0,
                    hmin = (tspan[end]-tspan[1])/1e9,
                    hmax = (tspan[end]-tspan[1]),
                    fout=(t,x)->nothing,)
    t, tend = tspan[1], tspan[end]

    fout(t, x)
    x = 1*x0
    xp, xs, k = allocate_memory(x0)

    F(t,x,k[1])
    h = (h0==0 ? initial_stepsize(F, x, k, k[2], k[3]) : h0)
    was_rejected = false
    while t < tfinal
        step(F, t, h, x, xp, xs, k)
        err = error_estimate(xp, xs, abstol, reltol)
        facmin = (was_rejected ? 1. : 5.)
        hnew = h*min(facmin, max(0.2, 0.9*(1./err)^(1./order)))
        hnew = min(hmax, hnew)
        if hnew<hmin
            error("Stepsize below hmin.")
        end
        if err>1
            was_rejected = true
            continue
        end
        was_rejected = false
        xp, x = x, xp
        fout(t, x)
    end
    return nothing
end


text:
"""
    LNR(;μ, σ, ϕ)

A lognormal race model object

# Fields

- `μ`: a vector of means in log-space
- `σ`: a standard deviation parameter in log-space
- `ϕ`: a encoding-response offset

# Example

```julia
using SequentialSamplingModels
dist = LNR(μ=[-2,-3], σ=1.0, ϕ=.3)
data = rand(dist, 10)
like = pdf.(dist, data)
loglike = logpdf.(dist, data)
```
# References

Rouder, J. N., Province, J. M., Morey, R. D., Gomez, P., & Heathcote, A. (2015).
The lognormal race: A cognitive-process model of choice and latency with desirable
psychometric properties. Psychometrika, 80(2), 491-513.
"""
struct LNR{T1,T2,T3} <: SequentialSamplingModel
    μ::T1
    σ::T2
    ϕ::T3
end

Broadcast.broadcastable(x::LNR) = Ref(x)

LNR(;μ, σ, ϕ) = LNR(μ, σ, ϕ)

function rand(dist::LNR)
    @unpack μ,σ,ϕ = dist
    x = @. rand(LogNormal(μ, σ)) + ϕ
    rt,resp = findmin(x)
    return resp,rt
end

rand(dist::LNR, N::Int) = [rand(dist) for i in 1:N]

function logpdf(d::LNR, r::Int, t::Float64)
    @unpack μ,σ,ϕ = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == r
            LL += logpdf(LogNormal(m, σ), t - ϕ)
        else
            LL += logccdf(LogNormal(m, σ), t - ϕ)
        end
    end
    return LL
end

function logpdf(d::LNR{T1,T2,Vector{T3}}, r::Int, t::Float64) where {T1,T2,T3}
    @unpack μ,σ,ϕ = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == r
            LL += logpdf(LogNormal(m, σ), t - ϕ[i])
        else
            LL += logccdf(LogNormal(m, σ), t - ϕ[i])
        end
    end
    return LL
end

logpdf(d::LNR, data::Tuple) = logpdf(d, data...)

pdf(d::LNR, data::Tuple) = pdf(d, data...)

function pdf(d::LNR, r::Int, t::Float64)
    @unpack μ,σ,ϕ = d
    density = 1.0
    for (i,m) in enumerate(μ)
        if i == r
            density *= pdf(LogNormal(m, σ), t - ϕ)
        else
            density *= (1 - cdf(LogNormal(m, σ), t - ϕ))
        end
    end
    return density
end