text:
export TimeSeriesLogger
"""

A struct for keeping track of numerous pieces of data over time throughout
a process. It's designed to be easy to use and relatively quickly. It
logs numeric types of any size (scalars, vectors, or matrices), as long
as the size is consistent from sample to sample.

Example:

    log = TimeSeriesLogger();    # Make a new logger.
    x   = [0; 0];
    for t = 1:100                # Make a random walk.
        x = x + randn(2, 1);     # Take a single step.
        add!(log, "walk", t, x); # Log a single sample, x, at time t.
    end

    tsl_plot(log)                # Plot everything.

    # We can also access specific logs by their names. In the above, we only
    #   have one log ('walk'). Let's get that log from the logger.

    t_log, x_log = get_log(log, "walk");    # Get a specific log.
    plot(x_log[:, 1], x_log[:, 2]);         # Do something with it.

    # We can make sure a log exists before trying to do anything with it:

    if contains(log,"walk")
        x_log = get_log(log, "walk");
        plot(x_log[:, 1], x_log[:, 2]);
    ;

    # If we want to log something but don't want it plotted when we call
    #   |plot|, then we can pass in false to the logger when we add the signal.

    add!(log, "var1", t, x, false);

    # To group items into different figures when they're plotted, give them
    #   "groups":

    add!(log, "var1", t, x, true, "group1");
    add!(log, "var2", t, y, true, "group2");
    add!(log, "foo",  t, bar, true, "group1");

    # All of the items with common group names will be added to the same
    #   figures as subplots.

    # Finally, we can clear out the logger with |initialize|. This deletes all
    #   data and returns the logger to its initial state.

    initialize!(log);

# Methods
-------

* initialize! Clear out all data.
* add!        Add a single data point to a time series.
* contains    See if a time series exists (by name).
* get_log     Return a log by name.
* tsl_plot    Plot the logged series.

# Notes

TimeSeriesLogger was created as the logging mechanism for odehybrid.
However, it's not dependent on that project and so can be used anywhere.

TimeSeriesLogger is not related to the 'timeseries' class in MATLAB.

Since this class never knows how much more data is going to be logged, it
can't preallocate the appropriate amount of space. However, increasing
the size of its data stores by one row every time a new sample is added
is very slow. To combat this, the data store starts off small and doubles
whenever the store is at capacity. Say we're logging a 4x1 signal. When
the first sample is added (say it's log k), data{k}{2} will be 1x4. When
the second signal is added, it becomes 2x4. For the third, it's 4x4, then
8x4, 16x4, etc. A separate counter stores how much of this allocated
space is currently used. This reduces the number of allocations from n to
log2(n). Practically, it saves a little time during logging without too
much complexity.

Online doc: http://www.anuncommonlab.com/doc/odehybrid/TimeSeriesLogger.html

Copyright 2014 An Uncommon Lab
"""
Base.@kwdef mutable struct TimeSeriesLogger

    # Properties
    names = []   # Unique identifier (within the group)
    times = []   # Time stamp corresponding to each data point
    data  = []   # Data, stored as {t, x}
    show  = []   # True iff the data should be plotted
    sizes = []   # Amount of allocated space
    count = []   # Amount of allocated space currently used
    group = []   # Plot group identifier

end;

export add!
"""
Add a new log or append to an existing log by name. Returns true iff a new log was created.

    add!(log, "var1", t, x);
    add!(log, "var1", t, x, true);  # Show log in plots (default)
    add!(log, "var1", t, x, false); # Don't show log in plots

    # Signals can be grouped together into figures by given them
    #   a common group argument. Here, both var1 and var2 logs will
    #   be plotted together.

    add!(log, "var1", t,  x, true, "group1");
    add!(log, "var2", t2, y, true, "group1");

NOTE that unlike in MATLAB, the signals are stored separately in times (t)
and data (x), (e.g., the kth pair would be log.times[k] and log.data[k])
"""
function add!(tsl::TimeSeriesLogger, name, t, x, show_it = true, group = "")

    # See if we've already added this variable
    index = isempty(tsl.names) ? [] : (tsl.names .== name)
    log_is_empty = (index == []) || !maximum(index)

    # NEW
    if !isa(x, Array)
        x = [x]
    end

    if log_is_empty
        # Initialize everything first time through

        new_times = Array{Float64, 1}()
        push!(new_times, t)
        push!(tsl.times, new_times)

        new_data = Array{Any, 1}()
        push!(new_data, x[:])
        push!(tsl.data, new_data)
        push!(tsl.names, name)

        push!(tsl.show, show_it)

        push!(tsl.sizes, 1)
        push!(tsl.count, 1)
        push!(tsl.group, group)

    else # otherwise, append data to existing set

        # If we have used up all the reserved space, double it
        c = tsl.count[index][1] + 1

        if c > tsl.sizes[index][1]
            push!(tsl.times[index][1], t)
            for i = 1:(tsl.sizes[index][1] - 1)
                push!(tsl.times[index][1], 0.0)
            end

            push!(tsl.data[index][1], x[:])
            for i = 1:(tsl.sizes[index][1] - 1)
                push!(tsl.data[index][1], zeros(maximum(size(x))))
            end

            tsl.sizes[index]   .= 2 * tsl.sizes[index][1]
        else # If space remains, add to it
            tsl.times[index][1][c]   = t
            tsl.data[index][1][c, :][1] .= x[:]
        end

        tsl.count[index] .= c

        # We should show this data if the user has ever asked to plot it
        tsl.show[index]  .= show_it || any(tsl.show[index])
        tsl.group[index] .= group

    end

    return log_is_empty
end;

export contains
"""
Return 'true' iff the TimeSeriesLogger contains this name
"""
function contains(tsl::TimeSeriesLogger, name)

    indices = findall(tsl.names .== name)
    return !isempty(indices)
end;

export get_log
"""
Return a specific log by name.

If one output is request, it returns the data. If two are
requested, it returns the time and the data. Returns empty if
the logs don't exist (use the 'contains' function to test for
this).

Example:

    x = get_log(log, "signal1");
    t, x = get_log(log, "signal1");

NOTE that |t| will be ns-by-1 and x will be ns-by-nx, where
nx is the number of elements in a single sample.
"""
function get_log(tsl::TimeSeriesLogger, name)

    index = (tsl.names .== name)

    if isempty(index)
        return ([], [])
    else
        t = tsl.times[index][1][1:tsl.count[index][1]]
        x = vec_to_mat(tsl.data[index][1][1:tsl.count[index][1]])
    end

    return (t, x)
end;

export initialize!
"""
Clears out all fields in a Time Series Logger (in-place)
"""
function initialize!(tsl::TimeSeriesLogger)

    tsl.names = [];
    tsl.times = [];
    tsl.data  = [];
    tsl.show  = [];
    tsl.sizes = [];
    tsl.count = [];
    tsl.group = [];

    return Nothing
end;

export tsl_plot
"""
Plot all of the signals, grouped appropriately into figures.
A custom x label can be added to figures as well
"""
function tsl_plot(tsl::TimeSeriesLogger, x_label)

    groups = unique(tsl.group)

    # For each group
    for g = 1:maximum(size(groups))

        # See what is to be shown and count them
        to_show = (tsl.group .== groups[g]) .& tsl.show[:]
        n = sum(to_show) #

        # Bail if there is no need to show anything
        if n == 0;
            continue;
        end

        plts = Array{Plots.Plot{Plots.GRBackend}, 1}(undef, n)

        ### TODO    ###############
        # Set up the figure
        # unique_figure("LoggerDefaultPlot$g", [groups[g] " Log"])

        # % Set up the figure.
        # unique_figure(sprintf('LoggerDefaultPlot%d', g), ...
        #               'Name', [groups{g} ' Log']);

        # clf();
        ###########################

        to_show = Int.(to_show)
        plt_idx = 1
        for k = 1:maximum(size(to_show))
            if (to_show[k] == 0)
                continue
            end

            N = tsl.count[k] #to_show[k]]
            plot_times = tsl.times[k][1:N]
            plot_data  = vec_to_mat(tsl.data[k][1:N, :])
            plot(plot_times, plot_data)

            ylabel!(tsl.names[k]) #to_show[k]])
            plts[plt_idx] = xlabel!(x_label)
            plt_idx += 1
        end
        display(plot(plts..., layout = (n, 1)))
    end
end;

function tsl_plot(tsl::TimeSeriesLogger)
    # Wrapper function that plots a TimeSeriesLogger without a label
    return tsl_plot(tsl, " ")
end;

# # NOT CURRENTLY FINISHED
# function unique_figure(id, name = " ", args = [])
#     """
#         h = unique_figure(id, varargin)

#         Creates a figure with the given ID (tag) or selects it if it already
#         exists, allowing one to easily reuse the same figure window identified
#         with text instead of handles. This is useful when, e.g., running a script
#         many times after clearing between runs without having to hard-code figure
#         numbers (which can become hard to keep track of).

#         Any additional arguments are passed along to the figure's |set| method.

#         Example:

#         h = unique_figure('trajectory');

#         Example with extra arguments:

#         h = unique_figure('trajectory', 'Name', 'Trajectory');
#     """
#     # See if there is a figure with this ID already
#     # h = findobj('Type', 'figure', 'Tag', id);
#     # if isempty(h)
#     h = plot("Tag", id, args);

#     # else

#     #     % If we found it, select it.
#     #     figure(h);

#     #     % If there were any additional arguments, pass them along. Note
#     #     % that if there weren't and we called set(h), then it would print h
#     #     % to the command window -- probably not what we want, hence the
#     #     % condition.
#     #     if nargin > 1
#     #         set(h, varargin{:});
#     #     end
#     # end

#     return h
# end;


text:
using PressureSensitiveMats
using Documenter

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

makedocs(;
    modules=[PressureSensitiveMats],
    authors="Carter Green <carter.green@icloud.com> and contributors",
    repo="https://github.com/carterjgreen/PressureSensitiveMats.jl/blob/{commit}{path}#{line}",
    sitename="PressureSensitiveMats.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://carterjgreen.github.io/PressureSensitiveMats.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "combiners.md",
        "movement.md",
        "occupancy.md",
        "utils.md"
    ],
)

deploydocs(;
    repo="github.com/carterjgreen/PressureSensitiveMats.jl",
    devbranch = "main",
)

ext: jl


text:
module Surrogates
using LinearAlgebra
using GLM
using Distributions
using Sobol
using LatinHypercubeSampling
using Stheno

abstract type AbstractSurrogate <: Function end
include("utils.jl")
include("Radials.jl")
include("Kriging.jl")
include("Sampling.jl")
include("Optimization.jl")
include("Lobachesky.jl")
include("LinearSurrogate.jl")
include("InverseDistanceSurrogate.jl")
include("SecondOrderPolynomialSurrogate.jl")
include("SthenoKriging.jl")
include("RandomForestSurrogate.jl")
include("NeuralSurrogate.jl")

export AbstractSurrogate, SamplingAlgorithm
export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
export RandomSample, KroneckerSample, GoldenSample
export SRBF,LCBS,EI,DYCORS,SOP,surrogate_optimize
export LobacheskySurrogate, lobachesky_integral, lobachesky_integrate_dimension
export LinearSurrogate
export RandomForestSurrogate
export SVMSurrogate
export NeuralSurrogate
export InverseDistanceSurrogate
export SecondOrderPolynomialSurrogate
export SthenoKriging

end


text:
function knot_set(directory::AbstractString,
                  knot_set::Matrix)
    open("$(directory)/knot.dat", "w") do file
        for i = 1:size(knot_set,1)
            write(file, string(knot_set[i,2]), " ==> ", knot_set[i,1])
            write(file, '\n')
        end
        write(file, "***\n")
    end
end

function knot_set(directory::AbstractString,
                  order::Integer,
                  num_splines::Integer,
                  Z::Real,
                  E_max::Real,
                  r_max::Real,
                  h_max::Real = 0,
                  h_int::Real = 0)
    h_max == 0 && (h_max = min(sqrt(1/E_max), r_max)) # Only valid for order == 8 ?
    n = floor(log(8Z)/log(2))
    h_int == 0 && (h_int = 1/2.0^n)
    knots = ["order of splines (ks)" order
             "number of splines (ns)" num_splines
             "nuclear charge (z)" Z
             "step size from 0 to 1 (h for z*r, = (1/2)^n)" h_int
             "maximum step size (hmax for r)" h_max
             "maximum r (rmax)" r_max]
    knot_set(directory, knots)
    knots
end

function load_knot_set(directory::AbstractString)
    open(joinpath(directory, "knot.dat")) do file
        read_parm = T -> parse(T, split(readline(file), "==>")[1])
        knot_set = Dict()
        knot_set[:ks] = read_parm(Int)
        knot_set[:ns] = read_parm(Int)
        knot_set[:Z] = read_parm(Float64)
        knot_set[:h] = read_parm(Float64)
        knot_set[:hmax] = read_parm(Float64)
        knot_set[:rmax] = read_parm(Float64)
        readline(file)
        knot_set[:t] = Vector{Float64}()
        while length(knot_set[:t]) < knot_set[:ns] + knot_set[:ks]
            tt = map(s -> parse(Float64, s), split(readline(file)))
            append!(knot_set[:t], tt)
        end
        readline(file)
        knot_set[:n_inner_knots] = read_parm(Int)
        knot_set[:n_exp_knots] = read_parm(Int)
        knot_set[:n_intervals] = read_parm(Int)
        knot_set[:max_k2_Ry] = read_parm(Float64)
        knot_set
    end
end


text:
abstract type AbstractMapMod <: ArrayModifier end

struct ReadonlyMapMod{F} <: AbstractMapMod
    f   ::  F
end
struct MapMod{F,Finv} <: AbstractMapMod
    f   ::  F
    finv::  Finv
end

ModStyle(::AbstractMapMod, ::Interface) = ModNothing()
ModStyle(::Type{<:AbstractMapMod}, ::IF_IndexStyle) = ModNothing()

ModStyle(::AbstractMapMod, ::IF_eltype) = ModRecursive()
mod_eltype_post(mod::AbstractMapMod, T) = Base.promote_op(mod.f,T)

ModStyle(::AbstractMapMod, ::IF_getindex) = ModRecursive()
mod_getindex_post(mod::AbstractMapMod, Z, I...) = mod.f(Z)

ModStyle(::MapMod, ::IF_setindex!) = ModRecursive()
mod_setindex!_pre(mod::MapMod, val, I...) = (mod.finv(val), I...)

ModStyle(::ReadonlyMapMod, ::IF_setindex!) = ModFinal()
mod_setindex_final(mod::ReadonlyMapMod, I...) = throw(InexactError("Readonly map is not writeable!"))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Location <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    subscriptionId::Any # spec type: Union{ Nothing, String } # spec name: subscriptionId
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    regionalDisplayName::Any # spec type: Union{ Nothing, String } # spec name: regionalDisplayName
    metadata::Any # spec type: Union{ Nothing, LocationMetadata } # spec name: metadata

    function Location(;id=nothing, subscriptionId=nothing, name=nothing, displayName=nothing, regionalDisplayName=nothing, metadata=nothing)
        o = new()
        validate_property(Location, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(Location, Symbol("subscriptionId"), subscriptionId)
        setfield!(o, Symbol("subscriptionId"), subscriptionId)
        validate_property(Location, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(Location, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(Location, Symbol("regionalDisplayName"), regionalDisplayName)
        setfield!(o, Symbol("regionalDisplayName"), regionalDisplayName)
        validate_property(Location, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type Location

const _property_map_Location = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("subscriptionId")=>Symbol("subscriptionId"), Symbol("name")=>Symbol("name"), Symbol("displayName")=>Symbol("displayName"), Symbol("regionalDisplayName")=>Symbol("regionalDisplayName"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_Location = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("subscriptionId")=>"String", Symbol("name")=>"String", Symbol("displayName")=>"String", Symbol("regionalDisplayName")=>"String", Symbol("metadata")=>"LocationMetadata")
Base.propertynames(::Type{ Location }) = collect(keys(_property_map_Location))
Swagger.property_type(::Type{ Location }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Location[name]))}
Swagger.field_name(::Type{ Location }, property_name::Symbol) =  _property_map_Location[property_name]

function check_required(o::Location)
    true
end

function validate_property(::Type{ Location }, name::Symbol, val)
end


text:
using Base64: base64encode
using SodiumSeal: KeyPair, seal, unseal
using Test

const K = KeyPair()
const M = rand(UInt8, 100)

@testset "SodiumSeal.jl" begin
    @testset "Bytes" begin
        ciphertext = seal(M, K)
        decrypted = unseal(ciphertext, K)
        @test decrypted == M
    end

    @testset "Strings" begin
        ciphertext = seal(base64encode(M), K)
        decrypted = unseal(ciphertext, K)
        @test decrypted == base64encode(M)
    end
end


text:
using Base.Test
using SparseVectors

let S = sprand(4, 8, 0.5)
    Sf = full(S)
    @assert isa(Sf, Matrix{Float64})

    # get a single column
    for j = 1:size(S,2)
        col = getcol(S, j)
        @test isa(col, SparseVector{Float64,Int})
        @test length(col) == size(S,1)
        @test full(col) == Sf[:,j]
    end

    # column views
    for j = 1:size(S,2)
        col = view(S, :, j)
        @test isa(col, SparseVectorView{Float64,Int})
        @test length(col) == size(S,1)
        @test full(col) == Sf[:,j]
    end

    # column-range views

    # non-empty range
    V = unsafe_colrange(S, 2:6)
    @test isa(V, SparseMatrixCSC{Float64,Int})
    @test size(V) == (4, 5)
    @test full(V) == Sf[:, 2:6]
    @test !isempty(V)

    # empty range
    V0 = unsafe_colrange(S, 2:1)
    @test isa(V0, SparseMatrixCSC{Float64,Int})
    @test size(V0) == (4, 0)
    @test isempty(V0)

end


text:
function snapshot!()
    ctx = Pkg.Types.Context();
    proj = ctx.env.project_file;
    man = ctx.env.manifest_file;

    cp(proj, joinpath(dirname(proj), "project_snapshot.toml"), force = true);
    cp(man, joinpath(dirname(proj), "manifest_snapshot.toml"), force = true);
end

function undo!()
    ctx = Pkg.Types.Context();
    snap_proj = joinpath(dirname(ctx.env.project_file), "project_snapshot.toml")
    snap_man = joinpath(dirname(ctx.env.manifest_file), "manifest_snapshot.toml")
    if isfile(snap_proj) && isfile(snap_man)
        cp(snap_proj, joinpath(dirname(snap_proj), "Project.toml"), force = true);
        cp(snap_man, joinpath(dirname(snap_proj), "Manifest.toml"), force = true);
        Pkg.activate(dirname(snap_proj))
    else
        @info "No valid snapshot found."
    end
end


text:
include("abm.jl")
using WebIO, CSSUtil, Interact, Mux
using Agents, DataFrames
using Plots

infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)

function generate_ui()

    function simulation_start(ob)

        stop_simulation = false
        println("starting simulation ...")

        params = COVID.create_params(C=8, max_travel_rate=0.01)
        model = COVID.model_initiation(;params...)
        N = sum(model.Ns) # Total initial population

        to_collect = [(:status, f) for f in (infected, recovered, length)]
        data = init_agent_dataframe(model, to_collect)

        nsteps = parse(Int, steps_widget[:steps_textbox][])

        @async for step in 1:1:nsteps
            try
                if stop_simulation
                    stop_simulation = false
                    println("stopping simulation")
                    break
                else
                    sleep(0.001)
                    #println("simulation step $step")
                    step!(model, COVID.agent_step!)
                    collect_agent_data!(data, model, to_collect, step)

                    # update plot
                    x = data.step
                    y_infected = log10.(data[:, Symbol("infected(status)")])
                    y_recovered = log10.(data[:, Symbol("recovered(status)")])
                    y_dead = log10.(N .- data[:, Symbol("length(status)")])

                    p = plot(x, [y_infected, y_recovered, y_dead], label=["infected" "recovered" "dead"])
                    xlabel!(p, "steps")
                    ylabel!(p, "log( count )")

                    if step == nsteps
                        simulation_complete = true
                        println("simulation complete")
                    end

                    plt[] = p
                end
            catch
                for (exc, bt) in Base.catch_stack()
                    showerror(stdout, exc, bt)
                    println()
                end

                stop_simulation = true
                break
            end
        end
        plt[] = plt[]
    end

    function simulation_stop(ob)
        stop_simulation = true
        return plt[]
    end

    function simulation_reset(ob)
        println("resetting simulation")
        plot()
    end

    # default parameters
    default_steps = 100

    # control
    stop_simulation = false
    simulation_complete = false

    # ui components
    html_title = HTML(string("
SIR Model based on the  Julia Agents.jl implementation
"))

    # TODO marco or function to generate labelled fields
    # steps field
    dict_steps_widget = OrderedDict(:steps_label => "Steps", :steps_textbox => textbox("Steps", value=string(default_steps)))
    #dict_steps_widget[:steps_textbox] = style(dict_steps_widget[:steps_textbox], Pair(:width, 5em))
    steps_widget = Widget{:mywidget}(dict_steps_widget)
    @layout! steps_widget hbox(:steps_label, CSSUtil.hskip(1em), :steps_textbox);

    start_button = button("Start")
    stop_button = button("Stop")
    reset_button = button("Reset")
    plt = Observable{Any}(plot())

    map!(simulation_start, plt, start_button)
    map!(simulation_stop, plt, stop_button)
    map!(simulation_reset, plt, reset_button)

    ui = vbox( # put things one on top of the other
        #node(:div, "hello"),
        html_title,
        steps_widget,
        hbox( # put things one next to the other
            pad(1em, start_button), # to allow some white space around the widget
            pad(1em, stop_button),
            pad(1em, reset_button),
        ),
        plt
    );

    function test(ob)
        try
            nsteps = parse(Int, steps_widget[:steps_textbox][])
            println("nsteps is $nsteps")
        catch
            for (exc, bt) in Base.catch_stack()
                showerror(stdout, exc, bt)
                println()
            end
        end
    end

    return ui
end

#display(ui)

# if we are not in Juno IDE then assume webserver environ
if @isdefined(Juno)
    display(ui)
else
    port = 8000
    println("server started on port $port ...")
    fetch(WebIO.webio_serve(page("/", req -> generate_ui()), port))
end


text:
function anagram(phrases)
    @show length(phrases)
    trues = []
    lower_nospace_phrases = copy(phrases)
    for i = 1 : length(phrases)
        lower_nospace_phrases[i] = join(sort(split(lowercase(replace(lower_nospace_phrases[i]," " => "")),"")))
    end

    i = 1
    j = i + 1

    while i <= length(phrases)-1
        while j <= length(phrases)
            if lower_nospace_phrases[i] == lower_nospace_phrases[j]
                if length(phrases)>1
                    deleteat!(phrases,j)
                    j = j - 1
                end
            end
            j = j + 1
        end
        i = i + 1
        j = i + 1
    end

    return phrases
end

@time(println(anagram(phrases)))

max_stars_repo_name: atyiadam/anagram-competition

max_issues_repo_name: atyiadam/anagram-competition

max_forks_repo_name: atyiadam/anagram-competition


text:
struct NaNChecker{F}
    fields :: F
end

"""
    NaNChecker(; fields)

Return a `NaNChecker`, which sets `sim.running=false` if a `NaN` is detected
in any member of `fields` when `NaNChecker(sim)` is called. `fields` should be
a container with key-value pairs like a dictionary or `NamedTuple`.
"""
NaNChecker(; fields) = NaNChecker(fields)

hasnan(field) = any(isnan, parent(field))

function (nc::NaNChecker)(simulation)
    for (name, field) in pairs(nc.fields)
        if hasnan(field)
            @info "NaN found in field $name. Stopping simulation."
            simulation.running = false
        end
    end
    return nothing
end


text:
"""
    SSIM([kernel], [(α, β, γ)]) <: FullReferenceIQI
    assess(iqi::SSIM, img, ref)
    ssim(img, ref)

Structural similarity (SSIM) index is an image quality assessment method based
on degradation of structural information.

The SSIM index is composed of three components: luminance, contrast, and
structure; `ssim = 𝐿ᵅ * 𝐶ᵝ * 𝑆ᵞ`, where `W := (α, β, γ)` controls relative
importance of each components. By default `W = (1.0, 1.0, 1.0)`.

In practice, a mean version SSIM is used. At each pixel, SSIM is calculated
locally with neighborhoods weighted by `kernel`, returning a ssim map;
`ssim` is actaully `mean(ssim_map)`.
By default `kernel = KernelFactors.gaussian(1.5, 11)`.

!!! info

    SSIM is defined only for gray images. RGB images are treated as 3d Gray
    images. General `Color3` images are converted to RGB images first, in which
    case, you could manually expand them using `channelview` if you don't want
    them converted to RGB first.

# Example

`ssim(img, ref)` should be sufficient to get a benchmark for algorithms. One
could also instantiate a customed SSIM, then pass it to `assess` or use it as a
function. For example:

```julia
iqi = SSIM(KernelFactors.gaussian(2.5, 17), (1.0, 1.0, 2.0))
assess(iqi, img, ref)
iqi(img, ref)
```

# Reference

[1] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2004). Image quality assessment: from error visibility to structural similarity. _IEEE transactions on image processing_, 13(4), 600-612.

[2] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2003). The SSIM Index for Image Quality Assessment. Retrived May 30, 2019, from http://www.cns.nyu.edu/~lcv/ssim/
"""
struct SSIM <: FullReferenceIQI
    kernel::AbstractArray{<:Real}
    W::NTuple{3}
    function SSIM(kernel, W)
        ndims(kernel) == 1 || throw(ArgumentError("only 1-d kernel is valid"))
        issymetric(kernel) || @warn "SSIM kernel is assumed to be symmetric"
        all(W .>= 0) || throw(ArgumentError("(α, β, γ) should be non-negative, instead it's $(W)"))
        new(kernel, W)
    end
end

# default values from [1]
const SSIM_KERNEL = KernelFactors.gaussian(1.5, 11) # kernel
const SSIM_W = (1.0, 1.0, 1.0) # (α, β, γ)
SSIM(kernel=SSIM_KERNEL) = SSIM(kernel, SSIM_W)

# api
# By default we don't crop the padding boundary to meet the ssim result from
# MATLAB Image Processing Toolbox, which is used more broadly than the original
# implementaion [2] (written in MATLAB as well).
# TODO: add keyword argument "crop=false" for compatibility
# -- Johnny Chen <johnnychen94@hotmail.com>
(iqi::SSIM)(x, ref) = mean(_ssim_map(iqi, x, ref))

@doc (@doc SSIM)
ssim(x, ref) = SSIM()(x, ref)

# Parameters `(K₁, K₂)` are used to avoid instability when denominator is very
# close to zero. Different from origianl implementation [2], we don't make it
# public since the ssim result is insensitive to these parameters according to
# [1].
# -- Johnny Chen <johnnychen94@hotmail.com>
const SSIM_K = (0.01, 0.03)

# SSIM is defined only for gray images,
# RGB images are treated as 3d gray images,
# other Color3 images are converted to RGB first.
function _ssim_map(iqi::SSIM, x::GenericGrayImage, ref::GenericGrayImage, peakval = 1.0, K = SSIM_K)
    if size(x) ≠ size(ref)
        err = ArgumentError("images should be the same size, instead they're $(size(x))-$(size(ref))")
        throw(err)
    end
    α, β, γ = iqi.W
    C₁, C₂ = @. (peakval * K)^2
    C₃ = C₂/2

    T = promote_type(float(eltype(ref)), float(eltype(x)))
    x = of_eltype(T, x)
    ref = of_eltype(T, ref)

    # calculate ssim in the neighborhood of each pixel, weighted by window
    window = kernelfactors(Tuple(repeated(iqi.kernel, ndims(ref))))

    μx = imfilter(x, window)   # equation (14) in [1]
    μy = imfilter(ref, window) # equation (14) in [1]
    μx² = μx .* μx
    μy² = μy .* μy
    μxy = μx .* μy
    σx² = imfilter(x.^2, window) .- μx²     # equation (15) in [1]
    σy² = imfilter(ref.^2, window) .- μy²   # equation (15) in [1]
    σxy = imfilter(x .* ref, window) .- μxy # equation (16) in [1]

    if [α, β, γ] ≈ [1.0, 1.0, 1.0]
        # equation (13) in [1]
        ssim_map = @. ((2μxy + C₁)*(2σxy + C₂))/((μx²+μy²+C₁)*(σx² + σy² + C₂))
    else
        σx_σy = @. sqrt(σx²*σy²)
        l = @. (2μxy + C₁)/(μx² + μy²) # equation (6) in [1]
        c = @. (2σx_σy + C₂)/(σx² + σy² + C₂) # equation (9) in [1]
        s = @. (σxy + C₃)/(σx_σy + C₃) # equation (10) in [1]

        ssim_map = @. l^α * c^β * s^γ # equation (12) in [1]
    end
    return ssim_map
end

_ssim_map(iqi::SSIM,
          x::AbstractArray{<:AbstractRGB},
          ref::AbstractArray{<:AbstractRGB},
          peakval = 1.0, K = SSIM_K) =
    _ssim_map(iqi, channelview(x), channelview(ref), peakval, K)

_ssim_map(iqi::SSIM,
          x::AbstractArray{<:Color3},
          ref::AbstractArray{<:Color3},
          peakval = 1.0, K = SSIM_K) =
    _ssim_map(iqi, of_eltype(RGB, x), of_eltype(RGB, ref), peakval, K)

# helpers
function issymetric(kernel)
    origin = first(axes(kernel, 1))
    center = (length(kernel)-1) ÷ 2 + origin
    kernel[origin:center] ≈ kernel[end:-1:center]
end


text:
#===============================================================================

The radical of n, rad(n), is the product of the distinct prime factors of n. For
example, 504 = 2^3 × 3^2 × 7, so rad(504) = 2 × 3 × 7 = 42.

If we calculate rad(n) for 1 ≤ n ≤ 10, then sort them on rad(n), and sorting on
n if the radical values are equal, we get:

 Unsorted |   Sorted  |
 n rad(n) |  n rad(n) |  k
 ---------+-----------+---
 1   1    |  1   1    |  1
 2   2    |  2   2    |  2
 3   3    |  4   2    |  3
 4   2    |  8   2    |  4
 5   5    |  3   3    |  5
 6   6    |  9   3    |  6
 7   7    |  5   5    |  7
 8   2    |  6   6    |  8
 9   3    |  7   7    |  9
10  10    | 10  10    | 10

Let E(k) be the kth element in the sorted n column; for example, E(4) = 8 and
E(6) = 9.

If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).

===============================================================================#

using Primes
using Base.Test

function compute_table(limit)
    res = Tuple{Int,Int}[]
    for n in 1:limit
        push!(res, (prod(keys(factor(n))), n))
    end

    sort!(res)
    map(pair -> pair[2], res)
end

test_table = compute_table(10)

@test test_table[4] == 8
@test test_table[6] == 9

actual_table = compute_table(100_000)

println("Euler 124: $(actual_table[10_000])")

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
using NearestNeighbors, StaticArrays
using Distances: Euclidean, Metric
import NearestNeighbors: KDTree

export AbstractNeighborhood
export FixedMassNeighborhood, FixedSizeNeighborhood
export neighborhood, KDTree

#####################################################################################
#                              Neighborhoods n stuff                                #
#####################################################################################
"""
    AbstractNeighborhood
Supertype of methods for deciding the neighborhood of points for a given point.

Concrete subtypes:
* `FixedMassNeighborhood(K::Int)` :
  The neighborhood of a point consists of the `K`
  nearest neighbors of the point.
* `FixedSizeNeighborhood(ε::Real)` :
  The neighborhood of a point consists of all
  neighbors that have distance < `ε` from the point.

See [`neighborhood`](@ref) for more.
"""
abstract type AbstractNeighborhood end

struct FixedMassNeighborhood <: AbstractNeighborhood
    K::Int
end
FixedMassNeighborhood() = FixedMassNeighborhood(1)

struct FixedSizeNeighborhood <: AbstractNeighborhood
    ε::Float64
end
FixedSizeNeighborhood() = FixedSizeNeighborhood(0.01)

"""
    neighborhood(point, tree, ntype)
    neighborhood(point, tree, ntype, n::Int, w::Int = 1)

Return a vector of indices which are the neighborhood of `point` in some
`data`, where the `tree` was created using `tree = KDTree(data [, metric])`.
The `ntype` is the type of neighborhood and can be any subtype
of [`AbstractNeighborhood`](@ref).

Use the second method when the `point` belongs in the data,
i.e. `point = data[n]`. Then `w` stands for the Theiler window (positive integer).
Only points that have index
`abs(i - n) ≥ w` are returned as a neighborhood, to exclude close temporal neighbors.
The default `w=1` is the case of excluding the `point` itself.

## References

`neighborhood` simply interfaces the functions
`knn` and `inrange` from
[NearestNeighbors.jl](https://github.com/KristofferC/NearestNeighbors.jl) by using
the argument `ntype`.
"""
function neighborhood(point::AbstractVector, tree,
                      ntype::FixedMassNeighborhood, n::Int, w::Int = 1)
    idxs, = knn(tree, point, ntype.K, false, i -> abs(i-n) < w)
    return idxs
end
function neighborhood(point::AbstractVector, tree, ntype::FixedMassNeighborhood)
    idxs, = knn(tree, point, ntype.K, false)
    return idxs
end

function neighborhood(point::AbstractVector, tree,
                      ntype::FixedSizeNeighborhood, n::Int, w::Int = 1)
    idxs = inrange(tree, point, ntype.ε)
    filter!((el) -> abs(el - n) ≥ w, idxs)
    return idxs
end
function neighborhood(point::AbstractVector, tree, ntype::FixedSizeNeighborhood)
    idxs = inrange(tree, point, ntype.ε)
    return idxs
end

KDTree(D::AbstractDataset, metric::Metric = Euclidean()) = KDTree(D.data, metric)


text:
@testset "Particles" begin
    T = transportoperator(grd, 100.0)
    @test T isa SparseMatrixCSC
    @testset "vᵀ T = 0 (T conserves mass)" begin
        @test norm(v) / norm(T' * v) > ustrip(upreferred(1u"Myr"))
    end

end


text:
using SnakeAI

const WIDTH = SnakeAI.WIDTH
const HEIGHT = SnakeAI.HEIGHT
const BACKGROUND = colorant"black"

game = SnakeAI.Game()
agent = SnakeAgent()

function draw()
    # Drawing the snake
    for p in game.snake.body
        outside = Rect(p.x, p.y, SnakeAI.BLOCK_SIZE, SnakeAI.BLOCK_SIZE)
        inside = Rect(p.x + 2, p.y + 2, SnakeAI.BLOCK_SIZE - 4, SnakeAI.BLOCK_SIZE - 4)
        draw(outside, colorant"#00262a", fill=true)
        draw(inside, colorant"#00717f", fill=true)

    end

    # Drawing the food
    f = Rect(game.food.x, game.food.y, SnakeAI.BLOCK_SIZE, SnakeAI.BLOCK_SIZE)
    draw(f, colorant"#7F0071", fill=true)
end

function update(g::GameZero.Game)
    global agent, game

    train!(agent, game)
end


text:
using ParallelAccelerator

@acc function kmeans()
    numCenter = 5
    iterNum = 30
    points = rand(20,10000)
    D,N = size(points) # number of features, instances
    centroids = rand(D, numCenter)

    for l in 1:iterNum
        dist :: Array{Array{Float64,1},1} = [ Float64[sqrt(sum((points[:,i].-centroids[:,j]).^2)) for j in 1:numCenter] for i in 1:N]
        labels :: Array{Int, 1} = [indmin(dist[i]) for i in 1:N]
        centroids :: Array{Float64,2} = [ sum(points[j,labels.==i])/sum(labels.==i) for j in 1:D, i in 1:numCenter]
    end
    return centroids
end

function main()

    tic()
    centroids_out = kmeans()
    time = toq()
    println("result = ", centroids_out)
    println("SELFTIMED ", time)

    #benchmark()

end

function benchmark()
  maxt = 0
  mint = Inf
  totalt = 0
  for i in 1:100
    tic()
    centroids_out = kmeans()
    t = toq()
    println("finished ",i)
    if t > maxt
      maxt = t
    end
    if t < mint
      mint = t
    end
    totalt += t
  end
  println("SELFTIMED AVERAGE ", totalt/100.0)
  println("SELFTIMED MINIMUM", mint)
  println("SELFTIMED MAXIMUM", maxt)
end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
"""
    minimumAssignmentProblem(C)

Solve minimum assignment problem using vector of costs `C`.
"""
function minimumAssignmentProblem(C)
    @assert ndims(C) == 2 && size(C, 1) == size(C, 2)

    m = size(C, 1)
    model = Model(solver = CplexSolver(CPXPARAM_ScreenOutput = getProperty("minimumAssignmentProblem.cplexLogging")))
    @variable(model, x[1:m, 1:m], Bin)
    @objective(model, Min, vecdot(C, x))

    for constraint in getAssignmentConstraints(m)
        JuMP.addconstraint(model, constraint(x))
	end

    status = solve(model)
    return map((x) -> Int(x), getvalue(x))
end

"""
    getAssignmentConstraints(m)

Return a list of constraints defining a set of feasible solutions of a minimum assignment problem.
Each constraint is function with one parameter, which is variable of a mathematical programming model.
"""
function getAssignmentConstraints(m)
    constraints = []
    for j in 1:m
        push!(constraints, x-> @LinearConstraint(sum(x[i,j] for i=1:m)==1))
    end
    for i in 1:m
        push!(constraints, x -> @LinearConstraint(sum(x[i,j] for j=1:m)==1))
    end
    constraints
end


text:
# Converging-diverging nozzles
export Nozzle, areas, positions, converging, diverging, throat, nozzleexit,
       machnumber,temperature, density, massflowrate, flow_quality, NozzleProcess

abstract type NozzleQuality end
abstract type SubsonicIsentropicNozzle <: NozzleQuality end
abstract type SupersonicWithShockNozzle <: NozzleQuality end
abstract type SupersonicIsentropicNozzle <: NozzleQuality end

struct Nozzle{T,S}
    Ai :: T
    At :: T
    Ae :: T
    x :: S
    areas :: Vector{T}
    throat_index :: Integer
end

function Nozzle(Ai::Area,At::Area,Ae::Area;xmin=-1.0,xmax=3.0,len=201)
    x = range(xmin,xmax,length=len)
    a = areas(x,Ai,At,Ae)
    Nozzle(Ai,At,Ae,x,a,_find_throat(a,At))
end

function areas(x,Ai::Area,At::Area,Ae::Area)
    a = Array{Area,1}(undef,length(x))
    araw = _nozzle_area(x,value(Ai),value(At),value(Ae))
    a .= Area.(araw)
end

function _find_throat(a::AbstractVector{T},At::Area) where T <: Area
    indx = findall(x -> (value(x) ≈ value(At)),a)
    isempty(indx) && error("Can't find the throat")
    return indx[1]
end

areas(n::Nozzle) = n.areas
positions(n::Nozzle) = n.x
throat(n::Nozzle) = n.At
nozzleexit(n::Nozzle) = n.Ae
inlet(n::Nozzle) = n.Ai
converging(n::Nozzle) = view(areas(n),1:n.throat_index)
diverging(n::Nozzle) = view(areas(n),n.throat_index+1:length(areas(n)))

function Base.show(io::IO, noz::Nozzle)
  println(io, "Converging-diverging nozzle")
  println(io, "   Inlet area (sq cm)= "*string(value(inlet(noz),u"cm^2")))
  println(io, "   Throat area (sq cm)= "*string(value(throat(noz),u"cm^2")))
  println(io, "   Exit area (sq cm)= "*string(value(nozzleexit(noz),u"cm^2")))
end

# Shaping of nozzle cross-section
_rad_conv(x,Rt,Ri) = Rt + (Ri-Rt)*x^2
_rad_div(x,Rt,Re,c) = Re + (Rt-Re)*exp(-c*x^2)
_radius_shape(x,Ri,Rt,Re,c) = x <= 0 ? _rad_conv(x,Rt,Ri) : _rad_div(x,Rt,Re,c)
function _nozzle_area(x,Ai,At,Ae;c=1)
    return π*_radius_shape.(x,sqrt(Ai/π),sqrt(At/π),sqrt(Ae/π),c).^2
end

struct NozzleProcess{process_type,T,S,G}
  noz :: Nozzle{T,S}
  p :: Vector{Pressure}
  M :: Vector{MachNumber}
  p0 :: StagnationPressure
  T0 :: StagnationTemperature
  pb_subcrit :: Pressure
  pb_supcrit_exitshock :: Pressure
  pb_supcrit :: Pressure
  pb :: Pressure
  As :: T
  gas :: G
  #NozzleProcess{process_type}(noz,p,M,p0,T0,pb,As,gas) = new{process_type,}
end

function NozzleProcess(noz::Nozzle{T,S},pb::Pressure{PU},p0::StagnationPressure{PSU},T0::StagnationTemperature{PTU};gas::PerfectGas=DefaultPerfectGas) where {T,S,PU,PSU,PTU}
  pb_subcrit = _pressure_subsonic_critical(noz,p0,gas)
  pb_supcrit_exitshock = _pressure_supersonic_shock_critical(noz,p0,gas)
  pb_supcrit = _pressure_supersonic_critical(noz,p0,gas)

  process_type = pb > pb_subcrit ? SubsonicIsentropicNozzle :
                (pb > pb_supcrit_exitshock ? SupersonicWithShockNozzle :
                                             SupersonicIsentropicNozzle)

  As = _shock_area(noz,pb,p0,process_type,gas)
  p = pressure_nozzle(noz,pb,p0,As,process_type,gas)
  M = machnumber_nozzle(noz,pb,p0,As,process_type,gas)

  NozzleProcess{process_type,T,S,typeof(gas)}(noz,p,M,p0,T0,
                  pb_subcrit,pb_supcrit_exitshock,pb_supcrit, pb,As,gas)
end

function Base.show(io::IO, nozproc::NozzleProcess)
  println(io, "Flow in a converging-diverging nozzle")
  println(io, "   Inlet area (sq cm)= "*string(value(inlet(nozproc.noz),u"cm^2")))
  println(io, "   Throat area (sq cm)= "*string(value(throat(nozproc.noz),u"cm^2")))
  println(io, "   Exit area (sq cm)= "*string(value(nozzleexit(nozproc.noz),u"cm^2")))
  println(io, "   Stagnation pressure (KPa) = "*string(value(nozproc.p0,u"kPa")))
  println(io, "   Stagnation temperature (K) = "*string(value(nozproc.T0,u"K")))
  println(io, "   Subsonic choked isentropic (KPa) = "*string(value(nozproc.pb_subcrit,u"kPa")))
  println(io, "   Supersonic choked isentropic (KPa) = "*string(value(nozproc.pb_supcrit,u"kPa")))
  println(io, "   Back pressure (KPa) = "*string(value(nozproc.pb,u"kPa")))
  println(io, "   Behavior --> "*flow_quality(nozproc))
end

machnumber(np::NozzleProcess) = np.M
pressure(np::NozzleProcess) = np.p

function temperature(np::NozzleProcess)
    Temp = Temperature[]
    for Mi in machnumber(np)
      T0_over_T = T0OverT(Mi,Isentropic,gas=np.gas)
      push!(Temp,Temperature(np.T0/T0_over_T))
    end
    Temp
end

function density(np::NozzleProcess)
    ρ0 = StagnationDensity(np.p0,np.T0,gas=np.gas)
    ρ = Density[]
    for (pi,Ti) in zip(pressure(np),temperature(np))
        push!(ρ,Density(pi,Ti,gas=np.gas))
    end
    ρ
end

function massflowrate(np::NozzleProcess)
    idx = np.noz.throat_index
    Mt = machnumber(np)[idx]
    pt = pressure(np)[idx]
    T0_over_T = T0OverT(Mt,Isentropic,gas=np.gas)
    Tt = Temperature(np.T0/T0_over_T)
    at = SoundSpeed(Tt,gas=np.gas)
    ut = Velocity(at,Mt)
    ρt = Density(pt,Tt,gas=np.gas)
    MassFlowRate(ρt,ut,throat(np.noz))
end

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

function machnumber_nozzle(noz::Nozzle,pb::Pressure,p01::StagnationPressure,As::Area,::Type{SupersonicWithShockNozzle},gas::PerfectGas)
    vcat(_machnumber_subsonic(converging(noz),converging(noz)[end],gas),
         _machnumber_diverging_nozzle_with_shock(diverging(noz),As,throat(noz),p01,gas))
end

function pressure_nozzle(noz::Nozzle,pb::Pressure,p01::StagnationPressure,As::Area,::Type{SupersonicWithShockNozzle},gas::PerfectGas)
    vcat(_pressure_subsonic(converging(noz),converging(noz)[end],p01,gas),
         _pressure_diverging_nozzle_with_shock(diverging(noz),As,throat(noz),p01,gas))
end

# ---

function machnumber_nozzle(noz::Nozzle,pb,p0,As,::Type{SupersonicIsentropicNozzle},gas::PerfectGas)
    vcat(_machnumber_subsonic(converging(noz),converging(noz)[end],gas),
         _machnumber_diverging_nozzle(diverging(noz),throat(noz),gas))
end

function pressure_nozzle(noz::Nozzle,pb,p0::StagnationPressure,As,::Type{SupersonicIsentropicNozzle},gas::PerfectGas)
    vcat(_pressure_subsonic(converging(noz),converging(noz)[end],p0,gas),
         _pressure_diverging_nozzle(diverging(noz),throat(noz),p0,gas))
end

# ---
machnumber_nozzle(noz::Nozzle,pb::Pressure,p0::StagnationPressure,As,::Type{SubsonicIsentropicNozzle},gas::PerfectGas) =
    _machnumber_subsonic(areas(noz),pb,p0,gas)

pressure_nozzle(noz::Nozzle,pb::Pressure,p0::StagnationPressure,As,::Type{SubsonicIsentropicNozzle},gas::PerfectGas) =
    _pressure_subsonic(areas(noz),pb,p0,gas)

# Subsonic nozzle processes

function _machnumber_subsonic(A::AbstractVector{T},pb::Pressure,p0::StagnationPressure,gas::PerfectGas) where T <: Area
  Mb = MachNumber(PressureRatio(pb/p0),Isentropic,gas=gas)
  Astar = AStar(A[end],Mb,gas=gas)
  _machnumber_subsonic(A,Astar,gas)
end

function _machnumber_subsonic(A::AbstractVector{T},Astar::Area,gas::PerfectGas) where T <: Area
  return SubsonicMachNumber(A,Astar,Isentropic,gas=gas)
end

function _pressure_subsonic(A::AbstractVector{T},pb::Pressure,p0::StagnationPressure,gas::PerfectGas) where T <: Area
  Mb = MachNumber(PressureRatio(pb/p0),Isentropic,gas=gas)
  Astar = AStar(A[end],Mb,gas=gas)
  _pressure_subsonic(A,Astar,p0,gas)
end

function _pressure_subsonic(A::AbstractVector{T},Astar::Area,p0::StagnationPressure,gas::PerfectGas) where {T <: Area}
  p_over_p0 = SubsonicPOverP0(A,Astar,Isentropic,gas=gas)
  p = Pressure[]
  for pr in p_over_p0
      push!(p,Pressure(pr*p0))
  end
  return p
end

function _pressure_subsonic_2(A::AbstractVector{T},Astar::Area,p0::StagnationPressure,gas::PerfectGas) where {T <: Area}
  p_over_p0 = SubsonicPOverP0(A,Astar,Isentropic,gas=gas)
  p = Pressure[]
  for pr in p_over_p0
      push!(p,Pressure(pr*p0))
  end
  return p
end

# Diverging nozzle processes

function _pressure_diverging_nozzle(A::AbstractVector{T},Astar::Area,p0::StagnationPressure,gas::PerfectGas) where {T <: Area}
    p = Pressure[]
    p_over_p0 = SupersonicPOverP0(A,Astar,Isentropic,gas=gas)
    for pr in p_over_p0
        push!(p,Pressure(pr*p0))
    end

    return p
end

function _machnumber_diverging_nozzle(A::AbstractVector{T},Astar::Area,gas::PerfectGas) where T <: Area
    return SupersonicMachNumber(A,Astar,Isentropic,gas=gas)
end

function _pressure_diverging_nozzle_with_shock(A::AbstractVector{T},As::Area,Astar1::Area,p01::StagnationPressure,gas::PerfectGas) where {T <: Area}
    i1 = _index_upstream_of_shock(A,As)

    p = Pressure[]
    pup_over_p01 = SupersonicPOverP0(A[1:i1],Astar1,Isentropic,gas=gas)
    for pr in pup_over_p01
        push!(p,Pressure(pr*p01))
    end

    M1, M2, Astar2, p02_over_p01, _ = _shock(As,Astar1,gas)

    p02 = StagnationPressure(p02_over_p01*p01)

    pdown_over_p02 = SubsonicPOverP0(A[i1+1:end],Astar2,Isentropic,gas=gas)
    for pr in pdown_over_p02
        push!(p,Pressure(pr*p02))
    end
    return p
end

function _machnumber_diverging_nozzle_with_shock(A::AbstractVector{T},As::Area,Astar1::Area,p01::StagnationPressure,gas::PerfectGas) where T <: Area
    i1 = _index_upstream_of_shock(A,As)

    M = SupersonicMachNumber(A[1:i1],Astar1,Isentropic,gas=gas)
    M1, M2, Astar2, _, _ = _shock(As,Astar1,gas)

    Mdown = SubsonicMachNumber(A[i1+1:end],Astar2,Isentropic,gas=gas)
    for Mi in Mdown
        push!(M,Mi)
    end
    return M
end

# Given an exit area Ae, shock area As, and Astar1, p01 in upstream nozzle
# return the pressure pe at exit
function _pressure_diverging_nozzle_with_shock(Ae::Area,As::Area,Astar1::Area,p01::StagnationPressure,gas::PerfectGas) where T <: Area
    value(Astar1) <= value(As) <= value(Ae) || error("Shock area outside of bounds")

    M1, M2, Astar2, p02_over_p01, _ = _shock(As,Astar1,gas)

    # conditions downstream of shock
    p02 = StagnationPressure(p02_over_p01*p01)

    pe_over_p02 = SubsonicPOverP0(Ae,Astar2,Isentropic,gas=gas)
    return Pressure(pe_over_p02*p02)
end

# Given an exit area Ae, shock area As, and Astar1, p01 in upstream nozzle
# return the Mach number Me at exit
function _machnumber_diverging_nozzle_with_shock(Ae::Area,As::Area,Astar1::Area,p01::StagnationPressure,gas::PerfectGas) where T <: Area
    value(Astar1) <= value(As) <= value(Ae) || error("Shock area outside of bounds")

    M1, M2, Astar2, _, _ = _shock(As,Astar1,gas)

    return SubsonicMachNumber(Ae,Astar2,Isentropic,gas=gas)
end

function _shock(As::Area,Astar1::Area,gas::PerfectGas)
    M1 = SupersonicMachNumber(AreaRatio(As/Astar1),Isentropic,gas=gas)
    M2 = MachNumber(M1,NormalShock,gas=gas)
    Astar2 = AStar(As,M2,Isentropic,gas=gas)
    p02_over_p01 = StagnationPressureRatio(M1,NormalShock,gas=gas)
    T2_over_T1 = TemperatureRatio(M1,NormalShock,gas=gas)
    return M1, M2, Astar2, p02_over_p01, T2_over_T1
end

# Find the index in array A corresponding to just before As.
# Assumes that A is diverging nozzle.
function _index_upstream_of_shock(A::AbstractVector{T},As::Area) where T <: Area
    value(A[1]) <= value(As) <= value(A[end]) || error("Shock area outside of bounds")
    indx = findall(x -> x >= 1.0,diff(sign.(value.(A) .- value(As))))
    return indx[end]
end

#=
Replace these in NozzleProcess...

function machnumber_nozzle(noz::Nozzle,As::Area,pb::Pressure,p01::StagnationPressure;gas::PerfectGas=DefaultPerfectGas)
  pb_subcrit = _pressure_subsonic_critical(noz,p01,gas=gas)
  pb_supcrit = _pressure_supersonic_shock_critical(noz,p01,gas=gas)

  M = pb > pb_subcrit ? machnumber_subsonic(areas(noz),pb,p01,gas=gas) :
     (pb > pb_supcrit ? machnumber_nozzle_with_shock(noz,As,pb,p01,gas=gas) :
                        machnumber_nozzle_supersonic(noz,gas=gas))

  return M
end

function pressure_nozzle(noz::Nozzle,As::Area,pb::Pressure,p01::StagnationPressure;gas::PerfectGas=DefaultPerfectGas)
  pb_subcrit = _pressure_subsonic_critical(noz,p01,gas=gas)
  pb_supcrit = _pressure_supersonic_shock_critical(noz,p01,gas=gas)
  p = pb > pb_subcrit ? pressure_subsonic(areas(noz),pb,p01,gas=gas) :
     (pb > pb_supcrit ? pressure_nozzle_with_shock(noz,As,pb,p01,gas=gas) :
                        pressure_nozzle_supersonic(noz,p01,gas=gas))

  # Set the final pressure to the back pressure
  p[end] = pb

  return p
end
=#

# Find area at which shock occurs
function _shock_area(noz::Nozzle,pb::Pressure,p01::StagnationPressure,::Type{SupersonicWithShockNozzle},gas::PerfectGas)
  pb_crit = _pressure_subsonic_critical(noz,p01,gas)
  pb < pb_crit || error("Back pressure too large for a shock")

  return Area(find_zero(x -> value(_pressure_diverging_nozzle_with_shock(nozzleexit(noz),Area(x),throat(noz),p01,gas))-value(pb),
    (value(diverging(noz)[1]),value(diverging(noz)[end-1])),order=8))
end

_shock_area(noz::Nozzle,pb,p01,::Type{SubsonicIsentropicNozzle},gas::PerfectGas) = throat(noz)
_shock_area(noz::Nozzle,pb,p01,::Type{SupersonicIsentropicNozzle},gas::PerfectGas) = throat(noz)

flow_quality(nozproc::NozzleProcess) = flow_quality(nozproc.noz,nozproc.pb,nozproc.p0,gas=nozproc.gas)

function flow_quality(noz::Nozzle,pb::Pressure,p0::StagnationPressure;gas::PerfectGas=DefaultPerfectGas)
  pb_subcrit = _pressure_subsonic_critical(noz,p0,gas)
  pb_supcrit_exitshock = _pressure_supersonic_shock_critical(noz,p0,gas)
  pb_supcrit = _pressure_supersonic_critical(noz,p0,gas)
  if pb > p0 || pb < Pressure(0.0)
    return error("Invalid back pressure")
  elseif pb > pb_subcrit
    return "Unchoked subsonic"
  elseif pb ≈ pb_subcrit
    return "Choked subsonic"
  elseif pb > pb_supcrit_exitshock
    return "Supersonic with normal shock"
  elseif pb > pb_supcrit
    return "Overexpanded supersonic"
  elseif pb ≈ pb_supcrit
    return "Perfectly expanded supersonic"
  else
    return "Underexpanded supersonic"
  end
  #println(pb_supcrit)
end

# Find critical back pressure for subsonic isentropic branch
function _pressure_subsonic_critical(noz::Nozzle,p0::StagnationPressure,gas::PerfectGas)
  pb_crit_over_p0 = SubsonicPOverP0(nozzleexit(noz),throat(noz),Isentropic,gas=gas)
  return Pressure(pb_crit_over_p0*p0)
end

# Find critical back pressure for supersonic isentropic branch with shock at exit
function _pressure_supersonic_shock_critical(noz::Nozzle,p0::StagnationPressure,gas::PerfectGas)
  pe_supcrit = _pressure_supersonic_critical(noz,p0,gas)
  Me_supcrit = SupersonicMachNumber(nozzleexit(noz),throat(noz),Isentropic,gas=gas)

  pb_over_pe_supcrit = PressureRatio(Me_supcrit,NormalShock,gas=gas)
  return Pressure(pb_over_pe_supcrit*pe_supcrit)
end

# Find critical back pressure for supersonic isentropic branch
function _pressure_supersonic_critical(noz::Nozzle,p0::StagnationPressure,gas::PerfectGas)
  pe_supcrit_over_p0 = SupersonicPOverP0(nozzleexit(noz),throat(noz),Isentropic,gas=gas)
  return Pressure(pe_supcrit_over_p0*p0)
end


text:
## Copyright (c) 2019 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD license.

# Set up pipes and start gnuplot process

# Async tasks to read/write to gnuplot's pipes.
const StartTasks = Condition()  # signal to start reading pipes

const ChanStdOut = Channel(10)
const ChanStdErr = Channel(10)

# This task reads all characters available from gnuplot's stdout.
@async begin
    wait(StartTasks)
    while true
        stdoutstr = String(readavailable(gstdout))
        process_running(gproc) || break
        put!(ChanStdOut, stdoutstr)
    end
end

# This task reads all characters available from gnuplot's stderr.
@async begin
    wait(StartTasks)
    while true
        stderrstr = String(readavailable(gstderr))
        process_running(gproc) || break
        put!(ChanStdErr, stderrstr)
    end
end

# initialize pipes and run gnuplot
gstdin = Pipe()
gstdout = Pipe()
gstderr = Pipe()
gproc = run(pipeline(`gnuplot`,
                     stdin = gstdin, stdout = gstdout, stderr = gstderr),
                     wait = false)
process_running(gproc) || error("There was a problem starting up gnuplot.")
close(gstdout.in)
close(gstderr.in)
close(gstdin.out)

# Start tasks to read and write gnuplot's pipes
NTasks = notify(StartTasks)


text:
# The license for this code is available at https://github.com/avik-pal/FastStyleTransfer.jl/blob/master/LICENSE.md

#-----------------------Utilities to load and save image---------------------------

im_mean = reshape([0.485, 0.458, 0.408], (1,1,3)) * 255
im_mean2 = reshape([0.485, 0.458, 0.408], (1,1,3,1)) * 255 |> gpu

# NOTE: The image returned is scaled to equal dimensions on both side
function load_image(filename; size_img::Int = -1, scale::Int = -1)
    img = load(filename)
    global original_size = size(img)
    if size_img != -1
        img = imresize(img, (size_img, size_img))
    elseif scale != -1
        dims = size(img, 1)
        img = imresize(img, (dims, dims))
    end
    img = Float32.(channelview(img)) * 255
    ndims(img) == 2 && return img
    permutedims(img, [3,2,1]) .- im_mean
end

function save_image(filename, img, display_img::Bool = false)
    img = reshape(img, (size(img, 2), size(img, 1), 3))
    img = (img .+ im_mean) / 255
    img = permutedims(img, [3,2,1])
    img -= minimum(img)
    img /= maximum(img)
    img = colorview(RGB{eltype(img)}, img)
    img = imresize(img, original_size)
    display_img && display(img)
    save(filename, img)
end

function load_dataset(path, batch, total)
    z = readdir(path)
    indices = randperm(length(z))[1:total]
    paths = [joinpath(path, i) for i in z[indices]]
    images = []
    for (counts, i) in enumerate(paths)
        img = load_image(i, size_img = 224)
        ndims(img) == 3 ? push!(images, img) : total -= 1 # Hack to avoid errors in case of MSCOCO
        counts % 100 == 0 && @info "$counts images have been loaded"
    end
    [cat(images[i]..., dims=4) for i in Iterators.partition(1:total, batch)]
end

# NOTE: The wrapper is quite slow
# function gram_matrix(x::CuArray)
#     w, h, ch, b = size(x)
#     local features = reshape(x, w*h, ch, b)
#     features = [features[:,:,i] for i in 1:b]
#     res = CuArrays.BLAS.gemm_batched('T', 'N', features, features)
#     cat(3, res...) / (w * h * ch)
# end

function gram_matrix(x)
    w, h, ch, b = size(x)
    local features = reshape(x, w*h, ch, b)
    features = [features[:,:,i] for i in 1:b]
    cat([features[i]' * features[i] for i in 1:b]..., dims=4) / Float32(w * h * ch)
end

##
##
## Deprecated?
##
##
#----------------------------Extension of certain functions------------------------------
#using Base.std

# Not as per the defination. Just as hack to get it working
#Base.std(x::TrackedArray, dim::Array; mean = Base.mean(x, dim)) =
#    sqrt.(sum((x .- mean).^2, dim) ./ (prod(size(x)[dim])-1))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SuiteSparse32_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SuiteSparse32")
JLLWrappers.@generate_main_file("SuiteSparse32", UUID("ca45d3f4-326b-53b0-9957-23b75aacb3f2"))
end  # module SuiteSparse32_jll


text:
module HDF5Utils

using Mmap, Random, Distributed
using Requires, ProgressMeter, Glob
using HDF5, FillArrays, DiskArrays
using FillArrays: Zeros, Fill

using HDF5: HDF5Dataset, HDF5Attribute, DatasetOrAttribute
using HDF5: H5T_COMPOUND, H5T_STRING, H5T_CSET_UTF8
using HDF5: H5P_FILE_ACCESS, H5P_FILE_CREATE, H5P_DEFAULT, H5P_DATASET_CREATE
using HDF5: H5S_ALL, H5S_SELECT_SET, H5S_SCALAR, H5F_CLOSE_STRONG, H5F_CLOSE_DEFAULT
using HDF5: h5t_create, h5t_insert, h5t_close, h5t_copy, h5t_get_size, h5t_set_size
using HDF5: h5t_get_cset, h5t_set_cset, h5t_get_class, h5t_get_nmembers
using HDF5: h5t_get_member_type, h5t_get_member_name, h5t_is_variable_str
using HDF5: h5s_create, h5s_close, h5s_select_hyperslab, h5d_read, h5d_write, h5f_close
using HDF5: hdf5_type_id, hdf5_to_julia, hdf5_to_julia_eltype, hyperslab
using HDF5: Herr, Hid, Hsize, isnull, libhdf5, writearray, readarray, checkvalid, checkprops
import DiskArrays: eachchunk, haschunks, readblock!, writeblock!, GridChunks, Chunked, Unchunked

export HDF5DiskArray
export d_zeros, copy_batch!, write_batch
export read_nonarray, write_nonarray
export h5load, h5loadv, h5save, h5readmmap, tryreadmmap
export h5concat, h5concat!
export MaxLenString, MLString
export d_create_virtual, VirtualLayout, VirtualSource
export enable_dag, disable_dag

include("diskarrays.jl")
include("util.jl")
include("batch.jl")
include("concat.jl")
include("io.jl")
include("alignment.jl")
include("conversion.jl")
include("mlstring.jl")
include("namedtuple.jl")
include("hyperslab.jl")
include("virtual.jl")
include("auto.jl")

@init @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" include("npystring.jl")

end # module


text:
using Memoize

@memoize fib(n) = n < 3 ? 1 : fib(n-1) + fib(n-2);
@time fib(40)
@time fib(40)
@time fib(39)

#=
julia> using Memoize

julia> @memoize fib(n) = n < 3 ? 1 : fib(n-1) + fib(n-2)
fib (generic function with 1 method)

julia> @time fib(40)
  0.027960 seconds (64.56 k allocations: 3.378 MiB)
102334155

julia> @time fib(40)
  0.000003 seconds (7 allocations: 208 bytes)
102334155

julia> @time fib(39)
  0.000003 seconds (7 allocations: 208 bytes)
63245986
=#

# Caching.jl

using Caching, CSV, DataFrames

@cache function read_csv(filename::AbstractString)
    println("Reading file: ", filename)
    @time df = CSV.File(filename) |> DataFrame
    return df
end

df = read_csv("Film_Permits.csv");

#=
julia> using Caching, CSV, DataFrames

julia> @cache function read_csv(filename::AbstractString)
           println("Reading file: ", filename)
           @time df = CSV.File(filename) |> DataFrame
           return df
       end

read_csv (cache with 0 entries, 0 in memory 0 on disk)

julia> df = read_csv("Film_Permits.csv")
Reading file: Film_Permits.csv
  7.853251 seconds (25.52 M allocations: 1010.414 MiB, 10.49% gc time)

julia> size(df)
(42576, 14)

julia> df_again = read_csv("Film_Permits.csv");

julia> df === df_again
true

julia> propertynames(read_csv)
(:name, :filename, :func, :cache, :offsets, :history, :max_size)

julia> read_csv.cache
Dict{UInt64,Any} with 1 entry:
  0x602537e6ffb4d3f0 => 42576×14 DataFrame. Omitted printing of 11 columns…

julia> read_csv.filename
"/Users/tomkwong/Downloads/_7d1d434983e10ae8_.bin"

shell> ls -la /Users/tomkwong/Downloads/_7d1d434983e10ae8_.bin
ls: /Users/tomkwong/Downloads/_7d1d434983e10ae8_.bin: No such file or directory

julia> @persist! read_csv
read_csv (cache with 1 entry, 1 in memory 1 on disk)

shell> ls -1s /Users/tomkwong/Downloads/_7d1d434983e10ae8_.bin
33760 /Users/tomkwong/Downloads/_7d1d434983e10ae8_.bin

julia> read_csv
read_csv (cache with 1 entry, 1 in memory 1 on disk)

julia> @empty! read_csv
read_csv (cache with 1 entry, 0 in memory 1 on disk)

julia> df = read_csv("Film_Permits.csv");

julia> size(df)
(42576, 14)

julia> read_csv
read_csv (cache with 1 entry, 0 in memory 1 on disk)

julia> @syncache! read_csv "disk"
read_csv (cache with 1 entry, 1 in memory 1 on disk)

=#


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for logging linear system equation solvers, and
## providing stopping rules.
##
## Contents
## - Abstract Types
## - `log_update!` function documentation
## - Loggers
## - Export Statements
##
############################################################################################

# Dependencies: LinearAlgebra

#############################################
# Abstract Types
#############################################

"""
    LinSysSolverLog

Abstract supertype for logging information about a randomized linear system solver.

Note, all descendants that are `mutable struct`s must have fields
- `iterations::Int64`, which tracks the number of iterations of the solver
- `converged::Bool`, which tracks whether the solver has converged to the solution (by some
    notion of convergence that is specified by the user).
"""
abstract type LinSysSolverLog end

#############################################
# `log_update!` function docstring
#############################################
"""
    log_update!(
        log::LinSysSolverLog,
        sampler::LinSysSampler,
        x::AbstractVector,
        samp:s:Tuple,
        iter::Int64,
        A,
        b)

A common interface for specifying different strategies for updating the `log` with supertype
    `LinSysSolverLog`. A `sampler` of supertype `LinSysSampler` can be used to provide
    specific multiple implementations for the same `log` type. `x` is the iterate value.
    `samp` is the output generated by the `sample` function (see `linears_samplers.jl`).
    `iter` is the current iteration counter. `A` and `b` specify the linear system.
"""
function log_update!(
    log::Nothing,
    sampler::LinSysSampler,
    x::AbstractVector,
    samp::Tuple,
    iter::Int64,
    A,
    b
)
    return nothing
end

#############################################
# Loggers
#############################################
include("linear_solver_logs/solve_log_oracle.jl")
include("linear_solver_logs/solve_log_full.jl")

#############################################
# Exports
#############################################
#export LinSysSolverLog
#export LSLogOracle, LSLogFull


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

module RangesTest

import Dates
using Test
using ExtendedDates

using InteractiveUtils: subtypes

let
    YD = YearDate
    SD = SemesterDate
    QD = QuarterDate
    MD = MonthDate
    WD = WeekDate
    DD = DayDate
    UD = UndatedDate
    for (f1, l1, f2, l2, f3, l3, f4, l4, period) in (
        (
            YD(2014),
            YD(2013),
            YD(2014),
            YD(2014),
            YD(-20),
            YD(20),
            typemin(YD),
            typemax(YD),
            (Dates.Year,),
        ),
        (
            SD(2014, 1),
            SD(2013, 1),
            SD(2014, 1),
            SD(2014, 1),
            SD(-2000, 1),
            SD(2000, 1),
            typemin(SD),
            typemax(SD),
            (Dates.Year, Semester),
        ),
        (
            QD(2014, 1),
            QD(2013, 1),
            QD(2014, 1),
            QD(2014, 1),
            QD(-2000, 1),
            QD(2000, 1),
            typemin(QD),
            typemax(QD),
            (Dates.Year, Semester, Dates.Quarter),
        ),
        (
            MD(2014, 1),
            MD(2013, 1),
            MD(2014, 1),
            MD(2014, 1),
            MD(-2000, 1),
            MD(2000, 1),
            typemin(MD),
            typemax(MD),
            (Dates.Year, Semester, Dates.Quarter, Dates.Month),
        ),
        (
            WD(2014, 1),
            WD(2013, 1),
            WD(2014, 1),
            WD(2014, 1),
            WD(1, 1),
            WD(2000, 1),
            typemin(WD),
            typemax(WD),
            (Dates.Week,),
        ),
        (
            DD(2014, 1, 1),
            DD(2013, 1, 1),
            DD(2014, 1, 1),
            DD(2014, 1, 1),
            DD(-2000, 1, 1),
            DD(2000, 1, 1),
            typemin(DD),
            typemax(DD),
            (Dates.Year, Semester, Dates.Quarter, Dates.Month, Dates.Week, Dates.Day),
        ),
        (
            UD(2014),
            UD(2013),
            UD(2014),
            UD(2014),
            UD(-2000),
            UD(2000),
            typemin(UD),
            typemax(UD),
            (Undated,),
        ),
    )

        for P in period
            for pos_step in (P(1), P(2), P(50), P(2048), P(10000))
                # empty range
                dr = f1:pos_step:l1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == f1
                @test last(dr) < f1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == typeof(f1)[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) < f1
                @test last(reverse(dr)) >= f1
                @test issorted(dr)
                @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                @test !(f1 in dr)
                @test !(l1 in dr)
                @test !(f1 - pos_step in dr)
                @test !(l1 + pos_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3))# , (f4, l4))
                    dr = f:pos_step:l
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == f
                    @test last(dr) <= l
                    @test minimum(dr) == first(dr)
                    @test maximum(dr) == last(dr)
                    @test dr[1] == f
                    @test dr[end] <= l
                    @test iterate(dr) == (first(dr), (length(dr), 1))
                    if len < 10#000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test last(reverse(dr)) == f
                    @test issorted(dr)
                    @test f in dr
                end
            end

            for neg_step in (P(-1), P(-2), P(-50), P(-2048), P(-10000))
                # empty range
                dr = l1:neg_step:f1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == l1
                @test last(dr) > l1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == typeof(f1)[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) > l1
                @test last(reverse(dr)) <= l1
                @test issorted(dr)
                @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                @test !(l1 in dr)
                @test !(l1 in dr)
                @test !(l1 - neg_step in dr)
                @test !(l1 + neg_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3))#, (f4, l4))
                    dr = l:neg_step:f
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == l
                    @test last(dr) >= f
                    @test minimum(dr) == last(dr)
                    @test maximum(dr) == first(dr)
                    @test dr[1] == l
                    @test dr[end] >= f
                    @test iterate(dr) == (first(dr), (length(dr), 1))

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test issorted(dr) == (len <= 1)
                    @test l in dr
                end
            end
        end
    end
end

# Dates are physical units, and ranges should require an explicit step.
# See #19896 and https://discourse.julialang.org/t/type-restriction-on-unitrange/6557/12
if VERSION >= v"1.7.0"
    @test_throws MethodError DayDate(2013, 1, 1):DayDate(2013, 2, 1)
end

# All the range representations we want to test
# Date ranges
dr = DayDate(2013, 1, 1):Day(1):DayDate(2013, 2, 1)
dr1 = DayDate(2013, 1, 1):Day(1):DayDate(2013, 1, 1)
dr2 = DayDate(2013, 1, 1):Day(1):DayDate(2012, 2, 1) # empty range
dr3 = DayDate(2013, 1, 1):Day(-1):DayDate(2012) # negative step
# Big ranges
dr4 = DayDate(0):Day(1):DayDate(20000, 1, 1)
dr5 = DayDate(0):Day(1):DayDate(200000, 1, 1)
dr6 = DayDate(0):Day(1):DayDate(2000000, 1, 1)
dr7 = DayDate(0):Day(1):DayDate(20000000, 1, 1)
dr8 = DayDate(0):Day(1):DayDate(200000000, 1, 1)
dr9 = typemin(DayDate):Day(1):typemax(DayDate)
# Other steps
dr10 = typemax(DayDate):Day(-1):typemin(DayDate)
dr11 = typemin(DayDate):Week(1):typemax(DayDate)

dr12 = typemin(DayDate):Month(1):typemax(DayDate)
dr13 = typemin(DayDate):Year(1):typemax(DayDate)

dr14 = typemin(DayDate):Week(10):typemax(DayDate)
dr15 = typemin(DayDate):Month(100):typemax(DayDate)
dr16 = typemin(DayDate):Year(1000):typemax(DayDate)
dr17 = typemax(DayDate):Week(-10000):typemin(DayDate)
dr18 = typemax(DayDate):Month(-100000):typemin(DayDate)
dr19 = typemax(DayDate):Year(-1000000):typemin(DayDate)
dr20 = typemin(DayDate):Day(2):typemax(DayDate)

drs = Any[
    dr,
    dr1,
    dr2,
    dr3,
    dr4,
    dr5,
    dr6,
    dr7,
    dr8,
    dr9,
    dr10,
    dr11,
    dr12,
    dr13,
    dr14,
    dr15,
    dr16,
    dr17,
    dr18,
    dr19,
    dr20,
]
drs2 = map(x -> DayDate(first(x)):step(x):DayDate(last(x)), drs)

@test map(length, drs) == map(x -> size(x)[1], drs)
@test map(length, drs) == map(x -> length(DayDate(first(x)):step(x):DayDate(last(x))), drs)
@test map(length, drs) == map(x -> length(reverse(x)), drs)
@test all(x -> findall(in(x), x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x -> reverse(x) == range(last(x), step = -step(x), length = length(x)), drs)
@test all(x -> minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x -> maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test all(drs[1:3]) do dd
    for (i, d) in enumerate(dd)
        @test d == (first(dd) + Day(i - 1))
    end
    true
end
@test_throws MethodError dr .+ 1
a = DayDate(2013, 1, 1)
b = DayDate(2013, 2, 1)
@test map!(x -> x + Day(1), Vector{DayDate}(undef, 32), dr) ==
      [(a+Day(1)):Day(1):(b+Day(1));]
@test map(x -> x + Day(1), dr) == [(a+Day(1)):Day(1):(b+Day(1));]

@test map(x -> a in x, drs[1:4]) == [true, true, false, true]
@test a in dr
@test b in dr
@test DayDate(2013, 1, 3) in dr
@test DayDate(2013, 1, 15) in dr
@test DayDate(2013, 1, 26) in dr
@test !(DayDate(2012, 1, 1) in dr)

@test all(x -> sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x -> step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test length(b:Day(-1):a) == 32
@test length(b:Day(1):a) == 0
@test length(b:Day(1):a) == 0
@test length(a:Day(2):b) == 16
@test last(a:Day(2):b) == DayDate(2013, 1, 31)
@test length(a:Day(7):b) == 5
@test last(a:Day(7):b) == DayDate(2013, 1, 29)
@test length(a:Day(32):b) == 1
@test last(a:Day(32):b) == DayDate(2013, 1, 1)
@test (a:Day(1):b)[1] == DayDate(2013, 1, 1)
@test (a:Day(1):b)[2] == DayDate(2013, 1, 2)
@test (a:Day(1):b)[7] == DayDate(2013, 1, 7)
@test (a:Day(1):b)[end] == b
@test first(a:Day(1):DayDate(20000, 1, 1)) == a
@test first(a:Day(1):DayDate(200000, 1, 1)) == a
@test first(a:Day(1):DayDate(2000000, 1, 1)) == a
@test first(a:Day(1):DayDate(20000000, 1, 1)) == a
@test first(a:Day(1):DayDate(200000000, 1, 1)) == a
@test first(a:Day(1):typemax(DayDate)) == a
@test first(typemin(DayDate):Day(1):typemax(DayDate)) == typemin(DayDate)

@test length(typemin(DayDate):Week(1):typemax(DayDate)) == 26351950414948059
# Big Month/Year ranges
@test length(typemin(DayDate):Month(1):typemax(DayDate)) == 6060531933867600
@test length(typemin(DayDate):Year(1):typemax(DayDate)) == 505044327822300

c = DayDate(2013, 6, 1)
@test length(a:Month(1):c) == 6
@test [a:Month(1):c;] == [a + Month(1) * i for i = 0:5]
@test [a:Month(2):DayDate(2013, 1, 2);] == [a]
@test [c:Month(-1):a;] == reverse([a:Month(1):c;])

@test length(range(Date(2000), step = Day(1), length = 366)) == 366
let n = 100000
    local a = DayDate(2000)
    for i = 1:n
        @test length(range(a, step = Day(1), length = i)) == i
    end
    return a + Day(n)
end

@test typeof(step(DayDate(2000):Day(1):DayDate(2001))) == Day

a = DayDate(2013, 1, 1)
b = DayDate(2013, 2, 1)
d = DayDate(2020, 1, 1)
@test length(a:Year(1):d) == 8
@test first(a:Year(1):d) == a
@test last(a:Year(1):d) == d
@test length(a:Month(12):d) == 8
@test first(a:Month(12):d) == a
@test last(a:Month(12):d) == d
@test length(a:Week(52):d) == 8
@test first(a:Week(52):d) == a
@test last(a:Week(52):d) == DayDate(2019, 12, 24)
@test length(a:Day(365):d) == 8
@test first(a:Day(365):d) == a
@test last(a:Day(365):d) == DayDate(2019, 12, 31)

a = DayDate(2013, 1, 1)
b = DayDate(2013, 2, 1)
@test length(a:Year(1):DayDate(2020, 2, 1)) == 8
@test length(a:Year(1):DayDate(2020, 6, 1)) == 8
@test length(a:Year(1):DayDate(2020, 11, 1)) == 8
@test length(a:Year(1):DayDate(2020, 12, 31)) == 8
@test length(a:Year(1):DayDate(2021, 1, 1)) == 9
@test length(DayDate(2000):Year(-10):DayDate(1900)) == 11
@test length(DayDate(2000, 6, 23):Year(-10):DayDate(1900, 2, 28)) == 11
@test length(DayDate(2000, 1, 1):Year(1):DayDate(2000, 2, 1)) == 1

let n = 100000
    local a, b
    a = b = DayDate(0)
    for i = 1:n
        @test length(a:Year(1):b) == i
        b += Year(1)
    end
end

let n = 10000,
    a = DayDate(1985, 12, 5),
    b = DayDate(1986, 12, 27),
    c = DayDate(1985, 12, 5),
    d = DayDate(1986, 12, 27)

    for i = 1:n
        @test length(a:Month(1):b) == 13
        @test length(a:Year(1):b) == 2
        @test length(c:Month(1):d) == 13
        @test length(c:Year(1):d) == 2
        a += Day(1)
        b += Day(1)
    end
end

let n = 100000
    local a, b
    a = b = DayDate(2000)
    for i = 1:n
        @test length(a:Month(1):b) == i
        b += Month(1)
    end
end

@test length(Year(1):Year(1):Year(10)) == 10
@test length(Year(10):Year(-1):Year(1)) == 10
@test length(Year(10):Year(-2):Year(1)) == 5
if VERSION >= v"1.7.0"
    @test length(typemin(Year):Year(1):typemax(Year)) == 0 # overflow
    @test_throws MethodError DayDate(0):DayDate(2000)
end
@test_throws MethodError DayDate(0):Year(10)
@test length(range(DayDate(2000), step = Day(1), length = 366)) == 366
@test last(range(DayDate(2000), step = Day(1), length = 366)) == DayDate(2000, 12, 31)
@test last(range(DayDate(2001), step = Day(1), length = 365)) == DayDate(2001, 12, 31)
@test last(range(DayDate(2000), step = Day(1), length = 367)) ==
      last(range(DayDate(2000), step = Month(12), length = 2)) ==
      last(range(DayDate(2000), step = Year(1), length = 2))

# Issue 5
lastdaysofmonth = [DayDate(2014, i, daysinmonth(2014, i)) for i = 1:12]
@test [DayDate(2014, 1, 31):Month(1):DayDate(2015);] == lastdaysofmonth

# Range addition/subtraction:
let d = Day(1)
    @test (DayDate(2000):d:DayDate(2001)) + d == (DayDate(2000)+d:d:DayDate(2001)+d)
    @test (DayDate(2000):d:DayDate(2001)) - d == (DayDate(2000)-d:d:DayDate(2001)-d)
end

end  # RangesTest module


text:
using Documenter, Plots
using HydrogenicAtoms

makedocs(
    sitename="HydrogenicAtoms",
    modules = [HydrogenicAtoms],
    pages = [
        "Home" => "index.md",
        "dirac.md"
    ]
)


text:
using LinearTimeVaryingModelsBase
using Test

# write your own tests here
t = Trajectory(randn(2,4), randn(2,4))
@test length(t) == 3
@test t.nx == t.nu == 2

for (x,u) in t
    println(x,u)
end


text:
mutable struct Core <: CoreInterfaces.ICore
    config::CoreInterfaces.IConfig
    logs::CoreInterfaces.ILog

    function Core()
        o = new()
        o.config = Config.SConfig()

        o.logs = Logs.SLog(o.config)
        Logs.configure(o.logs)
        o
    end
end


text:
using DynamicHMCModels, MCMCChains, Random

#Random.seed!(1233)

ProjDir = @__DIR__
cd(ProjDir)

include(joinpath(@__DIR__, "../Benchmarks/LBA/LBA_functions.jl"))

Base.@kwdef struct LBAModel{T}
   data::T
   N::Int
   Nc::Int
end

function make_transformation(model::LBAModel)
    as((v=as(Array,asℝ₊,Nc), A=asℝ₊, k=asℝ₊, tau=asℝ₊))
end

for i = 1:3
  N = 10
  v = [1.0, 1.5]
  Nc = length(v)
  data=simulateLBA(;Nd=N,v=v,A=.8,k=.2,tau=.4)

  model = LBAModel(; data=data, N=N, Nc=Nc)

  function (model::LBAModel)(θ)
      @unpack data=model
      @unpack v,A,k,tau=θ
      d=LBA(ν=v,A=A,k=k,τ=tau)
      minRT = minimum(x->x[2],data)
      logpdf(d,data)+sum(logpdf.(TruncatedNormal(0,3,0,Inf),v)) +
      logpdf(TruncatedNormal(.8,.4,0,Inf),A)+logpdf(TruncatedNormal(.2,.3,0,Inf),k)+
      logpdf(TruncatedNormal(.4,.1,0,minRT),tau)
  end

  d = [(c,r) for (c,r) in zip(data.choice,data.rt)]

  p = LBAModel(d,N,Nc)
  A = rand(Normal(.8, 1.0), 1)[1]
  k = rand(Normal(.2, 1.0), 1)[1]
  tau = rand(Normal(.4, 1.0), 1)[1]
  v=rand(Normal(0.0, 1.0),Nc)
  p((v=v,A=A,k=k,tau=tau))
  display(d)

end


text:
using Test
using Symbolics
using ForwardDiff
using LinearAlgebra
using SparseArrays
using DirectTrajectoryOptimization
const DTO = DirectTrajectoryOptimization

include("objective.jl")
include("dynamics.jl")
include("constraints.jl")
include("hessian_lagrangian.jl")
include("solve.jl")


text:
import Pkg; Pkg.update("ComplexRegions")
using Documenter, ComplexRegions, Plots

ENV["GKSwstype"] = "100"
DocMeta.setdocmeta!(ComplexRegions, :DocTestSetup, :(using ComplexRegions); recursive=true)
makedocs(sitename="ComplexRegions.jl",
    format = Documenter.HTML(),
    authors = "Toby Driscoll",
    pages = [
        "Introduction" => "index.md",
        "Curves" => "curves.md",
        "Paths" => "paths.md",
        "Polygons" => "polygons.md",
        "Intersections" => "intersections.md",
        "Regions" => "regions.md",
        "Möbius" => "mobius.md",
        "API Reference" => "api.md"
		],
	modules = [ComplexRegions],
    doctest = true
    )

deploydocs(
    repo = "github.com/complexvariables/ComplexRegions.jl.git"
#    versions = ["v#.#"],
#    make = nothing
    )


text:
using Weber
using Base.Test
include("find_timing.jl")

const check_timing = get(ENV,"WEBER_TIMING_TESTS","Yes") != "No"
const moment_eps = 2.5e-3

# warm up JIT...
find_timing() do
  addtrial(moment(50ms,() -> record(:a)),
           moment(50ms,() -> record(:b)) >> moment(100ms,() -> record(:d)),
           moment(100ms,() -> record(:c)),
           moment(100ms,() -> record(:e)))
end

comp_events,comp_times,_ = find_timing() do
  addtrial(moment(50ms,() -> record(:a)),
           moment(50ms,() -> record(:b)) >> moment(100ms,() -> record(:d)),
           moment(100ms,() -> record(:c)),
           moment(100ms,() -> record(:e)))
end

@testset "Compound Moments" begin
  @test comp_events == [:a,:b,:c,:d,:e]
  if check_timing
    comp_diff = maximum(diff(comp_times) - 0.05)
    @test comp_diff < moment_eps
  end
end


text:
struct SpecParams{T<:AF}
    lines::AA{T,1}
    depths::AA{T,1}
    conv_blueshifts::AA{T,1}
    variability::AA{Bool,1}
    coverage::Tuple{T,T}
    resolution::T
    lambdas::AA{T,1}
    soldata::SolarData{T}
end

"""
    SpecParams(; lines=[], depths=[], variability=[], resolution=7e5)

Construct a `SpecParams` composite type instance. If `variability` is not specified, all lines are variable by default.

# Arguments
- `lines::AbstractArray{Float64,1}=[]`: List of line centers (in angstroms)
- `depths::AbstractArray{Float64,1}=[]`: List of line depths
- `variability::AbstractArray{Bool,1}=[]`: Array of booleans controlling whether corresponding line has variability.
- `resolution::Float64=7e8`: Spectral resolution of spectrum
"""
function SpecParams(;lines=[], depths=[], variability=[], resolution=7e5,
                    buffer=0.75, fixed_width=false, fixed_bisector=false,
                    extrapolate=true, contiguous_only=false)
    @assert !isempty(lines)
    @assert !isempty(depths)
    @assert length(lines) == length(depths)
    @assert buffer >= 0.75

    # assign depths if needed, otherwise check same number of centers + depths
    if any(isnan.(depths))
        depths = rand(length(lines))
    else
        @assert !any(depths .> 1.0)
        @assert !any(depths .< 0.0)
    end

    # read in convective blueshifts
    df = CSV.read(soldir * "../convective_blueshift.dat", DataFrame,
                  header=1, delim=",", type=Float64)

    # assign convective blueshifts
    blueshifts = similar(lines)
    for i in eachindex(depths)
        ind = searchsortednearest(df.depth, depths[i])
        blueshifts[i] = df.blueshift[ind]
    end

    # assign fixed_width booleans
    if isempty(variability)
        variability = trues(length(lines))
    else
        @assert length(lines) == length(variability)
    end

    # calculate coverage
    # TODO: constant delta v, not wavelength?
    coverage = (minimum(lines) - buffer, maximum(lines) + buffer)

    # get generator for lambdas
    dλ = coverage[1]/resolution
    lambdas = coverage[1]:dλ:coverage[2]

    # get observational data
    soldata = SolarData(relative=true, fixed_width=fixed_width,
                        fixed_bisector=fixed_bisector, extrapolate=extrapolate,
                        contiguous_only=contiguous_only)
    return SpecParams(lines, depths, blueshifts, variability, coverage, resolution, lambdas, soldata)
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
### Aetherlang type system

## Type-specific tools
include("AetherlangExceptions.jl")
include("AetherlangPhantomCollection.jl")
include("AetherlangArrays.jl")
include("AetherlangStrings.jl")

## A named tuple for type references
const AETH_BUILTIN_TYPES = (
    unknown="unknown", # this is a type for system types that are not builtin
    number="number",
    string="string",
    symbol="symbol",
    timepoint="timepoint",
    void="void",
    callable="function",
    collection="array", # is only used to refer to a linear collection, such as tuple or array. define other types for other collections
    phantom="phantomarray" # is only used for phantoms on linear collections of type AETH_BUILTIN_TYPES.collection. define others types for phantoms on your structures
)

## AetherlangVoid type
struct AetherlangVoid
end
function Base.show(io::IO, o::AetherlangVoid)
    print(io, "void")
end

##
struct AetherlangObject{T}
    """A container type for Aetherlang. Can contain any Aetherlang type"""
    dataref::T
    type::Ref{String}
end
function Base.show(io::IO, o::AetherlangObject)
    try
        aetherlang_show(io, o.dataref)
    catch
        print(io, o.dataref)
    end
end

AetherlangObject() = AetherlangObject{AetherlangVoid}(AetherlangVoid(), Ref(AETH_BUILTIN_TYPES.void))
AetherlangObject(n::Nothing) = AetherlangObject{AetherlangVoid}(AetherlangVoid(), Ref(AETH_BUILTIN_TYPES.void))
AetherlangObject(s::String, issymbol::Bool) = AetherlangObject{String}(s, issymbol ? Ref(AETH_BUILTIN_TYPES.symbol) : Ref(AETH_BUILTIN_TYPES.string))
AetherlangObject(s::String) = AetherlangObject(s, false)
AetherlangObject(c::T where T <: AbstractVector) = AetherlangObject{typeof(c)}(c, Ref(AETH_BUILTIN_TYPES.collection))
#AetherlangObject(c::Tuple) = AetherlangObject{typeof(c)}(c, Ref(AETH_BUILTIN_TYPES.collection))
AetherlangObject(p::AetherlangPhantomCollection) = AetherlangObject{AetherlangPhantomCollection}(p, Ref(AETH_BUILTIN_TYPES.phantom))
AetherlangObject(a::T where T<:Any) = AetherlangObject{typeof(a)}(a, Ref(AETH_BUILTIN_TYPES.unknown))
function AetherlangObject(n::T where T<:Number)
    if typeof(n) <: Int
        if n <= 127 && n >= -128
            return AetherlangObject{Int8}(Int8(n), Ref(AETH_BUILTIN_TYPES.number))
        end
    end
    AetherlangObject{typeof(n)}(n, Ref(AETH_BUILTIN_TYPES.number))
end
AetherlangObject(b::Bool) = b ? AetherlangObject{Int8}(Int8(1), Ref(AETH_BUILTIN_TYPES.number)) : AetherlangObject()
AetherlangObject(c::Char) = AetherlangObject{String}(string(c), Ref(AETH_BUILTIN_TYPES.string))

## Namespace definition
Namespace = Dict{String, AetherlangObject{T} where T}

## Timepoint type and the corresponding constructor for AetherlangObject
struct AetherlangTimepoint
    names::Set{String} # a set of names that exist at this timepoint in given dimension
    line::UInt
end
AetherlangObject(t::AetherlangTimepoint) = AetherlangObject{AetherlangTimepoint}(t, Ref(AETH_BUILTIN_TYPES.timepoint))
function aetherlang_show(io::IO, t::AetherlangTimepoint)
    print(io, "(timepoint at line $(t.line))")
end

## More includes
include("AetherlangFunctions.jl")

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text: function f(x) x.a end  # aucun champ a

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# Orchestra types.
module Types

export Transformer,
       Learner,
       TestLearner,
       fit!,
       transform!

# All transformer types must have implementations
# of function `fit!` and `transform!`.
abstract Transformer

# Learner abstract type which all machine learners implement.
abstract Learner <: Transformer

# Test learner.
# Used to separate production learners from test.
abstract TestLearner <: Learner

# Trains transformer on provided instances and labels.
#
# @param transformer Target transformer.
# @param instances Training instances.
# @param labels Training labels.
function fit!(transformer::Transformer, instances::Matrix, labels::Vector)
  error(typeof(transformer), " does not implement fit!")
end

# Trains transformer on provided instances and labels.
#
# @param transformer Target transformer.
# @param instances Original instances.
# @return Transformed instances.
function transform!(transformer::Transformer, instances::Matrix)
  error(typeof(transformer), " does not implement transform!")
end

end # module