text:
function uniformCross(parent1::GAGenotype, parent2::GAGenotype, rng::Random.AbstractRNG)

    parent1Rep = parent1._representation
    parent2Rep = parent2._representation
    genLen = length(parent1Rep)
    genRepType = typeof(parent1Rep[1])
    genType = typeof(parent1)

    child1 = Array{genRepType}(undef, genLen)
    child2 = Array{genRepType}(undef, genLen)

    for i=1:genLen
        prob = rand(rng)
        if prob > 0.5
            child1[i] = parent1Rep[i]
            child2[i] = parent2Rep[i]
        else
            child1[i] = parent2Rep[i]
            child2[i] = parent1Rep[i]
        end
    end
    return [genType(child1), genType(child2)]
end

function kPointCross(parent1::GAGenotype, parent2::GAGenotype, rng::Random.AbstractRNG, k::Integer=1)

    if k < 1
        error("k_points ($k) should be greater than 0")
    end

    parent1Rep = parent1._representation
    parent2Rep = parent2._representation
    genLen = length(parent1Rep)
    if k > genLen
        error("k_points ($k) should be minor than the size of the genoytpe ($genLen)")
    end

    genRepType = typeof(parent1Rep[1])
    genType = typeof(parent1)

    child1 = Array{genRepType}(undef, genLen)
    child2 = Array{genRepType}(undef, genLen)

    points = Array{Integer}(undef, k+2)
    points[1] = 1
    points[end] = genLen

    for i=1:k
        points[i+1] = (rand(rng, UInt16) % genLen)+1
    end

    points = sort(points)

    for i=1:k
        indexes = points[i]:points[i+1]-1
        if isodd(i)
            child1[indexes] = parent1Rep[indexes]
            child2[indexes] = parent2Rep[indexes]
        else
            child1[indexes] = parent2Rep[indexes]
            child2[indexes] = parent1Rep[indexes]
        end
    end

    indexes = points[k+1]:points[k+2]
    if isodd(k)
        child1[indexes] = parent2Rep[indexes]
        child2[indexes] = parent1Rep[indexes]
    else
        child1[indexes] = parent1Rep[indexes]
        child2[indexes] = parent2Rep[indexes]
    end

    return [genType(child1), genType(child2)]
end

function singlePointCross(parent1::GAGenotype, parent2::GAGenotype, rng::Random.AbstractRNG)
    kPointCross(parent1, parent2, rng, 1)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# EXCLUDE FROM TESTING
using VectorEngine

function kernel(addr::Int64)
   @veshow(addr)
   return
end

vefunc = VectorEngine.vefunction(kernel, Tuple{Int64})

vefunc(1234)

synchronize()

ext: jl


text:
"""
## validate

Replacement for validateattributes
Numeric arrays only

All checks go through when there are NaN values

IN
    typeStr  ::  DataType
        e.g. Float64
        matches x being scalar or array of that type
  lb, ub
    can be scalars
"""
function validate(x, typeStr = Float64; sz = [], lb = [], ub = [], nonNan = true)
    if isa(typeStr, DataType)
        if !isa(x, Array{typeStr})
            throw(ErrorException("Not Array of $typeStr"))
        end
    end

    check_bounds(x, lb, ub);

    # Works when NaN present
    # if !isreal(x)
    #   throw(ErrorException("Complex values encountered"))
    # end
    if any(isinf.(x))
    throw(ErrorException("Inf encountered"))
    end

    if nonNan
        if any(isnan.(x))
          throw(ErrorException("NaN encountered"))
        end
    end

    # Size of a scalar is ()
    if length(sz) > 0
        if !isequal(size(x), sz)
            throw(DimensionMismatch("size(x) = $(size(x)),  sz = $sz"))
        end
    end
    return nothing
end

"""
## Validate numerical scalar
"""
function validate_scalar(x, typeStr = Float64; lb = [], ub = [], nonNan = true)
    if isa(typeStr, DataType)
        if !isa(x, typeStr)
            throw(ErrorException("Not a $typeStr"))
        end
    end
    check_bounds(x, lb, ub);
    if any(isinf.(x))
        throw(ErrorException("Inf encountered"))
    end
    if nonNan
        if any(isnan.(x))
            throw(ErrorException("NaN encountered"))
        end
    end
    return nothing
end

"""
## Check bounds
"""
function check_bounds(x, lb, ub)
    # Can check bounds including Nan values
    if !isempty(lb)
      if any(x .< lb)
        throw(ErrorException("Min x: $(minimum(x))"))
      end
    end
    if !isempty(ub)
      if any(x .> ub)
        throw(ErrorException("Max x: $(maximum(x))"))
      end
    end
    return nothing
end

# -----------


text:
"""
    SequenceSpecifics

The `SequenceSpecifics` module is intended only for implementing an API for
other formats to interact and exchange scanner sequence specific information.
It is intended to be robust and sufficient, as opposed to comprehensive
coverage of sequence design elements. All methods are exported and no
structures are defined.
"""
module SequenceSpecifics

export pulse_sequence,
       scanning_sequence,
       sequence_variant,
       scan_options

# TODO Pulse sequence type could probably have specific type parameters but for
# now this is just going to return String
"""
    pulse_sequence(x) -> String

A general description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).                                                                                |
"""
pulse_sequence(x::Any) = ""

"""
    scanning_sequence

Description of the type of data acquired. Corresponds to DICOM Tag 0018, 0020
`Scanning Sequence`.                                                                                |
"""
scanning_sequence(x::Any) = ""

"""
    sequence_variant

Variant of the ScanningSequence. Corresponds to DICOM Tag 0018, 0021 `Sequence
Variant`.                                                                                |
"""
function sequence_variant end

"""
    scan_options

Parameters of ScanningSequence. Corresponds to DICOM Tag 0018, 0022
`Scan Options`.                                                                                |
"""
function scan_options end

"""
    sequence_name

Manufacturer’s designation of the sequence name. Corresponds to DICOM Tag 0018,
0024 `Sequence Name`.                                                                                |
"""
function sequence_name end

"""
    pulse_sequence_details

Information beyond pulse sequence type that identifies the specific pulse
sequence used (i.e. "Standard Siemens Sequence distributed with the VB17
software," "Siemens WIP ### version #.##," or "Sequence written by X using a
version compiled on MM/DD/YYYY"). |
"""
function pulse_sequence_details end

"""
    nonlinear_gradient_correction(x) -> Bool

Stating if the image saved has been corrected for gradient nonlinearities by
the scanner sequence.                                                                                |
"""
function nonlinear_gradient_correction(x) end

# TODO does this mean JSON is in seconds and DICOM is in ms
"""
    echotime(x) -> Millisecond

The echo time (TE) for the acquisition.

This parameter is REQUIRED if corresponding fieldmap data is present or the
data comes from a multi echo sequence.

please note that the DICOM term is in milliseconds not seconds
"""
echotime(x::Any) = Millisecond(1)
dicomtag(::typeof(echotime)) = (0x0018, 0x0081)
jsontag(::typeof(echotime)) = "EchoTime"

"""
    inversiontime

The inversion time (TI) for the acquisition, specified in seconds. Inversion
time is the time after the middle of inverting RF pulse to middle of excitation
pulse to detect the amount of longitudinal magnetization. Corresponds to DICOM
Tag 0018, 0082 `Inversion Time` (please note that the DICOM term is in
milliseconds not seconds).                                                                                |
"""
inversiontime
dicomtag(::typeof(inversiontime)) = (0x0018, 0x0082)
jsontag(::typeof(inversiontime)) = "InversionTime"

"""
    slicetime

The time at which each slice was acquired within each volume (frame) of the
acquisition. Slice timing is not slice order -- rather, it is a list of times
(in JSON format) containing the time (in seconds) of each slice acquisition in
relation to the beginning of volume acquisition. The list goes through the
slices along the slice axis in the slice encoding dimension (see below). Note
that to ensure the proper interpretation of the `SliceTiming` field, it is
important to check if the OPTIONAL `SliceEncodingDirection` exists. In
particular, if `SliceEncodingDirection` is negative, the entries in
`SliceTiming` are defined in reverse order with respect to the slice axis
(i.e., the final entry in the `SliceTiming` list is the time of acquisition of
slice 0). This parameter is REQUIRED for sparse sequences that do not have the
`DelayTime` field set. In addition without this parameter slice time correction
will not be possible.
"""
function slicetiming end

"""
    slicedirection

Possible values: `i`, `j`, `k`, `i-`, `j-`, `k-` (the axis of the NIfTI data
along which slices were acquired, and the direction in which `SliceTiming` is
defined with respect to). `i`, `j`, `k` identifiers correspond to the first,
second and third axis of the data in the NIfTI file. A `-` sign indicates that
the contents of `SliceTiming` are defined in reverse order - that is, the first
entry corresponds to the slice with the largest index, and the final entry
corresponds to slice index zero. When present, the axis defined by
`SliceEncodingDirection` needs to be consistent with the ‘slice_dim’ field in
the NIfTI header. When absent, the entries in `SliceTiming` must be in the
order of increasing slice index as defined by the NIfTI header.                                                                                |
"""
function slicedirection end

"""
    dwelltime -> Nanosecond

Actual dwell time (in seconds) of the receiver per point in the readout
direction, including any oversampling. For Siemens, this corresponds to DICOM
field (0019,1018) (in ns). This value is necessary for the optional readout
distortion correction of anatomicals in the HCP Pipelines. It also usefully
provides a handle on the readout bandwidth, which isn’t captured in the other
metadata tags. Not to be confused with `EffectiveEchoSpacing`, and the frequent
mislabeling of echo spacing (which is spacing in the phase encoding direction)
as "dwell time" (which is spacing in the readout direction).                                                                                |
"""
dwelltime(x::Any) = Nanosecond(1)
dicomtag(::typeof(dwilltime)) = (0x0019, 0x1018)
jsontag(::typeof(dwilltime)) = "DwellTime"

"""
    flipangle

Flip angle for the acquisition
"""
function flipangle end
dicomtag(::typeof(flipangle)) = (0x0018,0x1314)
jsontag(::typeof(flipangle)) = "FlipAngle"

"""
    multiband_acceleration_factor

The multiband factor, for multiband acquisitions.
"""
function multiband_acceleration_factor end
jsontag(::typeof(multiband_acceleration_factor)) = "MultibandAccelerationFactor"

"""
    negative_contrast -> Bool
"""
function negative_contrast end
jsontag(::typeof(negative_contrast)) = "NegativeContrast"

end


text:
# Load data to use for tests
path = dirname(@__FILE__)
fred = CSV.read("$path/../reference/fred_160812.csv")

# Specify vintage and dates
custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_conditional_end, DSGE.quartertodate("2016-Q3")),
    :n_mon_anticipated_shocks => Setting(:n_mon_anticipated_shocks, 6))
m = Model990(testing = true, custom_settings = custom_settings)
m <= Setting(:rate_expectations_source, :ois)

# Test full, semi, and none conditional data
exp_data, exp_cond_data, exp_semicond_data =
    JLD2.jldopen("$path/../reference/load_data_out.jld2", "r") do file
        read(file, "data"), read(file, "cond_data"), read(file, "semi_cond_data")
    end
levels, semi_levels, full_levels =
    JLD2.jldopen("$path/../reference/transform_data_inputs.jld2", "r") do file
        read(file, "none"), read(file, "semi"), read(file, "full")
end
exp_data_rev_transforms =
    JLD2.jldopen("$path/../reference/transform_data_out.jld2", "r") do file
        read(file, "rev_transform")
    end
exp_hp_pop_hist, exp_hp_pop_forecast =
    JLD2.jldopen("$path/../reference/hp_data_out.jld2", "r") do file
        read(file, "exp_hp_pop_hist"), read(file, "exp_hp_pop_forecast")
    end
pop_forecast = DataFrame(CNP16OV = cumsum(fred[!,:CNP16OV][end] .* (ones(8) .* 1.03)))
pop_forecast[!,:date] = DSGE.get_quarter_ends(DSGE.next_quarter(fred[!,:date][end]), fred[!,:date][end] + Year(2))

@testset "Check helper functions for transform_data" begin
    @test !isempty(collect_data_transforms(m))
    @test collect_data_transforms(m, direction = :rev) == exp_data_rev_transforms
    hp_hist, hp_forecast = DSGE.transform_population_data(fred, pop_forecast, :CNP16OV,
                                                     pad_forecast_start = true)
    @test hp_hist[!,:filtered_population_recorded] == exp_hp_pop_hist
    @test hp_forecast[!,:filtered_population_forecast] == exp_hp_pop_forecast
    hp_hist, hp_forecast = DSGE.transform_population_data(fred, pop_forecast, :CNP16OV)
    @test hp_hist[!,:filtered_population_recorded] == exp_hp_pop_hist
    @test hp_forecast[!,:filtered_population_forecast] == exp_hp_pop_forecast

    hist, forecast = DSGE.transform_population_data(fred, pop_forecast, :CNP16OV,
                                    pad_forecast_start = true, use_hpfilter = false)
    @test @test_matrix_approx_eq hist[!,:dlpopulation_recorded] difflog(Float64.(fred[!,:CNP16OV]))
    @test @test_matrix_approx_eq forecast[!,:dlpopulation_forecast] difflog(Float64.(pop_forecast[!,:CNP16OV]))
    hist, forecast = DSGE.transform_population_data(fred, pop_forecast, :CNP16OV,
                                    use_hpfilter = false)
    @test @test_matrix_approx_eq hist[!,:dlpopulation_recorded] difflog(Float64.(fred[!,:CNP16OV]))
    @test @test_matrix_approx_eq forecast[!,:dlpopulation_forecast] difflog(Float64.(pop_forecast[!,:CNP16OV]))
end

@testset "Test transform_data for different conditional data" begin
    @info "The following warnings are expected."
    df = transform_data(m, levels; cond_type = :none, verbose = :none)
    start_date = date_presample_start(m)
    end_date = date_mainsample_end(m)
    df = df[start_date .<= df[!,:date] .<= end_date, :]
    DSGE.missing_cond_vars!(m, df, cond_type = :none, check_empty_columns = false)
    @test @test_matrix_approx_eq df_to_matrix(m,df) exp_data
    df = transform_data(m, semi_levels; cond_type = :semi, verbose = :none)
    end_date = date_conditional_end(m)
    df = df[start_date .<= df[!,:date] .<= end_date, :]
    DSGE.missing_cond_vars!(m, df, cond_type = :semi, check_empty_columns = false)
    @test @test_matrix_approx_eq df_to_matrix(m,df, cond_type = :semi) exp_semicond_data
    df = transform_data(m, full_levels; cond_type = :full, verbose = :none)
    end_date = date_conditional_end(m)
    df = df[start_date .<= df[!,:date] .<= end_date, :]
    DSGE.missing_cond_vars!(m, df, cond_type = :full, check_empty_columns = false)
    @test @test_matrix_approx_eq df_to_matrix(m,df, cond_type = :full) exp_cond_data
end

nothing


text:
ms_default=8
an_default = Plots.font(10,:left,:top,0.0,:black)

###########################################
"""
    plt = plotsgaugelocation(gauge::VisClaw.Gauge; offset=(0,0), font::Plots.Font, annotation_str=@sprintf(" %s",gauge.id), kwargs...)

    plt = plotsgaugelocation(gauges::Vector{VisClaw.Gauge}; offset=(0,0), font::Plots.Font, annotation_str=" ", kwargs...)

    plotsgaugelocation!(plt::Plots.Plot, gauge::VisClaw.Gauge; offset=(0,0), font::Plots.Font, annotation_str=@sprintf(" %s",gauge.id), kwargs...)

    plotsgaugelocation!(plt, gauges::Vector{VisClaw.Gauge}; offset=(0,0), font::Plots.Font, annotation_str=" ", kwargs...)

Function: plot a gauge location (with scatter plot) using Plots
"""
function plotsgaugelocation!(plt, gauge::VisClaw.Gauge; offset=(0,0), label="", font::Plots.Font=an_default, annotation_str=@sprintf(" %s",gauge.id), kwargs...)
    # plot
    plt = Plots.scatter!(plt, [gauge.loc[1]], [gauge.loc[2]]; label=label, ann=(gauge.loc[1]+offset[1], gauge.loc[2]+offset[2], Plots.text(annotation_str, font)), kwargs...)

    # return
    return plt

end
###########################################
"""
$(@doc plotsgaugelocation!)
"""
plotsgaugelocation(gauge::VisClaw.Gauge; offset=(0,0), label="", font::Plots.Font=an_default, annotation_str=@sprintf(" %s",gauge.id), kwargs...) =
plotsgaugelocation!(Plots.plot(), gauge; offset=(0,0), label=label, font=font, annotation_str=annotation_str, kwargs...)
###########################################
function plotsgaugelocation!(plt, gauges::Vector{VisClaw.Gauge}; offset=(0,0), label="", font::Plots.Font=an_default, annotation_str=" ", kwargs...)
    # get values in all gauges
    ngauges = length(gauges)
    loc_all = getfield.(gauges, :loc)
    loc = zeros(Float64, ngauges,2)
    for i=1:ngauges; loc[i,:] = loc_all[i]; end

    if (isa(annotation_str, Vector{String})) || (annotation_str==" ")
        if length(annotation_str) !== ngauges
            id_all = getfield.(gauges, :id)
            # annotation
            annotation_str = map(x -> @sprintf(" %s",x), id_all)
        end
    end
    if isa(annotation_str, String)
        annotation_str = fill(annotation_str, ngauges)
    end

    annotation_arg = Vector{Tuple}(undef,ngauges)
    for i=1:ngauges
        annotation_arg[i] = (loc_all[i][1]+offset[1], loc_all[i][2]+offset[2],
                             Plots.text(annotation_str[i], font))
    end
    # plot
    plt = Plots.scatter!(plt, loc[:,1], loc[:,2]; kwargs..., label=label, ann=annotation_arg)

    # return
    return plt
end
###########################################
plotsgaugelocation(gauges::Vector{VisClaw.Gauge}; offset=(0,0), label="", font::Plots.Font=an_default, annotation_str=" ", kwargs...) =
plotsgaugelocation!(Plots.plot(), gauges::Vector{VisClaw.Gauge}; offset=offset, label=label, font=font, annotation_str=annotation_str, kwargs...)


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

using Markdown
using InteractiveUtils

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

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

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

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	alphas = [0.5, 1.5, 1]
	lam = 2

	lambda(dist::Weibull) = shape(dist)*scale(dist)^(-shape(dist))
	theta(lam,alpha) = (alpha/lam)^(1/alpha)

	dists = [Weibull.(a,theta(lam,a)) for a in alphas]

	hA(dist,x) = pdf(dist,x)/ccdf(dist,x)
	hB(dist,x) = lambda(dist)*x^(shape(dist)-1)
end

# ╔═╡ 884176a0-21c0-11eb-06f8-55964a60d533
begin
	xGrid = 0.01:0.01:10
	hazardsA = [hA.(d,xGrid) for d in dists]
	hazardsB = [hB.(d,xGrid) for d in dists]

	Text("Maximum difference between implementations: $(maximum(maximum.(hazardsA-hazardsB)))")
end

# ╔═╡ 8841a7b0-21c0-11eb-313a-bfa4eb49f351
begin
	Cl = [:blue :red :green]
	Lb = [L"\lambda=" * string(lambda(d)) * ",   " * L"\alpha =" * string(shape(d))
			for d in dists]

	plot(xGrid, hazardsA, c=Cl, label=reshape(Lb, 1,:), xlabel="x",
		ylabel="Instantaneous failure rate", xlims=(0,10), ylims=(0,10))
end

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

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═884176a0-21c0-11eb-06f8-55964a60d533
# ╠═8841a7b0-21c0-11eb-313a-bfa4eb49f351
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Random
using Distributions
using CSV
using DataFrames
using Pipe
using Dates
using DataStructures

function map_day_to_index_Z(day)
    """
    Map the day to the appropriate index for the infection array Z noting the 35 day
    padding at the start.
    """
    return day + 31

end

function map_day_to_index_p(day)
    """
    Map the day to the appropriate index for the actual case data. This just wraps the UD
    mapping for clarity in the code.
    """
    return map_day_to_index_Z(day)

end

function map_day_to_index_UD(day)
    """
    Map the day to the appropriate index for the "detection" arrays U and D. This accounts
    for the fact that the first index corresponds to day 0.
    """
    # a branchless if statement for mapping between day and index
    return (day <= 0) * 1 + (day > 0) * (day + 1)

end

function map_day_to_index_cases(day)
    """
    Map the day to the appropriate index for the actual case data. This just wraps the UD
    mapping for clarity in the code.
    """
    return map_day_to_index_UD(day)

end

function sample_inf_time_delta()
    # these are the calculated cdf from doing
    # p = [pdf(Gamma(2.7, 1.0), x) for x in 1:21]
    # c = cumsum(p / sum(p))
    # we store c as a tuple as this treats it as a constant
    c = (
        0.23123495324177792, 0.5173637000503444, 0.7313700767804615, 0.8616192638327465, 0.9324257531818554, 0.968263917871601, 0.9855305708537918, 0.9935547309702418, 0.9971823528990054, 0.9987870872347043, 0.9994845898851373, 0.9997833904140334, 0.9999098409183975, 0.9999628010773516, 0.9999847842534525, 0.9999938383729761, 0.9999975420011928, 0.999999047825051, 0.9999996567609942, 0.9999999018144428, 1.0
    )
    rn = rand()

    return findfirst(rn <= cᵢ for cᵢ in c)

end

function sample_inf_time_omicron()

    c = (
        0.4068944,  0.66971566, 0.82150389, 0.90524088, 0.95029074,
        0.97414539, 0.98663916, 0.99313045, 0.99648247, 0.99820507, 0.99908686,
        0.99953678, 0.99976574, 0.99988197, 0.99994087, 0.99997066, 0.9999857,
        0.99999329, 0.99999711, 0.99999903, 1.0
    )
    rn = rand()

    return findfirst(rn <= cᵢ for cᵢ in c)

end

function sample_det_time_delta()
    c = (
        0.005481604976502295, 0.058918643630252246, 0.18959485142323212, 0.37100454410185607, 0.5556720050758871, 0.7101611913610963, 0.8230354111406564, 0.8977224959544793, 0.9435387039488634, 0.9700150992568369, 0.98459373725705, 0.9923072112983735, 0.9962539186531401, 0.998216511248062, 0.9991687521425704, 0.9996209873403942, 0.9998317584234413, 0.9999283682832782, 0.9999719976126707, 0.99999143977769, 1.0
    )
    rn = rand()

    return findfirst(rn <= cᵢ for cᵢ in c)

end

function sample_det_time_omicron()
    c = (
        0.05465972865344886, 0.20227289941119453, 0.39199410142203933, 0.5719151467164512, 0.7163612824630651, 0.8207303835305828, 0.890853438104743, 0.9355256297540415, 0.9628506906875296, 0.9790374698285156, 0.9883806383087265, 0.9936593224349681, 0.9965885012866724, 0.9981892154091149, 0.9990524882440441, 0.999512733615265, 0.9997556426337206, 0.9998827037093614, 0.9999486387768106, 0.9999826099893431, 1.0
    )
    rn = rand()

    return findfirst(rn <= cᵢ for cᵢ in c)

end

function sample_inf_time(omicron)
    """
    Sample infection times for num individuals based on the generation
    interval distribution, Gamma(shape_gen, scale_gen).
    """
    omicron == false && return sample_inf_time_delta()
    omicron == true && return sample_inf_time_omicron()

end

function sample_onset_time(omicron)
    """
    Sample incubation times for num individuals based on incubation period
    distribution, Gamma(shape_inc, scale_inc).
    """
    omicron == false && return sample_det_time_delta()
    omicron == true && return sample_det_time_omicron()

end

function sample_times(t, omicron)
    """
    Assuming an initial day t, sample the time until infection and onset for an individual
    and only take the ceiling once we have assigned them. This helps reduce the effect of
    ceiling errors induced by adding ceil numbers.
    """
    infection_time = t + sample_inf_time(omicron)
    onset_time = infection_time + sample_onset_time(omicron)

    return (infection_time, onset_time)

end

function set_simulation_constants(forecast_start_date, forecast_end_date, state)
    """
    Creates structures of the constants and type mappings. . This includes all the dynamical
    constants:
    - k = heterogeneity parameter
    - p_symp = probability of symptoms
    - γ = relative infectiousness of asymptomatic individuals
    - p_symp_given_detect = probability of symptoms given detection
    - p_asymp_given_detect = probability of being asymptomatic given detection
    - consistency_multiplier = chosen such that
        sim_cases < consistency_multiplier * actual cases
        results in cases being injected into the simulation. This is used to account
        for superspreading events after periods of low incidence.
    """
    # simulation_constants = Constants(forecast_start_date, forecast_end_date, state)
    simulation_constants = Constants(forecast_start_date, forecast_end_date)
    individual_type_map = IndividualTypeMap()

    return (simulation_constants, individual_type_map)

end

function NegativeBinomial2(μ, k)
    """
    Function for parameterisation of the negative Binomial in terms of mean and dispersion.
    """
    p = 1 / (1 + μ / k)

    return NegativeBinomial(k, p)

end

function sample_negative_binomial_limit(μ, k; approx_limit = 1000)
    """
    Samples from a NegBin(s, p) distribution. This uses a normal approximation
    when mu is large (i.e. s > approx_limit) to get a 10x runtime improvement.
    """
    X = zero(Int)

    X = rand(NegativeBinomial2(μ, k))
    # mean of NegBin(s, p) => this will boil down to TP
    # if μ <= approx_limit
    #     X = rand(NegativeBinomial2(μ, ϕ))
    # else
    #     σ = sqrt(ϕ)
    #     # use the standard normal transform X = μ + σZ
    #     X = ceil(Int, μ + σ * randn())
    # end

    return X

end

function sample_binomial_limit(n, p; approx_limit = 1000)
    """
    Samples from a Bin(n, p) distribution. This uses a normal approximation
    for np > approx_limit or n(1-p) > approx_limit to acheive a 10x runtime
    improvement.
    """
    X = zero(Int)

    X = rand(Binomial(n, p))
    # if n*p <= approx_limit || n*(1-p) <= approx_limit
    #     X = rand(Binomial(n, p))
    # else
    #     μ = n*p
    #     σ = sqrt(n*p*(1-p))
    #     # use standard normal transform
    #     X = ceil(Int, μ + σ * randn())
    # end

    return X

end

function read_in_TP(
    date,
    state;
    strain = "Delta",
    adjust_TP = false,
)

    # read in the reff file
    TP_file_name = "results/" * date * "/"

    if adjust_TP
        TP_file_name = TP_file_name * "soc_mob_R_adjusted_" * strain * date * ".csv"
    else
        TP_file_name = TP_file_name * "soc_mob_R_" * strain * date * ".csv"
    end

    # drop the first column
    df = CSV.read(TP_file_name, DataFrame, drop = [1])
    # extract the unique states
    unique_states = unique(df[!, "state"])

    TP_dict_local = Dict{String, Array}()
    # get the dates
    TP_dict_local["date"] = @pipe df |>
        filter(:state => ==(state), _) |>
        filter(:type => ==("R_L"), _) |>
        select(_, :date) |>
        Matrix(_)[:]
    # create a vector of the indices
    TP_dict_local["t_indices"] = collect(1:length(TP_dict_local["date"]))

    TP_dict_import = Dict{String, Array}()
    # get the dates
    TP_dict_import["date"] = @pipe df |>
        filter(:state => ==(state), _) |>
        filter(:type => ==("R_L"), _) |>
        select(_, :date) |>
        Matrix(_)[:]
    # create a vector of the indices
    TP_dict_import["t_indices"] = collect(1:length(TP_dict_import["date"]))

    # read in the posterior draws
    samples = CSV.read("results/" * date * "/posterior_sample_" * date * ".csv", DataFrame)
    R_I = samples.R_I
    R_I_omicron = samples.R_I
    voc_effect = zeros(size(R_I))

    if strain == "Delta"
        voc_effect .= samples.voc_effect_delta
    else
        voc_effect .= samples.voc_effect_omicron
    end

    R_I_omicron = R_I_omicron .* voc_effect

    # initialise an import matrix to use for each loop
    TP_import_matrix = zeros(length(TP_dict_import["t_indices"]), 2000)

    omicron_start_date = Dates.Date("2021-11-15")

    for state in unique_states
        # filter by state and RL and keep only TP trajectories
        df_state_matrix = @pipe df |>
            filter(:state => ==(state), _) |>
            filter(:type => ==("R_L"), _) |>
            select(_, Not(1:8)) |>
            Matrix(_)
        TP_dict_local[state] = df_state_matrix

        # fill the pre-Omicron R_I
        TP_import_matrix[TP_dict_import["date"] .< omicron_start_date, :] .= R_I[1:2000]'
        # fill the post-Omicron R_I
        TP_import_matrix[TP_dict_import["date"] .>= omicron_start_date, :] .= R_I_omicron[1:2000]'
        # add to the TP dict
        TP_dict_import[state] = TP_import_matrix
    end

    return (TP_dict_local, TP_dict_import)

end

function get_single_state_TP(
    forecast_start_date,
    date,
    state;
    strain = "Delta",
    adjust_TP = false,
)
    """
    Read in the TP for a given state from csv with filedate, date. Adjust the indices
    according to the forecast_start_date and then create matrices of local and import TP's.
    The indices vector is used to allow for pre-forecast window infection dates.
    """

    # read in the TP
    (TP_dict_local, TP_dict_import) = read_in_TP(
        date,
        state,
        strain = strain,
        adjust_TP = adjust_TP
    )
    # adjust the indices so that the forecast start date is t = 0
    TP_dates = TP_dict_local["date"]
    TP_indices = convert.(Int, Dates.value.(TP_dict_local["date"] - forecast_start_date))
    # make a matrix for a given state
    TP_local = TP_dict_local[state]
    TP_import = TP_dict_import[state]

    return (TP_indices, TP_local, TP_import)

end

function read_in_susceptible_depletion(file_date)
    """
    Read in the posterior drawn susceptible_depletion factors. This will be sorted/sampled
    in the same order as the posterior predictive TP's to ensure we use the appropriate
    posterior draws.
    """

    # read in the reff file
    file_name = "results/" *
        file_date *
        "/" *
        "sampled_susceptible_depletion_" *
        file_date *
        ".csv"

    # susceptible_depletion = Vector(CSV.read(file_name, DataFrame, drop = [1],)[:, 1])
    susceptible_depletion = Matrix(CSV.read(file_name, DataFrame, drop = [1],))

    return susceptible_depletion

end

function read_in_prop_omicron(file_date, state)
    """
    Read in the posterior drawn susceptible_depletion factors. This will be sorted/sampled
    in the same order as the posterior predictive TP's to ensure we use the appropriate
    posterior draws.
    """

    # prefix and suffix for the file name
    file_name_prefix = "results/" *
        file_date *
        "/"

    file_name_suffix = "_" *
        state *
        file_date *
        ".csv"

    # these vars are used to calculate the omicron proportion
    prop_vars = ("m0", "m1", "r", "tau")
    # prop_vars = ("m0", "r", "tau")
    # create vector of vectors to store the various parameter samples
    res = Dict{String, Vector{Float64}}()

    for v in prop_vars
        file_name = file_name_prefix * v * file_name_suffix
        res[v] = Vector(CSV.read(file_name, DataFrame, drop = [1],)[:, 1])
    end

    return res

end

function read_in_cases(
	date,
	rng;
    apply_delay = true,
    use_mean = false,
)

	# read in the reff file
	case_file_name = "data/interim_linelist_" * date * ".csv"
	# drop the first column
	df = CSV.read(case_file_name, DataFrame)
	# indicator for the NA dates
	is_confirmation = df.date_onset .== "NA"
	# get the confirmation dates
	confirm_dates = convert(Vector, df.date_confirmation)
	confirm_dates = Date.(confirm_dates)

	if apply_delay
        # shift them appropriately
        shape_rd = 2.33
        scale_rd = 1.35
        # sample from delay distribtuion
        rd = zeros(length(confirm_dates))
        if use_mean
            rd .= mean(Gamma(shape_rd, scale_rd))
        else
            # use the discretised pmf to sample the reporting delays
            rd = round.(rand(Gamma(shape_rd, scale_rd), length(confirm_dates)))
        end

        confirm_dates = confirm_dates - rd * Dates.Day(1)
    end

	# initialise array for complete_onset_dates
	complete_dates = deepcopy(confirm_dates)
	# fill the array with the most informative date
	complete_dates[.!is_confirmation] = Date.(
		convert(Vector, df.date_onset[.!is_confirmation])
	)
	complete_dates[is_confirmation] = confirm_dates[is_confirmation]

	# other useful fields
	state = df.state
	import_status = convert.(Int, df.import_status .== "imported")

	# create vector of all dates from the first (good) reported data in 2020 to the most
    # recent data in the linelist.
    # Note that in previous weeks, this was an issue as the data being provided
    # had onset dates which were imputed (and shouldn't have been)
	dates_since_start = range(
		Date("2020-03-01"),
		# Dates.Date(date),
		maximum(complete_dates),
		step = Day(1),
	)

	complete_dates_local = complete_dates[import_status .== 0]
	complete_dates_import = complete_dates[import_status .== 1]
	state_local = state[import_status .== 0]
	state_import = state[import_status .== 1]

	# construct df to hold all the linelists
	new_df = DataFrame("date" => dates_since_start)

	# construct df to hold all the linelists
    date_vec = dates_since_start
	local_case_dict = Dict{String, Vector{Int}}()
	import_case_dict = Dict{String, Vector{Int}}()

	# initialise arrays in the dictionary to copy into
	for s in unique(state)
		# filter by state
		local_case_dict[s] = zeros(Int, length(dates_since_start))
		import_case_dict[s] = zeros(Int, length(dates_since_start))
	end

	# initialise arrays to hold daily cases in
	local_cases = zeros(Int, length(dates_since_start))
	import_cases = zeros(Int, length(dates_since_start))

	# loop over states and then sum the number of cases on that day
	for s in unique(state)
		# filter by state
		complete_dates_local_tmp = complete_dates_local[state_local .== s]
		complete_dates_import_tmp = complete_dates_import[state_import .== s]
		# get cases on each day
		local_cases .= sum.(
			@views complete_dates_local_tmp .== dss for dss in dates_since_start
		)
		import_cases .= sum.(
			@views complete_dates_import_tmp .== dss for dss in dates_since_start
		)
		# append to the df with a deepcopy to avoid a 0 d
		local_case_dict[s] .= Int.(local_cases)
		import_case_dict[s] .= Int.(import_cases)
		# reset the case vectors
		local_cases .= 0
		import_cases .= 0
	end

    return (date_vec, local_case_dict, import_case_dict)

end


text:
using BenchmarkTools
using Plots
using Printf
using DelimitedFiles
using LinearAlgebra

function Y1(y₂)
    return y₂
end

function Y2(y₃)
    return y₃
end

function Y3(y₁,y₃,y₄,y₅,cμ,T∞)
    return -y₃*((y₅/(2*(y₄)))-(y₅/(y₄+cμ/T∞)))-y₁*y₃*((y₄+cμ/T∞)/(sqrt(y₄)*(1+cμ/T∞)));
end

function Y4(y₅)
    return y₅
end

function Y5(y₁,y₃,y₄,y₅,cμ,T∞,M∞,Pr,γ)
    return -y₅^2*((0.5/y₄)-(1/(y₄+cμ/T∞)))-Pr*y₁*y₅/sqrt(y₄)*(y₄+cμ/T∞)/(1+cμ/T∞)-(γ-1)*Pr*M∞^2*y₃^2;
end

function RK(N,h,y1,y2,y3,y4,y5,cμ,T∞,Pr,γ,M∞)
	for i=1:N

		k11 = Y1(y2[i]);
		k21 = Y2(y3[i]);
		k31 = Y3(y1[i], y3[i], y4[i], y5[i],cμ,T∞);
		k41 = Y4(y5[i]);
		k51 = Y5(y1[i], y3[i], y4[i], y5[i],cμ,T∞,M∞,Pr,γ);

		k12 = Y1(y2[i]+0.5*h*k21);
		k22 = Y2(y3[i]+0.5*h*k31);
		k32 = Y3(y1[i]+0.5*h*k11, y3[i]+0.5*h*k31, y4[i]+0.5*h*k41, y5[i]+0.5*h*k51,cμ,T∞);
		k42 = Y4(y5[i]+0.5*h*k51);
		k52 = Y5(y1[i]+0.5*h*k11, y3[i]+0.5*h*k31, y4[i]+0.5*h*k41, y5[i]+0.5*h*k51,cμ,T∞,M∞,Pr,γ);

		k13 = Y1(y2[i]+0.5*h*k22);
		k23 = Y2(y3[i]+0.5*h*k32);
		k33 = Y3(y1[i]+0.5*h*k12, y3[i]+0.5*h*k32, y4[i]+0.5*h*k42, y5[i]+0.5*h*k52,cμ,T∞);
		k43 = Y4(y5[i]+0.5*h*k52);
		k53 = Y5(y1[i]+0.5*h*k12, y3[i]+0.5*h*k32, y4[i]+0.5*h*k42, y5[i]+0.5*h*k52,cμ,T∞,M∞,Pr,γ);

		k14 = Y1(y2[i]+h*k23);
		k24 = Y2(y3[i]+h*k33);
		k34 = Y3(y1[i]+h*k13, y3[i]+h*k33, y4[i]+h*k43, y5[i]+h*k53,cμ,T∞);
		k44 = Y4(y5[i]+h*k53);
		k54 = Y5(y1[i]+h*k13, y3[i]+h*k33, y4[i]+h*k43, y5[i]+h*k53,cμ,T∞,M∞,Pr,γ);

		y5[i+1] = y5[i] + (1/6)*(k51 + 2*k52 + 2*k53 + k54)*h;
		y4[i+1] = y4[i] + (1/6)*(k41 + 2*k42 + 2*k43 + k44)*h;
		y3[i+1] = y3[i] + (1/6)*(k31 + 2*k32 + 2*k33 + k34)*h;
		y2[i+1] = y2[i] + (1/6)*(k21 + 2*k22 + 2*k23 + k24)*h;
		y1[i+1] = y1[i] + (1/6)*(k11 + 2*k12 + 2*k13 + k14)*h;
	end
	return y1,y2,y3,y4,y5
end

function selfsimilar(M∞=1.0, T∞=300, Tw=2, ηmax=10, N=101, ϵ=1e-9)
    Δη = ηmax/N
    γ  = 1.4   # Ratio of specific heat
    cμ = 110.4 # Sutherland law coefficient for [Kelvin]
    Pr = 0.72  # Prandtl Number
	adi = 1

    Δ = 1e-10     # Small number for shooting method (Decrease when ϵ is decreased)

    # Initializing the solution vectors
    y1 = zeros(N+1)   # f
    y2 = zeros(N+1)   # f'
    y3 = zeros(N+1)   # f''
    y4 = zeros(N+1)   # ρ(η)
    y5 = zeros(N+1)   # ρ(η)'
    η  = [i*Δη for i=0:N]

	if adi==1
		y1[1] = 0.0
    	y2[1] = 0.0
    	y5[1] = 0.0

    	α₀ = 0.1 	  # Initial Guess
    	β₀ = 2.3118      # Initial Guess
	elseif adi==0
		y1[1] = 0
    	y2[1] = 0
    	y4[1] = Tw

    	α₀ = 0.1 	  # Initial Guess
    	β₀ = 3.0      # Initial Guess
	end

    y2n = 1
    y4n = 1
    for i=1:100000
        if adi==1
    		y1[1] = 0.0
        	y2[1] = 0.0
        	y5[1] = 0.0

            y3[1] = α₀;    # Initial Guess
            y4[1] = β₀;    # Initial Guess
    	elseif adi==0
    		y1[1] = 0
        	y2[1] = 0
        	y4[1] = Tw

            y3[1] = α₀;    # Initial Guess
            y5[1] = β₀;    # Initial Guess
    	end

		# First solution for Newton's iteration
        y1,y2,y3,y4,y5 = RK(N,Δη,y1,y2,y3,y4,y5,cμ,T∞,Pr,γ,M∞)

        # Storing the freestream values for Newton's iteration method
        y₂ₒ = y2[N+1];
        y₄ₒ = y4[N+1];

        # Small number addition for Newton's iteration method
        if adi==1
    		y1[1] = 0.0
        	y2[1] = 0.0
        	y5[1] = 0.0

            y3[1] = α₀+Δ;   # Initial Guess + Small number
            y4[1] = β₀;    # Initial Guess
    	elseif adi==0
    		y1[1] = 0
        	y2[1] = 0
        	y4[1] = Tw

            y3[1] = α₀+Δ;    # Initial Guess + Small number
            y5[1] = β₀;    # Initial Guess
    	end

		# Second solution for Newton's iteration
        y1,y2,y3,y4,y5 = RK(N,Δη,y1,y2,y3,y4,y5,cμ,T∞,Pr,γ,M∞)

        # Storing the freestream values for Newton's iteration method
        y₂ₙ₁ = y2[N+1];
        y₄ₙ₁ = y4[N+1];

        # Small number addition for Newton's iteration method
        # Small number addition for Newton's iteration method
        if adi==1
    		y1[1] = 0.0
        	y2[1] = 0.0
        	y5[1] = 0.0

            y3[1] = α₀;   # Initial Guess
            y4[1] = β₀+Δ;    # Initial Guess + Small number
    	elseif adi==0
    		y1[1] = 0
        	y2[1] = 0
        	y4[1] = Tw

            y3[1] = α₀;    # Initial Guess
            y5[1] = β₀+Δ;    # Initial Guess + Small number
    	end

		# Third solution for Newton's iteration
        y1,y2,y3,y4,y5 = RK(N,Δη,y1,y2,y3,y4,y5,cμ,T∞,Pr,γ,M∞)

        # Storing the freestream values for Newton's iteration method
        y₂ₙ₂ = y2[N+1];
        y₄ₙ₂ = y4[N+1];

        # Calculation of the next initial guess with Newton's iteration method
        p₁₁ = (y₂ₙ₁-y₂ₒ)/Δ;
        p₂₁ = (y₄ₙ₁-y₄ₒ)/Δ;
        p₁₂ = (y₂ₙ₂-y₂ₒ)/Δ;
        p₂₂ = (y₄ₙ₂-y₄ₒ)/Δ;
        r₁ = 1-y₂ₒ;
        r₂ = 1-y₄ₒ;
        Δα = (p₂₂*r₁-p₁₂*r₂)/(p₁₁*p₂₂-p₁₂*p₂₁);
        Δβ = (p₁₁*r₂-p₂₁*r₁)/(p₁₁*p₂₂-p₁₂*p₂₁);
        α₀ = α₀ + Δα;
        β₀ = β₀ + Δβ;

        #@show(norm(y4))
		if (abs(y2[N+1]-1.0)<ϵ) && (abs(y4[N+1]-1.0)<ϵ) && (abs(y2n-norm(y2))<ϵ) && (abs(y4n-norm(y4))<ϵ)
        	break
    	end
        y2n = norm(y2)
        y4n = norm(y4)
    end

	# Copying values for logical names
    U = y2
    T = y4

    # Integration for η --> y transformation
    y = zeros(N+1);
    for i=2:N+1
       y[i] = y[i-1] + y4[i]*(η[i]-η[i-1]);
    end
    y = y*sqrt(2);

    return η,y,U,T,N
end

# ╔═╡ f044ac30-60d8-4ebe-8048-463ba4b8bebd
ηtest, ytest, Utest, Ttest, Ntest = selfsimilar();
result = @benchmark selfsimilar()

plot([Utest,Ttest],[ηtest,ηtest],
        title = "Compressible Blasius Profiles",
        label = ["U" "T"],
        legend = :topleft,
        xlabel = "U,T",
        ylabel = "\\eta",
        linewidth = 2,
        linecolor = :black,
        markershape = :circle,
        markercolor = :auto,
	)


text:
prob = copy(Problems.car_3obs_static)
n,m,N = size(prob)
U0 = deepcopy(controls(prob))
dt = prob.Z[1].dt

# Generate an initial state trajectory
N = prob.N
x0 = prob.x0
xf = prob.xf
x_guess = SVector{N}(range(x0[1], stop=xf[1], length=N))
y_guess = SVector{N}(range(x0[2], stop=xf[2], length=N))
X0 = [SVector(x_guess[k], y_guess[k], x0[3]) for k = 1:N]
Z0 = Traj(X0,U0,dt*ones(N))

# Create infeasible model
inf_model = InfeasibleModel(prob.model)
inf_con = InfeasibleConstraint(inf_model)

Z = infeasible_trajectory(inf_model, Z0)
@test states(Z) ≈ states(Z0)
rollout!(inf_model, Z, x0)
@test states(Z) ≈ states(Z0)

# Test constraints
model = prob.model
x,u0 = rand(model)
ui = @SVector rand(n)
u = [u0; ui]

z0 = KnotPoint(x,u0,dt)
z = KnotPoint(x,u,dt)

function test_con_allocs(con,z)
    allocs = @allocated evaluate(con, z)
    allocs += @allocated jacobian(con, z)
end

con = prob.constraints[1].con
idx_con = IndexedConstraint(n,n+m,con)
@test evaluate(idx_con, z) == evaluate(con, z0)
@test jacobian(idx_con, z) == jacobian(con, z0)
@test test_con_allocs(idx_con, z) == 0

bnd = BoundConstraint(n,m, u_min=[0,-2], u_max=[Inf,2])
idx_bnd = IndexedConstraint(n,n+m,bnd)

@test test_con_allocs(idx_bnd, z) == 0

evaluate(bnd, z0)
evaluate(idx_bnd, z)
@test evaluate(bnd, z0) == evaluate(idx_bnd, z)
@test [jacobian(bnd, z0) zeros(3,3)] == jacobian(idx_bnd, z)

conSet0 = copy(get_constraints(prob))
add_constraint!(conSet0, bnd, 1:N-1, 1)
conSet = change_dimension(conSet0, n, n+m)
evaluate!(conSet, Z)
jacobian!(conSet, Z)
@test size(conSet[1].∇c[1]) == (3,8)
@test size(conSet0[1].∇c[1]) == (3,5)

# Cost functions
cost0 = prob.obj[1]
ix = 1:n
iu = 1:m
idx_cost = IndexedCost(cost0,ix,iu)
@test stage_cost(cost0, x, u0) ≈ stage_cost(idx_cost, x, u)
@test stage_cost(cost0, x) ≈ stage_cost(idx_cost, x)
Qx0, Qu0 = gradient(cost0, x, u0)
Qx, Qu = gradient(idx_cost, x, u)
@test Qx0 == Qx
@test [Qu0; zeros(n)] == Qu

Qxx0, Quu0, Qux0 = hessian(cost0, x, u0)
Qxx, Quu, Qux = hessian(idx_cost, x, u)
@test Qxx0 == Qxx
@test Qux == zeros(n+m,n)
@test Quu == Diagonal([diag(Quu0); zeros(n)])

costfun = change_dimension(cost0, n, n+m)
@test stage_cost(costfun, x, u) == stage_cost(idx_cost, x, u)
@test gradient(costfun, x, u) == gradient(idx_cost, x, u)
@test hessian(costfun, x, u) == hessian(idx_cost, x, u)
@test costfun.Q isa Diagonal{Float64,<:SVector}
@test costfun.R isa Diagonal{Float64,<:SVector}
@test costfun.H isa SMatrix
@test costfun.q isa SVector
@test costfun.r isa SVector

costs = infeasible_objective(prob.obj, 2.3)
@test costs[1].Q == cost0.Q
@test costs[1].R == Diagonal([diag(cost0.R); fill(2.3,n)])
@test costs[1].q == cost0.q
@test costs[1].r == [cost0.r; zeros(n)]
@test sum(costs[1].H) == sum(cost0.H)

# Test the solve
sprob = copy(Problems.car_3obs_static)
sal = AugmentedLagrangianSolver(sprob)
solve!(sal)
max_violation(sal)

# Test Infeasible solve
sprob = copy(Problems.car_3obs_static)
initial_trajectory!(sprob, Z0)
sprob_inf = InfeasibleProblem(sprob, Z0, 0.01/sprob.Z[1].dt)
sal = AugmentedLagrangianSolver(sprob_inf)
Z_init = copy(sprob_inf.Z)

initial_trajectory!(sal, Z_init)
solve!(sal)
max_violation(sal)

initial_trajectory!(sal, Z_init)
@test (@allocated solve!(sal)) == 0


text:
function amelioration(sol::Solution, prob::Problem)::Solution
    chgt::Bool = true

    chgt21::Bool = true
	chgt12::Bool = true
    chgt11::Bool = true
    chgt01::Bool = true

	nbChgt::Int64 = 0

	nbChgt21::Int64 = 0
	nbChgt12::Int64 = 0
	nbChgt11::Int64 = 0
	nbChgt01::Int64 = 0

    #while chgt

		chgt = false
		nbChgt += 1

		nbChgt21 = 0
		nbChgt12 = 0
		nbChgt11 = 0
		nbChgt01 = 0

		@time while chgt21
			sol, chgt21 = exchange21(sol, prob)
			nbChgt21 += 1
		end
		# println("nbChgt21 -> ", nbChgt21)
        @time while chgt12
			sol, chgt12 = exchange12(sol, prob)
			nbChgt12 += 1
			if chgt12
				chgt = true
			end
		end
		# println("nbChgt12 -> ", nbChgt12)
		@time if !prob.equalWeights
			while chgt11
				sol, chgt11 = exchange11(sol, prob)
				nbChgt11 += 1
				if chgt11
					chgt = true
				end
			end
		end
		# println("nbChgt11 -> ", nbChgt11)
		@time if sol.nbVarLibre != 0
			while chgt01
				sol, chgt01 = exchange01(sol, prob)
				nbChgt01 += 1
				if chgt01
					chgt = true
				end
			end
		end
		# println("nbChgt01 -> ", nbChgt01)
	# end
	return sol
end


text:
module test_julia_environment_variables

using Test

julia_bin = joinpath(Sys.BINDIR, Base.julia_exename())

load_path = read(`$julia_bin -e 'println(Base.LOAD_PATH)'`, String)
@test occursin("""["@",""", load_path)

empty_load_path = read(setenv(`$julia_bin -e 'println(Base.LOAD_PATH)'`, Dict("JULIA_LOAD_PATH"=>"")), String)
@test endswith(empty_load_path, "String[]\n")

depot_path = read(`$julia_bin -e 'println(Base.DEPOT_PATH)'`, String)
@test endswith(depot_path, string(repr(Base.DEPOT_PATH), "\n"))

end # module test_julia_environment_variables


text:
#=
hdf5utils.jl
Utilities for reading HDF5 data files for the 2019 fastMRI challenge.
Code by Steven Whitaker
Documentation by Jeff Fessler
2020-01-24 update to HDF5 v0.15.2
=#

export h5_get_keys, h5_get_attributes, h5_get_ismrmrd
export h5_get_ESC, h5_get_RSS, h5_get_kspace

using HDF5: h5open, h5read, attributes, eltype, keys
import HDF5 # readmmap

"""
    keys = h5_get_keys(filename::String)
Get `keys` from file using `keys()`.
Returns an `Array{String}`
"""
function h5_get_keys(filename::String)
    return h5open(filename, "r") do file
        keys(file)
    end
end

"""
    attr = h5_get_attributes(filename::String)
Get attributes `attr` from file.  Returns a `Dict`
"""
function h5_get_attributes(filename::String)
    return h5open(filename, "r") do file
        a = attributes(file)
        attr = Dict{String,Any}()
        for s in keys(a)
            attr[s] = read(a[s])
        end
        attr
    end
end

"""
    hdr = h5_get_ismrmrd(filename::String)
Get ISMRM header data from file.  Returns a `?`
"""
function h5_get_ismrmrd(filename::String)
    return h5read(filename, "ismrmrd_header")
end

"""
    esc = h5_get_ESC(filename::String; T::DataType = ComplexF32)
Return `Array` of ESC (emulated single coil) data from file.

HDF5.jl reads the data differently than Python's h5py.
This is significant because the fastMRI paper
says that the datasets have certain dimensionality,
but the dimensions are permuted in Julia compared to Python,
hence the calls to `permutedims` in the following functions.
"""
function h5_get_ESC(filename::String; T::DataType = ComplexF32)
    data = T.(h5read(filename, "reconstruction_esc"))
    return permutedims(data, ndims(data):-1:1)
end

"""
    rss = h5_get_RSS(filename::String; T::DataType = ComplexF32)
Return `Array` of RSS (root sum of squares) data from file.
"""
function h5_get_RSS(filename::String; T::DataType = ComplexF32)
    data = T.(h5read(filename, "reconstruction_rss"))
    return permutedims(data, ndims(data):-1:1)
end

"""
    kspace = h5_get_kspace(filename::String; T::DataType = ComplexF32))
Return `Array` of kspace data from file.
"""
function h5_get_kspace(filename::String; T::DataType = ComplexF32)
    data = h5open(filename, "r") do file
        T = eltype(file["kspace"])
        HDF5.readmmap(file["kspace"], T)
    end
    return permutedims(T.(data), ndims(data):-1:1)
end

#=
Copied (basically) from
https://github.com/MagneticParticleImaging/MPIFiles.jl/blob/79711bf7af389f9e2dd4b0370e64040e5da1e193/src/Utils.jl#L33

superceded by eltype()
function h5_getcomplextype(dataset)
    T = HDF5.get_jl_type(
            HDF5.Datatype(
                HDF5.h5t_get_member_type(
                    datatype(dataset).id,
                    0,
                )
            )
        )
    return Complex{T}
end
=#


text:
module SQLite

using DataFrames
using DataArrays

export sqlitedb, readdlmsql, query, createtable, droptable

import Base: show, close

include("SQLite_consts.jl")
include("SQLite_api.jl")

type SQLiteDB
	file::String
	handle::Ptr{Void}
	resultset::Any
end
function show(io::IO,db::SQLiteDB)
    if db.handle == C_NULL
        print(io,"Null sqlite connection")
    else
        println(io,"sqlite connection")
        println(io,"-----------------")
        println(io,"File: $(db.file)")
        println(io,"Connection Handle: $(db.handle)")
        if isequal(db.resultset,null_resultset)
			print("Contains resultset? No")
		else
			print("Contains resultset(s)? Yes (access by referencing the resultset field (e.g. db.resultset))")
		end
    end
end

typealias TableInput Union(DataFrame,String)

const null_resultset = DataFrame(0)
const null_SQLiteDB = SQLiteDB("",C_NULL,null_resultset)
sqlitedb = null_SQLiteDB #Create default connection = null
const INTrx = r"^\d+$"
const STRINGrx = r"[^eE0-9\.\-\+]"i
const FLOATrx = r"^[+-]?([0-9]+(\.[0-9]*)?|\.[0-9]+)([eE][+-]?[0-9]+)?$"

#Core Functions
function connect(file::String)
	global sqlitedb
	handle = Array(Ptr{Void},1)
	if @FAILED sqlite3_open(file,handle)
		error("[sqlite]: Error opening $file; $(bytestring(sqlite3_errmsg(conn.handle)))")
	else
		return (sqlitedb = SQLiteDB(file,handle[1],null_resultset))
	end
end
function close(conn::SQLiteDB)
	# if is fine to close when conn.handle is NULL (as stated in sqlite3's document)
	if @FAILED sqlite3_close(conn.handle)
		error("[sqlite]: Error closing $(conn.file); $(bytestring(sqlite3_errmsg(conn.handle)))")
	else
		conn.file = ""
		conn.handle = C_NULL
		conn.resultset = null_resultset
	end
end
function internal_query(conn::SQLiteDB,q::String,finalize::Bool=true,stepped::Bool=true)
	stmt = Array(Ptr{Void},1)
	if @FAILED sqlite3_prepare_v2(conn.handle,utf8(q),stmt,[C_NULL])
        ret = bytestring(sqlite3_errmsg(conn.handle))
        internal_query(conn,"COMMIT")
        internal_query(conn,"PRAGMA synchronous = ON")
		error("[sqlite]: $ret")
	end
	stmt = stmt[1]
    r = 0
    if stepped
	   r = sqlite3_step(stmt)
    end
	if finalize
		sqlite3_finalize(stmt)
		return C_NULL, 0
	else
		return stmt, r
	end
end
function query(q::String,conn::SQLiteDB=sqlitedb)
	conn == null_SQLiteDB && error("[sqlite]: A valid SQLiteDB was not specified (and no valid default SQLiteDB exists)")
	stmt, status = SQLite.internal_query(conn,q,false)
	#get resultset metadata: column count, column types, and column names
	ncols = SQLite.sqlite3_column_count(stmt)
	colnames = Array(UTF8String,ncols)
	resultset = Array(Any,ncols)
	check = 0
	for i = 1:ncols
		colnames[i] = bytestring(SQLite.sqlite3_column_name(stmt,i-1))
		t = SQLite.sqlite3_column_type(stmt,i-1)
		if t == SQLite.SQLITE3_TEXT
			resultset[i] = DataArray(String,0)
			check += 1
		elseif t == SQLite.SQLITE_FLOAT
			resultset[i] = DataArray(Float64,0)
			check += 1
		elseif t == SQLite.SQLITE_INTEGER
			resultset[i] = DataArray(WORD_SIZE == 64 ? Int64 : Int32,0)
			check += 1
		else
			resultset[i] = DataArray(Any,0)
		end
	end
	#retrieve resultset
	while status != SQLITE_DONE
		for i = 1:ncols
			t = SQLite.sqlite3_column_type(stmt,i-1)
			if t == SQLITE3_TEXT
				r = bytestring( sqlite3_column_text(stmt,i-1) )
			elseif t == SQLITE_FLOAT
				r = sqlite3_column_double(stmt,i-1)
			elseif t == SQLITE_INTEGER
				r = WORD_SIZE == 64 ? sqlite3_column_int64(stmt,i-1) : sqlite3_column_int(stmt,i-1)
			else
				r = NA
			end
			push!(resultset[i],r)
		end
		status = sqlite3_step(stmt)
	end
    #this is for columns we couldn't get the type for earlier (NULL in row 1); should be the exception
	if check != ncols
        nrows = length(resultset[1])
		for i = 1:ncols
			if nrows > 0 && isna(resultset[i][1])
                d = resultset[i]
				for j = 2:nrows
                    if !isna(d[j])
                        t = typeof(d[j])
                        da = DataArray(t,nrows)
                        for k = 1:nrows
                            da[k] = d[k]
                        end
                        resultset[i] = da
                        break
                    end
                end
			end
		end
	end
	sqlite3_finalize(stmt)
	return (conn.resultset = DataFrame(resultset,Index(colnames)))
end
function createtable(input::TableInput,conn::SQLiteDB=sqlitedb;name::String="",delim::Char='\0',header::Bool=true,types::Array{DataType,1}=DataType[],infer::Bool=true)
	conn == null_SQLiteDB && error("[sqlite]: A valid SQLiteDB was not specified (and no valid default SQLiteDB exists)")
    #these 2 calls are for performance
    internal_query(conn,"PRAGMA synchronous = OFF")

    if typeof(input) == DataFrame
        r = df2table(input,conn,name)
    else
        r = dlm2table(input,conn,name,delim,header,types,infer)
    end
    internal_query(conn,"PRAGMA synchronous = ON")
    return r
end
function df2table(df::DataFrame,conn::SQLiteDB,name::String)
    #get df name for table name if not provided
    dfname = name
    if dfname == ""
        for sym in names(Main)
            if string(:(df)) == string(sym)
                dfname = string(sym)
            end
        end
    end
    # build column specifications
    ncols = length(df)
    colnames = copy(df.colindex.names)
    for col = 1 : ncols
        t = eltype(df[col])
        if t <: FloatingPoint
            colnames[col] *= " REAL"
        elseif t <: Integer
            colnames[col] *= " INT"
        elseif t <: String
            colnames[col] *= " TEXT"
        end
    end
    colspec = join(colnames, ",")
    # create table
    internal_query(conn,"create table $dfname ($colspec)")
    internal_query(conn,"BEGIN TRANSACTION")
	#prepare insert table with parameters for column values
	params = chop(repeat("?,",ncols))
	stmt, r = internal_query(conn,"insert into $dfname values ($params)",false,false)
	#bind, step, reset loop for inserting values
	for row = 1:nrow(df)
		for col = 1:ncols
            d = df[row,col]
            t = typeof(d)
            if t <: FloatingPoint
                SQLite.sqlite3_bind_double(stmt,col,d)
            elseif t <: Integer
                WORD_SIZE == 64 ? sqlite3_bind_int64(stmt,col,d) : sqlite3_bind_int(stmt,col,d)
            elseif t <: NAtype
                sqlite3_bind_null(stmt,col)
            else
                sqlite3_bind_text(stmt,col,string(d),length(string(d)),C_NULL)
            end
		end
		sqlite3_step(stmt)
		sqlite3_reset(stmt)
	end
	sqlite3_finalize(stmt)
    internal_query(conn,"COMMIT")
	return
end
function droptable(table::String,conn::SQLiteDB=sqlitedb)
	conn == null_SQLiteDB && error("[sqlite]: A valid SQLiteDB was not specified (and no valid default SQLiteDB exists)")
	internal_query(conn,"DROP TABLE $table")
	internal_query(conn,"VACUUM")
	return
end
#read raw file direct to sqlite table
function dlm2table(file::String,conn::SQLiteDB,name::String,delim::Char,header::Bool,types::Array{DataType,1},infer::Bool)
    #determine tablename and delimiter
    tablename = name
    if tablename == ""
        tablename = match(r"\w+(?=\.)",file).match
    end
    delimiter = delim
    if delimiter == '\0'
        delimiter = ismatch(r"csv$", file) ? ',' : ismatch(r"tsv$", file) ? '\t' : ismatch(r"wsv$", file) ? ' ' : error("Unable to determine separator used in $file")
    end
    #get column names/types: colnames, ncols, coltypes
    f = open(file)
    firstrow = split(chomp(readline(f)),delimiter)
    ncols = length(firstrow)
    if header
        colnames = firstrow
    else
        colnames = String["x$i" for i = 1:ncols]
        seekstart(f)
    end
    if infer
	    coltypes = Array(DataType,ncols)
	    check = falses(ncols)
	    for r in eachline(f)
			row = split_quoted(chomp(r),delimiter)
	        for i = 1:ncols
	        	if !check[i]
		        	if row[i] == "" #null/missing value
		        		continue
		        	elseif ismatch(INTrx,row[i]) #match a plain integer first
		        		colnames[i] *= " INT"; check[i] = true
		        	elseif ismatch(STRINGrx,row[i]) #then check if it's stringy
						colnames[i] *= " TEXT"; check[i] = true
		        	elseif ismatch(FLOATrx,row[i]) #if it's not integer or string, check if it's a float
						colnames[i] *= " REAL"; check[i] = true
		        	else #if it's still not a float, just make it a string
						colnames[i] *= " TEXT"; check[i] = true
		        	end
		        end
	    	end
	    	sum(check) == ncols && break
	    end
	    if sum(check) < ncols
	    	for i = 1:ncols
	    		if !coltypes[i]
	    			coltypes[i] = String
	    		end
	    	end
    	end
    	seekstart(f)
    	header && readline(f)
	elseif length(types) > 0
		if eltype(types) <: String
			for i = 1:ncols
				colnames[i] *= " " * types[i]
			end
		else
			for i = 1:ncols
				colnames[i] *= types[i] <: Integer ? " INT" : types[i] <: FloatingPoint ? " REAL" : " TEXT"
			end
		end
	end
	colnames = join(colnames,',')
    internal_query(conn,"create table $tablename ($colnames)")
    internal_query(conn,"BEGIN TRANSACTION")
    #prepare insert table with parameters for column values
    params = chop(repeat("?,",ncols))
    stmt, r = internal_query(conn,"insert into $tablename values ($params)",false,false)
    #bind, step, reset loop for inserting values
    for r in eachline(f)
		row = SQLite.split_quoted(chomp(r),delimiter)
	    for col = 1:ncols
	    	d = row[col]
			SQLite.sqlite3_bind_text(stmt,col,d,length(d),C_NULL)
	    end
    	SQLite.sqlite3_step(stmt)
    	SQLite.sqlite3_reset(stmt)
	end
    sqlite3_finalize(stmt)
    internal_query(conn,"COMMIT")
    close(f)
    return
end
#read raw file to sqlite table (call dlm2table), then run sql statement on table to return df (call to query)
function readdlmsql(input::String,conn::SQLiteDB=sqlitedb;sql::String="select * from file",name::String="file",delim::Char='\0',header::Bool=true,types::Array{DataType,1}=DataType[],infer::Bool=true)
	if conn == null_SQLiteDB
		handle = Array(Ptr{Void},1)
		file = tempname()
		SQLite.sqlite3_open(file,handle)
		conn = SQLite.SQLiteDB(file,handle[1],SQLite.null_resultset)
	end
	createtable(input,conn;name=name,delim=delim,header=header,types=types,infer=infer)
	return query(sql,conn)
end
function search_quoted(s::String, c::Char, i::Integer)
    if isempty(c)
        return 1 <= i <= endof(s)+1 ? i :
               i == endof(s)+2      ? 0 :
               error(BoundsError)
    end
    if i < 1 error(BoundsError) end
    i = nextind(s,i-1)
    while !done(s,i)
        d, j = next(s,i)
        if d == '"'
        	i = j
        	d, j = next(s,i)
        	while d != '"'
				i = j
        		d, j = next(s,i)
        	end
        end
        if d in c
            return i
        end
        i = j
    end
    return 0
end
search_quoted(s::String, c::Char) = search_quoted(s,c,start(s))
function split_quoted(str::String, splitter, limit::Integer, keep_empty::Bool)
    strs = String[]
    i = start(str)
    n = endof(str)
    r = search_quoted(str,splitter,i)
    j, k = first(r), last(r)+1
    while 0 < j <= n && length(strs) != limit-1
        if i < k
            if keep_empty || i < j
                push!(strs, str[i:j-1])
            end
            i = k
        end
        if k <= j; k = nextind(str,j) end
        r = search_quoted(str,splitter,k)
        j, k = first(r), last(r)+1
    end
    if keep_empty || !done(str,i)
        push!(strs, str[i:])
    end
    return strs
end
split_quoted(s::String, spl, n::Integer) = split_quoted(s, spl, n, true)
split_quoted(s::String, spl, keep::Bool) = split_quoted(s, spl, 0, keep)
split_quoted(s::String, spl)             = split_quoted(s, spl, 0, true)
end #sqlite module

function sqldf(q::String)
	handle = Array(Ptr{Void},1)
	file = tempname()
	SQLite.sqlite3_open(file,handle)
	conn = SQLite.SQLiteDB(file,handle[1],SQLite.null_resultset)
	#todo: do we need to change column names at all?
	tables = ref(String)
	#find tablenames
	for i in eachmatch(r"(?<=\bfrom\b\s)\w+",q)
		push!(tables,i.match)
	end
	#find join tables, if any
	for i in eachmatch(r"(?<=\bjoin\b\s)\w+",q)
		push!(tables,i.match)
	end
	#make sure all dfs exist
	check = 0
	for df in tables, sym in names(Main)
		if string(sym) == df
			check += 1
		end
	end
	check != length(tables) && error("[sqlite]: DataFrames specified in query were not found")
	for df in tables
		d = eval(symbol(df))
		SQLite.createtable(d,conn;name=df)
	end
	result = SQLite.query(q,conn)
	for df in tables
		SQLite.droptable(df,conn)
	end
	SQLite.sqlite3_close_v2(conn.handle)
	return result
end


text:
using Manopt, Manifolds, Test
@testset "Stochastic Gradient Descent" begin
    M = Sphere(2)
    # 5 point mean
    p = [0.0, 0.0, 1.0]
    s = 1.0
    pts = [
        exp(M, p, X) for
        X in [zeros(3), [s, 0.0, 0.0], [-s, 0.0, 0.0], [0.0, s, 0.0], [0.0, -s, 0.0]]
    ]
    F(y) = 1 / 2 * sum([distance(M, y, x)^2 for x in pts])
    gradF(M, y) = sum([-log(M, y, x) for x in pts])
    sgradF1(M, y) = [-log(M, y, x) for x in pts]
    function sgradF1!(M, X, y)
        for i in 1:length(pts)
            log!(M, X[i], y, pts[i])
            X[i] .*= -1
        end
        return X
    end
    sgradF2 = [((M, y) -> -log(M, y, x)) for x in pts]
    sgradF2! = [function f!(M, X, y)
        log!(M, X, y, x)
        X .*= -1
        return X
    end for x in pts]

    @testset "Constructors" begin
        p1 = StochasticGradientProblem(M, sgradF1)
        p1e = StochasticGradientProblem(M, sgradF1!; evaluation=MutatingEvaluation())
        p2 = StochasticGradientProblem(M, sgradF2)
        p2m = StochasticGradientProblem(M, sgradF2!; evaluation=MutatingEvaluation())
        @test get_gradient(p1, 1, p) == zeros(3)
        @test get_gradient(p2, 1, p) == zeros(3)
        for pr in [p1, p2, p2m]
            X = zero_vector(M, p)
            get_gradient!(pr, X, 1, p)
            @test X == get_gradient(pr, 1, p)
        end
        @test get_gradients(p1, p) == get_gradients(p2, p)
        @test get_gradients(p2m, p) == get_gradients(p2, p)
        Z = get_gradients(p2m, p)
        Z2 = similar.(Z)
        get_gradients!(p2m, Z2, p)
        @test Z == Z2
        Z3 = similar.(Z)
        get_gradients!(p1, Z3, p)
        @test Z == Z3
        Z4 = similar.(Z)
        get_gradients!(p2, Z4, p)
        @test Z == Z4
        Z5 = similar.(Z)
        get_gradients!(p1e, Z5, p)
        @test Z == Z5
        @test_throws ErrorException get_gradients(p1e, p)
        @test_throws ErrorException get_gradient!(p1e, Z4, 1, p)
        o = StochasticGradientDescentOptions(
            M, deepcopy(p), zero_vector(M, p); direction=StochasticGradient(deepcopy(p))
        )
        o.order = collect(1:5)
        o.order_type = :Linear
        initialize_solver!(p1, o)
        o.order_type = :Linear
        step_solver!(p1, o, 1)
        @test o.x == exp(M, p, get_gradient(p1, 1, p))
    end
    @testset "Momentum and Average Processor Constructors" begin
        p1 = StochasticGradientProblem(M, sgradF1)
        p2 = GradientProblem(M, F, gradF)
        m1 = MomentumGradient(M, p, StochasticGradient(zero_vector(M, p)))
        m2 = MomentumGradient(p1, p)
        @test typeof(m1.direction) == typeof(m2.direction) #both use StochasticGradient
        m3 = MomentumGradient(M, p)
        m4 = MomentumGradient(p2, p)
        @test m3.direction == m4.direction #both use Gradient
        a1 = AverageGradient(M, p, 10, StochasticGradient(zero_vector(M, p)))
        a2 = AverageGradient(p1, p, 10)
        @test typeof(a1.direction) == typeof(a2.direction) #both use StochasticGradient
        a3 = AverageGradient(M, p, 10)
        a4 = AverageGradient(p2, p, 10)
        @test a3.direction == a4.direction #both use Gradient
    end
    @testset "Comparing Stochastic Methods" begin
        x1 = stochastic_gradient_descent(M, sgradF1, p; order_type=:Linear)
        @test norm(x1) ≈ 1
        o1 = stochastic_gradient_descent(
            M, sgradF1, p; order_type=:Linear, return_options=true
        )
        x1a = get_solver_result(o1)
        @test x1 == x1a
        x2 = stochastic_gradient_descent(M, sgradF1, p; order_type=:FixedRandom)
        @test norm(x2) ≈ 1
        x3 = stochastic_gradient_descent(M, sgradF1, p; order_type=:Random)
        @test norm(x3) ≈ 1
        x4 = stochastic_gradient_descent(
            M,
            sgradF1,
            p;
            order_type=:Random,
            direction=AverageGradient(M, p, 10, StochasticGradient(zero_vector(M, p))),
        )
        @test norm(x4) ≈ 1
        x5 = stochastic_gradient_descent(
            M,
            sgradF1,
            p;
            order_type=:Random,
            direction=MomentumGradient(M, p, StochasticGradient(zero_vector(M, p))),
        )
        @test norm(x5) ≈ 1
    end
end


text:
using Documenter, ReadStat

makedocs(
	modules = [ReadStat],
	sitename = "ReadStat.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/ReadStat.jl.git"
)


text:
"""
Trains SeqModel given a config dict.

Returns dict containing results.

Differs from easy_sample.jl in that it uses masks and masked gibbs sampling.
Used for doing cross-validation on speckled hold-out data.
"""
function easy_sample_masked!(
        model::Union{SeqModel,DistributedSeqModel},
        spikes::Vector{Spike},
        masks::Vector{Mask},
        initial_assignments::Vector{Int64},
        config::Dict
    )

    # Save copy of initial assignments.
    _inits = copy(initial_assignments)

    # Draw annealed masked Gibbs samples.
    (
        unmasked_assignments,
        anneal_assignment_hist,
        anneal_train_log_p_hist,
        anneal_test_log_p_hist,
        anneal_latent_event_hist,
        anneal_globals_hist
    ) =
    annealed_masked_gibbs!(
        model,
        spikes,
        masks,
        initial_assignments,
        config[:num_anneals],
        config[:max_temperature],
        config[:num_spike_resamples_per_anneal],
        config[:samples_per_resample],
        config[:split_merge_moves_during_anneal],
        config[:split_merge_window],
        config[:save_every_during_anneal];
        verbose=true
    )

    # Draw regular masked Gibbs samples.
    (
        unmasked_assignments,
        assignment_hist,
        train_log_p_hist,
        test_log_p_hist,
        latent_event_hist,
        globals_hist
    ) =
    masked_gibbs!(
        model,
        spikes,
        masks,
        unmasked_assignments,
        config[:num_spike_resamples],
        config[:samples_per_resample],
        config[:split_merge_moves_after_anneal],
        config[:split_merge_window],
        config[:save_every_after_anneal];
        verbose=true
    )

    # return the results
    return Dict(

        # Initial assignment variables.
        :initial_assignments => _inits,

        # Results during annealing.
        :anneal_assignment_hist => anneal_assignment_hist,
        :anneal_train_log_p_hist => anneal_train_log_p_hist,
        :anneal_test_log_p_hist => anneal_test_log_p_hist,
        :anneal_latent_event_hist => anneal_latent_event_hist,
        :anneal_globals_hist => anneal_globals_hist,

        # Results after annealing.
        :final_assignments => unmasked_assignments,
        :assignment_hist => assignment_hist,
        :latent_event_hist => latent_event_hist,
        :globals_hist => globals_hist,
        :train_log_p_hist => train_log_p_hist,
        :test_log_p_hist => test_log_p_hist
        )
end


text:
function cleanaddress(;wb=missing,ws=missing,rng=missing,teststate=false)
    x = DataFrame(XLSX.readtable(wb,ws)...)
    regAve = r" Avenue"
    regSt = r" Street"
    regCt = r" Court"
    regPkwy = r" Parkway"
    regBlvd = r" Boulevard"
    regSte = r" Suite"
    regDr = r" Drive"
    regLn = r" Lane"
    regHwy = r" Highway"
    regSymbs = r"[,;:. -!@#$%^&()=<>//]" # + and * cannot be filtered,[] array
    function regexreplace(x)
        x = replace.(x,regAve => " Ave")
        x = replace.(x,regSt => " St")
        x = replace.(x,regCt => " Ct")
        x = replace.(x,regPkwy => " Pkwy")
        x = replace.(x,regBlvd => " Blvd")
        x = replace.(x,regSte => " Ste")
        x = replace.(x,regDr => " Dr")
        x = replace.(x,regLn => " Ln")
        x = replace.(x,regHwy => " Hwy")
        x = replace.(x,regSymbs => "")
    end
    x.address = regexreplace(x.address)
    reg2 = r","
    cleanstring(arr) = map(arr) do x
        s = string(x)
        s = Unicode.normalize(x, casefold=true)
    end
    x.address = cleanstring(x.address)
    outcell(; cell = split(rng,":")[1] |> String) = cell #default out location
    if teststate == false
        XLSX.writetable(wb, x, overwrite=true, sheetname=ws,anchor_cell=outcell()) #anchor_cell is the starting cell for pasting data to excel
    elseif teststate == true
        x
    end
end


text:
export splitTwoRotationsTranslations
##
##This function assumes that mAll is [m,rotation_translation_1,rotation_translation_2]
##
function splitTwoRotationsTranslations(mAll::Array{Float64,1},numMoves1::Int64,numMoves2::Int64,choose::Int64)
numElementsMove1 = numMoves1*5;
numElementsMove2 = numMoves2*5;
m_len = length(mAll) - numElementsMove1 - numElementsMove2;

Ind = 0;

J = sparse(1.0I,length(mAll),length(mAll));
if choose==1
	Ind = 1:(m_len + numElementsMove1)
else
	Ind = [1:m_len;(m_len + numElementsMove1+1):length(mAll)];
end
return mAll[Ind],J[Ind,:];
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int, k::Int)::Float64 = round(7(1 - prod(1 .- k ./ (6n ÷ 7 + 1:n))), digits=9)

compute(70, 20)

@benchmark compute(70, 20)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# original version with operations changed to support Julia 1.1,
# + some minor changes to variable naming, data manipulation, etc.

using LinearAlgebra
using MLDatasets
using Random

struct network_v1
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
    zs::Array{Array{Float64,1},1}
    activations::Array{Array{Float64,1},1}
    ∇_b::Array{Array{Float64,1},1}
    ∇_w::Array{Array{Float64,2},1}
    δ_∇_b::Array{Array{Float64,1},1}
    δ_∇_w::Array{Array{Float64,2},1}
    δs::Array{Array{Float64,1},1}
end

σ(z) = 1/(1+exp(-z))
σ_grad(z) = σ(z)*(1-σ(z))

function (net::network_v1)(a)
    for (w, b) in zip(net.weights, net.biases)
        a = σ.(w*a + b)
    end
    return a
end

function network_v1(sizes)
    num_layers = length(sizes)
    sizearr = sizes
    biases = [randn(y) for y in sizes[2:end]]
    weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    zs = [zeros(y) for y in sizes[2:end]]
    activations = [zeros(y) for y in sizes[1:end]]
    ∇_b = [zeros(y) for y in sizes[2:end]]
    ∇_w = [zeros(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    δ_∇_b = [zeros(y) for y in sizes[2:end]]
    δ_∇_w = [zeros(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    δs = [zeros(y) for y in sizes[2:end]]
    network_v1(num_layers, sizearr, biases, weights, zs, activations,∇_b,∇_w,δ_∇_b,δ_∇_w,δs)
end

function update_batch(net::network_v1, batch, η)
    ∇_b = net.∇_b
    ∇_w = net.∇_w

    for i in 1:length(∇_b)
      ∇_b[i] .= 0.0
    end

    for i in 1:length(∇_w)
      ∇_w[i] .= 0.0
    end

    δ_∇_b = net.δ_∇_b
    δ_∇_w = net.δ_∇_w

    for (x, y) in batch
        backprop!(net, x, y)
        for i in 1:length(∇_b)
            ∇_b[i] .+= δ_∇_b[i]
        end
        for i in 1:length(∇_w)
            ∇_w[i] .+= δ_∇_w[i]
        end
    end

	coef = (η/length(batch))
    for i in 1:length(∇_b)
        net.biases[i] .-= coef.*∇_b[i]
    end

    for i in 1:length(∇_w)
        net.weights[i] .-= coef.*∇_w[i]
    end
end

function backprop!(net::network_v1, x, y)
    ∇_b = net.δ_∇_b
    ∇_w = net.δ_∇_w

    len = net.num_layers - 1
    activations = net.activations
    activations[1] .= x
    zs = net.zs
    δs = net.δs

    for i in 1:len
        b, w, z = net.biases[i], net.weights[i], zs[i]
        mul!(z,w,activations[i]) # z = w * inp
        z .+= b
        activations[i+1] .= σ.(z)
    end

    δ = δs[end]
    δ .= (activations[end] .- y) .* σ_grad.(zs[end])
    ∇_b[end] .= δ

    mul!(∇_w[end], δ, activations[end-1]') # ∇_w[end] = δ * activations[end-1]'

    for l in 1:len-1
        z = zs[end-l]
        mul!(δs[end-l], net.weights[end-l+1]', δ) # δs[end-l] = net.weights[end-l+1]' * δ
        δ = δs[end-l]
        δ .*= σ_grad.(z)
        ∇_b[end-l] .= δ
        mul!(∇_w[end-l], δ, activations[end-l-1]') # ∇_w[end-l] = δ * activations[end-1-l]'
    end
    return nothing
end

function SGDtrain(net::network_v1, traindata, epochs, batch_size, η, testdata=nothing)
    n_test = testdata != nothing ? length(testdata) : nothing
    n = length(traindata)

    traindata = shuffle(traindata) # one time shuffle for performance, then only take random batch index

    println("========")
    for j in 1:epochs
        idx = randperm(n ÷ batch_size) .* batch_size

        batches = [traindata[k-batch_size+1 : k] for k in idx]

        @time for batch in batches
            update_batch(net, batch, η)
        end

        if testdata != nothing
            println("Epoch ", j,": ", evaluate(net, testdata), "/", n_test)
        else
            println("Epoch ", j," complete.")
        end
    end
end

function evaluate(net::network_v1, testdata)
    test_results = [(findmax(net(x))[2] - 1, y) for (x, y) in testdata]
    return sum(Int(x == y) for (x, y) in test_results)
end

function loaddata(rng = 1:60000)
    train_x, train_y = MNIST.traindata(Float64, Vector(rng))
    train_x = [train_x[:,:,x][:] for x in 1:size(train_x, 3)]
    train_y = [vectorize(x) for x in train_y]
    traindata = [(x, y) for (x, y) in zip(train_x, train_y)]

    test_x, test_y = MNIST.testdata(Float64)
    test_x = [test_x[:,:,x][:] for x in 1:size(test_x, 3)]
    testdata = [(x, y) for (x, y) in zip(test_x, test_y)]
    return traindata, testdata
end

function vectorize(n)
    ev = zeros(10)
    ev[n+1] = 1
    return ev
end

function main_v1()
    epochs = 10
    batch_size = 10
    η = 1.25
    net = network_v1([784, 30, 10])
    traindata, testdata = loaddata()
    SGDtrain(net, traindata, epochs, batch_size, η, testdata)
	# @profiler SGDtrain(net, traindata, 1, batch_size, η, testdata)
end


text:
function [DG,Param]=InitMethodDG(OrdPoly,Param,Trans)
Param.OrdPoly=OrdPoly;
DG.OrdPoly=OrdPoly;

[DG.w,DG.xw]=GaussLobattoQuad(DG.OrdPoly);
[DG.DW,DG.DS,DG.DV,DG.DVT]=DerivativeMatrixSingle(DG);
[Param.DW,Param.DS,Param.DV,Param.DVT]=DerivativeMatrixSingle(DG);

Param.i1(1).Ind=1:Param.OrdPoly+1;
Param.i1(2).Ind=Param.OrdPoly+1;
Param.i1(3).Ind=1:Param.OrdPoly+1;
Param.i1(4).Ind=1;

Param.i2(1).Ind=1;
Param.i2(2).Ind=1:Param.OrdPoly+1;
Param.i2(3).Ind=Param.OrdPoly+1;
Param.i2(4).Ind=1:Param.OrdPoly+1;

Param.V(1)=1;
Param.V(2)=1;
Param.V(3)=-1;
Param.V(4)=-1;

[Param.w,Param.xw]=GetQuadMeth(Param.OrdPoly+1,'GaussLobattoQuad','Edge');
%[Param.DW,Param.DS,Param.DV,Param.DVT]=DerivativeMatrix(Param);
[Param.DW,Param.DS]=DerivativeMatrixSingle(Param);

Param.MX=Param.w;
Param.MZ=Param.w;

Param.X=zeros(Param.OrdPoly+1,Param.OrdPoly+1,2,Param.Grid.NumFaces);
Param.dXdx=zeros(Param.OrdPoly+1,Param.OrdPoly+1,2,2,Param.Grid.NumFaces);
Param.J=zeros(Param.OrdPoly+1,Param.OrdPoly+1,Param.Grid.NumFaces);

for iF=1:Param.Grid.NumFaces
  [Param.X(:,:,:,iF),Param.J(:,:,iF),Param.dXdx(:,:,:,:,iF)]=...
    JacobiDG1(DG,Param.Grid.Faces(iF),Param.Grid,Trans,Param);
end

Param.N=zeros(Param.OrdPoly+1,2,Param.Grid.NumEdges);
Param.T1=zeros(Param.OrdPoly+1,2,Param.Grid.NumEdges);
Param.VolSurf=zeros(Param.OrdPoly+1,Param.Grid.NumEdges);
for iE=1:Param.Grid.NumEdges
  iF=Param.Grid.Edges(iE).F(1);
  FE=Param.Grid.Edges(iE).FE(1);
  if FE==1
    Param.T1(:,:,iE)=Param.dXdx(:,1,:,1,iF)...
      ./sqrt(Param.dXdx(:,1,1,1,iF).^2+Param.dXdx(:,1,2,1,iF).^2);
    Param.VolSurf(:,iE)=sqrt(Param.dXdx(:,1,1,1,iF).^2 ...
                            +Param.dXdx(:,1,2,1,iF).^2);
  elseif FE==2
    Param.T1(:,:,iE)=Param.dXdx(Param.OrdPoly+1,:,:,2,iF)...
      ./sqrt(Param.dXdx(Param.OrdPoly+1,:,1,2,iF).^2+Param.dXdx(Param.OrdPoly+1,:,2,2,iF).^2);
    Param.VolSurf(:,iE)=sqrt(Param.dXdx(Param.OrdPoly+1,:,1,2,iF).^2 ...
                            +Param.dXdx(Param.OrdPoly+1,:,2,2,iF).^2);

  elseif FE==3
    Param.T1(:,:,iE)=Param.dXdx(:,Param.OrdPoly+1,:,1,iF)...
      ./sqrt(Param.dXdx(:,Param.OrdPoly+1,1,1,iF).^2+Param.dXdx(:,Param.OrdPoly+1,2,1,iF).^2);
    Param.VolSurf(:,iE)=sqrt(Param.dXdx(:,Param.OrdPoly+1,1,1,iF).^2 ...
                            +Param.dXdx(:,Param.OrdPoly+1,2,1,iF).^2);
  elseif FE==4
    Param.T1(:,:,iE)=Param.dXdx(1,:,:,2,iF)...
      ./sqrt(Param.dXdx(1,:,1,2,iF).^2+Param.dXdx(1,:,2,2,iF).^2);
    Param.VolSurf(:,iE)=sqrt(Param.dXdx(1,:,1,2,iF).^2 ...
                            +Param.dXdx(1,:,2,2,iF).^2);
  end
  Param.N(:,1,iE)=Param.T1(:,2,iE);
  Param.N(:,2,iE)=-Param.T1(:,1,iE);
end
end


text:
# deprecated but left here "just in case"

export write_ppm

clamp(f, b, t) = f < b ? b : f > t ? t : f

function write_ppm(io, scanlines)
    h = length(scanlines)
    if h == 0
        return false
    end
    w = length(scanlines[1])

    function write(io::IO, c)
        val(rgb) = round(Int, 255clamp(sqrt(rgb), 0, 1))
        println(io, "$(val(c[1])) $(val(c[2])) $(val(c[3]))")
    end

	println(io, "P3\n$w $h\n255")
	foreach(scanline->foreach(p->write(io, p), scanline), scanlines[end:-1:1])
end


text:
# These functions are based on C functions in the chp.c program by Scott Aaronson
# These are pure Julia implementations, and do not link to the chp.c library
# chp.c is Copyright (c), Scott Aaronson, and is not allowed for commercial use

###########################
# Stablizer Circuit
# Circuits with only Hadamard, C-NOT, R (Phase)
# followed by only one bits measurement
###########################

export CHPCircuit,HadamardUnit,CNOTUnit,PhaseUnit

type CHPCircuit{N} <: AbstractQuCircuit{N}
    gates::Array{GateUnit,1}
end

CHPCircuit(num::Integer,gates::Array{GateUnit,1}) = CHPCircuit{num}(gates)
CHPCircuit(num::Integer) = CHPCircuit{num}(Array(GateUnit,0))

type HadamardUnit<:AbstractGateUnit{1}
    pos::Int
end

function process!{N}(unit::HadamardUnit,input::CHPState{N})
  a = unit.pos
  # swap x_{ia} with z_{ia}
  vec_temp = Array(Int,Int(2*N))
  vec_temp[:] = input.X[:,a]
  input.X[:,a] = input.Z[:,a]
  input.Z[:,a] = vec_temp[:]

  for i = 1:2N
      input.R[i] $= input.X[i,a]&input.Z[i,a]
  end

  return input
end

type CNOTUnit<:AbstractGateUnit{1}
    pos::Int
    ctrl::Int
end

function process!{N}(unit::CNOTUnit,input::CHPState{N})
    a = unit.pos
    b = unit.ctrl

    r = input.R
    x = input.X
    z = input.Z

    for i=1:2N
        r[i] $= x[i,a] & z[i,b] & (x[i,b]$z[i,a]$1)
        x[i,b] $= x[i,a]
        z[i,a] $= z[i,b]
    end
    return input
end

type PhaseUnit<:AbstractGateUnit{1}
    pos::Int
end

function process!{N}(unit::PhaseUnit,input::CHPState{N})
    a = unit.pos
    for i=1:2N
        input.R[i] $= input.X[i,a] & input.Z[i,a]
        input.Z[i,a] $= input.X[i,a]
    end
    return input
end


text:
"""
    sympy_matrix_to_sj(spmat::[ a sympy matrix object ])

convert spmat to a Symata `List` of `Lists`.
"""
function sympy_matrix_to_sj(spmat)
    jlist = pytosj(spmat.tolist())
    for i in 1:length(jlist)
        jlist[i] = pytosj(jlist[i])
    end
    mat = newargs()
    (nx,ny) = size(jlist)
    for i in 1:nx
        push!(mat, MList(jlist[i,:]))
    end
    MList(mat)
end

### Inverse

@mkapprule Inverse :nargs => 1

## FIXME: something changed somewhere.
##  so that sympy_matrix_to_sj(pyres) does not return a Symata expression. Maybe in PyCall ?

@sjdoc Inverse """
    Inverse(m)

compute the matrix inverse of `m`.
"""
@doap function Inverse(expr)
    arg1 = _sjtopy(mx[1])
    mat = sympy.Matrix(arg1)
    pyres = mat.inv()
    sympy_matrix_to_sj(pyres)
end

### Dot

@mkapprule Dot

@sjdoc Dot """
    Dot(a,b)
    a ⋅ b
    ⋅(a,b)

Compute the matrix product or inner product for `a` and `b` matrices or vectors.
The complex conjugate is applied to neither `a` nor `b`.

The infix operator `⋅` can usually be entered by typing `\\cdot` followed by `TAB`
"""
@doap function Dot(u::ListT,v::ListT)
    symlength(u) == 0 || symlength(v) == 0 && return mx
    vv = vectorq(v)
    if vectorq(u) && vv
        # Using MPlus is much faster than mplus. This means mplus should be implemented differently
        return  MPlus( [mmul(x[1], x[2]) for x in zip(margs(u), (margs(v)))]...)
#        return  mplus( [mmul(x[1], x[2]) for x in zip(margs(u), (margs(v)))]...)
    end
    mqu = matrixq(u)
    if mqu && vv && symlength(u[1]) == symlength(v)
       return matmulvec(u,v)
    end
    mqv = matrixq(v)
    if mqu && mqv
        (nu,mu) = matrixdims(u)
        (nv,mv) = matrixdims(v)
        mu != nv && return mx
        return matmulmat(u,v)
    end
    mx
end

"""
    matmulvec(m,v)

Matrix multiply matrix `m` and vector `v`.
"""
function matmulvec(m,v)
    n1 = symlength(m)
    n2 = symlength(m[1])
    a = newargs(n1)
    for i=1:n1
        c = 0
        r = m[i]
        for j=1:n2
            c += mmul(r[j], v[j])
        end
        a[i] = c
    end
    MList(a)
end

function matmulmat(a,b)
    (ma,n) = matrixdims(a)
    (n,mb) = matrixdims(b)
    c = zeromatrix(ma,mb)
    for i in 1:ma
        for j in 1:mb
            acc = 0 #  much faster
            for k in 1:n
                acc += a[i,k]*b[k,j]
            end
            c[i,j] = acc
        end
    end
    c
end

### IdentityMatrix

@mkapprule IdentityMatrix

@sjdoc IdentityMatrix """
    IdentityMatrix(n)

return the `n x n` identity matrix.
"""
@doap function IdentityMatrix(n::Integer)
    r = newargs(n)
    for i=1:n
        c = newargs(n)
        fill!(c,0)
        r[i] = MList(c)
    end
    mr = MList(r)
    for i=1:n
        mr[i,i] = 1
    end
    mr
end

### DiagonalMatrix

### ZeroMatrix

@mkapprule ZeroMatrix

@sjdoc ZeroMatrix """
    ZeroMatrix(n1,n2)

return an `n1 x n2` matrix of zeros.
"""

@doap  ZeroMatrix(n1,n2) = zeromatrix(n1,n2)

function zeromatrix(n2,n1)
    m = newargs(n2)
    for i=1:n2
        r = newargs(n1)
        fill!(r,0)
        m[i] = MList(r)
    end
    MList(m)
end

### Transpose

@mkapprule Transpose

@sjdoc Transpose """
    Transpose(m)

return the transpose of the matrix `m`.
"""

@doap function Transpose(x::Mxpr{:List})
    matrixq(x) || return mx
    a = margs(x)
    n = length(a)
    m = symlength(a[1])
    a0 = newargs(m)
    for i in 1:m
        a1 = newargs(n)
        for j in 1:n
#            a1[j] = margs(margs(x)[j])[i]
             a1[j] = x[j,i]  # is there a more efficient way ?
        end
        a0[i] = MList(a1)
    end
    MList(a0)
end

### Outer

### Tr

# assumes we have a matrix
matrixdims(m::ListT) = (length(m), length(m[1]))

@mkapprule Tr

@sjdoc Tr """
    Tr(m)

Compute the trace of matrix `m`.
"""
@doap function Tr(m::ListT)
    matrixq(m) || return mx
    n = min(matrixdims(m)...)
    a = newargs(n)
    for i=1:n
        a[i] = m[i,i]
    end
    MPlus(a)
end

### Eigenvalues

@mkapprule Eigenvalues :nargs => 1

@sjdoc Eigenvalues """
    Eigenvalues(m)

Compute the eigenvalues of matrix `m`.
"""
@doap function Eigenvalues(expr)
    arg1 = _sjtopy(mx[1])
    mat = sympy.Matrix(arg1)
    pyres = mat.eigenvals()
    eigdict = pyres |> pytosj
    a = newargs()
    for (k,v) in eigdict
        push!(a, MList(k,v))
    end
    MList(a)
end

### Eigenvectors

@mkapprule Eigenvectors :nargs => 1

@sjdoc Eigenvectors """
    Eigenvectors(m)

Compute the eigenvectors of matrix `m`.
"""
@doap function Eigenvectors(expr)
    arg1 = _sjtopy(mx[1])
    mat = sympy.Matrix(arg1)
    pyres = mat.eigenvects()
    args = newargs()
    for x in pyres
        push!(args, sympy_matrix_to_sj(x[3][1]))
    end
    MList(args)
end


text:
function map_constants(obuf)
    function getCustomType(camelName)
        const custom_type = [
            ( r"PacketType*", :UInt8 )
        ]

        match_idx = findfirst(x -> ismatch(x[1], camelName), custom_type)

        if match_idx > 0
            return custom_type[match_idx][2]
        else
            return nothing
        end
    end

    function toCamel(str)
        stream = IOBuffer()

        bumps = [string(bump) for bump in split(str, '_')]

        for bump in bumps
            print(stream, bump[1])
            print(stream, lowercase(bump[2:end]))
        end

        return takebuf_string(stream)
    end

    push!(obuf, "
# Convenience Constants
# with HSA_ prefix removed and
# in CamelCase
    ")

    for exu in obuf
        if isa(exu, Expr) && exu.head == :const && length(exu.args) >= 1
            assign = exu.args[1]

            if isa(assign, Expr) && assign.head == :(=)
                name_sym = assign.args[1]
                name = string(name_sym)

                # is a HSA_* constant? -> add a camelcase version without prefix
                if startswith(name, "HSA_") && isupper(replace(name, '_', ""))
                    shortName = replace(name, "HSA_", "")
                    camelName = toCamel(shortName)
                    camelName_sym = symbol(camelName)

                    customType = getCustomType(camelName)
                    if customType == nothing
                        push!(obuf, :(const $camelName_sym = $name_sym))
                    else
                        push!(obuf, :(const $camelName_sym = $customType($name_sym)))
                    end
                end
            end
        end
    end
end


text:
struct IdSet{T} <: AbstractSet{T}
  dict::IdDict{T,Nothing}
  IdSet{T}() where T = new(IdDict{T,Nothing}())
end

Base.eltype(::IdSet{T}) where T = T

IdSet() = IdSet{Any}()

Base.push!(s::IdSet{T}, x::T) where T = (s.dict[x] = nothing; s)
Base.delete!(s::IdSet{T}, x::T) where T = (delete!(s.dict, x); s)
Base.in(x, s::IdSet) = haskey(s.dict, x)

IdSet{T}(xs) where T = push!(IdSet{T}(), xs...)

IdSet(xs) = IdSet{eltype(xs)}(xs)

Base.collect(s::IdSet) = Base.collect(keys(s.dict))
Base.similar(s::IdSet, T::Type) = IdSet{T}()

@forward IdSet.dict Base.length

Base.iterate(s::IdSet, st...) = iterate(keys(s.dict), st...)


text:
for N in [Float64, Rational{Int}, Float32]
    # constructor from empty set
    E = EmptySet{N}()
    @test SymmetricIntervalHull(E) == E

    # normal constructor
    h = SymmetricIntervalHull(Ball1(N[2, 3], N(4)))

    # constructor attempt from an unbounded set
    @test_throws AssertionError SymmetricIntervalHull(HalfSpace(N[1], N(0)))

    # dimension
    @test dim(h) == 2

    # support vector
    d = N[1, 0]
    @test σ(d, h)[1] == N(6) && σ(-d, h)[1] == N(-6)
    d = N[0, 1]
    @test σ(d, h)[2] == N(7) && σ(-d, h)[2] == N(-7)

    # boundedness
    @test isbounded(h)

    # isempty
    @test !isempty(h)

    # an_element function
    @test an_element(h) == N[0, 0]

    # radius_hyperrectangle function (uncached and cached)
    h = SymmetricIntervalHull(Ball1(N[2, 3], N(4)))
    for i in 1:2
        @test radius_hyperrectangle(h, 1) == 6
        @test radius_hyperrectangle(h) == N[6, 7]
    end

    # high and low
    h = SymmetricIntervalHull(Singleton(N[1, 2]))
    @test high(h) == N[1, 2]
    @test low(h) == N[-1, -2]
end


text:
using GeometricMultigrid: mg!,Vcycle!,norm
using Optim

# count average V-cycles to get reltol 1e-3
itcount(d,s;kw...) = mg!(state(d...;smooth! = s,kw...);reltol=1e-3,abstol=0,smooth! = s, mxiter=32, kw...)
avecount(data,s;kw...) = sum(itcount(d,s;kw...) for d ∈ data)/length(data)

# measure average residual reduction after it V-cycles
function Δresid(st;it=1,kw...)
    r₀ = norm(st.r)
    for _ in 1:it; Vcycle!(st;kw...); end
    return log10(norm(st.r)/r₀)/it
end
loss(data;p=p₀,kw...) = sum(Δresid(state(d...;p,kw...);p,kw...) for d ∈ data)/length(data)

# optimize parameters to minimize loss
fit(data,p₀;smooth! = pseudo!,kw...) = Optim.minimizer(
    optimize(p->loss(data;p,smooth!,kw...),p₀, Newton(),
    Optim.Options(time_limit=120,show_trace=true,f_tol=1e-4); autodiff = :forward))