text:
# Autogenerated wrapper script for nghttp3_jll for x86_64-apple-darwin
export libnghttp3

JLLWrappers.@generate_wrapper_header("nghttp3")
JLLWrappers.@declare_library_product(libnghttp3, "@rpath/libnghttp3.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnghttp3,
        "lib/libnghttp3.0.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
`FWF.parsefield{T}(source::FWF.Source, ::Type{T}, row::Int, col::Int)` => `Nullable{T}`
`FWF.parsecol{T}(source::FWF.Source, ::Type{T}, col::Int)` => `Nullable{T}`

`source` is the source to read from
`T` is the type of the column / field
`row` is the row that is being read, if applicable
`col` is the column that is being read.

whitespace is ignored for numerial types, string trimming is configurable by options.
If `checkfornulls` is set in `opt` fields will be compared to the null list in `opt`, any values found will result in a missing value
Parsing happens for Integrer, Float64, Date
The parameter list without a row results in a whole column being streamed from the file.

"""
function parsefield end

missingon(source::FWF.Source) = (source.options.usemissings)
checkmissing(key::String, d::Set{String}) = key in d

function get_format(source::FWF.Source, col::Int)
    !haskey(source.options.dateformats, col) && return nothing
    return source.options.dateformats[col]
end

# This main dispatch will get called for field based parsing and dispatch to other methods
# We are going to pre-load whole lines and cache the results in the Source.
# Assume that every col==1 means load the next line
function parsefield(source::FWF.Source, ::Type{T}, row::Int, col::Int) where {T}
    # Assume that every col==1 means load the next line
    if !((col == source.lastcol+1) || (col == 1 && Data.size(source.schema)[2] == source.lastcol))
        throw(FWF.ParsingException("Out of order access col=$col lastcol=$(source.lastcol)"))
    end
    source.lastcol = col
    if col == 1
        readsplitline!(source.currentline, source)
    end

    tmp = source.currentline[col]
    ismissing(tmp) && return missing
    # We know ref is a string
    ref::String = tmp
    if missingon(source) && checkmissing(ref, source.options.missingvals)
        return missing
    end

    return parsefield(T, source.options.usemissings, ref, get_format(source, col))
end

# Batch of simple parsers to convert strings
usemissing_or_val(b, v) = b ? missing : v
if VERSION >= v"0.7.0-DEV"
    null_to_missing(x, b, v) = x == nothing ? usemissing_or_val(b, v) : x
else
    null_to_missing(x, b, v) = isnull(x) ? usemissing_or_val(b, v)  : unsafe_get(x)
end

parsefield(::Type{Int}, usemissing::Bool, string::String, format) =
    null_to_missing(tryparse(Int, string), usemissing, 0)
parsefield(::Type{Float64}, usemissing::Bool, string::String, format) =
    null_to_missing(tryparse(Float64, string), usemissing, 0.)
parsefield(::Type{String}, usemissing::Bool, string::String, format) = string
parsefield(::Type{Date}, usemissing::Bool, string::String, format) =
    null_to_missing(tryparse(Date, string, format), usemissing, Date())
parsefield(::Type{Union{Missing, T}}, usemissing::Bool, string::String, format) where {T} =
    (parsefield(T, usemissing, string, format))
parsefield(::Type{Missing}, usemissing::Bool, string::String, format) = (missing)

# Generic fallback
function parsefield(T, usemissing::Bool, string::String, format)
    return null_to_missing(tryparse(T, string), usemissing, nothing)
end

ext: jl


text:
## --- Utility functions for plotting

    function plot_rankorder_errorbar(data,uncert;
            scale=1,
            seriestype=:scatter,
            ylabel="",
            label="",
            xticks=[],
            xlabel="",
            color=:auto,
            markersize=2
        )
        h = plot()
        i0 = 0
        plot_rankorder_errorbar!(h,data,uncert,i0;
                scale=scale,
                seriestype=seriestype,
                label=label,
                ylabel=ylabel,
                xlabel=xlabel,
                xticks=xticks,
                color=color,
                markersize=markersize
        )
        return h
    end

    function plot_rankorder_errorbar!(h::Plots.Plot,data,uncert,i0;
            scale=1,
            seriestype=:scatter,
            label="",
            ylabel="",
            xlabel="",
            xticks=[],
            color=:auto,
            markersize=2
        )
        sI = sortperm(data)
        plot!(h, i0 .+ (1:length(sI))*scale,data[sI],
            yerror=uncert[sI],
            seriestype=seriestype,
            markerstrokecolor=color,
            color=color,
            label=label,
            ylabel=ylabel,
            xlabel=xlabel,
            xticks=xticks,
            markersize=markersize
        )
    end

## --- Remove outliers

    function screen_outliers(smpl::ChronAgeData; maxgap=100)
        # Variables from struct
        Name = collect(smpl.Name)::Array{String,1}
        Path = smpl.Path::String
        Age_Unit = smpl.Age_Unit::String

        # Create and populate directory of screened output
        screenedpath = joinpath(Path, "screened/")
        mkpath(screenedpath)
        for i=1:length(Name)
            # With screening
            # Maximum offset before cutoff
            data = readdlm(joinpath(Path, Name[i]*".csv"), ',', Float64)::Array{Float64,2}
            data = sortslices(data,dims=1) # Sort ages in ascending order
            nAnalyses = size(data,1)
            maxdt_sigma = maxgap*norm_width(nAnalyses)/nAnalyses

            # Rank-order plot of all individual ages for comparison
            hdl = plot(1:nAnalyses,data[:,1],yerror=data[:,2]*2/smpl.inputSigmaLevel, seriestype=:scatter, color=:red, markerstrokecolor=:red,label="rejected",legend=:topleft,framestyle=:box,fg_color_legend=:white)

            # Filter data to exclude outliers
            sigma_mutual = nanmean(data[:,2]) / smpl.inputSigmaLevel * sqrt(2)
            for j=nAnalyses:-1:2
                dt_sigma = abs(data[j,1]-data[j-1,1]) / sigma_mutual # Time gap divided by relative sigma

                # If we exceed the maximum allowed dt/sigma, delete any points
                # below (older than) the gap
                if dt_sigma>maxdt_sigma && j>2
                    data=data[1:j-1,:]
                end
            end
            plot!(hdl, 1:size(data,1),data[:,1],yerror=data[:,2]*2/smpl.inputSigmaLevel, seriestype=:scatter, color=:blue,markerstrokecolor=:blue,label="included",xlabel="N",ylabel="Age ($(Age_Unit))")
            savefig(hdl, joinpath(screenedpath, Name[i]*"_screening.pdf"))
            writedlm(joinpath(screenedpath, Name[i]*".csv"), data, ',')
        end
        smpl.Path = screenedpath
        return smpl
    end

## --- Bootstrap prior distribution shape

    """
    ```julia
    BootstrapCrystDistributionKDE(smpl::ChronAgeData)
    ```
    Bootstrap an estimate of the pre-eruptive (or pre-deposition) mineral
    crystallization distribution shape from a Chron.ChronAgeData object containing
    data for several samples, using a kernel density estimate of stacked sample data.
    """
    function BootstrapCrystDistributionKDE(smpl::ChronAgeData; cutoff::Number=-0.05)
        # Extact variables froms struct
        Name = collect(smpl.Name)::Array{String,1}
        Path = smpl.Path::String
        DistType = smpl.Age_DistType::Array{Float64,1}

        # Load all data points and scale from 0 to 1
        allscaled = Array{Float64}([])
        for i=1:length(Name)
            if DistType[i]==0
                data = readdlm(joinpath(Path, Name[i]*".csv"), ',', Float64)::Array{Float64,2}

                # Maximum extent of expected analytical tail (beyond eruption/deposition)
                maxTailLength = mean(data[:,2]) ./ smpl.inputSigmaLevel .* norm_quantile(1 - 1/(1+size(data,1)))
                included = (data[:,1] .- minimum(data[:,1])) .>= maxTailLength
                included .|= data[:,1] .> nanmedian(data[:,1]) # Don't exclude more than half (could only happen in underdispersed datasets)
                included .&= .~isnan.(data[:,1]) # Exclude NaNs

                # Include and scale only those data not within the expected analytical tail
                if sum(included)>0
                    scaled = data[included,1] .- minimum(data[included,1])
                    if maximum(scaled) > 0
                        scaled = scaled ./ maximum(scaled)
                    end
                    append!(allscaled, scaled)
                end
            end
        end

        # Calculate kernel density estimate, truncated at 0
        kd = kde(allscaled,npoints=2^7)
        t = kd.x .> cutoff # Ensure sharp cutoff at eruption / deposition
        return kd.density[t]
    end

    """
    ```julia
    BootstrapCrystDistributionKDE(data::AbstractArray, [sigma::AbstractArray])
    ```
    Bootstrap an estimate of thq;e pre-eruptive (or pre-deposition) mineral
    crystallization distribution shape from a 2-d array of sample ages (one row per
    sample, one column per datum, padded with NaNs as needed) and an equivalent-size
    array of one-sigma uncertainties,
    using a kernel density estimate of stacked sample data.
    """
    function BootstrapCrystDistributionKDE(data::AbstractArray{<:Number}; cutoff::Number=-0.05)
        # Load all data points and scale from 0 to 1
        allscaled = Array{float(eltype(data)),1}([])
        for i=1:size(data,2)
            scaled = data[:,i] .- minimum(data[:,i])
            if maximum(scaled) > 0
                scaled = scaled ./ maximum(scaled)
            end
            append!(allscaled, scaled)
        end

        # Calculate kernel density estimate, truncated at 0
        kd = kde(allscaled,npoints=2^7)
        t = kd.x .> cutoff
        return kd.density[t]
    end

    function BootstrapCrystDistributionKDE(data::AbstractArray{<:Number}, sigma::AbstractArray{<:Number}; cutoff::Number=-0.05)
        # Array to hold stacked, scaled data
        allscaled = Array{float(eltype(data)),1}([])

        # For each row of data
        for i=1:size(data,2)

            # Maximum extent of expected analytical tail (beyond eruption/deposition/cutoff)
            maxTailLength = nanmean(sigma[:,i]) .* norm_quantile(1 - 1/(1+count(.!isnan.(data[:,i]))))
            included = (data[:,i] .- nanminimum(data[:,i])) .>= maxTailLength
            included .|= data[:,i] .> nanmedian(data[:,i]) # Don't exclude more than half (could only happen in underdispersed datasets)
            included .&= .~isnan.(data[:,i]) # Exclude NaNs

            # Include and scale only those data not within the expected analytical tail
            if sum(included)>0
                scaled = data[included,i] .- minimum(data[included,i])
                if maximum(scaled) > 0
                    scaled = scaled ./ maximum(scaled)
                end
                append!(allscaled, scaled)
            end
        end

        # Calculate kernel density estimate, truncated at 0
        kd = kde(allscaled,npoints=2^7)
        t = kd.x .> cutoff
        return kd.density[t]
    end

## --- Fit and plot results from stationary distribution of depostion/eruption age distribution model

    """
    ```julia
    smpl = tMinDistMetropolis(smpl::ChronAgeData,nsteps::Int,burnin::Int,dist::Array{Float64})
    ```
    Calculate the minimum limiting (eruption/deposition) age of each sample in `smpl`
    using the `metropolis_min` function, assuming mineral ages for each sample are
    drawn from the source distribution `dist`. Fits a `bilinearexponential` function
    to the resulting stationary distribution for each sample.
    """
    function tMinDistMetropolis(smpl::ChronAgeData,nsteps::Int,burnin::Int,dist::Array{Float64})
        # Extract variables from struct
        Name = collect(smpl.Name)::Array{String,1}
        Path = smpl.Path::String
        Age_Unit = smpl.Age_Unit::String
        DistType = smpl.Age_DistType::Array{Float64,1}

        # Estimate the eruption/deposition distribution for each sample
        print("Estimating eruption/deposition age distributions...\n")
        for i=1:length(Name)
            if DistType[i] == 0 # A distribution to fit properly
                # Load data for each sample
                data = readdlm(joinpath(Path, Name[i]*".csv"), ',', Float64)::Array{Float64,2}
                print(i, ": ", Name[i], "\n") # Display progress

                # Run MCMC to estimate saturation and eruption/deposition age distributions
                tminDist = metropolis_min(nsteps,dist,data[:,1],data[:,2]/smpl.inputSigmaLevel; burnin=burnin) # Since we don't end up using any of the other distributions
                # (tminDist, tmaxDist, llDist, acceptanceDist) = metropolis_minmax(nsteps,dist,data[:,1],data[:,2]/smpl.inputSigmaLevel, burnin=burnin)

                # Fill in the strat sample object with our new results
                smpl.Age[i] = mean(tminDist)
                smpl.Age_sigma[i] = std(tminDist)
                smpl.Age_025CI[i] = percentile(tminDist,2.5)
                smpl.Age_975CI[i] = percentile(tminDist,97.5)
                smpl.Age_Distribution[i] = tminDist

                # Fit custom many-parametric distribution function to histogram
                edges = range(minimum(tminDist),maximum(tminDist),length=101) # Vector of bin edges
                hobj = fit(Histogram,tminDist,edges,closed=:left) # Fit histogram object

                t = hobj.weights.>0 # Only look at bins with one or more results
                N = hobj.weights[t] ./ length(tminDist) .* length(t) # Normalized number of MCMC steps per bin
                bincenters = cntr(hobj.edges[1])[t] # Vector of bin centers

                # Initial guess for parameters
                p = ones(5)
                p[1] = log(maximum(N))
                p[2] = mean(tminDist)
                p[3] = std(tminDist)

                # Fit nonlinear model
                fobj = curve_fit(bilinear_exponential,bincenters,N,p)
                smpl.Params[:,i] = fobj.param

                # Rank-order plot of analyses and eruption/deposition age range
                nAnalyses = length(data[:,1])
                h1 = plot_rankorder_errorbar(data[:,1],2*data[:,2]/smpl.inputSigmaLevel,ylabel="Age ($(Age_Unit))",label="Data (observed ages)")
                m = ones(nAnalyses).*smpl.Age[i]
                l = ones(nAnalyses).*smpl.Age_025CI[i]
                u = ones(nAnalyses).*smpl.Age_975CI[i]
                plot!(h1,1:nAnalyses,l,fillto=u,fillalpha=0.6,linealpha=0, label="Model ($(round(m[1],digits=3)) +$(round(u[1]-m[1],digits=3))/-$(round(m[1]-l[1],digits=3)) $(Age_Unit))")
                plot!(h1,1:nAnalyses,m,linecolor=:black,linestyle=:dot,label="",legend=:topleft,fg_color_legend=:white,framestyle=:box)
                savefig(h1,joinpath(Path, Name[i]*"_rankorder.pdf"))
                savefig(h1,joinpath(Path, Name[i]*"_rankorder.svg"))

                # Plot model fit to histogram
                h2 = plot(bincenters,N,label="Histogram",fg_color_legend=:white,framestyle=:box)
                plot!(h2,bincenters, bilinear_exponential(bincenters,smpl.Params[:,i]), label="Curve fit")
                plot!(h2,legend=:topleft, xlabel="Age ($(Age_Unit))", ylabel="Probability density")
                savefig(h2,joinpath(Path, Name[i]*"_distribution.pdf"))
                savefig(h2,joinpath(Path, Name[i]*"_distribution.pdf"))

            elseif DistType[i] == 1 # A single Gaussian
                # Load data for each sample
                data = readdlm(joinpath(Path, Name[i]*".csv"), ',', Float64)::Array{Float64,2}
                print(i, ": ", Name[i], "\n") # Display progress
                μ = data[1,1]
                σ = data[1,2]

                # Fill in the strat sample object with our new results
                smpl.Age[i] = μ
                smpl.Age_sigma[i] = σ
                smpl.Age_025CI[i] = μ - 1.95996398454*σ
                smpl.Age_975CI[i] = μ + 1.95996398454*σ

                # Initial guess for parameters
                p = ones(5)
                p[1] = log(normpdf(0,σ,0))
                p[2] = μ
                p[3] = σ

                # Fit nonlinear model
                x = μ .+ (-10σ:σ/10:10σ)
                fobj = curve_fit(bilinear_exponential, x, normpdf(μ, σ, x), p)
                fobj.param[5] = 1 # Must be symmetrical
                smpl.Params[:,i] = fobj.param
            end
        end

        # Save results as csv
        results = vcat(["Sample" "Age" "2.5% CI" "97.5% CI" "sigma"], hcat(Name,smpl.Age,smpl.Age_025CI,smpl.Age_975CI,smpl.Age_sigma))::Array{Any,2}
        writedlm(joinpath(Path, "distresults.csv"), results, ',')

        return smpl
    end

## ---


text:
using GeoStatsDevTools
using Test

# EMPTY


text:
module effective_interaction_rotated

# Optimized functions for calculation of Omega and Gamma

function Omega(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    cosξdivξ = cos(ξ)/ξ
    return 3.0/4*((1.0-3.0*cosθpow2) * (sin(ξ)+cosξdivξ)/ξ^2 - (1.0-cosθpow2)*cosξdivξ)
end

function Gamma(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    if ξ<1e-3
        return 3/2*(2/3 + 1/15*ξ^2 *(-2 + cosθpow2) + ξ^4*(1/140-cosθpow2/210))
    else
        sinξ = sin(ξ)/ξ
        return 3/2*((1-cosθpow2)*sinξ + (1-3*cosθpow2) * (cos(ξ)-sinξ)/ξ^2)
    end
end

function OmegaGamma(a, cosθpow2)
    #cosθpow2 = cos(θ)^2
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((1.0-3.0*cosθpow2) * (sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - (1.0-cosθpow2)*cosξ/(2*pi*a))
    gamma = 1.5*((1-cosθpow2)*sinξ/(2*pi*a) + (1-3*cosθpow2) * (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function OmegaGamma_orthogonal(a)
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - cosξ/(2*pi*a))
    gamma = 1.5*(sinξ/(2*pi*a) + (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function Omega_orthogonal(a)
    ξ = 2*pi*a
    cosξdξ = cos(ξ)/ξ
    return 3.0/4.0*(-cosξdξ + (sin(ξ)+cosξdξ)/ξ^2)
end

function Gamma_orthogonal(a)
    ξ = 2*a
    sincξ = sinc(ξ)
    return 3.0/2.0*(sincξ + (cospi(ξ)-sincξ)/(pi*ξ)^2)
end

function map_effectiveinteraction(a::Vector{Float64}, f::Function)
    omega_list = Float64[]
    gamma_list = Float64[]
    for ai=a
        omega, gamma = f(ai, e_dipole)
        push!(omega_list, omega)
        push!(gamma_list, gamma)
    end
    return omega_list, gamma_list
end

# Finite symmetric systems

# function triangle_orthogonal(a::Float64)
#     omega_eff = 2*Omega_orthogonal(a)
#     gamma_eff = 2*Gamma_orthogonal(a)
#     return omega_eff, gamma_eff
# end

"""
    effective_interaction_rotated.square_orthogonal(a, Nδ)

Effective Omega and Gamma for a square.

The polarization axis is orthogonal to the square plane.

# Arguments
* `a`: Edge length.
* `Nδ`: Phase shift (Number of atoms in 2π).
"""
function square_orthogonal(a::Float64, Nδ::Int)
    @assert 0 <= Nδ < 4
    δ = 2.0*pi*Nδ/4
    Ωcos = 2*Omega_orthogonal(a)*cos(δ) + Omega_orthogonal(sqrt(2.)*a)*cos(2*δ)
    Γcos = 2*Gamma_orthogonal(a)*cos(δ) + Gamma_orthogonal(sqrt(2.)*a)*cos(2*δ)
    return Ωcos, Γcos
end

# function rectangle_orthogonal(a::Float64, b::Float64)
#     d = sqrt(a^2 + b^2)
#     omega_eff, gamma_eff = OmegaGamma_orthogonal(a)
#     om, ga = OmegaGamma_orthogonal(b)
#     omega_eff += om; gamma_eff+=ga
#     om, ga = OmegaGamma_orthogonal(d)
#     omega_eff += om; gamma_eff+=ga
#     return omega_eff, gamma_eff
# end

"""
    effective_interaction_rotated.cube_orthogonal(a, dϕ)

Effective Omega and Gamma for a cube.

The polarization axis is orthogonal to the xy faces.

# Arguments
* `a`: edge length.
* `dϕ`: Phase shift.
"""
function cube_orthogonal(a::Float64, dϕ)
    omega_eff, gamma_eff = square_orthogonal(a, 0)
    sqrt2 = sqrt(2.)
    sqrt3 = sqrt(3.)
    Θdiag = atan(sqrt2, 1.)
    omega_eff += (Omega(a, 0.) + 2*Omega(sqrt(2.)*a, pi/4.) + Omega(sqrt3*a, Θdiag))*cos(dϕ)
    gamma_eff += (Gamma(a, 0.) + 2*Gamma(sqrt(2.)*a, pi/4.) + Gamma(sqrt3*a, Θdiag))*cos(dϕ)
    return omega_eff, gamma_eff
end

# function box_orthogonal(a::Float64, b::Float64, c::Float64)
#     omega_eff, gamma_eff = rectangle_orthogonal(a, b)

#     a2 = a^2
#     b2 = b^2
#     c2 = c^2

#     om, ga = OmegaGamma(c, 1.)
#     omega_eff += om; gamma_eff+=ga

#     dxz_pow2 = a2+c2
#     om, ga = OmegaGamma(sqrt(dxz_pow2), c2/dxz_pow2)
#     omega_eff += om; gamma_eff+=ga

#     dyz_pow2 = b2+c2
#     om, ga = OmegaGamma(sqrt(dyz_pow2), c2/dyz_pow2)
#     omega_eff += om; gamma_eff+=ga

#     dxyz_pow2 = a2+b2+c2
#     om, ga = OmegaGamma(sqrt(dxyz_pow2), c2/dxyz_pow2)
#     omega_eff += om; gamma_eff+=ga
#     return omega_eff, gamma_eff
# end

# Infinite 1D symmetric systems

# function chain(a::Float64, Θ, N::Int)
#     omega_eff::Float64 = 0.
#     gamma_eff::Float64 = 0.
#     d = 0.
#     for j=1:N
#         d += a
#         omega_eff += Omega(d, Θ)
#         gamma_eff += Gamma(d, Θ)
#     end
#     return 2*omega_eff, 2*gamma_eff
# end

"""
    effective_interaction_rotated.chain_orthogonal(a, N, dϕ)

Effective Omega and Gamma for an infinite chain.

The polarization axis is orthogonal to the chain and the calculation is
done by adding N spins left and N spins right of a central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
* `dϕ`: Phase shift between neighboring spins.
"""
function chain_orthogonal(a::Float64, N::Int, dϕ::Float64)
    omega_eff::Float64 = 0.
    gamma_eff::Float64 = 0.
    d::Float64 = 0.
    ϕ::Float64 = 0.
    for j=1:N
        d += a
        ϕ += dϕ
        omega_eff += Omega_orthogonal(d)*cos(ϕ)
        gamma_eff += Gamma_orthogonal(d)*cos(ϕ)
    end
    return 2*omega_eff, 2*gamma_eff
end

# Infinite 2D symmetric systems

# function squarelattice_orthogonal(a::Float64, N::Int)
#     omega_eff::Float64 = 0.
#     gamma_eff::Float64 = 0.
#     for j=1:N
#         d = a*j
#         omega_eff += 4*Omega_orthogonal(d)
#         gamma_eff += 4*Gamma_orthogonal(d)
#         dx2 = j^2
#         for i=1:j-1
#             dy2 = i^2
#             d = a*sqrt(dx2+dy2)
#             omega_eff += 8*Omega_orthogonal(d)
#             gamma_eff += 8*Gamma_orthogonal(d)
#         end
#         d = a*sqrt(2)*j
#         omega_eff += 4*Omega_orthogonal(d)
#         gamma_eff += 4*Gamma_orthogonal(d)
#     end
#     return omega_eff, gamma_eff
# end

# function hexagonallattice_orthogonal(a::Float64, N::Int)
#     omega_eff::Float64 = 0.
#     gamma_eff::Float64 = 0.
#     for n=1:2:N
#         nv1x = -0.5*n
#         nv1y_square = 3./4*n^2
#         d = a*sqrt(nv1y_square + nv1x^2)
#         omega_eff += 6*Omega_orthogonal(d)
#         gamma_eff += 6*Gamma_orthogonal(d)
#         for i=1:div(n-1,2)
#             d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
#             omega_eff += 12*Omega_orthogonal(d)
#             gamma_eff += 12*Gamma_orthogonal(d)
#         end
#     end
#     for n=2:2:N
#         nv1x = -0.5*n
#         nv1y_square = 3./4*n^2
#         d = a*sqrt(nv1y_square + nv1x^2)
#         omega_eff += 6*Omega_orthogonal(d)
#         gamma_eff += 6*Gamma_orthogonal(d)
#         for i=1:div(n-2,2)
#             d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
#             omega_eff += 12*Omega_orthogonal(d)
#             gamma_eff += 12*Gamma_orthogonal(d)
#         end
#         i = div(n-2,2)+1
#         d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
#         omega_eff += 6*Omega_orthogonal(d)
#         gamma_eff += 6*Gamma_orthogonal(d)
#     end
#     return omega_eff, gamma_eff
# end

# Infinite 3D symmetric systems

# function cubiclattice_orthogonal(a::Float64, N::Int)
#     omega_eff, gamma_eff = squarelattice_orthogonal(a, N)
#     for nx=1:N
#         for ny=0:nx
#             multiplicity = ((ny==0 || ny==nx) ? 8 : 16)
#             nrpow2 = nx^2 + ny^2
#             nr = sqrt(nrpow2)
#             for nz=1:N
#                 nzpow2 = nz^2
#                 d = a*sqrt(nrpow2+nzpow2)
#                 cosθpow2 = nzpow2/(nrpow2+nzpow2)
#                 om, ga = OmegaGamma(d, cosθpow2)
#                 omega_eff += multiplicity*om
#                 gamma_eff += multiplicity*ga
#             end
#         end
#     end
#     d::Float64 = 0.
#     for nz=1:N
#         d += a
#         om, ga = OmegaGamma(d, 1.)
#         omega_eff += 2*om
#         gamma_eff += 2*ga
#     end
#     return omega_eff, gamma_eff
# end

# function tetragonallattice_orthogonal(a::Float64, b::Float64, N::Int)
#     omega_eff, gamma_eff = squarelattice_orthogonal(a, N)
#     for nx=1:N
#         for ny=0:nx
#             multiplicity = ((ny==0 || ny==nx) ? 8 : 16)
#             nrpow2 = nx^2 + ny^2
#             nr = sqrt(nrpow2)
#             for nz=1:N
#                 nzpow2 = nz^2
#                 d = sqrt(a^2*nrpow2 + b^2*nzpow2)
#                 θ = atan2(a*nr, b*nz)
#                 omega_eff += multiplicity*Omega(d, θ)
#                 gamma_eff += multiplicity*Gamma(d, θ)
#             end
#         end
#     end
#     d::Float64 = 0.
#     for nz=1:N
#         d += b
#         omega_eff += 2*Omega(d, 0.)
#         gamma_eff += 2*Gamma(d, 0.)
#     end
#     return omega_eff, gamma_eff
# end

end # module


text:
include("../src/classification.jl")
using Main.Classification
using Random
using PyCall
using Statistics

metrics = pyimport("sklearn.metrics");
Random.seed!(5);
function traintestsplit(X, y, splitratio, portion)
    if portion > size(X, 1)
        @warn "Data set too small; using all"
        portion = size(X, 1)
    end
    M = shuffle(1:length(y))
    traintest = splitratio*portion |> x -> round(Int, x)
    Xtrain = X[M, :][1:traintest, :]
    Ytrain = y[M][1:traintest]
    Xtest  = X[M, :][traintest+1:portion, :]
    Ytest  = y[M][traintest+1:portion];
    (Xtrain, Ytrain), (Xtest, Ytest)
end
datasets = pyimport("sklearn.datasets")
cancer = datasets.load_breast_cancer()
wine = datasets.load_wine()
set = cancer
(Xtrain, Ytrain), (Xtest, Ytest) = traintestsplit(set["data"], set["target"], 0.8, 1000)

nagd = NAGDContext(learningrate=0.1, tolerance=1e-7, maxiterations=100000)
clf = LogisticClassifier(NesterovGradientDescent(nagd))
add_designmatrix!(clf, Xtrain)
addvalidationset!(clf, Xtest, Ytest)
@time fit!(clf, Ytrain)
@show clf.optimizer.converged
@show mean(Ytest .== predict(clf, Xtest))


text:
function ImageTransformations.restrict(f::AbstractPICDataStructure)
    data = resize_data(unwrapdata(f))
    grid = restrict(getdomain(f))
    newstruct(f, data, grid)
end

function resize_data(data)
    if hasunit(data)
        # use imresize directly on the underlying array
        apply_reinterpret(restrict, data)
    else
        restrict(data)
    end
end

function resize_data(data::StructArray, target_size...)
    StructArray(map(c->resize_data(c), components(data)))
end

function ImageTransformations.restrict(g::AbstractAxisGrid)
    grid_axes = collect(g)
    resized_axes = (resize_data(collect(g)) for g in grid_axes)
    AxisGrid(resized_axes..., names=propertynames(g))
end

# TODO: Investigate if sampling can be a better solution for particle data
function ImageTransformations.restrict(g::ParticlePositions)
    grid_axes = collect(g)
    resized_axes = (resize_data(g) for g in grid_axes)
    ParticlePositions(resized_axes..., names=propertynames(g))
end

function approx_target_size(::Type{T}) where T
    approx_target_size(domain_discretization(T), scalarness(T))
end
approx_target_size(::ParticleGrid, ::ScalarQuantity) = 7*10^5
approx_target_size(::ParticleGrid, ::VectorQuantity) = 70
approx_target_size(::LatticeGrid{1}, ::ScalarQuantity) = 600
approx_target_size(::LatticeGrid{2}, ::ScalarQuantity) = 400
approx_target_size(::LatticeGrid{3}, ::ScalarQuantity) = 200
approx_target_size(::LatticeGrid{2}, ::VectorQuantity) = 25
approx_target_size(::LatticeGrid{3}, ::VectorQuantity) = 15

"""
    downsample(f, target_size...)

Downsample the given input `f` to `target_size`. If no size is
provided, one is computed based on the type.
"""
function downsample(f, target_size...)
    all(size(f) .> target_size) ? downsample(restrict(f), target_size...) : f
end

function downsample(f; approx_size=nothing)
    T = typeof(f)
    if isnothing(approx_size)
        sz = approx_target_size(T)
    else
        sz = approx_size
    end
    largest_dim = maximum(size(f))
    factor = largest_dim ÷ sz
    factor == 0 && return f
    target_size = map(i->i÷factor, size(f))
    @debug "Resizing to size $target_size from $(size(f))"
    downsample(f, target_size...)
end


text:
"""

    PeriodicData{FT<:AbstractFloat}

Store information relevant to spatial periodicity in a system. `numdims` should be 0, 1, or
2, representing the number of periodic dimensions. `latticevecs` and `reciprocalvecs` should
each by 3-by-2 matrices, with each column representing a 3-element lattice or corresponding
reciprocal vector, with the vector being zeros along a dimension with no periodicity.

"""
struct PeriodicData{FT<:AbstractFloat}

    numdims::Integer
    latticevecs::Array{FT, 2}
    reciprocalvecs::Array{FT, 2}

end

"""

    readPeriodicData(FT::DataType=Float64, periodicfilename::AbstractString="")

Construct PeriodicData from an appropriate floating-point data type `FT` and a file of name
`periodicfilename` containing information about spatial periodicity; if no file is
specified, produce a default set of periodic data (no spatial periodicity at all).

"""
function readPeriodicData(FT::DataType=Float64, periodicfilename::AbstractString="")

    if (!(FT<:AbstractFloat) || isabstracttype(FT))
        error("readPeriodicData(): FT must be a concrete subtype of AbstractFloat");
    end

    if (periodicfilename == "")
        return PeriodicData(zero(Integer), zeros(FT, 3, 2), zeros(FT, 3, 2));
    end
    latticedata = vec(readdlm(periodicfilename, ',', FT));
    if (length(latticedata) == 3)
        return PeriodicData(one(Integer), hcat(latticedata, zero(latticedata)),
                            hcat(2*pi .* latticedata/dot(latticedata, latticedata),
                                 zero(latticedata)));
    elseif (length(latticedata) == 6)
        latticevecs = reshape(latticedata, 3, 2);
        eperp = cross(latticevecs[:,1], latticevecs[:,2]);
        Auc = norm(eperp);
        eperp = eperp ./ Auc;
        return PeriodicData(2, latticevecs,
                            (2*pi/Auc).*hcat(cross(latticevecs[:,2], eperp),
                                                   cross(eperp, latticevecs[:,1])));
    else
      error(string("readPeriodicData(): ", periodicfilename, " should have exactly 3 (for",
                   " 1D periodicity) or 6 (for 2D periodicity) real elements, each on a",
                   " separate line"));
    end
end

"""

    constructlatticefromreciprocal!(myPeriodicData::PeriodicData{<:AbstractFloat})

Reset the lattice vectors of `myPeriodicData` from its reciprocal vectors.

"""
function constructlatticefromreciprocal!(myPeriodicData::PeriodicData{<:AbstractFloat})

    if (myPeriodicData.numdims == 0)
        fill!(myPeriodicData.latticevecs, zero(eltype(myPeriodicData.latticevecs)));
    elseif (myPeriodicData.numdims == 1)
        myPeriodicData.latticevecs[:, 1] .= 2*pi*myPeriodicData.reciprocalvecs[:, 1] /
        dot(myPeriodicData.reciprocalvecs[:, 1], myPeriodicData.reciprocalvecs[:, 1]);
        myPeriodicData.latticevecs[:, 2] .= zeros(eltype(myPeriodicData.latticevecs), 3);
    elseif (myPeriodicData.numdims == 2)
        eperp = cross(myPeriodicData.reciprocalvecs[:, 1],
                      myPeriodicData.reciprocalvecs[:, 2]);
        ABZ = norm(eperp);
        eperp = eperp ./ ABZ;
        myPeriodicData.latticevecs .= (1/ABZ) .*
        hcat(cross(myPeriodicData.reciprocalvecs[:,2], eperp),
             cross(eperp, myPeriodicData.reciprocalvecs[:,1]));
    else
        error("PeriodicData: numdims must be 0, 1, or 2");
    end
end

"""

    constructreciprocalfromlattice!(myPeriodicData::PeriodicData{<:AbstractFloat})

Reset the reciprocal vectors of `myPeriodicData` from its lattice vectors.

"""
function constructreciprocalfromlattice!(myPeriodicData::PeriodicData{<:AbstractFloat})

    if (myPeriodicData.numdims == 0)
        fill!(myPeriodicData.reciprocalvecs, zero(eltype(myPeriodicData.latticevecs)));
    elseif (myPeriodicData.numdims == 1)
        myPeriodicData.reciprocalvecs[:, 1] .= 2*pi*myPeriodicData.latticevecs[:, 1] /
        dot(myPeriodicData.latticevecs[:, 1], myPeriodicData.latticevecs[:, 1]);
        myPeriodicData.reciprocalvecs[:, 2] .= zeros(eltype(myPeriodicData.reciprocalvecs),
                                                     3);
    elseif (myPeriodicData.numdims == 2)
        eperp = cross(myPeriodicData.latticevecs[:, 1], myPeriodicData.latticevecs[:, 2]);
        Auc = norm(eperp);
        eperp = eperp ./ Auc;
        myPeriodicData.reciprocalvecs .= (2*pi/Auc) .*
        hcat(cross(myPeriodicData.latticevecs[:, 2], eperp),
             cross(eperp, myPeriodicData.latticevecs[:, 1]));
    else
        error("PeriodicData: numdims must be 0, 1, or 2");
    end
end


text:
# Elimination data
struct KelData{S<:Number,E<:Number}
    s::Vector{S}
    e::Vector{E}
    a::Vector{Float64}
    b::Vector{Float64}
    r::Vector{Float64}
    ar::Vector{Float64}
    function KelData(s::Vector{S}, e::Vector{E}, a, b, r, ar)::KelData where S <: Number where E <: Number
        new{S, E}(s, e, a, b, r, ar)::KelData
    end
    function KelData()::KelData
        KelData(Float64[], Float64[], Float64[], Float64[], Float64[], Float64[])
    end
end

function Base.push!(keldata::KelData, s, e, a, b, r, ar)
    push!(keldata.s, s)
    push!(keldata.e, e)
    push!(keldata.a, a)
    push!(keldata.b, b)
    push!(keldata.r, r)
    push!(keldata.ar, ar)
end

function  Base.length(keldata::KelData)
    return length(keldata.a)
end

# Elimination settings
"""
    ElimRange(kelstart::Int, kelend::Int, kelexcl::Vector{Int})::ElimRange

Elimination settings for PK subject.

* `kelstart` - start point;
* `kelend` - end point;
* `kelexcl` - excluded points.
"""
mutable struct ElimRange{Symbol}
    kelstart::Int
    kelend::Int
    kelexcl::Vector{Int}
    function ElimRange(kelstart::Int, kelend::Int, kelexcl::Vector{Int}; time = false)::ElimRange
        if kelstart > kelend throw(ArgumentError("Kel start > kel end")) end
        if kelstart < 0 throw(ArgumentError("Kel start point < 0")) end
        if kelend   < 0 throw(ArgumentError("Kel endpoint < 0")) end
        if any(x -> x < 0, kelexcl) throw(ArgumentError("Exclude point < 0")) end
        if kelstart in kelexcl || kelend in kelexcl throw(ArgumentError("Kel start or kel end in exclusion")) end
        new{:point}(kelstart, kelend, kelexcl)::ElimRange
    end
    #=
    function ElimRange(kelstart, kelend)
        ElimRange(kelstart, kelend, Vector{Int}(undef, 0))
    end
    =#
    function ElimRange(;kelstart = 0, kelend = 0, kelexcl = Int[])
        ElimRange(kelstart, kelend, kelexcl)
    end
end

# Dose settings
"""
    DoseTime(dose::D, time::T, tau::TAU) where D <: Number where T <: Number where TAU <: Number

Dose settings.

* `dose` - dose;
* `time` - dose time;
* `tau` - tau (τ);

Dose time set 0 by default.
"""
struct DoseTime{D <: Number, T <: Number, TAU <: Number}
    dose::D
    time::T
    tau::TAU
    function DoseTime(dose::D, time::T, tau::TAU) where D <: Number where T <: Number where TAU <: Number
        if time < zero(T) throw(ArgumentError("Dose time can't be less zero!")) end
        new{D, T, TAU}(dose, time, tau)::DoseTime
    end
    function DoseTime(;dose = NaN, time = 0, tau = NaN)
        DoseTime(dose, time, tau)
    end
    #=
    function DoseTime(dose)
        DoseTime(dose, 0, NaN)
    end
    function DoseTime(dose, time)
        DoseTime(dose, time, NaN)
    end
    =#
end

# PK subject
mutable struct PKSubject{T <: Number, O <: Union{Number, Missing}, V <: Any} <: AbstractSubject
    time::Vector{T}
    obs::Vector{O}
    kelauto::Bool
    kelrange::ElimRange
    dosetime::DoseTime
    keldata::KelData
    id::Dict{Symbol, V}
    function PKSubject(time::Vector{T}, conc::Vector{O}, kelauto::Bool, kelrange::ElimRange, dosetime::DoseTime, keldata::KelData, sort::Dict{Symbol, V} = Dict{Symbol, Any}()) where T <: Number where O <: Union{Number, Missing} where V
        new{T, O, V}(time, conc, kelauto, kelrange, dosetime, keldata, sort)::PKSubject
    end
    function PKSubject(time::Vector, conc::Vector, kelauto::Bool, kelrange::ElimRange, dosetime::DoseTime, sort::Dict{Symbol, V}) where V
        PKSubject(time, conc, kelauto, kelrange, dosetime, KelData(), sort)
    end
    #=
    function PKSubject(time::Vector, conc::Vector, sort::Dict)
        PKSubject(time, conc, true, ElimRange(), DoseTime(NaN, 0), KelData(), sort)
    end
    function PKSubject(time::Vector, conc::Vector, kelauto::Bool, kelrange, dosetime; sort = Dict())
        PKSubject(time, conc, kelauto, kelrange, dosetime, KelData(), sort)
    end
    function PKSubject(time::Vector, conc::Vector; sort = Dict())
        PKSubject(time, conc, true, ElimRange(), DoseTime(NaN, 0), KelData(), sort)
    end
    =#
end

function Base.length(obj::T) where T <: AbstractSubject
    length(obj.time)
end

struct NCAResult{T} <: AbstractSubjectResult{T}
    data::T
    options::Dict{Symbol}
    result::Dict{Symbol, Float64}
    function NCAResult(subject::T, options, result) where T <: AbstractSubject
        new{T}(subject, options, result)
    end
    #=
    function NCAResult(subject::T, method, result) where T <: AbstractSubject
        NCAResult(subject, method, result, Dict())
    end
    =#
end

"""
    LimitRule(lloq::T, btmax, atmax, nan, rm::Bool) where T <: Real

    LimitRule(;lloq = NaN, btmax = NaN, atmax = NaN, nan = NaN, rm::Bool = false)

* `lloq` - LLOQ - low limit of quantification;
* `btmax` - value for points before Tmax;
* `atmat` - values for points after Tmax;
* `nan` - values for replacing `NaN`;
* `rm` - if `true`, removee all `NaN` points.

Rule for PK subject.

* STEP 1 (NaN step): replace all `NaN` and `missing` values with nan keyword value (if `nan` not NaN);
* STEP 2 (LLOQ step): replace values below `lloq` with `btmax` value if this value befor Tmax or with atmax if this value after Tmax (if `lloq` not NaN);
* STEP 3 (remove NaN): `rm` == true, then remove all `NaN` and `missing` values.
"""
struct LimitRule{T<:Real}
    lloq::T
    btmax::Float64
    atmax::Float64
    nan::Float64
    rm::Bool
    function LimitRule(lloq::T, btmax, atmax, nan, rm::Bool) where T <: Real
        new{T}(lloq, btmax, atmax, nan, rm)::LimitRule
    end
    function LimitRule(;lloq = NaN, btmax = NaN, atmax = NaN, nan = NaN, rm::Bool = false)
        LimitRule(lloq, btmax, atmax, nan, rm)
    end
end

function isapplicable(lr::LimitRule)
    !isnan(lr.lloq) || !isnan(lr.nan) || lr.rm ? true : false
end

#Urine PK subject
mutable struct UPKSubject{T <: Tuple{Number, Number}, O <: Union{Number, Missing}, VOL <: Union{Number, Missing}, V <: Any} <: AbstractSubject
    time::Vector{T}
    obs::Vector{O}
    vol::Vector{VOL}
    kelauto::Bool
    kelrange::ElimRange
    dosetime::DoseTime
    keldata::KelData
    id::Dict{Symbol, V}
    function UPKSubject(time::Vector{T}, conc::Vector{O}, vol::Vector{VOL}, kelauto::Bool, kelrange::ElimRange, dosetime::DoseTime, keldata::KelData, id::Dict{Symbol, V} = Dict{Symbol, Any}()) where T <: Tuple{Number, Number} where O <: Union{Number, Missing} where VOL <: Union{Number, Missing} where V
        new{T, O, VOL, V}(time, conc, vol, kelauto, kelrange, dosetime, keldata, id)
    end
    function UPKSubject(time::Vector, conc::Vector, vol::Vector, kelauto::Bool, kelrange::ElimRange, dosetime::DoseTime, id::Dict{Symbol, V}) where V
        UPKSubject(time, conc, vol, kelauto, kelrange, dosetime, KelData(), id)
    end
end

# PK subject
mutable struct PDSubject{T <: Number, O <: Union{Number, Missing}, V <: Any} <: AbstractSubject
    time::Vector{T}
    obs::Vector{O}
    bl::Float64
    th::Float64
    id::Dict{Symbol, V}
    function PDSubject(time::Vector{T}, conc::Vector{O}, bl, th, sort::Dict{Symbol, V} = Dict{Symbol, Any}()) where T <: Number where O <: Union{Number, Missing} where V
        new{T, O, V}(time, conc, bl, th, sort)::PDSubject
    end
    function PDSubject(time::Vector, conc::Vector, bl, th, sort::Dict{Symbol, V}) where V
        PDSubject(time, conc, bl, th, sort)
    end
end

struct NCAOptions{LR <: Union{Nothing, LimitRule}}
    adm::Symbol
    calcm::Symbol
    intpm::Symbol
    limitrule::LR
    verbose::Int
    warn::Bool
    io::IO
    modify!::Function
end

#=
function gettime(subj::PKSubject)
end
function getobs(subj::PKSubject)
end
function gettime(subj::UPKSubject)
end
function getobs(subj::UPKSubject)
end
=#


text:
import Dolo

@testset "features" begin

    path = Dolo.pkg_path

    examples_path = joinpath(path, "examples", "models")
    readdir(examples_path)

    model = Dolo.yaml_import(joinpath(examples_path, "rbc_iid.yaml"))
    ff = Dolo.features(model)
    @test ff[:one_dimensional] == false
    @test ff[:one_state] == false
    @test ff[:one_control] == false
    @test ff[:nonstochastic_transitions] == false
    @test ff[:bounds_are_constant] == true

    model = Dolo.yaml_import(joinpath(examples_path, "neoclassical.yaml"))
    ff = Dolo.features(model)
    @test ff[:one_dimensional] == true
    @test ff[:one_state] == true
    @test ff[:one_control] == true
    @test ff[:nonstochastic_transitions] == true
    @test ff[:bounds_are_constant] == true

    model = Dolo.yaml_import(joinpath(examples_path, "rbc_mc.yaml"))
    ff = Dolo.features(model)
    @test ff[:one_dimensional] == true
    @test ff[:one_state] == true
    @test ff[:one_control] == false
    @test ff[:nonstochastic_transitions] == true
    @test ff[:bounds_are_constant] == true

end


text:
module CUTENSOR

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cudaDataType
using ..CUDA: libcutensor,  @retry_reclaim

using CEnum

using Memoization

using DataStructures

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,Base.RefValue{cutensorHandle_t}}()

function handle()
    ctx = context()
    get!(task_local_storage(), (:CUTENSOR, ctx)) do
        handle = pop!(idle_handles, ctx) do
            handle = Ref{cutensorHandle_t}()
            cutensorInit(handle)
            handle
        end

        finalizer(current_task()) do task
            push!(idle_handles, ctx, handle) do
                # CUTENSOR doesn't need to actively destroy its handle
            end
        end

        handle
    end::Base.RefValue{cutensorHandle_t}
end

end


text:
import Base.length, Base.getindex, Base.iterate, Base.view

function view(o::GBIFRecords)
    view(o.occurrences, o.show)
end

function length(o::GBIFRecords)
    length(view(o))
end

function getindex(o::GBIFRecords, i::Int64)
    view(o)[i]
end

function getindex(o::GBIFRecords, r::UnitRange{Int64})
    view(o)[r]
end

function iterate(o::GBIFRecords)
    iterate(collect(view(o)))
end

function iterate(o::GBIFRecords, t::Union{Int64,Nothing})
    iterate(collect(view(o)), t)
end


text:
# Types for module InterestRates

"""
The type `DayCountConvention` sets the convention on how to count the number of days between dates, and also how to convert that number of days into a year fraction.

Given an initial date `D1` and a final date `D2`, here's how the distance between `D1` and `D2` are mapped into a year fraction for each supported day count convention:

* *Actual360* : `(D2 - D1) / 360`
* *Actual365* : `(D2 - D1) / 365`
* *BDays252* : `bdays(D1, D2) / 252`, where `bdays` is the business days between `D1` and `D2` from [BusinessDays.jl package](https://github.com/felipenoris/BusinessDays.jl).
"""
abstract type DayCountConvention end
struct Actual360 <: DayCountConvention end
struct Actual365 <: DayCountConvention end

mutable struct BDays252 <: DayCountConvention
    hc::HolidayCalendar
end

"""
The type `CompoundingType` sets the convention on how to convert a yield into an Effective Rate Factor.

Given a yield `r` and a maturity year fraction `t`, here's how each supported compounding type maps the yield to Effective Rate Factors:

* *ContinuousCompounding* : `exp(r*t)`
* *SimpleCompounding* : `(1+r*t)`
* *ExponentialCompounding* : `(1+r)^t`
"""
abstract type CompoundingType end
struct ContinuousCompounding <: CompoundingType end   # exp(r*t)
struct SimpleCompounding <: CompoundingType end       # (1+r*t)
struct ExponentialCompounding <: CompoundingType end  # (1+r)^t

"""
This package provides the following curve methods.

**Interpolation Methods**

* **Linear**: provides Linear Interpolation on rates.
* **FlatForward**: provides Flat Forward interpolation, which is implemented as a Linear Interpolation on the *log* of discount factors.
* **StepFunction**: creates a step function around given data points.
* **CubicSplineOnRates**: provides *natural cubic spline* interpolation on rates.
* **CubicSplineOnDiscountFactors**: provides *natural cubic spline* interpolation on discount factors.
* **CompositeInterpolation**: provides support for different interpolation methods for: (1) extrapolation before first data point (`before_first`), (2) interpolation between the first and last point (`inner`), (3) extrapolation after last data point (`after_last`).

For *Interpolation Methods*, the field `dtm` holds the number of days between `date` and the maturity of the observed yield, following the curve's day count convention, which must be given in advance, when creating an instance of the curve. The field `zero_rates` holds the yield values for each maturity provided in `dtm`. All yields must be anual based, and must also be given in advance, when creating the instance of the curve.

**Term Structure Models**

* **NelsonSiegel**: term structure model based on *Nelson, C.R., and A.F. Siegel (1987), Parsimonious Modeling of Yield Curve, The Journal of Business, 60, 473-489*.
* **Svensson**: term structure model based on *Svensson, L.E. (1994), Estimating and Interpreting Forward Interest Rates: Sweden 1992-1994, IMF Working Paper, WP/94/114*.

For *Term Structure Models*, the field `parameters` holds the constants defined by each model, as described below. They must be given in advance, when creating the instance of the curve.

For **NelsonSiegel** method, the array `parameters` holds the following parameters from the model:
* **beta1** = parameters[1]
* **beta2** = parameters[2]
* **beta3** = parameters[3]
* **lambda** = parameters[4]

For **Svensson** method, the array `parameters` hold the following parameters from the model:
* **beta1** = parameters[1]
* **beta2** = parameters[2]
* **beta3** = parameters[3]
* **beta4** = parameters[4]
* **lambda1** = parameters[5]
* **lambda2** = parameters[6]

### Methods hierarchy

As a summary, curve methods are organized by the following hierarchy.

* `<<CurveMethod>>`
    * `<<Interpolation>>`
        * `<<DiscountFactorInterpolation>>`
            * `CubicSplineOnDiscountFactors`
            * `FlatForward`
        * `<<RateInterpolation>>`
            * `CubicSplineOnRates`
            * `Linear`
            * `StepFunction`
        * `CompositeInterpolation`
    * `<<Parametric>>`
        * `NelsonSiegel`
        * `Svensson`
"""
abstract type CurveMethod end
abstract type Parametric <: CurveMethod end
abstract type Interpolation <: CurveMethod end

abstract type DiscountFactorInterpolation <: Interpolation end
abstract type RateInterpolation <: Interpolation end

struct CubicSplineOnRates <: RateInterpolation end
struct CubicSplineOnDiscountFactors <: DiscountFactorInterpolation end
struct FlatForward <: DiscountFactorInterpolation end
struct Linear <: RateInterpolation end
struct NelsonSiegel <: Parametric end
struct Svensson <: Parametric end
struct StepFunction <: RateInterpolation end

mutable struct CompositeInterpolation <: Interpolation
    before_first::Interpolation # Interpolation method to be applied before the first point
    inner::Interpolation
    after_last::Interpolation # Interpolation method to be applied after the last point
end

# Helper functions to identify cubic spline method
is_cubic_spline_on_rates(m::CurveMethod) = false
is_cubic_spline_on_rates(m::CubicSplineOnRates) = true
is_cubic_spline_on_rates(m::CompositeInterpolation) = is_cubic_spline_on_rates(m.before_first) || is_cubic_spline_on_rates(m.inner) || is_cubic_spline_on_rates(m.after_last)
is_cubic_spline_on_discount_factors(m::CurveMethod) = false
is_cubic_spline_on_discount_factors(m::CubicSplineOnDiscountFactors) = true
is_cubic_spline_on_discount_factors(m::CompositeInterpolation) = is_cubic_spline_on_discount_factors(m.before_first) || is_cubic_spline_on_discount_factors(m.inner) || is_cubic_spline_on_discount_factors(m.after_last)

abstract type AbstractIRCurve end

# Helper function to create splinefit results for method CubicSplineOnDiscountFactors
function _splinefit_discountfactors(curve::AbstractIRCurve)
    dtm_vec = curve_get_dtm(curve)
    curve_rates_vec = curve_get_zero_rates(curve)
    l = length(dtm_vec)
    yf_vec = Vector{Float64}(undef, l)
    discount_vec = Vector{Float64}(undef, l)

    for i = 1:l
        @inbounds yf_vec[i] = dtm_vec[i] / daysperyear(curve_get_daycount(curve))
        @inbounds discount_vec[i] = discountfactor(curve_get_compounding(curve), curve_rates_vec[i], yf_vec[i])
    end

    return splinefit(yf_vec, discount_vec)
end

mutable struct IRCurve <: AbstractIRCurve
    name::String
    daycount::DayCountConvention
    compounding::CompoundingType
    method::CurveMethod
    date::Date
    dtm::Vector{Int} # for interpolation methods, stores days_to_maturity on curve's daycount convention.
    zero_rates::Vector{Float64} # for interpolation methods, parameters[i] stores yield for maturity dtm[i].
    parameters::Vector{Float64} # for parametric methods, parameters stores model's constant parameters.
    dict::Dict{Symbol, Any}     # holds pre-calculated values for optimization, or additional parameters.

    # Constructor for Interpolation methods
    function IRCurve(name::AbstractString, _daycount::DayCountConvention,
        compounding::CompoundingType, method::M,
        date::Date, dtm::Vector{Int},
        zero_rates::Vector{Float64}, parameters = Vector{Float64}(), dict = Dict{Symbol, Any}()) where {M<:Interpolation}

        @assert !isempty(dtm) "Empty days-to-maturity vector"
        @assert !isempty(zero_rates) "Empty zero_rates vector"
        @assert length(dtm) == length(zero_rates) "dtm and zero_rates must have the same length"
        @assert issorted(dtm) "dtm should be sorted before creating IRCurve instance"

        new_curve = new(String(name), _daycount, compounding, method, date, dtm, zero_rates, parameters, dict)

        # Stores splinefit results for Cubic Spline methods
        if is_cubic_spline_on_rates(method)
            sp = splinefit(curve_get_dtm(new_curve), curve_get_zero_rates(new_curve))
            curve_set_dict_parameter!(new_curve, :spline_fit_on_rates, sp)
        end

        if is_cubic_spline_on_discount_factors(method)
            sp = _splinefit_discountfactors(new_curve)
            curve_set_dict_parameter!(new_curve, :spline_fit_on_discount_factors, sp)
        end

        return new_curve
    end

    # Constructor for Parametric methods
    function IRCurve(name::AbstractString, _daycount::DayCountConvention,
        compounding::CompoundingType, method::M,
        date::Date,
        parameters::Vector{Float64},
        dict = Dict{Symbol, Any}()) where {M<:Parametric}

        @assert !isempty(parameters) "Empty yields vector"
        new(String(name), _daycount, compounding, method, date, Vector{Int}(), Vector{Float64}(), parameters, dict)
    end
end

# Interface for curve types
"""
    curve_get_name(curve::AbstractIRCurve) → String

Returns the name of the curve.
"""
curve_get_name(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_daycount(curve::AbstractIRCurve) → DayCountConvention

Returns the DayCountConvention used by the curve. See DayCountConvention documentation.
"""
curve_get_daycount(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_compounding(curve::AbstractIRCurve) → CompoundingType

Returns the CompoundingType used by the curve. See CompoundingType documentation.
"""
curve_get_compounding(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_method(curve::AbstractIRCurve) → CurveMethod

Returns the CurveMethod used by the curve. See CurveMethod documentation.
"""
curve_get_method(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_date(curve::AbstractIRCurve) → Date

Returns the date when the curve is observed. All zero rate calculation will be performed based on this date.
"""
curve_get_date(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_dtm(curve::AbstractIRCurve) → Vector{Int}

Used for interpolation methods, returns `days_to_maturity` on curve's daycount convention.
"""
curve_get_dtm(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_zero_rates(curve::AbstractIRCurve) → Vector{Float64}

Used for interpolation methods, parameters[i] returns yield for maturity dtm[i].
"""
curve_get_zero_rates(curve::AbstractIRCurve) = error("method not defined")

"""
    curve_get_model_parameters(curve::AbstractIRCurve) → Vector{Float64}

Used for parametric methods, returns model's constant parameters.
"""
curve_get_model_parameters(curve::AbstractIRCurve) = error("method not defined")

# AbstractIRCurve interface implementation for IRCurve type
curve_get_name(curve::IRCurve) = curve.name
curve_get_daycount(curve::IRCurve) = curve.daycount
curve_get_compounding(curve::IRCurve) = curve.compounding
curve_get_method(curve::IRCurve) = curve.method
curve_get_date(curve::IRCurve) = curve.date
curve_get_dtm(curve::IRCurve) = curve.dtm
curve_get_zero_rates(curve::IRCurve) = curve.zero_rates
curve_get_model_parameters(curve::IRCurve) = curve.parameters
curve_get_dict_parameter(curve::IRCurve, sym::Symbol) = curve.dict[sym]

function curve_set_dict_parameter!(curve::IRCurve, sym::Symbol, value)
    curve.dict[sym] = value
end

function curve_get_spline_fit_on_rates(curve::IRCurve) :: Spline
    @assert is_cubic_spline_on_rates(curve_get_method(curve)) "Curve $(curve_get_name(curve)) with method $(curve_get_method(curve)) does not hold a spline_fit_on_rates result."
    return curve_get_dict_parameter(curve, :spline_fit_on_rates)
end

function curve_get_spline_fit_on_discount_factors(curve::IRCurve) :: Spline
    @assert is_cubic_spline_on_discount_factors(curve_get_method(curve)) "Curve $(curve_get_name(curve)) with method $(curve_get_method(curve)) does not hold a spline_fit_on_discount_factors result."
    return curve_get_dict_parameter(curve, :spline_fit_on_discount_factors)
end


text:
# # [Symbolic Regression of Nonlinear Time Continuous Systems](@id symbolic_regression_simple)
#
# !!! note
#
#   Symbolic regression is using regularized evolution, simulated annealing, and gradient-free optimization to find suitable equations.
#   Hence, the performance might differ and depends strongly on the hyperparameters of the optimization.
#   This example might not recover the groundtruth, but is showing off the use within `DataDrivenDiffEq.jl`.
#
# DataDrivenDiffEq offers an interface to [`SymbolicRegression.jl`](https://github.com/MilesCranmer/SymbolicRegression.jl) to infer more complex functions. To
# use it, simply load a sufficient version of `SymbolicRegression` (currently we support version >= 0.9).

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using SymbolicRegression
#md using Plots

A = [-0.9 0.2; 0.0 -0.5]
B = [0.0; 1.0]
u0 = [10.0; -10.0]
tspan = (0.0, 20.0)

f(u,p,t) = A*u .+ B .* sin(0.5*t)

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.01);

# We will use the data provided by our problem, but add the control signal `U = sin(0.5*t)` to it. Instead of using a function, like in [another example](@ref linear_continuous_controls)
X = Array(sol)
t = sol.t
U = permutedims(sin.(0.5*t))
prob = ContinuousDataDrivenProblem(X, t, U = U)

# And plot the problems data.

#md plot(prob)

# To solve our problem, we will use [`EQSearch`](@ref), which provides a wrapper for the symbolic regression interface.
# By default, it takes in a `Vector` of `Functions` and additional [keyworded arguments](https://astroautomata.com/SymbolicRegression.jl/v0.6/api/#Options). We will stick to simple operations
# like subtraction and multiplication, use a `L1DistLoss` , limit the maximum size and punish complex equations while fitting our equations on minibatches.

alg = EQSearch([-, *], loss = L1DistLoss(), verbosity = 0, maxsize = 9, batching = true, batchSize = 50, parsimony = 0.01f0)

# Again, we `solve` the problem to obtain a [`DataDrivenSolution`](@ref). Note that any additional keyworded arguments are passed onto
# symbolic regressions [`EquationSearch`](https://astroautomata.com/SymbolicRegression.jl/v0.6/api/#EquationSearch) with the exception of `niterations` which
# is `max_iter`

res = solve(prob, alg, max_iter = 300, numprocs = 0, multithreading = false)
#md println(res)

# We see that the system has been recovered correctly, indicated by the small error. A closer look at the equations r

system = result(res)
#md
println(system)

# Shows that while not obvious, the representation
# And also plot the prediction of the recovered dynamics

#md plot(res)

# To convert the result into an `ODESystem`, we substitute the control signal

u = controls(system)
t = get_iv(system)

subs_control = (u[1] => sin(0.5*t))

eqs = map(equations(system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

@named sys = ODESystem(
    eqs,
    get_iv(system),
    states(system),
    []
    );

# And simulated using `OrdinaryDiffEq.jl` using the (known) initial conditions and the parameter mapping of the estimation.
# Since the parameters are *hard numerical values* we do not need to include those.

x = states(system)
x0 = [x[1] => u0[1], x[2] => u0[2]]

ode_prob = ODEProblem(sys, x0, tspan)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

# And look at the result
#md plot(sol, color = :black)
#md plot!(estimate, color = :red, linestyle = :dash)

#md # ## [Copy-Pasteable Code](@id symbolic_regression_simple_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

@test all(l2error(res) .<= 5e-1) #src
@test Array(sol) ≈ Array(estimate) rtol = 5e-2 #src


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function airtovac(wave_air::AbstractFloat)
    if wave_air >= 2000
        wave_vac = wave_air
        for iter= 1:2
            sigma2 = (1e4/wave_vac)^2 # Convert to wavenumber squared
            # Computer conversion factor.
            fact = 1 + 5.792105e-2/(238.0185 - sigma2) +
                1.67917e-3/(57.362 - sigma2)
            wave_vac = wave_air*fact # Convert Wavelength
        end
        return wave_vac
    else
        return wave_air
    end
end

"""
    airtovac(wave_air) -> wave_vacuum

### Purpose ###

Converts air wavelengths to vacuum wavelengths.

### Explanation ###

Wavelengths are corrected for the index of refraction of air under standard conditions.
Wavelength values below \$2000 Å\$ will *not* be altered, take care within \$[1 Å, 2000
Å]\$.  Uses relation of Ciddor (1996).

### Arguments ###

* `wave_air`: the wavelength in air.

### Output ###

Vacuum wavelength in angstroms.

### Method ###

Uses relation of Ciddor (1996), Applied Optics 62, 958.

### Example ###

If the air wavelength is `w = 6056.125` (a Krypton line), then `airtovac(w)`
yields a vacuum wavelength of `6057.8019`.

```jldoctest
julia> using AstroLib

julia> airtovac(6056.125)
6057.801930991426
```

### Notes ###

`vactoair` converts vacuum wavelengths to air wavelengths.

Code of this function is based on IDL Astronomy User's Library.
"""
airtovac(wave_air::Real) = airtovac(float(wave_air))


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Test

# =========================================================================================
# Test tess_run_pipeline(Callback) /w no renderers.
@testset "tess_pipeline(Callback) /w no renderers" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local err      = "No renderers attached to the pipeline."

    @test (@test_logs (:error, err) tess_run_pipeline(pipeline) do add
        add(pix_with([Line("This is image one.")]), 72)
    end) == false
end

# =========================================================================================
# Test tess_pipeline(Callback) delete instance before run.
@testset "tess_pipeline(Callback) delete instance before run" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)

    @test tess_pipeline_text(pipeline) != nothing

    tess_delete!(inst)

    @suppress begin
        @test (tess_run_pipeline(pipeline) do add
            @test add(pix_with([Line("This is image one.")]), 72) == false
            return true
        end) == false
    end
end

# =========================================================================================
# Test tess_pipeline(Callback) delete instance inside pipeline.
@testset "tess_pipeline(Callback) delete instance inside pipeline" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)

    @test tess_pipeline_text(pipeline) != nothing

    @suppress begin
        @test (tess_run_pipeline(pipeline) do add
            @test add(pix_with([Line("This is image one.")]), 72) == true
            tess_delete!(inst)
            @test add(pix_with([Line("This is image two.")]), 72) == false
        end) == false
    end
end

# =========================================================================================
# Test tess_pipeline(Callback) delete pipeline inside pipeline.
@testset "tess_pipeline(Callback) delete pipeline inside pipeline" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)

    @test tess_pipeline_text(pipeline) != nothing

    @suppress begin
        @test (tess_run_pipeline(pipeline) do add
            @test add(pix_with([Line("This is image one.")]), 72) == true
            tess_delete!(pipeline)
            @test add(pix_with([Line("This is image two.")]), 72) == false
        end) == false
    end
end

# =========================================================================================
# Test tess_pipeline(Callback) /w failure.
@testset "tess_pipeline(Callback) /w failure" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)

    @test tess_pipeline_text(pipeline) != nothing

    @test (tess_run_pipeline(pipeline) do add
        @test add(pix_with([Line("This is image one.")]), 72) == true
        @test add(pix_with([Line("This is image two.")]), 72) == true
        return false
    end) == false
end


text:
using Blog

Blog.make()

max_stars_repo_name: piever/Blog

max_issues_repo_name: piever/Blog

max_forks_repo_name: piever/Blog


text:
# by Lazaro Alonso
using CairoMakie
let
    x = 0:0.05:4π
    fig = Figure(resolution = (600,400), font = "CMU Serif") # probably you need to install this font in your system
    ax = Axis(fig, xlabel = L"x", ylabel = L"f (x)", ylabelsize = 22,
        xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
        xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10)
    lines!(x, x -> sin(3x) / (cos(x) + 2)/x, label = L"\frac{\sin(3x)}{x(\cos(x) + 2)}")
    lines!(x, x-> cos(x)/x, label = L"\cos(x)/x")
    lines!(x, x-> exp(-x), label = L"e^{-x}")
    ylims!(-0.6,1.05)
    xlims!(-0.5,12)
    axislegend(L"f(x)", position = :rt, bgcolor = (:grey90, 0.25))
    fig[1,1] = ax
    #save("filename.png", fig, px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "latexSample.png"), fig, px_per_unit = 2.0) # HIDE
    display(fig)
end
using Pkg # HIDE
Pkg.status() # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Documenter, MLDatasets

makedocs(
    modules = [MLDatasets],
    clean = false,
    format = Documenter.HTML(
                prettyurls = haskey(ENV, "CI"),
                assets = [joinpath("assets", "favicon.ico")]
            ),
    sitename = "MLDatasets.jl",
    authors = "Hiroyuki Shindo, Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        "Available Datasets" => Any[
            "Image Classification" => Any[
                "MNIST handwritten digits" => "datasets/MNIST.md",
                "Fashion MNIST" => "datasets/FashionMNIST.md",
                "CIFAR-10" => "datasets/CIFAR10.md",
                "CIFAR-100" => "datasets/CIFAR100.md",
                "SVHN format 2" => "datasets/SVHN2.md",
            ],
        ],
        hide("Indices" => "indices.md"),
        "LICENSE.md",
    ]
)

deploydocs(repo = "github.com/JuliaML/MLDatasets.jl.git")


text:
mutable struct Counters
	"""linear counter for Separable_ExogStates"""
	i_exogstate::Int64
	"""linear counter for Separable_States"""
	i_state::Int64
	"""linear counter for Separable"""
	i_statechoice::Int64
	"""linear counter for choice"""
	i_choice::Int64
end

initialize_counter() = Counters(0,0,0,0)
getcounter(cnt::Counters, intdim::Type{Separable_ExogStates}) = cnt.i_exogstate
getcounter(cnt::Counters, intdim::Type{Separable_States}) = cnt.i_state
getcounter(cnt::Counters, intdim::Type{Separable}) = cnt.i_statechoice
function reset!(cnt::Counters)
	cnt.i_exogstate = 0
	cnt.i_state = 0
	cnt.i_statechoice = 0
	cnt.i_choice = 0
end

getreward(rewardcall::Type{pre_partial}, rewardfunc, mReward, tStateVectors, vChoices,
	cnt, ix, j, jprime) =
	rewardfunc(mReward[cnt.i_state], getindex.(tStateVectors, (j, ix.I...)), vChoices[jprime])

getreward(rewardcall::Type{jit}, rewardfunc, mReward, tStateVectors, vChoices,
	cnt, ix, j, jprime) =
	rewardfunc(getindex.(tStateVectors, (j, ix.I...)), vChoices[jprime])

getreward(rewardcall::Type{pre}, rewardfunc, mReward, tStateVectors, vChoices,
	cnt, ix, j, jprime) =  mReward[jprime, cnt.i_state] # nChoices x nStates

# function getreward(rewardcall, rewardfunc, mReward, tStateVectors, vChoices,
# 	cnt, ix, j, jprime)
# 	if rewardcall == :pre_partial
# 		return rewardfunc(mReward[cnt.i_state], getindex.(tStateVectors, (j, ix.I...)), vChoices[jprime])
# 	elseif rewardcall == :jit
# 		# need to be VERY careful with order of state vars here..
# 		return rewardfunc(getindex.(tStateVectors, (j, ix.I...)), vChoices[jprime])
# 	elseif rewardcall == :pre
# 		return mReward[jprime, cnt.i_state] # nChoices x nStates
# 	end
# end

########################
##### Solver ######
########################
_solve(p::DDP{NS,1},
	mTransition::Array{Float64,2}, mReward::Union{Array{Float64}, Nothing},
	opts::SolverOptions) where NS =
	_solve(p, mTransition, mReward, opts, p.tStateVectors,
		getchoicevars(p.tStateVectors, p.tChoiceVectors)[1],
		getnonchoicevars(p.tStateVectors, p.tChoiceVectors),)

# split _solve for function barrier for tuples
# Note: really hard to create functions for the innner loops without getting
# a performance hit. Possible to write a general maximum function that does
# concavity, but it's slower if outside the loop somehow. Didn't try gathering
# inputs into a struct though.
function _solve(p::DDP{NS,1},
				mTransition::Array{Float64,2}, mReward::Union{Array{Float64}, Nothing},
				opts::SolverOptions,
				tStateVectors,
				vChoices::Vector{Float64},
				tOtherStateVectors) where NS

	@unpack β, rewardfunc = p
	@unpack disp, disp_each_iter, max_iter, epsilon, rewardcall,
		monotonicity, concavity, intdim = opts

	nChoices = length(vChoices)

	nStates = prod(length.(tStateVectors))
	nOtherStates = prod(length.(tOtherStateVectors))

	mValFun    = zeros((nChoices, nOtherStates)) # matrix so that can multiply with mTransition
	mValFunNew = zeros(nStates) # vector so that can index easily

	mPolFunInd = zeros(Int16, nStates)

	mβEV = zeros(nChoices, size(mTransition, 1)) # depends on intdim

	# VFI initialization
	maxDifference::Float64 = Inf # to initialize
	iteration::Int64 = 0 # initialize counter
	if β == 0.0
	    tolerance = Inf
	elseif β == 1.0
	    throw(ArgumentError("method invalid for beta = 1"))
	else
	    tolerance = epsilon * (1-β) / (2*β)
	end

	# initialize for less memory allocation
	mValFunDiff = zeros(nStates)
	mValFunDiffAbs = zeros(nStates)

	# liquidationvalue::Float64 = -10000.
	# inactionvalue::Float64 = -10000.
	reward::Float64 = -Inf

	valueHighSoFar::Float64 = -Inf
	valueProvisional::Float64 = -Inf
	iChoiceStart::Int16 = 1
	iChoice::Int16 = 1

	cnt = initialize_counter()

	# will need beta times transpose of transition matrix
	mTransition_βT = β * transpose(mTransition)

	# VFI
	while maxDifference > tolerance

	mul!(mβEV, mValFun, mTransition_βT)

	reset!(cnt)

		# @inbounds
		for ix in CartesianIndices(length.(tOtherStateVectors)) # other states

			cnt.i_exogstate += 1

		    if monotonicity
		       iChoiceStart = 1
		    end

		    for j = 1:nChoices # first state

		        valueHighSoFar = -Inf
		        iChoice  = 1

				cnt.i_state += 1

				cnt.i_statechoice += iChoiceStart-1 # compensate if start later

		        for jprime = iChoiceStart:nChoices

					cnt.i_statechoice += 1

					reward = getreward(rewardcall, rewardfunc, mReward, tStateVectors, vChoices,
						cnt, ix, j, jprime)

		            valueProvisional = reward + mβEV[jprime, getcounter(cnt, intdim)]

		            if valueProvisional>valueHighSoFar
		            	valueHighSoFar = valueProvisional
		            	iChoice = jprime
		            elseif concavity
						cnt.i_statechoice += nChoices - jprime # adjust counter if finish early
		        	    break # break bc value will only be lower from now on
					end

		        end #jprime

		        mValFunNew[cnt.i_state] = valueHighSoFar
		        mPolFunInd[cnt.i_state] = iChoice

				if monotonicity
                    iChoiceStart = iChoice
				end

		    end #j

		end #i

	mValFunDiff .= mValFunNew .- @view(mValFun[:])
	mValFunDiffAbs .= abs.(mValFunDiff)
	maxDifference  = maximum(mValFunDiffAbs)

	# mcqueen!(mValFunNew, mValFun, p.β)

	copyto!(mValFun, mValFunNew)

	# howards improvement?

	if disp
		display_iter(iteration, disp_each_iter, maxDifference)
	end

	iteration += 1

	if iteration > max_iter
	    println("WARNING: $max_iter maximum iterations exceeded")
	    # println("Parameters used:")
		# println(p.params)
	    break
	end

	end #while

	# println("Iteration ", string(iteration), ": Sup change is = ", string(maxDifference))

	mPolFun = vChoices[mPolFunInd]

	nNodes = length.(tStateVectors)
	meshPolFun = reshape(mPolFun, tuple(nNodes...))
	meshValFun = reshape(mValFun, tuple(nNodes...))
	# meshExit   = reshape(mExit, tuple(p.nNodes...))

	meshValFun, (meshPolFun,)

end # _solve

# function mcqueen!(v, vold, beta)
#     b_l = beta /(1-beta) * minimum(v-vold);
#     b_u = beta /(1-beta) * maximum(v-vold);
#     # @show (b_l+b_u)/2
#     v[:] = v[:] + (b_l+b_u)/2;
#     nothing
# end


text:
# illustrates effects of collinearity
# try setting rho=0.1 or rho=0.9, and observe
# the effect on standard dev of coefficients
using Econometrics, Plots, LinearAlgebra
function wrapper(rho, n)
  V = [1.0 rho; rho 1.0]
  P = cholesky(V).U
  x = randn(n,2)*P
  x = [ones(n,1) x]
  b = ones(3)
  y = x*b + randn(n,1)
  b1 = x\y

  z = x[:,1:2]
  b2 = z\y

  R = [0.0 1.0 -1.0]
  r = 0.0
  b3 = ols(y, x, R=R, r=r, silent=true)[1]
  results = [b1' b2' b3']
end

function main()
reps = 1000
rho = 0.9
n = 30
data = zeros(reps,8)
for rep = 1:reps
    data[rep,:] = wrapper(rho,n)
end
println("correlation between x2 and x3: ",rho)
println("descriptive statistics for 1000 OLS replications")
dstats(data[:,1:3])
println("descriptive statistics for 1000 OLS replications, dropping x3")
dstats(data[:,4:5])
println("descriptive statistics for 1000 Restricted OLS replications, b2=b3")
dstats(data[:,6:8])
p1 = histogram(data[:,2],title="OLS")
p2 = histogram(data[:,5],title="Drop X3")
p3 = histogram(data[:,7],title="Restricted LS")
plot(p1,p2,p3, layout=(3,1), show=true, legend=false)
gui()
savefig("collin.png")
return
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
immutable NoncentralChisq <: ContinuousUnivariateDistribution
    ν::Float64
    λ::Float64
    function NoncentralChisq(ν::Real, λ::Real)
        @check_args(NoncentralChisq, ν > zero(ν))
        @check_args(NoncentralChisq, λ >= zero(λ))
    	new(ν, λ)
    end
end

@distr_support NoncentralChisq 0.0 Inf

### Parameters

params(d::NoncentralChisq) = (d.ν, d.λ)

### Statistics

mean(d::NoncentralChisq) = d.ν + d.λ
var(d::NoncentralChisq) = 2.0*(d.ν + 2.0*d.λ)
skewness(d::NoncentralChisq) = 2.0*sqrt2*(d.ν + 3.0*d.λ)/sqrt(d.ν + 2.0*d.λ)^3
kurtosis(d::NoncentralChisq) = 12.0*(d.ν + 4.0*d.λ)/(d.ν + 2.0*d.λ)^2

function mgf(d::NoncentralChisq, t::Float64)
    exp(d.λ * t/(1.0 - 2.0 * t))*(1.0 - 2.0 * t)^(-d.ν / 2.0)
end

function cf(d::NoncentralChisq, t::Float64)
    cis(d.λ * t/(1.0 - 2.0 * im * t))*(1.0 - 2.0 * im * t)^(-d.ν / 2.0)
end

### Evaluation & Sampling

@_delegate_statsfuns NoncentralChisq nchisq ν λ

rand(d::NoncentralChisq) = StatsFuns.Rmath.nchisqrand(d.ν, d.λ)


text:
#!/usr/bin/env julia

using Downloads, JSON3, Base.BinaryPlatforms, Scratch, SHA, Pkg
using Base: UUID

# Download versions.json, start iterating over Julia versions
versions_json_url = "https://julialang-s3.julialang.org/bin/versions.json"
num_concurrent_downloads = 8

@info("Downloading versions.json...")
json_buff = IOBuffer()
Downloads.download(versions_json_url, json_buff)
versions_json = JSON3.read(String(take!(json_buff)))

# Collect all versions that are >= 1.0.0, and are a stable release
versions = filter(versions_json) do (v, d)
    if VersionNumber(string(v)) < v"1.0.0"
        return false
    end
    if !d["stable"]
        return false
    end
    return true
end

# Build download URLs for each one, and then tack on the next release as well
function select_url_hash(data, host = HostPlatform())
    d = first(filter(f -> platforms_match(f.triplet, host), data.files))
    return (d.url, d.sha256)
end
version_urls = sort(select_url_hash.(values(versions)), by = pair -> pair[1])

function generate_nightly_url(major, minor, host = HostPlatform())
    # Map arch
    arch_str = Dict("x86_64" => "x64", "i686" => "x86", "aarch64" => "aarch64", "armv7l" => "armv7l", "ppc64le" => "ppc64le")[arch(host)]
    # Map OS name
    os_str = Dict("linux" => "linux", "windows" => "winnt", "macos" => "macos", "freebsd" => "freebsd")[os(host)]
    # Map wordsize tag
    wordsize_str = Dict("x86_64" => "64", "i686" => "32", "aarch64" => "aarch64", "armv7l" => "armv7l", "ppc64le" => "ppc64")[arch(host)]

    return string(
        "https://julialangnightlies-s3.julialang.org/bin/",
        # linux/
        os_str, "/",
        # x64/
        arch_str, "/",
        # 1.6/
        string(major), ".", string(minor), "/",
        "julia-latest-",
        # linux64
        os(host), wordsize_str,
        ".tar.gz",
    )
end
highest_release = maximum(VersionNumber.(string.(keys(versions))))
next_release_url = generate_nightly_url(highest_release.major, highest_release.minor + 1)
push!(version_urls, (next_release_url, ""))
@info("Identified $(length(version_urls)) versions to try...")

# Next, we're going to download each of these to a scratch space
scratch_dir = @get_scratch!("julia_installers")

# Ensure we always download the nightly
rm(joinpath(scratch_dir, basename(last(version_urls)[1])); force=true)

# Helper function to print out stdlibs from a Julia installation
function get_stdlibs(scratch_dir, julia_installer_name)
    installer_path = joinpath(scratch_dir, julia_installer_name)
    mktempdir() do dir
        @info("Extracting $(julia_installer_name)")
        run(`tar -C $(dir) --strip-components=1 -zxf $(installer_path)`)

        jlexe = joinpath(dir, "bin", @static Sys.iswindows() ? "julia.exe" : "julia")
        jlflags = ["--startup-file=no", "-O0"]
        jlvers = VersionNumber(readchomp(`$(jlexe) $(jlflags) -e 'print(VERSION)'`))
        jlvers = VersionNumber(jlvers.major, jlvers.minor, jlvers.patch)
        @info("Auto-detected Julia version $(jlvers)")

        if jlvers < v"1.1"
            stdlibs_str = readchomp(`$(jlexe) $(jlflags) -e 'import Pkg; println(repr(Pkg.Types.gather_stdlib_uuids()))'`)
        else
            stdlibs_str = readchomp(`$(jlexe) $(jlflags) -e 'import Pkg; println(repr(Pkg.Types.load_stdlib()))'`)
        end
        return (jlvers, stdlibs_str)
    end
end

jobs = Channel()
versions_dict = Dict()

@sync begin
    # Feeder task
    @async begin
        for (url, hash) in version_urls
            put!(jobs, (url, hash))
        end
        close(jobs)
    end

    # Consumer tasks
    for _ in 1:num_concurrent_downloads
        @async begin
            for (url, hash) in jobs
                fname = joinpath(scratch_dir, basename(url))
                if !isfile(fname)
                    @info("Downloading $(url)")
                    Downloads.download(url, fname)
                end

                if !isempty(hash)
                    calc_hash = bytes2hex(open(io -> sha256(io), fname, "r"))
                    if calc_hash != hash
                        @error("Hash mismatch on $(fname)")
                    end
                end

                version, stdlibs = get_stdlibs(scratch_dir, basename(url))
                versions_dict[version] = eval(Meta.parse(stdlibs))
            end
        end
    end
end

# Next, drop versions that are the same as the one "before" them:
sorted_versions = sort(collect(keys(versions_dict)))
versions_to_drop = VersionNumber[]
for idx in 2:length(sorted_versions)
    if versions_dict[sorted_versions[idx-1]] == versions_dict[sorted_versions[idx]]
        push!(versions_to_drop, sorted_versions[idx])
    end
end
for v in versions_to_drop
    delete!(versions_dict, v)
end

# Next, figure out which stdlibs are actually unresolvable, because they've never been registered
all_stdlibs = Dict{UUID,String}()
for (julia_ver, stdlibs) in versions_dict
    merge!(all_stdlibs, stdlibs)
end

registries = Pkg.Registry.reachable_registries()
unregistered_stdlibs = filter(all_stdlibs) do (uuid, name)
    return !any(haskey(reg.pkgs, uuid) for reg in registries)
end

# Helper function for getting these printed out in a nicely-sorted order
function print_sorted(io::IO, d::Dict)
    println(io, "Dict(")
    for pair in sort(collect(d), by = kv-> kv[2])
        println(io, "        ", pair, ",")
    end
    println(io, "    ),")
end

output_fname = joinpath(dirname(dirname(@__DIR__)), "src", "HistoricalStdlibs.jl")
@info("Outputting to $(output_fname)")
sorted_versions = sort(collect(keys(versions_dict)))
open(output_fname, "w") do io
    print(io, """
    ## This file autogenerated by ext/HistoricalStdlibGenerator/generate_historical_stdlibs.jl

    using Base: UUID

    # Julia standard libraries with duplicate entries removed so as to store only the
    # first release in a set of releases that all contain the same set of stdlibs.
    const STDLIBS_BY_VERSION = [
    """)
    for v in sorted_versions
        print(io, "    $(repr(v)) => ")
        print_sorted(io, versions_dict[v])
    end
    println(io, "]")

    print(io, """
    # Next, we also embed a list of stdlibs that must _always_ be treated as stdlibs,
    # because they cannot be resolved in the registry; they have only ever existed within
    # the Julia stdlib source tree, and because of that, trying to resolve them will fail.
    const UNREGISTERED_STDLIBS = Dict(
    """)
    for (uuid, name) in unregistered_stdlibs
        println(io, "    $(repr(uuid)) => $(repr(name)),")
    end
    println(io, ")")
end


text:
using Flux

versioninfo()
include("bench_utils.jl")

@info "Benchmark Dense"
include("dense.jl")

@info "Benchmark Conv"
include("conv.jl")

@info "Benchmark VGG"
include("vgg.jl")


text:
import Base.eltypeof, Base.promote_eltypeof, LinearAlgebra.lmul!, LinearAlgebra.rmul!

export opEye, opOnes, opZeros, opDiagonal, opRestriction, opExtension, BlockDiagonalOperator

"""`opEye()`

Identity operator.
```
opI = opEye()
v = rand(5)
@assert opI * v === v
```
"""
struct opEye <: AbstractLinearOperator{Any} end

lmul!(op :: opEye, x::AbstractArray{T, 1}) where {T} = nothing
rmul!(x::AbstractArray{T, 1}, op :: opEye) where {T} = nothing
lmul!(op :: opEye, A::AbstractArray{T, 2}) where {T} = nothing
rmul!(A::AbstractArray{T, 2}, op :: opEye) where {T} = nothing
lmul!(op :: opEye, T::AbstractLinearOperator) = nothing
rmul!(T::AbstractLinearOperator, op :: opEye) = nothing

*(::opEye, x::AbstractArray{T, 1} where {T}) = x
*(x::AbstractArray{T, 1} where {T}, ::opEye) = x
*(::opEye, A::AbstractArray{T, 2} where {T}) = A
*(A::AbstractArray{T, 2} where {T}, ::opEye) = A
*(::opEye, T::AbstractLinearOperator) = T
*(T::AbstractLinearOperator, ::opEye) = T
*(::opEye, T::opEye) = T

function show(io::IO, op::opEye)
  println(io, "Identity operator")
end

function mulOpEye!(res, v, α, β::T, n_min) where T
  if β != zero(T)
    res[1:n_min] .= @views α .* v[1:n_min] .+ β .* res[1:n_min]
    res[n_min+1:end] .= β
  else
    res[1:n_min] .= @views α .* v[1:n_min]
    res[n_min+1:end] .= 0
  end
end

"""
    opEye(Mv)
    opEye(T, n)
    opEye(n)

Identity operator of order `n` and of data type `T` (defaults to `Float64`) with preallocated vector `Mv`.
"""
function opEye(Mv::AbstractVector{T}) where T
  n = length(Mv)
  prod! = @closure (res, v, α, β) -> mulOpEye!(res, v, α, β, n)
  LinearOperator{T}(n, n, true, true, prod!, prod!, prod!, Mv, Mv, Mv)
end

opEye(T::DataType, n::Int) = opEye(zeros(T, n))
opEye(n::Int) = opEye(Float64, n)

# TODO: not type stable
"""
    opEye(T, nrow, ncol)
    opEye(nrow, ncol)

Rectangular identity operator of size `nrow`x`ncol` and of data type `T`
(defaults to `Float64`).
"""
function opEye(Mv::AbstractVector{T}, nrow::I, ncol::I) where {T,I<:Integer}
  length(Mv) == nrow || throw(LinearOperatorException("shape mismatch"))
  if nrow == ncol
    return opEye(T, nrow)
  end
  prod! = @closure (res, v, α, β) -> mulOpEye!(res, v, α, β, min(nrow, ncol))
  return LinearOperator{T}(nrow, ncol, false, false, prod!, prod!, prod!, Mv, Mv, Mv)
end

opEye(T::DataType, nrow::I, ncol::I) where {I<:Integer} =  opEye(zeros(T, nrow), nrow, ncol)
opEye(nrow::I, ncol::I) where {I<:Integer} = opEye(Float64, nrow, ncol)

function mulOpOnes!(res, v, α, β::T) where T
  if β != zero(T)
    res .= (α * sum(v)) .+ β .* res
  else
    res .= (α * sum(v))
  end
end

"""
    opOnes(Mv, nrow, ncol)
    opOnes(T, nrow, ncol)
    opOnes(nrow, ncol)

Operator of all ones of size `nrow`-by-`ncol` of data type `T` (defaults to
`Float64`) and storage vector Mv.
"""
function opOnes(Mv::AbstractVector{T}, nrow::I, ncol::I) where {T,I<:Integer}
  length(Mv) == nrow || throw(LinearOperatorException("shape mismatch"))
  prod! = @closure (res, v, α, β) -> mulOpOnes!(res, v, α, β)
  LinearOperator{T}(nrow, ncol, nrow == ncol, nrow == ncol, prod!, prod!, prod!, Mv, Mv, Mv)
end

opOnes(T::DataType, nrow::I, ncol::I) where {I<:Integer} = opOnes(zeros(T, nrow), nrow, ncol)
opOnes(nrow::I, ncol::I) where {I<:Integer} = opOnes(Float64, nrow, ncol)

function mulOpZeros!(res, v, α, β::T) where T
  if β != zero(T)
    res .*= β
  end
end

"""
    opZeros(Mv, nrow, ncol)
    opZeros(T, nrow, ncol)
    opZeros(nrow, ncol)

Zero operator of size `nrow`-by-`ncol`, of data type `T` (defaults to
`Float64`) and storage vector Mv.
"""
function opZeros(Mv::AbstractVector{T}, nrow::I, ncol::I) where {T,I<:Integer}
  length(Mv) == nrow || throw(LinearOperatorException("shape mismatch"))
  prod! = @closure (res, v, α, β) -> mulOpZeros!(res, v, α, β)
  LinearOperator{T}(nrow, ncol, nrow == ncol, nrow == ncol, prod!, prod!, prod!, Mv, Mv, Mv)
end

opZeros(T::DataType, nrow::I, ncol::I) where {I<:Integer} = opZeros(zeros(T, nrow), nrow, ncol)
opZeros(nrow::I, ncol::I) where {I<:Integer} = opZeros(Float64, nrow, ncol)

function mulSquareOpDiagonal!(res, d, v, α, β::T) where T
  if β != zero(T)
    res .= α .* d .* v .+ β .* res
  else
    res .= α .* d .* v
  end
end

"""
    opDiagonal(Mv, Maw, d)
    opDiagonal(d)

Diagonal operator with the vector `d` on its main diagonal and storage vectors `Mv` and `Maw`.
"""
function opDiagonal(Mv::AbstractVector{T}, Maw::AbstractVector{T}, d::AbstractVector{T}) where {T}
  prod! = @closure (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
  ctprod! = @closure (res, w, α, β) -> mulSquareOpDiagonal!(res, conj.(d), w, α, β)
  LinearOperator{T}(length(d), length(d), true, isreal(d), prod!, prod!, ctprod!, Mv, Mv, Maw)
end

function opDiagonal(d::AbstractVector{T}) where {T}
  nrow = length(d)
  Mv = zeros(T, nrow)
  Maw = isreal(d) ? Mv : zeros(T, nrow)
  return opDiagonal(Mv, Maw, d)
end

#TODO: not type stable
function mulOpDiagonal!(res, d, v, α, β, n_min)
  res[1:n_min] .= @views α .* d[1:n_min] .* v[1:n_min] .+ β .* res[1:n_min]
  res[n_min+1:end] .= 0
end
"""
    opDiagonal(Mv, Mtu, Maw, nrow, ncol, d)
    opDiagonal(nrow, ncol, d)

Rectangular diagonal operator of size `nrow`-by-`ncol` with the vector `d` on
its main diagonal and storage vectors `Mv`, `Mtu`, `Maw`.
"""
function opDiagonal(Mv::AbstractVector{T}, Mtu::AbstractVector{T}, Maw::AbstractVector{T},
                    nrow::I, ncol::I, d::AbstractVector{T}) where {T,I<:Integer}
  nrow == ncol <= length(d) && (return opDiagonal(d[1:nrow]))
  n_min = min(nrow, ncol)
  prod! = @closure (res, v, α, β) -> mulOpDiagonal!(res, d, v, α, β, n_min)
  tprod! = @closure (res, u, α, β) -> mulOpDiagonal!(res, d, u, α, β, n_min)
  ctprod! = @closure (res, w, α, β) -> mulOpDiagonal!(res, conj.(d), w, α, β, n_min)
  LinearOperator{T}(nrow, ncol, false, false, prod!, tprod!, ctprod!, Mv, Mtu, Maw)
end

function opDiagonal(nrow::I, ncol::I, d::AbstractVector{T}) where {T,I<:Integer}
  Mv = zeros(T, nrow)
  Mtu = zeros(T, ncol)
  Maw = isreal(d) ? Mtu : zeros(T, ncol)
  return opDiagonal(Mv, Mtu, Maw, nrow, ncol, d)
end

function mulRestrict!(res, I, v, α, β)
  res .= v[I]
end

function multRestrict!(res, I, u, α, β)
  res .= 0
  res[I] = u
end

"""
    Z = opRestriction(Mv, Mtu, I, ncol)
    Z = opRestriction(I, ncol)
    Z = opRestriction(:, ncol)

Creates a LinearOperator restricting a `ncol`-sized vector to indices `I` with storage Vectors `Mv`, `Mtu`.
The operation `Z * v` is equivalent to `v[I]`. `I` can be `:`.

    Z = opRestriction(k, ncol)

Alias for `opRestriction([k], ncol)`.
"""
function opRestriction(Mv::AbstractVector{T}, Mtu::AbstractVector{T}, Idx::LinearOperatorIndexType{I}, ncol::I) where {T,I<:Integer}
  all(1 .≤ Idx .≤ ncol) || throw(LinearOperatorException("indices should be between 1 and $ncol"))
  nrow = length(Idx)
  prod! = @closure (res, v, α, β) -> mulRestrict!(res, Idx, v, α, β)
  tprod! = @closure (res, u, α, β) -> multRestrict!(res, Idx, u, α, β)
  return LinearOperator{I}(nrow, ncol, false, false, prod!, tprod!, tprod!, Mv, Mtu, Mtu)
end

opRestriction(Idx::LinearOperatorIndexType{I}, ncol::I) where {I<:Integer} = opRestriction(zeros(length(Idx)), zeros(ncol), Idx, ncol)

opRestriction(::Colon, ncol::I) where {I<:Integer} = opEye(I, ncol)

opRestriction(k::I, ncol::I) where {I<:Integer} = opRestriction([k], ncol)

"""
    Z = opExtension(I, ncol)
    Z = opExtension(:, ncol)

Creates a LinearOperator extending a vector of size `length(I)` to size `ncol`,
where the position of the elements on the new vector are given by the indices
`I`.
The operation `w = Z * v` is equivalent to `w = zeros(ncol); w[I] = v`.

    Z = opExtension(k, ncol)

Alias for `opExtension([k], ncol)`.
"""
opExtension(Idx::LinearOperatorIndexType{I}, ncol::I) where {I<:Integer} = opRestriction(Idx, ncol)'

opExtension(::Colon, ncol::I) where {I<:Integer} = opEye(I, ncol)

opExtension(k::I, ncol::I) where {I<:Integer} = opExtension([k], ncol)

# indexing for linear operators
import Base.getindex
function getindex(
  op::AbstractLinearOperator,
  rows::Union{LinearOperatorIndexType{I}, I, Colon},
  cols::Union{LinearOperatorIndexType{I}, I, Colon},
) where {I<:Integer}
  R = opRestriction(rows, size(op, 1))
  E = opExtension(cols, size(op, 2))
  return R * op * E
end

eltypeof(op::AbstractLinearOperator) = eltype(op)  # need this for promote_eltypeof

"""
    BlockDiagonalOperator(Mv, Mtu, Maw, M1, M2, ..., Mn)
    BlockDiagonalOperator(M1, M2, ..., Mn)

Creates a block-diagonal linear operator:

    [ M1           ]
    [    M2        ]
    [       ...    ]
    [           Mn ]
"""
function BlockDiagonalOperator(Mv::AbstractVector{S}, Mtu::AbstractVector{S}, Maw::AbstractVector{S}, ops...) where {S}
  nrow = ncol = 0
  for op ∈ ops
    m, n = size(op)
    nrow += m
    ncol += n
  end
  T = promote_eltypeof(ops...)

  function prod!(y, x, α, β)
    # y = zeros(T, nrow)
    k = 0
    j = 0
    for op ∈ ops
      m, n = size(op)
      # y[(k + 1):(k + m)] .= op * x[(j + 1):(j + n)]
      @views mul!(y[(k + 1):(k + m)] , op, x[(j + 1):(j + n)], α, β)
      k += m
      j += n
    end
    # y
  end

  function tprod!(y, x, α, β)
    # y = zeros(T, ncol)
    k = 0
    j = 0
    for op ∈ ops
      m, n = size(op)
      # y[(k + 1):(k + n)] .= transpose(op) * x[(j + 1):(j + m)]
      @views mul!(y[(k + 1):(k + n)], transpose(op), x[(j + 1):(j + m)], α, β)
      k += n
      j += m
    end
    # y
  end

  function ctprod!(y, x, α, β)
    # y = zeros(T, ncol)
    k = 0
    j = 0
    for op ∈ ops
      m, n = size(op)
      # y[(k + 1):(k + n)] .= op' * x[(j + 1):(j + m)]
      @views mul!(y[(k + 1):(k + n)], adjoint(op), x[(j + 1):(j + m)], α, β)
      k += n
      j += m
    end
    # y
  end

  symm = all((issymmetric(op) for op ∈ ops))
  herm = all((ishermitian(op) for op ∈ ops))
  LinearOperator{T}(nrow, ncol, symm, herm, prod!, tprod!, ctprod!, Mv, Mtu, Maw)
end

function BlockDiagonalOperator(ops...)
  nrow = ncol = 0
  for op ∈ ops
    m, n = size(op)
    nrow += m
    ncol += n
  end
  T = promote_eltypeof(ops...)
  symm = all((issymmetric(op) for op ∈ ops))
  herm = all((ishermitian(op) for op ∈ ops))
  Mv = zeros(T, nrow)
  Mtu = symm ? Mv : zeros(T, ncol)
  Maw = herm ? Mv : zeros(T, ncol)
  return BlockDiagonalOperator(Mv, Mtu, Maw, ops...)
end


text:
using MLJTSVDInterface # substitute for correct interface pkg name
using Test
using TSVD
using MLJBase
using SparseArrays
using StableRNGs # for RNGs stable across all julia versions

const rng = StableRNGs.StableRNG(123)

@testset "tsvd transformer" begin
    n = 10
    p = 20
    prob_nonzero = 0.5

    # test with a sparse matrix
    X_sparse = sprand(rng, n, p, prob_nonzero)

    # use defaults - transform into an n x 2 dense matrix
    model = TSVDTransformer(rng=42)

    mach = machine(model, X_sparse)
    fit!(mach, verbosity=0)
    X_transformed = transform(mach, X_sparse)

    # also do the raw transformation with TSVD library
    U, s, V = tsvd(X_sparse, 2)

    @test size(X_transformed) == (10, 2)
    @test isapprox(s, fitted_params(mach).singular_values)
    @test size(V) == size(fitted_params(mach).components)

    # test with a dense matrix
    X_dense = rand(rng, n, p)

    mach = machine(model, X_dense)
    fit!(mach, verbosity=0)
    X_transformed = transform(mach, X_dense)

    # also do the raw transformation with TSVD library
    U, s, V = tsvd(X_dense, 2)

    @test size(X_transformed) == (10, 2)
    @test isapprox(s, fitted_params(mach).singular_values)
    @test size(V) == size(fitted_params(mach).components)

    # test tables
    X, _ = make_regression(100, 5)
    mach = machine(model, X)
    fit!(mach, verbosity=0)
    X_transformed = transform(mach, X)

    @test length(keys(X_transformed)) == 2

    # test with default RNG
    model = TSVDTransformer()

    mach = machine(model, X_sparse)
    fit!(mach, verbosity=0)
    X_transformed = transform(mach, X_sparse)

    # also do the raw transformation with TSVD library
    U, s, V = tsvd(X_sparse, 2)

    @test size(X_transformed) == (10, 2)
    @test isapprox(s, fitted_params(mach).singular_values)
    @test size(V) == size(fitted_params(mach).components)
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis
using ConstrainedControl

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1

p2 = [0.0;0.0;length1 / 2] # joint connection point

# Initial orientation
ϕ1 = 0
q1 = UnitQuaternion(RotX(ϕ1))
ϕ2 = 0
q2 = UnitQuaternion(RotX(ϕ2))

# Links
origin = Origin{Float64}()
link1 = Box(width, depth, length1, length1, color = RGBA(1., 1., 0.))
link2 = deepcopy(link1)

# Constraints
joint1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=p2))
joint2 = EqualityConstraint(Revolute(link1, link2, joint_axis; p1=-p2, p2=p2))

links = [link1;link2]
constraints = [joint1;joint2]

mech = Mechanism(origin, links, constraints)
setPosition!(origin,link1,p2 = p2,Δq = q1)
setPosition!(link1,link2,p1=-p2,p2 = p2,Δq = q2)

pid = PID(mech, getfield.(constraints,:id), [pi/2;-pi/4], P = [10.;10.], I = [10.;10], D = [5.;5.])

storage = simulate!(mech,10.,pid,record = true)
visualize(mech,storage)