text:
"""
    LMPTools

Convenience and utility functions for LongwaveModePropagator.jl.
"""
module LMPTools

using PyCall, Conda
using Rotations, StaticArrays
using HDF5, Dates
using GeographicLib, SatelliteToolbox
using LongwaveModePropagator

project_path(parts...) = normpath(@__DIR__, "..", parts...)

const GROUND_DATA = h5read(project_path("data", "conductivity_data.h5"), "/")

const TRANSMITTER = Dict(
    :NAA => Transmitter("NAA", 44.646394, -67.281069, 24e3),
    :NAU => Transmitter("NAU", 18.39875, -67.177433, 40.75e3),
    :NML => Transmitter("NML", 46.366014, -98.335544, 25.2e3),
    :NLK => Transmitter("NLK", 48.203611, -121.917056, 24.8e3),
    :NPM => Transmitter("NPM", 21.420, -158.151, 21.4e3)
)

const CHAOS = PyNULL()
const CHAOS_MODEL = PyNULL()

export TRANSMITTER
export range
export get_ground, get_groundcode, get_epsilon, get_sigma, groundsegments
export igrf, chaos, load_CHAOS_matfile
export zenithangle, isday, ferguson, flatlinearterminator, smoothterminator, fourierperturbation

function __init__()
    cp = try
        pyimport("chaosmagpy")
    catch
        # Install chaosmagpy
        Conda.add(["numpy", "scipy", "pandas", "cython", "cartopy", "matplotlib", "h5py"])

        Conda.pip_interop(true)
        Conda.pip("install", ["cdflib", "hdf5storage", "lxml"])
        Conda.pip("install", "chaosmagpy")

        pyimport("chaosmagpy")
    end

    chaosfile = newestchaos()
    @info "Loading $chaosfile"

    model = cp.load_CHAOS_matfile(joinpath(project_path("data"), chaosfile))
    copy!(CHAOS, cp)
    copy!(CHAOS_MODEL, model)
end

function newestchaos()
    datafiles = readdir(project_path("data"); sort=true)
    newestfile = ""
    for f in datafiles
        if occursin("CHAOS", f) && f > newestfile
            newestfile = f
        end
    end
    return newestfile
end

"""
    load_CHAOS_matfile(file)

Load a CHAOS model coefficient MATLAB `.mat` file for use by [`chaos`](@ref).

These files can be found at https://www.spacecenter.dk/files/magnetic-models/CHAOS-7/.

By default, LMPTools loads $(newestchaos()).
"""
function load_CHAOS_matfile(file)
    cp = pyimport("chaosmagpy")
    @info "Loading $file"
    model = cp.load_CHAOS_matfile(joinpath(project_path("data"), file))
    copy!(CHAOS, cp)
    copy!(CHAOS_MODEL, model)

    return nothing
end

###
# Extend library functions

function GeographicLib.GeodesicLine(tx::Transmitter, rx::Receiver)
    return GeodesicLine(tx.longitude, tx.latitude; lon2=rx.longitude, lat2=rx.latitude)
end

"""
    range(tx::Transmitter, rx::Receiver)

Return the range (great-circle distance, but on the WGS-84 ellipsoid) in meters between
`tx` and `rx`.
"""
function Base.range(tx::Transmitter, rx::Receiver)
    return inverse(tx.longitude, tx.latitude, rx.longitude, rx.latitude).dist
end

###
# Ground

function searchsortednearest(a, x)
    @assert issorted(a) "`a` must be sorted!"

    idx = searchsortedfirst(a, x)

    idx == 1 && return idx
    idx > length(a) && return length(a)
    a[idx] == x && return idx
    abs(a[idx]-x) < abs(a[idx-1]-x) ? (return idx) : (return idx-1)
end

"""
    get_ground(param, lat, lon)

Return the ground `param` "codemap", "epsilonmap", or "sigmamap" at position `lat`, `lon`
in degrees.

Convenience functions are available for the most common fields with arguments `lat`, `lon`
only:

    * [`get_groundcode`](@ref)
    * [`get_sigma`](@ref)
    * [`get_epsilon`](@ref)
"""
function get_ground(param, lat, lon)
    latidx = searchsortednearest(GROUND_DATA["lat"], lat)
    lonidx = searchsortednearest(GROUND_DATA["lon"], lon)

    return GROUND_DATA[param][latidx, lonidx]
end
get_groundcode(lat, lon) = get_ground("codemap", lat, lon)
get_sigma(lat, lon) = get_ground("sigmamap", lat, lon)
get_epsilon(lat, lon) = get_ground("epsilonmap", lat, lon)

"""
    groundsegments(tx::Transmitter, rx::Receiver; resolution=20e3, require_permittivity_change=false)

Return a `Vector` of `Ground`'s and a `Vector` of distances at which the ground changes
between `tx` and `rx` using path `resolution` in meters.

If `require_permittivity_change == true`, then if the relative permittivity doesn't change
we don't consider it a ground change. This helps reduce the number of short small changes
that occur during a long path which have the same permittivity but a small change in
conductivity that has a small influence on the electric field in the EIWG but comes a the
cost of additional (potentially nearly identical) segments in the waveguide model.
"""
function groundsegments(tx::Transmitter, rx::Receiver; resolution=20e3, require_permittivity_change=false)

    line = GeodesicLine(tx, rx)
    pts = waypoints(line, dist=resolution)

    lat, lon, dist = first(pts).lat, first(pts).lon, first(pts).dist
    lastground = GROUND[get_groundcode(lat, lon)]

    # initialize
    grounds = [lastground]
    distances = [dist]  # == 0.0

    for pt in pts
        lat, lon, dist = pt.lat, pt.lon, pt.dist
        ground = GROUND[get_groundcode(lat, lon)]

        if require_permittivity_change
            if ground.ϵᵣ != lastground.ϵᵣ
                push!(grounds, ground)
                push!(distances, dist)
                lastground = ground
            end
        elseif ground != lastground
            push!(grounds, ground)
            push!(distances, dist)
            lastground = ground
        end
    end

    return grounds, distances
end

###
# Magnetic field

"""
    igrf(geoaz, lat, lon, year; alt=60e3)

Return a `BField` from IGRF-13 for position (`lat`, `lon`)°  at fractional
`year`. `geoaz`° is the geodetic azimuth of the path from transmitter to receiver.
By default, the magnetic field at an altitude of 60,000 meters is returned,
but this can be overridden with the `alt` keyword argument.
"""
function igrf(geoaz, lat, lon, year; alt=60e3)
    n, e, v, t = igrf13syn(0, year, 1, alt/1000, 90-lat, mod(lon, 360))  # in nT!, @60km altitude

    # Rotate the nev frame to the propagation path xyz frame
    # negate az to correct rotation direction for downward pointing v
    R = RotXZ(π, -deg2rad(geoaz))
    Rd = R*SVector(n,e,v)

    return BField(t*1e-9, Rd[1]/t, Rd[2]/t, Rd[3]/t)
end

"""
    igrf(tx::Transmitter, rx::Receiver, year, dists; alt=60e3)

Return a `Vector{BField}` at each distance in `dists` in meters along the path from `tx`
to `rx` in fractional `year`.

If applying `igrf` to a single distance, it is recommended to use the `geoaz`, `lat`, `lon`
form of [`igrf`](@ref).
"""
function igrf(tx::Transmitter, rx::Receiver, year, dists; alt=60e3)
    line = GeodesicLine(tx, rx)
    geoaz = inverse(tx.longitude, tx.latitude, rx.longitude, rx.latitude).azi

    bfields = Vector{BField}(undef, length(dists))
    for (i, d) in enumerate(dists)
        lo, la, _ = forward(line, d)
        bfields[i] = igrf(geoaz, la, lo, year; alt=alt)
    end
    return bfields
end

"""
    chaos(geoaz, lat::Real, lon:Real, year; alt=60e3)

Return a `BField` from CHAOS-7 using internal and external sources for position
(`lat`, `lon`)°  at fractional `year`.

`geoaz`° is the geodetic azimuth of the path from transmitter to receiver.
By default, the magnetic field at an altitude of 60,000 meters is returned,
but this can be overridden with the `alt` keyword argument.

# References

[^1]: C. C. Finlay et al., “The CHAOS-7 geomagnetic field model and observed changes in the
    South Atlantic Anomaly,” Earth Planets Space, vol. 72, no. 1, Art. no. 1, Dec. 2020,
    doi: 10.1186/s40623-020-01252-9.
"""
function chaos(geoaz, lat::Number, lon::Number, year; alt=60e3)
    Re = 6371.2
    r = Re + alt/1000
    t = CHAOS.data_utils.dyear_to_mjd(year)  # MJD2000

    theta = 90 - lat  # geocentric co-lat (deg)
    phi = lon  # geocentric lon (deg)

    Br, Bt, Bp = CHAOS_MODEL(t, r, theta, phi)
    u, s, e = only(Br), only(Bt), only(Bp)

    # Rotate the "use" frame to the propagation path xyz frame
    # negate az to correct rotation direction for downward pointing v
    R = RotXZ(π, -deg2rad(geoaz))
    Rd = R*SVector(-s,e,-u)

    mag = hypot(u, s, e)
    return BField(mag*1e-9, Rd[1]/mag, Rd[2]/mag, Rd[3]/mag)
end

"""
    chaos(tx::Transmitter, rx::Receiver, year, dists; alt=60e3)

Return a `Vector{BField}` at each distance in `dists` in meters along the path from `tx`
to `rx` in fractional `year`.

If applying `chaos` to a single distance, it is recommended to use the `geoaz`, `lat`, `lon`
form of [`chaos`](@ref).
"""
function chaos(tx::Transmitter, rx::Receiver, year, dists; alt=60e3)
    line = GeodesicLine(tx, rx)
    geoaz = inverse(tx.longitude, tx.latitude, rx.longitude, rx.latitude).azi

    lats = Vector{Float64}(undef, length(dists))
    lons = similar(lats)
    for (i, d) in enumerate(dists)
        lo, la, _ = forward(line, d)
        lats[i] = la
        lons[i] = lo
    end

    bfields = chaos(geoaz, lats, lons, year; alt)

    return bfields
end

"""
    chaos(geoaz, lats, lons, year; alt=60e3)

Return a `Vector{BField}` at each latitude-longitude pair formed from each element of `lats`
and `lons` in fractional `year`.
"""
function chaos(geoaz, lats, lons, year; alt=60e3)
    length(lats) == length(lons) || throw(ArgumentError("length(lats) must equal length(lons)"))

    Re = 6371.2
    r = Re + alt/1000
    t = CHAOS.data_utils.dyear_to_mjd(year)  # MJD2000

    theta = 90 .- lats  # geocentric co-lat (deg)
    phi = lons  # geocentric lon (deg)

    u, s, e = CHAOS_MODEL(t, r, theta, phi)

    # Rotate the "use" frame to the propagation path xyz frame
    # negate az to correct rotation direction for downward pointing v
    R = RotXZ(π, -deg2rad(geoaz))

    bfields = Vector{BField}(undef, length(lats))
    for i in eachindex(bfields)
        Rd = R*SVector(-s[i],e[i],-u[i])

        mag = hypot(u[i], s[i], e[i])
        bfields[i] = BField(mag*1e-9, Rd[1]/mag, Rd[2]/mag, Rd[3]/mag)
    end

    return bfields
end

###
# Sun position

"""
    zenithangle(lat, lon, y, m, d, h, Δτ=nothing)

Return the solar zenith angle in degrees for year `y`, month `m`, day `d`, and hour `h` in
universal time (UT1). `Δτ` is the difference between terrestrial time and universal time.
The position in `lat`, `lon` is in degrees.

If `Δτ` is `nothing`, it will be computed by a linear extrapolation. Although this is not
accurate because `Δτ` is irregular, an error up to 30 seconds does not significantly affect
the result.

This function uses the calculation-optimized formulation of Algorithm 5 from

    Grena, R., “Five new algorithms for the computation of sun position from 2010 to 2110,”
    Solar Energy, 86, 2012, pp. 1323-1337. doi: 10.1016/j.solener.2012.01.024

which has a maximum error of 0.0027° (~10 arcsec) between year 2010 and 2110.

Refraction is not considered in this implementation.
"""
function zenithangle(lat, lon, yr::Int, mo::Int, dy::Int, hr::Real, Δτ=nothing)
    θ, ϕ = deg2rad(lon), deg2rad(lat)

    # Time scale computation
    if mo <= 2
        mo += 12
        yr -= 1
    end

    t = trunc(Int, 365.25*(yr - 2000)) + trunc(Int, 30.6001*(mo + 1)) -
        trunc(Int, 0.01*yr) + dy + 0.0416667*hr - 21958

    if isnothing(Δτ)
        Δτ = 96.4 + 0.00158*t
    end

    te = t + 1.1574e-5*Δτ

    # Heliocentric ecliptic longitude L
    ωa = 0.0172019715
    a = SVector(3.33024e-2, 3.512e-4, 5.2e-6)
    b = SVector(-2.0582e-3, -4.07e-5, -9e-7)
    s1, c1 = sincos(ωa*te)
    s2, c2 = 2*s1*c1, (c1 + s1)*(c1 - s1)
    s3, c3 = s2*c1 + c2*s1, c2*c1 - s2*s1
    s, c = SVector(s1, s2, s2), SVector(c1, c2, c3)

    β = 2.92e-5
    dβ = -8.23e-5
    ω = SVector(1.49e-3, 4.31e-3, 1.076e-2, 1.575e-2, 2.152e-2, 3.152e-2, 2.1277e-1)
    d = SVector(1.27e-5, 1.21e-5, 2.33e-5, 3.49e-5, 2.67e-5, 1.28e-5, 3.14e-5)
    φ = SVector(-2.337, 3.065, -1.533, -2.358, 0.074, 1.547, -0.488)

    L = 1.7527901 + 1.7202792159e-2*te
    for k = 1:3
        L += a[k]*s[k] + b[k]*c[k]
    end
    L += dβ*s1*sin(β*te)
    for i = 1:7
        L += d[i]*sin(ω[i]*te + φ[i])
    end

    # Nutation correction
    ωn = 9.282e-4
    ν = ωn*te - 0.8
    sν, cν = sincos(ν)
    Δλ = 8.34e-5*sν
    λ = L + π + Δλ
    ϵ = 4.089567e-1 - 6.19e-9*te + 4.46e-5*cν
    sλ, cλ = sincos(λ)
    sϵ, cϵ = sincos(ϵ)

    α = atan(sλ*cϵ, cλ)
    α < 0 && (α += 2π)  # [-π/2, π/2] → [0, 2π]

    δ = asin(sλ*sϵ)
    H = 1.7528311 + 6.300388099*t + θ - α + 0.92*Δλ
    H = mod2pi(H + π) - π  # hour angle

    # Calculate zenith and azimuth angle with parallax correction for topocentric coords
    sϕ, cϕ = sincos(ϕ)
    sδ, cδ = sincos(δ)
    cH = cos(H)
    # sH, cH = sincos(H)

    se0 = sϕ*sδ + cϕ*cδ*cH
    ep = asin(se0) - 4.26e-5*sqrt(1 - se0^2)

    # Γ = atan(sH, cH*sϕ - sδ*cϕ/cδ)  # azimuth

    # Refraction correction to elevation
    # Δre = 0.08422*P/(273 + T)/tan(ep + 0.003138/(ep + 0.08919))

    z = π/2 - ep  # - Δre  # zenith angle

    return rad2deg(z)
end

"""
    zenithangle(lat, lon, dt::DateTime, Δτ=nothing)

Return the solar zenith angle in degrees for universal `dt` with `Δτ` difference in seconds
between terrestrial time and universal time. The position `lat`, `lon` should be in degrees.
"""
function zenithangle(lat, lon, dt::DateTime, Δτ=nothing)
    y, m, d = yearmonthday(dt)
    h = hour(dt) + minute(dt)/60 + second(dt)/3600 + millisecond(dt)/3600_000

    return zenithangle(lat, lon, y, m, d, h, Δτ)
end

"""
    isday(sza)

Return `true` if solar zenith angle `sza` is less than 98°, otherwise return `false`.
"""
function isday(sza)
    abs(sza) < 98 ? true : false
end

"""
    ferguson(lat, sza, month_num::Int) → (h′, β)
    ferguson(lat, sza, dt::DateTime)

Ferguson ionosphere for geographic latitude `lat` in degrees, solar zenith angle `sza` in
degrees, and month number `month_num` or `DateTime` `dt`.

If the month term of the ionosphere fit is not desired, specify `month_num = 0`.

The higher order terms in [^Ferguson1980] (sunspot number and magnetic absorption) are ignored.

# References

[Ferguson1980]: J. A. Ferguson, "Ionospheric profiles for predicting nighttime VLF/LF
    propagation," NOSC/TR-530, 1980.
"""
function ferguson(lat, sza, month_num::Integer)
    lat = deg2rad(lat)
    sza = deg2rad(sza)

    cos_lat = cos(lat)
    cos_sza = cos(sza)

    month_scalar = cospi(2*(month_num-0.5)/12)

    hprime = 74.37 - 8.097*cos_sza + 5.779*cos_lat - 1.213*month_scalar
    beta = 0.3849 - 0.1658*cos_sza - 0.08584*cos_lat + 0.1296*month_scalar

    return hprime, beta
end
ferguson(lat, sza, dt::DateTime) = ferguson(lat, sza, month(dt))

"""
    flatlinearterminator(sza, hp_day=74, hp_night=86, b_day=0.3, b_night=0.5) → (h′, β)

A simple terminator model with flat day and night ionospheres and a linear transition
between `sza_min` and `sza_max`. Solar zenith angle `sza` is in degrees.
"""
function flatlinearterminator(sza, hp_day=74, hp_night=86, b_day=0.3, b_night=0.5;
    sza_min=90, sza_max=100)

    if sza_min < sza < sza_max
        # Linear trend between day and night
        m_hprime = (hp_night - hp_day) / (sza_max - sza_min)
        m_beta = (b_night - b_day) / (sza_max - sza_min)

        hprime = m_hprime*(sza - sza_min) + hp_day
        beta = m_beta*(sza - sza_min) + b_day
    elseif sza < sza_min
        hprime = hp_day
        beta = b_day
    else  # sza > sza_max
        hprime = hp_night
        beta = b_night
    end

    return hprime, beta
end

"""
    flatlinearterminator(sza::AbstractArray; hp_day=74, hp_night=86, b_day=0.3, b_night=0.5) → (h′, β)

Given an array of `sza`, return the corresponding arrays of h′s and βs.
"""
function flatlinearterminator(sza::AbstractArray, hp_day=74, hp_night=86, b_day=0.3, b_night=0.5;
    sza_min=90, sza_max=100)

    # Linear trend between day and night
    m_hprime = (hp_night - hp_day) / (sza_max - sza_min)
    m_beta = (b_night - b_day) / (sza_max - sza_min)

    hprimes = similar(sza, Float64)
    betas = similar(hprimes)
    for i in eachindex(sza)
        if sza_min < sza[i] < sza_max
            hprimes[i] = m_hprime*(sza[i] - sza_min) + hp_day
            betas[i] = m_beta*(sza[i] - sza_min) + b_day
        elseif sza[i] < sza_min
            hprimes[i] = hp_day
            betas[i] = b_day
        else  # sza > sza_max
            hprimes[i] = hp_night
            betas[i] = b_night
        end
    end

    return hprimes, betas
end

"""
    smoothterminator(sza, hp_day=74, hp_night=86, b_day=0.3, b_night=0.5;
        sza_min=90, sza_max=100, steepness=1) → (h′, β)

A smooth terminator model using a logistic curve from day to night ionospheres as a function
of solar zenith angle `sza`.

The `steepness` parameter `1` with the default values of `sza_min` and `sza_max`
approximately asymptotes to the day and night values by `sza_min` and `sza_max` and is
centered on their midpoint.
"""
function smoothterminator(sza, hp_day=74, hp_night=86, b_day=0.3, b_night=0.5;
    sza_min=90, sza_max=100, steepness=1)

    hp_night > hp_day || throw(ArgumentError("`hp_night` should be greater than `hp_day`."))
    b_night > b_day || throw(ArgumentError("`b_night` should be greater than `b_day`."))
    sza_max > sza_min || throw(ArgumentError("`sza_max` should be greater than `sza_min`."))

    hprime_range = hp_night - hp_day
    beta_range = b_night - b_day
    sza_mid = (sza_max - sza_min)/2 + sza_min

    # logistic curve
    hprime = hprime_range/(1 + exp(-steepness*(sza - sza_mid))) + hp_day
    beta = beta_range/(1 + exp(-steepness*(sza - sza_mid))) + b_day

    return hprime, beta
end

"""
    fourierperturbation(sza, coeff; a=deg2rad(45)/2)

Return the perturbation produced by the Fourier series:
```
pert = coeff[1] + coeff[2]*cos(sza*π/a) + coeff[3]*cos(sza*2π/a) + coeff[4]*sin(sza*π/a) + coeff[5]*sin(sza*2π/a)
```
for solar zenith angle `sza` in degrees.

This can be used to produce simulated ionospheres in conjunction with [`ferguson`](@ref).

# Examples
```julia-repl
julia> lat, lon = 45, -105
julia> dt = DateTime(2020, 1, 15);
julia> sza = zenithangle(lat, lon, dt);
julia> h, b = ferguson(lat, sza, dt);

julia> coeffH = (1.345, 0.668, -0.177, -0.248, 0.096);
julia> hpert = fourierperturbation(sza, coeffH)
1.7991484480844766

julia> h + hpert
79.5769471288202
```
"""
function fourierperturbation(sza, coeff; a=deg2rad(45)/2)
    sza = deg2rad(sza)

    pert = coeff[1] +
        coeff[2]*cospi(sza/a) +
        coeff[4]*sinpi(sza/a) +
        coeff[3]*cospi(sza*2/a) +
        coeff[5]*sinpi(sza*2/a)

    return pert
end

end  # module


text:
@discord_object struct Webhook
    id::Snowflake
    type::WebhookType.WebhookTypeEnum
    guild_id::Snowflake
    channel_id::Snowflake
    user::User
    name::String
    avatar::String
    token::String
end

ext: jl


text:
include("binarytree.jl")

function optimaltree(network, optdata::Dict; verbose::Bool = false)
    numtensors = length(network)
    allindices = unique(vcat(network...))
    numindices = length(allindices)
    costtype = valtype(optdata)
    allcosts = [get(optdata, i, one(costtype)) for i in allindices]
    maxcost = addcost(mulcost(reduce(mulcost, allcosts; init=one(costtype)), maximum(allcosts)), zero(costtype)) # add zero for type stability: Power -> Poly
    tensorcosts = Vector{costtype}(undef, numtensors)
    for k = 1:numtensors
        tensorcosts[k] = mapreduce(i->get(optdata, i, one(costtype)), mulcost, network[k], init=one(costtype))
    end
    initialcost = min(maxcost, addcost(mulcost(maximum(tensorcosts),minimum(tensorcosts)),zero(costtype))) # just some arbitrary guess

    if numindices <= 32
        return _optimaltree(UInt32, network, allindices, allcosts, initialcost, maxcost; verbose = verbose)
    elseif numindices <= 64
        return _optimaltree(UInt64, network, allindices, allcosts, initialcost, maxcost; verbose = verbose)
    elseif numindices <= 128 && !(@static Int == Int32 && Sys.iswindows() ? true : false)
        return _optimaltree(UInt128, network, allindices, allcosts, initialcost, maxcost; verbose = verbose)
    else
        return _optimaltree(BitVector, network, allindices, allcosts, initialcost, maxcost; verbose = verbose)
    end
end

storeset(::Type{BitSet}, ints, maxint) = sizehint!(BitSet(ints), maxint)
function storeset(::Type{BitVector}, ints, maxint)
    set = falses(maxint)
    for i in ints
        set[i] = true
    end
    return set
end
function storeset(::Type{T}, ints, maxint) where {T<:Unsigned}
    set = zero(T)
    u = one(T)
    for i in ints
        set |= (u<<(i-1))
    end
    return set
end
_intersect(s1::T, s2::T) where {T<:Unsigned} = s1 & s2
_intersect(s1::BitVector, s2::BitVector) = s1 .& s2
_intersect(s1::BitSet, s2::BitSet) = intersect(s1, s2)
_union(s1::T, s2::T) where {T<:Unsigned} = s1 | s2
_union(s1::BitVector, s2::BitVector) = s1 .| s2
_union(s1::BitSet, s2::BitSet) = union(s1, s2)
_setdiff(s1::T, s2::T) where {T<:Unsigned} = s1 & (~s2)
_setdiff(s1::BitVector, s2::BitVector) = s1 .& (.~s2)
_setdiff(s1::BitSet, s2::BitSet) = setdiff(s1, s2)
_isemptyset(s::Unsigned) = iszero(s)
_isemptyset(s::BitVector) = !any(s)
_isemptyset(s::BitSet) = isempty(s)

addcost(cost1::Integer, cost2::Integer, costs::Vararg{<:Integer}) = Base.Checked.checked_add(cost1, cost2, costs...)
mulcost(cost1::Integer, cost2::Integer) = Base.Checked.checked_mul(cost1, cost2)
addcost(cost1, cost2, costs...) = +(cost1, cost2, costs...)
mulcost(cost1, cost2) = cost1*cost2
function computecost(allcosts, ind1::T, ind2::T) where {T<:Unsigned}
    cost = one(eltype(allcosts))
    ind = _union(ind1, ind2)
    n = 1
    @inbounds while !iszero(ind)
        if isodd(ind)
            cost = mulcost(cost, allcosts[n])
        end
        ind = ind>>1
        n += 1
    end
    return cost
end
function computecost(allcosts, ind1::BitVector, ind2::BitVector)
    cost = one(eltype(allcosts))
    ind = _union(ind1, ind2)
    @inbounds for n in findall(ind)
        cost = mulcost(cost, allcosts[n])
    end
    return cost
end
function computecost(allcosts, ind1::BitSet, ind2::BitSet)
    cost = one(eltype(allcosts))
    ind = _union(ind1, ind2)
    @inbounds for n in ind
        cost = mulcost(cost, allcosts[n])
    end
    return cost
end

function _optimaltree(::Type{T}, network, allindices, allcosts::Vector{S}, initialcost::C, maxcost::C; verbose::Bool = false) where {T,S,C}
    numindices = length(allindices)
    numtensors = length(network)
    indexsets = Array{T}(undef, numtensors)

    tabletensor = zeros(Int, (numindices,2))
    tableindex = zeros(Int, (numindices,2))

    adjacencymatrix = falses(numtensors,numtensors)
    costfac = maximum(allcosts)

    @inbounds for n = 1:numtensors
        indn = findall(in(network[n]), allindices)
        indexsets[n] = storeset(T, indn, numindices)
        for i in indn
            if tabletensor[i,1] == 0
                tabletensor[i,1] = n
                tableindex[i,1] = _findfirst(isequal(allindices[i]), network[n])
            elseif tabletensor[i,2] == 0
                tabletensor[i,2] = n
                tableindex[i,2] = _findfirst(isequal(allindices[i]), network[n])
                n1 = tabletensor[i,1]
                adjacencymatrix[n1,n] = true
                adjacencymatrix[n,n1] = true
            else
                error("no index should appear more than two times: $i")
            end
        end
    end
    componentlist = connectedcomponents(adjacencymatrix)
    numcomponent = length(componentlist)

    # generate output structures
    costlist = Vector{C}(undef, numcomponent)
    treelist = Vector{Any}(undef, numcomponent)
    indexlist = Vector{T}(undef, numcomponent)

    # run over components
    @inbounds for c=1:numcomponent
        # find optimal contraction for every component
        component = componentlist[c]
        componentsize = length(component)
        costdict = Vector{Dict{T, C}}(undef, componentsize)
        treedict = Vector{Dict{T, Any}}(undef, componentsize)
        indexdict = Vector{Dict{T, T}}(undef, componentsize)

        for k=1:componentsize
            costdict[k] = Dict{T, C}()
            treedict[k] = Dict{T, Any}()
            indexdict[k] = Dict{T, T}()
        end

        for i in component
            s = storeset(T, [i], numtensors)
            costdict[1][s] = zero(C)
            treedict[1][s] = i
            indexdict[1][s] = indexsets[i]
        end

        # run over currentcost
        currentcost = initialcost
        previouscost = zero(initialcost)
        while currentcost <= maxcost
            nextcost = maxcost
            # construct all subsets of n tensors that can be constructed with cost <= currentcost
            for n=2:componentsize
                verbose && println("Constructing subsets of size $n with cost $currentcost")
                # construct subsets by combining two smaller subsets
                for k = 1:div(n-1,2)
                    for s1 in keys(costdict[k]), s2 in keys(costdict[n-k])
                        if _isemptyset(_intersect(s1, s2)) && get(costdict[n], _union(s1, s2), currentcost) > previouscost
                            ind1 = indexdict[k][s1]
                            ind2 = indexdict[n-k][s2]
                            cind = _intersect(ind1, ind2)
                            if !_isemptyset(cind)
                                s = _union(s1, s2)
                                cost = addcost(costdict[k][s1], costdict[n-k][s2], computecost(allcosts, ind1, ind2))
                                if cost <= get(costdict[n], s, currentcost)
                                    costdict[n][s] = cost
                                    indexdict[n][s] = _setdiff(_union(ind1,ind2), cind)
                                    treedict[n][s] = BinaryTreeNode(treedict[k][s1], treedict[n-k][s2])
                                elseif currentcost < cost < nextcost
                                    nextcost = cost
                                end
                            end
                        end
                    end
                end
                if iseven(n) # treat the case k = n/2 special
                    k = div(n,2)
                    it = keys(costdict[k])
                    elstate1 = iterate(it)
                    while elstate1 !== nothing
                        s1, nextstate1 = elstate1
                        elstate2 = iterate(it, nextstate1)
                        while elstate2 !== nothing
                            s2, nextstate2 = elstate2
                            if _isemptyset(_intersect(s1, s2)) && get(costdict[n], _union(s1, s2), currentcost) > previouscost
                                ind1 = indexdict[k][s1]
                                ind2 = indexdict[k][s2]
                                cind = _intersect(ind1, ind2)
                                if !_isemptyset(cind)
                                    s = _union(s1, s2)
                                    cost = addcost(costdict[k][s1], costdict[k][s2], computecost(allcosts, ind1, ind2))
                                    if cost <= get(costdict[n], s, currentcost)
                                        costdict[n][s] = cost
                                        indexdict[n][s] = _setdiff(_union(ind1,ind2), cind)
                                        treedict[n][s] = BinaryTreeNode(treedict[k][s1], treedict[k][s2])
                                    elseif currentcost < cost < nextcost
                                        nextcost = cost
                                    end
                                end
                            end
                            elstate2 = iterate(it, nextstate2)
                        end
                        elstate1 = iterate(it, nextstate1)
                    end
                end
            end
            currentbiggestset = _findlast(x->!isempty(x), costdict)
            verbose && println("Finished at cost $currentcost: maximum subset has size $currentbiggestset")
            if !isempty(costdict[componentsize])
                break
            end
            previouscost = currentcost
            currentcost = min(maxcost, nextcost*costfac)
        end
        if isempty(costdict[componentsize])
            error("Maxcost $maxcost reached without finding solution") # should be impossible
        end
        s = storeset(T, component, numtensors)
        costlist[c] = costdict[componentsize][s]
        treelist[c] = treedict[componentsize][s]
        indexlist[c] = indexdict[componentsize][s]
    end
    p  = sortperm(costlist)

    tree = treelist[p[1]]
    cost = costlist[p[1]]
    ind = indexlist[p[1]]
    for c = 2:numcomponent
        tree = BinaryTreeNode(tree, treelist[p[c]])
        cost = addcost(cost, costlist[p[c]], computecost(allcosts, ind, indexlist[p[c]]))
        ind = _union(ind, indexlist[p[c]])
    end
    return tree, cost
end

# For a given adjacency matrix of size n x n, connectedcomponents returns
# a list componentlist that contains integer vectors, where every integer
# vector groups the indices of the vertices of a connected component of the
# graph encoded by A. The number of connected components is given by
# length(componentlist).
#
# Used as auxiliary function to analyze contraction graph in contract.
function connectedcomponents(A::AbstractMatrix{Bool})
    n = size(A,1)
    @assert size(A,2) == n

    componentlist = Vector{Vector{Int}}(undef, 0)
    assignedlist = falses((n,))

    for i = 1:n
        if !assignedlist[i]
            assignedlist[i] = true
            checklist = [i]
            currentcomponent = [i]
            while !isempty(checklist)
                j = pop!(checklist)
                for k = findall(A[j,:])
                    if !assignedlist[k]
                        push!(currentcomponent, k)
                        push!(checklist,k)
                        assignedlist[k] = true;
                    end
                end
            end
            push!(componentlist, currentcomponent)
        end
    end
    return componentlist
end


text:
"""
    SemigroupMethod(epsilon, delta, max_iter)

Semigroup method from Algorithm 1 in [1].

[1] Taghvaei, A., & Mehta, P. G. (2016). Gain function approximation in the feedback particle filter. In 2016 IEEE 55th Conference on Decision and Control (CDC) (pp. 5446–5452). IEEE. https://doi.org/10.1109/CDC.2016.7799105

    SemigroupMethod(epsilon, delta, max_iter, lambda)

Semigroup method with regularization parameter `lambda`.
"""
struct SemigroupMethod <: GainEstimationMethod
    epsilon::Float64
    delta::Float64
    max_iter::Int
    lambda::Float64
    SemigroupMethod(eps, delta, max_iter, lambda=0.) = new(eps, delta, max_iter, lambda)
end

SemigroupMethod(epsilon::Float64, delta::Float64) = SemigroupMethod(epsilon, delta, 100)

function solve!(eq::PoissonEquation, method::SemigroupMethod)
    N = length(eq.positions)
    H = Htilde(eq)
    broadcast!(*, H, H, method.epsilon)

    # compute T operator
    T = zeros(eltype(eq.positions), N, N)
    n = state_dim(eq)
    for i in 1:N
        T[i,i] = one(eltype(eq.positions))
        for j in i+1:N
            sq = zero(eltype(eq.positions))
            for l in 1:n
                sq += (eq.positions[l,i] - eq.positions[l,j])^2
            end
            T[i,j] = exp(-sq/(4*method.epsilon))
            T[j,i] = T[i,j]
        end
    end
    broadcast!(/, T, T, sqrt.(sum(T,dims=1) .* sum(T,dims=2)))
    broadcast!(/, T, T, sum(T,dims=2))

    # add noise to regularize T
    if method.lambda > 0.
        for i in 1:N
            T[i,i] -=  method.lambda * rand(Distributions.Uniform(0.9,1))
        end
    end

    # solve fixed-point equation for potential
    newpotential = copy(eq.potential)::Array{Float64,2}
    fluctuation = 1.
    n = 1
    while fluctuation > method.delta
        LinearAlgebra.mul!(newpotential, eq.potential, T')
        broadcast!(+, newpotential, newpotential, H)
        broadcast!(-, newpotential, newpotential, Statistics.mean(newpotential, dims=2))
        fluctuation = maximum(abs.(newpotential - eq.potential))
        eq.potential .= newpotential
        n += 1
        if n == method.max_iter
            print("!")
            break
        end
    end

    # compute gain from potential
    gainhelper_semigroup!(eq, T)
    broadcast!(/, eq.gain, eq.gain, 2*method.epsilon)
end

function gainhelper_semigroup!(eq::PoissonEquation, T::AbstractMatrix)
    pos     = eq.positions
    pot     = eq.potential
    gain    = eq.gain

    n       = state_dim(eq)
    m       = obs_dim(eq)
    N       = no_of_particles(eq)

    @inbounds for l in 1:n, k in 1:m, i in 1:N
        Tpotpos = zero(eltype(gain))
        Tpot    = zero(eltype(gain))
        Tpos    = zero(eltype(gain))
        for j in 1:N
            Tpot        += T[i,j] * pot[k,j]
            Tpos        += T[i,j] * pos[l,j]
            Tpotpos     += T[i,j] * pot[k,j] * pos[l,j]
        end
        gain[l,i,k] = Tpotpos - Tpot * Tpos
    end
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
module Utilities

# import some common definitions from MAGMA
import MAGMA.libmagma

# Initialization and Finalization
export magmaInit, magmaFinalize

include("init.jl")

end  # module Utilities


text:
import OpenCL
const cl = OpenCL

const mult_kernel = """
   __kernel void mult(__global const float *a,
                      __global const float *b,
                      __global float *c)
    {
      int gid = get_global_id(0);
      c[gid] = a[gid] * b[gid];
    }
"""

const N = 100_000;
a = rand(Float32, N);
b = rand(Float32, N);

device, ctx, queue = cl.create_compute_context();

a_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=a);
b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b);
c_buff = cl.Buffer(Float32, ctx, :w, length(a));

p = cl.Program(ctx, source=mult_kernel) |> cl.build!;
k = cl.Kernel(p, "mult");

tm = @elapsed cl.call(queue, k, size(a), nothing, a_buff, b_buff, c_buff);
@printf "Time taken: %f\a sec.\n" tm

max_stars_repo_path: Module 3/Chapter10/cl-mmul.jl

max_issues_repo_path: Module 3/Chapter10/cl-mmul.jl

max_forks_repo_path: Module 3/Chapter10/cl-mmul.jl


text:
using PTools

shmcfg = [ShmCfg(:svar1, Int32, 64*1024), ShmCfg(:svar2, Uint8, (100,100))]
h = start_stasks(shmcfg)
ntasks = count_stasks(h)

svar1 = task_local_storage(:svar1)

for i in 1:64*1024
    svar1[i]=0
end
println(svar1)

offset_list = [i for i in 1:ntasks]
ntasks_list = [ntasks for i in 1:ntasks]

resp = pmap_stasks(h, (offset, ntasks) -> begin
                        svar1 = task_local_storage(:svar1)
                        svar2 = task_local_storage(:svar2)

                        mypid = myid()
                        for x in offset:ntasks:64*1024
                            svar1[x] = mypid
                        end

                        true

                    end,

                    offset_list, ntasks_list)

println(task_local_storage(:svar1))

stop_stasks(h, shmcfg)


text:
# An implementation of MINRES for the solution of the
# linear system Ax = b, or the linear least-squares problem
#
#  minimize ‖Ax - b‖
#
# where A is square and symmetric.
#
# MINRES is formally equivalent to applying the conjugate residuals method
# to Ax = b when A is positive definite, but is more general and also applies
# to the case where A is indefinite.
#
# This implementation follows the original implementation by
# Michael Saunders described in
#
# C. C. Paige and M. A. Saunders, Solution of Sparse Indefinite Systems
# of Linear Equations, SIAM Journal on Numerical Analysis, 12(4), 617-629, 1975.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Brussels, Belgium, June 2015.
# Montreal, August 2015.

export minres

"""Solve the shifted linear least-squares problem

  minimize ‖b - (A + λ I)x‖₂²

or the shifted linear system

  (A + λ I) x = b

using the MINRES method, where λ ≥ 0 is a shift parameter,
where A is square and symmetric.

MINRES is formally equivalent to applying CR to Ax=b when A is positive
definite, but is typically more stable and also applies to the case where
A is indefinite.

MINRES produces monotonic residuals ‖r‖₂ and optimality residuals ‖A'r‖₂.

A preconditioner M may be provided in the form of a linear operator and is
assumed to be symmetric and positive definite.
"""
function minres(A :: AbstractLinearOperator, b :: AbstractVector{T};
                M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
                atol :: Float64=1.0e-12, rtol :: Float64=1.0e-12, etol :: Float64=1.0e-8,
                window :: Int=5, itmax :: Int=0, conlim :: Float64=1.0e+8,
                verbose :: Bool=false) where T <: Number

  m, n = size(A)
  m == n || error("System must be square")
  size(b, 1) == m || error("Inconsistent problem size")
  verbose && @printf("MINRES: system of size %d\n", n)

  ϵM = eps(T)
  x = zeros(T, n)
  ctol = conlim > 0.0 ? 1 ./ conlim : 0.0;

  # Initialize Lanczos process.
  # β₁ M v₁ = b.
  r1 = copy(b)
  y = M * r1
  β₁ = @kdot(m, r1, y)
  β₁ < 0.0 && error("Preconditioner is not positive definite")
  β₁ == 0.0 && return (x, SimpleStats(true, true, [0.0], [0.0], "x = 0 is a zero-residual solution"))
  β₁ = sqrt(β₁)
  β = β₁

  oldβ = 0.0
  δbar = 0.0
  ϵ = 0.0
  rNorm = β₁
  rNorms = [β₁]
  ϕbar = β₁
  rhs1 = β₁
  rhs2 = 0.0
  γmax = 0.0
  γmin = Inf
  cs = -1.0
  sn = 0.0
  v = zeros(T, n)
  w = zeros(T, n)
  w1 = zeros(T, n)
  w2 = zeros(T, n)
  r2 = copy(r1)

  ANorm² = 0.0
  ANorm = 0.0
  Acond = 0.0
  ArNorm = 0.0
  ArNorms = [ArNorm]
  xNorm = 0.0

  xENorm² = 0.0
  err_lbnd = 0.0
  err_vec = zeros(T, window)

  verbose && @printf("%5s  %7s  %7s  %7s  %8s  %8s  %7s  %7s\n",
                     "Aprod", "‖r‖", "‖A'r‖", "β", "cos", "sin", "‖A‖", "κ(A)")
  verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e\n",
                     0, rNorm, ArNorm, β, cs, sn, ANorm, Acond)

  iter = 0
  itmax == 0 && (itmax = n)

  tol = atol + rtol * β₁
  status = "unknown"
  solved = solved_mach = solved_lim = (rNorm ≤ rtol)
  tired  = iter ≥ itmax
  ill_cond = ill_cond_mach = ill_cond_lim = false
  zero_resid = zero_resid_mach = zero_resid_lim = (rNorm ≤ tol)
  fwd_err = false

  while ! (solved || tired || ill_cond)
    iter = iter + 1

    # Generate next Lanczos vector.
    v = copy(y)
    @kscal!(n, 1 ./ β, v)
    y = A * v
    λ != 0.0 && @kaxpy!(n, -λ, v, y)          # (y = y - λ * v)
    iter ≥ 2 && @kaxpy!(n, -β / oldβ, r1, y)  # (y = y - β / oldβ * r1)

    α = dot(v, y)
    @kaxpy!(n, -α / β, r2, y)  # y = y - α / β * r2
    r1 = copy(r2)
    r2 = copy(y)
    y = M * r2
    oldβ = β
    β = @kdot(n, r2, y)
    β < 0.0 && error("Preconditioner is not positive definite")
    β = sqrt(β)
    ANorm² = ANorm² + α * α + oldβ * oldβ + β * β

    # Apply rotation to obtain
    #  [ δₖ    ϵₖ₊₁    ] = [ cs  sn ] [ δbarₖ  0    ]
    #  [ γbar  δbarₖ₊₁ ]   [ sn -cs ] [ αₖ     βₖ₊₁ ]
    oldϵ = ϵ
    δ = cs * δbar + sn * α
    γbar = sn * δbar - cs * α
    ϵ = sn * β
    δbar = -cs * β
    root = sqrt(γbar * γbar + δbar * δbar)
    ArNorm = ϕbar * root  # = ‖A'rₖ₋₁‖
    push!(ArNorms, ArNorm)

    # Compute the next plane rotation.
    γ = sqrt(γbar * γbar + β * β)
    γ = max(γ, ϵM)
    cs = γbar / γ
    sn = β / γ
    ϕ = cs * ϕbar
    ϕbar = sn * ϕbar

    # Update x.
    w1 = copy(w2)
    w2 = copy(w)
    w = (v - oldϵ * w1 - δ * w2) / γ
    @kaxpy!(n, ϕ, w, x)  # x = x + ϕ * w
    xENorm² = xENorm² + ϕ * ϕ

    # Compute lower bound on forward error.
    err_vec[mod(iter, window) + 1] = ϕ
    iter ≥ window && (err_lbnd = norm(err_vec))

    γmax = max(γmax, γ)
    γmin = min(γmin, γ)
    ζ = rhs1 / γ
    rhs1 = rhs2 - δ * ζ
    rhs2 = -ϵ * ζ

    # Estimate various norms.
    ANorm = sqrt(ANorm²)
    xNorm = @knrm2(n, x)
    ϵA = ANorm * ϵM
    ϵx = ANorm * xNorm * ϵM
    ϵr = ANorm * xNorm * rtol
    d = γbar
    d == 0.0 && (d = ϵA)

    rNorm = ϕbar

    test1 = rNorm / (ANorm * xNorm)
    test2 = root / ANorm
    push!(rNorms, rNorm)

    Acond = γmax / γmin

    verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e\n",
                       iter, test1, test2, β, cs, sn, ANorm, Acond)

    if iter == 1
      # A'b = 0 so x = 0 is a minimum least-squares solution
      β / β₁ ≤ 10 * ϵM && return (x, SimpleStats(true, true, [β₁], [0.0], "x is a minimum least-squares solution"))
    end

    # Stopping conditions that do not depend on user input.
    # This is to guard against tolerances that are unreasonably small.
    ill_cond_mach = (1.0 + 1.0 / Acond ≤ 1.0)
    solved_mach = (1.0 + test2 ≤ 1.0)
    zero_resid_mach = (1.0 + test1 ≤ 1.0)
    # solved_mach = (ϵx ≥ β₁)

    # Stopping conditions based on user-provided tolerances.
    tired = iter ≥ itmax
    ill_cond_lim = (1.0 / Acond ≤ ctol)
    solved_lim = (test2 ≤ tol)
    zero_resid_lim = (test1 ≤ tol)
    iter ≥ window && (fwd_err = err_lbnd ≤ etol * sqrt(xENorm²))

    zero_resid = zero_resid_mach | zero_resid_lim
    ill_cond = ill_cond_mach | ill_cond_lim
    solved = solved_mach | solved_lim | zero_resid | fwd_err
  end

  tired         && (status = "maximum number of iterations exceeded")
  ill_cond_mach && (status = "condition number seems too large for this machine")
  ill_cond_lim  && (status = "condition number exceeds tolerance")
  solved        && (status = "found approximate minimum least-squares solution")
  zero_resid    && (status = "found approximate zero-residual solution")
  fwd_err       && (status = "truncated forward error small enough")

  stats = SimpleStats(solved, !zero_resid, rNorms, ArNorms, status)
  return (x, stats)
end


text:
module TLA

using Pkg.Artifacts
artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
tla2tools_hash = artifact_hash("tla2tools", artifact_toml)

if tla2tools_hash == nothing || !artifact_exists(tla2tools_hash)
    tla2tools_hash = create_artifact() do artifact_dir
        tla2tools_url_base = "https://tla.msr-inria.inria.fr/tlatoolbox/dist/tla2tools.jar"
        download(
            "$(tla2tools_url_base)",
            joinpath(artifact_dir, "tla2tools.jar"),
        )
    end

    bind_artifact!(artifact_toml, "tla2tools", tla2tools_hash)
end

tla2tools_path = artifact"tla2tools"
tla2tools = tla2tools_path * "/tla2tools.jar"

function pcal(str)
    io = IOBuffer(str)
    run(pipeline(`java -cp $tla2tools pcal.trans $io`))
end

function pcal(str)
    try
        rm("test.tla")
    catch
    end
    io = open("test.tla", "w+")
    write(io, str)
    close(io)
    run(pipeline(`java -cp $tla2tools pcal.trans test.tla`))
    rm("test.tla")
end
pcal2(test)

tlc(str; default_ops = "-XX:+IgnoreUnrecognizedVMOptions - XX:+UseParallelGC") =
    run(`java $default_ops -cp $tla2tools tlc2.TLC $str`)
sany(str) = run(`java -cp $tla2tools tla2sany.SANY $str`)
tlatex(str) = run(`java -cp $tla2tools tla2tex.TLA $str`)

# use \bigvee<TAB> == ⋁ and \bigwedge<TAB> == ⋀
const unicode_dict = Dict(
    "∈" => raw"\in",
    "∀" => raw"\A",
    "∃" => raw"\E",
    "⋀" => raw"/\\",
    "⋁" => raw"\/",
)
#remove_unicode(str) = replace(str, r"∈|∀|∃|⋀|⋁" => unicode_dict[])
function replace_unicode(str)
    str = replace(str, "∈" => raw"\in")
    str = replace(str, '∀' => raw"\A")
    str = replace(str, '∃' => raw"\E")
    str = replace(str, '⋀' => raw"/\\")
    str = replace(str, '⋁' => raw"\/")
    str = replace(str, '↦' => "|->")
    str = replace(str, '⋃' => raw"\union")
    str = replace(str, '⋂' => raw"\intersect")
    str = replace(str, '∉' => raw"\notin")
    str = replace(str, '⟨' => "<<")
    str = replace(str, '⟩' => ">>")
    str = replace(str, '¬' => '~')
    str = replace(str, '⊆' => raw"\subseteq")
    str = replace(str, '≠' => "/=")
    str
end

replace_unicode("yow")

test = """
EXTENDS Integers
(*--algorithm wire
    variables
        people = {"alice", "bob"},
        acc = [alice |-> 5, bob |-> 5];
begin
    skip;
end algorithm;*)
"""

replace_unicode(test) |> pcal

end


text:
I_trues =
[0.44129802 0.39444956 0.30091208 0.25465866 0.30253774 0.39726529 0.44454934 ;
 0.39250948 0.35223311 0.27210130 0.23339588 0.27639321 0.35966691 0.40109329 ;
 0.35082997 0.31575536 0.24622966 0.21348201 0.25258702 0.32676664 0.36354469 ;
 0.29863696 0.26973242 0.21277028 0.18702656 0.22140580 0.28468958 0.31590800 ;
 0.26939388 0.24383765 0.19368088 0.17169133 0.20342523 0.26071536 0.28888259 ;
 0.22153000 0.20137817 0.16217962 0.14615212 0.17344750 0.22089470 0.24406575 ;
 0.20185562 0.18392101 0.14920591 0.13557847 0.16097330 0.20430272 0.22539038 ;
 0.16889020 0.15469709 0.12752450 0.11786073 0.13989890 0.17613018 0.19363900 ;
 0.13097041 0.12117362 0.10280748 0.09759154 0.11532259 0.14285044 0.15600064 ;
 0.10231626 0.09595756 0.08440503 0.08239375 0.09628240 0.11652976 0.12607100 ;
 0.08683737 0.08239718 0.07459712 0.07418876 0.08561228 0.10147599 0.10886768 ;
 0.06783499 0.06583006 0.06269131 0.06393739 0.07157324 0.08121402 0.08559886 ;
 0.05781217 0.05713933 0.05643322 0.05814731 0.06295443 0.06843439 0.07085458 ;
 0.05385796 0.05372482 0.05391877 0.05550325 0.05861141 0.06185271 0.06324324 ;
 0.05240320 0.05247332 0.05294178 0.05423678 0.05628819 0.05826948 0.05909602 ;
 0.05300496 0.05300496 0.05300496 0.05300496 0.05300496 0.05300496 0.05300496 ]

Q_trues =
[ -0.01753141 -0.06485313 -0.15965601 -0.20757277 -0.16128167 -0.06766886 -0.02078273;
  -0.01772108 -0.05944266 -0.14330674 -0.18659926 -0.14759865 -0.06687645 -0.02630489;
  -0.01524525 -0.05274909 -0.12838028 -0.16821107 -0.13473764 -0.06376036 -0.02795997;
  -0.00990697 -0.04270497 -0.10914792 -0.14510675 -0.11778344 -0.05766213 -0.02717801;
  -0.00608288 -0.03648858 -0.09825566 -0.13222805 -0.10800002 -0.05336629 -0.02557158;
   0.00131430 -0.02554220 -0.08036031 -0.11134114 -0.09162819 -0.04505874 -0.02122145;
   0.00477642 -0.02076298 -0.07299588 -0.10284244 -0.08476326 -0.04114469 -0.01875835;
   0.01119511 -0.01235218 -0.06066038 -0.08873702 -0.07303478 -0.03378526 -0.01355369;
   0.01968930 -0.00196538 -0.04650218 -0.07273772 -0.05901730 -0.02364220 -0.00534093;
   0.02701775  0.00644566 -0.03586341 -0.06078292 -0.04774078 -0.01412655  0.00326301;
   0.03129776  0.01117009 -0.03016559 -0.05432509 -0.04118074 -0.00790872  0.00926745;
   0.03668441  0.01698497 -0.02328501 -0.04623546 -0.03216694  0.00160101  0.01892054;
   0.03923450  0.01977043 -0.01979730 -0.04164830 -0.02631851  0.00847536  0.02619209;
   0.03983604  0.02053460 -0.01852852 -0.03954758 -0.02322116  0.01240671  0.03045076;
   0.03969919  0.02052770 -0.01814346 -0.03853982 -0.02148987  0.01473155  0.03300636;
   0.03755859  0.01877930 -0.01877930 -0.03755859 -0.01877930  0.01877930  0.03755859]

U_trues =
[ 0 0.04390364 0.07365528 0.08128298  0.06713098  0.03737934 0 ;
  0 0.04428736 0.07046980 0.07153176  0.05342685  0.02724440 0 ;
  0 0.04435776 0.06762732 0.06357361  0.04248541  0.01921585 0 ;
  0 0.04408551 0.06384064 0.05397200  0.02964161  0.00988649 0 ;
  0 0.04374666 0.06158007 0.04872177  0.02280850  0.00497510 0 ;
  0 0.04291904 0.05764878 0.04024242  0.01205313 -0.00267662 0 ;
  0 0.04248227 0.05594215 0.03677307  0.00775068 -0.00570920 0 ;
  0 0.04161531 0.05293867 0.03093600  0.00064404 -0.01067932 0 ;
  0 0.04036105 0.04917038 0.02406753 -0.00748419 -0.01629352 0 ;
  0 0.03913232 0.04592516 0.01855839 -0.01378109 -0.02057393 0 ;
  0 0.03828742 0.04388718 0.01529883 -0.01738884 -0.02298860 0 ;
  0 0.03686504 0.04073530 0.01057373 -0.02242105 -0.02629131 0 ;
  0 0.03563204 0.03822653 0.00708827 -0.02594929 -0.02854377 0 ;
  0 0.03475930 0.03654849 0.00488817 -0.02808194 -0.02987113 0 ;
  0 0.03413018 0.03538005 0.00341471 -0.02946560 -0.03071547 0 ;
  0 0.03252669 0.03252669 0.00000000 -0.03252669 -0.03252669 0 ; ]


text:
ConstraintCardType = CardType{:constraint}()
const ConstraintCard = NastranCard{ConstraintCardType}

type MPC <: ConstraintCard
    id::Int64
end
function convert(::Type{MPC},card)
    id = card[2]::Int64
    MPC(id)
end


text:
@testset "Bitpattern entropy" begin
    for N in [100,1000,10000,100000]
        # every bitpattern is only hit once, hence entropy = log2(N)
        @test isapprox(log2(N),bitpattern_entropy(rand(Float32,N)),atol=1e-1)
        @test isapprox(log2(N),bitpattern_entropy(rand(Float64,N)),atol=1e-1)
    end

    N = 1000_000   # more bitpattern than there are in 8 or 16-bit
    @test isapprox(16.0,bitpattern_entropy(rand(UInt16,N)),atol=1e-1)
    @test isapprox(16.0,bitpattern_entropy(rand(Int16,N)),atol=1e-1)

    @test isapprox(8.0,bitpattern_entropy(rand(UInt8,N)),atol=1e-1)
    @test isapprox(8.0,bitpattern_entropy(rand(Int8,N)),atol=1e-1)
end

@testset "Bitcount" begin
    @test bitcount(UInt8[1,2,4,8,16,32,64,128]) == ones(8)
    @test bitcount(collect(0x0000:0xffff)) == 2^15*ones(16)

    N = 100_000
    c = bitcount(rand(N))
    @test c[1] == 0         # sign always 0
    @test c[2] == 0         # first expbit always 0, i.e. U(0,1) < 1
    @test c[3] == N         # second expont always 1

    @test all(isapprox.(c[15:50],N/2,rtol=1e-1))
end

@testset "Bitcountentropy" begin

    # test the PRNG on uniformity
    N = 100_000
    H = bitcount_entropy(rand(UInt8,N))
    @test all(isapprox.(H,ones(8),rtol=5e-4))

    H = bitcount_entropy(rand(UInt16,N))
    @test all(isapprox.(H,ones(16),rtol=5e-4))

    H = bitcount_entropy(rand(UInt32,N))
    @test all(isapprox.(H,ones(32),rtol=5e-4))

    H = bitcount_entropy(rand(UInt64,N))
    @test all(isapprox.(H,ones(64),rtol=5e-4))

    # also for random floats
    H = bitcount_entropy(rand(N))
    @test H[1:5] == zeros(5)    # first bits never change
    @test all(isapprox.(H[16:55],ones(40),rtol=1e-4))
end

@testset "Bitinformation random" begin
    N = 100_000
    A = rand(UInt64,N)
    @test all(bitinformation(A) .< 1e-3)

    A = rand(Float32,N)
    bi = bitinformation(A)
    @test all(bi[1:4] .== 0.0) # the first bits are always 0011
    sort!(A)
    @test sum(bitinformation(A)) > 12
end

@testset "Bitinformation in chunks" begin
    N = 1000
    Nhalf = N ÷ 2
    A = rand(UInt32,N)
    n11,npair1,N1 = bitcount(A[1:Nhalf-1]),bitpaircount(A[1:Nhalf]),Nhalf
    n12,npair2,N2 = bitcount(A[Nhalf:end-1]),bitpaircount(A[Nhalf:end]),Nhalf-1
    @test bitinformation(n11+n12,npair1+npair2,N1+N2) == bitinformation(A)
end

@testset "Bitinformation dimensions" begin
    A = rand(Float32,30,40,50)
    @test bitinformation(A) == bitinformation(A,dims=1)

    # check that the :all_dimensions flag is
    bi1 = bitinformation(A,dims=1)
    bi2 = bitinformation(A,dims=2)
    bi3 = bitinformation(A,dims=3)
    @test bitinformation(A,:all_dimensions) == ((bi1 .+ bi2 .+ bi3)/3)

    # check that there is indeed more information in the sorted dimensions
    sort!(A,dims=2)
    @test sum(bitinformation(A,dims=1)) < sum(bitinformation(A,dims=2))
end

@testset "Mutual information" begin
    N = 10_000

    # equal probabilities for 00|01|10|11
    @test 0.0 == mutual_information([0.25 0.25;0.25 0.25])

    # every 0 or 1 in A is also a 0 or 1 in B
    @test 1.0 == mutual_information([0.5 0.0;0.0 0.5])

    # as before but more 1s means a lower entropy
    @test entropy([0.25,0.75],2) == mutual_information([0.25 0.0;0.0 0.75])

    # every bit is inverted
    @test 1.0 == mutual_information([0.0 0.5;0.5 0.0])

    # two independent arrays
    for T in [UInt8,UInt16,UInt32,UInt64,Float16,Float32,Float64]
        mutinf = bitinformation(rand(T,N),rand(T,N))
        for m in mutinf
            @test isapprox(0,m,atol=1e-3)
        end
    end

    # mutual information of identical arrays
    # for 0,1 occuring exactly 50/50 this is 1bit
    # but so in practice slightly lower for rand(UInt),
    # or clearly lower for low entropy bits in Float16/32/64
    # but identical to the bitcount_entropy (up to rounding errors)
    for T in [UInt8,UInt16,UInt32,UInt64,Float16,Float32,Float64]
        R = rand(T,N)
        mutinf = bitinformation(R,R)
        @test mutinf ≈ bitcount_entropy(R)
    end
end

@testset "Redundancy" begin
    N = 100_000

    # No redundancy in independent arrays
    for T in [UInt8,UInt16,UInt32,UInt64,Float16,Float32,Float64]
        redun = redundancy(rand(T,N),rand(T,N))
        for r in redun
            @test isapprox(0,r,atol=2e-3)
        end
    end

    # Full redundancy in identical arrays
    for T in [UInt8,UInt16,UInt32,UInt64,Float16,Float32,Float64]
        A = rand(T,N)
        H = bitcount_entropy(A)
        R = redundancy(A,A)
        for (r,h) in zip(R,H)
            if iszero(h)
                @test iszero(r)
            else
                @test isapprox(1,r,atol=2e-3)
            end
        end
    end
end

@testset "Mutual information with round to nearest" begin
    N = 100_000

    # compare shaving to round to nearest
    # for round to nearest take m more bits into account for the
    # joint probability
    m = 8

    for T in [Float32,Float64]
        R = rand(T,N)
        for keepbit in [5,10,15]
            mutinf_shave = bitinformation(R,shave(R,keepbit))
            mutinf_round = bitinformation(R,round(R,keepbit),m)
            for (s,r) in zip(mutinf_shave,mutinf_round)
                @test isapprox(s,r,atol=1e-2)
            end
        end
    end

    # shaving shouldn't change
    for T in [Float32,Float64]
        R = rand(T,N)
        for keepbit in [5,10,15]
            mutinf_shave = bitinformation(R,shave(R,keepbit))
            mutinf_round = bitinformation(R,shave(R,keepbit),m)
            for (s,r) in zip(mutinf_shave,mutinf_round)
                @test isapprox(s,r,atol=1e-2)
            end
        end
    end
end

@testset "Information of random set to zero" begin
    for T in (UInt32,UInt64,Float32,Float64)
        for N in [10_000,100_000]
            A = rand(T,N)
            # increase confidence here from the default 0.99 to avoid test failures
            # from false positives
            b = bitinformation(A,confidence=0.9999)
            for ib in b
                @test 0 == ib
            end

            m = bitinformation(A[1:end-1],A[2:end],confidence=0.9999)
            for im in m
                @test 0 == im
            end
        end
    end
end


text:
# Gutzwiller projected slater determinant wavefunction
# TODO: Does Gutzwiller need to carry the states?
mutable struct GutzwillerSlater{T}
    states :: Matrix{T}
    configuration :: Vector{Int}

    uplist :: Vector{Int}
    dnlist :: Vector{Int}

    upslater :: DetMatrix
    dnslater :: DetMatrix
end

"""
    random_gutzwiller_half(states)

returns a Gutzwiller wavefunction with randomly chosen half up spins
and half down spins.  The input `states` is a matrix with spatial
orbitals as columns.  The number of orbitals(columns) has to be
exactly half of the number of sites(rows).

"""
function random_gutzwiller_half(states::Matrix)
    num_sites = size(states, 1)
    num_states = size(states, 2)

    # even lattice size and half-filling
    @assert 2*num_states == num_sites

    ### TODO: can these be simplified in one function?
    configuration = random_ones(num_sites, num_states)
    uplist = findall(x->x==1, configuration)
    dnlist = findall(x->x==0, configuration)

    ### TODO: should be possible to slice them from states, directly!
    ## NOTE: while the states(orbitals) are originally given on cols,
    ## here we put states on rows and positions on columns. Notice the
    ## order of the for expression!

    upslater = DetMatrix([states[x,k] for k=1:num_states, x=uplist])
    dnslater = DetMatrix([states[x,k] for k=1:num_states, x=dnlist])

    return GutzwillerSlater(states, configuration,
                            uplist, dnlist,
                            upslater, dnslater)
end

"""
    propose_step!(wf)

propose and, if accepted, update a step for the Gutzwiller
wavefunction. The process is to choose randomly one up and down spin
fermion and swap their positions. The acceptance probability is the
minimum of 1 and square of the determinant ratio.

returns `true` if update the `wf`, and `false` otherwise.
"""
function propose_step!(wf::GutzwillerSlater)
    ### TODO: split these operations into multiple functions.
    # randomly choose an up and down
    ## NOTE: `up/dn_index` are random indeces in `up/dnlist` pointing to
    ## positions, namely, `up/dn_site` in the wf.configuration!
    up_index = rand(1:length(wf.uplist))
    dn_index = rand(1:length(wf.dnlist))
    up_site = wf.uplist[up_index]
    dn_site = wf.dnlist[dn_index]

    # new columns, i.e. just swap up and dn
    new_up_col = wf.dnslater.matrix[:, dn_index]
    new_dn_col = wf.upslater.matrix[:, up_index]

    # find the probability of the change
    up_detrat = det_ratio_1col(wf.upslater, new_up_col, up_index)
    dn_detrat = det_ratio_1col(wf.dnslater, new_dn_col, dn_index)

    # up_detrat = exact_det_ratio_1col(wf.upslater, new_up_col, up_index)
    # dn_detrat = exact_det_ratio_1col(wf.dnslater, new_dn_col, dn_index)

    if rand(Float64) < min(1, norm(up_detrat * dn_detrat)^2)
        # update the Gutzwiller wavefunction
        wf.configuration[up_site] = 0
        wf.configuration[dn_site] = 1

        wf.uplist[up_index], wf.dnlist[dn_index] = dn_site, up_site

        update_detmatrix_1col!(wf.upslater,
                               new_up_col, up_index, up_detrat)
        update_detmatrix_1col!(wf.dnslater,
                               new_dn_col, dn_index, dn_detrat)

        # exact_update_detmatrix_1col!(wf.upslater,
        #                        new_up_col, up_index, up_detrat)
        # exact_update_detmatrix_1col!(wf.dnslater,
        #                        new_dn_col, dn_index, dn_detrat)

        return true
    else
        return false
    end
end

function check_and_update_gutzwiller!(wf::GutzwillerSlater)
   check_and_update_detmatrix!(wf.upslater)
   check_and_update_detmatrix!(wf.dnslater)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
import Base.+

function +(buffer::Array{String,1}, content::String;
    prefix::Union{String,Nothing}=nothing,suffix::Union{String,Nothing}=nothing)

    # create a new content line -
    new_line = content

    # prefix -
    if (prefix !== nothing)
        new_line = prefix*new_line
    end

    # suffix -
    if (suffix !== nothing)
        new_line = new_line*suffix
    end

    # cache -
    push!(buffer,new_line)
end

function +(buffer::Array{String,1}, content_array::Array{String,1}; padding::Union{String,Nothing}=nothing)
    for line in content_array

        if (padding === nothing)
            push!(buffer, line)
        else
            push!(buffer,"$(padding)$(line)")
        end
    end
end

function contains(string,token)
    return occursin(token,string)
end

function include_function(path_to_src_file::String)::Array{String,1}

    # create src_buffer -
    src_buffer::Array{String,1} = String[]

    # read -
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    # return the raw buffer -
    return src_buffer
end

function include_function(path_to_src_file::String, prefix_pad_string::String)::String

    # create src_buffer -
    src_buffer::Array{String,1} = String[]

    # read -
    open(path_to_src_file,"r") do src_file
        for line in eachline(src_file)

            new_line_with_line_ending = line*"\n"
            push!(src_buffer,new_line_with_line_ending)
        end
    end

    string_value = ""
    for line in src_buffer
        string_value *= prefix_pad_string*line
    end

    return src_buffer
end

function _request_user_input(prompt::String="")::String
    print(prompt)
    return chomp(readline())
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
using Documenter, MRIsim

makedocs(sitename="MRIsim.jl")

deploydocs(
    repo = "github.com/cncastillo/MRIsim.jl.git",
)

#makedocs(format = LaTeX())
#julia --color=yes make.jl


text:
using Actors
import Actors: newLink, spawn

inc!(a::Ref{Int}, c=0) = a[] += c

function dinc!(dict::Dict, key, by)
	count = get!(dict, key, 0)
	dict[key] = count + by
end
dinc!(dict) = copy(dict)

function link(lk, n)
	send(lk, n + 1)
	stop()
end

function chain(n::Int)
	start = newLink(1)
	lk = start
	for i in 1:n
		lk = spawn(link, lk)
	end
	send(lk, 0)
	receive(start)
end

##
## The following is a basic implementation of a
## task chain (without actor overhead)
##
counter(prev) = ch -> put!(prev, take!(ch) + 1)

function task_chain(n)
	start = Channel{Int}(1)
	ch = start
	for i in 1:n
		ch = Channel{Int}(counter(ch), spawn=true)
	end
	put!(ch, 0)
	take!(start)
end


text:
module Wikipedia

using HTTP, Gumbo, Cascadia
import Cascadia: matchFirst

include("Articles.jl")
using .Articles

const PROTOCOL = "https://"
const DOMAIN_NAME = "en.m.wikipedia.org"
const RANDOM_PAGE_URL = PROTOCOL * DOMAIN_NAME * "/wiki/Special:Random"

export fetchrandom, fetchpage, articleinfo, persistedarticle

# function fetchpage(url)
#   url = startswith(url, "/") ? buildurl(url) : url
#
#   response = HTTP.get(url)
#
#   if response.status == 200 && length(response.body) > 0
#     String(response.body)
#   else
#     ""
#   end
# end

function fetchpage(url)
  url = startswith(url, "/") ? buildurl(url) : url
  response = HTTP.get(url)
  content = if response.status == 200 && length(response.body) > 0
              String(response.body)
            else
              ""
            end
  relative_url = collect(eachmatch(r"/wiki/(.*)$", (response.request.parent == nothing ? url : Dict(response.request.parent.headers)["Location"])))[1].match

  content, relative_url
end

function extractlinks(elem)
  map(eachmatch(Selector("a[href^='/wiki/']:not(a[href*=':'])"), elem)) do e
    e.attributes["href"]
  end |> unique
end

function extracttitle(elem)
  matchFirst(Selector("#section_0"), elem) |> nodeText
end

function extractimage(elem)
  e = matchFirst(Selector(".content a.image img"), elem)
  isa(e, Nothing) ? "" : e.attributes["src"]
end

function extractcontent(elem)
  matchFirst(Selector("#bodyContent"), elem) |> string
end

function fetchrandom()
  fetchpage(RANDOM_PAGE_URL)
end

function articledom(content)
  if ! isempty(content)
    return Gumbo.parsehtml(content)
  end

  error("Article content can not be parsed into DOM")
end

# function articleinfo(content)
#   dom = articledom(content)
#
#   Dict( :content => content,
#         :links => extractlinks(dom.root),
#         :title => extracttitle(dom.root),
#         :image => extractimage(dom.root)
#   )
# end

# function articleinfo(content)
#   dom = articledom(content)
#   Article(content,
#           extractlinks(dom.root),
#           extracttitle(dom.root),
#           extractimage(dom.root))
# end

# function articleinfo(content)
#   dom = articledom(content)
#   (content, extractlinks(dom.root), extracttitle(dom.root), extractimage(dom.root))
# end

function articleinfo(content)
  dom = articledom(content)
  (extractcontent(dom.root), extractlinks(dom.root), extracttitle(dom.root), extractimage(dom.root))
end

function persistedarticle(article_content, url)
  article = Article(articleinfo(article_content)..., url)
  save(article)

  article
end

function buildurl(article_url)
  PROTOCOL * DOMAIN_NAME * article_url
end

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
mutable struct CWorldVis
    w::CWorld
    s::Union{Vec2, Nothing}
    f::Union{Function, Nothing}
    g::Union{AbstractGrid, Nothing}
    title::Union{String, Nothing}
end

function CWorldVis(w::CWorld;
                   s=nothing,
                   f=nothing,
                   g=nothing,
                   title=nothing)
    return CWorldVis(w, s, f, g, title)
end

@recipe function f(v::CWorldVis)
    xlims --> v.w.xlim
    ylims --> v.w.ylim
    aspect_ratio --> 1
    title --> something(v.title, "Continuum World")
    if v.f !== nothing
        @series begin
            f = v.f
            width = v.w.xlim[2]-v.w.xlim[1]
            height = v.w.ylim[2]-v.w.ylim[1]
            n = 200 # number of pixels
            nx = round(Int, sqrt(n^2*width/height))
            ny = round(Int, sqrt(n^2*height/width))
            xs = range(v.w.xlim[1], stop=v.w.xlim[2], length=nx)
            ys = range(v.w.ylim[1], stop=v.w.ylim[2], length=ny)
            zg = Array{Float64}(undef, nx, ny)
            for i in 1:nx
                for j in 1:ny
                    zg[j,i] = f(Vec2(xs[i], ys[j]))
                end
            end
            seriescolor --> cgrad([:red, :white, :green])
            seriestype := :heatmap
            xs, ys, zg
        end
    end
    if v.g !== nothing
        @series begin
            g = v.g
            xs = collect(ind2x(g, i)[1] for i in 1:length(g))
            ys = collect(ind2x(g, i)[2] for i in 1:length(g))
            label --> "Grid"
            marker --> :+
            markercolor --> :blue
            seriestype := :scatter
            xs, ys
        end
    end
end

Base.show(io::IO, m::MIME, v::CWorldVis) = show(io, m, plot(v))
Base.show(io::IO, m::MIME"text/plain", v::CWorldVis) = println(io, v)


text:
using GraphLasso, DataFrames, Distributions, RDatasets
using LinearAlgebra
using Test

@testset "Test Graph Lasso" begin
    S_iris = cov(Matrix(dataset("datasets", "iris")[:,1:4]))
    α = 1.0
    W_R = [1.68569351 0.000000 0.2743155 0.01969989
           0.00000000 1.189979 0.0000000 0.00000000
           0.27431545 0.000000 4.1162779 0.29560940
           0.01969989 0.000000 0.2956094 1.58100626]
    Θ_R = [ 0.59973155 0.0000000 -0.03996709  0.00000000
            0.00000000 0.8403507  0.00000000  0.00000000
           -0.03996708 0.0000000  0.24890787 -0.04604166
            0.00000000 0.0000000 -0.04604166  0.64111722]
    W_iris, Θ_iris = graphlasso(S_iris, α; tol=1e-6, penalize_diag=true)
    @test norm(W_R-W_iris) < 1e-5
    @test norm(Θ_R-Θ_iris) < 1e-5
end


text:
using LinearAlgebra, Printf, Arpack
function branchAndBound(prob, #problem object
		K; # target sparsity
		outputFlag = 3, # 1, 2, or 3 depending on level of detail sought in output
		timeCap = 3600, #maximum time in seconds
		warmStart = zeros(1), #warm start vector

		maxDepth = 1000, #maximum number of active nodes in the tree
		localSearchSteps = 2, # number of steps to take to find good feasible solutions at each node
		ksegmentsTrigger = 0, # use alternative branching when the number of indices left to set falls below this number
		dimSelectSteps = 10, # number of steps to take in choosing the dimension to branch on
		branchDimSelectMethod = 0,
		UOE=5000, # buffer to add when resizing arrays containing nodes and bounds
		gap = .000001, # optimality gap for the algorithm
		eigGap = .00001, # convergence criterion for determining eigenavlues by power method
		eigCycles = 500, # max number of iterations of the power method
		valtobeat = 0) # if non-zero, the algorithm will terminate once it has a solution with objective value larger than this

	# computes the vector b that maximizes b'Qb
	# subject to the support of b being <= abs(y)
	# terminates convergence early once the eigenvalue crosses refineCap (if not 0)
	# refineCap is used when generating multiple upper bounds for a subproblem
	# since we want to take the least upper bound, we terminate this function once any
	# more refinement would no longer impact the least upper bound
	function bbMyeigmax(y, refineCap)
		yKeep = .!(y.==0)
		thism=0
		if highDim
			thisA = copy(A[:, yKeep])
			thism, ~ = size(thisA)
			Q = thisA*thisA'/(thism-1)
			b = startingEig[yKeep]
			b = thisA*(b / LinearAlgebra.norm(b))
			normb = 0
			newnorm = LinearAlgebra.norm(b)
		else
			Q = copy(Sigma[yKeep,yKeep])
			beta0 = startingEig[yKeep]
			normb = LinearAlgebra.norm(beta0)
			b = Q*beta0
			newnorm = LinearAlgebra.norm(b)
		end

		cycles = 0
		while abs(normb - newnorm) > eigGap*normb && cycles < eigCycles
			normb = newnorm
			b = Q*(b / normb)
			newnorm = LinearAlgebra.norm(b)
			if refineCap>0
				if newnorm > refineCap
					return refineCap, copy(b/normb)
				end
			end
			cycles = cycles + 1
		end

		if highDim
			b = thisA'*b
			b = thisA'*(thisA*b/LinearAlgebra.norm(b))/(thism-1)
			normb = LinearAlgebra.norm(b)
		end

		b = b/normb
		expandedb = zeros(length(y))
		expandedb[yKeep]=b

		return normb, expandedb
	end

	# Returns the upper and lower bounds for the problem at node y
	# oldub is the upper bound of the parent node, which provides a maximum value
	# that the upper bound at y could have.  Feeding oldub into eigen_bound helps
	# reduce time spent in that function
	function return_bounds(y, oldub)
		if sum(max.(y,0)) == K
			val, ~ = bbMyeigmax(max.(y,0), 0)
			return [val val]
		elseif sum(abs.(y)) == K
			val, ~ = bbMyeigmax(abs.(y), 0)
			return [val val]
		else
			eb = eigen_bound(y, oldub)
			true_upper, ~ = bbMyeigmax(y, eb)
			lower_val = YuanSubset(y)[1]
			return [lower_val true_upper]
		end
	end

	# Computes two upper bounds for the problem at node y
	# The first is based on gershgorin circle theorem
	# The second is based on the trace
	# The function returns the lesser of the two bounds
	# oldub (the upper bound from the parent node) provides a maximum value on
	# the upper bound at y.  It is used here to terminate the gershgorin calculation early
	# if the gershgorin bound is higher than oldub
	function eigen_bound(y, oldub)
		ypositive = (y.==1)
		numpositive = sum(ypositive)
		stillneed = K-numpositive

		#Uses maximum absolute column sums to provide an upper bound on eigenvalues
		#Inspired by gershgorin circle theorem
		#Hard to scale well
		stillneed = K-numpositive
		startingsums = sum(absSigma[:, ypositive],dims=2)

		eb1=0
		cutoff = oldub*(1-1e-6)

		for i=1:length(y) # scanning over all columns
			if y[i]==-1 || y[i]==1 # which columns to consider
			    newsum = startingsums[i] # must include these rows
			    added = 0
			    if y[i]==-1
			    	# if we choose column i, we must choose row i
			    	newsum = newsum + absSigma[i,i]
			    	added = 1
			    end
			    j=1
			    while added < stillneed
			        candidateIndex = permMat[i,j]
			        if y[candidateIndex]==-1 && candidateIndex != i
			            newsum = newsum + absSigma[i,candidateIndex]
			            added = added + 1
			        end
			        j=j+1
			    end
			    if newsum>eb1
			        eb1=newsum
			        if newsum>cutoff
			    		break
			    	end
			    end
			end
		end

		#based on how the trace is a bound on the eigenvalues since they're all positive
		startingsums = sum(diagSigma[ypositive])
		indicesLeft = sortedOrder[y.==-1]
		eb2 = startingsums+sum(selectsorted(diagSigma[indicesLeft],stillneed))

		return min(eb1,eb2)
	end

	# Returns true if y represents a terminal node (only one k-sparse support is feasible)
	# Returns false otherwise
	function isTerminal(y)
		if sum(max.(y,0)) >= K return true end
		if sum(abs.(y)) <= K return true end
		return false
	end

	# Selects the dimension to branch on from node y
	# Finds the dimension among those with y_i==-1 that corresponds to the largest
	# loading of the "eigenvector" b that maximizes b'Qb subject to the support of b
	# being where y_i == -1.
	function branchDimSelect(y)
		#"dumb", "smart", "random", "semirandom"
		negindices = findall(s -> s == -1, y)

		if branchDimSelectMethod == 0  # adaptive

			# dimSelectLookup is a dictionary
			# It starts empty, but once dim is determined for some set of negindices
			# it is saved in dimSelectLookup, to prevent having to recompute it later
			dim = get(dimSelectLookup, negindices, 0)
			if dim > 0
				return dim
			else
				Q = Sigma[negindices,negindices]
				b = startingEig[negindices]
				for i = 1:dimSelectSteps
					b = Q*(b / LinearAlgebra.norm(b))
				end
				dim = negindices[findmax(abs.(b))[2]]
				dimSelectLookup[negindices] = dim
			end
		elseif branchDimSelectMethod == 1 # random
			dim = rand(negindices)
		elseif branchDimSelectMethod == 2 # sort by Q_ii
			dim = intersect(sortedOrder, negindices)[1]
		else # sort by loading on first eigenvector
			dim = intersect(firsteigOrder, negindices)[1]
		end
		return dim
	end

	# Applies the Yuan and Zhang first order heuristic for localSearchSteps steps
	# In order to generate a lower bound for the problem at y
	function YuanSubset(y)
		beta = Hk(startingEig, K, y)

			for i = 1:localSearchSteps
				beta = Hk(sparseTimes(Sigma,beta), K, y)
			end
			bestObj = (transpose(beta)*sparseTimes(Sigma,beta))[1]

		return bestObj, beta
	end

	# Selects which node to branch on
	# Depends primarily on the parameter maxDepth and the size of the tree
	# If the size of the tree is small, will select the node with the highest upper bound
	# If the tree is large, it wil pick the most recently created node (depth first search)
	function select_node()
		nn = sum(.!death[1:num_nodes])
		if nn == 0
			return 0
		end
		if nn < maxDepth
			return max_ub_select()
		else
			return depth_select()
		end
	end

	# Function to select the node with the highest upper bound
	# This will necessarily not be terminal as long as the algorithm hasn't converged
	function max_ub_select()
		fm = findmax(upper_bounds)
		if fm[1] > 0
			return fm[2]
		else
			return 0
		end
	end

	# Function to select the most recently created node that is not terminal
	function depth_select()
		selected_node = num_nodes
		while (isTerminal(nodes[:, selected_node]) | death[selected_node])
			selected_node = selected_node-1
			if selected_node == 0
				return 0
			end
		end
		return selected_node
	end

	# Beginning the algorithm
	start = time();

	Sigma = prob.Sigma
	absSigma = abs.(Sigma)
	diagSigma = LinearAlgebra.diag(Sigma)
	sortedOrder = sortperm(-diagSigma) #order from largest to smallest

	permMat = round.(Int64,zeros(size(absSigma)))
	for i=1:size(permMat)[2]
	    permMat[i,:]=sortperm(vec(absSigma[i,:]),rev=true)
	end

	dimSelectLookup = Dict()

	A = prob.data
	m, n = size(A)

	# detects if high dimenisonal problem, in which case
	# methods should use A'A not Sigma for eigenvalue computations
	highDim = false
	if m < n/2
		highDim = true
	end

	# Uses the Yuan algorithm to generate a warm start if none was provided
	if length(warmStart)==1
		~, warmStart = subset(prob, K, timeLimit = max(20,timeCap/100))
	end

	#initializing variables
	nodes = zeros(n, UOE)
	nodes[:,1] = round.(Integer,-1*ones(n,1))
	upper_bounds = zeros(UOE)

	eigGapSave = eigGap
	eigGap = .000001
	if highDim
		upper_bounds[1], startingEig = myeigmax(A, ones(n), ones(n), true, 0)
	else
		upper_bounds[1], startingEig = myeigmax(Sigma, ones(n), ones(n), false, 0)
	end

	firsteigOrder = sortperm(-abs.(startingEig))

	eigGap = eigGapSave
	upper = maximum(upper_bounds)
	lower = (transpose(warmStart)*Sigma*warmStart)[1]
	death = falses(UOE)
	num_nodes = 1
	bestBound = 0
	BestBoundtime = time()

	# if K=n problem is trivial
	if K == n
		lower = upper
	end

	explored = 0
	lower_revised = 0
	best_node = (warmStart.!=0)*1

	#Initializes output
	println(" Nodes,   Left,  Objective,  Incumbent,       Gap(%),   Runtime(s)")
	toPrint=[Printf.@sprintf("%6d, %6d, %10f, %10f, %10.3f %%, %10.3f s \n", num_nodes, explored, upper, lower, (upper-lower)/(1e-10+upper)*100, time()-start)]
	printtime = time()

	while (upper - lower)/upper > gap  && time()-start < timeCap
		explored = explored + 1

		# Occasionally print updates and risize arrays
		if lower_revised == 1 || size(nodes)[2]-num_nodes < K+5 || explored < 100
			nodesToKeep = .!(upper_bounds.< lower) .& .!death
			nodesToKeep = findall(s->s,nodesToKeep)
			nodesToKeep = nodesToKeep[nodesToKeep.< num_nodes+1]
			upper_bounds = [upper_bounds[nodesToKeep]; zeros(UOE)]
			nodes = [nodes[:, nodesToKeep] zeros(n,UOE)]
			death = [death[nodesToKeep]; falses(UOE)]
			num_nodes = length(nodesToKeep)
		end

		upper = maximum(upper_bounds[.!death]);

		if lower_revised == 1 || explored < 100 || time() - printtime > 1
			printtime = time()
			if outputFlag == 1
				#Printf.@printf("%.4e, %.4e, %.4e \n", time()-start, upper, lower)
			elseif outputFlag == 2
				Printf.@printf("upper:  %.4e  \t lower: %.4e \n", upper, lower)
				Printf.@printf("%d nodes left  \t %d in \t %.4e seconds \n", num_nodes, explored, time()-start)
			elseif outputFlag == 3 && num_nodes>0
				s=Printf.@sprintf("%6d, %6d, %10f, %10f, %10.3f %%, %10.3f s\n", num_nodes, explored, upper, lower, (upper-lower)/(1e-10+upper)*100, time()-start)
				print(s)
				push!(toPrint, s)
			end
		end

		#selecting node to split on
		selected_node = select_node()
		if selected_node == 0
			break;
		end
		old_y = nodes[:, selected_node]

		# Determines whether to split in 2 on one dimension or use multiple segments
		if (K-sum(max.(old_y,0)) < ksegmentsTrigger) & (randn()>1)
			numOnes = sum(max.(old_y,0))
			negOnes = findall(s->s==-1, old_y)
			if length(negOnes) >= (K-numOnes+1)*2
				numBranches = round.(Integer,K-numOnes+1)
				newNodes = transpose(repeat(transpose(old_y), numBranches))
				partition = yourpart(negOnes,numBranches)
				for i = 1:numBranches
					newNodes[partition[i],i] = 0
				end
			else
				numBranches = 2
				newNodes = transpose(repeat(transpose(old_y), numBranches))
				branch_dim = branchDimSelect(old_y)
				newNodes[branch_dim ,1] = 1
				newNodes[branch_dim ,2] = 0
			end
		else
			numBranches = 2
			newNodes = transpose(repeat(transpose(old_y), numBranches))
			branch_dim = branchDimSelect(old_y)
			newNodes[branch_dim ,1] = 1
			newNodes[branch_dim ,2] = 0
		end

		# Computes bounds at new nodes and saves them if not dominated or terminal
		lower_revised = 0
		oldub = upper_bounds[selected_node]
		for i = 1:numBranches
			lb, ub = return_bounds(newNodes[:,i], oldub)
			if ub*(1-gap) > lower
				if lb > lower
					lower = lb
					best_node = copy(newNodes[:,i])
					lower_revised = 1
				end
				if .!(isTerminal(newNodes[:,i]))
					num_nodes = num_nodes + 1
					nodes[:,num_nodes] = copy(newNodes[:,i])
					upper_bounds[num_nodes] = ub
				end
			end
		end

		# removes old node
		death[selected_node] = true
		upper_bounds[selected_node] = 0

		# updates best feasible solution
		if lower > bestBound
			bestBound = lower
			BestBoundtime = time()
		end

		# terminates early if the goal was to find a superior solution than some candidate valtobeat
		if (valtobeat>0) & (lower > valtobeat)
			break;
		end

	end

	# records time taken
	timeToBestBound = BestBoundtime - start

	s=Printf.@sprintf("%6d, %6d, %10f, %10f, %10.3f %%, %10.3f s \n", num_nodes, explored, upper, lower, (upper-lower)/(1e-10+upper)*100, time()-start)
	push!(toPrint, s)

	# records whether the algorithm timed out
	timeOut = false
	if time()-start > timeCap
		timeOut = true
	end

	# saves final solutions
	y = (YuanSubset(best_node)[2].!=0)*1

	yKeep = .!(y.==0)
	Q = copy(Sigma[yKeep,yKeep])
	lambdas=Float64[]
	betas=Float64[]
	if size(Q,1)>1
		lambdas, betas, =Arpack.eigs(Q, which=:LR, nev=1)
	else
		lambdas=[1.0]
		betas=[1.0]
	end
	eig_soln = betas[:,1]
	obj = lambdas[1]
	xVal = zeros(length(y))
	xVal[yKeep] = eig_soln

	final_gap = max(0,(upper-lower)/(1e-10+upper)*100)

	return obj, # objective value
			xVal, # best feasible solution
			timeToBestBound, # time elapsed from the state to the identification of the best feasible solution
			time()-start, # total time
			timeOut, # whether the algorithm timed out
			explored, # how many nodes were searched
			toPrint, #complete output
			final_gap # final optimality gap

end


text:
# Cross validation

## cross validation generators

abstract CrossValGenerator

# K-fold

immutable Kfold <: CrossValGenerator
    permseq::Vector{Int}
    k::Int
    coeff::Float64

    function Kfold(n::Int, k::Int)
        2 <= k <= n || error("The value of k must be in [2, length(a)].")
        new(randperm(n), k, n / k)
    end
end

length(c::Kfold) = c.k

immutable KfoldState
    i::Int      # the i-th of the subset
    s::Int      # starting index
    e::Int      # ending index
end

start(c::Kfold) = KfoldState(1, 1, iround(c.coeff))
next(c::Kfold, s::KfoldState) =
    (i = s.i+1; (sort!(c.permseq[s.s:s.e]), KfoldState(i, s.e+1, iround(c.coeff * i))))
done(c::Kfold, s::KfoldState) = (s.i > c.k)

# LOOCV (Leave-one-out cross-validation)

function leave_one_out(n::Int, i::Int)
    @assert 1 <= i <= n
    x = Array(Int, n-1)
    for j = 1:i-1
        x[j] = j
    end
    for j = i+1:n
        x[j-1] = j
    end
    return x
end

immutable LOOCV <: CrossValGenerator
    n::Int
end

length(c::LOOCV) = c.n

start(c::LOOCV) = 1
next(c::LOOCV, s::Int) = (leave_one_out(c.n, s), s+1)
done(c::LOOCV, s::Int) = (s > c.n)

# Repeated random sub-sampling

immutable RandomSub <: CrossValGenerator
    n::Int    # total length
    sn::Int   # length of each subset
    k::Int    # number of subsets
end

length(c::RandomSub) = c.k

start(c::RandomSub) = 1
next(c::RandomSub, s::Int) = (sort!(sample(1:c.n, c.sn; replace=false)), s+1)
done(c::RandomSub, s::Int) = (s > c.k)

## Cross validation algorithm
#
#  estfun: model estimation function
#
#          model = estfun(train_inds)
#
#          it takes as input the indices of
#          the samples for training, and returns
#          a trained model.
#
#  evalfun: model evaluation function
#
#           v = evalfun(model, test_inds)
#
#           it applies a trained model to a subset
#           of samples whose indices are given in
#           test_inds, and returns an overall score.
#
#  n:   the number of samples in the whole data set
#
#  gen: an iterable object (e.g. an instance of CrossValGenerator)
#       where each element is a vector of indices
#
#  ord: ordering of the score
#
#  This function returns a tuple (best_model, best_score, best_traininds)
#
function cross_validate(estfun::Function, evalfun::Function, n::Integer, gen, ord::Ordering)
    best_model = nothing
    best_score = NaN
    best_inds = Int[]
    first = true

    for test_inds in gen
        train_inds = setdiff(1:n, test_inds)
        model = estfun(train_inds)
        score = evalfun(model, test_inds)
        if first || lt(ord, best_score, score)
            best_model = model
            best_score = score
            best_inds = train_inds
            first = false
        end
    end

    return (best_model, best_score, best_inds)
end

cross_validate(estfun::Function, evalfun::Function, n::Integer, gen) =
    cross_validate(estfun, evalfun, n, gen, Forward)


text:
_open_video(open_func, file; gray) = VsStream(
    video = if gray
        open_func(file, target_format=VideoIO.AV_PIX_FMT_GRAY8)
    else
        open_func(file)
    end
)

open_video(file; gray = false) = _open_video(VideoIO.openvideo, file; gray = gray)

open_camera(
    device = VideoIO.DEFAULT_CAMERA_DEVICE[];
    gray = false
) = _open_video(VideoIO.opencamera, device; gray = gray)

read_frame(stream::VsStream) = VsFrame(
    image = read(stream.video),
    time = position(stream.video)
)

function load_image(file; gray = false)
    img = load(file)
    if gray
        img = to_gray_image(img)
    end
    img
end

const LITERAL_TYPES = Union{Integer, AbstractFloat, Bool, Dates.DateTime, Dates.Time, Dates.Date, Symbol}

enum_prefix(::Type{T}) where T <: Enum = ""
function enum_to_string(x::T) where T <: Enum
    prefix = enum_prefix(T)
    s = string(x)
    if prefix == ""
        s
    else
        split(s, prefix, limit = 2)[2]
    end |> lowercase
end
function enum_from_string(s::String, ::Type{T}) where T <: Enum
    values = instances(T)
    s = lowercase(enum_prefix(T)) * to_snake_case(s)
    values[findfirst(values) do x
        lowercase(string(x)) == s
    end]
end

_to_generator(d::AbstractDict) = (Symbol(key) => value for (key, value) in d)
function to_snake_case(s::AbstractString)
    s = replace(s, r"\W+" => '_', r"([a-z])([A-Z])" => s"\1_\2", '-' => '_')
    lowercase(s)
end
function to_kebab_case(s::AbstractString)
    s = replace(s, r"\W+" => '-', r"([a-z])([A-Z])" => s"\1-\2", '_' => '-')
    lowercase(s)
end

serialize(object) = VsYAML.yaml(object)

_parse(val, ::Type{Any}; kwargs...) = val
_parse(val, T::Type; kwargs...) = convert(T, val)
_parse(int::Integer, ::Type{T}; kwargs...) where T <: Enum = T(int)
_parse(s, ::Type{T}; kwargs...) where T <: Enum = enum_from_string(string(s), T)
_parse(arr::AbstractArray, ::Type{<:AbstractArray{T}}; kwargs...) where T = [_parse(x, T; kwargs...) for x in arr]
_parse(arr::AbstractArray, ::Type{T}; kwargs...) where T <: Tuple =
    tuple(_parse(x, eltype(TE); kwargs...) for (x, TE) in zip(arr, T.types))
_parse(dict::AbstractDict, T::Type{<:AbstractDict{TKey, TValue}}; kwargs...) where {TKey, TValue} = T(
    _parse(key, TKey; kwargs...) => _parse(value, TValue; kwargs...)
    for (key, value) in dict
)
function _parse(dict::AbstractDict, T::Type; other_key = nothing, kwargs...)
    params = Dict{Symbol, Any}()
    fnames = fieldnames(T)
    other_key = if !isnothing(other_key)
        other_key = Symbol(other_key)
        if other_key ∈ fnames
            params[other_key] = Dict{Symbol, Any}()
            other_key
        end
    end
    for (key, value) in dict
        skey = to_snake_case(string(key))
        skey = Symbol(skey)
        if skey ∉ fnames
            if !isnothing(other_key)
                params[other_key][skey] = value
            end
            continue
        end
        TF = fieldtype(T, skey)
        params[skey] = _parse(value, TF; other_key = other_key, kwargs...)
    end
    T(; _to_generator(params)...)
end

"""
    deserialize(yaml, T)

Note that `T` must be a type that can be constructed by keywords, e.g., defined by `Base.@kwdef`
"""
function deserialize(yaml, T::Type; dicttype = Dict, other_key = nothing, kwargs...)
    dict = YAML.load(yaml; dicttype = dicttype)
    _parse(dict, T; other_key = other_key, kwargs...)
end


text:
module cost_functions

using code_structs
using channels

export hamming_cost
export exact_cost
export euclidean_cost

function hamming_cost(p::Float64, curr_seq::Array{Float64,1}, node::TreeNode, transition::Transition)
	return node.cost + sum(abs.(transition.output - curr_seq))
end

function exact_cost(p::Float64, curr_seq::Array{Float64,1}, node::TreeNode, transition::Transition)
	difference = sum(abs.(transition.output - curr_seq))
	return node.cost - (difference * log(p) + (length(transition.output) - difference) * log(1-p))
end

function euclidean_cost(p::Float64, curr_seq::Array{Float64,1}, node::TreeNode, transition::Transition)
	return node.cost + sum((BPSK(transition.output) - curr_seq).^2)
end

end


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@reexport using DataFrames

export sourceRand

export makeRand
function makeRand(m :: Model)
    fpre = sourceRand(m)
    f(;kwargs...) = Base.invokelatest(fpre; kwargs...)
end

export rand
rand(m::Model; kwargs...) = makeRand(m)(;kwargs...)

function rand(m::Model, n::Int64; kwargs...)
    r = makeRand(m)
    [r(;kwargs...) for j in 1:n] # |> DataFrame

end

struct FuncRand{AST} end

@generated function (::FuncRand{AST})(args...; kwargs...) where AST
    interpret(AST)
end

function sourceRand(m::Model)
    m = canonical(m)
    proc(m, st::Let)     = :($(st.x) = $(st.rhs))
    proc(m, st::Follows) = :($(st.x) = rand($(st.rhs)))
    proc(m, st::Return)  = :(return $(st.rhs))
    proc(m, st::LineNumber) = nothing

    body = buildSource(m, proc) |> striplines

    argsExpr = Expr(:tuple,freeVariables(m)...)

    stochExpr = begin
        vals = map(variables(m)) do x Expr(:(=), x,x) end
        Expr(:tuple, vals...)
    end

    ast = flatten(@q (
        begin
            @unpack $argsExpr = kwargs
            $body
            $stochExpr
        end
    ))

    FuncRand{expr2typelevel(ast)}()

end


text:
"""
Load estimates from known DGPs and check that the parameter uncertainty is estimated sufficiently well
"""

#region import and models

#%%
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"

using DrWatson
using DataFrames
using ScoreDrivenERGM
import ScoreDrivenERGM:DynNets, Utilities
using PyPlot
pygui(true)
using Statistics
using Logging
using SharedArrays
#endregion

import ScoreDrivenERGM:DynNets

@elapsed df = collect_results( datadir("sims", "dgp&FIl_est"))
@elapsed df = collect_results( datadir("sims", "dgp&FIl_conf"))
df["modelTag"] = string.(df["model"])

# select results for one set of parameters
begin
T = 300
N = 100
nSample = 50
indQuant = 1
modelTag = string(DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType = "FISH_D"))
modelTag = string(DynNets.SdErgmDirBin0Rec0_pmle(scoreScalingType = "FISH_D"))
parUncMethod = "WHITE-MLE"
parUncMethod = "NPB-MVN"

dgpSetting = DynNets.list_example_dgp_settings(DynNets.SdErgmDirBin0Rec0_mle()).dgpSetARlowlow

# # visual inspection of filter, dgp and conf bands
res = filter([:modelTag, :T, :N, :dgp, :S, :m] => (mTag,t,n, d, s, m) -> all((mTag==modelTag, t==T, n==N, d == dgpSetting, s == nSample, m == parUncMethod)), df)

nrow(res) != 1 ? error("$N, $T,  $(size(res))") : res = res[1,:]
    for nPlot = 1:10
        if !any(res.errInds[:,nPlot])
            DynNets.plot_filtered_and_conf_bands(res.model, N, res.allfVecT_filt[:,1:end,nPlot], res.allConfBandsFiltPar[:,1:end,:,:,nPlot] ; parDgpTIn=res.allParDgpT[:,1:T,nPlot], confBands2In = res.allConfBandsPar[:,1:end,:,:,nPlot] , offset = 1, indBand = indQuant)
        end
    end

using Statistics
using SharedArrays
# unrestrict static SD estimates

if true

names(res)

allCovWhite = res.allmvSDUnParEstCovWhite
covParBoot = cov(vEstSdUnParBootDist')
covParBoot = cov(Utilities.drop_bad_un_estimates(vEstSdUnParBootDist)')

parNames = [ "B_theta", "A_theta",  "B_eta", "A_eta"]
fig, ax = subplots(2,2)
for i = 3:6
    j=i
    ax[i-2].hist(log10.(allCovWhite[i,j,:]), range=quantile(log10.(allCovWhite[i,j,:]), [0.05, 0.95]), 20)
    ylim = ax[i-2].get_ylim()
    ax[i-2].vlines(log10(covParBoot[i,j]), ylim[1], ylim[2], color = "r" )
    ax[i-2].set_title("Var( $(parNames[i-2]))")
end
suptitle("Empirical Distribution of White Estimators \n $(res.modelTag) \n  nSample = $nSample , T = $T, N = $N, A=$(res.dgp.opt.A[1])")
tight_layout()

end

res = filter([:modelTag, :T, :N, :dgpSettings, :nSample] => (mTag,t,n, d, s) -> all((mTag==modelTag, t==T, n==N, d == dgpSetting, s == nSample)), df)

nrow(res) != 1 ? error("$N, $T,  $(size(res))") : res = res[1,:]

vEstSdResParBootDist =  res.allvEstSdResPar
vEstSdUnParBootDist = mapslices(x -> DynNets.unrestrict_all_par(res.model, x), vEstSdResParBootDist, dims=1)

parNames = ["w_theta", "B_theta", "A_theta", "w_eta", "B_eta", "A_eta"]
if false
sample1 = vEstSdResParBootDist
sample1 = vEstSdUnParBootDist
sample1 = Utilities.drop_bad_un_estimates(vEstSdUnParBootDist)
fig, ax = subplots(3,2)
for i = 1:6
    ax[i].hist(sample1[i,:], range=quantile(sample1[i,:], [0.01, 0.99]), 20, alpha = 0.4)
    ax[i].set_title(parNames[i])
    ylim = ax[2,1].get_ylim()
    # ax[i].vlines(res.vEstSdResPar[i], ylim[1], ylim[2], color = "r" )

end
end
if false
fig, ax = subplots(2,2)
ax[1,1].hist(res.allvEstSdResPar[2, :], 30)
ylim = ax[1,1].get_ylim()
ax[1,1].vlines(res.dgpSettings.opt.B, ylim[1], ylim[2], color = "r" )
ax[2,1].hist(res.allvEstSdResPar[5, :], 30)
ylim = ax[2,1].get_ylim()
ax[2,1].vlines(res.dgpSettings.opt.B, ylim[1], ylim[2], color = "r" )
ax[1,2].hist(res.allvEstSdResPar[3, :], 30)
ylim = ax[1,2].get_ylim()
ax[1,2].vlines(res.dgpSettings.opt.A, ylim[1], ylim[2], color = "r" )
ax[2,2].hist(res.allvEstSdResPar[6, :], 30)
ax[2,2].vlines(res.dgpSettings.opt.A, ylim[1], ylim[2], color = "r" )
ylim = ax[2,2].get_ylim()
# ax[2,2].vlines(res.dgpSettings.opt.A, ylim[1], ylim[2], color = "r" )
suptitle("Estimators' distributions for Static SD parameters  nSample = $nSample , T = $T, N = $N")
tight_layout()

end
end

# Domande
# 1. C'é un bias di campione finito. Lo stimatore, MLE o White, della covarianza degli stimatori dei parametri é ok?
# 2. Visto che posso costruire la distribuzione dei parametri con parametric bootstrap, perché non usare quella?
# 3. Che relazione c'é tra la distribuzione da parametric bootstrap e la mia idea di usare non parametric bootstrap campionando le osservazioni (che portano in dote la storia fino a quel punto, ma considerando la likelihood della sola osservazione)?

#endregion

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
module datasets
import ..datasets_dir
include(joinpath(datasets_dir, "voc2007", "voc.jl"))
end #module

using BinaryProvider
"""
    download_dataset(;name="all")

Download supported datasets. If folder exists, deletes first.
"""
function download_dataset(name::String = "all")
    !isdir(datasets_dir) && mkdir(datasets_dir)
    if any(name .== ["voc2007","VOC2007","voc", "VOC", "all"])
        voc_dir = joinpath(datasets_dir, "voc2007")
        voc_root = joinpath(voc_dir, "VOCdevkit")
        isdir(voc_root) && rm(voc_root, force = true, recursive = true)
        @info "Downloading dataset..."
        tmploc = joinpath(voc_dir,"VOCtrainval_06-Nov-2007.tar.gz")
        download("https://github.com/ianshmean/YOLO.jl/releases/download/VOC2007/VOCtrainval_06-Nov-2007.tar.gz",tmploc)
        @info "Extracting..."
        BinaryProvider.unpack(tmploc, voc_dir)
        rm(tmploc, force = true)
    end
end