text:
using .Networks

using Dates
using DataStructures
using Printf

"""
    Determines whether a decision should be :
    - already decided : DECIDED
    - to decide firmly (setting a common value for all scenarios) : TO_DECIDE
    - to decide freely (possibly setting different values for different scenarios): FREE
    The decision is based on the characteristic time period `delta` (delta can represent the DMO or DP)
"""
function compute_firmness(ech::Dates.DateTime, next_ech::Union{Nothing,Dates.DateTime},
                        ts::Dates.DateTime, delta::Dates.Period)
    if ( !isnothing(next_ech) && (next_ech < ech) )
        throw( error("next_ech (", next_ech, ") must be later than ech (", ech,").") )
    end

    final_decision_time = ts - delta

    if final_decision_time < ech
        return DECIDED
    elseif ( isnothing(next_ech) || (final_decision_time < next_ech) )
        return TO_DECIDE
    else
        return FREE
    end
end

function compute_firmness(ech::Dates.DateTime, next_ech::Union{Nothing,Dates.DateTime},
                    TS::Vector{Dates.DateTime}, generators::Vector{Networks.Generator})
    firmness = Firmness()
    for generator in generators
        gen_id = Networks.get_id(generator)
        dmo = Networks.get_dmo(generator)
        dp = Networks.get_dp(generator)

        for ts in TS
            if Networks.get_p_min(generator) > eps()
                #commitment
                commitment_firmness = compute_firmness(ech, next_ech, ts, dmo)
                set_commitment_firmness!(firmness, gen_id, ts, commitment_firmness)
            end

            #power level
            power_level_firmness = compute_firmness(ech, next_ech, ts, dp)
            set_power_level_firmness!(firmness, gen_id, ts, power_level_firmness)
        end
    end
    return firmness
end

#########################################
#       Schedule Firmness
#########################################

function verify_firmness(firmness::Firmness, schedule::Schedule;
                        excluded_ids::Union{Set{String},Vector{String}}=Set{String}())
    excluded_ids = Set{String}(excluded_ids)
    return ( verify_commitment_firmness(get_commitment_firmness(firmness), schedule, excluded_ids=excluded_ids)
            && verify_production_firmness(get_power_level_firmness(firmness), schedule, excluded_ids=excluded_ids) )
end

function verify_commitment_firmness(firmness::SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} },
                                    schedule::Schedule;
                                    excluded_ids::Set{String}=Set{String}())
    for (gen_id, generator_firmness) in firmness
        if gen_id in excluded_ids
            continue
        end
        if !verify_firmness(generator_firmness, get_sub_schedule(schedule, gen_id).commitment)
            @warn(@sprintf("commitment schedule of generator %s violates firmness", gen_id))
            return false
        end
    end
    return true
end

function verify_production_firmness(firmness::SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} },
                                    schedule::Schedule;
                                    excluded_ids::Set{String}=Set{String}())
    for (gen_id, generator_firmness) in firmness
        if gen_id in excluded_ids
            continue
        end
        if !verify_firmness(generator_firmness, get_sub_schedule(schedule, gen_id).production)
            @warn(@sprintf("production schedule of generator %s violates firmness", gen_id))
            return false
        end
    end
    return true
end

"""
    Verfies if a given schedule respects the firmness constraints at a given timestep
"""
function verify_firmness(generator_firmness::SortedDict{Dates.DateTime, DecisionFirmness},
                        scheduled_values::SortedDict{Dates.DateTime, UncertainValue{T}}) where T

    for (ts, decision_firmness) in generator_firmness
        #if firmness is defined, a schedule value must exist
        if !haskey(scheduled_values, ts)
            @warn(@sprintf("firmness violation : missing required schedule value at timestep %s", ts))
            return false
        end

        uncertain_value = scheduled_values[ts]
        if ((decision_firmness == TO_DECIDE) || (decision_firmness == DECIDED))
            if !is_definitive(uncertain_value)
                @warn(@sprintf("firmness violation : non-definitive value at timestep %s while a firm value is required", ts))
                return false
            end

        elseif (decision_firmness == FREE)
            # FIXME: are definitive values accepted as FREE ? (might be useful for EnergyMarketAtFO, which considers one scenario)
            # if (is_definitive(uncertain_value))
            #     @warn(@sprintf("firmness violation : definitive value at timestep %s while a non-firm value is required", ts))
            #     return false
            # end
            if (is_missing_values(uncertain_value))
                @warn(@sprintf("firmness violation : missing value at timestep %s while a non-firm value is required", ts))
                return false
            end

        else
            throw( error("Non-handled DecisionFirmness value : ", decision_firmness) )
        end
    end

    return true
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using Fluxion
using Base.Test

@testset "Helpers" begin
    @testset "Sequence" begin
        @testset "correctness" begin

            seq = Fluxion.event_sequence(Dict(
                :a => [x / 2 for x in 0:10],
                :b => 0:4,
                :c => [1.5, 2.5, 4]
                ))

            res = []

            t = 0
            while true
                events, seq = Fluxion.event_sequence_next_until(t, seq)
                push!(res, (t, events))
                t += 1
                if Fluxion.event_sequence_done(seq)
                    break
                end
            end

            res_ref = [
                (0, [(0.0, [:a, :b])]),
                (1, [(0.5, [:a]), (1.0, [:a, :b])]),
                (2, [(1.5, [:a, :c]), (2.0, [:a, :b])]),
                (3, [(2.5, [:a, :c]), (3.0, [:a, :b])]),
                (4, [(3.5, [:a]), (4.0, [:a, :b, :c])]),
                (5, [(4.5, [:a]), (5.0, [:a])])]

            @test res == res_ref

        end
    end
end

ext: jl


text:
function main()

    io = open("xenon.txt", "w");
    io_p = open("xenon_pressure.txt", "w");

    B1 = 103701.61e-8
    B2 = 31228.61e-8
    C1 = 12.75e-3*(1e-6)^2
    C2 = 0.561e-3*(1e-6)^2

    p0 = 1
    k0(λ) = 2π/λ
    a = 14.9e-6
    T = 293
    T0 = 273
    λ0 = 800e-9
    function δ(λ, T=T)
        return T0 / T * (B1 * λ^2/(λ^2 - C1) + B2 * λ^2 / (λ^2 - C2))
    end

    function n(λ, p, u, T=T)
        return sqrt(1 + δ(λ, T)*p/p0 - (u / k0(λ) / a)^2 )
    end

    e11 = 2.405
    e12 = 5.52
    e13 = 8.6573
    e23 = 10.1735
    e33 = 11.6198
    e43 = 13.0152

    ps = range(0, 4, length=1000)
    λs = range(400e-9, 1000e-9, length=1000)

    for p in ps
        res13 = (n(λ0, p, e11) - n(λ0/3, p, e13)) * 10^4
        res23 = (n(λ0, p, e11) - n(λ0/3, p, e23)) * 10^4
        res33 = (n(λ0, p, e11) - n(λ0/3, p, e33)) * 10^4
        write(io_p, "$p \t $res13 \t $res23 \t $res33 \n")
    end

    p = 0.96
    for x in λs
        i = x * 10^9
        res11 = (n(x, p, e11) - n(x/3, p, e11)) * 10^4
        res12 = (n(x, p, e11) - n(x/3, p, e12)) * 10^4
        res13 = (n(x, p, e11) - n(x/3, p, e13)) * 10^4
        res23 = (n(x, p, e11) - n(x/3, p, e23)) * 10^4
        res33 = (n(x, p, e11) - n(x/3, p, e33)) * 10^4
        write(io, "$i \t $res11 \t $res12 \t $res13 \t $res23 \t $res33 \n")
    end

    close(io)
    close(io_p)
end
main()


text:
module TestGAT

using Base.Meta: ParseError
using Test
using DataStructures: OrderedDict

import Catlab: GAT
using Catlab.GAT

# GAT expressions
#################

# Raw expressions
@test GAT.parse_raw_expr(:(Ob)) == :Ob
@test GAT.parse_raw_expr(:(Hom(X,Y))) == :(Hom(X,Y))
@test_throws ParseError GAT.parse_raw_expr(:("Ob"))
@test_throws ParseError GAT.parse_raw_expr(:(Hom(X,0)))

# Contexts
@test (GAT.parse_context(:((X::Ob, Y::Ob))) ==
       GAT.Context((:X => :Ob, :Y => :Ob)))
@test (GAT.parse_context(:((X::Ob, Y::Ob, f::Hom(X,Y)))) ==
       GAT.Context((:X => :Ob, :Y => :Ob, :f => :(Hom(X,Y)))))
@test GAT.parse_context(:(())) == GAT.Context()
@test_throws ParseError GAT.parse_context(:((X::Ob, X::Ob))) # Repeat variables

# Type constructor
expr = :(Ob::TYPE)
cons = GAT.TypeConstructor(:Ob, [], GAT.Context())
@test GAT.parse_constructor(expr) == cons

expr = (quote "Object" Ob::TYPE end).args[2]
cons = GAT.TypeConstructor(:Ob, [], GAT.Context(), "Object")
@test GAT.parse_constructor(expr) == cons

expr = :(Hom(X,Y)::TYPE <= (X::Ob, Y::Ob))
context = GAT.Context((:X => :Ob, :Y => :Ob))
cons = GAT.TypeConstructor(:Hom, [:X,:Y], context)
@test GAT.parse_constructor(expr) == cons

# Term constructor
expr = :(unit()::Ob)
cons = GAT.TermConstructor(:unit, [], :Ob, GAT.Context())
@test GAT.parse_constructor(expr) == cons

expr = (quote "Monoidal unit" munit()::Ob end).args[2]
cons = GAT.TermConstructor(:munit, [], :Ob, GAT.Context(), "Monoidal unit")
@test GAT.parse_constructor(expr) == cons

cons = GAT.TermConstructor(:id, [:X], :(Hom(X,X)), GAT.Context(:X => :Ob))
@test GAT.parse_constructor(:(id(X)::Hom(X,X) <= (X::Ob))) == cons
@test GAT.parse_constructor(:(id(X::Ob)::Hom(X,X))) == cons

expr = :(compose(f,g)::Hom(X,Z) <= (X::Ob, Y::Ob, Z::Ob, f::Hom(X,Y), g::Hom(Y,Z)))
context = GAT.Context((:X => :Ob, :Y => :Ob, :Z => :Ob,
                       :f => :(Hom(X,Y)), :g => :(Hom(Y,Z))))
cons = GAT.TermConstructor(:compose, [:f,:g], :(Hom(X,Z)), context)
@test GAT.parse_constructor(expr) == cons
expr = :(compose(f::Hom(X,Y), g::Hom(Y,Z))::Hom(X,Z) <= (X::Ob, Y::Ob, Z::Ob))
@test GAT.parse_constructor(expr) == cons

# Type transformations
bindings = Dict((:Ob => :Obj, :Hom => :Mor))
cons = GAT.TypeConstructor(:Hom, [:X,:Y],
  GAT.Context((:X => :Ob, :Y => :Ob)))
target = GAT.TypeConstructor(:Mor, [:X,:Y],
  GAT.Context((:X => :Obj, :Y => :Obj)))
@test GAT.replace_types(bindings, cons) == target

cons = GAT.TermConstructor(:compose, [:f,:g], :(Hom(X,Z)),
  GAT.Context((:X => :Obj, :Y => :Obj, :Z => :Obj,
               :f => :(Hom(X,Y)), :g => :(Hom(Y,Z)))))
target = GAT.TermConstructor(:compose, [:f,:g], :(Mor(X,Z)),
  GAT.Context((:X => :Obj, :Y => :Obj, :Z => :Obj,
               :f => :(Mor(X,Y)), :g => :(Mor(Y,Z)))))
@test GAT.replace_types(bindings, cons) == target

@test GAT.strip_type(:Ob) == :Ob
@test GAT.strip_type(:(Hom(X,Y))) == :Hom
@test GAT.strip_type(:(Hom(dual(X),dual(Y)))) == :Hom

# Signatures
############

""" Signature of the theory of categories
"""
@signature Category(Ob,Hom) begin
  Ob::TYPE
  Hom(dom, codom)::TYPE <= (dom::Ob, codom::Ob)

  id(X)::Hom(X,X) <= (X::Ob)
  compose(f,g)::Hom(X,Z) <= (X::Ob, Y::Ob, Z::Ob, f::Hom(X,Y), g::Hom(Y,Z))
end

@test isa(Category, Module)
@test occursin("theory of categories", string(Docs.doc(Category)))
@test sort(names(Category)) == sort([:Category, :Ob, :Hom])
@test isa(Category.Ob, Type) && isa(Category.Hom, Type)
@test isa(dom, Function) && isa(codom, Function)
@test isa(id, Function) && isa(compose, Function)

# Manually constructed signature of theory of categories
types = [
  GAT.TypeConstructor(:Ob, [], GAT.Context()),
  GAT.TypeConstructor(:Hom, [:dom,:codom],
    GAT.Context((:dom => :Ob, :codom => :Ob))),
]
terms = [
  GAT.TermConstructor(:id, [:X], :(Hom(X,X)), GAT.Context(:X => :Ob)),
  GAT.TermConstructor(:compose, [:f,:g], :(Hom(X,Z)),
    GAT.Context((:X => :Ob, :Y => :Ob, :Z => :Ob,
                 :f => :(Hom(X,Y)), :g => :(Hom(Y,Z))))),
]
category_signature = GAT.Signature(types, terms)

@test Category.class().signature == category_signature

""" Equivalent shorthand definition of Category signature
"""
@signature CategoryAbbrev(Ob,Hom) begin
  Ob::TYPE
  Hom(dom::Ob, codom::Ob)::TYPE

  id(X::Ob)::Hom(X,X)
  compose(f::Hom(X,Y),g::Hom(Y,Z))::Hom(X,Z) <= (X::Ob, Y::Ob, Z::Ob)
end

@test CategoryAbbrev.class().signature == category_signature

# Methods for signature
accessors = [ GAT.JuliaFunction(:(dom(::Hom)), :Ob),
              GAT.JuliaFunction(:(codom(::Hom)), :Ob) ]
constructors = [ GAT.JuliaFunction(:(id(X::Ob)), :Hom),
                 GAT.JuliaFunction(:(compose(f::Hom, g::Hom)), :Hom) ]
@test GAT.accessors(Category.class().signature) == accessors
@test GAT.constructors(Category.class().signature) == constructors
@test GAT.interface(Category.class()) == [accessors; constructors]

# Signature extension
@signature Semigroup(S) begin
  S::TYPE
  times(x::S,y::S)::S
end

@signature Semigroup(M) => MonoidExt(M) begin
  munit()::M
end

@test isa(Semigroup, Module) && isa(MonoidExt, Module)
@test length(methods(times)) == 2 # Semigroup.S, MonoidExt.M
@test length(methods(munit)) == 1 # MonoidExt.M

signature = GAT.Signature(
  [ GAT.TypeConstructor(:M, [], GAT.Context()) ],
  [ GAT.TermConstructor(:times, [:x,:y], :M,
      GAT.Context((:x => :M, :y => :M))),
    GAT.TermConstructor(:munit, [], :M, GAT.Context()) ],
)
@test MonoidExt.class().signature == signature

# GAT expressions in a signature
################################

sig = Category.class().signature
context = GAT.Context((:X => :Ob, :Y => :Ob, :Z => :Ob,
                       :f => :(Hom(X,Y)), :g => :(Hom(Y,Z))))
@test GAT.expand_in_context(:X, [:f,:g], context, sig) == :(dom(f))
@test (GAT.expand_in_context(:(Hom(X,Z)), [:f,:g], context, sig) ==
       :(Hom(dom(f),codom(g))))

context = GAT.Context((:X => :Ob, :Y => :Ob, :Z => :Ob, :f => :(Hom(X,Y))))
@test_throws ErrorException GAT.expand_in_context(:W, [:f], context, sig)
@test_throws ErrorException GAT.expand_in_context(:Z, [:f], context, sig)

context = GAT.Context((:X => :Ob, :Y => :Ob, :f => :(Hom(X,Y))))
@test GAT.equations(context, sig) == [ :(dom(f)) => :X, :(codom(f)) => :Y ]
@test GAT.equations([:f], context, sig) == []

context = GAT.Context((:X => :Ob, :Y => :Ob, :Z => :Ob,
                       :f => :(Hom(X,Y)), :g => :(Hom(Y,Z))))
@test (GAT.equations(context, sig) ==
       [ :(dom(f)) => :X, :(codom(f)) => :Y,
         :(dom(g)) => :Y, :(codom(g)) => :Z ])
@test GAT.equations([:f,:g], context, sig) == [ :(dom(g)) => :(codom(f)) ]

# Instances
###########

""" Vectors as an instance of the theory of semigroups
"""
@instance Semigroup(Vector) begin
  times(x::Vector, y::Vector) = [x; y]
end

@test times([1,2],[3,4]) == [1,2,3,4]

@signature Monoid(M) begin
  M::TYPE
  munit()::M
  times(x::M,y::M)::M

  times(xs::Vararg{M}) = foldl(times, xs)
end

# Incomplete instance of Monoid
@test_throws ErrorException @instance Monoid(String) begin
  times(x::AbsStringtractString, y::String) = string(x,y)
end

# Complete instance of Monoid
@instance Monoid(String) begin
  munit(::Type{String}) = ""
  times(x::String, y::String) = string(x,y)
end

@test length(methods(munit)) == 3 # Monoid, MonoidExt, String
@test munit(String) == ""
@test times("a", "b") == "ab"
@test times("a", "b", "c") == "abc"

# Reflection
@test invoke_term(Monoid, (String,), :munit) == ""
@test invoke_term(Monoid, (String,), :times, "a", "b") == "ab"
@test invoke_term(Monoid, (String,), :times, "a", "b", "c") == "abc"

end


text:
module reproSim
using Pkg
Pkg.add("Distributions")
Pkg.add("StatsBase")
Pkg.add("CSV")
Pkg.add("DataFrames")
using Distributions
using StatsBase
using CSV
using DataFrames
function getPop(n)
    """return array of length n. with index, x, y chomosome values"""
return hcat(1:n,rand(Normal(0,0.5),n),rand(Normal(0,0.5),n))
end

function fselect(pop,n)
    """"returns population with additional column for fecundity"""
        # get samples
    fitList = sample(pop[:,1], Weights(pop[:, 4]), n, replace=true)
    #add fecundity column to pop.
 pop = hcat(pop,[sum(i .== fitList) for i in pop[:,1]])
 #return pop. without non-reproducers
return pop[pop[:,5] .!= 0.6,:]
end

function sselect(pop,n)
    """" return half of pop. weighted by fitness"""

return pop[Int.(sample(Int.(1:n), Weights(pop[:, 5]), Int(n/2), replace=false)),:]
end

function ggselect(pop,n)
    """" return half of pop. weighted by fitness"""

return pop[Int.(sample(Int.(1:n/2), Weights(pop[:, 1]), Int(n/2), replace=true)),:]
end

function clone(pop,n,mut)
    x = []
    y = []

    for mother in 1:size(pop)[1] # for each mother

        for fecundity in 1:pop[mother,7] #according to fecundity
            push!(x,pop[mother,2])
            push!(y,pop[mother,3])

        end
    end

    return hcat(1:n,x .+ rand(Normal(0,mut),n),y .+ rand(Normal(0,mut),n))

end

function sex(pop,fathers,n,mut)
    x = []
    y = []

    for mother in 1:size(pop)[1] # for each mother
        father = fathers[mother,:] # pick a father

        for fecundity in 1:pop[mother,7] #according to fecundity
            push!(x,pop[mother,rand(2:3)])
            push!(y,father[rand(2:3)])
        end
    end

    return hcat(1:n,x .+ rand(Normal(0,mut),n), y .+ rand(Normal(0,mut),n))
end

function fitPop(pop,f,s,gg,theta)
    z = (pop[:,2] .+ pop[:,3]) ./2

    for c in [1-f,1-s,1-gg]
        pop= hcat(pop, c .+ ((1 .-c) .*exp.((-1/2) .*(z .- theta) .^2)))
    end

    return pop
end

function getW(pop,theta)
    z = pop[:,2] .+ pop[:,3] ./2

    return [mean(exp.((-1/2) .*(z .- theta) .^2)) var(z)]
end

function sim(vartheta,lambda,rMode,mut,n,generations,f,s,gg)
    theta = 0
    noise = vartheta*(1-(lambda^2))
    pop = getPop(n)
    pop = fitPop(pop,f,s,gg, theta)
    dat = Array{Float64}(undef, 0, 2)
    r = Int(n/2)

    for t in 1:generations

        pop = sselect(pop,n)

        if rMode == "clone"
            pop = fselect(pop,n)
            pop = clone(pop,n,mut)
        else
            fathers = ggselect(pop,n)
            pop = fselect(pop,n)
            pop = sex(pop,fathers,n,mut)
        end

        theta = theta*lambda +rand(Normal(0,noise))
        pop = fitPop(pop,f,s,gg, theta)
        dat = vcat(dat,getW(pop,theta))
    end

    wg = geomean(dat[:,1])
    wa = mean(dat[:,1])
    meanv = mean(dat[:,2])
    vmean = var(dat[:,1])

    return [wg wa meanv vmean vartheta lambda rMode mut n generations f s gg]
end

end # module


text:
using Flux: onehot, onehotbatch, onecold

@show onehot(1,[1,2,3])
@show onehotbatch([1,2], [1,2,3])
@show onecold([0.3,0.2,0.5],[1,2,3])

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Math function mappings to Metal intrinsics

using Base: FastMath

# Override Base functions with device counterparts
# Or if no Base equivalent, create device function
# The lack of metaprogramming here is for your benefit, dear user

# TODO: Add support for vector types

### Floating Point Intrinsics

## Metal only supports single and half-precision floating-point types (and their vector counterparts)
## For single precision types, there are precise and fast variants

# TODO: clamp, mix, saturate, sign, smoothstep, step

@device_override FastMath.abs_fast(x::Float32) = ccall("extern air.fast_fabs.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.abs(x::Float32) = ccall("extern air.fabs.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.abs(x::Float16) = ccall("extern air.fabs.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.acos_fast(x::Float32) = ccall("extern air.fast_acos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.acos(x::Float32) = ccall("extern air.acos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.acos(x::Float16) = ccall("extern air.acos.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.acosh_fast(x::Float32) = ccall("extern air.fast_acosh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.acosh(x::Float32) = ccall("extern air.acosh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.acosh(x::Float16) = ccall("extern air.acosh.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.asin_fast(x::Float32) = ccall("extern air.fast_asin.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.asin(x::Float32) = ccall("extern air.asin.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.asin(x::Float16) = ccall("extern air.asin.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.asinh_fast(x::Float32) = ccall("extern air.fast_asinh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.asinh(x::Float32) = ccall("extern air.asinh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.asinh(x::Float16) = ccall("extern air.asinh.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.atan_fast(x::Float32) = ccall("extern air.fast_atan.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.atan(x::Float32) = ccall("extern air.atan.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.atan(x::Float16) = ccall("extern air.atan.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.atanh_fast(x::Float32) = ccall("extern air.fast_atanh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.atanh(x::Float32) = ccall("extern air.atanh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.atanh(x::Float16) = ccall("extern air.atanh.f16", llvmcall, Float16, (Float16,), x)

@device_function ceil_fast(x::Float32) = ccall("extern air.fast_ceil.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.ceil(x::Float32) = ccall("extern air.ceil.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.ceil(x::Float16) = ccall("extern air.ceil.f16", llvmcall, Float16, (Float16,), x)

@device_function copysign_fast(x::Float32) = ccall("extern air.fast_copysign.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.copysign(x::Float32) = ccall("extern air.copysign.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.copysign(x::Float16) = ccall("extern air.copysign.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.cos_fast(x::Float32) = ccall("extern air.fast_cos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cos(x::Float32) = ccall("extern air.cos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cos(x::Float16) = ccall("extern air.cos.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.cosh_fast(x::Float32) = ccall("extern air.fast_cosh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cosh(x::Float32) = ccall("extern air.cosh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cosh(x::Float16) = ccall("extern air.cosh.f16", llvmcall, Float16, (Float16,), x)

@device_function cospi_fast(x::Float32) = ccall("extern air.fast_cospi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cospi(x::Float32) = ccall("extern air.cospi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.cospi(x::Float16) = ccall("extern air.cospi.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.exp_fast(x::Float32) = ccall("extern air.fast_exp.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp(x::Float32) = ccall("extern air.exp.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp(x::Float16) = ccall("extern air.exp.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.exp2_fast(x::Float32) = ccall("extern air.fast_exp2.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp2(x::Float32) = ccall("extern air.exp2.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp2(x::Float16) = ccall("extern air.exp2.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.exp10_fast(x::Float32) = ccall("extern air.fast_exp10.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp10(x::Float32) = ccall("extern air.exp10.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.exp10(x::Float16) = ccall("extern air.exp10.f16", llvmcall, Float16, (Float16,), x)

@device_function floor_fast(x::Float32) = ccall("extern air.fast_floor.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.floor(x::Float32) = ccall("extern air.floor.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.floor(x::Float16) = ccall("extern air.floor.f16", llvmcall, Float16, (Float16,), x)

# half/regular -> air.fma.f16
# half/(precise or fast) -> air.fma.f32
@device_override Base.fma(a::Float32, b::Float32, c::Float32) = ccall("extern air.fma.f32", llvmcall, Float32, (Float32,Float32,Float32,), a,b,c)
@device_override Base.fma(a::Float16, b::Float16, c::Float16) = ccall("extern air.fma.f32", llvmcall, Float16, (Float16,Float16,Float16,), a,b,c)

@device_function fract_fast(x::Float32) = ccall("extern air.fast_fract.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function fract(x::Float32) = ccall("extern air.fract.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function fract(x::Float16) = ccall("extern air.fract.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.log_fast(x::Float32) = ccall("extern air.fast_log.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log(x::Float32) = ccall("extern air.log.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log(x::Float16) = ccall("extern air.log.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.log2_fast(x::Float32) = ccall("extern air.fast_log2.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log2(x::Float32) = ccall("extern air.log2.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log2(x::Float16) = ccall("extern air.log2.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.log10_fast(x::Float32) = ccall("extern air.fast_log10.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log10(x::Float32) = ccall("extern air.log10.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.log10(x::Float16) = ccall("extern air.log10.f16", llvmcall, Float16, (Float16,), x)

@device_function rint_fast(x::Float32) = ccall("extern air.fast_rint.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function rint(x::Float32) = ccall("extern air.rint.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function rint(x::Float16) = ccall("extern air.rint.f16", llvmcall, Float16, (Float16,), x)

@device_function round_fast(x::Float32) = ccall("extern air.fast_round.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.round(x::Float32) = ccall("extern air.round.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.round(x::Float16) = ccall("extern air.round.f16", llvmcall, Float16, (Float16,), x)

@device_function rsqrt_fast(x::Float32) = ccall("extern air.fast_rsqrt.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function rsqrt(x::Float32) = ccall("extern air.rsqrt.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function rsqrt(x::Float16) = ccall("extern air.rsqrt.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.sin_fast(x::Float32) = ccall("extern air.fast_sin.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sin(x::Float32) = ccall("extern air.sin.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sin(x::Float16) = ccall("extern air.sin.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.sincos_fast(x::Float32) = ccall("extern air.fast_sincos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sincos(x::Float32) = ccall("extern air.sincos.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sincos(x::Float16) = ccall("extern air.sincos.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.sinh_fast(x::Float32) = ccall("extern air.fast_sinh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sinh(x::Float32) = ccall("extern air.sinh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sinh(x::Float16) = ccall("extern air.sinh.f16", llvmcall, Float16, (Float16,), x)

@device_function sinpi_fast(x::Float32) = ccall("extern air.fast_sinpi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sinpi(x::Float32) = ccall("extern air.sinpi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sinpi(x::Float16) = ccall("extern air.sinpi.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.sqrt_fast(x::Float32) = ccall("extern air.fast_sqrt.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sqrt(x::Float32) = ccall("extern air.sqrt.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.sqrt(x::Float16) = ccall("extern air.sqrt.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.tan_fast(x::Float32) = ccall("extern air.fast_tan.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.tan(x::Float32) = ccall("extern air.tan.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.tan(x::Float16) = ccall("extern air.tan.f16", llvmcall, Float16, (Float16,), x)

@device_override FastMath.tanh_fast(x::Float32) = ccall("extern air.fast_tanh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.tanh(x::Float32) = ccall("extern air.tanh.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.tanh(x::Float16) = ccall("extern air.tanh.f16", llvmcall, Float16, (Float16,), x)

@device_function tanpi_fast(x::Float32) = ccall("extern air.fast_tanpi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function tanpi(x::Float32) = ccall("extern air.tanpi.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_function tanpi(x::Float16) = ccall("extern air.tanpi.f16", llvmcall, Float16, (Float16,), x)

@device_function trunc_fast(x::Float32) = ccall("extern air.fast_trunc.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.trunc(x::Float32) = ccall("extern air.trunc.f32", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.trunc(x::Float16) = ccall("extern air.trunc.f16", llvmcall, Float16, (Float16,), x)

### Integer Intrinsics

# TODO: absdiff, addsat, clamp, extract_bits, hadd, insert_bits, mad24, madhi, madsat,
#       max3, median3, min, min3, mul24, mulhi, rhadd, rotate, subsat

@device_override Base.abs(x::Int64)   = ccall("extern air.abs.s.i64", llvmcall, Int64, (Int64,), x)
@device_override Base.abs(x::UInt64)  = ccall("extern air.abs.u.i64", llvmcall, UInt64, (UInt64,), x)
@device_override Base.abs(x::Int32)   = ccall("extern air.abs.s.i32", llvmcall, Int32, (Int32,), x)
@device_override Base.abs(x::UInt32)  = ccall("extern air.abs.u.i32", llvmcall, UInt32, (UInt32,), x)
@device_override Base.abs(x::Int16)   = ccall("extern air.abs.s.i16", llvmcall, Int16, (Int16,), x)
@device_override Base.abs(x::UInt16)  = ccall("extern air.abs.u.i16", llvmcall, UInt16, (UInt16,), x)
@device_override Base.abs(x::Int8)    = ccall("extern air.abs.s.i8", llvmcall, Int8, (Int8,), x)
@device_override Base.abs(x::UInt8)   = ccall("extern air.abs.u.i8", llvmcall, UInt8, (UInt8,), x)

@device_function clz(x::Int64)   = ccall("extern air.clz.i64", llvmcall, Int64, (Int64,), x)
@device_function clz(x::UInt64)  = ccall("extern air.clz.i64", llvmcall, UInt64, (UInt64,), x)
@device_function clz(x::Int32)   = ccall("extern air.clz.i32", llvmcall, Int32, (Int32,), x)
@device_function clz(x::UInt32)  = ccall("extern air.clz.i32", llvmcall, UInt32, (UInt32,), x)
@device_function clz(x::Int16)   = ccall("extern air.clz.i16", llvmcall, Int16, (Int16,), x)
@device_function clz(x::UInt16)  = ccall("extern air.clz.i16", llvmcall, UInt16, (UInt16,), x)
@device_function clz(x::Int8)    = ccall("extern air.clz.i8", llvmcall, Int8, (Int8,), x)
@device_function clz(x::UInt8)   = ccall("extern air.clz.i8", llvmcall, UInt8, (UInt8,), x)

@device_function ctz(x::Int64)   = ccall("extern air.ctz.i64", llvmcall, Int64, (Int64,), x)
@device_function ctz(x::UInt64)  = ccall("extern air.ctz.i64", llvmcall, UInt64, (UInt64,), x)
@device_function ctz(x::Int32)   = ccall("extern air.ctz.i32", llvmcall, Int32, (Int32,), x)
@device_function ctz(x::UInt32)  = ccall("extern air.ctz.i32", llvmcall, UInt32, (UInt32,), x)
@device_function ctz(x::Int16)   = ccall("extern air.ctz.i16", llvmcall, Int16, (Int16,), x)
@device_function ctz(x::UInt16)  = ccall("extern air.ctz.i16", llvmcall, UInt16, (UInt16,), x)
@device_function ctz(x::Int8)    = ccall("extern air.ctz.i8", llvmcall, Int8, (Int8,), x)
@device_function ctz(x::UInt8)   = ccall("extern air.ctz.i8", llvmcall, UInt8, (UInt8,), x)

@device_function popcount(x::Int64)   = ccall("extern air.popcount.i64", llvmcall, Int64, (Int64,), x)
@device_function popcount(x::UInt64)  = ccall("extern air.popcount.i64", llvmcall, UInt64, (UInt64,), x)
@device_function popcount(x::Int32)   = ccall("extern air.popcount.i32", llvmcall, Int32, (Int32,), x)
@device_function popcount(x::UInt32)  = ccall("extern air.popcount.i32", llvmcall, UInt32, (UInt32,), x)
@device_function popcount(x::Int16)   = ccall("extern air.popcount.i16", llvmcall, Int16, (Int16,), x)
@device_function popcount(x::UInt16)  = ccall("extern air.popcount.i16", llvmcall, UInt16, (UInt16,), x)
@device_function popcount(x::Int8)    = ccall("extern air.popcount.i8", llvmcall, Int8, (Int8,), x)
@device_function popcount(x::UInt8)   = ccall("extern air.popcount.i8", llvmcall, UInt8, (UInt8,), x)

@device_function reverse_bits(x::Int64)   = ccall("extern air.reverse_bits.i64", llvmcall, Int64, (Int64,), x)
@device_function reverse_bits(x::UInt64)  = ccall("extern air.reverse_bits.i64", llvmcall, UInt64, (UInt64,), x)
@device_function reverse_bits(x::Int32)   = ccall("extern air.reverse_bits.i32", llvmcall, Int32, (Int32,), x)
@device_function reverse_bits(x::UInt32)  = ccall("extern air.reverse_bits.i32", llvmcall, UInt32, (UInt32,), x)
@device_function reverse_bits(x::Int16)   = ccall("extern air.reverse_bits.i16", llvmcall, Int16, (Int16,), x)
@device_function reverse_bits(x::UInt16)  = ccall("extern air.reverse_bits.i16", llvmcall, UInt16, (UInt16,), x)
@device_function reverse_bits(x::Int8)    = ccall("extern air.reverse_bits.i8", llvmcall, Int8, (Int8,), x)
@device_function reverse_bits(x::UInt8)   = ccall("extern air.reverse_bits.i8", llvmcall, UInt8, (UInt8,), x)


text:
function max_divisor(number::Integer) ::Integer
    return Integer(floor(sqrt(number)))
end

function is_prime(number::Integer) ::Bool
    max::Integer = max_divisor(number)
    # False if any number evenly divides it.
    for x=2:max
        if number % x == 0
            return false
        end
    end
    # true otherwise.
    return true
end

function next_prime(number::Integer) ::Integer
    to_check = number
    while true
        to_check += 1
        if is_prime(to_check)
            return to_check
        end
    end
end

function first_n_primes(number::Integer) ::Vector{<:Integer}
    last::Integer = 0
    primes = Vector{Integer}()
    for n=1:number
        last = next_prime(last)
        push!(primes, last)
    end
    return primes
end

function first_n_twin_primes(number::Integer) ::Vector{<:Integer}
    last::Integer = 0
    primes = Vector{Integer}()
    for n=3:number
        while true
            last = next_prime(last)
            if is_prime(last + 2)
                push!(primes, last)
                push!(primes, last+2)
                last += 2
                break
            end
        end
    end
    return primes
end

sum(x^2 for x in first_n_twin_primes(10000))


text:
function get_random_weights(weights::Params)::Params
    params([Float32.(randn(size(w))) for w in weights])
end

function normalize_direction!(direction::Array{Float32}, weight::Array{Float32})
    if ndims(direction) == 1
        direction[:] .= weight
    else
        @inbounds @views for i in 1:size(direction, 1)
            direction[i, :] .*= (
                norm(weight[i, :]) / (norm(direction[i, :]) + 1e-6)
            )
        end
    end
end

function normalize_directions!(directions::Params, weights::Params)
    for (direction, weight) in zip(directions, weights)
        normalize_direction!(direction, weight)
    end
end

function create_directions(weights::Params)::Params
    directions = get_random_weights(weights)
    normalize_directions!(directions, weights)
    directions
end

function create_directions(weights1::Params, weights2::Params)::Params
    params([w2 - w1 for (w1, w2) in zip(weights1, weights2)])
end

function shift_weights(
    weights::Params, directions::Tuple{Params, Params},
    step::Tuple{Float32, Float32},
)::Params
    direction_x, direction_y = directions
    step_x, step_y = step
    params([
        @. weight + dx * step_x + dy * step_y
        for (weight, dx, dy) in zip(weights, direction_x, direction_y)
   ])
end


text:
"""
"""
struct CellDof{DS} <: CellDatum
  cell_dof::AbstractArray{<:Union{Dof,AbstractArray{<:Dof}}}
  trian::Triangulation
  domain_style::DS
end

get_data(f::CellDof) = f.cell_dof
get_triangulation(f::CellDof) = f.trian
DomainStyle(::Type{CellDof{DS}}) where DS = DS()

function change_domain(a::CellDof,::ReferenceDomain,::PhysicalDomain)
  @notimplemented
end

function change_domain(a::CellDof,::PhysicalDomain,::ReferenceDomain)
  @notimplemented
end

function change_domain(a::CellDof,ttrian::Triangulation,target_domain::DomainStyle)
  change_domain(a,DomainStyle(a),ttrian,target_domain)
end

function change_domain(a::CellDof,::ReferenceDomain,ttrian::Triangulation,::ReferenceDomain)
  msg = """\n
  We cannot move the given CellDof to the reference domain of the requested triangulation.
  Make sure that the given triangulation is either the same as the triangulation on which the
  CellDof is defined, or that the latter triangulation is the background of the former.
  """
  strian = get_triangulation(a)
  if strian === ttrian
    return a
  end
  @assert is_change_possible(strian,ttrian) msg
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  change_domain_ref_ref(a,ttrian,sglue,tglue)
end

function change_domain(a::CellDof,::PhysicalDomain,ttrian::Triangulation,::PhysicalDomain)
  msg = """\n
  We cannot move the given CellDof to the physical domain of the requested triangulation.
  Make sure that the given triangulation is either the same as the triangulation on which the
  CellDof is defined, or that the latter triangulation is the background of the former.
  """
  strian = get_triangulation(a)
  if strian === ttrian
    return a
  end
  @assert is_change_possible(strian,ttrian) msg
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  change_domain_phys_phys(a,ttrian,sglue,tglue)
end

function change_domain(a::CellDof,::PhysicalDomain,trian::Triangulation,::ReferenceDomain)
  a_trian = change_domain(a,trian,PhysicalDomain())
  change_domain(a_trian,ReferenceDomain())
end

function change_domain(a::CellDof,::ReferenceDomain,trian::Triangulation,::PhysicalDomain)
  a_phys = change_domain(a,PhysicalDomain())
  change_domain(a_phys,trian,PhysicalDomain())
end

function change_domain_ref_ref(
  a::CellDof,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::FaceToFaceGlue)
  sface_to_dof = get_data(a)
  mface_to_sface = sglue.mface_to_tface
  tface_to_mface = tglue.tface_to_mface
  tface_to_mface_map = tglue.tface_to_mface_map
  mface_to_dof = extend(sface_to_dof,mface_to_sface)
  tface_to_dof = lazy_map(Reindex(mface_to_dof),tface_to_mface)
  @notimplementedif ! isa(tface_to_mface_map,AbstractArray{<:GenericField{typeof(identity)}})
  CellDof(tface_to_dof,ttrian,ReferenceDomain())
end

function change_domain_phys_phys(
  a::CellDof,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::FaceToFaceGlue)
  sface_to_dof = get_data(a)
  mface_to_sface = sglue.mface_to_tface
  tface_to_mface = tglue.tface_to_mface
  mface_to_dof = extend(sface_to_dof,mface_to_sface)
  tface_to_dof = lazy_map(Reindex(mface_to_dof),tface_to_mface)
  CellDof(tface_to_dof,ttrian,PhysicalDomain())
end

function get_cell_points(dofs::CellDof)
  cell_to_x = lazy_map(get_nodes,get_data(dofs))
  CellPoint(cell_to_x,get_triangulation(dofs),DomainStyle(dofs))
end

# Evaluation of CellDof

(a::CellDof)(f) = evaluate(a,f)

function evaluate!(cache,s::CellDof,f::CellField)

  trian_f = get_triangulation(f)
  trian_s = get_triangulation(s)

  if trian_f !== trian_s
    @unreachable """\n
    A CellDof can only be evaluated on a CellField defined on the same Triangulation.
    """
  end

  b = change_domain(f,s.domain_style)
  lazy_map(evaluate,get_data(s),get_data(b))
end

function evaluate!(cache, ::CellField, ::CellDof)
  @unreachable """\n
  CellField (f) objects cannot be evaluated at CellDof (s) objects.
  However, CellDofs objects can be evaluated at CellField objects.
  Did you mean evaluate(f,s) instead of evaluate(s,f), i.e.
  f(s) instead of s(f)?
  """
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Replace <: Edit
        col::Symbol
        from::Any
        to::Any
    end

Replace values in `col` `from` -> `to`.

# Arguments
- `col::Symbol`: name of column containing values to be replaced
- `from::Any`: value to replace
- `to::Any`: new value
"""
mutable struct Replace <: Edit
    "name of column containing values to be replaced"
    col::Symbol
    "value to replace"
    from::Any
    "new value"
    to::Any
end

function Replace(; col, from, to, )
    Replace(col, from, to, )
end

"""Get [`Replace`](@ref) `col`."""
get_col(value::Replace) = value.col
"""Get [`Replace`](@ref) `from`."""
get_from(value::Replace) = value.from
"""Get [`Replace`](@ref) `to`."""
get_to(value::Replace) = value.to

"""Set [`Replace`](@ref) `col`."""
set_col!(value::Replace, val) = value.col = val
"""Set [`Replace`](@ref) `from`."""
set_from!(value::Replace, val) = value.from = val
"""Set [`Replace`](@ref) `to`."""
set_to!(value::Replace, val) = value.to = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
export TimeBasedSampleModel

import StatsBase: sample

"""
    TimeBasedSampleModel(nactions::Int, κ::Float64 = 1e-4)
"""
mutable struct TimeBasedSampleModel <: AbstractEnvironmentModel
    experiences::Dict{
        Any,
        Dict{Any,NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}}},
    }
    nactions::Int
    κ::Float64
    t::Int
    last_visit::Dict{Tuple{Any,Any},Int}
    TimeBasedSampleModel(nactions::Int, κ::Float64 = 1e-4) =
        new(
            Dict{
                Any,
                Dict{
                    Any,
                    NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}},
                },
            }(),
            nactions,
            κ,
            0,
            Dict{Tuple{Any,Any},Int}(),
        )
end

function extract_experience(t::AbstractTrajectory, m::TimeBasedSampleModel)
    if length(t) > 0
        get_trace(t, :state)[end],
        get_trace(t, :action)[end],
        get_trace(t, :reward)[end],
        get_trace(t, :terminal)[end],
        get_trace(t, :next_state)[end]
    else
        nothing
    end
end

function RLBase.update!(m::TimeBasedSampleModel, transition::Tuple)
    s, a, r, d, s′ = transition
    if haskey(m.experiences, s)
        m.experiences[s][a] = (reward = r, terminal = d, nextstate = s′)
    else
        m.experiences[s] = Dict{
            Any,
            NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}},
        }(a => (reward = r, terminal = d, nextstate = s′))
    end
    m.t += 1
    m.last_visit[(s, a)] = m.t
end

function sample(m::TimeBasedSampleModel)
    s = rand(keys(m.experiences))
    a = rand(1:m.nactions)
    r, d, s′ = get(m.experiences[s], a, (0.0, false, s))
    r += m.κ * sqrt(m.t - get(m.last_visit, (s, a), 0))
    s, a, r, d, s′
end


text:
##################################################
#  More tests to improve coverage
##################################################

# using Base.Test
# reload("ReverseDiffSource")
# m = ReverseDiffSource

#########  rdiff ###########

m.@deriv_rule .+(x::Real   , y::AbstractArray)    x     (a=0.; for x in ds ; a+=x ; end; a)

ex = m.rdiff( :( sum( x .+ [1., 2.]) ), x=Float64)
ex = m.rdiff( :( sum( x .+ [1., 2.]) ), x=Float64, allorders=false)

######### plot  ############

ex = quote
	a = 0.
	for i in 1:10
		a += i
	end
	a
end

m.plot( m.tograph(ex) )

########## show() for nodes and graphs  ##########

g = m.tograph(ex)
println(g.nodes[4])
println(g)

########## 'ignore' argument  ##########

ex = :(p^3+y)

@test m.rdiff( ex , p=Float64, y=Float64, ignore=:y) ==
        striplinenumbers(:(begin ; (p ^ 3 + y,3 * p ^ 2) ; end))
@test m.rdiff( ex , p=Float64, y=Float64, ignore=:p) ==
        striplinenumbers(:(begin ; (p ^ 3 + y,1.0) ; end))
@test m.rdiff( ex , p=Float64, y=Float64, ignore=[:p;]) ==
        striplinenumbers(:(begin ; (p ^ 3 + y,1.0) ; end))

@test m.rdiff( ex , p=Float64, y=Float64, ignore=:p, order=3) ==
        striplinenumbers(:(begin ; (p ^ 3 + y,1.0, 0.0, 0.0) ; end))

@test_throws ErrorException m.rdiff( ex , p=Float64, y=Float64, order=3)
@test_throws ErrorException m.rdiff( ex , p=Float64, y=Float64, ignore=(:y,:p), order=3)

ex = :(dot(p,p) ^ sum(y))

m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:p;])
m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:y;])
m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:y;])
m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:p;])

m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:y;], order=2)
m.rdiff( ex , p=Vector{Float64}, y=Vector{Float64}, ignore=[:p;], order=2)

####### error conditions  #############

@test_throws UndefVarError m.rdiff( :( x * abcd ), x=1.)    # undefined external
@test_throws ErrorException m.rdiff( :( log(x) ), x=-1.)    # unevaluable function
@test_throws ErrorException m.rdiff( :( [1] > [2] ), x=-1.) # unevaluable comparison

@test_throws ErrorException m.tograph(:(log(a) = 1,2))   # incorrect LHS

###### allorders rdiff flags  ################

m.rdiff( :( log(x) ), x=Float64, allorders=false)

###### BitArray type  ############

m.rdiff( :( sum(x .* falses(2)) ), x=Float64  )

####### Issue #25 (splinary not trigered) #############
# (happens when function are prefixed with a module)

module B; end

x = 2.
ex = quote
  B.max(1., x, 3.)
end

m.rdiff(ex, x=Float64)

function foo(x,y,z)
    return x + y + z
end
f = m.rdiff(foo,(Float64,Float64,Float64))

###### Issue #32  (allorders removing some legitimate variables)  ######

ex = m.rdiff( :(y*x^3+y^5) , x=Float64, y=Float64, order=1)
@test length(eval(:(x=2.;y=1.;$ex))) == 3

ex = m.rdiff( :(y*x^3+y^5) , x=Float64, y=Float64, order=1, allorders=false)
@test length(eval(:(x=2.;y=1.;$ex))) == 2

###### Issue #48  (order 0 regression)  ######

@test m.rdiff(:(dot(z, z)), z=Vector{Float64}, order=0) ==
		striplinenumbers(:(begin ; dot(z,z) ; end))


text:
# Define and register smooth functions
# These are "smooth" aka differentiable and avoid Gibbs effect
# These follow: `offset` + `smooth_wave` * `smooth_step` with zero output for `t < start_time`
function smooth_cos(x, δ, f, amplitude, ϕ, offset, start_time)
    offset + amplitude * cos(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
end

function smooth_damped_sin(x, δ, f, amplitude, damping, ϕ, offset, start_time)
    offset + exp((start_time - x)*damping)*amplitude*sin(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
end

function smooth_ramp(x, δ, height, duration, offset, start_time)
    offset + height/(duration) * (smooth_xH(x, δ, start_time) - smooth_xH(x, δ, start_time+duration))
end

function smooth_sin(x, δ, f, amplitude, ϕ, offset, start_time)
    offset + amplitude * sin(2*pi*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
end

function smooth_square(x, δ, f, amplitude, offset, start_time)
    offset + amplitude*2atan(sin(2π*(x - start_time)*f)/δ)/π * smooth_step(x, δ, one(x), zero(x), start_time)
end

function smooth_step(x, δ, height, offset, start_time)
    offset + height*(atan((x - start_time)/δ)/π + 0.5)
end

function smooth_triangular(x, δ, f, amplitude, offset, start_time)
    offset + amplitude * (1-2acos((1 - δ)sin(2π*(x - start_time)*f))/π) * smooth_step(x, δ, one(x), zero(x), start_time)
end

function smooth_xH(x, δ, tₒ)
    0.5*(x-tₒ) * (1+((x-tₒ)/sqrt((x-tₒ)^2+δ^2)))
end

function square(x, f, amplitude, offset, start_time)
    offset + (x > start_time) * (amplitude * (4*floor(f*(x - start_time)) - 2*floor(2*(x - start_time)*f) + 1))
end

function triangular(x, f, amplitude, offset, start_time)
    p = 1/f # period
    offset + (x > start_time) * (4 * amplitude * f * abs(abs((x - p/4 - start_time) % p) - p/2) - amplitude)
end

"""
Generate constant signal.

# Parameters:
- `k`: Constant output value

# Connectors:
- `output`
"""
function Constant(;name, k=1)
    @named output = RealOutput()
    pars = @parameters k=k
    eqs = [
        output.u ~ k
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate sine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function Sine(;name,
    frequency,
    amplitude=1,
    phase=0,
    offset=0,
    start_time=0,
    smooth=false)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
    equation = if smooth == false
        offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
    else
        δ = 1e-5
        smooth_sin(t, δ, frequency, amplitude, phase, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate cosine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""

function Cosine(;name,
    frequency,
    amplitude=1,
    phase=0,
    offset=0,
    start_time=0,
    smooth=false)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
    equation = if smooth == false
        offset + ifelse(t < start_time, zero(t), amplitude* cos(2*pi*frequency*(t - start_time) + phase))
    else
        δ = 1e-5
        smooth_cos(t, δ, frequency, amplitude, phase, offset, start_time)
    end
    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate current time signal.

# Parameters:
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`

# Connectors:
- `output`
"""
function ContinuousClock(;name, offset=0, start_time=0)
    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time
    eqs = [
        output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time)
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate ramp signal.

# Parameters:
- `height`: Height of ramp
- `duration`: [s] Duration of ramp (= 0.0 gives a Step)
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function Ramp(;name,
    height=1,
    duration=1,
    offset=0,
    start_time=0,
    smooth=false)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time height=height duration=duration
    equation = if smooth == false
        offset + ifelse(t < start_time, 0,
            ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))
    else
        δ = 1e-5
        smooth_ramp(t, δ, height, duration, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate smooth square signal.

# Parameters:
- `frequency`: [Hz] Frequency of square wave
- `amplitude`: Amplitude of square wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function Square(; name, frequency=1.0, amplitude=1.0,
    offset=0.0, start_time=0.0, smooth=false)
    δ = 1e-5

    @named output = RealOutput()
    pars = @parameters begin
        frequency=frequency
        amplitude=amplitude
        offset=offset
        start_time=start_time
    end

    equation = if smooth == false
        square(t, frequency, amplitude, offset, start_time)
    else
        δ = 1e-5
        smooth_square(t, δ, frequency, amplitude, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate step signal.

# Parameters:
- `height`: Height of step
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function Step(;name, height=1, offset=0, start_time=0, smooth=true)
    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time height=height
    equation = if smooth == false
        offset + ifelse(t < start_time, zero(t), height)
    else
        δ = 1e-5
        smooth_step(t, δ, height, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate exponentially damped sine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `damping`: [1/s] Damping coefficient of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function ExpSine(; name,
    frequency,
    amplitude=1,
    damping=0.1,
    phase=0,
    offset=0,
    start_time=0,
    smooth=false)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping

    equation = if smooth == false
        offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
    else
        δ = 1e-5
        smooth_damped_sin(t, δ, frequency, amplitude, damping, phase, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate smooth triangular signal for frequencies less than or equal to 25 Hz

# Parameters:
- `frequency`: [Hz] Frequency of square wave
- `amplitude`: Amplitude of square wave
- `offset`: Offset of output signal.
- `start_time`: [s] Output `y = offset` for `t < start_time`
- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
             It uses a smoothing factor of `δ=1e-5`

# Connectors:
- `output`
"""
function Triangular(; name, amplitude=1.0, frequency=1.0,
    offset=0.0, start_time=0.0, smooth=false)

    @named output = RealOutput()
    pars = @parameters begin
        amplitude=amplitude
        frequency=frequency
        offset=offset
        start_time=start_time
    end

    equation = if smooth == false
        triangular(t, frequency, amplitude, offset, start_time)
    else
        δ = 1e-5
        smooth_triangular(t, δ, frequency, amplitude, offset, start_time)
    end

    eqs = [
        output.u ~ equation
    ]

    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

# TODO:
# - Exponentials    Generate a rising and falling exponential signal
# - Pulse   Generate pulse signal of type Real
# - SawTooth    Generate saw tooth signal
# - Trapezoid   Generate trapezoidal signal of type Real


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 265ed188-5bf6-11eb-0d29-51ea74e6c7d1
using Pkg, DrWatson

# ╔═╡ ac6f2dfc-5bf6-11eb-3dc0-b5368b466ff7
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ a4042486-5bf5-11eb-0183-33fd00d868e4
md" ## Clip-07-16-18.1s.jl"

# ╔═╡ ac8227b2-5bf6-11eb-27dd-09924b99e6c2
stan7_9 = "
data {
    int<lower=1> K;
    int<lower=0> N;
    matrix[N, K] x;
    vector[N] y;
}
parameters {
    real a;
    vector[K] b;
    real<lower=0> sigma;
}
transformed parameters {
    vector[N] mu;
    mu = a + x * b;
}
model {
	a ~ normal(0, 100);
	b ~ normal(0, 10);
	sigma ~ exponential(1);
    y ~ normal(mu, sigma);          // observed model
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
";

# ╔═╡ 00ee9dc0-5c23-11eb-10c3-1dfb4f486e6f
begin
	tmpdir = joinpath(projectdir(), "tmp")
	m7_9s = SampleModel("m7.9s", stan7_9; tmpdir)
end;

# ╔═╡ 5d5b8982-64b2-11eb-0ca9-b94e3197ff31
begin
	Ns = [20, 100] 				# Number of observations
	rho = [0.15, -0.4]			# Covariance between x1 and x2
	L = 50 					    # Number of simulations
	K = 5 						# Number of slopes
	dev_is = zeros(L, K)
	dev_os = zeros(L, K)
	res = Vector{NamedTuple}(undef, length(Ns))
end;

# ╔═╡ 51b7be24-6fad-11eb-3c67-d53dcaf1ada6
md"
!!! note
	This might take some time ...
"

# ╔═╡ f253787c-64c3-11eb-2cd4-9322707100b7
begin
	for (ind, N) in enumerate(Ns)
		for i in 1:L
			for j = 1:K
				println("N = $(Ns[ind]), run = $i, no of b parms = $j")
				y, x_train, x_test = sim_train_test(;N, K, rho)
				data = (N = size(x_train, 1), K = size(x_train, 2),
					y = y, x = x_train,
					N_new = size(x_test, 1), x_new = x_test)
				rc7_9s = stan_sample(m7_9s; data)
				if success(rc7_9s)

					# use `logprob()`

					post7_9s_df = read_samples(m7_9s, :dataframe)
					lp_train = logprob(post7_9s_df, x_train, y, j)
					dev_is[i, j] = -2sum(lppd(lp_train))
					lp_test = logprob(post7_9s_df, x_test, y, j)
					dev_os[i, j] = -2sum(lppd(lp_test))

				end
			end
		end
		res[ind] = (
			mean_dev_is = mean(dev_is, dims=1),
			std_dev_is = std(dev_is, dims=1),
			mean_dev_os = mean(dev_os, dims=1),
			std_dev_os =std(dev_os, dims=1)
		)
	end
end;

# ╔═╡ 7eb8b028-657f-11eb-0abf-37c3819aea78
begin
	fig = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
	for i in 1:2
		xcoord = collect(1:K) .- 0.05
		ycoord = res[i].mean_dev_is[1,:]
		ylims = i == 1 ? (0,100) : (200, 350)
		fig[i] = scatter(xcoord, ycoord, xlab = "No of b parameters",
			ylab = "Deviance", ylims = ylims, leg=:bottomleft,
			lab = "In sample", markersize=2)
		for j in 1:K
			plot!([xcoord[j], xcoord[j]],
				[ycoord[j]-res[i].std_dev_is[j], ycoord[j]+res[i].std_dev_is[j]],
				lab=false, color=:darkblue)
		end
		title!("N = $(Ns[i]), L = $(L)")
		xcoord = collect(1:K) .+ 0.05
		ycoord = res[i].mean_dev_os[1,:]
		scatter!(xcoord, ycoord, lab = "Out of sample", markersize=2)
		for j in 1:K
			plot!([xcoord[j], xcoord[j]],
				[ycoord[j]-res[i].std_dev_os[j], ycoord[j]+res[i].std_dev_os[j]],
				lab=false, color=:red)
		end
		title!("N = $(Ns[i]), L = $(L)")
	end
	plot(fig..., layout=(1,2))
end

# ╔═╡ 918392e8-6888-11eb-344d-8d15802cfcac
res

# ╔═╡ 3f775c58-6fa1-11eb-3c36-dd192855ec6b
md"
!!! note
	Compare waic and psis estimates to in-sample estimate
"

# ╔═╡ 73872126-6f9d-11eb-027e-e3e81ada62b7
begin
	devs = Vector{Float64}(undef, 2)
	y, x_train, x_test = sim_train_test(;N=20, K=5, rho)
	data = (N = size(x_train, 1), K = size(x_train, 2),
		y = y, x = x_train,
		N_new = size(x_test, 1), x_new = x_test)
	rc7_9s = stan_sample(m7_9s; data)
	if success(rc7_9s)

		# use `logprob()`

		post7_9s_df = read_samples(m7_9s, :dataframe)
		lp_train = logprob(post7_9s_df, x_train, y, 2)
		devs[1] = -2sum(lppd(lp_train))
		lp_test = logprob(post7_9s_df, x_test, y, 2)
		devs[2] = -2sum(lppd(lp_test))
	end
end;

# ╔═╡ 28ce4ce2-6f9e-11eb-037c-85da29a882e3
devs

# ╔═╡ 5e6d444e-6f9f-11eb-2832-558bbc4303dc
if success(rc7_9s)
	waic(m7_9s)
end

# ╔═╡ 66ceb86a-6fa2-11eb-224d-9f617f07019d
if success(rc7_9s)
	loo7_9s, loos7_9s, pk7_9s = psisloo(m7_9s)
	-2loo7_9s
end

# ╔═╡ 72523564-6f9f-11eb-1481-1f7df6e36dfb
# Use `generated quantities`

if success(rc7_9s)
	nt7_9s = read_samples(m7_9s)
	-2sum(lppd(nt7_9s.log_lik'))
end

# ╔═╡ 74d4016e-6563-11eb-0c3b-b7c4f81baca6
md" ## End of clip-07-16-18.1s.jl"

# ╔═╡ Cell order:
# ╟─a4042486-5bf5-11eb-0183-33fd00d868e4
# ╠═265ed188-5bf6-11eb-0d29-51ea74e6c7d1
# ╠═ac6f2dfc-5bf6-11eb-3dc0-b5368b466ff7
# ╠═ac8227b2-5bf6-11eb-27dd-09924b99e6c2
# ╠═00ee9dc0-5c23-11eb-10c3-1dfb4f486e6f
# ╠═5d5b8982-64b2-11eb-0ca9-b94e3197ff31
# ╟─51b7be24-6fad-11eb-3c67-d53dcaf1ada6
# ╠═f253787c-64c3-11eb-2cd4-9322707100b7
# ╠═7eb8b028-657f-11eb-0abf-37c3819aea78
# ╠═918392e8-6888-11eb-344d-8d15802cfcac
# ╟─3f775c58-6fa1-11eb-3c36-dd192855ec6b
# ╠═73872126-6f9d-11eb-027e-e3e81ada62b7
# ╠═28ce4ce2-6f9e-11eb-037c-85da29a882e3
# ╠═5e6d444e-6f9f-11eb-2832-558bbc4303dc
# ╠═66ceb86a-6fa2-11eb-224d-9f617f07019d
# ╠═72523564-6f9f-11eb-1481-1f7df6e36dfb
# ╟─74d4016e-6563-11eb-0c3b-b7c4f81baca6


text:
using Terming

function main()
    # set term size
    Terming.displaysize(20, 75)
    # switch to alternate screen
    Terming.alt_screen(true)

    Terming.cmove(1, 1)
    Terming.println("Terminal is now switched to the alternate screen mode.")
    Terming.println("Press ENTER to switch back."); readline()

    # switch back from alternate screen
    Terming.alt_screen(false)

    return
end

main()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "gdk_pixbuf"
version = v"2.42.6"

# Collection of sources required to build gdk-pixbuf
sources = [
    ArchiveSource("https://gitlab.gnome.org/GNOME/gdk-pixbuf/-/archive/$(version)/gdk-pixbuf-$(version).tar.bz2",
                  "8a76cffe6a85f2602cf246c1c974eb475aea41c363a932f0c34695fa968f01fd"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gdk-pixbuf-*/
mkdir build && cd build

FLAGS=()
if [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(-Dx11=false)
fi

meson .. \
    -Dgir=false \
    -Dman=false \
    -Dinstalled_tests=false \
    -Dgio_sniffing=false \
    "${FLAGS[@]}" \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install

# Cleanup `loaders.cache` file, we're going to generate a new one on the user's machine
find ${prefix}/lib -name loaders.cache -delete
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgdk_pixbuf-2", "libgdk_pixbuf-2.0"], :libgdkpixbuf),
    ExecutableProduct("gdk-pixbuf-query-loaders", :gdk_pixbuf_query_loaders),
    FileProduct("lib/gdk-pixbuf-2.0/2.10.0/loaders", :gdk_pixbuf_loaders_dir),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Need a host gettext for msgfmt
    HostBuildDependency("Gettext_jll"),
    # Need a host glib for glib-compile-resources
    HostBuildDependency("Glib_jll"),
    Dependency("Glib_jll"; compat="2.68.1"),
    Dependency("JpegTurbo_jll"),
    Dependency("libpng_jll"),
    Dependency("Libtiff_jll"; compat="4.3.0"),
    Dependency("Xorg_libX11_jll"),
    BuildDependency("Xorg_xproto_jll"),
    BuildDependency("Xorg_kbproto_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
solver = POMCPSolver(estimate_value=RolloutEstimator(FeedWhenCrying()),
                    eps=0.01,
                    c=10.0,
                    tree_queries=50,
                    rng=MersenneTwister(2))

problem = BabyPOMDP()
policy = solve(solver, problem)
up = updater(policy)
belief = initialize_belief(up, initial_state_distribution(problem))
action(policy, belief)

dummy = IOBuffer()
show(dummy, MIME("text/html"), POMCPTreeVisualizer(belief))


text:
using Empirikos
using Test
using Hypatia
using Random
Zs_raw = NormalSample.(randn(10), rand(10))
Zs_compound = compound(Zs_raw)

Zs_standard = StandardNormalSample.(randn(100))
Zs_standard_compound = compound(Zs_standard)

@test Zs_standard_compound[1].vec == [StandardNormalSample()]

@test cdf.(Normal(0,2), Zs_standard) == cdf.(Normal(0,2), Zs_standard)
@test pdf.(Normal(0,2), Zs_standard) == pdf.(Normal(0,2), Zs_standard)

dkw_loc_standard_compound = fit(DvoretzkyKieferWolfowitz(α=0.05*exp(1)), Zs_standard_compound)
dkw_loc_standard = fit(DvoretzkyKieferWolfowitz(α=0.05), Zs_standard)

@test dkw_loc_standard.homoskedastic
@test !dkw_loc_standard_compound.homoskedastic

@test dkw_loc_standard.band ≈ dkw_loc_standard_compound.band

target = PosteriorMean(NormalSample(2.0, 2.0))

gcal = DiscretePriorClass(-2:0.2:2)
floc_method_standard = FLocalizationInterval(flocalization=DvoretzkyKieferWolfowitz(0.05),
                             convexclass=gcal,
                             solver=Hypatia.Optimizer)
floc_method_compound = FLocalizationInterval(flocalization=DvoretzkyKieferWolfowitz(0.05*exp(1)),
                             convexclass=gcal,
                             solver=Hypatia.Optimizer)

ci1 = confint(floc_method_standard, target, Zs_standard)
ci2 = confint(floc_method_compound, target, Zs_standard_compound)
@test ci1.lower ≈ ci2.lower
@test ci1.upper ≈ ci2.upper
@test ci1.α ≈ 0.05
@test ci2.α ≈ 0.05*exp(1)

# Same with actually heteroskedastic samples

Random.seed!(1)
σs = rand(1:10,1000)
μs = randn(1000)
Zs = NormalSample.(μs, σs)
Zs_compound = compound(Zs)
@test length(Zs_compound[1].vec) == 10

dkw_loc_not_compound = fit(DvoretzkyKieferWolfowitz(α=0.05), Zs)
dkw_loc_compound = fit(DvoretzkyKieferWolfowitz(α=0.05), Zs_compound)

@test !dkw_loc_not_compound.homoskedastic
@test !dkw_loc_compound.homoskedastic

@test dkw_loc_not_compound.band ≈ dkw_loc_compound.band
@test length(dkw_loc_not_compound.summary) == length(dkw_loc_compound.summary)

Random.seed!(1)
Zs_binom_heterosk = BinomialSample.( sample(1:4, 100), sample(20:22,100))
Zs_binom_summary = summarize(Zs_binom_heterosk)
@test nobs(Zs_binom_summary) == 100
@test length(keys(Zs_binom_summary)) == 3*4
comp_binom = compound(Zs_binom_heterosk)
Zs_binom_comp_summary = summarize(comp_binom)
@test nobs(Zs_binom_comp_summary) == 100
@test length(keys(Zs_binom_comp_summary)) == 4

Zs_binom_summary_comp = compound(Zs_binom_summary)
@test nobs(Zs_binom_summary_comp) == 100

Empirikos.multiplicity(Zs_binom_comp_summary) == Empirikos.multiplicity(Zs_binom_summary_comp)


text:
module Graphics

using Reexport
using Requires

include("WiringDiagramLayouts.jl")
include("Graphviz.jl")
include("GraphvizGraphs.jl")
include("GraphvizWiringDiagrams.jl")
include("ComposeWiringDiagrams.jl")
include("TikZ.jl")
include("TikZWiringDiagrams.jl")
include("YFilesWiringDiagrams.jl")

@reexport using .WiringDiagramLayouts
@reexport using .GraphvizGraphs
@reexport using .GraphvizWiringDiagrams
@reexport using .ComposeWiringDiagrams
@reexport using .TikZWiringDiagrams
@reexport using .YFilesWiringDiagrams

function __init__()
  @require Convex="f65535da-76fb-5f13-bab9-19810c17039a" begin
    @require SCS="c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" begin
      include("ConvexExternal.jl")
    end
  end
  @require TikzPictures="37f6aa50-8035-52d0-81c2-5a1d08754b2d" begin
    include("TikZExternal.jl")
  end
end

end


text:
# ---
# title: 970. Powerful Integers
# id: problem970
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/powerful-integers/description/>
# hidden: true
# ---
#
# Given two positive integers `x` and `y`, an integer is _powerful_  if it is
# equal to `x^i + y^j` for some integers `i >= 0` and `j >= 0`.
#
# Return a list of all _powerful_ integers that have value less than or equal to
# `bound`.
#
# You may return the answer in any order.  In your answer, each value should
# occur at most once.
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 2, y = 3, bound = 10
#     Output: [2,3,4,5,7,9,10]
#     Explanation:
#     2 = 2^0 + 3^0
#     3 = 2^1 + 3^0
#     4 = 2^0 + 3^1
#     5 = 2^1 + 3^1
#     7 = 2^2 + 3^1
#     9 = 2^3 + 3^0
#     10 = 2^0 + 3^2
#
#
# **Example 2:**
#
#
#
#     Input: x = 3, y = 5, bound = 15
#     Output: [2,4,6,8,10,14]
#
#
#
#
# **Note:**
#
#   * `1 <= x <= 100`
#   * `1 <= y <= 100`
#   * `0 <= bound <= 10^6`
#
#
## @lc code=start
using LeetCode

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


text:
using Documenter, StateSpaceReconstruction
push!(LOAD_PATH,"../src/")
ENV["GKSwstype"] = "100"

PAGES = [
    "Overview" => "index.md",
    "Embedding" => "embedding.md",
    "Rectangular partitions" => [
        "Types of partitions" => "partition/binningtypes.md",
        "Finding axis minima and step sizes" =>  "partition/minima_stepsizes.md",
        "Coordinate representation" => "partition/coordinate_representation.md",
        "Finding axis minima and step sizes" => "partition/minima_stepsizes.md",
        "Marginal visitation frequences" =>  "partition/marginal_visitation_frequency.md"
    ],
    "Simplex partitions" => [
        "Overview" => "simplexpartition/simplex.md"
    ]
]

makedocs(
    modules = [StateSpaceReconstruction],
    format = :markdown,
    sitename = "StateSpaceReconstruction.jl",
    authors = "Kristian Agasøster Haaga",
    pages = PAGES,
    clean = true
    # Use clean URLs, unless built as a "local" build
    #html_prettyurls = !("local" in ARGS),
    #html_canonical = "https://kahaaga.github.io/StateSpaceReconstruction.jl/latest/"
)

# deploydocs(;
#     Deps.pip("pygments", "mkdocs", "python-markdown-math"),
#     repo   = "github.com/kahaaga/StateSpaceReconstruction.jl.git",
#     julia  = "0.6",
#     target = "build",
#     deps = nothing,
#     make = nothing,
#     osname = "linux"
# )


text:
# An implementation of Nesterov's Accelerated method based on a reformulation
# that makes it similar to ordinary (Stochastic) Gradient Descent with Momentum
#
# Yoshua Bengio, Nicolas Boulanger-Lewandowski, Razvan Pascanu. Advances in
# Optimizing Recurrent Networks. arXiv:1212.0901 [cs.LG]

immutable Nesterov <: SolverMethod
end
type NesterovSolverState <: InternalSolverState
    learning_rate :: Float64
    momentum      :: Float64
    param_states  :: Vector{LayerState}
    param_history :: Vector{Vector{Blob}}
    last_momentum :: Float64
end

type NesterovSolverSnapshot <: SolverStateSnapshot
    iteration     :: Int
    obj_val       :: Float64
    learning_rate :: Float64
    momentum      :: Float64
end

#typealias NesterovSolverState SGDSolverState
#typealias NesterovSolverSnapashot SGDSolverSnapshot

function update(solver::Solver{Nesterov}, net::Net, state::SolverState{NesterovSolverState})
  # TODO check if this is identical with sgd.jl update()?
  for i = 1:length(state.internal.param_states)
    layer_state   = state.internal.param_states[i]
    history = state.internal.param_history[i]
    for j = 1:length(layer_state.parameters)
      hist_blob = history[j]
      gradient  = layer_state.parameters[j].gradient
      data_type = eltype(hist_blob)

      update_parameters!(net, solver.method, layer_state.parameters[j].learning_rate * state.internal.learning_rate,
          state.internal.last_momentum, state.internal.momentum, layer_state.parameters[j].blob, hist_blob, gradient, data_type)
    end
  end

  state.internal.last_momentum = state.internal.momentum
end

function shutdown(state::NesterovSolverState)
  map(x -> map(destroy, x), state.param_history)
end

function update_parameters!(net::Net{CPUBackend}, method::Nesterov, learning_rate,
    last_momentum, momentum, param_blob, hist_blob, gradient, data_type)

  # param_blob += -last_momentum* hist_blob (update with vt-1)
  BLAS.axpy!(length(hist_blob), convert(data_type, -last_momentum), hist_blob.data, 1, param_blob.data, 1)

  # hist_blob = last_momentum * hist_blob
  BLAS.scal!(length(hist_blob), convert(data_type, last_momentum), hist_blob.data, 1)

  # hist_blob = -learning_rate * gradient + hist_blob (calc vt)
  BLAS.axpy!(length(hist_blob), convert(data_type, -learning_rate), gradient.data, 1, hist_blob.data, 1)

  # param_blob += (1+momentum) * hist_blob (update with vt)
  BLAS.axpy!(length(hist_blob), convert(data_type, 1 + momentum), hist_blob.data, 1, param_blob.data, 1)
end


text:
include("load.jl")
EventTrackerTests.runtests()


text:
@testset "Defaults" begin
    @test Shuffle.DEFAULTS.shuffle === RandomShuffle()

    let mt1 = MersenneTwister(1234), mt2 = MersenneTwister(1234), x = rand(1000)
        @test shuffle(mt1, x) == shuffle(mt2, x, RandomShuffle())
    end

    let mt1 = MersenneTwister(123), mt2 = MersenneTwister(123), x = rand(1000)
        @test shuffle!(mt1, x) == shuffle!(mt2, x, RandomShuffle())
    end

    let mt1 = MersenneTwister(12), mt2 = MersenneTwister(12), x = rand(1000)
        @test nshuffle(mt1, x, 5) == nshuffle(mt2, x, 5, RandomShuffle())
    end

    let mt1 = MersenneTwister(1), mt2 = MersenneTwister(1), x = rand(1000)
        @test nshuffle!(mt1, x, 5) == nshuffle!(mt2, x, 5, RandomShuffle())
    end

    Shuffle.DEFAULTS.shuffle = Faro{:in}()
    @test Shuffle.DEFAULTS.shuffle === Faro{:in}()

    @test shuffle(collect(1:100)) == shuffle(collect(1:100), Faro{:in}())

    @test shuffle!(collect(1:100)) == shuffle!(collect(1:100), Faro{:in}())

    @test nshuffle(collect(1:100), 5) == nshuffle(collect(1:100), 5, Faro{:in}())

    @test nshuffle!(collect(1:100), 5) == nshuffle!(collect(1:100), 5, Faro{:in}())

    Shuffle.DEFAULTS.shuffle = RandomShuffle()
    @test Shuffle.DEFAULTS.shuffle === RandomShuffle()
end

@testset "Default RNG" begin
    rng = copy(default_rng())
    @test shuffle(collect(1:100), RandomShuffle()) ==
          shuffle(rng, collect(1:100), RandomShuffle())

    rng = copy(default_rng())
    @test shuffle!(collect(1:100), RandomShuffle()) ==
          shuffle!(rng, collect(1:100), RandomShuffle())

    rng = copy(default_rng())
    @test nshuffle(collect(1:100), 5, RandomShuffle()) ==
          nshuffle(rng, collect(1:100), 5, RandomShuffle())

    rng = copy(default_rng())
    @test nshuffle!(collect(1:100), 6, RandomShuffle()) ==
          nshuffle!(rng, collect(1:100), 6, RandomShuffle())
end


text:
@generated function harm(x::Integer...)
  l = length(x)
  if(l == 0)
    return 0
  end
  ex = :(1/x[1])
  for i = 2:l
    ex = :($ex + 1/x[$i])
  end
  ex = :($l/$ex)
  println(ex)
  return ex
end

function harm2(x::Integer...)
  l = length(x)
  if(l == 0)
    return 0
  end
  ex = :(1/x[1])
  for i = 2:l
    ex = :($ex + 1/x[$i])
  end
  ex = :($l/$ex)
  return ex
end


text:
nu = 0.01
buoyance_coef = 1.0

function u_exact(x,y)
    x*(1-x)*y*(1-y)
end

function v_exact(x,y)
    x*(1-x)*y*(1-y)
end

function p_exact(x,y)
    x*(1-x)*y*(1-y)
end

function t_exact(x,y)
    x*(1-x)*y*(1-y)
end

function ffunc_(x, y)
    x*y*(1 - x)*(1 - y)*(-x*y*(1 - x) + x*(1 - x)*(1 - y)) + x*y*(1 - x)*(1 - y)*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) - x*y*(1 - y) + 0.02*x*(1 - x) + y*(1 - x)*(1 - y) + 0.02*y*(1 - y)
end

function gfunc_(x, y)
    x*y*(1 - x)*(1 - y)*(-x*y*(1 - x) + x*(1 - x)*(1 - y)) + x*y*(1 - x)*(1 - y)*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) - buoyance_coef * x*y*(1 - x)*(1 - y) - x*y*(1 - x) + x*(1 - x)*(1 - y) + 0.02*x*(1 - x) + 0.02*y*(1 - y)
end

function hfunc_(x,y)
    -x*y*(1 - x) - x*y*(1 - y) + x*(1 - x)*(1 - y) + y*(1 - x)*(1 - y)
end

function heat_source_func(x, y)
    (2*x^2*y*(x - 1)*(2*y - 1) + x*y*(x - 1)*(y - 1)*(y^2 + 1)^2*(x*(x - 1)*(2*y - 1) + y*(2*x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
end

m = 20
n = 20
h = 1/n

mesh = Mesh(m, n, h, degree=2)
nnode = mesh.nnode
nedge = mesh.nedge
ndof = mesh.ndof
nelem = mesh.nelem
ngauss = get_ngauss(mesh)

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc_, mesh), mesh)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc_, mesh), mesh)
H = get_area(mesh) .* eval_f_on_fvm_pts(hfunc_, mesh)
B = constant(compute_interaction_matrix(mesh))

# compute F
Laplace = compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh)
heat_source = eval_f_on_gauss_pts(heat_source_func, mesh)
heat_source = constant(compute_fem_source_term1(heat_source, mesh))

kgauss = eval_f_on_gauss_pts(k_func, mesh); kgauss=stack(kgauss)
LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, mesh))

# xy = fem_nodes(mesh)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = dof_to_gauss_points(k, mesh)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, mesh)

function compute_residual(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]

    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    interaction = compute_interaction_term(p, mesh) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, mesh)
    f2 = compute_fem_source_term1(vgauss.*uy, mesh)
    f3 = -interaction[1:ndof]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, mesh)
    g2 = compute_fem_source_term1(vgauss.*vy, mesh)
    g3 = -interaction[ndof+1:end]
    g4 = Laplace*v
    g5 = -F2

    T_gauss = dof_to_gauss_points(T, mesh)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, mesh)

    G = g1 + g2 + g3 + g4 + g5 + buoyance_term

    H0 = -B * [u;v] + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, mesh) * T - heat_source
    R = [F;G;H0;T0]
    return R
    # return [F; G; constant(zeros(nelem)); constant(zeros(ndof)) ]
end

function compute_jacobian(S)
    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]

    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    M1 = constant(compute_fem_mass_matrix1(ux, mesh))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, mesh))

    N1 = constant(compute_fem_mass_matrix1(vy, mesh))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, mesh))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, mesh))

    gradT = eval_grad_on_gauss_pts1(T, mesh)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, mesh))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, mesh))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = constant(compute_fem_mass_matrix1(-buoyance_coef * constant(ones(ngauss)), mesh))
    T_mat = [SparseTensor(spzeros(ndof, ndof)); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(nelem, ndof))]

    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros(ndof, nelem)) M]]

    # return [J[1:2ndof, 1:2ndof] SparseTensor(spzeros(2ndof, nelem+ndof))
        # SparseTensor(spzeros(nelem+ndof, 2ndof))  spdiag(ones(nelem+ndof)) ]
    end

NT = 8    # number of iterations for Newton's method

bd = Array{Int64, 1}([])
eps = 1e-6
for j = 1:nnode
    nodex, nodey = mesh.nodes[j, 1], mesh.nodes[j, 2]
    if abs(nodex-0.0) <= eps || abs(nodex-1.0) <= eps || abs(nodey-0.0) <= eps || abs(nodey-1.0) <= eps
        global bd = [bd; j]
    end
end
for j = 1:nedge
    edgex, edgey = (mesh.nodes[mesh.edges[j, 1], :] .+ mesh.nodes[mesh.edges[j, 2], :]) ./ 2
    if abs(edgex-0.0) <= eps || abs(edgex-1.0) <= eps || abs(edgey-0.0) <= eps || abs(edgey-1.0) <= eps
        global bd = [bd; nnode + j]
    end
end

fvm_bd = Array{Int64, 1}([])
elemxy = fvm_nodes(mesh)
for j = 1:nelem
    elemx, elemy = elemxy[j, 1], elemxy[j, 2]
    if abs(elemx-0.0) <= 2*h && abs(elemy-0.0) <= h
        global fvm_bd = [fvm_bd; j]
    end
    if abs(elemx-1.0) <= 2*h && abs(elemy-0.0) <= h
        global fvm_bd = [fvm_bd; j]
    end
end

bd = [bd; bd .+ ndof; fvm_bd .+ 2*ndof; bd .+ (2*ndof+nelem)]
# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, residual = impose_Dirichlet_boundary_conditions(J, residual, bd, zeros(length(bd)))

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT
end

function body(i, S_arr)
    S = read(S_arr, i)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i+1, S_new)
    return i+1, S_arr
end

xy = fem_nodes(mesh)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)
t0 = @. t_exact(x,y)

xy = fvm_nodes(mesh)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(nelem+3*ndof))

i = constant(1, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, nelem+3*ndof))


text:
using Setfield

"""
    hyperoptim(f, candidates, algorithm = Optim.NelderMead(), opts = Optim.Options(); R = 30, η = 3, inner = RandomSampler(), threads = true, verbose = false, Rmin = 0.1, kwargs...)

Multi-start optimization with Optim and Hyperband. The interface is similar to that of `Optim.optimize`, but the initial guess is replaced by a vector of candidate vectors. The resources are used to set the allowed time for optimization by `Optim.optimize`.

# Arguments:
- `f` a function from `Rⁿ -> R`, i.e., it takes a vector of parameters.
- `candidates`: A vector of vectors
- `algorithm`: One of the Optim algorithms
- `opts`: `Optim.Options`
- `inner`: The inner sampler
- `threads`: Use threading to evaluate the cost function
- `verbose`: DESCRIPTION
- `Rmin`: The minimum amount of resources to use.
- `kwargs`: Are sent to `Optim.optimize`.

See `?hyperband` for the remaining arguments.
"""
function hyperoptim(f, candidates, algorithm = Optim.NelderMead(), opts = Optim.Options(); R=30, η=3, inner = RandomSampler(), threads=true, verbose = false, Rmin = 0.1, kwargs...)

    fun = function(r, pars)
        r = max(r, Rmin)
        @set! opts.time_limit = r
        res = Optim.optimize(f, pars, algorithm, opts; kwargs...)
        verbose && println("Resources: $(round(r, digits=3)), value: $(round(minimum(res), digits=6))")
        Optim.minimum(res), Optim.minimizer(res)
    end

    ho = hyperband(fun, candidates; R, η, inner, threads)
    ho
end


text:
function chebyshev_evaluate(weights::Array{T,1},x::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  polynomial_1 = chebyshev_polynomial(order[1],[x1])

  evaluated_polynomial = zero(T)

  for i = 1:order[1]+1

    evaluated_polynomial += polynomial_1[i]*weights[i]

  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,2},x::Array{T,1},order::Array{S,1},domain=[ones(1,2);-ones(1,2)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  polynomial_1 = chebyshev_polynomial(order[1],[x1])
  polynomial_2 = chebyshev_polynomial(order[2],[x2])

  evaluated_polynomial = zero(T)

  for j = 1:order[2]+1
    for i = 1:order[1]+1

      evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*weights[i,j]

    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,3},x::Array{T,1},order::Array{S,1},domain=[ones(1,3);-ones(1,3)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

  polynomial_1 = chebyshev_polynomial(order[1],[x1])
  polynomial_2 = chebyshev_polynomial(order[2],[x2])
  polynomial_3 = chebyshev_polynomial(order[3],[x3])

  evaluated_polynomial = zero(T)

  for k = 1:order[3]+1
    for j = 1:order[2]+1
      for i = 1:order[1]+1

        evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*weights[i,j,k]

      end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,4},x::Array{T,1},order::Array{S,1},domain=[ones(1,4);-ones(1,4)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

  polynomial_1 = chebyshev_polynomial(order[1],[x1])
  polynomial_2 = chebyshev_polynomial(order[2],[x2])
  polynomial_3 = chebyshev_polynomial(order[3],[x3])
  polynomial_4 = chebyshev_polynomial(order[4],[x4])

  evaluated_polynomial = zero(T)

  for l = 1:order[4]+1
    for k = 1:order[3]+1
      for j = 1:order[2]+1
        for i = 1:order[1]+1

          evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*weights[i,j,k,l]

  			 end
			end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,5},x::Array{T,1},order::Array{S,1},domain=[ones(1,5);-ones(1,5)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  polynomial_1 = chebyshev_polynomial(order[1],[x1])
  polynomial_2 = chebyshev_polynomial(order[2],[x2])
  polynomial_3 = chebyshev_polynomial(order[3],[x3])
  polynomial_4 = chebyshev_polynomial(order[4],[x4])
  polynomial_5 = chebyshev_polynomial(order[5],[x5])

  evaluated_polynomial = zero(T)

	for m = 1:order[5]+1
	  for l = 1:order[4]+1
      for k = 1:order[3]+1
        for j = 1:order[2]+1
          for i = 1:order[1]+1

            evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*polynomial_5[m]*weights[i,j,k,l,m]

					end
				end
      end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,6},x::Array{T,1},order::Array{S,1},domain=[ones(1,6);-ones(1,6)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  polynomial_1 = chebyshev_polynomial(order[1],[x1])
  polynomial_2 = chebyshev_polynomial(order[2],[x2])
  polynomial_3 = chebyshev_polynomial(order[3],[x3])
  polynomial_4 = chebyshev_polynomial(order[4],[x4])
  polynomial_5 = chebyshev_polynomial(order[5],[x5])
  polynomial_6 = chebyshev_polynomial(order[6],[x6])

  evaluated_polynomial = zero(T)

	for n = 1:order[6]+1
	  for m = 1:order[5]+1
	    for l = 1:order[4]+1
        for k = 1:order[3]+1
          for j = 1:order[2]+1
            for i = 1:order[1]+1

              evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*polynomial_5[m]*polynomial_6[n]*weights[i,j,k,l,m,n]

            end
          end
				end
      end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,1},x::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  polynomial_1 = chebyshev_polynomial(order,[x1])

  evaluated_polynomial = zero(T)

  for i = 1:order+1

    evaluated_polynomial += polynomial_1[i]*weights[i]

  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,2},x::Array{T,1},order::S,domain=[ones(1,2);-ones(1,2)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  polynomial_1 = chebyshev_polynomial(order,[x1])
  polynomial_2 = chebyshev_polynomial(order,[x2])

  evaluated_polynomial = zero(T)

  for j = 1:order+1
    for i = 1:order+1

      if (i+j <= order+2)

        evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*weights[i,j]

      end

    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,3},x::Array{T,1},order::S,domain=[ones(1,3);-ones(1,3)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

  polynomial_1 = chebyshev_polynomial(order,[x1])
  polynomial_2 = chebyshev_polynomial(order,[x2])
  polynomial_3 = chebyshev_polynomial(order,[x3])

  evaluated_polynomial = zero(T)

  for k = 1:order+1
    for j = 1:order+1
      for i = 1:order+1

        if (i+j+k <= order+3)

          evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*weights[i,j,k]

        end

      end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,4},x::Array{T,1},order::S,domain=[ones(1,4);-ones(1,4)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

  polynomial_1 = chebyshev_polynomial(order,[x1])
  polynomial_2 = chebyshev_polynomial(order,[x2])
  polynomial_3 = chebyshev_polynomial(order,[x3])
  polynomial_4 = chebyshev_polynomial(order,[x4])

  evaluated_polynomial = zero(T)

  for l = 1:order+1
    for k = 1:order+1
      for j = 1:order+1
        for i = 1:order+1

          if (i+j+k+l <= order+4)

            evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*weights[i,j,k,l]

          end

  			end
			end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,5},x::Array{T,1},order::S,domain=[ones(1,5);-ones(1,5)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  polynomial_1 = chebyshev_polynomial(order,[x1])
  polynomial_2 = chebyshev_polynomial(order,[x2])
  polynomial_3 = chebyshev_polynomial(order,[x3])
  polynomial_4 = chebyshev_polynomial(order,[x4])
  polynomial_5 = chebyshev_polynomial(order,[x5])

  evaluated_polynomial = zero(T)

  for m = 1:order+1
    for l = 1:order+1
      for k = 1:order+1
	      for j = 1:order+1
          for i = 1:order+1

            if (i+j+k+l+m <= order+5)

              evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*polynomial_5[m]*weights[i,j,k,l,m]

            end

					end
				end
      end
    end
  end

  return evaluated_polynomial

end

function chebyshev_evaluate(weights::Array{T,6},x::Array{T,1},order::S,domain=[ones(1,6);-ones(1,6)]) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  polynomial_1 = chebyshev_polynomial(order,[x1])
  polynomial_2 = chebyshev_polynomial(order,[x2])
  polynomial_3 = chebyshev_polynomial(order,[x3])
  polynomial_4 = chebyshev_polynomial(order,[x4])
  polynomial_5 = chebyshev_polynomial(order,[x5])
  polynomial_6 = chebyshev_polynomial(order,[x6])

  evaluated_polynomial = zero(T)

  for n = 1:order+1
    for m = 1:order+1
      for l = 1:order+1
        for k = 1:order+1
	        for j = 1:order+1
            for i = 1:order+1

              if (i+j+k+l+m+n <= order+6)

                evaluated_polynomial += polynomial_1[i]*polynomial_2[j]*polynomial_3[k]*polynomial_4[l]*polynomial_5[m]*polynomial_6[n]*weights[i,j,k,l,m,n]

              end

            end
					end
				end
      end
    end
  end

  return evaluated_polynomial

end

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
@doc """
A ZeroOneRollCall object represents a roll call data-set in a format
appropriate for binary classification analysis: 'yea' votes are encoded as 1's;
'nay' votes are encoded as 0's; all other vote values are encoded as `NaN`.
""" ->
immutable ZeroOneRollCall
    votes::Matrix{Float64}
end

@doc """
Translate a single integer-valued vote from a RollCall object that uses the
provided `encoding` into the zero-one encoding used by the ZeroOneRollCall
type.
""" ->
function transcode(vote::Integer, encoding::RollCallEncoding)
    ifelse(
        vote in encoding.nay,
        0.0,
        ifelse(
            vote in encoding.yea,
            1.0,
            NaN
        )
    )
end

@doc """
Convert a RollCall object into a ZeroOneRollCall object.
""" ->
function Base.convert(::Type{ZeroOneRollCall}, roll_call::RollCall)
    encoding, old_votes = roll_call.encoding, roll_call.votes
    n_legislators, n_bills = size(old_votes)
    votes = Array(Float64, n_legislators, n_bills)
    for j in 1:n_bills
        for i in 1:n_legislators
            votes[i, j] = transcode(old_votes[i, j], encoding)
        end
    end
    ZeroOneRollCall(votes)
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
module HHModel

using DifferentialEquations

include("./Dynamics.jl")
include("./ChannelZoo.jl")
include("./Misc.jl")
end


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 1, 12, :y, Triangle(15, 2), false),
  :properties => [1.0e5 0.2;],
  :x_coords => [0.0,  1.0,  6.0],
  :y_coords => [0.0,  -2.0],
  :support => [
    (1, [0 1]),
    (2, [0 1]),
    (3, [0 1]),
    (4, [0 1]),
    (5, [0 0]),
    (10, [0 0]),
    (15, [0 0]),
    (20, [0 0]),
    (25, [0 0]),
    (30, [0 0]),
    (35, [0 0]),
    (40, [0 0]),
    (41, [0 1]),
    (42, [0 1]),
    (43, [0 1]),
    (44, [0 1]),
    (45, [0 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -7/90]),
    (6, [0.0 -16/45]),
    (11, [0.0 -2/15]),
    (16, [0.0 -16/45]),
    (21, [0.0 -7/90])
    ]
)

@time m, dis_df, fm_df = p51(data)

@test m.loads[0:9] ≈ [0.0,-1.5912198206764464e-5,-1.1580746207864741e-5,
  -7.225530531847905e-6,-3.3533531045918466e-6,-9.31933942793864e-7,
    -1.559031671692096e-5,1.4927888685769456e-7,-1.1275746358774959e-5,
    4.539223317246215e-7]


text:
module SpecialFunctions

using OpenSpecFun_jll

export
    airyai,
    airyaiprime,
    airybi,
    airybiprime,
    airyaix,
    airyaiprimex,
    airybix,
    airybiprimex,
    besselh,
    besselhx,
    besseli,
    besselix,
    besselj,
    sphericalbesselj,
    besselj0,
    besselj1,
    besseljx,
    besselk,
    besselkx,
    bessely,
    sphericalbessely,
    bessely0,
    bessely1,
    besselyx,
    dawson,
    ellipk,
    ellipe,
    erf,
    erfc,
    erfcinv,
    erfcx,
    erfi,
    erfinv,
    logerfc,
    logerfcx,
    eta,
    digamma,
    invdigamma,
    polygamma,
    trigamma,
    gamma_inc,
    beta_inc,
    beta_inc_inv,
    gamma_inc_inv,
    ncbeta,
    ncF,
    hankelh1,
    hankelh1x,
    hankelh2,
    hankelh2x,
    zeta,
    sinint,
    cosint,
    lbinomial

include("bessel.jl")
include("erf.jl")
include("ellip.jl")
include("sincosint.jl")
include("gamma.jl")
include("gamma_inc.jl")
include("betanc.jl")
include("beta_inc.jl")
include("deprecated.jl")

for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, :logerfc, :logerfcx,
          :eta, :gamma, :invdigamma, :logfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
    @eval $(f)(::Missing) = missing
end
for f in (:beta, :lbeta)
    @eval $(f)(::Number, ::Missing) = missing
    @eval $(f)(::Missing, ::Number) = missing
    @eval $(f)(::Missing, ::Missing) = missing
end
polygamma(m::Integer, x::Missing) = missing

end # module