text:
@testset "SMatrix" begin
    @testset "Inner Constructors" begin
        @test SMatrix{1,1,Int,1}((1,)).data === (1,)
        @test SMatrix{1,1,Float64,1}((1,)).data === (1.0,)
        @test SMatrix{2,2,Float64,4}((1, 1.0, 1, 1)).data === (1.0, 1.0, 1.0, 1.0)

        # Bad input
        @test_throws Exception SMatrix{1,1,Int,1}()
        @test_throws Exception SMatrix{2,1,Int,2}((1,))
        @test_throws Exception SMatrix{1,1,Int,1}(())

        # Bad parameters
        @test_throws Exception SMatrix{1,1,Int,2}((1,))
        @test_throws Exception SMatrix{1,1,1,1}((1,))
        @test_throws Exception SMatrix{1,2,Int,1}((1,))
        @test_throws Exception SMatrix{2,1,Int,1}((1,))
    end

    @testset "Outer constructors and macro" begin
        @test SMatrix{1,1,Int}((1,)).data === (1,)
        @test SMatrix{1,1}((1,)).data === (1,)
        @test SMatrix{1}((1,)).data === (1,)

        @test SMatrix{2,2,Int}((1,2,3,4)).data === (1,2,3,4)
        @test SMatrix{2,2}((1,2,3,4)).data === (1,2,3,4)
        @test SMatrix{2}((1,2,3,4)).data === (1,2,3,4)

        @test ((@SMatrix [1.0])::SMatrix{1,1}).data === (1.0,)
        @test ((@SMatrix [1 2])::SMatrix{1,2}).data === (1, 2)
        @test ((@SMatrix [1 ; 2])::SMatrix{2,1}).data === (1, 2)
        @test ((@SMatrix [1 2 ; 3 4])::SMatrix{2,2}).data === (1, 3, 2, 4)

        @test ((@SMatrix Int[1.0])::SMatrix{1,1}).data === (1,)
        @test ((@SMatrix Float64[1 2])::SMatrix{1,2}).data === (1.0, 2.0)
        @test ((@SMatrix Float64[1 ; 2])::SMatrix{2,1}).data === (1.0, 2.0)
        @test ((@SMatrix Float64[1 2 ; 3 4])::SMatrix{2,2}).data === (1.0, 3.0, 2.0, 4.0)

        @test ((@SMatrix [i*j for i = 1:2, j=2:3])::SMatrix{2,2}).data === (2, 4, 3, 6)
        @test ((@SMatrix Float64[i*j for i = 1:2, j=2:3])::SMatrix{2,2}).data === (2.0, 4.0, 3.0, 6.0)

        @test (ex = macroexpand(:(@SMatrix [1 2; 3])); isa(ex, Expr) && ex.head == :error)

        @test ((@SMatrix zeros(2,2))::SMatrix{2, 2, Float64}).data === (0.0, 0.0, 0.0, 0.0)
        @test ((@SMatrix ones(2,2))::SMatrix{2, 2, Float64}).data === (1.0, 1.0, 1.0, 1.0)
        @test ((@SMatrix eye(2))::SMatrix{2, 2, Float64}).data === (1.0, 0.0, 0.0, 1.0)
        @test ((@SMatrix eye(2,2))::SMatrix{2, 2, Float64}).data === (1.0, 0.0, 0.0, 1.0)
        @test isa(@SMatrix(rand(2,2)), SMatrix{2, 2, Float64})
        @test isa(@SMatrix(randn(2,2)), SMatrix{2, 2, Float64})

        @test ((@SMatrix zeros(Float32, 2, 2))::SMatrix{2,2,Float32}).data === (0.0f0, 0.0f0, 0.0f0, 0.0f0)
        @test ((@SMatrix ones(Float32, 2, 2))::SMatrix{2,2,Float32}).data === (1.0f0, 1.0f0, 1.0f0, 1.0f0)
        @test ((@SMatrix eye(Float32, 2))::SMatrix{2, 2, Float32}).data === (1.0f0, 0.0f0, 0.0f0, 1.0f0)
        @test ((@SMatrix eye(Float32, 2, 2))::SMatrix{2, 2, Float32}).data === (1.0f0, 0.0f0, 0.0f0, 1.0f0)
        @test isa(@SMatrix(rand(Float32, 2, 2)), SMatrix{2, 2, Float32})
        @test isa(@SMatrix(randn(Float32, 2, 2)), SMatrix{2, 2, Float32})
    end

    @testset "Methods" begin
        m = @SMatrix [11 13; 12 14]

        @test isimmutable(m) == true

        @test m[1] === 11
        @test m[2] === 12
        @test m[3] === 13
        @test m[4] === 14

        @test Tuple(m) === (11, 12, 13, 14)

        @test size(m) === (2, 2)
        @test size(typeof(m)) === (2, 2)
        @test size(SMatrix{2,2}) === (2, 2)

        @test size(m, 1) === 2
        @test size(m, 2) === 2
        @test size(typeof(m), 1) === 2
        @test size(typeof(m), 2) === 2

        @test length(m) === 4

        @test_throws Exception m[1] = 1
    end
end


text:
push!(LOAD_PATH, dirname(@__FILE__))
@everywhere log(1.0) # Workaround for https://github.com/JuliaLang/julia/issues/18465
reload("Analysis")
Analysis.compute_permpval()
Analysis.compute_permpval_eachisi()
Analysis.compute_latency()
Analysis.compute_alignedinfoisi()
Analysis.compute_permpval_isiwin()
Analysis.compute_permpval_maintenance()
Analysis.compute_modulation_maintenance()
Analysis.compute_modulation_presentation()
Analysis.compute_maintenance_topstim_holdout()
Analysis.compute_maintenance_topstim_correct()
Analysis.compute_alignedinfoisi(pos1=true)

ext: jl

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
@testset "Test How to Stop II" begin
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
# We illustrate here the features of Stopping when the algorithm is used a
# subStopping.
#
###############################################################################
#using Test, Stopping

#Assume we want to solve "pb" starting from "x0" and solving at each step
#of the algorithm the subproblem "subpb".
# We can use this additional info to improve the stopping criterion.
x0 = ones(2)
pb = nothing
subpb = nothing
subsubpb = nothing

###############################################################################
#Initialize a Stopping for the main pb
main_stop = GenericStopping(pb, x0)
#We can then, initialize another stopping to the subproblem, and providing
#the main_stop as a keyword argument:
sub_stop = GenericStopping(subpb, x0, main_stp = main_stop, tol_check = (atol, rtol, opt0) -> atol)
#Note that by default main_stp is void
@test main_stop.main_stp == VoidStopping()

#The only difference appears in the event of a call to stop!, which now also
#check the time and resources of the main_pb.
OK = start!(sub_stop)
@test OK == false #no reason to stop just yet.

#Assume time is exhausted for the main_stop
main_stop.meta.start_time = 0.0 #force a timing failure in the main problem
stop!(sub_stop)

@test status(sub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
@test sub_stop.meta.tired == false
@test sub_stop.meta.main_pb == true

#The same applies if there is now a third subproblem
reinit!(main_stop)
reinit!(sub_stop)
subsub_stop = GenericStopping(subsubpb, x0, main_stp = sub_stop, tol_check = (atol, rtol, opt0) -> atol)
main_stop.meta.start_time = 0.0 #force a timing failure in the main problem
stop!(subsub_stop)

@test status(subsub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
@test subsub_stop.meta.tired   == false
@test subsub_stop.meta.main_pb == true
@test status(sub_stop, list = true) == [:ResourcesOfMainProblemExhausted]
@test sub_stop.meta.tired   == false
@test sub_stop.meta.main_pb == true
@test status(main_stop, list = true) == [:TimeLimit]
@test main_stop.meta.tired   == true
@test main_stop.meta.main_pb == false

end


text:
using Documenter, HypothesisTests

makedocs(
    format = :html,
    sitename = "HypothesisTests.jl",
    modules = [HypothesisTests],
    pages = [
        "index.md",
        "methods.md",
        "parametric.md",
        "nonparametric.md",
        "time_series.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaStats/HypothesisTests.jl.git",
    target = "build",
    julia  = "1.0",
    deps = nothing,
    make = nothing
)


text:
"""
    AggregationStack{T <: Tuple{Vararg{AbstractAggregation}}} <: AbstractAggregation

A container that implements a concatenation of one or more `AbstractAggregation`s.

Construct with e.g. `mean_aggregation([t::Type, ]d)`, `max_aggregation([t::Type, ]d)` for single
operators and with e.g. `pnormlse_aggregation([t::Type, ]d)` for concatenations. With these calls
all parameters inside operators are initialized randomly as `Float32` arrays, unless type `t` is
further specified. It is also possible to call the constructor directly, see Examples.

Intended to be used as a functor:

    (a::AggregationStack)(x, bags[, w])

where `x` is either `Missing`, `AbstractMatrix` or [`ArrayNode`](@ref),
`bags` is [`AbstractBags`](@ref) structure and optionally `w` is an `AbstractVector` of weights.

TODO mention that it can be constructed by vcatting operators as well.
mention flattening

# Examples
```jldoctest
julia> a = AggregationStack(SegmentedMean(2), SegmentedMax(2))
AggregationStack:
 SegmentedMean(ψ = Float32[0.0, 0.0])
 SegmentedMax(ψ = Float32[0.0, 0.0])

julia> a(Float32[0 1 2; 3 4 5], bags([1:1, 2:3]))
4×2 Matrix{Float32}:
 0.0  1.5
 3.0  4.5
 0.0  2.0
 3.0  5.0

julia> a = AggregationStack(SegmentedMean(2), AggregationStack(SegmentedMax(2)))
AggregationStack:
 SegmentedMean(ψ = Float32[0.0, 0.0])
 SegmentedMax(ψ = Float32[0.0, 0.0])

julia> a = SegmentedMeanMax(Float32, 2)
AggregationStack:
 SegmentedMean(ψ = Float32[0.0, 0.0])
 SegmentedMax(ψ = Float32[0.0, 0.0])

julia> vcat(SegmentedMean(2), SegmentedMax(2))
AggregationStack:
 SegmentedMean(ψ = Float32[0.0, 0.0])
 SegmentedMax(ψ = Float32[0.0, 0.0])

```

See also: [`AbstractAggregation`](@ref), [`SegmentedSum`](@ref), [`SegmentedMax`](@ref),
    [`SegmentedMean`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
struct AggregationStack{T <: Tuple{Vararg{AbstractAggregation}}} <: AbstractAggregation
    fs::T
    function AggregationStack(fs::Tuple{Vararg{AbstractAggregation}})
        ffs = _flatten_agg(fs)
        new{typeof(ffs)}(ffs)
    end
end

_flatten_agg(t) = tuple(vcat(map(_flatten_agg, t)...)...)
_flatten_agg(a::AggregationStack) = vcat(map(_flatten_agg, a.fs)...)
_flatten_agg(a::AbstractAggregation) = [a]

AggregationStack(fs::AbstractAggregation...) = AggregationStack(fs)

Flux.@functor AggregationStack

function (a::AggregationStack)(x::Union{AbstractArray, Missing}, bags::AbstractBags, args...)
    reduce(vcat, (f(x, bags, args...) for f in a.fs))
end

Flux.@forward AggregationStack.fs Base.getindex, Base.firstindex, Base.lastindex, Base.first,
    Base.last, Base.iterate, Base.eltype

Base.length(a::AggregationStack) = sum(length.(a.fs))
Base.size(a::AggregationStack) = tuple(sum(only, size.(a.fs)))

function Base.show(io::IO, ::MIME"text/plain", a::AggregationStack)
    print(io, "AggregationStack:\n")
    Base.print_array(io, a.fs |> collect)
end

function Base.show(io::IO, a::AggregationStack)
    if get(io, :compact, false)
        print(io, "AggregationStack")
    else
        s = join(repr.(a.fs; context=io), "; ")
        print(io, "[", s, "]")
    end
end


text:
function default_algorithm(prob::DiffEqBase.AbstractDDEProblem{uType,tType,lType,isinplace};kwargs...) where {uType,tType,lType,isinplace}
  o = Dict{Symbol,Any}(kwargs)
  extra_kwargs = Any[]; alg = MethodOfSteps(AutoTsit5(Rosenbrock23(autodiff=false))) # Standard default
  uEltype = eltype(prob.u0)

  alg_hints = get_alg_hints(o)
  tol_level = get_tolerance_level(o)

  if tol_level == :extreme_tol || tol_level == :low_tol
    stiff_alg = Rodas5(autodiff=false)
  else
    stiff_alg = Rosenbrock23(autodiff=false)
  end

  if :stiff ∈ alg_hints
    alg = MethodOfSteps(stiff_alg)
  elseif :nonstiff ∈ alg_hints || length(prob.u0) > 10000
    alg = MethodOfSteps(Tsit5())
  else # :auto
    if tol_level == :extreme_tol || tol_level == :low_tol
      alg = MethodOfSteps(AutoVern7(stiff_alg))
    else
      alg = MethodOfSteps(AutoTsit5(stiff_alg))
    end
  end

  # If adaptivity is not set and the tType is not a float, turn off adaptivity
  # Bad interaction with ForwardDiff
  #!(tType <: AbstractFloat) && (:adaptive ∉ keys(o)) && push!(extra_kwargs,:adaptive=>false)

  alg,extra_kwargs
end

max_stars_repo_name: ChrisRackauckas/DiffEq


text:
################################################################################
#                          TEST: Euler Angle and Axis
################################################################################

# Test constructions.
for k = 1:samples
    # Sample a random vector and angle.
    a = rand()
    v = rand(3)

    angleaxis = EulerAngleAxis(a,v)

    # Check if the representation was created correctly.
    @test angleaxis.a == a
    @test angleaxis.v == v

    # Test multiplication.
    a1  = -π + 2*π*rand()
    v1  = randn(3)
    v1 /= norm(v1)

    a2  = -π + 2*π*rand()
    v2  = randn(3)
    v2 /= norm(v2)

    ea1 = EulerAngleAxis(a1, v1)
    ea2 = EulerAngleAxis(a2, v2)
    ear = ea2*ea1

    q1  = angleaxis_to_quat(ea1)
    q2  = angleaxis_to_quat(ea2)
    qr  = q1*q2
    eaq = quat_to_angleaxis(qr)

    @test    ear.a ≈ eaq.a
    @test ear.v[1] ≈ eaq.v[1]
    @test ear.v[2] ≈ eaq.v[2]
    @test ear.v[3] ≈ eaq.v[3]

    # Test inversion.
    a   = -π + 2*π*rand()
    v   = randn(3)
    v  /= norm(v)
    ea  = EulerAngleAxis(a,v)
    ear = ea*inv_rotation(ea)

    @test norm(ear.v) ≈ 0 atol = 1e-7
    @test       ear.a ≈ 0 atol = 1e-7
end

# Test types.
@test typeof(EulerAngleAxis(1, [1;1;1])) <: EulerAngleAxis{Int}
@test typeof(EulerAngleAxis(1.f0, Float32[1;1;1])) <: EulerAngleAxis{Float32}
@test typeof(EulerAngleAxis(1.f0, [1;1;1])) <: EulerAngleAxis{Float32}

# Test multiplication that leads to unitary rotation.
let
    ea1 = EulerAngleAxis(pi, [1.0,0.0,0.0])
    ear = ea1*ea1

    @test norm(ear.v) ≈ 0
    @test       ear.a ≈ 0
end

# Test `Float16` and `Float32` multiplications.
let
    v   = SVector{3,Float32}(1,0,0)
    a   = Float32(1)
    ea1 = EulerAngleAxis(a,v)
    ear = ea1*ea1

    @test eltype(ear.a) == eltype(ear.v) == Float32

    v   = SVector{3,Float16}(1,0,0)
    a   = Float16(1)
    ea1 = EulerAngleAxis(a,v)
    ear = ea1*ea1

    @test eltype(ear.a) == eltype(ear.v) == Float16
end

# Test exceptions.
@test_throws Exception EulerAngleAxis(randn(), [1;1;1;1])
@test_throws Exception EulerAngleAxis(randn(), [1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{4,Int}[1;1;1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{4,Float32}[1;1;1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{4,Float64}[1;1;1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{2,Int}[1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{2,Float32}[1;1])
@test_throws Exception EulerAngleAxis(randn(), SVector{2,Float64}[1;1])


text:
#MCHAMMER Custom Distributions
#Distributions not currently implemented in Julia
#by Eric Torkia, August 2020

"""
    RiskEvent(Prob, Distribution, Trials; seed=0)

Risk Events are defined as conditional distributions that will inflate the 0.

The RiskEvent() allows you to conditionally sample any distribution for as many trials as defined.
*Prob* is the conditional probability of sampling the impact Distribution
*Distribution* is any univariate distribution from Distributions.jl
*Trials* is the number of total iterations.
*seed* allows you to set a seed for the RiskEvent. Left blank or set to 0, the seed is set to random.

"""
function RiskEvent(Prob, Distribution, Trials; seed=0)
    if seed != 0 Random.seed!(seed) end
    C_dist = rand(Bernoulli(Prob),Trials) .* rand(Distribution,Trials)
    return C_dist
end

## Metalogs
#Ported from https://github.com/colsmit/PyMetalog/blob/master/pymetalog/support.py Metalogs

function mlprobs(Array, Increment)
    ArrayLength = size(Array,1)
    SA = sort(Array,1,rev=false)
    ML_Probs = []
end


text:
function define_time(machine)
    define(machine, "now", op_now)
end

function op_now(machine, _...)
    push!(machine.stack, time() |> round |> Int64)
    return 1
end

op_time =
"""
: ago ;
: mins  60 * ;
: hours 3600 * ;
: days  86400 * ;
: weeks 604800 * ;
"""


text:
export get_variables, get_variabletype, has_variable, add_variable!, add_variables!
export get_objective, set_objective!
export has_constraint, get_constraints, get_constraint, add_constraint!, rm_constraint!
export get_slacks, get_relative_slacks, get_minslack

###############################################################################
#### Variables
###############################################################################
get_variables(pb::Problem) = pb.variables

function get_variabletype(pb::Problem, varName::String)
  if !haskey(pb.variables, varName)
    error(LOGGER, "get_variable(): Current Problem has no variable named ", varName)
  end
  pb.variables[varName]
end

has_variable(pb::Problem, var::Variable) = haskey(pb.variables, var.name) && pb.variables[var.name] == var.kind

function add_variable!(pb::Problem, x::Pair{String, T}) where T
  if haskey(pb.variables, x[1]) && !(x[2] <: pb.variables[x[1]])
    error(LOGGER, "add_variable!(): Attempting to add variable ", x.name, " (", x.kind, ") when ", x, " (", get_variabletype(pb, x.name), ") already exists.")
  end
  pb.variables[x[1]] = x[2]
end

"""
  add_variable!(pb, x::Variable)

  Add the `x` variable to `pb` set of variables.
"""
function add_variable!(pb::Problem, x::Variable)
  add_variable!(pb, x.name => x.kind)
end

"""
  add_variables!(pb, expo::Exponent)

  Add variables from `expo` exponent to `pb` set of variables.
"""
function add_variables!(pb::Problem, expo::Exponent)
  for var in keys(expo)
    add_variable!(pb, var)
  end
end

"""
  add_variables!(pb, p::Polynomial)

  Add variables from `p` polynomial to `pb` set of variables.
"""
function add_variables!(pb::Problem, p::Polynomial)
  vars = Set{Variable}()
  for expo in keys(p)
    for var in keys(expo)
      push!(vars, var)
    end
  end
  for var in vars
    add_variable!(pb, var)
  end
  return
end

"""
  add_variables!(pb, ctr::Constraint)

  Add variables from `ctr` constraint to `pb` set of variables.
"""
function add_variables!(pb::Problem, ctr::Constraint)
  add_variables!(pb, ctr.p)
  return
end

###############################################################################
#### Objective
###############################################################################
get_objective(pb::Problem) = pb.objective

"""
get_objective(pb::Problem, pt::Point)

Return the polynomial objective evaluated at `pt`.
"""
get_objective(pb::Problem, pt::Point) = evaluate(pb.objective, pt)

function set_objective!(pb::Problem, p::Polynomial)
  add_variables!(pb, p)
  pb.objective = p
  return
end

###############################################################################
#### Constraints
###############################################################################
has_constraint(pb::Problem, cstrName::String) = haskey(pb.constraints, cstrName)

get_constraints(pb::Problem) = pb.constraints

"""
  cstr = get_constraint(pb::Problem, cstrname::String)

Return the `cstrname` constraint from `pb`.
"""
get_constraint(pb::Problem, cstrname::String) = pb.constraints[cstrname]

"""
  cstr = get_constraint(pb::Problem, pt::Point)

Return the point (dict) of constraints names to their body's value when
evaluated at `pt`.
"""
function get_constraints(pb::Problem, pt::Point)
  img = Point()
  for (cstrName, cstr) in pb.constraints
    img[Variable(cstrName, Complex)] = evaluate(cstr.p, pt)
  end
  return img
end

"""
  add_constraint!(pb::Problem, cstrName::String, cstr::Constraint)

Add the constraint `cstr` under the `cstrname` name in the `pb` problem.
"""
function add_constraint!(pb::Problem, cstrName::String, cstr::Constraint)
  if haskey(pb.constraints, cstrName)
    warn(LOGGER, "add_constraint!(): A constraint with that name already exists ($cstrName)")
  end
  if real(cstr.ub) < real(cstr.lb) || imag(cstr.ub) < imag(cstr.lb)
    warn(LOGGER, "add_constraint!(): ", cstrName, " Lower bound is higher than upper bound ($(cstr.lb) - $(cstr.ub))")
  end
  add_variables!(pb, cstr)
  pb.constraints[cstrName] = cstr
end

rm_constraint!(pb::Problem, cstrName::String) = pop!(pb.constraints, cstrName)

###############################################################################
#### Solution problem intel evaluation
###############################################################################
"""
  pt = get_slacks(pb::Problem, pt::Point)

Return a point associating a constraint name to its slack, i.e. the minimum
algebraic distance between the real and imaginary parts of the body's value at
`pt` and its bounds.
"""
function get_slacks(pb::Problem, pt::Point)
  var_arr = Variable[]
  val_arr = Complex[]
  for (cstrName, cstr) in pb.constraints
    val = evaluate(cstr.p, pt)
    isa(val, Number) || error(LOGGER, "get_slacks(): constraint $cstrName not fully evaluated at provided point.\nEvaluated value is $val.")
    push!(var_arr, Variable(cstrName, Complex))
    push!(val_arr, min(real(val-cstr.lb), real(cstr.ub-val)) + min(imag(val-cstr.lb), imag(cstr.ub-val))*im)
  end
  return Point(var_arr, val_arr)
end

function get_relative_slacks(pb::Problem, pt::Point)
  var_arr = Variable[]
  val_arr = Complex[]
  for (cstrName, cstr) in pb.constraints
    val = evaluate(cstr.p, pt)
    lb = cstr.lb
    ub = cstr.ub
    isa(val, Number) || error(LOGGER, "get_slacks(): constraint $cstrName not fully evaluated at provided point.\nEvaluated value is $val.")
    push!(var_arr, Variable(cstrName, Complex))
    if lb==ub
      infeas = abs(real(val-ub))/max(abs(real(ub)),1) + im *abs(imag(val-ub))/max(abs(imag(ub)),1)
    else
      infeas_ub = max(0,real(val-ub))/max(abs(real(ub)),1) + im *max(0,imag(val-ub))/max(abs(imag(ub)),1)
      infeas_lb = max(0,real(lb-val))/max(abs(real(ub)),1) + im *max(0,imag(lb-val))/max(abs(imag(ub)),1)
      infeas = max(real(infeas_ub), real(infeas_lb)) + im * max(imag(infeas_ub), imag(infeas_lb))
    end
    push!(val_arr, infeas)
  end
  return Point(var_arr, val_arr)
end

function get_minslack(pb::Problem, pt::Point)
  minSlack = +Inf
  minCstrName = ""
  slacks = get_slacks(pb, pt)
  for (cstrName, slack) in slacks
    if minSlack > min(real(slack), imag(slack))
      minSlack = min(real(slack), imag(slack))
      minCstrName = cstrName
    end
  end
  return minSlack, minCstrName
end

function get_relativemaxslack(pb::Problem, pt::Point)
  maxRelSlack = - Inf
  maxCstrName = ""
  slacks = get_relative_slacks(pb, pt)
  for (cstrName, slack) in slacks
    max_slack = max(real(slack), imag(slack))
    if maxRelSlack < max_slack
      maxRelSlack = max_slack
      maxCstrName = cstrName
    end
  end
  return maxRelSlack, maxCstrName
end

function get_nb_infeasible_ctr_by_ctrtype(pb::Problem, pt::Point, epsilon::Float64)
  slacks = get_relative_slacks(pb, pt)
  nb_infeas_per_ctrtype = Dict{Tuple{String,String},Int64}()
  for (cstrName, slack) in slacks
    if real(slack) > epsilon || imag(slack) > epsilon
      ctr_elemid = split(cstrName.name, '_')
      scenario = ctr_elemid[1]
      bus_or_link = ctr_elemid[2]
      im_or_re = String(ctr_elemid[end])
      ctrtype = String(ctr_elemid[end-1])
      if ismatch(r"BALANCE", ctrtype)
        ctrtype = "BALANCE"
      end
      if !haskey(nb_infeas_per_ctrtype, (ctrtype,im_or_re))
        nb_infeas_per_ctrtype[(ctrtype,im_or_re)] = 0
      end
      nb_infeas_per_ctrtype[(ctrtype, im_or_re)] +=1
    end
  end
  return nb_infeas_per_ctrtype
end


text:
module HHPC

# Stdlib Requirements
using Logging
using Printf
using Random
using Statistics

# Non-stdlib essentials
using StaticArrays

# Open-loop requirements
using Graphs

# POMDP Stuff
using POMDPs
using POMDPModelTools
using POMDPModels
using POMDPPolicies
using POMDPSimulators
using LocalFunctionApproximation
using LocalApproximationValueIteration

# For saving and loading
using JLD2
using FileIO

# CMSSPs submodule
using CMSSPs.CMSSPModel

# Open-loop exports
export
    GraphTracker,
    metadatatype,
    bookkeepingtype,
    open_loop_plan!

# Closed-loop exports
export
    ModalHorizonPolicy,
    ModalFinInfHorPolicy,
    ModalStateAugmented,
    ModalAction,
    ModalFinHorMDP,
    ModalInfHorMDP,
    set_horizon_limit!,
    compute_terminalcost_localapprox!,
    finite_horizon_VI_localapprox!,
    infinite_horizon_VI_localapprox,
    compute_min_value_per_horizon_localapprox!,
    horizon_weighted_value,
    inf_hor_value,
    horizon_weighted_actionvalue,
    get_best_intramodal_action

# HHPC framework
# Requirements
# ensure that mode-switch MDP available?
export
    HHPCSolver,
    set_start_state!,
    set_open_loop_samples!

# Utils
export
    save_localapproxvi_policy_to_jld2,
    load_localapproxvi_policy_from_jld2,
    load_modal_horizon_policy_localapprox,
    save_modal_horizon_policy_localapprox,
    log2space_symmetric,
    polyspace_symmetric

# General functions that need to be implemented
export
    # Local
    expected_reward,
    get_relative_state,
    get_mdp,
    # Global
    update_vertices_with_context!,
    generate_goal_vertex_set!,
    generate_next_valid_modes,
    generate_bridge_vertex_set!,
    update_next_target!,
    # HHPC
    simulate_cmssp!,
    update_context_set!,
    display_context_future

# Local reqs
function get_relative_state end
function expected_reward end
function get_mdp end
# Global reqs
function update_vertices_with_context! end
function generate_goal_vertex_set! end
function generate_next_valid_modes end
function generate_bridge_vertex_set! end
function update_next_target! end
# HHPC reqs
function simulate_cmssp! end
function update_context_set! end
function get_bridging_action end
function display_context_future end

include("global_layer.jl")
include("local_layer.jl")
include("hhpc_framework.jl")
include("utils.jl")

end


text:
__precompile__()

module Equilibrium

using EFIT
using LinearAlgebra
using Interpolations
using ForwardDiff
using StaticArrays
using Trapz

import PolygonOps
import Contour
import Optim

const mu0 = 4*pi*1e-7

abstract type AbstractEquilibrium end
export AbstractEquilibrium

# Fallbacks
_not_implemented(M) = error("$(typeof(M)) has not implemented this functionality")
(M::AbstractEquilibrium)(x,y) = _not_implemented(M)
magnetic_axis(M::AbstractEquilibrium,r,z) = _not_implemented(M)
limits(M::AbstractEquilibrium) = _not_implemented(M)
psi_limits(M::AbstractEquilibrium) = _not_implemented(M)
psi_gradient(M::AbstractEquilibrium) = _not_implemented(M)
electric_potential(M::AbstractEquilibrium, psi) = _not_implemented(M)
electric_potential_gradient(M::AbstractEquilibrium) = _not_implemented(M)
pressure(M::AbstractEquilibrium, psi) = _not_implemented(M)
poloidal_current(M::AbstractEquilibrium, psi) = _not_implemented(M)
pressure_gradient(M::AbstractEquilibrium, psi) = _not_implemented(M)
poloidal_current_gradient(M::AbstractEquilibrium, psi) = _not_implemented(M)
cocos(M::AbstractEquilibrium) = _not_implemented(M)
B0Ip_sign(M::AbstractEquilibrium) = _not_implemented(M)

# Equilibrium API
export magnetic_axis, limits, psi_limits, psi_gradient, electric_potential, electric_potential_gradient
export pressure, poloidal_current, pressure_gradient, poloidal_current_gradient, safety_factor, cocos, B0Ip_sign
export plasma_current, beta_n

include("cocos.jl")
export COCOS, cocos, check_cocos, identify_cocos, transform_cocos
export cylindrical_cocos, poloidal_cocos, cylindrical_cocos_indices, poloidal_cocos_indices

include("boundary.jl")
export Boundary, PlasmaBoundary, FluxSurface, Wall, in_boundary, in_plasma, in_vessel
export boundary, flux_surface, circumference, average, area, area_average, volume, volume_average

include("solovev.jl")
export SolovevEquilibrium, solovev, clear_cache

include("fields.jl")
export Bfield, Efield, Jfield, EMFields, fields, gradB, curlB, poloidal_Bfield, poloidal_Jfield

include("efit.jl")
export EFITEquilibrium, efit

include("efit_io.jl")
export read_geqdsk

include("transp_io.jl")
export transp_potential!

end # module


text:
module NetworkEpidemics

using ProgressMeter
using LinearAlgebra
using OrdinaryDiffEq
using Graphs
using Graphs.SimpleGraphs
using StatsBase
using DataStructures
using FillArrays

include("utils/categorical_tree.jl")
include("utils/rand.jl")
include("utils/unzip.jl")
include("utils/laplacian.jl")
include("utils/arrays.jl")
include("utils/NotImplementedError.jl")

export normalized_laplacian

include("models/AbstractEpidemicModel.jl")

export AbstractEpidemicModel

include("models/models/metapopulation.jl")
include("models/models/contact_process.jl")
include("models/models/metaplex.jl")
include("models/models/heterogeneous_metaplex.jl")

include("stochastic/gillespie.jl")
include("stochastic/average.jl")
include("deterministic/meanfield.jl")

export gillespie, average, meanfield

include("models/dynamics/AbstractCompartimentalModel.jl")
include("models/dynamics/SI.jl")
include("models/dynamics/SIS.jl")
include("models/dynamics/SIR.jl")

export Metapopulation, ContactProcess, Metaplex
export HeterogeneousMetaplex, HeterogeneousMetapopulation
export AbstractCompartimentalModel
export SI, SIS, SIR

end # module


text:
# Standard isentropic vortex test case.  For a more complete description of
# the setup see for Example 3 of:
#
# @article{ZHOU2003159,
#   author = {Y.C. Zhou and G.ρw. Wei},
#   title = {High resolution conjugate filters for the simulation of flows},
#   journal = {Journal of Computational Physics},
#   volume = {189},
#   number = {1},
#   pages = {159--179},
#   year = {2003},
#   doi = {10.1016/S0021-9991(03)00206-7},
#   url = {https://doi.org/10.1016/S0021-9991(03)00206-7},
# }
#
# This version runs the isentropic vortex as a stand alone test (no dependence
# on CLIMA moist thermodynamics)

using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.SpaceMethods
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
using CLIMA.LinearSolvers
using CLIMA.GeneralizedConjugateResidualSolver

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const _nstate = 5
const _δρ, _δρu, _δρv, _δρw, _δρe = 1:_nstate
const _δρu⃗ = SVector(_δρu, _δρv, _δρw)
const stateid = (ρid = _δρ, ρuid = _δρu, ρvid = _δρv, ρwid = _δρw, ρeid = _δρe)
const statenames = ("δρ", "δρu", "δρv", "δρw", "δρe")
const _nauxstate = 5
const _a_ρ_ref, _a_ρu_ref, _a_ρv_ref, _a_ρw_ref, _a_ρe_ref = 1:_nauxstate
const γ_exact = 7 // 5

const _Tinf = 1
const _uinf = 2
const _vinf = 1

@inline function velocity_pressure(Q, aux)
    γ::eltype(Q) = γ_exact
    @inbounds begin
        δρ, δρu, δρv, δρw, δρe = Q[_δρ], Q[_δρu], Q[_δρv], Q[_δρw], Q[_δρe]

        ρ_ref, ρe_ref = aux[_a_ρ_ref], aux[_a_ρe_ref]
        ρu_ref, ρv_ref, ρw_ref = aux[_a_ρu_ref], aux[_a_ρv_ref], aux[_a_ρw_ref]

        ρ = ρ_ref + δρ
        ρu = ρu_ref + δρu
        ρv = ρv_ref + δρv
        ρw = ρw_ref + δρw
        ρe = ρe_ref + δρe

        ρinv = 1 / ρ
        u⃗ = SVector(ρinv * ρu, ρinv * ρv, ρinv * ρw)
        ((γ - 1) * (ρe - ρinv * (ρu^2 + ρv^2 + ρw^2) / 2), u⃗, ρinv)
    end
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t)
    P, u⃗, ρinv = velocity_pressure(Q, aux)
    n⃗ = SVector(n)
    γ::eltype(Q) = γ_exact
    abs(n⃗' * u⃗) + sqrt(ρinv * γ * P)
end
@inline function wavespeed_linear(n, Q, aux, t)
    ρ_ref, ρe_ref = aux[_a_ρ_ref], aux[_a_ρe_ref]
    ρu⃗_ref = SVector(aux[_a_ρu_ref], aux[_a_ρv_ref], aux[_a_ρw_ref])
    ρinv_ref = 1 / ρ_ref
    n⃗ = SVector(n)
    γ::eltype(Q) = γ_exact
    u⃗_ref = ρinv_ref * ρu⃗_ref
    P_ref = (γ - 1) * (ρe_ref - u⃗_ref' * ρu⃗_ref / 2)
    abs(n⃗' * u⃗_ref) + sqrt(ρinv_ref * γ * P_ref)
end

# physical flux function
@inline function eulerflux!(F, Q, QV, aux, t)
    P, u⃗, ρinv = velocity_pressure(Q, aux)
    @inbounds begin
        δρ, δρe = Q[_δρ], Q[_δρe]
        δρu⃗ = SVector(Q[_δρu], Q[_δρv], Q[_δρw])

        ρ_ref, ρe_ref = aux[_a_ρ_ref], aux[_a_ρe_ref]
        ρu⃗_ref = SVector(aux[_a_ρu_ref], aux[_a_ρv_ref], aux[_a_ρw_ref])

        ρ = ρ_ref + δρ
        ρe = ρe_ref + δρe

        ρu⃗ = ρu⃗_ref + δρu⃗

        F[:, _δρ] = ρu⃗
        F[:, _δρu⃗] = u⃗ * ρu⃗' + P * I
        F[:, _δρe] = u⃗ * (ρe + P)
    end
end

@inline function linearized_eulerflux!(F, Q, QV, aux, t)
    @inbounds begin
        FT = eltype(Q)
        γ::FT = γ_exact

        δρ, δρe = Q[_δρ], Q[_δρe]
        δρu⃗ = SVector(Q[_δρu], Q[_δρv], Q[_δρw])

        ρ_ref, ρe_ref = aux[_a_ρ_ref], aux[_a_ρe_ref]
        ρu⃗_ref = SVector(aux[_a_ρu_ref], aux[_a_ρv_ref], aux[_a_ρw_ref])

        ρinv_ref = 1 / ρ_ref
        u⃗_ref = ρinv_ref * ρu⃗_ref
        e_ref = ρinv_ref * ρe_ref

        P_ref = (γ - 1) * (ρe_ref - u⃗_ref' * ρu⃗_ref / 2)
        δP = (γ - 1) * (δρe - u⃗_ref' * δρu⃗ + u⃗_ref' * u⃗_ref * δρ)

        p_ref = ρinv_ref * P_ref

        # FIXME: Not sure these pure reference terms are right...
        #=
        F[:, _δρ ] = (ρu⃗_ref + δρu⃗)
        F[:, _δρu⃗] = (ρu⃗_ref * u⃗_ref' - (u⃗_ref * u⃗_ref') * δρ
                      + δρu⃗ * u⃗_ref' + u⃗_ref * δρu⃗'
                      + (P_ref + δP) * I)
        F[:, _δρe] = ((e_ref + p_ref) * ρu⃗_ref + (e_ref + p_ref) * δρu⃗ +
                      (δρe + δP) * u⃗_ref - (e_ref + p_ref) * u⃗_ref * δρ)
        =#
        F[:, _δρ] = δρu⃗
        F[:, _δρu⃗] = (
            +δρu⃗ * u⃗_ref' + u⃗_ref * δρu⃗' - (u⃗_ref * u⃗_ref') * δρ + δP * I
        )
        F[:, _δρe] = (
            (e_ref + p_ref) * δρu⃗ + u⃗_ref * (δρe + δP) -
            (e_ref + p_ref) * u⃗_ref * δρ
        )
    end
end

# initial condition
@inline function auxiliary_state_initialization!(aux, x, y, z)
    @inbounds begin
        FT = eltype(aux)
        γ::FT = γ_exact

        uinf::FT = _uinf
        vinf::FT = _vinf
        Tinf::FT = _Tinf

        ρ_ref = (Tinf)^(1 / (γ - 1))
        P_ref = ρ_ref^γ

        u_ref = uinf
        v_ref = vinf
        w_ref = zero(FT)

        ρu_ref = ρ_ref * u_ref
        ρv_ref = ρ_ref * v_ref
        ρw_ref = ρ_ref * w_ref

        ρe_ref =
            P_ref / (γ - 1) + (1 // 2) * ρ_ref * (u_ref^2 + v_ref^2 + w_ref^2)

        aux[_a_ρ_ref] = ρ_ref
        aux[_a_ρu_ref] = ρu_ref
        aux[_a_ρv_ref] = ρv_ref
        aux[_a_ρw_ref] = ρw_ref
        aux[_a_ρe_ref] = ρe_ref
    end
end

const halfperiod = 5
function isentropicvortex!(Q, t, x, y, z, aux)
    @inbounds begin
        ρ_ref, ρe_ref = aux[_a_ρ_ref], aux[_a_ρe_ref]
        ρu_ref, ρv_ref, ρw_ref = aux[_a_ρu_ref], aux[_a_ρv_ref], aux[_a_ρw_ref]

        FT = eltype(Q)

        γ::FT = γ_exact
        uinf::FT = _uinf
        vinf::FT = _vinf
        Tinf::FT = _Tinf
        λ::FT = 5

        xs = x - uinf * t
        ys = y - vinf * t

        # make the function periodic
        xtn = floor((xs + halfperiod) / (2halfperiod))
        ytn = floor((ys + halfperiod) / (2halfperiod))
        xp = xs - xtn * 2 * halfperiod
        yp = ys - ytn * 2 * halfperiod

        rsq = xp^2 + yp^2

        u = uinf - λ * (1 // 2) * exp(1 - rsq) * yp / π
        v = vinf + λ * (1 // 2) * exp(1 - rsq) * xp / π
        w = zero(FT)

        ρ =
            (
                Tinf - ((γ - 1) * λ^2 * exp(2 * (1 - rsq)) / (γ * 16 * π * π))
            )^(1 / (γ - 1))
        P = ρ^γ
        ρu = ρ * u
        ρv = ρ * v
        ρw = ρ * w
        ρe = P / (γ - 1) + (1 // 2) * ρ * (u^2 + v^2 + w^2)

        δρ = ρ - ρ_ref
        δρu = ρu - ρu_ref
        δρv = ρv - ρv_ref
        δρw = ρw - ρw_ref
        δρe = ρe - ρe_ref

        Q[_δρ], Q[_δρu], Q[_δρv], Q[_δρw], Q[_δρe] = δρ, δρu, δρv, δρw, δρe
    end
end

function main(
    mpicomm,
    FT,
    topl::AbstractTopology{dim},
    N,
    timeend,
    ArrayType,
    dt,
) where {dim}

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    # spacedisc = data needed for evaluating the right-hand side function
    nonlin_spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = _nstate,
        flux! = eulerflux!,
        numerical_flux! = (x...) ->
            NumericalFluxes.rusanov!(x..., eulerflux!, wavespeed),
        auxiliary_state_length = _nauxstate,
        auxiliary_state_initialization! = auxiliary_state_initialization!,
    )

    lin_spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = _nstate,
        flux! = linearized_eulerflux!,
        numerical_flux! = (x...) -> NumericalFluxes.rusanov!(
            x...,
            linearized_eulerflux!,
            (_...) -> 0,
        ), # central
        # wavespeed_linear), # rusanov!
        auxiliary_state_length = _nauxstate,
        auxiliary_state_initialization! = auxiliary_state_initialization!,
    )

    # This is a actual state/function that lives on the grid
    initialcondition(Q, x...) = isentropicvortex!(Q, FT(0), x...)
    Q = MPIStateArray(nonlin_spacedisc, initialcondition)
    dQ = similar(Q)

    linearsolver = GeneralizedConjugateResidual(3, Q, rtol = 1e-10)
    ark = ARK548L2SA2KennedyCarpenter(
        nonlin_spacedisc,
        lin_spacedisc,
        linearsolver,
        Q;
        dt = dt,
        t0 = 0,
    )

    eng0 = norm(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf """Update
        simtime = %.16e
        runtime = %s
        norm(Q) = %.16e""" ODESolvers.gettime(ark) Dates.format(
                convert(Dates.DateTime, Dates.now() - starttime[]),
                Dates.dateformat"HH:MM:SS",
            ) energy
        end
    end

    #= Paraview calculators:
    P = (0.4) * (ρe  - (ρu^2 + ρv^2 + ρw^2) / (2*ρ) - 9.81 * ρ * coordsZ)
    theta = (100000/287.0024093890231) * (P / 100000)^(1/1.4) / ρ
    =#
    step = [0]
    mkpath("vtk")
    cbvtk = GenericCallbacks.EveryXSimulationSteps(100) do (init = false)
        outprefix = @sprintf(
            "vtk/isentropicvortex_%dD_mpirank%04d_step%04d",
            dim,
            MPI.Comm_rank(mpicomm),
            step[1]
        )
        @debug "doing VTK output" outprefix
        writevtk(outprefix, Q, nonlin_spacedisc, statenames)
        pvtuprefix = @sprintf("isentropicvortex_%dD_step%04d", dim, step[1])
        prefixes = ntuple(
            i -> @sprintf(
                "vtk/isentropicvortex_%dD_mpirank%04d_step%04d",
                dim,
                i - 1,
                step[1]
            ),
            MPI.Comm_size(mpicomm),
        )
        writepvtu(pvtuprefix, prefixes, statenames)
        step[1] += 1
        nothing
    end

    solve!(Q, ark; timeend = timeend, callbacks = (cbinfo, cbvtk))

    # Print some end of the simulation information
    engf = norm(Q)
    Qe = MPIStateArray(
        nonlin_spacedisc,
        (Q, x...) -> isentropicvortex!(Q, FT(timeend), x...),
    )
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

function run(mpicomm, ArrayType, dim, Ne, N, timeend, FT, dt)
    brickrange = ntuple(
        j -> range(FT(-halfperiod); length = Ne[j] + 1, stop = halfperiod),
        dim,
    )
    topl =
        BrickTopology(mpicomm, brickrange, periodicity = ntuple(j -> true, dim))
    main(mpicomm, FT, topl, N, timeend, ArrayType, dt)
end

using Test
let
    CLIMA.init()
    ArrayTypes = (CLIMA.array_type(),)

    mpicomm = MPI.COMM_WORLD

    timeend = 1
    numelem = (5, 5, 1)

    polynomialorder = 4

    expected_error = Array{Float64}(undef, 2, 3) # dim-1, lvl
    expected_error[1, 1] = 5.7107915113422314e-01
    expected_error[1, 2] = 6.9428317829549141e-02
    expected_error[1, 3] = 3.2924199180417437e-03
    expected_error[2, 1] = 1.8059108422662815e+00
    expected_error[2, 2] = 2.1955161886731539e-01
    expected_error[2, 3] = 1.0411547208036219e-02
    lvls = integration_testing ? size(expected_error, 2) : 1

    @testset "$(@__FILE__)" for ArrayType in ArrayTypes
        for FT in (Float64,) #Float32)
            for dim in 2:3
                err = zeros(FT, lvls)
                for l in 1:lvls
                    Ne = ntuple(j -> 2^(l - 1) * numelem[j], dim)
                    dt = 3e-1 / Ne[1]
                    nsteps = ceil(Int64, timeend / dt)
                    dt = timeend / nsteps
                    @info (ArrayType, FT, dim)
                    err[l] = run(
                        mpicomm,
                        ArrayType,
                        dim,
                        Ne,
                        polynomialorder,
                        timeend,
                        FT,
                        dt,
                    )
                    @test err[l] ≈ FT(expected_error[dim - 1, l])
                end
                @info begin
                    msg = ""
                    for l in 1:(lvls - 1)
                        rate = log2(err[l]) - log2(err[l + 1])
                        msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                    end
                    msg
                end
            end
        end
    end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using CUDA.CUDNN:
    cudnnOpTensor,
    cudnnOpTensor!,
    cudnnOpTensorDescriptor,
        cudnnOpTensorDescriptor_t,
        cudnnCreateOpTensorDescriptor,
        cudnnSetOpTensorDescriptor,
        cudnnGetOpTensorDescriptor,
        cudnnDestroyOpTensorDescriptor,
    cudnnOpTensorOp_t,
        CUDNN_OP_TENSOR_ADD,  # 0,
        CUDNN_OP_TENSOR_MUL,  # 1,
        CUDNN_OP_TENSOR_MIN,  # 2,
        CUDNN_OP_TENSOR_MAX,  # 3,
        CUDNN_OP_TENSOR_SQRT, # 4, performed only on first arg
        CUDNN_OP_TENSOR_NOT,  # 5, performed only on first arg
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnDataType

@testset "cudnn/optensor" begin
    @test cudnnOpTensorDescriptor(C_NULL) isa cudnnOpTensorDescriptor
    @test Base.unsafe_convert(Ptr, cudnnOpTensorDescriptor(C_NULL)) isa Ptr
    @test cudnnOpTensorDescriptor(CUDNN_OP_TENSOR_ADD,cudnnDataType(Float32),CUDNN_NOT_PROPAGATE_NAN) isa cudnnOpTensorDescriptor

    (ax1,ax2,ay) = rand.((10,10,10))
    (cx1,cx2,cy) = CuArray.((ax1,ax2,ay))

    function optensortest(;
        op=CUDNN_OP_TENSOR_ADD,
        nanOpt=CUDNN_NOT_PROPAGATE_NAN,
        compType=(eltype(ax1) <: Float64 ? Float64 : Float32),
        alpha1=1,
        alpha2=1,
        beta=0,
    )
        f1 = (op === CUDNN_OP_TENSOR_ADD ? alpha1*ax1 .+ alpha2*ax2 :
              op === CUDNN_OP_TENSOR_MUL ? (alpha1*ax1) .* (alpha2*ax2) :
              op === CUDNN_OP_TENSOR_MIN ? min.(alpha1*ax1, alpha2*ax2) :
              op === CUDNN_OP_TENSOR_MAX ? max.(alpha1*ax1, alpha2*ax2) :
              op === CUDNN_OP_TENSOR_SQRT ? sqrt.(alpha1*ax1) :
              op === CUDNN_OP_TENSOR_NOT ? 1 .- ax1 :
              error("Unknown optensor"))
        f2 = f1 .+ beta * ay
        d = cudnnOpTensorDescriptor(op,cudnnDataType(compType),nanOpt)
        @test f1 ≈ cudnnOpTensor(cx1, cx2; op, compType, nanOpt, alpha1, alpha2) |> Array
        @test f1 ≈ cudnnOpTensor(cx1, cx2, d; alpha1, alpha2) |> Array
        @test f2 ≈ cudnnOpTensor!(copy(cy), cx1, cx2; op, compType, nanOpt, alpha1, alpha2, beta) |> Array
        @test f2 ≈ cudnnOpTensor!(copy(cy), cx1, cx2, d; alpha1, alpha2, beta) |> Array
    end

    optensortest(op = CUDNN_OP_TENSOR_ADD)
    optensortest(op = CUDNN_OP_TENSOR_MUL)
    optensortest(op = CUDNN_OP_TENSOR_MIN)
    optensortest(op = CUDNN_OP_TENSOR_MAX)
    optensortest(op = CUDNN_OP_TENSOR_SQRT)
    optensortest(op = CUDNN_OP_TENSOR_NOT)
    optensortest(nanOpt = CUDNN_PROPAGATE_NAN)
    optensortest(alpha1 = 2)
    optensortest(alpha2 = 2)
    optensortest(beta = 2)
end


text:
##### Beginning of file

Test.@test( Base.VERSION >= VersionNumber("1.0") )

Test.@test( RemoveLFS.version() > VersionNumber(0) )

Test.@test(
    RemoveLFS.version() ==
        RemoveLFS.version(RemoveLFS)
    )

Test.@test(
    RemoveLFS.version() ==
        RemoveLFS.version(first(methods(RemoveLFS.eval)))
    )

Test.@test(
    RemoveLFS.version() ==
        RemoveLFS.version(RemoveLFS.eval)
    )

Test.@test(
    RemoveLFS.version() ==
        RemoveLFS.version(RemoveLFS.eval, (Any,))
    )

Test.@test( RemoveLFS.version(TestModuleA) == VersionNumber("1.2.3") )

Test.@test( RemoveLFS.version(TestModuleB) == VersionNumber("4.5.6") )

Test.@test_throws(
    ErrorException,
    RemoveLFS._TomlFile(joinpath(mktempdir(),"1","2","3","4")),
    )

##### End of file


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
"""
This function check the presence of columns with an odd number of elements equal to 1 on the coboundary 1 matrix
"""
function check_odd_elements_columns_CSC(A::Lar.SparseMatrixCSC)#test HOMEWORK2
	sleep_time = 5 #seconds
	temp = 0
	result = []
	for i in range(1,length=size(A,2))
		push!(result, temp)
		temp = 0
		for j in range(1, length=size(A,1))
	        temp+=abs(A[j,i])
		end
	#println(temp)#show column sums
	end
	count = 0
	for k in result
		if isodd(k)
			count+=1
		end
	end
	if count!=0
		println("//////////////////////////////////////////////////////////////////////////////////////////////////////////")
		println("CHECKING COLUMNS WITH AN ODD NUMBER OF ELEMENTS")
		println("WARNING, THE MATRIX CONTAINS $count COLUMNS WITH AN ODD NUMBER OF ELEMENTS")
		println("the execution will sleep for $sleep_time seconds")
		println("//////////////////////////////////////////////////////////////////////////////////////////////////////////")
		sleep(sleep_time)
	end
end


text:
struct ParseFactorCovariate{S<:Unsigned, T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    basefactor::AbstractFactor{S}
    transform::Function
end

Base.length(var::ParseFactorCovariate) = length(var.basefactor)

function ParseFactorCovariate(name::String, basefactor::AbstractFactor{S}, transform::Function) where {S<:Unsigned}
    ParseFactorCovariate{S, Float32}(name, basefactor, transform)
end

function ParseFactorCovariate(name::String, basefactor::AbstractFactor{S}, transform::Function, ::Type{T}) where {S<:Unsigned} where {T<:AbstractFloat}
    ParseFactorCovariate{S, T}(name, basefactor, transform)
end

function covariate(basefactor::AbstractFactor{S}, transform::Function) where {S<:Unsigned}
    ParseFactorCovariate(getname(basefactor), basefactor, transform)
end

function slice(cov::ParseFactorCovariate{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:Unsigned} where {T<:AbstractFloat}
    levels = getlevels(cov.basefactor)
    parsed = Vector{T}(undef, length(levels) + 1)
    for (index, level) in enumerate(levels)
        parsed[index + 1] = convert(T, cov.transform(level))
    end
    parsed[1] = convert(T, cov.transform(MISSINGLEVEL))
    f = (i -> parsed[i + 1])
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(cov.basefactor, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, T)
end

function Base.map(covariate::ParseFactorCovariate, dataframe::AbstractDataFrame)
    ParseFactorCovariate(covariate.name, map(covariate.basefactor, dataframe), covariate.transform)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
"""
    fritzhugh_nagumo(parameters, UA)

Return the FritzHugh-Nagumo[1,2] system of ordinary differential equations as an in-place
function real_de!(du, u, p, t) to update the gradient vector du,
and a partial FritzHugh-Nagumo system as a universal differential equations[3], where
the linear term involved in a saddle-node bifurcation in the second equation must be
learned from data (see below).

Parameters: Real[] with parameters [b0, b1, I, E] of the FritzHugh-Nagumo system
UA: Differentiable universal function approximator in the UDE sense: UA(u, p)

The FritzHugh-Nagumo system:
dx/dt = x - 1/3 x^3 - y + I
dy/dt = Eb₀ + Eb₁x - Ey

The partial FritzHugh-Nagumo system as a UDE:
dx/dt = x - 1/3 x^3 - y + I
dy/dt = Eb₀ + Uₚ([x, y]) - Ey

where UAₚ is the universal function approximator with weights p.

We are interested in the system near a saddle-node bifurcation with
b0, b1, I, E = parameters = [0.9, 0.5, 1.2, 1.25] -> one stable state
b0, b1, I, E = parameters = [0.9, 0.5, 1.0, 1.25] -> two stable states, one unstable state

[1] FitzHugh, Richard. "Impulses and physiological states in theoretical models of nerve membrane."
Biophysical journal 1.6 (1961): 445-466.
[2] Nagumo, Jinichi, Suguru Arimoto, and Shuji Yoshizawa. "An active pulse transmission line
simulating nerve axon." Proceedings of the IRE 50.10 (1962): 2061-2070.
[3] Rackauckas, Christopher, et al. "Universal differential equations for scientific machine
learning." arXiv preprint arXiv:2001.04385 (2020).
"""
function fritzhugh_nagumo(parameters, UA)
    b0, b1, I, E = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] =  u[1] - 1/3 * u[1]^3 - u[2] + I
        du[2] = E*b0 + E*b1*u[1] - E *u[2]
    end

    univ_de = (u, p, t) -> begin
        x,y = u
        z = UA(u,p)
        [x - 1/3 * x.^3 - y + I,
        E*b0 + z[1] - E * y]
    end

    return (real_de!, univ_de)
end

"""
    genetic_toggle_switch(parameters, UA)

Return a genetic toggle switch (GTS) system of ordinary differential equations due to Gardner et al [1]
as an in-place function real_de!(du, u, p, t) to update the gradient vector du,
and a partial genetic toggle switch system as a universal differential equations[3], where
the non-linear term in the first equation involved in a pitchfork bifurcation must be
learned from data (see below).

Parameters: Real[] with parameters [a1, a2, k, s] of the GTS system
UA: Differentiable universal function approximator in the UDE sense: UA(u, p)

The GTS system:
dx/dt = a₁/(1+yᵏ) - x
dy/dt = a₂/(1+xˢ) - y

The partial GTS system as a UDE:
dx/dt = UAₚ([x,y]) - x
dy/dt = a₂/(1+xˢ) - y

where UAₚ is the universal function approximator with weights p.

We are interested in the system near a pitchfork bifurcation with, parameterization due to
Boshe and Ghosh[3]
a1, a2, k, s = parameters = [1.5, 1.5, 2.0, 2.0] -> one stable state
a1, a2, k, s = parameters = [3.5, 3.5, 2.0, 2.0] -> two stable states, one unstable state

[1] Gardner, Timothy S., Charles R. Cantor, and James J. Collins. "Construction of a genetic toggle
switch in Escherichia coli." Nature 403.6767 (2000): 339-342.
[2] Rackauckas, Christopher, et al. "Universal differential equations for scientific machine
learning." arXiv preprint arXiv:2001.04385 (2020).
[3] Bose, Indrani, and Sayantari Ghosh. "Bifurcation and criticality." Journal of Statistical
Mechanics: Theory and Experiment 2019.4 (2019): 043403.
"""
function genetic_toggle_switch(parameters, UA)
    a1, a2, k, s = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] = a1 /(1 + u[2]^k) - u[1]
        du[2] = a2 /(1 + u[1]^s) - u[2]
    end

    univ_de = (u, p, t) -> begin
        x,y = u
        z = UA(u,p)
        [z[1] - x,
        a2 /(1 + (x^s)) - y]
    end

    return (real_de!, univ_de)
end

"""
    truscott_brindley(parameters, UA)

Return the Truscott-Brindley[1] system of ordinary differential equations as an in-place
function real_de!(du, u, p, t) to update the gradient vector du,
and a partial Truscott-Brindley system as a universal differential equations[2], where
a non-linear term in the first equation and a linear term in the second equation involved
in a Hopf bifurcation must be learned from data (see below).

!! The Truscott Brindley system is not well behaved when solved numerically!

Parameters: Real[] with parameters [a, b, c, d] of the Truscott-Brindley system
UA: Differentiable universal function approximator in the UDE sense: UA(u, p)

The Truscott Brindley system:
dx/dt = bx(1-x) - y(x²/(a²+x²))
dy/dt = dy(x²/(a²+x²)) - cdy

The partial Truscott Brindley system as a UDE:
dx/dt = bx(1-x) + UAₚ([x, y])₁
dy/dt = dy(x²/(a²+x²)) + UAₚ([x,y])₂

where UAₚ is the universal function approximator with weights p, and UA₁ is the first element
of the vector valued output.

We are interested in the system near a Hopf bifurcation with
a, b, c, d = parameters = [0.053, 0.43, 0.024/(0.05*0.7), 0.05] -> oscillatory solution
a, b, c, d = parameters = [0.053, 0.43, 0.34, 0.05] -> steady state solution (Kirchfall)

[1] Truscott, J. E., and J. Brindley. "Ocean plankton populations as excitable media." Bulletin of
Mathematical Biology 56.5 (1994): 981-998.
[2] Rackauckas, Christopher, et al. "Universal differential equations for scientific machine
learning." arXiv preprint arXiv:2001.04385 (2020).
"""
function truscott_brindley(parameters, UA)
    a,b,c,d = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] = b*u[1]*(1-u[1]) - u[2]*((u[1]^2)/(a^2+u[1]^2))
        du[2] = d*u[2]*((u[1]^2)/(a^2+u[1]^2)) - c*d*u[2]
    end

    univ_de = (u, p, t) -> begin
        P,Z = u
        z = UA(u,p)
        [b*P*(1-P) + z[1],                 #z[1]=- Z*((P^2)/(p_[1]^2+P^2)),
         d*Z*((P^2)/(a^2+P^2)) + z[2]]  #z[2]=- p_[3]*p_[4]*Z]
    end

    return (real_de!, univ_de)
end

"""
    roessler(parameters, UA)

Return the Roessler[1] system of ordinary differential equations as an in-place
function real_de!(du, u, p, t) to update the gradient vector du,
and a partial Roessler system as a universal differential equations[2], where
a linear term in the second equation involved in a period doubling bifurcation must be
learned from data (see below).

Parameters: Real[] with parameters [a, b, c] of the Roessler system
UA: Differentiable universal function approximator in the UDE sense: UA(u, p)

The Roessler system:
dx/dt = - y - z
dy/dt = x + ay
dz/dt = b + xz - cz

The partial Roessler system as a UDE:
dx/dt = - y - z
dy/dt = x + UAₚ([x,y,z])
dz/dt = b + xz - cz

where UAₚ is the universal function approximator with weights p.

We are interested in the system near a period doubling bifurcation with
a, b, c = parameters = [0.1, 0.1, 4.0]
a, b, c = parameters = [0.1, 0.1, 6.0] -> period doubled

[1] Rössler, Otto E. "Chaotic behavior in simple reaction systems." Zeitschrift für Naturforschung
A 31.3-4 (1976): 259-264.
[2] Rackauckas, Christopher, et al. "Universal differential equations for scientific machine
learning." arXiv preprint arXiv:2001.04385 (2020).
"""
function roessler(parameters, UA)
    a,b,c = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] = - u[2] - u[3]
        du[2] = u[1] + a*u[2]
        du[3] = b + u[3]*u[1] - c*u[3]
    end

    univ_de = (u, p, t) -> begin
        x,y,z = u
        zf = UA(u,p)
        [-y - z,
        x + zf[1],
        b + z*x - c*z]
    end

    return (real_de!, univ_de)
end

""" Selkov Model
parameters_oscillation = Float32[0.1, 0.6]
parameters_steadystate = Float32[0.1, 0.15]

u₀ = Float32[1.0, 1.0]

tspan = (0.0f0, 15.0f)
s = 0.1
loss_weights = Float32[0.2, 0.8]

UA = regression_model(
    2, # input_d
    1, # output_d
    16, # hidden_d
    1, # hidden_layers
    tanh, # nonlinearity
    Flux.glorot_normal # initialization
)
"""
function selkov(parameters, UA)
    a, b = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] = -u[1] + a*u[2] + (u[1]^2)*u[2]
        du[2] = b - a*u[2] - (u[1]^2)*u[2]
    end

    univ_de = (u, p, t) -> begin
        x,y = u
        z = UA(u,p)
        [-x + z[1] + (x^2)*y,
         b - a*y - (x^2)*y]
    end

    return (real_de!, univ_de)
end

""" NPZ Model, the short documentation!
parameters_oscillation = Float32[0.2, 0.2, 0.4, 1.5, 0.03, 0.05, 0.15, 0.04, 0.6, 0.25, 0.33, 0.5, 0.6, 0.035]
parameters_steadystate = Float32[0.2, 0.2, 0.4, 1.0, 0.03, 0.05, 0.15, 0.04, 0.6, 0.25, 0.33, 0.5, 0.6, 0.035]

u₀ = Float32[0.4, 0.1, 0.05]

tspan = (0.0f0, 80.0)
s = 1.0
loss_weights = Float32[0.4, 0.6]
lr = 1e-3

UA = regression_model(
    3, # input_d
    1, # output_d
    32, # hidden_d
    1, # hidden_layers
    tanh, # nonlinearity
    Flux.glorot_normal # initialization
)

θ₀ = 1e-4 .* initial_params(UA)
"""
function NPZ(parameters, UA)
    a, b, c, d, e, k, r, s, N0, al, be, ga, la, nu = parameters

    real_de! = (du, u, p, t) -> begin
        du[1] = -1.0*(u[1]/(e + u[1]))*(a/(b + c*u[2]))*u[2] +
                r*u[2] +
                ((be*la*u[2]^2*u[3])/(nu^2+u[2]^2)) +
                ga*d*u[3]^2 +
                k*(N0-u[1])
        du[2] = (u[1]/(e + u[1]))*(a/(b + c*u[2]))*u[2] -
                r.*u[2] -
                ((la* u[2]^2*u[3])/(nu^2 + u[2]^2)) -
                (s + k)*u[2]
        du[3] = ((al*la*u[2]^2*u[3])/(nu^2+u[2]^2)) -
                d*u[3]^2
    end

    univ_de = (u, p, t) -> begin
        N,P,Z = u
        z = UA(u,p)
        [-1.0*(N/(e + N))*(a/(b + c*P))*P + r*P + ((be*la* P^2*Z)/(nu^2 + P^2)) + ga*d*Z^2 + k*(N0-N),
        (N/(e + N)) * (a/(b + c*P))*P - r*P - ((la* P^2 *Z)/(nu^2 + P^2)) - (s+k)*P,
        ((al*la*P^2*Z)/(nu^2 + P^2)) + z[1]]
    end

    return (real_de!, univ_de)
end

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add literature references to AdaptiveMHTuning docstring.

"""
    AdaptiveMHTuning(...) <: MHProposalDistTuning

Adaptive MCMC tuning strategy for Metropolis-Hastings samplers.

Adapts the proposal function based on the acceptance ratio and covariance
of the previous samples.

Fields:

* `λ`: Controls the weight given to new covariance information in adapting
  the proposal distribution. Defaults to `0.5`.

* `α`: Metropolis-Hastings acceptance ratio target, tuning will try to
  adapt the proposal distribution to bring the acceptance ratio inside this
  interval. Defaults to `IntervalSets.ClosedInterval(0.15, 0.35)`

* `β`: Controls how much the spread of the proposal distribution is
  widened/narrowed depending on the current MH acceptance ratio.

* `c`: Interval for allowed scale/spread of the proposal distribution.
  Defaults to `ClosedInterval(1e-4, 1e2)`.

* `r`: Reweighting factor. Take accumulated sample statistics of previous
  tuning cycles into account with a relative weight of `r`. Set to `0` to
  completely reset sample statistics between each tuning cycle.

Constructors:

```julia
AdaptiveMHTuning(
    λ::Real,
    α::IntervalSets.ClosedInterval{<:Real},
    β::Real,
    c::IntervalSets.ClosedInterval{<:Real},
    r::Real
)
```
"""
@with_kw struct AdaptiveMHTuning <: MHProposalDistTuning
    λ::Float64 = 0.5
    α::IntervalSets.ClosedInterval{Float64} = ClosedInterval(0.15, 0.35)
    β::Float64 = 1.5
    c::IntervalSets.ClosedInterval{Float64} = ClosedInterval(1e-4, 1e2)
    r::Real = 0.5
end

export AdaptiveMHTuning

mutable struct ProposalCovTuner{
    S<:MCMCBasicStats
} <: AbstractMCMCTunerInstance
    config::AdaptiveMHTuning
    stats::S
    iteration::Int
    scale::Float64
end

(tuning::AdaptiveMHTuning)(chain::MHIterator) = ProposalCovTuner(tuning, chain)

function ProposalCovTuner(tuning::AdaptiveMHTuning, chain::MHIterator)
    m = totalndof(getdensity(chain))
    scale = 2.38^2 / m
    ProposalCovTuner(tuning, MCMCBasicStats(chain), 1, scale)
end

_approx_cov(target::Distribution) = cov(unshaped(target))
_approx_cov(target::DistLikeDensity) = cov(target)
_approx_cov(target::AbstractPosteriorDensity) = cov(getprior(target))
_approx_cov(target::BAT.TransformedDensity{<:Any,<:BAT.DistributionTransform}) =
    BAT._approx_cov(target.trafo.target_dist)

function tuning_init!(tuner::ProposalCovTuner, chain::MHIterator)
    Σ_unscaled = _approx_cov(getdensity(chain))
    Σ = Σ_unscaled * tuner.scale

    next_cycle!(chain) # ToDo: This would be better placed in the burn-in algorithm
    chain.proposaldist = set_cov(chain.proposaldist, Σ)

    nothing
end

function tuning_update!(tuner::ProposalCovTuner, chain::MHIterator, samples::DensitySampleVector)
    stats = tuner.stats
    stats_reweight_factor = tuner.config.r
    reweight_relative!(stats, stats_reweight_factor)
    # empty!.(stats)
    append!(stats, samples)

    config = tuner.config

    α_min = minimum(config.α)
    α_max = maximum(config.α)

    c_min = minimum(config.c)
    c_max = maximum(config.c)

    β = config.β

    t = tuner.iteration
    λ = config.λ
    c = tuner.scale
    Σ_old = Matrix(get_cov(chain.proposaldist))

    S = convert(Array, stats.param_stats.cov)
    a_t = 1 / t^λ
    new_Σ_unscal = (1 - a_t) * (Σ_old/c) + a_t * S

    α = eff_acceptance_ratio(chain)

    max_log_posterior = stats.logtf_stats.maximum

    if α_min <= α <= α_max
        chain.info = MCMCIteratorInfo(chain.info, tuned = true)
        @debug "MCMC chain $(chain.info.id) tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"
    else
        chain.info = MCMCIteratorInfo(chain.info, tuned = false)
        @debug "MCMC chain $(chain.info.id) *not* tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"

        if α > α_max && c < c_max
            tuner.scale = c * β
        elseif α < α_min && c > c_min
            tuner.scale = c / β
        end
    end

    Σ_new = new_Σ_unscal * tuner.scale

    chain.proposaldist = set_cov(chain.proposaldist, Σ_new)
    tuner.iteration += 1

    nothing
end


text:
function run_intcode(program, x, y)::Int
    program[2] = x
    program[3] = y
    pointer = 1
    while pointer <= length(program)
        opcode = program[pointer]
        if opcode == 1
            pt1 = program[pointer+1]+1
            pt2 = program[pointer+2]+1
            pt3 = program[pointer+3]+1
            program[pt3] = program[pt1] + program[pt2]
        elseif opcode == 2
            pt1 = program[pointer+1]+1
            pt2 = program[pointer+2]+1
            pt3 = program[pointer+3]+1
            program[pt3] = program[pt1] * program[pt2]
        elseif opcode == 99
            break
        else
            throw(ErrorException("Invalid opcode"))
        end
        pointer = pointer+4
    end

    return program[1]
end

program = open("intcode.csv") do f
    parse.(Int, split(readlines(f)[1], ","))
end
println(run_intcode(deepcopy(program), 12, 2))
for i in 0:99
    for j in 0:99
        res = try
            run_intcode(deepcopy(program), i, j)
        catch e
            -1
        end
        if res == 19690720
            println(100*i+j)
            break
        end
    end
end


text:
"""
Overview: collections of routines and utilities to
perform matrix assembly
"""

"""
> assembly_forcing_vector!(Cp⁺,mesh,U,M,ϕ,κ_modes,κ_list,κ_phase)
It assemblies the non-autonomous forcing vector

\$ \\mathbf{F} = ∑_{i=1}^{N} κᵢMϕᵢ \$

- Cp⁺ : no-autonomous parametrisation data structure
- mesh : mesh data structure
- U : displacement field
- M : mass matrix
- ϕ : eigenmodes
- κ_modes : list of forced modes
- κ_list : forcing amplitude κᵢ
- κ_phase : phase of the applied forcing
"""
function assembly_forcing_vector!(Cp⁺::Parametrisation,mesh::Grid,U::Field,
                                  M,ϕ,κ_modes,κ_list,κ_phase)
  #
  for (i,κ) in enumerate(κ_modes)
    if (κ_phase[i]=='c')
      Cp⁺.nlr -= (1.0/2.0)*κ_list[i]*M*ϕ[:,κ]
    elseif (κ_phase[i]=='s')
      Cp⁺.nlr -= (1.0/2.0/im)*κ_list[i]*M*ϕ[:,κ]
    end
  end
  #
  return nothing
  #
end

"""
> assembly_H_nl!(Cp, entry,Ψ₁, Ψ₂, Ψ₃,mesh, U, mult = 1.0)
It assemblies cubic nonlinearities operator

\$ G(Ψ₁,Ψ₂,Ψ₃) = \\frac{1}{6} ∫_{Ω} γ(Ψ₁,Ψ₂):\\mathcal{A}:γ(Ψ₃,w) + γ(Ψ₁,Ψ₃):\\mathcal{A}:γ(Ψ₂,w) + γ(Ψ₃,Ψ₂):\\mathcal{A}:γ(Ψ₁,w) dΩ \$

- Cp : parametrisation data structure
- entry : entrance of the reference array
- Ψ₁ : mapping
- Ψ₂ : mapping
- Ψ₃ : mapping
- mesh : mesh data structure
- U : displacement field
- mult : integral multiplier
"""
function assembly_H_nl!(Cp::Parametrisation, entry::Int64,
                        Ψ₁, Ψ₂, Ψ₃,
                        mesh::Grid, U::Field, mult = 1.0)
  #
  X = zeros(Float64,nne_max*dim)
  dofs = zeros(Int64,nne_max*dim)
  Fₑ = zeros(ComplexF64,nne_max*dim)
  #
  Ψ₁ₑ = zeros(ComplexF64,nne_max*dim)
  Ψ₂ₑ = zeros(ComplexF64,nne_max*dim)
  Ψ₃ₑ = zeros(ComplexF64,nne_max*dim)
  #
  N     = zeros(Float64,nne_max)
  ∂N∂a  = zeros(Float64,(nne_max,dim))
  ∂N∂x  = zeros(Float64,(nne_max,dim))
  Jac   = zeros(Float64,(dim,dim))
  Jac⁻¹ = zeros(Float64,(dim,dim))
  #
  ∇U₁ = zeros(ComplexF64,(dim,dim))
  ∇U₂ = zeros(ComplexF64,(dim,dim))
  ∇U₃ = zeros(ComplexF64,(dim,dim))
  #
  e₁₂ = zeros(ComplexF64,(dim,dim))
  e₂₃ = zeros(ComplexF64,(dim,dim))
  e₁₃ = zeros(ComplexF64,(dim,dim))
  #
  eᵛ₁₂ = zeros(ComplexF64,dim*(dim-1))
  eᵛ₂₃ = zeros(ComplexF64,dim*(dim-1))
  eᵛ₁₃ = zeros(ComplexF64,dim*(dim-1))
  #
  σᵛ₁₂ = zeros(ComplexF64,dim*(dim-1))
  σᵛ₂₃ = zeros(ComplexF64,dim*(dim-1))
  σᵛ₁₃ = zeros(ComplexF64,dim*(dim-1))
  #
  sym∇ⁿˡ₁ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  sym∇ⁿˡ₂ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  sym∇ⁿˡ₃ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  #
  for iΩ ∈ mesh.Ω
    Dᵢⱼₖₗ = iΩ.mat.Dᵢⱼₖₗ
    for set = 1:iΩ.Sen
      etype = iΩ.Set[set]
      qr = select_quadrature_rule(etype)
      nn = iΩ.Senn[set]
      skip = iΩ.eskip[set]
      for e = 1:iΩ.ne[set]
        conn = @view iΩ.e2n[skip+1+(e-1)*nn:skip+e*nn]
        get_coor!(mesh, conn, X, nn)
        dofs!(U,nn,conn,dofs)
        #
        @inbounds for i = 1:nn*dim
          if (dofs[i]>0)
            Ψ₁ₑ[i] = Ψ₁[dofs[i]]
            Ψ₂ₑ[i] = Ψ₂[dofs[i]]
            Ψ₃ₑ[i] = Ψ₃[dofs[i]]
          else
            Ψ₁ₑ[i] = 0.0
            Ψ₂ₑ[i] = 0.0
            Ψ₃ₑ[i] = 0.0
          end
        end
        #
        integrate_H!(Fₑ,X,Dᵢⱼₖₗ,Ψ₁ₑ,Ψ₂ₑ,Ψ₃ₑ,
                    N,∂N∂a,∂N∂x,Jac,Jac⁻¹,
                    ∇U₁,∇U₂,∇U₃,
                    e₁₂,e₂₃,e₁₃,
                    eᵛ₁₂,eᵛ₂₃,eᵛ₁₃,
                    σᵛ₁₂,σᵛ₂₃,σᵛ₁₃,
                    sym∇ⁿˡ₁,sym∇ⁿˡ₂,sym∇ⁿˡ₃,
                    nn,etype,qr)
        #
        for i = 1:nn*dim
          if (dofs[i]>0)
            @inbounds Cp.nlr[dofs[i],entry] += Fₑ[i]*mult
          end
        end
        #
      end
    end
  end
  #
  return nothing
  #
end

"""
> assembly_G_nl!(Cp, entry,Ψ₁, Ψ₂,mesh, U, mult = 1.0)
It assemblies quadratic nonlinearities operator

\$ G(Ψ₁,Ψ₂) = \\frac{1}{2} ∫_{Ω} γ(Ψ₁,Ψ₂):\\mathcal{A}:ε(w) + γ(Ψ₁,w):\\mathcal{A}:ε(Ψ₂) + γ(w,Ψ₂):\\mathcal{A}:ε(Ψ₁) dΩ \$

- Cp : parametrisation data structure
- entry : entrance of the reference array
- Ψ₁ : mapping
- Ψ₂ : mapping
- mesh : mesh data structure
- U : displacement field
- mult : integral multiplier
"""
function assembly_G_nl!(Cp::Parametrisation, entry::Int64,
                        Ψ₁, Ψ₂,
                        mesh::Grid, U::Field, mult = 1.0)
  #
  X = zeros(Float64,nne_max*dim)
  dofs = zeros(Int64,nne_max*dim)
  Fₑ = zeros(ComplexF64,nne_max*dim)
  #
  Ψ₁ₑ = zeros(ComplexF64,nne_max*dim)
  Ψ₂ₑ = zeros(ComplexF64,nne_max*dim)
  #
  N     = zeros(Float64,nne_max)
  ∂N∂a  = zeros(Float64,(nne_max,dim))
  ∂N∂x  = zeros(Float64,(nne_max,dim))
  Jac   = zeros(Float64,(dim,dim))
  Jac⁻¹ = zeros(Float64,(dim,dim))
  #
  ∇U₁    = zeros(ComplexF64,(dim,dim))
  ∇U₂    = zeros(ComplexF64,(dim,dim))
  sym∇U₁ = zeros(ComplexF64,(dim,dim))
  sym∇U₂ = zeros(ComplexF64,(dim,dim))
  e₁₂    = zeros(ComplexF64,(dim,dim))
  #
  εᵛ₁ = zeros(ComplexF64,dim*(dim-1))
  εᵛ₂ = zeros(ComplexF64,dim*(dim-1))
  eᵛ₁₂ = zeros(ComplexF64,dim*(dim-1))
  #
  σᵛ₁ = zeros(ComplexF64,dim*(dim-1))
  σᵛ₂ = zeros(ComplexF64,dim*(dim-1))
  σᵛ₁₂ = zeros(ComplexF64,dim*(dim-1))
  #
  sym∇ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  sym∇ⁿˡ₁ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  sym∇ⁿˡ₂ = zeros(ComplexF64,(nne_max*dim,dim*(dim-1)))
  #
  for iΩ ∈ mesh.Ω
    Dᵢⱼₖₗ = iΩ.mat.Dᵢⱼₖₗ
    for set = 1:iΩ.Sen
      etype = iΩ.Set[set]
      qr = select_quadrature_rule(etype)
      nn = iΩ.Senn[set]
      skip = iΩ.eskip[set]
      for e = 1:iΩ.ne[set]
        conn = @view iΩ.e2n[skip+1+(e-1)*nn:skip+e*nn]
        get_coor!(mesh, conn, X, nn)
        dofs!(U,nn,conn,dofs)
        #
        @inbounds for i = 1:nn*dim
          if (dofs[i]>0)
            Ψ₁ₑ[i] = Ψ₁[dofs[i]]
            Ψ₂ₑ[i] = Ψ₂[dofs[i]]
          else
            Ψ₁ₑ[i] = 0.0
            Ψ₂ₑ[i] = 0.0
          end
        end
        #
        integrate_G!(Fₑ,X,Dᵢⱼₖₗ,Ψ₁ₑ,Ψ₂ₑ,
                     N,∂N∂a,∂N∂x,Jac,Jac⁻¹,
                     ∇U₁,∇U₂,sym∇U₁,sym∇U₂,e₁₂,
                     εᵛ₁,εᵛ₂,eᵛ₁₂,σᵛ₁,σᵛ₂,σᵛ₁₂,
                     sym∇,sym∇ⁿˡ₁,sym∇ⁿˡ₂,
                     nn,etype,qr)
        #
        for i = 1:nn*dim
          if (dofs[i]>0)
            @inbounds Cp.nlr[dofs[i],entry] += Fₑ[i]*mult
          end
        end
        #
      end
    end
  end
  #
  return nothing
  #
end

"""
> assembly_MCK!(mesh::Grid, U::Field, M::SparseMatrixCSC, C::SparseMatrixCSC, K::SparseMatrixCSC, α::Float64, β::Float64)
It assemblies mass, damping, and stiffness matrices
- mesh : Grid
- U : displacement field
- M : mass matrix
- C : damping matrix
- K : stiffness matrix
- α : Rayleigh damping mass proportional coefficient
- β : Rayleigh damping stiffness proportional coefficient
"""
function assembly_MCK!(mesh::Grid, U::Field,
                       M::SparseMatrixCSC,C::SparseMatrixCSC,K::SparseMatrixCSC,
                       α::Float64,β::Float64)
  #
  fill!(K.nzval,0.0)
  fill!(M.nzval,0.0)
  #
  X = zeros(Float64,nne_max*dim)
  dofs = zeros(Int64,nne_max*dim)
  Kₑ = zeros(Float64,(nne_max*dim,nne_max*dim))
  Mₑ = zeros(Float64,(nne_max*dim,nne_max*dim))
  #
  N = zeros(Float64,nne_max)
  ∂N∂a = zeros(Float64,(nne_max,dim))
  ∂N∂x = zeros(Float64,(nne_max,dim))
  sym∇ = zeros(Float64,(nne_max*dim,dim*(dim-1)))
  Jac = zeros(Float64,(dim,dim))
  Jac⁻¹ = zeros(Float64,(dim,dim))
  #
  for iΩ ∈ mesh.Ω
    Dᵢⱼₖₗ = iΩ.mat.Dᵢⱼₖₗ
    ρ = iΩ.mat.ρ
    for set = 1:iΩ.Sen
      etype = iΩ.Set[set]
      qr = select_quadrature_rule(etype)
      nn = iΩ.Senn[set]
      skip = iΩ.eskip[set]
      for e = 1:iΩ.ne[set]
        conn = @view iΩ.e2n[skip+1+(e-1)*nn:skip+e*nn]
        get_coor!(mesh, conn, X, nn)
        dofs!(U,nn,conn,dofs)
        integrate_MK!(Mₑ,Kₑ,X,ρ,Dᵢⱼₖₗ,N,∂N∂a,∂N∂x,sym∇,Jac,Jac⁻¹,nn,etype,qr)
        assembly_mat_sym!(K,Kₑ,dofs,nn*dim)
        assembly_mat_sym!(M,Mₑ,dofs,nn*dim)
      end
    end
  end
  #
  for i = 1:size(K.nzval)[1]
    @inbounds C.nzval[i] = α*M.nzval[i] + β*K.nzval[i]
  end
  #
  return nothing
end

"""
> assembly_mat_sym!(K::SparseMatrixCSC,Kₑ::Matrix{Float64},dofs::Array{Int64},rl::Int64)
It assembles elemental matrices into sparse matrices. Is assumes CSC format with lower triangular.
- K : global matrix
- Kₑ : elemental matrix
- dofs : degrees of freedom
- rl : row length
"""
function assembly_mat_sym!(K::SparseMatrixCSC,Kₑ::Matrix{Float64},
                           dofs::Array{Int64},rl::Int64)
  for irow = 1:rl
    if (dofs[irow]>0)
      for jcol = 1:rl
        if (dofs[jcol]>0 && dofs[irow]>=dofs[jcol])
          assembly!(K,Kₑ[irow,jcol],dofs[irow],dofs[jcol])
        end
      end
    end
  end
  return nothing
end

"""
> assembly!(K::SparseMatrixCSC, val::Float64, ir::Int64, jc::Int64)
Assembly value in sparse matrix using binary search
- K : global matrix
- val : value to assemble
- ir : row entry
- jc : column entry
"""
function assembly!(K::SparseMatrixCSC,val::Float64,ir::Int64,jc::Int64)
  @inbounds l = K.colptr[jc]-1
  @inbounds r = K.colptr[jc+1]
  while (true)
    m = (l+r)÷2
    if (K.rowval[m]>ir)
      r = m
    elseif (K.rowval[m]<ir)
      l = m
    else
      @inbounds K.nzval[m] += val
      return nothing
    end
  end
  #
  return nothing
end


text:
function ALrange(userinput; verbose=true)
    inputpass = string(replace(userinput, r"^(range)[\\:\\-\\ ]*"=>"")) |> strip
    (inputpass == "") && (inputpass = join(string.(replset.keys),","))

    for v in strip.(split(inputpass, r",| "))
        !(Symbol(v) ∈ replset.keys) && return replthrow("\n $v not found in replset!")
        verbose && println(v * " = {" * join(range(replset)[Symbol(v)], ",") * "}")
    end
end


text:
using Test
using ReverseDiff
include("../../src/noise_src_jl/smooth_max.jl")

# Inputs
k_smooth = 50
level_int = [1., 2., 3., 4., 5., 4., 3., 2., 1.]

# --- Test 1: evaluate smooth_max
@test smooth_max(k_smooth, level_int) == 5

# --- Test 2: evaluate derivatives of smooth_max
J_sol = [1.3838965267367376e-87, 7.17509597316441e-66, 3.7200759760208366e-44, 1.9287498479639178e-22, 1.0, 1.9287498479639178e-22, 3.7200759760208366e-44, 7.17509597316441e-66, 0.0]

f_smooth_max = (x)->smooth_max(k_smooth, x)
J = zeros(9)

@test ReverseDiff.gradient!(J, f_smooth_max, level_int) == J_sol

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
#----------------------------------------------------------------------------------------------#
#                                       abstract.test.jl                                       #
#----------------------------------------------------------------------------------------------#

@testset "abstract.test.jl:                                                       " begin
    # @test ...
end


text:
using Base
using Base: invpermute!!

"""Modifies perm vector in Statevector object to flip qubit at position `wire`"""
function flipqubit!(ψ::Statevector, wire::Int)
    @simd for i in 1:length(ψ)
        @inbounds ψ.perm[i] = (ψ.perm[i] - 1) ⊻ (2^(wire-1)) + 1
    end
end

"""Permutes perm vector in Statevector object, shifting qubits in "wires" to the fastest running qubit"""
function orderbit(n::Int, wires::NTuple{K,Int}, N::Int) where {K}
    result = 0
    l = N-length(wires)-1

    for i in N:-1:1
        if !(i in wires)
            bit = n >> (i-1) & 1
            result = result | (bit << l)
            l -= 1
        end
    end
    l = N-length(wires)
    for i in wires
        bit = n >> (i-1) & 1
        result = result | (bit << l)
        l += 1
    end
    return result
end

"""Permutes perm vector in Statevector object, shifting qubits in "wires" to the fastest running qubit"""
function orderqubit!(ψ::Statevector, wires::NTuple{K,Int}) where {K}
    @simd for i in 1:length(ψ)
        @inbounds ψ.perm[orderbit(i-1, wires, ψ.N) + 1] = i
    end
end

"""Swaps bit at position `p1` and `p2` integer `n`"""
function swapbit(n::Int, p1::Int, p2::Int)
    # Move p1'th to rightmost side
    bit1 = (n >> p1) & 1

    # Move p2'th to rightmost side
    bit2 = (n >> p2) & 1

    # XOR the two bits
    x = (bit1 ⊻ bit2)

    # Put the xor bit back to their original positions
    x = (x << p1) | (x << p2)

    # XOR 'x' with the original number so that the
    # two sets are swapped
    result = n ⊻ x
    return result
end

"""Permutes perm vector in Statevector object, swapping qubit at `wire1` and `wire2`"""
function swapqubit!(ψ::Statevector, wire1::Int, wire2::Int)
    if wire1 == wire2
        return
    end
    @simd for i in 1:length(ψ)
        @inbounds ψ.perm[i] = swapbit(ψ.perm[i]-1, wire1-1, wire2-1) + 1
    end
end

function resetpermute!(ψ::Statevector)
    @inbounds ψ.perm .= 1:2^ψ.N
    return
end

"""Tailored apply for XGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,XGate})
    wire = cg.iwire[1]
    flipqubit!(ψ, wire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    @inbounds ψ.state .= ψ.vec
    resetpermute!(ψ)
    return
end

"""Tailored apply for YGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,YGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)

    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = -im*ψ[x]
        else
            ψ[x] = im*ψ[x]
        end
    end

    flipqubit!(ψ, wire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    @inbounds ψ.state .= ψ.vec
    resetpermute!(ψ)
    return
end

"""Tailored apply for ZGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,ZGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)

    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = -ψ[x]
        end
    end
    return
end

"""Tailored apply for HadamardGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,HadamardGate})
    mid = 2^(ψ.N-1)
    ψ.state .*= 1/sqrt(2)

    orderqubit!(ψ, cg.iwire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    @views begin
        @inbounds ψ.state .= ψ.vec

        @inbounds ψ.state[mid+1:2mid] .*= -1

        ψ.state[1:mid] .+= ψ.vec[mid+1:end]
        ψ.state[mid+1:2mid] .+= ψ.vec[1:mid]

        @inbounds invpermute!!(ψ.state, ψ.perm)
        resetpermute!(ψ)
    end

    return
end

"""Tailored apply for SGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,SGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)

    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = im*ψ[x]
        end
    end
    return
end

"""Tailored apply for SdagGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,SdagGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)

    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = -im*ψ[x]
        end
    end
    return
end

"""Tailored apply for TGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,TGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)
    eπ4 = Base.exp(im*π/4)
    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = eπ4*ψ[x]
        end
    end
    return
end

"""Tailored apply for TdagGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,TdagGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)
    eπ4 = Base.exp(-im*π/4)
    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = eπ4*ψ[x]
        end
    end
    return
end

"""Tailored apply for RzGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,RzGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)
    θ = cg.gate.θ[]
    neπ2 = Base.exp(-im*θ/2)
    eπ2 = Base.exp(im*θ/2)
    for x in 1:2^ψ.N
        if (x-1) & bin != bin
            ψ[x] = neπ2*ψ[x]
        else
            ψ[x] = eπ2*ψ[x]
        end
    end
    return
end

"""Tailored apply for PhaseShiftGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,PhaseShiftGate})
    wire = cg.iwire[1]
    bin = 2^(wire-1)
    eϕ = Base.exp(im*cg.gate.ϕ[])
    for x in 1:2^ψ.N
        if (x-1) & bin == bin
            ψ[x] = eϕ*ψ[x]
        end
    end
    return
end

"""Tailored apply for a general single qubit gate"""
function _apply!(ψ::Statevector, cg::CircuitGate{1,T}) where {T<:AbstractGate}
    U = matrix(cg.gate)::Array{ComplexF64,2}
    orderqubit!(ψ, cg.iwire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    mid = 2^(ψ.N-length(cg.iwire))
    @views begin
        @inbounds mul!(ψ.state[1:mid], U[1,1], ψ.vec[1:mid])
        @inbounds mul!(ψ.state[mid+1:end], U[2,1], ψ.vec[1:mid])
        @inbounds mul!(ψ.vec[1:mid], U[1,2], ψ.vec[mid+1:end])
        @inbounds ψ.state[1:mid] .+= ψ.vec[1:mid]
        @inbounds mul!(ψ.vec[1:mid], U[2,2], ψ.vec[mid+1:end])
        @inbounds ψ.state[mid+1:mid+mid] .+= ψ.vec[1:mid]
        @inbounds invpermute!!(ψ.state, ψ.perm)
        resetpermute!(ψ)
    end
    return
end

"""Tailored apply for SwapGate"""
function _apply!(ψ::Statevector, cg::CircuitGate{2,SwapGate})
    i, j = cg.iwire
    i, j = i < j ? (i, j) : (j, i) # sort them
    swapqubit!(ψ, i, j)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    @inbounds ψ.state .= ψ.vec
    resetpermute!(ψ)
    return
end

"""Tailored apply for a general ControlledGate"""
function _apply(ψ::Statevector, cg::CircuitGate{M,ControlledGate{G}}) where {M,G}
    T = target_wires(cg.gate)
    C = control_wires(cg.gate)
    U = matrix(cg.gate.U)::Array{ComplexF64,2}
    orderqubit!(ψ, cg.iwire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)

    start = 2^ψ.N*(1-2^(-C))
    step = 2^(ψ.N-T-C)
    @views begin
        for i in 1:2^T
            @inbounds mul!(ψ.state[start+(i-1)*step+1:start+i*step], U[i,1], ψ.vec[start+1:start+step])
        end
        for j in 2:2^T
            for i in 1:2^T
                @inbounds mul!(ψ.vec[start+1:start+step], U[i,j], ψ.vec[start+(j-1)*step+1:start+j*step])
                @inbounds ψ.state[start+(i-1)*step+1:start+i*step] .+= ψ.vec[start+1:start+step]
            end
        end
        @inbounds invpermute!!(ψ.state, ψ.perm)
        resetpermute!(ψ)
    end
    return
end

"""
    apply(ψ::Statevector, cg::CircuitGate{M,G}) where {M,G}

Apply a [`CircuitGate`](@ref) to a quantum state vector `ψ`.
# Examples

```jldoctest
julia> cg = circuit_gate(1, HadamardGate());
julia> ψ = Statevector(ComplexF64[1 0]);
julia> apply!(ψ, cg)
julia> ψ.state
2-element Array{Complex{Float64},1}:
 0.7071067811865475 + 0.0im
 0.7071067811865475 + 0.0im
```
"""
function apply!(ψ::Statevector, cg::CircuitGate{M,G}) where {M,G}
    req_wires(cg) <= ψ.N || error("CircuitGate requires a minimum of $(req_wires(cg)) qubits, input vector `ψ` has $ψ.N qubits")
    _apply!(ψ, cg)
end

function _apply!(ψ::Statevector, cg::CircuitGate{M,G}) where {M,G}
    U = matrix(cg.gate)::Array{ComplexF64,2}
    orderqubit!(ψ, cg.iwire)
    @inbounds ψ.vec .= getindex.(Ref(ψ.state),ψ.perm)
    step = 2^(ψ.N-length(cg.iwire))
    @views begin
        for i in 1:2^M
            @inbounds mul!(ψ.state[(i-1)*step+1:i*step], U[i,1], ψ.vec[1:step])
        end
        for j in 2:2^M
            for i in 1:2^M
                @inbounds mul!(ψ.vec[1:step], U[i,j], ψ.vec[(j-1)*step+1:j*step])
                @inbounds ψ.state[(i-1)*step+1:i*step] .+= ψ.vec[1:step]
            end
        end
        @inbounds invpermute!!(ψ.state, ψ.perm)
        resetpermute!(ψ)
    end
end

"""
    apply(ψ::Statevector, cgs::Vector{<:CircuitGate})

Apply a sequence of [`CircuitGate`](@ref)(s) to a quantum state vector `ψ`.

# Examples

```jldoctest
julia> cgs = [circuit_gate(1, HadamardGate()),
                circuit_gate(1, X),
                circuit_gate(1, Y)];
julia> ψ = Statevector(ComplexF64[1 0]);
julia> apply!(ψ, cgs)
julia> ψ.state
2-element Array{Complex{Float64},1}:
 0.0 - 0.7071067811865475im
 0.0 + 0.7071067811865475im
```
"""
function apply!(ψ::Statevector, cgs::Vector{<:CircuitGate})
    req <= N || error("CircuitGates require a minimum of $req qubits, input vector `ψ` has $N qubits")
    for cg in cgs
        _apply!(ψ, cg)
    end
end

"""
    apply(ψ::Statevector, m::Moment)

returns state vector of `N` qubits after applying a `Moment{N}` object to a quantum state vector of `N` qubits `ψ`
"""
function apply!(ψ::Statevector, m::Moment)
    Nmoment = req_wires(m)
    Nmoment <= ψ.N || error("Moment affecting $Nmoment qubits applied to $N qubits")
    _apply!(ψ, m)
end

function _apply!(ψ::Statevector, m::Moment)
    for gate in m
        _apply!(ψ, gate)
    end
end

function apply!(ψ::Statevector, m::Vector{Moment})
    length(m) != 0 || error("Vector of length 0 cannot be applied")
    Nmoment = maximum(req_wires.(m))
    Nmoment <= ψ.N || error("Moment affecting $Nmoment qubits applied to $N qubits")
    for moment in m
        _apply!(ψ, moment)
    end
end

"""
    apply(ψ::Statevector, c::Circuit{N}) where {N}

returns list of expectation values from measurement operators in `c.meas` after applying circuit gates in `c.cgc` on state vector of `N` qubits `ψ`
"""
function apply!(ψ::Statevector, c::Circuit{N}) where {N}
    N == ψ.N || error("Size of vector `ψ` must match Circuit size of $(2^N)")
    length(c) != 0 || error("Circuit does not contain any gates")
    length(c.meas) != 0 || error("Circuit does not contain any measurement operators")
    for moment in c.moments
        _apply!(ψ, moment)
    end
    ψr = apply.((ψ.state,), c.meas)
    return real.(dot.((ψ.state,), ψr))
end


text:
using Documenter
using StateSpaceModels

# Set up to run docstrings with jldoctest
DocMeta.setdocmeta!(
    StateSpaceModels, :DocTestSetup, :(using StateSpaceModels); recursive=true
)

makedocs(;
    modules=[StateSpaceModels],
    doctest=true,
    clean=true,
    format=Documenter.HTML(mathengine=Documenter.MathJax2()),
    sitename="StateSpaceModels.jl",
    authors="Raphael Saavedra, Guilherme Bodin, and Mario Souto",
    pages=[
        "Home" => "index.md",
        "manual.md",
        "examples.md"
    ],
)

deploydocs(
        repo="github.com/LAMPSPUC/StateSpaceModels.jl.git",
        push_preview = true
    )


text:
#
# *** Marching Tetrahedra ***
#
# Marching Tetrahedra is an algorithm for extracting a triangular
# mesh representation of an isosurface of a scalar volumetric
# function sampled on a rectangular grid.
#
# We divide the cube into six tetrahedra. [It is possible to divide
# a cube into five tetrahedra, but not in a way that a translated
# version of the division would share face diagonals. (It reqires a
# reflection.)]
#
# Voxel corner and edge indexing conventions
#
#        Z
#        |
#
#        5------5------6              Extra edges not drawn
#       /|            /|              -----------
#      8 |           6 |              - face diagonals
#     /  9          /  10                - 13: 1 to 3
#    8------7------7   |                 - 14: 1 to 8
#    |   |         |   |                 - 15: 1 to 6
#    |   1------1--|---2  -- Y           - 16: 5 to 7
#    12 /          11 /                  - 17: 2 to 7
#    | 4           | 2                   - 18: 4 to 7
#    |/            |/                 - body diagonal
#    4------3------3                     - 19: 1 to 7
#
#  /
# X

# (X,Y,Z)-coordinates for each voxel corner ID
const voxCrnrPos = [[0 0 0],
                    [0 1 0],
                    [1 1 0],
                    [1 0 0],
                    [0 0 1],
                    [0 1 1],
                    [1 1 1],
                    [1 0 1]]'

# the voxel IDs at either end of the tetrahedra edges, by edge ID
const voxEdgeCrnrs = [[1 2],
                      [2 3],
                      [4 3],
                      [1 4],
                      [5 6],
                      [6 7],
                      [8 7],
                      [5 8],
                      [1 5],
                      [2 6],
                      [3 7],
                      [4 8],
                      [1 3],
                      [1 8],
                      [1 6],
                      [5 7],
                      [2 7],
                      [4 7],
                      [1 7]]'

# direction codes:
# 0 => +x, 1 => +y, 2 => +z,
# 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz
const voxEdgeDir = [1,0,1,0,1,0,1,0,2,2,2,2,3,4,5,3,4,5,6]

# For a pair of corner IDs, the edge ID joining them
# 0 denotes a pair with no edge
const voxEdgeIx = [[ 0  1 13  4  9 15 19 14],
                   [ 1  0  2  0  0 10 17  0],
                   [13  2  0  3  0  0 11  0],
                   [ 4  0  3  0  0  0 18 12],
                   [ 9  0  0  0  0  5 16  8],
                   [15 10  0  0  5  0  6  0],
                   [19 17 11 18 16  6  0  7],
                   [14  0  0 12  8  0  7  0]]

# voxel corners that comprise each of the six tetrahedra
const subTets = [[1 3 2 7],
                 [1 8 4 7],
                 [1 4 3 7],
                 [1 2 6 7],
                 [1 5 8 7],
                 [1 6 5 7]]'

# tetrahedron corners for each edge (indices 1-4)
const tetEdgeCrnrs = [[1 2],
                      [2 3],
                      [1 3],
                      [1 4],
                      [2 4],
                      [3 4]]'

# triangle cases for a given tetrahedron edge code
const tetTri = [[0 0 0 0 0 0],
                [1 3 4 0 0 0],
                [1 5 2 0 0 0],
                [3 5 2 3 4 5],
                [2 6 3 0 0 0],
                [1 6 4 1 2 6],
                [1 5 6 1 6 3],
                [4 5 6 0 0 0],
                [4 6 5 0 0 0],
                [1 6 5 1 3 6],
                [1 4 6 1 6 2],
                [2 3 6 0 0 0],
                [3 2 5 3 5 4],
                [1 2 5 0 0 0],
                [1 4 3 0 0 0],
                [0 0 0 0 0 0]]'

# Checks if a voxel has faces. Should be false for most voxels.
# This function should be made as fast as possible.
function hasFaces{T<:Real}(vals::Vector{T}, iso::T)
    hasFcs = false
    if vals[1] < iso
        for i = 2:8
            if vals[i] >= iso
                hasFcs = true
                break
            end
        end
    else
        for i = 2:8
            if vals[i] < iso
                hasFcs = true
                break
            end
        end
    end
    hasFcs
end

# Determines which case in the triangle table we are dealing with
function tetIx{T<:Real}(tIx::Int64, vals::Vector{T}, iso::T)
    crnrs = subTets[:,tIx]
    (vals[crnrs[1]] < iso ? 1 : 0) +
    (vals[crnrs[2]] < iso ? 2 : 0) +
    (vals[crnrs[3]] < iso ? 4 : 0) +
    (vals[crnrs[4]] < iso ? 8 : 0) + 1
end

# Determines a unique integer ID associated with the edge. This is used
# as a key in the vertex dictionary. It needs to be both unambiguous (no
# two edges get the same index) and unique (every edge gets the same ID
# regardless of which of its neighboring voxels is asking for it) in order
# for vertex sharing to be implemented properly.
function vertId(e::Int64, x::Int64, y::Int64, z::Int64,
                nx::Int64, ny::Int64)
    dx = voxCrnrPos[:,voxEdgeCrnrs[1,e]]
    voxEdgeDir[e]+7*(x-1+dx[1]+nx*(y-1+dx[2]+ny*(z-1+dx[3])))
end

# Assuming an edge crossing, determines the point in space at which it
# occurs.
function vertPos{T<:Real}(e::Int64, x::Int64, y::Int64, z::Int64,
                          vals::Vector{T}, iso::T)
    ixs = voxEdgeCrnrs[:,e]
    srcVal = float(vals[ixs[1]])
    tgtVal = float(vals[ixs[2]])
    a = (float(iso)-srcVal)/(tgtVal-srcVal)
    b = 1.0-a
    org = float([x,y,z])
    src = org+float(voxCrnrPos[:,ixs[1]])
    tgt = org+float(voxCrnrPos[:,ixs[2]])
    vrt = a*tgt+b*src
    Vertex(vrt[1],vrt[2],vrt[3])
end

# Gets the vertex ID, adding it to the vertex dictionary if not already
# present.
function getVertId{T<:Real}(e::Int64, x::Int64, y::Int64, z::Int64,
                            nx::Int64, ny::Int64,
                            vals::Vector{T}, iso::T,
                            vts::Dict{Int64,Vertex})
    vId = vertId(e,x,y,z,nx,ny)
    if !haskey(vts,vId)
        vts[vId] = vertPos(e,x,y,z,vals,iso)
    end
    vId
end

# Given a sub-tetrahedron case and a tetrahedron edge ID, determines the
# corresponding voxel edge ID.
function voxEdgeId(subTetIx::Int64, tetEdgeIx::Int64)
    tetCrnrs = tetEdgeCrnrs[:,tetEdgeIx]
    srcVoxCrnr = subTets[tetCrnrs[1],subTetIx]
    tgtVoxCrnr = subTets[tetCrnrs[2],subTetIx]
    voxEdgeIx[srcVoxCrnr,tgtVoxCrnr]
end

# Processes a voxel, adding any new vertices and faces to the given
# containers as necessary.
function procVox{T<:Real}(vals::Vector{T}, iso::T,
                          x::Int64, y::Int64, z::Int64,
                          nx::Int64, ny::Int64,
                          vts::Dict{Int64,Vertex}, fcs::Vector{Face})

    # check each sub-tetrahedron in the voxel
    for i = 1:6
        tIx = tetIx(i,vals,iso)

        for j = Range(1,3,2)
            e1 = tetTri[j,tIx]
            e2 = tetTri[j+1,tIx]
            e3 = tetTri[j+2,tIx]

            # bail if there are no more faces
            if e1 == 0 break end

            # add the face to the list
            vId(e) = getVertId(voxEdgeId(i,e),x,y,z,nx,ny,vals,iso,vts)
            push!(fcs,Face(vId(e1),vId(e2),vId(e3)))
        end
    end
end

# Given a 3D array and an isovalue, extracts a mesh represention of the
# an approximate isosurface by the method of marching tetrahedra.
function marchingTetrahedra{T<:Real}(lsf::AbstractArray{T,3},iso::T)
    vts = Dict{Int64,Vertex}()
    fcs = Array(Face,0)

    # process each voxel
    (nx,ny,nz) = size(lsf)
    for k = 1:nz-1, j = 1:ny-1, i = 1:nx-1
        vals = T[lsf[i+voxCrnrPos[1,l],j+voxCrnrPos[2,l],k+voxCrnrPos[3,l]]
                 for l = 1:8]
        if hasFaces(vals,iso)
            procVox(vals,iso,i,j,k,nx,ny,vts,fcs)
        end
    end

    (vts,fcs)
end

function isosurface(lsf,isoval)
    # get marching tetrahedra version of the mesh
    (vts,fcs) = marchingTetrahedra(lsf,isoval)

    # normalize the mesh representation
    prs = collect(vts)
    nV = size(prs,1)
    vtD = Dict{Int64,Int64}()
    for k = 1:nV
        vtD[prs[k][1]] = k
    end
    nF = size(fcs,1)
    newFace(f) = Face(vtD[f.v1],vtD[f.v2],vtD[f.v3])
    fcAry = Face[ newFace(fcs[i]) for i = 1:nF ]
    vtAry = Vertex[prs[i][2] for i = 1:nV]

    Mesh(vtAry,fcAry)
end
export isosurface


text:
using MAT,Images
export matconvnet, imgdata, make_image_grid

_mcnurl = "http://www.vlfeat.org/matconvnet/models"
_mcndir = joinpath(@__DIR__, "imagenet")

function matconvnet(name)
    global _mcncache
    if !@isdefined(_mcncache); _mcncache=Dict(); end
    if !haskey(_mcncache,name)
        matfile = "$name.mat"
        @info("Loading $matfile...")
        path = joinpath(_mcndir,matfile)
        if !isfile(path)
            println("Should I download $matfile?")
            readline()[1] == 'y' || error(:ok)
            isdir(_mcndir) || mkpath(_mcndir)
            download("$_mcnurl/$matfile",path)
        end
        _mcncache[name] = matread(path)
    end
    return _mcncache[name]
end

function imgdata(img, averageImage)
    global _imgcache
    if !@isdefined(_imgcache); _imgcache = Dict(); end
    if !haskey(_imgcache,img)
        if occursin("://",img)
            @info("Downloading $img")
            a0 = load(download(img))
        else
            a0 = load(img)
        end
        new_size = ntuple(i->div(size(a0,i)*224,minimum(size(a0))),2)
        a1 = Images.imresize(a0, new_size)
        i1 = div(size(a1,1)-224,2)
        j1 = div(size(a1,2)-224,2)
        b1 = a1[i1+1:i1+224,j1+1:j1+224]
        # ad-hoc solution for Mac-OS image
        macfix = convert(Array{FixedPointNumbers.Normed{UInt8,8},3}, channelview(b1))
        c1 = permutedims(macfix, (3,2,1))
        d1 = convert(Array{Float32}, c1)
        e1 = reshape(d1[:,:,1:3], (224,224,3,1))
        f1 = (255 * e1 .- averageImage)
        g1 = permutedims(f1, [2,1,3,4])
        _imgcache[img] = g1
    end
    return _imgcache[img]
end

function make_image_grid(images; gridsize=(8,8), scale=2.0, height=28, width=28)
    shape = (height, width)
    nchannels = size(first(images))[end]
    @assert nchannels == 1 || nchannels == 3
    shp = map(x->Int(round(x*scale)), shape)
    y = map(x->Images.imresize(x,shp), images)
    gridx, gridy = gridsize
    outdims = (gridx*shp[1]+gridx+1,gridy*shp[2]+gridy+1)
    out = zeros(outdims..., nchannels)
    for k = 1:gridx+1; out[(k-1)*(shp[1]+1)+1,:,:] .= 1.0; end
    for k = 1:gridy+1; out[:,(k-1)*(shp[2]+1)+1,:] .= 1.0; end

    x0 = y0 = 2
    for k = 1:length(y)
        x1 = x0+shp[1]-1
        y1 = y0+shp[2]-1
        out[x0:x1,y0:y1,:] .= y[k]

        y0 = y1+2
        if k % gridy == 0
            x0 = x1+2
            y0 = 2
        else
            y0 = y1+2
        end
    end

    out = convert(Array{Float64}, map(x->isnan(x) ? 0 : x, out))
    if nchannels == 1
        out = reshape(out, (size(out,1),size(out,2)))
        out = permutedims(out, (2,1))
    else
        out = permutedims(out, (3,1,2))
    end
    return out
end


text:
using Base: @deprecate
using TimeZones: TimeZones

# BEGIN LaxZonedDateTimes 1.0 deprecations

if isdefined(TimeZones, :localtime) && isdefined(TimeZones, :utc)
    import TimeZones: localtime, utc
end

@deprecate localtime(lzdt::LaxZonedDateTime) DateTime(lzdt) true
@deprecate utc(lzdt::LaxZonedDateTime) DateTime(lzdt, UTC) true

@deprecate DateTime(lzdt::LaxZonedDateTime, ::Type{Local}) DateTime(lzdt)
@deprecate Date(lzdt::LaxZonedDateTime, ::Type{Local}) Date(lzdt)
@deprecate Time(lzdt::LaxZonedDateTime, ::Type{Local}) Time(lzdt)

# END LaxZonedDateTimes 1.0 deprecations


text:
#Contains all the evaluation stuff
using JLD2
using PyPlot
using FileIO
using Statistics

d  = load("Results/10_29_evaluate_seg/1-2_all.jld2")["RESULT"]
e  = load("Results/10_29_evaluate_seg/1-2_ppos.jld2")["RESULT"]

a= load("Results/10_29_evaluate_seg/3_all.jld2")["result"]
b= load("Results/10_29_evaluate_seg/4_all.jld2")["result"]
c= load("Results/10_29_evaluate_seg/5_all.jld2")["result"]

a_all = 5
iters = 500

RE  = zeros(a_all,iters)
A2C = zeros(a_all,iters)
TNPG= zeros(a_all,iters)
TRPO= zeros(a_all,iters)
PPO = zeros(a_all,iters)
TD3 = zeros(a_all,iters)

RE[1,:]   = a[1,:]
A2C[1,:]  = a[2,:]
TNPG[1,:] = a[3,:]
TRPO[1,:] = a[4,:]
PPO[1,:]  = a[5,:]
TD3[1,:]  = a[6,:]

RE[2,:]   = b[1,:]
A2C[2,:]  = b[2,:]
TNPG[2,:] = b[3,:]
TRPO[2,:] = b[4,:]
PPO[2,:]  = b[5,:]
TD3[2,:]  = b[6,:]

RE[3,:]   = c[1,:]
A2C[3,:]  = c[2,:]
TNPG[3,:] = c[3,:]
TRPO[3,:] = c[4,:]
PPO[3,:]  = c[5,:]
TD3[3,:]  = c[6,:]

RE[4,:]   = d[1:6:3000]
A2C[4,:]  = d[2:6:3000]
TNPG[4,:] = d[3:6:3000]
TRPO[4,:] = d[4:6:3000]
PPO[4,:]  = e[5:6:3000]
TD3[4,:]  = e[6:6:3000]

RE[5,:]   = d[3001:6:end]
A2C[5,:]  = d[3002:6:end]
TNPG[5,:] = d[3003:6:end]
TRPO[5,:] = d[3004:6:end]
PPO[5,:]  = e[3005:6:end]
TD3[5,:]  = e[3006:6:end]

scale = 1
vis   = 0.5

function f_filt(x)
    l = length(x)
    tau = 0.1
    y = zeros(l)
    y[1] = x[1]

    for i = 2:l
        y[i] = (1-tau)*y[i-1] + tau *x[i]
    end

    return y
end

offset = 0

re_mean = f_filt(vec(mean(RE,dims = 1)) .+ offset)
re_std  = f_filt(scale*vec(std(RE,dims=1)))

a2c_mean = f_filt(vec(mean(A2C,dims = 1)) .+ offset)
a2c_std  = f_filt(scale*vec(std(A2C,dims=1)))

tnpg_mean = f_filt(vec(mean(TNPG,dims = 1)) .+ offset)
tnpg_std  = f_filt(scale*vec(std(TNPG,dims=1)))

trpo_mean = f_filt(vec(mean(TRPO,dims = 1)) .+ offset)
trpo_std  = f_filt(scale*vec(std(TRPO,dims=1)))

ppo_mean = f_filt(vec(mean(PPO,dims = 1)) .+ offset)
ppo_std  = f_filt(scale*vec(std(PPO,dims=1)))

td3_mean = f_filt(vec(mean(TD3,dims = 1)) .+ offset)
td3_std  = f_filt(scale*vec(std(TD3,dims=1)))

figure("seg")
clf()
grid(alpha=0.25)
fill_between(1:iters,re_mean .+ re_std, re_mean .- re_std, alpha=vis)
plot(re_mean, label="REINFORCE")

#fill_between(1:iters,a2c_mean .+ a2c_std, a2c_mean .- a2c_std, alpha=vis)
#plot(a2c_mean, label="A2C")
plot(0,0)
fill_between(1:iters,0,0,0)
#fill_between(1:iters,tnpg_mean .+ tnpg_std, tnpg_mean .- tnpg_std, alpha=vis)
#plot(tnpg_mean, label="TNPG")
plot(0,0)
fill_between(1:iters,0,0,0)

fill_between(1:iters,trpo_mean .+ trpo_std, trpo_mean .- trpo_std, alpha=vis)
plot(trpo_mean, label="TRPO")
fill_between(1:iters,ppo_mean .+ ppo_std, ppo_mean .- ppo_std, alpha=vis)
plot(ppo_mean, label="PPO")
fill_between(1:iters,td3_mean .+ td3_std, td3_mean .- td3_std, alpha=vis)
plot(td3_mean, label="PPO-soft")
xlabel("Iterations", fontsize = 16)
#ylabel("Cost", fontsize = 16)
#legend()
ylim([0,15])

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA