text:
module Mandrill
using Requests
using JSON
using Compat

if VERSION < v"0.4.0"
	Base.split(str, spl; limit=0, keep=true) = split(str, spl, limit, keep)
end

global const api = "https://mandrillapp.com/api/1.0"
# package code goes here
function __init__()
  global mandrill_key = get(ENV, "MANDRILL_KEY", "")
end

function key()
  if !isdefined(Mandrill, :mandrill_key) || isempty(mandrill_key)
      error("""API Key not provided
                   Please define MANDRILL_KEY as an environment variable
                   Or call Mandrill.key(k) before using any other function""")
  else
    return mandrill_key::ASCIIString
  end
end

key(k::ASCIIString) = global mandrill_key = k

macro mandrill_api(endpoint)
  splt = split(string(endpoint), '/'; keep=false)
  @assert size(splt, 1) == 2
  fn = symbol(string(splt[1], "_", splt[2]))
  d = quote
    function $(esc(fn))(req=Dict{ASCIIString, Any}())
      mandrill_api_request($(string(api, endpoint)), req)
    end
  end
end

function mandrill_api_request{T <: AbstractString, S <:Any}(endpoint, req::Dict{T,S})
  req["key"] = key()
  res = Requests.post(URI(endpoint), JSON.json(req))
  if res.status == 200
    return Requests.json(res)
  else
    mandrill_error(Requests.json(res))
  end
end

function mandrill_error(res::Dict)
  error("""Mandrill API $(res["status"])
  $(res["code"]): $(res["name"])
  $(res["message"])""")
end

@mandrill_api "/users/ping2"
@mandrill_api "/users/info"
@mandrill_api "/users/senders"
@mandrill_api "/messages/send"
@mandrill_api "/messages/send-template"
@mandrill_api "/messages/search"
@mandrill_api "/messages/search-time-series"
@mandrill_api "/messages/info"
@mandrill_api "/messages/content"
@mandrill_api "/messages/parse"
@mandrill_api "/messages/send-raw"
@mandrill_api "/messages/list-scheduled"
@mandrill_api "/messages/cancel-scheduled"
@mandrill_api "/messages/reschedule"
@mandrill_api "/tags/list"
@mandrill_api "/tags/delete"
@mandrill_api "/tags/info"
@mandrill_api "/tags/time-series"
@mandrill_api "/tags/all-time-series"
@mandrill_api "/rejects/add"
@mandrill_api "/rejects/list"
@mandrill_api "/rejects/delete"
@mandrill_api "/whitelists/add"
@mandrill_api "/whitelists/list"
@mandrill_api "/whitelists/delete"
@mandrill_api "/senders/list"
@mandrill_api "/senders/domains"
@mandrill_api "/senders/add-domain"
@mandrill_api "/senders/check-domain"
@mandrill_api "/senders/verify-domain"
@mandrill_api "/senders/info"
@mandrill_api "/senders/time-series"
@mandrill_api "/urls/list"
@mandrill_api "/urls/search"
@mandrill_api "/urls/time-series"
@mandrill_api "/urls/tracking-domains"
@mandrill_api "/urls/add-tracking-domain"
@mandrill_api "/urls/check-tracking-domain"
@mandrill_api "/templates/add"
@mandrill_api "/templates/info"
@mandrill_api "/templates/update"
@mandrill_api "/templates/publish"
@mandrill_api "/templates/delete"
@mandrill_api "/templates/list"
@mandrill_api "/templates/time-series"
@mandrill_api "/templates/render"
@mandrill_api "/webhooks/list"
@mandrill_api "/webhooks/add"
@mandrill_api "/webhooks/info"
@mandrill_api "/webhooks/update"
@mandrill_api "/webhooks/delete"
@mandrill_api "/subaccounts/list"
@mandrill_api "/subaccounts/add"
@mandrill_api "/subaccounts/info"
@mandrill_api "/subaccounts/update"
@mandrill_api "/subaccounts/delete"
@mandrill_api "/subaccounts/pause"
@mandrill_api "/subaccounts/resume"
@mandrill_api "/inbound/domains"
@mandrill_api "/inbound/add-domain"
@mandrill_api "/inbound/check-domain"
@mandrill_api "/inbound/delete-domain"
@mandrill_api "/inbound/routes"
@mandrill_api "/inbound/add-route"
@mandrill_api "/inbound/update-route"
@mandrill_api "/inbound/delete-route"
@mandrill_api "/inbound/send-raw"
@mandrill_api "/exports/info"
@mandrill_api "/exports/list"
@mandrill_api "/exports/rejects"
@mandrill_api "/exports/whitelist"
@mandrill_api "/exports/activity"
@mandrill_api "/ips/list"
@mandrill_api "/ips/info"
@mandrill_api "/ips/provision"
@mandrill_api "/ips/start-warmup"
@mandrill_api "/ips/cancel-warmup"
@mandrill_api "/ips/set-pool"
@mandrill_api "/ips/delete"
@mandrill_api "/ips/list-pools"
@mandrill_api "/ips/pool-info"
@mandrill_api "/ips/create-pool"
@mandrill_api "/ips/delete-pool"
@mandrill_api "/ips/check-custom-dns"
@mandrill_api "/ips/set-custom-dns"
@mandrill_api "/metadata/list"
@mandrill_api "/metadata/add"
@mandrill_api "/metadata/update"
@mandrill_api "/metadata/delete"

end # module


text:
module CKANQuery

using HTTP, JSON

export package_list, package_show, resource_list, resource_show
export CKANConnection, CKANDataset, CKANResource

const web_root = "https://open.canada.ca/data/en/"
const api_root = web_root * "api/3/"

include(joinpath(".", "types.jl"))
include(joinpath(".", "show.jl"))
include(joinpath(".", "url.jl"))
include(joinpath(".", "query.jl"))
include(joinpath(".", "actions.jl"))

end

ext: jl


text:
export ActionSetWrapper, MaxAndSkipWrapper

# --------------------------------
# ActionSetWrapper
# --------------------------------
struct ActionSetWrapper <: AbstractGymWrapper
    env::AbstractGymEnv
    action_space::Vector{Int}
end
Reinforce.actions(wrpenv::ActionSetWrapper) = wrpenv.action_space
# TODO maybe put this into reinforce (generalise the GymEnv->AbstractEnvironment)
Reinforce.actions(env::GymEnv) = Reinforce.actions(env, env.state)

"""
MaxAndSkipWrapper(env::AbstractGymEnv, nskip=4)

Take the agent's chosen action over the next `nskip` steps. Return the observation
as the pixelwise max of the last two frames (the `nskip-1`th, and `nskip`th frame).

N.b. using `NoFrameskip` versions of gym environments is encouraged if using this
wrapper, as this wrapper does the skipping instead of the gym.
"""
#=
This wrapper is both an observation wrapper and an action wrapper, but we'll
put it in this file because the observation wrapper file is pretty long already
=#
mutable struct MaxAndSkipWrapper{T, N} <: AbstractGymWrapper
    env::AbstractGymEnv
    nskip::Int
    done::Bool
    recent_frame::Ref{Array{T, N}}
end

function MaxAndSkipWrapper(env::AbstractGymEnv, nskip = 4)
    s = state(env)
    MaxAndSkipWrapper{eltype(s), ndims(s)}(env, nskip, false, s)
end

function Reinforce.step!{T, N}(wrpenv::MaxAndSkipWrapper{T, N}, action)
    """Repeat action, sum reward, and max over last observations."""
    total_reward = 0.0
    done = false
    local s
    for i in 1:wrpenv.nskip
        reward, s = step!(wrpenv.env, action)
        done = Reinforce.finished(wrpenv.env)
        i == wrpenv.nskip - 1 && (wrpenv.recent_frame[] = s)
        total_reward += reward
        done && break
    end
    # Note that the observation on the done=True frame doesn't matter
    wrpenv.recent_frame[] = max.(wrpenv.recent_frame[], s)
    wrpenv.done = done
    return wrpenv.recent_frame[], total_reward
end

Reinforce.state{T, N}(wrpenv::MaxAndSkipWrapper{T, N}) = wrpenv.recent_frame[]

Reinforce.finished(wrpenv::MaxAndSkipWrapper) = wrpenv.done

function Reinforce.reset!(wrpenv::MaxAndSkipWrapper)
    wrpenv.done = false
    Reinforce.reset!(wrpenv.env)
    wrpenv.recent_frame[] = state(wrpenv.env)
end


text:
include(joinpath("src", "LitteringModel.jl"))

@time begin
    include(joinpath("config", "cfg_garbage_collection.jl"))
    LitteringModel.run_batch(
        configs=configs, replicates=5, batchname="garbage_collection", compress=true
    )
end


text:
#lowpass filter a command vector
function apply_lp(mod_samps::AbstractVector, sampr::HasInverseTimeUnits, cutoff::HasInverseTimeUnits)
    if isinf(cutoff)
        return mod_samps
    end
    mod_samps3 = repeat(mod_samps, 11)
    fs = ustrip(uconvert(Unitful.Hz, sampr))
    cutoff = ustrip(uconvert(Unitful.Hz, cutoff))
    responsetype = Lowpass(cutoff; fs=fs)
    designmethod = Butterworth(4)
    filtd = filtfilt(digitalfilter(responsetype, designmethod), ustrip.(mod_samps3)).*unit(mod_samps[1])
    return filtd[5*length(mod_samps)+1:(6*length(mod_samps))]
end

#Use this whenever we need to calculate sample counts so that we are consistent with rounding
function calc_num_samps(duration::HasTimeUnits, sample_rate::HasInverseTimeUnits; should_warn = true)
    unrounded = duration * sample_rate
    rounded = round(Int, unrounded)
    rounding_err = abs(rounded - unrounded)/unrounded
    if rounding_err > 0.01
        @warn "The requested duration cannot be achieved to within 1% accuracy with the current sampling rate.  Consider increasing the sampling rate."
    end
    return rounded
end

#Returns a set of positions _centered_ within the span of pman and pmax and spaced by slice_spacing
#TODO: use this wherever such calculations are done
function slice_positions(pstart::HasLengthUnits, pstop::HasLengthUnits, slice_spacing::HasLengthUnits)
    pstart = uconvert(Unitful.μm, pstart)
    pstop = uconvert(Unitful.μm, pstop)
    slice_spacing = uconvert(Unitful.μm, slice_spacing)
    prng = abs(pstop - pstart)
    nslices = floor(Int, ustrip(prng/slice_spacing)) + 1
    leftover = (prng - ((nslices-1) * slice_spacing))/2
    if pstart > pstop
        leftover = -leftover
    end
    return [range(pstart + leftover, stop=pstop - leftover, length=nslices)...]
end

#This version first adjusts pmin and pmax to respect padding
function slice_positions(pmin::HasLengthUnits, pmax::HasLengthUnits, slice_spacing::HasLengthUnits, slice_pad::HasLengthUnits)
    @assert pmin <= pmax
    pstart = uconvert(Unitful.μm, pmin+slice_pad)
    pstop = uconvert(Unitful.μm, pmax-slice_pad)
    return slice_positions(pstart, pstop, slice_spacing)
end

#The sweep covers an interval that is closed on the left and open on the right.
#i.e. a sweep from 0 to 10 will include a sample at 0 but not at 10
function gen_sweep(pmin::HasLengthUnits, pmax::HasLengthUnits, tsweep::HasTimeUnits, sample_rate::HasInverseTimeUnits)
    tsweeps = uconvert(Unitful.s, tsweep)
    pminum = uconvert(Unitful.μm, pmin)
    pmaxum = uconvert(Unitful.μm, pmax)
    nsamps = calc_num_samps(tsweep, sample_rate)
    increment = (pmaxum - pminum) / nsamps
    return range(pminum, stop=pmaxum-increment, length=nsamps)
end

#moves at 90% of maximum safe speed from pmin to pmax
function gen_safe_sweep(pmin::HasLengthUnits, pmax::HasLengthUnits, sample_rate::HasInverseTimeUnits, rig::AbstractString)
    safe_tsweep = abs(pmax - pmin) / (.90 * PIEZO_MAX_SPEED[rig])
    return gen_sweep(pmin, pmax, safe_tsweep, sample_rate)
end

#Generate sets of samples describing the motion of the positioner during one stack
function gen_sawtooth(pmin::HasLengthUnits, pmax::HasLengthUnits, tfwd::HasTimeUnits, treset::HasTimeUnits, sample_rate::HasInverseTimeUnits)
    fwd_linspace = gen_sweep(pmin, pmax, tfwd, sample_rate)
    reset_linspace = gen_sweep(pmax, pmin, treset, sample_rate)
    return fwd_linspace, reset_linspace
end

gen_bidi_pos(pmin::HasLengthUnits, pmax::HasLengthUnits, tsweep::HasTimeUnits, sample_rate::HasInverseTimeUnits; lp_cutoff = Inf*Hz) = fit_smoothed_bidi(pmin, pmax, upreferred(1/(2*tsweep)), sample_rate; cutoff=lp_cutoff)

#start with triangle and smooth it:
#iteratively adjust limits of a waveform until the filtered version matches the desired pmin and pmax
function fit_smoothed_bidi(pmin, pmax, f, sr; cutoff = 150.0Hz)
	if cutoff == Inf*Hz
		fwd, bck =  ImagineInterface.gen_sawtooth(pmin, pmax, upreferred(1/(2*f)), upreferred(1/(2*f)), sr)
		return fwd,bck
	end
	thresh = 0.001μm
	v = 100.0μm
	wvform = fwd = bck = -1
	dxtra = 0.0μm
	while v > thresh
		fwd, bck =  ImagineInterface.gen_sawtooth(pmin-dxtra, pmax+dxtra, upreferred(1/(2*f)), upreferred(1/(2*f)), sr)
		wvform = apply_lp(vcat(fwd,bck), sr, cutoff)
		mx = maximum(wvform)
		mn = minimum(wvform)
		v = max(abs(pmin - mn), abs(pmax - mx))
		dxtra += v/2
	end
	if isnan(dxtra)
		error("Smoothing seems to have failed")
	end
    return wvform[1:div(length(wvform),2)], wvform[div(length(wvform),2)+1:end]
end

#returns a vector of sample-index intervals with starts separated by `interval_spacing`.
#The first interval is offset from the first sample of the sampled region by the `offset` keyword arg
#The `z_pad` kwarg, specified in non-temporal units, will prevent placement of intervals at the extremes of the sample space
#The `alignment` kwarg determines whether the first interval begins at the first valid sample (:start) or the last interval ends at the last valid sample (:stop)
#Thus :start and :stop will produce equal results for certain well-dividing sample counts, but most of the time they are different
function spaced_intervals(samples_space::StepRangeLen{TS}, interval_spacing::TS, interval_duration::TT, sample_rate::TTI;
                                delay::TT=uconvert(unit(TT), 0.0*Unitful.s), z_pad::TS = uconvert(unit(TS), 1.0*Unitful.μm),
                                alignment=:start, rig="ocpi-2") where{TS<:HasLengthUnits, TT<:HasTimeUnits, TTI<:HasInverseTimeUnits}
    if !in(alignment, (:start, :stop))
        error("Only :start and :stop alignment is supported")
    end
    samp_length = abs(step(samples_space))
    pad_samps = round(Int, z_pad/samp_length)
    delay_samps = calc_num_samps(delay, sample_rate)
    offset = pad_samps + delay_samps
    nsamps = length(samples_space) - offset - pad_samps
    dur_samps = calc_num_samps(interval_duration, sample_rate)
    spacing_samps = interval_spacing/samp_length
    if abs(spacing_samps-round(Int, spacing_samps))/spacing_samps > 0.01
        @warn "The requested spacing cannot be achieved to within 1% accuracy with the current sampling rate.  Consider increasing the sampling rate."
    end
    spacing_samps = round(Int, spacing_samps)
    if spacing_samps <= dur_samps
        error("The requested stack spacing results in overlapping camera exposure and/or laser pulse intervals.  Increase z-slice spacing, decrease exposure time, or change sampling rate.")
    end
    inter_samps = spacing_samps - dur_samps #number of samples between end of one interval and start of the next
    if inter_samps < ceil(Int, CAMERA_OFF_TIME[rig] * sample_rate)
        error("The requested spacing results in intervals which are too close in time for the jitter specification of the camera.  Increase interval_spacing, decrease interval_duration, or change sampling rate")
    end
    cycle_samps = inter_samps + dur_samps #number of samples in one whole cycle
    nintervals = div(nsamps, cycle_samps)
    extra = mod(nsamps-dur_samps, cycle_samps) #first cycle is partial
    if extra >= dur_samps
        nintervals+=1
        extra-=dur_samps
    end
    if nintervals == 0
        error("The requested interval_duration is longer than the sampled region.")
    end
    output = Array{ClosedInterval{Int}}(undef, nintervals)
    if alignment == :stop
        offset+=extra
    end
    output[1] = ClosedInterval(offset+1, offset+dur_samps)
    for i = 2:nintervals
        curstart = offset + (i-1)*cycle_samps+1
        output[i] = ClosedInterval(curstart, curstart+dur_samps-1) #offset+(i-1)*cycle_samps+1, offset+i*cycle_samps)
    end
    return output
end

function gen_pulses!(output::Vector{Bool}, pulsevec::Vector{ClosedInterval{Int}}; fillval=true)
    for p in pulsevec
        output[minimum(p):maximum(p)] .= fillval
    end
    return output
end
gen_pulses(nsamps::Int, pulsevec::Vector{ClosedInterval{Int}}; fillval=true) = gen_pulses!(fill(!fillval, nsamps), pulsevec; fillval=fillval)

function center(input::ClosedInterval{Int})
    w = width(input)
    halfw = div(w,2)
    return minimum(input) + halfw
end

#Return a new "scaled" interval centered on the center of the input and with width equal to frac * width(input)
function scale(input::ClosedInterval{Int}, frac::Float64)
    halfw = div(width(input),2)
    ctr = center(input)
    halfw_new = round(Int, frac*halfw)
    return ClosedInterval(ctr - halfw_new, ctr + halfw_new)
end

function gen_bidirectional_stack(pmin::TL, pmax::TL, z_spacing::TL,
                                 stack_time::TT, exp_time::TT, sample_rate::TTI, flash_frac::Real;
                                 z_pad::TL = 1.0*Unitful.μm, alternate_cameras = false,
                                 rig="ocpi-2") where{TL<:HasLengthUnits, TT<:HasTimeUnits, TTI<:HasInverseTimeUnits}
    if pmin == pmax
        error("Use the gen_2d_timeseries function instead of setting pmin and pmax to the same value")
    end
    flash = true
    if flash_frac >= 1.0
        @warn "flash_frac was set greater than 1.0, so keeping laser on throughout the stack"
        flash = false
    elseif flash_frac <= 0
        error("flash_frac must be positive")
    end

    pmin = uconvert(Unitful.μm, pmin)
    pmax = uconvert(Unitful.μm, pmax)
    z_spacing = uconvert(Unitful.μm, z_spacing)
    z_pad = uconvert(Unitful.μm, z_pad)

    stack_time = uconvert(Unitful.s, stack_time)
    exp_time = uconvert(Unitful.s, exp_time)
    sample_rate = uconvert(inv(Unitful.s), sample_rate)

    nsamps_stack = calc_num_samps(stack_time, sample_rate)
    posfwd, posback = gen_bidi_pos(pmin, pmax, stack_time, sample_rate)
    #offset by one sample going forward so that we don't use the end points of the triangle
    delay1samp = 1/sample_rate

    #We want the cameras to reach the global exposure state at the same time for each pulse in forward and reverse stacks, so we
    #don't line up the pulses exactly.  Assuming that the exposure is one las_frac's duration longer than the minimum exposure time for the given
    #should really put some more work into making this modular and adaptable to different exposure times / ROI sizes / line times /cameras
    #For now we deliver the flash at the end of the exposure pulse, and we align the forward and reverse stacks to the flash (not to the exposure)
    #Therefor it's up to the user to make sure the exposure time and ROI size combination is long enough that the flash fits within the global
    #exposure period (see PCO.Edge manual for more details)

    #distance traveled in one exposure time.  We need to adjust padding by that much to make sure that exposures can fit.
    exp_dist = exp_time * abs(pmax - pmin) / stack_time
    if exp_dist > z_pad
        @warn "Increasing z padding by $(exp_dist-z_pad) to leave space for bidirectional exposures\n"
    end
    exp_intervals_fwd = spaced_intervals(posfwd, z_spacing, exp_time, sample_rate; delay=delay1samp, z_pad = max(z_pad,exp_dist), alignment=:start, rig=rig)
    flash_nsamps = calc_num_samps(min(flash_frac,1.0) * exp_time, sample_rate)
    exp_nsamps = width(exp_intervals_fwd[1]) + 1
    offset_nsamps = exp_nsamps - flash_nsamps
    #The extra -1 is needed since the forward direction does not sample the last point (the first index in the reverse direction)
    exp_intervals_back = map(x-> ClosedInterval(length(posfwd)-(maximum(x)+offset_nsamps-2), length(posfwd)-(maximum(x)-flash_nsamps-1)), exp_intervals_fwd)
#    exp_intervals_back = spaced_intervals(posback, z_spacing, exp_time, sample_rate; delay=0.0*Unitful.s, z_pad = z_pad, alignment=:stop, rig=rig)

    if flash
        las_intervals_fwd = map(x-> ClosedInterval(maximum(x)-flash_nsamps+1, maximum(x)), exp_intervals_fwd)
        lasfwd = gen_pulses(nsamps_stack, las_intervals_fwd)
        las_intervals_back = map(x-> ClosedInterval(maximum(x)-flash_nsamps+1, maximum(x)), exp_intervals_back)
        lasback = gen_pulses(nsamps_stack, las_intervals_back)
    else
        lasfwd = fill(true, length(posfwd))
        lasback = fill(true, length(posback))
    end

    camfwd = gen_pulses(nsamps_stack, exp_intervals_fwd)
    camback = gen_pulses(nsamps_stack, exp_intervals_back)
    if alternate_cameras
        output = Dict("positioner" => vcat(posfwd, posback),
                    "camera_fwd" => vcat(camfwd, fill(false, length(camback))),
                    "camera_back" => vcat(fill(false, length(camfwd)), camback),
                    "laser_fwd" => vcat(lasfwd, fill(false, length(lasback))),
                    "laser_back" => vcat(fill(false, length(lasfwd)), lasback),
                    "nframes_per_cam" => length(exp_intervals_fwd))
        return output
    else
        output = Dict("positioner" => vcat(posfwd, posback), "camera" => vcat(camfwd, camback), "laser" => vcat(lasfwd, lasback), "nframes" => length(exp_intervals_fwd)*2)
        return output
    end
end

function gen_stepped_stack(pmin::TL, pmax::TL, z_spacing::TL, pause_time::TT,
                           reset_time::TT, exp_time::TT, sample_rate::TTI, flash_frac::Real;
                           rig="ocpi-2") where{TL<:HasLengthUnits, TT<:HasTimeUnits, TTI<:HasInverseTimeUnits}
    if exp_time > pause_time
        error("Pause duration must be greater than or equal to exposure duration")
    end
    if pmin == pmax
        error("Use the gen_2d_timeseries function instead of setting pmin and pmax to the same value")
    end
    flash = true
    if flash_frac > 1
        @warn "las_frac was set greater than 1, so defaulting to keeping laser on throughout the stack"
        flash = false
    elseif flash_frac <= 0
        error("las_frac must be positive")
    end

    pmin = uconvert(Unitful.μm, pmin)
    pmax = uconvert(Unitful.μm, pmax)
    z_spacing = uconvert(Unitful.μm, z_spacing)

    pause_time = uconvert(Unitful.s, pause_time)
    exp_time = uconvert(Unitful.s, exp_time)
    sample_rate = uconvert(inv(Unitful.s), sample_rate)

    curmin = curmax = pmin
    prng = pmax - pmin
    nslices = floor(Int, prng / z_spacing)+1
    pause_starts = Int[]
    pause_stops = Int[]
    posfwd = typeof(curmin)[]
    for i = 1:(nslices-1)
        push!(pause_starts, length(posfwd)+1)
        cur_pause = gen_sweep(curmax, curmax, pause_time, sample_rate)
        append!(posfwd, cur_pause)
        curmin = curmax
        curmax = curmax + z_spacing
        push!(pause_stops, length(posfwd))
        cur_sweep = gen_safe_sweep(curmin, curmax, sample_rate,rig)
        append!(posfwd, cur_sweep)
    end
    #final pause
    push!(pause_starts, length(posfwd)+1)
    append!(posfwd, gen_sweep(curmax, curmax, pause_time, sample_rate))
    nsamps_stack = length(posfwd)
    push!(pause_stops, nsamps_stack)
    #camera pulse intervals
    exp_intervals = Array{ClosedInterval{Int}}(undef, nslices)
    pause_nsamps = pause_stops[1] - pause_starts[1]
    exp_nsamps = calc_num_samps(exp_time, sample_rate)
    samp_offset = floor(Int, (pause_nsamps - exp_nsamps) / 2) #center the exposure within the positioner pause
    for i = 1:nslices
        exp_start = pause_starts[i] + samp_offset
        exp_intervals[i] = ClosedInterval(exp_start, exp_start+exp_nsamps-1)
    end
    camfwd = gen_pulses(nsamps_stack, exp_intervals)
    #laser pulse intervals
    if flash
        pulse_nsamps = calc_num_samps(flash_frac * exp_time, sample_rate)
        las_intervals = map(x-> ClosedInterval(maximum(x)-pulse_nsamps+1, maximum(x)), exp_intervals)
        lasfwd = gen_pulses(nsamps_stack, las_intervals)
    else
        lasfwd = fill(true, nsamps_stack)
    end
    #reset
    posreset = gen_sweep(curmax, pmin, reset_time, sample_rate)
    reset_digi = fill(false, length(posreset))
    output = Dict("positioner" => vcat(posfwd, posreset), "camera" => vcat(camfwd, reset_digi), "laser" => vcat(lasfwd, reset_digi), "nframes" => length(exp_intervals))
    return output
end

function gen_unidirectional_stack(pmin::TL, pmax::TL, z_spacing::TL, stack_time::TT, reset_time::TT,
                                  exp_time::TT, sample_rate::TTI, flash_frac::Real;
                                  z_pad::TL = 1.0*Unitful.μm,
                                  rig="ocpi-2") where{TL<:HasLengthUnits, TT<:HasTimeUnits, TTI<:HasInverseTimeUnits}
    if pmin == pmax
        error("Use the gen_2d_timeseries function instead of setting pmin and pmax to the same value")
    end
    flash = true
    if flash_frac > 1
        @warn "las_frac was set greater than 1, so defaulting to keeping laser on throughout the stack"
        flash = false
    elseif flash_frac <= 0
        error("las_frac must be positive")
    end

    pmin = uconvert(Unitful.μm, pmin)
    pmax = uconvert(Unitful.μm, pmax)
    z_spacing = uconvert(Unitful.μm, z_spacing)
    z_pad = uconvert(Unitful.μm, z_pad)

    stack_time = uconvert(Unitful.s, stack_time)
    exp_time = uconvert(Unitful.s, exp_time)
    sample_rate = uconvert(inv(Unitful.s), sample_rate)

    nsamps_stack = calc_num_samps(stack_time, sample_rate)
    posfwd, posreset = gen_sawtooth(pmin, pmax, stack_time, reset_time, sample_rate)

    exp_intervals = spaced_intervals(posfwd, z_spacing, exp_time, sample_rate; z_pad = z_pad, alignment=:start, rig=rig)

    if flash
        pulse_nsamps = calc_num_samps(flash_frac * exp_time, sample_rate)
        las_intervals = map(x-> ClosedInterval(maximum(x)-pulse_nsamps+1, maximum(x)), exp_intervals)
        lasfwd = gen_pulses(nsamps_stack, las_intervals)
    else
        lasfwd = fill(true, length(posfwd))
    end

    camfwd = gen_pulses(nsamps_stack, exp_intervals)
    reset_digi = fill(false, length(posreset))

    output = Dict("positioner" => vcat(posfwd, posreset), "camera" => vcat(camfwd, reset_digi), "laser" => vcat(lasfwd, reset_digi), "nframes" => length(exp_intervals))

    return output
end

function gen_2d_timeseries(position::TL, nframes::Int, exp_time::TT, inter_exp_time::TT,
                           sample_rate::TTI, flash_frac::Real) where{TL<:HasLengthUnits, TT<:HasTimeUnits, TTI<:HasInverseTimeUnits}
    cycle_time = exp_time + inter_exp_time
    exp_samps = round(Int, exp_time * sample_rate)
    cycle_samps = round(Int, cycle_time * sample_rate)
    inter_exp_samps = round(Int, inter_exp_time * sample_rate)
    nsamps = round(Int, nframes * cycle_time * sample_rate)
    pos = fill(position, nsamps)
    exp_intervals = Array{ClosedInterval{Int}}(undef, nframes)
    curstart = div(inter_exp_samps,2) #offset by half the interframe time
    for i = 1:nframes
        exp_intervals[i] = ClosedInterval(curstart, curstart+exp_samps-1)
        curstart += cycle_samps
    end
    pulse_nsamps = calc_num_samps(flash_frac * exp_time, sample_rate)
    las_intervals = map(x-> ClosedInterval(maximum(x)-pulse_nsamps+1, maximum(x)), exp_intervals)
    cam = gen_pulses(nsamps, exp_intervals)
    las = gen_pulses(nsamps, las_intervals)
    output = Dict("positioner" => pos, "camera" => cam, "laser" => las)
    return output
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
# Boyd & Vandenberghe, "Convex Optimization"
# Joëlle Skaf - 08/16/05
#
# Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14
#
# The goal is to find the largest Euclidean ball (i.e. its center and
# radius) that lies in a polyhedron described by affine inequalites in this
# fashion: P = {x : a_i'*x <= b_i, i=1,...,m} where x is in R^2

using Convex

# Generate the input data
a1 = [ 2;  1];
a2 = [ 2; -1];
a3 = [-1;  2];
a4 = [-1; -2];
b = ones(4, 1);

# Create and solve the model
r = Variable(1)
x_c = Variable(2)
p = maximize(r)
p.constraints += a1' * x_c + r * Base.norm(a1, 2) <= b[1];
p.constraints += a2' * x_c + r * Base.norm(a2, 2) <= b[2];
p.constraints += a3' * x_c + r * Base.norm(a3, 2) <= b[3];
p.constraints += a4' * x_c + r * Base.norm(a4, 2) <= b[4];
solve!(p)
p.optval

# Generate the figure
x = linspace(-1.5, 1.5);
theta = 0:pi/100:2*pi;
using Gaston
Gaston.set_terminal("x11")
plot(x, -x * a1[1] / a1[2] + b[1] / a1[2], "color", "black",
     x, -x * a2[1]/ a2[2] + b[2] / a2[2], "color", "red",
     x, -x * a3[1]/ a3[2] + b[3] / a3[2], "color", "green",
     x, -x * a4[1]/ a4[2] + b[4] / a4[2], "color", "blue",
     x_c.value[1] + r.value * cos(theta), x_c.value[2] + r.value * sin(theta), "color", "red",
     "title", "Largest Euclidean ball lying in a 2D polyhedron",
    );


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:PrioritizedDQN},
    ::Val{:CartPole},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_PrioritizedDQN_CartPole_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = StableRNG(seed)

    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    agent = Agent(
        policy = QBasedPolicy(
            learner = PrioritizedDQNLearner(
                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(),
                ),
                target_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(),
                ),
                loss_func = (ŷ, y) -> huber_loss(ŷ, y; agg = identity),
                stack_size = nothing,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 100,
                update_freq = 1,
                target_update_freq = 100,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArrayPSARTTrajectory(
            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
            if agent.policy.learner.update_step % agent.policy.learner.update_freq == 0
                with_logger(lg) do
                    @info "training" loss = agent.policy.learner.loss
                end
            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 the `PrioritizedDQNLearner` method with three dense layers to approximate the Q value.
    The testing environment is CartPoleEnv.
    """

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


text:
using NESSie.Format
using NESSie.Format: readhmo_nodes, readhmo_elements, readhmo_charges

testfiles = ((mktemp()..., (0,0,0)), # empty file
             (mktemp()..., (3,2,2))) # dummy file
write(testfiles[2][2], """
# stuffing
stuffing

BEG_NODL_DATA
	stuffing
    x 0.0 0.0 0.0
	x 1.0 2.0 3.0
	x -2.00001 1.337 42.0
END_NODL_DATA
stuffing

BEG_ELEM_DATA
	stuffing
	x x x 1 2 3
	x x x 2 3 1
END_ELEM_DATA
stuffing

BEG_CHARGE_DATA
	stuffing
	x 1.0 2.0 3.0 999.0
	x -2.00001 1.337 42.0 -999.9
END_CHARGE_DATA

""")

@testset "readhmo" begin
    try
        for (fname, fh, len) in testfiles, fn in ("single", "bulk"), T in testtypes
            seekstart(fh)
            nodes = (); elements = (); charges = ()
            if fn == "single"
                # check single methods (atstart=false)
                nodes = readhmo_nodes(fh, T)
                elements = readhmo_elements(fh, nodes)
                charges = readhmo_charges(fh, T)
            elseif fn == "bulk"
                # check bulk method (implicit single methods with atstart=true)
                model = readhmo(fh, T)
                nodes, elements, charges = (model.nodes, model.elements, model.charges)
            end
            # check return types
            @test typeof(nodes) == Vector{Vector{T}}
            @test typeof(elements) == Vector{Triangle{T}}
            @test typeof(charges) == Vector{Charge{T}}
            # check lengths
            @test length(nodes) == len[1]
            @test length(elements) == len[2]
            @test length(charges) == len[3]

            # check content
            if fname == testfiles[2][1]
                @test nodes[1] == T[0, 0, 0]
                @test nodes[2] == T[1, 2, 3]
                @test nodes[3] == T[-2.00001, 1.337, 42]
                @test elements[1].v1 === nodes[1]
                @test elements[1].v2 === nodes[2]
                @test elements[1].v3 === nodes[3]
                @test elements[2].v1 === nodes[2]
                @test elements[2].v2 === nodes[3]
                @test elements[2].v3 === nodes[1]
                @test charges[1].pos == T[1, 2, 3]
                @test charges[1].val == T(999)
                @test charges[2].pos == T[-2.00001, 1.337, 42]
                @test charges[2].val == T(-999.9)
            end
        end
    finally
        [begin close(fh); rm(fn); end for (fn, fh, _) in testfiles]
    end
end

@testset "writehmo" begin
    @test_skip writehmo
end


text:
mutable struct TimeStepWizard{FT, C, D}
                         cfl :: FT
               diffusive_cfl :: FT
                  max_change :: FT
                  min_change :: FT
                      max_Δt :: FT
                      min_Δt :: FT
    cell_advection_timescale :: C
    cell_diffusion_timescale :: D
end

infinite_diffusion_timescale(args...) = Inf # its not very limiting

"""
    TimeStepWizard(cfl=0.1, max_change=2.0, min_change=0.5, max_Δt=Inf, min_Δt=0.0)

Callback for adapting simulation time-steps `Δt` to maintain the advective
Courant-Freidrichs-Lewy (`cfl`) number, the `diffusive_cfl`, while maintaining
`max_Δt`, `min_Δt`, and satisfying `max_change` and `min_change` criteria
so `Δt` is not adapted "too quickly".

For more information on `cfl`, see
https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition.

Example
=======

To use `TimeStepWizard`, adapt in a [`Callback`](@ref) and add it to a `Simulation`:

```julia
julia> simulation = Simulation(model, Δt=0.9, stop_iteration=100)

julia> wizard = TimeStepWizard(cfl=0.2)

julia> simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))
```

Then when `run!(simulation)` is invoked, the time-step `simulation.Δt` will be updated every 4 iterations.
Note that the name `:wizard` is unimportant.
"""
function TimeStepWizard(FT=Float64; cfl = 0.1,
                                    diffusive_cfl = Inf,
                                    max_change = 2.0,
                                    min_change = 0.5,
                                    max_Δt = Inf,
                                    min_Δt = 0.0,
                                    cell_advection_timescale = cell_advection_timescale,
                                    cell_diffusion_timescale = infinite_diffusion_timescale)

    isfinite(diffusive_cfl) && # user wants to limit by diffusive CFL
    !(cell_diffusion_timescale === infinite_diffusion_timescale) && # user did not provide custom timescale
        (cell_diffusion_timescale = Oceananigans.TurbulenceClosures.cell_diffusion_timescale)

    C = typeof(cell_advection_timescale)
    D = typeof(cell_diffusion_timescale)

    return TimeStepWizard{FT, C, D}(cfl, diffusive_cfl, max_change, min_change, max_Δt, min_Δt,
                                    cell_advection_timescale, cell_diffusion_timescale)
end

using Oceananigans.Grids: topology

"""
     new_time_step(old_Δt, wizard, model)

Return a new time_step given `model.velocities` and model diffusivites,
and the parameters of the `TimeStepWizard` `wizard`.
"""
function new_time_step(old_Δt, wizard, model)

    advective_Δt = wizard.cfl * wizard.cell_advection_timescale(model)
    diffusive_Δt = wizard.diffusive_cfl * wizard.cell_diffusion_timescale(model)

    new_Δt = min(advective_Δt, diffusive_Δt)

    # Put the kibosh on if needed
    new_Δt = min(wizard.max_change * old_Δt, new_Δt)
    new_Δt = max(wizard.min_change * old_Δt, new_Δt)
    new_Δt = clamp(new_Δt, wizard.min_Δt, wizard.max_Δt)

    return new_Δt
end

(wizard::TimeStepWizard)(simulation) =
    simulation.Δt = new_time_step(simulation.Δt, wizard, simulation.model)


text:
#=
   Can one not have two models in the same file?
   Yes, it don't seems to be any problem...

=#
using Turing

# Define a simple Normal model with unknown mean and variance.
a = []
@model test_model(n) = begin
        d = TArray{Int}(undef, n)
        # d = Vector{Int}(undef, n)
        for i in 1:n
            d[i] ~ DiscreteUniform(0,9)
        end
        s = sum(d)
        push!(a, s)
end

#  Run sampler, collect results
chn = sample(test_model(10), MH(), 1000)
println("mean(a): $(mean(a))")

a = []
@model test_model2(n) = begin
        d = TArray{Int}(undef, n)
        # d = Vector{Int}(undef, n)
        for i in 1:n
            d[i] ~ DiscreteUniform(0,9)
        end
        s = sum(d)
        push!(a, s)
end

#  Run sampler, collect results
chn = sample(test_model2(10), MH(), 1000)
println("mean(a) 2: $(mean(a))")

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Test creation and evaluation of SPNs

@testset "Loading SPNs from file" begin

    @testset "Simple Categorical SPN" begin
        io = IOBuffer("""#
1 + 2 0.2 3 0.5 4 0.3
2 * 5 7
3 * 5 8
4 * 6 8
# Leaves
5 categorical 1 0.6 0.4
6 categorical 1 0.1 0.9
7 categorical 2 0.3 0.7
8 categorical 2 0.8 0.2
#""")
        S = SumProductNetwork(io)

        @testset "Saving to file then reading it" begin
            SPN = SumProductNetwork([
                SumNode([2,3,4],[0.2,0.5,0.3]),        # 1
                ProductNode([5,7]),                    # 2
                ProductNode([5,8]),                    # 3
                ProductNode([6,8]),                    # 4
                CategoricalDistribution(1, [0.6,0.4]), # 5
                CategoricalDistribution(1, [0.1,0.9]), # 6
                CategoricalDistribution(2, [0.3,0.7]), # 7
                CategoricalDistribution(2, [0.8,0.2])  # 8
            ])
            # @info "## Saving network to disk"
            SPNetworks.save(SPN,normpath("$(@__DIR__)/../assets/test.spn"))
            # @info "## Loading network from file"
            SPN2 = SumProductNetwork(normpath("$(@__DIR__)/../assets/test.spn"))
            # @info "## Comparing sizes"
            @test length(SPN) == length(SPN2) == length(S)
            @test nparams(SPN) == nparams(SPN2) == nparams(S)
            @testset "Evaluating simple SPNs at $a,$b" for a in 1:2, b in 1:2
                pred = SPN2([a,b])
                ref = SPN([a,b])
                ref2 = S([a,b])
                @test pred == ref == ref2
            end
        end
    end

    # @info "## Printing graph structure"
    # io = IOBuffer()
    # todot(io,SPN)
    # open("$(@__DIR__)/assets/spn.dot", "w") do f
    #     write(f, String(take!(io)))
    # end
    # println("Call dot -Tpdf assets/spn.dot -o assets/spn.pdf to generate image")

    @testset "XOR SPN" begin
        XOR = SumProductNetwork(normpath("$(@__DIR__)/../assets/xor.spn"))
        @testset "Evaluation XOR at $a,$b,$c" for a=1:2, b=1:2, c=1:2
            v = XOR([a,b,c])
            @test ((a + b + c - 3) % 2 == 1 && v ≈ 0.0) || ((a + b + c - 3) % 2 == 0 && v ≈ 0.25)
        end
    @test logpdf(XOR,fill(NaN,3)) ≈ 0.0
    end

    @testset "HMM SPN"  begin
        HMM = SumProductNetwork(normpath("$(@__DIR__)/../assets/hmm.spn"))
        results = [0.11989139999999997,0.06615860000000003,0.29298060000000004,
                0.1709694,0.0708666,0.03658340000000001,0.1561014,0.08644860000000001]
        @testset "Evaluating HMM at $a,$b,$c" for a=1:2, b=1:2, c=1:2
            v = HMM([a,b,c])
            @test v ≈ results[4*(a-1) + 2*(b-1) + c]
        end
        @test logpdf(HMM,fill(NaN,3)) ≈ 0.0
    end

    @testset "Reading Small SPN in pyspn's format" begin
        PySPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/example.pyspn.spn"); offset = 1)
        @test length(PySPN) == 13
        @test length(scope(PySPN)) == 2
        @test nparams(PySPN) == 12
        @test logpdf(PySPN,[NaN,NaN]) ≈ 0.0 atol=1e-6
        @test PySPN(2,1) ≈ 0.4
        @test logpdf(PySPN, [2,NaN]) ≈ log(0.55)
    end

    @testset "Reading Breast-Cancer SPN in pyspn's format" begin
        PySPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/breast-cancer.spn"); offset = 1)
        @test length(scope(PySPN)) == 10
        @test logpdf(PySPN,fill(NaN,10)) ≈ 0.0 atol=1e-6
    end
end


text:
using Combinatorics
using SparseArrays
using ArnoldiMethod

function CASCI{T}(Alg::SparseHamiltonian) where T <: AbstractFloat

    @output "Getting molecule...\n"
    molecule = Molecule()
    @output "Computing AO Integrals...\n"
    #aoint = ConventionalAOIntegrals()

    @output "Calling RHF module...\n"
    refwfn = Fermi.HartreeFock.RHF(molecule)
    ints = refwfn.ints

    @output "Transforming Integrals for CAS computation...\n"
    # Read options
    frozen = Fermi.CurrentOptions["cas_frozen"]

    nmo = refwfn.ndocc + refwfn.nvir

    act_elec = 2*(refwfn.ndocc - frozen)

    if act_elec < 0
        error("\nInvalid number of frozen orbitals ($frozen) for $(2*refwfn.ndocc) electrons.")
    end

    # Active = -1 means FCI, with frozen
    if Fermi.CurrentOptions["cas_active"] == -1
        active = nmo - frozen
    else
        active = Fermi.CurrentOptions["cas_active"]
    end

    if active ≤ act_elec/2
        error("\nNumber of active orbitals ($active) too small for $(act_elec) active electrons")
    end

    if active+frozen > nmo
        error("\nNumber of active ($active) and frozen orbitals ($frozen) greater than number of orbitals ($nmo)")
    end

    s = 1:(frozen+active)
    h = T.(Fermi.Integrals.transform_fock(ints["T"] + ints["V"], ints.orbs["FU"][s], ints.orbs["FU"][s]))
    V = T.(Fermi.Integrals.transform_eri(ints["μ"], ints.orbs["FU"][s], ints.orbs["FU"][s], ints.orbs["FU"][s], ints.orbs["FU"][s]))

    aoint = nothing
    CASCI{T}(refwfn, h, V, frozen, act_elec, active, Alg)
end

function CASCI{T}(refwfn::Fermi.HartreeFock.RHF, h::Array{T,2}, V::Array{T,4}, frozen::Int, act_elec::Int, active::Int, Alg::SparseHamiltonian) where T <: AbstractFloat

    # Print intro
    Fermi.ConfigurationInteraction.print_header()
    @output "\n    • Computing FCI with the SparseMatrix algorithm.\n\n"

    nroot = Fermi.CurrentOptions["cas_nroot"]

    @output "\n →  ACTIVE SPACE\n"
    @output "Frozen Orbitals:  {:3d}\n" frozen
    @output "Active Electrons: {:3d}\n" act_elec
    @output "Active Orbitals:  {:3d}\n" active

    @output "\nGenerating determinants ..."
    t = @elapsed begin
        dets = get_determinants(act_elec, active, frozen)
        Ndets = length(dets)
    end

    @output " done in {} seconds.\n" t
    @output "\nNumber of Determinants: {:10d}\n" Ndets

    @output "\nBuilding Sparse Hamiltonian..."

    t = @elapsed H = get_sparse_hamiltonian_matrix(dets, h, V, Fermi.CurrentOptions["cas_cutoff"])
    @output " done in {:5.5f} seconds.\n" t
    @output "Hamiltonian Matrix size: {:10.3f} Mb\n" Base.summarysize(H)/10^6

    @output "Diagonalizing Hamiltonian for {:3d} eigenvalues..." nroot
    t = @elapsed begin
        decomp, history = partialschur(H, nev=nroot, tol=10^-12, which=LM())
        λ, ϕ = partialeigen(decomp)
    end
    @output " done in {:5.5f} seconds.\n" t
    @output "\n Final FCI Energy: {:15.10f}\n" λ[1]+refwfn.molecule.Vnuc

    # Sort dets by importance
    C = ϕ[:,1]
    sp = sortperm(abs.(ϕ[:,1]), rev=true)
    C = ϕ[:,1][sp]
    dets = dets[sp]

    @output "\n • Most important determinants:\n\n"

    for i in 1:50
        @output "{:15.5f}      {}\n" C[i]  detstring(dets[i], frozen+active)
    end
    @output "\n"
    return CASCI{T}(refwfn, λ[1]+T(refwfn.molecule.Vnuc), dets, C)
end

function get_determinants(Ne::Int, No::Int, nfrozen::Int)

    Nae = Int(Ne/2)
    occ_string = repeat('1', nfrozen)

    zeroth = repeat('1', Nae)*repeat('0', No-Nae)

    perms = multiset_permutations(zeroth, length(zeroth))

    dets = [Determinant[] for i = 1:Threads.nthreads()]
    @sync for αstring in perms
        Threads.@spawn begin
        for βstring in perms
            α = occ_string*join(αstring)
            β = occ_string*join(βstring)
            _det = Determinant(α, β)
            push!(dets[Threads.threadid()], _det)
        end
        end #Threads.@spawn
    end

    dets = vcat(dets...)
    # Determinant list is sorted by its excitation level w.r.t the first determinant (normally HF)
    #sort!(dets, by=d->excitation_level(dets[1], d))

    return dets
end

function get_sparse_hamiltonian_matrix(dets::Array{Determinant,1}, h::Array{T,2}, V::Array{T,4}, tol::Float64) where T <: AbstractFloat

    Ndets = length(dets)
    Nα = sum(αlist(dets[1]))
    Nβ = sum(αlist(dets[1]))

    αind = [Array{Int64,1}(undef,Nα) for i = 1:Threads.nthreads()]
    βind = [Array{Int64,1}(undef,Nβ) for i = 1:Threads.nthreads()]
    vals = [T[] for i = 1:Threads.nthreads()]
    ivals = [Int64[] for i = 1:Threads.nthreads()]
    jvals = [Int64[] for i = 1:Threads.nthreads()]

    @sync for i in 1:Ndets
        Threads.@spawn begin
        D1 = dets[i]
        αindex!(D1, αind[Threads.threadid()])
        βindex!(D1, βind[Threads.threadid()])
        elem = 0.0
        for j in i:Ndets
            D2 = dets[j]
            αexc = αexcitation_level(D1,D2)
            βexc = βexcitation_level(D1,D2)
            el = αexc + βexc
            if el > 2
                continue
            elseif el == 2
                elem = Hd2(D1, D2, V, αexc)
            elseif el == 1
                elem = Hd1(αind[Threads.threadid()], βind[Threads.threadid()], D1, D2, h, V, αexc)
            else
                elem = Hd0(αind[Threads.threadid()], βind[Threads.threadid()], h, V)
            end
            if abs(elem) > tol
                push!(vals[Threads.threadid()], elem)
                push!(ivals[Threads.threadid()], i)
                push!(jvals[Threads.threadid()], j)
            end
        end
    end #Threads.@spawn
    end

    ivals = vcat(ivals...)
    jvals = vcat(jvals...)
    vals  = vcat(vals...)
    return Symmetric(sparse(ivals, jvals, vals))
end


text:
type BatesProcess<: AbstractHestonProcess
  hestonProcess::HestonProcess
  lambda::Float64
  nu::Float64
  delta::Float64
  m::Float64
end

function BatesProcess(riskFreeRate::YieldTermStructure, dividendYield::YieldTermStructure, s0::Quote, v0::Float64,
                      kappa::Float64, theta::Float64, sigma::Float64, rho::Float64, lambda::Float64, nu::Float64,
                      delta::Float64, d::AbstractHestonDiscretization = FullTruncation())

  hp = HestonProcess(riskFreeRate, dividendYield, s0, v0, kappa, theta, sigma, rho, d)
  return BatesProcess(hp, lambda, nu, delta, expm1(nu + 0.5 * delta * delta))
end

get_risk_free_rate(bp::BatesProcess) = bp.hestonProcess.riskFreeRate
get_dividend_yield(bp::BatesProcess) = bp.hestonProcess.dividendYield
get_s0(bp::BatesProcess) = bp.hestonProcess.s0
get_v0(bp::BatesProcess) = bp.hestonProcess.v0
get_kappa(bp::BatesProcess) = bp.hestonProcess.kappa
get_theta(bp::BatesProcess) = bp.hestonProcess.theta
get_sigma(bp::BatesProcess) = bp.hestonProcess.sigma
get_rho(bp::BatesProcess) = bp.hestonProcess.rho
get_nu(bp::BatesProcess) = bp.nu
get_lambda(bp::BatesProcess) = bp.lambda
get_delta(bp::BatesProcess) = bp.delta
get_m(bp::BatesProcess) = bp.m


text:
#= License
Copyright 2019, 2020 (c) Yossi Bokor Katharine Turner

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

__precompile__()

module PersistentHomologyTransfer

#### Requirements ####
using CSV
using Hungarian
using DataFrames
using LinearAlgebra
using SparseArrays
using Eirene

#### Exports ####

export 	PHT,
		Recenter,
		Direction_Filtration,
		Evaluate_Barcode,
		Total_Rank_Exact,
		Total_Rank_Grid,
		Total_Rank_Auto,
		Average_Rank_Grid,
		Create_Heat_Map,
		Set_Mean_Zero,
		Weighted_Inner_Product,
		Weighted_Inner_Product_Matrix,
		Principal_Component_Scores,
		Average_Discretised_Rank,
		unittest

#### First some functions to recenter the curves ####
function Find_Center(points)
	n_p = size(points,1)

	c_x = Float64(0)
	c_y = Float64(0)

	for i in 1:n_p
		c_x += points[i,1]
		c_y += points[i,2]
	end

	return Float64(c_x/n_p), Float64(c_y/n_p)
end

function Recenter(points)

	points = convert(Array{Float64}, points)
	center = Find_Center(points)

	for i in 1:size(points)[1]
		points[i,1] = points[i,1] - center[1]
		points[i,2] = points[i,2] - center[2]
	end

	return points
end

function Evaluate_Rank(barcode, point)

	n = size(barcode)[1]
	count = 0

	if point[2] < point[1]
		#return 0
	else
		for i in 1:n
			if barcode[i,1] <= point[1]
				if barcode[i,2] >= point[2]
					count +=1
				end
			end
		end
		return count
	end
end

function Total_Rank_Exact(barcode)
	@assert size(barcode,2) == 2

	rks = []

	n = size(barcode,1)

	b = copy(barcode)
	reshape(b, 2,n)
	for i in 1:n
		for j in 1:n
			if barcode[i,1] < barcode[j,1]
				if barcode[i,2] < barcode[j,2]
					b = vcat(b, [barcode[j,1] barcode[i,2]])
					m += 1
				end
			end
		end
	end

	for i in 1:n
		append!(rks, Evaluate_Rank(barcode, b[i,:]))
	end
	return b, rks

end

function Total_Rank_Grid(barcode, x_g, y_g) #the grid should be an array, with 0s in all entries below the second diagonal.

		@assert size(x_g) == size(y_g) #I should maybe change this as you don't REALLY need to use the same size grid.....

		n_g = size(x_g,1)
		rks = zeros(n_g,n_g)
		n_p = size(barcode,1)

		for i in 1:n_p
			point = barcode[i,:]
			x_i = findfirst(>=(point[1]), x_g)
			y_i = findfirst(<=(point[2]), y_g)
			for j in x_i:n_g-y_i+1
				for k in j:n_g-y_i+1
					rks[n_g-k+1,j] += 1
				end
			end
		end

	return rks

end

function Average_Rank_Grid(list_of_barcodes, x_g, y_g)

	rks = zeros(length(x_g),length(y_g))

	n_b = length(list_of_barcodes)

	for i in 1:n_b
		rk_i = Total_Rank_Grid(list_of_barcodes[i], x_g,y_g)
		rks = rks .+ rk_i
	end

	rks = rks/n_b

	return rks
end

function Average_Rank_Point(list_of_barcodes, x,y)

	rk = 0
	n_b = length(list_of_barcodes)
	if y >= x
	for i in 1:n_b
			rk += Evaluate_Rank(list_of_barcodes[i], [x,y])
		end
		return rk/n_b
	else
		 return 0
	end
end

function Create_Heat_Map(barcode, x_g, y_g)

	f(x,y) =  begin
					if x > y
						return 0
					else
						return Evaluate_Rank(barcode,[x,y])
					end
				end

	#Z = map(f, X, Y)

	p1 = contour(x_g, y_g, f, fill=true)

	return p1
end

# Let us do PCA for the rank functions using Kate and Vanessa's paper.
# So, I first need to calculate the pointwise norm

function Set_Mean_Zero(discretised_ranks)
	n_r = length(discretised_ranks)
	println(n_r)
	grid_size = size(discretised_ranks[1])

	for i in 1:n_r
		@assert size(discretised_ranks[i]) == grid_size
	end

	mu = zeros(grid_size)

	for i in 1:n_r
		mu = mu .+ discretised_ranks[i]
	end

	mu = mu./n_r

	normalised = []
	for i in 1:n_r
		append!(normalised, [discretised_ranks[i] .- mu])
	end
	return normalised
end

function Weighted_Inner_Product(disc_rank_1, disc_rank_2, weights)

	wip = sum((disc_rank_1.*disc_rank_2).*weights)

	return wip
end

function Weighted_Inner_Product_Matrix(discretised_ranks, weights)
	n_r = length(discretised_ranks)
	D = Array{Float64}(undef, n_r, n_r)

	for i in 1:n_r
		for j in i:n_r
			wip = Weighted_Inner_Product(discretised_ranks[i], discretised_ranks[j], weights)
			D[i,j] = wip
			D[j,i] = wip
		end
	end

	return D
end

function Principal_Component_Scores(inner_prod_matrix, dimension)
	F = LinearAlgebra.eigen(inner_prod_matrix, permute = false, scale=false) # this sorts the eigenvectors in ascending order
	n_r = size(inner_prod_matrix,1)
	lambda = Array{Float64}(undef, 1,dimension)
	w = Array{Float64}(undef, size(F.vectors)[1],dimension)
	n_v = length(F.values)

	for i in 1:dimension
		lambda[i] = F.values[n_v-i+1]
		w[:,i] = F.vectors[:,n_v-i+1]
	end

	s = Array{Float64}(undef, n_r,dimension)

	for i in 1:size(inner_prod_matrix,1)
		for j in 1:dimension

			den = sqrt(sum([w[k,j]*sum(w[l,j]*inner_prod_matrix[k,l] for l in 1:n_r) for k in 1:n_r]))
			numerator = sum(w[:,j].*inner_prod_matrix[:,i])
			s[i,j] = numerator/den
		end
	end
	return s
end

function Average_Discretised_Rank(list_of_disc_ranks)
	average = Array{Float64}(undef, size(list_of_disc_ranks[1]))
	n_r = length(list_of_disc_ranks)

	for i in n_r
		average = average .+ list_of_disc_ranks[i]
	end

	return average/n_r
end

function Direction_Filtration(ordered_points, direction; out = "barcode")
	number_of_points = length(ordered_points[:,1]) #number of points
	heights = zeros(number_of_points) #empty array to be changed to heights for filtration
	fv = zeros(2*number_of_points) #blank fv Eirene
	for i in 1:number_of_points
		heights[i]= ordered_points[i,1]*direction[1] + ordered_points[i,2]*direction[2] #calculate heights in specificed direction
	end

	for i in 1:number_of_points
		fv[i]= heights[i] # for a point the filtration step is the height
	end

	for i in 1:(number_of_points-1)
		fv[(i+number_of_points)]=maximum([heights[i], heights[i+1]]) # for an edge between two adjacent points it enters when the 2nd of the two points does
	end

	fv[2*number_of_points] = maximum([heights[1] , heights[number_of_points]]) #last one is a special snowflake
	dv = [] # template dv for Eirene

	for i in 1:number_of_points
		append!(dv,0) # every point is 0 dimensional
	end

	for i in (1+number_of_points):(2*number_of_points)
		append!(dv,1) # edges are 1 dimensional
	end

	D = zeros((2*number_of_points, 2*number_of_points))

	for i in 1:number_of_points
		D[i,(i+number_of_points)]=1 # create boundary matrix and put in entries
	end

	for i in 2:(number_of_points)
		D[i, (i+number_of_points-1)]=1 # put in entries for boundary matrix
	end

	D[1, (2*number_of_points)]=1

	ev = [number_of_points, number_of_points] # template ev for Eirene

	S  = sparse(D) 	# converting as required for Eirene
	rv = S.rowval 	# converting as required for Eirene
	cp = S.colptr	# converting as required for Eirene
	C = Eirene.eirene(rv=rv,cp=cp,ev=ev,fv=fv) # put it all into Eirene

	if out == "barcode"
		return barcode(C, dim=0)
	elseif out == "one_cycle"
		return barcode(C, dim=0), maximum(heights)
	else
		return C
	end
end

#### Wrapper for the PHT function ####

function PHT(curve_points, directions; one_cycle = "n") ##accepts an ARRAY of points

	if typeof(directions) ==  Int64
		println("auto generating directions")
		dirs = Array{Float64}(undef, directions,2)
		for n in 1:directions
			dirs[n,1] = cos(n*pi/(directions/2))
			dirs[n,2] = sin(n*pi/(directions/2))
		end
		println("Directions are:")
		println(dirs)
	else
		println("using directions provided")
		dirs = copy(directions)
	end
	pht = []

	if one_cycle == "y"
		c_1 = []
	end

	for i in 1:size(dirs,1)

		if one_cycle == "y"
			pd,c_1 = Direction_Filtration(curve_points, dirs[i,:], out ="one_cycle")
			append!(cycle_1, c_1)
			pht = vcat(pht, [pd])
		else
			pd = Direction_Filtration(curve_points, dirs[i,:])
			pht = vcat(pht, [pd])
		end
	end

	if one_cycle == "y"
		return pht, cycle_1

	else
		return pht
	end
end

#### Wrapper for PCA ####
function PCA(ranks, dimension, weights)

	normalised = Set_Mean_Zero(ranks)

	D = Weighted_InnerProd_Matrix(normalised, weights)

	return Principal_Component_Scores(D, dimension)
end

#### Unittests ####
function test_1()
	return PHT([0,0,0],0)
end

function test_2()
	pht = PHT([1 1; 5 5], 1)
	if pht == [0.9999999999999998 4.999999999999999]
		return []
	else
		println("Error: test_2, pht = ")
		return pht
	end
end

function unittest()

	x = Array{Any}(undef, 2)

	x[1] = test_1()
	x[2] = test_2()

	for p 	= 	1:length(x)
		if !isempty(x[p])
			println(p)
			return x
		end
	end
	return []
end

end# module


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs
import QuantumOpticsBase

@testset "initial state" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, Energy)
    mol2 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    vibindices = getVibIndices(agg, aggInds)
    aggIndsLen = length(aggInds)
    basis = GenericBasis([aggIndsLen])
    FCFact = getFranckCondonFactors(agg, aggInds)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    T = 300.0
    mu_array = [[1, 2]]
    W0_ref = [
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25023985364402346 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24999994249781035 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24999994249781035 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2497602613603558
    ]
    # resulting in NaNs
    #=
    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_composite_ref, W0)
    =#

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-3 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-15 > D(W0_composite_ref, W0)

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(
        T,
        [0.0, 0.8, 0.2],
        Ham,
        vibindices,
        aggInds;
        diagonalize = false,
        diagonal = true,
    )
    @test 1e-15 > D(W0_composite_ref, W0)

    T = 1e-10
    W0_ref = [
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
    ]
    # resulting in NaNs
    #=
    println(Ham)
    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false)
    println(W0)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_composite_ref, W0)

    # W0 = thermal_state(T, mu_array, Ham, aggInds; diagonalize=true)
    # @test 1e-15 > D(W0_ref, W0.data)
    # println(W0.data)
    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(
        T,
        [0.0, 0.8, 0.2],
        Ham,
        vibindices,
        aggInds;
        diagonalize = false,
        diagonal = true,
    )
    @test 1e-15 > D(W0_composite_ref, W0)
    =#

    T = 300.0
    mu_array = [[1, 2], [2, 1]]
    W0_ref = [
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25023985364402346 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24999994249781035 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24999994249781035 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24976026136035578
    ]
    # resulting in NaNs
    #=
    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_composite_ref, W0)
    =#

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-3 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-15 > D(W0_composite_ref, W0)

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(
        T,
        [0.0, 0.8, 0.2],
        Ham,
        vibindices,
        aggInds;
        diagonalize = false,
        diagonal = true,
    )
    @test 1e-15 > D(W0_composite_ref, W0)

    aggInds = getIndices(agg)
    FCFact = getFranckCondonFactors(agg, aggInds)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)
    T = 300.0
    mu_array = [[1, 1]]
    W0_ref = [
        0.25023985364402346 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.24999994249781035 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.24999994249781035 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.24976026136035578 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
    ]
    #=
    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = false)
    @test 1e-15 > D(W0_composite_ref, W0)
    =#

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-3 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, vibindices, aggInds; diagonalize = true)
    @test 1e-15 > D(W0_composite_ref, W0)

    W0 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    @test 1e-15 > D(W0_ref, W0.data)
    W01 = thermal_state(T, [[2, 1]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W02 = thermal_state(T, [[1, 2]], Ham, vibindices, aggInds; diagonalize = false, diagonal = true)
    W0_composite_ref = 0.8 * W01 + 0.2 * W02
    normalize!(W0_composite_ref)
    W0 = thermal_state_composite(
        T,
        [0.0, 0.8, 0.2],
        Ham,
        vibindices,
        aggInds;
        diagonalize = false,
        diagonal = true,
    )
    @test 1e-15 > D(W0_composite_ref, W0)

end # testset


text:
"""
    center_predictors(X::AbstractMatrix)

Centers the columns of a matrix `X` of predictors to mean 0.

Returns a tuple with:
1. `μ_X`: 1xK `Matrix` of `Float64`s of the means of the K columns in the original `X`
matrix.
2. `X_centered`: A `Matrix` of `Float64`s with the same dimensions as the original matrix
`X` with the columns centered on mean μ=0.

# Arguments
- `X::AbstractMatrix`: a matrix of predictors where rows are observations and columns are
variables.
"""
function center_predictors(X::AbstractMatrix)
    μ_X = mapslices(mean, X; dims=1)
    X_centered = similar(X, Float64)
    for (idx, val) in enumerate(eachcol(X))
        X_centered[:, idx] = val .- μ_X[1, idx]
    end
    return vec(μ_X), X_centered
end

"""
    standardize_predictors(X::AbstractMatrix)

Standardizes the columns of a matrix `X` of predictors to mean 0 and standard deviation 1.

Returns a tuple with:
1. `μ_X`: 1xK `Matrix` of `Float64`s of the means of the K columns in the original `X`
matrix.
2. `σ_X`: 1xK `Matrix` of `Float64`s of the standard deviations of the K columns in the
original `X` matrix.
3. `X_std`: A `Matrix` of `Float64`s with the same dimensions as the original matrix
`X` with the columns centered on mean μ=0 and standard deviation σ=1.

# Arguments
- `X::AbstractMatrix`: a matrix of predictors where rows are observations and columns are
variables.
"""
function standardize_predictors(X::AbstractMatrix)
    μ_X = mapslices(mean, X; dims=1)
    σ_X = mapslices(std, X; dims=1)
    X_std = similar(X, Float64)
    for (idx, val) in enumerate(eachcol(X))
        X_std[:, idx] = (val .- μ_X[1, idx]) ./ σ_X[1, idx]
    end
    return vec(μ_X), vec(σ_X), X_std
end

"""
    standardize_predictors(x::AbstractVector)

Standardizes the vector `x` to mean 0 and standard deviation 1.

Returns a tuple with:
1. `μ_X`: `Float64`s of the mean of the original vector `x`.
2. `σ_X`: `Float64`s of the standard deviations of the original vector `x`.
3. `x_std`: A `Vector` of `Float64`s with the same length as the original vector
`x` with the values centered on mean μ=0 and standard deviation σ=1.

# Arguments
- `x::AbstractVector`: a vector.
"""
function standardize_predictors(x::AbstractVector)
    μ_x = mean(x)
    σ_x = std(x)
    x_std = (x .- μ_x) ./ σ_x
    return μ_x, σ_x, x_std
end

"""
    tuple_length(::NTuple{N, Any}) where {N} = Int(N)

This is a hack to get the length of any tuple.
"""
tuple_length(::NTuple{N,Any}) where {N} = Int(N)

"""
    convert_str_to_indices(v::AbstractVector)

Converts a vector `v` to a vector of indices, i.e. a vector where all the entries are
integers. Returns a tuple with the first element as the converted vector and the
second element a `Dict` specifying which string is which integer.

This function is especially useful for random-effects varying-intercept hierarchical models.
Normally `v` would be a vector of group membership with values such as `"group_1"`,
`"group_2"` etc. For random-effect models with varying-intercepts, Turing needs the group
membership values to be passed as `Int`s.
"""
function convert_str_to_indices(v::AbstractVector)
    d = Dict{eltype(v),Int}()
    v_int = Int[]
    for i in v
        n = get!(d, i, length(d) + 1)
        push!(v_int, n)
    end
    return v_int, d
end


text:
using HDF5

macro gcvalid(args...)
	Expr(:block, quote
		gc_enable()
		gc()
		gc_disable()
	end,
	[:(@assert HDF5.isvalid($x)) for x in args]...)
end

macro closederror(x)
	quote
		try
			$x
		catch e
			isa(e, ErrorException) || rethrow(e)
			e.msg == "File or object has been closed" || error("Attempt to access closed object did not throw")
		end
	end
end

gc_disable()
fn = joinpath(tempdir(),"test.h5")
for i = 1:10
	file = h5open(fn, "w")
    memtype_id = HDF5.h5t_create(HDF5.H5T_COMPOUND, 2*sizeof(Float64))
    HDF5.h5t_insert(memtype_id, "real", 0, HDF5.hdf5_type_id(Float64))
    HDF5.h5t_insert(memtype_id, "imag", sizeof(Float64), HDF5.hdf5_type_id(Float64))
    dt = HDF5Datatype(memtype_id)
    t_commit(file, "dt", dt)
    ds = dataspace((2,))
	d = d_create(file, "d", dt, ds)
	g = g_create(file, "g")
	a = a_create(file, "a", dt, ds)
	@gcvalid dt ds d g a
	close(file)

	@closederror read(d)
	for obj in (d, g)
		@closederror a_read(obj, "a")
		@closederror a_write(obj, "a", 1)
	end
	for obj in (g, file)
		@closederror d_open(obj, "d")
		@closederror d_read(obj, "d")
		@closederror d_write(obj, "d", 1)
		@closederror read(obj, "x")
		@closederror write(obj, "x", "y")
	end
end
for i = 1:10
	file = h5open(fn, "r")
	dt = file["dt"]
	d = file["d"]
	ds = dataspace(d)
	g = file["g"]
	a = attrs(file)["a"]
	@gcvalid dt ds d g a
	close(file)
end


text:
module FrechetDiff

using LinearAlgebra
using ForwardDiff

include("graph.jl")
include("program.jl")
include("emit_code.jl")
include("differentiate.jl")
include("rules.jl")
include("materialize.jl")

export record_function, differentiate, emit_code, materialize_hessian

end # module


text:
using ContinuumArrays, LinearAlgebra, LazyArrays, Base64, FillArrays, QuasiArrays, BandedMatrices, Test
using QuasiArrays: ApplyQuasiArray, ApplyStyle, MemoryLayout, mul, MulQuasiMatrix, Vec
import LazyArrays: MulStyle, LdivStyle, arguments
import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, SubBasisLayout, AdjointMappedBasisLayout, MappedBasisLayout, coefficients

@testset "Splines" begin
    @testset "HeavisideSpline" begin
        H = HeavisideSpline([1,2,3])

        @test axes(H) ≡ (axes(H,1),axes(H,2)) ≡ (Inclusion(1..3), Base.OneTo(2))
        @test size(H) ≡ (size(H,1),size(H,2)) ≡ (ℵ₁, 2)

        @test_throws BoundsError H[0.1, 1]
        @test H[1.1,1] ≡ H'[1,1.1] ≡ transpose(H)[1,1.1] ≡ 1.0
        @test H[2.1,1] ≡ H'[1,2.1] ≡ transpose(H)[1,2.1] ≡ 0.0
        @test H[1.1,2] ≡ 0.0
        @test H[2.1,2] ≡ 1.0
        @test_throws BoundsError H[2.1,3]
        @test_throws BoundsError H'[3,2.1]
        @test_throws BoundsError transpose(H)[3,2.1]
        @test_throws BoundsError H[3.1,2]

        @test all(H[[1.1,2.1], 1] .=== H'[1,[1.1,2.1]] .=== transpose(H)[1,[1.1,2.1]] .=== [1.0,0.0])
        @test all(H[1.1,1:2] .=== H[1.1,:] .=== [1.0,0.0])
        @test all(H[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])

        @test_throws BoundsError H[[0.1,2.1], 1]
        @test MemoryLayout(typeof(H)) == BasisLayout()
        @test ApplyStyle(*, typeof(H), typeof([1,2])) isa MulStyle

        f = H*[1,2]
        @test f isa ApplyQuasiArray
        @test axes(f) == (Inclusion(1.0..3.0),)
        @test f[1.1] ≈ 1
        @test f[2.1] ≈ 2

        @test @inferred(H'H) == @inferred(materialize(applied(*,H',H))) == Eye(2)
        @test summary(f) == stringmime("text/plain", f) == "HeavisideSpline([1, 2, 3]) * [1, 2]"

        @testset "sum/cumsum" begin
            H = HeavisideSpline(range(0,1;length=1000));
            x = axes(H,1)
            @test sum(H/H \ exp.(x)) ≈ ℯ-1 atol=1E-5
            @test last(cumsum(H/H \ exp.(x))) ≈ sum(H/H\exp.(x))
        end
    end

    @testset "LinearSpline" begin
        @testset "Evaluation" begin
            L = LinearSpline([1,2,3])
            @test size(L) == (ℵ₁, 3)

            @test_throws BoundsError L[0.1, 1]
            @test L[1.1,1] == L'[1,1.1] == transpose(L)[1,1.1] ≈ 0.9
            @test L[2.1,1] === L'[1,2.1] === transpose(L)[1,2.1] === 0.0
            @test L[1.1,2] ≈ 0.1
            @test L[2.1,2] ≈ 0.9
            @test L[2.1,3] == L'[3,2.1] == transpose(L)[3,2.1] ≈ 0.1
            @test_throws BoundsError L[3.1,2]

            @test L[[1.1,2.1], 1] == L'[1,[1.1,2.1]] == transpose(L)[1,[1.1,2.1]] ≈ [0.9,0.0]
            @test L[1.1,1:2] ≈ [0.9,0.1]
            @test L[[1.1,2.1], 1:2] ≈ [0.9 0.1; 0.0 0.9]

            @test_throws BoundsError L[[0.1,2.1], 1]
        end

        @testset "Expansion" begin
            L = LinearSpline([1,2,3])
            f = L*[1,2,4]

            @test basis(f) == L
            @test coefficients(f) == [1,2,4]
            @test axes(f) == (Inclusion(1.0..3.0),)
            @test f[1.1] ≈ 1.1
            @test f[2.1] ≈ 2.2

            δ = DiracDelta(1.2,1..3)
            L = LinearSpline([1,2,3])
            @test @inferred(δ'L) ≈ [0.8, 0.2, 0.0]
            @test δ'f == mul(δ',f) == dot(δ,f) == f[1.2]

            @test @inferred(L'L) == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])

            @testset "Algebra" begin
                @test 2f == f*2 == 2 .* f == f .* 2
                @test 2\f == f/2 == 2 .\ f == f ./ 2
                @test sum(f) ≈ 4.5
                @test L \ BroadcastQuasiArray(*,2,f) ≈ [2,4,8]
            end
        end

        @testset "==" begin
            L = LinearSpline([1,2,3])
            H = HeavisideSpline([1,2,3])
            @test L == L
            @test L ≠ H
            H = HeavisideSpline([1,1.5,2.5,3])
            @test_throws ArgumentError L == H
        end

        @testset "Adjoint layout" begin
            L = LinearSpline([1,2,3])
            @test MemoryLayout(L') isa AdjointBasisLayout
            @test @inferred([3,4,5]'*L') isa ApplyQuasiArray
        end

        @testset "+" begin
            L = LinearSpline([1,2,3])
            b = @inferred(L*[3,4,5] + L*[1.,2,3])
            @test ApplyStyle(\, typeof(L), typeof(b)) == LdivStyle()
            @test (L\b) == [4,6,8]
            B = BroadcastQuasiArray(+, L, L)
            @test L\B == 2Eye(3)

            @test L*[3,4,5] + L*[1.,2,3] + L*[4,5,6] == @inferred(broadcast(+, L*[3,4,5], L*[1.,2,3], L*[4,5,6])) == L*[8,11,14]

            b = L*[3,4,5] - L*[1.,2,3]
            @test (L\b) == [2,2,2]
            B = BroadcastQuasiArray(-, L, L)
            @test L\B == 0Eye(3)

            @testset "sub" begin
                v = ApplyQuasiArray(*, L[:,2:end], [1,2])
                f = L * [1,2,3]
                @test v + f == f + v == L*[1,3,5]
                @test v + v == L*[0,2,4]
            end
        end
    end

    @testset "Algebra" begin
        L = LinearSpline([1,2,3])
        f = L*[1,2,4]
        g = L*[5,6,7]

        @test MemoryLayout(f) isa ExpansionLayout
        @test MemoryLayout(2f) isa ExpansionLayout
        @test MemoryLayout(f*2) isa ExpansionLayout
        @test MemoryLayout(2\f) isa ExpansionLayout
        @test MemoryLayout(f/2) isa ExpansionLayout
        @test MemoryLayout(f+g) isa ExpansionLayout
        @test MemoryLayout(f-g) isa ExpansionLayout
        @test f[1.2] == 1.2
        @test (2f)[1.2] == (f*2)[1.2] == 2.4
        @test (2\f)[1.2] == (f/2)[1.2] == 0.6
        @test (f+g)[1.2] ≈ f[1.2] + g[1.2]
        @test (f-g)[1.2] ≈ f[1.2] - g[1.2]
    end

    @testset "Derivative" begin
        L = LinearSpline([1,2,3])
        f = L*[1,2,4]

        D = Derivative(axes(L,1))
        @test copy(D) == D

        @test D*L isa MulQuasiMatrix
        @test length((D*L).args) == 2
        @test eltype(D*L) == Float64

        M = applied(*, (D*L).args..., [1,2,4])
        @test eltype(materialize(M)) == Float64

        M = applied(*, D, L, [1,2,4])
        @test materialize(M) isa ApplyQuasiArray

        fp = D*L*[1,2,4]

        @test eltype(fp) == Float64

        @test fp isa Vec
        @test length(fp.args) == 2
        @test fp[1.1] ≈ 1
        @test fp[2.2] ≈ 2

        fp = D*f
        @test length(fp.args) == 2
        @test fp[1.1] ≈ 1
        @test fp[2.2] ≈ 2

        @testset "View derivatives" begin
            L = LinearSpline(1:5)
            H = HeavisideSpline(1:5)
            x = axes(L,1)
            D = Derivative(x)

            @test view(D, x, x) ≡ D
            @test_throws BoundsError view(D, Inclusion(0..1), x)
            jr = [1,3,4]
            @test H \ (D*L[:,jr]) == H\(L[:,jr]'D')' == (H \ (D*L))[:,jr]

            @test D*L[:,jr]*[1,2,3] == D * L * [1,0,2,3,0]
            @test L[:,jr]'D'D*L[:,jr] == (L'D'D*L)[jr,jr]

            a = affine(0..1, 1..5)
            D̃ = Derivative(axes(a,1))
            @test H[a,:] \ (D̃ * L[a,:]) == H[a,:] \ (L[a,:]'D̃')' ==  4*(H\(D*L))
            @test_throws DimensionMismatch D * L[a,:]
            @test_throws DimensionMismatch D̃ * L[:,jr]

            @test ContinuumArrays.simplifiable(*, D, L[:,jr]) isa Val{true}
            @test ContinuumArrays.simplifiable(*, L[:,jr]', D') isa Val{true}
            @test ContinuumArrays.simplifiable(*, D̃, L[a,jr]) isa Val{true}
            @test ContinuumArrays.simplifiable(*, L[a,jr]', D̃') isa Val{true}
        end

        @testset "maps broadcasted" begin
            L = LinearSpline(1:5)
            a = affine(0..1, 1..5)
            M = L[a,:]
            M̃ = L[a,1:5]
            x = axes(M,1)

            @test M == M
            @test M == M̃
            @test M̃ == M

            @test (x .* M)[0.25,:] ≈ (x .* M̃)[0.25,:] ≈ 0.25 * M[0.25,:]
            @test (exp.(x) .* M)[0.25,:] ≈ exp(0.25) * M[0.25,:]

            @test M \ [ones(x) x] ≈ [M\ones(x) M\x]
        end
    end

    @testset "Weak Laplacian" begin
        H = HeavisideSpline(0:2)
        L = LinearSpline(0:2)
        D = Derivative(axes(L,1))

        @test apply(*,L',D') isa MulQuasiMatrix
        @test MemoryLayout(typeof(L')) isa AdjointBasisLayout
        @test (L'D') isa MulQuasiMatrix

        A = (L'D') * (D*L)
        @test A isa BandedMatrix
        @test A == (D*L)'*(D*L) == [1.0 -1 0; -1.0 2.0 -1.0; 0.0 -1.0 1.0]
        @test bandwidths(A) == (1,1)
    end

    @testset "Views" begin
        L = LinearSpline(0:2)
        @test view(L,0.1,1)[1] == L[0.1,1]

        L = LinearSpline(0:2)
        B1 = view(L,:,1)
        @test B1 isa SubQuasiArray{Float64,1}
        @test size(B1) == (ℵ₁,)
        @test B1[0.1] == L[0.1,1]
        @test_throws BoundsError B1[2.2]

        B = view(L,:,1:2)
        @test B isa SubQuasiArray{Float64,2}
        @test B[0.1,:] == L[0.1,1:2]

        B = @view L[:,2:end-1]
        @test B[0.1,:] == [0.1]

        L = LinearSpline([1,2,3,4])
        @test L[:,2:3] isa SubQuasiArray
        @test axes(L[:,2:3]) ≡ (Inclusion(1..4), Base.OneTo(2))
        @test L[:,2:3][1.1,1] == L[1.1,2]
        @test_throws BoundsError L[0.1,1]
        @test_throws BoundsError L[1.1,0]

        @test MemoryLayout(typeof(L[:,2:3])) isa SubBasisLayout
        @test L\L[:,2:3] isa BandedMatrix
        @test L\L[:,2:3] == [0 0; 1 0; 0 1.0; 0 0]
        @test L[:,2:3]\L isa BandedMatrix
        @test L[:,2:3]\L == [0 1 0 0; 0 0 1 0]

        @testset "Subindex of splines" begin
            L = LinearSpline(range(0,stop=1,length=10))
            @test L[:,2:end-1] isa SubQuasiArray
            @test L[:,2:end-1][0.1,1] == L[0.1,2]
            v = randn(8)
            f = L[:,2:end-1] * v
            @test f[0.1] ≈ (L*[0; v; 0])[0.1]
        end

        @testset "sub-of-sub" begin
            L = LinearSpline([1,2,3])
            V = view(L,:,1:2)
            @test V == V
            V2 = view(V,1.1:0.1:2,:)
            @test V2 == L[1.1:0.1:2,1:2]
        end

        @testset "sub-colon" begin
            L = LinearSpline([1,2,3])
            @test L[:,1][1.1] == (L')[1,:][1.1] == L[1.1,1]
            @test L[:,1:2][1.1,:] == (L')[1:2,:][:,1.1] == L[1.1,1:2]
        end

        @testset "transform" begin
            L = LinearSpline([1,2,3])
            x = axes(L,1)
            @test (L \ x) == pinv(L)x == [1,2,3]
            @test factorize(L[:,2:end-1]) isa ContinuumArrays.ProjectionFactorization
            @test L[:,1:2] \ x == [1,2]
            @test L \ [x one(x)] ≈ [L\x L\one(x)]
            @test factorize(L) \ QuasiOnes(x, Base.OneTo(3)) ≈ L \ QuasiOnes(x, Base.OneTo(3)) ≈ ones(3,3)

            L = LinearSpline(range(0,1; length=10_000))
            x = axes(L,1)
            @test L[0.123,:]'* (L \ exp.(x)) ≈ exp(0.123) atol=1E-9
            @test L[0.123,2:end-1]'* (L[:,2:end-1] \ exp.(x)) ≈ exp(0.123) atol=1E-9

            @test L \ zeros(x) ≡ Zeros(10_000)
        end
    end

    @testset "sum" begin
        H = HeavisideSpline([1,2,3,6])
        @test sum(H; dims=1) * [1,1,1] == [5]
        x = Inclusion(0..1)
        B = H[5x .+ 1,:]
        @test sum(B; dims=1) * [1,1,1] == [1]
        @test sum(H[:,1:2]; dims=1) * [1,1] == [2]
        @test sum(H'; dims=2) == permutedims(sum(H; dims=1))

        u = H * randn(3)
        @test sum(u[5x .+ 1]) ≈ sum(view(u,5x .+ 1)) ≈ sum(u)/5

        L = LinearSpline([1,2,3,6])
        D = Derivative(axes(L,1))
        @test sum(D*L; dims=1) ≈ sum((D*L)'; dims=2)' ≈ [-1 zeros(1,2) 1]
    end

    @testset "Poisson" begin
        L = LinearSpline(range(0,stop=1,length=10))
        B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
        @test B'B == (L'L)[2:end-1,2:end-1]
        D = Derivative(axes(L,1))
        @test apply(*,D,B) isa MulQuasiMatrix
        @test D*B isa MulQuasiMatrix
        @test apply(*,D,B)[0.1,1] == (D*B)[0.1,1] == 9
        @test apply(*,D,B)[0.2,1] == (D*B)[0.2,1] == -9
        Δ = -((D*B)'*(D*B)) # Weak Laplacian
        @test Δ isa BandedMatrix

        @test B'D' isa MulQuasiMatrix
        @test length((B'D').args) == 2

        @test *(B',D',D,B) isa BandedMatrix

        @test Δ == -(*(B',D',D,B))
        @test Δ == -(B'D'D*B)
        @test Δ == -((B'D')*(D*B))
        @test_broken Δ == -B'*(D'D)*B
        @test Δ == -(B'*(D'D)*B)

        f = L*exp.(L.points) # project exp(x)
        u = B * (Δ \ (B'f))

        @test u[0.1] ≈ -0.06612902692412974
    end

    @testset "Helmholtz" begin
        L = LinearSpline(range(0,stop=1,length=10))
        B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
        D = Derivative(axes(L,1))

        A = -((B'D')*(D*B)) + 100^2*B'B # Weak Laplacian

        f = L*exp.(L.points) # project exp(x)
        u = B * (A \ (B'f))

        @test u[0.1] ≈ 0.00012678835289369413
    end

    @testset "Change-of-variables" begin
        x = Inclusion(0..1)
        y = 2x .- 1
        L = LinearSpline(range(-1,stop=1,length=10))
        @test L[y,:][0.1,:] == L[2*0.1-1,:]
        @test L[y,:]'L[y,:] isa SymTridiagonal
        @test L[y,:]'L[y,:] == 1/2*(L'L)

        D = Derivative(axes(L,1))
        H = HeavisideSpline(L.points)
        @test H\((D*L) * 2) ≈ (H\(D*L))*2 ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

        @test MemoryLayout(L[y,:]) isa MappedBasisLayout
        @test MemoryLayout(L[y,:]') isa AdjointMappedBasisLayout
        a,b = arguments((D*L)[y,:])
        @test H[y,:]\a == Eye(9)
        @test H[y,:] \ (D*L)[y,:] isa BandedMatrix
        @test @inferred(grid(L[y,:])) ≈ (grid(L) .+ 1) ./ 2

        D = Derivative(x)
        @test (D*L[y,:])[0.1,1] ≈ -9
        @test H[y,:] \ (D*L[y,:]) isa BandedMatrix
        @test H[y,:] \ (D*L[y,:]) ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

        B = L[y,2:end-1]
        @test MemoryLayout(typeof(B)) isa MappedBasisLayout
        @test B[0.1,1] == L[2*0.1-1,2]
        @test B\B == Eye(8)
        @test L[y,:] \ B == Eye(10)[:,2:end-1]
        @test B'B == (1/2)*(L'L)[2:end-1,2:end-1]

        @testset "algebra" begin
            f = L[y,:] * randn(10)
            g = L[y,:] * randn(10)
            @test MemoryLayout(f + g) isa ExpansionLayout
            @test (f+g)[0.1] ≈ f[0.1] + g[0.1]
        end

        @testset "vec demap" begin
            @test L[y,:] \ exp.(axes(L,1))[y] ≈ L[y,:] \ exp.(y) ≈  factorize(L[y,:]) \ exp.(y)
            @test ContinuumArrays.demap(view(axes(L,1),y)) == axes(L,1)

            @test L[y,:] \ (y .* exp.(y)) ≈ L[y,:] \ BroadcastQuasiVector(y -> y*exp(y), y)
            @test L[y,:] \ (y .* L[y,1:3]) ≈ [L[y,:]\(y .* L[y,1]) L[y,:]\(y .* L[y,2]) L[y,:]\(y .* L[y,3])]

            c = randn(size(L,2))
            @test L[y,:] \ (L[y,:] * c) ≈ c
            @test ContinuumArrays.demap(L[y,:] * c) == L*c
        end

        @testset "Mapped and BroadcastLayout{typeof(+)}" begin
            @test L[y,:] \ (y .+ y) ≈ L[y,:] \ (2y)
            @test L[y,:] \ (y .- y) ≈ zeros(10)
        end
    end

    @testset "diff" begin
        L = LinearSpline(range(-1,stop=1,length=10))
        f = L * randn(size(L,2))
        h = 0.0001;
        @test diff(f)[0.1] ≈ (f[0.1+h]-f[0.1])/h
    end

    @testset "show" begin
        x = Inclusion(0..1)
        H = HeavisideSpline([1,2,3,6])
        B = H[5x .+ 1,:]
        u = H * [1,2,3]
        @test stringmime("text/plain", B) == "HeavisideSpline([1, 2, 3, 6]) affine mapped to 0..1"
    end

    @testset "A \\ ( c .* B) == c .* (A\\B) #101" begin
        L = LinearSpline(0:5)
        @test L \ (2L) == 2(L\L)
    end

    @testset "H\\D*L" begin
        L = LinearSpline(0:5)
        H = HeavisideSpline(L)
        x = axes(L,1)
        D = Derivative(x)
        @test H\D*L == H\(D*L)
    end
end


text:
module NastyGenerator

using Random

pn(i) = "P$i"
randvers(k::Int) = VersionNumber(rand(0:k), rand(0:k), rand(0:k))

function randvspec(k::Int)
    lb, ub = randvers(k), randvers(k)
    if lb > ub
        lb, ub = ub, lb
    end

    slb = rand() < 0.1 ? "0" : string(lb)
    sub = rand() < 0.1 ? "*" : string(ub)

    return "$slb-$sub"
end

"""
Generates a random graph with 2 planted solutions (or quasi-solutions if sat==false).
With the right parameters, this is quite hard to solve for large sizes.
For added fun, the two planted solutions are cycles.
We want the solver to find the best of the two, of course.
This is an extremely unrealistic scenario and it's only intended to stress-test the solver.
Note that the "problematic" output assumes that all non-planted versions will be
uninstallable, which is only the case for some regimes of the parameters (e.g. large
enough d).
"""
function generate_nasty(n::Int,             # size of planted solutions
                        m::Int;             # size of the graph
                        k::Int = 10,        # version number limit
                        q::Int = 10,        # versions per package (upper bound)
                        d::Int = 10,        # neighbors per package
                        seed::Integer = 32524,
                        sat::Bool = true    # create a satisfiable problem?
                       )
    @assert m ≥ n
    d ≤ m-1 || @warn "d=$d, should be ≤ m-1=$(m-1)"

    srand(seed)

    allvers = [sort(unique(randvers(k) for j = 1:q)) for i = 1:m]

    planted1 = [rand(2:length(allvers[i])) for i = 1:n]

    planted2 = [rand(1:(planted1[i]-1)) for i = 1:n]

    deps = []
    problematic = []

    # random dependencies
    for i = 1:m, j = 1:length(allvers[i])
        if i ≤ n && (planted1[i] == j || planted2[i] == j)
            if j == planted1[i]
                if i < n
                    push!(deps, [pn(i), allvers[i][j], pn(i+1), "$(allvers[i+1][planted1[i+1]])-*"])
                else
                    if !sat
                        push!(deps, [pn(i), allvers[i][j], pn(1), "0-$(allvers[1][planted2[1]])"])
                    else
                        push!(deps, [pn(i), allvers[i][j], pn(1), "0-*"])
                    end
                end
            else # j == planted2[i]
                if i < n
                    push!(deps, [pn(i), allvers[i][j], pn(i+1), "0-$(allvers[i+1][planted2[i+1]])"])
                else
                    if !sat
                        push!(deps, [pn(i), allvers[i][j], pn(1), "$(allvers[1][planted1[1]])-*"])
                    else
                        push!(deps, [pn(i), allvers[i][j], pn(1), "0-*"])
                    end
                end
            end
            sat || push!(problematic, [pn(i), allvers[i][j]])
            continue
        end

        s = shuffle([1:(i-1); (i+1):m])[1:min(d,m-1)]
        for a in s
            push!(deps, [pn(i), allvers[i][j], pn(a), randvspec(k)])
        end
        push!(problematic, [pn(i), allvers[i][j]])
    end

    reqs = [[pn(1), "*"]]
    # reqs = [[pn(1), string(allvers[1][planted1[1]])]]

    # info("SOLUTION: $([(i,planted1[i]) for i = 1:n])")
    # info("REST: $([(i,length(allvers[i])+1) for i = (n+1):m])")

    want = Dict(pn(i) => allvers[i][planted1[i]] for i = 1:n)

    return deps, reqs, want, problematic
end

end # module

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
"""
  AbstractStorage

Module defining interface to use for creating Storage objects
"""
module AbstractStorage

"""
Define an abstract Container Type to be used for storage like objects
"""
abstract type Container end

"""
  Container(::Type{C}, d::Dict)

Create a container object
"""
Container(::Type{C}, d::Dict) where {C<:Container} = nothing

"""
  minimaldict(container::C)

Create a dictionary from `container` containing the minimum
information required to locate the container.
"""
minimaldict(container::C) where {C<:Container} = Dict(n=>getfield(container, n) for n ∈ fieldnames(C))

"""
  session(container::C)

Create a authenticated interaction with a container
"""
session(container::C) where {C<:Container} = nothing

"""
  scrubsession!(container::C)

Remove sensitive information from `session` (e.g. token, client secret)
"""
scrubsession(container::C) where {C<:Container} = copy(container)

"""
  scrubsession!(container::C)

Remove sensitive information from `session` (e.g. token, client secret)
"""
scrubsession!(container::C) where {C<:Container} = container

export Container, minimaldict, scrubsession, scrubsession!, session

end


text:
module krpc
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module schema
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("krpc_pb.jl")
  end
end


text:
#= Component types =#

abstract type PowerConsumer end
abstract type PowerGenerator end

abstract type Actuator <: PowerConsumer end
abstract type Sensor <: PowerConsumer end

abstract type Antenna <: PowerConsumer end
abstract type AttitudeActuator <: Actuator end
abstract type AttitudeSensor <: Sensor end
abstract type OnBoardComputer <: PowerConsumer end
abstract type Payload <: PowerConsumer end

moments(::AttitudeActuator) = error("Not implemented by user")
power_generated(::PowerGenerator) = error("Not implemented by user")
power_consumption(::PowerConsumer) = error("Not implemented by user")
frequency(::Antenna) = error("Not implemented by user")
bandwidth(::Antenna) = error("Not implemented by user")

#= Perturbation types =#

abstract type PerturbationModel

"Aerodynamic model, containing forces and moments."
abstract type AerodynamicModel <: PerturbationModel end
abstract type MagneticModel <: PerturbationModel end
abstract type GravityModel <: PerturbationModel end

forces(::AerodynamicModel) = error("Not implemented by user")
moments(::AerodynamicModel) = error("Not implemented by user")


text:
module EMD

using Dierckx, ElasticArrays
include("utils.jl")
include("bEMD.jl")
include("sEMD.jl")
export sEMD, bEMD

end


text:
using Random

## Leer datos
include("../../utils/preprocesado_datos.jl")
include("../../utils/validation.jl")
include("../../learner/euclidean-1-NN.jl")

file_path= "src/Instancias_APC/"
struct DataFile
    route::String
    class_atributte::String
end

files = [
    DataFile(file_path*"ionosphere.arff", "class"),
    DataFile(file_path*"parkinsons.arff", "Class"),
    DataFile(file_path*"spectf-heart.arff", "Class"),
]
csv_file_path  = "src/resultados/1NN/"

process_name = [
    DataFile(csv_file_path*"1NN-ionosphere.result.csv", "Datos Iosfera"),
    DataFile(csv_file_path*"1NN-parkinsons.result.csv", "Datos parkinson "),
    DataFile(csv_file_path*"1NN-spectf-heart.result.csv", "Datos ataques corazón")
]
Random.seed!(0)
for i in 1:length(files)
    file = files[i]
    data , labels = DataLabelArff(file.route, file.class_atributte)
    ## Entrenar para obtener el peso
    println("==============================================")
    println("$(process_name[i].class_atributte)")
    w = ones(size(data)[2])

    clasificador(data,labels)=LearnerEuclideanOneNN(data, labels), 0, w
    ## Evaluar resultado con cross validation
    VerboseCrossValidation(data, labels, 5, clasificador, process_name[i].route)
end


text:
struct ParticlePositions{N,T,V<:AbstractVector{T},Nm1,Names} <: AbstractGrid{N,T,Names}
    grid::NamedTuple{Names,Tuple{V,Vararg{V,Nm1}}}
    minvals::MVector{N,T}
    maxvals::MVector{N,T}
end

# We need to ensure that we have at least one argument
function ParticlePositions(firstarg::T, args::Vararg{T,Nm1};
    names=:auto, mins=:auto, maxs=:auto) where {Nm1, T<:AbstractVector}

    all_args = (firstarg, args...)
    ParticlePositions(all_args; names, mins, maxs)
end

function ParticlePositions(args::Tuple{T,Vararg{T,Nm1}};
    names=:auto, mins=:auto, maxs=:auto) where {Nm1, T<:AbstractVector}

    N = Nm1 + 1
    names = replace_default_names(names, N)
    grid = NamedTuple{names}(args)

    mins = mins == :auto ? MVector(map(minimum, args)) : mins
    maxs = maxs == :auto ? MVector(map(maximum, args)) : maxs

    ParticlePositions(grid, mins, maxs)
end

function ParticlePositions{N,T}(::UndefInitializer) where {N,T}
    grid = ntuple(N) do i
        T[]
    end
    mins = MVector{N,T}(undef)
    maxs = MVector{N,T}(undef)

    ParticlePositions(grid...; mins, maxs)
end

function Base.empty!(grid::ParticlePositions{N,T}) where {N,T}
    for grid_dir in grid
        empty!(grid_dir)
    end
    mins = minimum(grid)
    maxs = maximum(grid)
    mins .= typemax(T)
    maxs .= typemin(T)

    return grid
end

function Base.empty(::ParticlePositions{N,T}) where {N,T}
    undefp = ParticlePositions{N,T}(undef)
    minimum(undefp) .= typemax(T)
    maximum(undefp) .= typemin(T)

    return undefp
end

function update_minmax!(grid::ParticlePositions, new_grid::ParticlePositions)
    mins = minimum(grid)
    maxs = maximum(grid)
    new_mins = minimum(new_grid)
    new_maxs = maximum(new_grid)
    for (m,nm) in zip(mins, new_mins)
        if nm > m
            m = nm
        end
    end
    for (m,nm) in zip(maxs, new_maxs)
        if nm > m
            m = nm
        end
    end
end

function update_minmax!(grid::ParticlePositions, new_grid)
    mins = minimum(grid)
    maxs = maximum(grid)
    new_mins = minimum.(new_grid)
    new_maxs = maximum.(new_grid)
    for (m,nm) in zip(mins, new_mins)
        if nm > m
            m = nm
        end
    end
    for (m,nm) in zip(maxs, new_maxs)
        if nm > m
            m = nm
        end
    end
end

function Base.append!(grid::ParticlePositions, new_grid)
    for (grid_dir, new_g) in zip(grid, new_grid)
        append!(grid_dir, new_g)
    end
    update_minmax!(grid, new_grid)

    return grid
end

function Base.selectdim(g::ParticlePositions, dir::Symbol, idx)
    grid = getdomain(g)
    # For ParticleGrids we always have the same number of elements on
    # each axis, so we have to slice on all dimensions
    for dir in propertynames(g)
        selected_axis = getproperty(g, dir)
        selection = selectdim(selected_axis, 1, idx)
        grid = setindex!!(grid, selection, dir)
    end

    return parameterless_type(g)(grid..., names=propertynames(grid))
end

# This makes size(field) == size(grid)
Base.size(g::ParticlePositions) = (length(first(g)), )

Base.minimum(g::ParticlePositions) = getfield(g, :minvals)
Base.maximum(g::ParticlePositions) = getfield(g, :maxvals)

Base.sort(g::ParticlePositions, dir) = ParticlePositions(sort(getproperty(g, dir)), g.minvals, g.maxvals)

function Base.sort!(f::T, dim) where T <: AbstractPICDataStructure
    sort!(domain_discretization(T), f, dim)
end

function Base.sort!(::ParticleGrid, f, dir)
    grid = getdomain(f)
    sort_idx = sortperm(getproperty(grid, dim))
    permute!.(grid, (sort_idx,))
    permute!(unwrapdata(f), sort_idx)

    return f
end