text:
mutable struct MapPoint
    id::Int64
    """
    Id of the KeyFrame, from which it was created.
    """
    kfid::Int64
    """
    Set of KeyFrame's ids that are visible from this MapPoint.
    """
    observer_keyframes_ids::OrderedSet{Int64}
    """
    Mean descriptor.
    """
    descriptor::BitVector
    keyframes_descriptors::Dict{Int64, BitVector}
    """
    Descriptor distances to each KeyFrame. kfid => hamming_distance.
    """
    descriptor_distances_map::Dict{Int64, Float64}
    """
    Position in world coordinate system.
    """
    position::Point3f
    """
    True if the MapPoint has been initialized.
    """
    is_3d::Bool
    """
    True if the MapPoint is visible in the current Frame.
    """
    is_observed::Bool

    mappoint_lock::ReentrantLock
end

function MapPoint(::Val{:invalid})
    MapPoint(
        -1, -1, OrderedSet{Int64}(), Point3f(0, 0, 0),
        0, false, false, ReentrantLock())
end

function MapPoint(id, kfid, descriptor, is_observed::Bool = true)
    observed_keyframes_ids = OrderedSet{Int64}(kfid)
    keyframes_descriptors = Dict{Int64, BitVector}(kfid => descriptor)
    descriptor_distances_map = Dict{Int64, Float64}(kfid => 0.0)
    position = Point3f(0, 0, 0)
    is_3d = false

    MapPoint(
        id, kfid, observed_keyframes_ids,
        descriptor, keyframes_descriptors, descriptor_distances_map,
        position, is_3d, is_observed, ReentrantLock())
end

@inline is_valid(m::MapPoint)::Bool = m.id != -1

function add_keyframe_observation!(m::MapPoint, kfid)
    lock(m.mappoint_lock) do
        push!(m.observer_keyframes_ids, kfid)
    end
end

function get_observers(m::MapPoint)
    lock(m.mappoint_lock) do
        deepcopy(m.observer_keyframes_ids)
    end
end

function get_observers_number(m::MapPoint)
    lock(m.mappoint_lock) do
        length(m.observer_keyframes_ids)
    end
end

function get_position(m::MapPoint)
    lock(m.mappoint_lock) do
        m.position
    end
end

function set_position!(m::MapPoint, position)
    lock(m.mappoint_lock) do
        m.position = position
        m.is_3d = true
    end
end

function remove_kf_observation!(m::MapPoint, kfid)
    lock(m.mappoint_lock) do
        kfid in keys(m.observer_keyframes_ids) || return
        pop!(m.observer_keyframes_ids, kfid)
        if isempty(m.observer_keyframes_ids)
            m.descriptor |> empty!
            m.keyframes_descriptors |> empty!
            m.descriptor_distances_map |> empty!
            return
        end
        # Set new anchor KeyFrame if removed previous anchor.
        kfid == m.kfid && (m.kfid = m.observer_keyframes_ids[1];)
        # Find the most representative descriptor among observer KeyFrames.
        mindist = length(m.descriptor) * 8.0
        minid = -1

        kfid in keys(m.keyframes_descriptors) || return
        kfid_desc = m.keyframes_descriptors[kfid]
        for (kfd, kfd_desc) in m.keyframes_descriptors
            kfd == kfid && continue
            dist = hamming_distance(kfid_desc, kfd_desc)
            m.descriptor_distances_map[kfd] -= dist

            desc_distance = m.descriptor_distances_map[kfd]
            if desc_distance < mindist
                mindist = desc_distance
                minid = kfd
            end
        end

        pop!(m.keyframes_descriptors, kfid)
        pop!(m.descriptor_distances_map, kfid)
        minid > 0 && (m.descriptor = m.keyframes_descriptors[minid])
    end
end

function add_descriptor!(m::MapPoint, kfid, descriptor::BitVector)
    lock(m.mappoint_lock) do
        kfid in keys(m.keyframes_descriptors) && return

        descriptor_distance = 0.0
        m.keyframes_descriptors[kfid] = descriptor
        m.descriptor_distances_map[kfid] = descriptor_distance
        length(m.keyframes_descriptors) == 1 && (m.descriptor = descriptor; return)
        # Find the most representative descriptor among observer KeyFrames.
        mindist = length(m.descriptor) * 8.0
        minid = -1
        for (kfd, kfd_desc) in m.keyframes_descriptors
            dist = hamming_distance(descriptor, kfd_desc)
            m.descriptor_distances_map[kfd] += dist
            dist < mindist && (mindist = dist; minid = kfd;)
            descriptor_distance += dist
        end
        # Get descriptor with minimal distance.
        descriptor_distance < mindist && (minid = kfid;)
        m.descriptor = m.keyframes_descriptors[minid]
        m.descriptor_distances_map[kfid] = descriptor_distance
    end
end

"""
Check if MapPoint is bad and update it if so.

Bad MapPoint is a 3D point that is observed by < 2 KeyFrames
and not observed by current Frame (e.g. m.observed = false).
If bad, set `is_3d` to `false`.
"""
function is_bad!(m::MapPoint)::Bool
    lock(m.mappoint_lock) do
        length(m.observer_keyframes_ids) < 2 && !m.is_observed && m.is_3d &&
            (m.is_3d = false; return true)
        isempty(m.observer_keyframes_ids) && !m.is_observed &&
            (m.is_3d = false; return true)
        false
    end
end

function mappoint_min_distance(m1::MapPoint, m2::MapPoint)
    min_distance = 1e6
    lock(m1.mappoint_lock) do
        for d1 in values(m1.keyframes_descriptors), d2 in values(m2.keyframes_descriptors)
            distance = hamming_distance(d1, d2)
            distance < min_distance && (min_distance = distance;)
        end
    end
    min_distance
end


text:
function repl_startdebugger_request(conn, params::NamedTuple{(:debugPipename,),Tuple{String}}, crashreporting_pipename)
    hideprompt() do
        debug_pipename = params.debugPipename
        try
            @debug "Trying to connect to debug adapter."
            socket = Sockets.connect(debug_pipename)
            try
                DebugAdapter.startdebug(socket, function (err, bt)
                    if is_disconnected_exception(err)
                        @debug "connection closed"
                    else
                        printstyled(stderr, "Error while running the debugger", color=:red, bold=true)
                        printstyled(stderr, " (consider adding a breakpoint for uncaught exceptions):\n", color=:red)
                        Base.display_error(stderr, err, bt)
                    end
                end)
            finally
                close(socket)
            end
        catch err
            global_err_handler(err, catch_backtrace(), crashreporting_pipename, "Debugger")
        end
    end
end

function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

function debugger_getdebugitems_request(conn, params)
    accessor = params.juliaAccessor
    out = DebugConfigTreeItem[]
    loaded_modules = Base.loaded_modules_array()
    if accessor == "#root"
        # root modules
        for mod in loaded_modules
            push!(out, DebugConfigTreeItem(string(mod), true, string(mod)))
        end
    else
        obj = get_obj_by_accessor(accessor)
        if obj isa Module
            for name in names(obj; all = true)
                isdefined(obj, name) || continue
                strname = string(name)
                startswith(strname, '#') && continue
                this = getfield(obj, name)
                this === obj && continue

                if this isa Base.Callable || this isa Module
                    push!(out, DebugConfigTreeItem(strname, this isa Module, string(accessor, ".", strname)))
                end
            end
        end
    end
    return sort!(out, lt = (x, y) -> x.hasChildren == y.hasChildren ? x.label < y.label : x.hasChildren)
end

function get_obj_by_accessor(accessor, super = nothing)
    parts = split(accessor, '.')
    @assert length(parts) > 0
    top = popfirst!(parts)
    if super === nothing
        # try getting module from loaded_modules_array first and then from Main:
        loaded_modules = Base.loaded_modules_array()
        ind = findfirst(==(top), string.(loaded_modules))
        if ind !== nothing
            root = loaded_modules[ind]
            if length(parts) > 0
                return get_obj_by_accessor(join(parts, '.'), root)
            end
            return root
        else
            return get_obj_by_accessor(accessor, Main)
        end
    else
        if isdefined(super, Symbol(top))
            this = getfield(super, Symbol(top))
            if length(parts) > 0
                if this isa Module
                    return get_obj_by_accessor(join(parts, '.'), this)
                end
            else
                return this
            end
        end
    end
    return nothing
end

macro enter(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/enter", (code = $(string(command)), filename = $(string(__source__.file)))))
end

macro run(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/run", (code = $(string(command)), filename = $(string(__source__.file)))))
end

ext: jl

max_stars_repo_name: novitk/julia-vscode

max_issues_repo_name: novitk/julia-vscode

max_forks_repo_name: novitk/julia-vscode


text:
"""
    Experiment

An experiment has the models to compare, a stopping rule to make decision.
"""
abstract type Experiment end

"""
    ExperimentABN{T,n} <: Experiment

An experiment with stopping rule of type `T` and `n` models. Models must have the same
`ProbabilisticModel` type.
"""
mutable struct ExperimentABN{T,n} <: Experiment
    models::Dict{String,ProbabilisticModel}
    modelnames::Vector{String}
    rule::T
    winner::Union{String,Nothing}

    function ExperimentABN{T,n}(models, rule::T; modelnames=nothing) where {T <: StoppingRule,n}
        length(models) >= 2 || throw(ArgumentError("Number of models needs to be greater than 1."))
        length(models) == n || throw(ArgumentError("Number of models needs to be equal to $n."))

        if modelnames === nothing
            modelnames = ["control"]
            for i in 1:(length(models) - 1)
                push!(modelnames, "variant $i")
            end
        end
        @assert(length(models) == length(modelnames),
                "Number of model names is different from number of models.")
        return new(Dict(zip(modelnames, models)), modelnames, rule, nothing)
    end
end

function ExperimentABN(models, rule::T; modelnames=nothing) where T <: StoppingRule
    return ExperimentABN{T,length(models)}(models, rule, modelnames=modelnames)
end

"""
    ExperimentAB(models, rule; modelnames=nothing) where T <: StoppingRule

An experiment of `ExperimentABN` with two models.
"""
function ExperimentAB(models, rule::T; modelnames=nothing) where T <: StoppingRule
    return ExperimentABN{T,2}(models, rule, modelnames=modelnames)
end

function Base.show(io::IO, experiment::ExperimentABN)
    modelinfo = ""
    for (i, modelname) in enumerate(experiment.modelnames)
        modelinfo = string(modelinfo, "  \"", modelname, "\": ", experiment.models[modelname])
        if i < length(experiment.modelnames)
            modelinfo = string(modelinfo, "\n")
        end
    end
    print(io, typeof(experiment), "\n", modelinfo)
end

function update!(experiment::ExperimentABN,
                 listofstats::Union{Vector{T},Vector{Vector{T}}}) where T <: ModelStatistics
    for (i, stats) in enumerate(listofstats)
        modelname = experiment.modelnames[i]
        update!(experiment, modelname, stats)
    end
    return nothing
end

function update!(experiment::ExperimentABN,
    modelname::String, stats::Union{T,Vector{T}}) where T <: ModelStatistics
    model = experiment.models[modelname]
    update!(model, stats)
    return nothing
end

function defaultparams(experiment::ExperimentABN)
    model1 = [model for (modelname, model) in experiment.models][1]
    return defaultparams(model1)
end

"""
    ExperimentBF{M} <: Experiment

`ExperimentBF` is an experiment using a Bayes Factor between the
null and alternative hypothesis as the stopping rule.

# Constructors

    ExperimentBF(kwargs...)

## Keywords

- `model::M`: Prior of effect size of alternative hypothesis
- `p0::Float64`: Probablity of null hypothesis
- `winner::Union{String, Nothing`: Decision to reject the null hypothesis or not
- `rule::BayesFactorThresh`: Stopping rule using Bayes Factor as the threshold
- `stats`: Statistics for calculating the bayes factor. Default is `nothing`.
- `modelnames`: Names of the hypotheses. Default is `["null", "alternative"]`.

"""
mutable struct ExperimentBF{M<:BayesFactorModel} <: Experiment
    model::M
    p0::Float64
    winner::Union{String, Nothing}
    rule::BayesFactorThresh
    stats::Union{NormalStatistics, TwoNormalStatistics, Nothing}
    modelnames::Vector{String}

    function ExperimentBF(;
        model, rule, p0=0.5, stats=nothing, modelnames=["null", "alternative"])
        return new{typeof(model)}(model, p0, nothing, rule, stats, modelnames)
    end
end

function update!(experiment::ExperimentBF,
    stats::Union{NormalStatistics, TwoNormalStatistics})
    if experiment.stats === nothing
        experiment.stats = stats
    else
        experiment.stats = update!(experiment.stats, stats)
    end
    return nothing
end

function bayesfactor(experiment::ExperimentBF{NormalEffectSize})
    experiment.stats !== nothing || error("Experiment statistics is not initialized.")
    return bayesfactor(experiment.model, experiment.stats)
end

function bayesfactor(experiment::ExperimentBF{StudentTEffectSize})
    experiment.stats !== nothing || error("Experiment statistics is not initialized.")
    tstats = StudentTStatistics(experiment.stats)
    return bayesfactor(experiment.model, tstats)
end

"""
    expectedloss(modelA, modelB; lossfunc, numsamples)

Approximating the expected loss for choosing model A over model B.

    expectedlosses(experiment::ExperimentABN; lossfunc, numsamples)

Approximating the expected loss for all models in the experiment.
For A/B/N experiment, the expected loss for each model is the maximum of the losses
by comparing that model to all the other models.
"""
function expectedloss(modelA::T, modelB::T, parameters;
    lossfunc=upliftloss, numsamples=10_000) where T <: ProbabilisticModel
    sampleA = samplepost(modelA, parameters, numsamples)
    sampleB = samplepost(modelB, parameters, numsamples)
    length(sampleA) == length(sampleB) ||
        throw(ArgumentError("Approximating expected loss requires equal length."))
    losschooseA = expectedloss(sampleA, sampleB, lossfunc, numsamples)
    return losschooseA
end

function expectedloss(modelA::T, modelB::T; lossfunc=upliftloss,
        numsamples=10_000) where T <: ProbabilisticModel
    parameters = defaultparams(modelA)
    return expectedloss(modelA, modelB, parameters, lossfunc=lossfunc, numsamples=numsamples)
end

function expectedloss(sampleA::Vector{T}, sampleB::Vector{T}, lossfunc, numsamples) where T<:Real
    return sum(lossfunc(sampleA, sampleB)) / numsamples
end

function expectedlosses(experiment::ExperimentABN{ExpectedLossThresh,n},
    parameters::Vector{Symbol}; lossfunc=upliftloss, numsamples=10_000) where n
    modelnames = experiment.modelnames
    expectedlosses = Vector{Float64}()
    # TODO: Need a test to make sure the max loss is selected.
    for (refmodelindex, refmodelname) in enumerate(modelnames)
        refmodel = experiment.models[refmodelname]
        maxloss = -Inf
        for compmodelname in modelnames
            if refmodelname == compmodelname
                continue
            end
            compmodel = experiment.models[compmodelname]
            loss = expectedloss(
                refmodel, compmodel, parameters; lossfunc=lossfunc, numsamples=numsamples)
            if loss > maxloss
                maxloss = loss
            end
        end
        push!(expectedlosses, maxloss)
    end
    return (modelnames, expectedlosses)
end

function expectedlosses(experiment::ExperimentABN{ExpectedLossThresh,n};
    lossfunc=upliftloss, numsamples=10_000) where n
    parameters = defaultparams(experiment)
    return expectedlosses(experiment, parameters, lossfunc=lossfunc, numsamples=numsamples)
end

function probbeatall(experiment::ExperimentABN{ProbabilityBeatAllThresh,n},
    parameters::Vector{Symbol}; numsamples=10_000) where n
    modelnames = experiment.modelnames
    models = [experiment.models[modelname] for modelname in modelnames]
    samples = samplepost(models, parameters, numsamples)
    probsbeatall = Vector{Float64}(undef, length(modelnames))
    for refmodelindex = 1:length(modelnames)
        refsample = samples[:, refmodelindex]
        numbeatsall = refsample .- samples |>
            x -> mapslices(x -> x .>= 0, x, dims=1) |>
            x -> sum(x, dims=2) |>
            x -> sum(x .== length(modelnames))
        probbeatsall = numbeatsall / numsamples
        probsbeatall[refmodelindex] = probbeatsall
    end
    return (modelnames, probsbeatall)
end

function probbeatall(experiment::ExperimentABN{ProbabilityBeatAllThresh,n}; numsamples=10_000) where n
    parameters = defaultparams(experiment)
    return probbeatall(experiment, parameters, numsamples=numsamples)
end

# Stopping rule check functions
function checkrule(stoppingrule::ExpectedLossThresh, decisionvalue)
    return decisionvalue < stoppingrule.threshold
end

function checkrule(stoppingrule::ProbabilityBeatAllThresh, decisionvalue)
    return decisionvalue > stoppingrule.threshold
end

"""
    metrics(experiment, parameters, numsamples)

Returns the winner's index and key metrics of the experiment.
"""
function metrics(experiment::ExperimentABN{ExpectedLossThresh,n},
    parameters::Vector{Symbol}; numsamples=10_000) where n
    modelnames, explosses = expectedlosses(experiment, parameters, numsamples=numsamples)
    minindex = argmin(explosses)
    return (minindex, explosses)
end

function metrics(experiment::ExperimentABN{ProbabilityBeatAllThresh,n},
    parameters::Vector{Symbol}; numsamples=10_000) where n
    modelnames, probabilities = probbeatall(experiment, parameters, numsamples=numsamples)
    maxindex = argmax(probabilities)
    return (maxindex, probabilities)
end

function metrics(experiment::ExperimentABN; numsamples=10_000)
    parameters = defaultparams(experiment)
    return metrics(experiment, parameters, numsamples=numsamples)
end

function metrics(experiment::ExperimentBF)
    (experiment.stats !== nothing) || error("The experiment has no data.")
    return bayesfactor(experiment)
end

"""
    decide!(experiment; numsamples=10_000)

Make decision based on an experiment result and its stopping rule.
"""
function decide!(experiment::ExperimentABN, parameters::Vector{Symbol}; numsamples=10_000)
    winnerindex, winnermetric = metrics(experiment, parameters::Vector{Symbol}, numsamples=numsamples)
    experiment.winner = nothing
    if checkrule(experiment.rule, winnermetric[winnerindex])
        experiment.winner = experiment.modelnames[winnerindex]
    end
    return (experiment.winner, winnermetric)
end

function decide!(experiment::ExperimentABN; numsamples=10_000)
    parameters = defaultparams(experiment)
    return decide!(experiment, parameters, numsamples=numsamples)
end

function decide!(experiment::ExperimentBF)
    modelnames = experiment.modelnames
    bayesfactor = metrics(experiment)
    rule = experiment.rule
    experiment.winner = _decide(rule, bayesfactor, modelnames)
    return (experiment.winner, bayesfactor)
end

function _decide(onesidethresh::OneSidedBFThresh, bayesfactor::Real, modelnames)
    threshold = onesidethresh.threshold
    if bayesfactor > threshold
        return modelnames[2]
    else
        return nothing
    end
end

function _decide(twosidethresh::TwoSidedBFThresh, bayesfactor::Real, modelnames)
    threshold = twosidethresh.threshold
    if bayesfactor > threshold
        return modelnames[2]
    elseif bayesfactor < 1/threshold
        return modelnames[1]
    else
        return nothing
    end
end


text:
using VarInfoX
using Test

# Ignore some objects (just enough for pass these particular test; maybe more
# things should be ignored).
ignore_some_rows(rows) = filter(((name,),) -> name ∉ ("Workqueue", "uvhandles"), rows)

@testset "VarInfoX.jl" begin
    a = VarInfoX.varinfo_seq(Base; recursive = true, all = true)
    b = VarInfoX.varinfo_dac_module(Base; recursive = true, all = true)
    c = VarInfoX.varinfo_parallel_names(Base; recursive = true, all = true)
    a = ignore_some_rows(a)
    b = ignore_some_rows(b)
    c = ignore_some_rows(c)
    @test Set(a) == Set(b) == Set(c)
end


text:
"""
```
plot_impulse_response(m, shock, var, class, input_type, cond_type;
    title = "", kwargs...)

plot_impulse_response(m, shock, vars, class, input_type, cond_type;
    forecast_string = "", plotroot = figurespath(m, \"forecast\"),
    titles = [], kwargs...)
```

Plot the responses of `var` to `shock`. By default, only 90% bands are plotted.

### Inputs

- `m::AbstractDSGEModel`
- `shock::Symbol`: e.g. `:g_sh`
- `var::Symbol` or `vars::Vector{Symbol}`: response variable(s), e.g. `:obs_gdp`
- `class::Symbol`
- `input_type::Symbol`
- `cond_type::Symbol`

### Keyword Arguments

- `forecast_string::String`
- `plotroot::String`: if nonempty, plots will be saved in that directory
- `title::String` or `titles::Vector{String}`
- `verbose::Symbol`

See `?irf` for additional keyword arguments, all of which can be passed
into `plot_history_and_forecast`.

### Output

- `p::Plot` or `plots::OrderedDict{Symbol, Plot}`
"""
function plot_impulse_response(m::AbstractDSGEModel, shock::Symbol, var::Symbol, class::Symbol,
                               input_type::Symbol, cond_type::Symbol;
                               title::String = "",
                               kwargs...)

    plots = plot_impulse_response(m, shock, [var], class, input_type, cond_type;
                                  titles = isempty(title) ? String[] : [title],
                                  kwargs...)
    return plots[var]
end

function plot_impulse_response(m::AbstractDSGEModel, shock::Symbol, vars::Vector{Symbol}, class::Symbol,
                               input_type::Symbol, cond_type::Symbol;
                               input_type2::Symbol = Symbol(),
                               forecast_string::String = "",
                               plotroot::String = figurespath(m, "forecast"),
                               titles::Vector{String} = String[],
                               verbose::Symbol = :low,
                               kwargs...)
    # Read in MeansBands
    mb = read_mb(m, input_type, cond_type, Symbol(:irf, class), forecast_string = forecast_string)
    if input_type2 != Symbol()
        mb2 = read_mb(m, input_type2, cond_type, Symbol(:irf, class), forecast_string = forecast_string)
    else
        mb2 = MeansBands()
    end
    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(m, var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title) in zip(vars, titles)
        # Call recipe
        plots[var] = irf(shock, var, mb, mb2; title = title, input_type = input_type, input_type2 = input_type2, kwargs...)

        # Save plot
        if !isempty(plotroot)
            output_file = get_forecast_filename(plotroot, filestring_base(m), input_type, cond_type,
                                                Symbol("irf_", detexify(shock), "_", detexify(var)),
                                                forecast_string = forecast_string,
                                                fileformat = plot_extension())
            save_plot(plots[var], output_file, verbose = verbose)
        end
    end
    return plots
end

@userplot Irf

"""
```
irf(shock, var, mb; flip_sign = false, label_mean_bands = false,
    mean_color = :black, bands_color = :blue, bands_pcts = [\"90.0%\"])
```

User recipe called by `plot_impulse_response`.

### Inputs

- `shock::Symbol`: e.g. `:g_sh`
- `var::Symbol`: e.g. `:obs_gdp`
- `mb::MeansBands`

### Keyword Arguments

- `flip_sign::Bool`: whether to flip the sign of the impulse response while
  plotting
- `label_mean_bands::Bool`
- `mean_color`
- `bands_color`
- `bands_pcts::Vector{String}`

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
irf
@recipe function f(irf::Irf;
                   flip_sign = false,
                   label_mean_bands = false,
                   mean_color = :black,
                   bands_color = :blue,
                   bands_pcts = which_density_bands(irf.args[3], uniquify = true),
                   input_type::Symbol = Symbol(),
                   input_type2::Symbol = Symbol())
    # Error checking
  #=  if length(irf.args) != 3 || typeof(irf.args[1]) != Symbol || typeof(irf.args[2]) != Symbol ||
        typeof(irf.args[3]) != MeansBands

        error("irf must be given two Symbols and a MeansBands. Got $(typeof(irf.args))")
    end=#
    shock, var, mb, mb2 = irf.args

    varshock = Symbol(var, "__", shock)
    sign = flip_sign ? -1 : 1

    quarters_ahead = collect(1:size(mb.means,1))
    # Bands
   for pct in bands_pcts
        @series begin
            fillcolor := bands_color
            fillalpha := 0.1
            linealpha := 0
            label     := label_mean_bands ? "$pct Bands" : ""
            fillrange := sign * mb.bands[varshock][!,Symbol(pct, " UB")]
            quarters_ahead, sign * mb.bands[varshock][!,Symbol(pct, " LB")]
        end
    end

    # Mean
    @series begin
        label     := label_mean_bands ? "Mean"*string(input_type) : ""
        linewidth := 2
        linecolor := mean_color
        quarters_ahead, sign * mb.means[!, varshock]
    end

   #= if input_type2 != Symbol()
        @series begin
            label     := label_mean_bands ? "Mean"*string(input_type2) : ""
            linewidth := 2
            linecolor := :blue
            @show typeof(sign * mb2.means[varshock])
            sign * mb2.means[varshock]
        end
    end =#
end


text:
using Test
using ExcelReaders
using DataFrames
using Mimi
using MimiDICE2016
using CSVFiles

using MimiDICE2016: getparams

@testset "MimiDICE2016" begin

#------------------------------------------------------------------------------
#   1. Run tests on the whole model
#------------------------------------------------------------------------------

@testset "MimiDICE2016-model" begin

m = MimiDICE2016.get_model();
run(m)

f = openxl(joinpath(@__DIR__, "../data/DICE2016R-090916ap-v2.xlsm"))

#Test Precision
Precision = 1.0e-10

#Time Periods
T=100

#TATM Test (temperature increase)
True_TATM = getparams(f, "B99:CW99", :all, "Base", T);
@test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol = Precision

#MAT Test (carbon concentration atmosphere)
True_MAT = getparams(f, "B87:CW87", :all, "Base", T);
@test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol = Precision

#DAMFRAC Test (damages fraction)
True_DAMFRAC = getparams(f, "B105:CW105", :all, "Base", T);
@test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol = Precision

#DAMAGES Test (damages $)
True_DAMAGES = getparams(f, "B106:CW106", :all, "Base", T);
@test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol = Precision

#E Test (emissions)
True_E = getparams(f, "B112:CW112", :all, "Base", T);
@test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol = Precision

#YGROSS Test (gross output)
True_YGROSS = getparams(f, "B104:CW104", :all, "Base", T);
@test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol = Precision

#AL test (total factor productivity)
True_AL = getparams(f, "B21:CW21", :all, "Base", T);
@test maximum(abs, m[:totalfactorproductivity, :AL] .- True_AL) ≈ 0. atol = Precision

#CPC Test (per capita consumption)
True_CPC = getparams(f, "B126:CW126", :all, "Base", T);
@test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol = Precision

#FORCOTH Test (exogenous forcing)
True_FORCOTH = getparams(f, "B73:CW73", :all, "Base", T);
@test maximum(abs, m[:radiativeforcing, :FORCOTH] .- True_FORCOTH) ≈ 0. atol = Precision

#FORC Test (radiative forcing)
True_FORC = getparams(f, "B100:CW100", :all, "Base", T);
@test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol = Precision

#Utility Test
True_UTILITY = getparams(f, "B129:B129", :single, "Base", T);
@test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol = Precision

end #MimiDICE2016-model testset

#------------------------------------------------------------------------------
#   2. Run tests on SCC
#------------------------------------------------------------------------------

@testset "Standard API" begin

m = MimiDICE2016.get_model()
run(m)

# Test the errors
@test_throws ErrorException MimiDICE2016.compute_scc()  # test that it errors if you don't specify a year
@test_throws ErrorException MimiDICE2016.compute_scc(year=2021)  # test that it errors if the year isn't in the time index
@test_throws ErrorException MimiDICE2016.compute_scc(last_year=2299)  # test that it errors if the last_year isn't in the time index
@test_throws ErrorException MimiDICE2016.compute_scc(year=2105, last_year=2100)  # test that it errors if the year is after last_year

# Test the SCC
scc1 = MimiDICE2016.compute_scc(year=2020)
@test scc1 isa Float64

# Test that it's smaller with a shorter horizon
scc2 = MimiDICE2016.compute_scc(year=2020, last_year=2200)
@test scc2 < scc1

# Test that it's smaller with a larger prtp
scc3 = MimiDICE2016.compute_scc(year=2020, last_year=2200, prtp=0.02)
@test scc3 < scc2

# Test with a modified model
m = MimiDICE2016.get_model()
update_param!(m, :t2xco2, 5)
scc4 = MimiDICE2016.compute_scc(m, year=2020)
@test scc4 > scc1   # Test that a higher value of climate sensitivty makes the SCC bigger

# Test compute_scc_mm
result = MimiDICE2016.compute_scc_mm(year=2030)
@test result.scc isa Float64
@test result.mm isa Mimi.MarginalModel
marginal_temp = result.mm[:climatedynamics, :TATM]
@test all(marginal_temp[1:findfirst(isequal(2030), MimiDICE2016.model_years)] .== 0.)
@test all(marginal_temp[findfirst(isequal(2035), MimiDICE2016.model_years):end] .!= 0.)

end

end #MimiDICE2016 testset

nothing


text:
module DDF

using Reexport

# The order of these include statements matters
include("Defs.jl")
include("Algorithms.jl")
include("ZeroOrOne.jl")
include("SparseOps.jl")
include("Meshing.jl")
include("Manifolds.jl")
include("ManifoldConstructors.jl")
include("Funs.jl")
include("Ops.jl")
include("ManifoldOps.jl")
include("Continuum.jl")
include("Plotting.jl")

@reexport using .Algorithms
@reexport using .Continuum
@reexport using .Defs
@reexport using .Funs
@reexport using .ManifoldConstructors
@reexport using .ManifoldOps
@reexport using .Manifolds
@reexport using .Meshing
@reexport using .Ops
@reexport using .Plotting
@reexport using .SparseOps
@reexport using .ZeroOrOne

end


text:
using RegressionAndOtherStories

hibbs = CSV.read(ros_data("ElectionsEconomy", "hibbs.dat"),
    DataFrame; delim=" ")
f = open(ros_datadir("ElectionsEconomy", "hibbs.csv"), "w")
CSV.write(f, hibbs)
close(f)

hibbs = CSV.read(ros_datadir("ElectionsEconomy", "hibbs.csv"),
    DataFrame)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

# %%
# This notebook is a semi top-down explanation. This cell needs to be
# executed first so that the operators and helper functions are defined
# All of this is explained in the later half of the notebook

using Compose, Interact
Compose.set_default_graphic_size(2inch, 2inch)

points_f = [
    (.1, .1),
    (.9, .1),
    (.9, .2),
    (.2, .2),
    (.2, .4),
    (.6, .4),
    (.6, .5),
    (.2, .5),
    (.2, .9),
    (.1, .9),
    (.1, .1)
]

f = compose(context(), stroke("black"), line(points_f))

rot(pic) = compose(context(rotation=Rotation(-deg2rad(90))), pic)
flip(pic) = compose(context(mirror=Mirror(deg2rad(90), 0.5w, 0.5h)), pic)
above(m, n, p, q) =
    compose(context(),
            (context(0, 0, 1, m/(m+n)), p),
            (context(0, m/(m+n), 1, n/(m+n)), q))

above(p, q) = above(1, 1, p, q)

beside(m, n, p, q) =
    compose(context(),
            (context(0, 0, m/(m+n), 1), p),
            (context(m/(m+n), 0, n/(m+n), 1), q))

beside(p, q) = beside(1, 1, p, q)

over(p, q) = compose(context(),
                (context(), p), (context(), q))

rot45(pic) =
    compose(context(0, 0, 1/sqrt(2), 1/sqrt(2),
        rotation=Rotation(-deg2rad(45), 0w, 0h)), pic)

# Utility function to zoom out and look at the context
zoomout(pic) = compose(context(),
                (context(0.2, 0.2, 0.6, 0.6), pic),
                (context(0.2, 0.2, 0.6, 0.6), fill(nothing), stroke("black"), strokedash([0.5mm, 0.5mm]),
                    polygon([(0, 0), (1, 0), (1, 1), (0, 1)])))

function read_path(p_str)
    tokens = [try parsefloat(x) catch symbol(x) end for x in split(p_str, r"[\s,]+")]
    path(tokens)
end

fish = compose(context(units=UnitBox(260, 260)), stroke("black"),
            read_path(strip(readall("fish.path"))))

rotatable(pic) = @manipulate for θ=0:0.001:2π
    compose(context(rotation=Rotation(θ)), pic)
end

blank = compose(context())

fliprot45(pic) = rot45(compose(context(mirror=Mirror(deg2rad(-45))),pic))

# Hide this cell.
display(MIME("text/html"), """<script>
var cell = \$(".container .cell").eq(0), ia = cell.find(".input_area")
if (cell.find(".toggle-button").length == 0) {
ia.after(
    \$('<button class="toggle-button">Toggle hidden code</button>').click(
        function (){ ia.toggle() }
        )
    )
ia.hide()
}
</script>""")

# %% [markdown]
# # Functional Geometry
# *Functional Geometry* is a paper by Peter Henderson ([original (1982)](users.ecs.soton.ac.uk/peter/funcgeo.pdf), [revisited (2002)](https://cs.au.dk/~hosc/local/HOSC-15-4-pp349-365.pdf)) which deconstructs the MC Escher woodcut *Square Limit*
#
# ![Square Limit](http://i.imgur.com/LjRzmNM.png)

# %% [markdown]
# > A picture is an example of a complex object that can be described in terms of its parts.
# Yet a picture needs to be rendered on a printer or a screen by a device that expects to
# be given a sequence of commands. Programming that sequence of commands directly is
# much harder than having an application generate the commands automatically from the
# simpler, denotational description.

# %% [markdown]
# A `picture` is a *denotation* of something to draw.
#
# e.g. The value of f here denotes the picture of the letter F

# %% [markdown]
# Original at http://nbviewer.jupyter.org/github/shashi/ijulia-notebooks/blob/master/funcgeo/Functional%20Geometry.ipynb

# %% [markdown]
# ## In conclusion
#
# We described Escher's *Square Limit* from the description of its smaller parts, which in turn were described in terms of their smaller parts.
#
# This seemed simple because we chose to talk in terms of an *algebra* to describe pictures. The primitives `rot`, `flip`, `fliprot45`, `above`, `beside` and `over` fit the job perfectly.
#
# We were able to describe these primitives in terms of `compose` `contexts`, which the Compose library knows how to render.
#
# Denotation can be an easy way to describe a system as well as a practical implementation method.
#
# [Abstraction barriers](https://mitpress.mit.edu/sicp/full-text/sicp/book/node29.html) are useful tools that can reduce the cognitive overhead on the programmer. It entails creating layers consisting of functions which only use functions in the same layer or layers below in their own implementation. The layers in our language were:
#
#     ------------------[ squarelimit ]------------------
#     -------------[ quartet, cycle, nonet ]-------------
#     ---[ rot, flip, fliprot45, above, beside, over ]---
#     -------[ compose, context, line, path,... ]--------
#
# Drawing this diagram out is a useful way to begin building any library.
#

max_stars_repo_name: st--/jupytext

max_issues_repo_name: st--/jupytext

max_forks_repo_name: st--/jupytext


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# Base.Array related interface

import Base: reshape

"""
    reshape(sym::SymbolicNode, dim; reverse=false, name)
    reshape(sym::SymbolicNode, dim...; reverse=false, name)

Reshape SymbolicNode operator

Some dimensions of the shape can take special values from the set
{0, -1, -2, -3, -4}.
The significance of each is explained below:

- `0`  copy this dimension from the input to the output shape.

  Example:

  - input shape = (2,3,4), shape = (4,0,2), output shape = (4,3,2)
  - input shape = (2,3,4), shape = (2,0,0), output shape = (2,3,4)

- `-1` infers the dimension of the output shape by using the remainder of the
  input dimensions keeping the size of the new array same as that of the input
  array. At most one dimension of shape can be -1.

  Example:

  - input shape = (2,3,4), shape = (6,1,-1), output shape = (6,1,4)
  - input shape = (2,3,4), shape = (3,-1,8), output shape = (3,1,8)
  - input shape = (2,3,4), shape=(-1,), output shape = (24,)

- `-2` copy all/remainder of the input dimensions to the output shape.

  Example:

  - input shape = (2,3,4), shape = (-2,), output shape = (2,3,4)
  - input shape = (2,3,4), shape = (2,-2), output shape = (2,3,4)
  - input shape = (2,3,4), shape = (-2,1,1), output shape = (2,3,4,1,1)

- `-3` use the product of two consecutive dimensions of the input shape as the
  output dimension.

  Example:

  - input shape = (2,3,4), shape = (-3,4), output shape = (6,4)
  - input shape = (2,3,4,5), shape = (-3,-3), output shape = (6,20)
  - input shape = (2,3,4), shape = (0,-3), output shape = (2,12)
  - input shape = (2,3,4), shape = (-3,-2), output shape = (6,4)

- `-4` split one dimension of the input into two dimensions passed subsequent
  to -4 in shape (can contain -1).

  Example:

  - input shape = (2,3,4), shape = (-4,1,2,-2), output shape = (1,2,3,4)
  - input shape = (2,3,4), shape = (2,-4,-1,3,-2), output shape = (2,1,3,4)

If the argument `reverse` is set to `1`, then the special values are inferred
from right to left.

  Example:

  - with `reverse=false`, for input shape = (10,5,4), shape = (-1,0),
    output shape would be (40,5)
  - with `reverse=true`, output shape will be (50,4).
"""
reshape(sym::SymbolicNode, dim::NTuple{N, Integer}; kwargs...) where {N} =
  _reshape(sym, dim; kwargs...)
reshape(sym::SymbolicNode, dim::Integer...; kwargs...) =
  _reshape(sym, dim; kwargs...)

@inline function _reshape(sym::SymbolicNode, dim::NTuple{N,Integer};
                          reverse::Bool=false, name::String="") where N
  op = _get_cached_libmx_op_handle("reshape")
  node = _create_atomic_symbol(op.value, ["shape", "reverse"],
                               [dump_mx_param(dim), dump_mx_param(!reverse)])
  name = get!(DEFAULT_NAME_MANAGER, name, "reshape")
  _compose!(node, name=name, data=sym)
end

################################################################################
# Base.getindex
################################################################################

"""
    getindex(self :: SymbolicNode, idx :: Union{Int, Base.Symbol, AbstractString})

Get a node representing the specified output of this node. The index could be
a symbol or string indicating the name of the output, or a 1-based integer
indicating the index, as in the list of [`list_outputs`](@ref).
"""
function Base.getindex(self :: SymbolicNode, idx :: Union{Base.Symbol, AbstractString})
  idx   = Symbol(idx)
  i_idx = findall(idx .== list_outputs(self))
  @assert(length(i_idx) > 0, "Cannot find output with name '$idx'")
  @assert(length(i_idx) < 2, "Found duplicated output with name '$idx'")
  Base.getindex(self, i_idx[1])
end
function Base.getindex(self :: SymbolicNode, idx :: Int)
  ref_hdr = Ref{MX_handle}(0)
  # note Julia is 1-based, while MXNet is 0-based
  @mxcall(:MXSymbolGetOutput, (MX_handle, MX_uint, Ref{MX_handle}), self, idx-1, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end


text:
using ITensorParallel
using Documenter

DocMeta.setdocmeta!(ITensorParallel, :DocTestSetup, :(using ITensorParallel); recursive=true)

makedocs(;
    modules=[ITensorParallel],
    authors="Matthew Fishman <mfishman@flatironinstitute.org> and contributors",
    repo="https://github.com/mtfishman/ITensorParallel.jl/blob/{commit}{path}#{line}",
    sitename="ITensorParallel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mtfishman.github.io/ITensorParallel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mtfishman/ITensorParallel.jl",
    devbranch="main",
)

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
using Statistics
using Plots

include("ga.jl")

using .GA

str_toint(s) = parse(Int, join(Int.(s)), base=2)

fitness(x::Integer) = 20 * (cos(x/3 + 1)/(x/3 + 1))
fitness(x::Vector{Bool}) = fitness(str_toint(x))

x = 0:63
y = fitness.(x)

@show x[argmin(y)]
@show y[argmin(y)]

ff_plt = plot(
    x,
    y,
    label=nothing,
    title="Fitness Function",
    size=(1024, 512),
)
scatter!(x, y, label=nothing)
savefig(ff_plt, "figs/f_f_fig.png")

let quz="uniform-1", n_pop=8, n_genes=ceil(Integer, log2(63)), uniform_p=1.0, n_repeat=100, stop_bybest=true
    ga01_his = SimpleGA[]
    ga05_his = SimpleGA[]
    ga10_his = SimpleGA[]
    ga50_his = SimpleGA[]

    for i in 1:n_repeat
        ga01 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.01,
            uniform_p=uniform_p,
        )
        optimize!(ga01, stop_bybest=stop_bybest)
        push!(ga01_his, ga01)

        ga05 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.05,
        )
        optimize!(ga05, stop_bybest=stop_bybest)
        push!(ga05_his, ga05)

        ga10 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.1,
        )
        optimize!(ga10, stop_bybest=stop_bybest)
        push!(ga10_his, ga10)

        ga50 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.5,
        )
        optimize!(ga50, stop_bybest=stop_bybest)
        push!(ga50_his, ga50)
    end

    mean_f(his) = mean(fitness.(broadcast(x -> x.best_chro, his)))
    var_f(his) = var(fitness.(broadcast(x -> x.best_chro, his)))

    mean_i(his) = mean(broadcast(x -> x.iter_n, his))
    var_i(his) = var(broadcast(x -> x.iter_n, his))

    fgm_c(his) = sum(broadcast(x -> str_toint(x.best_chro) == 5, his))

    mut_plt_gm = plot(
        [1, 5, 10, 50],
        fgm_c.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Global Minimum Count",
        size=(1024, 512),
    )
    savefig(mut_plt_gm, "figs/$(quz)_mut_gm_fig.png")

    mut_plt_f_m = plot(
        [1, 5, 10, 50],
        mean_f.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Fitness Mean",
        size=(1024, 512),
    )
    savefig(mut_plt_f_m, "figs/$(quz)_mut_f_m_fig.png")
    mut_plt_f_v = plot(
        [1, 5, 10, 50],
        var_f.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Fitness Variance",
        size=(1024, 512),
    )
    savefig(mut_plt_f_v, "figs/$(quz)_mut_f_v_fig.png")

    mut_plt_i_m = plot(
        [1, 5, 10, 50],
        mean_i.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - N.Iterations Mean",
        size=(1024, 512),
    )
    savefig(mut_plt_i_m, "figs/$(quz)_mut_i_m_fig.png")
    mut_plt_i_v = plot(
        [1, 5, 10, 50],
        var_i.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - N.Iterations Variance",
        size=(1024, 512),
    )
    savefig(mut_plt_i_v, "figs/$(quz)_mut_i_v_fig.png")
end

let quz="uniform-0", n_pop=8, n_genes=ceil(Integer, log2(63)), uniform_p=0.0, n_repeat=100, stop_bybest=true
    ga01_his = SimpleGA[]
    ga05_his = SimpleGA[]
    ga10_his = SimpleGA[]
    ga50_his = SimpleGA[]

    for i in 1:n_repeat
        ga01 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.01,
            uniform_p=uniform_p,
        )
        optimize!(ga01, stop_bybest=stop_bybest)
        push!(ga01_his, ga01)

        ga05 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.05,
        )
        optimize!(ga05, stop_bybest=stop_bybest)
        push!(ga05_his, ga05)

        ga10 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.1,
        )
        optimize!(ga10, stop_bybest=stop_bybest)
        push!(ga10_his, ga10)

        ga50 = SimpleGA(
            n_pop=n_pop,
            n_genes=n_genes,
            fitness=fitness,
            mutation_p=0.5,
        )
        optimize!(ga50, stop_bybest=stop_bybest)
        push!(ga50_his, ga50)
    end

    mean_f(his) = mean(fitness.(broadcast(x -> x.best_chro, his)))
    var_f(his) = var(fitness.(broadcast(x -> x.best_chro, his)))

    mean_i(his) = mean(broadcast(x -> x.iter_n, his))
    var_i(his) = var(broadcast(x -> x.iter_n, his))

    fgm_c(his) = sum(broadcast(x -> str_toint(x.best_chro) == 5, his))

    mut_plt_gm = plot(
        [1, 5, 10, 50],
        fgm_c.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Global Minimum Count",
        size=(1024, 512),
    )
    savefig(mut_plt_gm, "figs/$(quz)_mut_gm_fig.png")

    mut_plt_f_m = plot(
        [1, 5, 10, 50],
        mean_f.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Fitness Mean",
        size=(1024, 512),
    )
    savefig(mut_plt_f_m, "figs/$(quz)_mut_f_m_fig.png")
    mut_plt_f_v = plot(
        [1, 5, 10, 50],
        var_f.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - Fitness Variance",
        size=(1024, 512),
    )
    savefig(mut_plt_f_v, "figs/$(quz)_mut_f_v_fig.png")

    mut_plt_i_m = plot(
        [1, 5, 10, 50],
        mean_i.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - N.Iterations Mean",
        size=(1024, 512),
    )
    savefig(mut_plt_i_m, "figs/$(quz)_mut_i_m_fig.png")
    mut_plt_i_v = plot(
        [1, 5, 10, 50],
        var_i.([ga01_his, ga05_his, ga10_his, ga50_his]),
        label=nothing,
        title="Mutation Change - N.Iterations Variance",
        size=(1024, 512),
    )
    savefig(mut_plt_i_v, "figs/$(quz)_mut_i_v_fig.png")
end


text:
using RandomMatrixDistributions
using Distributions
using Test

@testset "MatrixSamplers" begin

    n = 1000
    p = 400
    γ = p/n

    # test white Wishart
    for beta in [1, 4, 8]
        #println("Testing white Wishart with β = ", beta)

        λs = randeigvals(SpikedWishart(beta, n, p))

        @test length(λs) == p
        @test eltype(λs) <: Real

        @test minimum(λs) > 0

        #println("Testing Jacobi with β = ", beta)
        λs = randeigvals(Jacobi(beta, n, n, p))

        @test length(λs) == p
        @test eltype(λs) <: Real

        @test minimum(λs) > 0
    end

    #println("Testing Spiked Wishart with β = 1")
    s = 8
    W = SpikedWishart(1, n, p, [0.1, 2, 4, s])
    λs = randeigvals(W)/n

    @test length(λs) == p
    @test eltype(λs) <: Real

    # Check if the maximum spike is reasonable
    # This tests can fail in principle, but only with negligible probability
    #μ = (1 + s) * (1 + γ/s)
    #σ = (1 + s) * sqrt(γ * 2 * (1 - γ/s^2))

    spikedist = covspikedist(W)
    μ = mean(spikedist)[end]
    σ = sqrt(var(spikedist)[end])

    @test abs(maximum(λs) - μ)/σ < 4

    @test all(randeigstat(W, eigmax, 10) .> 0)

end


text:
using MatrixAlgebra
using Test

@testset "arithmetic" begin
    a = RatApp(1.0, 0, 1, 0)
    @test a == one(RatApp{Float64})
end


text:
"""
    GeneralFieldModule

Module for general fields.
"""
module GeneralFieldModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FieldModule.Field
import FinEtools.FieldModule.@add_Field_fields

"""
    GeneralField{T<:Number} <: Field

General field.
"""
mutable struct GeneralField{T<:Number} <: Field
    @add_Field_fields()
end

"""
    GeneralField(data::FMat{T}=[]) where {T<:Number}

Constructor of general field.
"""
function GeneralField(data::FMat{T}=[]) where {T<:Number}
    values = deepcopy(data)
    dofnums = 0*similar(values,FInt)
    is_fixed = similar(values,Bool)
    fill!(is_fixed, 0)
    fixed_values = zeros(T,size(values))
    nfreedofs = 0
    return GeneralField(values, dofnums, is_fixed, fixed_values, nfreedofs)
end

function GeneralField(data::FVec{T}) where {T<:Number}
    return GeneralField(reshape(data, length(data), 1))
end

end


text:
module Emulate
using Distributions
using Random
using LinearAlgebra
using Flux
using FFTW
using Zygote

using Tullio
using CUDA
using CUDAKernels
using KernelAbstractions
using ChainRulesCore

import Distributions: rand
import Distributions: mean, cov

include("gaussian_processes.jl")

greet() = println("Hello World!")

end # module

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


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

using Markdown
using InteractiveUtils

# ╔═╡ feadb782-9905-4e05-a0c9-e1b1b54fb3b0
begin
	using Distributions
	n_samples = 100
	p_1 = 0.2
	x_data = rand(Binomial(1, p_1), n_samples)
end

# ╔═╡ c54da392-7d0b-11eb-278e-b3a7ec910fc7
begin
	using Gaston; set(showable="svg")
	p0 = range(0, 1;length=100)
	p1 = 1 .- p0
	plot(p0, p1, Axes(xrange = (-1, 1), yrange = (-1, 1), zeroaxis="", xtics="axis", ytics="axis"))
	plot!([0.1,1], [1,0.2], supp=["p₀","p₁"], w="labels")
end

# ╔═╡ e42579e4-5e98-4cca-a425-a3b356dd03ba
begin
	using LinearAlgebra
	p = [0.8, 0.2]
	norm(p, 1)
end

# ╔═╡ 6de3cf91-1739-49ac-aac5-0dcd72871461
begin
	using Yao
	q = zero_state(1)
end

# ╔═╡ 0ef307d4-84d3-11eb-10e7-07cc58485fba
include("./essentials.jl"); import .essentials

# ╔═╡ 4f69ca95-0111-40d3-b74b-34e999f4f185
md"""
Probability theory is a cornerstone for machine learning. We can think of quantum states as probability distributions with certain properties that make them different from our classical notion of probabilities. Contrasting these properties is an easy and straightforward introduction to the most basic concepts we need in quantum computing.

Apart from probability theory, linear algebra is also critical for many learning protocols. As we will see, geometry and probabilities are intrinsically linked in quantum computing, but geometric notions are also familiar in dealing with classical probability distributions. This notebook first talks about classical probabilities and stochastic vectors, and introduces quantum states as a natural generalization.

Throughout this course, we will assume finite probability distributions and finite dimensional spaces. This significantly simplifies notation and most quantum computers operate over finite dimensional spaces, so we do not lose much in generality.

# Classical probability distributions

Let us toss a biased coin. Without getting too technical, we can associate a random variable $X$ with the output: it takes the value 0 for heads and the value 1 for tails. We get heads with probability $P(X=0) = p_0$ and tails with $P(X=1) = p_1$ for each toss of the coin. In classical, Kolmogorovian probability theory, $p_i\geq 0$ for all $i$, and the probabilities sum to one: $\sum_i p_i = 1$. Let's sample this distribution
"""

# ╔═╡ c3898a67-d4b2-46a0-b5cc-69420762cd62
md"""
We naturally expect that the empirically observed frequencies also sum to one:
"""

# ╔═╡ 02a05ce8-35ad-4192-a9d6-e7d3cc29e6ba
begin
	frequency_of_zeros, frequency_of_ones = 0, 0
	for x in x_data
	    if x==1
	        frequency_of_ones += 1/n_samples
	    else
    	    frequency_of_zeros += 1/n_samples
		end
	end
	frequency_of_ones+frequency_of_zeros
end

# ╔═╡ a776d9b8-9adb-4985-9830-165acb8d805b
md"""
Since $p_0$ and $p_1$ must be non-negative, all possible probability distributions are restricted to the positive orthant. The normalization constraint puts every possible distribution on a straight line. This plot describes all possible probability distributions by biased and unbiased coins.
"""

# ╔═╡ 8c7034ce-8286-4d06-9fe1-17581695219c
md"""
We may also arrange the probabilities in a vector $\vec{p} = \begin{bmatrix} p_0 \\ p_1 \end{bmatrix}$. Here, for notational convenience, we put an arrow above the variable representing the vector, to distinguish it from scalars. You will see that quantum states also have a standard notation that provides convenience, but goes much further in usefulness than the humble arrow here.

A vector representing a probability distribution is called a *stochastic vector*. The normalization constraint essentially says that the norm of the vector is restricted to one in the $l_1$ norm. In other words, $||\vec{p}||_1 = \sum_i |p_i| = 1$. This would be the unit circle in the $l_1$ norm, but since $p_i\geq 0$, we are restricted to a quarter of the unit circle, just as we plotted above. We can easily verify this with numpy's norm function:
"""

# ╔═╡ 13533ef2-3be4-4aea-981d-0bc77ec5db57
md"""
We know that the probability of heads is just the first element in the $\vec{p}$, but since it is a vector, we could use linear algebra to extract it. Geometrically, it means that we project the vector to the first axis. This projection is described by the matrix $\begin{bmatrix} 1 & 0\\0 & 0\end{bmatrix}$. The length in the $l_1$ norm gives the sought probability:
"""

# ╔═╡ 0d931337-ad9b-45d1-85fa-2b89f052099d
begin
	Π₀ = [1 0;0 0]
	norm(Π₀ * p)
end

# ╔═╡ 29c91dc7-cfda-42d4-bd0f-f479c611339e
md"""
We can repeat the process to get the probability of tails:
"""

# ╔═╡ c1ec3e28-93fc-42d5-8ef7-454cdb7f7c33
begin
	Π₁ = [0 0; 0 1]
	norm(Π₁ * p)
end

# ╔═╡ de543889-d698-4926-a2a2-9aea48578fd8
md"""
The two projections play an equivalent role to the values 0 and 1 when we defined the probability distribution. In fact, we could define a new random variable called $\Pi$ that can take the projections $\Pi_0$ and $\Pi_1$ as values and we would end up with an identical probability distribution. This may sound convoluted and unnatural, but the measurement in quantum mechanics is essentially a random variable that takes operator values, such as projections.
"""

# ╔═╡ caa23a10-d43e-4a8c-b987-af69b55ed16e
md"""
What happens when we want to transform a probability distribution to another one? For instance, to change the bias of a coin, or to describe the transition of a Markov chain. Since the probability distribution is also a stochastic vector, we can apply a matrix on the vector, where the matrix has to fulfill certain conditions. A left *stochastic matrix* will map stochastic vectors to stochastic vectors when multiplied from the left: its columns add up to one. In other words, it maps probability distributions to probability distributions. For example, starting with a unbiased coin, the map $M$ will transform the distribution to a biased coin:
"""

# ╔═╡ cac5a2c0-6913-4997-b29d-154250be44dd
let
	p = [0.5, 0.5]
	M = [0.7 0.6; 0.3 0.4]
	norm(M * p)
end

# ╔═╡ 1dee601c-2f9f-46d8-be37-99d31e589ebe
md"""
One last concept that will come handy is entropy. A probability distribution's entropy is defined as $H(p) = - \sum_i p_i \log_2 p_i$. Let us plot it over all possible probability distributions of coin tosses:
"""

# ╔═╡ 15416de9-5fe1-4cb7-a5bb-b4227fd9ad4a
let
	ϵ = 10^-10
	p_0 = range(ϵ, 1-ϵ;length = 100)
	p_1 = 1 .- p_0
	H = -1 .* (p_0 .* log.(2,p_0) .+ p_1 .* log.(2, p_1))
	plot(p_0, H, Axes(ylabel="'H'", xlabel="'p₀'", xrange = (0,1), yrange = (0,1)))
	plot!([0.5, 0.5], [0,1], ls = :dash)
end

# ╔═╡ d803bf96-6cd4-4874-a3f7-ed7c0bf16267
md"""
Here we can see that the entropy is maximal for the unbiased coin. This is true in general: the entropy peaks for the uniform distribution. In a sense, this is the most unpredictable distribution: if we get heads with probability 0.2, betting tails is a great idea. On the other hand, if the coin is unbiased, then a deterministic strategy is of little help in winning. Entropy quantifies this notion of surprise and unpredictability.
"""

# ╔═╡ 4cf8cd38-ac26-48f5-8507-9ac56fb90867
md"""
# Quantum states

A classical coin is a two-level system: it is either heads or tails. At a first look a quantum state is a probability distribution, and the simplest case is a two-level state, which we call a qubit. Just like the way we can write the probability distribution as a column vector, we can write a quantum state as a column vector. For notational convenience that will become apparent later, we write the label of a quantum state in what is called a ket in the Dirac notation. So for instance, for some qubit, we can write

``|\psi\rangle = \begin{bmatrix} a_0 \\ a_1 \\ \end{bmatrix}.``

In other words, a ket is just a column vector, exactly like the stochastic vector in the classical case. Instead of putting an arrow over the name of the variable to express that it is a vector, we use the ket to say that it is a column vector that represents a quantum state. There's more to this notation, as we will see.

The key difference to classical probability distributions and stochastic vectors is the normalization constraint. The square sum of their absolute values adds up to 1:

$\sqrt{|a_0|^2+|a_1|^2}=1,$

where $a_0, a_1\in \mathbb{C}$. In other words, we are normalizing in the $l_2$ norm instead of the $l_1$ norm. Furthermore, we are no longer restricted to the positive orthant: the components of the quantum state vector, which we call *probability amplitudes*, are complex valued.

Let us introduce two special qubits, corresponding to the canonical basis vectors in two dimensions: $|0\rangle$ and $|1\rangle$.

$|0\rangle = \begin{bmatrix} 1 \\ 0 \\ \end{bmatrix}, \,\,\, |1\rangle = \begin{bmatrix} 0 \\ 1 \\ \end{bmatrix}.$

This basis is also called the computational basis in quantum computing.

We can expand an arbitrary qubit state in this basis:

$|\psi\rangle = \begin{bmatrix} a_0 \\ a_1 \\ \end{bmatrix}=a_0\begin{bmatrix} 1 \\ 0 \\ \end{bmatrix} + a_1 \begin{bmatrix} 0 \\ 1 \\ \end{bmatrix}= a_0|0\rangle+a_1|1\rangle.$

This expansion in a basis is called a superposition. If we sample the qubit state, we obtain the outcome 0 with probability $|a_0|^2$, and 1 with probability $|a_1|^2$. This is known as the Born rule; you will learn more about measurements and this rule in a subsequent notebook.

For now, let's take a look at how we can simulate classical coin tossing on a quantum computer. Let's start with a completely biased case where we get heads with probability 1. This means that our qubit $|\psi\rangle=|0\rangle$. We create a circuit of a single qubit and a single classical register where the results of the sampling (measurements) go.
"""

# ╔═╡ 1b623852-7074-4211-bbcd-a65ac55c3070
md"""
Any qubit is initialized in $|0\rangle$, so if we measure it right away, we should get our maximally biased coin.
"""

# ╔═╡ 9e2a3016-ac1a-447c-a223-3c4b80c62025
measure(q)

# ╔═╡ 3e692d90-58ce-4562-a3fa-4d73e24b7b62
md"""
Let us execute it a hundred times and study the result
"""

# ╔═╡ ee624636-c6fa-434a-bd6a-7048e755dffa
result = measure(q, nshots=100)

# ╔═╡ 0ebecb30-7db3-11eb-0dab-6ba6b1b2522a
begin
	using StatsBase
	[fit(Histogram, Int.(result), 0:2).weights string.(0:1)]
end

# ╔═╡ d6fc78fa-ba2e-45bf-af00-45fee2fd017e
md"""
As expected, all of our outcomes are $0$.
To understand the possible quantum states, we use the Bloch sphere visualization. Since the probability amplitudes are complex and there are two of them for a single qubit, this would require a four-dimensional space. Now since the vectors are normalized, this removes a degree of freedom, allowing a three-dimensional representation with an appropriate embedding. This embedding is the Bloch sphere. It is slightly different than an ordinary sphere in three dimensions: we identify the north pole with the state $|0\rangle$, and the south pole with $|1\rangle$. In other words, two orthogonal vectors appear as if they were on the same axis -- the axis Z. The computational basis is just one basis: the axes X and Y represent two other bases. Any point on the surface of this sphere is a valid quantum state. This is also true the other way around: every pure quantum state is a point on the Bloch sphere. Here it 'pure' is an important technical term and it essentially means that the state is described by a ket (column vector). Later in the course we will see other states called mix states that are not described by a ket (you will see later that these are inside the Bloch sphere).

To make it less abstract, let's plot our $|0\rangle$ on the Bloch sphere:
"""

# ╔═╡ ba3d9b32-88dd-11eb-3b2a-0f96968e5bcb
essentials.plot_blochsph(state(zero_state(1)))

# ╔═╡ 63946ad1-d28d-4759-8898-085f189e57db
md"""
Compare this sphere with the straight line in the positive orthant that describes all classical probability distributions of coin tosses. You can already see that there is a much richer structure in the quantum probability space.

Let us pick another point on the Bloch sphere, that is, another distribution. Let's transform the state $|0\rangle$ to $\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$. This corresponds to the unbiased coin, since we will get 0 with probability $|\frac{1}{\sqrt{2}}|^2=1/2$, and the other way around. There are many ways to do this transformation. We pick a rotation around the Y axis by $\pi/2$, which corresponds to the matrix $\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & -1\\1 & 1\end{bmatrix}$.
"""

# ╔═╡ a9bfa0e3-b128-43b1-af54-579d2bdcb3b7
essentials.plot_hist(measure(zero_state(1) |> Ry(π/2), nshots=1024))

# ╔═╡ 0ee44bd4-e41f-4977-b71c-9bc04bba0385
md"""
To get an intuition why it is called a rotation around the Y axis, let's plot it on the Bloch sphere:
"""

# ╔═╡ 770ae76e-6521-4613-b781-3627847ffae9
essentials.plot_blochsph(state(zero_state(1) |> Ry(π/2)))

# ╔═╡ 85bea7c1-71ba-4b6b-9341-affd291c0a7c
md"""
It does exactly what it says: it rotates from the north pole of the Bloch sphere.

Why is interesting to have complex probability amplitudes instead of non-negative real numbers? To get some insight, take a look what happens if we apply the same rotation to $|1\rangle$. To achieve this, first we flip $|0\rangle$ to $|1\rangle$ by applying a NOT gate (denoted by X in quantum computing) and then the rotation.
"""

# ╔═╡ a2eaec0e-63a1-4771-899e-240289173b90
essentials.plot_blochsph(state(zero_state(1) |> chain(1, put(1=>X), put(1=>Ry(π/2)))))

# ╔═╡ d054d2e6-da5d-4070-8b93-c81f07af756d
md"""
You can verify that the result is $\frac{1}{\sqrt{2}}(-|0\rangle + |1\rangle)$. That is, the exact same state as before, except that the first term got a minus sign: it is a negative probability amplitude. Note that the difference cannot be observed from the statistics:
"""

# ╔═╡ 0c6bdade-da51-4442-90d4-f3c36aa06e4c
begin
	essentials.plot_hist(measure(zero_state(1) |> chain(1, put(1=>X), put(1=>Ry(π/2))), nshots=1024))
end

# ╔═╡ ca2b7743-b2b3-4335-8c4c-2647492072ed
md"""
It still looks like an approximately unbiased coin. Yet, that negative sign -- or any complex value -- is what models *interference*, a critically important phenomenon where probability amplitudes can interact in a constructive or a destructive way. To see this, if we apply the rotation twice in a row on $|0\rangle$, we get another deterministic output, $|1\rangle$, although in between the two, it was some superposition.
"""

# ╔═╡ 21a24eb3-6104-4f6a-bae5-d11e504aad62
essentials.plot_hist(measure(zero_state(1) |> chain(1, put(1=>Ry(π/2)), put(1=>Ry(π/2))), nshots=1024))

# ╔═╡ 5846d249-d73f-46d7-a1b9-09efc42ef331
md"""
Many quantum algorithms exploit interference, for instance, the seminal [Deutsch-Josza algorithm](https://en.wikipedia.org/wiki/Deutsch–Jozsa_algorithm), which is among the simplest to understand its significance.
"""

# ╔═╡ 2994d86d-c2ba-417b-aa0e-b0c7d0956008
md"""
# More qubits and entanglement

We have already seen that quantum states are probability distributions normed to 1 in the $l_2$ norm and we got a first peek at interference. If we introduce more qubits, we see another crucial quantum effect emerging. To do that, we first have to define how we write down the column vector for describing two qubits. We use a tensor product, which, in the case of qubits, is equivalent to the Kronecker product. Given two qubits, $|\psi\rangle=\begin{bmatrix}a_0\\a_1\end{bmatrix}$ and $|\psi'\rangle=\begin{bmatrix}b_0\\b_1\end{bmatrix}$, their product is $|\psi\rangle\otimes|\psi'\rangle=\begin{bmatrix}a_0b_0\\ a_0b_1\\ a_1b_0\\ a_1b_1\end{bmatrix}$. Imagine that you have two registers $q_0$ and $q_1$, each can hold a qubit, and both qubits are in the state $|0\rangle$. Then this composite state would be described by according to this product rule as follows:
"""

# ╔═╡ f98f70cf-aaf2-4710-9ef7-33c401b10c06
let
	q0 = [1, 0]
	q1 = [1, 0]
	kron(q0, q1)
end

# ╔═╡ 2f1a3452-a122-41f2-8ab7-df4e1cc01b02
md"""
This is the $|0\rangle\otimes|0\rangle$ state, which we often abbreviate as $|00\rangle$. The states $|01\rangle$, $|10\rangle$, and $|11\rangle$ are defined analogously, and the four of them give the canonical basis of the four dimensional complex space, $\mathbb{C}^2\otimes\mathbb{C}^2$.

Now comes the interesting and counter-intuitive part. In machine learning, we also work with high-dimensional spaces, but we never construct it as a tensor product: it is typically $\mathbb{R}^d$ for some dimension $d$. The interesting part of writing the high-dimensional space as a tensor product is that not all vectors in can be written as a product of vectors in the component space.

Take the following state: $|\phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$. This vector is clearly in $\mathbb{C}^2\otimes\mathbb{C}^2$, since it is a linear combination of two of the basis vector in this space. Yet, it cannot be written as $|\psi\rangle\otimes|\psi'\rangle$ for some $|\psi\rangle$, $|\psi'\rangle\in\mathbb{C}^2$.

To see this, assume that it can be written in this form. Then

$|\phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle+|11\rangle) = \begin{bmatrix}a_0b_0\\ a_0b_1\\ a_1b_0\\ a_1b_1\end{bmatrix} = a_0b_0|00\rangle + a_0b_1|01\rangle + a_1b_0|10\rangle + a_1b_1|11\rangle.$

$|01\rangle$ and $|10\rangle$ do not appear on the left-hand side, so their coefficients must be zero: $a_1b_0=0$ and $a_0b_1=0$. This leads to a contradiction, since $a_1$ cannot be zero ($a_1b_1=1$), so $b_0$ must be zero, but $a_0b_0=1$. Therefore $|\phi^+\rangle$ cannot be written as a product.

States that cannot be written as a product are called entangled states. This is the mathematical form of describing a phenomenon of strong correlations between random variables that exceed what is possible classically. Entanglement plays a central role in countless quantum algorithms. A simple example is [quantum teleportation](https://en.wikipedia.org/wiki/Quantum_teleportation). We will also see its applications in quantum machine learning protocols.

We will have a closer look at entanglement in a subsequent notebook on measurements, but as a teaser, let us look at the measurement statistics of the $|\phi^+\rangle$ state. The explanation of the circuit preparing it will also come in a subsequent notebook.
"""

# ╔═╡ a7cdf304-48d2-4485-b473-dd2248bb7ce6
essentials.plot_hist(measure(zero_state(2) |> chain(2, put(1=>H), control(1, 2=>X)), nshots=1024))

# ╔═╡ 5a2944be-209c-4fbe-b1b1-19a5f9df4316
md"""
Notice that 01 or 10 never appear in the measurement statistics.
"""

# ╔═╡ ec3512fa-22de-4123-bb1f-440b899ac901
md"""
# Further reading

Chapter 9 in Quantum Computing since Democritus by Scott Aaronson describes a similar approach to understanding quantum states -- in fact, the interference example was lifted from there.
"""

# ╔═╡ Cell order:
# ╟─4f69ca95-0111-40d3-b74b-34e999f4f185
# ╠═feadb782-9905-4e05-a0c9-e1b1b54fb3b0
# ╟─c3898a67-d4b2-46a0-b5cc-69420762cd62
# ╠═02a05ce8-35ad-4192-a9d6-e7d3cc29e6ba
# ╟─a776d9b8-9adb-4985-9830-165acb8d805b
# ╠═c54da392-7d0b-11eb-278e-b3a7ec910fc7
# ╟─8c7034ce-8286-4d06-9fe1-17581695219c
# ╠═e42579e4-5e98-4cca-a425-a3b356dd03ba
# ╟─13533ef2-3be4-4aea-981d-0bc77ec5db57
# ╠═0d931337-ad9b-45d1-85fa-2b89f052099d
# ╟─29c91dc7-cfda-42d4-bd0f-f479c611339e
# ╠═c1ec3e28-93fc-42d5-8ef7-454cdb7f7c33
# ╟─de543889-d698-4926-a2a2-9aea48578fd8
# ╟─caa23a10-d43e-4a8c-b987-af69b55ed16e
# ╠═cac5a2c0-6913-4997-b29d-154250be44dd
# ╟─1dee601c-2f9f-46d8-be37-99d31e589ebe
# ╠═15416de9-5fe1-4cb7-a5bb-b4227fd9ad4a
# ╟─d803bf96-6cd4-4874-a3f7-ed7c0bf16267
# ╟─4cf8cd38-ac26-48f5-8507-9ac56fb90867
# ╠═6de3cf91-1739-49ac-aac5-0dcd72871461
# ╟─1b623852-7074-4211-bbcd-a65ac55c3070
# ╠═9e2a3016-ac1a-447c-a223-3c4b80c62025
# ╟─3e692d90-58ce-4562-a3fa-4d73e24b7b62
# ╠═ee624636-c6fa-434a-bd6a-7048e755dffa
# ╠═0ebecb30-7db3-11eb-0dab-6ba6b1b2522a
# ╟─d6fc78fa-ba2e-45bf-af00-45fee2fd017e
# ╠═0ef307d4-84d3-11eb-10e7-07cc58485fba
# ╠═ba3d9b32-88dd-11eb-3b2a-0f96968e5bcb
# ╟─63946ad1-d28d-4759-8898-085f189e57db
# ╠═a9bfa0e3-b128-43b1-af54-579d2bdcb3b7
# ╟─0ee44bd4-e41f-4977-b71c-9bc04bba0385
# ╠═770ae76e-6521-4613-b781-3627847ffae9
# ╟─85bea7c1-71ba-4b6b-9341-affd291c0a7c
# ╠═a2eaec0e-63a1-4771-899e-240289173b90
# ╟─d054d2e6-da5d-4070-8b93-c81f07af756d
# ╠═0c6bdade-da51-4442-90d4-f3c36aa06e4c
# ╟─ca2b7743-b2b3-4335-8c4c-2647492072ed
# ╠═21a24eb3-6104-4f6a-bae5-d11e504aad62
# ╟─5846d249-d73f-46d7-a1b9-09efc42ef331
# ╟─2994d86d-c2ba-417b-aa0e-b0c7d0956008
# ╠═f98f70cf-aaf2-4710-9ef7-33c401b10c06
# ╟─2f1a3452-a122-41f2-8ab7-df4e1cc01b02
# ╠═a7cdf304-48d2-4485-b473-dd2248bb7ce6
# ╟─5a2944be-209c-4fbe-b1b1-19a5f9df4316
# ╟─ec3512fa-22de-4123-bb1f-440b899ac901

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc


text:
nullable_returntype(::Type{T}) where {T} = isconcretetype(T) ? T : Union{}
# An element type satisfying for all A:
# unsafe_get(A)::unsafe_get_eltype(A)
_unsafe_get_eltype(::Type{T}) where T = Type{T}
_unsafe_get_eltype(x) = typeof(x)

# Copied from julia 0.6
maptoTuple(f) = Tuple{}
maptoTuple(f, a, b...) = Tuple{f(a),maptoTuple(f, b...).types...}

_nullable_eltype(f, A, As...) =
    Base._return_type(f, maptoTuple(_unsafe_get_eltype, A, As...))

function Dates.DateTime(dt::DataValue{T}, format::AbstractString; locale::Dates.Locale=Dates.ENGLISH) where {T <: AbstractString}
    isna(dt) ? DataValue{DateTime}() : DataValue{DateTime}(DateTime(unsafe_get(dt), format, locale=locale))
end

function Dates.Date(dt::DataValue{T}, format::AbstractString; locale::Dates.Locale=Dates.ENGLISH) where {T <: AbstractString}
    isna(dt) ? DataValue{Date}() : DataValue{Date}(Date(unsafe_get(dt), format, locale=locale))
end

for f in (:(Base.abs), :(Base.abs2), :(Base.conj), :(Base.sign))
    @eval begin
        function $f(a::DataValue{T}) where {T}
            if isna(a)
                DataValue{T}()
            else
                DataValue($f(unsafe_get(a)))
            end
        end
    end
end

for f in (:(Base.acos), :(Base.acosh), :(Base.asin), :(Base.asinh),
        :(Base.atan), :(Base.atanh), :(Base.sin), :(Base.sinh), :(Base.cos),
        :(Base.cosh), :(Base.tan), :(Base.tanh), :(Base.exp), :(Base.exp2),
        :(Base.expm1), :(Base.log), :(Base.log10), :(Base.log1p),
        :(Base.log2), :(Base.exponent), :(Base.sqrt), :(Dates.value))
    @eval begin
        function $f(a::DataValue{T}) where {T}
            if isna(a)
                DataValue{Float64}()
            else
                DataValue{Float64}($f(unsafe_get(a)))
            end
        end
    end
end

for op in (:+, :-, :*, :/, :%, :&, :|, :^, :<<, :>>, :div, :mod, :fld,
        :min, :max)
    @eval begin
        import Base.$(op)
        function $op(a::DataValue{T1}, b::DataValue{T2}) where {T1,T2}
            nonnull = hasvalue(a) && hasvalue(b)
            S = _nullable_eltype($op, a, b)
            if nonnull
                return DataValue($op(get(a), get(b)))
            else
                return DataValue{nullable_returntype(S)}()
            end
        end

        function $op(a::DataValue{T1}, b::T2) where {T1,T2}
            nonnull = hasvalue(a)
            S = _nullable_eltype($op, a, b)
            if nonnull
                return DataValue($op(get(a), b))
            else
                return DataValue{nullable_returntype(S)}()
            end
        end

        function $op(a::T1, b::DataValue{T2}) where {T1,T2}
            nonnull = hasvalue(b)
            S = _nullable_eltype($op, a, b)
            if nonnull
                return DataValue($op(a, get(b)))
            else
                return DataValue{nullable_returntype(S)}()
            end
        end
    end
end


text:
# Tools for implementing device functionality

# which Julia types map to a given LLVM type
const llvmtypes = Dict{Type,Symbol}(
    Nothing => :void,
    Int8    => :i8,
    Int16   => :i16,
    Int32   => :i32,
    Int64   => :i64,
    Float32 => :float,
    Float64 => :double,
)

# which LLVM types map to a given Julia type
const jltypes = Dict{Symbol,Type}(
    :void   => Nothing,
    :i8     => Int8,
    :i16    => Int16,
    :i32    => Int32,
    :i64    => Int64,
    :float  => Float32,
    :double => Float64
)

# which Julia types map to a given ROC lib type
const fntypes = Dict{Type,Symbol}(
    Nothing => :void,
    Int8    => :i8,
    Int16   => :i16,
    Int32   => :i32,
    Int64   => :i64,
    Float16 => :f16,
    Float32 => :f32,
    Float64 => :f64
)

# Decode an expression of the form:
#
#    function(arg::arg_type, arg::arg_type, ... arg::arg_type)::return_type
#
# Returns a tuple containing the function name, a vector of argument, a vector of argument
# types and the return type (all in symbolic form).
function decode_call(e)
    @assert e.head == :(::)

    # decode the return type expression: single symbol (the LLVM type), or a tuple of 2
    # symbols (the LLVM and corresponding Julia type)
    retspec = e.args[2]
    if isa(retspec, Symbol)
        rettype = retspec
    else
        @assert retspec.head == :tuple
        @assert length(retspec.args) == 2
        rettype = (retspec.args[1], retspec.args[2])
    end

    call = e.args[1]
    @assert call.head == :call

    fn = Symbol(call.args[1])
    args = Symbol[arg.args[1] for arg in call.args[2:end]]
    argtypes = Symbol[arg.args[2] for arg in call.args[2:end]]

    return fn, args, argtypes, rettype
end

# Generate a `llvmcall` statement calling an intrinsic specified as follows:
#
#     intrinsic(arg::arg_type, arg::arg_type, ... arg::arg_type)::return_type [attr]
#
# The argument types should be valid LLVM type identifiers (eg. i32, float, double).
# Conversions to the corresponding Julia type are automatically generated; make sure the
# actual arguments are of the same type to make these conversions no-ops. The optional
# argument `attr` indicates which LLVM function attributes (such as `readnone` or `nounwind`)
# to add to the intrinsic declaration.

# For example, the following call:
#     `@wrap __some_intrinsic(x::float, y::double)::float`
#
# will yield the following `llvmcall`:
# ```
#     Base.llvmcall(("declare float @__somme__intr(float, double)",
#                    "%3 = call float @__somme__intr(float %0, double %1)
#                     ret float %3"),
#                   Float32, Tuple{Float32,Float64},
#                   convert(Float32,x), convert(Float64,y))
# ```
macro wrap(call, attrs="")
    intrinsic, args, argtypes, rettype = decode_call(call)

    # decide on intrinsic return type
    if isa(rettype, Symbol)
        # only LLVM return type specified, match against known LLVM/Julia type combinations
        llvm_ret_typ = rettype
        julia_ret_typ = jltypes[rettype]
    else
        # both specified (for when there is a mismatch, eg. i32 -> UInt32)
        llvm_ret_typ = rettype[1]
        julia_ret_typ = rettype[2]
    end

    llvm_args = String["%$i" for i in 0:length(argtypes)]
    if llvm_ret_typ == :void
        llvm_ret_asgn = ""
        llvm_ret = "void"
    else
        llvm_ret_var = "%$(length(argtypes)+1)"
        llvm_ret_asgn = "$llvm_ret_var = "
        llvm_ret = "$llvm_ret_typ $llvm_ret_var"
    end
    llvm_declargs = join(argtypes, ", ")
    llvm_defargs = join(("$t $arg" for (t,arg) in zip(argtypes, llvm_args)), ", ")

    julia_argtypes = (jltypes[t] for t in argtypes)
    julia_args = (:(convert($argtype, $(esc(arg)))) for (arg, argtype) in zip(args, julia_argtypes))

    dest = ("""declare $llvm_ret_typ @$intrinsic($llvm_declargs)""",
            """$llvm_ret_asgn call $llvm_ret_typ @$intrinsic($llvm_defargs)
                ret $llvm_ret""")

    return quote
        Base.llvmcall($dest, $julia_ret_typ, Tuple{$(julia_argtypes...)}, $(julia_args...))
    end
end

# julia.h: jl_datatype_align
Base.@pure function datatype_align(::Type{T}) where {T}
    # typedef struct {
    #     uint32_t nfields;
    #     uint32_t alignment : 9;
    #     uint32_t haspadding : 1;
    #     uint32_t npointers : 20;
    #     uint32_t fielddesc_type : 2;
    # } jl_datatype_layout_t;
    field = T.layout + sizeof(UInt32)
    unsafe_load(convert(Ptr{UInt16}, field)) & convert(Int16, 2^9-1)
end

# generalization of word-based primitives

# extract bits from a larger value
@inline function extract_word(val, ::Val{i}) where {i}
    extract_value(val, UInt32, Val(32*(i-1)))
end

@generated function extract_value(val, ::Type{sub}, ::Val{offset}) where {sub, offset}
    T_val = convert(LLVMType, val)
    T_sub = convert(LLVMType, sub)
    bytes = Core.sizeof(val)
    T_int = LLVM.IntType(8*bytes, JuliaContext())

    # create function
    llvm_f, _ = create_function(T_sub, [T_val])
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)
        equiv = bitcast!(builder, parameters(llvm_f)[1], T_int)
        shifted = lshr!(builder, equiv, LLVM.ConstantInt(T_int, offset))
        # extracted = and!(builder, shifted, 2^32-1)
        extracted = trunc!(builder, shifted, T_sub)
        ret!(builder, extracted)
    end

    call_function(llvm_f, UInt32, Tuple{val}, :( (val,) ))
end

# insert bits into a larger value
@inline function insert_word(val, word::UInt32, ::Val{i}) where {i}
    insert_value(val, word, Val(32*(i-1)))
end

@generated function insert_value(val, sub, ::Val{offset}) where {offset}
    T_val = convert(LLVMType, val)
    T_sub = convert(LLVMType, sub)
    bytes = Core.sizeof(val)
    T_out_int = LLVM.IntType(8*bytes, JuliaContext())

    # create function
    llvm_f, _ = create_function(T_val, [T_val, T_sub])
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)
        equiv = bitcast!(builder, parameters(llvm_f)[1], T_out_int)
        ext = zext!(builder, parameters(llvm_f)[2], T_out_int)
        shifted = shl!(builder, ext, LLVM.ConstantInt(T_out_int, offset))
        inserted = or!(builder, equiv, shifted)
        orig = bitcast!(builder, inserted, T_val)
        ret!(builder, orig)
    end

    call_function(llvm_f, val, Tuple{val, sub}, :( (val, sub) ))
end

# split the invocation of a function `op` on a value `val` with non-struct eltype
# into multiple smaller invocations on byte-sized partial values.
@generated function split_value_invocation(op::Function, val, args...)
    # TODO: control of lower-limit

    ex = quote
        Base.@_inline_meta
    end

    # disassemble into words
    words = Symbol[]
    for i in 1:Core.sizeof(val)÷4
        word = Symbol("word$i")
        push!(ex.args, :( $word = extract_word(val, Val($i)) ))
        push!(words, word)
    end

    # perform the operation
    for word in words
        push!(ex.args, :( $word = op($word, args...)) )
    end

    # reassemble
    push!(ex.args, :( out = zero(val) ))
    for (i,word) in enumerate(words)
        push!(ex.args, :( out = insert_word(out, $word, Val($i)) ))
    end

    push!(ex.args, :( out ))

    return ex
end

# split the invocation of a function `op` on a value `val`
# by invoking the function on each of its fields
@generated function recurse_value_invocation(op::Function, val, args...)
    ex = quote
        Base.@_inline_meta
    end

    fields = fieldnames(val)
    if isempty(fields)
        push!(ex.args, :( split_value_invocation(op, val, args...) ))
    else
        ctor = Expr(:new, val)
        for field in fields
            push!(ctor.args, :(
                recurse_value_invocation(op, getfield(val, $(QuoteNode(field))), args...) ))
        end
        push!(ex.args, ctor)
    end

    return ex
end

# split the invocation of a function `op` on a pointer `ptr` with non-struct eltype
# into multiple smaller invocations on any supported pointer as listed in `supported_ptrs`.
@generated function split_pointer_invocation(op::Function, ptr, ::Type{supported_ptrs},
                                             args...) where {supported_ptrs}
    T = eltype(ptr)
    elsize(x) = Core.sizeof(eltype(x))
    supported_ptrs = reverse(Base.uniontypes(supported_ptrs))

    ex = quote
        Base.@_inline_meta
    end

    # disassemble
    vals = Tuple{Symbol,Int,Type}[]
    offset = 0
    while offset < Core.sizeof(T)
        val = Symbol("value.$(length(vals)+1)")

        # greedy selection of next pointer type
        remaining = Core.sizeof(T)-offset
        valid = filter(ptr->elsize(ptr)<=remaining, supported_ptrs)
        if isempty(valid)
            error("Cannot partition $T into values of $supported_typs")
        end
        ptr = first(sort(collect(valid); by=elsize, rev=true))

        push!(vals, (val, offset, ptr))
        offset += elsize(ptr)
    end

    # perform the operation
    for (val, offset, ptr) in vals
        subptr = :(convert($ptr, ptr+$offset))
        push!(ex.args, :( $val = op($subptr, args...)) )
    end

    # reassemble
    push!(ex.args, :( out = zero($T) ))
    for (val, offset, ptr) in vals
        push!(ex.args, :( out = insert_value(out, $val, Val($offset)) ))
    end

    push!(ex.args, :( out ))

    return ex
end

# split the invocation of a function `op` on a pointer `ptr`
# by invoking the function on a pointer to each of its fields
@generated function recurse_pointer_invocation(op::Function, ptr, ::Type{supported_ptrs},
                                               args...) where {supported_ptrs}
    T = eltype(ptr)
    ex = quote
        Base.@_inline_meta
    end
    fields = fieldnames(T)

    if isempty(fields)
        push!(ex.args, :( split_pointer_invocation(op, ptr, supported_ptrs, args...) ))
    else
        ctor = Expr(:new, T)
        for (i,field) in enumerate(fields)
            field_typ = fieldtype(T, i)
            field_offset = fieldoffset(T, i)
            field_ptr_typ = :($(ptr.name.wrapper){$field_typ})
            # NOTE: this ctor is a leap of faith
            subptr = :(convert($field_ptr_typ, ptr+$field_offset))
            push!(ctor.args, :(
                recurse_pointer_invocation(op, $subptr, supported_ptrs, args...) ))
        end
        push!(ex.args, ctor)
    end

    return ex
end

# calculate the size of an LLVM type
llvmsize(::LLVM.LLVMHalf) = sizeof(Float16)
llvmsize(::LLVM.LLVMFloat) = sizeof(Float32)
llvmsize(::LLVM.LLVMDouble) = sizeof(Float64)
llvmsize(::LLVM.IntegerType) = div(Int(intwidth(GenericValue(LLVM.Int128Type(), -1))), 8)
llvmsize(ty::LLVM.ArrayType) = length*llvmsize(eltype(ty))
# TODO: VectorType, StructType, PointerType
llvmsize(ty) = error("Unknown size for type: $ty, typeof: $(typeof(ty))")


text:
@doc raw"""
    Oblique{N,M,𝔽} <: AbstractPowerManifold{𝔽}

The oblique manifold $\mathcal{OB}(n,m)$ is the set of 𝔽-valued matrices with unit norm
column endowed with the metric from the embedding. This yields exactly the same metric as
considering the product metric of the unit norm vectors, i.e. [`PowerManifold`](@ref) of the
$(n-1)$-dimensional [`Sphere`](@ref).

The [`Sphere`](@ref) is stored internally within `M.manifold`, such that all functions of
[`AbstractPowerManifold`](@ref) can be used directly.

# Constructor

    Oblique(n,m)

Generate the manifold of matrices $\mathbb R^{n × m}$ such that the $m$ columns are unit
vectors, i.e. from the [`Sphere`](@ref)`(n-1)`.
"""
struct Oblique{N,M,𝔽,S} <:
       AbstractPowerManifold{𝔽,Sphere{S,𝔽},ArrayPowerRepresentation} where {N,M}
    manifold::Sphere{S,𝔽}
end

function Oblique(n::Int, m::Int, field::AbstractNumbers = ℝ)
    return Oblique{n,m,field,n - 1}(Sphere(n - 1, field))
end

Base.:^(M::Sphere{N,𝔽}, m::Int) where {N,𝔽} = Oblique{manifold_dimension(M) + 1,m,𝔽,N}(M)

@doc raw"""
    check_manifold_point(M::Oblique{n,m},p)

Checks whether `p` is a valid point on the [`Oblique`](@ref)`{m,n}` `M`, i.e. is a matrix
of `m` unit columns from $\mathbb R^{n}$, i.e. each column is a point from
[`Sphere`](@ref)`(n-1)`.
"""
check_manifold_point(::Oblique, ::Any)
function check_manifold_point(M::Oblique{n,m}, p; kwargs...) where {n,m}
    if size(p) != (n, m)
        return DomainError(
            length(p),
            "The matrix in `p` ($(size(p))) does not match the dimension of $(M).",
        )
    end
    return check_manifold_point(PowerManifold(M.manifold, m), p; kwargs...)
end
@doc raw"""
    check_tangent_vector(M::Oblique p, X; check_base_point = true, kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`Oblique`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and columnswise
a tangent vector to the columns of `p` on the [`Sphere`](@ref).
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
"""
function check_tangent_vector(
    M::Oblique{n,m},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {n,m}
    if check_base_point && size(p) != (n, m)
        return DomainError(
            length(p),
            "The matrix `p` ($(size(p))) does not match the dimension of $(M).",
        )
    end
    if size(X) != (n, m)
        return DomainError(
            length(X),
            "The matrix `X` ($(size(X))) does not match the dimension of $(M).",
        )
    end
    return check_tangent_vector(
        PowerManifold(M.manifold, m),
        p,
        X;
        check_base_point = check_base_point,
        kwargs...,
    )
end

get_iterator(M::Oblique{n,m}) where {n,m} = Base.OneTo(m)

@generated function manifold_dimension(::Oblique{n,m,𝔽}) where {n,m,𝔽}
    return (n * real_dimension(𝔽) - 1) * m
end

@generated representation_size(::Oblique{n,m}) where {n,m} = (n, m)

function Base.show(io::IO, ::Oblique{n,m,𝔽}) where {n,m,𝔽}
    return print(io, "Oblique($(n),$(m); field = $(𝔽))")
end


text:
const MR_UP_GOOD = [1,0,0,0,0,0,0,-1,-1]
const COST_UP_GOOD = [1,1,1,1,-1,-1,-1,-1,-1,-1,-1]

function format_and_class( change :: Real ) :: Tuple
    gnum = format( abs(change), commas=true, precision=2 )
    glclass = "";
    glstr = ""
    if change > 20.0
        glstr = "positive_strong"
        glclass = "text-success"
    elseif change > 10.0
        glstr = "positive_med"
        glclass = "text-success"
    elseif change > 0.01
        glstr = "positive_weak"
        glclass = "text-success"
    elseif change < -20.0
        glstr = "negative_strong"
        glclass = "text-danger"
    elseif change < -10
        glstr = "negative_med"
        glclass = "text-danger"
    elseif change < -0.01
        glstr = "negative_weak"
        glclass = "text-danger"
    else
        glstr = "nonsig"
        glclass = "text-body"
        gnum = "";
    end
    ( gnum, glclass, glstr )
end


text:
println("Running sdss tests.")

using Base.Test
import SloanDigitalSkySurvey: SDSS, PSF

const field_dir =
  joinpath(Pkg.dir("SloanDigitalSkySurvey"), "dat", "sample_field")
const run_num = "003900"
const camcol_num = "6"
const field_num = "0269"

function test_least_squares_psf()
  raw_psf_comp =
    SDSS.load_psf_data(field_dir, run_num, camcol_num, field_num, 1);
  psf = PSF.get_psf_at_point(500.0, 500.0, raw_psf_comp);

  opt_result, mu_vec, sigma_vec, weight_vec =
    PSF.fit_psf_gaussians_least_squares(
      psf, K=2, optim_method=:nelder_mead);

  x_mat = PSF.get_x_matrix_from_psf(psf);
  psf_fit = PSF.render_psf(opt_result.minimum, x_mat);

  @test_approx_eq sum((psf_fit - psf) .^ 2) opt_result.f_minimum
  @test 0 < opt_result.f_minimum < 1e-4

end

function test_psf_transforms()

  mu_vec = Vector{Float64}[ Float64[1, 2], Float64[-1, -2], Float64[1, -1] ]
  sigma_vec = Array(Matrix{Float64}, 3)
  sigma_vec[1] = Float64[ 1 0.1; 0.1 1]
  sigma_vec[2] = Float64[ 1 0.3; 0.3 2]
  sigma_vec[3] = Float64[ 0.5 0.2; 0.2 0.5]
  weight_vec = Float64[0.4, 0.6, 0.1]

  par = wrap_parameters(mu_vec, sigma_vec, weight_vec)
  mu_vec_test, sigma_vec_test, weight_vec_test = unwrap_parameters(par)

  using Base.Test

  for k=1:3
    @test_approx_eq mu_vec_test[k] mu_vec[k]
    @test_approx_eq sigma_vec_test[k] sigma_vec[k]
    @test_approx_eq weight_vec_test[k] weight_vec[k]
  end

end

function test_load_psf()
  raw_psf_comp =
    SDSS.load_psf_data(field_dir, run_num, camcol_num, field_num, 1);
  raw_psf = PSF.get_psf_at_point(10.0, 10.0, raw_psf_comp);

  # Check that the raw psf is approximately symmetric
  center_point = (size(raw_psf)[1] + 1) / 2
  @test_approx_eq_eps(sum(collect(1:size(raw_psf)[1]) .* raw_psf) / sum(raw_psf),
                      center_point, 0.5)
  @test_approx_eq_eps(sum(collect(1:size(raw_psf)[1])' .* raw_psf) / sum(raw_psf),
                      center_point, 0.5)

end

function test_load_field()
  b = 1
  band_gain, band_dark_variance =
    SDSS.load_photo_field(field_dir, run_num, camcol_num, field_num)
  nelec, calib_col, sky_grid, sky_x, sky_y, sky_image, wcs =
    SDSS.load_raw_field(field_dir, run_num, camcol_num, field_num,
                        b, band_gain[b])

  @test size(sky_image) == size(nelec)
  @test length(calib_col) == size(nelec)[1]
  @test sky_image[10, 5] == sky_grid[sky_x[10], sky_y[5]]

  nelec_original = deepcopy(nelec)
  SDSS.mask_image!(nelec, field_dir, run_num, camcol_num, field_num, b);

  # Test that at least some pixels are bad
  @test sum(isnan(nelec)) > 0

  # Test that no more than 1% are bad.  Not sure this is a good test
  # in general but it works for the current files.
  @test sum(isnan(nelec)) / prod(size(nelec)) < 0.01
end

function test_load_catalog()
  # Just check that the catalog loads.
  cat_df = SDSS.load_catalog_df(field_dir, run_num, camcol_num, field_num);
  @test size(cat_df, 1) == 805
end

test_load_psf()
test_load_field()
test_load_catalog()
test_least_squares_psf()


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

import JSON

"""
    gen( re1, im1, re2, im2, name )

Generate fixture data and write to file.

# Arguments

* `re1`: real components for first complex number
* `im1`: imaginary components for first complex number
* `re2`: real components for second complex number
* `im2`: imaginary components for second complex number
* `name::AbstractString`: output filename

# Examples

``` julia
julia> re1 = rand( 1000 );
julia> im1 = rand( 1000 );
julia> re2 = rand( 1000 );
julia> im2 = rand( 1000 );
julia> gen( re1, im1, re2, im2, \"data.json\" );
```
"""
function gen( re1, im1, re2, im2, name )
	zre = Array{Float64}( undef, length(re1) );
	zim = Array{Float64}( undef, length(re1) );
	for i in eachindex(re1)
		z = Complex{Float64}( re1[i], im1[i] ) / Complex{Float64}( re2[i], im2[i] );
		zre[ i ] = real( z );
		zim[ i ] = imag( z );
	end

	# Store data to be written to file as a collection:
	data = Dict([
		("re1", re1),
		("im1", im1),
		("re2", re2),
		("im2", im2),
		("qre", zre),
		("qim", zim),
	]);

	# Based on the script directory, create an output filepath:
	filepath = joinpath( dir, name );

	# Write the data to the output filepath as JSON:
	outfile = open( filepath, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# Get the filename:
file = @__FILE__;

# Extract the directory in which this file resides:
dir = dirname( file );

# Large positive real components:
re1 = rand( 500 ) .* 1.0e300;
re2 = rand( 500 ) .* 1.0e300;
im1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
gen( re1, re2, im1, im2, "large_positive_real_components.json" );

# Large negative real components:
re1 = -rand( 500 ) .* 1.0e300;
re2 = -rand( 500 ) .* 1.0e300;
im1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
gen( re1, im1, re2,  im2, "large_negative_real_components.json" );

# Large positive imaginary components:
re1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
re2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im1 = rand( 500 ) .* 1.0e300;
im2 = rand( 500 ) .* 1.0e300;
gen( re1, im1, re2, im2, "large_positive_imaginary_components.json" );

# Large negative imaginary components:
re1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
re2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im1 = -rand( 500 ) .* 1.0e300;
im2 = -rand( 500 ) .* 1.0e300;
gen( re1, im1, re2, im2, "large_negative_imaginary_components.json" );

# Tiny positive real components:
re1 = rand( 500 ) .* 1.0e-324;
re2 = rand( 500 ) .* 1.0e-324;
im1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
gen( re1, im1, re2, im2, "tiny_positive_real_components.json" );

# Tiny negative real components:
re1 = -rand( 500 ) .* 1.0e-324;
re2 = -rand( 500 ) .* 1.0e-324;
im1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
gen( re1, im1, re2, im2, "tiny_negative_real_components.json" );

# Tiny positive imaginary components:
re1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
re2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im1 = rand( 500 ) .* 1.0e-324;
im2 = rand( 500 ) .* 1.0e-324;
gen( re1, im1, re2, im2, "tiny_positive_imaginary_components.json" );

# Tiny negative imaginary components:
re1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
re2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im1 = -rand( 500 ) .* 1.0e-324;
im2 = -rand( 500 ) .* 1.0e-324;
gen( re1, im1, re2, im2, "tiny_negative_imaginary_components.json" );

# Real components different scales:
re1 = rand( 500 ) .* 1.0e200;
re2 = rand( 500 ) .* 1.0e-100;
im1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
gen( re1, im1, re2, im2, "real_component_scales.json" );

# Imaginary components different scales:
re1 = ( rand( 500 ) .* 20.0 ) .- 10.0;
re2 = ( rand( 500 ) .* 20.0 ) .- 10.0;
im1 = rand( 500 ) .* 1.0e200;
im2 = rand( 500 ) .* 1.0e-100;
gen( re1, im1, re2, im2, "imaginary_component_scales.json" );

# Components different scales:
re1 = rand( 500 ) .* 1.0e200;
re2 = rand( 500 ) .* 1.0e200;
im1 = ( rand( 500 ) .* 1.0e-200 ) .- 0.5e-200;
im2 = ( rand( 500 ) .* 1.0e-200 ) .- 0.5e-200;
gen( re1, im1, re2, im2, "component_scales1.json" );

# Components different scales:
re1 = ( rand( 500 ) .* 1.0e-200 ) .- 0.5e-200;
re2 = ( rand( 500 ) .* 1.0e-200 ) .- 0.5e-200;
im1 = rand( 500 ) .* 1.0e200;
im2 = rand( 500 ) .* 1.0e200;
gen( re1, im1, re2, im2, "component_scales2.json" );

# Normal:
re1 = ( rand( 500 ) .* 100.0 ) .- 50.0;
re2 = ( rand( 500 ) .* 100.0 ) .- 50.0;
im1 = ( rand( 500 ) .* 100.0 ) .- 50.0;
im2 = ( rand( 500 ) .* 100.0 ) .- 50.0;
gen( re1, im1, re2, im2, "data.json" );

max_stars_repo_name: ghalimi/stdlib

max_issues_repo_name: ghalimi/stdlib

max_forks_repo_name: sthagen/stdlib


text:
export Hexagonal

"""
    Hexagonal{Shape, BC} <: BoundedLattice{2, BC}

Hexagonal lattice.

Ref: https://en.wikipedia.org/wiki/Hexagonal_lattice
"""
struct Hexagonal{Shape <: Tuple, BC} <: BoundedLattice{2, BC} end

Hexagonal(x, y; boundary=Periodic) = Hexagonal{Tuple{x, y}, boundary}()
nameof(::Hexagonal) = "Hexagonal Lattice"


text:
module PrimeFactorsKata

function factorsOf(x)
    factors = Int[]
    divisor = 2
    while x > 1
        while x % divisor == 0
            append!(factors, divisor)
            x = convert(Int, x / divisor)
        end
        divisor += 1
    end
    return factors
end

export factorsOf

end # module


text:
# construct the values of 1/Jacobian at each LGL point for a given element
# Note: This is agnostic as to whether or not the elemnt is curved
# TODO: this could be adapted to reuse the routine already persent in CurvedMesh
function calc_inverse_jacobian!(inverse_jacobian::AbstractArray{<:Any, 3}, element, X_xi, X_eta,
                                Y_xi, Y_eta)
  @. @views inverse_jacobian[:, :, element] = inv(X_xi[:, :, element] * Y_eta[:, :, element] -
                                                  X_eta[:, :, element] * Y_xi[:, :, element])
  return inverse_jacobian
end

include("mappings_geometry_curved_2d.jl")
include("mappings_geometry_straight_2d.jl")