text:
# Bond guessing algorithm derived from from Chemfiles, itself from VMD

# Data from Blue Obelisk's data repository, corrected for H, C, O, N, S and F
# (see https://github.com/chemfiles/chemfiles/issues/301#issuecomment-574100048)
const vdwradii = Float32[1.0, 1.4, 2.2, 1.9, 1.8, 1.5, 1.4, 1.3, 1.2, 1.54, 2.4,
                         2.2, 2.1, 2.1, 1.95, 1.9, 1.8, 1.88, 2.8, 2.4, 2.3,
                         2.15, 2.05, 2.05, 2.05, 2.05, 2.0, 2.0, 2.0, 2.1, 2.1,
                         2.1, 2.05, 1.9, 1.9, 2.02, 2.9, 2.55, 2.4, 2.3, 2.15,
                         2.1, 2.05, 2.05, 2.0, 2.05, 2.1, 2.2, 2.2, 2.25, 2.2,
                         2.1, 2.1, 2.16, 3.0, 2.7, 2.5, 2.48, 2.47, 2.45, 2.43,
                         2.42, 2.4, 2.38, 2.37, 2.35, 2.33, 2.32, 2.3, 2.28,
                         2.27, 2.25, 2.2, 2.1, 2.05, 2.0, 2.0, 2.05, 2.1, 2.05,
                         2.2, 2.3, 2.3, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.4, 2.0,
                         2.3, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
                         2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
                         2.0, 2.0, 2.0, 2.0, 2.0, 2.0]

# Data from PeriodicTable.jl
const ismetal = Bool[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1,
                     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1,
                     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,
                     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                     1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1]

const ismetalloid = Bool[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

"""
    guess_bonds(pos, types, mat, options)

Return the bonds guessed from the positions, types and cell matrix, given as a
`Vector{Vector{Tuple{Int,Float32}}}`.

The `i`-th entry of the list is a list, whose entries are of the form `(j, dist)` which
indicates that the representatives of vertices `i` and `j` distant of at most `dist` are
bonded together.
"""
function guess_bonds(pos, types, mat, options)
    # Algorithm from chemfiles, itself from VMD
    @ifwarn if options.bonding != Bonding.Guess
        @warn "Guessing bonds with custom algorithm (from Chemfiles and VMD). This may take a while for big structures and may be inexact."
        @info "To avoid guessing bonds, use a file format that contains the bonds."
    end
    n = length(pos)
    bonds = [Tuple{Int,Float32}[] for _ in 1:n]
    @toggleassert n == length(types)
    radii = Vector{Float32}(undef, n)
    for (i, typ) in enumerate(types)
        t = get(atomic_numbers, typ, nothing)
        if t isa Int
            radii[i] = vdwradii[t]*(1 + options.wider_metallic_bonds*(ismetal[t]|ismetalloid[t])*0.5)
        else
            radii[i] = 0.0
            @ifwarn @warn lazy"Unrecognized atom type \"$t\" will be considered a dummy atom."
        end
    end
    cutoff = 3*(options.cutoff_coeff^3.1) * max(maximum(radii), 0.833)
    cutoff2 = 13*options.cutoff_coeff/15
    buffer, ortho, safemin = prepare_periodic_distance_computations(mat)
    for i in 1:n
        radius_i = radii[i]
        iszero(radius_i) && continue
        posi = pos[i]
        typi = types[i]
        #ignoreifmetallic = typi === :C
        #ignoreifC = ismetal[atomic_numbers[typi]]
        skiphomoatomic = typi ∈ options.ignore_homoatomic_bonds ||
                         (options.ignore_homometallic_bonds && ismetal[atomic_numbers[typi]])
        acceptonlyO = options.structure === StructureType.Zeolite && typi !== :O
        acceptallbutO = options.structure === StructureType.Zeolite && typi === :O
        for j in (i+1):n
            typj = types[j]
            skiphomoatomic && typj === typi && continue
            acceptonlyO && typj !== :O && continue
            acceptallbutO && typj === :O && continue
            #(ignoreifC & (typj === :C)) && continue
            #ignoreifmetallic && ismetal[atomic_numbers[typj]] && continue
            radius_j = radii[j]
            iszero(radius_j) && continue
            posj = pos[j]
            buffer .= posi .- posj
            d1 = periodic_distance!(buffer, mat, ortho, safemin)
            maxdist = cutoff2*(radius_i + radius_j)
            if d1 < cutoff && 0.5 < d1 < maxdist
                push!(bonds[i], (j, maxdist))
                push!(bonds[j], (i, maxdist))
            end
        end
    end
    return bonds
end

ext: jl


text:
using DataFrames, XLSX, Statistics
using Query
using Gadfly, Cairo, Fontconfig

table = XLSX.readtable(
    "./results/SSITmethods_2minmax_16core_mcgonagall.xlsx",
    "method_problem_results")
df = DataFrame(table...)

theme() = Theme(
    panel_stroke=colorant"black",
    background_color=colorant"white",
    # point_size=.07cm,
    grid_color=colorant"grey")

function save_plot(p, name)
    img = PNG("analysis/graphs/$name.png", 30cm, 35cm, dpi=300)
    draw(img, p)
end

feasible_df = df |> @filter(_.infeasibility == 0 && _.dataset > 4) |> DataFrame
infeasible_df = df |> @filter(_.infeasibility > 0 && _.dataset > 4) |> DataFrame

names(feasible_df)
summary_stats = combine(
    groupby(feasible_df, :method_name),
    :lowest_gap=>mean,
    :total_time=>mean)

xticks = [.005, .05, .25, .5, .75]
p = plot(feasible_df, x="lowest_gap", y="total_time",
    ygroup="case", xgroup="dataset",
    color="method_name",
    # Scale.x_log10(),
    Geom.subplot_grid(Geom.point),
    theme())
save_plot(p, "infeasible_scatter")

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
# Somewhat scrappy reference MNIST implementation.

using WGANGP
using Flux
using CUDA
using Base.Iterators: partition
using Plots, Images
using Printf
using Dates
using ProgressMeter
using TensorBoardLogger
using Logging
using BSON: @save, @load
using Parameters: @with_kw
using MLDatasets

@with_kw struct Hyperparameters
    batch_size::Int = 16
    critic_train_factor::Int = 5
    epochs::Int = 100
    verbose_freq::Int = 20
    sample_freq::Int = 1000
    gpu::Bool = true
end

function create_montage(genr, z, n, m; skip_batchnorm = true, kws...)
    skip_batchnorm && @eval Flux.istraining() = false
    out_tensor = cpu(genr(z))
    skip_batchnorm && @eval Flux.istraining() = true
    img_tensor = Gray.(clamp.(out_tensor, 0.0f0, 1.0f0))
    imgs = [img_tensor[:, :, 1, i] for i in 1:n*m]
    return Plots.plot(
        mosaicview(imgs; nrow = n, fillvalue = 1, npad = 5);
        title = "Generated images", size = (600, 600),
        xaxis = false, yaxis = false, axis = nothing,
        kws...
    )
end

# TODO/WARN: Actually doesn't mutate genr and crit if they live on CPU when
# this function is called, so make sure to catch them when they come home from the gym.
function train!(genr, crit, opt_genr, opt_crit, data_tensor, hps = Hyperparameters(); run_prefix = nothing, runid = nothing)

    if hps.gpu == true
        device = gpu
    else
        error("CPU training not implemented")
    end

    genr = genr |> device
    crit = crit |> device

    # Infer latent_dim from generator network architecture
    latent_dim = size(genr.layers[1].weight, 2)

    # Create some folders
    runid = isnothing(runid) ? "mnist" : runid
    path = mkpath("output/$runid")
    imspath = mkpath("$path/ims")

    # Tensorboard logger
    tblog = TBLogger("logs/$runid")

    # For output images
    fixed_noise = randn(Float32, latent_dim, 25) |> device

    train_steps = 0
    data = [data_tensor[:, :, :, r] |> device for
            r in partition(1:size(data_tensor, 4), hps.batch_size)]
    p = Progress(hps.epochs * length(data), "Training...")

    with_logger(tblog) do
        for epoch in 1:hps.epochs
            for i = 1:length(data)

                z = CUDA.randn(latent_dim, hps.batch_size)
                x_generated = genr(z)
                x_true = data[i]

                loss_crit = step_critic!(opt_crit, crit, x_true, x_generated, hps.batch_size)

                if i % hps.critic_train_factor == 1 || hps.critic_train_factor == 1
                    z = CUDA.randn(latent_dim, hps.batch_size)
                    loss_genr = step_generator!(opt_genr, genr, crit, z)
                end

                if train_steps % hps.verbose_freq == 0
                    @info "Losses" critic = loss_crit generator = loss_genr
                    ProgressMeter.update!(p, train_steps;
                        showvalues = [
                            ("Epoch", epoch),
                            ("Step", train_steps),
                            ("Critic loss", loss_crit),
                            ("Generator loss", loss_genr)]
                    )
                end
                if train_steps % hps.sample_freq == 0
                    # Save generated images
                    static_output_montage = create_montage(genr, fixed_noise, 5, 5; size = (800, 800))
                    png(static_output_montage, @sprintf("%s/static_step_%06d.png", imspath, train_steps))
                    sample_output_montage = create_montage(genr, CUDA.randn(latent_dim, 25), 5, 5; size = (800, 800))
                    png(sample_output_montage, @sprintf("%s/sample_step_%06d.png", imspath, train_steps))
                    @info "Generated" static = static_output_montage sampled = sample_output_montage
                end
                train_steps += 1
            end
        end
    end

    return genr, crit
end

function create_mnist_models(latent_dim = 100)
    genr = Chain(
        Dense(latent_dim, 7 * 7 * 256, relu),
        x -> reshape(x, 7, 7, 256, :),
        ConvTranspose((5, 5), 256 => 128; stride = 1, pad = 2),
        BatchNorm(128, relu),
        ConvTranspose((4, 4), 128 => 64; stride = 2, pad = 1),
        BatchNorm(64, relu),
        ConvTranspose((4, 4), 64 => 1, sigmoid; stride = 2, pad = 1)
    ) |> gpu

    crit = Chain(
        Conv((5, 5), 1 => 64, selu; stride = 2, pad = 1),
        Conv((5, 5), 64 => 128, selu; stride = 2, pad = 1),
        Conv((5, 5), 128 => 256, selu; stride = 2, pad = 1),
        x -> reshape(x, 1024, :),
        Dense(1024, 1)
    ) |> gpu

    return (genr, crit)
end

function prepare_mnist_data()
    ims = MNIST.traintensor(Float32)
    return Flux.unsqueeze(ims, 3)
end

genr, crit = create_mnist_models()
opt_genr = ADAM(1e-4, (0.9, 0.99))
opt_crit = ADAM(1e-4, (0.9, 0.99))
data = prepare_mnist_data()
hps = Hyperparameters()
train!(genr, crit, opt_genr, opt_crit, data, hps)


text:
export Map_From_VectFunc

@doc raw"""

An [`Abstract_Map`](@ref) instance that can be easily created by
providing a Julia function

# Constructor
```julia
Map_From_VectFunc(n_θ,f)
```
where ``(X,θ)->Y=f(X,θ)`` is the wrapped function.

# Example

```jldoctest
using NLS_Fit

f(X,θ) = θ[1]*X

map = Map_From_VectFunc(1,f)

X = Float64[1,2,3]

θ = Float64[4]

eval_map(map,X,θ)

# output
3-element Vector{Float64}:
  4.0
  8.0
 12.0

```

"""
struct Map_From_VectFunc{F<:Function} <: Abstract_Map
    _n_θ::Int
    _f::F
end

parameter_size(m::Map_From_VectFunc) = m._n_θ

function eval_map(m::Map_From_VectFunc,
                  X::AbstractVector,
                  θ::AbstractVector)
    @assert length(θ) == parameter_size(m)

    m._f(X,θ)
end


text:
using Test

using MultivariatePolynomials
const MP = MultivariatePolynomials

using SemialgebraicSets

struct DummyAlgo <: SemialgebraicSets.AbstractGröbnerBasisAlgorithm end
function SemialgebraicSets.gröbnerbasis!(::AbstractVector{<:MP.APL}, ::DummyAlgo)
    error("Dummy algo")
end

struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
function SemialgebraicSets.solvealgebraicequations(
    ::SemialgebraicSets.AlgebraicSet,
    ::DummySolver,
)
    error("Dummy solver")
end

function algo_and_solver_test(func)
    Mod.@polyvar x y
    p = [y^2 - x * y, x^2 + x * y]
    V = func(p, DummyAlgo())
    @test_throws ErrorException("Dummy algo") p[1] in V
    @test_throws ErrorException("Dummy algo") rem(p[1], V.I)
    V = func(p, DummySolver())
    @test_throws ErrorException("Dummy solver") collect(V)
    V = func(p, DummyAlgo(), DummySolver())
    @test_throws ErrorException("Dummy solver") collect(V)
    @test_throws ErrorException("Dummy algo") rem(p[1], V.I)
end

@testset "Algo and solver" begin
    algo_and_solver_test(algebraicset)
    algo_and_solver_test(projectivealgebraicset)
end


text:
"""
Code for finding a line or wavelength range in spectra.

Author: Eric Ford
Created: August 2020
Contact: https://github.com/eford/
"""

""" Return list of all order indices that contain a pixel with wavelength lambda """
function find_orders_with_line(goal::Real,lambda::AbstractArray{T,2}) where T<:Real
   order_min(i) = lambda[1,i]
   order_max(i) = lambda[end,i]
#=
    for i in 1:10
       println("# i= ",i," order_min= ",order_min(i)," order_max= ",order_max(i), "   goal= ",goal)
   end
   flush(stdout)
=#
   findall(i->order_min(i)<=goal<=order_max(i), 1:size(lambda,2) )
end

""" Return list of all order indices that contain all pixels with wavelengths between goal_lo and goal_hi """
function find_orders_with_line(goal_lo::Real,goal_hi::Real,lambda::AbstractArray{T,2}) where T<:Real
   order_min(i) = lambda[1,i]
   order_max(i) = lambda[end,i]
   findall(i->order_min(i)<=goal_lo && goal_hi<=order_max(i), 1:size(lambda,2) )
end

""" Return list of all order indices that include any wavelengths between goal_lo and goal_hi """
function find_orders_in_range(goal_lo::Real,goal_hi::Real,lambda::AbstractArray{T,2}) where T<:Real
   order_min(i) = lambda[1,i]
   order_max(i) = lambda[end,i]
   findall(i-> (goal_lo<=order_min(i)<=goal_hi) || (goal_lo<=order_max(i)<=goal_hi), 1:size(lambda,2) )
end

# Find indicies for pixels around lines
const Δλoλ_fit_line_default = 5*(1.8*1000/speed_of_light_mps)
const Δλoλ_edge_pad_default = 0*(1.8*1000/speed_of_light_mps)

""" Return a range of columns indices with wavelengths within Δ of line_center """
function find_cols_to_fit(wavelengths::AbstractArray{T,1}, line_center::Real; Δ::Real = Δλoλ_fit_line_default) where T<:Real
    @assert Δ >= zero(Δ)
    first = findfirst(x->x>=line_center*(1-Δ),wavelengths)
    last = findlast(x->x<=line_center*(1+Δ),wavelengths)
    if isnothing(first) || isnothing(last)   return 0:0   end
    if last<first    return last:first    end
    return first:last
end

""" Return a range of columns indices with wavelengths between line_lo and line_hi """
function find_cols_to_fit(wavelengths::AbstractArray{T,1}, line_lo::Real, line_hi::Real; Δ::Real = Δλoλ_edge_pad_default) where T<:Real
    @assert line_lo < line_hi
    first = findfirst(x->x>=line_lo*(1-Δ),wavelengths)
    last = findlast(x->x<=line_hi*(1+Δ),wavelengths)
    if isnothing(first) || isnothing(last)   return 0:0   end
    if last<first    return last:first    end
    return first:last
end

""" Return list of (pixels, order_idx) pairs that contain pixels with desireed wavelengths.
    Excludes locations that contain any pixels with var == NaN.
"""
function findall_line end

function findall_line(goal::Real,lambda::AbstractArray{T1,2},var::AbstractArray{T2,2}; Δ::Real = Δλoλ_fit_line_default) where {T1<:Real, T2<:Real}
    @assert lambda[1,1] <= goal <= lambda[end,end]
    @assert size(lambda) == size(var)
    @assert Δ >= zero(Δ)
    orders = find_orders_with_line(goal,lambda)
    @assert length(orders) >= 1
    locs = map(o->(pixels=find_cols_to_fit(lambda[:,o],goal,Δ=Δ),order=o), orders)
    locs_good_idx = findall(t->!any(isnan.(var[t[1],t[2]])),locs)
    #locs_good_idx = findall(t-> !(first(t.pixels)==0 || last(t.pixels)==0 || t.order==0) && (!any(isnan.(var[t.pixels,t.order]))) ,locs)
    if length(locs) != length(locs_good_idx)
        locs = locs[locs_good_idx]
    end
    return locs
end

function findall_line(goal_lo::Real,goal_hi::Real, lambda::AbstractArray{T1,2},var::AbstractArray{T2,2}; Δ::Real = Δλoλ_edge_pad_default, verbose::Bool = false) where {T1<:Real, T2<:Real}
    @assert lambda[1,1] <= goal_lo < goal_hi <= lambda[end,end]
    orders = find_orders_with_line(goal_lo,goal_hi,lambda)
    #if ! (length(orders) >= 1) return end
    #=
    if verbose
        for i in 1:5
        println("# i= ",i," min(order)= ",minimum(lambda[:,i])," max(order)= ",maximum(lambda[:,i]), "   goal_lo= ",goal_lo, " goal_hi = ",goal_hi)
        end
    end
    flush(stdout)
    =#
    @assert length(orders) >= 1
    locs = map(o->(pixels=find_cols_to_fit(lambda[:,o],goal_lo, goal_hi,Δ=Δ),order=o), orders)
    #locs_good_idx = findall(t->!any(isnan.(var[t[1],t[2]])),locs)
    locs_good_idx = findall(t-> !(first(t.pixels)==0 || last(t.pixels)==0 || t.order==0) && (!any(isnan.(var[t.pixels,t.order]))) ,locs)
    if length(locs) != length(locs_good_idx)
        locs = locs[locs_good_idx]
    end
    return locs
end

function findall_line(goal::Real,lambda::AbstractArray{T1,1},var::AbstractArray{T2,1}; Δ::Real = Δλoλ_fit_line_default) where {T1<:Real, T2<:Real}
    @assert lambda[1] <= goal <= lambda[end]
    @assert size(lambda) == size(var)
    @assert Δ >= zero(Δ)
    locs = find_cols_to_fit(lambda,goal,Δ=Δ)
    locs_good_idx = findall(t->!any(isnan.(var[t])),locs)
    #locs_good_idx = findall(t-> !(first(t.pixels)==0 || last(t.pixels)==0 || t.order==0) && (!any(isnan.(var[t.pixels,t.order]))) ,locs)
    if length(locs) != length(locs_good_idx)
        locs = locs[locs_good_idx]
    end
    return locs
end

function findall_line(goal_lo::Real,goal_hi::Real, lambda::AbstractArray{T1,1},var::AbstractArray{T2,1}; Δ::Real = Δλoλ_edge_pad_default, verbose::Bool = false) where {T1<:Real, T2<:Real}
    @assert lambda[1] <= goal_lo < goal_hi <= lambda[end]
#=    if verbose
        for i in 1:5
        println("# i= ",i," min(order)= ",minimum(lambda[:,i])," max(order)= ",maximum(lambda[:,i]), "   goal_lo= ",goal_lo, " goal_hi = ",goal_hi)
        end
    end
    flush(stdout)
    =#
    locs = find_cols_to_fit(lambda,goal_lo, goal_hi,Δ=Δ)
    locs_good_idx = findall(t-> !any(isnan.(var[t])) ,locs)
    #locs_good_idx = findall(t-> !(first(t.pixels)==0 || last(t.pixels)==0 || t.order==0) && (!any(isnan.(var[t.pixels,t.order]))) ,locs)
    if length(locs) != length(locs_good_idx)
        locs = locs[locs_good_idx]
    end
    return locs
end

function findall_line(goal::Real,spectra::AS; Δ::Real = Δλoλ_fit_line_default) where {AS<:AbstractSpectra}
    findall_line(goal,spectra.λ,spectra.var, Δ=Δ)
end

function findall_line(goal_lo::Real,goal_hi::Real,spectra::AS; Δ::Real = Δλoλ_edge_pad_default) where {AS<:AbstractSpectra}
    findall_line(goal_lo,goal_hi,spectra.λ,spectra.var, Δ=Δ)
end

""" Return (pixels, order_idx) pair that contain "best" region of spectra, based on highest SNR. """
function find_line_best end

function find_line_best(goal::Real,lambda::AbstractArray{T1,2},flux::AbstractArray{T2,2},var::AbstractArray{T3,2}; Δ::Real = Δλoλ_fit_line_default) where {T1<:Real, T2<:Real, T3<:Real}
    locs = findall_line(goal,lambda,var,Δ=Δ)
    if length(locs) == 0   return  missing end
    #scores = map( t->sum( flux[t[1],t[2]] ./ var[t[1],t[2]])/sum( 1.0 ./ var[t[1],t[2]]), locs)
    scores = map( t->calc_snr(flux[t[1],t[2]],var[t[1],t[2]]), locs)
    idx_best = findmax(scores)
    locs[idx_best[2]]
end

function find_line_best(goal_lo::Real,goal_hi::Real, lambda::AbstractArray{T1,2},flux::AbstractArray{T2,2},var::AbstractArray{T3,2}; Δ::Real = Δλoλ_edge_pad_default) where {T1<:Real, T2<:Real, T3<:Real}
    locs = findall_line(goal_lo,goal_hi,lambda,var,Δ=Δ)
    if length(locs) == 0
        println("=>(",goal_lo, ", ",goal_hi, ")  Δ=",Δ)
        return  missing
    end
    #scores = map( t->sum( flux[t[1],t[2]] ./ var[t[1],t[2]])/sum( 1.0 ./ var[t[1],t[2]]), locs)
    scores = map( t->calc_snr(flux[t[1],t[2]],var[t[1],t[2]]), locs)
    idx_best = findmax(scores)
    locs[idx_best[2]]
end

function find_line_best(goal::Real,lambda::AbstractArray{T1,1},flux::AbstractArray{T2,1},var::AbstractArray{T3,1}; Δ::Real = Δλoλ_fit_line_default) where {T1<:Real, T2<:Real, T3<:Real}
    cols = find_cols_to_fit(lambda,goal, Δ=Δ)
    @assert( ( first(cols)==0 && last(cols)==0)  || !any(isnan.(var[cols])) )
    return cols
    #=
    locs = findall_line(goal,lambda,var,Δ=Δ)
    if length(locs) == 0   return  missing end
    #scores = map( t->sum( flux[t[1],t[2]] ./ var[t[1],t[2]])/sum( 1.0 ./ var[t[1],t[2]]), locs)
    return locs
    scores = map( t->calc_snr(flux[t[1],t[2]],var[t[1],t[2]]), locs)
    idx_best = findmax(scores)
    locs[idx_best[2]]
    =#
end

function find_line_best(goal_lo::Real,goal_hi::Real, lambda::AbstractArray{T1,1},flux::AbstractArray{T2,1},var::AbstractArray{T3,1}; Δ::Real = Δλoλ_edge_pad_default) where {T1<:Real, T2<:Real, T3<:Real}
    cols = find_cols_to_fit(lambda,goal_lo, goal_hi, Δ=Δ)
    @assert( ( first(cols)==0 && last(cols)==0)  || !any(isnan.(var[cols])) )
    return cols
    #=
    locs = findall_line(goal_lo,goal_hi,lambda,var,Δ=Δ)
    if length(locs) == 0
        println("=>(",goal_lo, ", ",goal_hi, ")  Δ=",Δ)
        return  missing
    end
    return locs
    #scores = map( t->sum( flux[t[1],t[2]] ./ var[t[1],t[2]])/sum( 1.0 ./ var[t[1],t[2]]), locs)
    scores = map( t->calc_snr(flux[t],var[t]), locs)
    idx_best = findmax(scores)
    locs[idx_best[2]]
    =#
end

function find_line_best(goal::Real,spectra::AS; Δ::Real = Δλoλ_fit_line_default) where {AS<:AbstractSpectra}
    find_line_best(goal,spectra.λ,spectra.flux,spectra.var, Δ=Δ)
end

function find_line_best(goal_lo::Real,goal_hi::Real,spectra::AS; Δ::Real = Δλoλ_edge_pad_default) where {AS<:AbstractSpectra}
    find_line_best(goal_lo,goal_hi,spectra.λ,spectra.flux,spectra.var, Δ=Δ)
end

""" Find pixels included in a range of wavelengths """
function find_pixels_for_line_in_chunk( chunk::AbstractChunkOfSpectrum, λ_min::Real, λ_max::Real )# ; plan::LineFinderPlan = LineFinderPlan() )
  idx_lo = searchsortedfirst(chunk.λ, λ_min, by=x->x>=λ_min)
  idx_tmp = searchsortedlast(chunk.λ[idx_lo:end], λ_max, by=x->x<=λ_max, rev=true)
  idx_hi = idx_lo + idx_tmp - 1
  return idx_lo:idx_hi
end

function find_pixels_for_line_in_chunklist( chunk_list::AbstractChunkList, λ_min::Real, λ_max::Real; verbose::Bool = true )
  ch_idx_all = findall(c-> (λ_min <= minimum(chunk_list.data[c].λ)) && (maximum(chunk_list.data[c].λ) <= λ_max) ,1:length(chunk_list))
  println("Hello")
  #map(c->(chunk_idx=c, pixels=find_pixels_for_line_in_chunk(chunk_list.data[c], λ_min, λ_max) ), ch_idx)
  ch_idx = 0
  if length(ch_idx_all) > 1
    snr_of_chunks_with_line = map(c->RvSpectMLBase.calc_snr(chunk_list.data[c].flux, chunk_list.data[c].var), ch_idx_all)
    ch_idx_to_keep = argmax(snr_of_chunks_with_line)
    ch_idx = ch_idx_all[ch_idx_to_keep]
    if verbose
      println(" Found λ=",λ_min,"-",λ_max," in chunks: ", ch_idx_all, " containing ", length.(ch_idx_all), " pixels.")
      println(" SNRs = ", snr_of_chunks_with_line)
      println(" Keeping chunk #",ch_idx)
    end
  elseif length(ch_idx_all) == 1
    ch_idx = first(ch_idx_all)
    if verbose
      println(" Found λ=",λ_min,"-",λ_max," in chunk: ", ch_idx, " containing ", length(ch_idx), " pixels.")
      snr_of_chunk_with_line = RvSpectMLBase.calc_snr(chunk_list.data[ch_idx].flux, chunk_list.data[ch_idx].var)
      println(" SNRs = ", snr_of_chunk_with_line)
    end

  end
  if ch_idx == 0
    error("Didn't find λ = " *string(λ_min)*" - " *string(λ_max)* " in chunklist.")

  end
  return (chunk_idx=ch_idx, pixels=find_pixels_for_line_in_chunk(chunk_list.data[ch_idx], λ_min, λ_max) )
end

function find_pixels_for_line_in_chunklist( chunk_list::AbstractChunkList, λ_min::Real, λ_max::Real, chunk_id::Integer)
  return (chunk_idx=chunk_id, pixels=find_pixels_for_line_in_chunk(chunk_list.data[chunk_id], λ_min, λ_max) )
end

""" `is_in_wavelength_range_list_any_order(λ; list )`
Return true if λ is between lambda_lo and lambda_hi for any row in list
"""
function is_in_wavelength_range_list_any_order(λ::Real; list::DataFrame  )
    @assert hasproperty(list, :lambda_lo)
    @assert hasproperty(list, :lambda_hi)
    idx =  searchsortedfirst(list[:,:lambda_hi], λ)
    return idx>size(list,1) || !(list[idx,:lambda_lo]<=λ<=list[idx,:lambda_hi]) ?  false : true
end

""" `is_in_wavelength_range_list(λ; order, list )`
Return true if λ is between lambda_lo and lambda_hi for any row in list
"""
function is_in_wavelength_range_list(λ::Real; order::Integer = 0, list::DataFrame  )
    @assert hasproperty(list, :lambda_lo)
    @assert hasproperty(list, :lambda_hi)
    if order == 0
        return is_in_wavelength_range_list_any_order(λ, list=list)
    else
        @assert hasproperty(list, :order)
        list_for_order = list |> @filter(_.order==order) |> DataFrame
        @assert issorted(list_for_order.lambda_hi)
        idx =  searchsortedfirst(list_for_order[:,:lambda_hi], λ)
        return idx>size(list_for_order,1) || !(list_for_order[idx,:lambda_lo]<=λ<=list_for_order[idx,:lambda_hi]) ?  false : true
    end
end

""" `is_in_wavelength_range_list(λ_lo, λ_hi; list )`
Return true if there is overlap between (λ_lo, λ_hi) and lambda_lo and lambda_hi for any row in list
# TODO: test
"""
function is_in_wavelength_range_list(λ_lo::Real, λ_hi::Real; list::DataFrame  )
    @assert λ_lo < λ_hi
    @assert hasproperty(list, :lambda_lo)
    @assert hasproperty(list, :lambda_hi)
    idx =  searchsortedfirst(list[:,:lambda_hi], λ_lo)
    if idx>size(list,1)    return false  end
    if λ_lo<=list[idx,:lambda_hi] &&  λ_hi>=list[idx,:lambda_lo]
        return true
    else
        return false
    end
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
#=
This is where I create denoising functions for the plug and play methods
=#

using BSON: @load
using MLDatasets: FashionMNIST
using NNlib
using Flux
using Zygote
using LinearAlgebra

#=
Here I copy Babhru's VAE denoiser
=#

# This function determines the architecture
function create_vae()
    # Define the encoder and decoder networks
    encoder_features = Chain(
        Dense(784,500, relu),
        Dense(500,500, relu)
    )
    encoder_μ = Chain(encoder_features, Dense(500, 20))
    encoder_logvar = Chain(encoder_features, Dense(500, 20))

    decoder = Chain(
        Dense(20, 500, relu, bias = false),
        Dense(500, 500, relu, bias = false),
        Dense(500, 784, sigmoid, bias = false)
    )
    return encoder_μ, encoder_logvar, decoder
end

# This function loads the parameters to the model
# For now I'm using a pre-trained model

function load_model(load_dir::String, epoch::Int)
    print("Loading model...")
    @load joinpath(load_dir, "model-$epoch.bson") encoder_μ encoder_logvar decoder
    println("Done")
    return encoder_μ, encoder_logvar, decoder
end

function reconstruct_images(encoder_μ, encoder_logvar, decoder, x)
    # Forward propagate through mean encoder and std encoders
    μ = encoder_μ(x)
    logvar = encoder_logvar(x)
    # Apply reparameterisation trick to sample latent
    z = μ + randn(Float32, size(logvar)) .* exp.(0.5f0 * logvar)
    # Reconstruct from latent sample
    x̂ = decoder(z)
    return clamp.(x̂, 0 ,1)
end

function decoder_loss_function(z_variable, decoder, y_noisy)
    return 0.5*norm(y_noisy - decoder(z_variable))^2
end

# TO DO: rewrite as iterator
function gradient_descent(loss_function, z0, decoder, y_noisy, maxit, stepsize_numerator)
    zk = z0
    for t in 0:maxit
        loss_k = loss_function(zk, decoder, y_noisy)
        #println("at iteration $t loss is $loss_k \n")
        #resize_and_save_single_MNIST_image(clamp.(decoder(zk), 0.0, 1.0), "./scripts/temp/run_denoising", "during_GD_$t")
        zk .= zk - (stepsize_numerator/sqrt(t+1))*Zygote.gradient(z -> loss_function(z, decoder, y_noisy), zk)[1]
    end
    return zk
end

struct gradient_descent_state
    z_curr
    loss_value
    gradient
end

mutable struct denoising_problem
    z # optimization variable
    b # fit D(z) to observed image vector b
    loss_function # z -> 1/2||b - D(z)||_{2}^2
    loss_value # 1/2||b - D(z_curr)||_{2}^2
    decoder # D
end

function decoder_denoiser!(x,y,encoder_μ, encoder_logvar, decoder,num_iter)
    #TO DO: load model outside of denoiser otherwise it's slow
    #model_epoch_number = 20
    #model_dir = "./src/utilities/saved_models/MNIST"
    #encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)
    #resize_and_save_single_MNIST_image(y_curr,"./scripts/temp/denoising_decoder","before_denoising")
    loss_function = decoder_loss_function

    # as an initial point for gradient descent, take the current iterate (y) and apply encoder
    μ = encoder_μ(y)
    logvar = encoder_logvar(y)
    z0 = μ + randn(Float32, size(logvar)) .* exp.(0.5f0 * logvar)

   # num_iter = 20 #number of iterations of gradient descent

    z_out = gradient_descent(loss_function,z0,decoder,y,num_iter,0.05)

    x .= decoder(z_out)
end


text:
module testExamples
using LinearAlgebra
using Test

srcDir = dirname(dirname(@__FILE__))*"\\src\\Examples\\"
include(srcDir *"Examples_sep.jl")
include(srcDir *"Examples_ent.jl")
include(srcDir *"Examples.jl")

using .Examples_sep
using .Examples_ent
using .Examples

@testset "Test the call examples" begin
    ρ_ex = Examples.get_example()
    @test ρ_ex["RANDℂb5"] == Examples.get_example("RANDℂb5")
    @test ρ_ex["RANDℂb5"] == Examples.get_example(5)
end

@testset "Test the call example meta data" begin
    ρ_ex = Examples.get_example()
    df = Examples.get_example_meta()
    @test df[df.Example .== "Eq14Hor97a",:] ==  Examples.get_example_meta("Eq14Hor97a")
    @test df[df.Example .== "Eq14Hor97a",:].Example[1] ==  Examples.get_example_meta(1).Example
    length(ρ_ex) == size(df)[1]
end

##
df = ex.get_example_meta()
entdf = df[df.isSeparable .== "ent",:]
ρ_dict = ex.get_example()

for n_ex in 1:23
    exa = entdf.Example[n_ex]
    ρ = ρ_dict[exa]
    d = entdf.Size[n_ex]
    println("$(entdf.Example[n_ex]) $(la.eigmin(ρ))")
    @assert us.isPPᵀ(ρ,d) == entdf.hasPPᵀ[n_ex]
end

##
df = ex.get_example_meta()
sepdf = df[df.isSeparable .== "sep",:]
ρ_dict = ex.get_example()

for n_ex in 1:36
    exa = sepdf.Example[n_ex]
    ρ = ρ_dict[exa]
    d = sepdf.Size[n_ex]
    println("$(sepdf.Example[n_ex]) $(la.eigvals(ρ)[1])")
    @assert us.isPPᵀ(ρ,d) == sepdf.hasPPᵀ[n_ex]
end

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
## Julia 1.1.0

# Initialize
using Plots
plotlyjs() # Backend PlotlyJS

angles = 8 * 2π

a = 4
b = 8

steps = 100
anim = @animate for i = 1:steps
    Θ = range(0, stop=i * angles / (steps), length=10*i)
    x = a .* cos.(Θ)
    y = a .* sin.(Θ)
    z = b .* Θ
    plot(x, y, z, size=(960, 960))
end

gif(anim, "./gif/circular-helix.gif", fps = 30)

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
include("kepler.jl")

#function kepler_init(time::Float64,mass::Float64,elements::Array{Float64,1})
function kepler_init(time::T,mass::T,elements::Array{T,1}) where {T <: Real}
# Takes orbital elements of a single Keplerian; returns positions & velocities.
# This is 3D), so 6 orbital elements specified, the code returns 3D.  For
# Inclination = pi/2, motion is in X-Z plane; sky plane is X-Y.
# Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
period = elements[1]
# Compute the semi-major axis in AU (or other units specified by GNEWT):
semi = cbrt(GNEWT*mass*period^2/4/pi^2)
# Convert to eccentricity & longitude of periastron:
ecc2=elements[3]^2+elements[4]^2
ecc=sqrt(ecc2)
#omega = atan2(elements[4],elements[3])
omega = atan(elements[4],elements[3])
# The true anomaly at the time of transit:
f1 = 1.5*pi-omega
# Compute the time of periastron passage:
sqrt1mecc2 = sqrt(1.0-ecc^2)
tp=(elements[2]+period*sqrt1mecc2/2.0/pi*(ecc*sin(f1)/(1.0+ecc*cos(f1))
#    -2.0/sqrt1mecc2*atan2(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
    -2.0/sqrt1mecc2*atan(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
# Compute the mean anomaly
n = 2pi/period
m=n*(time-tp)
# Kepler solver:
if ecc > 0.0
  f = kepler(m,ecc)
else
  f = m
end
ecosfp1 = 1.0+ecc*cos(f)
fdot = n*ecosfp1^2/sqrt1mecc2^3
# Compute the radial distance:
r=semi*(1.0-ecc^2)/ecosfp1
rdot = semi*n/sqrt1mecc2*ecc*sin(f)
# For now assume plane-parallel:
#inc = pi/2
#capomega = pi
inc = elements[5]
capomega = elements[6]
# Now, compute the positions
x = zeros(eltype(elements),3)
v = zeros(eltype(elements),3)
if abs(capomega-pi) > 1e-15
  coscapomega = cos(capomega) ; sincapomega = sin(capomega)
else
  coscapomega = -1.0 ; sincapomega = 0.0
end
cosomegapf = cos(omega+f) ; sinomegapf = sin(omega+f)
if abs(inc-pi/2) > 1e-15
  cosinc = cos(inc) ; sininc = sin(inc)
else
  cosinc = 0.0 ; sininc = 1.0
end
x[1]=r*(coscapomega*cosomegapf-sincapomega*sinomegapf*cosinc)
x[2]=r*(sincapomega*cosomegapf+coscapomega*sinomegapf*cosinc)
x[3]=r*sinomegapf*sininc
rdotonr = rdot/r
# Compute the velocities:
rfdot = r*fdot
v[1]=x[1]*rdotonr+rfdot*(-coscapomega*sinomegapf-sincapomega*cosomegapf*cosinc)
v[2]=x[2]*rdotonr+rfdot*(-sincapomega*sinomegapf+coscapomega*cosomegapf*cosinc)
v[3]=x[3]*rdotonr+rfdot*cosomegapf*sininc
return x,v
end

function kepler_init(time::T,mass::T,elements::Array{T,1},jac_init::Array{T,2}) where {T <: Real}
# Takes orbital elements of a single Keplerian; returns positions & velocities.
# This is 3D), so 6 orbital elements specified, the code returns 3D.  For
# Inclination = pi/2, motion is in X-Z plane; sky plane is X-Y.
# Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
# Returns the position & velocity of Keplerian.
# jac_init is the derivative of (x,v,m) with respect to (elements,m).
period = elements[1]
n = 2pi/period
t0 = elements[2]
# Compute the semi-major axis in AU (or other units specified by GNEWT):
semi = cbrt(GNEWT*mass*period^2/4/pi^2)
dsemidp = 2third*semi/period
dsemidm = third*semi/mass
# Convert to eccentricity & longitude of periastron:
ecosomega = elements[3]
esinomega = elements[4]
ecc=sqrt(esinomega^2+ecosomega^2)
deccdecos = ecosomega/ecc
deccdesin = esinomega/ecc
#omega = atan2(esinomega,ecosomega)
# The true anomaly at the time of transit:
#f1 = 1.5*pi-omega
# Compute the time of periastron passage:
sqrt1mecc2 = sqrt(1.0-ecc^2)
#tp=(t0+period*sqrt1mecc2/2.0/pi*(ecc*sin(f1)/(1.0+ecc*cos(f1))
#    -2.0/sqrt1mecc2*atan2(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
den1 = esinomega-ecosomega-ecc
tp = (t0 - sqrt1mecc2/n*ecosomega/(1.0-esinomega)-2/n*
#     atan2(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1))
     atan(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1))
dtpdp = (tp-t0)/period
fac = sqrt((1.0-ecc)/(1.0+ecc))
den2 = 1.0/den1^2
theta = fac*(esinomega+ecosomega+ecc)/den1
#println("theta: ",theta," tan(atan2()): ",tan(atan2(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1)))
dthetadecc = ((ecc+ecosomega)^2+2*(1.0-ecc^2)*esinomega-esinomega^2)/(sqrt1mecc2*(1.0+ecc))*den2
dthetadecos = 2fac*esinomega*den2
dthetadesin = -2fac*(ecosomega+ecc)*den2
dtpdecc = ecc/sqrt1mecc2/n*ecosomega/(1.0-esinomega)-2/n/(1.0+theta^2)*dthetadecc
dtpdecos = dtpdecc*deccdecos -sqrt1mecc2/n/(1.0-esinomega)-2/n/(1.0+theta^2)*dthetadecos
dtpdesin = dtpdecc*deccdesin -sqrt1mecc2/n*ecosomega/(1.0-esinomega)^2-2/n/(1.0+theta^2)*dthetadesin
dtpdt0 = 1.0
# Compute the mean anomaly
m=n*(time-tp)
dmdp = -m/period
dmdtp = -n
# Kepler solver: instead of true anomaly, return eccentric anomaly:
ekep=ekepler(m,ecc)
cosekep = cos(ekep); sinekep = sin(ekep)
# Compute the radial distance:
r=semi*(1.0-ecc*cosekep)
drdekep = semi*ecc*sinekep
drdecc = -semi*cosekep
denom = semi/r
dekepdecos = sinekep*denom*deccdecos
dekepdesin = sinekep*denom*deccdesin
dekepdm   = denom
inc = elements[5]
capomega = elements[6]
# Now, compute the positions
coscapomega = cos(capomega) ; sincapomega = sin(capomega)
cosomega = ecosomega/ecc; sinomega = esinomega/ecc
cosinc = cos(inc) ; sininc = sin(inc)
# Define rotation matrices (M&D 2.119-2.120):
P1 = [cosomega -sinomega 0.0; sinomega cosomega 0.0; 0.0 0.0 1.0]
P2 = [1.0 0.0 0.0; 0.0 cosinc -sininc; 0.0 sininc cosinc]
P3 = [coscapomega -sincapomega 0.0; sincapomega coscapomega 0.0; 0.0 0.0 1.0]
P321 = P3*P2*P1
# Express x & v in terms of eccentric anomaly (2.121, 2.41, 2.68):
xplane = semi*[cosekep-ecc;  sqrt1mecc2*sinekep;  0.0] # position vector in orbital plane
vplane = [-sinekep; sqrt1mecc2*cosekep; 0.0]
x = P321*xplane
# Now derivatives:
dxda = x/semi
dxdekep = P321*semi*vplane
P32 = P3*P2
#dxdecc  = -x/ecc + P321*semi*[-1.0; -ecc/sqrt1mecc2*sinekep; 0.0]
# Get rid of cancellations in the prior expression:
dxdecc  = -P321*semi/ecc*[cosekep; sinekep/sqrt1mecc2; zero(T)]
# These may need to be rewritten to flag ecc = 0.0 case:
dxdecos = dxdecc*deccdecos + P32/ecc*xplane
dxdesin = dxdecc*deccdesin + P32/ecc*[-xplane[2]; xplane[1]; 0.0]
dxdinc  = P3*[0.0 0.0 0.0; 0.0 -sininc -cosinc; 0.0 cosinc -sininc]*P1*xplane
dxdcom  = [-sincapomega -coscapomega 0.0; coscapomega -sincapomega 0.0; 0.0 0.0 0.0]*P2*P1*xplane
# Compute the velocities:
v = P321*n*semi*denom*vplane
#decc = 1e-5
#dvp = P321*ecc/(ecc+decc)*n*semi/(1.0-(ecc+decc)*cosekep)*[-sinekep; sqrt(1.0-(ecc+decc)^2)*cosekep; 0.0]
#dvm = P321*ecc/(ecc-decc)*n*semi/(1.0-(ecc-decc)*cosekep)*[-sinekep; sqrt(1.0-(ecc-decc)^2)*cosekep; 0.0]
dvda    = v/semi
dvdp    = -v/period
dvdekep = -v*ecc*sinekep*denom+P321*n*semi*denom*[-cosekep; -sqrt1mecc2*sinekep; 0.0]
dvdecc  = -v/ecc + v*cosekep*denom + P321*n*semi*denom*[0.0; -ecc/sqrt1mecc2*cosekep; 0.0]
#dvdecc_num = .5*(dvp-dvm)/decc
#println("dvdecc: ",dvdecc," dvdecc_num: ",dvdecc_num)
#dvdecc = dvdecc_num
dvdecos = dvdecc*deccdecos + P32*n*semi*denom/ecc*vplane
dvdesin = dvdecc*deccdesin + P32*n*semi*denom/ecc*[-vplane[2]; vplane[1]; 0.0]
dvdinc  = P3*[0.0 0.0 0.0; 0.0 -sininc -cosinc; 0.0 cosinc -sininc]*P1*n*semi*denom*vplane
dvdcom  = [-sincapomega -coscapomega 0.0; coscapomega -sincapomega 0.0; 0.0 0.0 0.0]*P2*P1*n*semi*denom*vplane
# Now, take derivatives (11/15/2017 notes):
# Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
fill!(jac_init,0.0)
jac_init[1:3,1] = dxda*dsemidp + dxdekep*dekepdm*(dmdp+dmdtp*dtpdp)
#if  T != BigFloat
#  println("position vs. period: term 1 ",dxda*dsemidp," term 2: ",dxdekep*dekepdm*dmdp," term 3: ",dxdekep*dekepdm*dmdtp*dtpdp," sum: ",jac_init[1:3,1])
#end
jac_init[1:3,2] = dxdekep*dekepdm*dmdtp*dtpdt0
jac_init[1:3,3] = dxdecos + dxdekep*(dekepdm*dmdtp*dtpdecos + dekepdecos)
#if  T != BigFloat
#  println("position vs. ecosom : term 1 ",dxdecos," term 2: ",dxdekep*dekepdm*dmdtp*dtpdecos," term 3: ",dxdekep*dekepdm*dekepdecos," sum: ",jac_init[1:3,3])
#end
jac_init[1:3,4] = dxdesin + dxdekep*(dekepdm*dmdtp*dtpdesin + dekepdesin)
#if  T != BigFloat
#  println("position vs. esinom : term 1 ",dxdesin," term 2: ",dxdekep*dekepdm*dmdtp*dtpdesin," term 3: ",dxdekep*dekepdm*dekepdesin," sum: ",jac_init[1:3,4])
#end
jac_init[1:3,5] = dxdinc
jac_init[1:3,6] = dxdcom
jac_init[1:3,7] = dxda*dsemidm
jac_init[4:6,1] = dvdp + dvda*dsemidp +  dvdekep*dekepdm*(dmdp+dmdtp*dtpdp)
#if  T != BigFloat
#  println("velocity vs. period: term 1 ",dvdp," term 2: ",dvda*dsemidp," term 3: ",dvdekep*dekepdm*dmdp," term 4: ",dvdekep*dekepdm*dmdtp*dtpdp," sum: ",jac_init[4:6,1])
#end
jac_init[4:6,2] = dvdekep*dekepdm*dmdtp*dtpdt0
jac_init[4:6,3] = dvdecos + dvdekep*(dekepdm*dmdtp*dtpdecos + dekepdecos)
#if  T != BigFloat
#  println("velocity vs. ecosom : term 1 ",dvdecos," term 2: ",dvdekep*dekepdm*dmdtp*dtpdecos," term 3: ",dvdekep*dekepdm*dekepdecos," sum: ",jac_init[4:6,3])
#end
jac_init[4:6,4] = dvdesin + dvdekep*(dekepdm*dmdtp*dtpdesin + dekepdesin)
#if  T != BigFloat
#  println("velocity vs. esinom : term 1 ",dvdesin," term 2: ",dvdekep*dekepdm*dmdtp*dtpdesin," term 3: ",dvdekep*dekepdm*dekepdesin," sum: ",jac_init[4:6,4])
#end
jac_init[4:6,5] = dvdinc
jac_init[4:6,6] = dvdcom
jac_init[4:6,7] = dvda*dsemidm
jac_init[7,7] = 1.0
#return x,v,tp,dtpdecos,dtpdesin
#return x,v,ekep,dekepdm*dmdtp*dtpdecos + dekepdecos,dekepdm*dmdtp*dtpdesin + dekepdesin
return x,v
end


text:
using Test, StateSpaceModels, Statistics, CSV, LinearAlgebra, Random, Distributions

const SSM = StateSpaceModels

# Run tests
include("test_utils.jl")
include("test_userdefined.jl")
include("test_locallevel.jl")
include("test_lineartrend.jl")
include("test_structural.jl")
include("test_regression.jl")
include("test_simulate.jl")
include("test_opt_methods.jl")
include("test_results.jl")

max_stars_repo_name: PaulMainwood/StateSpaceCPG

max_issues_repo_name: PaulMainwood/StateSpaceCPG

max_forks_repo_name: PaulMainwood/StateSpaceCPG


text:
export State
export is_density_matrix, is_pure, is_mixed

"""
    is_density_matrix( ρ :: Matrix; atol=ATOL :: Float64 ) :: Bool

Returns `true` if `ρ` is a valid density matrix. The following constraints must be
satisfied for all density matrices:
* Hermitian: `ρ' == ρ`
* Positive Semidefinite: `eigmin(ρ) ≥ 0`
* Trace-one: `tr(ρ) == 1`
"""
function is_density_matrix(ρ::Matrix; atol=ATOL :: Float64)::Bool
    if !is_hermitian(ρ, atol=atol)
        return false
    elseif !isapprox(tr(ρ), 1, atol=atol)
        return false
    elseif !is_positive_semidefinite(ρ, atol=atol)
        return false
    end

    return true
end

"""
    _density_matrix_error(ρ; atol=ATOL :: Float64)

Throws a `DomainError` indicating why density matrix `ρ` is not a density matrix
the error messages are self-explanatory and are listed as follows:
* `"Density matrix `ρ` is not hermitian."`
* `"Density matrix `ρ` is not trace-one."`
* `"Density matrix `ρ` is not positive semi-definite."`
"""
function _density_matrix_error(ρ; atol=ATOL :: Float64)
    if !is_hermitian(ρ, atol=atol)
        throw(DomainError(ρ, "Density matrix `ρ` is not hermitian."))
    elseif !isapprox(tr(ρ), 1, atol=atol)
        throw(DomainError(ρ, "Density matrix `ρ` is not trace-one."))
    elseif !is_positive_semidefinite(ρ, atol=atol)
        throw(DomainError(ρ, "Density matrix `ρ` is not positive semi-definite."))
    end
end

"""
    State( ρ :: Matrix{T<:Number} ) <: Operator{T}

The density matrix representation of a quantum state. The constructor, `State(ρ)`
throws a `DomainError` if [`is_density_matrix`](@ref) evaluates to `false`.
"""
struct State{T<:Number} <: Operator{T}
    M :: Matrix{T}
    atol :: Float64
    State(
        ρ :: Matrix{<:Number}; atol=ATOL :: Float64
    ) = is_density_matrix(ρ, atol=atol) ? new{eltype(ρ)}(ρ, atol) : _density_matrix_error(ρ, atol=atol)
    State(
        ρ :: Adjoint{T,Matrix{T}}; atol=ATOL :: Float64
    ) where T <: Number = is_density_matrix(ρ[:,:], atol=atol) ? new{eltype(ρ)}(ρ[:,:], atol) : _density_matrix_error(ρ[:,:], atol=atol)
end
is_density_matrix(::State) = true

"""
    kron(states :: Vararg{State}; atol=ATOL :: Float64) :: State

Performs the kronecker product on the supplied `states` to construct a new `State`: ``\\rho_A \\otimes \\rho_B = \\rho_{AB}.``

This method extends [`Base.kron`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.kron).
"""
kron(states :: Vararg{State}; atol=ATOL :: Float64) :: State = State(
    kron(map(ρ -> ρ.M, states)...),
atol=atol)

"""
    partial_trace(ρ::State, system::Vector{Int64}, id::Int64) :: State

Takes the partial_trace of a `State` `ρ` to construct a new `State`: ``\\text{Tr}_B[\\rho_{AB}] = \\rho_A``.
"""
partial_trace(ρ::State, system::Vector{Int64}, id::Int64; atol=ATOL :: Float64) :: State = begin
    State(partial_trace(ρ.M, system, id), atol=atol)
end

"""
    eigvals(ρ :: State) :: Vector{Float64}

Computes the eigenvalues of `State` `ρ`. The eigenvalues of a density matrix
are real. A `DomainError` is thrown if the imaginary component of any eigenvalues
exceeds the absolute tolerance `ρ.atol`. This method extends [`LinearAlgebra.eigvals`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigvals).
"""
function eigvals(ρ :: State) :: Vector{Float64}
    λs = eigvals(ρ.M)
    for λ in λs
        if !isapprox(imag(λ), 0, atol=ρ.atol)
            throw(DomainError(λ), "eigenvalue λ has a imaginary component outside the absolute tolerance.")
        end
    end

    real.(λs)
end

"""
    is_pure(ρ :: State) :: Bool

Returns `true` if `ρ` is pure, i.e. `rank(ρ) == 1` within the absolute tolerance `ρ.atol`.
This method can also be used on `Matrix` types.

```julia
is_pure(ρ :: Matrix; atol=ATOL :: Float64) :: Bool
```
"""
function is_pure(ρ :: State) :: Bool
    rank(ρ) == 1
end
is_pure(ρ :: Matrix; atol=ATOL :: Float64) :: Bool = is_pure(State(ρ, atol=atol))

"""
    is_mixed(ρ :: State) :: Bool

Returns `true` if `ρ` is mixed, i.e. `rank(ρ) > 1` within the absolute tolerance `ρ.atol`.
This method also can be used on `Matrix` types.

```julia
is_mixed(ρ :: Matrix; atol=ATOL :: Float64) :: Bool
```
"""
function is_mixed(ρ :: State) :: Bool
    rank(ρ) > 1
end
is_mixed(ρ :: Matrix; atol=ATOL :: Float64) :: Bool = is_mixed(State(ρ, atol=atol))


text:
using BenchmarkTools, Plots, SNOPT7, Ipopt

T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-3,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3,
    opts_al=opts_al,
    opts_pn=opts_pn)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol)

# iLQR
prob_ilqr = copy(Problems.pendulum)
@time p0, s0 = solve(prob_ilqr, opts_ilqr)

# ALTRO w/ Newton
prob_altro = copy(Problems.pendulum)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
max_violation_direct(p1)
plot(p1.X,title="Pendulum state (ALTRO)")
plot(p1.U,title="Pendulum control (ALTRO)")

# DIRCOL w/ Ipopt
prob_ipopt = copy(Problems.pendulum)
rollout!(prob_ipopt)
prob_ipopt = update_problem(prob_ipopt,model=Dynamics.pendulum) # get continuous time model
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
max_violation_direct(p2)
plot(p2.X,title="Pendulum state (Ipopt)")
plot(p2.U,title="Pendulum control (Ipopt)")

# DIRCOL w/ SNOPT
prob_snopt = copy(Problems.pendulum)
rollout!(prob_snopt)
prob_snopt = update_problem(prob_snopt,model=Dynamics.pendulum) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
max_violation_direct(p3)
plot(p3.X,title="Pendulum state (SNOPT)")
plot(p3.U,title="Pendulum control (SNOPT)")

## Minimum Time
dt_max = 0.15
dt_min = 1.0e-3
dt = 0.15/2
R_minimum_time = 15.0

opts_altro_mt = ALTROSolverOptions{T}(verbose=verbose,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3,
    opts_al=opts_al,
    opts_pn=opts_pn,
    R_minimum_time=R_minimum_time,
    dt_max=dt_max,
    dt_min=dt_min)

opts_ipopt_mt = DIRCOLSolverMTOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    opts=Dict(:print_level=>3,
        :tol=>max_con_viol,
        :constr_viol_tol=>max_con_viol),
    R_min_time=R_minimum_time,
    h_max=dt_max,
    h_min=dt_min)

opts_snopt_mt = DIRCOLSolverMTOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    opts=Dict(:Major_print_level=>0,
        :Minor_print_level=>0,
        :Major_optimality_tolerance=>max_con_viol,
        :Major_feasibility_tolerance=>max_con_viol,
        :Minor_feasibility_tolerance=>max_con_viol),
    R_min_time=R_minimum_time,
    h_max=dt_max,
    h_min=dt_min)

# ALTRO w/ Newton
prob_altro_mt = update_problem(copy(Problems.pendulum),dt=dt,tf=0.)
@time p4, s4 = solve(prob_altro_mt, opts_altro_mt)
@benchmark p4, s4 = solve($prob_altro_mt, $opts_altro_mt)
max_violation_direct(p4)
total_time(p4)
plot(p4.X,title="Pendulum state (Min. Time) (ALTRO)")
plot(p4.U,title="Pendulum control (Min. Time) (ALTRO)")

# DIRCOL w/ Ipopt
prob_ipopt_mt = copy(Problems.pendulum)
rollout!(prob_ipopt_mt)
prob_ipopt_mt = update_problem(prob_ipopt_mt,model=Dynamics.pendulum,dt=dt,tf=0.) # get continuous time model
@time p5, s5 = solve(prob_ipopt_mt, opts_ipopt_mt)
@benchmark p5, s5 = solve($prob_ipopt_mt, $opts_ipopt_mt)
max_violation_direct(p5)
total_time(p5)
U5 = to_array([p5.U[k][1:p5.model.m] for k = 1:p5.N])
plot(p5.X,title="Pendulum state (Min. Time) (Ipopt)")
plot(U5',title="Pendulum control (Min. Time) (Ipopt)")

# DIRCOL w/ SNOPT
prob_snopt_mt = copy(Problems.pendulum)
rollout!(prob_snopt_mt)
prob_snopt_mt = update_problem(prob_snopt_mt,model=Dynamics.pendulum,dt=dt,tf=0.) # get continuous time model
@time p6, s6 = solve(prob_snopt_mt, opts_snopt_mt)
@benchmark p6, s6 = solve($prob_snopt_mt, $opts_snopt_mt)
max_violation_direct(p6)
total_time(p6)
U6 = to_array([p6.U[k][1:p6.model.m] for k = 1:p6.N])
plot(p6.X,title="Pendulum state (Min. Time) (SNOPT)")
plot(U6',title="Pendulum control (Min. Time) (SNOPT)")

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using EyeOfRa
using Test
using Revise: entr, watched_files
using Plots

@testset "All tests" begin
    function LongTest(n)
        @show rand(n)
    end

    # The below comments are what I would like to have for a test with a timeout...
    # however that version ends up blocking and defeats the point. I don't
    # understand what is going on.

    # @testset "Redirect tests" begin
    #     # Test for blocking - need to implement a dodgy timeout here.
    #     ch = Channel()
    #     @async begin
    #         EyeOfRa.TestFunction(LongTest, (), ())
    #         put!(ch, :success)
    #     end
    #     @async begin
    #         sleep(5)
    #         put!(ch, :failure)
    #     end
    #     @test take!(ch) == :success
    # end
    @testset "Redirect tests" begin
        @test_nowarn EyeOfRa.TestFunction(LongTest, (2000,), ())
        @test_nowarn EyeOfRa.TestFunction(LongTest, (4000,), ())
    end

    function BasicSuccess(x)
        println("one")
        println(stderr, "two")
        x * 2
    end
    function BasicFailure(x)
        error("Failed")
    end
    function BasicInference()
        if rand(Bool)
            1
        else
            1.0
        end
    end

    @testset "TestFunction" begin
        result = EyeOfRa.TestFunction(BasicSuccess, (2,), ())
        @test result.result == 4
        @test result.stdout == "one\n"
        @test result.stderr == "two\n"

        result = EyeOfRa.TestFunction(BasicFailure, (2,), ())
        @test result.result isa Exception

        result = EyeOfRa.TestFunction(BasicInference, (), ())
        @test result.result == 1
        @test result.inferred_type == Union{Float64,Int}
    end

    @testset "Revise interaction" begin
        # This is for an independent check on Revise upgrades, so that they don't
        # break compatibility. It is not a direct check of this package's code.
        mktempdir() do dir
            filename = joinpath(dir, "EyeOfRaTesting.jl")
            write(filename, """
module EyeOfRaTesting
end
""")

            pushfirst!(LOAD_PATH, dir)
            @eval using EyeOfRaTesting
            popfirst!(LOAD_PATH)

            ch = Channel{Symbol}(2)

            task = @async entr([], [EyeOfRaTesting], pause=0) do
                if !isdefined(EyeOfRaTesting, :testfunc)
                    push!(ch, :before)
                else
                    @assert EyeOfRaTesting.testfunc(2) == 2
                    push!(ch, :after)
                    throw(InterruptException())
                end
            end

            @test take!(ch) == :before

            write(filename, """
module EyeOfRaTesting
    testfunc(x) = x
end
""")
            # Wait a couple of seconds
            for i = 1:5
                sleep(0.2)
                istaskdone(task) && break
            end
            @test istaskdone(task)
            # InterruptExcpetion is not eaten by entr in the inner function.
            # Need to check this from the CompositeException object.
            istaskfailed(task) && task.exception.exceptions[].task.exception != InterruptException() && throw(TaskFailedException(task))
            if istaskdone(task) && !istaskfailed(task)
                @test take!(ch) == :after
            end
        end
    end

    function CaptureStdout(func)
        mktemp() do path,io
            redirect_stdout(func, io)
            close(io)
            read(path, String)
        end
    end

    done = false
    @testset "Emacs images" begin
        pushdisplay(EyeOfRa.EMACS_DISPLAY())

        try
            output = CaptureStdout() do
                display(EyeOfRa.CLEAR())
            end
            @test occursin("<emacs-clear></emacs-clear>", output)

            output = CaptureStdout() do
                display(plot(1:10))
            end
            @test occursin(r"<emacs-svg>.*</emacs-svg>", output)
        finally
            popdisplay(EyeOfRa.EMACS_DISPLAY())
        end
    end
end


text:
# comparison times for the zeta function
using SpecialFunctions
m = 100000
b = 1000.0
b = 8.0
b = 1.0

x = 2.0*b .* rand(m) .- b
y = 2.0*b .* rand(m) .- b
z = x + im*y
Z = Array{Complex{Float64}}(undef,m)

zeta(1.0)
val, t, bytes, gctime, memallocs = @timed begin
    for i=1:m
        Z[i] = zeta(z[i])
    end
end
@printf("  finished %d calculations in time t=%.2f seconds, i.e., %.1f microsec each\n", m, t, (t/m)*1.0e6)

# Mathematica timing results
# z in [-1,1], [-8,8], [-16,16]
time_s = [16.06, 17.9, 18.9]
time_micros = 1.0e6 .* (time_s ./ 10000) #


text:
function glGenOne(glGenFn)
	id = GLuint[0]
	glGenFn(1, id)
	glCheckError("generating a buffer, array, or texture")
	id[]
end

glGenBuffer() = glGenOne(glGenBuffers)

glGenVertexArray() = glGenOne(glGenVertexArrays)

glGenTexture() = glGenOne(glGenTextures)

function getInfoLog(obj::GLuint)
	# Return the info log for obj, whether it be a shader or a program.
	isShader = glIsShader(obj)
	getiv = isShader == GL_TRUE ? glGetShaderiv : glGetProgramiv
	getInfo = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog
	# Get the maximum possible length for the descriptive error message
	len = GLint[0]
	getiv(obj, GL_INFO_LOG_LENGTH, len)
	maxlength = len[]
	# TODO: Create a macro that turns the following into the above:
	# maxlength = @glPointer getiv(obj, GL_INFO_LOG_LENGTH, GLint)
	# Return the text of the message if there is any
	if maxlength > 0
		buffer = zeros(GLchar, maxlength)
		sizei = GLsizei[0]
		getInfo(obj, maxlength, sizei, buffer)
		len = sizei[]
		unsafe_string(pointer(buffer), len)
	else
		""
	end
end

function validateShader(shader)
	success = GLint[0]
	glGetShaderiv(shader, GL_COMPILE_STATUS, success)
	success[] == GL_TRUE
end

function glErrorMessage()
# Return a string representing the current OpenGL error flag, or the empty string if there's no error.
	err = glGetError()
	err == GL_NO_ERROR ? "" :
	err == GL_INVALID_ENUM ? "GL_INVALID_ENUM: An unacceptable value is specified for an enumerated argument. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_VALUE ? "GL_INVALID_VALUE: A numeric argument is out of range. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_OPERATION ? "GL_INVALID_OPERATION: The specified operation is not allowed in the current state. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_FRAMEBUFFER_OPERATION ? "GL_INVALID_FRAMEBUFFER_OPERATION: The framebuffer object is not complete. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_OUT_OF_MEMORY ? "GL_OUT_OF_MEMORY: There is not enough memory left to execute the command. The state of the GL is undefined, except for the state of the error flags, after this error is recorded." : "Unknown OpenGL error with error code $err."
end

function glCheckError(actionName="")
	message = glErrorMessage()
	if length(message) > 0
		if length(actionName) > 0
			error("Error ", actionName, ": ", message)
		else
			error("Error: ", message)
		end
	end
end

function createShader(source, typ)
# Create the shader
	shader = glCreateShader(typ)::GLuint
	if shader == 0
		error("Error creating shader: ", glErrorMessage())
	end
	# Compile the shader
	glShaderSource(shader, 1, convert(Ptr{UInt8}, pointer([convert(Ptr{GLchar}, pointer(source))])), C_NULL)
	glCompileShader(shader)
	# Check for errors
	!validateShader(shader) && error("Shader creation error: ", getInfoLog(shader))
	shader
end

function createShaderProgram(f, vertexShader, fragmentShader)
	# Create, link then return a shader program for the given shaders.
	# Create the shader program
	prog = glCreateProgram()
	if prog == 0
		error("Error creating shader program: ", glErrorMessage())
	end
	# Attach the vertex shader
	glAttachShader(prog, vertexShader)
	glCheckError("attaching vertex shader")
	# Attach the fragment shader
	glAttachShader(prog, fragmentShader)
	glCheckError("attaching fragment shader")
	f(prog)
	# Finally, link the program and check for errors.
	glLinkProgram(prog)
	status = GLint[0]
	glGetProgramiv(prog, GL_LINK_STATUS, status)
	if status[] == GL_FALSE then
		glDeleteProgram(prog)
		error("Error linking shader: ", glGetInfoLog(prog))
	end

	glDeleteShader(vertexShader)
	glDeleteShader(fragmentShader)

	prog
end

createShaderProgram(vertexShader, fragmentShader) = createShaderProgram(prog->0, vertexShader, fragmentShader)

global GLSL_VERSION = ""

function createcontextinfo()
	global GLSL_VERSION
	glsl = split(unsafe_string(glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
	if length(glsl) >= 2
		glsl = VersionNumber(parse(Int, glsl[1]), parse(Int, glsl[2]))
		GLSL_VERSION = string(glsl.major) * rpad(string(glsl.minor),2,"0")
	else
		error("Unexpected version number string. Please report this bug! GLSL version string: $(glsl)")
	end

	glv = split(unsafe_string(glGetString(GL_VERSION)), ['.', ' '])
	if length(glv) >= 2
		glv = VersionNumber(parse(Int, glv[1]), parse(Int, glv[2]))
	else
		error("Unexpected version number string. Please report this bug! OpenGL version string: $(glv)")
	end
	dict = Dict{Symbol,Any}(
	    :glsl_version   => glsl,
	    :gl_version     => glv,
	    :gl_vendor	    => unsafe_string(glGetString(GL_VENDOR)),
	    :gl_renderer	=> unsafe_string(glGetString(GL_RENDERER)),
	    #:gl_extensions => split(unsafe_string(glGetString(GL_EXTENSIONS))),
	)
end

function get_glsl_version_string()
	if isempty(GLSL_VERSION)
		error("couldn't get GLSL version, GLUTils not initialized, or context not created?")
	end
	return "#version $(GLSL_VERSION)\n"
end

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
@testset "LabeledValue" begin
    lbls1 = Dict(1=>"a", 2=>"b")
    lbls2 = Dict{Union{Float64, Missing}, String}(1.0=>"b")
    v1 = LabeledValue(1, lbls1)
    v2 = LabeledValue(2, lbls1)
    v3 = LabeledValue(1.0, lbls2)
    v4 = LabeledValue(missing, lbls2)

    @test v1 == v3
    @test isequal(v4 == v1, missing)
    @test isequal(v1, v3)
    @test isequal(v4, v4)
    @test v1 < v2
    @test isless(v1, v2)
    @test isapprox(v1, v3)

    @test v1 == 1
    @test 1 == v1
    @test v1 == "a"
    @test "a" == v1
    @test ismissing(v1 == missing)
    @test ismissing(missing == v1)
    @test isequal(v1, 1)
    @test isequal(1, v1)
    @test isequal(v4, missing)
    @test isequal(missing, v4)
    @test v1 < 2
    @test 0 < v1
    @test ismissing(v1 < missing)
    @test ismissing(missing < v1)
    @test isless(v1, 2)
    @test isless(0, v1)
    @test isless(v1, missing)
    @test !isless(missing, v1)
    @test isapprox(v3, 1.0)
    @test isapprox(1.0, v3)

    @test hash(v1) == hash(v1.value)
    d = Dict{LabeledValue, Int}(v1 => 1)
    @test haskey(d, v3)

    @test unwrap(v1) === 1
    @test labels(v1) == "a"
    @test labels(v4) == "missing"

    @test sprint(show, v1) == "a"
    @test sprint(show, v4) == "missing"
    @test sprint(show, MIME("text/plain"), v1) == "1 => a"

    @test convert(String, v1) == "a"
end

@testset "LabeledArray" begin
    vals = repeat(1:3, inner=2)
    lbls = Dict(i=>string(i) for i in 1:3)
    x = LabeledArray{Int, 1}(vals, lbls)
    @test eltype(x) == LabeledValue{Int}
    @test size(x) == (6,)
    @test IndexStyle(typeof(x)) == IndexLinear()

    @test x[1] === LabeledValue(1, lbls)
    @test x[2:3] == [1, 2]
    @test x[isodd.(1:6)] == [1, 2, 3]
    @test_throws ArgumentError LabeledArray(string.(vals), Dict{String,String}())

    vals1 = [1, 2, missing]
    x1 = LabeledArray(vals1, convert(Dict{Union{Int,Missing},String}, lbls))
    @test isequal(x1[3], missing)
    @test x1[3] == "missing"

    @test refarray(x) === x.values
    @test refarray(view(x, [1, 3])) == 1:2
    @test refarray(reshape(x, 3, 2)) == reshape(x.values, 3, 2)
    @test collect(x) == x

    x2 = LabeledArray([1:3 1:3], lbls)
    ra = refarray(collect(x2))
    @test size(ra) == (3, 2)
    @test ra isa Matrix{Int64}

    lbs = labels(x)
    @test size(lbs) == size(x)
    @test IndexStyle(typeof(lbs)) == IndexLinear()
    @test lbs[1] == "1"
    @test eltype(typeof(lbs)) == String
    @test length(lbs) == length(x)

    v = view(x, 1:3)
    lbs = labels(v)
    @test size(lbs) == size(v)
    @test lbs[3] == "2"

    lbs1 = labels(x1)
    @test lbs1[3] == "missing"

    y = LabeledArray(copy(vals), copy(lbls))
    @test x == y
    y.labels[4] = "4"
    @test x == y
    y.values[1] = 2
    @test x != y

    @test x == vals
    @test vals == x
    @test x == string.(vals)
    @test string.(vals) == x

    @test copy(x) == x
    @test typeof(copy(x)) == typeof(x)

    @test convert(Vector{Int}, x) === x.values
    @test convert(Vector{Int16}, x) == x.values

    c = CategoricalArray(labels(x))
    @test c == string.(vals)
end


text:
function tvcred_parameterize!(m::AbstractDSGEModel)
    # Parameterize
    for i in 1:5
        adj = (i == 2 || i == 3) ? 0.25 : 1.
        ModelConstructors.set_regime_val!(m[:σ_g], i, adj * m[:σ_g].value)
        ModelConstructors.set_regime_val!(m[:σ_b], i, adj * m[:σ_b].value)
        ModelConstructors.set_regime_val!(m[:σ_μ], i, adj * m[:σ_μ].value)
        ModelConstructors.set_regime_val!(m[:σ_ztil], i, adj * m[:σ_ztil].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_f], i, adj * m[:σ_λ_f].value)
        ModelConstructors.set_regime_val!(m[:σ_λ_w], i, adj * m[:σ_λ_w].value)
        ModelConstructors.set_regime_val!(m[:σ_r_m], i, adj * m[:σ_r_m].value)
        ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, adj * m[:σ_σ_ω].value)
        ModelConstructors.set_regime_val!(m[:σ_μ_e], i, adj * m[:σ_μ_e].value)
        ModelConstructors.set_regime_val!(m[:σ_γ], i, adj * m[:σ_γ].value)
        ModelConstructors.set_regime_val!(m[:σ_π_star], i, adj * m[:σ_π_star].value)
        ModelConstructors.set_regime_val!(m[:σ_lr], i, adj * m[:σ_lr].value)
        ModelConstructors.set_regime_val!(m[:σ_z_p], i, adj * m[:σ_z_p].value)
        ModelConstructors.set_regime_val!(m[:σ_tfp], i, adj * m[:σ_tfp].value)
        err_adj = (i in 2:4) ? 10. : adj
        ModelConstructors.set_regime_val!(m[:σ_gdpdef], i, err_adj * m[:σ_gdpdef].value)
        ModelConstructors.set_regime_val!(m[:σ_corepce], i, err_adj * m[:σ_corepce].value)
        ModelConstructors.set_regime_val!(m[:σ_gdp], i, adj * m[:σ_gdp].value)
        ModelConstructors.set_regime_val!(m[:σ_gdi], i, adj * m[:σ_gdi].value)

        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, adj * m[Symbol("σ_r_m$(j)")])
        end
    end

    ModelConstructors.set_regime_val!(m[:σ_z_p], 2, 0.)
    ModelConstructors.set_regime_val!(m[:σ_z_p], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_r_m], 2, 10. * m[:σ_r_m].value)
    ModelConstructors.set_regime_val!(m[:σ_r_m], 3, 10. * m[:σ_r_m].value)
    for i = 2:3
        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, m[Symbol("σ_r_m$(j)")].value)
        end
    end

    # Set non-standard shocks
    ModelConstructors.set_regime_val!(m[:σ_φ], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 2, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 3, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 4, 400.)
    ModelConstructors.set_regime_val!(m[:σ_φ], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_φ_prop], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 2, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 3, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 4, 4.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc_prop], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ygap], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ygap], 2, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ygap], 3, 20.)
    ModelConstructors.set_regime_val!(m[:σ_ygap], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ygap], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_pgap], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_pgap], 2, 0.)
    ModelConstructors.set_regime_val!(m[:σ_pgap], 3, 20.)
    ModelConstructors.set_regime_val!(m[:σ_pgap], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_pgap], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 2, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 3, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 4, 5.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], 5, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 2, 2.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 4, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid_prop], 5, 0.)
    m[:σ_biidc1].prior = m[:σ_biidc].prior
    m[:σ_biidc1].valuebounds = m[:σ_biidc].valuebounds
    ModelConstructors.set_regime_val!(m[:σ_biidc1], 1, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc1], 2, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc1], 3, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc1], 4, ModelConstructors.regime_val(m[:σ_biidc], 2))
    ModelConstructors.set_regime_val!(m[:σ_biidc1], 5, 0.)
end

function set_regime_vals_fnct(m, n)
    if n > 6
        for p in m.parameters
            regime_switching = !isempty(p.regimes)
            if regime_switching
                if haskey(p.regimes, :value)
                    for i in 6:n
                        ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 1))
                    end
                end
            end
        end
    end
end


text:
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using Distributions

include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))
using CompositeDistributions
import Compat.view

# Set Prior for Population Parameters
d1 = Rayleigh(0.3)
d2 = MvNormal([0.0, 1.0], ones(2))
param_prior = CompositeDist( ContinuousDistribution[d1,d2] )

# Code to generate simulated data given array of model parameters
num_data_default = 100
num_outputs = 2
function gen_data(theta::Array, n::Integer = num_data_default)
  data = Array{Float64}(3,n)
  data[1,:] = rand(Rayleigh(theta[1]),n)
  data[2:3,:] = rand(MvNormal(theta[2:3],ones(2)),n)
  return data
end

# Function to adjust originally proposed model parameters, so that they will be valid
function normalize_theta13_pos!(theta::Array)
 theta[1] = abs(theta[1])
 theta[3] = abs(theta[3])
 return theta
end

# Function to test if the proposed model parameter are valid
function is_valid_theta13_pos(theta::Array)
  theta[1]>0.0  && theta[3]>0.0
end

# True Population Parameters
theta_true = [0.3, 0.0, 1.0]

# Tell ABC what it needs to know for a simulation
abc_plan = abc_pmc_plan_type(gen_data,ABC.calc_summary_stats_mean_var,ABC.calc_dist_max, param_prior; is_valid=is_valid_theta13_pos,num_max_attempt=10000);

# Generate "true/observed data"
data_true = gen_data(theta_true)   # Draw "real" data from same model as for analysis
ss_true = abc_plan.calc_summary_stats(data_true)
#println("theta= ",theta_true," ss= ",ss_true, " d= ", 0.)

# Run ABC simulation
@time pop_out = run_abc(abc_plan,ss_true);

#=
# Optional plotting of results
using PyPlot

hist(pop_out.weights*length(pop_out.weights));
hist(pop_out.dist);

num_grid_x = 100
num_grid_y = 100
limit = 1.0
x = collect(linspace(theta_true[1]-limit,theta_true[1]+limit,num_grid_x));
y = collect(linspace(theta_true[2]-limit,theta_true[2]+limit,num_grid_y));
z = zeros(Float64,(num_param,length(x),length(y)))
for i in 1:length(x), j in 1:length(y)
    z[1,i,j] = x[i]
    z[2,i,j] = y[j]
end
z = reshape(z,(num_param,length(x)*length(y)))
zz = [ ABC.pdf(ABC.GaussianMixtureModelCommonCovar(pop_out.theta,pop_out.weights,ABC.cov_weighted(pop_out.theta',pop_out.weights)),vec(z[:,i])) for i in 1:size(z,2) ]
zz = reshape(zz ,(length(x),length(y)));
levels = [exp(-0.5*i^2)/sqrt(2pi^num_param) for i in 0:5];
PyPlot.contour(x,y,zz',levels);
plot(pop_out.theta[1,:],pop_out.theta[2,:],".");
=#


text:
# Check if two arrays of edges are the same, regardless of directions
function isSame(primsEdges, primsMST)
    for edge in primsEdges
        if(!(in(Edge(edge.dst, edge.src), primsMST)) && !(in(edge, primsMST)))
            return false
        end
    end
    true
end

# Benchmarking PrimsI
@benchmark parallelPrims(graph, numberOfVertices) setup=(graph = randomWeightedGraph(numberOfVertices, connections)) samples=samplesToRun seconds=maxRuntime

# Benchmarking PrimsII
@benchmark prims(graph, numberOfThreads) setup=(graph = randomWeightedGraph(numberOfVertices, connections)) samples=samplesToRun seconds=maxRuntime

# Benchmarking A*
@benchmark astar(graph, graphWidth, startV, endV) setup=(graph = makeMaze(makeGrid(graphWidth))) samples=samplesToRun seconds=maxRuntime


text:
mutable struct Transit_Struct{T}
  # Structure to hold arrays and other quantities for computing transit:
  r       :: T           # radius ratio
  b       :: T           # impact parameter
  u_n     :: Array{T,1}  # limb-darkening coefficients
  n       :: Int64       # number of limb-darkening coefficients
  v_max   :: Int64       # maximum coefficient in I_v and J_v
  c_n     :: Array{T,1}  # Green's basis coefficients
  sn      :: Array{T,1}  # Green's basis terms
  Iv      ::  Array{T,1} # integral over sin(phi)^{2v}
  Jv      ::  Array{T,1} # integral over sin(phi)^{2v} (1-sin(phi)^2/k^2)^{3/2}
  Kv      ::  Array{T,1} # integral over cos(phi)^{2v}
  Lv      ::  Array{T,1} # integral over cos(phi)^{2v} (1-sin(phi)^2/k^2)^{3/2}
  grad    :: Bool        # true for gradient; false for no gradient
  dIvdk   ::  Array{T,1} # derivative of I_v with respect to k
  dJvdk   ::  Array{T,1} # derivative of J_v with respect to k
  dKvdk   ::  Array{T,1} # derivative of K_v with respect to k
  dLvdk   ::  Array{T,1} # derivative of L_v with respect to k
  s2_grad :: Array{T,1}  # gradient of s2 with respect to (r,b) - this is S_1 = s_2 from starry
  dsndr   ::  Array{T,1} # derivatives of Green's basis with respect to r
  dsndb   ::  Array{T,1} # derivatives of Green's basis with respect to b
  dcdu    :: Array{T,2}  # derivatives c_n with respect to u_n
  dfdrbc  :: Array{T,1}  # derivative of flux with respect to c_n
  dfdrbu  :: Array{T,1}  # derivative of flux with respect to u_n
  nmax    :: Int64       # maximum number of terms in series expansions of I_v and J_v
  Iv_coeff :: Array{T,1} # coefficients for series expansion of I_v
  Jv_coeff :: Array{T,3} # coefficients for series expansion of J_v
  dJvdk_coeff :: Array{T,3} # coefficients for series expansion of dJ_v/dk
  k2      :: T           # k^2 = (1-(r-b)^2)/(4br)
  k       :: T           # k = sqrt(k^2)
  kc      :: T           # k_c = sqrt(1-k^2) (unless k > 1, then it is sqrt(1-1/k^2))
  kck     :: T           # k_c * k
  kap0    :: T           # kappa = sin^{-1}(k) (=kappa_0 in M&A)
  pimkap1 :: T           # pi-kappa_1
  kite_area2 :: T        # 2*A_kite
  Eofk    :: T           # E(k^2) is complete elliptic integral of first kind
  Em1mKdm :: T           # (E(m)-(1-m)K(m))/m is complete elliptic integral with m=k^2
  onembmr2:: T           # 1-(b-r)^2
  onembpr2:: T           # 1-(b+r)^2
  onembmr2inv :: T       # 1/(1-(b-r)^2)
  sqonembmr2  :: T       # sqrt(1-(b-r)^2)
  fourbr  :: T           # 4*b*r
  sqbr    :: T           # sqrt(b*r)
  sqbrinv :: T           # 1/sqrt(b*r)
  fourbrinv  :: T           # 1/(4*b*r)
  k2inv   :: T           # 1/k^2
  kc2     :: T           # 1-k^2 (or 1-1/k^2 if k > 1)
  sqrt1mr2:: T           # sqrt(1-r^2)
  den     :: T           # 1/(c_1 + c_2*2/3)
  third   :: T           # 1/3
  twothird:: T           # 2/3
  sqr1mr  :: T           # sqrt(r*(1-r)) if r < 1
  bincoeff:: Array{T,2}  # binomial(n,i)
end

using SpecialFunctions
include("compute_c_n_struct.jl")
include("IJv_coeff.jl")

function transit_init(r::T,b::T,u_n::Array{T,1},grad::Bool) where {T <: Real}
# Initializs a transit structure.
n = length(u_n)
if iseven(n)
  v_max = round(Int64,n/2)+2
else
  v_max = round(Int64,(n-1)/2)+2
end
nmax = 50
trans = Transit_Struct{T}(r,b,u_n,n,v_max,
  zeros(T,n+1),    # c_n
  zeros(T,n+1),    # sn
  zeros(T,v_max+1),# Iv
  zeros(T,v_max+1),# Jv
  zeros(T,v_max+1),# Kv
  zeros(T,v_max+1),# Lv
  grad,            # grad
  zeros(T,v_max+1),# dIvdk
  zeros(T,v_max+1),# dJvdk
  zeros(T,v_max+1),# dKvdk
  zeros(T,v_max+1),# dLvdk
  zeros(T,2),      # s2_grad
  zeros(T,n+1),    # dsndr
  zeros(T,n+1),    # dsndb
  zeros(T,n+1,n),  # dcdu
  zeros(T,n+3),    # dfdrbc
  zeros(T,n+2),    # dfdrbu
  nmax,		   # nmax
  zeros(T,nmax),   # Iv_coeff
  zeros(T,2,2,nmax),   # Jv_coeff for k^2 < 1 & k^2 > 1; v_max & v_max-1; series coefficients
  zeros(T,2,2,nmax),   # dJvdk_coeff for k^2 < 1 & k^2 > 1; v_max & v_max-1; series coefficients
  zero(T),         # k^2
  zero(T),         # k
  zero(T),         # k_c
  zero(T),         # k_c*k
  zero(T),         # kappa/kap0
  zero(T),         # pimkap1/pi-kappa_1
  zero(T),         # kite_area2/2*A_kite
  zero(T),         # E(k^2)
  zero(T),         # (E(m)-(1-m)K(m))/m
  zero(T),         # 1-(b-r)^2
  zero(T),         # 1-(b+r)^2
  zero(T),         # 1/(1-(b-r)^2)
  zero(T),         # sqrt(1-(b-r)^2)
  zero(T),         # 4*b*r
  zero(T),	   # sqrt(b*r)
  zero(T),         # 1/sqrt(b*r)
  zero(T),         # 1/(4*b*r)
  zero(T),         # 1/k^2
  zero(T),	   # 1-k^2 (or 1-1/k^2 if k > 1)
  zero(T),         # sqrt(1-r^2)
  zero(T),         # den = 1/(pi*(c[1] + c[2]*2/3))
  one(T)/3,        # 1/3
  convert(T,2)/3,  # 2/3
  zero(T),         # sqrt(r*(1-r))
  zeros(v_max+1,v_max) # binomial(n,i)
)
# Initialize the series coefficients for I_{v_max}:
Iv_series_coeff!(trans)
# Initialize the series coefficients for J_{v_max} and J_{v_max-1}, and if
# t.grad is true, will also compute coefficients for dJ/dk_{v_max} and _{v_max-1}:
dJvdk_series_coeff!(trans)
# Compute binomial coefficients:
for n=1:trans.v_max
  trans.bincoeff[1,n]=one(T)
  trans.bincoeff[n+1,n]=one(T)
  for i=1:n-1
    if n == 1
      trans.bincoeff[i+1,n]=binomial(n,i)
    else
      trans.bincoeff[i+1,n]=trans.bincoeff[i,n-1]+trans.bincoeff[i+1,n-1]
    end
  end
end
if grad
  compute_c_n_grad!(trans)
else
  compute_c_n!(trans)
end
trans.den = inv(pi*(trans.c_n[1]+trans.twothird*trans.c_n[2]))  # for c_2 and above, the flux is zero.
return trans
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# Autogenerated wrapper script for cilkrts_jll for x86_64-linux-gnu
export libcilkrts

JLLWrappers.@generate_wrapper_header("cilkrts")
JLLWrappers.@declare_library_product(libcilkrts, "libcilkrts.so.5")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcilkrts,
        "lib/libcilkrts.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
### Histograms of recurrence structures

# Add one item to position `p` in the histogram `h` that has precalculated length `n`
# - update the histogram and return its new length
@inline function extendhistogram!(h::Vector{Int}, n::Int, p::Int)
    if p > n
        append!(h, zeros(p-n))
        n = p
    end
    h[p] += 1
    return n
end

@inline function extend_skeleton_histogram!(h::Matrix{Int}, r::Int, c::Int, p::Int)
    h[1,end] = p
    h = hcat(h, [0, c, r]) # revert order for c and r, due to the initial transpose
    return h
end

# Calculate the histograms of segments and distances between segments
# from the indices of rows and columns/diagonals of the matrix
# `theiler` is used for histograms of vertical structures
# `distances` is used to simplify calculations of the distances are not wanted
function _linehistograms(rows::T, cols::T, lmin::Integer, theiler::Integer,
     distances::Bool=true) where {T<:AbstractVector{Int}}

    # check bounds
    n = length(rows)
    if length(cols) != n
        throw(ErrorException("mismatch between number of row and column indices"))
    end
    # histogram for segments
    bins = [0]
    nbins = 1
    # histogram for distances between segments
    bins_d = [0]
    nbins_d = 1
    # find the first point outside the Theiler window
    firstindex = 1
    while (firstindex ≤ n) && (abs(rows[firstindex]-cols[firstindex]) < theiler)
        firstindex += 1
    end
    if firstindex > n
        return ([0],[0])
    end
    # Iterate over columns
    cprev = cols[firstindex]
    r1 = rows[firstindex]
    rprev = r1 - 1 # coerce that (a) is not hit in the first iteration
    dist = 0 # placeholder for distances between segments
    @inbounds for i in firstindex:n
        r = rows[i]
        c = cols[i]
        if abs(r-c) ≥ theiler
            # Search the second and later segments in the column
            if c == cprev
                if r-rprev != 1 # (a): there is a separation between rprev and r
                    # update histogram of segments
                    current_vert = rprev-r1+1
                    if current_vert ≥ lmin
                        nbins = extendhistogram!(bins, nbins, current_vert)
                    end
                    # update histogram of distances if it there were at least
                    # two previous segments in the column
                    if distances
                        halfline = div(current_vert, 2)
                        if dist != 0
                            nbins_d = extendhistogram!(bins_d, nbins_d, dist+halfline)
                        end
                        # update the distance
                        dist = r-rprev+halfline-1
                    end
                    r1 = r # update the start of the next segment
                end
                rprev = r  # update the previous position
            else # hit in the first point of a new column
                # process the last fragment of the previous column
                current_vert = rprev-r1+1
                if current_vert ≥ lmin
                    nbins = extendhistogram!(bins, nbins, current_vert)
                end
                if  distances
                    if dist!= 0
                        halfline = div(current_vert, 2)
                        nbins_d = extendhistogram!(bins_d, nbins_d, dist+halfline)
                    end
                    dist = 0
                end
                # initialize values for searching new fragments
                cprev = c
                r1 = r
                rprev = r
            end
        end
    end
    # Process the latest fragment
    current_vert = rprev-r1+1
    if current_vert ≥ lmin
        nbins = extendhistogram!(bins, nbins, current_vert)
    end
    if  distances && dist != 0
        halfline = div(current_vert, 2)
        nbins_d = extendhistogram!(bins_d, nbins_d, dist+halfline)
    end
    return (bins, bins_d)
end

deftheiler(x::Union{RecurrenceMatrix,JointRecurrenceMatrix}) = 1
deftheiler(x::CrossRecurrenceMatrix) = 0
deftheiler(x::SparseMatrixCSC) = 1
function diagonalhistogram(x::Union{ARM,SparseMatrixCSC}; lmin::Integer=2, theiler::Integer=deftheiler(x), kwargs...)
    (theiler < 0) && throw(ErrorException(
        "Theiler window length must be greater than or equal to 0"))
    (lmin < 1) && throw(ErrorException("lmin must be 1 or greater"))
    m,n = size(x)
    rv = rowvals(x)[:]
    dv = colvals(x) .- rowvals(x)
    loi_hist = Int[]
    if issymmetric(x)
        # If theiler==0, the LOI is counted separately to avoid duplication
        if theiler == 0
            if all(isequal(1),x.nzval)==false #CrossRecurrenceMatrix won't work if all values aren't equal to 1 and theiler=0
                xtmp = deepcopy(x)
                xtmp.nzval .= 1
                loi_hist = verticalhistograms(CrossRecurrenceMatrix(hcat(diag(xtmp,0))),lmin=lmin)[1]
            else
                loi_hist = verticalhistograms(CrossRecurrenceMatrix(hcat(diag(x,0))),lmin=lmin)[1]
            end
        end
        inside = (dv .< max(theiler,1))
        f = 2
    else
        inside = (abs.(dv) .< theiler)
        f = 1
    end
    # remove Theiler window and short by diagonals
    deleteat!(rv, inside)
    deleteat!(dv, inside)
    rv = rv[sortperm(dv)]
    dv = sort!(dv)
    dh = f.*_linehistograms(rv, dv, lmin, 0, false)[1]
    # Add frequencies of LOI if suitable
    if (nbins_loi = length(loi_hist)) > 0
        nbins = length(dh)
        if nbins_loi > nbins
            loi_hist[1:nbins] .+= dh
            dh = loi_hist
        else
            dh[1:nbins_loi] .+= loi_hist
        end
    end
    return dh
end

function verticalhistograms(x::Union{ARM,SparseMatrixCSC};
    lmin::Integer=2, theiler::Integer=deftheiler(x), distances=true, kwargs...)
    (theiler < 0) && throw(ErrorException(
        "Theiler window length must be greater than or equal to 0"))
    (lmin < 1) && throw(ErrorException("lmin must be 1 or greater"))
    m,n=size(x)
    rv = rowvals(x)
    cv = colvals(x)
    return _linehistograms(rv,cv,lmin,theiler,distances)
end

vl_histogram(x; kwargs...) = verticalhistograms(x; kwargs...)[1]
rt_histogram(x; kwargs...) = verticalhistograms(x; kwargs...)[2]

"""
    recurrencestructures(x::AbstractRecurrenceMatrix;
                             diagonal=true,
                             vertical=true,
                             recurrencetimes=true,
                             kwargs...)

Return a dictionary with the
histograms of the recurrence structures contained in the recurrence matrix `x`,
with the keys `"diagonal"`, `"vertical"` or
`"recurrencetimes"`, depending on what keyword arguments are given as `true`.

## Description
Each item of the dictionary is a vector of integers, such that the `i`-th
element of the vector is the number of lines of length `i` contained in `x`.

* `"diagonal"` counts the diagonal lines, i.e. the recurrent trajectories.
* `"vertical"` counts the vertical lines, i.e. the laminar states.
* `"recurrencetimes"` counts the vertical distances between recurrent states,
    i.e. the recurrence times.

All the points of the matrix are counted by default. The keyword argument
`theiler` can be passed to rule out the lines around the main
diagonal. See the arguments of the function [`rqa`](@ref) for further details.

"Empty" histograms are represented always as `[0]`.

*Notice*: There is not a unique operational definition of "recurrence times". In the
analysis of recurrence plots, usually the  "second type" of recurrence times as
defined by Gao and Cai [1] are considered, i.e. the distance between
consecutive (but separated) recurrent structures in the vertical direction of
the matrix. But that distance is not uniquely defined when the vertical recurrent
structures are longer than one point. The recurrence times calculated here are
the distance between the midpoints of consecutive lines, which is a balanced
estimator of the Poincaré recurrence times [2].

## References

[1] J. Gao & H. Cai. "On the structures and quantification of recurrence plots".
[*Physics Letters A*, 270(1-2), 75–87 (2000)](https://www.sciencedirect.com/science/article/pii/S0375960100003042?via%3Dihub).

[2] N. Marwan & C.L. Webber, "Mathematical and computational foundations of
recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence
Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).
"""
function recurrencestructures(x::Union{ARM,SparseMatrixCSC};
    diagonal=true, vertical=true, recurrencetimes=true,
    theiler=deftheiler(x), kwargs...)

    # Parse arguments for diagonal and vertical structures
    histograms = Dict{String,Vector{Int}}()
    if diagonal
        kw_d = Dict(kwargs)
        haskey(kw_d, :theilerdiag) && (kw_d[:theiler] = kw_d[:theilerdiag])
        histograms["diagonal"] = diagonalhistogram(x; theiler=theiler, kw_d..., lmin=1)
    end
    if vertical || recurrencetimes
        kw_v = Dict(kwargs)
        haskey(kw_v, :theilervert) && (kw_v[:theiler] = kw_v[:theilervert])
        vhist = verticalhistograms(x; theiler=theiler, kw_v..., lmin=1)
        vertical && (histograms["vertical"] = vhist[1])
        recurrencetimes && (histograms["recurrencetimes"] = vhist[2])
    end
    return histograms
end


text:
filename = "math-assignment.txt"

file = open(filename, "w")

for _ in 1:100
    x = rand(5:600)
    y = rand(10:602)
    println(file, x, " + ", y, " = ")
end

close(file)


text:
"""
    function ct_kenya(t,minpoint_y;schoolseffect = 0.,
                    measures_x = 24.,minpoint_x = 55.,returnpoint_x = 254.,schoolsreopen = 319. )

Defines piece-wise linear contact rates.
"""
function ct_kenya(t,minpoint_y;schoolseffect = 0.,
                    measures_x = 24.,minpoint_x = 55.,returnpoint_x = 254.,schoolsreopen = 319. )
    ct = 1.0*(t<measures_x)
    ct += (1.0 - schoolseffect  + ((t-measures_x)*(minpoint_y  - 1.0)/(minpoint_x - measures_x)))*(t ≥ measures_x && t < minpoint_x)
    ct += (minpoint_y - schoolseffect + ((t-minpoint_x)*(1.0   - minpoint_y)/(returnpoint_x-minpoint_x) ))*(t ≥minpoint_x && t < returnpoint_x)
    ct += (1.0 - schoolseffect )*(t ≥ returnpoint_x && t<schoolsreopen)
    ct += 1.0*(t ≥ schoolsreopen)
end

"""
    function twoSEStransmission(du,u,p,t)

In-place creation of vector field for the transmission model
"""
function twoSEStransmission(du,u,p,t)
    #Get state
	S₁,E₁,I₁,R₁,W₁,C₁,S₂,E₂,I₂,R₂,W₂,C₂,C_sero₁,C_sero₂  = u
    #Get parameters
	R₀,α,γ,ϵ,σ,N₁,N₂,ω,schooleffect,ct_min1,ct_min2 = p
    #Adjust for schools being open or shut
    Reff = R₀*(t < 24 || t ≥ 319) + R₀*schooleffect*(t≥ 24 && t< 319)
    #Rate of visiting POIs for lower and higher SES groups
    ct1 = ct_kenya(t,ct_min1)
    ct2 = ct_kenya(t,ct_min2)
    #Force of infection on lower and higher SES individuals visiting POIs
    λ₁ = γ*Reff*(ct1*ϵ*I₁ + ct2*(1-ϵ)*I₂)
    λ₂ = γ*Reff*(ct1*(1-ϵ)*I₁ + ct2*ϵ*I₂)
    #Dynamics
	du[1] = -ct1*(S₁/N₁)*λ₁
	du[2] = ct1*((S₁+σ*W₁)/N₁)*λ₁ - α*E₁
	du[3] = α*E₁ - γ*I₁
	du[4] = γ*I₁ - ω*R₁
	du[5] = ω*R₁ - σ*ct1*(W₁/N₁)*λ₁
	du[6] = ct1*((S₁+σ*W₁)/N₁)*λ₁
	du[7] = -ct2*(S₂/N₂)*λ₂
	du[8] = ct2*((S₂+σ*W₂)/N₂)*λ₂ - α*E₂
	du[9] = α*E₂ - γ*I₂
	du[10] = γ*I₂ - ω*R₂
	du[11] = ω*R₂ - σ*ct2*(W₂/N₂)*λ₂
	du[12] = ct2*((S₂+σ*W₂)/N₂)*λ₂
    du[13] = ct1*(S₁/N₁)*λ₁
    du[14] = ct2*(S₂/N₂)*λ₂
end

"""
    function twoSEStransmission_with_lower_transmissibility(du,u,p,t)

In-place creation of vector field for the transmission model. Model where previously infected are less transmissible.
"""
function twoSEStransmission_with_lower_transmissibility(du,u,p,t)
    #Get state
	S₁,E₁,I₁,R₁,W₁,C₁,S₂,E₂,I₂,R₂,W₂,C₂,C_sero₁,C_sero₂,EW₁,IW₁,EW₂,IW₂ = u
    #Get parameters
	R₀,α,γ,ϵ,σ,N₁,N₂,ω,schooleffect,ct_min1,ct_min2,ι = p
    #Adjust for schools being open or shut
    Reff = R₀*(t < 24 || t ≥ 319) + R₀*schooleffect*(t≥ 24 && t< 319)
    #Rate of visiting POIs for lower and higher SES groups
    ct1 = ct_kenya(t,ct_min1)
    ct2 = ct_kenya(t,ct_min2)
    #Force of infection on lower and higher SES individuals visiting POIs
    λ₁ = γ*Reff*(ct1*ϵ*(I₁ + ι*IW₁) + ct2*(1-ϵ)*(I₂ + ι*IW₂))
    λ₂ = γ*Reff*(ct1*(1-ϵ)*(I₁ + ι*IW₁) + ct2*ϵ*(I₂ + ι*IW₂))
    #Dynamics
	du[1] = -ct1*(S₁/N₁)*λ₁
	du[2] = ct1*(S₁/N₁)*λ₁ - α*E₁
	du[3] = α*E₁ - γ*I₁
	du[4] = γ*I₁ + γ*IW₁ - ω*R₁
	du[5] = ω*R₁ - σ*ct1*(W₁/N₁)*λ₁
	du[6] = ct1*((S₁+σ*W₁)/N₁)*λ₁
	du[7] = -ct2*(S₂/N₂)*λ₂
	du[8] = ct2*(S₂/N₂)*λ₂ - α*E₂
	du[9] = α*E₂ - γ*I₂
	du[10] = γ*I₂ + γ*IW₂ - ω*R₂
	du[11] = ω*R₂ - σ*ct2*(W₂/N₂)*λ₂
	du[12] = ct2*((S₂+σ*W₂)/N₂)*λ₂
    du[13] = ct1*(S₁/N₁)*λ₁
    du[14] = ct2*(S₂/N₂)*λ₂
    du[15] = ct1*σ*(W₁/N₁)*λ₁ - α*EW₁
    du[16] = α*EW₁ - γ*IW₁
    du[17] = ct2*σ*(W₂/N₂)*λ₂ - α*EW₂
    du[18] = α*EW₂ - γ*IW₂
end

function build_prob()
    u0 = [4.3e6*0.7 - 100.,100.,0.,0.,0.,0.,
            4.3e6*0.3 - 100.,100.,0.,0.,0.,0.,
            0.0,0.0]
    f = ODEFunction(twoSEStransmission;syms = [:S₁,:E₁,:I₁,:R₁,:W₁,:C₁,:S₂,:E₂,:I₂,:R₂,:W₂,:C₂,:C_sero₁,:C_sero₂])
    return ODEProblem(f,u0,(0,365))
end

function build_prob(::Val{:lower_trans})
    u0 = [4.3e6*0.7 - 100.,100.,0.,0.,0.,0.,
            4.3e6*0.3 - 100.,100.,0.,0.,0.,0.,
            0.0,0.0,
            0.0,0.0,0.0,0.0]
    f = ODEFunction(twoSEStransmission_with_lower_transmissibility;syms = [:S₁,:E₁,:I₁,:R₁,:W₁,:C₁,:S₂,:E₂,:I₂,:R₂,:W₂,:C₂,:C_sero₁,:C_sero₂,:EW₁,:IW₁,:EW₂,:IW₂])
    return ODEProblem(f,u0,(0,365))
end

"""
function simple_conv(incidence,kern)

Direct implementation of convolution [incidence ⋆ kern](t)
"""
function simple_conv(incidence,kern)
    z = similar(incidence)
    n = length(incidence)
    m = length(kern)
    for t = 1:n
        z[t] = 0.
        for s = 1:(t-1)
            if t-s <= m && t-s >= 1
                z[t] += kern[t-s]*incidence[s]
            end
        end
    end
    return z
end


text:
using Pkg
Pkg.add("Documenter")
#push!(LOAD_PATH,"../src/")
using Documenter, Solaris

makedocs(
    sitename = "Solaris.jl",
)

deploydocs(
    repo = "github.com/ethemidori/Solaris.jl.git",
)


text:
"""
    forwardsub(L,b)

Solve the lower-triangular linear system with matrix `L` and right-hand side
vector `b`.
"""
function forwardsub(L,b)

n = size(L,1)
x = zeros(n)
for i = 1:n
    s = i > 1 ? sum( L[i,j]*x[j] for j=1:i-1 ) : 0
    x[i] = ( b[i] - s ) / L[i,i]
end

return x
end

"""
    backsub(U,b)

Solve the upper-triangular linear system with matrix `U` and right-hand side
vector `b`.
"""
function backsub(U,b)

n = size(U,1)
x = zeros(n)
for i = n:-1:1
    s = i < n ? sum( U[i,j]*x[j] for j=i+1:n ) : 0
    x[i] = ( b[i] - s ) / U[i,i]
end

return x
end

"""
    lufact(A)

Compute the LU factorization of square matrix `A`, returning the factors.
"""
function lufact(A)

n = size(A,1)
L = Matrix(Diagonal(ones(n)))    # puts ones on diagonal
U = copy(A)

# Gaussian elimination
for j = 1:n-1
  for i = j+1:n
    L[i,j] = U[i,j] / U[j,j]   # row multiplier
    U[i,j:n] = U[i,j:n] - L[i,j]*U[j,j:n]
  end
end

return L,triu(U)
end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
using OpenStreetMapGraphs, Test
# execute test with:
# using Pkg; Pkg.test("OpenStreetMapGraphs")

function tests()
    @testset "load_save" begin
        query = overpassquery(53.67,10.16,53.,9.72)
        @test typeof(query) == String
        @test length(query) > 0
    end

    @testset "parse tag" begin
        word = OpenStreetMapGraphs.rulewords("PH off")
        @test typeof(word) == Array{AbstractString, 1}
        @test word == ["off"]
        word2 = OpenStreetMapGraphs.rulewords("PH")
        @test typeof(word2) == Array{AbstractString, 1}
        @test word2 == []
    end

    @testset "maxspeed tags" begin
        @testset "maxspeed" begin
            way = Dict(
                "type" => "way",
                "tags" => Dict(
                    "maxspeed"=> "12",
                    "maxspeed:forward" => "123",
                    "maxspeed:backward" => "1234"
                )
            )
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test maxspeed[:maxspeed] === 12
            @test typeof(maxspeed[:maxspeed]) <: Real

            @test maxspeed[:uvmaxspeed] === 123
            @test typeof(maxspeed[:uvmaxspeed]) <: Real

            @test maxspeed[:vumaxspeed] === 1234
            @test typeof(maxspeed[:vumaxspeed]) <: Real
        end
        @testset "maxspeed:conditional" begin
            way = Dict(
                "type" => "way",
                "tags" => Dict(
                    "maxspeed"=> "12",
                    "maxspeed:conditional" => "30@(22:00-06:00)",
                )
            )
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test length(maxspeed[:maxspeed]) === 168
            @test typeof(maxspeed[:maxspeed]) == Array{Union{Missing,Real},1}
            @test maxspeed[:maxspeed][6] === 30
            @test maxspeed[:maxspeed][7] === 12
            @test maxspeed[:maxspeed][22] === 12
            @test maxspeed[:maxspeed][23] === 30

            way["tags"]["maxspeed:conditional"] = "30 @ (06:00-22:00; Su,PH off)"
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test length(maxspeed[:maxspeed]) === 168
            @test typeof(maxspeed[:maxspeed]) == Array{Union{Missing,Real},1}
            @test maxspeed[:maxspeed][6] === 12
            @test maxspeed[:maxspeed][7] === 30
            @test maxspeed[:maxspeed][22] === 30
            @test maxspeed[:maxspeed][23] === 12
            @test maxspeed[:maxspeed][165] === 12
            @test maxspeed[:maxspeed][166] === 12

            way["tags"]["maxspeed:conditional"] = "30 @ (Mo-Fr 00:00-06:30,22:00-24:00; Sa,Su 00:00-24:00)"
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test length(maxspeed[:maxspeed]) === 168
            @test typeof(maxspeed[:maxspeed]) == Array{Union{Missing,Real},1}
            @test maxspeed[:maxspeed][6] === 30
            @test maxspeed[:maxspeed][7] === 12
            @test maxspeed[:maxspeed][22] === 12
            @test maxspeed[:maxspeed][23] === 30
            @test maxspeed[:maxspeed][165] === 30
            @test maxspeed[:maxspeed][166] === 30

            way["tags"]["maxspeed:conditional"] = "30 @ (Mo-Fr 00:00-06:00;Mo-Sa 22:00-24:00;Su 00:00-24:00)"
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test length(maxspeed[:maxspeed]) === 168
            @test typeof(maxspeed[:maxspeed]) == Array{Union{Missing,Real},1}
            @test maxspeed[:maxspeed][6] === 30
            @test maxspeed[:maxspeed][7] === 12
            @test maxspeed[:maxspeed][22] === 12
            @test maxspeed[:maxspeed][23] === 30
            @test maxspeed[:maxspeed][126] === 12
            @test maxspeed[:maxspeed][126] === 12
            @test maxspeed[:maxspeed][142] === 12
            @test maxspeed[:maxspeed][143] === 30
            @test maxspeed[:maxspeed][165] === 30
            @test maxspeed[:maxspeed][166] === 30

            way["tags"]["maxspeed:conditional"] = "30 @ (Mo-Sa 06:00-22:00);50 @ We"
            maxspeed =  OpenStreetMapGraphs.maxspeed(way)
            @test length(maxspeed[:maxspeed]) === 168
            @test typeof(maxspeed[:maxspeed]) == Array{Union{Missing,Real},1}
            @test maxspeed[:maxspeed][6] === 12
            @test maxspeed[:maxspeed][7] === 30
            @test maxspeed[:maxspeed][22] === 30
            @test maxspeed[:maxspeed][23] === 12
            @test maxspeed[:maxspeed][70] === 50
            @test maxspeed[:maxspeed][166] === 12
        end
    end
end

tests()