text:
a=["126.88669" "17:55:28.559";
"112.820335" "17:55:29.037";
"91.86494" "17:55:29.316";
"91.08089" "17:55:29.552";
"90.91018" "17:55:29.881";
"89.58628" "17:55:30.213";
"89.219666" "17:55:30.527";
"88.87194" "17:55:30.966";
"87.85196" "17:55:31.366";
"86.70168" "17:55:31.805";
"86.45967" "17:55:32.129";
"87.81331" "17:55:32.667";
"88.95815" "17:55:32.997";
"89.598076" "17:55:33.398";
"90.88459" "17:55:33.73";
"92.66952" "17:55:34.177";
"93.01398" "17:55:34.476";
"93.92483" "17:55:34.783";
"95.02025" "17:55:35.098";
"94.21429" "17:55:35.347";
"91.82439" "17:55:35.64";
"86.28223" "17:55:35.905";
"83.442795" "17:55:36.205";
"78.69276" "17:55:36.501";
"76.92684" "17:55:36.842";
"73.44886" "17:55:37.186";
"56.237144" "17:55:37.458";
"59.54601" "17:55:37.821";
"60.77915" "17:55:38.202";
"54.70414" "17:55:38.553";
"49.85453" "17:55:38.918";
"51.20992" "17:55:39.266";
"52.65848" "17:55:39.659";
"50.036953" "17:55:40.035";
"49.78403" "17:55:40.397";
"52.24312" "17:55:40.789";
"52.3873" "17:55:41.158";
"50.091232" "17:55:41.571";
"50.03581" "17:55:41.991";
"51.66831" "17:55:42.324";
"51.65034" "17:55:42.743";
"53.212296" "17:55:43.124";
"58.219124" "17:55:43.57";
"77.74377" "17:55:44.087";
"65.11755" "17:55:44.739";
"145.0894" "17:55:45.367";
"201.49815" "17:55:46.348";
"281.0915" "17:55:51.522";
"166.84415" "17:55:51.929";
"98.68957" "17:55:55.685";
"73.52383" "17:56:05.847";
"63.965332" "17:56:06.656";
"59.05876" "17:56:09.628";
"55.19937" "17:56:13.793";
"55.72312" "17:56:14.791";
"50.995964" "17:56:15.191";
"49.618893" "17:56:18.446";
"51.554836" "17:56:20.064";
"412.13364" "17:56:23.149";
"230.6789" "17:56:51.45";
"62.255943" "17:57:29.169";
"84.60381" "17:57:29.716";
"99.69191" "17:57:34.064";
"104.43569" "17:57:34.537";
"106.3232" "17:57:34.98";
"107.02868" "17:57:35.473";
"107.223816" "17:57:35.858";
"107.169876" "17:57:36.348";
"106.96233" "17:57:36.83";
"106.71725" "17:57:37.233";
"106.404045" "17:57:37.756";
"106.04085" "17:57:38.305";
"105.63985" "17:57:38.861";
"105.244644" "17:57:49.982";
"105.01264" "17:57:50.363";
"104.75897" "17:57:50.645";
"104.4583" "17:57:51.193";
"104.13869" "17:57:51.562";
"103.805084" "17:57:51.85";
"103.45662" "17:57:52.132";
"103.069786" "17:57:52.727";
"102.64776" "17:57:53.306";
"102.19315" "17:57:53.894";
"101.70799" "17:57:54.506";
"101.19382" "17:57:55.063";
"100.651634" "17:57:55.628";
"100.08208" "17:57:56.188";
"99.485275" "17:57:56.758";
"98.86389" "17:57:57.059";
"98.21697" "17:57:57.329";
"97.53991" "17:57:57.886";
"96.83273" "17:57:58.477";
"96.09504" "17:57:59.095";
"95.32621" "17:57:59.672";
"94.515686" "17:58:00.009";
"93.751785" "17:58:00.36";
"92.94315" "17:58:00.913";
"92.09068" "17:58:01.422";
"91.1969" "17:58:01.981";
"90.2673" "17:58:02.528";
"89.348305" "17:58:02.869";
"88.42868" "17:58:03.232";
"87.60162" "17:58:03.601";
"87.05892" "17:58:03.981";
"87.197044" "17:58:04.341";
"89.49103" "17:58:04.731";
"95.34069" "17:58:05.074";
"99.3908" "17:58:05.349";
"101.72479" "17:58:05.655";
"103.11622" "17:58:06.103";
"103.46161" "17:58:06.641";
"103.82472" "17:58:07";
"103.66056" "17:58:07.501";
"103.02916" "17:58:08.018";
"102.25988" "17:58:08.436";
"100.82844" "17:58:08.966";
"97.86634" "17:58:09.453";
"92.90056" "17:58:09.782";
"88.56462" "17:58:10.151";
"88.00033" "17:58:10.458";
"88.12315" "17:58:10.764";
"88.340904" "17:58:11.098";
"88.674355" "17:58:11.452";
"88.94319" "17:58:11.82";
"89.26047" "17:58:12.225";
"89.93408" "17:58:12.552";
"91.648834" "17:58:12.898";
"94.64191" "17:58:13.276";
"97.46883" "17:58:13.591";
"99.5846" "17:58:13.909";
"101.18376" "17:58:14.22";
"101.57606" "17:58:14.757";
"101.15649" "17:58:15.283";
"99.71409" "17:58:15.766";
"96.73128" "17:58:16.177";
"93.31119" "17:58:16.505";
"90.87826" "17:58:16.825";
"90.0079" "17:58:17.126";
"90.17752" "17:58:17.49";
"90.52873" "17:58:17.809";
"91.2447" "17:58:18.182";
"94.39628" "17:58:18.557";
"98.48616" "17:58:18.911";
"100.82646" "17:58:19.234";
"100.638115" "17:58:19.801";
"98.30174" "17:58:20.286";
"93.44149" "17:58:20.672";
"90.5696" "17:58:21.027";
"90.30207" "17:58:21.34";
"90.41641" "17:58:21.694";
"90.858246" "17:58:22.031";
"93.03901" "17:58:22.376";
"98.3454" "17:58:22.727";
"101.62711" "17:58:23.053";
"101.9411" "17:58:23.536";
"100.78429" "17:58:24.007";
"97.33588" "17:58:24.399";
"92.27823" "17:58:24.833";
"91.83477" "17:58:25.186";
"91.997795" "17:58:25.55";
"92.14125" "17:58:25.913";
"92.51152" "17:58:26.304";
"94.99264" "17:58:26.685";
"100.58009" "17:58:27.042";
"103.47047" "17:58:27.343";
"103.37792" "17:58:27.86";
"101.62857" "17:58:28.387";
"97.59933" "17:58:28.79";
"94.43933" "17:58:29.194";
"94.014786" "17:58:29.536";
"94.10075" "17:58:29.915";
"94.24824" "17:58:30.295";
"94.55415" "17:58:30.739";
"94.68708" "17:58:31.339";
"94.72834" "17:58:31.709";
"94.749275" "17:58:32.335";
"94.596924" "17:58:32.935";
"94.3649" "17:58:33.274";
"94.06945" "17:58:33.618";
"93.73108" "17:58:33.969";
"93.30565" "17:58:34.591";
"92.91356" "17:58:34.958";
"92.39469" "17:58:35.286";
"91.80461" "17:58:35.642";
"91.29023" "17:58:35.988";
"91.417625" "17:58:36.366";
"94.31918" "17:58:36.741";
"100.01713" "17:58:37.143";
"103.471405" "17:58:37.462";
"105.63808" "17:58:37.76";
"106.504715" "17:58:38.339";
"106.47781" "17:58:39.452";
"105.93881" "17:58:40.078";
"104.882454" "17:58:40.721";
"103.059715" "17:58:41.28";
"97.86982" "17:58:41.878";
"92.07422" "17:58:42.257";
"92.28134" "17:58:42.646";
"93.13538" "17:58:43.041";
"94.14432" "17:58:43.46"]
using Dates , Plots
function parse_dates(x)
    t=split(x,(':','.'))
    w=map(x->parse(Float64,x),t)
    Dates.Time(w...)
end
function norm_time(x)
    abs_times = x
    start_time = abs_times[1]
    norm_times = Vector(undef,length(abs_times))
    for i in 1:length(abs_times)
        norm_times[i] = abs_times[i]-start_time
    end
    norm_times
end

loss_values = map(x->parse(Float64,x),a[:,1])
time_values = map(x->parse_dates(x),a[:,2])
norm_time_vals = norm_time(time_values)

plot(time_values[1:20],loss_values[1:20],linecolor="blue")


text:
using OVERT

ϵ = 0.1

function tora_dynamics(x::Array{T, 1} where {T <: Real},
                       u::Array{T, 1} where {T <: Real})
    dx1 = x[2]
    dx2 = ϵ * sin(x[3]) - x[1]
    dx3 = x[4]
    dx4 = u[1]
    return [dx1, dx2, dx3, dx4]
end

tora_dim2 = :($ϵ * sin(x3) - x1)

function tora_dynamics_overt(range_dict::Dict{Symbol, Array{T, 1}} where {T <: Real},
                             N_OVERT::Int,
					         t_idx::Union{Int, Nothing}=nothing)
	if isnothing(t_idx)
		v1 = tora_dim2
	else
    	v1 = "$ϵ * sin(x3_$t_idx) - x1_$t_idx"
    	v1 = Meta.parse(v1)
	end

    oA_out = overapprox(v1, range_dict; N=N_OVERT)
    return oA_out, [oA_out.output]
end

function tora_update_rule(input_vars, control_vars, overt_output_vars)
    integration_map = Dict(input_vars[1] => input_vars[2],
                           input_vars[2] => overt_output_vars[1],
                           input_vars[3] => input_vars[4],
                           input_vars[4] => control_vars[1])
    return integration_map
end

tora_input_vars = [:x1, :x2, :x3, :x4]
tora_control_vars = [:u1]

Tora = OvertProblem(
	tora_dynamics,
	tora_dynamics_overt,
	tora_update_rule,
	tora_input_vars,
	tora_control_vars
)

ext: jl


text:
export someprimes, allprimes

type PrimeIt{T<:Integer}
    n1::T
    n2::T
end

type PrimeInfIt{T<:Integer}
    n1::T
end

someprimes(n::Integer) = PrimeIt(convert(typeof(n),2),n)
someprimes(n1::Integer, n2::Integer) = PrimeIt(promote(n1,n2)...)
allprimes(n1::Integer) = PrimeInfIt(n1)
allprimes() = PrimeInfIt(2)

Base.start(pri::PrimeIt) = pri.n1
Base.start(pri::PrimeInfIt) = pri.n1
Base.next(pri::PrimeIt, state)  = (state, nextprime(state))
Base.next(pri::PrimeInfIt, state) = (state, nextprime(state))
Base.done(pri::PrimeIt, state) = state > pri.n2
Base.done(pri::PrimeInfIt, state) = false


text:
abstract type AbstractNormalSample{T} <: ContinuousEBayesSample{T} end

"""
    NormalSample(Z,σ)

An observed sample ``Z`` drawn from a Normal distribution with known variance ``\\sigma^2 > 0``.

```math
Z \\sim \\mathcal{N}(\\mu, \\sigma^2)
```

``\\mu`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``\\mu``.

```jldoctest
julia> NormalSample(0.5, 1.0)          #Z=0.5, σ=1
Z=     0.5 | σ=1.0
```
"""
struct NormalSample{T,S} <: AbstractNormalSample{T}
    Z::T
    σ::S
end

# TODO: Should not need this eventually.
function NormalSample(Z::P, σ::S) where {T,P<:EBInterval{T},S}
    NormalSample{EBInterval{T},S}(Z, σ)
end

function NormalSample(σ::S) where {S}
    NormalSample(missing, σ)
end

"""
    StandardNormalSample(Z)

An observed sample ``Z`` drawn from a Normal distribution with known variance ``\\sigma^2 =1``.

```math
Z \\sim \\mathcal{N}(\\mu, 1)
```

``\\mu`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``\\mu``.

```jldoctest
julia> StandardNormalSample(0.5)          #Z=0.5
Z=     0.5 | σ=1.0
```
"""
struct StandardNormalSample{T} <: AbstractNormalSample{T}
    Z::T
end

StandardNormalSample() = StandardNormalSample(missing)

eltype(Z::AbstractNormalSample{T}) where {T} = T
support(Z::AbstractNormalSample) = RealInterval(-Inf, +Inf)

response(Z::AbstractNormalSample) = Z.Z
var(Z::AbstractNormalSample) = Z.σ^2
var(Z::StandardNormalSample) = one(eltype(response(Z)))
var(Z::StandardNormalSample{Missing}) = 1.0

std(Z::AbstractNormalSample) = Z.σ
std(Z::StandardNormalSample) = one(eltype(response(Z)))
std(Z::StandardNormalSample{Missing}) = 1.0

nuisance_parameter(Z::AbstractNormalSample) = std(Z)

likelihood_distribution(Z::AbstractNormalSample, μ) = Normal(μ, std(Z))

function Base.show(io::IO, Z::AbstractNormalSample{<:Real})
    Zz = response(Z)
    print(io, "Z=")
    print(io, lpad(round(Zz, digits = 4),8))
    print(io, " | ", "σ=")
    print(io, rpad(round(std(Z), digits = 3),5))
end

function Base.show(io::IO, Z::AbstractNormalSample{<:Interval})
    Zz = response(Z)
    print(io, "Z ∈ ")
    show(IOContext(io, :compact => true), Zz)
    print(io, " | ", "σ=")
    print(io, rpad(round(std(Z), digits = 3),5))
end

# Targets

# TODO: Note this is not correct for intervals.
function cf(target::MarginalDensity{<:AbstractNormalSample}, t)
    error_dbn = likelihood_distribution(location(target))
    cf(error_dbn, t)
end

# Conjugate computations
function default_target_computation(::BasicPosteriorTarget,
    ::AbstractNormalSample,
    ::Normal
)
    Conjugate()
end

function marginalize(Z::AbstractNormalSample, prior::Normal)
    prior_var = var(prior)
    prior_μ = mean(prior)
    likelihood_var = var(Z)
    marginal_σ = sqrt(likelihood_var + prior_var)
    Normal(prior_μ, marginal_σ)
end

function posterior(Z::AbstractNormalSample, prior::Normal)
    z = response(Z)
    sigma_squared = var(Z)
    prior_mu = mean(prior)
    prior_A = var(prior)

    post_mean =
        (prior_A) / (prior_A + sigma_squared) * z +
        sigma_squared / (prior_A + sigma_squared) * prior_mu
    post_var = prior_A * sigma_squared / (prior_A + sigma_squared)
    Normal(post_mean, sqrt(post_var))
end

function _set_defaults(
    convexclass::DiscretePriorClass,
    Zs::VectorOrSummary{<:AbstractNormalSample};  #TODO for MultinomialSummary
    hints,
)
    eps = get(hints, :eps, 1e-4)
    prior_grid_length = get(hints, :prior_grid_length, 300)::Integer
    _sample_min, _sample_max = extrema(response.(Zs))
    # This won't handle infinity correctly. TODOO
    # Also TODO once switch from Intervals.jl to IntervalSets.jl
    _sample_min = isa(_sample_min, Interval) ? first(_sample_min) : _sample_min
    _sample_max = isa(_sample_max, Interval) ? last(_sample_max) : _sample_max
    _grid = range(_sample_min - eps; stop = _sample_max + eps, length = prior_grid_length)
    DiscretePriorClass(_grid)
end

function _set_defaults(
    convexclass::GaussianScaleMixtureClass,
    Zs::AbstractVector{<:AbstractNormalSample};  #TODO for MultinomialSummary
    hints,
)
    grid_scaling = get(hints, :grid_scaling, sqrt(2))

    σ_min =  minimum(std.(Zs))./ 10
    _max = maximum(response.(Zs).^2 .-  var.(Zs))
    _σ_max = _max > 0.0 ? 2*sqrt(_max) : 8*σ_min
    σ_max = get(hints, :σ_max, _σ_max) #somewhat redundant computations above.

    npoint = ceil(Int, log2(σ_max/σ_min)/log2(grid_scaling))
    σ_grid = σ_min*grid_scaling.^(0:npoint)

    GaussianScaleMixtureClass(σ_grid)
end

# Target specifics
function Base.extrema(density::MarginalDensity{<:AbstractNormalSample{<:Real}})
    (0.0, 1 / sqrt(2π * var(location(density))))
end


text:
# generate data
if Pkg.installed("JLD") == nothing
    Pkg.add("JLD")
end

using JLD

function gendata(n)

	# covariance matrix
	covar = (arg = randn(n,n); arg*arg'/1000)
    # expected return
    u = rand(Float64,n)/100

#return(covar)

save("randdata.jld","covar", covar, "u",u)

end


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using MeasureTheory
using ParameterHandling
using GalacticOptim
using AdvancedHMC
using AdvancedMH
import Distributions as Dists
using HypercubeTransform
using Optim
using Plots
using ForwardDiff
using JuKeBOX
const Jk = JuKeBOX
using BlackBoxOptim
using GCMAES

load_ehtim()
obs = ehtim.obsdata.load_uvfits("/home/ptiede/Research/ngEHT/ngEHTChallenge2/data/Challenge_2/synthetic_data/M87_GRMHD/eht2022/230/M87_GRMHD_230_122_eht2022_04_allnoise_10s_230GHz_seed-1_netcal.uvfits")
obs.add_scans()
obsflg = obs.flag_uvdist(uv_min=0.9e9)
obsavg = obsflg.avg_coherent(0.0, scan_avg=true)
obsnn = obsavg.add_fractional_noise(0.02)
dlcamp = extract_lcamp(obsnn; count="min")
dcphase = extract_cphase(obsnn; count="min-cut0bl")

lklhd = RadioLikelihood(dlcamp, dcphase)

function make_model(θ)
    (;spin, inc, rpeak, width, β, χ, ι, θg, ξ) = θ
    #(;spin, inc, θg, ξ) = θ
    o = Jk.Observer(1, inc)
    g, acc = Jk.bam(1, spin, 1.0, rpeak, width, β, χ, ι)
    macc = SingleStokes(Jk.SimpleModel(acc, g, o), :I)
    mimg = modelimage(macc; fovx=30.0, fovy=30.0, nx=64, ny=64)
    #fm = flux(mimg)
    return rotated(stretched(mimg, θg, θg), ξ)
end

function test()
    params = transform(tf, rand(dimension(tf)))
    model = make_model(params)
    #heatmap(intensitymap(model, 150.0, 150.0, 256, 256))
end

prior = (
         spin = Dists.Uniform(-0.98, 0.98),
         inc = Dists.Uniform(0.01, π/2),
         rpeak = Dists.Uniform(2.1, 10.0),
         width = Dists.Uniform(1.0, 10.0),
         β = Dists.Uniform(0.0, 0.98),
         χ = Dists.Uniform(0.0, π/2),
         ι = Dists.Uniform(0.0, π/2),
         θg = Dists.Uniform( μas2rad(0.1), μas2rad(10.0)),
         ξ = Dists.Uniform(-π, π),
         #f = Dists.Uniform(0.0, 1.0)
        )

function make_lpost(lklhd, prior, make_model)
    tr = asflat(prior)
    function lpost(x)
        y, logjac = transform_and_logjac(tr, x)
        return logdensity(lklhd, make_model(y))  +
               logjac
    end
end

lp = make_lpost(lklhd, prior, make_model)
nd = dimension(lp.tr)
p = make_model(transform(lp.tr, rand(nd)))
lp(rand(nd))
@btime lp($rand(nd))
∇lpost(x) = ForwardDiff.gradient(lp, x)
∇lpost(rand(nd))
@time ∇lpost(rand(nd));
f = OptimizationFunction((x,p)->-lp(x))
prob = GalacticOptim.OptimizationProblem(f, rand(nd), nothing, lb = fill(0.01, nd), ub = fill(0.99, nd))
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters=10_000, maxtime=1000.0)
f = OptimizationFunction((x,p)->-lp(x), GalacticOptim.AutoForwardDiff())
prob = GalacticOptim.OptimizationProblem(f, sol.u, nothing, lb = fill(0.01, nd), ub = fill(0.99, nd))
sol = solve(prob, GCMAESOpt())
@show sol.minimum
xopt = transform(lp.tr, sol.u)

mopt = make_model(xopt)
dvis = extract_vis(obsavg)
ac = arrayconfig(dvis)
u,v = getuv(ac)
scatter(hypot.(u,v), abs.(visibility.(dvis.data)))
scatter!(hypot.(u,v), amplitude.(Ref(mopt), u, v))
img = intensitymap(mopt, μas2rad(100.0), μas2rad(100.0), 64, 64)

x,y = imagepixels(img)
heatmap(rad2μas(x), rad2μas(y), img, xflip=true, aspect_ratio=:equal, size=(600,500))


text:
module QPControlTests

using Test
using LinearAlgebra
using Random
using QPControl
using RigidBodyDynamics
using RigidBodyDynamics.Contact
using ValkyrieRobot.BipedControlUtil
using ValkyrieRobot
using StaticArrays
using Rotations
using MathOptInterface
using OSQP
using OSQP.MathOptInterfaceOSQP: OSQPSettings
using Parametron
using StaticUnivariatePolynomials

import Parametron: mock_model, setdirty!

const QPC = QPControl
const RBD = RigidBodyDynamics
const MOI = MathOptInterface

macro test_noalloc(expr)
    quote
        let
            $expr
            allocs = @allocated $expr
            @test allocs == 0
        end
    end |> esc
end

function defaultoptimizer()
    optimizer = OSQP.Optimizer()
    MOI.set(optimizer, OSQPSettings.Verbose(), false)
    MOI.set(optimizer, OSQPSettings.EpsAbs(), 1e-8)
    MOI.set(optimizer, OSQPSettings.EpsRel(), 1e-16)
    MOI.set(optimizer, OSQPSettings.MaxIter(), 20000)
    MOI.set(optimizer, OSQPSettings.AdaptiveRhoInterval(), 25) # required for deterministic behavior
    optimizer
end

include("trajectories.jl")

using QPControl.Trajectories
using RigidBodyDynamics.PDControl
include("tasks.jl")
include("controller.jl")
include("notebooks.jl")

end # module


text:
export rule

@rule NormalMeanPrecision(:τ, Marginalisation) (q_out::Any, q_μ::Any) = begin
    θ = 2 / (var(q_out) + var(q_μ) + abs2(mean(q_out) - mean(q_μ)))
    α = convert(typeof(θ), 1.5)
    return Gamma(α, θ)
end

@rule NormalMeanPrecision(:τ, Marginalisation) (q_out_μ::Any, ) = begin
    m, V = mean_cov(q_out_μ)
    θ = 2 / (V[1,1] - V[1,2] - V[2,1] + V[2,2] + abs2(m[1] - m[2]))
    α = convert(typeof(θ), 1.5)
    return Gamma(α, θ)
end


text:
""" Graphviz support for Catlab's graph types.
"""
module GraphvizGraphs
export parse_graphviz, to_graphviz, to_graphviz_property_graph

using StaticArrays: StaticVector, SVector

using ...Theories
using ...Graphs, ...CategoricalAlgebra.Subobjects
import ..Graphviz

# Property graphs
#################

""" Parse Graphviz output in JSON format.

Returns a property graph with graph layout and other metadata. Each node has a
position and size.

All units are in points. Note that Graphviz has 72 points per inch.
"""
function parse_graphviz(doc::AbstractDict)::AbstractPropertyGraph
  graph = doc["directed"] ? PropertyGraph{Any}() : SymmetricPropertyGraph{Any}()
  nsubgraphs = doc["_subgraph_cnt"] # Subgraphs are ignored.

  # Graph-level layout: bounds and padding.
  # It seems, but is not documented, that the first two numbers in the Graphviz
  # bounding box are always zero.
  set_gprops!(graph,
    bounds = SVector{2}(parse_vector(doc["bb"])[3:4]),
    pad = 72*parse_point(get(doc, "pad", "0,0")),
    rankdir = get(doc, "rankdir", "TB"),
  )

  # Add vertex for each Graphviz node.
  node_keys = ("id", "name", "comment", "label", "shape", "style")
  for node in doc["objects"][nsubgraphs+1:end]
    props = Dict{Symbol,Any}(
      Symbol(k) => node[k] for k in node_keys if haskey(node, k))
    props[:position] = parse_point(node["pos"])
    props[:size] = 72*SVector(parse(Float64, node["width"]),
                              parse(Float64, node["height"]))
    add_vertex!(graph, props)
  end

  # Add edge for each Graphviz edge.
  edge_keys = ("id", "comment", "label", "xlabel", "headlabel", "taillabel",
               "headport", "tailport")
  for edge in doc["edges"]
    if get(edge, "style", nothing) == "invis"
      # Omit invisible edges, which are used to tweak the layout in Graphviz.
      continue
    end
    props = Dict{Symbol,Any}(
      Symbol(k) => edge[k] for k in edge_keys if haskey(edge, k))
    props[:spline] = parse_spline(edge["pos"])
    src = edge["tail"] - nsubgraphs + 1
    tgt = edge["head"] - nsubgraphs + 1
    add_edge!(graph, src, tgt, props)
  end

  graph
end

""" Parse an array of floats in Graphviz's comma-separated format.
"""
parse_vector(s::AbstractString) = [ parse(Float64, x) for x in split(s, ",") ]

""" Parse Graphviz point.

http://www.graphviz.org/doc/info/attrs.html#k:point
"""
parse_point(s::AbstractString) = SVector{2}(parse_vector(s))

""" Parse Graphviz spline.

In Graphviz, a "spline" is a cubic B-spline of overlapping cubic Bezier curves.
It consists of 3n+1 points, where n is the number of Bezier curves.

http://www.graphviz.org/doc/info/attrs.html#k:splineType
https://web.archive.org/web/20170418034924/http://www.graphviz.org/content/how-convert-b-spline-bezier
"""
function parse_spline(spline::AbstractString)
  points = StaticVector{2,Float64}[]
  start, stop = nothing, nothing
  for s in split(spline, " ")
    if startswith(s, "s,"); start = parse_point(s[3:end])
    elseif startswith(s, "e,"); stop = parse_point(s[3:end])
    else push!(points, parse_point(s)) end
  end
  # Prefer explicit start or end points to the spline start and end points.
  # Thus, endpoints may pass into the node shape but should not fall short.
  if !isnothing(start); points[0] = start end
  if !isnothing(stop); points[end] = stop end
  points
end

""" Convert a property graph to a Graphviz graph.

This method is usually more convenient than direct AST manipulation for creating
simple Graphviz graphs. For more advanced features, like nested subgraphs, you
must use the Graphviz AST.
"""
function to_graphviz(g::AbstractPropertyGraph)::Graphviz.Graph
  gv_name(v::Int) = "n$v"
  gv_path(e::Int) = [gv_name(src(g,e)), gv_name(tgt(g,e))]

  # Add Graphviz node for each vertex in property graph.
  stmts = Graphviz.Statement[]
  for v in vertices(g)
    push!(stmts, Graphviz.Node(gv_name(v), vprops(g, v)))
  end

  # Add Graphviz edge for each edge in property graph.
  is_directed = !(g isa SymmetricPropertyGraph)
  for e in edges(g)
    # In undirected case, only include one edge from each pair.
    if is_directed || (e <= inv(g,e))
      push!(stmts, Graphviz.Edge(gv_path(e), eprops(g, e)))
    end
  end

  attrs = gprops(g)
  Graphviz.Graph(
    name = get(attrs, :name, "G"),
    directed = is_directed,
    prog = get(attrs, :prog, is_directed ? "dot" : "neato"),
    stmts = stmts,
    graph_attrs = Graphviz.as_attributes(get(attrs, :graph, Dict())),
    node_attrs = Graphviz.as_attributes(get(attrs, :node, Dict())),
    edge_attrs = Graphviz.as_attributes(get(attrs, :edge, Dict())),
  )
end

# Graphs
########

""" Convert a graph to a Graphviz graph.

A simple default style is applied. For more control over the visual appearance,
first convert the graph to a property graph, define the Graphviz attributes as
needed, and finally convert the property graph to a Graphviz graph.
"""
function to_graphviz(g::HasGraph; kw...)
  to_graphviz(to_graphviz_property_graph(g; kw...))
end

""" Convert graph or other structure to a property graph suitable for Graphviz.

This function is an intermediate step in many methods of the generic function
[`to_graphviz`](@ref), but can be useful in its own right for customizing the
Graphviz graph beyond whatever options are supported by [`to_graphviz`](@ref).
"""
function to_graphviz_property_graph(g::AbstractGraph;
    prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Union{Symbol,Bool}=false, edge_labels::Union{Symbol,Bool}=false)
  PropertyGraph{Any}(g, v -> node_label(g, node_labels, v),
                     e -> edge_label(g, edge_labels, e);
    prog = prog,
    graph = merge!(Dict(:rankdir => "LR"), graph_attrs),
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
  )
end

node_label(g, name::Symbol, v::Int) = Dict(:label => string(g[v, name]))
node_label(g, labels::Bool, v::Int) = Dict(:label => labels ? string(v) : "")

edge_label(g, name::Symbol, e::Int) = Dict(:label => string(g[e, name]))
edge_label(g, labels::Bool, e::Int) =
  labels ? Dict(:label => string(e)) : Dict{Symbol,String}()

function default_node_attrs(labels::Union{Symbol,Bool})
  shape = labels isa Symbol ? "ellipse" : (labels ? "circle" : "point")
  Dict(:shape => shape, :width => "0.05", :height => "0.05", :margin => "0")
end

# Symmetric graphs
##################

function to_graphviz_property_graph(g::AbstractSymmetricGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Union{Symbol,Bool}=false, edge_labels::Union{Symbol,Bool}=false)
  SymmetricPropertyGraph{Any}(g, v -> node_label(g, node_labels, v),
                              e -> symmetric_edge_label(g, edge_labels, e);
    prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
end

symmetric_edge_label(g, name::Symbol, e::Int) = edge_label(g, name, e)

function symmetric_edge_label(g, labels::Bool, e::Int)
  if labels
    e′ = inv(g,e)
    Dict(:label => (e == e′ ? string(e) : "($(min(e,e′)),$(max(e,e′)))"))
  else
    Dict{Symbol,String}()
  end
end

# Reflexive graphs
##################

function to_graphviz_property_graph(g::AbstractReflexiveGraph;
    prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Union{Symbol,Bool}=false, edge_labels::Union{Symbol,Bool}=false,
    show_reflexive::Bool=false)
  pg = PropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
  )
  add_vertices!(pg, nv(g))
  for v in vertices(g)
    set_vprops!(pg, v, node_label(g, node_labels, v))
  end
  reflexive_edges = Set(refl(g))
  for e in edges(g)
    is_reflexive = e ∈ reflexive_edges
    if show_reflexive || !is_reflexive
      e′ = add_edge!(pg, src(g,e), tgt(g,e))
      is_reflexive && set_eprop!(pg, e′, :style, "dashed")
      set_eprops!(pg, e′, edge_label(g, edge_labels, e))
    end
  end
  pg
end

# Symmetric reflexive graphs
############################

function to_graphviz_property_graph(g::AbstractSymmetricReflexiveGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Union{Symbol,Bool}=false, edge_labels::Union{Symbol,Bool}=false,
    show_reflexive::Bool=false)
  pg = SymmetricPropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
  add_vertices!(pg, nv(g))
  for v in vertices(g)
    set_vprops!(pg, v, node_label(g, node_labels, v))
  end
  reflexive_edges = Set(refl(g))
  for e in edges(g)
    is_reflexive = e ∈ reflexive_edges
    if (is_reflexive ? show_reflexive : e <= inv(g,e))
      e′ = first(add_edge!(pg, src(g,e), tgt(g,e)))
      is_reflexive && set_eprop!(pg, e′, :style, "dashed")
      set_eprops!(pg, e′, symmetric_edge_label(g, edge_labels, e))
    end
  end
  pg
end

# Half-edge graphs
##################

to_graphviz(g::AbstractHalfEdgeGraph; kw...) =
  to_graphviz(to_graphviz_property_graph(g; kw...))

function to_graphviz_property_graph(g::AbstractHalfEdgeGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Union{Symbol,Bool}=false, edge_labels::Union{Symbol,Bool}=false)
  pg = SymmetricPropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
  for v in vertices(g)
    add_vertex!(pg, label=node_labels ? string(v) : "")
  end
  for e in half_edges(g)
    if e == inv(g,e)
      # Dangling half-edge: add an invisible vertex.
      v = add_vertex!(pg, style="invis", shape="none", label="")
      e′ = first(add_edge!(pg, vertex(g,e), v))
      set_eprops!(pg, e′, symmetric_edge_label(g, edge_labels, e))
    elseif e < inv(g,e)
      # Pair of distict half-edges: add a standard edge.
      e′ = first(add_edge!(pg, vertex(g,e), vertex(g,inv(g,e))))
      set_eprops!(pg, e′, symmetric_edge_label(g, edge_labels, e))
    end
  end
  pg
end

# Subgraphs
###########

to_graphviz(subgraph::Subobject{<:HasGraph}; kw...) =
  to_graphviz(to_graphviz_property_graph(subgraph; kw...))

function to_graphviz_property_graph(
    subgraph::Subobject{<:Union{AbstractGraph,AbstractSymmetricGraph}};
    subgraph_node_attrs::AbstractDict=Dict(:color => "cornflowerblue"),
    subgraph_edge_attrs::AbstractDict=Dict(:color => "cornflowerblue"), kw...)
  pg = to_graphviz_property_graph(ob(subgraph); kw...)
  f = hom(subgraph)
  for v in vertices(dom(f))
    set_vprops!(pg, f[:V](v), subgraph_node_attrs)
  end
  for e in edges(dom(f))
    set_eprops!(pg, f[:E](e), subgraph_edge_attrs)
  end
  pg
end

end

max_forks_repo_name: epatters/CompCat


text:
export SR, sr_none, sr_multiplicative, sr_shift
export sr_diag, sr_cholesky, sr_div, sr_minres, sr_lsq, sr_cg, sr_qlp

@enum SRPreconditionType sr_none=1 sr_shift=2 sr_multiplicative=3
@enum SRAlgorithm sr_diag=1 sr_cholesky=2 sr_pivcholesky=3 sr_div=4 sr_qlp=11 sr_minres=12 sr_lsq=14 sr_cg=15
is_iterative(alg::SRAlgorithm) = Int(alg) >=10

################################################################################
######   Preconditioning algorithm definition (structure holding params)  ######
################################################################################
"""
    SR([use_iterative=true, ϵ=0.001, λ0=100, b=0.95, λmin=1e-4, [precondition_type=sr_shift, algorithm=sr_qlp, precision=1e-4])

Stochastic Reconfiguration preconditioner which corrects the gradient according
to the natural gradient computed as S^-1 ∇C. Using this algorithm will lead to
the computation of the S matrix together with the gradient of the cost function
∇C. To compute the natural gradient S^-1∇C an iterative scheme (Minres-QLP) or
a direct inversion is used.

The linear system x = S^-1 ∇C is by default solved with `minres_qlp` iterative
solver. Alternatively you can use `sr_minres`, `sr_cg` or `sr_lsq` to use respectively
the minres, conjugate gradient and least square solvers from IterativeSolvers.jl.
For small systems you can also solve it by computing the pseudo-inverse (`sr_diag`),
the cholesky factorisation (`sr_cholesky`) the pivoted-cholesky factorisation (`sr_pivcholesky`),
and using the automatic julia solver, usually involving qr decomposition (`sr_div`).
Those non-iterative methods are all from Base.LinearAlgebra.

If `use_iterative=true` the inverse matrix `S^-1` is not computed, and an iterative
MINRES-QLP algorithm is used to compute the product S^-1*F

If `precondition_type=sr_shift` then a diagonal uniform shift is added to S
S --> S+ϵ*identity

If `precondition_type=sr_multiplicative` then a diagonal multiplicative shift is added to S
S --> S + max(λ0*b^n,λmin)*Diagonal(diag(S)) where n is the number of the iteration.
"""
mutable struct SR{T1,T2,T3,TP,TT} <: Algorithm
    sr_diag_shift::T1  # cutoff
    sr_diag_mult::T2
    sr_precision::T3   #cutoff
    precondition_type::SRPreconditionType
    algorithm::SRAlgorithm
    use_fullmat::Bool
    λ0::TP
    b::TP
    λmin::TP
    time_evo::TT
end

function SR(T::Type=STD_REAL_PREC; ϵ=0.001, precision=10e-5,
            precondition_type=sr_shift, full_matrix=false, algorithm=sr_cholesky,
            λ0=100.0, b=0.95, λmin=1e-4, time_evo=false)
    return SR(T(ϵ), T(1.0), T(precision), precondition_type,
              algorithm, full_matrix,
              T(λ0), T(b), T(λmin),
              Val(time_evo))
end

is_iterative(alg::SR) = is_iterative(algorithm(alg))
algorithm(alg::SR) = alg.algorithm
# -------------- Base.show extension for nice printing -------------- #
Base.show(io::IO, mm::MIME"text/plain", sr::SR) = begin
    print(io, "SR{$(typeof(sr.sr_diag_shift)), $(typeof(sr.λ0))}\n")
    print(io, "\tActive preconditioner: $(sr.precondition_type)\n")
    if sr.precondition_type == sr_multiplicative
        print(io, "\t - λ0   : $(sr.λ0)\n")
        print(io, "\t - b    : $(sr.b)\n")
        print(io, "\t - λmin : $(sr.λmin)")
    else
        print(io, "\t - ϵ : $(sr.sr_diag_shift)")
    end
end

function algorithm_cache(alg::SR, prob, net, par_cache)
    if is_iterative(alg)
        _sr_iterative_cache(alg, prob, net, par_cache)
    else
        _sr_direct_cache(alg, prob, net, par_cache)
    end
end


text:
module MicrostructureNoise
lgamma(x) = (logabsgamma(x))[1]
using SpecialFunctions, Distributions

using Statistics

include("microstructure.jl")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Select(col₁, col₂, ..., colₙ)
    Select([col₁, col₂, ..., colₙ])

The transform that selects columns `col₁`, `col₂`, ..., `colₙ`.
"""
struct Select{N} <: Stateless
  cols::NTuple{N,Symbol}
end

Select(cols::NTuple{N,AbstractString}) where {N} =
  Select(Symbol.(cols))

Select(cols::AbstractVector) = Select(Tuple(cols))

Select(cols...) = Select(cols)

isrevertible(::Type{<:Select}) = true

function apply(transform::Select, table)
  # retrieve relevant column names
  allcols = collect(Tables.columnnames(table))
  select  = collect(transform.cols)
  reject  = setdiff(allcols, select)

  # keep track of indices to revert later
  sinds = indexin(select, allcols)
  rinds = indexin(reject, allcols)

  # sort indices to facilitate reinsertion
  sperm  = sortperm(sinds)
  sorted = sortperm(rinds)
  reject = reject[sorted]
  rinds  = rinds[sorted]

  # original columns
  cols = Tables.columns(table)

  # selected columns
  scols = [Tables.getcolumn(cols, name) for name in select]

  # rejected columns
  rcols = [Tables.getcolumn(cols, name) for name in reject]

  # table with selected columns
  𝒯 = (; zip(select, scols)...)
  stable = 𝒯 |> Tables.materializer(table)

  stable, (reject, rcols, sperm, rinds)
end

function revert(::Select, newtable, cache)
  # selected columns
  cols   = Tables.columns(newtable)
  select = Tables.columnnames(newtable)
  scols  = [Tables.getcolumn(cols, name) for name in select]

  # rejected columns
  reject, rcols, sperm, rinds = cache

  # restore rejected columns
  anames = collect(select[sperm])
  acols  = collect(scols[sperm])
  for (i, rind) in enumerate(rinds)
    insert!(anames, rind, reject[i])
    insert!(acols, rind, rcols[i])
  end
  𝒯 = (; zip(anames, acols)...)
  𝒯 |> Tables.materializer(newtable)
end

"""
    Reject(col₁, col₂, ..., colₙ)
    Reject([col₁, col₂, ..., colₙ])

The transform that discards columns `col₁`, `col₂`, ..., `colₙ`.
"""
struct Reject{N} <: Stateless
  cols::NTuple{N,Symbol}
end

Reject(cols::NTuple{N,AbstractString}) where {N} =
  Reject(Symbol.(cols))

Reject(cols::AbstractVector) = Reject(Tuple(cols))

Reject(cols...) = Reject(cols)

isrevertible(::Type{<:Reject}) = true

function apply(transform::Reject, table)
  allcols = Tables.columnnames(table)
  reject  = collect(transform.cols)
  select  = Tuple(setdiff(allcols, reject))
  strans  = Select(select)
  newtable, scache = apply(strans, table)
  newtable, (strans, scache)
end

function revert(::Reject, newtable, cache)
  strans, scache = cache
  revert(strans, newtable, scache)
end


text:
mutable struct Packet
    seq :: Int
    src :: UInt16
    dst :: UInt16
    ttl :: UInt8

    time_start :: Int
    time_done :: Int
    done :: Bool
end

Base.@kwdef mutable struct Node <: AbstractAgent
    id::Int

    # (x, y) coordinates
    pos :: Tuple{Int, Int}

    # role = :relay | :sink
    role :: Symbol

    # channel = 37 | 38 | 39
    channel :: UInt8 = 37

    # state = :scanning | :advertising | :sleeping
    state :: Symbol = :scanning

    # the time of the event's (current state) begging
    event_start :: Int = rand(0:50)

    # whether the device is currently sending data (this removes the need for redundant allocations)
    transmitting :: Bool = false

    # the length of the interval between transmissions on different channels (ms)
    t_interpdu :: UInt = 10

    # the length of the interval between scanning on different channels (ms)
    t_scan_interval :: UInt = 20

    # the length of the interval before advertising the received packet
    t_back_off_delay :: UInt = 50

    # the number of extra retransmissions of the received packet
    n_retx_transmit_count :: UInt = 0

    # the length of the delay between bonus transmissions of the received packet (ms)
    t_retx_transmit_delay :: UInt = 20

    # additional random component of delay between retransmissions (min:max ms)
    rt_retx_random_delay :: UnitRange{Int} = 20:50

    # the current number of extra transmissions left
    n_transmit_left :: UInt = 0

    # buffer containing seq of the received packets
    received_packets :: CircularBuffer{Int} = CircularBuffer{Int}(20)

    # transmitting power
    tx_power :: Float64 = 1e-3

    # the reference to the current holding packet
    packet_seq :: Int = 0

    # own ttl for the current packet
    packet_ttl :: UInt8 = 0

    pocket = CircularBuffer{Tuple{Int, Int, UInt8}}(20)
end

Base.@kwdef mutable struct Source <: AbstractAgent
    id :: Int
    pos :: Tuple{Int, Int}

    role :: Symbol = :source
    state :: Symbol = :advertising

    channel :: UInt8 = 37
    event_start :: UInt = 0
    transmitting :: Bool = false

    t_interpdu :: UInt = 1

    # the number of extra transmissions for a packet originating from the source
    n_og_transmit_count :: UInt = 2

    # the length of the delay between additional transmissions of the original packet (ms)
    t_og_transmit_delay :: UInt = 10

    # the range of the additional random component to delay between additional original transmissions
    rt_og_random_delay :: UnitRange{Int} = 0:20

    # the current number of extra transmissions left
    n_transmit_left :: UInt = 0

    # transmitting power
    tx_power :: Float64 = 1e-3

    packet_seq :: Int = 0
    packet_ttl :: UInt8 = 0
end


text:
using GeometricIntegrators.BasisFunctions
using GeometricIntegrators.CommonFunctions
using Test

include("vandermonde_tests.jl")
include("lagrange_basis_tests.jl")
include("lagrange_polynomials_tests.jl")
include("legendre_basis_tests.jl")


text:
# load all with naive unify
include("config.jl")
include("misc.jl")
include("utils.jl")
include("types.jl")
include("primitives.jl")
include("subst.jl")
include("unifybase.jl")
include("naiveunify.jl")
include("reso.jl")
include("merge.jl")
include("newcore.jl")
include("coreprint.jl")
include("dvcreso.jl")
include("analyzer.jl")
include("repl.jl")
include("parser.jl")

println("... naivelogic is loaded")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
const FASTAI_DATA_RECIPES = Dict{String,Vector{DatasetRecipe}}(
    # Vision datasets
    [name => [ImageFolders()] for name in (
        "imagenette", "imagenette-160", "imagenette-320",
        "imagenette2", "imagenette2-160", "imagenette2-320",
        "imagewoof", "imagewoof-160", "imagewoof-320",
        "imagewoof2", "imagewoof2-160", "imagewoof2-320",
        "cifar10", "cifar100", "caltech_101", "mnist_png",
        "mnist_sample", "CUB_200_2011"
    )]...,
    [name => [ImageFolders(filefilterfn=f -> !(occursin("unsup", f)))]
        for name in ("imagewang-160", "imagewang-320", "imagewang")]...,
    "camvid" => [ImageSegmentationFolders()],
    "camvid_tiny" => [ImageSegmentationFolders()],
    "pascal_2007" => [ImageTableMultiLabel()],
    "mnist_tiny" => [ImageFolders(filefilterfn=f -> !occursin("test", f))],
    "mnist_var_size_tiny" => [ImageFolders(filefilterfn=f -> !occursin("test", f))],

    # Tabular datasets
    "adult_sample" => [
        TableDatasetRecipe(file="adult.csv"),
        TableClassificationRecipe(TableDatasetRecipe(file="adult.csv"), :salary),
        TableRegressionRecipe(TableDatasetRecipe(file="adult.csv"), :age),
    ],
)

"""
    const FASTAI_DATA_REGISTRY

The default `DataRegistry` containing every dataset in
the fastai dataset collection.
"""
const FASTAI_DATA_REGISTRY = DatasetRegistry(
    Dict(d => () -> datasetpath(d) for d in DATASETS),
    FASTAI_DATA_RECIPES,
)


text:
using LRUCache

@testset "SDF" begin
    include("sdf_header.jl")
    fn = "0002.sdf"
    blocks = file_summary(fn)

    @testset "Units" begin
        for (key, block) in pairs(blocks)
            if hasfield(typeof(block), :units)
                @test block.units == units[string(key)]
            end
        end
    end

    unsupported = Symbol.(["cpu/proton", "cpu/electron", "run_info", "cpu_rank"])
    @testset "Data" begin
        open(fn, "r") do f
            @testset "$key" for (key, block) in pairs(blocks)
                if key ∉ unsupported
                    @test read!(f, block) == data[string(key)]
                end
            end
        end
    end

    @testset "Utilities" begin
        @test nameof(blocks[:ex]) == "Electric Field/Ex"
        @test labels(blocks[:grid]) == ("X", "Y", "Z")
        @test labels(blocks[Symbol("grid/electron")]) == ("X", "Y", "Z")
    end

    @testset "Cache" begin
        cache = LRU{Tuple{String,AbstractBlockHeader},AbstractArray}(maxsize=2)
        open(fn, "r") do f
            t1 = 0.
            t2 = 0.
            @testset "store" begin
                ex = @timed cached_read(f, blocks[:ex], cache)
                @test read(f, blocks[:ex]) == ex.value
                t1 = ex.time
                @test read(f, blocks[:ey]) == cached_read(f, blocks[:ey], cache)
            end
            @testset "load" begin
                ex = @timed cached_read(f, blocks[:ex], cache)
                @test read(f, blocks[:ex]) == ex.value
                @test read(f, blocks[:ey]) == cached_read(f, blocks[:ey], cache)
                t2 = ex.time
                # test that the cache was used
                @test t2 < 100t1
                @test length(cache) == 2
            end
        end
    end
end


text:
function getSeasStat(series,ilon,ilat,NY,NpY,msc,stdmsc,n)
    fill!(msc,0.0)
    fill!(stdmsc,0.0)
    fill!(n,0.0)
    for iyear=0:(NY-1)
        for iday = 1:NpY
            v=series[ilon,ilat,iyear*NpY+iday]
            if !isnan(v)
                stdmsc[iday]=stdmsc[iday]+v*v
                msc[iday]=msc[iday]+v
                n[iday]=n[iday]+1
            end
        end
    end
    for iday = 1:NpY
        if n[iday]>0
            stdmsc[iday]=sqrt(max(stdmsc[iday]/n[iday]-msc[iday]*msc[iday]/(n[iday]*n[iday]),zero(eltype(series))))
            msc[iday]=msc[iday]/n[iday]
        else
            msc[iday]=nan(eltype(series))
            stdmsc[iday]=nan(eltype(series))
        end
    end
    msc,stdmsc
end

function smooth_circular!(xin::Vector,xout::Vector,wl::Integer=9)
    mod(wl,2)==1 || error("Window length must be uneven")
    offs = div(wl,2)
    fwl=float64(wl)
    n=length(xin)
    awl=0
    s=0.0
    for i=(n-offs+1):n
        if !isnan(xin[i])
            s=s+xin[i]
            awl=awl+1
        end
    end
    for i=1:(offs+1)
        if !isnan(xin[i])
            s=s+xin[i];
            awl=awl+1
        end
    end
    #Start main loop
    for i=1:offs
        xout[i]=s/awl
        if !isnan(xin[n-offs+i])
            s=s-xin[n-offs+i]
            awl=awl-1
        end
        if !isnan(xin[offs+1+i])
            s=s+xin[offs+1+i]
            awl=awl+1
        end
    end
    for i=(offs+1):(n-offs-1)
        xout[i]=s/awl
        if !isnan(xin[i-offs])
            s=s-xin[i-offs]
            awl=awl-1
        end
        if !isnan(xin[offs+1+i])
            s=s+xin[offs+1+i]
            awl=awl+1
        end
    end
    for i=(n-offs):n
        xout[i]=s/awl
        if !isnan(xin[i-offs])
            s=s-xin[i-offs]
            awl=awl-1
        end
        if !isnan(xin[i-n+offs+1])
            s=s+xin[i-n+offs+1]
            awl=awl+1
        end
    end
end

function anomalies(series::Array,ilon,ilat,msc,stdmsc,smsc,sstdmsc,n,series_anomaly;NpY::Int=46,filt=9)
    # length of time series
    N    = size(series,3)
    @assert mod(N,NpY) == 0
    # divide length of time series by samples/year -> nr years
    NY=ifloor(N/NpY)
    # mean seasonal cycle
    getSeasStat(series,ilon,ilat,NY,NpY,msc,stdmsc,n)
    # smooth MAC
    if filt > 1
        smooth_circular!(msc, smsc, filt)
    else
        copy!(smsc,msc)
    end
    #deseasoanlized time series, difference of time series and mac
    for iyear=0:(NY-1),iday=1:NpY
        series_anomaly[ilon,ilat,iyear*NpY+iday]  = series[ilon,ilat,iyear*NpY+iday] - smsc[iday]
    end
end

function get_anomalies(x,NpY;filt=9)
    anomar = similar(x)
    get_anomalies!(x,NpY,anomar,filt=filt)
    return(anomar)
end

function get_anomalies!(x,NpY;filt=9)
    get_anomalies!(x,NpY,x,filt=filt)
    return(x)
end

function get_anomalies!(x,NpY,anomar;filt=9)
    nlon=size(x,1)
    nlat=size(x,2)
    msc     = Array(Float64,NpY) #Allocate once
    stdmsc  = Array(Float64,NpY) #Allocate once
    smsc    = Array(Float64,NpY) #Allocate once
    sstdmsc = Array(Float64,NpY) #Allocate once
    n       = Array(Int64,NpY)
    for lon in 1:nlon, lat in 1:nlat
        anomalies(x,lon,lat,msc,stdmsc,smsc,sstdmsc,n,anomar,NpY=NpY,filt=filt)
        #println("$lon $lat")
    end
    return(anomar)
end

max_stars_repo_name: CAB-LAB/ExtremeStats

max_issues_repo_name: CAB-LAB/ExtremeStats

max_forks_repo_name: CAB-LAB/ExtremeStats


text:
@testset "BOX algorithm" begin

    # one-dimensional
    prob, tspan = exponential_1d()
    sol = solve(prob, tspan=tspan, BOX(δ=0.01))
    @test isa(sol.alg, BOX)
    @test setrep(sol) <: Hyperrectangle
    @test setrep(sol) == Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}
    @test dim(sol) == 1

    # higher-dimensional homogeneous
    prob, tspan = linear5D_homog()
    sol = solve(prob, tspan=tspan, BOX(δ=0.01))
    @test setrep(sol) == Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}
    @test dim(sol) == 5

    # static option
    sol = solve(prob, tspan=tspan, BOX(δ=0.01), static=false) # TODO add static = true
    #HS = Hyperrectangle{Float64,StaticArrays.SArray{Tuple{5},Float64,1,5},StaticArrays.SArray{Tuple{5},Float64,1,5}} # TODO add
    #@test setrep(sol) == HS # TODO add

    # TODO higher-dimensional with input
    #prob, tspan = linear5D()
    #sol = solve(prob, tspan=tspan, BOX(δ=0.01))
end


text:
using Libdl

function _find_gcclibdir()

  gcc_bin = Sys.which("gcc")

  if  "/usr/bin" == gcc_bin[1:8]
    return _find_gcclibdir_in_system()
  else
    return _find_gcclibdir_custom()
  end

end

function _find_gcclibdir_in_system()
  for v in [5,6,7,8,9]
    gcclibdir = "/usr/lib/gcc/x86_64-linux-gnu/$v"
    if isdir(gcclibdir)
      l = Libdl.find_library(joinpath(gcclibdir,"libgomp"))
      if l != ""
        return gcclibdir
      end
    end
  end
  s = """
  libgomp could not be found in system. Try
  \$ sudo apt-get install libgomp1
  """
  error(s)
  ""
end

function _find_gcclibdir_custom()
  gcc_bin = Sys.which("gcc")
  gcc_root = gcc_bin[1:(end-8)]
  gcclibdir = joinpath(gcc_root,"lib64")
  if ! isdir(gcclibdir)
    error("lib64 directory not gound in gcc instalation: $gcclibdir")
  end
  gcclibdir
end

deps_jl = "deps.jl"

if isfile(deps_jl)
  rm(deps_jl)
end

if !haskey(ENV,"MKLROOT")
  s = """
  Environment variable MKLROOT not found.

  Please install intel mkl math lbrary and rebuild.
  """
  error(s)
end

mklroot = ENV["MKLROOT"]

@info "MKLROOT found at: $mklroot"

mkllibdir = joinpath(mklroot,"lib/intel64")

if ! isdir(mkllibdir)
  s = """
  MKL lib directory not found: $mkllibdir
  """
  error(s)
end

@info "MKL libraries found at: $mkllibdir"

gcclibdir = _find_gcclibdir()

if ! isdir(gcclibdir)
  s = """
  gcc lib directory not found: $gcclibdir
  """
  error(s)
end

@info "gcc libraries found at: $gcclibdir"

open(deps_jl,"w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const mklroot = $mklroot))
  println(f, :(const mkllibdir = $mkllibdir))
  println(f, :(const gcclibdir = $gcclibdir))
end


text:
module StateVectorM

export StateVector, GateIndices
export init_state!, getindex, setindex!, getIndicesAfterExclusion, generateBitPatterns
export applyPauliX!, applyHadamard!, applyRX!, applyRY!

import Base.getindex
import Base.setindex!

mutable struct StateVector
    num_qubits::Int64
    data::Vector{ComplexF64}
end

Base.@ccallable function create_StateVector(num_qubits::Int64)::StateVector
    return StateVector(num_qubits, zeros(2^num_qubits).+0im)
end

mutable struct StateVectorPtr
    num_qubits::Int64
    data::Ptr{ComplexF64}
end
mutable struct StateVectorRef
    num_qubits::Int64
    data::Ref{ComplexF64}
end

function init_state!(sv::StateVector)
    sv.data = zeros(2^sv.num_qubits).+0im
    sv[1] = 1 + 0im
    ;
end

struct GateIndices
    internal::Vector{Int64}
    external::Vector{Int64}
    GateIndices(wires, num_qubits) = new(generateBitPatterns(wires, num_qubits), generateBitPatterns(getIndicesAfterExclusion(wires, num_qubits), num_qubits))
end

function getindex(sv::StateVector, index::Int64)
    return getindex(sv.data, index)
end

function setindex!(sv::StateVector, value::ComplexF64, index::Int64)
    sv.data[index] = value
end

Base.@ccallable function getIndicesAfterExclusion(indicesToExclude::Vector{Int64}, num_qubits::Int64)::Vector{Int64}
    indices = Set{Int64}(1:1:num_qubits)
    for excl_idx = indicesToExclude
        pop!(indices, excl_idx);
    end
    collect(indices)
end

Base.@ccallable function generateBitPatterns(qubitIndices::Vector{Int64}, num_qubits::Int64)::Vector{Int64}
    indices = zeros(1)
    for index_it = Iterators.reverse(qubitIndices)
        value = exp2(num_qubits - index_it)
        current_size = length(indices)
        for i = 1:1:current_size
            push!(indices, indices[i] + value)
        end
    end

    return indices
end

Base.@ccallable function applyPauliX!(sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    inverse::Bool = false)::Cvoid
    for ext_idx = externalIndices
        tmp1 = sv[1 + ext_idx + indices[1]]
        sv[1 + ext_idx + indices[1]] = sv[1 + ext_idx + indices[2]]
        sv[1 + ext_idx + indices[2]] = tmp1
    end
    ;
end

Base.@ccallable function applyPauliX!(sv::StateVector, wire::Int64, inverse::Bool = false)::Cvoid
    gi = GateIndices([wire], sv.num_qubits);
    applyPauliX!(sv::StateVector, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyHadamard!(sv::StateVector,
                        indices::Vector{Int64},
                        externalIndices::Vector{Int64},
                        inverse::Bool = false)::Cvoid
    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = (1/sqrt(2))*(v1 + v2)
        sv[1 + ext_idx + indices[2]] = (1/sqrt(2))*(v1 - v2)
    end
    ;
end

Base.@ccallable function applyHadamard!(sv::StateVector, wire::Int64, inverse::Bool = false)::Cvoid
    gi = GateIndices([wire], sv.num_qubits);
    applyHadamard!(sv::StateVector, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyRX!(sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    param::Float64, inverse::Bool = false)::Cvoid
    c = cos(param / 2);
    js = (inverse == true) ?  -1im*sin(-param / 2) : 1im*sin(-param / 2);

    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = c * v1 + js * v2;
        sv[1 + ext_idx + indices[2]] = js * v1 + c * v2;
    end
    ;
end

Base.@ccallable function applyRX!(sv::StateVector, wire::Int64, inverse::Bool = false)::Cvoid
    gi = GateIndices([wire], sv.num_qubits);
    applyRX!(sv::StateVector, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyRY!(  sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    param::Float64, inverse::Bool = false)::Cvoid
    c = cos(param / 2);
    s = (inverse == true) ? -sin(param / 2) : sin(param / 2);

    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = c * v1 - s * v2;
        sv[1 + ext_idx + indices[2]] = s * v1 + c * v2;
    end
    ;
end

Base.@ccallable function applyRY!(sv::StateVector, wire::Int64, inverse::Bool)::Cvoid
    gi = GateIndices([wire], sv.num_qubits);
    applyRY!(sv::StateVector, gi.internal, gi.external, inverse);
end

end

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
# ---
# title: 241. Different Ways to Add Parentheses
# id: problem241
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Divide and Conquer
# link: <https://leetcode.com/problems/different-ways-to-add-parentheses/description/>
# hidden: true
# ---
#
# Given a string of numbers and operators, return all possible results from
# computing all the different possible ways to group numbers and operators. The
# valid operators are `+`, `-` and `*`.
#
# **Example 1:**
#
#
#
#     Input: "2-1-1"
#     Output: [0, 2]
#     Explanation:
#     ((2-1)-1) = 0
#     (2-(1-1)) = 2
#
# **Example 2:**
#
#
#
#     Input:"2*3-4*5"
#     Output: [-34, -14, -10, -10, 10]
#     Explanation: (2*(3-(4*5))) = -34
#     ((2*3)-(4*5)) = -14
#     ((2*(3-4))*5) = -10
#     (2*((3-4)*5)) = -10
#     (((2*3)-4)*5) = 10 ****
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CareData
    using StatsBase
    using Statistics

    using FileIO: load
    using DataFrames
    using Query
    using IterableTables
    using CSVFiles
    using CSV
    using Logging

    using IteratorInterfaceExtensions
    import TableTraits: isiterabletable

    using GlobalDecls
    using ONSCodes
    using Utils
    import LAModelData: get_staying_put_rates, get_18s_level_from_doe

    export load_all
    export get_yp, make_yp_frame, make_carer_frame
    export Carer, YP
    export create_base_datasets, add_carer_to_frame!
    export DataSettings, default_data_settings
    export CPIINDEX, AFC_SURVEY_YEAR, THIS_YEAR, uprate
    export CarerOutcomes, addcareroutcomestoframe!, makecareroutcomesframe
    export fillcomponents!

    THIS_YEAR=2020
    AFC_SURVEY_YEAR=2018

    MINIMUM_WAGE_2019 = 6.15*40 # 18-20
 	JSA_IS_2019 = 57.90

    MINIMUM_WAGE_2020 = 6.45*40 # 18-20
 	JSA_IS_2020 = 57.90 # FIXME unchanged 2019/20 rate

    MINIMUM_WAGE=MINIMUM_WAGE_2020
 	JSA_IS=JSA_IS_2020

    @enum EmploymentStatus TrainingOrEmployment HigherEd OtherEd NEET Unknown
    @enum DataPublisher OFSTED DFE AFC

    export EmploymentStatus, TrainingOrEmployment, HigherEd, OtherEd, NEET, Unknown
    export DataPublisher, OFSTED, DFE, AFC

    function efroms( s :: AbstractString ) :: EmploymentStatus
        t :: EmploymentStatus = Unknown
        if s == "TrainingOrEmployment"
            t = TrainingOrEmployment
        elseif s == "HigherEd"
            t = HigherEd
        elseif s == "OtherEd"
            t = OtherEd
        elseif s == "NEET"
            t = NEET
        else
            t = Unknown
        end
        t
    end

    function efroms( s :: EmploymentStatus)::EmploymentStatus
        s
    end

    #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 whose activity is training or employment
    # from DFe CareLeavers_Act_InTouch19to212018_amended.csv
    CL_All_19to21_2018  = 28_510 # all
    ED_PCS_2018 = (
        #  Care leavers in the year ending 31 March 2018 aged 19 to 21 who were looked after for a total of at least 13 weeks after their 14th birthday including some time after their 16th birthday [1]
        :CL_Act_TE19to21_pc => 25, # 25

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 whose activity is higher education i.e. studies beyond A level
        :CL_Act_HE19to21_pc => 6, # HE 6

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 whose activity is education other than higher education
        :CL_Act_OE19to21_pc => 20, # FE 20

        # Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 who are not in education, employment or training
        :CL_Act_NEET19to21_pc => 39, # NEET 39

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 for whom local authority does not have activity information
        :CL_Act_NoInf19to21_pc => 10 )
         #No inf 10
    # 2019 version
    CL_All_19to21_2019  = 29_930 # all
    ED_PCS_2019 = (
        #  Care leavers in the year ending 31 March 2018 aged 19 to 21 who were looked after for a total of at least 13 weeks after their 14th birthday including some time after their 16th birthday [1]
        :CL_Act_TE19to21_pc => 25.0, # 25

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 whose activity is higher education i.e. studies beyond A level
        :CL_Act_HE19to21_pc => 6.0, # HE 6

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 whose activity is education other than higher education
        :CL_Act_OE19to21_pc => 21.0, # FE 20

        # Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 who are not in education, employment or training
        :CL_Act_NEET19to21_pc => 39.0, # NEET 39

        #  Percentage of care leavers in the year ending 31 March 2018 aged 19 to 21 for whom local authority does not have activity information
        :CL_Act_NoInf19to21_pc => 9.0 ) #No inf 10

    ED_PCS = ED_PCS_2019 # FIXME move this to a function
    CL_All_19to21 = CL_All_19to21_2019

    function educstatus( ccode :: AbstractString ) :: EmploymentStatus
        # FIXME should use LA level stuff and a search
        # FIXME 2020 needed
        e ::EmploymentStatus = NEET
        r = rand(1:90)
        if( r <= 25)
            e = TrainingOrEmployment
        elseif( r <= 31 )
            e = HigherEd
        elseif( r <= 51 )
            e = OtherEd
        end
        e
    end

    # from ofsted ss 2019, addn table 1
    # we'll use this is proportions in each skill class
    # Approved for one type of care	20,260
    # Approved for two types of care	12,735
    # Approved for three types of care	4,905
    # Approved for four or more types of care	5,575
    # then (Main Table, col1 4920 new approved - assume level 0 )
    SKILLS_LEVELS_2018 = [5000, 20060.0-5000,12735.0,4905.0,5575.0]
    SKILLS_PROPS_2018 = cumsum(SKILLS_LEVELS_2018)/sum(SKILLS_LEVELS_2018)

    # 5_110 is 8400 approved - 3290 friends and family
    SKILLS_LEVELS_2019 = [5_110, 20_755.0-5000,13_510.0,5_200.0,4_985.0]
    SKILLS_PROPS_2019 = cumsum(SKILLS_LEVELS_2019)/sum(SKILLS_LEVELS_2019)

    SKILLS_LEVELS=SKILLS_LEVELS_2019
    SKILLS_PROPS=SKILLS_PROPS_2019

    function skilllevel()
        r = rand()
        nsk = size( SKILLS_PROPS )[1]
        for i in 1:nsk
            if r <= SKILLS_PROPS[i]
                return i
            end
        end
    end

    # actual and forecast CPI index rebased to 2019=100; from ONS chart 3.18
    # FIXME unchanged for update
    CPIINDEX = Dict(
        2010 => 0.816463685631,
        2011 => 0.843392062033,
        2012 => 0.881038496092,
        2013 => 0.905956397748,
        2014 => 0.929192356761,
        2015 => 0.942768148942,
        2016 => 0.943146193059,
        2017 => 0.949367808928,
        2018 => 0.974840594465,
        2019 => 1.00,
        2020 => 1.020495466904,
        2021 => 1.040612242049,
        2022 => 1.062698654012,
        2023 => 1.084820802953,
        2024 => 1.106897548551,
        2025 => 1.129423568993,
        2026 => 1.152408007287,
        2027 => 1.175860192508,
        2028 => 1.199789643583,
        2029 => 1.224206073155,
        2030 => 1.249119391524,
        2031 => 1.274539710671,
        2032 => 1.300030504884,
        2033 => 1.326031114982

    )

    function uprate( x :: Real, newyear :: Integer, oldyear :: Integer = -1 ) :: Real
        if oldyear == -1
            oldyear = newyear-1
        end
        x * (CPIINDEX[newyear]/CPIINDEX[oldyear])
    end

    # from OFSTED
    total_in_placements = [50600,51315,51805,51805,52005,53040]
    total_reaching_18 = [3160,3335,3435,4025,3430,3435]
    age18_and_remain = [1685, 1750, 1790, 2190, 1570, 1695]

    ny = size( total_in_placements )[1]
    # @info  "size $ny"
    growth = zeros( ny-1 )

    for i in 1:ny-1
        growth[i]=1+(total_in_placements[i+1]-total_in_placements[i])/total_in_placements[i]
    end

    function make_yp_frame( n::Integer )
        DataFrame(
            year          = zeros( Integer, n ),
            ccode         = Vector{AbstractString}(undef,n),
            carer         = zeros( Integer, n ),
            age           = zeros( Integer, n ),
            housing_costs = zeros( n ),
            earnings      = zeros( n ),
            benefits      = zeros( n ),
            employment_status = Vector{EmploymentStatus}(undef, n ))
    end

    function make_carer_frame( n::Integer )
        DataFrame(
            year              = zeros( Integer, n ),
            ccode             = Vector{AbstractString}(undef,n),
            id                = zeros( Integer, n ),
            age               = zeros( Integer, n ),
            housing_costs     = zeros( n ),
            earnings          = zeros( n ),
            benefits          = zeros( n ),
            employment_status = zeros( Integer, n ),
            skill_level       = zeros( Integer, n ))
     end

     mutable struct Carer
        id                :: Integer
        age               :: Integer
        housing_costs     :: Real
        earnings          :: Real
        benefits          :: Real
        employment_status :: Integer
        skill_level       :: Integer
     end

     function carerfromrow( i ) :: Carer
         Carer( i.id, i.age, i.housing_costs, i.earnings, i.benefits, i.employment_status, i.skill_level )
     end

     function ypfromrow( q ) :: YP
        YP( q.age, q.housing_costs, q.earnings, q.benefits, efroms( q.employment_status ))
     end

     function get_carers(
        carer_dataset,
        ccode :: AbstractString,
        year  :: Integer ) :: Array{Carer}
        @assert isiterabletable( carer_dataset ) "data needs to implement IterableTables"
        q = @from i in carer_dataset begin
            @where (i.ccode == ccode) && (i.year == year)
            @select i
            @collect
        end
        lq = length( q )[1]
        # @info  "lq = $lq"
        carers = Array{Carer}( undef, lq )
        j = 0
        for i in q
            j += 1
            carers[j] = Carer( i.id, i.age, i.housing_costs, i.earnings, i.benefits, i.employment_status, i.skill_level )
        end
        carers
     end # get_carers

    function add_carer_to_frame!( yp_dataset, year ::Integer, ccode :: AbstractString, carer :: Carer )
        @assert isiterabletable( yp_dataset ) "data needs to implement IterableTables; is "*typeof(yp_dataset)
        d = [ year, ccode, carer.id, carer.age, carer.housing_costs,
             carer.earnings, carer.benefits, carer.employment_status, carer.skill_level]
        push!( yp_dataset, d )
    end

    mutable struct CarerOutcomes
        payments_from_la      :: Real
        income_recieved       :: Real
        contributions_from_yp :: Real
        housing_cont          :: Real
        other_cont            :: Real
        fee                   :: Real
        allowance             :: Real
    end

    function CarerOutcomes(; fee :: Real, allowance::Real, housing_cont :: Real, other_cont :: Real ) :: CarerOutcomes
        c = CarerOutcomes( 0.0, 0.0, 0.0, housing_cont, other_cont, fee, allowance )
        fillcomponents!( c )
        c
    end

    function fillcomponents!( c :: CarerOutcomes )
        c.contributions_from_yp =
            c.housing_cont + c.other_cont
        c.income_recieved = c.fee + c.allowance
        c.payments_from_la = max( 0.0, c.income_recieved - c.contributions_from_yp)
    end

    import Base.≈

    function ≈(left :: CarerOutcomes, right::CarerOutcomes )::Bool
        ( left.payments_from_la ≈ right.payments_from_la ) &&
        ( left.contributions_from_yp ≈ right.contributions_from_yp ) &&
        ( left.housing_cont ≈ right.housing_cont ) &&
        ( left.income_recieved ≈ right.income_recieved )
    end

    function makecareroutcomesframe( n :: Integer )
        DataFrame(
            sysno         = zeros( Integer, n ),
            iteration     = zeros( Integer, n ),
            year          = zeros( Integer, n ),
            ccode         = Vector{AbstractString}(undef,n),
            rcode         = Vector{AbstractString}(undef,n),
            carer         = zeros( Integer, n ),
            yp_age        = zeros( Integer, n ),
            payments_from_la      = zeros( n ),
            contributions_from_yp = zeros( n ),
            income_recieved       = zeros( n )
            )
    end

    function addcareroutcomestoframe!(
        out_dataset,
        sysno :: Integer,
        iteration :: Integer,
        year  :: Integer,
        ccode :: AbstractString,
        rcode :: AbstractString,
        carerid::Integer,
        yp_age :: Integer,
        outcomes :: CarerOutcomes )
        @assert isiterabletable( out_dataset ) "data needs to implement IterableTables; is "*typeof(out_dataset)
        d = [ sysno, iteration, year, ccode, rcode, carerid, yp_age, outcomes.payments_from_la, outcomes.contributions_from_yp, outcomes.income_recieved]
        push!( out_dataset, d )
    end

    function make_council_frame( n::Integer )
        DataFrame(
            name  = Vector{AbstractString}(undef,n),
            ccode = Vector{AbstractString}(undef,n),
            rcode = Vector{AbstractString}(undef,n))
    end

    mutable struct YP
        age :: Integer
        housing_costs :: Real
        earnings      :: Real
        benefits      :: Real
        employment_status :: EmploymentStatus
    end

    function makelaoutcomesframe( n::Integer )
        DataFrame(
            year      = zeros( Integer, n ),
            iteration = zeros( Integer, n ),
            council   = Vector{AbstractString}(undef,n),
            ccode     = Vector{AbstractString}(undef,n),
            rcode     = Vector{AbstractString}(undef,n),
            incomes   = zeros(n),
            payments  = zeros(n),
            grant_allocations = zeros(n))
    end

    function ageyp!( yp :: YP, year :: Integer )
        yp.age += 1
        yp.housing_costs = uprate(yp.housing_costs, year)
        yp.earnings = uprate(yp.earnings, year)
        yp.benefits = uprate(yp.benefits, year)
    end

    function agecarer!( carer :: Carer, year :: Integer )
        carer.age += 1
        carer.housing_costs = uprate(carer.housing_costs, year)
        carer.earnings = uprate(carer.earnings, year)
        carer.benefits = uprate(carer.benefits, year)
    end

    function get_yp(
        yp_dataset,
        year       :: Integer,
        carer      :: Carer,
        yearmatch  = nothing ) :: Union{YP,Nothing}
        @assert isiterabletable( yp_dataset ) "data needs to implement IterableTables"
        hits=[]
        if yearmatch !== nothing
            hits = (yp_dataset.carer .== carer.id) .& yearmatch
        else
            hits = (yp_dataset.carer .== carer.id) .& (yp_dataset.year .== year)
        end
        n = 0
        yp = nothing
        for q in eachrow(yp_dataset[hits,:])
            yp = YP( q.age, q.housing_costs, q.earnings, q.benefits, q.employment_status )
            n += 1
        end
        @assert n <= 1
        yp
    end

    """
    looks up a random cat A LHA allowance for this LA. (Random since you can have > 1 BHA for each LA)
    if missing returns the mean of all values in CAT A for 2019  LHA Tables 2019/20
    see: https://www.gov.uk/government/publications/local-housing-allowance-lha-rates-applicable-from-april-2019-to-march-2020
    see: http://sticerd.lse.ac.uk/case/ for what we really want to do here
    """
    function get_shared_accommodation_rate( ccode :: AbstractString )::Real
        v = ONSCodes.get_brma_value( ccode )
        if v < 0.0
            v = 69.16
        end
        v
    end

    function get_yp2( yp_dataset, year :: Integer, carer :: Carer ) :: Union{YP,Nothing}
        @assert isiterabletable( yp_dataset ) "data needs to implement IterableTables"
        q = @from i in yp_dataset begin
            @where i.carer == carer.id && i.year == year
            @select i
            @collect
        end
        n = size( q )[1]
        @assert n <= 1
        if n == 1
            q = q[1]  # handle as a dataframe row
            YP( q.age, q.housing_costs, q.earnings, q.benefits, q.employment_status )
        end
    end

    function addyptoframe!( yp_dataset, year ::Integer, ccode :: AbstractString, carerid::Integer, q :: YP )
        @assert isiterabletable( yp_dataset ) "data needs to implement IterableTables"
        d = [ year, ccode, carerid, q.age, q.housing_costs, q.earnings, q.benefits, q.employment_status]
        push!( yp_dataset, d )
    end

    function newperson()::YP
        YP( 18,0.0,0.0, Unknown )
    end

    mutable struct DataSettings
        name      :: AbstractString
        dataset   :: AbstractString
        description :: AbstractString
        agglevel  :: AggLevel
        poolyears :: Bool
        avggrowth :: Real
        prp_reach_18 :: Real
        startyear :: Integer
        endyear   :: Integer
        annual_sp_increment :: Real
        num_iterations :: Integer
        targets   :: Array{AbstractString}
        reaching_18s_source :: DataPublisher
        datayear  :: Integer
    end

    function default_data_settings()::DataSettings
        DataSettings(
            "",
            "",
            "Description of the run",
            local_authority,
            false,
            mean(growth),
            (total_reaching_18./total_in_placements)[ny],
            2017,
            2026,
            0.0,
            1,
            [],
            OFSTED,
            THIS_YEAR )
    end

    function add_aged_data!(
        ;
        yp_dataset,
        carer_dataset,
        pid       :: Integer,
        ccode     :: AbstractString,
        ladata    :: Dict,
        stayingrates :: Vector{<:Real},
        startyear :: Integer
        )
        @assert isiterabletable( yp_dataset ) "yp_dataset needs to implement IterableTables"
        @assert isiterabletable( carer_dataset ) "carer_dataset needs to implement IterableTables"

        employment_status = educstatus( ccode )
        housing = get_shared_accommodation_rate( ccode )*CPIINDEX[startyear]

        benefits = 0.0
        earnings = 0.0
        if employment_status == TrainingOrEmployment
            earnings = MINIMUM_WAGE*CPIINDEX[startyear]
        end
        benefits = 0.0
        if employment_status == NEET
            benefits = JSA_IS*CPIINDEX[startyear]
        end
        yp = YP( 18, housing, earnings, benefits, employment_status )
        earnings = 0.0
        employment_status = 0
        skill_level = skilllevel()
        housing = 0.0
        benefits = 0.0
        carer = Carer( pid, 40, housing, earnings, benefits, employment_status, skill_level )
        thisyear = startyear
        @info  "stayingrates=$stayingrates"
        for p in 1:3 ### FIXME parameterise 3
            r = rand()
            stayr = stayingrates[p]
            if r > stayr
                @info  "exiting; r=$r stayrate $stayr "
                break
            end
            addyptoframe!( yp_dataset, thisyear, ccode, carer.id, yp )
            add_carer_to_frame!( carer_dataset, thisyear, ccode, carer )
            ageyp!( yp, thisyear )
            agecarer!( carer, thisyear )
            thisyear += 1
        end
    end # add_aged_data!

    function include_this( ccode :: AbstractString; settings :: DataSettings ) :: Boolean
        if ONSCodes.is_aggregate( ccode )
            return false;
        end
        if ccode in SKIPLIST
            return false;
        end
        if size(settings.targets)[1] == 0
            return true
        end
        return ccode in settings.targets
    end

    function create_base_datasets(
        ofdata   :: DataFrame,
        settings :: DataSettings )
        if size(settings.targets)[1] > 0
            targets  = ofdata.ccode .|> [x->x in settings.targets]
            t = settings.targets
            @info  "settings.targets = $t"
            @info  "targets = $targets"
            ofdata = ofdata[targets,:]
        end
        @info  "settings $settings"
        ofdata.number_of_children_or_young_people_in_placements_at_31_march =
            Missings.replace( ofdata[!,:number_of_children_or_young_people_in_placements_at_31_march],0)
        ncouncils = size( ofdata )[1]
        councils = make_council_frame( ncouncils )
        carer_data = make_carer_frame(0)
        yp_data = make_yp_frame(0)
        nc = 0
        pid = 0
        r = 0
        for ofdat in eachrow(ofdata)
            r += 1
            ccode = ofdat.ccode
            if include_this( ccode, settings )
                nc += 1
                if nc > 2000 # test break if needed
                    break
                end
                councils[:name][nc] = ofdat[:council]
                councils[:ccode][nc] = ccode
                councils[:rcode][nc] = ofdat[:rcode]
                base_staying_rates = get_staying_put_rates(
                    ccode,
                    settings.agglevel,
                    settings.poolyears )
                total_fostered = ofdat[:number_of_children_or_young_people_in_placements_at_31_march]
                if settings.reaching_18s_source == OFSTED
                    new_reached_18_base = settings.prp_reach_18 * total_fostered
                else
                    new_reached_18_base = get_18s_level_from_doe( ccode )
                    if ismissing( new_reached_18_base )
                        # tf = typeof(total_fostered )
                        @warn  "new_reached_18 missing for ccode $ccode"
                        new_reached_18_base = settings.prp_reach_18 * total_fostered # fallback
                    end
                end
                g = 1.0
                for year in (settings.startyear-2):settings.endyear # FIXME hack
                    y = year-2019 # fixme make 2019 some sort of parameter
                    stayingrates = base_staying_rates .+ ( settings.annual_sp_increment*y )
                    @info  "at year $year stayingrates = $stayingrates base_staying_rates $base_staying_rates"
                    g = settings.avggrowth^y
                    new_reached_18 = round(Integer,new_reached_18_base*g)
                    @info  "adding $new_reached_18 people for council $ccode total fostered $total_fostered y=$y g=$g new_reached_18_base=$new_reached_18_base"
                    for i in 1:new_reached_18
                        pid += 1
                        add_aged_data!(
                            yp_dataset = yp_data,
                            carer_dataset = carer_data,
                            pid=pid,
                            ccode=ccode,
                            ladata = Dict(),
                            stayingrates=stayingrates,
                            startyear=year
                        )
                    end # people
                end # year
            end # a valid council
        end # councils
        year_matches = Dict()
        for year in (settings.startyear-2):settings.endyear+2
            year_matches[year] = yp_data.year .== year
        end
        (
            councils     = councils,
            carer_data   = carer_data,
            yp_data      = yp_data,
            year_matches = year_matches )
    end

    function load_all( year :: Integer )::NamedTuple
        ofdata = CSV.File( DATADIR*"edited/$(year)/OFDATA.csv", missingstrings=["x","","-", "c",".", "", "..", "*"],
        types=make_type_block(14:1000)) |> DataFrame
        lcnames = Symbol.(Utils.basiccensor.(string.(names(ofdata))))
        rename!( ofdata, lcnames )
        col = 0
        ofs = size( ofdata )[1]
        for name in lcnames
            col+=1
            if( col >13 )
                for row in 1:ofs
                    if ofdata[row,name] === missing
                        ofdata[row,name] = 0
                    end
                end
            end
        end
        # loadtable( DATADIR*"edited/OFDATA.csv", indexcols=[:ccode] )
        grantdata = load( DATADIR*"edited/$(year)/GRANTS_$(year).csv" ) |> DataFrame
        # loadtable(  DATADIR*"edited/GRANTS_2019.csv", indexcols=[:ccode] )
        lcnames = Symbol.(Utils.basiccensor.(string.(names(grantdata))))
        rename!(grantdata, lcnames )
        (ofdata=ofdata, grantdata=grantdata )
    end # load_all

end


text:
module TestQuadratic

using Test
using TreeParzen

objective(params) = (params[:x] - 3)^2

best = fmin(objective, Dict(:x => HP.Uniform(:x, -5.0, 5.0)), 50)
@test abs(best[:x] - 3) < .25

end
true


text:
export IR
module IR
using MLStyle
using MLFS.HM: HMT
import MLFS.HM
import MLFS: InstResolCtx, ExprImpl
using MLFS.TypedIO: @typedIO
export applyImplicits, applyExplicits
using Setfield: @set

@as_record struct Expr
    ln :: LineNumberNode
    ty :: Union{Nothing, HMT}
    expr :: ExprImpl
end

@typedIO Expr(LineNumberNode, Union{Nothing, HMT}, ExprImpl)

@data Decl begin
    Perform(impl::Expr)
    Assign(sym::Symbol, ty::HMT, impl::Expr)
end

@typedIO Decl = begin
    Perform(Expr)
    Assign(Symbol, HMT, Expr)
end

@data ExprImpl′ <: ExprImpl begin
    ETypeVal(HMT)
    EExt(Any)
    EVar(Symbol)
    ELet(Vector{Decl}, Expr)
    EITE(Expr, Expr, Expr)
    EFun(Symbol, Expr)
    EApp(Expr, Expr)
    ETup(Vector{Expr})
    EInt(Int64, UInt8)
    EFloat(Float64, UInt8)
    EStr(String)
    EChar(Char)
    EBool(Bool)
    EIm(Expr, HMT, InstResolCtx)
end

@typedIO ExprImpl = begin
    ETypeVal(HMT)
    EExt(Any)
    EVar(Symbol)
    ELet(Vector{Decl}, Expr)
    EITE(Expr, Expr, Expr)
    EFun(Symbol, Expr)
    EApp(Expr, Expr)
    ETup(Vector{Expr})
    EInt(Int64, UInt8)
    EFloat(Float64, UInt8)
    EStr(String)
    EChar(Char)
    EBool(Bool)
    EIm(Expr, HMT, InstResolCtx)
end

function applyImplicits(e::ExprImpl, implicits::Vector{<:HMT}, finalty::HMT, ln::LineNumberNode, localImplicits::InstResolCtx)
    for im in implicits
        e = EIm(Expr(ln, nothing, e), im, localImplicits)
    end
    Expr(ln, finalty, e)
end

function applyExplicits(e::ExprImpl, explicits::Vector{Expr}, targety::HMT, ln::LineNumberNode)
    for explicit in explicits
        exp = Expr(ln, nothing, e)
        e = EApp(exp, explicit)
    end
    Expr(ln, targety, e)
end

function gTrans(self::Function, root::Decl)
    @match root begin
        Perform(impl) => Perform(self(impl))
        Assign(sym, t, impl) => Assign(sym, t, self(impl))
    end
end

function gTrans(self::Function, root::Expr)
    self(root.expr)
end

function gTrans(self::Function, root::ExprImpl)
    @match root begin
        ETypeVal() || EExt() || EVar() || EChar() ||
        EStr() || EBool() || EFloat() || EInt() =>
           root

        ELet(decls, expr) =>
            ELet(IR.Decl[self(decl) for decl in decls], self(expr))
        EITE(a1, a2, a3) =>
            EITE(self(a1), self(a2), self(a3))
        EFun(s, e) => EFun(s, self(e))
        EApp(f, a) => EApp(self(f), self(a))
        ETup(xs) => ETup(IR.Expr[self(x) for x in xs])
        EIm(expr, t, insts) => EIm(self(expr), t, insts)
    end
end

function gTransCtx(self::Function, ctx::Ctx, root::Decl) where Ctx
    !(root) = self(ctx, root)
    @match root begin
        Perform(impl) => Perform(!impl)
        Assign(sym, t, impl) => Assign(sym, t, !impl)
    end

end

function gTransCtx(self::Function, ctx::Ctx, root::Expr) where Ctx
    self(ctx, root.expr)
end

function gTransCtx(self::Function, ctx::Ctx, root::ExprImpl) where Ctx
    !(root) = self(ctx, root)
    @match root begin
        ETypeVal() || EExt() || EVar() || EChar() ||
        EStr() || EBool() || EFloat() || EInt() =>
           root

        ELet(decls, expr) =>
            ELet(IR.Decl[!decl for decl in decls], !expr)
        EITE(a1, a2, a3) =>
            EITE(!a1, !a2, !a3)
        EFun(s, e) => EFun(s, !e)
        EApp(f, a) => EApp(!f, !a)
        ETup(xs) => ETup(IR.Expr[!x for x in xs])
        EIm(expr, t, insts) => EIm(!expr, t, insts)
    end
end

end # module


text:
export lsomp
# Based on the algorithm developed by
# Baptista, R., Stolbunov, V., & Nair, P. B. (2019).
# Some greedy algorithms for sparse polynomial chaos expansions.
# Journal of Computational Physics, 387, 303-325.

# Solve ψc = u using the orthogonal matching pursuit with least-square

function lsomp(ψ::AbstractMatrix{T}, u::AbstractVector{T}; invert::Bool=true, verbose::Bool = true, ϵrel::Float64 = 1e-1, maxterms::Int64=typemax(Int64)) where {T}
    m, n = size(ψ)

    residue = copy(u)
    idxset = Int64[]
    dict = collect(1:n)
    ϵu = norm(u)
    ϵhist = [ϵu]
    ϵrel *= ϵu
    ϕ = zeros(T, n)

    for j=1:n
        for i=1:m
        ϕ[j] += ψ[i,j]^2
        end
    end

    ek = zeros(T, n)
    q = zeros(T, m)
    new_idx = 1

    @inbounds for k=1:n
        # Step 6 & 7 update denominators and choose best candidate

        entry = 0.0
        for (j, idx) in enumerate(dict)
            ψj = view(ψ,:,idx)

            ratio = residue'*ψj
            ratio = ratio^2
            ratio /= ϕ[idx]

            if ratio > entry
                entry = copy(ratio)
                new_idx = copy(idx)
            end
        end

        # Step 8 Update set of selected basis
        push!(idxset, new_idx)

        # Step 9 Update candidate dictionary
        filter!(x-> x != new_idx, dict)
        # setdiff!(dict, new_idx)

        # Step 10 Compute residual
        c = view(ψ,:,idxset)\u
        residue = u - view(ψ,:,idxset)*c

        # Step 12 Calculate stopping critera
        ϵ = norm(residue)
        if verbose == true
            push!(ϵhist, copy(ϵ))
        end
        if ϵ < ϵrel || k == maxterms
            break
        end
    end

    # Solve the system with the set of indices
    if verbose == true
        if invert == true
            c = zeros(length(idxset))
            c .= view(ψ,:, idxset)\u
            return idxset, c, ϵhist
        else
            return idxset, ϵhist
        end
    else
        if invert == true
            c = zeros(length(idxset))
            c .= view(ψ,:, idxset)\u
            return idxset, c
        else
            return idxset
        end
    end
end

max_forks_repo_name: mleprovost/QROMP


text:
using LeapSeconds

export insideleap, getleap

insideleap(ep::Epoch{S}) where {S} = false
getleap(ep::Epoch{S}) where {S} = 0.0

struct TAIOffset
    date::TAIEpoch{Float64}
    jd::Float64
    start::TAIEpoch{Float64}
    reference::TAIEpoch{Float64}
    leap::Float64
    offset::Float64
    slope_utc::Float64
    slope_tai::Float64
end

@inline function getoffset(t::TAIOffset, ep::Epoch)
    t.slope_tai == 0.0 && return t.offset

    t.offset + value(ep - t.reference) * t.slope_tai
end

@inline function getoffset(to::TAIOffset, d::Date, t::Time)
    days = AstroDates.julian(d) - to.jd
    fraction = secondinday(t)
    to.offset + days * (to.slope_utc * SECONDS_PER_DAY) + fraction * to.slope_utc
end

Base.isless(t::TAIOffset, ep::Epoch) = isless(ep, t.date)
Base.isless(ep::Epoch, t::TAIOffset) = isless(ep, t.date)

Base.isless(t::TAIOffset, jd::Float64) = isless(jd, t.jd)
Base.isless(jd::Float64, t::TAIOffset) = isless(jd, t.jd)

const TAI_OFFSETS = TAIOffset[]

@inline function findoffset(ep)
    idx = searchsortedlast(TAI_OFFSETS, ep)
    idx == 0 && return nothing

    TAI_OFFSETS[idx]
end

const EPOCHS = [LeapSeconds.EPOCHS; LeapSeconds.LS_EPOCHS]
const OFFSETS = [LeapSeconds.OFFSETS; LeapSeconds.LEAP_SECONDS]
const DRIFT_EPOCHS = [LeapSeconds.DRIFT_EPOCHS;
                      zeros(length(LeapSeconds.LS_EPOCHS))]
const DRIFT_RATES = [LeapSeconds.DRIFT_RATES;
                     zeros(length(LeapSeconds.LS_EPOCHS))]

for (ep, offset, dep, rate) in zip(EPOCHS, OFFSETS, DRIFT_EPOCHS, DRIFT_RATES)
    tai = TAIEpoch(ep * days, origin=:modified_julian)
    previous = isempty(TAI_OFFSETS) ? 0.0 : getoffset(last(TAI_OFFSETS), date(tai), AstroDates.H00)
    ref = TAIEpoch(TAIEpoch(dep * days, origin=:modified_julian), offset)
    start = TAIEpoch(tai, previous)
    start_offset = offset + (ep - dep) * rate
    stop = TAIEpoch(tai, start_offset)
    slope = rate / SECONDS_PER_DAY
    leap = value(stop - start) / (1 + slope)
    o = TAIOffset(start, ep + LeapSeconds.MJD_EPOCH,
                  TAIEpoch(start, leap),
                  ref,
                  leap, offset, slope, slope / (1 + slope))
    push!(TAI_OFFSETS, o)
end

@inline function insideleap(ep::UTCEpoch)
    ep = TAIEpoch(ep)
    offset = findoffset(ep)
    offset === nothing && return false

    ep < offset.start
end

@inline function getleap(ep::UTCEpoch)
    ep = TAIEpoch(ep)
    offset = findoffset(ep)
    offset === nothing && return 0.0

    offset.leap
end


text:
using Test
using Dash
using Dash: HandlerState, main_registry, start_reload_poll, enable_dev_tools!
import HTTP
import JSON3
@testset "reload state" begin
    rm("assets/test2.png", force = true)
    app = dash()
    enable_dev_tools!(app, dev_tools_hot_reload_watch_interval = 1.)
    app.layout = html_div()
    state = HandlerState(app, main_registry())
    start_reload_poll(state)
    initial_hash = state.reload.hash
    sleep(1)
    @test length(state.reload.changed_assets) == 0
    write("assets/test2.png", "")
    sleep(2)
    @test state.reload.hash != initial_hash
    @test length(state.reload.changed_assets) == 1
    @test length(state.reload.changed_assets) == 1
    @test state.reload.changed_assets[1].url == "/assets/test2.png"
    rm("assets/test2.png", force = true)

end
@testset "reload handler" begin
    rm("assets/test2.css", force = true)
    app = dash()
    enable_dev_tools!(app, dev_tools_hot_reload = true, dev_tools_hot_reload_watch_interval = 1.)
    app.layout = html_div()
    handler = make_handler(app)
    write("assets/test2.css", "")
    sleep(2)
    response = HTTP.handle(handler, HTTP.Request("GET", "/_reload-hash"))
    data = JSON3.read(response.body)
    @test length(data.files) == 1
    @test data.files[1].url == "/assets/test2.css"
    @test data.files[1].is_css == true
    rm("assets/test2.css", force = true)
end


text:
using MAT, Test

function check(filename, result)
    matfile = matopen(filename)
    for (k, v) in result
        @test exists(matfile, k)
        got = read(matfile, k)
        if !isequal(got, v) || (typeof(got) != typeof(v) && (!isa(got, String) || !(isa(v, String))))
            close(matfile)
            error("""
                Data mismatch reading $k from $filename ($format)

                Got $(typeof(got)):

                $(repr(got))

                Expected $(typeof(v)):

                $(repr(v))
                """)
        end
    end
    @test union!(Set(), names(matfile)) == union!(Set(), keys(result))
    close(matfile)

    mat = matread(filename)
    if !isequal(mat, result)
        error("""
            Data mismatch reading $filename ($format)

            Got:

            $(repr(mat))

            Expected:

            $(repr(result))
            """)
        close(matfile)
        return false
    end

    return true
end
cd(dirname(@__FILE__))
for filename in readdir("v4")
    #println("testing $filename")
    d = matread("v4/$filename")
    matwrite4("v4/tmp.mat", d)
    check("v4/tmp.mat", d)
    rm("v4/tmp.mat")
end


text:
import Base.iterate
export iterate

function iterate(x::Variable, s=0)
    return s >= length(x) ? nothing : (x[s+1], s+1)
end


text:
using LinearAlgebra, Statistics, Plots
gr(fmt=:png);

x_vals = [0 0 0; 2 -3 -4]
y_vals = [0 0 0; 4 3 -3.5]

plot(x_vals, y_vals, arrow = true, color = :blue,
     legend = :none, xlims = (-5, 5), ylims = (-5, 5),
     annotations = [(2.2, 4.4, "[2, 4]"),
                    (-3.3, 3.3, "[-3, 3]"),
                    (-4.4, -3.85, "[-4, -3.5]")],
     xticks = -5:1:5, yticks = -5:1:5,
     framestyle = :origin)

x = ones(3)
y = [2,4,6]

dot(x,y)
norm(x)

# fixed linear function, to generate a plane
f(x, y) = 0.2x + 0.1y

# lines to vectors
x_vec = [0 0; 3 3]
y_vec = [0 0; 4 -4]
z_vec = [0 0; f(3, 4) f(3, -4)]

# draw the plane
n = 20
grid1 = range(-5, 5, length = n)
z2 = [ f(grid1[row], grid1[col]) for row in 1:n, col in 1:n ]
wireframe(grid1, grid1, z2, fill = :blues, gridalpha =1 )
plot!(x_vec, y_vec, z_vec, color = [:blue :green], linewidth = 3, labels = "",
      colorbar = false)

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
using ScTenifold
X0=rand(100,1000);
X1=copy(X0)
X1[4,:].=0.0

d,fc,p,adjp=ScTenifold.sctenifoldnet(X0,X1,donorm=false)

#@show Threads.nthreads()
#@time Z0=ScTenifoldNet.tenrnet(X0, donorm=false)
#@time Z1=ScTenifoldNet.tenrnet(X1, donorm=false)
#Z0=0.5*(Z0+Z0')
#Z1=0.5*(Z1+Z1')
#@time d,aln0,aln1=ScTenifoldNet.manialn(Z0,Z1)
#fc,p,adjp=ScTenifoldNet.drgenes(d)

#using StatsPlots, Distributions
#x=rand(Chisq(1), length(fc))
#qqplot(x, fc)

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
# declare some types -
struct Result{T}
    value::T
end

struct AVKError
    message::String
end

function _is_path_valid(path_to_file::String)::Result{Bool}

    # the config file should be a json file, and should be reachable -
    # TODO: fill me in ...
    return Result{Bool}(true)
end

function _is_string_empty(raw_string::String)::Result{Bool}

    # if we have an empty string - return true
    if isempty(raw_string) == true
        return Result{Bool}(true)
    end

    # default Return is flase -
    return Result{Bool}(false)
end

function _check_missing_api_key(user_model::AVKUserModel)::(Union{T, Nothing} where T<:Any)

    # do we have the alpha_vantage_api_key -
    if (hasfield(AVKUserModel, :alphavantage_api_key) == false)
        # throw -
        return Result{AVKError}(AVKError("user model is missing API key information"))
    end

    # get the key -
    api_key = user_model.alphavantage_api_key

    # check -
    if (isempty(api_key) == true)

        # formulate an error message -
        error_message = "the API key is empty in the user model"

        # throw -
        return Result{AVKError}(AVKError(error_message))
    end

    #
    return nothing
end

function _check_missing_symbol(stock_symbol::String)::(Union{T, Nothing} where T<:Any)

    if (isempty(stock_symbol) == true)

        # formulate an error message -
        error_message = "missing stock symbol"

        # throw -
        return Result{AVKError}(AVKError(error_message))
    end

    # return nothing -
    return nothing
end

function _check_json_api_return_data(api_call_raw_data::String)::(Union{T, Nothing} where T<:Any)

    # well formed JSON?
    if _is_valid_json(api_call_raw_data).value == false
        return AVKError("invalid JSON $(api_call_raw_data)")
    end

    # need to check to see if legit data is coming back from the service -
    api_data_dictionary = JSON.parse(api_call_raw_data)
    if (haskey(api_data_dictionary,"Error Message") == true)

        # grab the error mesage -
        error_message = api_data_dictionary["Error Message"]

        # throw -
        return Result{AVKError}(AVKError(error_message))
    end

    # need to check - are we hitting the API call limit?
    if (haskey(api_data_dictionary,"Note") == true)

        # grab the error mesage -
        error_message = api_data_dictionary["Note"]

        # throw -
        return Result{AVKError}(AVKError(error_message))
    end

    # default -
    return nothing
end

function _check_user_model(user_model::AVKUserModel)::(Union{T, Nothing} where T<:Any)
    return Result{Bool}(true)
end

"""
    _is_valid_json(raw_string)->Bool

Checks to see if the input `string` is a valid JSON structure.
Returns `true` indicating valid JSON, `false` otherwise.
"""
function _is_valid_json(raw_string::String)::Result{Bool}

    # check: do we have an empty string?
    if (_is_string_empty(raw_string).value == true)
        return Result{Bool}(false)
    end

    # otherwise, to check to see if the string is valid JSON, try to
    # parse it.
    try
        JSON.parse(raw_string)
        return Result{Bool}(true)
    catch
        return Result{Bool}(false)
    end
end

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI