text:
#########################################################################
#
#   Bi-directional Dict()
#
#     - enforces unicity of keys and values
#     - provides lookup by values
#
#########################################################################

import Base: setindex!, getindex, haskey, delete!,
  keys, values, start, next, done

type BiDict{K,V}
  kv::Dict{K,V}
  vk::Dict{V,K}

  BiDict() = new(Dict{K,V}(), Dict{V,K}())

  function BiDict(ks, vs)
    n = length(ks)
    length(unique(ks)) != n && error("Duplicate keys")
    length(unique(vs)) != n && error("Duplicate values")
    h = BiDict{K,V}()
    for i=1:n
      h.kv[ks[i]] = vs[i]
      h.vk[vs[i]] = ks[i]
    end
    return h
  end

  function BiDict(d)
    n = length(d)
    vs = values(d)
    length(unique(vs)) != n && error("Duplicate values")
    h = BiDict{K,V}()
    for (k,v) in d
      h.kv[k] = v
      h.vk[v] = k
    end
    return h
  end
end

# BiDict()                                                = BiDict{Any,Any}()
# BiDict{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V}) = BiDict{K,V}(ks,vs)
# BiDict(ks, vs)                                          = BiDict{Any,Any}(ks, vs)
# BiDict{K,V}(d::Dict{K,V})                               = BiDict{K,V}(d)
# BiDict(d)                                               = BiDict{Any,Any}(d)

# use first dict (kv) by default for base functions
function setindex!(bd::BiDict, v, k)
  k2 = get(bd.vk, v, nothing)  # existing key for v ?
  if k2 != nothing
    delete!(bd.kv, k2)
    delete!(bd.vk, v)
  end

  v2 = get(bd.kv, k, nothing)  # existing value for k ?
  if v2 != nothing
    delete!(bd.kv, k)
    delete!(bd.vk, v2)
  end

  bd.kv[k] = v
  bd.vk[v] = k
end

function delete!(bd::BiDict, k)
  !haskey(bd.kv, k) && error("unknown key $k")
  v = bd.kv[k]
  delete!(bd.kv, k)
  delete!(bd.vk, v)
end

length(bd::BiDict)      = length(bd.kv)
haskey(bd::BiDict, k)   = haskey(bd.kv, k)
keys(bd::BiDict)        = keys(bd.kv)
values(bd::BiDict)      = values(bd.kv)

start(bd::BiDict)       = start(bd.kv)
next(bd::BiDict, i)     = next(bd.kv, i)
done(bd::BiDict, i)     = done(bd.kv, i)

getindex(bd::BiDict, k) = bd.kv[k]


text:
module SimradEK60

using SimradRaw

export Sv, TS, pings, power, powerdb,  athwartshipangle, alongshipangle, R, filetime, power2db, octet2deg

# Float32 is okay for most practical purposes and within 10^-5 of EchoView
# Float64 gets to 10^-7 of Echoview
FLOAT_TYPE = Float32

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

# We split the RAW file into a series of EK60 specific pings.

# Whilst it is tempting to include latitude and longitude in the ping,
# this requires NMEA parsing and decisions about interpolation. Some
# people use external GPSs. Since the processing approach could be a
# user preference and potential source of errors, we choose not to
# include them. A caller can easily access NMEA datagrams produced by
# SimradRaw.jl and call PyNmea.

struct EK60Ping
    power::Vector{Int16}
    athwartshipangle::Vector{Int8}
    alongshipangle::Vector{Int8}
    filetime::UInt64
    frequency::FLOAT_TYPE
    soundvelocity::FLOAT_TYPE
    sampleinterval::FLOAT_TYPE
    absorptioncoefficient::FLOAT_TYPE
    transmitpower::FLOAT_TYPE
    pulselength::FLOAT_TYPE
    gain::FLOAT_TYPE
    equivalentbeamangle::FLOAT_TYPE
    sacorrection::FLOAT_TYPE
    sacorrectiontable::Vector{FLOAT_TYPE}
    pulselengthtable::Vector{FLOAT_TYPE}
end

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

function load(filename::AbstractString)
    ps = collect(pings(filename))
end

"""
    R(r, s, T)

TVG range correction for Ex60

R is the corrected range (m)

r the uncorrected range (m)

s is the TvgRangeCorrectionOffset value, See Echoview documentation,
Simrad Time Varied Gain (TVG) range correction, http://bit.ly/2pqzS2D

T is the sample thickness (m)
"""
R(r, s, T) = max.(0f0, r .- s*T)

"""
    R(ping::EK60Ping; soundvelocity = nothing, rangecorrectionoffset=2)

Returns the corrected range (depth) of samples in `ping`.
"""
function R(ping::EK60Ping; soundvelocity = nothing, rangecorrectionoffset=2)

    if soundvelocity == nothing
        soundvelocity = ping.soundvelocity
    end

    p = ping.power

    samplethickness = soundvelocity .* ping.sampleinterval / 2 # in metres of range

    r = [x* samplethickness for x in 0:length(p)-1]
    R(r, rangecorrectionoffset, samplethickness)

end

function R(pings::Vector{EK60Ping}; soundvelocity = nothing, rangecorrectionoffset=2)

    r = [R(ping,
           soundvelocity=soundvelocity,
           rangecorrectionoffset=rangecorrectionoffset) for ping in pings]

    myhcat(r)

end

R(pings::Channel{EK60Ping}; soundvelocity = nothing, rangecorrectionoffset=2) =
    R(collect(pings),
      soundvelocity=soundvelocity,
      rangecorrectionoffset=rangecorrectionoffset)

###

"""
    filetime(ping::EK60Ping)

Returns the FILETIME timestamp for a ping
"""
filetime(ping::EK60Ping) = ping.filetime

function filetime(pings::Vector{EK60Ping})
    [ping.filetime for ping in pings]
end

filetime(pings::Channel{EK60Ping}) = filetime(collect(pings))

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

# We must apply the SONAR equation and instrument corrections as
# described in http://bit.ly/2o1oOrq

## % An earlier MATLAB implementation
## % Sv = recvPower + 20 log10(Range) + (2 *  alpha * Range) - (10 * ...
## %           log10((xmitPower * (10^(gain/10))^2 * lambda^2 * ...
## %           c * tau * 10^(psi/10)) / (32 * pi^2)) - (2 * SaCorrection)

"""
    Sv(Pr, λ, G, Ψ, c, α, Pt, τ, Sa, R)
where:

R = the corrected range (m) - see TVG Range Correction

Pr = received power (dB re 1 W) - see Simrad EK numbers to Power

Pt = transmitted power (W)

α = absorption coefficient (dB/m)

G0 = transducer peak gain (non-dimensional) calculated as 10^(G/10)

G is the Transducer gain (dB re 1)

λ = wavelength (m) = c/f

f = frequency (Hz)

c = sound speed (m/s)

τ = transmit pulse duration (s) - also known as the
TransmittedPulseLength

ψ = Equivalent Two-way beam angle (Steradians) calculated as 10^(Ψ/10)

Ψ is the Two-way beam angle (dB re 1 Steradian)

Sa = Simrad correction factor (dB re 1m−1) determined during
calibration. This represents the correction required to the Sv
constant to harmonize the TS and NASC measurements.

"""
function Sv(Pr, λ, G, Ψ, c, α, Pt, τ, Sa, R)

    # TVG is applied to samples with ranges greater than 1 meter.
    tvg =  20log10.(max.(1, R))

    csv = 10log10.((Pt * (10^(G/10))^2 *  λ^2 * c * τ * 10^(Ψ/10)) /
                   (32 * FLOAT_TYPE(pi)^2))

    # Ignore absorption for R < 1m
    r =  [(x < 1) ? 0 : x for x in R]

    Pr .+ tvg + (2 * α * r) .- csv .- 2Sa
end

"""
     TS(Pr, λ, G, α, Pt, R)

Target strength

Pr = received power (dB re 1 W) - see Simrad EK numbers to Power

G is the Transducer gain (dB re 1)

α = absorption coefficient (dB/m)

Pt = transmitted power (W)

R = the corrected range (m) - see TVG Range Correction

"""
function TS(Pr, λ, G, α, Pt, R)

    # TVG is applied to samples with ranges greater than 1 meter.
    tvg =  40log10.(max.(1, R))

    csv = 10log10.((Pt * (10^(G/10))^2 *  λ^2) /
                   (16 * FLOAT_TYPE(pi)^2))

    # Ignore absorption for R < 1m
    r =  [(x < 1) ? 0 : x for x in R]

    Pr .+ tvg + (2 * α * r) .- csv
end

"""
    pings(datagrams::Channel{SimradRaw.Datagram})

Returns a `Channel` of `Ping`s from a `Channel` of `Datagram`s.

"""
function pings(datagrams::Channel{SimradRaw.Datagram})

    function _it(chn1)
        config = nothing
        for datagram in datagrams
            if datagram.dgheader.datagramtype == "RAW0"

                transducer = config.configurationtransducers[datagram.channel]

                if transducer.frequency != datagram.frequency
                    @error("$(config.configurationtransducer)  $datagram")
                end

                idx = findfirst(datagram.pulselength .== transducer.pulselengthtable)

                sacorrection = transducer.sacorrectiontable[idx]
                gain = transducer.gaintable[idx]

                ping = EK60Ping(datagram.power,
                                datagram.athwartshipangle,
                                datagram.alongshipangle,
                                SimradRaw.filetime(datagram.dgheader.datetime),
                                datagram.frequency,
                                datagram.soundvelocity,
                                datagram.sampleinterval,
                                datagram.absorptioncoefficient,
                                datagram.transmitpower,
                                datagram.pulselength,
                                gain,
                                transducer.equivalentbeamangle,
                                sacorrection,
                                transducer.sacorrectiontable,
                                transducer.pulselengthtable)

                put!(chn1, ping)

            elseif datagram.dgheader.datagramtype == "CON0"
                config = datagram
            end
        end
    end

    return Channel(_it, ctype=EK60Ping)
end

function pings(datagrams::Vector{SimradRaw.Datagram})
    # Is there a neater method here?
    chnl = Channel(c->foreach(i->put!(c,i), datagrams), ctype=SimradRaw.Datagram)
    pings(chnl)
end

"""
    pings(filename::AbstractString)

Returns a `Channel` over `Ping`s from the RAW file designated by
`filename`.

"""
function pings(filename::AbstractString)
    pings([filename])
end

"""
    pings(filenames::Vector{AbstractString}})

Returns a `Channel` of `Pings` from the RAW files designated by
`filenames`.

"""
function pings(filenames::Vector{String})
    pings(datagrams(filenames))

end

"""
    Sv(ping::EK60Ping;
            frequency=nothing,
            gain=nothing,
            equivalentbeamangle=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            pulselength=nothing,
            sacorrection=nothing,
            rangecorrectionoffset=2)

Returns a `Vector` of Sv, the (Mean) Volume backscattering strength (MVBS) in (dB re
1 m-1) for a given `ping`.

The function accepts a number of optional arguments which, if
specified, override the ping's own settings. This facilitates external
calibration.

"""
function Sv(ping::EK60Ping;
            frequency=nothing,
            gain=nothing,
            equivalentbeamangle=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            pulselength=nothing,
            sacorrection=nothing,
            rangecorrectionoffset=2)

    if frequency == nothing
        frequency = ping.frequency
    end

    if gain == nothing
        gain = ping.gain
    end

    if equivalentbeamangle == nothing
        equivalentbeamangle = ping.equivalentbeamangle
    end

    if soundvelocity == nothing
        soundvelocity = ping.soundvelocity
    end

    if absorptioncoefficient == nothing
        absorptioncoefficient = ping.absorptioncoefficient
    end

    if transmitpower == nothing
        transmitpower= ping.transmitpower
    end

    if pulselength == nothing
        pulselength = ping.pulselength
    end

    if sacorrection == nothing
        sacorrection = ping.sacorrection
    end

    pdb = powerdb(ping)

    rangecorrected = R(ping,
                       soundvelocity = soundvelocity,
                       rangecorrectionoffset=rangecorrectionoffset)

    λ =  soundvelocity / frequency # calculate wavelength

    Sv(pdb, λ, gain, equivalentbeamangle,
            soundvelocity, absorptioncoefficient,
            transmitpower, pulselength, sacorrection,
            rangecorrected)

end

"""
    TS(ping::EK60Ping;
            frequency=nothing,
            gain=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            rangecorrectionoffset=0)

Target strength

"""
function TS(ping::EK60Ping;
            frequency=nothing,
            gain=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            rangecorrectionoffset=0)

    if frequency == nothing
        frequency = ping.frequency
    end

    if gain == nothing
        gain = ping.gain
    end

    if soundvelocity == nothing
        soundvelocity = ping.soundvelocity
    end

    if absorptioncoefficient == nothing
        absorptioncoefficient = ping.absorptioncoefficient
    end

    if transmitpower == nothing
        transmitpower= ping.transmitpower
    end

    pdb = powerdb(ping)

    rangecorrected = R(ping,
                       soundvelocity = soundvelocity,
                       rangecorrectionoffset = rangecorrectionoffset)

    λ =  soundvelocity / frequency # calculate wavelength

    TS(pdb, λ, gain, absorptioncoefficient, transmitpower, rangecorrected)

end

function myhcat(s; missingvalue=NaN32)
    ls = map(length,s)
    minlength =minimum(ls)
    maxlength =maximum(ls)
    if minlength == maxlength
        return hcat(s...)
    end

    etype = typeof(missingvalue)

    array = Array{etype}(undef, maxlength, length(s))

    array .= missingvalue
    for i in 1:length(s)
        r = 1:length(s[i])
        array[r, i] .= s[i][r]
    end
    return array
end

"""
    Sv(pings::Vector{EK60Ping};
        frequency=nothing,
        gain=nothing,
        equivalentbeamangle=nothing,
        soundvelocity=nothing,
        absorptioncoefficient=nothing,
        transmitpower=nothing,
        pulselength=nothing,
        sacorrection=nothing,
        rangecorrectionoffset=2)

Returns an `Array` of Sv, the Volume backscattering strength in (dB re
1 m-1) for a set of contiguous `pings`.

The function accepts a number of optional arguments which, if
specified, override the pings' own settings. This facilitates external
calibration.

"""
function Sv(pings::Vector{EK60Ping};
            frequency=nothing,
            gain=nothing,
            equivalentbeamangle=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            pulselength=nothing,
            sacorrection=nothing,
            rangecorrectionoffset=2)
    s = [Sv(ping,
            frequency=frequency,
            gain=gain,
            equivalentbeamangle=equivalentbeamangle,
            soundvelocity=soundvelocity,
            absorptioncoefficient=absorptioncoefficient,
            transmitpower=transmitpower,
            pulselength=pulselength,
            sacorrection=sacorrection,
            rangecorrectionoffset=rangecorrectionoffset) for ping in pings]

    myhcat(s)
end

"""
    TS(pings::Vector{EK60Ping};
            frequency=nothing,
            gain=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            rangecorrectionoffset=0)

Target strength
"""
function TS(pings::Vector{EK60Ping};
            frequency=nothing,
            gain=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            rangecorrectionoffset=0)
    s = [TS(ping,
            frequency=frequency,
            gain=gain,
            soundvelocity=soundvelocity,
            absorptioncoefficient=absorptioncoefficient,
            transmitpower=transmitpower,
            rangecorrectionoffset=rangecorrectionoffset) for ping in pings]

    myhcat(s)
end

function Sv(pings::Channel{EK60Ping};
            frequency=nothing,
            gain=nothing,
            equivalentbeamangle=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            pulselength=nothing,
            sacorrection=nothing,
            rangecorrectionoffset=2)
    Sv(collect(pings),
       frequency=frequency,
       gain=gain,
       equivalentbeamangle=equivalentbeamangle,
       soundvelocity=soundvelocity,
       absorptioncoefficient=absorptioncoefficient,
       transmitpower=transmitpower,
       pulselength=pulselength,
       sacorrection=sacorrection,
       rangecorrectionoffset=rangecorrectionoffset)
end

function TS(pings::Channel{EK60Ping};
            frequency=nothing,
            gain=nothing,
            soundvelocity=nothing,
            absorptioncoefficient=nothing,
            transmitpower=nothing,
            rangecorrectionoffset=0)
    TS(collect(pings),
       frequency=frequency,
       gain=gain,
       soundvelocity=soundvelocity,
       absorptioncoefficient=absorptioncoefficient,
       transmitpower=transmitpower,
       rangecorrectionoffset=rangecorrectionoffset)
end

#

"""
    athwartshipangle(ping::EK60Ping)

Returns the athwartship phase angle difference vector for the given
ping.

"""
athwartshipangle(ping::EK60Ping) = ping.athwartshipangle

function athwartshipangle(pings::Vector{EK60Ping})
    s = [ping.athwartshipangle for ping in pings]
    myhcat(s, missingvalue=Int8(0))
end

athwartshipangle(pings::Channel{EK60Ping}) = athwartshipangle(collect(pings))

@deprecate(athwartshipanglematrix, athwartshipangle)

#

"""
    alongshipangle(ping::EK60Ping)

Returns the alongship phase angle difference vector for the given
ping.

"""
alongshipangle(ping::EK60Ping) = ping.alongshipangle

function alongshipangle(pings::Vector{EK60Ping})
    s = [ping.alongshipangle for ping in pings]
    myhcat(s, missingvalue=Int8(0))
end

alongshipangle(pings::Channel{EK60Ping}) = alongshipangle(collect(pings))

@deprecate(alongshipanglematrix, alongshipangle)

#

"""
    power(ping::EK60Ping)

Returns power data for a ping in manufacturer units. For decibels,
see `powerdb`.

"""
power(ping::EK60Ping) = ping.power

function power(pings::Vector{EK60Ping})
    s = [ping.power for ping in pings]
    myhcat(s, missingvalue=Int16(0))
end

power(pings::Channel{EK60Ping}) = power(collect(pings))

#

const POWER_MULTIPLIER = FLOAT_TYPE(10 * log10(2) / 256)

"""
    power2db(x)

EK60 power data is stored in a "compressed", manufacturer specific
format. This function converts to decibels.

"""
power2db(x) = x * POWER_MULTIPLIER

"""
    powerdb(x)

Returns power data in decibels. x is one or more `Ping`s.

"""
powerdb(x) = power2db(power(x))

"""
    octet2deg(x)

EK60 split beam angles are stored as signed octets with a resolution
of 180/128.  Convert from octet to degrees.
"""
octet2deg(x) = x * 180/128

end # module

ext: jl


text:
module Bingomatic
    using Plots, StatsBase, TextWrap
    export make_card, sample_words, savefig

    """
    * `word_pool`: a list of words
    * `n`: the dimensions of the sampled words
    * `free_space`: sets center cell to "Free Space" if true
    """
    function sample_words(word_pool, n...=(5,5)...; free_space=true)
        n_cells = prod(n)
        center = Int(ceil(n_cells/2))
        length(word_pool) < n_cells ? error("list must have more elements than cells") : nothing
        iseven(n_cells) ? error("Must have an odd number of cells") : nothing
        cells = sample(word_pool, n, replace=false)
        cells[center] = free_space ? "Free Space" : cells[center]
        return cells
    end

    """
    * `labels`: an array containing bingo words
    * `n_rows`: the number of rows and columns in each bingo card
    * `word_size`: the size of words in points
    * `line_width`: the number of characters per line for each word or phrase
    * `break_words`: break long words when wrapping words
    * `kwargs`: variable keyword arguments to pass to the `plot` function
    """
    function make_card(list; word_size=12, line_width=8, break_words=false, kwargs...)
        min_x = 0
        max_x = 2
        n_rows = size(list, 1)
        n_cells = n_rows^2
        offset = 1/n_rows
        iseven(n_rows) ? error("n_rows must be odd") : nothing
        length(list) < n_cells ? error("labels must have more elements than cells") : nothing
        ticks = range(min_x, max_x, length=n_rows+1)
        scale = range(min_x+offset, max_x-offset, length=n_rows)
        x = [[i j] for i in scale for j in scale] |> x -> vcat(x...)
        w_list = map(x -> TextWrap.wrap(x, width=line_width, break_long_words=break_words), list)
        return scatter(x[:,1], x[:,2], markersize=0, xlims=(min_x,max_x), ylims=(min_x,max_x),
            series_annotations = text.(w_list[:], :center, word_size), label = false, grid=true,
            gridstyle=:solid, gridlinewidth=1.5, gridalpha=1, xticks=ticks, yticks=ticks,
            xaxis=font(0), yaxis=font(0), xlabel="", ylabel="", frame_style=:box, size=(400,400);
            kwargs...)
    end
end


text:
@defcomp ch4forcing begin
    c_N2Oconcentration=Parameter(index=[time],unit="ppbv")
    c_CH4concentration=Parameter(index=[time],unit="ppbv")
    f0_CH4baseforcing=Parameter(unit="W/m2", default=0.550)
    fslope_CH4forcingslope=Parameter(unit="W/m2", default=0.036)
    c0_baseN2Oconc=Parameter(unit="ppbv", default=322.)
    c0_baseCH4conc=Parameter(unit="ppbv", default=1860.)
    f_CH4forcing=Variable(index=[time],unit="W/m2")
    over_baseoverlap=Variable(unit="W/m2")
    over=Variable(index=[time],unit="W/m2")

    function run_timestep(p, v, d, t)

        #from p.16 in Hope 2009
        if is_first(t)
            #calculate baseline forcing overlap in first time period
            v.over_baseoverlap=-0.47*log(1+2.0e-5*(p.c0_baseN2Oconc*p.c0_baseCH4conc)^0.75+5.3e-15*p.c0_baseCH4conc*(p.c0_baseCH4conc*p.c0_baseN2Oconc)^1.52)
        end

        v.over[t]=-0.47*log(1+2.0e-5*(p.c_CH4concentration[t]*p.c0_baseN2Oconc)^0.75+5.3e-15*p.c_CH4concentration[t]*(p.c0_baseN2Oconc*p.c_CH4concentration[t])^1.52)
        v.f_CH4forcing[t]=p.f0_CH4baseforcing+p.fslope_CH4forcingslope*(sqrt(p.c_CH4concentration[t])-sqrt(p.c0_baseCH4conc))+v.over[t]-v.over_baseoverlap
    end
end


text:
using Base.Test

# check to see what's on our path and what directory we are in
println(ENV["PATH"])
println("Current working directory ", pwd())
println("Files in this directory ")
run(`ls`)

##################################
### Visibility interpolation
##################################

include("interpolation.jl")

###################################
### Basic unit tests
###################################

# Do more simple unit tests here
include("ensemble_sampler.jl")

###################################
### Functional, integration tests
###################################

######################################
### Standard model, fixed distance

# Just make sure the directory is clean
run(`bash clean.sh`)

# Make a fake model so we can test off of it.
println("Initializing directory with a new standard project.")
run(`DJ_initialize.jl --new-project=standard`)

# Here test the fix parameters reading routines
include("parameters_fixed.jl")

# Set up the config file and pos0.npy
println("Editing config.yaml to use the standard model with fixed distance.")
run(`python edit_config.py --fix-distance --model=standard`)

# Run the initialization and plotting scripts
run(`make all`)

# The tricky thing is that this sampling will go on for a long time. How can we get the essential test of venus.jl without needing to wait for a full iteration across 24 walkers?  I mean, is it possible to run with fewer walkers than parameters?

# Run venus.jl and see if we sample properly.
run(`venus.jl --test`)

# Clean up before moving on to the next test
run(`make clean`)

######################################
### Vertical model, fixed distance

# Just make sure the directory is clean
run(`bash clean.sh`)

# Make a fake model so we can test off of it.
println("Initializing directory with a new vertical project.")
run(`DJ_initialize.jl --new-project=vertical`)

# Set up the config file and pos0.npy
println("Editing config.yaml to use the vertical model with fixed distance.")
run(`python edit_config.py --fix-distance --model=vertical`)

# Run the initialization and plotting scripts
run(`make all`)

# Run venus.jl and see if we sample properly. --test flag limits the sampling to a few walkers.
run(`venus.jl --test`)

# Clean up before moving on to the next test
run(`make clean`)

######################################
### Cavity model, fixed distance

# Just make sure the directory is clean
run(`bash clean.sh`)

# Make a fake model so we can test off of it.
println("Initializing directory with a new vertical project.")
run(`DJ_initialize.jl --new-project=cavity`)

# Set up the config file and pos0.npy
println("Editing config.yaml to use the cavity model with fixed distance.")
run(`python edit_config.py --fix-distance --model=cavity`)

# Run the initialization and plotting scripts
run(`make all`)

# The tricky thing is that this sampling will go on for a long time. How can we get the essential test of venus.jl without needing to wait for a full iteration across 24 walkers?  I mean, is it possible to run with fewer walkers than parameters?

# Run venus.jl and see if we sample properly.
run(`venus.jl --test`)

# Clean up before moving on to the next test
run(`make clean`)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Serializer: serialize_cycle, deserialize_cycle, writetag,
                      __deserialized_types__, serialize_typename, deserialize_typename,
                      TYPENAME_TAG, reset_state, serialize_type

import Base.Serializer: object_number, lookup_object_number, remember_object

mutable struct ClusterSerializer{I<:IO} <: AbstractSerializer
    io::I
    counter::Int
    table::ObjectIdDict
    pending_refs::Vector{Int}

    pid::Int                                     # Worker we are connected to.
    tn_obj_sent::Set{UInt64}                     # TypeName objects sent
    glbs_sent::Dict{UInt64, UInt64}              # (key,value) -> (object_id, hash_value)
    glbs_in_tnobj::Dict{UInt64, Vector{Symbol}}  # Track globals referenced in
                                                 # anonymous functions.
    anonfunc_id::UInt64

    function ClusterSerializer{I}(io::I) where I<:IO
        new(io, 0, ObjectIdDict(), Int[], worker_id_from_socket(io),
            Set{UInt64}(), Dict{UInt64, UInt64}(), Dict{UInt64, Vector{Symbol}}(), 0)
    end
end
ClusterSerializer(io::IO) = ClusterSerializer{typeof(io)}(io)

const object_numbers = WeakKeyDict()
const obj_number_salt = Ref(0)
function object_number(s::ClusterSerializer, @nospecialize(l))
    global obj_number_salt, object_numbers
    if haskey(object_numbers, l)
        return object_numbers[l]
    end
    # a hash function that always gives the same number to the same
    # object on the same machine, and is unique over all machines.
    ln = obj_number_salt[]+(UInt64(myid())<<44)
    obj_number_salt[] += 1
    object_numbers[l] = ln
    return ln::UInt64
end

const known_object_data = Dict{UInt64,Any}()

function lookup_object_number(s::ClusterSerializer, n::UInt64)
    return get(known_object_data, n, nothing)
end

function remember_object(s::ClusterSerializer, @nospecialize(o), n::UInt64)
    known_object_data[n] = o
    if isa(o, TypeName) && !haskey(object_numbers, o)
        # set up reverse mapping for serialize
        object_numbers[o] = n
    end
    return nothing
end

function deserialize(s::ClusterSerializer, ::Type{TypeName})
    full_body_sent = deserialize(s)
    number = read(s.io, UInt64)
    if !full_body_sent
        tn = lookup_object_number(s, number)::TypeName
        remember_object(s, tn, number)
        deserialize_cycle(s, tn)
    else
        tn = deserialize_typename(s, number)
    end

    # retrieve arrays of global syms sent if any and deserialize them all.
    foreach(sym->deserialize_global_from_main(s, sym), deserialize(s))
    return tn
end

function serialize(s::ClusterSerializer, t::TypeName)
    serialize_cycle(s, t) && return
    writetag(s.io, TYPENAME_TAG)

    identifier = object_number(s, t)
    send_whole = !(identifier in s.tn_obj_sent)
    serialize(s, send_whole)
    write(s.io, identifier)
    if send_whole
        # Track globals referenced in this anonymous function.
        # This information is used to resend modified globals when we
        # only send the identifier.
        prev = s.anonfunc_id
        s.anonfunc_id = identifier
        serialize_typename(s, t)
        s.anonfunc_id = prev
        push!(s.tn_obj_sent, identifier)
        finalizer(t) do x
            cleanup_tname_glbs(s, identifier)
        end
    end

    # Send global refs if required.
    syms = syms_2b_sent(s, identifier)
    serialize(s, syms)
    foreach(sym->serialize_global_from_main(s, sym), syms)
    nothing
end

function serialize(s::ClusterSerializer, g::GlobalRef)
    # Record if required and then invoke the default GlobalRef serializer.
    sym = g.name
    if g.mod === Main && isdefined(g.mod, sym)
        if (binding_module(Main, sym) === Main) && (s.anonfunc_id != 0) &&
            !startswith(string(sym), "#") # Anonymous functions are handled via FULL_GLOBALREF_TAG

            push!(get!(s.glbs_in_tnobj, s.anonfunc_id, []), sym)
        end
    end

    invoke(serialize, Tuple{AbstractSerializer, GlobalRef}, s, g)
end

# Send/resend a global object if
# a) has not been sent previously, i.e., we are seeing this object_id for the first time, or,
# b) hash value has changed or
# c) is a bits type
function syms_2b_sent(s::ClusterSerializer, identifier)
    lst = Symbol[]
    check_syms = get(s.glbs_in_tnobj, identifier, [])
    for sym in check_syms
        v = getfield(Main, sym)

        if isbits(v)
            push!(lst, sym)
        else
            oid = object_id(v)
            if haskey(s.glbs_sent, oid)
                # We have sent this object before, see if it has changed.
                s.glbs_sent[oid] != hash(sym, hash(v)) && push!(lst, sym)
            else
                push!(lst, sym)
            end
        end
    end
    return unique(lst)
end

function serialize_global_from_main(s::ClusterSerializer, sym)
    v = getfield(Main, sym)

    oid = object_id(v)
    record_v = true
    if isbits(v)
        record_v = false
    elseif !haskey(s.glbs_sent, oid)
        # set up a finalizer the first time this object is sent
        try
            finalizer(v) do x
                delete_global_tracker(s,x)
            end
        catch ex
            # Do not track objects that cannot be finalized.
            if isa(ex, ErrorException)
                record_v = false
            else
                rethrow(ex)
            end
        end
    end
    record_v && (s.glbs_sent[oid] = hash(sym, hash(v)))

    serialize(s, isconst(Main, sym))
    serialize(s, v)
end

function deserialize_global_from_main(s::ClusterSerializer, sym)
    sym_isconst = deserialize(s)
    v = deserialize(s)
    if sym_isconst
        @eval Main const $sym = $v
    else
        @eval Main $sym = $v
    end
end

function delete_global_tracker(s::ClusterSerializer, v)
    oid = object_id(v)
    if haskey(s.glbs_sent, oid)
        delete!(s.glbs_sent, oid)
    end

    # TODO: A global binding is released and gc'ed here but it continues
    # to occupy memory on the remote node. Would be nice to release memory
    # if possible.
end

function cleanup_tname_glbs(s::ClusterSerializer, identifier)
    delete!(s.glbs_in_tnobj, identifier)
end

# TODO: cleanup from s.tn_obj_sent

# Specialized serialize-deserialize implementations for CapturedException to partially
# recover from any deserialization errors in `CapturedException.ex`

function serialize(s::ClusterSerializer, ex::CapturedException)
    serialize_type(s, typeof(ex))
    serialize(s, string(typeof(ex.ex))) # String type should not result in a deser error
    serialize(s, ex.processed_bt)       # Currently should not result in a deser error
    serialize(s, ex.ex)                 # can result in a UndefVarError on the remote node
                                        # if a type used in ex.ex is undefined on the remote node.
end

function original_ex(s::ClusterSerializer, ex_str, remote_stktrace)
    local pid_str = ""
    try
        pid_str = string(" from worker ", worker_id_from_socket(s.io))
    end

    stk_str = remote_stktrace ? "Remote" : "Local"
    ErrorException(string("Error deserializing a remote exception", pid_str, "\n",
                          "Remote(original) exception of type ", ex_str, "\n",
                          stk_str,  " stacktrace : "))
end

function deserialize(s::ClusterSerializer, t::Type{<:CapturedException})
    ex_str = deserialize(s)
    local bt
    local capex
    try
        bt = deserialize(s)
    catch e
        throw(CompositeException([
            original_ex(s, ex_str, false),
            CapturedException(e, catch_backtrace())
        ]))
    end

    try
        capex = deserialize(s)
    catch e
        throw(CompositeException([
            CapturedException(original_ex(s, ex_str, true), bt),
            CapturedException(e, catch_backtrace())
        ]))
    end

    return CapturedException(capex, bt)
end

"""
    clear!(syms, pids=workers(); mod=Main)

Clears global bindings in modules by initializing them to `nothing`.
`syms` should be of type `Symbol` or a collection of `Symbol`s . `pids` and `mod`
identify the processes and the module in which global variables are to be
reinitialized. Only those names found to be defined under `mod` are cleared.

An exception is raised if a global constant is requested to be cleared.
"""
function clear!(syms, pids=workers(); mod=Main)
    @sync for p in pids
        @async remotecall_wait(clear_impl!, p, syms, mod)
    end
end
clear!(sym::Symbol, pid::Int; mod=Main) = clear!([sym], [pid]; mod=mod)
clear!(sym::Symbol, pids=workers(); mod=Main) = clear!([sym], pids; mod=mod)
clear!(syms, pid::Int; mod=Main) = clear!(syms, [pid]; mod=mod)

clear_impl!(syms, mod::Module) = foreach(x->clear_impl!(x,mod), syms)
clear_impl!(sym::Symbol, mod::Module) = isdefined(mod, sym) && @eval(mod, global $sym = nothing)

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using SQLREPL
using Test

@testset "SQLREPL.jl" begin
    # Write your tests here.
end


text:
module Cairn
include("rpc.jl")
include("client.jl")
include("server.jl")
end


text:
@testset "tree creation" begin
	eM = eMesh_sphere()

	@test_throws ErrorException eMesh_to_tree(eM)
	num_tri = length(eM.tri)
	num_tet = length(eM.tet)

	tree = eMesh_to_tree(as_tri_eMesh(eM))
	@test leafNumber(tree) == num_tri
	@test Set(extractData(tree)) == Set(1:num_tri)
	@test treeDepth(tree) < (log2(num_tri) * 1.3)

	tree = eMesh_to_tree(as_tet_eMesh(eM))
	@test leafNumber(tree) == num_tet
	@test Set(extractData(tree)) == Set(1:num_tet)
	@test treeDepth(tree) < (log2(num_tet) * 1.3)
end


text:
if !isdir("downloads")
    mkdir("downloads")
end
cd("downloads")

if !isdir("BSR")
    println("Downloading BSR code...")
    run(`wget http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/BSR/BSR_code.tgz`)
    println("Extracting...")
    run(`tar xvzf BSR_code.tgz`)
    rm("BSR_code.tgz")
end


text:
using Documenter, KernelDensityDistributionEsimates

makedocs(;
    modules=[KernelDensityDistributionEsimates],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chriselrod/KernelDensityDistributionEsimates.jl/blob/{commit}{path}#L{line}",
    sitename="KernelDensityDistributionEsimates.jl",
    authors="Chris Elrod",
    assets=[],
)

deploydocs(;
    repo="github.com/chriselrod/KernelDensityDistributionEsimates.jl",
)


text:
struct DrWatson32{f} <: AbstractWatson
    val::Float32
end

# conversions back from DrWatsons
Base.Float64(x::DrWatson32) = Float64(x.val)
Base.Float32(x::DrWatson32) = x.val
Base.Float16(x::DrWatson32) = Float16(x.val)
Base.Int64(x::DrWatson32) = Int64(x.val)
Base.Int32(x::DrWatson32) = Int32(x.val)
Base.Int16(x::DrWatson32) = Int16(x.val)

# conversion among DrWatsons
DrWatson32{f}(x::DrWatson32) where f = DrWatson32{f}(x.val)
DrWatson32(x::DrWatson32{f}) where {f} = DrWatson32{f}(x.val)

Base.promote_rule(::Type{Int64},::Type{DrWatson32{f}}) where f = DrWatson32{f}
Base.promote_rule(::Type{Int32},::Type{DrWatson32{f}}) where f = DrWatson32{f}
Base.promote_rule(::Type{Int16},::Type{DrWatson32{f}}) where f = DrWatson32{f}

Base.promote_rule(::Type{Float64},::Type{DrWatson32{f}}) where f = DrWatson32{f}
Base.promote_rule(::Type{Float32},::Type{DrWatson32{f}}) where f = DrWatson32{f}
Base.promote_rule(::Type{Float16},::Type{DrWatson32{f}}) where f = DrWatson32{f}

Base.bitstring(x::DrWatson32) = bitstring(x.val)

Base.eps(::Type{DrWatson32{f}}) where f = eps(Float32)
Base.eps(x::DrWatson32) = eps(x.val)

Base.typemin(::DrWatson32{f}) where f = DrWatson32{f}(typemin(Float32))
Base.typemax(::DrWatson32{f}) where f = DrWatson32{f}(typemax(Float32))
Base.floatmin(::DrWatson32{f}) where f = DrWatson32{f}(floatmin(Float32))
Base.floatmax(::DrWatson32{f}) where f = DrWatson32{f}(floatmax(Float32))
Base.precision(::DrWatson32{f}) where f = DrWatson32{f}(precision(Float32))

function +(x::DrWatson32{f},y::DrWatson32{f}) where f
    r = x.val + y.val
    scent(f,r)
    return DrWatson32{f}(r)
end

function -(x::DrWatson32{f},y::DrWatson32{f}) where f
    r = x.val - y.val
    scent(f,r)
    return DrWatson32{f}(r)
end

function *(x::DrWatson32{f},y::DrWatson32{f}) where f
    r = x.val * y.val
    scent(f,r)
    return DrWatson32{f}(r)
end

function /(x::DrWatson32{f},y::DrWatson32{f}) where f
    r = x.val / y.val
    scent(f,r)
    return DrWatson32{f}(r)
end

function ^(x::DrWatson32{f},y::DrWatson32{f}) where f
    r = x.val ^ y.val
    scent(f,r)
    return DrWatson32{f}(r)
end

for O in ( :(-), :(+),
           :sign,
           :prevfloat, :nextfloat,
           :round, :trunc, :ceil, :floor,
           :inv, :abs, :sqrt, :cbrt,
           :exp, :expm1, :exp2, :exp10,
           :log, :log1p, :log2, :log10,
           :rad2deg, :deg2rad, :mod2pi, :rem2pi,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
           :sinc, :sinpi, :cospi,
           :sind, :cosd, :tand, :cscd, :secd, :cotd,
           :asind, :acosd, :atand, :acscd, :asecd, :acotd
          )
    @eval begin
        function Base.$O(x::DrWatson32{f}) where f
            r = $O(x.val)
            scent(f,r)
            return DrWatson32{f}(r)
        end
    end
end

for O in ( :(<), :(<=))     # do not trigger scent for those
    @eval begin
        Base.$O(x::DrWatson32, y::DrWatson32) = $O(x.val, y.val)
    end
end


text:
macro load_soss_fmi()
    expr = quote
        @info "Loading foreign model interface to \u001b[3m\u001b[34;1mSoss.jl\u001b[0m\n\n          \u001b[34;1mhttps://github.com/cscherrer/Soss.jl\n\nThis interface is a work in progress! Please post issues to:\n\n          \u001b[34;1mhttps://github.com/femtomc/Jaynes.jl/issues\n\n"

        using PyCall
        using Conda

        # Install SymPy.
        if PyCall.conda
            Conda.add("pip")
            pip = joinpath(Conda.BINDIR, "pip")
            run(`$pip install sympy`)
        else
            try
                pyimport("sympy")
            catch ee
                typeof(ee) <: PyCall.PyError || rethrow(ee)
                warn("""
                     Python Dependencies not installed
                     Please either:
                     - Rebuild PyCall to use Conda, by running in the julia REPL:
                     - `ENV[PYTHON]=""; Pkg.build("PyCall"); Pkg.build("Jaynes")`
                     - Or install the dependencies: `pip install sympy`
                     """)
            end
        end

        import Jaynes: has_top, get_top, has_sub, get_sub, get_score, collect!, add_call!, haskey
        import Jaynes: simulate, propose, generate, update, regenerate, score
        using Soss
        import Soss: logpdf

        # ------------ Trace ------------ #

        mutable struct SossTrace{NT <: NamedTuple} <: Jaynes.Trace
            choices::NT
            SossTrace(nt::NT) where NT = new{NT}(nt)
            SossTrace{NT}() where NT = new{NT}()
        end
        add_call!(st::SossTrace, addr, choices) = st.choices = choices

        # ------------ Call site ------------ #

        mutable struct SossModelCallSite{M, A} <: Jaynes.CallSite
            trace::SossTrace
            score::Float64
            model::M
            args::A
        end

        get_score(gfcs::SossModelCallSite) = gfcs.score
        haskey(cs::SossModelCallSite, addr) = haskey(trace, addr)
        getindex(cs::SossModelCallSite, addrs...) = getindex(trace, addr)
        get_ret(cs::SossModelCallSite) = cs.trace

        # ------------ Pretty printing ------------ #

        function collect!(par::P, addrs::Vector{Any}, chd::Dict{Any, Any}, tr::ST, meta) where {P <: Tuple, ST <: SossTrace}
            for (k, v) in pairs(tr.choices)
                push!(addrs, (par..., k))
                chd[(par..., k)] = v
                meta[(par..., k)] = "(Soss)"
            end
        end

        function collect!(addrs::Vector{Any}, chd::Dict{Any, Any}, tr::ST, meta) where ST <: SossTrace
            for (k, v) in pairs(tr.choices)
                push!(addrs, (k, ))
                chd[(k, )] = v
            end
        end

        # ------------ Contexts ------------ #

        function (ctx::Jaynes.SimulateContext)(c::typeof(foreign),
                                               addr::T,
                                               model::M,
                                               args...) where {T <: Jaynes.Address, M <: Soss.Model}
            choices = rand(model(args...))
            score = Soss.logpdf(m(args...), choices)
            Jaynes.add_call!(ctx, addr, SossModelCallSite(SossTrace(choices), score, model, args))
            return choices
        end

        # Convenience.
        function simulate(model::M, args...) where M <: Soss.Model
            ctx = Jaynes.Simulate()
            addr = gensym()
            ret = ctx(foreign, addr, model, args...)
            return ret, get_sub(ctx.tr, addr)
        end

        function (ctx::Jaynes.ProposeContext)(c::typeof(foreign),
                                              addr::T,
                                              model::M,
                                              args...) where {T <: Jaynes.Address, M <: Soss.Model}
            choices = rand(model(args...))
            score = Soss.logpdf(m(args...), choices)
            Jaynes.add_call!(ctx, addr, SossModelCallSite(SossTrace(choices), score, model, args))
            increment!(ctx, score)
            return choices
        end

        # Convenience.
        function propose(model::M, args...) where M <: Soss.Model
            ctx = Propose()
            addr = gensym()
            ret = ctx(foreign, addr, model, args...)
            return ret, get_top(ctx.tr, addr), ctx.score
        end

        function (ctx::Jaynes.GenerateContext)(c::typeof(foreign),
                                               addr::T,
                                               model::M,
                                               args...) where {T <: Jaynes.Address, M <: Soss.Model}
            data = Jaynes.get_top(ctx.select, addr)
            w, choices = Soss.weightedSample(model(args...), data)
            score = Soss.logpdf(m(args...), choices)
            Jaynes.add_call!(ctx, addr, SossModelCallSite(SossTrace(choices), score, model, args))
            Jaynes.increment!(ctx, w)
            return choices
        end

        # Convenience.
        function generate(sel::L, model::M, args...) where {L <: NamedTuple, M <: Soss.Model}
            addr = gensym()
            v_sel = selection([(addr, ) => sel])
            ctx = Generate(v_sel)
            ret = ctx(foreign, addr, model, args...)
            return ret, get_sub(ctx.tr, addr), ctx.weight
        end

        function (ctx::Jaynes.UpdateContext)(c::typeof(foreign),
                                             addr::T,
                                             model::M,
                                             args...) where {T <: Jaynes.Address, M <: Soss.Model}
            prev = Jaynes.get_prev(ctx, addr)
            kvs = Jaynes.get_top(ctx.select, addr)
            data = namedtuple(Dict{Symbol, Any}(kvs))
            w, choices = Soss.weightedSample(model(args...), data)
            score = Soss.logpdf(m(args...), choices)
            Jaynes.add_call!(ctx, addr, SossModelCallSite(SossTrace(choices), score, model, args))
            Jaynes.increment!(ctx, w - get_score(prev))
            return choices
        end

        # Convenience.
        function update(sel::L, soss_cl::C) where {L <: NamedTuple, C <: SossModelCallSite}
            addr = gensym()
            v_sel = selection([(addr, ) => sel])
            ctx = Jaynes.UpdateContext(cl, v_sel, NoChange())
            ret = ctx(foreign, addr, model, args...)
            return ret, get_sub(ctx.tr, addr), ctx.weight, Jaynes.UndefinedChange(), nothing
        end

        function (ctx::Jaynes.RegenerateContext{C, T})(c::typeof(foreign),
                                                       addr::A,
                                                       model::M,
                                                       args...) where {A <: Jaynes.Address, M <: Soss.Model, C <: SossModelCallSite, T <: SossTrace}
            targeted = dump_queries(Jaynes.get_sub(ctx.select, addr))
            prev = ctx.prev
            kvs = Dict{Symbol, Any}()
            for (k, v) in pairs(prev.trace.choices)
                haskey(ctx.select, (addr, k)) && continue
                kvs[k] = v
            end
            data = Soss.namedtuple(kvs)
            w, choices = Soss.weightedSample(model(args...), data)
            score = Soss.logpdf(m(args...), choices)
            ctx.tr.choices = choices
            ctx.score += score
            Jaynes.increment!(ctx, score - get_score(prev))
            return choices
        end

        function (ctx::Jaynes.RegenerateContext)(c::typeof(foreign),
                                                 addr::T,
                                                 model::M,
                                                 args...) where {T <: Jaynes.Address, M <: Soss.Model}
            targeted = dump_queries(Jaynes.get_sub(ctx.select, addr))
            prev = Jaynes.get_prev(ctx, addr)
            kvs = Dict{Symbol, Any}()
            for (k, v) in pairs(prev.trace.choices)
                k in targeted && continue
                kvs[k] = v
            end
            data = Soss.namedtuple(kvs)
            w, choices = Soss.weightedSample(model(args...), data)
            score = Soss.logpdf(m(args...), choices)
            Jaynes.add_call!(ctx, addr, SossModelCallSite(SossTrace(choices), score, model, args))
            Jaynes.increment!(ctx, score - get_score(prev))
            return choices
        end

        # Convenience.
        function regenerate(sel::L, soss_cl::C) where {L <: Jaynes.UnconstrainedSelection, C <: SossModelCallSite}
            addr = gensym()
            v_sel = selection(addr => sel)
            ctx = Jaynes.Regenerate(soss_cl, v_sel, NoChange())
            ret = ctx(foreign, addr, soss_cl.model, soss_cl.args...)
            return ret, SossModelCallSite(ctx.tr, ctx.score, soss_cl.model, soss_cl.args), ctx.weight, Jaynes.UndefinedChange(), nothing
        end

        function (ctx::Jaynes.ScoreContext)(c::typeof(foreign),
                                            addr::T,
                                            model::M,
                                            args...) where {T <: Jaynes.Address, M <: Soss.Model}
            kvs = Jaynes.get_sub(ctx.select, addr)
            choices = namedtuple(Dict{Symbol, Any}(kvs))
            score = Soss.logpdf(m(args...), choices)
            Jaynes.increment!(ctx, score)
            return choices
        end

        # Convenience.
        function score(sel::L, model::M, args...) where {L <: Jaynes.UnconstrainedSelection, M <: SossModelCallSite}
            addr = gensym()
            v_sel = selection(addr => sel)
            ctx = Jaynes.Score(v_sel)
            ret = ctx(foreign, addr, model, args...)
            return ret, ctx.weight
        end
    end

    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end


text:
using Flux
using Flux: onehotbatch
using MLDatasets
using MLDataPattern: splitobs

get_test_data(args) = @eval $(Symbol(args.name, :_testdata))(;b_size=$args.batchsize, splitr_=$args.splitr_)
get_train_data(args) = @eval $(Symbol(args.name, :_traindata))(;b_size=$args.batchsize, splitr_=$args.splitr_)

function CIFAR10_traindata(;b_size=32, splitr_=.1)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    x, y = CIFAR10.traindata(Float32)

    (train_x, train_y), (val_x, val_y) = splitobs((x, y), at=1-splitr_)

    train_y = onehotbatch(train_y, 0:9)
    val_y = onehotbatch(val_y, 0:9)

    train = Flux.Data.DataLoader((train_x, train_y), batchsize=b_size, shuffle=true)
    val = Flux.Data.DataLoader((val_x, val_y), batchsize=b_size)

    return train, val
end

function CIFAR10_testdata(;b_size=32, splitr_=.1)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
    test_x, test_y = CIFAR10.testdata(Float32)
    test_y = onehotbatch(test_y, 0:9)

    return Flux.Data.DataLoader((test_x, test_y), batchsize=b_size)
end

function MNIST_traindata(;b_size=32, splitr_=.1)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    x, y = MNIST.traindata(Float32)
    x = reshape(x, 28, 28, 1, :)
    (train_x, train_y), (val_x, val_y) = splitobs((x, y), at=1-splitr_)

    train_y = onehotbatch(train_y, 0:9)
    val_y = onehotbatch(val_y, 0:9)

    train = Flux.Data.DataLoader((train_x, train_y), batchsize=b_size, shuffle=true)
    val = Flux.Data.DataLoader((val_x, val_y), batchsize=b_size)

    return train, val
end

function MNIST_testdata(;b_size=32, splitr_=.1)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
    test_x, test_y = MNIST.testdata(Float32)
    test_x = reshape(test_x, 28, 28, 1, :)
    test_y = onehotbatch(test_y, 0:9)

    return Flux.Data.DataLoader((test_x, test_y), batchsize=b_size)
end

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
using ColorTypes
using ColorTypes.FixedPointNumbers
using Test

@test isempty(detect_ambiguities(ColorTypes, Base, Core))

using Documenter
doctest(ColorTypes, manual = false)

# if the test below fails, please extend the list of types at the call to
# make_alpha in types.jl (this is the price of making that list explicit)
ctypes = union(setdiff(ColorTypes.parametric3, (XRGB, RGBX)), (Gray,))
@test Set(ctypes) ==
    Set([DIN99d, DIN99o, DIN99, HSI, HSL, HSV, LCHab, LCHuv,
         LMS, Lab, Luv, XYZ, YCbCr, YIQ, xyY, BGR, RGB, Gray])

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :register_error_hint)
    @testset "error_hints" begin
        # ColorVectorSpace, if needed, should not be imported before this
        include("error_hints.jl")
    end
end

@testset "conversions" begin
    include("conversions.jl")
end

@testset "operations" begin
    include("operations.jl")
end

@testset "show" begin
    include("show.jl")
end

@testset "traits" begin
    include("traits.jl")
end

@testset "types" begin
    include("types.jl")
end

 @testset "misc." begin
     # This tests the same thing as the last, but in a user-observable way
     let a = Array{Gray}(undef, 1)
         a[1] = Gray(0)
         a[1] = 1
         @test a[1] === Gray(1)
     end

     # broadcast
     @test eltype(broadcast(RGB, [BGR(1,0,0)])) === RGB{N0f8}
     addred(x1::AbstractRGB, x2::AbstractRGB) = red(x1) + red(x2)
     @test addred.([RGB(1,0,0)], RGB(1.0,0,0)) == [2.0]
 end

@testset "type-centric" begin
    @testset "$C" for C in ColorTypes.parametric3
        @test eltype(C{Float32}) === Float32
        et = (C <: AbstractRGB) ? N0f8 : Float32
        @test color_type(C(1,0,0)) === C{et}
        @test color_type(C) === C
        @test color_type(C{Float32}) === C{Float32}
        @test eltype(C{Float32}(1,0,0)) === Float32
    end

    for C in setdiff(ColorTypes.parametric3, [XRGB,RGBX])
        @testset "$A" for A in (alphacolor(C), coloralpha(C))
            @test eltype(A{Float32}) === Float32
            @test color_type(A) === C
            @test color_type(A{Float32}) === C{Float32}
            @test eltype(A(1,0.8,0.6,0.4)) === Float64
            c = A{Float64}(1,0.8,0.6,0.4)
            @test color_type(c) === C{Float64}
            cc = color(c)
            @test cc === C{Float64}(1,0.8,0.6)
            @test A(cc) === A{Float64}(1,0.8,0.6,1)
            @test A(cc, 0.4)  === c
            @test A(cc, 0x01) === A{Float64}(1,0.8,0.6,1)
            @test A{Float32}(cc, 0x01) === A{Float32}(1,0.8,0.6,1)
            @test C(c         ) === C{Float64}(1,0.8,0.6)
            @test C{Float32}(c) === C{Float32}(1,0.8,0.6)
            @test convert(A, cc) === A{Float64}(1,0.8,0.6,1)
            @test A(cc) === A{Float64}(1,0.8,0.6,1)
            @test A{Float64}(cc) === A{Float64}(1,0.8,0.6,1)
            @test convert(A, cc, 0.4)  === c
            @test convert(A, cc, 0x01) === A{Float64}(1,0.8,0.6,1)
            @test convert(A{Float32}, cc, 0x01) === A{Float32}(1,0.8,0.6,1)
            @test convert(C,          c) === C{Float64}(1,0.8,0.6)
            @test convert(C{Float32}, c) === C{Float32}(1,0.8,0.6)
            @test C{Float32}(c) === C{Float32}(1,0.8,0.6)
        end
    end
end


text:
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, -1<=x[1:2]<=10)

@constraint(m, con1[i=1:2], x[i]<=[3,4][i])
@constraint(m, con2, 0<=x[1])
@constraint(m, con3, 0<=x[2])
@constraint(m, con4[i=1:2], [5,4][i]<=x[i])
@constraint(m, con5, x[1] <= 9)
@constraint(m, con6, x[2] <= 6)

M = (10,(10,10,10))
@disjunction(m,(con1,con2,con3),(con4,con5,con6),
                reformulation=:BMR,
                M=M)


text:
# wrappers

function cudnnSoftmaxForward(src::CuArray{T,4}, dest::CuArray{T,4}=src;
                             algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                             mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                             alpha=1.0, beta=0.0) where T
    cudnnSoftmaxForward(handle(), algorithm, mode,
                        cptr(alpha, src), TensorDesc(src), src,
                        cptr(beta, dest), TensorDesc(dest), dest)
    return dest
end

function cudnnSoftmaxBackward(src::CuArray{T,4}, srcDiff::CuArray{T,4}, destDiff::CuArray=srcDiff;
                              algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                              mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                              alpha=1.0, beta=0.0) where T
    cudnnSoftmaxBackward(handle(), algorithm, mode,
                         cptr(alpha, src), TensorDesc(src), src,
                         TensorDesc(srcDiff), srcDiff,
                         cptr(beta, destDiff), TensorDesc(destDiff), destDiff)
    return destDiff
end


text:
module puMPS

using LinearAlgebra
using TensorOperations
using LinearMaps
using KrylovKit
using Optim

include("MPS.jl")
using .MPS

export puMPState, rand_puMPState, mps_tensor, num_sites, bond_dim, phys_dim, set_mps_tensor!,
       apply_blockTM_l, blockTM_dense,
       expect_nn, expect,
       canonicalize_left!,
       minimize_energy_local!, vumps_opt!

export MPO_PBC_uniform, MPO_open_uniform, MPO_PBC_uniform_split, MPO_PBC_split

include("states.jl")

export puMPSTvec, tvec_tensor, momentum, mps_tensors, excitations!, excitations,
       tangent_space_metric_and_MPO, tangent_space_MPO,
       Hn_in_basis, overlap, fidelity

include("tangentspace.jl")

export ising_local_MPO, ising_PBC_MPO, ising_PBC_MPO_split, ising_OBC_MPO, ising_Hn_MPO_split,
       heis_local_MPO,
       ANNNI_local_MPO, ANNNI_OBC_MPO, ANNNI_PBC_MPO_split, ANNNI_Hn_MPO_split,
       OBF_local_MPO, OBF_OBC_MPO, OBF_PBC_MPO_split, OBF_Hn_MPO_split,
       potts3_local_MPO, potts3_OBC_MPO, potts3_PBC_MPO_split, potts3_Hn_MPO_split

include("models.jl")

end

max_stars_repo_name: FuTen/puMPS

max_issues_repo_name: FuTen/puMPS

max_forks_repo_name: FuTen/puMPS


text:
abstract type Authentication
end

struct BaseURL
    uri::HTTP.URI
end
_get_http_uri(base_url::BaseURL) = base_url.uri
function _get_http_uri_string(uri::HTTP.URI)::String
    return Base.string(uri)
end

BaseURL(base_url::AbstractString) = BaseURL(HTTP.URI(base_url))

struct Client{A <: Authentication}
    base_url::BaseURL
    auth::A
end
get_auth(client::Client) = client.auth
get_base_url(client::Client) = client.base_url
function Base.shred!(client::Client)::Nothing
    Base.shred!(client.auth)
    return nothing
end

struct Credential
    secret_buffer::Base.SecretBuffer
end
Credential() = Credential(Base.SecretBuffer())

struct AnonymousAuth <: Authentication end

struct JWTAuth <: Authentication
    jwt_cred::Credential
end
JWTAuth() = JWTAuth(Credential())

struct OAuth2 <: Authentication
    oauth2_cred::Credential
end
OAuth2() = OAuth2(Credential())

struct UsernamePassAuth <: Authentication
    username::String
    password_cred::Credential
end
@inline function UsernamePassAuth(username::AbstractString)
    return UsernamePassAuth(username, Credential())
end
@inline function UsernamePassAuth(; username::AbstractString)
    return UsernamePassAuth(username, Credential())
end


text:
mutable struct MoominModel
	data::Data
	genes::Genes
	reactions::Reactions

	function MoominModel(data::Data, genes::Genes, reactions::Reactions)
		new(data, genes, reactions)
	end
end

function nReactions(model::MoominModel)
	return nReactions(model.reactions)
end

"""
Read metabolic model and DE data from files
"""
function readModel(pathToData, pathToModel;
	delimiter='\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC",
	modelStr="model",
	pThresh=0.9, alpha=3, precision=7)

	data = readData(pathToData; delimiter=delimiter, geneIDhead=geneIDhead,
		PPDEhead=PPDEhead, logFChead=logFChead)

	modelFile = matopen(pathToModel)
	vars = matread(pathToModel)
	network = vars[modelStr]

	S = network["S"]
	nReac = size(S, 2)

	genes = Genes(getField(network, "genes", nReac))

	reactionIDs = getField(network, "rxns", nReac)
	reactionNames = getField(network, "rxnNames", nReac)
	reversible = [i<0 for i in network["lb"][:]]
	grRules = getField(network, "grRules", nReac)
	subsystems = getField(network, "subSystems", nReac)
  	assGenes = getAssGenes(getField(network, "rules", nReac))
	reactions = Reactions(reactionIDs, reactionNames, S, reversible, grRules, subsystems,
	 						assGenes)

	mapDataToGenes!(genes, data)
	colour!(genes, pThresh, alpha, precision)
	colour!(reactions, genes, pThresh, alpha, precision)

	close(modelFile)
	return MoominModel(data, genes, reactions)
end

function getField(network, fieldName, nReactions)
	if haskey(network, fieldName)
		if typeof(network[fieldName][1]) == Array{Any,2}
			return convert(Array{String,1}, [row[1] for row in network[fieldName]][:])
		else
			return convert(Array{String,1}, network[fieldName][:])
		end
	else
		@warn "Model has no entry for $fieldName"
		return fill("NA", nReactions)
	end
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
"""
The abstract base for types that manage the objective function evaluation.
`P` is the optimization problem it is used for.
"""
abstract type Evaluator{P <: OptimizationProblem} end

fitness_scheme(e::Evaluator) = fitness_scheme(problem(e))
fitness_type(e::Evaluator) = fitness_type(fitness_scheme(e))
worst_fitness(e::Evaluator) = worst_fitness(fitness_scheme(e))
numdims(e::Evaluator) = numdims(problem(e))
search_space(e::Evaluator) = search_space(problem(e))
describe(e::Evaluator) = "Problem: $(name(e.problem)) (dimensions = $(numdims(e)))"
problem_summary(e::Evaluator) = "$(name(e.problem))_$(numdims(e))d"

shutdown!(e::Evaluator) = e # do nothing

"""
Default implementation of the `Evaluator`.

`FP` is the original problem's fitness type
`FA` is the fitness type actually stored by the archive.
"""
mutable struct ProblemEvaluator{FP, FA, A<:Archive, P<:OptimizationProblem} <: Evaluator{P}
    # FIXME F is the fitness type of the problem, but with current Julia it's
    # not possible to get and use it at declaration type
    problem::P
    archive::A
    num_evals::Int
    last_fitness::FP

    ProblemEvaluator(problem::P, archive::A) where {P<:OptimizationProblem, A<:Archive} =
        new{fitness_type(fitness_scheme(problem)),fitness_type(archive),A,P}(
                problem, archive,
                0, nafitness(fitness_scheme(problem)))
end

ProblemEvaluator(problem::OptimizationProblem; archiveCapacity::Int = 10) =
    ProblemEvaluator(problem, TopListArchive(fitness_scheme(problem), numdims(problem), archiveCapacity))

problem(e::Evaluator) = e.problem
num_evals(e::ProblemEvaluator) = e.num_evals

"""
    fitness(params::Individual, e::ProblemEvaluator, tag::Int=0)

Evaluate the fitness and implicitly update the archive with the provided
parameters and calculated fitness.

Returns the fitness in the archived format.
"""
function fitness(params::Individual, e::ProblemEvaluator, tag::Int=0)
    e.last_fitness = fit = fitness(params, e.problem)
    e.num_evals += 1
    fita = archived_fitness(fit, e.archive)
    candi = add_candidate!(e.archive, fita, params, tag, e.num_evals)
    return fita
end

"""
    last_fitness(e::Evaluator)

Get the fitness of the last evaluated candidate.

Leads to nicer code if we can first test if it is better or worse than existing candidates
and only want the fitness itself if the criteria fulfilled.
"""
last_fitness(e::Evaluator) = e.last_fitness

is_better(f1::F, f2::F, e::Evaluator) where {F} = is_better(f1, f2, fitness_scheme(e))

is_better(candidate, f, e::Evaluator) = is_better(fitness(candidate, e), f, fitness_scheme(e))

function best_of(candidate1::Individual, candidate2::Individual, e::Evaluator)
    f1 = fitness(candidate1, e)
    f2 = fitness(candidate2, e)
    if is_better(f1, f2, e)
        return candidate1, f1
    else
        return candidate2, f2
    end
end

function update_fitness!(e::ProblemEvaluator{FP,FA}, candidate::Candidate{FA}) where {FP,FA}
    # evaluate fitness if not known yet
    if isnafitness(candidate.fitness, fitness_scheme(e.archive))
        candidate.fitness = fitness(candidate.params, e, candidate.tag)
    end
    return candidate
end

function update_fitness!(e::ProblemEvaluator{FP,FA}, candidates::AbstractVector{Candidate{FA}}) where {FP,FA}
    @inbounds for candidate in candidates
        update_fitness!(e, candidate)
    end
    return candidates
end

function rank_by_fitness!(e::Evaluator, candidates::AbstractVector{<:Candidate})
    fs = fitness_scheme(e)
    sort!(update_fitness!(e, candidates);
          # FIXME use lt=fitness_scheme(a) when v0.5 #14919 would be fixed
          by=fitness, lt=(x, y) -> is_better(x, y, fs))
end

# called by check_stop_condition(OptRunController)
check_stop_condition(e::Evaluator, ctrl) = check_stop_condition(e.archive, e.problem, ctrl)


text:
"""
    DoubleIntegrator{N,M}

A simple point mass in `M` dimensions with direct control over acceleration. Has `N` state
and `M` controls.

# Constructor
    DoubleIntegrator([D=1])

where `D` is the dimensionality of the space in which the particle moves. For example,
for a particle in 3D space, `D = 3`.
"""
@autodiff struct DoubleIntegrator{N,M} <: ContinuousDynamics
    pos::SVector{M,Int}
    vel::SVector{M,Int}
    gravity::SVector{M,Float64}
end

function DoubleIntegrator(D=1; gravity=zeros(D))
    pos = SVector{D,Int}(1:D)
    vel = SVector{D,Int}(D .+ (1:D))
    DoubleIntegrator{2D,D}(pos,vel, gravity)
end

RobotDynamics.state_dim(::DoubleIntegrator{N,M}) where {N,M} = N
RobotDynamics.control_dim(::DoubleIntegrator{N,M}) where {N,M} = M

@generated function dynamics(di::DoubleIntegrator{N,M}, x, u) where {N,M}
    vel = [:(x[$i]) for i = M+1:N]
    us = [:(u[$i] + di.gravity[$i]) for i = 1:M]
    :(SVector{$N}($(vel...),$(us...)))
end
@generated function dynamics!(di::DoubleIntegrator{N,M}, xdot, x, u) where {N,M}
    vel = [:(xdot[$i-M] = x[$i]) for i = M+1:N]
    us = [:(xdot[$(i+M)] = u[$i] + di.gravity[$i]) for i = 1:M]
    quote
        $(vel...)
        $(us...)
        return nothing
    end
end

Base.position(::DoubleIntegrator{<:Any,2}, x) = @SVector [x[1], x[2], 0]
orientation(::DoubleIntegrator, x) = UnitQuaternion(I)


text:
module HWfuncapp

using FastGaussQuadrature # to get chebyshevnodes

# you dont' have to use PyPlot. I did much of it in PyPlot, hence
# you will see me often qualify code with `PyPlot.plot` or so.
using PyPlot
import ApproXD: getBasis, BSpline
using Distributions
using ApproxFun
using Plots

ChebyT(x,deg) = cos(acos(x)*deg)
unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1	#[a,b] -> [-1,1]

function q1(n)
	f(x) = x .+ 2x.^2 - exp.(-x)

	# without using PyPlot, just erase the `PyPlot.` part
	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png"))
	return Dict(:error=>maximum(abs,err))
end

function q2(b::Number)
	@assert b > 0
	# use ApproxFun.jl to do the same:

	Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q2.png"))
end

function q3(b::Number)

	# p is your plot
	return (p,integral)
end

# optinal
function q4()

	return fig
end

# I found having those useful for q5
mutable struct ChebyType
	f::Function # fuction to approximate
	nodes::Union{Vector,LinRange} # evaluation points
	basis::Matrix # basis evaluated at nodes
	coefs::Vector # estimated coefficients

	deg::Int 	# degree of chebypolynomial
	lb::Float64 # bounds
	ub::Float64

	# constructor
	function ChebyType(_nodes::Union{Vector,LinRange},_deg,_lb,_ub,_f::Function)
		n = length(_nodes)
		y = _f(_nodes)
		_basis = Float64[ChebyT(unitmap(_nodes[i],_lb,_ub),j) for i=1:n,j=0:_deg]
		_coefs = _basis \ y  # type `?\` to find out more about the backslash operator. depending the args given, it performs a different operation
		# create a ChebyComparer with those values
		new(_f,_nodes,_basis,_coefs,_deg,_lb,_ub)
	end
end

# function to predict points using info stored in ChebyType
function predict(Ch::ChebyType,x_new)

	true_new = Ch.f(x_new)
	basis_new = Float64[ChebyT(unitmap(x_new[i],Ch.lb,Ch.ub),j) for i=1:length(x_new),j=0:Ch.deg]
	basis_nodes = Float64[ChebyT(unitmap(Ch.nodes[i],Ch.lb,Ch.ub),j) for i=1:length(Ch.nodes),j=0:Ch.deg]
	preds = basis_new * Ch.coefs
	preds_nodes = basis_nodes * Ch.coefs

	return Dict("x"=> x_new,"truth"=>true_new, "preds"=>preds, "preds_nodes" => preds_nodes)
end

function q5(deg=(5,9,15),lb=-1.0,ub=1.0)

	runge(x) = 1.0 ./ (1 .+ 25 .* x.^2)

	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png"))

end

function q6()

	# compare 2 knot vectors with runge's function

	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q6.png"))

end

function q7()

	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q7.png"))
end

	# function to run all questions
function runall()
	@info("running all questions of HW-funcapprox:")
	q1(15)
	q2(3)
	q3(10)
	q4()
	q5()
	q6()
	q7()
end

end # module