text:
function momentum_grid(kmin, kmax, nk, nθ; spacing=:momentum)
    kmag = if spacing == :momentum
        range(kmin, stop=kmax, length=nk)
    elseif spacing == :energy
        .√(2range(kmin^2/2, stop=kmax^2/2, length=nk))
    else
        throw(ArgumentError("Unknown spacing $(spacing)"))
    end
    nθ == 1 && return (kmag,kmag,nothing)

    x = range(0, stop=1, length=nθ)
    k = Array{SVector{3,Float64}}(undef, length(kmag), nθ)
    for j = 1:nθ
        s,c = sincospi(x[j])
        k̂ = SVector{3}(s, 0, c)
        for (i,km) in enumerate(kmag)
            k[i,j] = km*k̂
        end
    end
    k,kmag,x*π
end

momentum_grid(Emin::Unitful.Energy, Emax::Unitful.Energy, args...; kwargs...) =
    momentum_grid(√(2austrip(Emin)), √(2austrip(Emax)), args...; kwargs...)

export momentum_grid


text:
struct NoiselessQCtrl
        drift_Ham
        ctrl_Hams
        n_steps
        transfer_time
        init_density_op
        target_density_op
        init_guess
end

struct GateTransfer
    drift_Ham
    ctrl_Hams
    n_steps
    transfer_time
    tar_gate
    init_guess
end

struct StateTransfer
    drift_Ham
    ctrl_Hams
    n_steps
    transfer_time
    init_state
    tar_state
    init_guess
end

export GateTransfer, StateTransfer

ext: jl


text:
module PhysicsTutorials
# inspiration taken from DiffEqTutorials.jl

using Weave, Pkg

repo_directory = joinpath(@__DIR__,"..")
cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")
htmlfile = joinpath(@__DIR__, "..", "templates", "julia_html.tpl")

function weave_file(folder,file,build_list=(:html,:pdf,:notebook); kwargs...)
  jmdfile = joinpath(repo_directory,"tutorials",folder,file)
  args = Dict{Symbol,String}(:folder=>folder,:file=>file)
  if :script ∈ build_list
    println("Building Script")
    tangle(jmdfile)
  end
  if :html ∈ build_list
    println("Building HTML")
    weave(jmdfile,doctype = "md2html",args=args; css=cssfile, template=htmlfile, kwargs...)
  end
  if :pdf ∈ build_list
    println("Building PDF")
    dir = joinpath(repo_directory,"pdf",folder)
    weave(jmdfile,doctype="md2pdf",args=args; template=latexfile, kwargs...)
  end
  if :github ∈ build_list
    println("Building Github Markdown")
    dir = joinpath(repo_directory,"markdown",folder)
    weave(jmdfile,doctype = "github",args=args; kwargs...)
  end
  if :notebook ∈ build_list
    println("Building Notebook")
    dir = joinpath(repo_directory,"notebook",folder)
    # Weave.convert_doc(jmdfile,joinpath(dir,file[1:end-4]*".ipynb"))
    Weave.notebook(jmdfile; out_path=dirname(jmdfile))
  end
end

function weave_all()
  for folder in readdir(joinpath(repo_directory,"tutorials"))
    weave_folder(folder)
  end
end

function weave_folder(folder)
  for file in readdir(joinpath(repo_directory,"tutorials",folder))
    println("Building $(joinpath(folder,file)))")
    try
      weave_file(folder,file)
    catch
    end
  end
end

function tutorial_footer(folder=nothing, file=nothing; remove_homedir=true)
    display("text/markdown", """
    ## Appendix

     This tutorial is part of the [PhysicsTutorials.jl](https://github.com/JuliaPhysics/PhysicsTutorials.jl) repository.
    """)
    # if folder !== nothing && file !== nothing
    #     display("text/markdown", """
    #     To locally run this tutorial, do the following commands:
    #     ```
    #     using PhysicsTutorials
    #     PhysicsTutorials.weave_file("$folder","$file")
    #     ```
    #     """)
    # end

    ctx = Pkg.API.Context()
    pkgs = Pkg.Display.status(Pkg.API.Context(), use_as_api=true);
    projfile = ctx.env.project_file
    remove_homedir && (projfile = replace(projfile, homedir() => "~"))

    display("text/markdown","""
    Package Information:
    """)

    md = ""
    md *= "```\nStatus `$(projfile)`\n"

    for pkg in pkgs
        if !isnothing(pkg.new.ver)
          md *= "[$(string(pkg.uuid))] $(string(pkg.name)) $(string(pkg.new.ver))\n"
        else
          md *= "[$(string(pkg.uuid))] $(string(pkg.name))\n"
        end
    end
    md *= "```"
    display("text/markdown", md)
end

function open_notebooks()
  Base.eval(Main, Meta.parse("import IJulia"))
  path = joinpath(repo_directory,"tutorials")
  IJulia.notebook(;dir=path)
end

end # module


text:
using BenchmarkTools

test_input = parse.(Int64, readlines("1/input.txt"));

function sum2(input, sum2val, is_sorted=false)
    sorted_exp = !is_sorted ? sort(input) : input;

    itr_exp = sorted_exp[sorted_exp .< sum2val / 2];
    chk_exp = sorted_exp[sum2val / 2 .<= sorted_exp .< sum2val]
    for expense in itr_exp
        matches = searchsorted(chk_exp, sum2val - expense)
        if matches.start <= matches.stop
            return expense * (sum2val - expense)
        end
    end
    return 0
end

function sum3(input, sum3val, is_sorted=false)
    sorted_exp = !is_sorted ? sort(input) : input;

    for expense in sorted_exp
        prod23 = sum2(sorted_exp[sorted_exp .> expense], sum3val - expense, true);
        if prod23 != 0
            return expense * prod23;
        end
    end
    return 0
end

# println("Part 1: $(sum2(test_input, 2020))")
# println("Part 2: $(sum3(test_input, 2020))")

@benchmark sum2(test_input, 2020)
@benchmark sum3(test_input, 2020)


text:
"""
# Get information about a time series in the FRED database.

    get_metadata(symbol)

## Arguments

    - `symbol` : String specifying the name (id) of the time series.

## Examples

```jldoctest
julia> get_metadata("GDPC1")
julia> get_metadata("FEDFUNDS")
julia> get_metadata("T10Y2Y")
```
"""
function get_metadata(symbol::String)

    url = "https://api.stlouisfed.org/fred/series"

    parameters = Dict("api_key" => ENV["API_FRED"],
                      "file_type" => "json",
                      "series_id" => symbol)

    response = HTTP.request("GET", url; query = parameters)
    results = JSON.parse(String(response.body))["seriess"][1]

    println("")
    println("Metadata for: " * symbol)
    println("Title: " * results["title"])
    println("Units: " * results["units"])
    println("Adjustment: " * results["seasonal_adjustment"])
    println("Frequency: " * results["frequency"])
    println("Notes: " * results["notes"])

    return(results)
end


text:
import ..RvSpectMLBase.InstrumentsCommon: get_λ_range
import ..RvSpectMLBase.InstrumentsCommon: choose_obs_idx_for_init_guess

function get_λ_range(data::ST) where { T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
                                       IT<:AnyTheoreticalInstrument, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT} }
   (λmin, λmax) = extrema(data.λ)
   return (min=λmin, max=λmax)
end

""" `calc_λs`
Generate an array of n wavelengths uniformly spaced in log(lambda)
"""
function calc_λs end

""" `calc_λs(λ_min, λ_max, n)`
Generate an array of n wavelengths uniformly spaced in log(lambda) from λ_min to λ_max.
"""
function calc_λs(λ_min::Real, λ_max::Real, n::Integer)
   r = range(log(λ_min), stop=log(λ_max), length=n)
   return exp.(r)
end

""" `calc_λs(inst)`
Generate an array of wavelengths based on a theoretical instrument's properties
"""
function calc_λs(inst::TheoreticalInstrument1D)
   Δlnλ = 1.0/inst.resolution
   logλ_min = log(λ_min(inst))
   logλ_max = log(λ_max(inst))
   n = ceil(Int,(logλ_max-logλ_min) / Δlnλ )
   @assert 1 <= n <= max_pixel(inst)
   r = range(logλ_min, stop=logλ_max, length=n)
   return exp.(r)
end

""" `calc_λs(inst)`
Compute the wavelength for one pixel of a theoretical 1d instrument
"""
function calc_λ(pixel, inst::TheoreticalInstrument1D)
   @assert min_pixel(inst) 1<= pixel <= max_pixel(inst)
   Δlnλ = 1.0/inst.resolution
   logλ_min = log(λ_min(inst))
   logλ_max = log(λ_max(inst))
   n = ceil(Int,(logλ_max-logλ_min) / Δlnλ )
   @assert 1 <= n <= max_pixel(inst)
   r = range(logλ_min, stop=logλ_max, length=n)
   @assert 1<= pixel <= length(r)
   exp(r[pixel])
end

""" `calc_λs(inst)`
Generate an array of wavelengths based on a theoretical instrument's properties
"""
function calc_λs(inst::TheoreticalInstrument2D)
   num_orders = length(inst.λ_min)
   num_pixels = inst.pixels_per_order
   λ = zeros(num_pixels,num_orders)
   for i in 1:num_orders
      #λ_max_order = exp(logλ_min + (logλ_max-logλ_min)*(i  )/num_orders)
      #λ_min_order = exp(logλ_min + (logλ_max-logλ_min)*(i-1)/num_orders)
      λ[:,i] .= calc_λs(inst.λ_min[i], inst.λ_max[i], num_pixels)
   end
   return λ
end

function choose_obs_idx_for_init_guess(df::DataFrame, inst::AnyTheoreticalInstrument)
   @assert size(df,1) >= 1
   return 1
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
"""
Table Schema main type
https://github.com/frictionlessdata/TableSchema.jl#schema
"""
mutable struct Schema
    errors::Array{SchemaError}
    descriptor::Dict
    primary_key::Array{String}
    foreign_keys::Array{Dict}
    missing_values::Array{String}
    fields::Array{Field}

    function Schema(d::Dict, strict::Bool=false)
        fls = haskey(d, "fields") ? [ Field(f) for f in d["fields"] ] : []
        pk = haskey(d, "primaryKey") ? d["primaryKey"] : []
        if isa(pk, String); pk = [pk]; end
        fks = haskey(d, "foreignKeys") ? d["foreignKeys"] : []
        if isa(fks, String); fks = [fks]; end
        mvs = haskey(d, "missingValues") ? d["missingValues"] : []
        if isa(mvs, String); mvs = [mvs]; end

        # Validate schema types
        if !(isa(pk, Array))
            s = new([SchemaError("primaryKey must be a string array")], Dict(),[],[],[],[])
        elseif !(isa(fks, Array))
            s = new([SchemaError("foreignKeys must be a dictionary array")], Dict(),[],[],[],[])
        elseif !(isa(mvs, Array))
            s = new([SchemaError("missingValues must be a string array")], Dict(),[],[],[],[])
        else
            s = new([], d, pk, fks, mvs, fls)
        end

        # Detailed validation
        validate(s, strict)
    end

    function Schema(a::Array, strict::Bool=false)
        s = new([SchemaError("Descriptor must be an object, not an array")], Dict(),[],[],[],[])
        validate(s, strict)
    end

    Schema(filename::String, strict::Bool=false) =
        Schema(fetch_json(filename), strict)

    Schema(strict::Bool=false) =
        Schema(Dict(), strict)
end

function fetch_json(filename::String)
    if match(r"^https?://", filename) !== nothing
        req = request("GET", filename)
        JSON.parse(String(req.body))
    else
        JSON.parsefile(filename)
    end
end

function validate(s::Schema, strict::Bool=false)
    if isempty(s.descriptor)
        push!(s.errors, SchemaError("Missing Descriptor"))
    end
    if length(s.fields) == 0
        push!(s.errors, SchemaError("No Fields specified"))
    end

    # Validate each Field
    for fld in s.fields
        try
            validate(fld)
        catch ex
            if isa(ex, FieldError)
                push!(s.errors, SchemaError(ex))
            else
                throw(ex)
            end
        end
    end

    # Validate primary keys
    for key in s.primary_key
        if !(has_field(s, key))
            push!(s.errors, SchemaError("Missing field as defined in primaryKey", key))
        end
    end
    for e in s.errors
        if occursin(e.message, "primaryKey") && e.key != ""
            deleteat!(s.primary_key, findall(in([e.key]), s.primary_key))
        end
    end

    # Validate foreign keys
    for key in s.foreign_keys
        if !(haskey(key, "fields"))
            push!(s.errors, SchemaError("'fields' is required on all foreignKeys"))
            key["fields"] = []
        end
        if !(haskey(key, "reference"))
            push!(s.errors, SchemaError("'reference' is required on all foreignKeys"))
            key["reference"] = Dict()
        end
        if !(haskey(key["reference"], "resource"))
            push!(s.errors, SchemaError("'resource' is required on each reference of a foreignKey"))
            key["reference"]["resource"] = ""
        end
        if !(haskey(key["reference"], "fields"))
            push!(s.errors, SchemaError("'fields' is required on each reference of a foreignKey"))
            key["reference"]["fields"] = []
        end
        if isa(key["fields"], String)
            key["fields"] = [key["fields"]]
        end
        if isa(key["reference"]["fields"], String)
            key["reference"]["fields"] = [key["reference"]["fields"]]
        end
        # Validate field references
        for f in key["fields"]
            if !(has_field(s, f))
                push!(s.errors, SchemaError("Missing field as defined in foreignKeys fields", f))
            end
        end
        if key["reference"]["resource"] == ""
            for f in key["reference"]["fields"]
                if !(has_field(s, f))
                    push!(s.errors, SchemaError("Missing field as defined in foreignKeys references", f))
                end
            end
        end
        if length(key["fields"]) != length(key["reference"]["fields"])
            push!(s.errors, SchemaError("Number of source and outer foreign key references must match"))
        end
        # Handle errors
        for e in s.errors
            if occursin(e.message, "foreignKeys fields") && e.key != ""
                deleteat!(key["fields"],
                    findall(in([e.key]), key["fields"]))
            end
            if occursin(e.message, "foreignKeys references") && e.key != ""
                deleteat!(key["reference"]["fields"],
                    findall(in([e.key]), key["reference"]["fields"]))
            end
        end
    end

    # Error handling
    if strict && length(s.errors)>0
        throw(s.errors[1])
    end

    # TODO: find out why returning something else breaks the module
    s
end

function guess_type(value)
    if isa(value, Int)
        return "integer"
    elseif isa(value, Number)
        return "number"
    elseif isa(value, Bool)
        return "boolean"
    elseif isa(value, AbstractString)
        try
            df = DateFormat("YYYY-MM-DDThh:mm:ssZ")
            dd = Date(value, df)
            return "datetime"
        catch
            # continue
        end
        try
            df = DateFormat("HH:MM:SS")
            dd = Date(value, df)
            return "time"
        catch
            # continue
        end
        try
            df = DateFormat("y-m-d")
            dd = Date(value, df)
            return "date"
        catch
            # continue
        end
        gp = split(value, ",")
        if length(gp) == 2
            try
                lon = parse(Float64, gp[1])
                lat = parse(Float64, gp[2])
                if lon <= 180 && lon >= -180 &&
                    lat <= 90 && lat >= -90
                        return "geopoint"
                end
            catch
                # continue
            end
        end
        try
            obj = JSON.parse(value)
            if isa(obj, Array)
                return "array"
            elseif isa(obj, Dict)
                return "object"
            end
        catch
            # continue
        end
        return "string"
    else
        return nothing
    end
end

function infer(s::Schema, rows::Array{Any}, headers::Array{String}, maxrows=MAX_ROWS_INFER)
    for (c, header) in enumerate(headers)
        if has_field(s, header)
            continue
        end
        type_match = Dict()
        col = view(rows, :, c)
        if maxrows > length(col)
            maxrows = length(col)
        end
        for (r, val) in enumerate(col[1:maxrows])
            guess = guess_type(val)
            if guess == nothing
                @warn "Could not guess type at ($(r), $(c))"
            else
                if !haskey(type_match, guess)
                    type_match[guess] = 0
                end
                type_match[guess] = type_match[guess] + 1
                @debug "Guessed $(guess) at ($(r), $(c))"
            end
            # print(val)
            # print("\n")
        end
        f = Field(header)
        if length(type_match)>0
            sorted = sort(collect(type_match), by=x->x[2], rev=true)
            f.typed = sorted[1][1]
        end
        add_field(s, f)
    end
end

function set_primary_key(s::Schema, name::String)
    if !(has_field(s, name))
        throw(SchemaError("Fields defined as Primary Key must exist"))
    end
    if !(name in s.primary_key)
        push!(s.primary_key, name)
    end
end

function set_foreign_key(s::Schema, source_fields::Array{String}, resource::String, outer_fields::Array{String})
    for key in source_fields
        if !(has_field(s, key))
            throw(SchemaError("Source fields must reference existing schema fields"))
        end
    end
    if resource == ""
        for key in outer_fields
            if !(has_field(s, key))
                throw(SchemaError("Outer fields must reference existing schema fields to set empty reference"))
            end
        end
    end
    push!(s.foreign_keys, Dict(
            "fields" => source_fields,
            "reference" => Dict(
                "resource" => resource,
                "fields" => outer_fields
            )
        )
    )
end

function cast_row(s::Schema, row::Array, fail_fast=false, check_constraints=true)
    result = []; errors = []

    lrow = length(row); lsf = length(s.fields)
    (lrow != lsf) && throw(CastError(
        "Row length $lrow does not match fields count $lsf"))

    for (i, field) in enumerate(s.fields)
        value = row[i]
        try
            push!(result, cast_value(field, value, check_constraints))
        catch e
            if isa(e, CastError) && !fail_fast
                push!(errors, e)
            else
                throw(e)
            end
        end
    end

    le = length(errors)
    (le>0) && throw(CastError(
        "There are $le cast errors (see exception.errors)", errors))

    result
end

function build(s::Schema)
    d = Dict()
    d["fields"] = [ build(f) for f in s.fields ]
    d["primaryKey"] = s.primary_key
    d["foreignKeys"] = s.foreign_keys
    d["missingValues"] = s.missing_values
    s.descriptor = d
    d
end

function commit(s::Schema ; strict=nothing)
    strict !== nothing &&
        throw(ErrorException("strict parameter not implemented"))
    @debug "Building schema descriptor"
    s.descriptor = build(s)
    s.errors = []
    validate(s)
    length(s.errors) > 0 &&
        throw(ErrorException("Invalid schema: check errors"))
    true
end

function save(s::Schema, target::String)
    commit(s)
	indent = 4
	@debug "Saving table schema to $target"
	open(target, "w") do io
		JSON.print(io, s.descriptor, indent)
	end
end

field_names(s::Schema) = [ f.descriptor.name for f in s.fields ]
get_field(s::Schema, name::String) = [ f for f in s.fields if f.name == name ][1]
get_field_index(s::Schema, name::String) = findall(in([get_field(s, name)]), s.fields)
has_field(s::Schema, name::String) = length([ true for f in s.fields if f.name == name ]) > 0
add_field(s::Schema, d::Dict) = push!(s.fields, Field(d))
add_field(s::Schema, f::Field) = push!(s.fields, f)
remove_field(s::Schema, name::String) = deleteat!(s.fields, get_field_index(s, name))

is_empty(s::Schema) = Base.isempty(s.fields)
is_valid(s::Schema) = (!Base.isempty(s.descriptor) && length(s.errors) == 0)

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
using Random
using Statistics
using StatsBase
using Plots
import Base.+
import Base.-

include("Point_Struct.jl")
include("fun.jl")
include("constants.jl")

# This script reproduce the turbulance distribution function (a gaussian)
# Figure 4 in Errard's article

p_0 = Point{Float64}(4000, 0.8, 0.8, undef, undef, undef)()
p_1 = [Point{Float64}(4000, 0.8, i, undef, undef, undef)() for i = (1:1000)*( (π/2) / 1000 )  ]

dr_vec = [ diff_cart(i,p_0) for i in p_1  ]

DR = [( i.x^2 + i.y^2 + i.z^2 )^0.5  for i in dr_vec ]

w = [ χ₁(p_0, i) for i in p_1 ]

plot(DR, w, seriestype=:scatter, xlims = (0,1000))

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
using BenchmarkTools, OMEinsum, CUDA
CUDA.allowscalar(false)

a = randn(Float32, 50, 50)
@benchmark ein"ij,jk->ik"($a,$a) seconds = 1
@benchmark ein"ij,ik,il->jkl"($a,$a,$a) seconds = 1
@benchmark a*a seconds = 1

a = a |> CuArray
y = y |> CuArray
@benchmark ein"ij,jk->ik"($a,$a) seconds = 1
@benchmark ein"ij,ik,il->jkl"($a,$a,$a) seconds = 1
@benchmark a*a seconds = 1


text:
# mesh tools

meshgrid(v::AbstractVector) = meshgrid(v, v)

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

# Construct mesh for quick reconstruction
function buildmesh!(dc::DepthCamera)
  H, W = dc.shape
  xs,ys = collect(1:W), collect(1:H)
  fxinv = 1.0 / dc.K[1,1];
  fyinv = 1.0 / dc.K[2,2];

  xs = (xs-dc.K[1,3]) * fxinv
  xs = xs[1:dc.skip:end]
  ys = (ys-dc.K[2,3]) * fyinv
  ys = ys[1:dc.skip:end]

  dc.xs, dc.ys = meshgrid(xs, ys);
  nothing
end


text:
function steiner(n)
    model = Model(optimizer_with_attributes(CS.Optimizer,
        "logging" => [],
        "seed" => 4,

        "traverse_strategy" => :BFS,
        # "traverse_strategy"=>:DFS,
        # "traverse_strategy"=>:DBFS,

        # "branch_split"=>:Smallest,
        # "branch_split"=>:Biggest,
        "branch_split" => :InHalf,

        # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
        # "branch_strategy" => :IMPS, # default
        "branch_strategy" => :ABS, # Activity Based Search
        "activity.decay" => 0.999, # default 0.999
        "activity.max_probes" => 1, # default, 10
        "activity.max_confidence_deviation" => 20, # default 20

        # "simplify"=>false,
        # "simplify"=>true, # default

        # "backtrack" => false, # default true
        # "backtrack_sorting" => false, # default true

        # "lp_optimizer" => cbc_optimizer,
        # "lp_optimizer" => glpk_optimizer,
        # "lp_optimizer" => ipopt_optimizer,
    ))

    @assert (n % 6 == 1 || n % 6 == 3)

    nb = round(Int, (n * (n - 1)) / 6) # number of sets

    @variable(model, x[1:nb,1:n], Bin)
    @constraint(model, x[1,1] == 1) # symmetry breaking

    # atmost 1 element in common
    for i in 1:nb
        @constraint(model,sum(x[i,:]) == 3)

        for j in i + 1:nb
            b = @variable(model, [1:n], Bin)
            for k in 1:n
                @constraint(model, b[k] := { x[i,k] && x[j,k] })
            end
            @constraint(model, sum(b) <= 1)
        end
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    @assert status == MOI.OPTIMAL
end


text:
# # Fitting and Predicting - Time Data
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/fitpredictTime.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/fitpredictTime.ipynb)

# A large majority of scientists and engineers who undertake rheological experiments fit their data with one or several viscoelastic models in order to classify materials, quantify their behaviour and predict their response to external perturbations.

# Standard linear viscoelastic models take the form of an ordinary differential equation between stress σ and strain ϵ. Fitting models and predicting their response in the time domain then requires computing viscoelastic hereditary integrals such as:

# ```math
# \sigma(t) = \int_{0}^t G(t - \tau) \frac{d \epsilon(\tau)}{d \tau} d \tau
# ```

# where G is the relaxation response (stress response to a step in strain) of the material, which can be analytically or numerically defined for each rheological model and its set of rheological parameters. A similar relation exists to calculate the strain from the stress history:

# ```math
# \epsilon(t) = \int_{0}^t J(t - \tau) \frac{d \sigma(\tau)}{d \tau} d \tau
# ```

# where J is the creep modulus of the material. Fitting and predicting behaviour then becomes non-trivial and standardised tools to process are needed.

using RHEOS

# ## Fitting

# #### Step 1: Loading in Experimental Data

# RHEOS has a convenient function to import data from CSV files (more information in [File I/O](@ref) section).

data = importcsv("assets/data_time.csv", t_col = 1, ϵ_col = 2, σ_col = 3)
rheotimedatatype(data)

# The time, strain and stress data are stored into a [`RheoTimeData`](@ref) struct.

# #### Step 2: Fitting Models

# Let's fit a Maxwell model via its relaxation modulus, G, as our data is from a stress relaxation test, thus the strain is imposed. The first argument is our data, the second argument tells RHEOS which model to fit ([`RheoModelClass`](@ref)) and the final argument tells RHEOS whether to fit the model using a relaxation modulus (*strain* imposed) or creep modulus (*stress* imposed). The function will return a [`RheoModel`](@ref) object, i.e. a [`RheoModelClass`](@ref) with fixed values for the parameters.

maxwell_model = modelfit(data, Maxwell, strain_imposed)

# For the fitting process RHEOS relies on the optimistion package [NLopt.jl](https://nlopt.readthedocs.io/en/latest/). RHEOS makes use of a local derivative free algorithm, specifically the Tom Rowan's "Subplex" algorithm originally introduced in his [PhD thesis at The University of Texas at Austin, 1990](https://dl.acm.org/doi/book/10.5555/100816).

# Next, we will fit a fractional Maxwell (spring) model (the only difference from the above model is that the dash-pot is replaced by a spring-pot). This time we will also add upper and lower bounds on the model parameters. This is highly recommended for fractional models in particular as values less than 0 or greater than 1 for the spring-pot parameter are unphysical and can cause errors in the Mittag-Leffler function used (more information about the fractional models in review TBA). To facilitate the convergence of the fitting algorithm, it is also suggested to define an initial set of parameters (the function uses 0.5 for all parameters not defined). For initial parameters, lower bounds and upper bounds, it is possible to only specify some of the parameters and the rest will be filled in automatically.

## Lower bounds
Lo = (cₐ = 0.0, a = 0.01, k = 0.0)
## Upper bounds
Hi = (cₐ = Inf, a = 0.99, k = Inf)
## Initial parameters
P0 = (cₐ = 0.1, a = 0.2, k = 2.0)

Fmaxwell_model = modelfit(data, FractS_Maxwell, strain_imposed, lo = Lo, hi = Hi, p0 = P0)

# Further optional parameters can be defined in the fitting function, such as the relative tolerance (default parameter of 1e-4) or the finite difference formula to use for the derivative (`"BD"` = Backward Difference (default), `"CD"` = Central Difference). If the user is interested in the model parameters on each optimisation iteration, they can be printed on the terminal by enabling `verbose = true`.

#md # !!! warning "Important"
#md #     Often an ideal step loading is assumed, which simplifies the equations in the first section to σ(t) = G(t)ϵ₀ and ϵ(t) = J(t)σ₀. This enables the modulus to be used directly for the fitting, which vastly reduces the computational burden. For this reason, RHEOS provides a dedicated function [`modelstepfit`](@ref) that requires the same parameters as [`modelfit`](@ref). If this assumption is appropriate for the data then fitting can be sped up greatly by use of this function.

#nb # Important Note: Often an ideal step loading is assumed, which simplifies the equations in the first section to σ(t) = G(t)ϵ₀ and ϵ(t) = J(t)σ₀. This enables the modulus to be used directly for the fitting, which vastly reduces the computational burden. For this reason, RHEOS provides a dedicated function `modelstepfit` that requires the same parameters as `modelfit`. If this assumption is appropriate for the data then fitting can be sped up greatly by use of this function.

# ## Predicting

# RHEOS allows the user to simulate the response of a model (with defined parameters, [`RheoModel`](@ref) struct) to an imposed loading.

# Given an incomplete data set (time data with only either stress or strain missing) and model with values substituted into parameters ([`RheoModel`](@ref)), return a new "complete" dataset based on the model with the simulated missing variable.

# #### Assess Quality of the Fit

# The ability of predicting model's response is first exploited to assess the quality of the fits above. The incomplete [`RheoTimeData`](@ref) variable is defined by extracting the time and the imposed variable (for the current example, the strain).

data_ext = onlystrain(data)
## alternatively: data_ext = extract(data, 1);

rheotimedatatype(data_ext)

# For the prediction, RHEOS' function [`modelpredict`](@ref) requires the incomplete data set and a model with fixed parameters ([`RheoModel`](@ref)). For the assessment of the fitting quality the [`RheoModel`](@ref) required is the output of the fitting function.

maxwell_predict = modelpredict(data_ext, maxwell_model)

## Now we can plot data and model together for comparison
using PyPlot

fig, ax = subplots(1, 1, figsize = (7, 5))
ax.plot(data.t, data.σ, "o", markersize = 5)
ax.plot(maxwell_predict.t, maxwell_predict.σ, color = "red")
ax.set_xlabel("Time")
ax.set_ylabel("Stress")
#!nb fig #hide

# #### Simulate different model behaviours

# The ability of predicting model's response can be exploited to simulate the behaviour of the model to other external loading conditions. For example, we will explore the response of the fitted model to a creep experiment.

# First we need to define an incomplete [`RheoTimeData`](@ref) struct, which can be achieved via the data generation functions provided in RHEOS (see section [Generating Data](@ref)).

## Define timeline
dσ = timeline(t_end = 10)
## and a step in stress
dσ = stressfunction(dσ, hstep())

## we can now predict the creep response of the Maxwell model
maxwell_creepPredict = modelsteppredict(dσ, maxwell_model)

## Init plot figure
fig, ax = subplots(1, 1, figsize = (7, 5))

## Visualisation of the imposed loading
ax.plot(maxwell_creepPredict.t, maxwell_creepPredict.σ, "--")
## Visualisation of the simulated response
ax.plot(maxwell_creepPredict.t, maxwell_creepPredict.ϵ)
ax.set_xlabel("Time")
ax.set_ylabel("Strain")
#!nb fig #hide

#md # !!! warning "Important"
#md #     As for the fitting procedure, if an ideal step loading is assumed, RHEOS provides a dedicated function [`modelsteppredict`](@ref) that requires the same parameters as [`modelpredict`](@ref) which significantly reduces the computation time.

#nb # Important Note: As for the fitting procedure, if an ideal step loading is assumed, RHEOS provides a dedicated function modelsteppredict that requires the same parameters as modelpredict which significantly reduces the computating time.

# Similarly, more complex behaviour can be simulated (such as a stairs loading).

## Define timeline
dσ = timeline(t_end = 10)
## and a step in stress
dσ = stressfunction(dσ, stairs())

## we can now predict the stairs response of the Maxwell model
maxwell_creepPredict = modelpredict(dσ, maxwell_model)
## Plotting
fig, ax = subplots(1, 1, figsize = (7, 5))
ax.plot(maxwell_creepPredict.t, maxwell_creepPredict.σ, "--")
ax.plot(maxwell_creepPredict.t, maxwell_creepPredict.ϵ)
ax.set_xlabel("Time")
ax.set_ylabel("Strain")
#!nb fig #hide


text:
#!/usr/bin/env julia

export pe6

pe6(n::Integer=100)=sum(1:n)^2-sum((1:n).^2)

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe6()|>print
end

max_stars_repo_path: 6. Sum square difference/6.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 6. Sum square difference/6.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 6. Sum square difference/6.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
export NoBoundaries

struct NoBoundaries{N} <: BoundaryComponent{N} end

Base.show(io::IO, ::Type{NoBoundaries{N}}) where {N} = print(io, "NoBoundaries")

get_positions(::NoBoundaries{N}) where {N} = SVectorF{N}[]
get_velocities(::NoBoundaries{N}) where {N} = SVectorF{N}[]
get_pressures(::NoBoundaries) = Float[]

function compute_mass_density(
    ::NoBoundaries,
    ::EnergyComponent,
    ::EquationOfState,
    ::Unsigned,
    ::Unsigned,
)
    nothing
end

function updateboundaries!(
    ::NoBoundaries,
    ::AbstractVector{SVectorF{N}},
    ::Velocities{N},
    ::MassComponent{MD,<:Kernel{N}},
    ::Pressures,
) where {N,MD} end

function evolveboundaries!(::NoBoundaries, ::Number) end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
@testset "examples" begin

function find_sources(path::String, sources=String[])
    if isdir(path)
        for entry in readdir(path)
            find_sources(joinpath(path, entry), sources)
        end
    elseif endswith(path, ".jl")
        push!(sources, path)
    end
    sources
end

examples_dir = joinpath(@__DIR__, "..", "examples")
examples = find_sources(examples_dir)
filter!(file -> readline(file) != "# EXCLUDE FROM TESTING", examples)

cd(examples_dir) do
    examples = relpath.(examples, Ref(examples_dir))
    @testset for example in examples
        # construct a command
        cmd = `$(Base.julia_cmd())`
        if Base.JLOptions().project != C_NULL
            # --project isn't preserved by julia_cmd()
            cmd = `$cmd --project=$(unsafe_string(Base.JLOptions().project))`
        end
        cmd = `$cmd $example`

        # run it and conditionally show output
        out = Pipe()
        rv = success(pipeline(cmd, stdout=out, stderr=out))
        close(out.in)
        rv || print(read(out, String))
        @test rv
    end
end

end


text:
#################################################
"""
```julia
randHermitian(d, n; diag, norm )

randHermitian(n; norm)
```
- `d` : entry distribution
- `n`  : dimensions
- `norm` : default `false`, if `norm` set to `true`, then the matrix will be normalized with ``n^{-1/2}``.
- `diag` : default `diag = d`, diagonal entry distribution.
    To use a different distribution (say Circular(2)) for digonal elements, set `diag = Circular(2)`.
    The diagonal entries will always be forced to have imgainary part `0`.
- See also [`GUE`](@ref)

# Examples
Generates a 2 by 2 random Hermitian matrix with off-diagonal entries from the Standard Complex Gaussian, and Standard Normal on the diagonal.
```julia
randHermitian(2)

2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
  0.382095+0.0im        -0.708469-0.0636734im
 -0.708469+0.0636734im   0.336952+0.0im
```
Generate a 3 by 3 Hermitian matrix, with off-diagonal entries [`Circular`](@ref)(1) and diagonal entries uniformly `-1` or `1`.
```julia
randHermitian(Circular(1),3,diag = (-1,1))

3×3 Hermitian{ComplexF64, Matrix{ComplexF64}}:
     1.0+0.0im       1.56259-0.676099im  1.39468-0.295073im
 1.56259+0.676099im     -1.0+0.0im       1.53369+0.296583im
 1.39468+0.295073im  1.53369-0.296583im     -1.0+0.0im
```
Generate a random 2 by 2 Symmetric Matrix with entries  `Poisson(2)` rvs.
This is also be done with `randSymmetric(Poisson(2),3)`
```julia
randHermitian(Poisson(2),3)

3×3 Hermitian{Int64, Matrix{Int64}}:
 3  1  0
 1  1  2
 0  2  1
```
Entries uniformly from ``\\{1,2,3,...,10\\}``
```julia
randHermitian(1:10,2)

2×2 Hermitian{Int64, Matrix{Int64}}:
 10  7
  7  6
```

Entries either -1 or pi with equal probability
```julia
randHermitian([-1,pi],3)

3×3 Hermitian{Float64, Matrix{Float64}}:
  3.14159  -1.0      -1.0
 -1.0      -1.0       3.14159
 -1.0       3.14159   3.14159
```
"""
function randHermitian(d::D, n::Int; diag=d::D, norm = false::Bool)  where D<:S

    M = randTriangular(d,n,diag=diag)

    if norm
        M/=sqrt(n)
    end

    return Hermitian(convert(Matrix,M))
end

function randHermitian(n::Int; norm = false::Bool)
    return randHermitian(ComplexNormal(), n, diag=Normal(), norm = norm)
end

##########################################################
"""
```julia
randSymmetric(d, n; Diag, norm)

randSymmetric(n; norm)
```
- `d` : entry distribution
- `n` : dimensions
- `norm` : default  `false`, if `norm` set to `true`, then the matrix will be normalized with ``n^{-1/2}``.
- `diag` : default `diag = d`, the distribution for diagonal entries.
    To use a different distribution (say Binomial) for digonal elements, set `diag = Binomial(1,0.5)`
- See also [`GOE`](@ref)

# Examples
Generates a 3 by 3 random Symmetric matrix with entries from the Standard Gaussian.
```julia
randSymmetric(3)

3×3 Symmetric{Float64, Matrix{Float64}}:
 -0.230698  -1.72846     0.306362
 -1.72846    0.0845915  -0.0116108
  0.306362  -0.0116108  -0.559046
```
"""
function  randSymmetric(d::D, n::Int; diag = d::D, norm = false::Bool)  where D<:S

    M = randTriangular(d,n,diag=diag)

    if norm
        M/=sqrt(n)
    end

    return Symmetric(convert(Matrix,M))
end

function randSymmetric(n::Int; norm = false::Bool)
    return randSymmetric(Normal(), n, norm = norm)
end

#########################################################

"""
```julia
GOE <: ContinuousMatrixDistribution
GOE(n)
```
- `n` : dimension
- The **Gaussian Orthogonal Ensemble** (`GOE`) is an ensemble of random ``n \\times n`` Symmetric matrices
    ``M_{n}`` in which the upper-triangular entries are iid with distribution ``N(0,1)_{\\mathbf{R}}``,
    and the diagonal entries are iid with distribution ``N(0,2)_{\\mathbf{R}}``, and independent of the upper-triangular ones

```julia
rand(M::GOE, norm::bool)
```
 -  `norm` : default `false`,  if `norm` set to `true`, then the matrix will be normlaized with ``\\operatorname{min}(n,m)^{-1/2}``.

# Examples

Generate a 3 by 3 random matrix from GOE(3)
```julia
rand(GOE(3))

3×3 Symmetric{Float64, Matrix{Float64}}:
 -1.62391   -0.451433    0.863883
 -0.451433   0.0271799  -0.524854
  0.863883  -0.524854   -0.00930624
```

```julia
rand(GOE(3),norm=true)

3×3 Symmetric{Float64, Matrix{Float64}}:
  0.302141   0.152634   -0.711679
  0.152634  -0.0629327   0.103075
 -0.711679   0.103075    1.51861
```

"""
struct GOE <: ContinuousMatrixDistribution
    n::Int
end

rand( M::GOE; norm = false::Bool) = randSymmetric(Normal(),M.n,diag=Normal(0,sqrt(2)),norm=norm)

"""
```julia
GUE <: ContinuousMatrixDistribution
GUE(n)
```
- `n` : dimension
- The **Gaussian Unitary Ensemble** (`GUE`) is an ensemble of random ``n \\times n`` Hermitian matrices
    ``M_{n}`` in which the upper-triangular entries are iid with distribution ``N(0,1)_{\\mathbf{C}}``,
    and the diagonal entries are iid with distribution ``N(0,1)_{\\mathbf{R}}``, and independent of the upper-triangular ones

```julia
rand(M::GUE, norm::bool)
```
- `norm` : default `false`,  if `norm` set to `true`, then the matrix will be normlaized with ``\\operatorname{min}(n,m)^{-1/2}``.

# Examples

Generate a 3 by 3 random matrix from GUE(3)
```julia
rand(GUE(3))

3×3 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 -0.883413+0.0im         1.09872+0.874884im     -0.1985-1.04778im
   1.09872-0.874884im    1.55483+0.0im        -0.488532+1.18694im
   -0.1985+1.04778im   -0.488532-1.18694im   -0.0823873+0.0im
```

```julia
rand(GUE(2),norm=true)
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 -0.457089+0.0im       -0.672713-0.102234im
 -0.672713+0.102234im   0.380126+0.0im
```
"""
struct GUE <: ContinuousMatrixDistribution
    n::Int
end

rand( M::GUE; norm=false::Bool) = randHermitian(ComplexNormal(),M.n,diag=Normal(),norm=norm)

max_stars_repo_path: src/MatrixModels/Hermitian Matrices.jl

max_issues_repo_path: src/MatrixModels/Hermitian Matrices.jl

max_forks_repo_path: src/MatrixModels/Hermitian Matrices.jl


text:
using Blosc2: Compressor, fill_available_compressors!, available_compressors, available_compressors_names, has_compressor, default_compressor
@testset "compressors" begin
    cprs = Dict{Symbol,Compressor}()
    fill_available_compressors!(cprs)
    @test haskey(cprs, :blosclz)
    @test haskey(available_compressors(), :blosclz)
    @test has_compressor(:blosclz)

    @test Blosc2.default_compressor_name() == :blosclz

    @test Blosc2.compressor_by_name(:blosclz) == Blosc2.default_compressor()
end

@testset "filters" begin
    @test has_filter(:shuffle)
    @test !has_filter(:shuffle2)

    f = Blosc2.make_filter(:nofilter)
    @test f.id == 0
    @test f.meta == 0

    f = Blosc2.make_filter(:delta, 10)
    @test f.id == Blosc2.Lib.BLOSC_DELTA
    @test f.meta == 10

    @test :shuffle in available_filter_names()

    @test filter_description(:nofilter) == "No filter."

    pipe = filter_pipeline(:nofilter, :trunc_prec=>22, :shuffle)
    @test pipe.filters[1].name == :nofilter
    @test pipe.filters[1].id == 0
    @test pipe.filters[1].meta == 0

    @test pipe.filters[2].name == :trunc_prec
    @test pipe.filters[2].id == Blosc2.Lib.BLOSC_TRUNC_PREC
    @test pipe.filters[2].meta == 22

    @test pipe.filters[3].name == :shuffle
    @test pipe.filters[3].id == Blosc2.Lib.BLOSC_SHUFFLE
    @test pipe.filters[3].meta == 0

    for i in 4:Blosc2.max_filters_count()
        @test pipe.filters[i].name == :nofilter
        @test pipe.filters[i].id == 0
        @test pipe.filters[i].meta == 0

    end
    pipe = filter_pipeline()
    for i in 1:Blosc2.max_filters_count()
        @test pipe.filters[i].name == :nofilter
        @test pipe.filters[i].id == 0
        @test pipe.filters[i].meta == 0

    end

end

@testset "params" begin
    params = Blosc2.CompressionParams()
    @test params.compressor.code == default_compressor().code
    params = Blosc2.CompressionParams(Int32)
    @test params.typesize == sizeof(Int32)
    @test params.compressor.name == :blosclz

    params = Blosc2.CompressionParams(Int8; compressor = :lz4, level = 9, nthreads = 3, blocksize = 10, splitmode = true)
    @test params.typesize == sizeof(Int8)
    @test params.compressor.name == :lz4
    @test params.level == 9
    @test params.nthreads == 3
    @test params.blocksize == 10
    @test params.splitmode == true

    cp = Blosc2.make_cparams(params)
    @test cp.typesize == sizeof(Int8)
    @test cp.compcode == 1
    @test cp.clevel == 9
    @test cp.nthreads == 3
    @test cp.blocksize == 10
    @test cp.splitmode == 1
    @test cp.filters == (Blosc2.Lib.BLOSC_SHUFFLE, 0, 0, 0, 0, 0)
    @test cp.filters_meta == (0, 0, 0, 0, 0, 0)

    params = Blosc2.CompressionParams(Int8; filter_pipeline = filter_pipeline(:trunc_prec=>22, :shuffle))
    cp = Blosc2.make_cparams(params)
    @test cp.filters == (Blosc2.Lib.BLOSC_TRUNC_PREC, Blosc2.Lib.BLOSC_SHUFFLE, 0, 0, 0, 0)
    @test cp.filters_meta == (22, 0, 0, 0, 0, 0)

    dparams = Blosc2.DecompressionParams(nthreads = 10)
    dp = Blosc2.make_dparams(dparams)
    @test dp.nthreads == 10
end

@testset "context" begin
    ctx = Blosc2.Context(Blosc2.CompressionParams(level = 9))
    @test ctx.params.nthreads == 1
    @test ctx.params.level == 9
    @test typeof(ctx) == Blosc2.CompressionContext
    @test typeof(ctx.params) == Blosc2.CompressionParams

    ctx = Blosc2.CompressionContext(Blosc2.CompressionParams(level = 8))
    @test ctx.params.nthreads == 1
    @test ctx.params.level == 8
    @test typeof(ctx) == Blosc2.CompressionContext
    @test typeof(ctx.params) == Blosc2.CompressionParams

    ctx = Blosc2.CompressionContext(level = 8, nthreads = 2)
    @test ctx.params.nthreads == 2
    @test ctx.params.level == 8
    @test typeof(ctx) == Blosc2.CompressionContext
    @test typeof(ctx.params) == Blosc2.CompressionParams

    ctx = Blosc2.CompressionContext()
    @test ctx.params == Blosc2.CompressionParams()
    @test typeof(ctx.params) == Blosc2.CompressionParams

    ctx = Blosc2.Context(Blosc2.DecompressionParams(nthreads = 2))
    @test ctx.params.nthreads == 2
    @test typeof(ctx) == Blosc2.DecompressionContext
    @test typeof(ctx.params) == Blosc2.DecompressionParams

    ctx = Blosc2.DecompressionContext(Blosc2.DecompressionParams(nthreads = 3))
    @test ctx.params.nthreads == 3
    @test typeof(ctx) == Blosc2.DecompressionContext
    @test typeof(ctx.params) == Blosc2.DecompressionParams

    ctx = Blosc2.DecompressionContext(nthreads = 3)
    @test ctx.params.nthreads == 3
    @test typeof(ctx) == Blosc2.DecompressionContext
    @test typeof(ctx.params) == Blosc2.DecompressionParams

    ctx = Blosc2.DecompressionContext()
    @test typeof(ctx.params) == Blosc2.DecompressionParams
    @test ctx.params == Blosc2.DecompressionParams()

    ctx = CompressionContext(Int32)
    @test ctx.params.typesize == 4
    ctx = CompressionContext(Int32, typesize = 2)
    @test ctx.params.typesize == 2
end

@testset "compress!/decompress!/info" begin
    n = 10000
    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    sz = Blosc2.compress!(Blosc2.CompressionContext(), buffer, data)
    @test Blosc2.compressed_sizeof(buffer) == sz
    @test Blosc2.uncompressed_sizeof(buffer) == n * sizeof(Int64)
    res = Vector{Int64}(undef, n)
    Blosc2.decompress!(Blosc2.DecompressionContext(),  res, buffer)
    @test res == data

    res = Vector{Int64}(undef, ceil(Int64, n/2))
    @test_throws BoundsError Blosc2.decompress!(Blosc2.DecompressionContext(),  res, buffer)

    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    sz = Blosc2.compress!(buffer, data)
    @test Blosc2.compressed_sizeof(buffer) == sz
    @test Blosc2.uncompressed_sizeof(buffer) == n * sizeof(Int64)
    res = Vector{Int64}(undef, n)
    Blosc2.decompress!(res, buffer)
    @test res == data

    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    sz = Blosc2.compress!(buffer, data; level = 0)
    @test Blosc2.compressed_sizeof(buffer) == sz
    @test sz == sizeof(buffer)
    @test Blosc2.uncompressed_sizeof(buffer) == n * sizeof(Int64)
    res = Vector{Int64}(undef, n)
    dsz = Blosc2.decompress!(res, buffer)
    @test dsz == n
    @test res == data
end

@testset "compress/decompress" begin
    n = 100000
    data = rand(1:1000, n)
    buffer = compress(CompressionContext(), data)
    @test sizeof(buffer) > 0
    @test sizeof(buffer) < max_compressed_size(data)
    res = decompress(DecompressionContext(), Int, buffer)
    @test res == data

    data = rand(Base.OneTo{Int16}(1000), n)
    buffer = compress(data, level = 9)
    @test sizeof(buffer) > 0
    @test sizeof(buffer) < max_compressed_size(data)
    res = decompress(Int16, buffer, nthreads = 2)
    @test res == data

    data = rand(Base.OneTo{Int16}(1000), n)
    buffer = compress(data, 51:60, level = 9)
    @test sizeof(buffer) > 0
    @test sizeof(buffer) < max_compressed_size(data)
    res = decompress(Int16, buffer, nthreads = 2)
    @test res == data[51:60]

end

@testset "unsafe compress/decompress" begin
    n = 10000
    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    sz = GC.@preserve data buffer begin
        unsafe_compress!(CompressionContext(), pointer(buffer), sizeof(buffer), pointer(data), sizeof(data))
    end
    @test sz > 0
    res = Vector{Int64}(undef, n)
    dsz = GC.@preserve data buffer begin
        unsafe_decompress!(DecompressionContext(), pointer(res), sizeof(res), pointer(buffer), sizeof(buffer))
    end
    @test dsz == n
    @test res == data

    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    sz = GC.@preserve data buffer begin
        unsafe_compress!(pointer(buffer), sizeof(buffer), pointer(data), sizeof(data), level = 1)
    end
    @test sz > 0
    res = Vector{Int64}(undef, n)
    dsz = GC.@preserve data buffer begin
        unsafe_decompress!(pointer(res), sizeof(res), pointer(buffer), sizeof(buffer), nthreads = 2)
    end
    @test dsz == n
    @test res == data
end

@testset "offsets compress/decompress" begin
    n = 10000
    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    pos = 1
    sz = compress!(CompressionContext(),buffer, data, pos, 1:500)
    @test sz > 0
    pos += sz
    sz = compress!(CompressionContext(),buffer, data, pos, 501:1000)
    @test sz > 0
    result = Vector{Int64}(undef, 1500)

    sz = decompress!(DecompressionContext(), result, buffer, 501, pos)
    @test sz == 500
    @test result[501:1000] == data[501:1000]

    sz = decompress!(DecompressionContext(), result, buffer, 1001, 1)
    @test sz == 500
    @test result[1001:1500] == data[1:500]

    n = 10000
    data = rand(Base.OneTo{Int16}(1000), n)
    buffer = Blosc2.make_compress_buffer(data)
    pos = 1
    sz = compress!(buffer, data, pos, 1:500)
    @test sz > 0
    pos += sz
    sz = compress!(buffer,  data, pos, 501:1000)
    @test sz > 0
    result = Vector{Int16}(undef, 1500)

    sz = decompress!(result, buffer, 501, pos)
    @test sz == 500
    @test result[501:1000] == data[501:1000]

    sz = decompress!(result, buffer, 1001, 1)
    @test sz == 500
    @test result[1001:1500] == data[1:500]
end

@testset "decompress_items!" begin
    n = 10000
    data = rand(1:1000, n)
    buffer = Blosc2.make_compress_buffer(data)
    pos = 1
    sz = compress!(CompressionContext(),buffer, data, pos, 1:500)
    @test sz > 0
    pos += sz
    sz = compress!(CompressionContext(),buffer, data, pos, 501:1000)
    @test sz > 0
    result = Vector{Int64}(undef, 1500)

    r = decompress_items!(DecompressionContext(), result,  buffer, 11:20, 50, pos)
    @test r == 10
    @test result[50:59] == data[511:520]

    r = decompress_items!(result,  buffer, 11:20)
    @test r == 10
    @test result[1:10] == data[11:20]

end

@testset "filters test" begin
    n = 10000
    data = rand(1:1000, n)
    buffer = compress(data)
    sz1 = sizeof(buffer)
    sdata = sort(data)
    sbuffer = compress(sdata, filter_pipeline = filter_pipeline(:delta, :shuffle))
    sz2 = sizeof(sbuffer)
    @test sz2 < sz1
    @test decompress(Int, sbuffer) == sdata
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

# messages to the supervisor
struct ChildInit{L,I}
    from::L
    init::I
end
struct Delete end
struct Terminate end
struct Which end
struct Strategy
    strategy::Symbol
end
struct NodeFailure{T}
    pids::T
end

const strategies = (:one_for_one, :one_for_all, :rest_for_one)
const restarts   = (:permanent, :temporary, :transient)
const sv_options = (:strategy, :max_restarts, :max_seconds, :spares)

isnormal(reason) = reason in (:normal, :shutdown, :done, :timed_out)

"""
```
Supervisor(; strategy=:one_for_one, max_restarts=3, max_seconds=5, kwargs...)
```
Supervisor functor with data and behavior.

# Fields (acquaintances)
- `option::Dict{Symbol,Any}`: supervisor option
- `childs::Array{Child,1}`: supervised childs,
- `rtime::Array{Float64,1}`: last restart times.

# Options
- `strategy::Symbol`: supervision strategy, can be
    either `:one_for_one`, `:one_for_all` or `:rest_for_one`.
- `max_restarts::Int`: maximum number of restarts allowed in a time frame,
- `max_seconds::Float64`: time frame in which `max_restarts` applies, defaults to 5,
- `kwargs...`: further option to extend supervisor behavior.
"""
struct Supervisor{O,C,T}
    option::O
    childs::C
    rtime::T
end
function Supervisor(; strategy=:one_for_one, max_restarts=3, max_seconds=5, kwargs...)
    s = Supervisor(Dict{Symbol,Any}(), Child[], Float64[])
    s.option[:strategy]     = strategy
    s.option[:max_restarts] = max_restarts
    s.option[:max_seconds]  = max_seconds
    merge!(s.option, kwargs)
    return s
end

function shutdown_child(c::Child)
    if c.lk isa Link
        try
            send(c.lk, Connect(self(), true))
            exit!(c.lk, :shutdown)
        catch exc
            log_error("shutdown_child", exc)
        end
    end
end

function restart_child!(c::Child, act::_ACT)
    if c.lk.chn isa RemoteChannel
        lk = !isnothing(c.start) ? c.start(act.bhv, pid=c.lk.pid) :
            !isnothing(act.init) ? spawn(act.init, pid=c.lk.pid) :
                spawn(act.bhv, pid=c.lk.pid)
    else
        lk = !isnothing(c.start) ? c.start(act.bhv) :
            !isnothing(act.init) ? spawn(act.init) :
                spawn(act.bhv)
    end
    c.lk.chn = lk.chn
    c.lk.pid = lk.pid
    send(lk, Connect(Super(self())))
    update!(lk, c.lk, s=:self)
end
function restart_child!(c::Child, ::Nothing)
    c.lk[] = Threads.@spawn c.start()
    _supervisetask(c.lk, self(), timeout=c.info.timeout, pollint=c.info.pollint)
end

function shutdown_restart_child!(c::Child)
    act = isnothing(c.start) ? diag(c.lk, :act) : nothing
    shutdown_child(c)
    restart_child!(c, act)
end

function must_restart(c::Child, reason)
    return c.info.restart == :permanent ? true :
        c.info.restart == :temporary ? false :
        isnormal(reason) ? false : true
end

function restart_limit!(s::Supervisor)
    length(s.rtime) > s.option[:max_restarts] && popfirst!(s.rtime)
    push!(s.rtime, time_ns() / 1e9)
    return length(s.rtime) > s.option[:max_restarts] &&
        s.rtime[end] - s.rtime[begin] ≤ s.option[:max_seconds]
end

function restart!(s::Supervisor, c::Child, msg::Exit)
    if s.option[:strategy] == :one_for_one
        log_warn("supervisor: restarting")
        restart_child!(c, msg.state)
    elseif s.option[:strategy] == :one_for_all
        log_warn("supervisor: restarting all")
        for child in s.childs
            child.lk == c.lk ?
                restart_child!(child, msg.state) :
                shutdown_restart_child!(child)
        end
    else
        log_warn("supervisor: restarting rest")
        ix = findfirst(x -> x.lk == c.lk, s.childs)
        for child in s.childs[ix:end]
            child.lk == c.lk ?
                restart_child!(child, msg.state) :
                shutdown_restart_child!(child)
        end
    end
end

#
# Supervisor behavior methods
#
function (s::Supervisor)(msg::ChildInit)
    ix = findfirst(c -> c.lk == msg.from, s.childs)
    isnothing(ix) && throw(AssertionError("child not found"))
    s.childs[ix].init = msg.init
end
function (s::Supervisor)(msg::Exit)
    ix = findfirst(c -> c.lk == msg.from, s.childs)
    isnothing(ix) && throw(AssertionError("child not found"))
    if must_restart(s.childs[ix], msg.reason)
        if restart_limit!(s)
            log_warn("supervisor: restart limit $(s.option[:max_restarts]) exceeded!")
            send(self(), Exit(:shutdown, fill(nothing, 3)...))
        else
            restart!(s, s.childs[ix], msg)
        end
    else
        ## todo: warn about temporary actor failure!
        act = task_local_storage("_ACT")
        filter!(c -> c.lk != msg.from, act.conn)
        filter!(c -> c.lk != msg.from, s.childs)
    end
end
function (s::Supervisor)(child::Child)
    act = task_local_storage("_ACT")
    ix = findfirst(c -> child.lk == c.lk, s.childs)
    isnothing(ix) && push!(s.childs, child)
    ix = findfirst(c -> child.lk == c.lk, act.conn)
    isnothing(ix) && push!(act.conn, child)
    if child.lk isa Link
        send(child.lk, Connect(Super(act.self)))
        myid() == child.lk.pid || rnfd_add(s, child.lk)
    end
end
function (s::Supervisor)(::Delete, lst)
    lk = lst isa Symbol ? whereis(lst) : lst
    act = task_local_storage("_ACT")
    filter!(c -> c.lk != lk, act.conn)
    filter!(c -> c.lk != lk, s.childs)
    if lk isa Link
        trysend(lk, Connect(act.self, true))
        lk.pid ≠ myid() &&  send(rnfd(s), Remove(lk))
    end
end
function (s::Supervisor)(::Terminate, ls)
    lk = ls isa Symbol ? whereis(ls) : ls
    s(Delete(), lk)
    try
        send(lk, Connect(self(), true))
        exit!(lk, :shutdown)
    catch
    end
end
(s::Supervisor)(::Which) = s.childs
(s::Supervisor)(msg::Strategy) = s.option[:strategy] = msg.strategy

#
# API functions
#
"""
```
supervisor( strategy=:one_for_one,
            max_restarts::Int=3,
            max_seconds::Real=5;
            name=nothing, kwargs...)
```
Spawn a supervisor actor with an empty child list and
return a link to it.

# Arguments
The following arguments are mandatory and go into a supervisor's
options:
- `strategy=:one_for_one`: [supervision strategy](@ref strategies), can be
    either `:one_for_one`, `:one_for_all` or `:rest_for_one`,
- `max_restarts::Int=3`: maximum number of restarts
    allowed in a time frame,
- `max_seconds::Real=5`: time frame in which
    `max_restarts` applies,

# Keyword arguments and further options
- `name=nothing`: name (Symbol) under which a supervisor should
    be registered, if `nothing` it doesn't get registered,
- `kwargs...`: further keyword arguments to the
    [`Supervisor`](@ref) and to [`spawn`](@ref). Keyword arguments
    not taken by `spawn` are supervisor options, used to extend a
    supervisor's functionality.

# Reserved options
- `strategy`, `max_restarts`, `max_seconds`, `name` are reserved,
    see above,
- `spares=[5,6,7]`: spare `pid`s, where you can give spare pids
    (e.g. `[5,6,7]`) to a supervisor used for actor restarts after
    node failures.
"""
function supervisor(strategy=:one_for_one, max_restarts::Int=3, max_seconds::Real=5; name=nothing, kwargs...)
    @assert strategy in strategies "Unknown strategy $strategy"
    sv_kw = merge((; strategy, max_restarts, max_seconds), structdiff(NamedTuple(kwargs), NamedTuple{spawn_kws}))
    lk = spawn(Supervisor(; sv_kw...); structdiff(NamedTuple(kwargs), sv_kw)...)
    !isnothing(name) && register(name, lk)
    send(lk, Update(:mode, :supervisor))
    return lk
end

"""
    count_children(sv)

Return a named tuple containing children counts for the given
supervisor `sv` (link or registered name).
"""
function count_children(sv)
    childs = call(sv, Which())
    d = Dict(:all => length(childs))
    ms = [c.lk.mode for c in childs if c.lk isa Link]
    tasks = length(childs) - length(ms)
    for m in Set(ms)
        d[m] = length(filter(==(m), ms))
    end
    tasks > 0 &&  (d[:tasks] = tasks)
    (; d...)
end

"""
```
start_actor(start, sv, cb=nothing, restart::Symbol=:transient;
            name::Union{Symbol,Nothing}=nothing, kwargs...)
```
Tell a supervisor `sv` to start an actor, to add
it to its childs list and to return a link to it.

# Arguments
- `start`: start behavior of the child, a callable object,
- `sv::Link`: link to a started supervisor,
- `cb=nothing`: callback (a callable object, gets the last
    actor behavior as argument and must return a [`Link`](@ref));
    if `nothing`, the actor gets restarted with its [`init!`](@ref)
    callback or with its last behavior,
- `restart::Symbol=:transient`: [restart option](@ref restart), one of
    `:permanent`, `:temporary`, `:transient`,
- `name::Union{Symbol,Nothing}=nothing`, name (Symbol) under which
    the actor should be registered,
- `kwargs...`: keyword arguments to [`spawn`](@ref).
"""
function start_actor(start, sv::Link, cb=nothing, restart::Symbol=:transient;
    name::Union{Symbol,Nothing}=nothing, kwargs...)
    @assert restart in restarts "Not a known restart strategy: $restart"
    lk = spawn(start; kwargs...)
    isnothing(name) || register(name, lk)
    send(sv, Child(lk, cb, lk.chn isa RemoteChannel ? start : nothing, name, (; restart)))
    return lk
end

function _supervisetask(r::Ref{Task}, sv; timeout::Real=5.0, pollint::Real=0.1)
    res = timedwait(() -> r[].state != :runnable, timeout; pollint)
    res == :ok ?
        r[].state == :done ?
            send(sv, Exit(:done, r, r[], nothing)) :
            send(sv, Exit(r[].exception, r, r[], nothing)) :
        send(sv, Exit(res, r, r[], nothing))
end

"""
```
start_task(start, sv, cb=nothing;
           timeout::Real=5.0, pollint::Real=0.1)
```
Spawn a task with a `start` behavior and tell the supervisor `sv` to supervise it (with
restart strategy `:transient`) and return a reference to it.

# Parameters
- `start`: must be a callable object (with no arguments),
- `sv`: link or registered name of a supervisor,
- `cb=nothing`: callback for restart (a callable object, must
    return a `Task`); if `cb=nothing`, the task gets restarted
    with its `start` function,
- `timeout::Real=5.0`: how long [seconds] should the task
    be supervised,
- `pollint::Real=0.1`: polling interval [seconds].
"""
function start_task(start, sv, cb=nothing; timeout::Real=5.0, pollint::Real=0.1)
    r = Ref(Threads.@spawn start())
    !isinf(timeout) && @async _supervisetask(r, sv; timeout, pollint)
    send(sv, Child(r, isnothing(cb) ? start : cb, start, nothing, (; restart=:transient, timeout, pollint)))
    return r
end

"""
    terminate_child(sv, ls)

Tell a supervisor `sv` to remove a child `ls` from its
childs list and to terminate it with reason `:shutdown`.

# Arguments
- `sv`: link or registered name of a supervisor,
- `ls`: link or registered name of a child.
"""
terminate_child(sv, ls) = call(sv, Terminate(), ls)

"""
    which_children(sv, info=false)

Tell a supervisor `sv` to return its childs list.

# Arguments
- `sv`: link or registered name of a supervisor,
- `info=false`: If `info=true`, return a list of named tuples
    with child information.
"""
function which_children(sv, info=false)
    function cinfo(c::Child)
        if c.lk isa Link
            i = Actors.info(c.lk)
            (actor = i.mode, bhv = i.bhvf, pid = i.pid, thrd = i.thrd, task = i.task, id = i.tid, name = i.name, restart = c.info.restart)
        else
            (task = c.lk[], restart = c.info.restart)
        end
    end
    childs = call(sv, Which())
    return info ? map(cinfo, childs) : childs
end

"""
```
supervise(sv, child=self(); cb=nothing, restart::Symbol=:transient)
```
Tell a supervisor `sv` to supervise the given `child` actor.

# Arguments
- `sv`: link or registered name of a supervisor,
- `child`: link or registered name of an actor to supervise.

# Keyword Arguments
- `cb=nothing`: callback (a callable object), takes the
    previous actor behavior as argument and must return
    a [`Link`](@ref) to a new actor; if `nothing`, the
    actor gets restarted with its [`init!`](@ref) callback
    or its previous behavior.
- `restart::Symbol=:transient`: [restart option](@ref restart), one of
    `:permanent`, `:temporary`, `:transient`,
"""
function supervise(sv, child=self(); cb=nothing, restart::Symbol=:transient)
    @assert restart in restarts "Not a known restart strategy: $restart"
    act = try
        a = task_local_storage("_ACT")
        a.self == child ? a : diag(child, :act)
    catch
        diag(child, :act)
    end
    send(sv, Child(child isa Symbol ? whereis(child) : child, cb, isnothing(act.init) ? act.bhv : act.init, act.name, (; restart)))
end

"""
    unsupervise(sv, lst=self())

Tell a supervisor `sv` to delete `ls` from its childs list.

# Arguments
- `sv`: link or registered name of a supervisor,
- `lst`: link, registered name of an actor or task.
"""
unsupervise(sv, lst=self()) = send(sv, Delete(), lst)

"""
    set_strategy(sv, strategy::Symbol)

Tell a supervisor `sv` to change its restart `strategy`.
"""
function set_strategy(sv, strategy::Symbol)
    @assert strategy in strategies "$strategy not known!"
    send(sv, Strategy(strategy))
end


text:
# Autogenerated wrapper script for mlpack_jll for x86_64-w64-mingw32-cxx11
export libmlpack, libmlpack_julia_adaboost, libmlpack_julia_approx_kfn, libmlpack_julia_bayesian_linear_regression, libmlpack_julia_cf, libmlpack_julia_dbscan, libmlpack_julia_decision_stump, libmlpack_julia_decision_tree, libmlpack_julia_det, libmlpack_julia_emst, libmlpack_julia_fastmks, libmlpack_julia_gmm_generate, libmlpack_julia_gmm_probability, libmlpack_julia_gmm_train, libmlpack_julia_hmm_generate, libmlpack_julia_hmm_loglik, libmlpack_julia_hmm_train, libmlpack_julia_hmm_viterbi, libmlpack_julia_hoeffding_tree, libmlpack_julia_image_converter, libmlpack_julia_kde, libmlpack_julia_kernel_pca, libmlpack_julia_kfn, libmlpack_julia_kmeans, libmlpack_julia_knn, libmlpack_julia_krann, libmlpack_julia_lars, libmlpack_julia_linear_regression, libmlpack_julia_linear_svm, libmlpack_julia_lmnn, libmlpack_julia_local_coordinate_coding, libmlpack_julia_logistic_regression, libmlpack_julia_lsh, libmlpack_julia_mean_shift, libmlpack_julia_nbc, libmlpack_julia_nca, libmlpack_julia_nmf, libmlpack_julia_pca, libmlpack_julia_perceptron, libmlpack_julia_preprocess_binarize, libmlpack_julia_preprocess_describe, libmlpack_julia_preprocess_one_hot_encoding, libmlpack_julia_preprocess_scale, libmlpack_julia_preprocess_split, libmlpack_julia_radical, libmlpack_julia_random_forest, libmlpack_julia_softmax_regression, libmlpack_julia_sparse_coding, libmlpack_julia_util

using boost_jll
using armadillo_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("mlpack")
JLLWrappers.@declare_library_product(libmlpack, "libmlpack.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_adaboost, "libmlpack_julia_adaboost.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_approx_kfn, "libmlpack_julia_approx_kfn.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_bayesian_linear_regression, "libmlpack_julia_bayesian_linear_regression.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_cf, "libmlpack_julia_cf.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_dbscan, "libmlpack_julia_dbscan.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_stump, "libmlpack_julia_decision_stump.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_tree, "libmlpack_julia_decision_tree.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_det, "libmlpack_julia_det.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_emst, "libmlpack_julia_emst.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_fastmks, "libmlpack_julia_fastmks.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_generate, "libmlpack_julia_gmm_generate.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_probability, "libmlpack_julia_gmm_probability.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_train, "libmlpack_julia_gmm_train.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_generate, "libmlpack_julia_hmm_generate.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_loglik, "libmlpack_julia_hmm_loglik.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_train, "libmlpack_julia_hmm_train.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_viterbi, "libmlpack_julia_hmm_viterbi.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_hoeffding_tree, "libmlpack_julia_hoeffding_tree.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_image_converter, "libmlpack_julia_image_converter.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_kde, "libmlpack_julia_kde.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_kernel_pca, "libmlpack_julia_kernel_pca.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_kfn, "libmlpack_julia_kfn.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_kmeans, "libmlpack_julia_kmeans.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_knn, "libmlpack_julia_knn.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_krann, "libmlpack_julia_krann.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_lars, "libmlpack_julia_lars.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_regression, "libmlpack_julia_linear_regression.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_svm, "libmlpack_julia_linear_svm.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_lmnn, "libmlpack_julia_lmnn.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_local_coordinate_coding, "libmlpack_julia_local_coordinate_coding.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_logistic_regression, "libmlpack_julia_logistic_regression.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_lsh, "libmlpack_julia_lsh.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_mean_shift, "libmlpack_julia_mean_shift.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_nbc, "libmlpack_julia_nbc.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_nca, "libmlpack_julia_nca.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_nmf, "libmlpack_julia_nmf.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_pca, "libmlpack_julia_pca.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_perceptron, "libmlpack_julia_perceptron.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_binarize, "libmlpack_julia_preprocess_binarize.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_describe, "libmlpack_julia_preprocess_describe.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_one_hot_encoding, "libmlpack_julia_preprocess_one_hot_encoding.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_scale, "libmlpack_julia_preprocess_scale.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_split, "libmlpack_julia_preprocess_split.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_radical, "libmlpack_julia_radical.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_random_forest, "libmlpack_julia_random_forest.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_softmax_regression, "libmlpack_julia_softmax_regression.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_sparse_coding, "libmlpack_julia_sparse_coding.dll")
JLLWrappers.@declare_library_product(libmlpack_julia_util, "libmlpack_julia_util.dll")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, armadillo_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libmlpack,
        "bin\\libmlpack.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_adaboost,
        "bin\\libmlpack_julia_adaboost.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_approx_kfn,
        "bin\\libmlpack_julia_approx_kfn.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_bayesian_linear_regression,
        "bin\\libmlpack_julia_bayesian_linear_regression.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_cf,
        "bin\\libmlpack_julia_cf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_dbscan,
        "bin\\libmlpack_julia_dbscan.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_decision_stump,
        "bin\\libmlpack_julia_decision_stump.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_decision_tree,
        "bin\\libmlpack_julia_decision_tree.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_det,
        "bin\\libmlpack_julia_det.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_emst,
        "bin\\libmlpack_julia_emst.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_fastmks,
        "bin\\libmlpack_julia_fastmks.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_generate,
        "bin\\libmlpack_julia_gmm_generate.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_probability,
        "bin\\libmlpack_julia_gmm_probability.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_train,
        "bin\\libmlpack_julia_gmm_train.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_generate,
        "bin\\libmlpack_julia_hmm_generate.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_loglik,
        "bin\\libmlpack_julia_hmm_loglik.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_train,
        "bin\\libmlpack_julia_hmm_train.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_viterbi,
        "bin\\libmlpack_julia_hmm_viterbi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hoeffding_tree,
        "bin\\libmlpack_julia_hoeffding_tree.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_image_converter,
        "bin\\libmlpack_julia_image_converter.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kde,
        "bin\\libmlpack_julia_kde.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kernel_pca,
        "bin\\libmlpack_julia_kernel_pca.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kfn,
        "bin\\libmlpack_julia_kfn.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kmeans,
        "bin\\libmlpack_julia_kmeans.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_knn,
        "bin\\libmlpack_julia_knn.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_krann,
        "bin\\libmlpack_julia_krann.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lars,
        "bin\\libmlpack_julia_lars.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_linear_regression,
        "bin\\libmlpack_julia_linear_regression.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_linear_svm,
        "bin\\libmlpack_julia_linear_svm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lmnn,
        "bin\\libmlpack_julia_lmnn.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_local_coordinate_coding,
        "bin\\libmlpack_julia_local_coordinate_coding.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_logistic_regression,
        "bin\\libmlpack_julia_logistic_regression.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lsh,
        "bin\\libmlpack_julia_lsh.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_mean_shift,
        "bin\\libmlpack_julia_mean_shift.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nbc,
        "bin\\libmlpack_julia_nbc.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nca,
        "bin\\libmlpack_julia_nca.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nmf,
        "bin\\libmlpack_julia_nmf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_pca,
        "bin\\libmlpack_julia_pca.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_perceptron,
        "bin\\libmlpack_julia_perceptron.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_binarize,
        "bin\\libmlpack_julia_preprocess_binarize.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_describe,
        "bin\\libmlpack_julia_preprocess_describe.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_one_hot_encoding,
        "bin\\libmlpack_julia_preprocess_one_hot_encoding.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_scale,
        "bin\\libmlpack_julia_preprocess_scale.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_split,
        "bin\\libmlpack_julia_preprocess_split.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_radical,
        "bin\\libmlpack_julia_radical.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_random_forest,
        "bin\\libmlpack_julia_random_forest.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_softmax_regression,
        "bin\\libmlpack_julia_softmax_regression.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_sparse_coding,
        "bin\\libmlpack_julia_sparse_coding.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_util,
        "bin\\libmlpack_julia_util.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    set(str,var::Symbol,val)

Substitutes the value a single field `var` with `val`. Returns the modified `str`.
"""
function set(str, var::Symbol, val)
    !any(i -> i == var, fieldnames(typeof(str))) && @warn "\n No field in struct $(typeof(str)) is named $(var)!\n Fields are $(fieldnames(typeof(str)))"
    dict = Dict{Symbol,Any}()
    for (i, name) in enumerate(fieldnames(typeof(str)))
        if typeof(getfield(str, name)) <: Parameter && name == var
            dict[name] = Parameter(val, getfield(getfield(str, name), :description))
        elseif typeof(getfield(str, name)) <: Parameter
            dict[name] = Parameter(getfield(getfield(str, name), :value), getfield(getfield(str, name), :description))
        elseif typeof(getfield(str, name)) <: BoundedParameter && name == var
            dict[name] = BoundedParameter(val, getfield(getfield(str, name), :lb), getfield(getfield(str, name), :ub), getfield(getfield(str, name), :description))
        elseif typeof(getfield(str, name)) <: BoundedParameter
            dict[name] = BoundedParameter(getfield(getfield(str, name), :value), getfield(getfield(str, name), :lb), getfield(getfield(str, name), :ub), getfield(getfield(str, name), :description))
        elseif name == var
            dict[name] = val
        else
            dict[name] = getfield(str, name)
        end
    end
    x = NamedTuple{Tuple(keys(dict))}(values(dict))
    str = typeof(str)(; x...)
end
"""
    (←)(str, tpl::Tuple)

Substitutes the value of a single or multiple fields of struct `str`. The tuple has the form (field::Symbol, new_value). Returns the modified `str`.
"""
function (←)(str, tpl::Tuple)
    if typeof(values(tpl)[1]) ≠ Symbol
        for value in values(tpl)
            str = set(str, value[1], value[2])
        end
    else
        str = set(str, tpl[1], tpl[2])
    end
    return str
end


text:
@testset "GeodesicComplex" begin

    @testset "Construction & data structure properties" begin
        for (pts, metric, r) in datasets
            lmks, cover = GeodesicComplexes.getcover(nothing, pts, r, KDTree(pts, metric))
            @test all(1 .≤ lmks .≤ length(pts))
            @test reduce(union, cover) == Set(1:length(pts))
            @test issorted(lmks)
            @test all(eachindex(lmks)) do i
                all(evaluate(metric, pts[lmks[i]], pts[j]) ≤ r for j in cover[i])
            end
            gc = GeodesicComplex(pts, r, metric = metric)

            # TODO: indexing
            lmrks = landmarks(gc)
            othrs = nonlandmarks(gc)
            @test allunique(lmrks)
            @test allunique(othrs)
            @test length(intersect(lmrks, othrs)) == 0
            @test length(lmrks) + length(othrs) == length(points(gc))
            @test Set(union(lmrks, othrs)) == Set(points(gc))
            nlandmarks(gc) == length(lmrks)

            @test all(weight(e) ≈ evaluate(metric, lmrks[src(e)], lmrks[dst(e)])
                      for e in edges(gc))
            @test all(r ≤ weight(e) ≤ 2r
                      for e in edges(gc))
        end
    end

    @testset "LightGraphs interface" begin
        for (pts, metric, r) in datasets
            gc = GeodesicComplex(pts, r, metric = metric)
            # Run some functions to see that the interface is defined properly
            @test betweenness_centrality(gc) isa Array
            @test dijkstra_shortest_paths(gc, 1) isa LightGraphs.DijkstraState
        end
    end

    @testset "Plotting" begin
        # Count number of series.
        d = Dict{Symbol, Any}()
        for (pts, metric, r) in datasets
            gc = GeodesicComplex(pts, r, metric = metric)
            # ¯\_(ツ)_/¯
            @test length(RecipesBase.apply_recipe(d, gc)) == 2
        end
    end
end


text:
# The following method was removed due to the fact that it is a type piracy from Distributions.jl
# Too bad;

#Distributions.convolve(d1::Distributions.Gamma,d2::Distributions.Gamma) = UnivariateGammaConvolution([d1.α,d2.α],[d1.θ,d2.θ])
Distributions.convolve(d1::Distributions.Gamma,d2::UnivariateGammaConvolution) = UnivariateGammaConvolution([d1.α,d2.α...], [d1.θ,d2.θ...])
Distributions.convolve(d1::UnivariateGammaConvolution,d2::Distributions.Gamma) = Distribution.convolve(d2,d1)
# Products of UnivariateGammaConvolution

"""
    Distributions.product_distribution(d1,d2)

Implements the special case of the product of two gamma or univariate gamma convolutions distributions, and
output as the result a multivariate gamma convolution with independent margins.
"""
function Distributions.product_distribution(dists::AbstractVector{<:UnivariateGammaConvolution})

    # number of gammas :
    T = Base.promote_eltype(dists)
    d = length(dists)
    ns = [length(d.α) for d in dists]
    n = sum(ns)
    α = zeros(T,n)
    θ = zeros(T,(n,d))
    current_n = 1
    for i in 1:d
        α[current_n:(current_n-1+ns[i])] = dists[i].α
        θ[current_n:(current_n-1+ns[i]),i] = dists[i].θ
        current_n +=ni
    end
    return MultivariateGammaConvolution(α,θ)
end


text:
using Pkg
Pkg.activate(@__DIR__)

using Documenter, DimensionalData, CoordinateTransformations, Dates

CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing

docsetup = quote
    using DimensionalData, Random, Dates
    Random.seed!(1234)
end
DocMeta.setdocmeta!(DimensionalData, :DocTestSetup, docsetup; recursive=true)

makedocs(
    modules = [DimensionalData],
    sitename = "DimensionalData.jl",
    format = Documenter.HTML(
        prettyurls = CI,
    ),
    pages = [
        "Introduction" => "index.md",
        "Crash course" => "course.md",
        "API" => "api.md",
        "For Developers" => "developer.md"
        ],
    strict=true,
)

if CI
    deploydocs(
        repo = "github.com/rafaqz/DimensionalData.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
"""
This Direct method assumes that there are a fixed total number
of transitions and that each clock is marked with :Index,
numbered from 1, as an extra argument to the AddTransition!.
"""
struct FixedDirect
    tree::PrefixSearchTree{Float64}
    N::Int
    clock_index::Dict{Int, Any}
    init::Bool
    FixedDirect(N::Int) = new(
        PrefixSearchTree(Float64, N), N,
        Dict{Int,Any}(), true
        )
end

function fd_indexof(kind::Array{Any,1})
    for (symbol, value) in kind
        if symbol==:index
            return value
        end
    end
    error("Need an index=<int> for AddTransition!")
end

function Next(propagator::FixedDirect, process, rng)
    if propagator.init
        hazards=Array{Tuple{Int, Float64}, 1}()
        Hazards(process, rng) do clock, now, enabled, rng2
            lambda=Parameters(Distribution(clock.intensity))[1]
            index=fd_indexof(clock.kind)
            propagator.clock_index[index]=clock
            push!(hazards, (index, lambda))
        end
        push!(propagator.tree, hazards)
        propagator.init=false
    end
    total=sum(propagator.tree)
    if total>eps(Float64)
        (index, value)=choose(propagator.tree, rand(rng)*total)
        clock=propagator.clock_index[index]
        return (Time(process)-log(rand(rng))/total, clock)
    else
        return (Inf, nothing)
    end
end

function Observer(propagator::FixedDirect)
    function fdobserve(clock, time, updated, rng)
        if updated != :Disabled && updated != :Fired
            index = fd_indexof(clock.kind)
            propagator.clock_index[index] = clock
            lambda = Parameters(Distribution(clock.intensity))[1]
            push!(propagator.tree, index, lambda)
        else
            index=fd_indexof(clock.kind)
            push!(propagator.tree, index, 0.0)
        end
    end
end


text:
# lnprob evaluation for Gaussian

using visibilities
using image
using gauss_model
using gridding
using constants

const dvarr = DataVis("data/V4046Sgr_fake.hdf5") #Read in the whole fake dataset
nlam = length(dvarr)

# Dimensions of the problem
const nx = 256
const ny = 256

const ra = fftspace(10., nx) # [arcsec]
const dec = fftspace(10., ny) # [arcsec]

# convert ra and dec in [arcsec] to radians, and then take the sin to convert to ll, mm
ll = sin(ra * arcsec) # direction cosines
mm = sin(dec * arcsec)

function f(p::Vector{Float64})
    # p is naturally in arcsec

    # Compute the array of images
    data = Array(Float64, (ny, nx, nlam))
    lams = Array(Float64, nlam)
    sigma_x, sigma_y = p[3:4]
    for k=1:nlam
        data[:, :, k] = imageGauss(ll, mm, Float64[0., 0., sigma_x, sigma_y], k)
        lams[k] = dvarr[k].lam # [microns]
    end

    skim = SkyImage(data, ra, dec, lams)

    mu_RA, mu_DEC = p[1:2]

    # Apply the gridding correction function before doing the FFT
    corrfun!(skim, 1.0, mu_RA, mu_DEC) #alpha = 1.0

    # Now calculate the sum of lnprobs in a loop
    lnp = Array(Float64, nlam)
    for i=1:nlam
        # FFT the appropriate image chanel
        vis_fft = transform(skim, i)

        # Interpolate the `vis_fft` at the same locations as the DataSet
        mvis = ModelVis(dvarr[i], vis_fft)

        # Do the phase shift for the offset
        phase_shift!(mvis, mu_RA, mu_DEC)

        # Calculate chi^2 between these two
        lnp[i] = lnprob(dvarr[i], mvis)
    end
    return sum(lnp)
end

# Necessary wrapper for NLopt requires gradient as an argument (even if it's not used)
function fgrad(p::Vector, grad::Vector)
    val = f(p)
    println(p, " : ", val)
    return val
end

function fp(p::Vector)
    val = f(p)
    println(p, " : ", val)
    return val
end

# # Now try optimizing the function using NLopt
using NLopt

starting_param = [0.0, 0.0, 1.0, 1.0]

nparam = length(starting_param)
opt = Opt(:LN_COBYLA, nparam)

max_objective!(opt, fgrad)
xtol_rel!(opt,1e-4)

(optf,optx,ret) = optimize(opt, starting_param)
println(optf, " ", optx, " ", ret)

quit()

using LittleMC

using Distributions
using PDMats

mc = MC(fp, 500, [1.2, 1.0, 1.0, 1.0], PDiagMat([0.03^2, 0.03^2, 0.01^2, 0.01^2]))

start(mc)

println(mean(mc.samples, 2))
println(std(mc.samples, 2))

runstats(mc)

write(mc, "mc.hdf5")

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey