text:
"""
    dediac(s::String)

Dediacritize the input `String` object.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> dediac(verses(crpsdata[1][1])[1])
"bsm {llh {lrHm`n {lrHym"
julia> dediac(arabic(verses(crpsdata[1][1])[1]))
"بسم ٱلله ٱلرحمٰن ٱلرحيم"
```
"""
function dediac(s::String)
    trans = Transliterator()
    isarabic = in(Symbol(s[1]), collect(keys(trans.encode))) ? true : false
    return !isarabic ? replace(s, trans.rx_diacs => s"") :
        replace(s, trans.rx_ardiacs => s"")
end


text:
################################################################################
# Different potentials for molecular dynamics.
################################################################################

module Potentials

type LennardJonesModel
  factor1::Float64
  factor2::Float64
end

"""
Contructor for the LennardJonesModel.
"""
function LennardJonesModel(sigma::Float64, epsilon::Float64)
  factor1::Float64 = epsilon * 24 * 2 * sigma^12
  factor2::Float64 = epsilon * 24 * sigma^6

  return LennardJonesModel(factor1, factor2)
end

"""
Calculating the Lennard-Jones potential force.

m: LennardJonesModel for particles.
r: Distance of particles.
"""
function lennardJonesPotential(m::LennardJonesModel, r::Float64)
  return ((m.factor1 / r^13) - (m.factor2 / r^7))::Float64
end

end

ext: jl

max_stars_repo_name: Christian-Kr/simLJP-julia

max_issues_repo_name: Christian-Kr/simLJP-julia

max_forks_repo_name: Christian-Kr/simLJP-julia


text:
## PART 3

println("\nPart 3");

# Estimate the covariance matrix for these quantities
println("\n\tCalculating the covariance matrices, Cmoldy_1069 and Cmoldy_1304");

# Calculate the mean values of each quantity
mdavg_1069 = zeros(3);
mdavg_1304 = zeros(3);
for n in 1:3
    mdavg_1069[n] = mean(md_1069[:,n]);
    mdavg_1304[n] = mean(md_1304[:,n]);
end

# Calculate the covariance for T = 1.069...
println("\n\tCovariance matrix for T = 1.069:\n")
D = md_1069 - ones(lengthmd) * transpose(mdavg_1069);
Cmoldy_1069 = transpose(D) * D ./ lengthmd
show(Cmoldy_1069)

# ...and for T = 1.304.
println("\n\n\tCovariance matrix for T = 1.304\n")
D = md_1304 - ones(lengthmd) * transpose(mdavg_1304);
Cmoldy_1304 = transpose(D) * D ./ lengthmd;
show(Cmoldy_1304)


text:
######################################################################
# Application Interface
# This is one of many possible ways that devs can use to integrate
# with the Pia framework.
# -----
# Licensed under Apache License 2.0
module Applications
using ExtraFun: Ident
using ..Pia
using ..Packets
import Base.Threads: @spawn

export Application
struct Application <: AbstractApplication
    appname::String
    version::VersionNumber
    protocols::ProtocolStack
end

export appname, appversion
appname(app::AbstractApplication) = app.appname
appversion(app::AbstractApplication) = app.version

export application
function application(name::AbstractString, version::VersionNumber;
        transport::TransportProtocol = PiaLink(),
        cryptography::Optional{CryptographyProtocol} = nothing,
        discovery::DiscoveryProtocol = PiaPortal(32),
        data::Optional{DataProtocol} = nothing,
        return_task::Bool = false,
    )
    app = Application(name, version, ProtocolStack(transport, cryptography, discovery, data))
    for protocol ∈ protocolstack(app)
        Pia.init(protocol, app)
        Pia.restore(protocol)
    end

    task = @spawn Pia.loop(app)
    if return_task
        return app, task
    else
        return app
    end
end
function application(cb, name::AbstractString, version::VersionNumber; kwargs...)
    app, task = application(name, version; kwargs..., return_task=true)
    try
        wait(@async(cb(app)))
        wait(task)
    finally
        for protocol ∈ protocolstack(app)
            Pia.store(protocol)
        end
        close(app)
    end
end

# ===== process_packet =====
# This special method leverages meta programming to achieve dynamic dispatch partially on value rather than type.
# Identified packets are handled. The default fallback forwards packets to waiting futures only if previously registered.
# This solution is the backbone of the secure multiplexed connection.

process_packet(ident::Ident, app::AbstractApplication, ch::PeerChannel, bytes) = process_packet(ident, app, ch.transport, ch.endpoint, bytes)
process_packet(ident::Ident, app::AbstractApplication, tp::TransportProtocol, phys, bytes) = throw(MethodError(process_packet, (ident, app, tp, phys, bytes)))

# --- Handshake Start ---
process_packet(::Ident{Packets.HandshakeStart}, app::AbstractApplication, ch::PeerChannel, bytes) = process_packet(Ident{Packets.HandshakeStart}(), app, ch.transport, ch.endpoint, bytes)

# --- Ping ---
function process_packet(::Ident{Packets.Ping}, ::AbstractApplication, ch::PeerChannel, _)
    if haskey(peers(ch.transport), ch.endpoint)
        @sync begin
            @async write(ch, packetbytes(Packets.Pong))
            @async refreshpeer(ch.transport, ch.endpoint)
        end
    end
end

# --- Goodbye ---
process_packet(::Ident{Packets.Goodbye}, app::AbstractApplication, ch::PeerChannel, bytes) = process_packet(Ident{Packets.Goodbye}(), app, ch.transport, ch.endpoint, bytes)

function Base.close(app::AbstractApplication)
    close(protocolstack(app).transport)
end

end # module Applications


text:
"""
Make all star inference functiouns: prior, loglike, log (unnormalized) posterior.
"""
function make_star_inference_functions(imgs::Vector,
          entry::CatalogEntry;
          pos_delta::Array{Float64, 1}=[2., 2.],
          patches::Array{ImagePatch, 1}=nothing,
          background_images::Array{Array{Float64, 2}, 1}=nothing)
    # position log prior --- same for both star and galaxy (constrains to a
    # small window around the existing catalog location.  Because the range
    # of RA,DEC values is so small, numerical underflow becomes an issue in
    # the slice sampler.  The generic parameter representation of the (ra,dec)
    # position of the sampler will be on [0, 1]^2.  We transform this value
    # before feeding it into the likelihood (but slice sampler machineary will
    # be on the scale of [0, 1]).
    pos_logprior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
        make_location_prior(imgs[1], entry.pos; pos_pixel_delta=pos_delta)

    # star parameters are [lnfluxes, upos]
    loglike = make_star_loglike(imgs;
                                patches=patches,
                                background_images=background_images,
                                pos_transform=uniform_to_deg)

    function logprior(th)
        lnfluxes, pos = th[1:5], uniform_to_deg(th[6:end])
        return logflux_logprior(lnfluxes; is_star=true) +
               pos_logprior(pos)
    end

    function sample_prior()
        lnfluxes = sample_logfluxes(; is_star=true)
        pos      = rand(2)
        return [lnfluxes..., pos...]
    end

    function logpost(th)
        llprior = logprior(th)
        if llprior < -1e100
            return llprior
        end
        return loglike(th) + llprior
    end

    # quick test
    return loglike, logprior, logpost, sample_prior,
           ra_lim, dec_lim, uniform_to_deg, deg_to_uniform
end

"""
Make all star inference functiouns: prior, loglike, log (unnormalized) posterior.
"""
function make_gal_inference_functions(imgs::Vector,
          entry::CatalogEntry;
          pos_delta::Array{Float64, 1}=[2., 2.],
          patches::Array{ImagePatch, 1}=nothing,
          background_images::Array{Array{Float64, 2}, 1}=nothing)
    # constrained location prior (same as above)
    pos_logprior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
        make_location_prior(imgs[1], entry.pos; pos_pixel_delta=pos_delta)

    # gal parameters are [lnfluxes, upos, gal_shape]
    loglike = make_gal_loglike(imgs;
                               patches=patches,
                               background_images=background_images,
                               pos_transform=uniform_to_deg)
    gal_logprior = make_gal_logprior()
    function logprior(th)
        _, pos, _ = th[1:5], uniform_to_deg(th[6:7]), th[8:end]
        return gal_logprior(th) + pos_logprior(pos)
    end

    #th_cat = parameters_from_catalog(entry; is_star=false)
    #cat_shape = th_cat[8:end]
    function sample_prior()
        lnfluxes = sample_logfluxes(; is_star=false)
        pos = rand(2)
        shape = sample_galaxy_shape()
        return [lnfluxes..., pos..., shape...]
    end

    function logpost(th)
        llprior = logprior(th)
        if llprior < -1e100
            return llprior
        end
        return loglike(th) + llprior
    end
    return loglike, logprior, logpost, sample_prior,
           ra_lim, dec_lim, uniform_to_deg, deg_to_uniform
end

"""
Star log likelihood maker.  Creates a star log prob function that is
parameterized by a flat vector

  - th : [log_fluxes; unconstrained_pos]

Args:
  imgs: Array of observed data Images with the .pixel field
  patches...
  background_images...
  pos_transform....
  use_raw_psf: use patch provided raw psf that interpolates on a grid as
    appearance model, not the MOG approximation
"""
function make_star_loglike(imgs::Vector;
                           patches::Array{ImagePatch, 1}=nothing,
                           background_images::Array{Array{Float64, 2}, 1}=nothing,
                           pos_transform::Function=nothing)
    # create background images --- sky noise and neighbors if there
    if background_images == nothing
        background_images = make_empty_background_images(imgs)
    end

    # patch offset (0 if no patches passed in)
    if patches == nothing
        offsets = zeros(2, length(imgs))
        active_bitmaps = [BitArray(ones(img.H, img.W)) for img in imgs]
    else
        offsets = hcat([convert(Array{Float64, 1}, p.bitmap_offset) - 1.
                        for p in patches]...)
        active_bitmaps = [p.active_pixel_bitmap for p in patches]
    end

    # cache the log(data!) term (lgamma in poisson lnpdf) --- this call is
    # a fixed constant wrt params, and only enters into the likelihood
    # as a sum, so we can compute it here and cache it
    lgamma_const = compute_lgamma_sum(imgs, active_bitmaps)

    # create function to return
    function star_loglike(th::Array{Float64, 1})

        # unpack log fluxes and location (ra, dec)
        lnfluxes, pos = th[1:5], th[6:end]
        if pos_transform != nothing
            pos = pos_transform(pos)
        end

        ll = 0.
        for ii in 1:length(imgs)
            img, active_bitmap = imgs[ii], active_bitmaps[ii]

            # sky pixel intensity (sky image)
            background = background_images[ii]

            # band-specific flux --- do bounds check
            bflux = exp(lnfluxes[img.b])
            if isinf(bflux)  # if bflux overflows, then return -Inf logprob
                return -Inf
            end

            # create source flux image
            src_pixels = zeros(Float32, img.H, img.W)
            Model.write_star_nmgy!(pos, bflux, patches[ii], src_pixels;
                                   write_to_patch=true)

            # add background, convert flux to ave elec count
            src_pixels += background
            src_pixels .*= img.nelec_per_nmgy

            # sum per-pixel likelihood contribution
            for h in 1:img.H, w in 1:img.W
                rate_hw    = src_pixels[h, w]
                if isinf(rate_hw)
                    return -Inf
                end
                pixel_data = img.pixels[h,w]
                is_active  = active_bitmap[h, w]
                if !isnan(pixel_data) && is_active
                    ll += (pixel_data*log(rate_hw) - rate_hw)
                    #if isnan(ll)
                    #    println("found nan in star loglike!")
                    #    println("rate: ", rate_hw)
                    #    println("pixel_data: ", pixel_data)
                    #    println("Bflux: ", bflux)
                    #    throw("NOOOOOO")
                    #end
                end
            end
        end
        return ll - lgamma_const
    end

    return star_loglike
end

"""
Galaxy log likelihood maker.  Creates a galaxy log prob function that is
parameterized by a flat vector

  - th : [log_fluxes, unconstrained_pos, unconstrained_gal_shape]

Args:
  imgs: Array of observed data Images with the .pixel field
  pos0: Initial location of the source in ra/dec
  pos_delta: (optional) determines how much ra/dec we allow the sampler
    to drift away from pos0

Note on Galaxy Shapes: the `gal_scale` parameter is in pixels.  The
  Celeste parameterization sigma^2_{cel} is slightly different from the
  Photo and SDSS parameterization for :half_light_radius_px (sigma^2_{sdss}).
  The two have the following relationship:

      sigma^2_{cel}  = sigma^2_{sdss} / sqrt(gal_ab)
   => sigma^2_{sdss} = sigma^2_{cel} * sqrt(gal_ab)

  The function `write_galaxy_unit_flux` is defined for sigma^2_{cel}, that is
  a scale parameter in pixels that indicates the length of the major axis
  (as given by the bivariate_normals.jl#get_bvn_cov function).

  Also note that the scale prior defined below is over sigma^2_{cel}.
"""
function make_gal_loglike(imgs::Vector;
                          patches::Array{ImagePatch, 1}=nothing,
                          background_images::Array{Array{Float64, 2}, 1}=nothing,
                          pos_transform::Function=nothing)

    # create background images --- sky noise and neighbors if there
    if background_images == nothing
        background_images = make_empty_background_images(imgs)
    end

    # patch offset (0 if no patches passed in)
    if patches == nothing
        offsets = zeros(2, length(imgs))
        active_bitmaps = [BitArray(ones(img.H, img.W)) for img in imgs]
    else
        offsets = hcat([convert(Array{Float64, 1}, p.bitmap_offset) - 1.
                        for p in patches]...)
        active_bitmaps = [p.active_pixel_bitmap for p in patches]
    end

    # cache the log(data!) term (lgamma in poisson lnpdf) --- this call is
    # a fixed constant wrt params, and only enters into the likelihood
    # as a sum, so we can compute it here and cache it
    lgamma_const = compute_lgamma_sum(imgs, active_bitmaps)

    function pretty_print_galaxy_params(th)
        lnfluxes, unc_pos, ushape = th[1:5], th[6:7], th[8:end]
        pos   = constrain_pos(unc_pos)
        gal_frac_dev, gal_ab, gal_angle, gal_scale = ushape
        Log.info(@sprintf "lnfluxes     = %s \n" string(lnfluxes))
        Log.info(@sprintf "pos (ra,dec) = %2.4f, %2.4f \n" pos[1] pos[2])
        Log.info(@sprintf "gal shape:\n")
        Log.info(@sprintf "  frac_dev   = %2.4f \n" gal_frac_dev)
        Log.info(@sprintf "  ab ratio   = %2.4f \n" gal_ab)
        Log.info(@sprintf "  angle      = %2.4f \n" gal_angle)
        Log.info(@sprintf "  scale      = %2.4f \n" gal_scale)
    end

    # make galaxy log like function
    function gal_loglike(th::Array{Float64, 1}; print_params=false)

        # unpack location and log fluxes (smushed)
        lnfluxes, pos, ushape = th[1:5], th[6:7], th[8:end]
        gal_frac_dev, gal_ab, gal_angle, gal_scale = ushape
        if pos_transform != nothing
            pos = pos_transform(pos)
        end

        if print_params
          pretty_print_galaxy_params(th)
        end

        ll = 0.
        for ii in 1:length(imgs)
            img, active_bitmap = imgs[ii], active_bitmaps[ii]

            # sky pixel intensity (sky image)
            background = background_images[ii]

            # image specific flux
            bflux = exp(lnfluxes[img.b])
            if isinf(bflux)
                return -Inf
            end

            # create and cache unit flux src image
            src_pixels = zeros(Float32, img.H, img.W)
            Model.write_galaxy_nmgy!(pos, bflux, gal_frac_dev, gal_ab,
                gal_angle, gal_scale, img.psf, [patches[ii]][:,:], src_pixels;
                write_to_patch=true)

            #src_pixels2 = zeros(img.H, img.W)
            #write_galaxy_unit_flux(pos, img.psf, img.wcs, 1.,
            #    gal_frac_dev, gal_ab, gal_angle, gal_scale, src_pixels2; flux=bflux)
            #println("gal is approx", isapprox(src_pixels, src_pixels2))
            #println("  ... rmse", mean( (src_pixels .- src_pixels2).^2 ))
            src_pixels += background
            src_pixels .*= img.nelec_per_nmgy

            # sum per-pixel likelihood contribution
            for h in 1:img.H, w in 1:img.W
                rate_hw = src_pixels[h, w]
                if isinf(rate_hw)
                    return -Inf
                end
                pixel_data = img.pixels[h,w]
                is_active  = active_bitmap[h,w]
                if !isnan(pixel_data) && is_active
                    ll += (pixel_data*log(rate_hw) - rate_hw)
                end
            end
        end
        return ll - lgamma_const
    end

    return gal_loglike
end

##########
# priors #
##########

"""
Create a uniform prior around a RA/DEC location --- specify size of box
by number of pixels
"""
function make_location_prior(img::Image,
                             pos0::Array{Float64, 1};
                             pos_pixel_delta::Array{Float64, 1} = [1., 1.])

    # figure out lower and upper bounds on RA, Dec
    pos0_pix = WCS.world_to_pix(img.wcs, pos0)
    pos0_pix_lower = pos0_pix - .5 * pos_pixel_delta
    pos0_pix_upper = pos0_pix + .5 * pos_pixel_delta
    pos0_world_lower = WCS.pix_to_world(img.wcs, pos0_pix_lower)
    pos0_world_upper = WCS.pix_to_world(img.wcs, pos0_pix_upper)

    # lower and upper bounds on the ra/dec
    ra_lo, ra_hi   = sort([pos0_world_lower[1], pos0_world_upper[1]])
    dec_lo, dec_hi = sort([pos0_world_lower[2], pos0_world_upper[2]])
    Log.info(@sprintf " ... limiting RA  to [%3.5f, %3.5f]" ra_lo ra_hi)
    Log.info(@sprintf " ... limiting DEC to [%3.5f, %3.5f]" dec_lo dec_hi)

    # corresponding uniform log likelihoods
    llra  = log(1./(ra_hi - ra_lo))
    lldec = log(1./(dec_hi - dec_lo))

    function pos_logprior(pos)
        if !inrange(pos[1], ra_lo, ra_hi)
            return -Inf
        end
        if !inrange(pos[2], dec_lo, dec_hi)
            return -Inf
        end
        return llra + lldec
    end

    function uniform_to_deg(u)
        ra  = (ra_hi  - ra_lo ) * u[1] + ra_lo
        dec = (dec_hi - dec_lo) * u[2] + dec_lo
        return [ra, dec]
    end

    function deg_to_uniform(radec)
        u = [ (radec[1] - ra_lo)  / (ra_hi  - ra_lo),
              (radec[2] - dec_lo) / (dec_hi - dec_lo) ]
        return u
    end

    return pos_logprior, [ra_lo, ra_hi], [dec_lo, dec_hi],
           uniform_to_deg, deg_to_uniform
end

function make_gal_logprior()
    # distributions over galaxy parameters
    prior = Model.construct_prior()

    function gal_logprior(th)
        lnfluxes, u, ushape = th[1:5], th[6:7], th[8:end]

        # first compute prior --- make sure all in bounds
        gal_frac_dev, gal_ab, gal_angle, gal_scale = ushape
        if !inrange(gal_frac_dev, 0., 1.)
            return -Inf
        end
        if !inrange(gal_ab, 0., 1.)
            return -Inf
        end
        if !inrange(gal_angle, 0., pi)
            return -Inf
        end
        if !inrange(gal_scale, 1e-5, Inf)
            return -Inf
        end

        # uniform over angle, ll log normal over scale
        llangle = -log(pi)
        llscale = logpdf(prior.galaxy.gal_radius_px, gal_scale)
        if isinf(llangle)
          throw(" angle bad!")
        end
        if isinf(llscale)
          throw(" scale bad!")
        end
        ll = MCMC.logflux_logprior(lnfluxes; is_star=false) + llangle + llscale
        return ll

    end
    return gal_logprior
end

param_prior = Model.construct_prior()

function sample_galaxy_shape()
    gal_frac_dev = rand()
    gal_ab       = rand()
    gal_angle    = rand() * pi
    gal_scale    = rand(param_prior.galaxy.gal_radius_px)
    return [gal_frac_dev, gal_ab, gal_angle, gal_scale]
end

function make_empty_background_images(imgs::Vector)
    background_images = []
    for img in imgs
        # sky pixel intensity (sky image)
        #epsilon    = img.sky[1, 1]
        #iota       = img.nelec_per_nmgy[1]
        #sky_pixels = [epsilon * iota for h=1:img.H, w=1:img.W]
        sky_pixels = img.sky
        Log.info("sky image size : ", size(sky_pixels))
        push!(background_images, sky_pixels)
    end
    return background_images
end

function compute_lgamma_sum(imgs::Vector, active_bitmaps)
    lgamma_const = 0.
    for ii in 1:length(imgs)
        img, active_bitmap = imgs[ii], active_bitmaps[ii]
        for h in 1:img.H, w in 1:img.W
            pixel_data = img.pixels[h,w]
            is_active  = active_bitmap[h,w]
            if !isnan(pixel_data) && is_active
                lgamma_const += lgamma(pixel_data+1.)
            end
        end
    end
    return lgamma_const
end

function make_position_transformations(pos0::Array{Float64, 1},
                                       pos_delta::Array{Float64, 1})

    # create a function that constrains the pixel location to be within
    # pos0 +- [pos_delta]
    pos_lo = pos0 .- pos_delta
    pos_hi = pos0 .+ pos_delta
    pos_range = pos_hi .- pos_lo

    function constrain_pos(upos::Array{Float64, 1})
        unit = 1 ./ (1. + exp.(-upos))
        cpos  = (unit.*pos_range) .+ pos_lo
        return cpos
    end

    function unconstrain_pos(cpos::Array{Float64, 1})
        unit = (cpos .- pos_lo) ./ pos_range
        upos = log.(unit) .- log.(1 .- unit)
        return upos
    end

    return constrain_pos, unconstrain_pos
end

#####################################
# simple log prob implementation    #
#####################################

"""
given a single image, create a function that returns log likelihood of image
as a function of the log flux = [lnu, lnr, lng, lni, lnz]
"""
function make_single_image_logflux_loglike(img0::Image,
                                           pos::Array{Float64, 1};
                                           is_star               = true,
                                           gal_frac_dev::Float64 = .5,
                                           gal_axis_ratio::Float64       = .7,
                                           gal_angle::Float64    = .79,
                                           gal_radius_px::Float64    = 4.)
    # sky pixel intensity
    epsilon    = img0.sky[1, 1]
    iota       = img0.nelec_per_nmgy[1]
    sky_pixels = [epsilon * iota for h=1:img0.H, w=1:img0.W]

    # create and cache unit flux src image
    src_pixels = [0. for h=1:img0.H, w=1:img0.W]
    if is_star
        write_star_unit_flux(img0, pos, src_pixels)
    else
        write_galaxy_unit_flux(img0, pos, gal_frac_dev,
                               gal_axis_ratio, gal_angle, gal_radius_px, src_pixels)
    end

    # create log like handle and return
    function lnflux_loglike(lnflux)

        # compute model image
        rates = sky_pixels .+ exp(lnflux)*src_pixels

        # sum per-pixel likelihood contribution
        ll = 0.
        H, W = size(rates)
        for h in 1:H, w in 1:W
            ll += poisson_lnpdf(img0.pixels[h,w], rates[h,w])
        end

        return ll
    end

    return lnflux_loglike
end

"""
create vectorized star log flux log like function with multiple images
"""
function make_logflux_loglike(imgs::Vector,
                              pos::Array{Float64, 1};
                              is_star = true,
                              gal_frac_dev::Float64 = .5,
                              gal_axis_ratio::Float64       = .7,
                              gal_angle::Float64    = .79,
                              gal_radius_px::Float64    = 4.)

    # create per image scalar log flux function (caches unit flux pixel image)
    if is_star
        scalar_funs = [make_single_image_logflux_loglike(img, pos;
                          is_star=true)
                       for img in imgs]
    else
        scalar_funs = [make_single_image_logflux_loglike(img, pos;
                          is_star=false,
                          gal_frac_dev=gal_frac_dev,
                          gal_axis_ratio=gal_axis_ratio,
                          gal_angle=gal_angle,
                          gal_radius_px=gal_radius_px)
                       for img in imgs]
    end

    function lnflux_loglike(lnflux::Vector{Float64})
        # for each of those image/function pairs, compute flux-specific ll
        ll = 0.
        for i in 1:length(imgs)
            b   = imgs[i].b
            ll += scalar_funs[i](lnflux[b])
        end
        return ll
    end

    return lnflux_loglike
end

#############################
# Priors and prior samplers #
#############################

const pp = Model.load_prior()

"""
prior over log fluxes (for stars and galaxies)
"""
function logflux_logprior(lnfluxes::Vector{Float64}; is_star::Bool=true)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2

    # convert to color space for prior calc
    lnr, colors = logfluxes_to_colors(lnfluxes)

    # compute brightness distribution
    llr = logpdf(Normal(pp.flux_mean[type_i], sqrt(pp.flux_var[type_i])), lnr)

    # compute color likelihoods (mixture model --- computes all component lls)
    nc, nk, ns = size(pp.color_mean)
    llk = [logpdf(MvNormal(pp.color_mean[:,k,type_i],
                           pp.color_cov[:,:,k,type_i]), colors)
           for k in 1:nk]
    lnpik = log.(pp.k[:, type_i])
    llc   = Model.logsumexp(llk .+ lnpik)

    # add and return
    return llr + llc
end

"""
More numerically stable converter from log fluxes to colors
"""
function logfluxes_to_colors(lnfluxes::Vector{Float64})
    lnr = lnfluxes[3]
    colors = Vector{Float64}(4)
    colors[1] = lnfluxes[2] - lnfluxes[1]
    colors[2] = lnfluxes[3] - lnfluxes[2]
    colors[3] = lnfluxes[4] - lnfluxes[3]
    colors[4] = lnfluxes[5] - lnfluxes[4]
    return lnr, colors
end

"""
sample fluxes from prior
"""
function sample_logfluxes(; is_star=true, lnr=nothing)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2

    # sample log r
    if lnr == nothing
        lnr = rand(Normal(pp.flux_mean[type_i], sqrt(pp.flux_var[type_i])))
    end

    # sample colors
    k = rand(Distributions.Categorical(pp.k[:, type_i]))
    c = rand(MvNormal(pp.color_mean[:,k,type_i], pp.color_cov[:,:,k,type_i]))

    # convert to fluxes
    lnfluxes = log.(Model.colors_to_fluxes(lnr, c))
    return lnfluxes
end

"""
sample colors from prior
"""
function sample_colors(; is_star=true)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2
    # sample colors
    k = rand(Distributions.Categorical(pp.k[:, type_i]))
    c = rand(MvNormal(pp.color_mean[:,k,type_i], pp.color_cov[:,:,k,type_i]))
    return c
end

function sample_logr(; is_star=true)
    type_i = is_star ? 1 : 2
    lnr = rand(Normal(pp.flux_mean[type_i], sqrt(pp.flux_var[type_i])))
    return lnr
end

"""
Convert catalog entry into unconstrained parameters for star or gal loglikes
"""
function parameters_from_catalog(entry::CatalogEntry;
                                 unconstrain_pos::Function= x-> x,
                                 is_star=true, epsilon=1e-5)
    if is_star
        return vcat([log.(entry.star_fluxes), unconstrain_pos(entry.pos)]...)
    else
        #ushape = Model.unconstrain_gal_shape([
        #  entry.gal_frac_dev, entry.gal_ab, entry.gal_angle, entry.gal_scale])
        ushape = [clamp(entry.gal_frac_dev, epsilon, 1-epsilon),
                  clamp(entry.gal_axis_ratio, epsilon, 1-epsilon),
                  clamp(entry.gal_angle, epsilon, pi-epsilon),
                  clamp(entry.gal_radius_px, epsilon, Inf)]
        return vcat([log.(entry.gal_fluxes), unconstrain_pos(entry.pos), ushape]...)
    end
end


text:
using OptimalTransport

using Distances
using PythonOT: PythonOT

using Random
using Test

const POT = PythonOT

Random.seed!(100)

@testset "quadratic.jl" begin
    @testset "quadreg" begin
        M = 250
        N = 200

        # create two uniform histograms
        μ = fill(1 / M, M)
        ν = fill(1 / N, N)

        # create random cost matrix
        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

        # compute optimal transport map (Julia implementation + POT)
        eps = 0.25
        γ = quadreg(μ, ν, C, eps, QuadraticOTNewton())

        γ_pot = POT.Smooth.smooth_ot_dual(μ, ν, C, eps; stopThr=1e-9)

        # need to use a larger tolerance here because of a quirk with the POT solver
        @test γ ≈ γ_pot rtol = 1e-1
    end
end


text:
import DataFrames.DataFrame

"""
    polynomial(df::DataFrame, degree::T = 2) where T <: Integer

Calculate [polynomial](https://en.wikipedia.org/wiki/Polynomial_regression) interaction
terms between columns in a `DataFrame`.

If you have a `DataFrame` with 3 columns: `x`, `y`, and `z`, you can get degree-2
polynomial interaction terms: `x*x`, `x*y`, `x*z`, `y*y`, `y*z`, and `z*z`.

# Examples

```julia-repl
julia> using DataFrames

julia> df = DataFrame(a=1:10,b=repeat(0:1,5))
10×2 DataFrame
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     1      0
   2 │     2      1
   3 │     3      0
   4 │     4      1
   5 │     5      0
   6 │     6      1
   7 │     7      0
   8 │     8      1
   9 │     9      0
  10 │    10      1

julia> polynomial(df,2)
10×5 DataFrame
 Row │ a      a_a    a_b    b      b_b
     │ Int64  Int64  Int64  Int64  Int64
─────┼───────────────────────────────────
   1 │     1      1      0      0      0
   2 │     2      4      2      1      1
   3 │     3      9      0      0      0
   4 │     4     16      4      1      1
   5 │     5     25      0      0      0
   6 │     6     36      6      1      1
   7 │     7     49      0      0      0
   8 │     8     64      8      1      1
   9 │     9     81      0      0      0
  10 │    10    100     10      1      1

julia> polynomial(df,3)
10×9 DataFrame
 Row │ a      a_a    a_a_a  a_a_b  a_b    a_b_b  b      b_b    b_b_b
     │ Int64  Int64  Int64  Int64  Int64  Int64  Int64  Int64  Int64
─────┼───────────────────────────────────────────────────────────────
   1 │     1      1      1      0      0      0      0      0  0
   2 │     2      4      8      4      2      2      1      1  1
   3 │     3      9     27      0      0      0      0      0  0
   4 │     4     16     64     16      4      4      1      1  1
   5 │     5     25    125      0      0      0      0      0  0
   6 │     6     36    216     36      6      6      1      1  1
   7 │     7     49    343      0      0      0      0      0  0
   8 │     8     64    512     64      8      8      1      1  1
   9 │     9     81    729      0      0      0      0      0  0
  10 │    10    100   1000    100     10     10      1      1  1
```
"""
function polynomial(df::DataFrame, degree::T = 2) where T <: Integer
    # Check inputs
    if degree < 1
        error("`degree` must be ≥ 1")
    end
    # Build new colnames based on degree
    colnames = names(df)
    last = [[c] for c = colnames]
    new_cols = last[:]
    for _ = 1:(degree-1)
        next = [[l...,c] for l = last for c = colnames]
        append!(new_cols,next)
        last = next
    end
    # Only keep unique columns
    new_cols = Set([sort(cs) for cs = new_cols])
    new_cols = sort([[cs...] for cs = new_cols])
    # Calculate new columns
    rename_cols = cols -> join(cols,"_")
    get_cols = col_names -> (df[!,c] for c = col_names)
    mult_cols = cols -> reduce((l,r) -> l .* r, cols)
    get_mult = col_names -> mult_cols(get_cols(col_names))
    # Return the resulting dataframe
    DataFrame([
        rename_cols(c) => get_mult(c)
        for c = new_cols
    ])
end


text:
"This code defines some simple finite difference approximations, applies them to
approximate u'(x), and compares the result to the exact solution."

# run file by typing "include("fd_simple.jl")", copying and pasting into the REPL, or
# selecting code and pressing Ctrl/Command-Enter

Dfwd(u,x,h) = (u(x+h) - u(x))/h
Dback(u,x,h) = (u(x) - u(x-h))/(h)
D2(u,x,h) = (u(x+h) - u(x-h))/(2*h)

u(x) = sin(1+pi*x)
dudx(x) = pi*cos(1+pi*x)

x̄ = .5
h = .001
@show Dfwd(u,x̄,h) - dudx(x̄)
@show Dback(u,x̄,h) - dudx(x̄)
@show D2(u,x̄,h) - dudx(x̄)


text:
@reexport module _Finance

using Base.Dates: year, month, day
using Base.Dates: yearmonthday, yearmonth, monthday
using Base.Dates: dayofyear, lastdayofmonth, isleapyear

using Roots: fzero

include("daycount.jl")
include("cashflow.jl")

end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

using Parameters

# initialise sim
testCityDataFolder = "data/cities/small/1"
global simDemandCoverage = initSim(joinpath(testCityDataFolder, "sim_config.xml"));

# initialise demand and demand coverage
initDemand!(simDemandCoverage, demandFilename = joinpath(testCityDataFolder, "demand", "demand.csv"))
coverTimes = Dict([p => 20/60/24 for p in priorities])
coverTimes[highPriority] = 12/60/24
initDemandCoverage!(simDemandCoverage, coverTimes = coverTimes, rasterCellNumRows = 2, rasterCellNumCols = 2)

@testset "demand coverage init" begin
	# check sim.demandCoverage after being initialised with function initDemandCoverage!

	# shorthand
	global simDemandCoverage
	sim = simDemandCoverage
	@unpack demand, demandCoverage, stations, numStations, travel = sim;
	@unpack nodesPoints, points, pointsCoverageModes, pointsCoverageModeLookup = demandCoverage;
	numPoints = length(points)
	numNodes = length(sim.net.fGraph.nodes)

	# check that demand of points adds up to demand of raster, for each raster
	# note that this is not the true demand, as it has not been multiplied by a DemandMode.rasterMultiplier value
	for i = 1:demand.numRasters
		@test isapprox(sum(demand.rasters[i].z), sum(demandCoverage.rastersPointDemands[i]), rtol = eps(Float))
	end

	# check demandCoverage.nodesPoints
	@assert(all(i -> points[i].index == i, 1:numPoints))
	@test vcat(nodesPoints...) |> sort == [1:numPoints;] # each point is included exactly once
	@test all(i -> all(j -> points[j].nearestNodeIndex == i, nodesPoints[i]), 1:numNodes) # nodesPoints matches data stored in point.nearestNodeIndex for all points

	# check demandCoverage.pointsCoverageModes
	@test allunique([[pcm.coverTime, pcm.travelMode.index] for pcm in pointsCoverageModes]) # each mode should have a unique combination of coverTime and travelMode
	for (i,pointsCoverageMode) in enumerate(pointsCoverageModes)
		@assert pointsCoverageMode.index == i
		@assert pointsCoverageMode.points == demandCoverage.points

		# will check fields of pointsCoverageMode: pointSets, stationSets, stationsCoverPointSets
		@unpack pointSets, stationSets, stationsCoverPointSets, travelMode = pointsCoverageMode

		@assert(length(pointSets) == length(stationSets))
		@test vcat(pointSets...) |> sort == [1:numPoints;] # each point should be included in exactly one point set
		@test allunique(stationSets) # each set of stations should be unique, otherwise any duplicates should be combined

		stationsCoverPoints = fill(false, numStations, numPoints) # stationsCoverPoints[i,j] is true if station i covers point j
		@assert(all(i -> stations[i].index == i, 1:numStations))
		for (i,station) in enumerate(stations)
			for (j,point) in enumerate(points)
				# calculate travel time from station to point
				time1 = offRoadTravelTime(travelMode, station.nearestNodeDist)
				pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, station.nearestNodeIndex, point.nearestNodeIndex)
				time2 = offRoadTravelTime(travelMode, point.nearestNodeDist)
				travelTime = time1 + pathTravelTime + time2

				# check if station covers point
				if travelTime <= pointsCoverageMode.coverTime
					stationsCoverPoints[i,j] = true
				end
			end
		end

		@test all(i -> all(j -> findall(stationsCoverPoints[:,j]) == stationSets[i], pointSets[i]), 1:length(pointSets)) # assume that both vectors (stationsCoverPoints[:,j] and stationSets[i]) are sorted the same
		@test all(i -> Set(vcat(pointSets[stationsCoverPointSets[i]]...)) == Set(findall(stationsCoverPoints[i,:])), 1:numStations) # check stationsCoverPointSets
	end

	# check pointsCoverageModeLookup
	for pcm in pointsCoverageModes
		@test JEMSS.getPointsCoverageMode(demandCoverage, pcm.travelMode, pcm.coverTime) == pcm
	end
end

@testset "demand coverage calc" begin
	# compare results of fast (using sim.demandCoverage) and slow calculation of demand point coverage

	# shorthand
	global simDemandCoverage
	sim = simDemandCoverage
	@unpack demand, demandCoverage, stations, numStations, travel = sim;
	@unpack points = demandCoverage;
	numPoints = length(points)

	stationsNumAmbs = ones(Int, numStations)
	pointCoverCounts = Vector{Int}(undef, numPoints)
	times = vcat(travel.setsStartTimes, demand.setsStartTimes) |> unique |> sort # simulating to these points in time will cause each combination of travel and demand states to be checked
	for t in times
		# simulateToTime!(sim, t) # probably unnecessary, and sim.time may stop before t
		sim.time = t # is this safe?
		demandsPointSetsCoverCounts = JEMSS.calcPointSetsCoverCounts!(sim, sim.time, stationsNumAmbs)
		for demandPriority in priorities
			travelMode = JEMSS.getTravelMode!(travel, sim.responseTravelPriorities[demandPriority], sim.time)
			coverTime = demandCoverage.coverTimes[demandPriority]
			pointCoverCounts[:] .= 0
			for (i,station) in enumerate(stations)
				for (j,point) in enumerate(points)
					# calculate travel time from station to point
					time1 = offRoadTravelTime(travelMode, station.nearestNodeDist)
					pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, station.nearestNodeIndex, point.nearestNodeIndex)
					time2 = offRoadTravelTime(travelMode, point.nearestNodeDist)
					travelTime = time1 + pathTravelTime + time2

					# check if station covers point
					if travelTime <= coverTime
						pointCoverCounts[j] += stationsNumAmbs[i]
					end
				end
			end

			# compare pointSetsCoverCounts with pointCoverCounts
			pointSetsCoverCounts = demandsPointSetsCoverCounts[demandPriority] # to be checked
			pointsCoverageMode = JEMSS.getPointsCoverageMode(demandCoverage, travelMode, coverTime)
			pointSets = pointsCoverageMode.pointSets # shorthand
			@test all(i -> all(j -> pointCoverCounts[j] == pointSetsCoverCounts[i], pointSets[i]), 1:length(pointSets)) # check that calculated point cover counts are correct

			# also check pointSetsDemands match sum of pointDemands for points in pointSets
			demandMode = JEMSS.getDemandMode!(demand, demandPriority, sim.time)
			rasterIndex = demandMode.rasterIndex # shorthand
			pointSetsDemands = demandCoverage.pointSetsDemands[pointsCoverageMode.index, rasterIndex]
			pointDemands = demandCoverage.rastersPointDemands[rasterIndex]
			@test all(i -> isapprox(sum(j -> pointDemands[j], pointSets[i]), pointSetsDemands[i], rtol = eps(Float)*length(pointSets[i])), 1:length(pointSets)) # check that pointSetsDemands correctly stores the total demand of each pointSets
		end
	end

	# reset sim state, in case sim is needed again
	sim.time = sim.startTime
	travel.recentSetsStartTimesIndex = 1
	demand.recentSetsStartTimesIndex = 1
end


text:
using MieSeries
using Base.Test

p = [1.0, 2.1, -3.7]
K = 4

Y = MieSeries.Ylm(K,p)
Y_from_matlab = [
    0.282094791773878 + 0.000000000000000im
    0.079054029942852 - 0.166013462879989im
    -0.413657340830090 + 0.000000000000000im
    -0.079054029942852 - 0.166013462879989im
    -0.068963090518276 - 0.084939876884680im
    -0.149655973558723 + 0.314277544473318im
    0.362782863277375 + 0.000000000000000im
    0.149655973558723 + 0.314277544473318im
    -0.068963090518276 + 0.084939876884680im
    -0.061128748571790 + 0.014799854825926im
    0.154472358060769 + 0.190259209341709im
    0.191065304223659 - 0.401237138869684im
    -0.184433850338810 + 0.000000000000000im
    -0.191065304223659 - 0.401237138869684im
    0.154472358060769 - 0.190259209341709im
    0.061128748571790 + 0.014799854825926im]

norm(Y-Y_from_matlab)
for i in eachindex(Y)
    @eval @test_approx_eq Y[$i] Y_from_matlab[$i]
end

Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz = MieSeries.Xlm(K,p)

Xx_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.206828670415045 + 0.000000000000000im
    0.000000000000000 - 0.166013462879989im
    -0.206828670415045 + 0.000000000000000im
    -0.061096795363054 + 0.128303270262413im
    0.153237367829497 - 0.034676559530382im
    0.000000000000000 + 0.314277544473318im
    0.153237367829497 + 0.034676559530382im
    0.061096795363054 + 0.128303270262413im
    0.054614225945323 + 0.067266788554357im
    0.065596704579624 - 0.177906321035805im
    -0.021710262651990 + 0.086841050607960im
    0.000000000000000 - 0.401237138869684im
    -0.021710262651990 - 0.086841050607960im
    -0.065596704579624 - 0.177906321035805im
    0.054614225945323 - 0.067266788554357im]

Xy_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.000000000000000 + 0.206828670415045im
    0.000000000000000 + 0.079054029942852im
    0.000000000000000 - 0.206828670415045im
    0.128303270262413 + 0.061096795363054im
    0.034676559530382 - 0.209545495447878im
    0.000000000000000 - 0.149655973558723im
    0.034676559530382 + 0.209545495447878im
    -0.128303270262413 + 0.061096795363054im
    0.067266788554357 - 0.054614225945323im
    -0.188371398743793 - 0.108821257220184im
    -0.086841050607960 + 0.162723587686820im
    0.000000000000000 + 0.191065304223659im
    -0.086841050607960 - 0.162723587686820im
    0.188371398743793 - 0.108821257220184im
    0.067266788554357 + 0.054614225945323im]

Xz_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.055899640652715 + 0.117389245370701im
    0.000000000000000 + 0.000000000000000im
    -0.055899640652715 - 0.117389245370701im
    0.056308127618382 + 0.069353119060764im
    0.061096795363054 - 0.128303270262413im
    0.000000000000000 + 0.000000000000000im
    0.061096795363054 + 0.128303270262413im
    -0.056308127618382 + 0.069353119060764im
    0.052939049164722 - 0.012817050251573im
    -0.089184657508741 - 0.109846205729241im
    -0.055155802413164 + 0.115827185067644im
    0.000000000000000 + 0.000000000000000im
    -0.055155802413164 - 0.115827185067644im
    0.089184657508741 - 0.109846205729241im
    0.052939049164722 + 0.012817050251573im]

@test norm(Xx - Xx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Xy - Xy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Xz - Xz_from_matlab) < sqrt(eps(eltype(p)))

Yx_from_matlab = [
    0.064547345156897 + 0.000000000000000im
    0.018088698925201 - 0.037986267742922im
    -0.094650748379100 + 0.000000000000000im
    -0.018088698925201 - 0.037986267742922im
    -0.015779746867279 - 0.019435465349727im
    -0.034243438949520 + 0.071911221793991im
    0.083009936290289 + 0.000000000000000im
    0.034243438949520 + 0.071911221793991im
    -0.015779746867279 + 0.019435465349727im
    -0.013987136764422 + 0.003386419620560im
    0.035345496987908 + 0.043534043210913im
    0.043718489312336 - 0.091808827555906im
    -0.042201117296689 + 0.000000000000000im
    -0.043718489312336 - 0.091808827555906im
    0.035345496987908 - 0.043534043210913im
    0.013987136764422 + 0.003386419620560im]

Yy_from_matlab = [
    0.135549424829483 + 0.000000000000000im
    0.037986267742922 - 0.079771162260136im
    -0.198766571596109 + 0.000000000000000im
    -0.037986267742922 - 0.079771162260136im
    -0.033137468421285 - 0.040814477234427im
    -0.071911221793991 + 0.151013565767382im
    0.174320866209607 + 0.000000000000000im
    0.071911221793991 + 0.151013565767382im
    -0.033137468421285 + 0.040814477234427im
    -0.029372987205287 + 0.007111481203177im
    0.074225543674607 + 0.091421490742917im
    0.091808827555906 - 0.192798537867402im
    -0.088622346323047 + 0.000000000000000im
    -0.091808827555906 - 0.192798537867402im
    0.074225543674607 - 0.091421490742917im
    0.029372987205287 + 0.007111481203177im]

Yz_from_matlab = [
    -0.238825177080517 + 0.000000000000000im
    -0.066928186023243 + 0.140549190648810im
     0.350207769002669 + 0.000000000000000im
     0.066928186023243 + 0.140549190648810im
     0.058385063408931 + 0.071911221793991im
     0.126700724113223 - 0.266071520637768im
    -0.307136764274069 + 0.000000000000000im
    -0.126700724113223 - 0.266071520637768im
     0.058385063408931 - 0.071911221793991im
     0.051752406028363 - 0.012529752596074im
    -0.130778338855260 - 0.161075959880378im
    -0.161758410455643 + 0.339692661956851im
     0.156144133997750 + 0.000000000000000im
     0.161758410455643 + 0.339692661956851im
    -0.130778338855260 + 0.161075959880378im
    -0.051752406028363 - 0.012529752596074im]

@test norm(Yx - Yx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Yy - Yy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Yz - Yz_from_matlab) < sqrt(eps(eltype(p)))

Zx_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.026860347512988 + 0.231510614278609im
    0.000000000000000 + 0.066928186023243im
    -0.026860347512988 - 0.231510614278609im
    0.135679873108069 + 0.085050235134019im
    0.058715266725129 - 0.239055014523741im
    0.000000000000000 - 0.126700724113223im
    0.058715266725129 + 0.239055014523741im
    -0.135679873108069 + 0.085050235134019im
    0.082386704863588 - 0.052395848498920im
    -0.202331852112205 - 0.144911736140107im
    -0.100023706233593 + 0.193420082720831im
    0.000000000000000 + 0.161758410455644im
    -0.100023706233593 - 0.193420082720831im
    0.202331852112205 - 0.144911736140107im
    0.082386704863588 + 0.052395848498920im]

Zy_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.187894526174186 - 0.026860347512988im
    0.000000000000000 + 0.140549190648810im
    0.187894526174186 + 0.026860347512988im
    0.038841244654719 - 0.124492234448281im
    -0.143712605222459 + 0.058715266725129im
    0.000000000000000 - 0.266071520637768im
    -0.143712605222459 - 0.058715266725129im
    -0.038841244654719 - 0.124492234448281im
    -0.058350340882630 - 0.054016226339494im
    -0.035128304403982 + 0.175752233030384im
    0.031000635931996 - 0.100023706233593im
    0.000000000000000 + 0.339692661956851im
    0.031000635931996 + 0.100023706233593im
    0.035128304403982 + 0.175752233030384im
    -0.058350340882630 + 0.054016226339494im]

Zz_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.099383285798055 + 0.047325374189550im
    0.000000000000000 + 0.097859861185336im
    0.099383285798055 - 0.047325374189550im
    0.058715266725129 - 0.047671204650641im
    -0.065697622768118 - 0.031284582270532im
    0.000000000000000 - 0.185257004716901im
    -0.065697622768118 + 0.031284582270532im
    -0.058715266725129 - 0.047671204650641im
    -0.010851084051334 - 0.044818898327529im
    -0.074621970637991 + 0.060585933303702im
    -0.009438478588217 - 0.004494513613437im
    0.000000000000000 + 0.236517027179738im
    -0.009438478588217 + 0.004494513613437im
    0.074621970637991 + 0.060585933303702im
    -0.010851084051334 + 0.044818898327529im]

@test norm(Zx - Zx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Zy - Zy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Zz - Zz_from_matlab) < sqrt(eps(eltype(p)))


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
## envelope
function get_env(t::TF, pulse::Pulse{TL,TF})::TF where {TL,TF}
    @unpack tL, tR, T, env_type, ncyc = pulse
    if t <= tL || tR <= t
        return 0
    end
    if env_type == :sin2
        return sin((π * t) / (ncyc * T))^2
    elseif env_type == :sin4
        return sin((π * t) / (ncyc * T))^4
    else
        env_type == :flat || env_type == :delta
        return 1
    end
end
function get_env(t::TVF, pulse::Pulse)::TVF where {TVF <: AbstractVector{<:AbstractFloat}}
    env = zeros(eltype(t), size(t))
    for it in eachindex(t)
        env[it] = get_env(t[it], pulse)
    end
    return env
end
## electric field
function getEtAt(t::TF, pulse::Pulse{TL,TF})::Tuple{TF,TF} where {TL,TF}
    return (getEt(t, pulse), getAt(t, pulse))
end
function getEt(t::TF, pulse::Pulse{Et,TF})::TF where {TF}
    @unpack E, ω, ϕ = pulse
    return E * get_env(t, pulse) * sin(ω * t + ϕ)
end
function getEt(t::TF, pulse::Pulse{At,TF})::TF where {TF}
    @unpack tL, tR, T = pulse
    if t <= tL || tR <= t
        return zero(TF)
    end
    dt = T * 1e-5
    Et =
        (1 / 60dt) * (getAt(t + 3dt, pulse) - getAt(t - 3dt, pulse)) -
        (3 / 20dt) * (getAt(t + 2dt, pulse) - getAt(t - 2dt, pulse)) +
        (3 / 4dt) * (getAt(t + dt, pulse) - getAt(t - dt, pulse))
    return -Et
end
function getEt(t::TVF, pulse::Pulse)::TVF where {TVF <: AbstractVector{<:AbstractFloat}}
    Et = zeros(eltype(t), length(t))
    for it in eachindex(t)
        Et[it] = getEt(t[it], pulse)
    end
    return Et
end
## vector potential
function getAt(t::TF, pulse::Pulse{Et,TF})::TF where {TF}
    @unpack tL, tR, T = pulse
    if t <= tL || tR <= t
        return 0
    end
    nt = 1
    tprev = pulse.tprev[1]
    At = pulse.Atprev[1]
    if tprev > t
        tprev = 0
        At = 0
    end
    while (t - tprev) / nt > T / 100000
        nt *= 10
    end
    dt = (t - tprev) / nt
    for it = 1:nt
        t0 = tprev + (it - 1) * dt
        t1 = t0 + dt / 2
        t2 = t0 + dt
        At += (-dt / 6) * (getEt(t0, pulse) + 4getEt(t1, pulse) + getEt(t2, pulse))
    end
    pulse.tprev[1] = t
    pulse.Atprev[1] = At
    return At
end
function getAt(t::TF, pulse::Pulse{At,TF})::TF where {TF}
    @unpack tL, tR, T, A, ω, ϕ = pulse
    if t <= tL || tR <= t
        return 0
    end
    return A * get_env(t, pulse) * sin(ω * t + ϕ)
end
function getAt(t::TVF, pulse::Pulse)::TVF where {TVF <: AbstractVector{<:AbstractFloat}}
    At = zeros(eltype(t), size(t))
    for it = eachindex(t)
        At[it] = getAt(t[it], pulse)
    end
    return At
end


text:
using Juno
using BenchmarkTools

f = (x,y) -> x^2 + y
x = linspace(0.0, 10.0, 500);
Auw = f.(x, x')
A = f.(x, x') .% (2π);
unwrap!(A, (false, false), 1234)
@btime unwrap!($A)
@code_warntype(unwrap!(A))
Juno.@profiler unwrap!(A)

a = 6rand()
b = 6rand()
@btime find_period($a, $b)

f = (x) -> x^2
x = linspace(0.0, 20.0, 1e6);
B = f.(x) .% (2π);
Juno.@profiler unwrap(B)
@btime unwrap!($B)
p1 = TestPixel(0.1)
p2 = TestPixel(0.2)
p3 = TestPixel(0.3)
p4 = TestPixel(0.4)
p5 = TestPixel(0.5)
p6 = TestPixel(0.6)
p7 = TestPixel(0.7)
p8 = TestPixel(0.8)
merge_pixels!(p1, p2, 1)
merge_pixels!(p1, p3, 0)
merge_pixels!(p4, p5, 0)
merge_pixels!(p6, p7, 0)
merge_into_group!(p5, p7, 0)
merge_into_group!(p2, p7, 0)
merge_pixels!(p1, p8, 0)

p8.head.last

max_forks_repo_name: platawiec/Unwrap


text:
using Documenter, Example

include("common.jl")

recursive_replace(@__DIR__, r"\.zh(?<ext>(\.md)?)$" => s"\g<ext>")

makedocs(
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-132847825-3",
        lang = "zh-CN",
        footer = "本站基于 [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl) 和 [Julia 编程语言](https://julialang.org/) 构建，所有内容默认遵循[CC BY 4.0](https://creativecommons.org/licenses/by/4.0/)协议。",
        assets = [
            "assets/custom.css",
        ]
    ),
    sitename="pilgrim",
    pages = [
        "👋 关于" => "index.md",
        "🔗 友链" => "blogroll.md",
        "🗃 研究" => [

        ],
        "💻 技术" => [

        ],
        "📚 杂谈" => [

        ]
    ]
)

max_stars_repo_name: findmyway/Pilgrimygy

max_issues_repo_name: findmyway/Pilgrimygy

max_forks_repo_name: findmyway/Pilgrimygy


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

using Test, SafeTestsets, Distributed

function redirect_devnull(f)
    open(@static(Sys.iswindows() ? "nul" : "/dev/null"), "w") do io
        redirect_stdout(io) do
            f()
        end
    end
end

length(procs()) == 1 && addprocs(1)

@safetestset "Basics"        begin include("test_basics.jl") end
@testset "Actors"            begin include("test_actors.jl") end
@safetestset "Communication" begin include("test_com.jl") end
@safetestset "API"           begin include("test_api.jl") end
@testset "Registry"          begin include("test_registry.jl") end

println("running examples, output suppressed!")
redirect_devnull() do
    @safetestset "Factorial"     begin include("../examples/factorial.jl") end
    @safetestset "Simple"        begin include("../examples/simple.jl") end
    @safetestset "Simple msg."   begin include("../examples/simple_msg.jl") end
    @testset     "Simple distr." begin include("../examples/simple_distr.jl") end
    @testset     "Registry"      begin include("../examples/registry_ex.jl") end
    @safetestset "Stack"         begin include("../examples/stack.jl") end
end


text:
using MolecularDynamicsIntegrators: move!, VerletIntegrator

@testset "Verlet" begin
    Random.seed!(2017012811)

    N = 3
    r = [
        Vector3(0, 0, 0),
        Vector3(0, 1, 0),
        Vector3(0, 1, 1)
    ]
    v = [rand(3) for _ in 1:N]
    dt = 0.001
    mass = i -> 1.0
    force = lj_potential_forces
    steps = 2000

    setup = VerletIntegrator(r, v, dt, mass, force)

    for i = 1:steps
        move!(setup)
    end

    @test vector_equal(setup.positions, [
        [0.012925251233403398, 0.12468008323252931, 0.42136795208804606],
        [0.6118280564741769, 2.5418601972010393, 1.7430189196629002],
        [1.136453017521721, 3.511798489296003, 1.9158143836134611]
    ])
    @test vector_equal(setup.velocities, [
        [-0.006794248644095888, -0.017369573618650585, 0.1533066209776989],
        [0.305914028237079, 0.7709300986005134, 0.8715094598314336],
        [0.581483383021656, 1.335608859882911, 0.5152845468730523]
    ])
end


text:
using DelimitedFiles
using PyPlot
using Distributions
using DataFrames
using CSV
using Random
using Roots

Random.seed!(123)
PyPlot.svg(true)

cpalette10 = [
"#1f77b4"
"#ff7f0e"
"#2ca02c"
"#d62728"
"#9467bd"
"#e377c2"
"#8c564b"
"#7f7f7f"
"#bcbd22"
"#17becf"]

cpalette20 = ["#1f77b4",
 "#aec7e8",
 "#ff7f0e",
 "#ffbb78",
 "#2ca02c",
 "#98df8a",
 "#d62728",
 "#ff9896",
 "#9467bd",
 "#c5b0d5",
 "#8c564b",
 "#c49c94",
 "#e377c2",
 "#f7b6d2",
 "#7f7f7f",
 "#c7c7c7",
 "#bcbd22",
 "#dbdb8d",
 "#17becf",
 "#9edae5"]

include("types/tariffs.jl")
include("types/consumers.jl")
include("types/financials.jl")
include("types/pvsystems.jl")

include("functions/printing_results.jl")
include("functions/balance_energy.jl")
include("functions/balance_energy_nmw.jl")
include("functions/consumer.jl")
include("functions/billing.jl")
include("functions/financials.jl")
include("functions/pvsystem.jl")
include("functions/depictions.jl")
include("functions/tariff.jl")
include("functions/tariffs_hist.jl")
include("functions/printing_tariffs.jl")
include("functions/optimal_pvsystem.jl")
include("functions/retrieve_nsrdb_data.jl")

include("data/tariff_data.jl")
include("data/tariff_data_alternative.jl")
include("data/pvsystem_data.jl")
include("data/consumer_data.jl");

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
using Plots
plotlyjs()

plot([0,6,15,29,53,107])
plot!([0, 6, 12, 18, 24, 30])

max_stars_repo_name: NeuroGameJam/homepage

max_issues_repo_name: NeuroGameJam/homepage

max_forks_repo_name: NeuroGameJam/homepage


text:
# First-order irreversible series reaction
id = :PerezGalvan2017a

pL = [4.5; 0.2]
pU = [0.2; 1.8]
PerezGalvan2017a_x0(p::Vector{T}) where T = T[one(T); zero(T)]
function PerezGalvan2017a_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    du[1] = -p[1]*u[1]
    du[2] = p[1]*u[1] - p[2]*u[2]
    return
end
tspan = (0.0, 1.0)

prob = ODERelaxProb(PerezGalvan2017a_f!, tspan, PerezGalvan2017a_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0098135417300923"
source = "PerezGalvan2017"
desc = "First-order irreversible series reaction"


text:
export RegGrid, gridtopology

"Specification for M by N grid with node degree D"
struct RegGrid{M,N,D}
    width::Float64
end

"Generate nodes of a regular grid"
function gridnodes(grid::RegGrid{M,N,D}) where {M,N,D}
    w = grid.width
    [Node(w*(j - 1), w*(M - i)) for j=1:N for i=1:M]
end

"2d array filled with indices to corresponding 1d array"
nodeidx(m,n) = reshape(1:m*n, m,n)

"generate coordinates for grid neighbors"
function candneighbors(i::Int, j::Int, grid::RegGrid{M,N,4}) where {M,N}
    [          (i-1, j),
     (i, j-1),           (i, j+1),
               (i+1, j)          ]
end
function candneighbors(i::Int, j::Int, grid::RegGrid{M,N,8}) where {M,N}
    [(i-1, j-1), (i-1, j), (i-1, j+1),
     (i, j-1),             (i, j+1),
     (i+1, j-1), (i+1, j), (i+1, j+1)]
end

"valid grid neighbors"
function gridneighbors(i::Int, j::Int, grid::RegGrid{M,N,D}) where {M,N,D}
    filter(t -> t[1] in 1:M && t[2] in 1:N,
           candneighbors(i, j, grid))
end

"Generate (anti-parallel) arcs of a regular grid"
function gridarcs(grid::RegGrid{M,N,D}) where {M,N,D}
    arcs = Arc[]
    idx = nodeidx(M,N)
    for j=1:N
        for i=1:M
            for n in gridneighbors(i, j, grid)
                push!(arcs, Arc(idx[i,j], idx[n...]))
            end
        end
    end
    arcs
end

function gridtopology(grid::RegGrid{M,N,D}) where {M,N,D}
    Topology(gridnodes(grid), gridarcs(grid))
end


text:
#=
    Contents: This file contains code for:
              - Setting the initial values of the biases and weights equal to the final values of a previous run.
	        This is helpful for re-estimating a model on updated training data, where the original and updated training data largely overlap.
	      - Changing the loss function (in our example from Accuracy to ACE)

    Notes:
    1. The model is a toy example with 4 outcomes (categories).
       The model is a poor fit to the data, but this is unimportant. The point of the example is to demonstrate the use of some non-default settings.
    2. For categorical outcomes, use 0-based categories! Some of the loss functions assume this, such as ACE.
    3. Incomplete batches are padded with repeated instances of an artificial observation.
       This is bad because the artificial data is over-represented and thus biases the results.
       The ideal solution is to distribute the observations from the incomplete batch among the complete batches.
       This would result in batches of variable but similar size, and thus the estimate of the gradient would not be significantly affected.
       But this doesn't happen.
       For simplicity we instead drop these extra observations, so that the number of observations in the data set is a multiple of the batch_size.
=#

using RDatasets
using MXNet

################################################################################
### Data: Exam scores discretised into 4 categories (use zero-based categories!).
df = dataset("mlmRev", "Gcsemv");    # 1905 x 5
complete_cases!(df)                  # 1523 x 5
n = nrow(df)
df[:written] = zeros(Int, n)
df[:course]  = zeros(Int, n)
for i = 1:n
    # Categorise :Written
    if df[i, :Written] <= 20.0
	df[i, :written] = 0
    elseif df[i, :Written] <= 40.0
	df[i, :written] = 1
    elseif df[i, :Written] <= 60.0
	df[i, :written] = 2
    else
	df[i, :written] = 3
    end

    # Categorise :Course
    if df[i, :Course] <= 25.0
	df[i, :course] = 0
    elseif df[i, :Course] <= 50.0
	df[i, :course] = 1
    elseif df[i, :Course] <= 75.0
	df[i, :course] = 2
    else
	df[i, :course] = 3
    end
end
df = df[1:1500, :]    # Ensure nrows is a multiple of batch_size (100 in our example, see below)

x = convert(Vector{Float64}, df[:course])
y = convert(Vector{Float64}, df[:written])

################################################################################
### Hyperparameters

# Architecture
mlp = @mx.chain mx.Variable(:data) =>
        mx.FullyConnected(name = :h1, num_hidden = 10) =>
	mx.Activation(name = :h1_out, act_type = :sigmoid) =>
        mx.FullyConnected(name = :out, num_hidden = 4) =>
	mx.SoftmaxOutput(name = :softmax)

# Hyperparameters
n_epoch    = 100
batch_size = 100
learn_rate = 0.1
mom        = 0.9
wt_decay   = 0.00001

# Connect data, network architecture and hyperparameters
train_prov = mx.ArrayDataProvider(x, y; batch_size = batch_size)
eval_prov  = mx.ArrayDataProvider(x, y; batch_size = batch_size)
opt        = mx.SGD(lr = learn_rate, momentum = mom, weight_decay = wt_decay)    # Optimizing algorithm

################################################################################
### Run 1: Basic run, storing initial and final state.

# Learn
mdl1 = mx.FeedForward(mlp, context = mx.cpu())                                               # Model targets the local CPU
cb = mx.do_checkpoint("first", frequency = n_epoch, save_epoch_0 = true)                     # Write initial and final states to disk
mx.fit(mdl1, opt, train_prov, n_epoch = n_epoch, eval_data = eval_prov, callbacks = [cb])    # Random initial biases and weights

################################################################################
### Run 2: Load the previously trained model and run it some more, starting where Run 1 finished.

# Load final state of 1st run from disk
arch, arg_params, aux_params = mx.load_checkpoint("first", 100)    # arch is the network structure, arg_params contains the weights and biases
mdl2 = mx.FeedForward(arch, context = mx.cpu())                    # Only populates the arch and ctx fields
mdl2.arg_params = arg_params                                       # Populate the arg_params fields
cb   = mx.do_checkpoint("second", frequency = n_epoch, save_epoch_0 = true)
mx.fit(mdl2, opt, train_prov, n_epoch = n_epoch, eval_data = eval_prov, callbacks = [cb])

# Test whether the final state of 1st run equals the initial state of 2nd run
run(`diff first-0100.params second-0000.params`)    # Throws error if not true, does nothing otherwise

#=
    # Other useful functions
    arch       = mx.load("first-symbol.json", mx.SymbolicNode)
    arg_params = mx.load("first-0100.params", mx.NDArray)
=#

################################################################################
### Run 3: Change the loss function from the default Accuracy to ACE

mdl3 = mx.FeedForward(mlp, context = mx.cpu())
mx.fit(mdl3, opt, train_prov, n_epoch = n_epoch, eval_data = eval_prov, eval_metric = mx.ACE())
#mx.fit(mdl3, opt, train_prov, n_epoch = n_epoch, eval_data = eval_prov, eval_metric = mx.Accuracy())    # Default eval_metric
#mx.fit(mdl3, opt, train_prov, n_epoch = n_epoch, eval_data = eval_prov, eval_metric = mx.MultiACE(4))

# Test manually
probs = mx.predict(mdl3, eval_prov)
LL    = 0.0
for i = 1:size(y, 1)
    LL += log(probs[Int(y[i]) + 1, i])
end
-LL / size(y, 1)    # Should equal the value of ACE from the final iteration of fit(mdl3, ...)

# EOF


text:
module test_bukdu_controllers

using Test # @test_throws
using Bukdu # routes get ApplicationController Routing
import Bukdu.Actions: index

@test_throws Routing.AbstractControllerError routes(() -> get("/", ApplicationController, index))

end # module test_bukdu_controllers

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BasicDQN},
    ::Val{:EmptyRoom},
    ::Nothing;
    seed = 123,
    save_dir = nothing,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_BasicDQN_EmptyRoom$(t)")
    end
    log_dir = joinpath(save_dir, "tb_log")
    lg = TBLogger(log_dir, min_level = Logging.Info)
    rng = StableRNG(seed)

    inner_env = GridWorlds.EmptyRoomDirected(rng = rng)
    action_space_mapping = x -> Base.OneTo(length(RLBase.action_space(inner_env)))
    action_mapping = i -> RLBase.action_space(inner_env)[i]
    env = RLEnvs.ActionTransformedEnv(
        inner_env,
        action_space_mapping = action_space_mapping,
        action_mapping = action_mapping,
    )
    env = RLEnvs.StateOverriddenEnv(env, x -> vec(Float32.(x)))
    env = RewardOverriddenEnv(env, x -> x - convert(typeof(x), 0.01))
    env = MaxTimeoutEnv(env, 240)

    ns, na = length(state(env)), length(action_space(env))
    agent = Agent(
        policy = QBasedPolicy(
            learner = BasicDQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = Chain(
                        Dense(ns, 128, relu; initW = glorot_uniform(rng)),
                        Dense(128, 128, relu; initW = glorot_uniform(rng)),
                        Dense(128, na; initW = glorot_uniform(rng)),
                    ) |> cpu,
                    optimizer = ADAM(),
                ),
                batch_size = 32,
                min_replay_history = 100,
                loss_func = huber_loss,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000)

    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNStep() do t, agent, env
            with_logger(lg) do
                @info "training" loss = agent.policy.learner.loss
            end
        end,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                @info "training" reward = total_reward_per_episode.rewards[end] log_step_increment =
                    0
            end
        end,
    )

    description = """
    This experiment uses three dense layers to approximate the Q value.
    The testing environment is EmptyRoom.

    You can view the runtime logs with `tensorboard --logdir $log_dir`.
    Some useful statistics are stored in the `hook` field of this experiment.
    """

    Experiment(agent, env, stop_condition, hook, description)
end


text:
# This file manages the tests for the `Exceptions.jl` package

# A flag to throw the first captured exception back
fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"

# A flag to hide additional LoadErrors
quiet = length(ARGS) > 0 && ARGS[1] == "-q"

# A flag to check if any error occurred
anyerrors = false

# Specify tests
tests_dir = joinpath(@__DIR__, "tests")
tests = ["exception.jl"]

# Print info
println("\e[1;32mRUNNING TESTS\e[0m for `Exceptions.jl`")

# Run tests
for test in joinpath.(tests_dir, tests)

    # Run a test
    try

        include(test)

        # If everything is fine, print about that
        println("\e[1;32mPASSED\e[0m: $(basename(test))")

    # If error occurred, note that
    catch e

        # Alter the global variable
        global anyerrors = true

        # Print about the fail
        println("\e[1;31mFAILED\e[0m: $(basename(test))")

        # Exit immediately is needed
        if fatalerrors
            rethrow(e)
        # Show the error if needed
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end

    end
end

# If any errors, throw an exception
anyerrors && throw("Some tests have failed.")


text:
include("../../src/preprocessing/one_hot_encoder.jl")
using Test, DataFrames

test_data = DataFrame(
    x = ["one", "two", "three", "four"],
    y = [1.0, 2.0, 2.0, 2.0]
)

@testset "one_hot_encoder" begin
    # string column
    @test isequal(one_hot_encoder(test_data, "x"), DataFrame(
        y = [1.0, 2.0, 2.0, 2.0],
        x_one = [1, 0, 0, 0],
        x_two = [0, 1, 0, 0],
        x_three = [0, 0, 1, 0],
        x_four = [0, 0, 0, 1]
    ))
    # symbol column
    @test isequal(one_hot_encoder(test_data, :x), DataFrame(
        y = [1.0, 2.0, 2.0, 2.0],
        x_one = [1, 0, 0, 0],
        x_two = [0, 1, 0, 0],
        x_three = [0, 0, 1, 0],
        x_four = [0, 0, 0, 1]
    ))
    # different prefix
    @test isequal(one_hot_encoder(test_data, "x", prefix="new"), DataFrame(
        y = [1.0, 2.0, 2.0, 2.0],
        new_one = [1, 0, 0, 0],
        new_two = [0, 1, 0, 0],
        new_three = [0, 0, 1, 0],
        new_four = [0, 0, 0, 1]
    ))
    # decimals in values
    @test isequal(
        one_hot_encoder(test_data, "y", prefix="new"),
        DataFrame(
            x = ["one", "two", "three", "four"],
            one = [1, 0, 0, 0],
            two = [0, 1, 1, 1]
        ) |> x -> rename!(x, [:x, Symbol("new_1.0"), Symbol("new_2.0")])
    )
    # remove original = false
    @test isequal(one_hot_encoder(test_data, "x", remove_original=false), DataFrame(
        x = ["one", "two", "three", "four"],
        y = [1.0, 2.0, 2.0, 2.0],
        x_one = [1, 0, 0, 0],
        x_two = [0, 1, 0, 0],
        x_three = [0, 0, 1, 0],
        x_four = [0, 0, 0, 1]
    ))
end


text:
push!(LOAD_PATH, "../src/")

using Documenter, CJKFrequencies

makedocs(
    sitename="CJKFrequencies.jl Documentation",
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    ),
    modules=[CJKFrequencies],
    pages=[
        "Home" => "index.md",
        "API Reference" => "api_reference.md",
        "Developer Docs" => "devdocs.md"
    ]
)

deploydocs(
    repo = "github.com/tmthyln/CJKFrequencies.jl.git",
    devbranch = "main",
    devurl="latest"
    )


text:
# Invertible conditional HINT layer from Kruse et al. (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export ConditionalLayerHINT

"""
    CH = ConditionalLayerHINT(nx, ny, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, permute=true) (2D)

    CH = ConditionalLayerHINT(nx, ny, nz, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, permute=true) (3D)

 Create a conditional HINT layer based on coupling blocks and 1 level recursion.

 *Input*:

 - `nx`, `ny`, `nz`: spatial dimensions of both `X` and `Y`.

 - `n_in`, `n_hidden`: number of input and hidden channels of both `X` and `Y`

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 - `permute`: bool to indicate whether to permute `X` and `Y`. Default is `true`

 *Output*:

 - `CH`: Conditional HINT coupling layer.

 *Usage:*

 - Forward mode: `Zx, Zy, logdet = CH.forward_X(X, Y)`

 - Inverse mode: `X, Y = CH.inverse(Zx, Zy)`

 - Backward mode: `ΔX, ΔY, X, Y = CH.backward(ΔZx, ΔZy, Zx, Zy)`

 - Forward mode Y: `Zy = CH.forward_Y(Y)`

 - Inverse mode Y: `Y = CH.inverse(Zy)`

 *Trainable parameters:*

 - None in `CH` itself

 - Trainable parameters in coupling layers `CH.CL_X`, `CH.CL_Y`, `CH.CL_YX` and in
   permutation layers `CH.C_X` and `CH.C_Y`.

 See also: [`CouplingLayerBasic`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct ConditionalLayerHINT <: NeuralNetLayer
    CL_X::CouplingLayerHINT
    CL_Y::CouplingLayerHINT
    CL_YX::CouplingLayerBasic
    C_X::Union{Conv1x1, Nothing}
    C_Y::Union{Conv1x1, Nothing}
    logdet::Bool
    is_reversed::Bool
end

@Flux.functor ConditionalLayerHINT

# 2D Constructor from input dimensions
function ConditionalLayerHINT(nx::Int64, ny::Int64, n_in::Int64, n_hidden::Int64, batchsize::Int64;
                              k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet=true, permute=true)

    # Create basic coupling layers
    CL_X = CouplingLayerHINT(nx, ny, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet, permute="none")
    CL_Y = CouplingLayerHINT(nx, ny, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet, permute="none")
    CL_YX = CouplingLayerBasic(nx, ny, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet)

    # Permutation using 1x1 convolution
    permute == true ? (C_X = Conv1x1(n_in)) : (C_X = nothing)
    permute == true ? (C_Y = Conv1x1(n_in)) : (C_Y = nothing)

    return ConditionalLayerHINT(CL_X, CL_Y, CL_YX, C_X, C_Y, logdet, false)
end

# 3D Constructor from input dimensions
function ConditionalLayerHINT(nx::Int64, ny::Int64, nz:: Int64, n_in::Int64, n_hidden::Int64, batchsize::Int64;
                              k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet=true, permute=true)

    # Create basic coupling layers
    CL_X = CouplingLayerHINT(nx, ny, nz, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet, permute="none")
    CL_Y = CouplingLayerHINT(nx, ny, nz, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet, permute="none")
    CL_YX = CouplingLayerBasic(nx, ny, nz, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet)

    # Permutation using 1x1 convolution
    permute == true ? (C_X = Conv1x1(n_in)) : (C_X = nothing)
    permute == true ? (C_Y = Conv1x1(n_in)) : (C_Y = nothing)

    return ConditionalLayerHINT(CL_X, CL_Y, CL_YX, C_X, C_Y, logdet, false)
end

function forward(X, Y, CH::ConditionalLayerHINT; logdet=nothing, x_lane=false)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    # Y-lane
    ~isnothing(CH.C_Y) ? (Yp = CH.C_Y.forward(Y)) : (Yp = copy(Y))
    logdet ? (Zy, logdet1) = CH.CL_Y.forward(Yp) : Zy = CH.CL_Y.forward(Yp)

    # X-lane: coupling layer
    ~isnothing(CH.C_X) ? (Xp = CH.C_X.forward(X)) : (Xp = copy(X))
    logdet ? (X, logdet2) = CH.CL_X.forward(Xp) : X = CH.CL_X.forward(Xp)

    # X-lane: conditional layer
    logdet ? (Zx, logdet3) = CH.CL_YX.forward(Yp, X)[2:3] : Zx = CH.CL_YX.forward(Yp, X)[2]

    if x_lane == false
        logdet ? (return Zx, Zy, logdet1 + logdet2 + logdet3) : (return Zx, Zy)
    else
        logdet ? (return Zx, Zy, logdet2 + logdet3) : (return Zx, Zy)
    end
end

function inverse(Zx, Zy, CH::ConditionalLayerHINT; logdet=nothing, x_lane=false)
    isnothing(logdet) ? logdet = (CH.logdet && CH.is_reversed) : logdet = logdet

    # Y-lane
    logdet ? (Yp, logdet1) = CH.CL_Y.inverse(Zy; logdet=true) : Yp = CH.CL_Y.inverse(Zy; logdet=false)
    ~isnothing(CH.C_Y) ? (Y = CH.C_Y.inverse(Yp)) : (Y = copy(Yp))

    # X-lane: conditional layer
    YZ = tensor_cat(Yp, Zx)
    logdet ? (X, logdet2) = CH.CL_YX.inverse(Yp, Zx; logdet=true)[2:3] : X = CH.CL_YX.inverse(Yp, Zx)[2]

    # X-lane: coupling layer
    logdet ? (Xp, logdet3) = CH.CL_X.inverse(X; logdet=true) : Xp = CH.CL_X.inverse(X; logdet=false)
    ~isnothing(CH.C_X) ? (X = CH.C_X.inverse(Xp)) : (X = copy(Xp))

    if x_lane == false
        logdet ? (return X, Y, logdet1 + logdet2 + logdet3) : (return X, Y)
    else
        logdet ? (return X, Y, logdet2 + logdet3) : (return X, Y)
    end
end

function backward(ΔZx, ΔZy, Zx, Zy, CH::ConditionalLayerHINT)

    # Y-lane
    ΔYp, Yp = CH.CL_Y.backward(ΔZy, Zy)

    # X-lane: conditional layer
    ΔYp_, ΔX, X = CH.CL_YX.backward(ΔYp.*0f0, ΔZx, Yp, Zx)[[1,2,4]]
    ΔYp += ΔYp_

    # X-lane: coupling layer
    ΔXp, Xp = CH.CL_X.backward(ΔX, X)

    # 1x1 Convolutions
    if isnothing(CH.C_X) || isnothing(CH.C_Y)
        ΔX = copy(ΔXp); X = copy(Xp)
        ΔY = copy(ΔYp); Y = copy(Yp)
    else
        ΔX, X = CH.C_X.inverse((ΔXp, Xp))
        ΔY, Y = CH.C_Y.inverse((ΔYp, Yp))
    end
    return ΔX, ΔY, X, Y
end

function backward_inv(ΔX, ΔY, X, Y, CH::ConditionalLayerHINT)

    # 1x1 Convolutions
    if isnothing(CH.C_X) || isnothing(CH.C_Y)
        ΔXp = copy(ΔX); Xp = copy(X)
        ΔYp = copy(ΔY); Yp = copy(Y)
    else
        ΔXp, Xp = CH.C_X.forward((ΔX, X))
        ΔYp, Yp = CH.C_Y.forward((ΔY, Y))
    end

    # X-lane: coupling layer
    ΔX, X = backward_inv(ΔXp, Xp, CH.CL_X)

    # X-lane: conditional layer
    ΔYp_, ΔZx, Zx = backward_inv(ΔYp.*0f0, ΔX, Yp, X, CH.CL_YX)[[1,2,4]]
    ΔYp += ΔYp_

    # Y-lane
    ΔZy, Zy = backward_inv(ΔYp, Yp, CH.CL_Y)

    return ΔZx, ΔZy, Zx, Zy
end

function forward_Y(Y, CH::ConditionalLayerHINT)
    ~isnothing(CH.C_Y) ? (Yp = CH.C_Y.forward(Y)) : (Yp = copy(Y))
    Zy = CH.CL_Y.forward(Yp; logdet=false)
    return Zy

end

function inverse_Y(Zy, CH::ConditionalLayerHINT)
    Yp = CH.CL_Y.inverse(Zy; logdet=false)
    ~isnothing(CH.C_Y) ? (Y = CH.C_Y.inverse(Yp)) : (Y = copy(Yp))
    return Y
end

# Clear gradients
function clear_grad!(CH::ConditionalLayerHINT)
    clear_grad!(CH.CL_X)
    clear_grad!(CH.CL_Y)
    clear_grad!(CH.CL_YX)
    ~isnothing(CH.C_X) && clear_grad!(CH.C_X)
    ~isnothing(CH.C_Y) && clear_grad!(CH.C_Y)
end

# Get parameters
function get_params(CH::ConditionalLayerHINT)
    p = get_params(CH.CL_X)
    p = cat(p, get_params(CH.CL_Y); dims=1)
    p = cat(p, get_params(CH.CL_YX); dims=1)
    ~isnothing(CH.C_X) && (p = cat(p, get_params(CH.C_X); dims=1))
    ~isnothing(CH.C_Y) && (p = cat(p, get_params(CH.C_Y); dims=1))
    return p
end

# Put parameters
function put_params!(CH::ConditionalLayerHINT, Params::Array{Any,1})
    idx_s = 1
    counter = 0

    nlayers_x = length(CH.CL_X.CL)
    isnothing(CH.CL_X.C) ? counter += 5*nlayers_x : counter += 5*nlayers_x + 3
    put_params!(CH.CL_X, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 0
    nlayers_y = length(CH.CL_Y.CL)
    isnothing(CH.CL_Y.C) ? counter += 5*nlayers_y : counter += 5*nlayers_y + 3
    put_params!(CH.CL_Y, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 5
    put_params!(CH.CL_YX, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 3
    ~isnothing(CH.C_X) && put_params!(CH.C_X, Params[idx_s:idx_s+counter-1]); idx_s += counter
    ~isnothing(CH.C_Y) && put_params!(CH.C_Y, Params[idx_s:idx_s+counter-1]); idx_s += counter
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(H::ConditionalLayerHINT, tag::Bool)
    H.is_reversed = tag
    tag_as_reversed!(H.CL_X, tag)
    tag_as_reversed!(H.CL_Y, tag)
    tag_as_reversed!(H.CL_YX, tag)
    return H
end


text:
"""
    Supplier(order_cost, onorder::NumDist, leadtime::Int = 0, capacity::number = Inf; name = "")
"""
mutable struct Supplier{F,L <:LeadTime}
    order_cost::F
    capacity::Float64
    leadtime::L
    pull_orders::IdDict{Any, Float64}
    batchsize_log::Vector{Float64}
    cost_log::Vector{Float64}
    name::String
end

function Supplier(order_cost; leadtime = leadtime::LeadTime = LeadTime(0, 0), capacity::Number = Inf, name = "supplier")
    @assert hasmethod(order_cost, Tuple{Supplier}) "order cost must have a method with `(::Supplier) signature`"
    Supplier(order_cost, Float64(capacity), leadtime, IdDict{Any,Float64}(), zeros(0), zeros(0), name)
end

function Supplier(fixed_order_cost::Number, variable_order_cost::Number; leadtime::LeadTime = LeadTime(0, 0), capacity::Number = Inf, name = "supplier")
    Supplier(FixedLinearOrderCost(fixed_order_cost,variable_order_cost), leadtime = leadtime, capacity = capacity, name = name)
end

RLBase.state(s::Supplier) = state(s.leadtime)
state_size(s::Supplier) = state_size(s.leadtime)
print_state(s::Supplier) = print_state(s.leadtime)

function pull!(sup::Supplier, quantity::Number, issuer) #add leadtime
    @assert 0 == length(sup.pull_orders) "supplier received multiple orders"
    push!(sup.pull_orders, issuer => min(quantity, sup.capacity))
    nothing
end

function dispatch!(sup::Supplier)
    destination, quantity = only(sup.pull_orders)
    push!(sup.batchsize_log, quantity)
    push!(sup.leadtime, quantity, destination)
    dispatch!(sup.leadtime)
    nothing
end

function reward!(sup::Supplier)
    cost = sup.order_cost(sup)
    push!(sup.cost_log, cost)
    empty!(sup.pull_orders)
    return -cost + reward!(sup.leadtime)
end

inventory_position(sup::Supplier) = inventory_position(sup.leadtime)

function RLBase.reset!(s::Supplier)
    empty!(s.cost_log)
    empty!(s.batchsize_log)
    reset!(s.leadtime)
end

children(sup::Supplier) = ()

mutable struct FixedLinearOrderCost{T1,T2}
    K::T1
    c::T2
end
function (f::FixedLinearOrderCost)(sup::Supplier)
    q = sum(values(sup.pull_orders))
    return (q > 0)*(f.K + f.c*q)
end

mutable struct LinearOrderCost{T2}
    c::T2
end
(f::LinearOrderCost)(sup::Supplier) = f.c*sum(values(sup.pull_orders))

Base.show(io::IO, sup::Supplier{F}) where {F} = print(io, "Supplier($(sup.name), $F, $(sup.leadtime.leadtime))")

max_stars_repo_path: src/elementary blocks/supplier.jl

max_issues_repo_path: src/elementary blocks/supplier.jl

max_forks_repo_path: src/elementary blocks/supplier.jl


text:
function foo()
    a = 10
    b = 20
    c1 = (a + b)
    c2 = (a - b)
    c3 = (a * b)
    c4 = (a / b)
    c5 = -(a)
    d = 2.0
    e1 = (a + d)
    e2 = (a - d)
    e3 = (a * d)
    e4 = (a / d)
    f = -3.0
    g = -(a)
end

function add1(x::Int8, y::Int8)::Int16
    return (x + y)
end

function add2(x::Int16, y::Int16)::Int32
    return (x + y)
end

function add3(x::Int32, y::Int32)::Int64
    return (x + y)
end

function add4(x::Int64, y::Int64)::Int64
    return (x + y)
end

function add5(x::UInt8, y::UInt8)::UInt16
    return (x + y)
end

function add6(x::UInt16, y::UInt16)::UInt32
    return (x + y)
end

function add7(x::UInt32, y::UInt32)::UInt64
    return (x + y)
end

function add8(x::UInt64, y::UInt64)::UInt64
    return (x + y)
end

function add9(x::Int8, y::UInt16)::UInt32
    return (x + y)
end

function sub(x::Int8, y::Int8)::Int8
    return (x - y)
end

function mul(x::Int8, y::Int8)::Int16
    return (x * y)
end

function fadd(x::Int8, y::Float64)::Float64
    return (x + y)
end


text:
## Linear programming in MatrixOptim.jl

"Solution for a linear programming problem."
struct SolutionLinear
    obj::Float64
    x::Array{Float64,2}
end

"Solution for a mixed integer programming problem."
struct SolutionMix
    obj::Float64
    x::Array{Float64,2}
    y::Array{Float64,2}
end

@doc raw"""
Constructor for mixed integer linear programming (MILP) problems.

Any MILP can be written in the following standard format:
`` \begin{aligned}
\min \quad & \mathbf{c}^T \mathbf{x} + \mathbf{f}^T \mathbf{y} \\
\text{s.t.} \quad & \mathbf{A} \mathbf{x} + \mathbf{D} \mathbf{y} \geq \mathbf{b} \\
& \mathbf{y} \in Y \\
& \mathbf{x} \geq 0
\end{aligned} ``

# Arguments

- `c`: coefficient vector for continuous variables in objective function.
- `a`: coefficient matrix for continuous variables in constraints.
- `b`: coefficient matrix in constraints.
- `f`: coefficient vector for integer variables in objective function.
- `d`: coefficient matrix for integer variables in constraints.
"""
struct ModelMilp
    c::Array{Float64,2}
    a::Array{Float64,2}
    b::Array{Float64,2}
    f::Union{Array{Float64,2}, Missing}
    d::Union{Array{Float64,2}, Missing}
    solution::Union{SolutionMix, Missing}

    function ModelMilp(c, a, b, f=missing, d=missing)
        check_col_vec(c, "c")
        check_col_vec(f, "f")
        check_col_vec(b, "b")
        new(c, a, f, d, b, missing)
    end
end

"Solve the model."
function solve!(model::Union{ModelLinear, ModelMix})::Union{SolutionLinear, SolutionMix}
    if isa(model, ModLinear)
        obj, x, u = solve_lp(model.c, model.b, model.a)
    else
        obj, x, u = solve_milp(model.c, model.b, model.a, model.f, model.b)
    end
    model.solution = Solution(obj, hcat(x), hcat(u))

    return solution
end

"""
    solve_lp(n_x, vec_c, vec_b, mat_aCap)

Solve a linear programming problem in its matrix form.

# Arguments
- vec_c: [column vector] coefficient vector in objective function
- vec_b: [column vector] right-hand-side coefficient vector in constraint
- mat_aCap: Coefficient Matrix
"""
function solve_lp(vec_c, vec_b, mat_aCap)
    n_x = length(vec_c)
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model, vec_x[1: n_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    @constraint(model, vec_cons, mat_aCap * vec_x .>= vec_b)
    optimize!(model)
    vec_result_u = dual_vec(vec_cons)
    obj = objective_value(model)
    vec_result_x =  value_vec(vec_x)
    return obj, vec_result_x, vec_result_u
end

function solve_milp(n_x, vec_c, vec_b, mat_aCap, vec_f, mat_b)
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model, vec_x[1: n_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    @constraint(model, vec_cons, mat_aCap * vec_x .>= vec_b)
    optimize!(model)
    vec_result_u = dual_vec(vec_cons)
    obj = objective_value(model)
    vec_result_x = value_vec(vec_x)
    return obj, vec_result_x, vec_result_u
end

"Append a constraint"
function append_cons!(
        mod::ModelMix,
        vec_a_new::Array{Int64,2},
        vec_b_new::Array{Int64,2},
        b_new::Float64
    )
    check_col_vec(vec_a_new, "vec_a_new")
    check_col_vec(vec_b_new, "vec_b_new")
    mod.mat_aa = vcat(mod.mat_aa, vec_a_new')
    mod.mat_bb = vcat(mod.mat_bb, vec_b_new')
    mod.vec_b = vcat(mode.vec_b, b_new)
end


text:
## A SortedSet is a wrapper around balancedTree with
## methods similiar to those of the julia Set.

type SortedSet{K, Ord <: Ordering}
    bt::BalancedTree23{K,Nothing,Ord}

## Zero-argument constructor, or possibly one argument to specify order.

    function SortedSet(o::Ord=Forward)
        bt1 = BalancedTree23{K,Nothing,Ord}(o)
        new(bt1)
    end

end

typealias SetSemiToken IntSemiToken

# The following definition was moved to tokens2.jl
#typealias SetToken Tuple{SortedSet, IntSemiToken}

## This constructor takes an iterable; order defaults
## to Forward

## This one takes an iterable; ordering type is optional.

SortedSet{Ord <: Ordering}(iter, o::Ord=Forward) =
sortedset_with_eltype(iter, eltype(iter), o)

function sortedset_with_eltype{K,Ord}(iter, ::Type{K}, o::Ord)
    h = SortedSet{K,Ord}(o)
    for k in iter
        insert!(h, k)
    end
    h
end

## This function looks up a key in the tree;
## if not found, then it returns a marker for the
## end of the tree.

@inline function find(m::SortedSet, k_)
    ll, exactfound = findkey(m.bt, convert(keytype(m),k_))
    IntSemiToken(exactfound? ll : 2)
end

## This function inserts an item into the tree.
## It returns a bool and a token.
## The bool is true if the inserted item is new.
## It is false if there was already an item
## with that key.
## The token points to the newly inserted item.

@inline function insert!(m::SortedSet, k_)
    b, i = insert!(m.bt, convert(keytype(m),k_), nothing, false)
    b, IntSemiToken(i)
end

## push! is similar to insert but returns the set

@inline function push!(m::SortedSet, k_)
    b, i = insert!(m.bt, convert(keytype(m),k_), nothing, false)
    m
end

## First and last return the first and last (key,data) pairs
## in the SortedDict.  It is an error to invoke them on an
## empty SortedDict.

@inline function first(m::SortedSet)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return m.bt.data[i].k
end

@inline function last(m::SortedSet)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return m.bt.data[i].k
end

@inline function in(k_, m::SortedSet)
    i, exactfound = findkey(m.bt, convert(keytype(m),k_))
    return exactfound
end

@inline eltype{K,Ord <: Ordering}(m::SortedSet{K,Ord}) = K
@inline eltype{K,Ord <: Ordering}(::Type{SortedSet{K,Ord}}) = K
@inline keytype{K,Ord <: Ordering}(m::SortedSet{K,Ord}) = K
@inline orderobject(m::SortedSet) = m.bt.ord

haskey(m::SortedSet, k_) = in(k_, m)

@inline function delete!(m::SortedSet, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    delete!(m.bt, i)
    m
end

@inline function pop!(m::SortedSet, k_)
    k = convert(keytype(m),k_)
    i, exactfound = findkey(m.bt, k)
    !exactfound && throw(KeyError(k_))
    d = m.bt.data[i].d
    delete!(m.bt, i)
    k
end

@inline function pop!(m::SortedSet)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    k = m.bt.data[i].k
    delete!(m.bt, i)
    k
end

## Check if two SortedSets are equal in the sense of containing
## the same K entries.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedSet, m2::SortedSet)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        throw(ArgumentError("Cannot use isequal for two SortedSets unless their element types and ordering objects are equal"))
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1 = deref((m1,p1))
        k2 = deref((m2,p2))
        if !eq(ord,k1,k2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

function union!{K, Ord <: Ordering}(m1::SortedSet{K,Ord}, iterable_item)
    for k in iterable_item
        push!(m1,convert(K,k))
    end
    m1
end

function union(m1::SortedSet, others...)
    mr = packcopy(m1)
    for m2 in others
        union!(mr, m2)
    end
    mr
end

function intersect2{K, Ord <: Ordering}(m1::SortedSet{K, Ord}, m2::SortedSet{K, Ord})
    ord = orderobject(m1)
    mi = SortedSet(K[], ord)
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1) || p2 == pastendsemitoken(m2)
            return mi
        end
        k1 = deref((m1,p1))
        k2 = deref((m2,p2))
        if lt(ord,k1,k2)
            p1 = advance((m1,p1))
        elseif lt(ord,k2,k1)
            p2 = advance((m2,p2))
        else
            push!(mi,k1)
            p1 = advance((m1,p1))
            p2 = advance((m2,p2))
        end
    end
end

function intersect{K, Ord <: Ordering}(m1::SortedSet{K,Ord}, others::SortedSet{K,Ord}...)
    ord = orderobject(m1)
    for s2 in others
        if !isequal(ord, orderobject(s2))
            throw(ArgumentError("Cannot intersect two SortedSets unless their ordering objects are equal"))
        end
    end
    if length(others) == 0
        return m1
    else
        mi = intersect2(m1, others[1])
        for s2 = others[2:end]
            mi = intersect2(mi, s2)
        end
        return mi
    end
end

function symdiff{K, Ord <: Ordering}(m1::SortedSet{K,Ord}, m2::SortedSet{K,Ord})
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2))
        throw(ArgumentError("Cannot apply symdiff to two SortedSets unless their ordering objects are equal"))
    end
    mi = SortedSet(K[], ord)
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        m1end = p1 == pastendsemitoken(m1)
        m2end = p2 == pastendsemitoken(m2)
        if m1end && m2end
            return mi
        elseif m1end
            push!(mi, deref((m2,p2)))
            p2 = advance((m2,p2))
        elseif m2end
            push!(mi, deref((m1,p1)))
            p1 = advance((m1,p1))
        else
            k1 = deref((m1,p1))
            k2 = deref((m2,p2))
            if lt(ord,k1,k2)
                push!(mi, k1)
                p1 = advance((m1,p1))
            elseif lt(ord,k2,k1)
                push!(mi, k2)
                p2 = advance((m2,p2))
            else
                p1 = advance((m1,p1))
                p2 = advance((m2,p2))
            end
        end
    end
end

function setdiff{K, Ord <: Ordering}(m1::SortedSet{K,Ord}, m2::SortedSet{K,Ord})
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2))
        throw(ArgumentError("Cannot apply setdiff to two SortedSets unless their ordering objects are equal"))
    end
    mi = SortedSet(K[], ord)
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return mi
        elseif p2 == pastendsemitoken(m2)
            push!(mi, deref((m1,p1)))
            p1 = advance((m1,p1))
        else
            k1 = deref((m1,p1))
            k2 = deref((m2,p2))
            if lt(ord,k1,k2)
                push!(mi, deref((m1,p1)))
                p1 = advance((m1,p1))
            elseif lt(ord,k2,k1)
                p2 = advance((m2,p2))
            else
                p1 = advance((m1,p1))
                p2 = advance((m2,p2))
            end
        end
    end
end

function setdiff!(m1::SortedSet, iterable)
    for p in iterable
        i = find(m1, p)
        if i != pastendsemitoken(m1)
            delete!((m1,i))
        end
    end
end

function issubset(iterable, m2::SortedSet)
    for k in iterable
        if !in(k, m2)
            return false
        end
    end
    return true
end

function packcopy{K,Ord <: Ordering}(m::SortedSet{K,Ord})
    w = SortedSet(K[], orderobject(m))
    for k in m
        push!(w, k)
    end
    w
end

function packdeepcopy{K, Ord <: Ordering}(m::SortedSet{K,Ord})
    w = SortedSet(K[], orderobject(m))
    for k in m
        newk = deepcopy(k)
        push!(w, newk)
    end
    w
end

function Base.show{K,Ord <: Ordering}(io::IO, m::SortedSet{K,Ord})
    print(io, "SortedSet(")
    keys = K[]
    for k in m
        push!(keys, k)
    end
    print(io, keys)
    println(io, ",")
    print(io, orderobject(m))
    print(io, ")")
end

similar{K,Ord<:Ordering}(m::SortedSet{K,Ord}) =
SortedSet{K,Ord}(orderobject(m))


text:
using Binning, Test

@testset "layer tests" begin
    # Write your tests here.

    # 2x + 3y
    x = range(-10, stop=10, length=100)
    testBin = DomainBinner(x, 20)
    @test length(testBin.b) == 19 # for 20 bins there are 19 interior edges

    # test that function behaves appropriately
    @test length(testBin(rand(100))) == 20

    # test that layer can be called on data batches
    # the output need to have the right shape for subsequent
    # dense layers to be able to matrix multiply.
    # out = testBin(rand(100, 2))
    # @test size(out) == (20, 2)

end


text:
using HDF5: attributes, create_group, h5open
using HDF5: File as h5file
using NeXLCore

"""
   Base.write(h5::HDF5.File, path::AbstractString, mats::AbstractArray{<:Material})

Write an AbstractArray{<:Material} to an HDF5 file.
"""
function Base.write(h5::h5file, path::AbstractString, mats::AbstractArray{<:Material})
    els=sort(collect(mapreduce(m->keys(m), union, mats)))
    data = zeros(Float64, length(els), size(mats)...)
    for (i, el) in enumerate(els), ci in CartesianIndices(mats)
        data[i, ci] = NeXLUncertainties.value(mats[ci][el])
    end
    g = create_group(h5, path)
    attributes(g)["elements"] = map(el->symbol(el), els)
    g["massfractions", compress=3] = data
    mats
end

"""
   Base.read(h5::HDF5.File, path::AbstractString, ::Type{Material})

Read an AbstractArray{<:Material} from an HDF5 file.
"""
function Base.read(h5::h5file, path::AbstractString, ::Type{Material})
    @assert haskey(h5, path)
    g = h5[path]
    els = map(el->parse(Element,el), read(attributes(g)["elements"]))
    tmp = g["massfractions"]
    data = zeros(eltype(tmp), size(tmp)...)
    copy!(data, tmp[axes(tmp)...])
    tmp = nothing
    @assert length(els) == size(data,1) "The number of elements must match the depth of the mass-fractions data."
    return map(CartesianIndices(size(data)[2:end])) do ci
        m = data[:, ci]
        material("M$(ci.I)", Dict(els[i] => m[i] for i in filter(j -> m[j] > 0, eachindex(m))))
    end
end


text:
erf(x::Double64) = Double64Float128(erf, x)
erfc(x::Double64) = Double64Float128(erfc, x)
gamma(x::Double64) = Double64Float128(gamma, x)
lgamma(x::Double64) = Double64Float128(lgamma, x)
loggamma(x::Double64) = Double64Float128(lgamma, x)

erf(x::Double32) = Double32Float128(erf, x)
erfc(x::Double32) = Double32Float128(erfc, x)
gamma(x::Double32) = Double32Float128(gamma, x)
lgamma(x::Double32) = Double32Float128(lgamma, x)
loggamma(x::Double32) = Double32Float128(lgamma, x)

erf(x::Double16) = Double16Float128(erf, x)
erfc(x::Double16) = Double16Float128(erfc, x)
gamma(x::Double16) = Double16Float128(gamma, x)
lgamma(x::Double16) = Double16Float128(lgamma, x)
loggamma(x::Double16) = Double16Float128(lgamma, x)


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, Measures
	using DataFrames, Plots, LaTeXStrings
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.31"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	delta = 0.01
	grid = 0:delta:1
	f(x,y) = 9/8*(4x+y)*sqrt((1-x)*(1-y))
	z = [f(x,y) for y in grid, x in grid]

	densityIntegral = sum(z)*delta^2
	Text("2-dimensional Riemann sum over density: $(densityIntegral)")
end

# ╔═╡ f3fb9a0e-2294-11eb-227d-135a8b801a20
begin
	probB = sum([sum([f(x,y)*delta for y in x:delta:1])*delta for x in grid])
	Text("2-dimensional Riemann sum to evaluate probability: $(probB)")
end

# ╔═╡ f3fbcc54-2294-11eb-1a96-b1676ddb7772
begin
	p1 = surface(grid, grid, z,
		c=cgrad([:blue, :red]), la=1, camera=(60,50),
		ylabel="y", zlabel=L"f(x,y)", legend=:none)
	p2 = contourf(grid, grid, z,
		c=cgrad([:blue, :red]))
	p2 = contour!(grid, grid, z,
		c=:black, xlims=(0,1), ylims=(0,1), ylabel="y", ratio=:equal)

	plot(p1, p2, size=(800, 400), xlabel="x", margin=5mm)
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.31"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═f3fb9a0e-2294-11eb-227d-135a8b801a20
# ╠═f3fbcc54-2294-11eb-1a96-b1676ddb7772
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3