text:
module Model

using Random
using Agents
using DataFrames
using ImageFiltering
using CairoMakie: RGBf0, RGBAf0
using InteractiveDynamics
using DelimitedFiles

name = "Rice-Brown Plant Hopper"

"""
	neighbors_at(n::Integer)

Return the directions on the grid.
"""
function neighbors_at(n::Integer)
    k = 2 * n + 1
    kern = zeros(Int, k, k)
    center = k ÷ 2 + 1
    dist = map(CartesianIndices(kern)) do I
        i, j = Tuple(I)
        center - (abs(i - center) + abs(j - center))
    end
    dist[center, center] = 0
    map(findall(dist .> 0)) do I
        (I[1] - center, I[2] - center)
    end
end

const default_parameters = (#
    energy_miss=0.025,
    age_init=168,
    age_reproduce=504,
    age_old=600,
    age_die=720,
    pr_reproduce=Dict(true => 0.188, false => 0.157),
    pr_egg_death=0.0025,
    pr_old_death=0.04,
    offspring_max=12,
    offspring_min=5,
    energy_max=1.0,
    energy_transfer=0.1,
    energy_consume=0.025,
    energy_move=0.2,
    energy_reproduce=0.8,
    move_directions=Dict(true => neighbors_at(1), false => neighbors_at(2)),
)

"""
	gencrop_3x3(::Type{T})

Generate a 3x3 rice map with type `T`.
"""
function gencrop_3x3(T::DataType=Float32)
    flower_position = [31:35; 61:65]
    nan = convert(T, NaN)
    one_ = one(T)
    return food = [
        begin
            if x in flower_position || y in flower_position
                nan
            else
                one_
            end
        end for (x, y) in Iterators.product(1:100, 1:100)
    ]
end

Base.@kwdef mutable struct BPH <: AbstractAgent
    id::Int
    pos::Dims{2}
    energy::Float16
    age::Int
    isfemale::Bool
    isshortwing::Bool
end

function init_model(envmap, nb_bph_init, init_position, pr_killed; seed, kwargs...)
    return init_model(; envmap, nb_bph_init, init_position, pr_killed, seed, kwargs...)
end
function init_model(; envmap, nb_bph_init::Int, init_position, pr_killed, seed, kwargs...)
    rng = MersenneTwister(seed)
    food = collect(transpose(init_envmap(envmap)))
    pr_killed = init_pr_killed(pr_killed, food)
    init_position = Symbol(init_position)

    # PROPERTIES

    params = merge(default_parameters, kwargs)
    properties = (#
        food=food,
        total_bph=nb_bph_init,
        death_natural=0,
        death_predator=0,
        pr_killed=pr_killed,
        pr_killed_positions=convert.(Tuple, findall(!iszero, pr_killed)),
        energy_full=1.0 - params.energy_transfer,
        params...,
    )

    # MODEL

    space = GridSpace(size(food); periodic=false)
    scheduler = Schedulers.by_id
    model = ABM(BPH, space; scheduler=scheduler, properties=properties, rng=rng)

    # AGENTS CREATION

    positions =
        let p = if init_position === :corner
                Iterators.product(1:5, 1:5)
            elseif init_position === :random_c1
                Iterators.product(1:(size(food, 1) ÷ 2), 1:size(food, 1))
            elseif init_position === :random_c2
                Iterators.product(1:(size(food, 1) ÷ 3), 1:size(food, 1))
            elseif init_position === :border
                Iterators.product(1:5, 1:size(food, 1))
            else
                @assert false "Postition not in [:corner, :random_c1, :random_c2, :border]"
            end
            filter(pos -> !isnan(food[pos...]), collect(p))
        end
    for _ in 1:nb_bph_init
        isshortwing = rand(model.rng, Bool)
        bph = BPH(; #
            id=nextid(model),
            pos=rand(model.rng, positions),
            energy=rand(model.rng, 0.4:0.01:0.6),
            age=rand(model.rng, 0:300),
            isfemale=rand(model.rng, Bool),
            isshortwing=isshortwing,
        )
        add_agent_pos!(bph, model)
    end

    # RETURN
    return model
end

"""
	init_envmap(filepath::abstractstring)

Return the food map
"""
function init_envmap(filepath::AbstractString)
    content = read(filepath, String)
    if occursin(",", content)
        readdlm(filepath, ',')
    else
        readdlm(filepath)
    end
end

"""
	init_envmap(envmap::AbstractMatrix)

Guard function
"""
function init_envmap(envmap::AbstractMatrix)
    return envmap
end

"""
	init_pr_killed(pr_killed::Real, food; gauss=2.5)

Return matrix of death pr. Obtained by using gauss kernel to filter the food matrix.
"""
function init_pr_killed(pr_killed::Real, food; gauss=2.5)
    return imfilter(isnan.(food) * pr_killed, Kernel.gaussian(gauss))
end
"""
	init_pr_killed(pr_killed::AbstractMatrix)

Identity, guard function
"""
function init_pr_killed(pr_killed::AbstractMatrix, args...)
    return pr_killed
end

# Agents behaviors
function agent_step!(agent, model)
    # position
    x, y = agent.pos

    # Older
    agent.age = agent.age + 1

    # Step wise energy loss
    agent.energy = agent.energy - (agent.age ≥ model.age_init) * model.energy_consume

    # Move conditionally
    if (agent.age ≥ model.age_init && agent.energy ≥ model.energy_move) && (
        agent.energy ≥ model.energy_full ||
        isnan(model.food[x, y]) ||
        rand(model.rng) > (model.food[x, y] * 0.5)
    )
        thres = rand(model.rng)
        directions = filter(model.move_directions[agent.isshortwing]) do (dx, dy)
            food = get(model.food, (x + dx, y + dy), -1.0)
            thres ≤ (isnan(food) / 2 + !isnan(food) * food)
        end
        if isempty(directions)
            walk!(agent, rand(model.rng, model.move_directions[agent.isshortwing]), model)
        else
            walk!(agent, rand(model.rng, directions), model)
        end
    end

    # Eat conditionally
    if model.food[x, y] > 0 && agent.age ≥ model.age_init
        transfer = min(#
            model.energy_transfer,
            model.food[x, y],
            model.energy_max - agent.energy,
        )
        model.food[x, y] -= transfer
        agent.energy += transfer
        #min(agent.energy + transfer, model.energy_max)
    end

    # Reproduce conditionally
    if (
        agent.isfemale && # is female
        agent.age ≥ model.age_reproduce && # Old enough
        agent.energy ≥ model.energy_reproduce && # Energy requirement
        rand(model.rng) ≤ model.pr_reproduce[agent.isshortwing] # Have RNG Jesus by your side
    )
        nb_offspring = rand(model.rng, (model.offspring_min):(model.offspring_max))
        isshortwing = rand(model.rng, Bool)
        for _ in 1:nb_offspring
            id = nextid(model)
            agent = BPH(;
                id=id,
                pos=agent.pos,
                energy=0.4,
                age=0,
                isfemale=rand(model.rng, Bool),
                isshortwing=isshortwing,
            )
            add_agent_pos!(agent, model)
        end
        agent.energy -= 0.1
    end

    # Die conditionally
    if (agent.energy ≤ 0) || # Exausted
       (agent.age ≥ model.age_die) || # Too old
       (
           model.age_die > agent.age ≥ model.age_old && # Old
           rand(model.rng) ≤ model.pr_old_death # And weak
       ) ||
       (
           agent.age < model.age_init && # Young
           rand(model.rng) ≤ model.pr_egg_death # And weak
       ) # then
        kill_agent!(agent, model)
        return nothing
    end

    # End of agent step
end

# Environment behaviors
function model_step!(model)
    # Rice getting more energy
    # Energy cap is 1.0
    # Dead rice receive no energy
    @. model.food = min(model.food + model.food * 0.008f0 * (model.food > 0), 1)

    # Randomly select a bph at every flower
    for pos in model.pr_killed_positions
        pos = Tuple(pos)
        if isempty(pos, model)
            continue
        elseif rand(model.rng) < model.pr_killed[pos...]
            agents_to_kill = collect(agents_in_position(pos, model))
            n = length(agents_to_kill)
            perm = randperm(model.rng, n)
            for i in Iterators.take(perm, rand(model.rng, 1:3))
                @inbounds kill_agent!(agents_to_kill[i], model)
            end
        end
    end
end

# Plotting ultilities
function ac(model)
    return function (agent)
        if isnan(model.food[agent.pos...])
            return RGBf0(0, 0, 0)
        elseif agent.age < model.age_init
            RGBf0(0.0f0, 0.0f0, 1.0f0)
        else
            RGBf0(1.0f0, 0.0f0, 0.0f0)
        end
    end
end

agent_markers = Dict(true => :circle, false => :utriangle)
function am(agent)
    return agent_markers[agent.isshortwing]
end

function heatarray(model)
    return model.food
end

function video(crop, nb_bph_init, position, pr_killed0; seed, kwargs...)
    return video(
        "BPH-$(nb_bph_init)-$(position)-$(pr_killed0)-$(seed).mp4",
        crop,
        nb_bph_init,
        position,
        pr_killed0;
        seed=seed,
        kwargs...,
    )
end

function video(
    videopath::String, crop, nb_bph_init, position, pr_killed0; seed, frames=2880, kwargs...
)
    @info "Video seed: $seed"
    model = init_model(crop, nb_bph_init, position, pr_killed0; seed=seed, kwargs...)
    return abm_video(
        videopath,
        model,
        agent_step!,
        model_step!;#
        frames=frames,
        framerate=24,
        ac=ac(model),
        am=am,
        heatarray=heatarray,
        heatkwargs=(
            nan_color=RGBAf0(1.0, 1.0, 0.0, 0.5),
            colormap=[RGBAf0(0, 1.0, 0, i) for i in 0:0.01:1],
            colorrange=(0, 1),
        ),
    )
end

# Data collection

adata, mdata = let food(model) = count(model.food .≥ 0.5), bph(agent) = agent.energy > 0
    [(bph, count)], [food]
end

post_process = function (adf, mdf)
    return rightjoin(adf, mdf; on=:step)
end

# END MODULE
end

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
module MLJModels

import MLJModelInterface

using ScientificTypes
const ScientificTypesBase = ScientificTypes.ScientificTypesBase

# still needed b/s UnivariateFinite type used in ThresholdPredictor:
using MLJBase

using Pkg, Pkg.TOML, OrderedCollections, Parameters
using Tables, CategoricalArrays, StatsBase, Statistics, Dates
import Distributions
import REPL # stdlib, needed for Term

# from loading.jl:
export load, @load, @iload, @loadcode, info

# from model_search:
export models, localmodels, matching

# from model/Constant
export ConstantRegressor, ConstantClassifier,
        DeterministicConstantRegressor, DeterministicConstantClassifier

# from model/ThresholdPredictors
export BinaryThresholdPredictor

# from model/Transformers
export FeatureSelector, StaticTransformer, UnivariateDiscretizer,
    UnivariateStandardizer, Standardizer, UnivariateBoxCoxTransformer,
    OneHotEncoder, ContinuousEncoder, FillImputer, UnivariateFillImputer,
    UnivariateTimeTypeToContinuous

const srcdir = dirname(@__FILE__) # the directory containing this file
const MMI = MLJModelInterface

if VERSION < v"1.3"
    nonmissingtype = ScientificTypes.ScientificTypesBase.nonmissing
end

nonmissing = nonmissingtype

include("utilities.jl")

Handle = NamedTuple{(:name, :pkg), Tuple{String,String}}
(::Type{Handle})(name,string) = NamedTuple{(:name, :pkg)}((name, string))

# load utilities for reading model metadata from file:
include("metadata.jl")

# read in the metadata:
metadata_file = joinpath(srcdir, "registry", "Metadata.toml")
const INFO_GIVEN_HANDLE = info_given_handle(metadata_file)
const PKGS_GIVEN_NAME = pkgs_given_name(INFO_GIVEN_HANDLE)
const AMBIGUOUS_NAMES = ambiguous_names(INFO_GIVEN_HANDLE)
const NAMES = model_names(INFO_GIVEN_HANDLE)
const MODEL_TRAITS_IN_REGISTRY = model_traits_in_registry(INFO_GIVEN_HANDLE)

# model search and registry code:
include("model_search.jl")
include("loading.jl")
include("registry/src/info_dict.jl")
include("registry/src/Registry.jl")
include("registry/src/check_registry.jl")
import .Registry.@update

# load built-in models:
include("builtins/Constant.jl")
include("builtins/Transformers.jl")
include("builtins/ThresholdPredictors.jl")

end # module

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""WebhookClientConfig contains the information to make a TLS connection with the webhook

    IoK8sApiAdmissionregistrationV1WebhookClientConfig(;
        caBundle=nothing,
        service=nothing,
        url=nothing,
    )

    - caBundle::Vector{UInt8} : `caBundle` is a PEM encoded CA bundle which will be used to validate the webhook's server certificate. If unspecified, system trust roots on the apiserver are used.
    - service::IoK8sApiAdmissionregistrationV1ServiceReference : `service` is a reference to the service for this webhook. Either `service` or `url` must be specified.  If the webhook is running within the cluster, then you should use `service`.
    - url::String : `url` gives the location of the webhook, in standard URL form (`scheme://host:port/path`). Exactly one of `url` or `service` must be specified.  The `host` should not refer to a service running in the cluster; use the `service` field instead. The host might be resolved via external DNS in some apiservers (e.g., `kube-apiserver` cannot resolve in-cluster DNS as that would be a layering violation). `host` may also be an IP address.  Please note that using `localhost` or `127.0.0.1` as a `host` is risky unless you take great care to run this webhook on all hosts which run an apiserver which might need to make calls to this webhook. Such installs are likely to be non-portable, i.e., not easy to turn up in a new cluster.  The scheme must be \"https\"; the URL must begin with \"https://\".  A path is optional, and if present may be any string permissible in a URL. You may use the path to pass an arbitrary string to the webhook, for example, a cluster identifier.  Attempting to use a user or basic auth e.g. \"user:password@\" is not allowed. Fragments (\"#...\") and query parameters (\"?...\") are not allowed, either.
"""
mutable struct IoK8sApiAdmissionregistrationV1WebhookClientConfig <: SwaggerModel
    caBundle::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: caBundle
    service::Any # spec type: Union{ Nothing, IoK8sApiAdmissionregistrationV1ServiceReference } # spec name: service
    url::Any # spec type: Union{ Nothing, String } # spec name: url

    function IoK8sApiAdmissionregistrationV1WebhookClientConfig(;caBundle=nothing, service=nothing, url=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1WebhookClientConfig, Symbol("caBundle"), caBundle)
        setfield!(o, Symbol("caBundle"), caBundle)
        validate_property(IoK8sApiAdmissionregistrationV1WebhookClientConfig, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(IoK8sApiAdmissionregistrationV1WebhookClientConfig, Symbol("url"), url)
        setfield!(o, Symbol("url"), url)
        o
    end
end # type IoK8sApiAdmissionregistrationV1WebhookClientConfig

const _property_map_IoK8sApiAdmissionregistrationV1WebhookClientConfig = Dict{Symbol,Symbol}(Symbol("caBundle")=>Symbol("caBundle"), Symbol("service")=>Symbol("service"), Symbol("url")=>Symbol("url"))
const _property_types_IoK8sApiAdmissionregistrationV1WebhookClientConfig = Dict{Symbol,String}(Symbol("caBundle")=>"Vector{UInt8}", Symbol("service")=>"IoK8sApiAdmissionregistrationV1ServiceReference", Symbol("url")=>"String")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1WebhookClientConfig }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1WebhookClientConfig))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1WebhookClientConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1WebhookClientConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1WebhookClientConfig }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1WebhookClientConfig[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1WebhookClientConfig)
    true
end

function validate_property(::Type{ IoK8sApiAdmissionregistrationV1WebhookClientConfig }, name::Symbol, val)
    if name === Symbol("caBundle")
    end
end


text:
"blue notebook p123 :
This is the derivation of KI with mohr-coulomb friction, allowing free parametrization of ψ
!!! wrong derivation !!!"
# function compute_N1(r::Rheology,D)
#     α = cosd(r.ψ)
#     1/(((D/r.D₀)^(1/3)-1)*α + r.β)^(3/2)
# end

# function compute_N2(r::Rheology,D)
#     α = cosd(r.ψ)
#     2*sqrt(α)*((D-r.D₀)/r.D₀)^(1/6)
# end

# function compute_N3(r::Rheology,D)
#     α = cosd(r.ψ)
#     r.D₀^(2/3)/(α^2*(1-D^(2/3)))
# end

# compute_N1N2N3(r,D) = (compute_N1(r,D), compute_N2(r,D), compute_N3(r,D))

# compute_M1(r) = 0.5*sind(r.ψ)*( r.μ*(1-cosd(2*r.ψ)) - sind(2*r.ψ) )
# compute_M3(r) = 0.5*sind(r.ψ)*( r.μ*(1+cosd(2*r.ψ)) + sind(2*r.ψ) )
# compute_M1M3(r) = (compute_M1(r), compute_M3(r))
# # eq 15 Bhat2012 (A1 : *c2*c3), Perol&Bhat2016 (A1 : ...*c2)*c3):
# # Perol&Bhat2016 is the corrected version, and the one implemented
# compute_A1(N1,N2,N3,M1) = (N1 + N2*N3)*M1
# compute_A3(N1,N2,N3,M3) = (N1 + N2*N3)*M3 + N2

# compute_A1A3(N1,N2,N3,M1,M3) = (compute_A1(N1,N2,N3,M1), compute_A3(N1,N2,N3,M3))
# compute_A1A3(r::Rheology,D) = compute_A1A3(compute_N1N2N3(r,D)...,compute_M1M3(r)...)

# compute_KI(r, σ₁, σ₃, A1, A3) = sqrt(r.a/π) * (A1*σ₁ + A3*σ₃)
# compute_KI(r, σ₁, σ₃, D) = compute_KI(r, σ₁, σ₃, compute_A1A3(r,D)...)

# Bhat 2011 with KI twisted to match extreme cases where l=0 and l>>1 :
function compute_c1(r,D)
    α = cosd(r.ψ)
    return 1/(π^2 * α^(3/2) * ((D/r.D₀)^(1/3) - 1 + r.β/α)^(3/2))
end

function compute_c2(r,D)
    α = cosd(r.ψ)
    return 2 * ((D/r.D₀)^(1/3) - 1)^(1/2) * (r.D₀^(2/3)/(1-D^(2/3))) / (π^2 * α^(3/2))
end

function compute_c3(r,D)
    α = cosd(r.ψ)
    return 2/π * sqrt(α) * ((D/r.D₀)^(1/3) - 1)^(1/2)
end
compute_c1c2c3(r::Rheology,D) = (compute_c1(r,D), compute_c2(r,D), compute_c3(r,D))

compute_A1(r) = π * sqrt(r.β/3) * (sqrt(1+r.μ^2) - r.μ) #TODO: Check if + or - μ makes sense
compute_A3(r) = compute_A1(r) * (sqrt(1+r.μ^2) + r.μ) / (sqrt(1+r.μ^2) - r.μ)
compute_A1A3(r::Rheology) = (compute_A1(r), compute_A3(r))

"principal stress formulation of KI from Bhat 2011"
compute_KI(r::Rheology, σ₁, σ₃, A1, A3, c1, c2, c3) = sqrt(π*r.a) * ((σ₃*A3 - σ₁*A1) * (c1 + c2) + σ₃*c3)
compute_KI(r::Rheology, σ₁, σ₃, D) = compute_KI(r, σ₁, σ₃, compute_A1A3(r)..., compute_c1c2c3(r,D)...)

compute_KI_no_interaction(r::Rheology, σ₁, σ₃, A1, A3, c1, c2, c3) = sqrt(π*r.a) * ((σ₃*A3 - σ₁*A1)*c1 + σ₃*c3)
compute_KI_no_interaction(r::Rheology, σ₁, σ₃, D) = compute_KI_no_interaction(r, σ₁, σ₃, compute_A1A3(r)..., compute_c1c2c3(r,D)...)
#################
# use those instead to drop the match with extreme cases, vizualizing the difference shows a massive underestime of KI,
# it should not be used.

# c1 when droping Beta in the wedging force term of KI
function compute_c1psi(r,D)
    α = cosd(r.ψ)
    return 1/(π^2 * α^(3/2) * (D/r.D₀)^(1/2))
end
compute_c1c2c3psi(r::Rheology,D) = (compute_c1psi(r,D), compute_c2(r,D), compute_c3(r,D))
# mohr coulomb formulation of A1 & A3 :
compute_A1psi(r) = 0.5*π * sind(r.ψ) * (r.μ*(cosd(2*r.ψ) - 1) + sind(2*r.ψ))
compute_A3psi(r) = 0.5*π * sind(r.ψ) * (r.μ*(cosd(2*r.ψ) + 1) + sind(2*r.ψ))
compute_A1A3psi(r::Rheology) = (compute_A1psi(r), compute_A3psi(r))
##################

max_stars_repo_name: Leooop/SimpleDamage

max_issues_repo_name: Leooop/SimpleDamage

max_forks_repo_name: Leooop/SimpleDamage


text:
# The random parser executes random valid moves to test the cost fn.
#
# The public interface for rparse takes the following arguments:
#
# pt::ParserType: ArcHybridR1, ArcEagerR1, ArcHybrid13, ArcEager13
# s::Sentence or c::Corpus: input sentence(s), a single parser is returned for s, a vector of parsers for c
# ncpu::Integer: (optional) performs parallel processing
# oracle: (kwarg) probability of following the oracle

function rparse{T<:Parser}(pt::Type{T}, s::Sentence; o...)
    p = pt(s)
    rparse(p, s; o...)
end

function rparse{T<:Parser}(pt::Type{T}, c::Corpus; o...)
    pa = map(pt, c)
    rparse(pa, c; o...)
end

function rparse(p::Parser, s::Sentence; oracle=0)
    c = Array(Cost, p.nmove)
    totalcost = 0
    while anyvalidmoves(p)
        movecosts(p, s.head, s.deprel, c)
        if rand() < oracle
            r = indmin(c)
        else
            f = find(c .< typemax(Cost))
            r = f[rand(1:length(f))]
        end
        totalcost += c[r]
        move!(p, r)
    end
    # truecost only makes sense if s has oracle parse
    if oracle > 0 && totalcost != truecost(p, s)
        error("Bad cost estimate")
    end
    return p
end

function rparse{T<:Parser}(pa::Vector{T}, c::Corpus; o...)
    for i=1:length(c); rparse(pa[i], c[i]; o...); end
    return pa
end

# Not tested
function rparse_dbg{T<:Parser}(pt::Type{T}, s::Sentence)
    dump(vcat([1:wcnt(s)]', s.head', s.deprel'))
    p = pt(s)
    c = Array(Position, p.nmove)
    cost = 0
    while anyvalidmoves(p)
        movecosts(p, s.head, s.deprel, c)
        if rand() < 0.1
            f = find(c .< typemax(Cost))
            r = f[rand(1:length(f))]
        else
            r = indmin(c)
        end
        cost += c[r]
        move!(p, r)
        println("$((r,Int(c[r]),cost)) $(Int(p.stack[1:p.sptr])) $(p.wptr) $(Int(c))")
    end
    @assert cost == truecost(p,s)
    return p
end

# Not tested
function rparse{T<:Parser}(pt::Type{T}, c::Corpus, ncpu::Integer; o...)
    (ncpu == 1) && return rparse(pt, c)
    @date Main.resetworkers(ncpu)
    sa = distribute(c)          # distributed sentence array
    pa = map(pt, sa)            # distributed parser array
    @sync for p in procs(sa)
        @spawnat p rparse(localpart(pa), localpart(sa), f; o...)
    end
    pa = convert(Vector{pt}, pa)
    @date Main.rmworkers()
    return pa
end


text:
using DelimitedFiles
using Plots; gr();
using LaTeXStrings
# export DelimitedFiles, readdlm, writedlm

l25 = readdlm("25ncells0.2t.txt", Float64)
l50 = readdlm("50ncells0.2t.txt", Float64)
l100 = readdlm("100ncells0.2t.txt", Float64)

nl25 = readdlm("25ncells0.2tnosl.txt", Float64)
nl50 = readdlm("50ncells0.2tnosl.txt", Float64)
nl100 = readdlm("100ncells0.2tnosl.txt", Float64)

as25 = readdlm("25ncells0.8t.txt", Float64)
as50 = readdlm("50ncells0.8t.txt", Float64)
as100 = readdlm("100ncells0.8t.txt", Float64)

l50r = zeros(25)
l100r = zeros(25)
nl50r = zeros(25)
nl100r = zeros(25)
as50r = zeros(25)
as100r = zeros(25)
for i in 1:length(l25)
    l50r[i] = l50[2i-1]
    l100r[i] = l100[4i-3]
    nl50r[i] = nl50[2i-1]
    nl100r[i] = nl100[4i-3]
    as50r[i] = as50[2i-1]
    as100r[i] = as100[4i-3]
end
xt = range(0,stop=0.96,length=25)

plot(xt,[l25,l50r,l100r],label=["25Cells" "50Cells" "100Cells"],title="Plot at 3 different resolution at time = 0.2 with the slope limiter")
savefig("3plotwithlimiter.png")

plot(xt,[nl25,nl50r,nl100r],label=["25Cells" "50Cells" "100Cells"],title="Plot at 3 different resolution at time = 0.2 without the slope limiter")
savefig("3plotwithlimiternolimiter.png")

plot(xt,[l50r-l25,2(l100r-l50r)],label=["dx/2-dx" "2(dx/4-dx/2)"],title="order of convergence at time = 0.2 with the slope limiter, a=2")
savefig("orderofconvergencelimiter.png")

plot(xt,[nl50r-nl25,2(nl100r-nl50r)],label=["dx/2-dx" "2(dx/4-dx/2)"],title="order of convergence at time = 0.2 without the slope limiter, a=2")
savefig("orderofconvergencenolimiter.png")

plot(xt,[as25,as50r,as100r],label=["25Cells" "50Cells" "100Cells"],title="Plot at 3 different resolution at time = 0.8")
savefig("3plotlatertime.png")

plot(xt,[as50r-as25,2(as100r-as50r)],label=["dx/2-dx" "2(dx/4-dx/2)"],title="order of convergence at time = 0.8, with a=2")
savefig("orderofconvergencelatertime.png")
# l


text:
include("getdigit_n_xc.jl")
include("simpson_2.jl")
include("simpson_noparallel.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n_xc
using .Simpson_2
using .Simpson_noparallel
using .MyPrintf

const DIGIT = 17
const MAXDH = BigFloat(1) / BigFloat(500)
const MAXN = 1000

function gauss(xc::BigFloat, dh::BigFloat)
    res = Simpson_2.simpson(x::BigFloat -> exp(- x * x), BigFloat(0), xc, dh)

    return BigFloat(4) * res * res
end

function gauss(xc::BigFloat, n::Int64)
    res = Simpson_noparallel.simpson(x::BigFloat -> exp(- x * x), BigFloat(0), xc, n)
    return BigFloat(4) * res * res
end

function main()
    xc = GetDigit_n_xc.getxc(gauss, MAXDH, DIGIT, true)
    n = GetDigit_n_xc.getn_even_linear(gauss, xc, MAXN, DIGIT, true)
    precision = GetDigit_n_xc.getdigit(gauss, xc, n, DIGIT, true)
    precision = GetDigit_n_xc.getdigit_2(gauss, xc, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, gauss(xc, n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("xc = %.14f, n = %d, 計算%d桁で、%d桁求まりました\n", xc, n, precision, DIGIT)
end

main()


text:
using DifferentialEquations
using Test
using MRIgeneralizedBloch
using MRIgeneralizedBloch:
    apply_hamiltonian_gbloch_superlorentzian!, apply_hamiltonian_freeprecession!, apply_hamiltonian_gbloch!

##
max_error = 5e-2

## set parameters
ω1 = π / 500e-6
ω0 = 200.0
B1 = 1.0
m0s = 0.15
R1f = 0.3
R1s = 2.0
R2f = 1 / 65e-3
T2s = 10e-6
Rx = 30.0
TRF = 500e-6
mfun(p, t; idxs = nothing) = typeof(idxs) <: Number ? 0.0 : zeros(30)
alg = MethodOfSteps(DP8())
N = Inf

# ApproxFun
G = interpolate_greens_function(greens_superlorentzian, 0, 100)
dGdT2s = interpolate_greens_function(dG_o_dT2s_x_T2s_superlorentzian, 0, 100)

## baseline IDE solution
m0 = [0.5 * (1 - m0s), 0.0, 0.5 * (1 - m0s), m0s, 1.0]
gBloch_sol = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0, TRF),
        (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, N),
    ),
    alg,
)

## Analytical gradients (using ApproxFun)
grad_list = [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]
m0 = zeros(5 * (length(grad_list) + 1), 1)
m0[1] = 0.5 * (1 - m0s)
m0[3] = 0.5 * (1 - m0s)
m0[4] = m0s
m0[5] = 1.0

gBloch_sol_grad = solve(
    DDEProblem(
        apply_hamiltonian_gbloch!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G, dGdT2s, grad_list),
    ),
    alg,
)

## FD derivative wrt. m0s
dm0s = 1e-9

m0 = m0[1:5]
gBloch_sol_dm0s = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, (m0s + dm0s), R1f, R2f, Rx, R1s, T2s, N),
    ),
    alg,
)

t = 0:1e-5:TRF
dxf = similar(t)
dyf = similar(t)
dzf = similar(t)
dzs = similar(t)
dxf_fd = similar(t)
dyf_fd = similar(t)
dzf_fd = similar(t)
dzs_fd = similar(t)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dm0s(t[i])[1] - gBloch_sol(t[i])[1]) / dm0s
    dyf_fd[i] = (gBloch_sol_dm0s(t[i])[2] - gBloch_sol(t[i])[2]) / dm0s
    dzf_fd[i] = (gBloch_sol_dm0s(t[i])[3] - gBloch_sol(t[i])[3]) / dm0s
    dzs_fd[i] = (gBloch_sol_dm0s(t[i])[4] - gBloch_sol(t[i])[4]) / dm0s

    dxf[i] = gBloch_sol_grad(t[i])[6]
    dyf[i] = gBloch_sol_grad(t[i])[7]
    dzf[i] = gBloch_sol_grad(t[i])[8]
    dzs[i] = gBloch_sol_grad(t[i])[9]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. R1f
dR1f = 1e-9

gBloch_sol_dR1f = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, (R1f + dR1f), R2f, Rx, R1s, T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dR1f(t[i])[1] - gBloch_sol(t[i])[1]) / dR1f
    dyf_fd[i] = (gBloch_sol_dR1f(t[i])[2] - gBloch_sol(t[i])[2]) / dR1f
    dzf_fd[i] = (gBloch_sol_dR1f(t[i])[3] - gBloch_sol(t[i])[3]) / dR1f
    dzs_fd[i] = (gBloch_sol_dR1f(t[i])[4] - gBloch_sol(t[i])[4]) / dR1f

    dxf[i] = gBloch_sol_grad(t[i])[11]
    dyf[i] = gBloch_sol_grad(t[i])[12]
    dzf[i] = gBloch_sol_grad(t[i])[13]
    dzs[i] = gBloch_sol_grad(t[i])[14]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test [dxf;dyf;dzf;dzs] ≈ [dxf_fd;dyf_fd;dzf_fd;dzs_fd] rtol = max_error

## test derivative wrt. R2f
dR2f = 1e-9

gBloch_sol_dR2f = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1f, (R2f + dR2f), Rx, R1s, T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dR2f(t[i])[1] - gBloch_sol(t[i])[1]) / dR2f
    dyf_fd[i] = (gBloch_sol_dR2f(t[i])[2] - gBloch_sol(t[i])[2]) / dR2f
    dzf_fd[i] = (gBloch_sol_dR2f(t[i])[3] - gBloch_sol(t[i])[3]) / dR2f
    dzs_fd[i] = (gBloch_sol_dR2f(t[i])[4] - gBloch_sol(t[i])[4]) / dR2f

    dxf[i] = gBloch_sol_grad(t[i])[16]
    dyf[i] = gBloch_sol_grad(t[i])[17]
    dzf[i] = gBloch_sol_grad(t[i])[18]
    dzs[i] = gBloch_sol_grad(t[i])[19]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error

## test derivative wrt. Rx
dRx = 1e-6

gBloch_sol_dRx = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1f, R2f, (Rx + dRx), R1s, T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dRx(t[i])[1] - gBloch_sol(t[i])[1]) / dRx
    dyf_fd[i] = (gBloch_sol_dRx(t[i])[2] - gBloch_sol(t[i])[2]) / dRx
    dzf_fd[i] = (gBloch_sol_dRx(t[i])[3] - gBloch_sol(t[i])[3]) / dRx
    dzs_fd[i] = (gBloch_sol_dRx(t[i])[4] - gBloch_sol(t[i])[4]) / dRx

    dxf[i] = gBloch_sol_grad(t[i])[21]
    dyf[i] = gBloch_sol_grad(t[i])[22]
    dzf[i] = gBloch_sol_grad(t[i])[23]
    dzs[i] = gBloch_sol_grad(t[i])[24]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. R1s
dR1s = 1e-5

gBloch_sol_dR1s = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1f, R2f, Rx, (R1s + dR1s), T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dR1s(t[i])[1] - gBloch_sol(t[i])[1]) / dR1s
    dyf_fd[i] = (gBloch_sol_dR1s(t[i])[2] - gBloch_sol(t[i])[2]) / dR1s
    dzf_fd[i] = (gBloch_sol_dR1s(t[i])[3] - gBloch_sol(t[i])[3]) / dR1s
    dzs_fd[i] = (gBloch_sol_dR1s(t[i])[4] - gBloch_sol(t[i])[4]) / dR1s

    dxf[i] = gBloch_sol_grad(t[i])[26]
    dyf[i] = gBloch_sol_grad(t[i])[27]
    dzf[i] = gBloch_sol_grad(t[i])[28]
    dzs[i] = gBloch_sol_grad(t[i])[29]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. T2s
dT2s = 1e-14

gBloch_sol_dT2s = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, (T2s + dT2s), N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dT2s(t[i])[1] - gBloch_sol(t[i])[1]) / dT2s
    dyf_fd[i] = (gBloch_sol_dT2s(t[i])[2] - gBloch_sol(t[i])[2]) / dT2s
    dzf_fd[i] = (gBloch_sol_dT2s(t[i])[3] - gBloch_sol(t[i])[3]) / dT2s
    dzs_fd[i] = (gBloch_sol_dT2s(t[i])[4] - gBloch_sol(t[i])[4]) / dT2s

    dxf[i] = gBloch_sol_grad(t[i])[31]
    dyf[i] = gBloch_sol_grad(t[i])[32]
    dzf[i] = gBloch_sol_grad(t[i])[33]
    dzs[i] = gBloch_sol_grad(t[i])[34]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

## test derivative wrt. ω0
dω0 = 1

gBloch_sol_dω0 = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, (ω0 + dω0), m0s, R1f, R2f, Rx, R1s, T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dω0(t[i])[1] - gBloch_sol(t[i])[1]) / dω0
    dyf_fd[i] = (gBloch_sol_dω0(t[i])[2] - gBloch_sol(t[i])[2]) / dω0
    dzf_fd[i] = (gBloch_sol_dω0(t[i])[3] - gBloch_sol(t[i])[3]) / dω0
    dzs_fd[i] = (gBloch_sol_dω0(t[i])[4] - gBloch_sol(t[i])[4]) / dω0

    dxf[i] = gBloch_sol_grad(t[i])[36]
    dyf[i] = gBloch_sol_grad(t[i])[37]
    dzf[i] = gBloch_sol_grad(t[i])[38]
    dzs[i] = gBloch_sol_grad(t[i])[39]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error

## test derivative wrt. B1
dB1 = 1e-9

gBloch_sol_dB1 = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, (B1 + dB1), ω0, m0s, R1f, R2f, Rx, R1s, T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dB1(t[i])[1] - gBloch_sol(t[i])[1]) / dB1
    dyf_fd[i] = (gBloch_sol_dB1(t[i])[2] - gBloch_sol(t[i])[2]) / dB1
    dzf_fd[i] = (gBloch_sol_dB1(t[i])[3] - gBloch_sol(t[i])[3]) / dB1
    dzs_fd[i] = (gBloch_sol_dB1(t[i])[4] - gBloch_sol(t[i])[4]) / dB1

    dxf[i] = gBloch_sol_grad(t[i])[41]
    dyf[i] = gBloch_sol_grad(t[i])[42]
    dzf[i] = gBloch_sol_grad(t[i])[43]
    dzs[i] = gBloch_sol_grad(t[i])[44]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error

######################################################################
## apparent R1
######################################################################
R1a = 1

## baseline IDE solution
m0 = [0.5 * (1 - m0s), 0.0, 0.5 * (1 - m0s), m0s, 1.0]
gBloch_sol = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0, TRF),
        (ω1, B1, ω0, m0s, R1a, R2f, Rx, R1a, T2s, N),
    ),
    alg,
)

## Analytical gradients (using ApproxFun)
grad_list = [grad_R1a()]
m0 = zeros(5 * (length(grad_list) + 1), 1)
m0[1] = 0.5 * (1 - m0s)
m0[3] = 0.5 * (1 - m0s)
m0[4] = m0s
m0[5] = 1.0

gBloch_sol_grad = solve(
    DDEProblem(
        apply_hamiltonian_gbloch!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, R1a, R2f, Rx, R1a, T2s, G, dGdT2s, grad_list),
    ),
    alg,
)

## test derivative wrt. R1a
dR1a = 1e-9
m0 = m0[1:5]

gBloch_sol_dR1a = solve(
    DDEProblem(
        apply_hamiltonian_gbloch_superlorentzian!,
        m0,
        mfun,
        (0.0, TRF),
        (ω1, B1, ω0, m0s, (R1a + dR1a), R2f, Rx, (R1a + dR1a), T2s, N),
    ),
    alg,
)

for i = 1:length(t)
    dxf_fd[i] = (gBloch_sol_dR1a(t[i])[1] - gBloch_sol(t[i])[1]) / dR1a
    dyf_fd[i] = (gBloch_sol_dR1a(t[i])[2] - gBloch_sol(t[i])[2]) / dR1a
    dzf_fd[i] = (gBloch_sol_dR1a(t[i])[3] - gBloch_sol(t[i])[3]) / dR1a
    dzs_fd[i] = (gBloch_sol_dR1a(t[i])[4] - gBloch_sol(t[i])[4]) / dR1a

    dxf[i] = gBloch_sol_grad(t[i])[6]
    dyf[i] = gBloch_sol_grad(t[i])[7]
    dzf[i] = gBloch_sol_grad(t[i])[8]
    dzs[i] = gBloch_sol_grad(t[i])[9]
end

@test dxf ≈ dxf_fd rtol = max_error
@test dyf ≈ dyf_fd rtol = max_error
@test dzf ≈ dzf_fd rtol = max_error
@test dzs ≈ dzs_fd rtol = max_error * 10

## plot the functions for debug purposes
# using Plots
# plotlyjs()
# plot(t, dxf, ticks=:native, label="xf", legend=:bottomleft)
# plot!(t, dxf_fd, label="xf_fd")
# plot!(t, dyf, ticks=:native, label="yf")
# plot!(t, dyf_fd, label="yf_fd")
# plot!(t, dzf, ticks=:native, label="zf")
# plot!(t, dzf_fd, label="zf_fd")
# plot!(t, dzs, label="zs")
# plot!(t, dzs_fd, label="zs_fd")


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

################################################################################
## function dftranspose(df)
##
## inputs: df - [DataFrame] a DataFrame object to be transposed
## results t - [DataFrame] a transposed version of the df DataFrame.
##		   withhead - [Boolean] if True, first column of df will become column
##		   names for t. Otherwise, first column will first row and column names
##		   will be generic (e.g. x1:xN)
##
## Note this function is necessary because no stock function to transpose
## DataFrames appears to exist.
################################################################################
"""
   df = dftranspose(df::DataFrame, withhead::Bool)

Returns a transpose of a Dataframe.\n
FIXME: This is for DataFrames@0.20.2, as used in GenX.
Versions 0.21+ could use stack and unstack to make further changes while retaining the order
"""
function dftranspose(df::DataFrame, withhead::Bool)
	if withhead
		colnames = cat(:Row, Symbol.(df[!,1]), dims=1)
		return DataFrame([[names(df)]; collect.(eachrow(df))], colnames)
	else
		return DataFrame([[names(df)]; collect.(eachrow(df))], [:Row; Symbol.("x",axes(df, 1))])
	end
end # End dftranpose()

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# f(x)=A*X
function FAX(A::Array,X::Array)
    if IsScalarArray(A)
        return (A[1].*X,nothing)
    elseif IsScalarArray(X)
        return (A.*X[1],nothing)
    else
        return (A*X,nothing)
    end
end

function FAX_inplace(value,auxvalue,A::Array,X::Array)
    if IsScalarArray(A)
        copy!(value,A[1]*X)
    elseif IsScalarArray(X)
        copy!(value,A*X[1])
    else
        BLAS.gemm!('N','N',1.0,A,X,0.0,value)
    end
end

function DAX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X)
    if derivativeIDX==1
        if IsScalarArray(A)
            axpy!(1.0,[sum(X.*grad_c)],grad_n)
        elseif IsScalarArray(X)
            axpy!(X[1],grad_c,grad_n)
        else
            BLAS.gemm!('N','T',1.0,grad_c,X,1.0,grad_n) # grad_c*X'
        end
    elseif derivativeIDX==2
        if IsScalarArray(A)
            axpy!(A[1],grad_c,grad_n)
        elseif IsScalarArray(X)
            axpy!(1.0,[sum(A.*grad_c)],grad_n)
        else
            BLAS.gemm!('T','N',1.0,A,grad_c,1.0,grad_n) # A'*grad_c
        end
    end
end

if PROC=="GPU"
    gemm!(T1::Char,T2::Char,alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},beta::Float64,C::CudaArray{Float64})=CUBLAS.gemm!(T1,T2,alpha,A,B,beta,C)
    gemm!(T1::Char,T2::Char,alpha,A::CudaArray{Float32},B::CudaArray{Float32},beta,C::CudaArray{Float32})=CUBLAS.gemm!(T1,T2,Float32(alpha),A,B,Float32(beta),C)
    export gemm!

    function FAX_inplace(value::CudaArray,auxvalue,A::CudaArray,X::CudaArray)
        if IsScalarArray(A)
            copy!(value,X); scale!(A,value) # nb argument converse of Base.scale!
        elseif IsScalarArray(X)
            copy!(value,A); scale!(X,value) # nb argument converse of Base.scale!
        else
            #CUBLAS.gemm!('N','N',1.0,A,X,0.0,value)
            gemm!('N','N',1.0,A,X,0.0,value)
        end
    end

    function DAX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray)
        if derivativeIDX==1
            if IsScalarArray(A)
                tmp=CudaArray(Float64,size(X))
                vmult!(1.0,X,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
                #axpy!(1.0,sum(tmp),grad_n);
                free(tmp)
            elseif IsScalarArray(X)
                alphaaxpy!(1.0,X,grad_c,grad_n)
            else
                #CUBLAS.gemm!('N','T',1.0,grad_c,X,1.0,grad_n) # grad_c*X'
                gemm!('N','T',1.0,grad_c,X,1.0,grad_n) # grad_c*X'
            end
        elseif derivativeIDX==2
            if IsScalarArray(A)
                alphaaxpy!(1.0,A,grad_c,grad_n)
            elseif IsScalarArray(X)
                tmp=CudaArray(Float64,size(A))
                vmult!(1.0,A,grad_c,tmp)
                axpy!(1.0,sum(tmp),grad_n); free(tmp)
            else
                #CUBLAS.gemm!('T','N',1.0,A,grad_c,1.0,grad_n) # A'*t
                gemm!('T','N',1.0,A,grad_c,1.0,grad_n) # A'*t
            end
        end
    end

    function DAX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float32},X::CudaArray{Float32})
        if derivativeIDX==1
            if IsScalarArray(A)
                tmp=CudaArray(Float32,size(X))
                vmult!(1.0,X,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
                #axpy!(1.0,sum(tmp),grad_n);
                free(tmp)
            elseif IsScalarArray(X)
                alphaaxpy!(1.0,X,grad_c,grad_n)
            else
                #CUBLAS.gemm!('N','T',1.0,grad_c,X,1.0,grad_n) # grad_c*X'
                gemm!('N','T',1.0,grad_c,X,1.0,grad_n) # grad_c*X'
            end
        elseif derivativeIDX==2
            if IsScalarArray(A)
                alphaaxpy!(1.0,A,grad_c,grad_n)
            elseif IsScalarArray(X)
                tmp=CudaArray(Float32,size(A))
                vmult!(1.0,A,grad_c,tmp)
                axpy!(1.0,sum(tmp),grad_n); free(tmp)
            else
                #CUBLAS.gemm!('T','N',1.0,A,grad_c,1.0,grad_n) # A'*t
                gemm!('T','N',1.0,A,grad_c,1.0,grad_n) # A'*t
            end
        end
    end

end

Derivative[FAX]=DAX
export FAX

Inplace[FAX]=FAX_inplace

Derivative[FAX]=DAX
import Base.*
*(A::ADnode,B::ADnode)=ADnode(FAX,[A B])

*(A::Real,B::ADnode)=ADnode(FAX,[ADconst(A) B])
*(A::ADnode,B::Real)=ADnode(FAX,[A ADconst(B)])

@gpu *(A::CudaArray,B::CudaArray)=CUBLAS.gemm('N','N',A,B)
export *

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Pathfinder, Random, Test, Turing

Random.seed!(0)

@model function regression_model(x, y)
    σ ~ truncated(Normal(); lower=0)
    α ~ Normal()
    β ~ filldist(Normal(), size(x, 2))
    y_hat = muladd(x, β, α)
    y .~ Normal.(y_hat, σ)
    return (; y)
end

@testset "Turing integration" begin
    x = 0:0.01:1
    y = sin.(x) .+ randn.() .* 0.2 .+ x
    X = [x x .^ 2 x .^ 3]
    model = regression_model(X, y)
    expected_param_names = Symbol.(["α", "β[1]", "β[2]", "β[3]", "σ"])

    result = pathfinder(model; ndraws=10_000)
    @test result isa PathfinderResult
    @test result.input === model
    @test size(result.draws) == (5, 10_000)
    @test result.draws_transformed isa MCMCChains.Chains
    @test result.draws_transformed.info.pathfinder_result isa PathfinderResult
    @test sort(names(result.draws_transformed)) == expected_param_names
    @test all(>(0), result.draws_transformed[:σ])
    init_params = Vector(result.draws_transformed.value[1, :, 1])
    chns = sample(model, NUTS(), 10_000; init_params, progress=false)
    @test mean(chns).nt.mean ≈ mean(result.draws_transformed).nt.mean rtol = 0.1

    result = multipathfinder(model, 10_000; nruns=4)
    @test result isa MultiPathfinderResult
    @test result.input === model
    @test size(result.draws) == (5, 10_000)
    @test length(result.pathfinder_results) == 4
    @test result.draws_transformed isa MCMCChains.Chains
    @test result.draws_transformed.info.pathfinder_result isa MultiPathfinderResult
    @test sort(names(result.draws_transformed)) == expected_param_names
    @test all(>(0), result.draws_transformed[:σ])
    init_params = Vector(result.draws_transformed.value[1, :, 1])
    chns = sample(model, NUTS(), 10_000; init_params, progress=false)
    @test mean(chns).nt.mean ≈ mean(result.draws_transformed).nt.mean rtol = 0.1
end


text:
module ALGCM1D

using Plots
using Plots.PlotMeasures
using LinearAlgebra
using Dates
using Dierckx
using Measures
using Parameters
using Statistics

include("constants.jl")
include("grid.jl")
include("parameters.jl")
include("time.jl")

include("heaviside.jl")
include("delta.jl")
include("angle_of_incidence.jl")
include("rad.jl")
include("surface.jl")
include("holtslag.jl")
include("turbulence.jl")

include("atm_parameters.jl")
include("hydraulic_variables.jl")
include("soil_parameters.jl")
include("Temperature_dist.jl")
include("metric_head.jl")
include("soil.jl")

include("initialize.jl")
include("atmosphere.jl")

p1 = plot(
    1,
    #xlim = (0, 360*24*2),
    #ylim = (0, 50),
    marker = 2,
    legend=false
)
p2 = plot(
    1,
    #xlim = (0, 360*24*2),
    #ylim = (0.01, 0.015),
    marker = 2,
    legend=false
)

function init()
  include("initialize.jl")
end

function run(stime)

  global mtime=stime # model clock
  global jd,t
  global plt

  Q=zeros(n)
  sw=zeros(n)
  lw=zeros(n)
  lh=zeros(n)
  sh=zeros(n)

  for i=1:n

      ########################################################################
      # time
      ########################################################################

      epochdays=Dates.date2epochdays(Date(mtime))
      t=(Dates.datetime2epochms(mtime)/1000-epochdays*24*60*60)/3600;
      #julian day
      jd=Dates.datetime2julian(mtime)-Dates.datetime2julian(DateTime(year(mtime)));

      println(i,"   ",mtime)
      ########################################################################
      # Radiation
      # Input: date, lat and lon
      # Output: longwave and shortwave radiation
      # Future input: cloud fraction and chemestry
      ########################################################################
      rad(mtime)
      #println("SW")
      #println(SW)
      #println("---------")

      if rem(i,360)==0;
        #p1=scatter([i],[TS[end]])
        #display(p1)
        p2=scatter!([Dates.format.(mtime,"dd:HH")],[hpbl],xlabel="time [DAY:HOUR]",ylabel="PBL height",xrot=60)
        display(p2)
      end
      ########################################################################
      # Surface
      # Input from land: surface volumetric water content and temperature
      # Input from atmosphere: surface specific humidity, temperature and wind
      # Output: latent heat, sensible heat and wind stress.
      ########################################################################
      surface(mtime)

      ########################################################################
      # turbulence
      # Input: temperature, specific humidity, air densiy, surface wind speed,
      # latent heat, sensible heat, wind stress
      # Output: vertical diffusion for momentum, heat and water
      ########################################################################
      turbulence(mtime)

      ########################################################################
      # Soil
      # Input: fluxes
      ########################################################################
      soil(mtime)

      ########################################################################
      # Atmosphere
      # Input: fluxes, albedo
      ########################################################################
      atmosphere(mtime)

      if rem(i,3600000)==0;#3600*6/dt
        title = plot(title = mtime, grid = false, axis = false,
        showaxis = false, bottom_margin = -10pt,yaxis=nothing,titlefont = 100)
        #=
        println("UA")
        println(UA)
        println("Θv")
        println(Θv)
        println("qA")
        println(qA)
        println("KAm")
        println(KAm)
        =#
        f1=plot(ZAC,UA,ylabel=("UA"));
        f2=plot(ZAC,Θv.-273.15,ylabel=("Θv")); #ylim([19 21]);...
        f3=plot(ZAC,qA,ylabel=("q"));
        f4=plot(ZAC,KAm,ylabel=("KAm"));
        f5=plot(soil_numerical_parameters.z_s,theta_l,ylabel=("theta_l"));
        f6=plot(soil_numerical_parameters.z_s,TS.-273.15,ylabel=("TS")); #ylim([19 21]);...
        plt=plot(title,f1,f2,f3,f4,f5,f6,layout=grid(7, 1, heights=[0.1 ,0.15, 0.15, 0.15, 0.15, 0.15, 0.15]),legend=false,titlefontsize=6,
        size = (700, 1200),margin=8mm)
        display(plt)
      end
      mtime=mtime+Second(ΔT)
      Q[i]=Qnet
      sw[i]=sum(SW)
      lw[i]=LW
      lh[i]=LH
      sh[i]=SH
  end

end

end


text:
### A Nbdev script file ###
### Autogenerated file. Don't modify. ###

module Export
#export
import Pluto: Notebook, Cell, load_notebook_nobackup

#export
begin
"""
> struct Nucleus--> This is the lowest entity of a code notebook. This type represents the code cell in a Pluto notebook.
> * Fields:
>   * code--> String which makes up a code cell
"""
Base.@kwdef mutable struct Nucleus
	code::String=""
	end

"""
> Nucleus(code)-->Getter for accessing the constructer of the Nucleus type
"""
Nucleus(code) = Nucleus(code = code)
end

#export
begin
"""
> struct Nb--> a notebook having nuclei made up of code cells only.
> * Fields:
>   * nuclei--> Array of Nucleus type
>   * name--> name of notebook
"""
mutable struct Nb
	nuclei::Array{Nucleus,1}
	name::AbstractString
end

"""
> nuclei(nuclei::Array)--> Helper to collect an array of Nucleus.
"""
nuclei(nuclei::Array) = nuclei

"""
> nuclei(nuclei::Array)--> Helper to collect an name of a notebook.
"""
name(name::AbstractString) = name
end

#export
begin
"""
Headers that would appear in a code file
"""
const header = "### A Nbdev script file ###"
const subheader = "### Autogenerated file. Don't modify. ###"
end

#export
"""
This had to be done because having the cell iddelimeter as a constant gives an error during parsing the notebook
"""
celliddelim=string("#"," ╔═╡ ")

#export
"""
marks the end of a cell
"""
const cellsuff = "\n\n"

#export
"""
> struct Scrubbednb--> Represents a notebook from which all but code to be exported are scrubbed off.
> * Fields-->
>   * cells--> Vector of strings.
"""
mutable struct Scrubbednb
	cells::Array{String, 1}
end

#export
"""
> collectnuclei(notebook::Scrubbednb, marker)--> Reads a notebook, collects code from code cells and returns an array of Nucelus.
"""
function collectnuclei(nb::Scrubbednb, marker)
	nbcells=nb.cells
	nuclei=[]

    for i in 1:length(nbcells)
    	raw_code=nbcells[i]
    	nucleus=Nucleus(raw_code)

    	if startswith(raw_code, marker)
            push!(nuclei,nucleus)
        end

    end
	 nuclei
end

#export
begin
"""
>loadscrubbednb(io, path)::Scrubbednb--> These are modification of the Pluto.load_notebook methods. Scrubs the notebook of all stuff but the code to be exported.
"""
function loadscrubbednb(io, path)::Scrubbednb
    collectedcells = []

	# ignore first bits of file
    readuntil(io,celliddelim)

    last_read = ""
    while !eof(io)
        cellidstr = String(readline(io))
        if cellidstr == "Cell order:"
            break
        else
            rawcode = String(readuntil(io, celliddelim))
            # change Windows line endings to Linux
            normalisedcode = replace(rawcode, "\r\n" => "\n")
            # remove the cell appendix
            code = normalisedcode[1:prevind(normalisedcode, end, length(cellsuff))]

            readcell = code
            push!(collectedcells, readcell)
        end
    end

	Scrubbednb(collectedcells)
end

function loadscrubbednb(path::String)::Scrubbednb
    local loaded
    open(path, "r") do io
        loaded = loadscrubbednb(io, path)
    end
    loaded
end
end

#export
begin
"""
> loadnb(filename::String, marker::String)--> High level function which calls _load_nb
"""
loadnb(filename::String, marker::String) = _loadnb(filename, marker)

"""
> _load_nb(filename::String, marker::String)--> creates a scrubbed notebook and returns a curated Nb type having code to be exported.
"""
function _loadnb(fn::String, marker::String)
		scrubbednb=loadscrubbednb(fn)
		collectednuclei=nuclei(collectnuclei(scrubbednb, marker))
		Nb(collectednuclei, basename(fn))
end

end

#export
"""
> strip --> Helper fucntion to replace a substring y in a string x with blank.
"""
strip=(x::String,y) -> replace(x, y=>"")

#export
begin
"""
> savenb(io, nb)--> Reads the supplied notebook and creates an io and writes stuffs like the module name and the content to the created io.
"""
function savenb(io, nb)
    println(io, header)
    println(io, subheader)
    println(io, "")
    println(io, "module $(uppercasefirst(strip(strip(nb.name, r"[0-9_]"), r".jl")))")
	for nucleus in nb.nuclei
			println(io, nucleus.code*"\n")
    end
	#notebook
	print(io, "end")
end

"""
> savenb(nb::Nb, path::String)--> Creates a file in the supplied path with the name in the NB type.
"""
function savenb(nb::Nb, srcdir::String)
	fn=uppercasefirst(strip(nb.name, r"[0-9_]"))
	open(joinpath(srcdir, fn), "w") do io
        savenb(io, nb)
    end
end
end

#export
begin
"""
> readfilenames(nbsdir::String)--> Reads files in the directory and subdirectories in the given path. Reads only the files with ".jl" extension
"""
function readfilenames(nbsdir::String)
	files=[]
	for file in readdir(nbsdir)
			if endswith(file, ".jl")
				push!(files,joinpath(nbsdir, file))
			end
	end
	files
end

end

#export
"""
> exportfile(f::String, srcdir::String, marker::String)--> Loads the file in the supplied path and reads the cells which are marked with "#export". Then saves the notebook in the given path
"""
function exportfile(f::String, srcdir::String, marker::String)
	nb=loadnb(f, marker)
	if !isempty(nb.nuclei)
	    savenb(nb, srcdir)
	end
end

#export
"""
> exportcontent(fns::AbstractVector, srcdir::String, marker::String)--> maps the `export_file` function to each files
"""
function exportcontent(fns::AbstractVector, srcdir::String, marker::String)
	map(f->exportfile(f, srcdir, marker), fns)
end

#export
"""
> getfileextension(fn)--> get the file extensions in the pwd
"""
function getfileextension(fn)
    return fn[findlast(isequal('.'),fn):end]
end

#export
"""
> notebook2script(nbsdir::String, srcdir::String)--> Export all the code from the provided notebook directory to the given source directory
"""
function notebook2script(nbsdir::String, srcdir::String)
	exportcontent(readfilenames(nbsdir), srcdir, "#export")
end

#export
export notebook2script

end


text:
using Test

using MultivariateBases

using DynamicPolynomials

function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
    @polyvar x[1:2]
    for basis in [
        maxdegree_basis(B, x, degree),
        basis_covering_monomials(B, monomials(x, 0:degree))
    ]
        n = binomial(2 + degree, 2)
        @test length(basis) == n
        @test typeof(copy(basis)) == typeof(basis)
        @test nvariables(basis) == 2
        @test variables(basis) == x
        @test monomialtype(typeof(basis)) == monomialtype(x[1])
        @test typeof(empty_basis(typeof(basis))) == typeof(basis)
        @test length(empty_basis(typeof(basis))) == 0
        @test polynomialtype(basis, Float64) == polynomialtype(x[1], Float64)
        @test polynomial(i -> 0.0, basis) isa polynomialtype(basis, Float64)
        @test polynomial(zeros(n, n), basis, Float64) isa polynomialtype(basis, Float64)
        @test polynomial(ones(n, n), basis, Float64) isa polynomialtype(basis, Float64)
    end
end

function orthogonal_test(B::Type{<:AbstractMultipleOrthogonalBasis}, univ::Function, even_odd_separated::Bool)
    @test MultivariateBases.even_odd_separated(B) == even_odd_separated
    @polyvar x y
    univariate_x = univ(x)
    univariate_y = univ(y)

    @testset "Univariate $var" for (var, univ) in [(x, univariate_x), (y, univariate_y)]
        basis = maxdegree_basis(B, (var,), 4)
        for i in 1:5
            @test basis.polynomials[length(basis) + 1 - i] == univ[i]
        end
    end

    @testset "basis_covering_monomials" begin
        monos = basis_covering_monomials(B, monovec([x^2 * y, y^2]))
        if even_odd_separated
            exps = [(2, 1), (0, 2), (0, 1), (0, 0)]
        else
            exps = [(2, 1), (2, 0), (1, 1), (0, 2), (1, 0), (0, 1), (0, 0)]
        end
        for i in 1:length(monos)
            @test monos.polynomials[i] == univariate_x[exps[i][1] + 1] * univariate_y[exps[i][2] + 1]
        end
        monos = basis_covering_monomials(B, monovec([x^4, x^2, x]))
        if even_odd_separated
            exps = [4, 2, 1, 0]
        else
            exps = [4, 3, 2, 1, 0]
        end
        for i in 1:length(monos)
            @test monos.polynomials[i] == univariate_x[exps[i] + 1]
        end
    end
end

@testset "Monomial" begin
    include("monomial.jl")
end
@testset "Scaled" begin
    include("scaled.jl")
end
@testset "Fixed" begin
    include("fixed.jl")
end
@testset "Hermite" begin
    include("hermite.jl")
end
@testset "Laguerre" begin
    include("laguerre.jl")
end
@testset "Legendre" begin
    include("legendre.jl")
end
@testset "Chebyshev" begin
    include("chebyshev.jl")
end

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
using Statistics

abstract type Forecaster
end

function train!(
    forecaster::T, train_dict::SeriesDict, info_dict::SeriesInfoDict) where T<:Forecaster
end

function predict(
    forecaster::T, train_dict::SeriesInfoDict, info_dict::SeriesInfoDict)::SeriesDict where T<:Forecaster
end

struct ForecastScore
    smape
    mase
end

const ForecastScoreDict = Dict{String, ForecastScore}

function calc_mae(y, y_hat)
    mean(abs.(y_hat .- y))
end

function calc_smape(y, y_hat)
    200 * mean(abs.(y_hat .- y) ./ (abs.(y_hat) .+ abs.(y)))
end

function calc_mase(y_train, y_test, y_hat, freq)
    y_train_naive = [y_train[i - freq] for i in freq + 1:length(y_train)]
    mae_train_naive = calc_mae(y_train[freq + 1:end], y_train_naive)
    mae_test = calc_mae(y_test, y_hat)
    mae_test / mae_train_naive
end

function split_series(series_dict::SeriesDict, info_dict::SeriesInfoDict)
    train_data = SeriesDict()
    test_data = SeriesDict()

    for (id, y) in series_dict
        info = info_dict[id]
        pivot = length(y) - info.horizon
        train_data[id] = y[1:pivot]
        test_data[id] = y[pivot + 1:end]
    end

    train_data, test_data
end

function evaluate(forecaster::T, train_dict, test_dict, info_dict) where T<:Forecaster
    train!(forecaster, train_dict, info_dict)
    forecast_dict = predict(forecaster, train_dict, info_dict)
    score_dict = ForecastScoreDict()

    for id in keys(test_dict)
        info = info_dict[id]
        y_train = train_dict[id]
        y_test = test_dict[id]
        y_hat = forecast_dict[id]

        smape = calc_smape(y_test, y_hat)
        mase = calc_mase(y_train, y_test, y_hat, info.freq)
        score_dict[id] = ForecastScore(smape, mase)
    end

    score_dict
end

function calc_stats(score_dict::ForecastScoreDict)
    scores = values(score_dict)
    smape = mean([s.smape for s in scores])
    mase = mean([s.mase for s in scores])
    ForecastScore(smape, mase)
end


text:
import FileIO: open, close, save, @format_str, Stream
import IndirectArrays: IndirectArray
import MappedArrays: mappedarray
import Colors: colormap, RGB
import ImageContainers: ImageContainer, storeimage

const CURRENT_COLORMAP = colormap("Blues", 256)
current_colormap() = CURRENT_COLORMAP

function set_colormap(cmap::AbstractVector)
    @assert length(cmap)==256 "Number of gradation must be 256"
    @. CURRENT_COLORMAP = RGB(cmap)
end

function matshow(A::AbstractMatrix; f::Function=x->clamp(1+floor(Int, 256*x), 1, 256))
    #ImageContainer{:jlc}(IndirectArray(mappedarray(f, A), CURRENT_COLORMAP))
    storeimage(:jlc, IndirectArray(mappedarray(f, A), CURRENT_COLORMAP))
end

function addframe(io, img::ImageContainer{:jlc})
    save(Stream(format"BMP", io), img.content)
end

function openanim(func::Function, fmt::Symbol=:gif)
    fin = Pipe()
    fout = IOBuffer()

    if fmt == :mp4
        fproc = run(
            pipeline(`ffmpeg -v 0 -i - -pix_fmt yuv420p -f matroska -`,
                stdin=fin, stdout=fout),
            wait=false)

        try
            func(fin)
        finally
            close(fin)
        end

        wait(fproc)
        close(fproc)
        return ImageContainer{:mp4}(take!(fout))

    elseif fmt == :gif
        palette = tempname() * ".bmp"
        save(palette, IndirectArray(reshape(1:256, 16, 16)', SimpleHeatmaps.CURRENT_COLORMAP))

        fproc = run(
            pipeline(`ffmpeg -v 0 -i - -i $palette -lavfi paletteuse=dither=sierra2_4a -f gif -`,
                stdin=fin, stdout=fout),
            wait=false)

        try
            func(fin)
        finally
            close(fin)
        end

        wait(fproc)
        close(fproc)
        return ImageContainer{:gif}(take!(fout))

    else
        error("Only support mp4 / gif: $fmt")
    end
end


text:
export loadSubset

immutable Idx
    numDims::Int
    dims
    data
end

toString(idx::Idx) = "numDims: $(idx.numDims)\ndims: $(idx.dims)"

Base.print(io::IOBuffer, idx::Idx) = Base.print(io, toString(idx))

Base.print(idx::Idx) = Base.print(toString(idx))

readIdx{T}(::Type{T}, path) = open(path) do f
    bytes = readbytes(f, 4)

    bytes[3] != 8 && error("Unsupported data type")

    numDims = bytes[4]
    dims = Array{Int}(numDims)
    size = 1
    for i in 1:numDims
        dim = ntoh(read(f, Int32))
        size *= dim
        dims[i] = dim
    end

    bytes = readbytes(f, size)
    data = Array{T}(size)
    for i in 1:size
        data[i] = convert(T, bytes[i])
    end

    return Idx(numDims, dims, data)
end

include("MNISTData.jl")

loadSubset(ctx, expandLabels = true, frac = 0.6f0) = scope!(ctx) do this
    frac = min(frac, 0.8f0)

    cd = dirname(@__FILE__)
    dataRoot = joinpath(cd, "data")
    imageData = readIdx(Float32, joinpath(dataRoot, "images-subset"))
    labelData = readIdx(UInt32, joinpath(dataRoot, "labels-subset"))

    rIDims = reverse(imageData.dims)
    images = array(ctx, imageData.data, rIDims...)

    r = randu(ctx, Float32, 10000)
    cond = r .< frac
    trainIndices = where(cond)
    testIndices = where(!cond)

    trainImages = lookup(images, trainIndices, 2) ./ 255.0f0
    testImages = lookup(images, testIndices, 2) ./ 255.0f0

    numClasses = 10
    numTrain = dims(trainImages, 2)
    numTest = dims(testImages, 2)

    local trainLabelsArr, testLabelsArr

    if expandLabels
        trainLabels = [0.0f0 for x = 1:numClasses, y = 1:numTrain]
        testLabels = [0.0f0 for x = 1:numClasses, y = 1:numTest]

        assert(jType(trainIndices) == UInt32)
        assert(jType(testIndices) === UInt32)

        hTrainIdx = host(trainIndices)
        hTestIdx = host(testIndices)

        process = (num, indices, labels) ->
            for i in 1:num
                idx = indices[i] + 1
                label = labelData.data[idx]
                assert(label >= 0 && label <= 9)
                labels[label + 1, i] = 1.0f0
            end

        process(numTrain, hTrainIdx, trainLabels)
        process(numTest, hTestIdx, testLabels)

        trainLabelsArr = array(ctx, trainLabels)
        testLabelsArr = array(ctx, testLabels)
    else
        labels = array(ctx, labelData.data)
        trainLabelsArr = labels[trainIndices]
        testLabelsArr = labels[testIndices]
    end

    featureSize = DimT(elements(trainImages) / numTrain)
    trainFeats = transpose(moddims(trainImages, featureSize, numTrain))
    testFeats = transpose(moddims(testImages, featureSize, numTest))
    trainTarget = transpose(trainLabelsArr)
    testTarget = transpose(testLabelsArr)

    MNISTData(
        this.result(trainFeats),
        this.result(testFeats),
        this.result(trainTarget),
        this.result(testTarget))
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
"""
A regular expression parser transforming a PCRE string to a CombinedParser equivalent to the regular expression.
"""
module Regexp
using ..CombinedParsers
import TextParse
import TextParse: AbstractToken
using AutoHashEquals

using LazyStrings
import LazyStrings: reversed, reverse_index

import ..CombinedParsers: LeafParser, WrappedParser, CombinedParser, ConstantParser, Either, SideeffectParser
import ..CombinedParsers: parser, prune_captures, deepmap_parser, _deepmap_parser, print_constructor
import ..CombinedParsers: iterate_state, iterate_state_constant
import ..CombinedParsers: regex_prefix, regex_suffix, regex_inner, _regex_string, regex_string, _log_names
import ..CombinedParsers: state_type, leftof, tuple_pos, tuple_state
import ..CombinedParsers: _prevind, _nextind, _leftof, _rightof
_indexed_captures(x,a...) = x

import Base: SubString, ==

include("pcre.jl")

"""
SequenceWithCaptures ensapsulates a sequence to be parsed, and parsed captures.

This struct will allow for captures a sequence-level state.
For next version, a match-level state passed as iterate_state argument is considered.

See also [`ParserWithCaptures`](@ref)
"""
@auto_hash_equals struct SequenceWithCaptures{S,T} <: StringWrapper
    x::S
    subroutines::Vector{CombinedParser}
    captures::Vector{Vector{UnitRange{Int}}}
    names::Dict{Symbol,Vector{Int}}
    state::T
    function SequenceWithCaptures(x,subroutines, captures, names, state)
        new{typeof(x),typeof(state)}(x,subroutines, captures, names, state)
    end
    SequenceWithCaptures(x,cs::SequenceWithCaptures,state) =
        let S=typeof(x)
            new{S,typeof(state)}(x,cs.subroutines,cs.captures,cs.names,state)
        end
    SequenceWithCaptures(x,cs::SequenceWithCaptures) =
        let S=typeof(x)
            new{S,typeof(cs.state)}(x,cs.subroutines,cs.captures,cs.names,cs.state)
        end
end
import Base: empty!
Base.empty!(sequence::SequenceWithCaptures) =
    for c in sequence.captures
        Base.empty!(c)
    end
copy_captures(x::SequenceWithCaptures,state) =
    SequenceWithCaptures(x.x,x.subroutines, [ copy(c) for c in x.captures ],x.names,state)
reversed(x::SequenceWithCaptures) = SequenceWithCaptures(reversed(x.x),x)
reverse_index(x::SequenceWithCaptures,a...) = reverse_index(x.x,a...)
with_options(flags::UInt32,x::SequenceWithCaptures) =
    SequenceWithCaptures(with_options(flags,x.x),x)
function Base.show(io::IO, x::SequenceWithCaptures)
    print(io, "SequenceWithCaptures ")
    show(io,x.x)
end

with_options(set_flags::UInt32, unset_flags::UInt32,x::SequenceWithCaptures) =
    SequenceWithCaptures(
        with_options(set_flags, unset_flags,x.x),
        x
    )

"""
    ==(pcre_m::RegexMatch,pc_m::ParseMatch)

equal iif values of `.match`, `.offset`, `.ncodeunits` and `.captures` are equal.
"""
function ==(pcre_m::RegexMatch,pc_m::ParseMatch)
    pcre_m.match==pc_m.match &&
        pcre_m.match.offset==pc_m.match.offset &&
        pcre_m.match.ncodeunits==pc_m.match.ncodeunits &&
        pcre_m.captures==pc_m.captures
end
==(pc_m::ParseMatch,pcre_m::RegexMatch) =
    pcre_m==pc_m

match_string(x::SubString,y::UnitRange{<:Integer}) =
    SubString(x.string,x.offset+y.start,x.offset+y.stop)

match_string(x::Tuple{<:AbstractString,UnitRange{<:Integer}},y::UnitRange{<:Integer}) =
    let rel = min(x[2].start,x[2].stop)-1
        SubString(x[1],y.start,y.stop)
    end

match_string(x::AbstractString,y::UnitRange{<:Integer}) =
    SubString(x,y.start,y.stop)

export Capture
"""
Capture a parser result, optionally with a name.
`index` field is initialized when calling `ParserWithCaptures` on the parser.

[`ParserWithCaptures`](@ref)
"""
@auto_hash_equals struct Capture{P,S,T} <: WrappedParser{P,S,T}
    parser::P
    name::Union{Nothing,Symbol}
    index::Int
    Capture(name::Union{Nothing,Symbol},x_,index=-1) =
        let x = parser(x_)
            new{typeof(x),state_type(x),result_type(x)}(x,name==Symbol("") ? nothing : name,index)
        end
    Capture(x::Capture,index) =
        new{typeof(x.parser),state_type(x),result_type(x)}(x.parser,x.name,index)
end
Capture(x,index=-1) =
    Capture(nothing,x,index)
function print_constructor(io::IO,x::Capture)
    print_constructor(io,x.parser)
    print(io, " |> Capture ", x.index )
end

regex_string(x::Capture) = _regex_prefix(x)*regex_string(x.parser)*")"
_regex_prefix(x::Capture) =
    let name = (x.name===nothing ? "" : "?<$(x.name)>")
        "($name"
    end
regex_prefix(x::Capture) =
    _regex_prefix(x::Capture)*regex_prefix(x.parser)
regex_suffix(x::Capture) = regex_suffix(x.parser)*")"

function _deepmap_parser(f::Function,mem::AbstractDict,x::Capture,a...;kw...)
    Capture(x.name,deepmap_parser(f,mem,x.parser,a...;kw...),x.index)
end

Base.get(x::Capture, sequence, till, after, i, state) =
    get(x.parser, sequence, till, after, i, state)

@inline function iterate_state(parser::Capture, sequence, till, posi, next_i, state)
    r = iterate_state(parser.parser, sequence, till, posi, next_i, state)
    if r !== nothing ## set only if found (e.g. if repeated capture capture last)
        set_capture(sequence,parser.index,posi,_prevind(sequence,tuple_pos(r)))
    elseif state !== nothing
        prune_captures(sequence, posi)
    end
    r
end

set_capture(sequence::AbstractString, index::Int, start, stop) = nothing
set_capture(sequence::StringWithOptions, index::Int, start,stop) =
    set_capture(sequence.x,index,start,stop)
set_capture(sequence::SequenceWithCaptures, index::Int, start,stop) =
    push!((@inbounds sequence.captures[index]), start:stop)
set_capture(sequence::SequenceWithCaptures{<:ReversedString}, index::Int, start,stop) =
    push!((@inbounds sequence.captures[index]),
          reverse_index(sequence.x,
                        stop):reverse_index(sequence.x,
                                            start))

function prune_captures(sequence::SequenceWithCaptures,after_i)
    @inbounds for i in 1:length(sequence.captures)
        @inbounds cv = sequence.captures[i]
        while !isempty(cv) && (@inbounds  cv[end].stop >= after_i)
            pop!(cv)
        end
    end
end

export Backreference
"""
    Backreference(f::Function,index::Integer)

    Backreference(f::Function,name::Union{Nothing,Symbol},index::Integer)

    Backreference(f::Function,name::AbstractString)

Parser matching previously captured sequence, optionally with a name.
`index` field is recursively set when calling 'ParserWithCaptures` on the parser.
"""
@auto_hash_equals struct Backreference <: LeafParser{Int,AbstractString}
    name::Union{Nothing,Symbol}
    index::Int
    fallback::Function
    Backreference(f::Function,index::Integer) =
        new(nothing,index,f)
    Backreference(f::Function,name::Union{Nothing,Symbol},index::Integer) =
        new(name,index,f)
    Backreference(f::Function,name::AbstractString) =
        new(Symbol(name),-1,f)
end

_regex_string(x::Backreference) =
   if x.name !== nothing
       string(x.name)
   else
       string(x.index)
   end

regex_inner(x::Backreference) =
    "\\g{"*_regex_string(x) *"}"

capture_index(name,delta,index,context) =
    if ( index<0 || delta!=Symbol("") )
        if name !== nothing
            index
        else
            length(context.subroutines)+index+1
        end
    else
        index
    end

function Base.get(x::Backreference, sequence, till, after, i, state)
    sequence[i:_prevind(sequence,i+state)]
end

@inline function _leftof(str,i,parser::Backreference,x)
    i-x
end

@inline function _rightof(str,i,parser::Backreference,x)
    i+x
end

function resolve_index(p::Backreference, sequence::SequenceWithCaptures)
    index = p.index
    if index < 0 && p.name !== nothing
        for i in sequence.names[p.name]
            if !isempty(sequence.captures[i])
                return i
            end
        end
    end
    ( index<0 || isempty(sequence.captures[index]) ) ? -1 : index
end

capture_substring(p::ParserOptions{<:Backreference}, sequence::SequenceWithCaptures) =
    with_options(p.set_flags, p.unset_flags,capture_substring(p.parser, sequence))

function capture_substring(p::Backreference, sequence::SequenceWithCaptures)
    index = resolve_index(p, sequence)
    index<0 && return nothing
    SubString(sequence.x, sequence.captures[index][end])
end

@inline function iterate_state(p::Union{Backreference,ParserOptions{<:Backreference}}, sequence::SequenceWithCaptures, till, posi, next_i, state::Nothing)
    r = iterate_state_constant(
        ConstantParser(capture_substring(p, sequence)),
        sequence, till, posi, next_i, state)
    r === nothing && return nothing
    tuple_pos(r), tuple_pos(r)-next_i
end

@inline iterate_state(p::Union{Backreference,ParserOptions{<:Backreference}}, sequence::SequenceWithCaptures, till, posi, next_i, state::Int) =
    return nothing

iterate_state_condition(p::Backreference, sequence, till, posi, next_i, state) =
    resolve_index(p, sequence)>0

export Subroutine
"""
Parser matching preceding capture, optionally with a name.
`index` field is recursively set when calling `ParserWithCaptures` on the parser.
"""
@auto_hash_equals struct Subroutine{S,T} <: CombinedParser{S,T}
    name::Union{Nothing,Symbol}
    delta::Symbol
    index::Int
    Subroutine{S,T}(name::Union{Nothing,Symbol},delta::Symbol,index::Integer) where {S,T} =
        new{S,T}(name,delta,index)
    Subroutine(name::Union{Nothing,Symbol},delta::Symbol,index::Integer) =
        new{Any,Any}(name,delta,index)
    Subroutine(index::Int) =
        new{Any,Any}(nothing,Symbol(""),index)
    Subroutine(name::AbstractString) =
        new{Any,Any}(Symbol(name),Symbol(""),-1)
    Subroutine() =
        new{Any,Any}(nothing,Symbol(""),-1)
end
children(x::Subroutine) = tuple()
function regex_prefix(x::Subroutine)
    "(?" *
        if x.name !== nothing
            "&$(x.name),$(x.index)"
        else
            if x.delta == Symbol("")
                string(x.delta)
            else
                ""
            end*string(x.index)
        end
end
regex_suffix(x::Subroutine) = ")"
regex_inner(x::Subroutine) = ""

_deepmap_parser(::Function,mem::AbstractDict,x::Subroutine) = x

function iterate_state_condition(cond::Subroutine, sequence, till, posi, next_i, state)
    sequence.state === nothing && return false
    if cond.name === nothing && cond.index < 0
        true
    elseif cond.index>0
        cond.index == sequence.state.index
    else
        @show cond, sequence.state
        error()
    end
end

@inline function _leftof(sequence,i,parser::Subroutine,x)
    _leftof(sequence,i,sequence.subroutines[index(parser,sequence)].parser,x)
end

@inline function _rightof(sequence,i,parser::Subroutine,x)
    _rightof(sequence,i,sequence.subroutines[index(parser,sequence)].parser,x)
end

"""
    index(parser::Subroutine,sequence)

Index of a subroutine.
["If you make a subroutine call to a non-unique named subpattern, the one that corresponds to the first occurrence of the name is used."](https://www.pcre.org/original/doc/html/pcrepattern.html#SEC16)
(what about "In the absence of duplicate numbers (see the previous section) this is the one with the lowest number."?)
"""
index(parser::Subroutine,sequence) =
    parser.index <= 0 ? first(sequence.names[parser.name]) : parser.index

@inline function iterate_state(parser::Subroutine, sequence::SequenceWithCaptures, till, posi, next_i, state)
    iterate_state(
        sequence.subroutines[index(parser,sequence)].parser,
        copy_captures(sequence,parser), till, posi, next_i, state)
end

export DupSubpatternNumbers
"""
Parser wrapper for `ParserWithCaptures`, setting reset_index=true in `deepmap_parser(::typeof(_indexed_captures),...)`.

```jldoctest
julia> p = re"(?|(a)|(b))\\1"
🗄 Sequence |> regular expression combinator with 1 capturing groups
├─ |🗄 Either |> DupSubpatternNumbers
│  ├─ (a)  |> Capture 1
│  └─ (b)  |> Capture 1
└─ \\g{1} Backreference
::Tuple{Char, AbstractString}

julia> match(p, "aa")
ParseMatch("aa", 1="a")

julia> match(p, "bb")
ParseMatch("bb", 1="b")

```

See also [pcre doc](https://www.pcre.org/original/doc/html/pcrepattern.html#dupsubpatternnumber)
"""
@auto_hash_equals struct DupSubpatternNumbers{P,S,T} <: WrappedParser{P,S,T}
    parser::P
    DupSubpatternNumbers(parser) =
        new{typeof(parser),state_type(parser),result_type(parser)}(parser)
end

_deepmap_parser(f::Function,mem::AbstractDict,x::DupSubpatternNumbers, a...;kw...) =
    DupSubpatternNumbers(deepmap_parser(f,mem,x.parser,a...;kw...))

export Conditional
"""
Conditional parser, `iterate_state` cycles conditionally on `iterate_state_condition` through matches in field `yes` and `no` respectively.
"""
@auto_hash_equals struct Conditional{C,Y,N,S,T} <: CombinedParser{S,T}
    condition::C
    yes::Y
    no::N
    Conditional(condition,yes,no) =
        new{typeof(condition),typeof(yes),typeof(no),
            Pair{Symbol,Union{state_type(yes),state_type(no)}},
            Union{result_type(yes),result_type(no)}}(condition,yes,no)
end

function regex_prefix(x::Conditional)
    "(?("*_regex_string(x.condition)*")"
end
function regex_suffix(x::Conditional)
    ")"
end
regex_inner(x::Conditional) =
    regex_string(x.yes)*(isa(x.no, Always) ? "" : ( "|" * regex_string(x.no)))

children(x::Conditional) = x.no isa Always ? tuple(x.yes) : tuple(x.yes,x.no)

function _deepmap_parser(f::Function,mem::AbstractDict,x::Conditional,a...;kw...)
    Conditional(deepmap_parser(f,mem,x.condition,a...;kw...),
                deepmap_parser(f,mem,x.yes,a...;kw...),
                deepmap_parser(f,mem,x.no,a...;kw...))
end

@inline Base.get(parser::Conditional, sequence, till, after, i, state) =
    get(state.first == :yes ? parser.yes : parser.no, sequence, till, after, i, state.second)

iterate_state_condition(cond::WrappedParser, sequence, till, posi, next_i, state) =
    iterate_state_condition(cond.parser, sequence, till, posi, next_i, state)
iterate_state_condition(cond, sequence, till, posi, next_i, state) =
    iterate_state(cond, sequence, till, posi, next_i, state) !== nothing

@inline function _leftof(str,i,parser::Conditional,state)
    leftof(str,i,state.first == :yes ? parser.yes : parser.no, state.second)
end

@inline function _rightof(str,i,parser::Conditional,state)
    rightof(str,i,state.first == :yes ? parser.yes : parser.no, state.second)
end

@inline function iterate_state(parser::Conditional, sequence, till, posi, next_i, state::Nothing)
    c = iterate_state_condition(parser.condition, sequence, till, posi, next_i, state)
    cparse = c ? parser.yes : parser.no
    s = iterate_state(cparse,
                 sequence, till, posi, next_i, state)
    s === nothing && return nothing
    tuple_pos(s), (c ? :yes : :no) => tuple_state(s)
end

@inline function iterate_state(parser::Conditional, sequence, till, posi, next_i, state::Pair)
    iterate_state(state.first == :yes ? parser.yes : parser.no, sequence, till, posi, next_i, state.second)
end

include("indexed_captures.jl")

include("re-parser.jl")

end


text:
#=
The benchmarks are defined alongside the tests. In order to generate the
benchmark, run the following in the REPL

    julia> result = PkgBenchmark.benchmarkpkg("Commutator")
    julia>  PkgBenchmark.export_markdown("benchmark.md", result)

=#

include(joinpath(dirname(@__FILE__), "..", "test", "runtests.jl"))


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, Pkg

name = "figtree"
version = v"0.9.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://sourceforge.net/projects/figtree/files/figtree-$(version).zip", "f0b39cef8a0ab56075cce0c656b1b4327d9bd1acb31aad927241974ac4d10d82"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/figtree-*/
atomic_patch -p1 ../patches/fix-dlext.patch
make -j${nproc}
mkdir -p ${libdir}
mv "lib/libann_figtree_version.${dlext}" "${libdir}/libann_figtree_version.${dlext}"
mv "lib/libfigtree.${dlext}" "${libdir}/libfigtree.${dlext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libfigtree", :libfigtree),
    LibraryProduct("libann_figtree_version", :libann_figtree_version)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
export ExactLinearAlgebra

# COV_EXCL_START
if VERSION < v"1.5-"
    function Base._round_rational(::Type{T}, x::Rational{Tr}, ::RoundingMode{:ToZero}) where {T,Tr}
        if denominator(x) == zero(Tr) && T <: Integer
            throw(DivideError())
        elseif denominator(x) == zero(Tr)
            return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
        end
        return convert(T,div(x.num,x.den))
    end

    function Base._round_rational(::Type{T}, x::Rational{Tr}, ::RoundingMode{:Down}) where {T,Tr}
        if denominator(x) == zero(Tr) && T <: Integer
            throw(DivideError())
        elseif denominator(x) == zero(Tr)
            return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
        end
        return convert(T,fld(x.num,x.den))
    end

    function Base._round_rational(::Type{T}, x::Rational{Tr}, ::RoundingMode{:Up}) where {T,Tr}
        if denominator(x) == zero(Tr) && T <: Integer
            throw(DivideError())
        elseif denominator(x) == zero(Tr)
            return convert(T, copysign(one(Tr)//zero(Tr), numerator(x)))
        end
        return convert(T,cld(x.num,x.den))
    end
end
# COV_EXCL_STOP

module ExactLinearAlgebra
    import LinearAlgebra

    ScalarType = Union{Integer, Complex{<:Integer}, <:Rational, Complex{<:Rational}}

    function get_cofactor_matrix_unsafe!(
        out::AbstractMatrix{I1},
        mat::AbstractMatrix{I2},
        row::Integer,
        col::Integer,
    ) where {I1<:ScalarType, I2<:ScalarType}
        out[1:row-1, 1:col-1] = mat[1:row-1, 1:col-1]
        out[1:row-1, col:end] = mat[1:row-1, col+1:end]
        out[row:end, 1:col-1] = mat[row+1:end, 1:col-1]
        out[row:end, col:end] = mat[row+1:end, col+1:end]
        return out
    end

    function determinant_unsafe(mat::AbstractMatrix{S}) where {S<:ScalarType}
        n = size(mat)[1]
        n == 1 && return mat[1]
        out = Matrix{S}(undef, (n-1, n-1))
        sign = 1
        D = zero(S)
        for f = 1:n
            get_cofactor_matrix_unsafe!(out, mat, 1, f)
            @inbounds D += sign * mat[1,f] * determinant_unsafe(out)
            sign = -sign
        end
        return D
    end

    function determinant(mat::AbstractMatrix{S}) where {S<:ScalarType}
        n, m = size(mat)
        n != m && throw(ArgumentError("matrix needs to be square"))
        n <= 0 && throw(ArgumentError("matrix empty"))
        return determinant_unsafe(mat)
    end

    determinant(mat::AbstractMatrix{<:AbstractFloat}) = LinearAlgebra.det(mat)

    function inverse(mat::AbstractMatrix{S}) where {S<:ScalarType}
        n, m = size(mat)
        n != m && throw(ArgumentError("matrix needs to be square"))
        n <= 0 && throw(ArgumentError("matrix empty"))
        n == 1 && return ones(S, (1,1)) .// mat[1]
        cofactor = Array{S}(undef, (n, n))
        temp = Array{S}(undef, (n-1, n-1))
        for r in 1:n, c in 1:n
            sign = 1 - 2 * mod(r+c, 2)
            get_cofactor_matrix_unsafe!(temp, mat, r, c)
            @inbounds cofactor[r,c] = sign * determinant_unsafe(temp)
        end
        D = sum(mat[1,:] .* cofactor[1,:])
        return transpose(cofactor) // D
    end

    inverse(mat::AbstractMatrix{<:AbstractFloat}) = LinearAlgebra.inv(mat)

end # module ExactLinearAlgebra

module ExprParser

export parse_expr

SYMBOL_DATABASE = Dict(
    :i => im,
    :im => im,
    :pi => pi,
    :π => π,
    :+ => (+),  :- => (-),  :* => (*),  :/ => (/),  :\ => (\),  :^ => (^), :// => (//),
    :cos => cos,  :sin => sin,  :tan => tan,
    :exp => exp,  :cis => cis,
    :cospi => cospi,  :sinpi => sinpi,
    :sqrt => sqrt,
    :log => log,
    :abs => abs,  :abs2 => abs2,  :sign => sign,
    :conj => conj, :real => real, :imag => imag,
    :angle => angle
)

myeval(expr::Number) = expr
function myeval(expr::Symbol)
    if haskey(SYMBOL_DATABASE, expr)
        return SYMBOL_DATABASE[expr]
    else
        error("unsupported symbol $expr")
    end
end

function myeval(expr::Expr)
    if expr.head == :call
        ftn, args = Iterators.peel(myeval.(expr.args))
        return ftn(args...)
    elseif expr.head == :vect
        return [(myeval.(expr.args))...]
    elseif expr.head == :tuple
        return tuple((myeval.(expr.args))...)
    elseif expr.head == :hcat
        return hcat((myeval.(expr.args))...)
    elseif expr.head == :row
        return hcat((myeval.(expr.args))...)
    elseif expr.head == :vcat
        return vcat((myeval.(expr.args))...)
    else
        error("unsupported expression $expr")
    end
end

parse_expr(expr::Number) = expr
parse_expr(expr::AbstractString) = myeval(Meta.parse(expr))
parse_expr(expr::AbstractArray) = [parse_expr(elem) for elem in expr]

end

using .ExprParser: parse_expr

#= # commented for now
function parse_table(s::AbstractString)
    s = strip(s)
    rows = Vector{Int}[]
    for line in split(s, "\n")
        line = strip(line)
        row = [parse(Int, x) for x in split(line)]
        !isempty(row) && push!(rows, row)
    end
    return transpose(hcat(rows...))
end
=#

function cleanup_number(x::AbstractFloat, tol::Real)
    units = [16, 16/sqrt(3), 16*sqrt(3)]
    for unit in units
        if isapprox(x*unit, round(x*unit); atol=tol)
            return round(x*unit)/unit
        end
    end
    return x
end

cleanup_number(x::Integer, tol::Real) = x
cleanup_number(x::Rational{<:Integer}, tol::Real) = x

function cleanup_number(x::Complex, tol::Real)
    return cleanup_number(real(x), tol) + im * cleanup_number(imag(x), tol)
end

cleanup_number(x::AbstractArray, tol::Real) = [cleanup_number(y, tol) for y in x]

function extended_gcd(a::Integer, b::Integer)
    sign_a, sign_b = sign(a), sign(b)
    s, old_s = 0, 1
    r, old_r = abs(b), abs(a)
    while r != 0
        quotient = old_r ÷ r
        (old_r, r) = (r, old_r - quotient * r)
        (old_s, s) = (s, old_s - quotient * s)
    end
    if b != 0
        bezout_t = (old_r - old_s * abs(a)) ÷ abs(b)
    else
        bezout_t = 0
    end
    return (old_r, [sign_a * old_s, sign_b * bezout_t])
end


text:
struct HS{bits} <: Encoding
    key::Vector{UInt8}
end

"HMAC SHA-256"
const HS256 = HS{256}

"HMAC using SHA-384"
const HS384 = HS{384}

"HMAC using SHA-512"
const HS512 = HS{512}

to_byte_array(str::AbstractString) = Vector{UInt8}(Base.CodeUnits(str))

HS{bits}(key::AbstractString) where {bits} = HS{bits}(to_byte_array(key))
alg(::HS{bits}) where {bits} = "HS$(bits)"

function has_valid_signature(encoding::HS,
        header_and_claims_encoded::AbstractString,
        signature_encoded::AbstractString) :: Bool

    return signature_encoded == sign(encoding, header_and_claims_encoded)
end

sign(encoding::HS256, data::AbstractString) = base64url_encode(SHA.hmac_sha2_256(encoding.key, data))
sign(encoding::HS384, data::AbstractString) = base64url_encode(SHA.hmac_sha2_384(encoding.key, data))
sign(encoding::HS512, data::AbstractString) = base64url_encode(SHA.hmac_sha2_512(encoding.key, data))