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

@doc raw"""HorizontalPodAutoscalerSpec describes the desired functionality of the HorizontalPodAutoscaler.

    IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec(;
        maxReplicas=nothing,
        metrics=nothing,
        minReplicas=nothing,
        scaleTargetRef=nothing,
    )

    - maxReplicas::Int32 : maxReplicas is the upper limit for the number of replicas to which the autoscaler can scale up. It cannot be less that minReplicas.
    - metrics::Vector{IoK8sApiAutoscalingV2beta2MetricSpec} : metrics contains the specifications for which to use to calculate the desired replica count (the maximum replica count across all metrics will be used).  The desired replica count is calculated multiplying the ratio between the target value and the current value by the current number of pods.  Ergo, metrics used must decrease as the pod count is increased, and vice-versa.  See the individual metric source types for more information about how each type of metric must respond. If not set, the default metric will be set to 80% average CPU utilization.
    - minReplicas::Int32 : minReplicas is the lower limit for the number of replicas to which the autoscaler can scale down.  It defaults to 1 pod.  minReplicas is allowed to be 0 if the alpha feature gate HPAScaleToZero is enabled and at least one Object or External metric is configured.  Scaling is active as long as at least one metric value is available.
    - scaleTargetRef::IoK8sApiAutoscalingV2beta2CrossVersionObjectReference : scaleTargetRef points to the target resource to scale, and is used to the pods for which metrics should be collected, as well as to actually change the replica count.
"""
mutable struct IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec <: SwaggerModel
    maxReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: maxReplicas
    metrics::Any # spec type: Union{ Nothing, Vector{IoK8sApiAutoscalingV2beta2MetricSpec} } # spec name: metrics
    minReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: minReplicas
    scaleTargetRef::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2CrossVersionObjectReference } # spec name: scaleTargetRef

    function IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec(;maxReplicas=nothing, metrics=nothing, minReplicas=nothing, scaleTargetRef=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec, Symbol("maxReplicas"), maxReplicas)
        setfield!(o, Symbol("maxReplicas"), maxReplicas)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec, Symbol("metrics"), metrics)
        setfield!(o, Symbol("metrics"), metrics)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec, Symbol("minReplicas"), minReplicas)
        setfield!(o, Symbol("minReplicas"), minReplicas)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec, Symbol("scaleTargetRef"), scaleTargetRef)
        setfield!(o, Symbol("scaleTargetRef"), scaleTargetRef)
        o
    end
end # type IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec

const _property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec = Dict{Symbol,Symbol}(Symbol("maxReplicas")=>Symbol("maxReplicas"), Symbol("metrics")=>Symbol("metrics"), Symbol("minReplicas")=>Symbol("minReplicas"), Symbol("scaleTargetRef")=>Symbol("scaleTargetRef"))
const _property_types_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec = Dict{Symbol,String}(Symbol("maxReplicas")=>"Int32", Symbol("metrics")=>"Vector{IoK8sApiAutoscalingV2beta2MetricSpec}", Symbol("minReplicas")=>"Int32", Symbol("scaleTargetRef")=>"IoK8sApiAutoscalingV2beta2CrossVersionObjectReference")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec)
    (getproperty(o, Symbol("maxReplicas")) === nothing) && (return false)
    (getproperty(o, Symbol("scaleTargetRef")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerSpec }, name::Symbol, val)
end


text:
module TestFlowchart

using FlowChart
using FF, Test

@testset "testset1" begin
  @test 1==1
end

end #module

ext: jl


text:
using LightGraphs
using DataFrames
using CSV
using GraphPlot
using Compose
using ColorSchemes
using Colors

using Plots
pyplot()

function create_id_nodename_map(node_names::AbstractArray{T}) where T
  id_name_map = Dict{Int8, T}()
  for (i, c) in enumerate(node_names)
      id_name_map[i] = c
  end
  id_name_map
end

function add_edges!(G, node_edge_map, name_id_map; verbose::Bool=false)
  for (x, y) in node_edge_map
    verbose && println("($(name_id_map[x]): $x) → ($(name_id_map[y]): $y)")
    add_edge!(G, name_id_map[x], name_id_map[y])
  end
end

create_df_map(df::DataFrame, a::T, b::T) where T <: Symbol = Dict(zip(df[a], df[b]))

function team_maps(df::DataFrame)
  column_pairs = [[:killer, :killer_team_id], [:victim, :victim_team_id]]
  [create_df_map(df, a, b) for (a, b) in column_pairs]
end

get_players(df::DataFrame) = union(df[:killer], df[:victim])

function player_id_maps(df::DataFrame)
  id_player_map = create_id_nodename_map(get_players(df))
  player_id_map = map(reverse, id_player_map)
  id_player_map, player_id_map
end

function build_graph(df::DataFrame)
  node_edge_map = zip(df[:killer], df[:victim])
  id_player_map, player_id_map = player_id_maps(df)

  G = DiGraph(length(id_player_map))
  add_edges!(G, node_edge_map, player_id_map)

  G
end

function print_properties(G::AbstractGraph)
  println("is_directed: $(is_directed(G))")
  println("n_vetrices: $(nv(G))")
  println("n_edges: $(ne(G))")
  println("has_self_loops: $(has_self_loops(G))")
end

function player_kill_counts(df)
  players = get_players(df)
  player_kills = Dict(zip(players, zeros(Int8, length(players))))

  for k in df[:killer]
    player_kills[k] += 1
  end

  player_kills
end

team_colour_map(n_teams) = get(ColorSchemes.rainbow, collect(linspace(0, 1, n_teams)))

function player_team_colour_map(players, player_team_map, team_colours)

  player_team_colour = Dict{String, RGB{Float64}}()
  for player in players
    team_id = player_team_map[player]
    if team_id == 0
      player_team_colour[player] = RGB{Float64}(128/256, 128/256, 128/256)
    else
      player_team_colour[player] = team_colours[player_team_map[player]]
    end
  end

  player_team_colour
end

function generate_graph_plot(G::AbstractGraph, player_kills, id_player_map, player_team_colour, output::String)

  srand(312132)
  loc_x, loc_y = spring_layout(G)

  node_sizes = [player_kills[id_player_map[i]] + 1 for i in 1:nv(G)]  # add one to offset zero sized nodes
  # node_sizes .^= 0.3  # zoom factor
  node_sizes = node_sizes .^ 0.5

  draw(
    PNG(output, (1920 - 1920 * 0.1)*2px, (1080 - 1080 * 0.1)*2px),
    gplot(
      G,
      nodelabel=[id_player_map[i] for i in 1:nv(G)],
      nodesize=node_sizes,
      # nodelabelsize=node_sizes,
      nodefillc=[player_team_colour[id_player_map[i]] for i in 1:nv(G)],
      arrowlengthfrac=0.02,
      loc_x, loc_y
    ),
  )

end

function main()
  # match_id = "64c6a84d-49b4-4c6c-943e-26bc8676b611"
  match_id = "1859feb8-e65e-46eb-9ef0-555082002695"
  filename = "pubg_kill_events_$(match_id).csv"

  path_kill_events = "c:/workspace/pubg-analytics/output/$filename"

  df = CSV.read(path_kill_events)
  println(tail(df))

  killer_team_map, victim_team_map = team_maps(df)
  player_team_map = merge(killer_team_map, victim_team_map)

  players = get_players(df)
  id_player_map, player_id_map = player_id_maps(df)

  G = build_graph(df)
  print_properties(G)

  player_kills = player_kill_counts(df)

  n_teams = maximum(values(player_team_map))  # + 1  # for non-killer caused deaths
  team_colours = team_colour_map(n_teams)

  player_team_colour = player_team_colour_map(players, player_team_map, team_colours)

  output_path = "c:/workspace/pubg-analytics/output/pubg_kill_graph_$match_id.png"
  generate_graph_plot(G, player_kills, id_player_map, player_team_colour, output_path)

end

main()

max_stars_repo_name: rokkuran/pubg-analytics

max_issues_repo_name: rokkuran/pubg-analytics

max_forks_repo_name: rokkuran/pubg-analytics


text:
"""
    pdi{T<:Number}(x::Array{T, 1})

This function will range the values of each row, so that the strongest link has
a value of one. This works for deterministic and quantitative networks.

#### References

Poisot, T., Bever, J.D., Nemri, A., Thrall, P.H., Hochberg, M.E., 2011. A
conceptual framework for the evolution of ecological specialisation. Ecol. Lett.
14, 841–851. https://doi.org/10.1111/j.1461-0248.2011.01645.x

Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
of ecological specialization estimators. Methods in Ecology and Evolution 3,
537–544. https://doi.org/10.1111/j.2041-210X.2011.00174.x
"""
function pdi(x::T) where {T<:AbstractArray{<:Number}}
    dp = (1.0 .- sort(x/maximum(x), rev=true))[2:end]
    return sum(dp)/(length(x)-1)
end

"""
    specificity(N::DeterministicNetwork)

Measure of specificity in a deterministic network. This returns a value between
0 and 1, where 1 indicates maximal specificity.

#### References

- Poisot, T., Bever, J.D., Nemri, A., Thrall, P.H., Hochberg, M.E., 2011. A
  conceptual framework for the evolution of ecological specialisation. Ecol.
  Lett. 14, 841–851. https://doi.org/10.1111/j.1461-0248.2011.01645.x

- Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
  of ecological specialization estimators. Methods in Ecology and Evolution 3,
  537–544. https://doi.org/10.1111/j.2041-210X.2011.00174.x
"""
function specificity(N::DeterministicNetwork)
    return Dict([s => pdi(N[i,:]) for (i,s) in enumerate(species(N; dims=1))])
end


text:
using MathOptInterface
const MOI = MathOptInterface
using BARON, CPLEX, Gurobi

struct FactoryConfig
    time_limit::Float64
    num_threads::Int
    cplex_path::String
end

const primal_factory(config::FactoryConfig) = MOI.OptimizerWithAttributes(
    Gurobi.Optimizer,
    MOI.TimeLimitSec() => config.time_limit,
    MOI.NumberOfThreads() => config.num_threads,
    MOI.RawParameter("NonConvex") => 2,
    MOI.RawParameter("MIPFocus") => 1,
)

# CutUp works for minimization (which we are doing)
function cplex_direct_factory(config::FactoryConfig, cutoff::Float64)
    return MOI.OptimizerWithAttributes(
        CPLEX.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("CPXPARAM_OptimalityTarget") => 3,
        MOI.RawParameter("CPX_PARAM_CUTUP") => cutoff,
    )
end
function cplex_mip_factory(config::FactoryConfig, cutoff::Float64)
    return MOI.OptimizerWithAttributes(
        CPLEX.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("CPX_PARAM_CUTUP") => cutoff,
    )
end
function gurobi_direct_factory(config::FactoryConfig, cutoff::Float64)
    return MOI.OptimizerWithAttributes(
        Gurobi.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("NonConvex") => 2,
        MOI.RawParameter("Cutoff") => cutoff,
    )
end
function gurobi_direct_mipfocus_factory(
    config::FactoryConfig,
    cutoff::Float64,
    mipfocus::Int,
)
    return MOI.OptimizerWithAttributes(
        Gurobi.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("NonConvex") => 2,
        MOI.RawParameter("MIPFocus") => mipfocus,
        MOI.RawParameter("Cutoff") => cutoff,
    )
end
function gurobi_direct_mipfocus_3_factory(config::FactoryConfig, cutoff::Float64)
    return gurobi_direct_mipfocus_factory(config, cutoff, 3)
end

function gurobi_mip_factory(config::FactoryConfig, cutoff::Float64)
    return MOI.OptimizerWithAttributes(
        Gurobi.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("Cutoff") => cutoff,
    )
end

function gurobi_mip_mipfocus_factory(config::FactoryConfig, cutoff::Float64, mipfocus::Int)
    return MOI.OptimizerWithAttributes(
        Gurobi.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.NumberOfThreads() => config.num_threads,
        MOI.RawParameter("MIPFocus") => mipfocus,
        MOI.RawParameter("Cutoff") => cutoff,
    )
end
function gurobi_mip_mipfocus_3_factory(config::FactoryConfig, cutoff::Float64)
    return gurobi_mip_mipfocus_factory(config, cutoff, 3)
end

function baron_direct_factory(config::FactoryConfig, cutoff::Float64)
    return MOI.OptimizerWithAttributes(
        BARON.Optimizer,
        MOI.TimeLimitSec() => config.time_limit,
        MOI.RawParameter("threads") => config.num_threads,
        MOI.RawParameter("CplexLibName") => config.cplex_path,
        MOI.RawParameter("CutOff") => cutoff,
    )
end


text:
# https://juliadocs.github.io/Documenter.jl/stable/man/guide/#Package-Guide
push!(LOAD_PATH, "../src/")

using Documenter, DocHosting

makedocs(sitename = "Blogging ",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         modules = [DocHosting],
         pages = [
        "Introduction" => "index.md",
        "Article 1" => "article1.md"
    ]
)

deploydocs(
    repo = "github.com/kunov/DocHosting.jl.git",
)


text:
"""
    layerdimensions()

Retrive dimensions of a given fixed-input-size layer.
"""
layerdimensions(::Any) = (1,1)
layerdimensions(l::Flux.Dense) = size(l.W)
layerdimensions(l::Flux.RNNCell) = (size(l.Wh)[2], size(l.Wi)[2])
layerdimensions(l::Flux.LSTMCell) = (size(l.Wh)[2], size(l.Wi)[2])
layerdimensions(l::Flux.GRUCell) = (size(l.Wh)[2], size(l.Wi)[2])
layerdimensions(r::Flux.Recur) = layerdimensions(r.cell)

"""
    FIXED_INPUT_DIM_LAYERS

List of layers with fixed-sized input data
"""
const FIXED_INPUT_DIM_LAYERS = (Flux.Dense, Flux.Recur, Flux.RNNCell, Flux.LSTMCell, Flux.GRUCell) # list of types of layers with fixed input dimensions

"""
    get_dimensions(m::Flux.Chain, input_data = nothing)

Return the dimensions of the input and of the output data of each hidden layer.

If `input_data` is not given, the first layer is required to be a layer
with fixed input dimensions, such as Flux.Dense or Flux.Recur,
otherwise the given data is used to infer the dimensions of each layer.
"""
function get_dimensions(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing)

    if (input_data === nothing) & (m.layers[1] isa Union{FIXED_INPUT_DIM_LAYERS...})
        input_data = rand(Float32, layerdimensions(m.layers[1])[2])
    elseif input_data === nothing
        throw(ArgumentError("An `input_data` is required when the first layer accepts variable-dimension input"))
    else
        input_data = convert.(Float32, input_data)
    end

    chain_dimensions = vcat(size(input_data), [size(m[1:nl](input_data)) for nl in 1:length(m.layers)])
    return chain_dimensions
end

function get_dimensions(m::Flux.Chain)
    if m.layers[1] isa Union{FIXED_INPUT_DIM_LAYERS...}
        input_data = rand(Float32, layerdimensions(m.layers[1])[2])
    else
        throw(ArgumentError("An `input_data` is required when the first layer accepts variable-dimension input"))
    end
    return get_dimensions(m, input_data)
end

function get_dimensions(m::Flux.Chain, ldim::Tuple)
    input_data = rand(Float32, ldim)
    return get_dimensions(m, input_data)
end

"""
    UnitVector{T}

Structure for unit vectors in a linear space

Used for generating a basis to infer the layer connection
"""
struct UnitVector{T} <: AbstractVector{T}
    idx::Int
    length::Int
end

Base.getindex(x::UnitVector{T}, i) where T = x.idx==i ? one(T) : zero(T)
Base.length(x::UnitVector) = x.length
Base.size(x::UnitVector) = (x.length,)

"""
    neuron_indices(ldim:Tuple) -> Vector{NTuple{N, Int}} where N

Return all possible indices for a given Tuple `ldim`.
"""
function neuron_indices(ldim::Tuple)
    return [Tuple(1+mod(div(i,prod(ldim[1:j-1])), ldim[j]) for j in 1:length(ldim)) for i in 0:prod(ldim)-1]
end

"""
    neuron_connections(m::Flux.Chain) -> Vector{Dict{Tuple, Vector{Tuple}}}

Return all the connections from every neuron in each layer to the corresponding neurons in the next layer.
"""
function neuron_connections(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing)
    chain_dimensions = get_dimensions(m, input_data)
    connections = Vector{Dict{Tuple, Vector{Tuple}}}()

    for (ln, l) in enumerate(m)
        ldim = chain_dimensions[ln]
        layer_connections = Dict{Tuple,Array{Tuple,1}}()
        foreach(1:prod(ldim)) do idx
            affected = Array{Tuple,1}()
            basis_element = reshape(UnitVector{Int}(idx, prod(ldim)),ldim...)
            for rv in convert.(Float32, rand(Int16,1000))
                union!(affected, Tuple.(findall(x -> abs(x) > eps(), l(rv .* basis_element))))
            end
            push!(layer_connections, Tuple(findfirst(x->x==1, basis_element)) => affected)
        end
        push!(connections, layer_connections)
    end

    return connections
end

"""
    get_max_width(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing)

Get the maximum display width for the chain.
"""
get_max_width(m::Flux.Chain, input_data::Union{Nothing,Array} = nothing) =
    mapreduce(x->x[1], max, get_dimensions(m,input_data))


text:
using CairoMakie
using TuringGLM

function plot_chains(chns)
    params = names(chns, :parameters)

    n_chains = length(chains(chns))
    n_samples = length(chns)

    fig = Figure(; resolution=(1_000, 800))

    for (i, param) in enumerate(params)
        ax = Axis(fig[i, 1]; ylabel=string(param))
        for chain in 1:n_chains
            values = chns[:, param, chain]
            lines!(ax, 1:n_samples, values; label=string(chain))
        end

        hideydecorations!(ax; label=false)
        if i < length(params)
            hidexdecorations!(ax; grid=false)
        else
            ax.xlabel = "Iteration"
        end
    end

    for (i, param) in enumerate(params)
        ax = Axis(fig[i, 2]; ylabel=string(param))
        for chain in 1:n_chains
            values = chns[:, param, chain]
            density!(ax, values; label=string(chain))
        end

        hideydecorations!(ax)
        if i == length(params)
            ax.xlabel = "Parameter estimate"
        end
    end

    axes = [only(contents(fig[i, 2])) for i in 1:length(params)]

    return fig
end


text:
function price_elasticity(inv_sigma, df, p_points; included = ones(size(eachcol(df[:, r"shares"]),1), size(eachcol(df[:, r"shares"]),1)),
    whichProducts = [1,1], trueS = true, bO = 2 .*ones(size(included,2),1), iv::Integer = 0, deltas = -1 .* convert(Array{Float64}, df[:,r"p"]), maxes = [])
# --------------------------------------------------------------------
# Code takes in estimates from inverse_demand and returns implied price
# elasticity of first good
# --------------------------------------------------------------------

# inv_sigma     -- estimates returned by inverse_demand.m
# s             -- market shares
# p_points      -- vector of prices corresponding to s
# deltas        -- index (excluding xi) at which to evaluate demand
# b0            -- order of Bernstein polynomials
# whichProducts -- [derivative of s_j, with respect to k]
# included      -- matrix indicating which products are substitutes
# trueS         -- evaluate derivatives at realized shares or calculate
#                   counterfactual shares
# maxes         -- object with maxes.maxs = max(s), maxes.mins = min(s),
#                   used if s is small enough to cause numerical issues and needs to be rescaled

s = convert(Array{Float64}, df[:, r"shares"]);
J = size(s,2);
numBadMarkets = 0;
if size(bO,1) !=J
    if size(bO,2)==J
        bO=bO';
    elseif length(bO) ==1
        bO = bO*ones(J,1);
    else
        throw("Error: Matrix of polynomial orders has wrong number of elements")
    end
end

if trueS == true
    if size(deltas,1) ==1
        throw("Error: When evaluating at realized prices, deltas must be the same length as the data")
    end
end

bernO = convert.(Integer, bO);
order = bernO;

sig = [];
for j = 1:1:J
    push!(sig, inv_sigma[j])
end
if !isempty(maxes)
    maxs = maxes[1];
    mins = maxes[2];
end

## Declare prices and shares to evaluate derivatives
svec = similar(deltas)
if trueS == false
    for p_i = 1:size(deltas,1)
        s!(sj) = solve_s_nested_flexible(sj, inv_sigma, deltas[p_i,:]', J, bernO, included, maxes,nothing);
        ans = nlsolve(s!, 1/(2*J) .* ones(J))
        svec[p_i,:] =  ans.zero;
    end
else
    svec = s;
end

# Share Jacobian
dsids = zeros(J,J,size(deltas,1)) # initialize matrix of \partial s^{-1} /\partial s
for j1 = 1:J
    if j1>1
        if j1<J
            perm_order = hcat(convert.(Integer,j1:J)', convert.(Integer, 1:j1-1)');
            perm_order = convert.(Integer, perm_order)
        else
            perm_order = convert.(Integer,1:J)
            perm_order = perm_order .- 1
            perm_order[1] = J
        end
    else
        perm_order = convert.(Integer, 1:J)
    end
    s_perm = similar(svec)
    for i = 1:J
        s_perm[:,i] = svec[:,perm_order[i]];
    end
    for j2 = 1:J
        if included[j1,j2]==1 # cross-derivatives are zero if not included in j1 demand
            tempmat_s = zeros(size(deltas,1),1)
            for j3 = 1:1:J
                stemp = s_perm[:,j3];
                if included[j1,perm_order[j3]]==1 # design matrix only uses included products
                    if j2==perm_order[j3]
                        tempmat_s = [tempmat_s dbern(stemp, bernO[j1,1])];
                    else
                        tempmat_s = [tempmat_s bern(stemp, bernO[j1,1])];
                    end
                end
            end

            # Market characteristics
            marketvars = nothing
            if marketvars!=nothing
                for i = 1:size(marketvars,2)
                    tempmat_s = [tempmat_s bern(marketvars[:,i], bernO[j1,1])];
                end
            end
            tempmat_s = tempmat_s[:,2:end]
            tempmat_s = fullInteraction(tempmat_s, convert.(Integer, sum(included[j1,:])), convert.(Integer, bernO[j1,1]));
            if !isempty(maxes)
                dsids[j1,j2,:] = tempmat_s*inv_sigma[j1] ./ (maxs[j2] - mins[j2])  ;
            else
                dsids[j1,j2,:] = tempmat_s*inv_sigma[j1];
            end
        else
            zerotemp = zeros(size(svec,1),1);
            dsids[j1,j2,:] = zerotemp;
        end
    end
end

Jmat = [];
J_sp = zeros(size(svec[:,1]));
all_own = zeros(size(svec,1),J);
if !isempty(maxes)
    svec2 = svec.*(maxs - mins) + mins;
else
    svec2 = svec;
end
for ii = 1:length(dsids[1,1,:])
    J_s = [];
    for j1 = 1:J
        J_sj1 = [];
        for j2 = 1:J
            if j2==1
                J_sj1 = dsids[j1,j2,ii]
            else
                J_sj1 = [J_sj1 dsids[j1,j2,ii]]
            end
        end
        if j1==1
            J_s = J_sj1;
        else
            J_s = vcat(J_s, J_sj1);
        end
    end
    temp = -1*inv(J_s);
    push!(Jmat, temp)
    ps = p_points[ii,:]./svec2[ii,:]
    all_own[ii,:] = -1*inv(J_s)*ps;
    J_sp[ii,1] = temp[whichProducts[1],whichProducts[2]];
    #print("Market $ii \n")
end

esep = J_sp.*p_points[:,1]./svec2[:,1]; # own-price varying
numBadMarkets = 0
# if trueS==1
#     print("There were $numBadMarkets bad markets")
# end
return esep, Jmat, svec, all_own
end


text:
x0 = [1., 0.]
P0 = Matrix(1.0I, 2, 2)

Φ = [0.8 0.5; -0.1 0.8]
b = zeros(2)
Q = [0.2 0.0; 0.0 1.0]

yshadow = [0.0]
H = [1.0 0.0]
R = Matrix(1.0I, 1, 1)

E = LinearEvolution(Φ, Gaussian(b, Q))
Obs = LinearObservationModel(H, R)
M = LinearStateSpaceModel(E, Obs)


text:
s = open("input.txt") do file
    readlines(file);
end

function EvalExpression(expr)
    # println("=====")
    # println(expr);
    newExp = expr;
    m = match(r"(\d+) \+ (\d+)", expr)
    while m != nothing
        val = parse(Int, m.captures[1]) + parse(Int, m.captures[2]);
        newExp = replace(newExp, string(m.match)=>"$(val)", count=1);

        m = match(r"(\d+) \+ (\d+)", newExp)
    end

    # println(newExp)

    nums = [];
    for m in eachmatch(r"\d+", newExp)
        push!(nums, parse(Int, m.match));
    end
    # println(nums);
    return *(nums...);
end

function DoTheMath(expression)
    println(expression);
    expressionStack = [];
    currentExpr = "";
    for i = 1:length(expression)
        if expression[i] == '('
            push!(expressionStack, currentExpr);
            currentExpr = "";
        elseif expression[i] == ')'
            val = EvalExpression(currentExpr);
            currentExpr = pop!(expressionStack) * "$(val)";
        else
            currentExpr *= string(expression[i])
        end
    end
    # println(currentExpr);

    val = EvalExpression(currentExpr);
    # println(val)
    # println("----------")
    return Int128(val);
end

eachAns = map(DoTheMath, s);

answer = sum(eachAns)


text:
@testset "reservoir simple" begin
    res = Wflow.SimpleReservoir{Float64}(
        Δt = 86400.0,
        demand = [52.523],
        maxrelease = [420.184],
        maxvolume = [25_000_000.0],
        volume = [1.925e7],
        totaloutflow = [0.0],
        inflow = [0.0],
        area = [1885665.353626924],
        targetfullfrac = [0.8],
        targetminfrac = [0.2425554726620697],
        precipitation = [4.2],
        evaporation = [1.5],
        outflow = [NaN],
        percfull = [NaN],
        demandrelease = [NaN],
    )

    Wflow.update(res, 1, 100.0, 86400.0)
    @test res.outflow[1] ≈ 91.3783714867453
    @test res.totaloutflow[1] ≈ 7.895091296454794e6
    @test res.volume[1] ≈ 2.0e7
    @test res.percfull[1] ≈ 0.80
    @test res.demandrelease[1] ≈ 52.5229994727611
    @test res.precipitation[1] ≈ 4.2
    @test res.evaporation[1] ≈ 1.5
end

@testset "natural lake" begin
    lake = Wflow.NaturalLake{Float64}(
        Δt = 86400.0,
        lowerlake_ind = [0],
        area = [180510409.0],
        threshold = [0.0],
        storfunc = [1],
        outflowfunc = [3],
        totaloutflow = [0.0],
        inflow = [0.0],
        b = [0.22],
        e = [2.0],
        sh = [missing],
        hq = [missing],
        storage = Wflow.initialize_storage([1], [180510409.0], [18.5], [missing]),
        waterlevel = [18.5],
        precipitation = [20.0],
        evaporation = [3.2],
        outflow = [NaN],
    )

    Wflow.update(lake, 1, 2500.0, 181, 86400.0)
    @test lake.outflow[1] ≈ 85.31903276150577
    @test lake.totaloutflow[1] ≈ 7.371564430594098e6
    @test lake.storage[1] ≈ 3.551103576940606e9
    @test lake.waterlevel[1] ≈ 19.672569557695734
    @test lake.precipitation[1] ≈ 20.0
    @test lake.evaporation[1] ≈ 3.2
end

datadir = joinpath(@__DIR__, "data")
sh = [
    Wflow.read_sh_csv(joinpath(datadir, "lake_sh_1.csv")),
    Wflow.read_sh_csv(joinpath(datadir, "lake_sh_2.csv")),
]
@testset "linked lakes (HBV)" begin
    @test keys(sh[1]) == (:H, :S)
    @test typeof(values(sh[1])) == Tuple{Vector{Float},Vector{Float}}

    lake = Wflow.NaturalLake{Float}(
        Δt = 86400.0,
        lowerlake_ind = [2, 0],
        area = [472461536.0, 60851088.0],
        threshold = [393.7, 0.0],
        storfunc = [2, 2],
        inflow = [0.0, 0.0],
        totaloutflow = [0.0, 0.0],
        outflowfunc = [2, 1],
        b = [140.0, 0.0],
        e = [1.5, 1.5],
        sh = sh,
        hq = [missing, Wflow.read_hq_csv(joinpath(datadir, "lake_hq_2.csv"))],
        waterlevel = [395.03027, 394.87833],
        precipitation = [10.0, 10.0],
        evaporation = [2.0, 2.0],
        outflow = [NaN, NaN],
        storage = Wflow.initialize_storage(
            [2, 2],
            [472461536.0, 60851088.0],
            [395.03027, 394.87833],
            sh,
        ),
    )

    Wflow.update(lake, 1, 500.0, 15, 86400.0)
    Wflow.update(lake, 2, 500.0, 15, 86400.0)
    @test lake.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2
    @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3
    @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4
    @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2
    lake.totaloutflow .= 0.0
    lake.inflow .= 0.0
    Wflow.update(lake, 1, 500.0, 15, 86400.0)
    Wflow.update(lake, 2, 500.0, 15, 86400.0)
    @test lake.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2
    @test lake.totaloutflow ≈ [-2.2446764487487033e7, 2.070723775102806e7] atol = 1e3
    @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4
    @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2

end


text:
#==========================================================================================#

# TYPE

mutable struct IV <: GMM

    method::String
    sample::Microdata
    β::Vector{Float64}
    V::Matrix{Float64}
    W::Matrix{Float64}

    IV() = new()
end

#==========================================================================================#

# CONSTRUCTOR

function IV(MD::Microdata, method::String)
    obj        = IV()
    obj.sample = MD
    obj.method = method
    return obj
end

#==========================================================================================#

# INTERFACE

function fit(::Type{IV}, MD::Microdata; novar::Bool = false, method::String = "TSLS")

    if method == "OLS"

        FSD                   = Microdata(MD)
        FSD.mapping[:control] = vcat(MD.mapping[:treatment], MD.mapping[:control])

        pop!(FSD.mapping, :treatment)
        pop!(FSD.mapping, :instrument)

        obj = OLS(FSD)

        _fit!(obj, getweights(obj))

    elseif method == "Reduced form"

        FSD                   = Microdata(MD)
        FSD.mapping[:control] = vcat(MD.mapping[:instrument], MD.mapping[:control])

        pop!(FSD.mapping, :treatment)
        pop!(FSD.mapping, :instrument)

        obj = OLS(FSD)

        _fit!(obj, getweights(obj))

    elseif method == "First stage"

        FSD                    = Microdata(MD)
        FSD.mapping[:control]  = vcat(MD.mapping[:instrument], MD.mapping[:control])
        FSD.mapping[:response] = MD.mapping[:treatment]

        pop!(FSD.mapping, :treatment)
        pop!(FSD.mapping, :instrument)

        obj = OLS(FSD)

        _fit!(obj, getweights(obj))

    elseif length(MD.mapping[:treatment]) == length(MD.mapping[:instrument])

        obj   = IV(MD, "Method of moments")
        k     = length(MD.mapping[:instrument]) + length(MD.mapping[:control])
        obj.W = Matrix{Float64}(I, k, k)

        _fit!(obj, getweights(obj))

    elseif (method == "TSLS") | (method == "2SLS")

        obj   = IV(MD, "Two-step GMM")
        obj.W = crossprod(getmatrix(obj, :instrument, :control), getweights(obj))

        _fit!(obj, getweights(obj))

    elseif (method == "Two-step GMM") | (method == "Optimal GMM")

        obj   = IV(MD, method) ; _fit!(obj, getweights(obj))
        obj.W = wmatrix(obj, getcorr(obj), getweights(obj))

        _fit!(obj, obj.W, getweights(obj))

    else
        throw("unknown method")
    end

    novar || _vcov!(obj, getcorr(obj), getweights(obj))

    return obj
end

#==========================================================================================#

# ESTIMATION

function _fit!(obj::IV, ::UnitWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)

    if obj.method == "Method of moments"
        obj.β = (z' * x) \ (z' * y)
    else
        zγ    = z * (z \ x)
        obj.β = (zγ' * x) \ (zγ' * y)
    end
end

function _fit!(obj::IV, W::Matrix{Float64}, ::UnitWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)

    zγ    = z * (W \ (z' * x))
    obj.β = (zγ' * x) \ (zγ' * y)
end

function _fit!(obj::IV, w::AbstractWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)
    v = Diagonal(w) * z

    if obj.method == "Method of moments"
        obj.β = (v' * x) \ (v' * y)
    else
        γ     = (v' * z) \ (v' * x)
        vγ    = v * γ
        obj.β = (vγ' * x) \ (vγ' * y)
    end
end

function _fit!(obj::IV, W::Matrix{Float64}, w::AbstractWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)
    v = Diagonal(w) * z

    vγ    = v * (W \ (v' * x))
    obj.β = (vγ' * x) \ (vγ' * y)
end

#==========================================================================================#

# SCORE (MOMENT CONDITIONS)

score(obj::IV) = Diagonal(residuals(obj)) * getmatrix(obj, :instrument, :control)

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

function jacobian(obj::IV, ::UnitWeights)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)
    return - z' * x
end

function jacobian(obj::IV, w::AbstractWeights)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)
    v = Diagonal(w) * z
    return - v' * x
end

# VARIANCE MATRIX

function _vcov!(obj::IV, corr::Homoscedastic, ::UnitWeights)

    x  = getmatrix(obj, :treatment, :control)
    z  = getmatrix(obj, :instrument, :control)
    r  = residuals(obj)
    σ² = sum(abs2, r) / dof_residual(obj)
    γ  = z \ x
    V  = x' * z * γ

    obj.V = lmul!(σ², inv(V))
end

function _vcov!(obj::IV, corr::Homoscedastic, w::AbstractWeights)

    x  = getmatrix(obj, :treatment, :control)
    z  = getmatrix(obj, :instrument, :control)
    v  = Diagonal(w) * z
    r  = residuals(obj)
    σ² = sum(abs2.(r), w) / dof_residual(obj)
    γ  = (v' * z) \ (v' * x)
    V  = x' * v * γ

    obj.V = lmul!(σ², inv(V))
end

#==========================================================================================#

# LINEAR PREDICTOR

function predict(obj::IV, MD::Microdata)
    if getnames(obj, :treatment, :control) != getnames(MD, :treatment, :control)
        throw("missing variables")
    end
    getmatrix(MD, :treatment, :control) * obj.β
end

# FITTED VALUES

fitted(obj::IV, MD::Microdata) = predict(obj, MD)

# DERIVATIVE OF FITTED VALUES

jacobexp(obj::IV) = getmatrix(obj, :treatment, :control)

#==========================================================================================#

# UTILITIES

coefnames(obj::IV) = getnames(obj, :treatment, :control)
mtitle(obj::IV)    =  "Linear GMM"

# COEFFICIENT OF DETERMINATION

adjr2(obj::IV) = 1.0 - (1.0 - r2(obj)) * (nobs(obj) - 1) / dof_residual(obj)
r2(obj::IV)    = _r2(obj, getweights(obj))

function _r2(obj::IV, ::UnitWeights)
    y   = response(obj)
    ŷ   = fitted(obj)
    rss = sum(abs2, y .- ŷ)
    tss = sum(abs2, y .- mean(y))
    return 1.0 - rss / tss
end

function _r2(obj::IV, w::AbstractWeights)
    y   = response(obj)
    μ   = mean(y, w)
    ŷ   = fitted(obj)
    rss = sum(abs2.(y .- ŷ), w)
    tss = sum(abs2.(y .- μ), w)
    return 1.0 - rss / tss
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, Example

# sets default format or read provided argument
if length(ARGS)==0
    outputFormat = "html"
else
    outputFormat = ARGS[1]
end

# sets different format and pages settings for each output format
if outputFormat == "pdf"

    makedocs(
      sitename="Documentation of Example Project",
      format = Documenter.LaTeX(),
      pages = [
        "index.md",
        "more.md",
        "thirdA.md",
        "thirdB.md"
      ],
      authors = "AuthorName"
    )

elseif outputFormat == "html"

  makedocs(
    sitename="Documentation of Example Project",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
      "Intro side title" => "index.md",
      "More side title" => "more.md",
      "Third A side title" => "thirdA.md",
      "Third B side title" => "thirdB.md"
    ]
  )

else
  error("Output format not considered yet. Please open an issue.")
end

if outputFormat == "html"
  deploydocs(
    repo = "github.com/jorgepz/Julia_Docs_Testbed.git",
    push_preview=true,
    devbranch = "main"
  )
end


text:
function fsort!(vec::CategoricalArray)
    cnts = zeros(UInt, length(vec.pool))

    @inbounds for r in vec.refs
        cnts[r] += 1
    end

    # vec.refs .= StatsBase.inverse_rle(1:length(vec.pool), cnts) # the below is 5x faster
    j = 0
    @inbounds for ref in 1:length(vec.pool)
        for k = 1:cnts[ref]
            j += 1
            vec.refs[j] = ref
        end
    end

    vec
end

fsort(vec::CategoricalArray) = fsort!(copy(vec))

function fsort2!(vec::CategoricalArray)
    Base.Sort.sort_int_range!(vec.refs, length(vec.pool), 0)
    vec
end

fsort2(vec::CategoricalArray) = fsort2!(copy(vec))


text:
using IndexFunArrays, LinearAlgebra, NDTools
export damp_edge_outside

"""
    damp_edge_outside(img::AbstractArray{T,N}, border, mykernel, usepixels)

Extrapolates the data by filling in blurred information outside the edges. This is a bit like DampEdge but using a normalized convolution with a kernel

# Arguments:
+ `img`: image to extrapolate from
+ `border`: percentage of border pixels to add
+ `mykernel`: A kernel can be provided. Default: 1/r^3 kernel of support `norm(border size) * sqrt(2)`
+ `usepixels`: number of border pixels in the image to use on each edge. Default: 2. Zero means use all pixels

# Example:
```julia
julia> using TestImages, FourierTools, Colors

julia> img = permutedims(Float32.(testimage("fabio_gray_512.png")),[2 1]);

julia> img_d = damp_edge_outside(img);

julia> Gray.(log.(1 .+abs.(ft(img)))) ./ 10

julia> Gray.(log.(1 .+abs.(ft(img_d)))) ./ 10
```
"""
function damp_edge_outside(img::AbstractArray{T, N}, border=0.1, mykernel=nothing, usepixels=2) where {T,N}
    kernelpower=3;
    rborder=ceil.(Int,border.*size(img));
    new_size=size(img).+ rborder;
    if isnothing(mykernel)
        mykernel = (one(T)./rr(new_size)) .- one(T)./norm(rborder.*sqrt(2.0));
        # clip at zero
        mykernel .= max.(mykernel, zero(T));
        mykernel=mykernel .^ kernelpower;
        midpos = size(mykernel) .÷2 .+ 1;
        # This does not matter as it only applies to points directly having data.
        mykernel[midpos...] = zero(T);
    end
    transfer=ft(mykernel);
    # weight image
    wimg=ones(eltype(img),size(img));
    nimg=copy(img);
    if usepixels > 0
        inner_nimg = NDTools.select_region_view(nimg, new_size = size(img) .- 2 .*usepixels)
        inner_wimg = NDTools.select_region_view(wimg, new_size = size(img) .- 2 .*usepixels)
        inner_nimg .= zero(T)
        inner_wimg .= zero(T)
    end
    nimg = NDTools.select_region(nimg,new_size=new_size);
    # just mark every pixel
    wimg = NDTools.select_region(wimg,new_size=new_size);

    # should this only work for 2D?
    nimg2 = real(ift(ft(nimg) .* transfer));
    wimg2 = real(ift(ft(wimg) .* transfer));

    nimg = nimg2 ./ wimg2;
    if usepixels > 0
        # replace the original in the middle
        NDTools.select_region!(img, nimg)
    end
    return nimg;

    # ToDo: Phase subpixel peak determination
    # ToDo: DampEdge (with Gaussian filter)
    # ToDo: make it work in ND
end


text:
"""
Features extraction from returns proposed by Betina and al.
"""
#KLT stands for the long-term period
#KST stands for the short-term period
#ϱ is the number of standard deviations above the means
function feature_indexes_sgns(asset_rtn_sgn, kst, klt, kmom, ϱ=2)
    #current time
    t = size(asset_rtn_sgn, 1)
    # sig_1 and sig_2 are
    # short-term and long-term simple moving averages
    sig_1 = (1 / kst) * sum(asset_rtn_sgn[t - d] for d in 0:(kst - 1))
    sig_2 = (1 / klt) * sum(asset_rtn_sgn[t - d] for d in 0:(klt - 1))

    # sig_3 and sig_4 are exponential moving average (short and long)
    sig_3 =
        (1.0 / sum(exp(-d / kst) for d in 1:kst)) *
        sum(exp(-d / kst) * asset_rtn_sgn[t - d] for d in 0:(kst - 1))
    sig_4 =
        (1.0 / sum(exp(-d / klt) for d in 1:klt)) *
        sum(exp(-d / klt) * asset_rtn_sgn[t - d] for d in 0:(klt - 1))
    # RSI
    R_pos = 0
    R_neg = 0
    num_pos = 0
    num_neg = 0
    for i in (t - kmom + 1):t
        if (asset_rtn_sgn[i] > 0)
            R_pos = R_pos + asset_rtn_sgn[i]
            num_pos = num_pos + 1
        else
            R_neg = R_neg + asset_rtn_sgn[i]
            num_neg = num_neg + 1
        end
    end
    med_pos = R_pos / num_pos
    if num_neg == 0
        med_neg = 1
    else
        med_neg = abs(R_neg / num_neg)
    end
    sig_5 = med_pos / med_neg
    # estimated shor-t and long-t standard deviations
    σ_st_6 = sqrt((1.0 / kst) * sum((asset_rtn_sgn[t - d] - sig_1)^2 for d in 0:(kst - 1)))
    σ_lt_7 = sqrt((1.0 / klt) * sum((asset_rtn_sgn[t - d] - sig_2)^2 for d in 0:(klt - 1)))
    # sig_6 and sig_7 are

    if sig_1 <= 0
        sig_6 = sig_1 - ϱ * σ_st_6
    else
        sig_6 = sig_1 + ϱ * σ_st_6
    end
    if sig_2 <= 0
        sig_7 = sig_2 - ϱ * σ_lt_7
    else
        sig_7 = sig_2 + ϱ * σ_lt_7
    end
    σ_st_8 = sqrt((1.0 / kst) * sum((asset_rtn_sgn[t - d] - sig_3)^2 for d in 0:(kst - 1)))
    σ_lt_9 = sqrt((1.0 / klt) * sum((asset_rtn_sgn[t - d] - sig_4)^2 for d in 0:(klt - 1)))
    if sig_3 <= 0
        sig_8 = sig_3 - ϱ * σ_st_8
    else
        sig_8 = sig_3 + ϱ * σ_st_8
    end
    if sig_4 <= 0
        sig_9 = sig_4 - ϱ * σ_lt_9
    else
        sig_9 = sig_4 + ϱ * σ_lt_9
    end

    return [sig_1; sig_2; sig_3; sig_4; sig_5; sig_6; sig_7; sig_8; sig_9]
end

"""
Mixed signals predictor proposed by Betina and al. Univariate
"""
function mixed_signals_predict_return(asset_rtn_sgn, num_t, kst_a, klt_a, kmom, ϱ=2)
    numD = size(asset_rtn_sgn, 1)
    num_train = min(numD - 1, num_t)

    ## optimize features weights (least squares based on past returns)
    # signal features extraction train (from klt+1 to klt+num_train)
    #kst = max(min(numD-num_train-10,kst_a),0)
    #klt = max(min(numD-num_train-10,klt_a),0)
    sig_features = zeros(num_train, 9) #each line is the i-element's featurs
    for i in 1:num_train
        #kst = max(min(numD-i-2,kst_a),0)
        klt = max(min(numD - i - 2, klt_a), 0)
        kst = Int64(floor(klt / 2.0))
        sig_features[num_train + 1 - i, :] = feature_indexes_sgns(
            asset_rtn_sgn[1:(numD - i)], kst, klt, kmom
        )
    end
    # optimize weights (least squares - pseudo inverse)
    try
        weights_opt = \(
            sig_features'sig_features,
            sig_features'asset_rtn_sgn[(numD - num_train + 1):numD],
        )
    catch
        weights_opt = pinv(sig_features, 1E-07) * asset_rtn_sgn[(numD - num_train + 1):numD]
    end
    #  normalize features weights
    weights_opt = weights_opt ./ sum(weights_opt)
    ## extract features current return
    kst = min(numD, kst_a)
    klt = min(numD, klt_a)
    features = feature_indexes_sgns(asset_rtn_sgn[1:numD], kst, klt, kmom)
    ## extrapolate return

    return (features'weights_opt)[1]
end


text:
module TestUtils

export random_complex_matrix,
    random_real_matrix,
    random_hermitian_matrix,
    random_complex_sparse_matrix,
    random_real_sparse_matrix,
    random_hermitian_sparse_matrix,
    random_state_vector,
    random_hermitian_real_matrix,
    random_hermitian_sparse_real_matrix
export dummy_control_problem
export test

using Logging
using Random
using Distributions
using LinearAlgebra
using SparseArrays
using Coverage
using LocalCoverage
using Printf

import ..Objective
import ..ControlProblem
import ..getcontrols
import ..discretize
import ..discretize_on_midpoints

"""Run a package test-suite in a subprocess.

```julia
test(
    file="test/runtests.jl";
    root=pwd(),
    project="test",
    code_coverage="user",
    show_coverage=(code_coverage == "user"),
    color=<inherit>,
    compiled_modules=<inherit>,
    startup_file=<inherit>,
    depwarn=<inherit>,
    inline=<inherit>,
    check_bounds="yes",
    track_allocation=<inherit>,
    threads=<inherit>,
    genhtml=false,
    covdir="coverage"
)
```

runs the test suite of the package located at `root` by running `include(file)`
inside a new julia process.

This is similar to what `Pkg.test()` does, but differs in the "sandboxing"
approach. While `Pkg.test()` creates a new temporary sandboxed environment,
`test()` uses an existing environment in `project` (the `test` subfolder by
default). This allows testing against the dev-versions of other packages. It
requires that the `test` folder contains both a `Project.toml` and a
`Manifest.toml` file.

The `test()` function also differs from directly including `test/runtests.jl`
in the REPL in that it can generate coverage data and reports (this is only
possible when running tests in a subprocess).

If `show_coverage` is passed as `true` (default), a coverage summary is shown.
Further, if `genhtml` is `true`, a full HTML coverage report will be generated
in `covdir` (relative to `root`). This requires the `genhtml` executable (part
of the [lcov](http://ltp.sourceforge.net/coverage/lcov.php) package). Instead
of `true`, it is also possible to pass the path to the `genhtml` exectuable.

All other keyword arguments correspond to the respective command line flag for
the `julia` executable that is run as the subprocess.

This function is intended to be exposed in a project's development-REPL.
"""
function test(
    file="test/runtests.jl";
    root=pwd(),
    project="test",
    code_coverage="user",
    show_coverage=(code_coverage == "user"),
    color=(Base.have_color === nothing ? "auto" : Base.have_color ? "yes" : "no"),
    compiled_modules=(Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no"),
    startup_file=(Base.JLOptions().startupfile == 1 ? "yes" : "no"),
    depwarn=(Base.JLOptions().depwarn == 2 ? "error" : "yes"),
    inline=(Bool(Base.JLOptions().can_inline) ? "yes" : "no"),
    track_allocation=(("none", "user", "all")[Base.JLOptions().malloc_log+1]),
    check_bounds="yes",
    threads=Threads.nthreads(),
    genhtml::Union{Bool,AbstractString}=false,
    covdir="coverage"
)
    julia = Base.julia_cmd().exec[1]
    cmd = [
        julia,
        "--project=$project",
        "--color=$color",
        "--compiled-modules=$compiled_modules",
        "--startup-file=$startup_file",
        "--code-coverage=$code_coverage",
        "--track-allocation=$track_allocation",
        "--depwarn=$depwarn",
        "--check-bounds=$check_bounds",
        "--threads=$threads",
        "--inline=$inline",
        "--eval",
        "include(\"$file\")"
    ]
    @info "Running '$(join(cmd, " "))' in subprocess"
    run(Cmd(Cmd(cmd), dir=root))
    tracefile = joinpath(root, "lcov.info")
    if show_coverage || genhtml
        logger = Logging.SimpleLogger(Logging.Error)
        local coverage
        Logging.with_logger(logger) do
            coverage = Coverage.process_folder(joinpath(root, "src"))
        end
        if show_coverage
            coverage_summary(coverage)
        end
        (genhtml === true) && (genhtml = "genhtml")
        (genhtml === false) && (genhtml = "")
        if !isempty(genhtml)
            covdir = joinpath(root, covdir)
            LocalCoverage.CoverageTools.LCOV.writefile(tracefile, coverage)
            branch = try
                strip(read(`git rev-parse --abbrev-ref HEAD`, String))
            catch e
                @warn "git branch could not be detected.\nError message: $(sprint(Base.showerror, e))"
            end
            title = isnothing(branch) ? "N/A" : "on branch $(branch)"
            try
                run(`$(genhtml) -t $(title) -o $(covdir) $(tracefile)`)
            catch e
                @error(
                    "Failed to run $(genhtml). Check that lcov is installed.\nError message: $(sprint(Base.showerror, e))"
                )
            end
            @info(
                "Generated coverage HTML. Serve with 'LiveServer.serve(dir=\"$(relpath(covdir, pwd()))\")'"
            )
        end
    end
end

"""Construct a random complex matrix of size N×N with spectral radius ρ.

```julia
random_complex_matrix(N, ρ)
```
"""
function random_complex_matrix(N, ρ)
    σ = 1 / √N
    d = Normal(0.0, σ)
    H = ρ * (rand(d, (N, N)) + rand(d, (N, N)) * 1im) / √2
end

"""Construct a random real-valued matrix of size N×N with spectral radius ρ.

```julia
random_real_matrix(N, ρ)
```
"""
function random_real_matrix(N, ρ)
    σ = 1 / √N
    d = Normal(0.0, σ)
    H = ρ * rand(d, (N, N))
end

"""Construct a random Hermitian matrix of size N×N with spectral radius ρ.

```julia
random_hermitian_matrix(N, ρ)
```
"""
function random_hermitian_matrix(N, ρ)
    σ = 1 / √N
    d = Normal(0.0, σ)
    X = (rand(d, (N, N)) + rand(d, (N, N)) * 1im) / √2
    H = ρ * (X + X') / (2 * √2)
end

"""Construct a random Hermitian real matrix of size N×N with spectral radius ρ.

```julia
random_hermitian_real_matrix(N, ρ)
```
"""
function random_hermitian_real_matrix(N, ρ)
    σ = 1 / √N
    d = Normal(0.0, σ)
    X = rand(d, (N, N))
    H = ρ * (X + X') / (2 * √2)
end

"""Construct a random sparse complex matrix.

```julia
random_complex_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_complex_sparse_matrix(N, ρ, sparsity)
    σ = 1 / √(sparsity * N)
    d = Normal(0.0, σ)
    Hre = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    Him = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    H = ρ * (Hre + Him * 1im) / √2
end

"""Construct a random sparse real-valued matrix.

```julia
random_real_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_real_sparse_matrix(N, ρ, sparsity)
    σ = 1 / √(sparsity * N)
    d = Normal(0.0, σ)
    H = ρ * sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
end

"""Construct a random sparse Hermitian matrix.

```julia
random_hermitian_sparse_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_hermitian_sparse_matrix(N, ρ, sparsity)
    σ = 1 / √(sparsity * N)
    d = Normal(0.0, σ)
    H1 = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    H2 = copy(H1)
    H2.nzval .= rand(d, length(H2.nzval))
    X = (H1 + H2 * 1im) / √2
    return 0.5ρ * (X + X') / √2
end

"""Construct a random sparse Hermitian real matrix.

```julia
random_hermitian_sparse_real_matrix(N, ρ, sparsity)
```

returns a matrix of size N×N with spectral radius ρ and the given sparsity
(number between zero and one that is the approximate fraction of non-zero
elements).
"""
function random_hermitian_sparse_real_matrix(N, ρ, sparsity)
    σ = 1 / √(sparsity * N)
    d = Normal(0.0, σ)
    H = sprand(N, N, sparsity, (dims...) -> rand(d, dims...))
    return 0.5ρ * (H + H') / √2
end

"""Return a random, normalized Hilbert space state vector of dimension `N`.

```julia
random_state_vector(N)
```
"""
function random_state_vector(N)
    Ψ = rand(N) .* exp.((2π * im) .* rand(N))
    Ψ ./= norm(Ψ)
    return Ψ
end

"""Set up a dummy control problem.

```julia
problem = dummy_control_problem(;
    N=10, n_objectives=1, n_controls=1, n_steps=50, dt=1.0, sparsity=0.5,
    complex_operators=true, hermitian=true, kwargs...)
```

Sets up a control problem with random (sparse) Hermitian matrices.

# Arguments

* `N`: The dimension of the Hilbert space
* `n_objectives`: The number of objectives in the optimization. All objectives
  will have the same Hamiltonian, but random initial and target states.
* `n_controls`: The number of controls, that is, the number of control terms in
  the control Hamiltonian. Each control is an array of random values,
  normalized on the intervals of the time grid.
* `n_steps`: The number of time steps (intervals of the time grid)
* `dt`: The time step
* `sparsity`: The sparsity of the Hamiltonians, as a number between 0.0 and
  1.0. For `sparsity=1.0`, the Hamiltonians will be dense matrices.
* `complex_operators`: Whether or not the drift/control operators will be
  complex-valued or real-valued.
* `hermitian`: Whether or not all drift/control operators will be Hermitian matrices.
* `kwargs`: All other keyword arguments are passed on to
  [`ControlProblem`](@ref)
"""
function dummy_control_problem(;
    N=10,
    n_objectives=1,
    n_controls=1,
    n_steps=50,
    dt=1.0,
    sparsity=0.5,
    complex_operators=true,
    hermitian=true,
    kwargs...
)

    tlist = collect(range(0; length=(n_steps + 1), step=dt))
    pulses = [rand(length(tlist) - 1) for l = 1:n_controls]
    for l = 1:n_controls
        # we normalize on the *intervals*, not on the time grid points
        pulses[l] ./= norm(pulses[l])
    end
    controls = [discretize(pulse, tlist) for pulse in pulses]

    function random_op(N, ρ, sparsity, complex_operators, hermitian)
        if sparsity < 1.0
            if hermitian
                if complex_operators
                    H = random_hermitian_sparse_matrix(N, ρ, sparsity)
                end
            else
                if complex_operators
                    H = random_complex_sparse_matrix(N, ρ, sparsity)
                end
            end
        else
            if hermitian
                if complex_operators
                    H = random_hermitian_matrix(N, ρ)
                else
                    H = random_hermitian_real_matrix(N, ρ)
                end
            else
                if complex_operators
                    H = random_complex_matrix(N, ρ)
                else
                    H = random_real_matrix(N, ρ)
                end
            end
        end
    end

    hamiltonian = []
    H_0 = random_op(N, 1.0, sparsity, complex_operators, hermitian)
    push!(hamiltonian, H_0)
    for control ∈ controls
        H_c = random_op(N, 1.0, sparsity, complex_operators, hermitian)
        push!(hamiltonian, (H_c, control))
    end

    objectives = [
        Objective(;
            initial_state=random_state_vector(N),
            generator=tuple(hamiltonian...),
            target_state=random_state_vector(N)
        ) for k = 1:n_objectives
    ]

    return ControlProblem(
        objectives=objectives,
        pulse_options=Dict(
            control => Dict(:lambda_a => 1.0, :update_shape => t -> 1.0) for
            control in controls
        ),
        tlist=tlist,
        kwargs...
    )
end

"""Result returned by [`optimize_with_dummy_method`](@ref)."""
mutable struct DummyOptimizationResult
    tlist::Vector{Float64}
    iter_start::Int64  # the starting iteration number
    iter_stop::Int64 # the maximum iteration number
    iter::Int64  # the current iteration number
    J_T::Float64  # the current value of the final-time functional J_T
    J_T_prev::Float64  # previous value of J_T
    guess_controls::Vector{Vector{Float64}}
    optimized_controls::Vector{Vector{Float64}}
    converged::Bool
    message::String

    function DummyOptimizationResult(problem)
        tlist = problem.tlist
        controls = getcontrols(problem.objectives)
        iter_start = get(problem.kwargs, :iter_start, 0)
        iter = iter_start
        iter_stop = get(problem.kwargs, :iter_stop, 20)
        guess_controls = [discretize(control, tlist) for control in controls]
        J_T = 0.0
        J_T_prev = 0.0
        optimized_controls = [copy(guess) for guess in guess_controls]
        converged = false
        message = "in progress"
        new(
            tlist,
            iter_start,
            iter_stop,
            iter,
            J_T,
            J_T_prev,
            guess_controls,
            optimized_controls,
            converged,
            message
        )
    end

end

struct DummyOptimizationWrk
    objectives
    adjoint_objectives
    kwargs
    controls
    pulses0::Vector{Vector{Float64}}
    pulses1::Vector{Vector{Float64}}
    result
end

function DummyOptimizationWrk(problem)
    objectives = [obj for obj in problem.objectives]
    adjoint_objectives = [adjoint(obj) for obj in problem.objectives]
    controls = getcontrols(objectives)
    kwargs = Dict(problem.kwargs)
    tlist = problem.tlist
    if haskey(kwargs, :continue_from)
        @info "Continuing previous optimization"
        result = kwargs[:continue_from]
        if !(result isa DummyOptimizationResult)
            result = convert(DummyOptimizationResult, result)
        end
        result.iter_stop = get(problem.kwargs, :iter_stop, 20)
        result.converged = false
        result.message = "in progress"
        pulses0 = [
            discretize_on_midpoints(control, tlist) for control in result.optimized_controls
        ]
    else
        result = DummyOptimizationResult(problem)
        pulses0 = [discretize_on_midpoints(control, tlist) for control in controls]
    end
    pulses1 = [copy(pulse) for pulse in pulses0]
    return DummyOptimizationWrk(
        objectives,
        adjoint_objectives,
        kwargs,
        controls,
        pulses0,
        pulses1,
        result,
    )
end

function update_result!(wrk::DummyOptimizationWrk, ϵ⁽ⁱ⁺¹⁾, i::Int64)
    res = wrk.result
    res.J_T_prev = res.J_T
    res.J_T = sum([norm(ϵ) for ϵ ∈ ϵ⁽ⁱ⁺¹⁾])
    (i > 0) && (res.iter = i)
    if i >= res.iter_stop
        res.converged = true
        res.message = "Reached maximum number of iterations"
    end
end

"""Run a dummy optimization.

```julia
result = optimize(problem, method=:dummymethod)
```

runs through and "optimization" of the given `problem` where in each iteration,
the amplitude of the guess pulses is diminished by 10%. The (summed) vector
norm of the the control serves as the value of the optimization functional.
"""
function optimize_with_dummy_method(problem)
    # This is connected to the main `optimize` method in the main
    # QuantumcontrolBase.jl
    iter_start = get(problem.kwargs, :iter_start, 0)
    check_convergence! = get(problem.kwargs, :check_convergence, res -> res)
    wrk = DummyOptimizationWrk(problem)
    ϵ⁽ⁱ⁾ = wrk.pulses0
    ϵ⁽ⁱ⁺¹⁾ = wrk.pulses1
    update_result!(wrk, ϵ⁽ⁱ⁺¹⁾, 0)
    println("# iter\tJ_T")
    @printf("%6d\t%.2e\n", 0, wrk.result.J_T)
    i = wrk.result.iter  # = 0, unless continuing from previous optimization
    while !wrk.result.converged
        i = i + 1
        for l = 1:length(ϵ⁽ⁱ⁺¹⁾)
            ϵ⁽ⁱ⁺¹⁾[l] .= 0.9 * ϵ⁽ⁱ⁾[l]
        end
        update_result!(wrk, ϵ⁽ⁱ⁺¹⁾, i)
        @printf("%6d\t%.2e\n", i, wrk.result.J_T)
        check_convergence!(wrk.result)
        ϵ⁽ⁱ⁾, ϵ⁽ⁱ⁺¹⁾ = ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾
    end
    for l = 1:length(ϵ⁽ⁱ⁾)
        wrk.result.optimized_controls[l] = discretize(ϵ⁽ⁱ⁾[l], problem.tlist)
    end
    return wrk.result
end

end


text:
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V, (VV,EV,FV,CV) = Lar.simplex(3, true)
tetra = V, EV,FV,CV
twotetra = Lar.Struct([ tetra, Lar.t(0.25,0.25,0.25), tetra ])
V,EV,FV,CV = Lar.struct2lar(twotetra)
GL.VIEW([ GL.GLGrid(V,FV, GL.COLORS[1],.2) ]);

GL.VIEW(GL.numbering(.5)((V,[[[v] for v=1:size(V,2)],EV,FV,CV]) ));

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.space_arrangement( W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)

# copCF = OK !!

EV = Lar.cop2lar(copEV)
CV = [collect(Set(cat([FV[f]  for  f in cell]))) for cell in cop2lar(copCF)]
FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
FV = [union(EV[e] for e in FE[f]) for f=1:length(FE)]
FV = convert(Lar.Cells, FV)
W = convert(Lar.Points, V') # V'
WW = [[k] for k=1:size(W,2)]

triangulated_faces = Lar.triangulate(V, [copEV, copFE])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
V = convert(Lar.Points, V')
GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99));

EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));


text:
# The following test cases and comments are adapted from http://brainfuck.org/tests.b

using Test
if !@isdefined(BrainFuck)
    include(joinpath(@__DIR__, "..", "src", "BrainFuck.jl"))
end
output = IOBuffer()

# Here are some little programs for testing brainfuck implementations.

@testset "IO test" begin
    # This is for testing i/o; give it a return followed by an EOF. (Try it both
    # with file input — a file consisting only of one blank line — and with
    # keyboard input, i.e. hit return and then ctrl-d (Unix) or ctrl-z (Windows).)
    # It should give two lines of output; the two lines should be identical, and
    # should be lined up one over the other. If that doesn't happen, ten is not
    # coming through as newline on output.
    # The content of the lines tells how input is being processed; each line
    # should be two uppercase letters.
    # Anything with O in it means newline is not coming through as ten on input.
    # LK means newline input is working fine, and EOF leaves the cell unchanged
    # (which I recommend).
    # LB means newline input is working fine, and EOF translates as 0.
    # LA means newline input is working fine, and EOF translates as -1.
    # Anything else is fairly unexpected.
    iotest = ">,>+++++++++,>+++++++++++[<++++++<++++++<+>>>-]<<.>.<<-.>.>.<<."

    BrainFuck.interpret(iotest, input = "\n", output = output)
    @test String(take!(output)) == "LK\nLK\n"

    BrainFuck.interpret(iotest, input = open(joinpath(@__DIR__, "empty.txt")), output = output)
    @test String(take!(output)) == "LK\nLK\n"
end

@testset "Memory size" begin
    memorysizetest = "++++[>++++++<-]>[>+++++>+++++++<<-]>>++++<[[>[[>>+<<-]<]>>>-]>-[>+>+<<-]>]+++++[>+++++++<<++>-]>.<<."
    # Goes to cell 30000 and reports from there with a #. (Verifies that the
    # array is big enough.)
    BrainFuck.interpret(memorysizetest, output = output)
    @test String(take!(output)) == "#\n"

    BrainFuck.compile(memorysizetest)(output = output)
    @test String(take!(output)) == "#\n"
end

@testset "Several obscure tests" begin
    # Tests for several obscure problems. Should output an H.
    obscuretest = """[]++++++++++[>>+>+>++++++[<<+<+++>>>-]<<<<-]
    "A*\$";?@![#>>+<<]>[>>]<<<<[>++<[-]]>.>."""
    @test BrainFuck.interpret(obscuretest, output = String) == "H\n"
    @test BrainFuck.compile(obscuretest)(output = String) == "H\n"
end

@testset "Unmatched brackets" begin
    # Should ideally give error message "unmatched [" or the like, and not give
    # any output. Not essential.
    unmatchedleft = "+++++[>+++++++>++<<-]>.>.["
    @test_throws ArgumentError BrainFuck.interpret(unmatchedleft)
    @test_throws ArgumentError BrainFuck.compile(unmatchedleft)

    # Should ideally give error message "unmatched ]" or the like, and not give
    # any output. Not essential.
    unmatchedright = "+++++[>+++++++>++<<-]>.>.]["
    @test_throws ArgumentError BrainFuck.interpret(unmatchedright)
    @test_throws ArgumentError BrainFuck.compile(unmatchedright)
end

@testset "Deep brackets" begin
    # [Daniel Cristofani's] pathological program rot13.b is good for testing the
    # response to deep brackets; the input "~mlk zyx" should produce the output "~zyx mlk".
    rot13 = """
        ,
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>++++++++++++++<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>>+++++[<----->-]<<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>++++++++++++++<-
        [>+<-[>+<-[>+<-[>+<-[>+<-
        [>++++++++++++++<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>>+++++[<----->-]<<-
        [>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-[>+<-
        [>++++++++++++++<-
        [>+<-]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
        ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]>.[-]<,]
        """

    @test BrainFuck.interpret(rot13, input = "~mlk zyx", output = String) == "~zyx mlk"
    @test BrainFuck.compile(rot13)(input = "~mlk zyx", output = String) == "~zyx mlk"
end

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
# Test reading and writing of SAC files
using Test
using SAC

@testset "IO" begin
    mktemp() do tempfile, tempio

        # Reading
        let s = SAC.sample()
            # Reals
            @test s.t[1:100] ≈ Float32[-0.09728,-0.09728,-0.09856,-0.09856,-0.09728,-0.096,-0.09472,-0.09344,-0.09344,-0.09344,-0.09344,-0.09344,-0.09472,-0.09472,-0.09344,-0.09344,-0.09216,-0.09216,-0.09216,-0.09216,-0.09088,-0.09088,-0.09216,-0.09344,-0.09472,-0.09472,-0.09472,-0.09472,-0.09472,-0.09472,-0.09344,-0.09344,-0.09216,-0.09088,-0.09088,-0.09216,-0.09216,-0.09216,-0.09344,-0.09472,-0.096,-0.09856,-0.09856,-0.09856,-0.09728,-0.09728,-0.09856,-0.09984,-0.09984,-0.09984,-0.09984,-0.09984,-0.10112,-0.10112,-0.10112,-0.10112,-0.1024,-0.1024,-0.10368,-0.1024,-0.10496,-0.10496,-0.10624,-0.10368,-0.10368,-0.1024,-0.10368,-0.10368,-0.10368,-0.10368,-0.10496,-0.10624,-0.10624,-0.10496,-0.10368,-0.10368,-0.10496,-0.10624,-0.10496,-0.10496,-0.10496,-0.10624,-0.1088,-0.1088,-0.1088,-0.1088,-0.1088,-0.10624,-0.10496,-0.10496,-0.10496,-0.10624,-0.10624,-0.10624,-0.10368,-0.10368,-0.10368,-0.10368,-0.1024,-0.10112]
            # Integers
            @test s.nzyear == 1981
            # Strings
            @test s.kevnm == "K8108838"
            @test s.kstnm == "CDV"
            @test s.khole == SAC.sac_cnull
            # Logicals
            @test s.leven == true
        end

        # Writing
        let s1 = SAC.sample()
            @test begin
                write([s1], [tempfile])
                s2 = SAC.read(tempfile)
                s1 == s2
            end
            @test begin
                s1 == SAC.sample()
                write(s1, tempfile)
                s2 = SAC.read(tempfile)
                s1 == s2
            end
        end

        # Test for npts matching file length
        let s = SAC.sample()
            write(s, tempfile)
            data = read(tempfile)
            write(tempfile, data[1:end-4]) # Remove last data point
            @test_throws ErrorException SAC.read(tempfile)
            s′ = SAC.read(tempfile, check_npts=false)
            s.t = s.t[1:end-1] # Remove last datum
            @test s.t == s′.t
        end

        let s = SAC.sample()
            s.kevnm = SAC.sac_cnull
            write(s, tempfile)
            local d = open(tempfile, "r") do f
                read(f)
            end
            @test String(d[449:(449+15)]) == "-12345  -12345  "
        end

        # Handling spaces in kevnm correctly
        let s = SACtr(0.1, 2), kevnm = "1234567 910"
            s.kevnm = kevnm
            write(s, tempfile)
            local s2 = SAC.read(tempfile)
            @test s2.kevnm == kevnm
        end

        # Cutting
        let
            write(SAC.sample(), tempfile)
            local s1 = SAC.read(tempfile)
            cut!(s1, 12, 13)
            local s2 = SAC.read_cut(tempfile, 12, 13)
            @test s1 == s2
        end

        # Printing
        let
            open(tempfile, "w") do f
                show(f, SAC.SACtr(1, 1))
            end
            @test read(tempfile, String) == "SAC.SACtr(delta=1.0, b=0.0, npts=1)"
        end
        let
            teststring = VERSION < v"0.6-" ?
                "SAC.SACtr[SAC.SACtr(delta=1.0, b=0.0, npts=1),SAC.SACtr(delta=1.0, b=0.0, npts=1)]" :
                VERSION >= v"0.7-" ?
                "SACtr[SAC.SACtr(delta=1.0, b=0.0, npts=1), SAC.SACtr(delta=1.0, b=0.0, npts=1)]" :
                "SAC.SACtr[SAC.SACtr(delta=1.0, b=0.0, npts=1), SAC.SACtr(delta=1.0, b=0.0, npts=1)]"
            open(tempfile, "w") do f
                show(f, [SAC.SACtr(1, 1), SAC.SACtr(1, 1)])
            end
            @test read(tempfile, String) == teststring
        end

    end
end


text:
###
### **acb.h** -- complex numbers
###

### Types, macros and constants

### Memory management
arbcall"void acb_init(acb_t x)"
arbcall"void acb_clear(acb_t x)"
#mo arbcall"acb_ptr _acb_vec_init(slong n)" # clashes with similar method for arb
arbcall"void _acb_vec_clear(acb_ptr v, slong n)"
arbcall"slong acb_allocated_bytes(const acb_t x)"
arbcall"slong _acb_vec_allocated_bytes(acb_srcptr vec, slong len)"
#mo arbcall"double _acb_vec_estimate_allocated_bytes(slong len, slong prec)" # clashes with similar method for arb

### Basic manipulation
arbcall"void acb_zero(acb_t z)"
arbcall"void acb_one(acb_t z)"
arbcall"void acb_onei(acb_t z)"
arbcall"void acb_set(acb_t z, const acb_t x)"
arbcall"void acb_set_ui(acb_t z, ulong x)"
arbcall"void acb_set_si(acb_t z, slong x)"
arbcall"void acb_set_d(acb_t z, double x)"
#ni arbcall"void acb_set_fmpz(acb_t z, const fmpz_t x)"
arbcall"void acb_set_arb(acb_t z, const arb_t c)"
arbcall"void acb_set_si_si(acb_t z, slong x, slong y)"
arbcall"void acb_set_d_d(acb_t z, double x, double y)"
#ni arbcall"void acb_set_fmpz_fmpz(acb_t z, const fmpz_t x, const fmpz_t y)"
arbcall"void acb_set_arb_arb(acb_t z, const arb_t x, const arb_t y)"
#ni arbcall"void acb_set_fmpq(acb_t z, const fmpq_t x, slong prec)"
arbcall"void acb_set_round(acb_t z, const acb_t x, slong prec)"
#ni arbcall"void acb_set_round_fmpz(acb_t z, const fmpz_t x, slong prec)"
arbcall"void acb_set_round_arb(acb_t z, const arb_t x, slong prec)"
arbcall"void acb_swap(acb_t z, acb_t x)"
arbcall"void acb_add_error_mag(acb_t x, const mag_t err)"
arbcall"void acb_get_mid(acb_t m, const acb_t x)"

### Input and output
#mo arbcall"void acb_print(const acb_t x)" # clashes with Base.print
#ns arbcall"void acb_fprint(FILE * file, const acb_t x)"
arbcall"void acb_printd(const acb_t x, slong digits)"
#ns arbcall"void acb_fprintd(FILE * file, const acb_t x, slong digits)"
arbcall"void acb_printn(const acb_t x, slong digits, ulong flags)"
#ns arbcall"void acb_fprintn(FILE * file, const acb_t x, slong digits, ulong flags)"

### Random number generation
#ns arbcall"void acb_randtest(acb_t z, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void acb_randtest_special(acb_t z, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void acb_randtest_precise(acb_t z, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void acb_randtest_param(acb_t z, flint_rand_t state, slong prec, slong mag_bits)"

### Precision and comparisons
arbcall"int acb_is_zero(const acb_t z)"
arbcall"int acb_is_one(const acb_t z)"
arbcall"int acb_is_finite(const acb_t z)"
arbcall"int acb_is_exact(const acb_t z)"
arbcall"int acb_is_int(const acb_t z)"
arbcall"int acb_is_int_2exp_si(const acb_t x, slong e)"
arbcall"int acb_equal(const acb_t x, const acb_t y)"
arbcall"int acb_equal_si(const acb_t x, slong y)"
arbcall"int acb_eq(const acb_t x, const acb_t y)"
arbcall"int acb_ne(const acb_t x, const acb_t y)"
arbcall"int acb_overlaps(const acb_t x, const acb_t y)"
arbcall"void acb_union(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_get_abs_ubound_arf(arf_t u, const acb_t z, slong prec)"
arbcall"void acb_get_abs_lbound_arf(arf_t u, const acb_t z, slong prec)"
arbcall"void acb_get_rad_ubound_arf(arf_t u, const acb_t z, slong prec)"
arbcall"void acb_get_mag(mag_t u, const acb_t x)"
arbcall"void acb_get_mag_lower(mag_t u, const acb_t x)"
#ni arbcall"int acb_contains_fmpq(const acb_t x, const fmpq_t y)"
#ni arbcall"int acb_contains_fmpz(const acb_t x, const fmpz_t y)"
arbcall"int acb_contains(const acb_t x, const acb_t y)"
arbcall"int acb_contains_zero(const acb_t x)"
arbcall"int acb_contains_int(const acb_t x)"
arbcall"int acb_contains_interior(const acb_t x, const acb_t y)"
arbcall"slong acb_rel_error_bits(const acb_t x)"
arbcall"slong acb_rel_accuracy_bits(const acb_t x)"
arbcall"slong acb_rel_one_accuracy_bits(const acb_t x)"
arbcall"slong acb_bits(const acb_t x)"
arbcall"void acb_indeterminate(acb_t x)"
arbcall"void acb_trim(acb_t y, const acb_t x)"
arbcall"int acb_is_real(const acb_t x)"
#ni arbcall"int acb_get_unique_fmpz(fmpz_t z, const acb_t x)"

### Complex parts
arbcall"void acb_get_real(arb_t re, const acb_t z)"
arbcall"void acb_get_imag(arb_t im, const acb_t z)"
arbcall"void acb_arg(arb_t r, const acb_t z, slong prec)"
arbcall"void acb_abs(arb_t r, const acb_t z, slong prec)"
arbcall"void acb_sgn(acb_t r, const acb_t z, slong prec)"
arbcall"void acb_csgn(arb_t r, const acb_t z)"

### Arithmetic
arbcall"void acb_neg(acb_t z, const acb_t x)"
arbcall"void acb_neg_round(acb_t z, const acb_t x, slong prec)"
arbcall"void acb_conj(acb_t z, const acb_t x)"
arbcall"void acb_add_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_add_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_add_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_add_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_add(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_sub_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_sub_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_sub_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_sub_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_sub(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_mul_onei(acb_t z, const acb_t x)"
arbcall"void acb_div_onei(acb_t z, const acb_t x)"
arbcall"void acb_mul_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_mul_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_mul_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_mul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_mul(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_mul_2exp_si(acb_t z, const acb_t x, slong e)"
#ni arbcall"void acb_mul_2exp_fmpz(acb_t z, const acb_t x, const fmpz_t e)"
arbcall"void acb_sqr(acb_t z, const acb_t x, slong prec)"
arbcall"void acb_cube(acb_t z, const acb_t x, slong prec)"
arbcall"void acb_addmul(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_addmul_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_addmul_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_addmul_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_addmul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_submul(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_submul_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_submul_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_submul_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_submul_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_inv(acb_t z, const acb_t x, slong prec)"
arbcall"void acb_div_ui(acb_t z, const acb_t x, ulong y, slong prec)"
arbcall"void acb_div_si(acb_t z, const acb_t x, slong y, slong prec)"
#ni arbcall"void acb_div_fmpz(acb_t z, const acb_t x, const fmpz_t y, slong prec)"
arbcall"void acb_div_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_div(acb_t z, const acb_t x, const acb_t y, slong prec)"

### Dot product
arbcall"void acb_dot_precise(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void acb_dot_simple(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void acb_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void acb_approx_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)"

### Mathematical constants
arbcall"void acb_const_pi(acb_t y, slong prec)"

### Powers and roots
arbcall"void acb_sqrt(acb_t r, const acb_t z, slong prec)"
arbcall"void acb_sqrt_analytic(acb_t r, const acb_t z, int analytic, slong prec)"
arbcall"void acb_rsqrt(acb_t r, const acb_t z, slong prec)"
arbcall"void acb_rsqrt_analytic(acb_t r, const acb_t z, int analytic, slong prec)"
#ni arbcall"void acb_quadratic_roots_fmpz(acb_t r1, acb_t r2, const fmpz_t a, const fmpz_t b, const fmpz_t c, slong prec)"
arbcall"void acb_root_ui(acb_t r, const acb_t z, ulong k, slong prec)"
#ni arbcall"void acb_pow_fmpz(acb_t y, const acb_t b, const fmpz_t e, slong prec)"
arbcall"void acb_pow_ui(acb_t y, const acb_t b, ulong e, slong prec)"
arbcall"void acb_pow_si(acb_t y, const acb_t b, slong e, slong prec)"
arbcall"void acb_pow_arb(acb_t z, const acb_t x, const arb_t y, slong prec)"
arbcall"void acb_pow(acb_t z, const acb_t x, const acb_t y, slong prec)"
arbcall"void acb_pow_analytic(acb_t r, const acb_t x, const acb_t y, int analytic, slong prec)"
arbcall"void acb_unit_root(acb_t res, ulong order, slong prec)"

### Exponentials and logarithms
arbcall"void acb_exp(acb_t y, const acb_t z, slong prec)"
arbcall"void acb_exp_pi_i(acb_t y, const acb_t z, slong prec)"
arbcall"void acb_exp_invexp(acb_t s, acb_t t, const acb_t z, slong prec)"
arbcall"void acb_expm1(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_log(acb_t y, const acb_t z, slong prec)"
arbcall"void acb_log_analytic(acb_t r, const acb_t z, int analytic, slong prec)"
arbcall"void acb_log1p(acb_t z, const acb_t x, slong prec)"

### Trigonometric functions
arbcall"void acb_sin(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_cos(acb_t c, const acb_t z, slong prec)"
arbcall"void acb_sin_cos(acb_t s, acb_t c, const acb_t z, slong prec)"
arbcall"void acb_tan(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_cot(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_sin_pi(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_cos_pi(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_sin_cos_pi(acb_t s, acb_t c, const acb_t z, slong prec)"
arbcall"void acb_tan_pi(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_cot_pi(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_sec(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_csc(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_csc_pi(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_sinc(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_sinc_pi(acb_t s, const acb_t z, slong prec)"

### Inverse trigonometric functions
arbcall"void acb_asin(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_acos(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_atan(acb_t res, const acb_t z, slong prec)"

### Hyperbolic functions
arbcall"void acb_sinh(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_cosh(acb_t c, const acb_t z, slong prec)"
arbcall"void acb_sinh_cosh(acb_t s, acb_t c, const acb_t z, slong prec)"
arbcall"void acb_tanh(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_coth(acb_t s, const acb_t z, slong prec)"
arbcall"void acb_sech(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_csch(acb_t res, const arb_t z, slong prec)"

### Inverse hyperbolic functions
arbcall"void acb_asinh(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_acosh(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_atanh(acb_t res, const acb_t z, slong prec)"

### Lambert W function
#ni arbcall"void acb_lambertw_asymp(acb_t res, const acb_t z, const fmpz_t k, slong L, slong M, slong prec)"
#ni arbcall"int acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)"
#ni arbcall"void acb_lambertw_bound_deriv(mag_t res, const acb_t z, const acb_t ez1, const fmpz_t k)"
#ni arbcall"void acb_lambertw(acb_t res, const acb_t z, const fmpz_t k, int flags, slong prec)"

### Rising factorials
arbcall"void acb_rising_ui_bs(acb_t z, const acb_t x, ulong n, slong prec)"
arbcall"void acb_rising_ui_rs(acb_t z, const acb_t x, ulong n, ulong step, slong prec)"
arbcall"void acb_rising_ui_rec(acb_t z, const acb_t x, ulong n, slong prec)"
arbcall"void acb_rising_ui(acb_t z, const acb_t x, ulong n, slong prec)"
arbcall"void acb_rising(acb_t z, const acb_t x, const acb_t n, slong prec)"

### Gamma function
arbcall"void acb_gamma(acb_t y, const acb_t x, slong prec)"
arbcall"void acb_rgamma(acb_t y, const acb_t x, slong prec)"
arbcall"void acb_lgamma(acb_t y, const acb_t x, slong prec)"
arbcall"void acb_digamma(acb_t y, const acb_t x, slong prec)"
arbcall"void acb_log_sin_pi(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_polygamma(acb_t res, const acb_t s, const acb_t z, slong prec)"
arbcall"void acb_barnes_g(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_log_barnes_g(acb_t res, const acb_t z, slong prec)"

### Zeta function
arbcall"void acb_zeta(acb_t z, const acb_t s, slong prec)"
arbcall"void acb_hurwitz_zeta(acb_t z, const acb_t s, const acb_t a, slong prec)"
arbcall"void acb_bernoulli_poly_ui(acb_t res, ulong n, const acb_t x, slong prec)"

### Polylogarithms
arbcall"void acb_polylog(acb_t w, const acb_t s, const acb_t z, slong prec)"
arbcall"void acb_polylog_si(acb_t w, slong s, const acb_t z, slong prec)"

### Arithmetic-geometric mean
arbcall"void acb_agm1(acb_t m, const acb_t z, slong prec)"
arbcall"void acb_agm1_cpx(acb_ptr m, const acb_t z, slong len, slong prec)"
arbcall"void acb_agm(acb_t m, const acb_t x, const acb_t y, slong prec)"

### Other special functions
arbcall"void acb_chebyshev_t_ui(acb_t a, ulong n, const acb_t x, slong prec)"
arbcall"void acb_chebyshev_u_ui(acb_t a, ulong n, const acb_t x, slong prec)"
arbcall"void acb_chebyshev_t2_ui(acb_t a, acb_t b, ulong n, const acb_t x, slong prec)"
arbcall"void acb_chebyshev_u2_ui(acb_t a, acb_t b, ulong n, const acb_t x, slong prec)"

### Piecewise real functions
arbcall"void acb_real_abs(acb_t res, const acb_t z, int analytic, slong prec)"
arbcall"void acb_real_sgn(acb_t res, const acb_t z, int analytic, slong prec)"
arbcall"void acb_real_heaviside(acb_t res, const acb_t z, int analytic, slong prec)"
arbcall"void acb_real_floor(acb_t res, const acb_t z, int analytic, slong prec)"
arbcall"void acb_real_ceil(acb_t res, const acb_t z, int analytic, slong prec)"
arbcall"void acb_real_max(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec)"
arbcall"void acb_real_min(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec)"
arbcall"void acb_real_sqrtpos(acb_t res, const acb_t z, int analytic, slong prec)"

### Vector functions
arbcall"void _acb_vec_zero(acb_ptr A, slong n)"
arbcall"int _acb_vec_is_zero(acb_srcptr vec, slong len)"
arbcall"int _acb_vec_is_real(acb_srcptr v, slong len)"
arbcall"void _acb_vec_set(acb_ptr res, acb_srcptr vec, slong len)"
arbcall"void _acb_vec_set_round(acb_ptr res, acb_srcptr vec, slong len, slong prec)"
arbcall"void _acb_vec_neg(acb_ptr res, acb_srcptr vec, slong len)"
arbcall"void _acb_vec_add(acb_ptr res, acb_srcptr vec1, acb_srcptr vec2, slong len, slong prec)"
arbcall"void _acb_vec_sub(acb_ptr res, acb_srcptr vec1, acb_srcptr vec2, slong len, slong prec)"
arbcall"void _acb_vec_scalar_submul(acb_ptr res, acb_srcptr vec, slong len, const acb_t c, slong prec)"
arbcall"void _acb_vec_scalar_addmul(acb_ptr res, acb_srcptr vec, slong len, const acb_t c, slong prec)"
arbcall"void _acb_vec_scalar_mul(acb_ptr res, acb_srcptr vec, slong len, const acb_t c, slong prec)"
arbcall"void _acb_vec_scalar_mul_ui(acb_ptr res, acb_srcptr vec, slong len, ulong c, slong prec)"
arbcall"void _acb_vec_scalar_mul_2exp_si(acb_ptr res, acb_srcptr vec, slong len, slong c)"
arbcall"void _acb_vec_scalar_mul_onei(acb_ptr res, acb_srcptr vec, slong len)"
arbcall"void _acb_vec_scalar_div_ui(acb_ptr res, acb_srcptr vec, slong len, ulong c, slong prec)"
arbcall"void _acb_vec_scalar_div(acb_ptr res, acb_srcptr vec, slong len, const acb_t c, slong prec)"
arbcall"void _acb_vec_scalar_mul_arb(acb_ptr res, acb_srcptr vec, slong len, const arb_t c, slong prec)"
arbcall"void _acb_vec_scalar_div_arb(acb_ptr res, acb_srcptr vec, slong len, const arb_t c, slong prec)"
#ni arbcall"void _acb_vec_scalar_mul_fmpz(acb_ptr res, acb_srcptr vec, slong len, const fmpz_t c, slong prec)"
#ni arbcall"void _acb_vec_scalar_div_fmpz(acb_ptr res, acb_srcptr vec, slong len, const fmpz_t c, slong prec)"
arbcall"slong _acb_vec_bits(acb_srcptr vec, slong len)"
arbcall"void _acb_vec_set_powers(acb_ptr xs, const acb_t x, slong len, slong prec)"
arbcall"void _acb_vec_unit_roots(acb_ptr z, slong order, slong len, slong prec)"
#ni arbcall"void _acb_vec_add_error_arf_vec(acb_ptr res, arf_srcptr err, slong len)"
#ni arbcall"void _acb_vec_add_error_mag_vec(acb_ptr res, mag_srcptr err, slong len)"
arbcall"void _acb_vec_indeterminate(acb_ptr vec, slong len)"
arbcall"void _acb_vec_trim(acb_ptr res, acb_srcptr vec, slong len)"
#ni arbcall"int _acb_vec_get_unique_fmpz_vec(fmpz * res,  acb_srcptr vec, slong len)"
arbcall"void _acb_vec_sort_pretty(acb_ptr vec, slong len)"


text:
include("init/resources.jl")
include("init/components.jl")


text:
function find_other_triangle_edge(v1::Integer,v2::Integer,skip::Integer,t::AbstractArray{T,2}) where {T<:Integer}
    for i in 1:size(t)[end]
        if in(v1,t[:,i]) & in(v2,t[:,i]) & !(i==skip)
            return i
        end
    end
    return -1
end

function find_triangle_vertex(v::Integer,t::AbstractArray{T,2}) where {T<:Integer}

    for i in 1:size(t,2)
        if in(v,t[:,i])
            return i
        end
    end
    return size(t,2) + 1 # -1
end

### Deprecaction Warnings
FaceRing(v,faces) = error("Deprecated: use FaceVRing")
FaceRingS(v,faces) = error("Deprecated: use DoubleVertexVRing")
VeretexRing(v,faces) = error("Deprecated: use VertexVRing")

struct FaceVRing
    v::Int
    faces::Array{Int,2}
end

start(iter::FaceVRing) = find_triangle_vertex(iter.v,iter.faces)
done(iter::FaceVRing,ti::Int) = ti<=size(iter.faces,2) ? false : true

function next(iter::FaceVRing,i::Int)

     v = iter.v
    nexti = find_triangle_vertex(iter.v,iter.faces[:,i+1:end]) + i  # possible botleneck here
    return i, nexti
end

Base.iterate(iter::FaceVRing) = next(iter,start(iter))
function Base.iterate(iter::FaceVRing,ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end

struct DoubleVertexVRing
    v::Int
    faces::Array{Int,2}
end

start(iter::DoubleVertexVRing) = find_triangle_vertex(iter.v,iter.faces)
done(iter::DoubleVertexVRing,ti::Int) = ti<=size(iter.faces,2) ? false : true

function next(iter::DoubleVertexVRing,i::Int)

    v = iter.v
    nexti = find_triangle_vertex(iter.v,iter.faces[:,i+1:end]) + i  # possible botleneck here
    face = iter.faces[:,i]
    w = face.==v
    cw = w[[3,1,2]]
    ccw = w[[2,3,1]]

    return (face[cw]...,face[ccw]...), nexti
end

Base.iterate(iter::DoubleVertexVRing) = next(iter,start(iter))
function Base.iterate(iter::DoubleVertexVRing,ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end

### Unordered Veretex ring

struct VertexVRing
    v::Int
    faces::Array{Int,2}
end

start(iter::VertexVRing) = find_triangle_vertex(iter.v,iter.faces)
done(iter::VertexVRing,ti::Int) = ti<=size(iter.faces,2) ? false : true

function next(iter::VertexVRing,ti::Int)

    v = iter.v
    faces = iter.faces
    face = faces[:,ti]
    cw = (face.==v)[[3,1,2]]
    vi, = face[cw]

    nexti = find_triangle_vertex(iter.v,iter.faces[:,ti+1:end]) + ti  # possible botleneck
    return vi, nexti
end

Base.iterate(iter::VertexVRing) = next(iter,start(iter))
function Base.iterate(iter::VertexVRing,ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end


text:
# Reference: Experience Study Calculations, 2016, Society of Actuaries
# https://www.soa.org/globalassets/assets/Files/Research/2016-10-experience-study-calculations.pdf
# examples on pages 41-45
@testset "Fractional Year and distribution of deaths" begin
    @testset "SOA Experience Study Calcuations distribution examples" begin
        soa_mort = UltimateMortality([0.12])

        methods = [Balducci(), Uniform(), Constant()]
        time_targets = Dict(
            1 / 12 => [0.9888, 0.9900, 0.9894],
            6 / 12 => [0.9362, 0.9400, 0.9381],
            7 / 12 => [0.9263, 0.9300, 0.9281],
            12 / 12 => [0.8800, 0.8800, 0.8800],
        )

        # test whole ages with assumption argument
        for method in methods
            @test survival(soa_mort, 0, 1, method) == 0.88
            @test survival(soa_mort, 1, method) == 0.88
            @test decrement(soa_mort, 0, 1, method) == 0.12
            @test decrement(soa_mort, 1, method) == 0.12

            # test floating point whole ages
            @test decrement(soa_mort, 0, 12 / 12, method) == 0.12
            @test decrement(soa_mort, 12 / 12, method) == 0.12
        end

        # test fractional ages
        for i = 1:length(methods)
            for (t, target) in time_targets
                @test round(survival(soa_mort, t, methods[i]), digits=4) ==
                      target[i]

                @test round(MortalityTables.decrement_partial_year(soa_mort, 0, t, methods[i]),
                    digits=4) ≈  1 - target[i]

                @test round(decrement(soa_mort, t, methods[i]), digits=4) ==
                      round(1 - target[i], digits=4)
            end
        end

        # test time zero when given distribution of deaths
        for m in methods
            @test decrement(soa_mort, 0, m) == 0.0
            @test survival(soa_mort, 0, m) == 1.0
        end
    end
    @testset "Multi-year examples" begin
        mort = UltimateMortality([0.20, 0.50])
        methods = [Balducci(), Uniform(), Constant()]

        # these sample values calculated manually
        time_targets = Dict(
            1 + 1 / 12 => [0.7385, 0.7667, 0.7551],
            1 + 6 / 12 => [0.5333, 0.6000, 0.5657],
            1 + 7 / 12 => [0.5053, 0.5667, 0.5339],
            1 + 12 / 12 => [0.4000, 0.4000, 0.4000],
        )

        # test whole ages with assumption argument
        for method in methods
            @test decrement(mort, 1, method) ≈ 0.20
            @test decrement(mort, 2, method) ≈ 0.20 + 0.8 * 0.5
            @test decrement(mort, 1, 2, method) ≈ 0.50
        end

        # test fractional ages
        for i = 1:length(methods)
            for (t, target) in time_targets
                @test round(survival(mort, t, methods[i]), digits=4) ==
                      target[i]
                @test round(decrement(mort, t, methods[i]), digits=4) ==
                      round(1 - target[i], digits=4)
            end
        end
    end

    @testset "Error when asking for a fractional without assumption" begin
        mort = UltimateMortality([0.20, 0.50])
        @test_throws MethodError decrement(mort, 0, 1, 1.5)
        @test_throws MethodError survival(mort, 0, 1, 1.5)
    end

    @testset "Issue #60 - starting with non-integer ages" begin
        m = UltimateMortality([0.5 for i in 1:8])

        @test survival(m, 1, 2) ≈ 0.5
        @test survival(m, 1.5, 2.5, Constant()) ≈ 0.5
        @test survival(m, 1.5, 3.5, Constant()) ≈ 0.25
        @test survival(m, 1.5, 1.5 + eps(), Constant()) ≈ 1.0
        @test survival(m, 1, 1 + eps(), Constant()) ≈ 1.0

    end

    @testset "Issue #88 - Stackoverflow" begin
        @test survival(MortalityTables.mortality_vector([0.5,0.5],start_age=50),50,50.5,MortalityTables.Uniform()) ≈ 1 - 0.5 * 0.5
    end

end


text:
# Adaptor for the API/ABI expected by the Julia runtime code.
function core_parser_hook(code, filename, lineno, offset, options)
    try
        # TODO: Check that we do all this input wrangling without copying the
        # code buffer
        if code isa Core.SimpleVector
            # The C entry points will pass us this form.
            (ptr,len) = code
            code = String(unsafe_wrap(Array, ptr, len))
        end
        io = IOBuffer(code)
        seek(io, offset)

        stream = ParseStream(io)
        rule = options === :all ? :toplevel : options
        if rule !== :toplevel
            # To copy the flisp parser driver, we ignore leading trivia when
            # parsing statements or atoms
            bump_trivia(stream)
        end
        JuliaSyntax.parse(stream; rule=rule)

        if any_error(stream)
            e = Expr(:error, ParseError(SourceFile(code), stream.diagnostics))
            ex = options === :all ? Expr(:toplevel, e) : e
        else
            ex = build_tree(Expr, stream, wrap_toplevel_as_kind=K"None")
            if Meta.isexpr(ex, :None)
                # The None wrapping is only to give somewhere for trivia to be
                # attached; unwrap!
                ex = only(ex.args)
            end
        end

        # Note the next byte in 1-based indexing is `last_byte(stream) + 1` but
        # the Core hook must return an offset (ie, it's 0-based) so the factors
        # of one cancel here.
        last_offset = last_byte(stream)

        # Rewrap result in an svec for use by the C code
        return Core.svec(ex, last_offset)
    catch exc
        @error("JuliaSyntax parser failed — falling back to flisp!",
               exception=(exc,catch_backtrace()),
               offset=offset,
               code=code)
    end
    return Core.Compiler.fl_parse(code, filename, offset, options)
end

# Core._parse gained a `lineno` argument in
# https://github.com/JuliaLang/julia/pull/43876
# Prior to this, the following signature was needed:
function core_parser_hook(code, filename, offset, options)
    core_parser_hook(code, filename, LineNumberNode(0), offset, options)
end

# Hack:
# Meta.parse() attempts to construct a ParseError from a string if it receives
# `Expr(:error)`.
Base.Meta.ParseError(e::JuliaSyntax.ParseError) = e

"""
Connect the JuliaSyntax parser to the Julia runtime so that it replaces the
flisp parser for all parsing work.

That is, JuliaSyntax will be used for `include()` `Meta.parse()`, the REPL, etc.
"""
function enable_in_core!()
    # TODO: Use invoke_in_world to freeze the world age at the time this was enabled.
    Base.eval(Core, :(_parse = $core_parser_hook))
    nothing
end

"""
Revert to the flisp parser for all parsing work.
"""
function disable_in_core!()
    Base.eval(Core, :(_parse = Core.Compiler.fl_parse))
    nothing
end