text:
using Flux
println("NNPDE_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
println("Starting Soon!")
using GalacticOptim
using Optim
using Quadrature,Cubature, Cuba
using QuasiMonteCarlo

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

## Example 1, 1D ode
@parameters θ
@variables u(..)
@derivatives Dθ'~θ

# 1D ODE
eq = Dθ(u(θ)) ~ θ^3 + 2*θ + (θ^2)*((1+3*(θ^2))/(1+θ+(θ^3))) - u(θ)*(θ + ((1+3*(θ^2))/(1+θ+θ^3)))

# Initial and boundary conditions
bcs = [u(0.) ~ 1.0]

# Space and time domains
domains = [θ ∈ IntervalDomain(0.0,1.0)]
# Discretization
dt = 0.1
# Neural network
chain = FastChain(FastDense(1,12,Flux.σ),FastDense(12,1))

discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             nothing; #init_params
                                             phi = nothing,
                                             autodiff=false,
                                             derivative = nothing,
                                             strategy = NeuralPDE.GridTraining(dx=dt))

pde_system = PDESystem(eq,bcs,domains,[θ],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1); cb = cb, maxiters=1000)
phi = discretization.phi

analytic_sol_func(t) = exp(-(t^2)/2)/(1+t+t^3) + t^2
ts = [domain.domain.lower:dt/10:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(phi(t,res.minimizer)) for t in ts]

@test u_predict ≈ u_real atol = 0.2

# using Plots
# t_plot = collect(ts)
# plot(t_plot ,u_real)
# plot!(t_plot ,u_predict)

## Example 2, 2D Poisson equation
function test_2d_poisson_equation(chain, strategy)
    @parameters x y
    @variables u(..)
    @derivatives Dxx''~x
    @derivatives Dyy''~y

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

    # Initial and boundary conditions
    bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
           u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    # Space and time domains
    domains = [x ∈ IntervalDomain(0.0,1.0),
               y ∈ IntervalDomain(0.0,1.0)]

    # chain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))
    discretization = NeuralPDE.PhysicsInformedNN(chain,
                                                 strategy = strategy)

    pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
    prob = NeuralPDE.discretize(pde_system,discretization)

    res = GalacticOptim.solve(prob, ADAM(0.1); cb = cb, maxiters=500)
    phi = discretization.phi

    dx = 0.1
    xs,ys = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains]
    analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)

    u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
    u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
    diff_u = abs.(u_predict .- u_real)

    @test u_predict ≈ u_real atol = 3.0

    # p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
    # p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
    # p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
    # plot(p1,p2,p3)
end

# Discretization
dx = 0.1
grid_strategy = NeuralPDE.GridTraining(dx=dx)
fastchain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))
fluxchain = Chain(Dense(2,12,Flux.σ),Dense(12,12,Flux.σ),Dense(12,1))
chains = [fluxchain, fastchain]
for chain in chains
    test_2d_poisson_equation(chain, grid_strategy)
end

stochastic_strategy = NeuralPDE.StochasticTraining(number_of_points = 100)
quadrature_strategy = NeuralPDE.QuadratureTraining(algorithm=HCubatureJL(),
                                                   reltol = 1e-2, abstol = 1e-2,
                                                   maxiters = 50)
quasirandom_strategy = NeuralPDE.QuasiRandomTraining(sampling_method = UniformSample(),
                                                     number_of_points = 100,
                                                     number_of_minibatch = 10)

strategies = [stochastic_strategy, quadrature_strategy,quasirandom_strategy]
for strategy in strategies
    chain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))
    test_2d_poisson_equation(chain, strategy)
end

function run_2d_poisson_equation(strategy)
    @parameters x y
    @variables u(..)
    @derivatives Dxx''~x
    @derivatives Dyy''~y

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

    # Initial and boundary conditions
    bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
           u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    # Space and time domains
    domains = [x ∈ IntervalDomain(0.0,1.0),
               y ∈ IntervalDomain(0.0,1.0)]
    # Discretization
    dx = 0.1

    chain = FastChain(FastDense(2,4,Flux.σ),FastDense(4,1))

    discretization = NeuralPDE.PhysicsInformedNN(chain,
                                                 strategy = strategy)

    pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
    prob = NeuralPDE.discretize(pde_system,discretization)

    res = GalacticOptim.solve(prob, ADAM(0.01); cb = cb,  maxiters=2)
end

algs = [CubaVegas(), CubaSUAVE(),HCubatureJL(), CubatureJLh(), CubatureJLp()]
#CubaDivonne(),CubaCuhre() doesn't work with dim = 2
for alg in algs
    strategy =  NeuralPDE.QuadratureTraining(algorithm = alg,reltol=1e-8,abstol=1e-8,maxiters=600)
    run_2d_poisson_equation(strategy)
end

## Example 3, high-order ode
@parameters x
@variables u(..)
@derivatives Dxxx'''~x
@derivatives Dx'~x

# ODE
eq = Dxxx(u(x)) ~ cos(pi*x)

# Initial and boundary conditions
bcs = [u(0.) ~ 0.0,
       u(1.) ~ cos(pi),
       Dx(u(1.)) ~ 1.0]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.05

# Neural network
chain = FastChain(FastDense(1,8,Flux.σ),FastDense(8,1))

discretization = NeuralPDE.PhysicsInformedNN(chain,strategy = NeuralPDE.GridTraining(dx=dx))
pde_system = PDESystem(eq,bcs,domains,[x],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.01); cb = cb, maxiters=2000)
phi = discretization.phi

analytic_sol_func(x) = (π*x*(-x+(π^2)*(2*x-3)+1)-sin(π*x))/(π^3)

xs = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(x) for x in xs]
u_predict  = [first(phi(x,res.minimizer)) for x in xs]

@test u_predict ≈ u_real atol = 1.0

# x_plot = collect(xs)
# plot(x_plot ,u_real)
# plot!(x_plot ,u_predict)

## Example 4, system of pde
@parameters x, y
@variables u1(..), u2(..)
@derivatives Dx'~x
@derivatives Dy'~y

# System of pde
eqs = [Dx(u1(x,y)) + 4*Dy(u2(x,y)) ~ 0,
      Dx(u2(x,y)) + 9*Dy(u1(x,y)) ~ 0,
      3*u1(x,0) ~ 2*u2(x,0)]

# Initial and boundary conditions
bcs = [u1(x,0) ~ 2x, u2(x,0) ~ 3x]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]

# Neural network
chain1 = FastChain(FastDense(2,8,Flux.σ),FastDense(8,1))
chain2 = FastChain(FastDense(2,8,Flux.σ),FastDense(8,1))

discretization = NeuralPDE.PhysicsInformedNN([chain1,chain2],strategy = NeuralPDE.GridTraining(dx=0.1))
pde_system = PDESystem(eqs,bcs,domains,[x,y],[u1,u2])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=300)
phi = discretization.phi

analytic_sol_func(x,y) =[1/3*(6x - y), 1/2*(6x - y)]
xs,ys = [domain.domain.lower:0.01:domain.domain.upper for domain in domains]
u_real  = [[analytic_sol_func(x,y)[i] for x in xs  for y in ys] for i in 1:2]

initθ = discretization.initθ
acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers = [res.minimizer[s] for s in sep]
u_predict  = [[phi[i]([x,y],minimizers[i])[1] for x in xs  for y in ys] for i in 1:2]

@test u_predict ≈ u_real atol = 10.0

# p1 =plot(xs, ys, u_predict, st=:surface);
# p2 = plot(xs, ys, u_real, st=:surface);
# plot(p1,p2)

## Example 5, 2d wave equation, neumann boundary condition
#here we use low level api for build solution
@parameters x, t
@variables u(..)
@derivatives Dxx''~x
@derivatives Dtt''~t
@derivatives Dt'~t

#2D PDE
C=1
eq  = Dtt(u(x,t)) ~ C^2*Dxx(u(x,t))

# Initial and boundary conditions
bcs = [u(0,t) ~ 0.,# for all t > 0
       u(1,t) ~ 0.,# for all t > 0
       u(x,0) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(x,0)) ~ 0. ] #for all  0 < x < 1]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           t ∈ IntervalDomain(0.0,1.0)]

# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

phi = NeuralPDE.get_phi(chain)
derivative = NeuralPDE.get_numeric_derivative()
initθ = DiffEqFlux.initial_params(chain)

indvars = [x,t]
depvars = [u]
dim = length(domains)

_pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,phi, derivative,initθ)

bc_indvars = NeuralPDE.get_bc_varibles(bcs,indvars,depvars)
_bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars, phi, derivative,initθ,
                                              bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

train_sets = NeuralPDE.generate_training_sets(domains,dx,bcs,indvars,depvars)
pde_train_set,bcs_train_set,train_set = train_sets
pde_bounds, bcs_bounds = NeuralPDE.get_bounds(domains,bcs,indvars,depvars)

quadrature_strategy = NeuralPDE.QuadratureTraining(algorithm=HCubatureJL(),reltol= 1e-3,abstol= 1e-3,maxiters=20)
pde_loss_function = NeuralPDE.get_loss_function(_pde_loss_function,
                                                pde_bounds,
                                                quadrature_strategy)

dx = 0.1
grid_strategy = NeuralPDE.GridTraining(dx=dx)
bc_loss_function = NeuralPDE.get_loss_function(_bc_loss_functions,
                                               bcs_train_set,
                                               grid_strategy)

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

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

res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=400)

xs,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
analytic_sol_func(x,t) =  sum([(8/(k^3*pi^3)) * sin(k*pi*x)*cos(C*k*pi*t) for k in 1:2:50000])

u_predict = reshape([first(phi([x,t],res.minimizer)) for x in xs for t in ts],(length(xs),length(ts)))
u_real = reshape([analytic_sol_func(x,t) for x in xs for t in ts], (length(xs),length(ts)))

@test u_predict ≈ u_real atol = 1.0

# diff_u = abs.(u_predict .- u_real)
# p1 = plot(xs, ts, u_real, linetype=:contourf,title = "analytic");
# p2 =plot(xs, ts, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(xs, ts, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)

## Example 6, pde with mixed derivative
@parameters x y
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y
@derivatives Dy'~y
@derivatives Dx'~x

eq = Dxx(u(x,y)) + Dx(Dy(u(x,y))) - 2*Dyy(u(x,y)) ~  -1.

# Initial and boundary conditions
bcs = [u(x,0) ~ x,
       Dy(u(x,0)) ~ x,
       u(x,0) ~ Dy(u(x,0))]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.1
quadrature_strategy = NeuralPDE.QuadratureTraining(algorithm=HCubatureJL(),reltol= 1e-2,abstol= 1e-2,maxiters=10)
# Neural network
chain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))
discretization = NeuralPDE.PhysicsInformedNN(chain, strategy = quadrature_strategy)
pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=200)
phi = discretization.phi

analytic_sol_func(x,y) = x + x*y +y^2/2
xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)

@test u_predict ≈ u_real atol = 1.0

# p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
# p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)


text:
export gradiente
export newton
export newton_rc

function gradiente(
  nlp :: AbstractNLPModel;
  max_time = 10.0,
  max_iter = 10_000,
  η₁ = 1e-2,
  atol = 1e-6,
  rtol = 1e-6,
)

t₀ = time()
Δt = time() - t₀
iter = 0

x = copy(nlp.meta.x0)
f(x) = obj(nlp, x)
∇f(x) = grad(nlp, x)
fx = f(x)
gx = ∇f(x)

x⁺ = similar(x)

ϵ = atol + rtol * norm(gx)

status = :unknown

# Cond. de parada
resolvido = norm(gx) < ϵ # true ou false
cansado = Δt > max_time || iter > max_iter
while !(resolvido || cansado) # Não tiver satisfeito as condições de parada
  # Calculo a direção e passo e xₖ₊₁
  d = -gx
  slope = dot(gx, d)

  α = 1.0
  x⁺ = x + α * d
  f⁺ = f(x⁺)
  while f⁺ ≥ fx + η₁ * α * slope
      α = α / 2
      x⁺ = x + α * d
      f⁺ = f(x⁺)
      if α < 1e-8
          status = :small_step
          break
      end
  end
  if status != :unknown
      break
  end
  x .= x⁺
  fx = f⁺
  gx = ∇f(x)

  resolvido = norm(gx) < ϵ
  Δt = time() - t₀
  iter += 1
  cansado = Δt > max_time || iter > max_iter
end

if resolvido
  status = :first_order
elseif cansado
  if Δt > max_time
      status = :max_time
  elseif iter > max_iter
      status = :max_iter
  end
end

return GenericExecutionStats(status, nlp, objective=fx,
  solution=x, dual_feas=norm(gx), iter=iter, elapsed_time=Δt)
end

##

function newton(
  nlp :: AbstractNLPModel;
  max_time = 10.0,
  max_iter = 10_000,
  η₁ = 1e-2,
  atol = 1e-6,
  rtol = 1e-6,
)

t₀ = time()
Δt = time() - t₀
iter = 0

x = copy(nlp.meta.x0)
f(x) = obj(nlp, x)
∇f(x) = grad(nlp, x)
H(x) = Symmetric(hess(nlp, x), :L)
fx = f(x)
gx = ∇f(x)

x⁺ = similar(x)

ϵ = atol + rtol * norm(gx)

status = :unknown

# Cond. de parada
resolvido = norm(gx) < ϵ # true ou false
cansado = Δt > max_time || iter > max_iter
while !(resolvido || cansado) # Não tiver satisfeito as condições de parada
  # Calculo a direção e passo e xₖ₊₁
  Hx = H(x)
  F = cholesky(Hx, check=false)
  if !issuccess(F)
      status = :not_desc
      break
  end
  d = -(F \ gx)
  slope = dot(gx, d)

  α = 1.0
  x⁺ = x + α * d
  f⁺ = f(x⁺)
  while f⁺ ≥ fx + η₁ * α * slope
      α = α / 2
      x⁺ = x + α * d
      f⁺ = f(x⁺)
      if α < 1e-8
          status = :small_step
          break
      end
  end
  if status != :unknown
      break
  end
  x .= x⁺
  fx = f⁺
  gx = ∇f(x)

  resolvido = norm(gx) < ϵ
  Δt = time() - t₀
  iter += 1
  cansado = Δt > max_time || iter > max_iter
end

if resolvido
  status = :first_order
elseif cansado
  if Δt > max_time
      status = :max_time
  elseif iter > max_iter
      status = :max_iter
  end
end

return GenericExecutionStats(status, nlp, objective=fx,
  solution=x, dual_feas=norm(gx), iter=iter, elapsed_time=Δt)
end

##

function newton_rc(
  nlp::AbstractNLPModel;
  atol::Real = 1e-6,
  rtol::Real = 1e-6,
  max_eval::Int = 1000,
  max_iter::Int = 0,
  max_time::Float64 = 10.0
)

  if !unconstrained(nlp)
    error("Problem is not unconstrained")
  end

  x = copy(nlp.meta.x0)

  f(x) = obj(nlp, x)
  ∇f(x) = grad(nlp, x)
  H(x) = Symmetric(hess(nlp, x), :L)

  fx = f(x)
  ∇fx = ∇f(x)

  ϵ = atol + rtol * norm(∇fx)
  t₀ = time()

  iter = 0
  Δt = time() - t₀
  solved = norm(∇fx) < ϵ # First order stationary
  tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

  # status must be one of a few options found in SolverTools.show_statuses()
  # A good default value is :unknown.
  status = :unknown

  # log_header is up for some rewrite in the future. For now, it simply prints the column names with some spacing
  @info log_header(
    [:iter, :fx, :ngx, :nf, :Δt],
    [Int, Float64, Float64, Int, Float64],
    hdr_override=Dict(:fx => "f(x)", :ngx => "‖∇f(x)‖", :nf => "#f")
  )
  # log_row uses the type information of each value, thus we use `Any` here.
  @info log_row(
    Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
  )

  # This template implements a simple steepest descent method without any hopes of working.
  # This is where most of your change will happen
  η₁ = 1e-2
  η₂ = 0.75
  Δ = 1.0
  iter = 0

  while !(solved || tired)
    fx = f(x)
    gx = ∇f(x)
    Hx = H(x)
    model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
    @variable(model, d[1:2])
    @objective(model, Min, fx + dot(d, gx) + dot(d, Hx * d) / 2)
    @NLconstraint(model, d[1]^2 + d[2]^2 ≤ Δ^2)
    optimize!(model)
    d = value.(d)

    Ared = f(x) - f(x + d)
    Pred = f(x) - (fx + dot(d, gx) + dot(d, Hx * d) / 2)
    ρ = Ared / Pred
    if ρ < η₁
        Δ = Δ / 2
    elseif ρ < η₂
        x = x + d
    else
        x = x + d
        Δ = 2Δ
    end

    iter += 1
    Δt = time() - t₀
    solved = norm(∇fx) < ϵ # First order stationary
    tired = neval_obj(nlp) ≥ max_eval > 0|| #Evaluations
            iter ≥ max_iter > 0 || #iterations
            Δt ≥ max_time > 0 # Excess time

    @info log_row(
      Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
    )
  end

  if solved
    status = :first_order
  elseif tired
    if neval_obj(nlp) ≥ max_eval > 0
      status = :max_eval
    elseif iter ≥ max_iter > 0
      status = :max_iter
    elseif Δt ≥ max_time > 0
      status = :max_time
    end
  end

  return GenericExecutionStats(
    status,
    nlp,
    solution=x,
    objective=f(x),
    dual_feas=norm(∇fx),
    elapsed_time=Δt,
    iter=iter
  )
end

ext: jl


text:
#############################################################################
# lovasz_extension.jl
# Handles the Lovasz extensions of set functions.
#############################################################################

export lovasz, LovaszExtAbsAtom
export sign, monotonicity, curvature, evaluate
export ConvexLovaszAbs

mutable struct LovaszExtAbsAtom <: AbstractExpr
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr, Variable}
  size::Tuple{Int, Int}
  variable::Variable
  func::SubmodFunc

  function LovaszExtAbsAtom(F::SubmodFunc, x::AbsAtom)
    if length(F.setvariables) == 1
      if length(x.children) == 1 && isa(x.children[1], Variable)
        if evaluate(F, [])[1] != 0
          error("A combinatorial function should be 0 at the empty set to derive its Lovasz extension.")
        else
          var = get_v(x)
          if length(var) != 1
            error("Functions with other than one variable are not supported.")
          else
            if var[1].size[1] == F.setvariables[1].cardinality
              children = (F, var[1])
              return new(:lovaszabs, hash(children), children, (1, 1), var[1], F)
            else
              error("The size of the continuous variable should be the same as the baseset of the combinatorial variable of the combinatorial function.")
            end
          end
        end
      else
        error("Only able to define Lovasz extesions on absolute values of a variable.")
      end
    else
      error("A combinatorial function should be sigle-variant to obtain its Lovasz extension.")
    end
  end
end

lovasz(F::SubmodFunc, x::AbsAtom) = LovaszExtAbsAtom(F, x)

function sign(F::LovaszExtAbsAtom)
  return sign(F.children[1])
end

function monotonicity(F::LovaszExtAbsAtom)
  return monotonicity(F.children[1])
end

function curvature(F::LovaszExtAbsAtom)
  modd = modularity(F.children[1])
  if modd == SubModularity()
    return ConvexVexity()
  elseif modd == SuperModularity()
    return ConcaveVexity()
  elseif modd == Modularity()
    return AffineVexity()
  elseif modd == ConstModularity()
    return ConstVexity()
  end
end

vexity(F::LovaszExtAbsAtom) = curvature(F)

function evaluate(f::LovaszExtAbsAtom)
  n = f.children[2].size[1]
  y = abs.(f.children[2].value[:])
  i = sortperm(y, rev = true)
  V = sort(f.children[1].setvariables[1].baseset)
  storage = zeros(n + 1)
  x = zeros(n)
  for ii = 2:n+1
    storage[ii] = evaluate(f.children[1], V[i[1:ii-1]])[1]
    x[i[ii - 1]] = storage[ii] - storage[ii - 1]
  end
  # Compatibility with Convex.jl
  val = zeros(1, 1)
  val[1] = dot(x, y)
  return val
end

function get_cv(F::LovaszExtAbsAtom)
  variable = []
  return push!(variable, F.variable)
end

function get_sv(F::LovaszExtAbsAtom)
  return []
end

function get_v(F::LovaszExtAbsAtom)
  variable = []
  return push!(variable, F.variable)
end


text:
# Copies Corner plots to figures directory:

run(`cp ../static/corner_photdyn.png  ..`)
run(`cp ../static/corner_ttv.png  ..`)


text:
using InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.8.0")

max_stars_repo_name: ykamit/lecture-source-jl

max_issues_repo_name: ykamit/lecture-source-jl

max_forks_repo_name: ykamit/lecture-source-jl


text:
module FlightSimulator

#using DifferentialEquations

export SSB, MercuryBarycenter, VenusBarycenter, EarthMoonBarycenter, MarsBarycenter, JupiterBarycenter
export SaturnBarycenter, NeptuneBarycenter, UranusBarycenter, PlutoBarycenter
export Sun, Mercury, Venus, Earth, Moon, Mars, Jupiter, Saturn, Neptune, Uranus, Pluto
export SunType, MercuryType, VenusType, EarthType, MoonType, MarsType, JupiterType, SaturnType, NeptuneType, UranusType, PlutoType
export abstractCelestialBody
export PlanetarySystem, GetInertialFramePos, GetGravAccel, LoadKernels, naif_id, parent

#include global Dictionaries
#include("solver\\GlobalDictionaries.jl")

#Include utility functions
include("utils/LoadUtils.jl")

#Include base function first
include("Environment/LoadEnvironment.jl")

#simulator Options
#include("solver\\SimulationOptions.jl")

#Model Properties Type
#include("ModelOptions.jl")

#Include the Aerodynamic database files
#include("Vehicle\\AeroDataBase.jl")
#include("Vehicle\\VehicleAeroDataBase.jl")

#Add Planets
#include("Planets.jl")

#Add Spacecrafts
#include("Vehicle\\Spacecraft.jl")

#Add Nozzle
#include("Vehicle\\Nozzle.jl")

#Add Actuators
#include("Vehicle\\Actuators.jl")

#include run simulation ability
#include("RunSimulation.jl")

end #module FlightSimulator


text:
using SampleChainsWeighted
using Test

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

max_stars_repo_name: cscherrer/SampleChainsWeighted

max_issues_repo_name: cscherrer/SampleChainsWeighted

max_forks_repo_name: cscherrer/SampleChainsWeighted


text:
import ..Utils: UnionFind, find, unite!

"""
    lowest_common_ancestors(tree, parent, pairs)

For a list of pairs `[(u1, v1), ..., (uk, vk)]` compute the
[lowest common ancestors](https://en.wikipedia.org/wiki/Lowest_common_ancestor)
`[a1, ..., ak]`, i.e. `ai` is the lowest node in the tree that is
ancestor (parent, grand-parent, grand-grand-parent, etc.) of both, `ui` and `vi`.

Runtime is `O(length(tree) + length(pairs))`.
More precisely, due to the use of the
[`UnionFind`](@ref) data structure, the inverse Ackermann function is also included;
however the impact is at most 4 for all practical instances (``n < 10^{100}``).

See CLRS (3rd edition) page 584 or
[`networkx`](https://github.com/networkx/networkx/blob/master/networkx/algorithms/lowest_common_ancestors.py).
"""
function lowest_common_ancestors(
    tree::ChildrenIndex,
    parent::Vector{Int},
    pairs::IncidenceIndex,
)::Vector{Int}
    @assert parent[root_node(tree)] == root_node(tree)

    local k = num_edges(pairs)      # number of pairs
    local n = num_nodes(tree)
    local uf = UnionFind(n)
    local colors = fill(false, n)
    local ancestors = Vector{Int}(undef, n)
    local stack = Vector{Int}()
    lcas = Vector{Int}(undef, k)
    sizehint!(stack, n)

    lcas .= -1
    ancestors .= 1:n
    dfs_walk(tree, stack) do v::Int
        if v >= 0   # v just finished
            for (u, ei) in pairs[v]
                if colors[u]
                    lcas[ei] = ancestors[uf[u]]
                end
            end

            colors[v] = true
            if v != parent[v]
                unite!(uf, uf[v], uf[parent[v]])
                ancestors[uf[parent[v]]] = parent[v]
            end
        end
    end

    return lcas
end

lowest_common_ancestors(tree, parent, pairs::Vector{Tuple{Int, Int}}) =
    lowest_common_ancestors(tree, parent, IncidenceIndex(length(tree), pairs))

lowest_common_ancestors(tree, parent, pairs::Vector{Pair{Int, Int}}) =
    lowest_common_ancestors(tree, parent, IncidenceIndex(length(tree), pairs))

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
@testset "test_transferfunction" begin
# Naming convention:
# ------------------
# {type}_{dims}
# type: C: Continuous, D: Discrete
# dims: "npnuny" (np = # poles)

# CONTINUOUS
C_111 = tf([1, 2], [1, 5])
C_211 = tf([1, 2, 3], [1, 8, 15])
C_212 = tf(vecarray(2, 1,[1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]))
C_221 = tf(vecarray(1, 2,[1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]))
C_222 = [C_221; C_221]
C_022 = tf(4*[1 0; 0 1])
s = tf("s")

# DISCRETE
D_111 = tf([1, 2], [1, -0.5], 0.005)
D_211 = tf([1, 2, 3], [1, -0.2, -0.15], 0.005)
D_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, -0.2, -0.15], [1, -0.2, -0.15]), 0.005)
D_222 = [D_221; D_221]
D_022 = tf(4*[1 0; 0 1], 0.005)
z = tf("z", 0.005)

# TESTS
# Constructors
@test s == tf([1, 0], [1])
@test z == tf([1, 0], [1], 0.005)
@test C_022 == tf(vecarray(2, 2, [4], [0], [0], [4]), vecarray(2, 2, [1], [1], [1], [1]))
@test D_022 == tf(vecarray(2, 2, [4], [0], [0], [4]), vecarray(2, 2, [1], [1], [1], [1]), 0.005)
@test C_022 == [tf(4) 0;0 4]
@test C_022 == tf([4 0;0 4])
@test D_022 == tf([4 0;0 4], 0.005)

# Test equality
@test tf([1,2], [2,3,4]) == tf(2*[1,2], 2*[2,3,4])
@test tf([1.0], [2.0,3.0]) == tf(π*[1.0], π*[2.0,3.0])
@test tf([1.0+2.0im], [2.0+im,3.0]) == tf((π+im)*[1+2.0im], (π+im)*[2.0+im,3.0])

# Test approximate equlity
# rtol should just be on the order of ϵ, no particular reason that exactly ϵ
# would work, but apparently it does
ϵ = 1e-14
@test tf([1,2], [2,3,4]) != tf(2*[1,2], (2+ϵ)*[2,3,4])
@test tf([1,2], [2,3,4]) ≈   tf(2*[1,2], (2+ϵ)*[2,3,4]) rtol=ϵ
@test tf([1.0], [2.0,3.0]) != tf((π+ϵ)*[1.0], π*[2.0,3.0])
@test tf([1.0], [2.0,3.0]) ≈ tf((π+ϵ)*[1.0], π*[2.0,3.0]) rtol=ϵ
@test tf([1.0+2.0im], [2.0+im,3.0]) != tf((π+(1+ϵ)im)*[1+2.0im], (π+ϵ+im)*[2.0+(1+ϵ*im),3.0])
@test tf([1.0+2.0im], [2.0+im,3.0]) ≈ tf((π+(1+ϵ)im)*[1+2.0im], (π+ϵ+im)*[2.0+(1+ϵ)*im,3.0]) rtol=ϵ

# Addition
@test C_111 + C_111 == tf([2,14,20], [1,10,25])
@test C_222 + C_222 == tf(vecarray(2, 2, [2,20,68,108,90], [0,2,20,62,60],
              [2,20,68,108,90], [0,2,20,62,60]),
  vecarray(2, 2, [1,16,94,240,225], [1,16,94,240,225],
              [1,16,94,240,225], [1,16,94,240,225]))
@test C_222 + 1 == tf(vecarray(2, 2, [2,10,18], [1,9,17], [2,10,18], [1,9,17]),
  vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15]))
@test D_111 + D_111 == tf([2,3,-2], [1,-1,0.25], 0.005)

# Subtraction
@test C_111 - C_211 == tf([0,3,18,15], [1,13,55,75])
@test 1 - C_222 == tf(vecarray(2, 2, [0,6,12], [1,7,13], [0,6,12], [1,7,13]), vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15]))
# We are not doing enough to identify zero numerator here
@test_broken D_111 - D_211 - tf([0,0.3,-2.55,1.2], [1,-0.7,-0.05,0.075], 0.005) == tf([0.0], [1], 0.005)

# Multiplication
@test C_111 * C_221 == tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]),
  vecarray(1, 2, [1,13,55,75], [1,13,55,75]))
@test C_212 * C_111 == tf(vecarray(2, 1, [1,4,7,6], [0,1,4,4]),
  vecarray(2, 1, [1,13,55,75], [1,13,55,75]))
@test 4*C_222 == tf(vecarray(2, 2, [4,8,12], [0,4,8], [4,8,12], [0,4,8]),
  vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15]))
 # We are not doing enough to identify zero numerator here
@test_broken D_111 * D_221 - tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]),
  vecarray(1, 2, [1,-0.7,-0.05,0.075], [1.0,-0.7,-0.05,0.075]), 0.005) ==
tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)

# Division
@test 1/C_111 == tf([1,5], [1,2])
@test C_212/C_111 == tf(vecarray(2, 1, [1,7,13,15], [0,1,7,10]),
  vecarray(2, 1, [1,10,31,30], [1,10,31,30]))
@test 1/D_111 == tf([1.0,-0.5], [1.0,2.0], 0.005)

# Indexing
@test size(C_222) == (2, 2)
@test size(C_212) == (2, 1)
@test C_222[1,1] == tf([1, 2, 3], [1, 8, 15])
@test C_222[1:1,1] == tf([1, 2, 3], [1, 8, 15])
@test C_222[1,1:2] == C_221
@test size(C_222[1,[]]) == (1,0)

# Errors
@test_throws ErrorException tf(vecarray(1, 1, [1,7,13,15]),
   vecarray(2, 1, [1,10,31,30], [1,10,31,30]))

# Printing
res = ("TransferFunction:\nInput 1 to Output 1\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 1 to Output 2\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 2 to Output 1\n    s + 2\n-------------\ns^2 + 8s + 15\n\nInput 2 to Output 2\n    s + 2\n-------------\ns^2 + 8s + 15\n\nContinuous-time transfer function model")
@test_broken sprint(show, C_222) == res
res = ("TransferFunction:\nInput 1 to Output 1\nz^2 + 2.0z + 3.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 1 to Output 2\nz^2 + 2.0z + 3.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 1\n     z + 2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 2\n     z + 2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nSample Time: 0.005 (seconds)\nDiscrete-time transfer function model")
@test_broken sprint(show, D_222) == res

@test tf(zpk([1.0 2; 3 4])) == tf([1 2; 3 4])

# Type stability Continuous-time
@test eltype(fill(tf("s"),2)) <: TransferFunction
@test eltype(fill(tf([1],[1,1]),2)) <: TransferFunction
@test eltype(fill(tf(1.0,[1,1]),2)) <: TransferFunction
@test eltype(fill(tf([1 2; 3 4]),2)) <: TransferFunction
@test eltype(fill(tf(1)+tf(2),2)) <: TransferFunction
@test eltype(fill(tf(1)/tf(2),2)) <: TransferFunction
@test eltype(fill(tf(1)+1,2)) <: TransferFunction

# Type stability Discrete-time
@test eltype(fill(tf("z",1.0),2)) <: TransferFunction
@test eltype(fill(tf([1],[1,1],1),2)) <: TransferFunction
@test eltype(fill(tf(1.0,[1,1],1),2)) <: TransferFunction
@test eltype(fill(tf([1 2; 3 4],1),2)) <: TransferFunction
@test eltype(fill(tf(1,1)+tf(2,1),2)) <: TransferFunction
@test eltype(fill(tf(1,1)/tf(2,1),2)) <: TransferFunction
@test eltype(fill(tf(1,1.0)+1,2)) <: TransferFunction

# Errors
@test_throws ErrorException C_111 + C_222             # Dimension mismatch
@test_throws ErrorException C_111 - C_222             # Dimension mismatch
@test_throws ErrorException C_111 * C_222             # Dimension mismatch
@test_throws ErrorException [s 0; 1]                  # Dimension mismatch
@test_throws ErrorException D_111 + C_111             # Sampling time mismatch
@test_throws ErrorException D_111 - C_111             # Sampling time mismatch
@test_throws ErrorException D_111 * C_111             # Sampling time mismatch
D_diffTs = tf([1], [2], 0.1)
@test_throws ErrorException D_111 + D_diffTs          # Sampling time mismatch
@test_throws ErrorException D_111 - D_diffTs          # Sampling time mismatch
@test_throws ErrorException D_111 * D_diffTs          # Sampling time mismatch
@test_throws ErrorException tf([1], [2], -0.1)        # Negative samping time
@test_throws ErrorException tf("s", 0.01)             # s creation can't be discrete
@test_throws ErrorException tf("z", 0)                # z creation can't be continuous
@test_throws ErrorException tf("z")                   # z creation can't be continuous
@test_throws ErrorException [z 0]                     # Sampling time mismatch (inferec could be implemented)

# Test polynomial inputs
Polynomial = ControlSystems.Polynomials.Polynomial
v1, v2, v3, v4 = [1,2,3], [4,5,6], [7,8], [9,10,11]
p1, p2, p3, p4 = Polynomial.(reverse.((v1,v2,v3,v4)))

S1v = tf(v1, v2)
S2v = tf(v3, v4)
Sv = [S1v S2v; 2S1v 3S2v]

# Test polynomial SISO constructor
S1p = tf(p1,p2)
# Test polynomial matrix MIMO constructor
Sp = tf([p1 p3; 2p1 3p3], [p2 p4; p2 p4])

@test S1v == S1p
@test Sv == Sp

# Test constructors
@test ControlSystems.SisoRational{Float64}(v1,v2) isa ControlSystems.SisoRational{Float64}
@test ControlSystems.SisoRational{Float64}(v1,1.0*v2) isa ControlSystems.SisoRational{Float64}
@test ControlSystems.SisoRational{Int64}(v1,1.0*v2) isa ControlSystems.SisoRational{Int64}

@test_throws InexactError ControlSystems.SisoRational{Int64}(v1,1.5*v2)

end


text:
export TracerProvider

"""
    TracerProvider(;kw...)

# Keyword Arguments

  - `sampler`=[`DEFAULT_ON`](@ref)
  - `resource`=[`Resource`](@ref)()
  - `span_processor`=[`CompositSpanProcessor`](@ref)()
  - `id_generator`=[`RandomIdGenerator`](@ref)()
  - `limit_info`=[`LimitInfo`](@ref)()
"""
Base.@kwdef struct TracerProvider{
    S<:AbstractSampler,
    R<:Resource,
    SP<:AbstractSpanProcessor,
    RNG,
} <: AbstractTracerProvider
    sampler::S = DEFAULT_ON
    resource::R = Resource()
    span_processor::SP = CompositSpanProcessor()
    id_generator::RNG = RandomIdGenerator()
    limit_info::LimitInfo = LimitInfo()
    is_shut_down::Ref{Bool} = Ref(false)
end

force_flush!(p::TracerProvider) = force_flush!(p.span_processor)

function shut_down!(p::TracerProvider)
    shut_down!(p.span_processor)
    p.is_shut_down[] = true
end

Base.push!(p::TracerProvider, sp::AbstractSpanProcessor) = push!(p.span_processor, sp)
"""
    Span(name::String, tracer::Tracer{TracerProvider};kw...)

# Keyword Arguments

  - `context=current_context()`,
  - `kind=SPAN_KIND_INTERNAL`,
  - `attributes=Dict{String, TAttrVal}()`,
  - `links=Link[]`,
  - `events=Event[]`,
  - `start_time::UInt=time()*10^9`, the nanoseconds.
  - `is_remote=false`
"""
function OpenTelemetryAPI.Span(
    name::String,
    tracer::Tracer{<:TracerProvider};
    context = current_context(),
    kind = SPAN_KIND_INTERNAL,
    attributes = Dict{String,TAttrVal}(),
    links = Link[],
    events = OpenTelemetryAPI.Event[],
    start_time = UInt(time() * 10^9),
    is_remote = false,
)
    provider = tracer.provider
    parent_span_ctx = context |> current_span |> span_context
    if isnothing(parent_span_ctx)
        trace_id = generate_trace_id(provider.id_generator)
    else
        trace_id = parent_span_ctx.trace_id
    end

    attributes = DynamicAttrs(
        attributes;
        count_limit = provider.limit_info.span_attribute_count_limit,
        value_length_limit = provider.limit_info.span_attribute_value_length_limit,
    )

    links = Limited(links; limit = provider.limit_info.span_link_count_limit)
    events = Limited(events; limit = provider.limit_info.span_event_count_limit)

    sampling_result = should_sample(
        provider.sampler,
        context,
        trace_id,
        name,
        kind,
        attributes,
        links,
        isnothing(parent_span_ctx) ? TraceState() : parent_span_ctx.trace_state,
    )

    span_ctx = SpanContext(
        trace_id = trace_id,
        span_id = generate_span_id(provider.id_generator),
        is_remote = is_remote,
        trace_flag = TraceFlag(sampled = is_sampled(sampling_result)),
        trace_state = sampling_result.trace_state,
    )

    is_no_op_span = provider.is_shut_down[] || !is_recording(sampling_result)
    s = Span(
        Ref(name),
        tracer,
        is_no_op_span ? INVALID_SPAN_CONTEXT : span_ctx,
        parent_span_ctx,
        kind,
        start_time,
        Ref{Union{Nothing,UInt64}}(is_no_op_span ? start_time : nothing),
        attributes,
        links,
        events,
        Ref(SpanStatus(SPAN_STATUS_UNSET)),
    )
    if !is_no_op_span
        on_start(tracer.provider.span_processor, s)
    end
    s
end

function OpenTelemetryAPI.end!(s::Span{<:Tracer{<:TracerProvider}}, t = UInt(time() * 10^9))
    if is_recording(s)
        s.end_time[] = t
        on_end(s.tracer.provider.span_processor, s)
    else
        @warn "the span is not recording."
    end
end


text:
using Halton

## Halton Sequence

# First 10 entries in Halton Sequence with base 3
H = Vector{Float64}(10)
B = 3
HaltonSeq!(H, B)

# Skip first 100 entries
HaltonSeq!(H, B, skip=100)

## Halton Draws

# 1000 Normal draws using a Halton sequence with base 17.
# Default: Skip the first 500 draws
H = Vector{Float64}(1000)
B = 17
HaltonDraws!(H, B)

# Check your work:
using Plots
gr()
histogram(H)

# 1000 Normal draws using a Halton sequence with base 17.
## Use Distributions.Uniform and skip no entries
## Requires Distributions.quantile(distr,x)
H = Vector{Float64}(1000)
B = 17
HaltonDraws!(H, B, distr=Distributions.Uniform(), skip=0)


text:
using OutlierDetectionBenchmark
using Test

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


text:
#=== Set everything up for FLOWFarm ===#
cd("/Users/nbaker/Documents/GitHub/FLOWFarm.jl/test/")
include("iea37_specific_functions.jl")
using FLOWFarm; const ff = FLOWFarm
import YAML
using PyPlot
using CSV

#--- Read in windfarm boundary data ---#
# Which case study we're doing. 'cs3' or 'cs4'
str_case = "4"
#- Rip the boundary coordinates from the .yaml file -#
file_dir = "./inputfiles/"
bnry_file_name_orig = "iea37-boundary-cs" * str_case * ".yaml"
bnry_file_name = string(file_dir,bnry_file_name_orig)
bndry_x, bndry_y = getBndryCs4YAML(bnry_file_name)
bndry_x_clsd, bndry_y_clsd = ff.closeBndryLists(bndry_x, bndry_y)

# Place the turbines for the first area
layout_file_name = "./inputfiles/iea37-ex-opt4.yaml"
~, ~, fname_turb, ~ = ff.get_turb_loc_YAML(layout_file_name)
turbine_file_name = string("./inputfiles/",fname_turb)
~, ~, ~, ~, turb_diam, ~ = ff.get_turb_atrbt_YAML(turbine_file_name)
region = 1

# Plot boundary
bndry_x = bndry_x_clsd[region]
bndry_y = bndry_y_clsd[region]

plot(bndry_x, bndry_y)
num_bndry_pts = length(bndry_x)
println(num_bndry_pts)
#scatter(bndry_x[20],bndry_y[20])#, turb_diam/2.0, fill=true,color="black")
# Plot boundary verticies
for i = 1:num_bndry_pts
     scatter(bndry_x[i],bndry_y[i])#, turb_diam/2.0, fill=true,color="black")
     plt.text(bndry_x[i]+turb_diam, bndry_y[i]+turb_diam, string(i))
end

# Formatting
axis("square")
axis("off")
plt.show()


text:
mutable struct LoadedState
    roomName::String
    filename::String
    room::Union{Maple.Room, Nothing}
    side::Union{Maple.Side, Nothing}
    map::Union{Maple.Map, Nothing}
    lastSavedHash::Unsigned
end

function LoadedState(roomName::String, filename::String)
    side = nothing
    map = nothing
    room = nothing

    if isfile(filename)
        side = loadSide(filename)
        map = side.map
        room = getRoomByName(map, roomName)

        if map !== nothing
            EntityIds.updateValidIds(map)
        end
    end

    return LoadedState(roomName, filename, isa(room, Maple.Room) ? room : nothing, side, map, hash(side))
end


text:
using Documenter
using Literate
using AutomotiveDrivingModels

# generate tutorials and how-to guides using Literate
src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "lit")
notebooks = joinpath(src, "notebooks")

for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = true)
    Literate.notebook(ipath, notebooks, execute = false)
end

makedocs(
	modules = [AutomotiveDrivingModels],
	format = Documenter.HTML(),
	sitename="AutomotiveDrivingModels.jl",
	pages=[
		"Home" => "index.md",
		"Tutorials" => [
			"tutorials/straight_roadway.md",
			"tutorials/stadium.md",
			"tutorials/intersection.md",
			"tutorials/crosswalk.md",
			"tutorials/sidewalk.md",
			"tutorials/feature_extraction.md"
		],
		"Manual" => [
			"Roadways.md",
			"actions.md",
			"states.md",
			"agent_definitions.md",
			"behaviors.md",
			"simulation.md",
			"collision_checkers.md",
			"feature_extraction.md"
		],
    "Internals" =>[
			 "Vec.md"
		]
	]
)

deploydocs(
    repo = "github.com/sisl/AutomotiveDrivingModels.jl.git",
)


text:
using TV-VAR
using Test

@testset "TV-VAR.jl" begin
    # Write your tests here.
end


text:
# Copyright 2017,2018 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using StatsBase
include("loaddata.jl")

@testset "encode by farthest points" begin
    using KernelMethods
    X, y = loadiris()
    kmodel = KernelClassifier(X, y, folds=3, ensemble_size=3, size=31, score=accuracy)
    yh = predict(kmodel, X)
    acc = mean(y .== yh)
    @info "===== KernelClassifier accuracy: $acc"
    @test acc > 0.9
end


text:
module tippling_examples

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
geostiffness = FEMMCorotBeamModule.geostiffness
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio, save_to_json

function tippling_1()
    # Parameters:
    E = 1000000.0;
    nu = 0.3;
    L =   30.0; # Length of the beam
    b =  0.5; # width
    h =  4.0; # height
    magn = -0.2056167583560*1e4/L^2;
    magn_scale = 10000
    neigvs = 4;

    # Reference buckling factors: note that there are pairs, positive and
    # negative. This is due to the symmetry of the problem.
    reffs = [48.5475, 124.1839]

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [-1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    n = 32;
    members = []
    push!(members, frame_member([0 0 0; 0 0 L], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, E, nu)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = selectnode(fens; box = Float64[0 0 0 0 0 0], tolerance = L/10000)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the global discrete system
    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    tipn = selectnode(fens; box = Float64[0 0 0 0 L L], tolerance = L/10000)
    loadbdry = FESetP1(reshape(tipn, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[0, -magn*b*h/magn_scale, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve the static problem
    U = K\F
    scattersysvec!(dchi, U[:])
    @show  dchi.values[tipn, :]

    # Compute the geometric stiffness
    u0.values = dchi.values[:,1:3]
    update_rotation_field!(Rfield0, dchi)
    Kg = geostiffness(femm, geom0, u0, Rfield0, dchi);

    # Solve the eigenvalue problem
    d,v,nconv = eigs(-Kg, K; nev=neigvs, which=:LM, explicittransform=:none)
    fs = 1.0 ./ (d) ./ magn_scale
    println("Buckling factors: $fs [ND]")
    println("Reference: $reffs [ND]")

    # Visualize vibration modes
    scattersysvec!(dchi, 10*v[:, 1])
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[-L/2 -L/2 0]; [L/2 L/2 1.1*L]]),
    plot_nodes(fens),
    plot_solid(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    save_to_json(pl, "plot.json")
    # @show p.plot.data
    true
end # tippling_1

function allrun()
    println("#####################################################")
    println("# tippling_1 ")
    tippling_1()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
immutable ClassicalRustModel <: StatisticalModel

	#Linear Rust Model u[k,a]=aa[k,a,itheta]*theta + bb[k,a]
	aa::Array{Float64,3}		  	#aa[k,a,itheta]
	bb::Array{Float64,2}		  	#bb[k,a]

	#Initial distribution
	mu::Array{Float64,2}			#mu[k,a]

	core::RustModelCore
end

rustmodel(aa,bb,mu,beta,pi)=ClassicalRustModel(aa,bb,mu,rustmodelcore(beta,pi))

rustmodel(aa,bb,mu,beta,pi,horizonindex::Int)=ClassicalRustModel(aa,bb,mu,rustmodelcore(beta,pi,horizonindex))

checkdp(model::ClassicalRustModel)=checkdp(model.core)

function coef!(model::ClassicalRustModel,theta)
	dk,da=size(model.bb)
	dtheta=length(theta)
	for ia=1:da
		for ik=1:dk
			model.core.u[ik,ia]=model.bb[ik,ia]
			for itheta=1:dtheta
				model.core.u[ik,ia]+=model.aa[ik,ia,itheta]*theta[itheta]
			end
		end
	end
	solvedp!(model.core)
end

#sample (x,y) from a matrix of joint probabilities.
function wsample2(mu)
	dx,dy=size(mu)
	ind2sub((dx,dy),wsample(1:dx*dy,vec(mu)))
end

#an alternative would be to compute the transition matrix for y=(k,a) and
#draw from the Markov chain
function rand(model::ClassicalRustModel,T::Int)
	dk,da=size(model.bb)
	data=zeros(Int,T,2)
	ka1=wsample2(model.mu)
	data[1,1]=ka1[1]
	data[1,2]=ka1[2]
	for t=2:T
		data[t,1]=wsample(1:dk,reshape(model.core.pi[data[t-1,2]][data[t-1,1],:],dk))
		data[t,2]=wsample(1:da,reshape(model.core.p[data[t,1],:],da))
	end
	data
end

function loglikelihood(model::ClassicalRustModel,data::Array{Int,2})
	T=size(data,1)
	llk=log(model.mu[data[1,1],data[1,2]])
	for t=2:T
		llk+=log(model.core.pi[data[t-1,2]][data[t-1,1],data[t,1]])
		llk+=log(model.core.p[data[t,1],data[t,2]])
	end
	llk/T
end

dim(model::ClassicalRustModel)=size(model.aa,3)


text:
"""
Check that all elements of a vector have the same length.
Return the common length or throw a DimensionMismatch error.
"""
function checklength(atoms::Tuple{Vararg{AbstractAtom}}, msg)
    n = length(atoms[1])
    if !all(a->length(a)==n, atoms)
        throw(DimensionMismatch(msg))
    end
    return n
end

function checklength(sets::Tuple{Vararg{AbstractAtomicSet}}, msg)
    n = length(sets[1])
    if !all(A->length(A)==n, sets)
        throw(DimensionMismatch(msg))
    end
    return n
end

function checklength(faces::Tuple{Vararg{AbstractFace}}, msg)
    n = length(faces[1])
    if !all(F->length(F)==n, faces)
        throw(DimensionMismatch(msg))
    end
    return n
end


text:
const codegen_ctxs = Dict{Symbol,Any}()
codegen_ctxs[:julia] = JuliaASTCodegenCtx

function device_m(__module__::Module, ex::Expr; target::Symbol = :julia)
    ir = YaoIR(__module__, ex) #= default parsing pass =#

    # simple validation
    ir.pure_quantum = is_pure_quantum(ir)
    ir.qasm_compatible = is_qasm_compatible(ir)

    # TODO: code optimization/transformation pass
    # TODO: switch compile target

    @timeit_debug to "codegen" code = codegen(codegen_ctxs[target](ir), ir)
    return code
end


text:
function scan!(op, xs)
    acc = xs[begin]
    for i in firstindex(xs)+1:lastindex(xs)
        acc = op(acc, @inbounds xs[i])
        @inbounds xs[i] = acc
    end
    return acc
end

function scan!(op, acc, xs)
    for i in eachindex(xs)
        acc = op(acc, @inbounds xs[i])
        @inbounds xs[i] = acc
    end
    return acc
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    precompile(Tuple{Type{LoopVectorization.LoopSet},Array{Symbol,1},Array{Int64,1},Array{LoopVectorization.Loop,1},Dict{Symbol,LoopVectorization.Operation},Array{LoopVectorization.Operation,1},Array{Int64,1},Array{Int64,1},LoopVectorization.LoopOrder,Expr,Expr,Array{Tuple{Int64,Symbol},1},Array{Tuple{Int64,Int64},1},Array{Tuple{Int64,Float64},1},Array{Int64,1},Array{Int64,1},Array{Symbol,1},Array{Symbol,1},Array{Symbol,1},Array{LoopVectorization.ArrayReferenceMeta,1},Array{Float64,2},Array{Float64,2},Array{Bool,1},Array{Bool,1},Base.RefValue{LoopVectorization.UnrollSpecification},Base.RefValue{Bool},Base.RefValue{LoopVectorization.LoopStartStopManager},Base.RefValue{Int64},Base.RefValue{Bool},Symbol})
    precompile(Tuple{Type{LoopVectorization.LoopSet},Expr,Symbol})
    precompile(Tuple{Type{LoopVectorization.Loop},LoopVectorization.LoopSet,Expr,Symbol,Type{VectorizationBase.StaticLowerUnitRange{1}}})
    precompile(Tuple{typeof(LoopVectorization._avx_loopset),Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Core.SimpleVector,Any})
    precompile(Tuple{typeof(LoopVectorization.add_compute!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,LoopVectorization.ArrayReferenceMetaPosition})
    precompile(Tuple{typeof(LoopVectorization.add_compute!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,Nothing})
    precompile(Tuple{typeof(LoopVectorization.add_constant!),LoopVectorization.LoopSet,Float64,Array{Symbol,1},Symbol,Int64})
    precompile(Tuple{typeof(LoopVectorization.add_constant_offset_load_elmination_cost!),SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},Base.RefValue{Bool},LoopVectorization.LoopSet,LoopVectorization.Operation,Float64,LoopVectorization.UnrollSymbols,Bool,Bool,Int64,Int64,Bool})
    precompile(Tuple{typeof(LoopVectorization.add_if!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,LoopVectorization.ArrayReferenceMetaPosition})
    precompile(Tuple{typeof(LoopVectorization.add_if!),LoopVectorization.LoopSet,Symbol,Expr,Int64,Int64,Nothing})
    precompile(Tuple{typeof(LoopVectorization.add_mref!),LoopVectorization.LoopSet,LoopVectorization.ArrayReferenceMeta,Int64,Type{VectorizationBase.PermutedDimsStridedPointer{(2, 1),(2, 1),Float64,VectorizationBase.PackedStridedPointer{Float64,1}}},Symbol})
    precompile(Tuple{typeof(LoopVectorization.add_pow!),LoopVectorization.LoopSet,Symbol,Symbol,Int64,Int64,Int64})
    precompile(Tuple{typeof(LoopVectorization.argmeta_and_consts_description),LoopVectorization.LoopSet,Array{Symbol,1}})
    precompile(Tuple{typeof(LoopVectorization.array_reference_meta!),LoopVectorization.LoopSet,Symbol,SubArray{Any,1,Array{Any,1},Tuple{UnitRange{Int64}},true},Int64,Nothing})
    precompile(Tuple{typeof(LoopVectorization.choose_tile),LoopVectorization.LoopSet})
    precompile(Tuple{typeof(LoopVectorization.cost_vec_buf),LoopVectorization.LoopSet})
    precompile(Tuple{typeof(LoopVectorization.determine_unroll_factor),LoopVectorization.LoopSet,Array{Symbol,1},Symbol,Symbol})
    precompile(Tuple{typeof(LoopVectorization.evaluate_cost_tile),LoopVectorization.LoopSet,Array{Symbol,1},LoopVectorization.UnrollSymbols})
    precompile(Tuple{typeof(LoopVectorization.evaluate_cost_unroll),LoopVectorization.LoopSet,Array{Symbol,1},Symbol,Float64})
    precompile(Tuple{typeof(LoopVectorization.generate_call),LoopVectorization.LoopSet,Tuple{Int8,Int8,Int8},Bool})
    precompile(Tuple{typeof(LoopVectorization.pushop!),LoopVectorization.LoopSet,LoopVectorization.Operation,Symbol})
    precompile(Tuple{typeof(LoopVectorization.stride_penalty),LoopVectorization.LoopSet,LoopVectorization.Operation,Array{Symbol,1},Array{Int64,1}})
    precompile(Tuple{typeof(LoopVectorization.substitute_broadcast),Expr,Symbol})
    precompile(Tuple{typeof(println),Base.GenericIOBuffer{Array{UInt8,1}},Array{LoopVectorization.Operation,1}})
    precompile(Tuple{typeof(push!),LoopVectorization.LoopSet,Expr,Int64,Int64})
end


text:
function generate_inputdata1d!(X, index_output, network::NNPP, f, data, element_id, mesh, c2e)
    a = data[1]
    b = data[end]
    a1 = a[1]
    b1 = b[1]
    nnodes = length(data)
    n,_ = gauss_lobatto_nodes_weights(nnodes)
    _, inverse_vandermonde_legendre = vandermonde_legendre(n)

    modal = zeros(1,nnodes)
    indicator= zeros(1,nnodes)

    for i in 1:nnodes
        indicator[1,i] = f((b1-a1)/2*n[i]+(a1+b1)/2)
    end

    multiply_dimensionwise!(modal, inverse_vandermonde_legendre, indicator)

    total_energy = 0.0
    for i in 1:nnodes
        total_energy += modal[1, i]^2
    end
    total_energy_clip1 = 0.0
    for i in 1:(nnodes-1)
        total_energy_clip1 += modal[1, i]^2
    end
    total_energy_clip2 = 0.0
    for i in 1:(nnodes-2)
        total_energy_clip2 += modal[1, i]^2
    end

    x1 = (total_energy - total_energy_clip1)/total_energy
    x2 = (total_energy_clip1 - total_energy_clip2)/total_energy_clip1

    # directly normalize
    norm = max(maximum(abs.([x1,x2])),1)
    X[1,index_output] = x1 /norm
    X[2,index_output] = x2 /norm

end

function generate_inputdata1d!(X, index_output, network::NNPPh, f, data, element_id, mesh, c2e)
    # Todo; unterschied von a und a1???
    a = data[1]
    b = data[end]
    a1 = a[1]
    b1 = b[1]
    nnodes = length(data)
    n,_ = gauss_lobatto_nodes_weights(nnodes)
    _, inverse_vandermonde_legendre = vandermonde_legendre(n)

    modal = zeros(1,nnodes)
    indicator= zeros(1,nnodes)

    for i in 1:nnodes
        indicator[1,i] = f((b1-a1)/2*n[i]+(a1+b1)/2)
    end

    multiply_dimensionwise!(modal, inverse_vandermonde_legendre, indicator)

    total_energy = 0.0
    for i in 1:nnodes
        total_energy += modal[1, i]^2
    end
    total_energy_clip1 = 0.0
    for i in 1:(nnodes-1)
        total_energy_clip1 += modal[1, i]^2
    end
    total_energy_clip2 = 0.0
    for i in 1:(nnodes-2)
        total_energy_clip2 += modal[1, i]^2
    end

    x1 = (total_energy - total_energy_clip1)/total_energy
    x2 = (total_energy_clip1 - total_energy_clip2)/total_energy_clip1

    # directly normalize
    norm = max(maximum(abs.([x1,x2,nnodes])),1)
    X[1,index_output] = x1 / norm
    X[2,index_output] = x2 / norm
    X[3,index_output] = nnodes / norm

end

function generate_inputdata1d!(X, index_output, network::NNRH, f, data, element_id, mesh, c2e)
    a = data[1]
    b = data[end]
    leaf_cell_ids = leaf_cells(mesh.tree)
    nnodes = length(data)
    indicator = zeros(nnodes)
    n,_ = gauss_lobatto_nodes_weights(nnodes)

    # calc the interface values for the central element
    for i in 1:nnodes
        indicator[i] = f((b-a)/2*n[i]+(a+b)/2)
    end

    x_a = (b-a)/2*n[1]+(a+b)/2
    t_a = x_a*2/(b-a)-(a+b)/(b-a)

    x_b = (b-a)/2*n[end]+(a+b)/2
    t_b = x_b*2/(b-a)-(a+b)/(b-a)
    interface_value_right = 0
    interface_value_left = 0
    for i = 1:nnodes
        phi_a,_= legendre_polynomial_and_derivative(i-1,t_a)./sqrt(i-0.5)
        phi_b,_= legendre_polynomial_and_derivative(i-1,t_b)./sqrt(i-0.5)
        interface_value_left += indicator[i]*phi_a
        interface_value_right += indicator[i]*phi_b
    end

    # now calc all mean values
    # get neighbors id
    cell_id = leaf_cell_ids[element_id]

    # Todo: Allokation hier ist schlecht
    neighbor_ids = Array{Int64}(undef, 2)
    c2e[leaf_cell_ids[element_id]] = element_id

    for direction in 1:n_directions(mesh.tree)
        if has_neighbor(mesh.tree, cell_id, direction)
            neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
            if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor
                if direction == 1
                    neighbor_ids[direction] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
                else
                    neighbor_ids[direction] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
                end
            else # Cell has same refinement level neighbor
                neighbor_ids[direction] = c2e[neighbor_cell_id]
            end
        else # Cell is small and has large neighbor
            parent_id = mesh.tree.parent_ids[cell_id]
            neighbor_cell_id = mesh.tree.neighbor_ids[direction, parent_id]
            neighbor_ids[direction] = c2e[neighbor_cell_id]
        end
    end

    u_mean = get_mean(f, data)

    get_data!(data, leaf_cell_ids, mesh, neighbor_ids[1], nnodes)
    ul_mean = get_mean(f, data)
    get_data!(data, leaf_cell_ids, mesh, neighbor_ids[2], nnodes)
    ur_mean = get_mean(f, data)

    # Write inputdata in X and directly normalize
    norm = max(maximum(abs.([ul_mean, u_mean, ur_mean, interface_value_left, interface_value_right ])),1)
    X[1,index_output] = ul_mean / norm
    X[2,index_output] = u_mean / norm
    X[3,index_output] = ur_mean / norm
    X[4,index_output] = interface_value_left / norm
    X[5,index_output] = interface_value_right / norm
end

function get_mean(f, data)
    a = data[1]
    b = data[end]

    nnodes = length(data)
    indicator = zeros(nnodes)
    nab = zeros(nnodes)
    n,_ = gauss_lobatto_nodes_weights(nnodes)
    for i in 1:nnodes
        nab[i] = (b-a)/2*n[i]+(a+b)/2
        indicator[i] = f(nab[i])
    end

    mean = 0
    for x in nab
        t = x*2/(b-a)-(a+b)/(b-a)
        for i = 1:nnodes
            phi,_= legendre_polynomial_and_derivative(i-1,t)./sqrt(i-0.5)
            mean += indicator[i]*phi
        end
    end

    return mean / nnodes
end

function get_data!(data, leaf_cell_ids, mesh, element_id, nnodes)
    cell_id = leaf_cell_ids[element_id]
    dx = length_at_cell(mesh.tree, cell_id)
    nodes,_ = gauss_lobatto_nodes_weights(nnodes)
    for i in 1:nnodes
        data[i] = (mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
    end
end


text:
"""
Author: Eric Ford
Adapted from: https://github.com/eford/RvSpectraKitLearn.jl/blob/master/src/deriv_spectra_gp.jl

GPs using Stheno.jl, TemporalGPs., KernelFunctions.jl, etc.
"""

"""
Module for interpolating via Gaussian Process Regression based on Stheno and TemporalGPs packages.
"""
module TemporalGPInterpolation
using LinearAlgebra
using PDMats
using StaticArrays
using TemporalGPs, KernelFunctions
using Dates
using RvSpectMLBase

#export make_kernel_data, make_kernel_obs_pred,
export gp_marginal
export predict_mean, predict_deriv, predict_deriv2, predict_mean_and_deriv, predict_mean_and_derivs

# import Stheno: AbstractGP
# instead of Stheno, import temporal GPs functions to replace Stheno functions
import TemporalGPs: GP as AbstractGP
import TemporalGPs: marginals
import TemporalGPs: GP
import Distributions: AbstractMvNormal

using Distributions
const Fx_PosteriorType = Distribution{Multivariate,Continuous}

import ..numerical_deriv, ..calc_dfluxdlnlambda, ..calc_dfluxdlnlambda!, ..calc_d2fluxdlnlambda2, ..calc_d2fluxdlnlambda2!
include("low_level.jl")
export  construct_gp_prior, construct_gp_posterior, log_pdf_gp_posterior

include("convenience.jl")
export interp_chunk_to_shifted_grid_gp_temporal, interp_chunk_to_shifted_grid_gp_temporal!
export interp_chunk_to_grid_gp_temporal, interp_chunk_to_grid_gp_temporal!

end # module


text:
using Distributions, OrdinaryDiffEq, GpABC, Distances, PyPlot

# GaussianProcesses/Bioinformatics paper/l2_norm_two_params/...

idx1 = 1
idx2 = 9
n_design_point = 200
n_plot = 100
contour_colors = ["white", "#FFE9EC", "#FFBBC5", "#FF8B9C", "#FF5D75", "#FF2F4E", "#D0001F", "#A20018", "#990017", "#800013"]

true_params =  [2.0, 2.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.),
            Uniform(75., 125.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)]
param_indices = [idx1, idx2]
Tspan = (0.0, 10.0)
x0 = [3.0, 2.0, 1.0]
solver = RK4()
saveat = 0.1

function ODE_3GeneReg(dx, x, pars, t)
    dx[1] = pars[1]/(1+pars[7]*x[3]) - pars[4]*x[1]
    dx[2] = pars[2]*pars[8]*x[1]./(1+pars[8]*x[1]) - pars[5]*x[2]
    dx[3] = pars[3]*pars[9]*x[1]*pars[10]*x[2]./(1+pars[9]*x[1])./(1+pars[10]*x[2]) - pars[6]*x[3]
end

function simulate_particle(param1::T, param2::T) where {T<:Real}
    params = copy(true_params)
    params[idx1] = param1
    params[idx2] = param2

    prob = ODEProblem(ODE_3GeneReg, x0, Tspan, params)
    Obs = solve(prob, solver, saveat=saveat)

    return Obs
end

function simulate_particles_product(param1::AbstractArray{T, 1}, param2::AbstractArray{T, 1}) where {T<:Real}
    map(x->simulate_particle(x...), Iterators.product(param1, param2))
end

function emulate_particles_product(param1::AbstractArray{T, 1}, param2::AbstractArray{T, 1}, gpem::GPModel) where {T<:Real}
    test_x = vcat(map(x->[x...]',
        reshape(collect(Iterators.product(param1, param2)),
            (length(param1) * length(param2),) ))...)
    ret = gp_regression(test_x, gpem)
    return reshape(ret[1], (length(param1), length(param2)))
end

function simulate_particles(params::AbstractArray{T, 2}) where {T<:Real}
    map(i->simulate_particle(params[i, :]...), range(1, size(params, 1)))
end

function build_distances(model_results, ref_data::AbstractArray{T,2},
    summary_stats_func::Function=GpABC.keep_all_summary_statistic,
    distance_func::Function=Distances.euclidean) where {T<:Real}
    ref_stats = summary_stats_func(ref_data)
    return map(x->distance_func(x, ref_stats), summary_stats_func.(model_results))
end

reference_data = simulate_particle(true_params[idx1], true_params[idx2])
contour_x = linspace(params(priors[idx1])..., n_plot)
contour_y = linspace(params(priors[idx2])..., n_plot)
sim_result = build_distances(
    simulate_particles_product(contour_x, contour_y),
    reference_data)

training_x = hcat(rand(priors[idx1], n_design_point), rand(priors[idx2], n_design_point))
training_y = build_distances(simulate_particles(training_x), reference_data)
gpm = GPModel(training_x=training_x, training_y=training_y, kernel=SquaredExponentialArdKernel())
gp_train(gpm)
em_result = emulate_particles_product(contour_x, contour_y, gpm)

subplot(121)
subplots_adjust(wspace  =  0.5)
contourf(contour_x, contour_y, sim_result, 8, colors=contour_colors)
# scatter3D(repmat(contour_x, 1, n_plot), repmat(contour_y', n_plot, 1), sim_result)
title("L2 norm - Simulation")
xlabel("Parameter $(idx1)")
ylabel("Parameter $(idx2)")
subplot(122)
contourf(contour_x, contour_y, em_result, 8, colors=contour_colors)
# scatter3D(repmat(contour_x, 1, n_plot), repmat(contour_y', n_plot, 1), em_result)
scatter(training_x[:, 1], training_x[:,2], marker="x", label="training points")
legend()
xlim(params(priors[idx1])...)
ylim(params(priors[idx2])...)
title("L2 norm - Emulation")
xlabel("Parameter $(idx1)")
ylabel("Parameter $(idx2)")
# x = simulate_prior_distances(5, 6, reference_data)


text:
using celerite
using Test
using LinearAlgebra
using Profile
using Random

include("test_solver.jl")
include("test_cholesky.jl")


text:
# ---
# title: 1106. Parsing A Boolean Expression
# id: problem1106
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/parsing-a-boolean-expression/description/>
# hidden: true
# ---
#
# Return the result of evaluating a given boolean `expression`, represented as a
# string.
#
# An expression can either be:
#
#   * `"t"`, evaluating to `True`;
#   * `"f"`, evaluating to `False`;
#   * `"!(expr)"`, evaluating to the logical NOT of the inner expression `expr`;
#   * `"&(expr1,expr2,...)"`, evaluating to the logical AND of 2 or more inner expressions `expr1, expr2, ...`;
#   * `"|(expr1,expr2,...)"`, evaluating to the logical OR of 2 or more inner expressions `expr1, expr2, ...`
#
#
#
# **Example 1:**
#
#
#
#     Input: expression = "!(f)"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: expression = "|(f,t)"
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: expression = "&(t,f)"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: expression = "|(&(t,f,t),!(t))"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= expression.length <= 20000`
#   * `expression[i]` consists of characters in `{'(', ')', '&', '|', '!', 't', 'f', ','}`.
#   * `expression` is a valid expression representing a boolean, as given in the description.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import HoTRG_lemon.LatticeModule:
	getHausdorffDim,
	getFractalDim,
	getLegextension,
	getTensorP,
	getTensorQ,
	setTensorP!,
	setTensorQ!

FractalSimulators = [
					 :Classical2dFractalSimulator,
					 :Quantum2dFractalSimulator,
					 :Quantum2dFractalInititerSimulator
					]

#
##### update the SimulatorModule

for FractalSimulator in FractalSimulators

	eval(quote

	  function setNorm!{T}(simulator::$(FractalSimulator){T}, norm::T, which::Int, iteration::Int)
			simulator.normalizationFactor[iteration, which] = norm
	  end

		function getNorm{T}(simulator::$(FractalSimulator){T}, which::Int, iteration::Int)
			return simulator.normalizationFactor[iteration, which]
		end

		function setCoefficient!{T}(simulator::$(FractalSimulator){T}, coefficient::T, which::Int, iteration::Int)
			simulator.coefficients[iteration, which] = coefficient
		end

		function getCoefficient{T}(simulator::$(FractalSimulator){T}, which::Int, iteration::Int)
			if which == 0
				coef = simulator.coefficients[iteration,:]
			else
				coef = simulator.coefficients[iteration, which]
			end
			return coef
		end

	end)

		# inherit from FractalLattice
		single_args = [
					   :getHausdorffDim,
					   :getFractalDim,
					   :getLegextension,
					   :getTensorQ
					   ]

		for f in single_args
			eval(quote

			  function $(f)(simulator::$(FractalSimulator))
			  return $(f)(simulator.lattice)
				end

			end)
		end

		getTensors = [
					 :getTensorP
					 ]

		for f in getTensors
			eval(quote

			  function $(f)(simulator::$(FractalSimulator), which::Int = 0)
				  return $(f)(simulator.lattice, which)
			  end

			 end)
		end

		setTensors = [
					  :setTensorP!
					  ]
		for f in setTensors
			eval(quote

			  function $(f)(simulator::$(FractalSimulator), which::Int, tensor::Array)
			  return $(f)(simulator.lattice, which, tensor)
			  end

			 end)
		end

		setTensors2 = [
					  :setTensorQ!
					 ]
		for f in setTensors2
			eval(quote

			  function $(f)(simulator::$(FractalSimulator), tensor::Array)
			  return $(f)(simulator.lattice, tensor)
			  end

			 end)
		end
end


text:
using BenchmarkTools
using Random
using OMEinsum

Random.seed!(0)

const SUITE = BenchmarkGroup()

SUITE["manyinds"] = BenchmarkGroup()
suite = SUITE["manyinds"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["abcdefghijklmnop,flnqrcipstujvgamdwxyz->bcdeghkmnopqrstuvwxyz"]
    dims = [2]
    for str in args, dim in dims
        code = ein"abcdefghijklmnop,flnqrcipstujvgamdwxyz->bcdeghkmnopqrstuvwxyz"
        arr1 = rand(T, map(i->2, OMEinsum.getixs(code)[1])...)
        arr2 = rand(T, map(i->2, OMEinsum.getixs(code)[2])...)
        suite[string(T)]["tiny"] = @benchmarkable code($arr1,$arr2)
    end
end

# Matrix multiplication
SUITE["matmul"] = BenchmarkGroup()
suite = SUITE["matmul"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["tiny"   => rand(T,2,2),
    "small"  => rand(T,10,10),
    "medium" => rand(T,10^2,10^2),
    "large"  => rand(T, 10^3,10^3)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ij,jk -> ik"($m,$m)
    end
end

# Matrix batch-multiplication
SUITE["batchmul"] = BenchmarkGroup()
suite = SUITE["batchmul"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["tiny"   => rand(T,2,2,3),
    "small"  => rand(T,10,10,3),
    "medium" => rand(T,10^2,10^2,3),
    "large"  => rand(T, 10^3,10^3,3)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijk,jlk -> ilk"($m,$m)
    end
end

#inner - reduction to scalar
SUITE["dot"] = BenchmarkGroup()
suite = SUITE["dot"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["tiny"   => rand(T,fill(2,3)...)
    "small"  => rand(T,fill(10,3)...)
    "medium" => rand(T,fill(20,3)...)
    "large"  => rand(T,fill(50,3)...)
    "huge"   => rand(T,fill(10^2,3)...)]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijk,ijk -> "($m, $m)
    end
end

#trace
SUITE["trace"] = BenchmarkGroup()
suite = SUITE["trace"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["tiny"   => rand(T,fill(2,2)...)
    "small"  => rand(T,fill(10,2)...)
    "medium" => rand(T,fill(10^2,2)...)
    "large"  => rand(T,fill(10^3,2)...)]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ii -> "($m)
    end
end

#partial trace
SUITE["ptrace"] = BenchmarkGroup()
suite = SUITE["ptrace"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args = ["tiny"   => rand(T,fill(2,3)...)
    "small"  => rand(T,fill(10,3)...)
    "medium" => rand(T,fill(20,3)...)
    "large"  => rand(T,fill(50,3)...)
    "huge"   => rand(T,fill(10^2,3)...)]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"iij -> j"($m)
    end
end

#diagonal
SUITE["diag"] = BenchmarkGroup()
suite = SUITE["diag"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,3)...)
    "small"  => rand(T,fill(10,3)...)
    "medium" => rand(T,fill(20,3)...)
    "large" => rand(T,fill(30,3)...)
    "huge" => rand(T,fill(100,3)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijj -> ij"($m)
    end
end

#permutation
SUITE["perm"] = BenchmarkGroup()
suite = SUITE["perm"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(1,4)...)
    "small"  => rand(T,fill(2,4)...)
    "medium" => rand(T,fill(10,4)...)
    "large"  => rand(T,fill(30,4)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijkl -> ljki"($m)
    end
end

# tensor contraction
SUITE["tcontract"] = BenchmarkGroup()
suite = SUITE["tcontract"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(1,3)...)
    "small"  => rand(T,fill(2,3)...)
    "medium" => rand(T,fill(10,3)...)
    "large"  => rand(T,fill(30,3)...)
    ]

    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijk, jlk -> il"($m,$m)
    end
end

# star contraction
SUITE["star"] = BenchmarkGroup()
suite = SUITE["star"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,2)...)
    "small"  => rand(T,fill(10,2)...)
    "medium" => rand(T,fill(30,2)...)
    "large"  => rand(T,fill(100,2)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ij,ik,il -> jkl"($m,$m,$m)
    end
end

#star and contract
SUITE["starandcontract"] = BenchmarkGroup()
suite = SUITE["starandcontract"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,2)...)
    "small"  => rand(T,fill(10,2)...)
    "medium" => rand(T,fill(30,2)...)
    "large"  => rand(T,fill(100,2)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ij,ik,ik -> j"($m,$m,$m)
    end
end

# index-sum
SUITE["indexsum"] = BenchmarkGroup()
suite = SUITE["indexsum"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,3)...)
    "small"  => rand(T,fill(10,3)...)
    "medium" => rand(T,fill(30,3)...)
    "large"  => rand(T,fill(100,3)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijk -> ik"($m)
    end
end
# Hadamard
SUITE["hadamard"] = BenchmarkGroup()
suite = SUITE["hadamard"]
for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,3)...)
    "small"  => rand(T,fill(10,3)...)
    "medium" => rand(T,fill(30,3)...)
    "large"  => rand(T,fill(100,3)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ijk,ijk -> ijk"($m,$m)
    end
end

# Outer
SUITE["outer"] = BenchmarkGroup()
suite = SUITE["outer"]

for T in (Float32, Float64, ComplexF32, ComplexF64)
    suite[string(T)] = BenchmarkGroup()
    args =  [
    "tiny"   => rand(T,fill(2,2)...)
    "small"  => rand(T,fill(10,2)...)
    "medium" => rand(T,fill(50,2)...)
    "large"  => rand(T,fill(100,2)...)
    ]
    for (k,m) in args
        suite[string(T)][k] = @benchmarkable ein"ij,kl -> ijkl"($m,$m)
    end
end


text:
#!/usr bin/env julia

using FASTX
using BioSequences
using CSV
using DataFrames

#function to subset a text AA sequence file (txt) based on repeat domain index txt file
#file: input text file with AA sequence
#index_file: txt file with start and end AA numbers for each domain, one per line
#output: output file to save fasta file of repeats
#prot_name: text of protein name to put into the fasta file as description
#alpha is default, identifies repeats with unique letter, false identifies with number
function subset_sequence(file::String, index_file::String, output::String, prot_name::String, alpha::Bool = true)
    #import sequence from text file
    prot_txt_string = read(file, String)

    #convert to BioSequence
    prot_aa_sequence = LongSequence{AminoAcidAlphabet}(prot_txt_string)

    #import the indexes of each repeat domain
    index_df = CSV.read(index_file, header = 0)

    n_repeat = nrow(index_df)
    println("Number of Repeats:")
    println(n_repeat)

    #create matrix of empty aa BioSequences, length of the number of repeats
    repeat_matrix = fill(aa"", n_repeat)

    #put each repeat into a new position in the array
    for i in 1:n_repeat
        global repeat_matrix[i,1] = prot_aa_sequence[Int(index_df[i,1]):Int(index_df[i,2])]
    end

    #count the length of each repeat
    length_array = fill(1, n_repeat)

    for i in 1:n_repeat
        global length_array[i] = length(repeat_matrix[i,1])
    end

    #create empty string matrix
    title_matrix = fill("", n_repeat)

    #fill title matrix with repeat numbers
    for i in 1:n_repeat
        global title_matrix[i] = string(prot_name, "_Repeat_Number", i)
    end

    #check whether identifier is number or letter, create alphabet matrix:
    if alpha == true
        alpha_matrix = letter_list(n_repeat)
    end

    #Write each aa repeat sequence to a record in a fasta file, either with alphabet or number
    if alpha == true
        open(FASTA.Writer, output) do writer
            for i in 1:n_repeat
                write(writer, FASTA.Record(string(alpha_matrix[i], " ", title_matrix[i]), string(length_array[i], " Residues"), repeat_matrix[i,1]))
            end
        end
    else
        open(FASTA.Writer, output) do writer
            for i in 1:n_repeat
                write(writer, FASTA.Record(string(i, " ", title_matrix[i]), string(length_array[i], " Residues"), repeat_matrix[i,1]))
            end
        end
    end

    println("Individual repeat sequences saved to fasta file:")
    println(output)
end

#function to convert fasta file to text file
function fasta_to_text(fasta_file::String, output_file::String)

    #read the fasta file
    r = open(FASTA.Reader, fasta_file)

    #store sequence as Biosequence variable
    for record in r
	    global seq = sequence(record)
    end

    #convert to string
    seq_str = convert(String, seq)

    open(output_file, "w") do f
	    write(f, seq_str)
    end
    println(string("saved sequence as ", output_file))
end