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

using Markdown
using InteractiveUtils

# ╔═╡ 49674784-f070-488d-b53b-9cebf0ce101e
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Agents", version="4"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="DataFrames", version="1"),
        Pkg.PackageSpec(name="LightGraphs", version="1"),
        Pkg.PackageSpec(name="Distributions", version="0.25"),
		Pkg.PackageSpec(name="GraphRecipes", version="0.5"),
    ])
    using Agents
	using PlutoUI
	using Plots
	using Random
	using DataFrames
	using LightGraphs
	using Distributions
	using GraphRecipes
	using LinearAlgebra: diagind
	PlutoUI.TableOfContents()
end

# ╔═╡ e5f4abe0-be15-11eb-038f-5da4c19a690d
md"""

# ABM Example 2: The spread of SARS-CoV-2 (Graph model)

[Source](https://juliadynamics.github.io/Agents.jl/stable/examples/sir/)

Here we add one more category of individuals: those who are infected, but do not know it. Transmission rate for infected and diagnosed individuals is lower than infected and undetected.
"""

# ╔═╡ 07c4c6aa-c7bb-4543-868d-975de9f37ba9
mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::Int            # Which city
    days_infected::Int  # number of days since is infected
    status::Symbol      # 1: S, 2: I, 3:R
end

# ╔═╡ 99c0bed6-a533-4058-bf85-c81fac0f5976
function make_SIRgraph(;
    Ns,
    migration_rates,
    β_und,
    β_det,
    infection_period = 30,
    reinfection_probability = 0.05,
    detection_time = 14,
    death_rate = 0.02,
    Is = [zeros(Int, length(Ns) - 1)..., 1],
    seed = 0,
)

    rng = MersenneTwister(seed)
    @assert length(Ns) ==
    length(Is) ==
    length(β_und) ==
    length(β_det) ==
    size(migration_rates, 1) "length of Ns, Is, and B, and number of rows/columns in migration_rates should be the same "
    @assert size(migration_rates, 1) == size(migration_rates, 2) "migration_rates rates should be a square matrix"

    C = length(Ns) # Number of cities

    # normalize migration_rates
    migration_rates_sum = sum(migration_rates, dims = 2)
    for c in 1:C
        migration_rates[c, :] ./= migration_rates_sum[c]
    end

    properties = (;
        Ns,
        Is,
        β_und,
        β_det,
        migration_rates,
        infection_period,
        reinfection_probability,
        detection_time,
        C,
        death_rate
    )

    space = GraphSpace(complete_digraph(C))
    model = ABM(PoorSoul, space; properties, rng)

    # Add initial individuals
    for city in 1:C, n in 1:Ns[city]
        ind = add_agent!(city, model, 0, :S) # Susceptible
    end
    # add infected individuals
    for city in 1:C
        inds = ids_in_position(city, model)
        for n in 1:Is[city]
            agent = model[inds[n]]
            agent.status = :I # Infected
            agent.days_infected = 1
        end
    end
    return model
end

# ╔═╡ 867f455d-39c0-458f-a374-72de9a95076e
function make_SIRgraphParams(;
	C,
    max_travel_rate,
    infection_period = 30,
    reinfection_probability = 0.05,
    detection_time = 14,
    death_rate = 0.02,
    Is = [zeros(Int, C - 1)..., 1],
    seed = 19,
)
	# For reproducibility
	Random.seed!(seed)

	# City population
    Ns = rand(50:5000, C)

	# Undetected transmission
    β_und = rand(0.3:0.02:0.6, C)

	# Detected transmission (10% of undetected)
    β_det = β_und ./ 10

	# Migrate from city i to city j
	# People in small cities tend to migrate to bigger cities
	migration_rates = zeros(C, C)
    for c in 1:C, c2 in 1:C
        migration_rates[c, c2] = (Ns[c] + Ns[c2]) / Ns[c]
    end

	# Normalize migration rates
	maxM = maximum(migration_rates)
    migration_rates = (migration_rates .* max_travel_rate) ./ maxM

	# Migrate to self = 1
    migration_rates[diagind(migration_rates)] .= 1.0

	return (; Ns,
        β_und,
        β_det,
        migration_rates,
        infection_period,
        reinfection_probability,
        detection_time,
        death_rate,
        Is)
end

# ╔═╡ a6471587-2890-4d4d-bb94-168993d021a1
SIRgraphparams = make_SIRgraphParams(C = 8, max_travel_rate = 0.01)

# ╔═╡ 340a98fc-05e0-4094-8944-9ef1d20bdc70
function migrate!(agent, model)
    pid = agent.pos
    d = DiscreteNonParametric(1:(model.C), model.migration_rates[pid, :])
    m = rand(model.rng, d)
    if m ≠ pid
        move_agent!(agent, m, model)
    end
end

# ╔═╡ 538e400e-5832-474b-aa32-7fc4d311577b
function transmit!(agent, model)
    agent.status == :S && return
    rate = if agent.days_infected < model.detection_time
        model.β_und[agent.pos]
    else
        model.β_det[agent.pos]
    end

    d = Poisson(rate)
    n = rand(model.rng, d)
    n == 0 && return

    for contactID in ids_in_position(agent, model)
        contact = model[contactID]
        if contact.status == :S ||
           (contact.status == :R && rand(model.rng) ≤ model.reinfection_probability)
            contact.status = :I
            n -= 1
            n == 0 && return
        end
    end
end

# ╔═╡ 1b2e2fb0-d5f1-4e47-a593-fddf9c10c8b8
update!(agent, model) = agent.status == :I && (agent.days_infected += 1)

# ╔═╡ 44745a79-2f35-4b2d-881c-8fbf07c918d7
function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand(model.rng) ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end

# ╔═╡ 557f060c-a9ac-4aee-85a5-b6ad48065489
function agent_step!(agent::PoorSoul, model)
    migrate!(agent, model)
    transmit!(agent, model)
    update!(agent, model)
    recover_or_die!(agent, model)
end

# ╔═╡ bcee6280-2b48-4bf5-a975-96637d1bfde2
total_infected(m) = count(a.status == :I for a in allagents(m))

# ╔═╡ 787286dd-f0bb-4278-8220-f81f7e87ab44
"Coloring infected node"
infected_fraction(xs) = cgrad(:inferno)[count(x.status == :I for x in xs) / length(xs)]

# ╔═╡ 86a11e20-3a6c-46ef-b062-f3de19d2db75
xpos = sinpi.(range(0, 2, length=9)[1:8])

# ╔═╡ 59383688-ea08-4a1b-ab9d-adc1bc1874ec
ypos = cospi.(range(0, 2, length=9)[1:8])

# ╔═╡ 0a003441-169a-4e2f-9c27-3c074d9606a2
begin
	model = make_SIRgraph(; SIRgraphparams...)
	anim = @animate for i in 1:40
		Agents.step!(model, agent_step!, 1)
		infected = total_infected(model)
		abm_plot_on_graph(model;
			size= (600, 600),
			ac = infected_fraction,
			curves=false,
			x=xpos, y=ypos,
			title = "Step $i: $infected infected")
	end

	mp4(anim, fps = 5)
end

# ╔═╡ 7ad1fc2b-d4d6-4431-a06a-e8fda045bee6
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 77180529-5223-4a0e-933b-447a3e9f56ae

# ╔═╡ Cell order:
# ╠═e5f4abe0-be15-11eb-038f-5da4c19a690d
# ╠═07c4c6aa-c7bb-4543-868d-975de9f37ba9
# ╠═99c0bed6-a533-4058-bf85-c81fac0f5976
# ╠═867f455d-39c0-458f-a374-72de9a95076e
# ╠═a6471587-2890-4d4d-bb94-168993d021a1
# ╠═557f060c-a9ac-4aee-85a5-b6ad48065489
# ╠═340a98fc-05e0-4094-8944-9ef1d20bdc70
# ╠═538e400e-5832-474b-aa32-7fc4d311577b
# ╠═1b2e2fb0-d5f1-4e47-a593-fddf9c10c8b8
# ╠═44745a79-2f35-4b2d-881c-8fbf07c918d7
# ╠═bcee6280-2b48-4bf5-a975-96637d1bfde2
# ╠═787286dd-f0bb-4278-8220-f81f7e87ab44
# ╠═86a11e20-3a6c-46ef-b062-f3de19d2db75
# ╠═59383688-ea08-4a1b-ab9d-adc1bc1874ec
# ╠═0a003441-169a-4e2f-9c27-3c074d9606a2
# ╠═7ad1fc2b-d4d6-4431-a06a-e8fda045bee6
# ╠═49674784-f070-488d-b53b-9cebf0ce101e
# ╠═77180529-5223-4a0e-933b-447a3e9f56ae

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
"""
    data = NodeData(sf::StackFrame, status::UInt8, span::UnitRange{Int})

Data associated with a single node in a flamegraph. `sf` is the stack frame
(see `?StackTraces.StackFrame`). `status` is a bitfield with information about this
node or any "suppressed" nodes immediately called by this one:

- `status & 0x01` is nonzero for runtime dispatch
- `status & 0x02` is nonzero for garbage collection

By default, C-language stackframes are omitted, but information about
their identity is accumulated into their caller's `status`.

`length(span)` is the number of times this stackframe was captured at this
depth and location in the flame graph. The starting index begins with the caller's
starting `span` but increments to ensure each child's `span` occupies a distinct
subset of the caller's `span`. Concretely, `span` is the range of indexes
that will be occupied by this stackframe when the flame graph is rendered.
"""
struct NodeData
    sf::StackFrame
    status::UInt8             # a bitfield, see below
    span::UnitRange{Int}
end

# status bitfield values
const runtime_dispatch = UInt8(1)
const gc_event         = UInt8(2)
const repl             = UInt8(4)

const defaultpruned = Tuple{Symbol,Symbol}[]

# This allows Revise to correct the location information in profiles
if VERSION >= v"1.5.0-DEV.9"
    # ref https://github.com/JuliaLang/julia/pull/34235
    lineinfodict(data::Vector{<:Unsigned}) = Profile.getdict(data)
else
    # Use the definition of Profile.getdict from Julia 1.5.0-DEV.9+
    function lineinfodict(data::Vector{<:Unsigned})
        # Lookup is expensive, so do it only once per ip.
        udata = unique(data)
        dict = Profile.LineInfoDict()
        for ip in udata
            st = Profile.lookup(convert(Ptr{Cvoid}, ip))
            # To correct line numbers for moving code, put it in the form expected by
            # Base.update_stackframes_callback[]
            stn = map(x->(x, 1), st)
            try Base.invokelatest(Base.update_stackframes_callback[], stn) catch end
            dict[UInt64(ip)] = map(first, stn)
        end
        return dict
    end
end

"""
    lidict = lineinfodict(uips)

Look up location information for each instruction pointer in `uips`.
This is a `Dict(UInt64=>Vector{StackFrame})`, where the `UInt64` is the instruction pointer.
The reason a single instruction pointer gives a `Vector{StackFrame}` is because
of inlining; the first entry corresponds to the instruction that actually ran,
and the later entries correspond to the call chain.

See also [`unique_ips`](@ref).
"""
lineinfodict(s::Set) = lineinfodict(collect(s))

"""
    g = flamegraph(data=Profile.fetch(); lidict=nothing, C=false, combine=true, recur=:off, norepl=true, pruned=[], filter=nothing)

Compute a graph representing profiling data. To compute it for the currently-collected
profiling information, omit both `data` and `lidict`; if you are computing it for saved
profiling data, supply both. (`data` and `lidict` must be a matched pair from `Profile.retrieve()`.)

You can control the strategy with the following keywords:

- `C`: if `true`, include stackframes collected from `ccall`ed code.
- `recur` (supported on Julia 1.4+): represent recursive calls as if they corresponded to
  iteration.
- `norepl`: if true, the portions of stacktraces deeper than `REPL.eval_user_input` are discarded.
- `pruned`: a list of `(funcname, filename)` pairs that trigger the termination of this branch
  of the flame graph. You can use this to prevent very "tall" graphs from deeply-recursive
  calls, e.g., `pruned = [("sort!", "sort.jl")]` would omit nodes corresponding to Julia's
  `sort!` function and anything called by it. See also `recur` for an alternative strategy.
- `combine`: if true, instruction pointers that correspond to the same line of code are
  combined into a single stackframe
- `filter`: drop all branches that do not satisfy the filter condition. `Condition` can be `string`, `regex`
  or any function, that can be applied to `NodeData`. For example, `filter = "mapslices"` or equivalently
  `filter = x -> (x.sf.func == :mapslices)` removes all branches that do not contain `mapslices` Node as a child
  or ancestor.

`g` can be inspected using [`AbstractTrees.jl`'s](https://github.com/JuliaCollections/AbstractTrees.jl)
`print_tree`.
"""
function flamegraph(data=Profile.fetch(); lidict::Union{Dict{UInt64,Vector{Base.StackTraces.StackFrame}},Dict{UInt64,Base.StackTraces.StackFrame},Nothing}=nothing, C::Bool=false, combine::Bool=true,
        recur::Symbol=:off, norepl::Bool=true, pruned=defaultpruned, filter=nothing)
    if lidict === nothing
        lidict = lineinfodict(unique(data))
    end
    data_u64 = convert(Vector{UInt64}, data)
    root = combine ? StackFrameTree{StackFrame}() : StackFrameTree{UInt64}()
    # Build the tree with C=true, regardless of user setting. This is because
    # we need the C frames to set the status flag. They will be omitted by `flamegraph!`
    # as needed.
    if VERSION >= v"1.8.0-DEV.460"
        root, _ = Profile.tree!(root, data_u64, lidict, #= C =# true, recur)
    elseif VERSION >= v"1.4.0-DEV.128"
        root = Profile.tree!(root, data_u64, lidict, #= C =# true, recur)
    else
        root = Profile.tree!(root, data_u64, lidict, #= C =# true)
    end
    if isempty(root.down)
        Profile.warning_empty()
        return nothing
    end
    root.count = sum(pr->pr.second.count, root.down)  # root count seems borked
    g = flamegraph!(Node(NodeData(root.frame, status(root, C), 1:root.count)), root; C=C, pruned=pruned)
    if norepl
        prunerepl!(g)
    end

    if filter != nothing
        filtergraph!(g, filter)
    end

    return g
end

function status(node, C::Bool)
    st = status(node.frame)
    C && return st
    # If we're suppressing C frames, check all C-frame children
    for child in values(node.down)
        child.frame.from_c || continue
        st |= status(child, C)
    end
    return st
end

function status(sf::StackFrame)
    st = UInt8(0)
    if sf.from_c && (sf.func === :jl_invoke || sf.func === :jl_apply_generic || sf.func === :ijl_apply_generic)
        st |= runtime_dispatch
    end
    if sf.from_c && startswith(String(sf.func), "jl_gc_")
        st |= gc_event
    end
    if !sf.from_c && sf.func === :eval_user_input && endswith(String(sf.file), "REPL.jl")
        st |= repl
    end
    return st
end

function nodematches(needle::Union{AbstractString,Regex,AbstractChar}, g::Node)
    return occursin(needle, string(g.data.sf))
end

function nodematches(needle::Function, g::Node)
    return needle(g.data)
end

function canfind(needle, g::Node, memo = Dict{Node, Bool}())
    haskey(memo, g) && return memo[g]
    if nodematches(needle, g)
        memo[g] = true
    else
        memo[g] = any(canfind(needle, child, memo) for child in g)
    end
    return memo[g]
end

function filtergraph!(g::Node, filter, memo = Dict{Node, Bool}())
    if !canfind(filter, g, memo)
        if isroot(g)
            @warn "The filter condition results in the root node pruning, so the filter is ignored"
        else
            prunebranch!(g)
        end
    else
        for child in g
            nodematches(filter, child) && continue
            filtergraph!(child, filter, memo)
        end
    end
end

function flamegraph!(graph, ptree; C=false, pruned=defaultpruned, hstart=first(graph.data.span))
    nexts = collect(values(ptree.down))
    lilist = collect(frame.frame for frame in nexts)
    p = Profile.liperm(lilist)
    for i in p
        down = nexts[i]
        frame, count = down.frame, down.count
        ispruned(frame, pruned) && continue
        if !C && frame.from_c
            flamegraph!(graph, down; C=C, pruned=pruned, hstart=hstart)
        else
            child = addchild(graph, NodeData(frame, status(down, C), hstart:hstart+count-1))
            flamegraph!(child, down; C=C, pruned=pruned)
        end
        hstart += count
    end
    return graph
end

function ispruned(frame, (fname, file)::Tuple{Any,Any})
    return frame.func == Symbol(fname) && frame.file == Symbol(file)
end

ispruned(frame, pruned) = any(t->ispruned(frame, t), pruned)

function prunerepl!(node)
    for c in node
        if c.data.status & repl != 0
            # Add children directly to the root node
            root = node
            replchild = c
            while !isroot(root)
                root = root.parent
                replchild = replchild.parent
            end
            graftchildren!(root, c)
            # Eliminate all nodes in between. This might include some that don't
            # call REPL code, but as this is also internal it seems OK.
            prunebranch!(replchild)
            return true
        else
            prunerepl!(c) && return true
        end
    end
    return false
end

ext: jl


text:
"""

	| F | ≤ 1

"""

struct AtMostOne <: Constraint

	F::Variable

end

### MDR


text:
using FFTW

"""
    plan_transforms(::PPN, A::Array; planner_flag=FFTW.PATIENT)

Return the transforms required to solve Poisson's equation with periodic
boundary conditions in the x,y dimensions and staggered Neumann boundary
conditions in the z-dimension.

Fast Fourier transforms (FFTs) are used in the x,y dimensions and real-to-real
discrete cosine transforms are used in the z dimension. Note that the DCT-II is
used for the DCT and the DCT-III for the IDCT which correspond to REDFT10 and
REDFT01 in FFTW.

They operatore on an array with the shape of `A`, which is needed to plan
efficient transforms. `A` will be mutated.

Note that the transforms returns operate on Arrays and so only work on CPUs.
"""
function plan_transforms(::PPN, A::Array, planner_flag=FFTW.PATIENT)
    FFT!  = FFTW.plan_fft!(A, [1, 2]; flags=planner_flag)
    IFFT! = FFTW.plan_ifft!(A, [1, 2]; flags=planner_flag)
    DCT!  = FFTW.plan_r2r!(A, FFTW.REDFT10, 3; flags=planner_flag)
    IDCT! = FFTW.plan_r2r!(A, FFTW.REDFT01, 3; flags=planner_flag)
    return FFT!, IFFT!, DCT!, IDCT!
end

"""
    plan_transforms(::PPN, A::Array; planer_flag=FFTW.PATIENT)

Similar to `plan_transforms(::PPN, ...)` but return the transforms required to
solve Poisson's equation with periodic boundary conditions in the x dimension
and staggered Neumann boundary conditions in the y,z dimensions.
"""
function plan_transforms(::PNN, A::Array, planner_flag=FFTW.PATIENT)
    FFT!  = FFTW.plan_fft!(A, 1; flags=planner_flag)
    IFFT! = FFTW.plan_ifft!(A, 1; flags=planner_flag)
    DCT!  = FFTW.plan_r2r!(A, FFTW.REDFT10, [2, 3]; flags=planner_flag)
    IDCT! = FFTW.plan_r2r!(A, FFTW.REDFT01, [2, 3]; flags=planner_flag)
    return FFT!, IFFT!, DCT!, IDCT!
end

struct PoissonSolverCPU{BC, AAR, AAC, FFTT, DCTT, IFFTT, IDCTT} <: AbstractPoissonSolver
        bcs :: BC
        kx² :: AAR
        ky² :: AAR
        kz² :: AAR
    storage :: AAC
       FFT! :: FFTT
       DCT! :: DCTT
      IFFT! :: IFFTT
      IDCT! :: IDCTT
end

function PoissonSolverCPU(pbcs::PoissonBCs, grid::AbstractGrid, planner_flag=FFTW.PATIENT)
    Nx, Ny, Nz, _ = unpack_grid(grid)

    # The eigenvalues of the discrete form of Poisson's equation correspond
    # to discrete wavenumbers so we call them k² to make the analogy with
    # solving Poisson's equation in spectral space.
    kx² = λi(grid, pbcs)
    ky² = λj(grid, pbcs)
    kz² = λk(grid, pbcs)

    # Storage for RHS and Fourier coefficients is hard-coded to be Float64
    # because of precision issues with Float32.
    # See https://github.com/climate-machine/Oceananigans.jl/issues/55
    storage = zeros(Complex{Float64}, grid.Nx, grid.Ny, grid.Nz)

    FFT!, IFFT!, DCT!, IDCT! = plan_transforms(pbcs, storage, planner_flag)

    PoissonSolverCPU(pbcs, kx², ky², kz², storage, FFT!, DCT!, IFFT!, IDCT!)
end

normalize_idct_output(::PPN, grid::AbstractGrid, ϕ::Array) = (@. ϕ = ϕ / (2grid.Nz))
normalize_idct_output(::PNN, grid::AbstractGrid, ϕ::Array) = (@. ϕ = ϕ / (4grid.Ny*grid.Nz))

"""
    solve_poisson_3d!(solver::PoissonSolverCPU, grid::RegularCartesianGrid)

Solve Poisson equation on a uniform staggered grid (Arakawa C-grid) with
appropriate boundary conditions as specified by `solver.bcs` using planned FFTs
and DCTs. The right-hand-side RHS is stored in solver.storage which the solver
mutates to produce the solution, so it will also be stored in solver.storage.

We should describe the algorithm in detail in the documentation.
"""
function solve_poisson_3d!(solver::PoissonSolverCPU, grid::RegularCartesianGrid)
    Nx, Ny, Nz, _ = unpack_grid(grid)
    kx², ky², kz² = solver.kx², solver.ky², solver.kz²

    # We can use the same storage for the RHS and the solution ϕ.
    RHS, ϕ = solver.storage, solver.storage

    solver.DCT! * RHS  # Calculate DCTᶻ(f) in place.
    solver.FFT! * RHS  # Calculate FFTˣʸ(f) in place.

    # Solve the discrete Poisson equation in spectral space. We are essentially
    # computing the Fourier coefficients of the solution from the Fourier
    # coefficients of the RHS.
    @. ϕ = -RHS / (kx² + ky² + kz²)

    # Setting DC component of the solution (the mean) to be zero. This is also
    # necessary because the source term to the Poisson equation has zero mean
    # and so the DC component comes out to be ∞.
    ϕ[1, 1, 1] = 0

    solver.IFFT! * ϕ  # Calculate IFFTˣʸ(ϕ) in place.
    solver.IDCT! * ϕ  # Calculate IDCTᶻ(ϕ) in place.

    # Must normalize by 2N for each dimension transformed via FFTW.REDFT.
    normalize_idct_output(solver.bcs, grid, ϕ)

    nothing
end


text:
module Quadprog
# written by Stefan Harmeling, 2014-11-12

export quadprog

using Ipopt  # a julia wrapper for Ipopt

# References for Ipopt:
# A. Wächter and L. T. Biegler,
# On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming,
# Mathematical Programming 106(1), pp. 25-57, 2006

function quadprog(Q, c;
                  A   =  Array(Float64, (0, length(c))),
                  b   =  Array(Float64, 0),
                  Aeq =  Array(Float64, (0, length(c))),
                  beq =  Array(Float64, 0),
                  lb  = -Inf*ones(length(c)),
                  ub  =  Inf*ones(length(c)),
                  x0  =  zeros(length(c)))

    # DESCRIPTION:
    # min  0.5 * x' * Q * x + c' * x
    # s.t. A   * x <= b
    #      Aeq * x == beq
    #      lb <= x <= ub
    # init x with x0

    if ~issym(Q)           error("Q must be symmetric") end
    n   = length(x0)       # number of variables
    m   = size(A,   1)     # number of inequality constraints
    meq = size(Aeq, 1)     # number of equality constraints
    if length(b)   != m    error("size missmatch") end
    if length(beq) != meq  error("size missmatch") end
    if size(A, 2) != size(Aeq, 2) != length(ub) != length(lb) != length(x0) != n
        error("size missmatch")
    end

    # setup the Ipopt problem
    g_L = [-Inf*ones(m), beq]   # lower constraint bounds
    g_U = [b, beq]              # upper constraint bounds
    eval_f(x) = (0.5*x'*Q*x + c'*x)[1]
    eval_grad_f(x, grad_f) = begin grad_f[:] = Q*x + c end
    eval_g(x, g) = begin g[:] = [A*x, Aeq*x] end
    eval_jac_g(x, mode, rows, cols, values) = begin
        # grad_g = [A; Aeq]
        if mode == :Structure
            rows[:] = (1:(m+meq)) * ones(Int64, n)'
            cols[:] = ones(Int64, m+meq) * (1:n)'
        else
            values[:] = [A; Aeq]
        end
    end
    eval_h(x, mode, rows, cols, obj_factor, lambda, values) = begin
        if mode == :Structure
            rows[:] = [1:n] * ones(Int64, n)'
            cols[:] = ones(Int64, n) * [1:n]'
        else
            values[:] = -obj_factor * Q
        end
    end
    prob = createProblem(n, lb, ub, m+meq, g_L, g_U, (m+meq)*n, n*n,
                         eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
    addOption(prob, "print_level", 0)
    addOption(prob, "tol", 1e-12)
    prob.x = x0
    status = solveProblem(prob)
    return (prob.x, prob.obj_val, status)
end

end # of module


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, SpecialFunctions, LinearAlgebra, Test

@testset "broadcast" begin
    @testset "Special functions" begin
        x = Fun()
        @test exp(x) ≈ exp.(x) atol=10eps()
        f = Fun(exp)
        @test exp.(f) ≈ exp(f)
        @test besselj.(1,f) ≈ besselj(1,f)

        @test atan.(f,1)(0.1) ≈ atan(f(0.1),1)
        @test atan.(f,f)(0.1) ≈ atan(f(0.1),f(0.1))
    end

    @testset  "triplet with mixed scalar-fun" begin
        f = Fun(exp)
        @test ((a,b,c) -> exp(a +b+c)).(f,2.0,f) ≈ exp(2f + 2.0)
    end

    @testset "Arrays" begin
        x = Fun()
        @test exp.( x .+ [1,2,3]) isa Vector
        @test exp.(x .+ [1,2]) ≈ [exp(x+1),exp(x+2)]

        @test maximum.( [x;x] .+ [x;x]) isa Fun
        @test maximum.( [x;x] .+ [x;x]) ≈ 2x

        @test exp.([x,x] .+ 2 .+ [x;x]) isa Vector
        @test exp.( [x,x] .+ 2 .+ [x;x]) ≈ [exp(2x+2),exp(2x+2)]

        # since array broadcasting takes presidence, the following
        # does not try to call a constructor:
        @test ApproxFunBase.tocanonical.(x,[0.1,0.2]) ≈ [0.1,0.2]
    end

    @testset "broadcast!" begin
        x = Fun()
        f = Fun()
        f .= exp.(x)
        @test ≈(exp(x),f;atol=10eps())

        f = Fun(Ultraspherical(1))
        f .= exp.(x)
        @test f(0.1) ≈ exp(0.1)
    end

    @testset "norm" begin
        f = Fun(x ->[exp(x),sin(x)])
        F = [Fun(exp),Fun(sin)]
        @test norm.(f) ≈ Fun(x ->norm([exp(x),sin(x)]))
        @test norm.(F) ≈ [norm(Fun(exp)),norm(Fun(sin))] # isa Vector{Float64}

        @test norm.(f,[1,2]) ≈ norm.(F,[1,2]) ≈ [norm(f[1],1),norm(f[2],2)] # isa Vector{Float64}
    end
end


text:
using ImageIO
using SubpixelRegistration
using Test
using TestImages
using StableRNGs

rng = StableRNG(121)

source = Float64.(testimage("cameraman"))

@testset "SubpixelRegistration.jl" begin

    @testset "pixel shift" for _ in 1:100
        shift = (rand(rng, -25:25), rand(rng, -25:25))
        shifted = @inferred fourier_shift(source, shift)
        result = @inferred phase_offset(source, shifted; upsample_factor = 1)
        @test all(result.shift .≈ -1 .* shift)

        registered = @inferred register(source, shifted; upsample_factor = 1)
        @test registered ≈ fourier_shift(shifted, result.shift, result.phasediff)
    end

    @testset "subpixel shift" for _ in 1:100
        shift = 10 .* randn(rng, 2)
        shifted = fourier_shift(source, shift)

        f = 10
        result = @inferred phase_offset(source, shifted; upsample_factor = f)
        @test all(isapprox.(result.shift, -1 .* shift, atol = inv(f)))

        registered = @inferred register(source, shifted; upsample_factor = f)
        @test registered ≈ fourier_shift(shifted, result.shift, result.phasediff)
    end

    @testset "coregister" begin
        shifts = 10 .* randn(rng, 19, 2)
        cube = cat(source, (fourier_shift(source, shift) for shift in eachrow(shifts))..., dims = 3)
        f = 100
        cube_shift = @inferred coregister(cube; dims = 3, upsample_factor = f)
        for slice in eachslice(cube_shift, dims = 3)
            @test all(isapprox.(slice, source; atol = 1e-2))
        end

        # alternate axis of iteration
        cubep = permutedims(cube, (3, 2, 1))
        cubep_shift = @inferred coregister(cubep, dims = 1, upsample_factor = f)
        for slice in eachslice(cubep_shift, dims = 1)
            @test all(isapprox.(slice', source; atol = 1e-2))
        end
    end

end


text:
using StochasticDiffEq

f(u,p,t,W) = 1.01u.+0.87u.*W
u0 = 1.00
tspan = (0.0,1.0)
prob = RODEProblem(f,u0,tspan)
sol = solve(prob,RandomEM(),dt=1/100)

f(u,p,t,W,du) = (du.=1.01u.+0.87u.*W)
u0 = ones(4)
prob = RODEProblem(f,u0,tspan)
sol = solve(prob,RandomEM(),dt=1/100)

f(u,p,t,W) = 2u*sin(W)
u0 = 1.00
tspan = (0.0,5.0)
prob = RODEProblem{false}(f,u0,tspan)
sol = solve(prob,RandomEM(),dt=1/100)

function f(du,u,p,t,W)
  du[1] = 2u[1]*sin(W[1] - W[2])
  du[2] = -2u[2]*cos(W[1] + W[2])
end
u0 = [1.00;1.00]
tspan = (0.0,5.0)
prob = RODEProblem(f,u0,tspan)
sol = solve(prob,RandomEM(),dt=1/100)

function f(du,u,p,t,W)
  du[1] = -2W[3]*u[1]*sin(W[1] - W[2])
  du[2] = -2u[2]*cos(W[1] + W[2])
end
u0 = [1.00;1.00]
tspan = (0.0,5.0)
prob = RODEProblem(f,u0,tspan,rand_prototype=zeros(3))
sol = solve(prob,RandomEM(),dt=1/100)


text:
"""
Main module for `LLLplus.jl` -- lattice reduction and related tools for Julia.

As an example of the functions in the package, see [`lll`](@ref), which does
Lenstra–Lenstra–Lovász lattice reduction of a matrix.

"""
module LLLplus

using LinearAlgebra
using Printf

export
    lll,
    cvp,
    svp,
    hkz,ishkzreduced,
    brun,
    gauss,
    sizereduction,
    seysen,
    vblast,
    subsetsum,lagariasodlyzko,mdsubsetsum,
    integerfeasibility,
    rationalapprox,spigotBBP,
    hardsphere, hard_sphere,
    issizereduced,islllreduced,orthogonalitydefect,hermitefactor,seysencond,
    gen_qary_b,gen_qary!,
    dataTypeForGram,intTypeGivenBitsRequired

include("lll.jl")          # lll, gauss, sizereduction
include("l2.jl")           # lll on Rational{Integer} with `l2`
include("cvp.jl")          # cvp, svp
include("hkz.jl")
include("bkz.jl")
include("brun.jl")
include("seysen.jl")
include("vblast.jl")
include("applications.jl") # subsetsum, integerfeasibility, rationalapprox
include("utilities.jl")
include("latticegen.jl")

include("hard_sphere.jl")  # may be deprecated in future

end # LLLplus module


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

mutable struct WinRMConfiguration <: SwaggerModel
    listeners::Any # spec type: Union{ Nothing, Vector{WinRMListener} } # spec name: listeners

    function WinRMConfiguration(;listeners=nothing)
        o = new()
        validate_property(WinRMConfiguration, Symbol("listeners"), listeners)
        setfield!(o, Symbol("listeners"), listeners)
        o
    end
end # type WinRMConfiguration

const _property_map_WinRMConfiguration = Dict{Symbol,Symbol}(Symbol("listeners")=>Symbol("listeners"))
const _property_types_WinRMConfiguration = Dict{Symbol,String}(Symbol("listeners")=>"Vector{WinRMListener}")
Base.propertynames(::Type{ WinRMConfiguration }) = collect(keys(_property_map_WinRMConfiguration))
Swagger.property_type(::Type{ WinRMConfiguration }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_WinRMConfiguration[name]))}
Swagger.field_name(::Type{ WinRMConfiguration }, property_name::Symbol) =  _property_map_WinRMConfiguration[property_name]

function check_required(o::WinRMConfiguration)
    true
end

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


text:
module TestOpenml

using Test
using MLJOpenML
using ScientificTypes

for M in [OpenML, MLJOpenML]
    @test length(M.list_tags()) > 200
    table = M.list_datasets(
        tag = "OpenML100",
        filter = "number_instances/100..1000/number_features/1..10")
    @test issubset([:id, :name, :status], schema(table).names)
    M.describe_dataset(61)
    table = M.load(61)
    @test :sepallength in schema(table).names
end

end
true


text:
x, xi  = 2.0, 0.5
y, yi  = 4.0, 0.25
z, zi  = x + y, 1.0 / ( x + y )

multiplier = (n1, n2) -> (m) -> n1 * n2 * m

numlist  = [x ,  y,  z]
numlisti = [xi, yi, zi]

@show collect(multiplier(n, invn)(0.5) for (n, invn) in zip(numlist, numlisti))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
const BroadMod = Broadcast.Broadcasted{Broadcast.ArrayStyle{ModelingStruct}}

Base.BroadcastStyle(::Type{<:ModelingStruct}) = Broadcast.ArrayStyle{ModelingStruct}()

# Goes through fields and reattaches views to new vector
function _review(v, ms::ModelingStruct{T,V,K,P}) where {T,V,K,P}
    tup = map(x -> _review(v, x), _structure(ms))
    return ModelingStruct(v, NamedTuple{K}(tup), Val(false))
end
_review(v, a::AbstractVector) = map(x -> _review(v, x), a)
_review(v::A, r::SubArray) where {A} = view(v, r.indices[1])

# Create uninitialized ModelingStruct at new memory location
Base.similar(ms::ModelingStruct) = _review(similar(_data(ms)), ms)
Base.similar(ms::ModelingStruct, ::Type{T}) where {T} = _review(similar(_data(ms),T), ms)
Base.similar(bc::BroadMod, ::Type{T}) where T = similar(bc.args[1], T)

# Needed for usage with LAPACK. Stolen from:
#   https://github.com/JuliaDiffEq/LabelledArrays.jl/blob/master/src/larray.jl
function Base.unsafe_convert(::Type{Ptr{T}}, ms::ModelingStruct{T}) where T
    return Base.unsafe_convert(Ptr{T}, _data(ms))
end

# Conversion to NamedTuple (note, does not preserve numeric types of original NamedTuple)
function _namedtuple(ms::ModelingStruct{T,V,K}) where {T,V,K}
    data = []
    for key in K
        val = getproperty(ms, key) |> _namedtuple
        push!(data, key => val)
    end
    return (; data...)
end
_namedtuple(v::AbstractVector) = _namedtuple.(v)
_namedtuple(x) = x

Base.convert(::Type{NamedTuple}, x::ModelingStruct) = _namedtuple(x)
Base.NamedTuple(x::ModelingStruct) = _namedtuple(x)

max_stars_repo_name: jonniedie/NamedViewVectors


text:
import Canvas, HTTP, JSON, Test

import Revise, Canvas
auth=Canvas.OAuth2(ENV["CANVAS_TOKEN"])
api = Canvas.CanvasAPI(ENV["CANVAS_API"]; auth=auth)
Canvas.setapi!(api)

cs, pd = Canvas.courses()
c = cs[1]
@assert c.id == 7609

c′ = Canvas.course(c)
c′′ = Canvas.course(7609)
c = Canvas.course(7609; query=Dict("include[]"=>["permissions"]))

fs, pd = Canvas.files(7609)


text:
# * Source terms

abstract type AbstractSourceTerm end

# ** Dipole

struct DipoleSourceTerm{Dipole<:Function,Field<:Union{ElectricFields.AbstractField,AbstractVector}} <: AbstractSourceTerm
    d::Dipole
    F::Field

    DipoleSourceTerm(d::Dipole, F::Field) where {Dipole<:Function, Field<:ElectricFields.AbstractField} =
        new{Dipole,Field}(d, F)

    function DipoleSourceTerm(d::Dipole, F::Field, fs::Number) where {Dipole<:Function, Field<:ElectricFields.AbstractField}
        t = timeaxis(F, fs)
        Fv = field_amplitude.(F, t)
        new{Dipole,typeof(Fv)}(d, Fv)
    end

    DipoleSourceTerm(d::Dipole, Fv::FieldAmplitude, Av::AbstractVector, t::AbstractRange) where {Dipole<:Function, FieldAmplitude<:AbstractVector} =
        new{Dipole,FieldAmplitude}(d, Fv)
end

dipole(F::Number, d::Number) = F*d
dipole(F::SVector{3}, d::SVector{3}) = dot(F, d)
dipole(F::Number, d::SVector{3}) = F*d[3]

source_term(d::DipoleSourceTerm{<:Any,<:ElectricFields.AbstractField}, t, p) =
    dipole(field_amplitude(d.F, t), d.d(p))

source_term(d::DipoleSourceTerm{<:Any,<:AbstractVector}, i::Integer, p) =
    dipole(d.F[i], d.d(p))

# * Ionizations channels

"""
    IonizationChannel(E, st)

Represents an ionization channel with energy `E` above the neutral
(which is taken to have energy `0 Ha`).

"""
struct IonizationChannel{T,SourceTerm<:AbstractSourceTerm}
    E::T
    st::SourceTerm
end

source_term(ic::IonizationChannel, args...) =
    source_term(ic.st, args...)

IonizationChannel(E, args...) =
    IonizationChannel(austrip(E), DipoleSourceTerm(Base.Fix2(d_hyd, hyd_α(austrip(E))), args...))

Base.show(io::IO, ic::IonizationChannel) =
    write(io, "IonizationChannel: Iₚ = $(ic.E) Ha = $(27.211ic.E) eV")

# * Couplings

abstract type AbstractCanonicalMomentumConservation end
struct CanonicalMomentumConservation <:  AbstractCanonicalMomentumConservation end
struct NoCanonicalMomentumConservation <:  AbstractCanonicalMomentumConservation end

abstract type AbstractCoupling end
Base.iszero(::AbstractCoupling) = false
canonical_momentum_conservation(::AbstractCoupling) = NoCanonicalMomentumConservation()
canonical_momentum_conservation(::Type{<:AbstractCoupling}) = NoCanonicalMomentumConservation()

struct NoCoupling <: AbstractCoupling end
Base.iszero(::NoCoupling) = true
Base.zero(::AbstractCoupling) = NoCoupling()

# ** Dipole couplings

struct DipoleCoupling{DipoleMoment<:Union{Number,SVector{3}},Field<:ElectricFields.AbstractField} <: AbstractCoupling
    d::DipoleMoment
    F::Field
end

canonical_momentum_conservation(::DipoleCoupling) = CanonicalMomentumConservation()
canonical_momentum_conservation(::Type{<:DipoleCoupling}) = CanonicalMomentumConservation()

# ** Coulomb couplings

momentum_transfer(k::Number, p::Number) = k-p
momentum_transfer(k::SVector{3}, p::SVector{3}) = k-p
momentum_transfer(k::SVector{3}, p::Number) = SVector{3}(k[1], k[2], k[3]-p)

struct CoulombCoupling{Coupling} <: AbstractCoupling
    coupling::Coupling
end

(cc::CoulombCoupling)(𝐤, 𝐩, _) = cc.coupling(𝐤, 𝐩)


text:
make_keys(ssh_keygen_file = settings("ssh_keygen_file")) = mktempdir() do temp
    cd(temp) do
        @info "Generating ssh key"
        filename = ".documenter"
        try
            run(`$ssh_keygen_file -f $filename -N ""`)
        catch x
            if isa(x, Base.UVError)
                error("Cannot find $ssh_keygen_file")
            else
                rethrow()
            end
        end
        read(string(filename, ".pub"), String),
            read(filename, String) |> chomp |> base64encode
    end
end


text:
"""
An infinite length data stream generator based on a random decision tree

# Arguments:
* `stream::Stream`: A container of data stream options
* `classification_function`: One of the 4 functions `sine_func_[0, 1, 2, 3]`.
    Determines how the label is generated from features.

"""
@enum NodeType begin
    numeric
    categoric
end
struct RTreeNode
    class_label::Union{Int, Missing}
    type::NodeType
    split_att_index::Int
    split_att_value::Float64
    children::Array{RTreeNode}
end

struct RandomTreeGenerator <: AbstractDataStream
    stream::Stream
    tree_random_seed::Int
    tree_random_state::AbstractRNG
    n_categories_per_cat_feature::Int
    max_leaf_depth::Int
    min_leaf_depth::Int
    fraction_leaves_per_level::Float64
    tree_root::RTreeNode
    cat_var_selector::Base.RefValue{Int}
end

"""
Initialize a RandomTreeGenerator using an Int to select a classification_function.
"""
function RandomTreeGenerator(;tree_random_state::Int=-1,
                        sample_random_state::Int=-1,
                        n_classes::Int=2,
                        n_cat_features=5,
                        n_num_features=5,
                        n_categories_per_cat_feature=5,
                        max_leaf_depth=5,
                        min_leaf_depth=3,
                        fraction_leaves_per_level=0.15
                        )
    if tree_random_state == -1
        tree_random_state = Random.rand(1:1000)
    end
    if sample_random_state == -1
        sample_random_state = Random.rand(1:1000)
    end
    rng = MersenneTwister(sample_random_state)
    tree_rng = MersenneTwister(tree_random_state)

    return RandomTreeGenerator(Stream(n_samples=-1,
                                    n_targets=1,
                                    n_features=(n_cat_features * n_categories_per_cat_feature)+n_num_features,
                                    n_num_features=n_num_features,
                                    n_cat_features=n_cat_features,
                                    n_classes=n_classes,
                                    n_cat_features_idx=[i for i in n_num_features+1:n_num_features + (n_cat_features * n_categories_per_cat_feature)],
                                    current_sample_x=Array{Float64, 1}(),
                                    current_sample_y=Array{Int, 1}(),
                                    sample_idx=0,
                                    feature_names=vcat(["f$(i)" for i in 1:n_num_features], ["c$(i)$(v)" for i in n_num_features+1:n_num_features + n_cat_features for v in 1:n_categories_per_cat_feature]),
                                    target_names=["label"],
                                    target_values=[0, 1],
                                    name="RandomTreeGenerator",
                                    random_state_seed=sample_random_state,
                                    random_state=rng),
                                    tree_random_state,
                                    tree_rng,
                                    n_categories_per_cat_feature,
                                    max_leaf_depth,
                                    min_leaf_depth,
                                    fraction_leaves_per_level,
                                    _generate_random_tree(
                                        tree_rng,
                                        n_classes,
                                        n_cat_features,
                                        n_num_features,
                                        n_categories_per_cat_feature,
                                        max_leaf_depth,
                                        min_leaf_depth,
                                    fraction_leaves_per_level),
                                        Ref(0)
                                )
end

function _generate_random_tree(tree_rng::AbstractRNG,
                                n_classes::Int,
                                n_cat_features::Int,
                                n_num_features::Int,
                                n_categories_per_cat_feature::Int,
                                max_leaf_depth::Int,
                                min_leaf_depth::Int,
                                fraction_leaves_per_level::Float64)
    _generate_random_tree_node(0, [i for i in 1:n_num_features+n_cat_features], [0.0 for i in 1:n_num_features], [1.0 for i in 1:n_num_features],
                            tree_rng, n_classes, n_cat_features, n_num_features, n_categories_per_cat_feature, max_leaf_depth, min_leaf_depth, fraction_leaves_per_level)
end

function  _generate_random_tree_node(current_depth::Int,
                                    nominal_candidates::Array{Int},
                                    min_numeric_values::Array{Float64},
                                    max_numeric_values::Array{Float64},
                                    tree_rng::AbstractRNG,
                                    n_classes::Int,
                                    n_cat_features::Int,
                                    n_num_features::Int,
                                    n_categories_per_cat_feature::Int,
                                    max_leaf_depth::Int,
                                    min_leaf_depth::Int,
                                    fraction_leaves_per_level::Float64)
    leaf_rng = rand(tree_rng) > fraction_leaves_per_level
    set_leaf = (current_depth >= max_leaf_depth) || (current_depth >= min_leaf_depth && leaf_rng)
    if set_leaf
        class_label = rand(tree_rng, 1:n_classes)
        leaf = RTreeNode(class_label, numeric, -1, -1, [])
        return leaf
    end

    split_attribute = rand(tree_rng, 1:length(nominal_candidates))
    if split_attribute <= n_num_features
        min_value = min_numeric_values[split_attribute]
        max_value = max_numeric_values[split_attribute]
        split_point = min_value + (max_value - min_value) * rand(tree_rng)
        left_child = _generate_random_tree_node(
            current_depth + 1,
            nominal_candidates,
            min_numeric_values,
            [max_numeric_values[1:split_attribute]..., split_point, max_numeric_values[split_attribute+2:length(max_numeric_values)]...],
            tree_rng,
            n_classes,
            n_cat_features,
            n_num_features,
            n_categories_per_cat_feature,
            max_leaf_depth,
            min_leaf_depth,
            fraction_leaves_per_level
        )
        right_child = _generate_random_tree_node(
            current_depth + 1,
            nominal_candidates,
            [min_numeric_values[1:split_attribute]..., split_point, min_numeric_values[split_attribute+2:length(min_numeric_values)]...],
            max_numeric_values,
            tree_rng,
            n_classes,
            n_cat_features,
            n_num_features,
            n_categories_per_cat_feature,
            max_leaf_depth,
            min_leaf_depth,
            fraction_leaves_per_level
        )
        return RTreeNode(missing, numeric, split_attribute, split_point, [left_child, right_child])
    else
        new_nominal_candidates = [nominal_candidates[1:split_attribute]..., nominal_candidates[split_attribute+2:length(nominal_candidates)]...]
        children = [_generate_random_tree_node(
            current_depth + 1,
            new_nominal_candidates,
            min_numeric_values,
            max_numeric_values,
            tree_rng,
            n_classes,
            n_cat_features,
            n_num_features,
            n_categories_per_cat_feature,
            max_leaf_depth,
            min_leaf_depth,
            fraction_leaves_per_level) for i in 1:n_categories_per_cat_feature]
            return RTreeNode(missing, categoric, split_attribute, 0, children)
    end
end

function _classify_instance(node, X, n_num_features, n_categories_per_cat_feature)
    if !ismissing(node.class_label)
        return node.class_label
    end
    attribute_value = X[node.split_att_index]
    if node.type == numeric
        child_id = attribute_value < node.split_att_value ? 1 : 2
        return _classify_instance(node.children[child_id], X, n_num_features, n_categories_per_cat_feature)
    else
        child_id = get_cat_value_from_onehot(X, node.split_att_index, n_num_features, n_categories_per_cat_feature)
        return _classify_instance(node.children[child_id], X, n_num_features, n_categories_per_cat_feature)
    end
end

"""
Each input vector X is laid out [N1, N2, ... C1_A1, C1_A2...]
I.e all numeric attributes are first then all categoric.
Each categoric attribute is one hot encoded, so covers
indexes [i, i+num_levels] in X.
We scan this range for a 1 to find the true value.
"""
function get_cat_value_from_onehot(X, attribute_index, n_num_features, n_categories_per_cat_feature)
    num_prev_cat_attrs = attribute_index - n_num_features - 1
    starting_index = (n_num_features + 1) + num_prev_cat_attrs * n_categories_per_cat_feature
    for i in starting_index:starting_index+n_categories_per_cat_feature - 1
        if X[i] == 1.0
            return 1 + (i - starting_index)
        end
    end
end

function classify_instance(s::RandomTreeGenerator, X)
    _classify_instance(s.tree_root, X, s.stream.n_num_features, s.n_categories_per_cat_feature)
end

"""
Each feature is uniform random in range [0, 1]
"""
function next_num_feature(s::RandomTreeGenerator, f_id::Int)
    rand(s.stream.random_state)
end

"""
Return 0 or 1 to build a random one hot encoded element for
each categoric variable.
"""
function next_cat_feature(s::RandomTreeGenerator, f_id::Int)
    cat_id = f_id - s.stream.n_num_features
    first_elem = f_id % s.n_categories_per_cat_feature == 1
    if first_elem
        s.cat_var_selector[] = rand(s.stream.random_state, 1:s.n_categories_per_cat_feature)
    end
    if cat_id % s.n_categories_per_cat_feature == s.cat_var_selector[] - 1
        return 1.0
    end
    return 0.0
end


text:
using Plots
using DifferentialEquations
plotlyjs()

# Rossler space equation
function rossler(t::Float64, x::Vector{Float64}, dx::Vector{Float64})
    dx[1] = -x[2] - x[3]
    dx[2] = x[1] + 0.2*x[2]
    dx[3] = 0.2 + x[3]*(x[1]-5.7)
    return dx
end

# Simulate Rossler
time_span = (0.0, 1000.0)
x0 = [1.0, 1.0, 2.0]
prob = ODEProblem(rossler, x0, time_span)
sol = solve(prob, RK4(), reltol=1e-8)

plot(sol, vars=(1, 2, 3))

max_stars_repo_name: MACSIN-UFMG/julia-rossler-example

max_issues_repo_name: MACSIN-UFMG/julia-rossler-example

max_forks_repo_name: MACSIN-UFMG/julia-rossler-example


text:
# Write function to automate plotting of scale-location plot
# - Plots scale-location plot for assessing linear model fits

function plot_scale_loc(model)

    # Extract predicted values from model object
    pred = GLM.predict(model)

    # Extract residuals from model object
    res = GLM.residuals(model)

    # Calculate mean of residuals
    res_mean = StatsBase.mean(res)

    # Calculate length (n) of residuals vector
    n = length(res)

    # Calculate standardised residuals
    stdRes = (res - res_mean * ones(n)) / (Statistics.std(res))

    # Calculate square-root of standardised residuals
    sqrtStdRes = sqrt.(abs.(stdRes))

    # Create plot
    p = Gadfly.plot(x = pred,
                    y = sqrtStdRes,
                    Gadfly.layer(Gadfly.Geom.point),
                    Gadfly.layer(Gadfly.Geom.smooth(method = :loess,
                                                    smoothing = 0.9)),
                    # Change plot aesthetics
         Gadfly.Guide.title("Scale-Location"),
         Gadfly.Guide.xlabel("Fitted values",
                             orientation=:horizontal),
         Gadfly.Guide.ylabel("(|Standardized residuals|)"))

end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
include("imageio.jl")
import imageio

# Split out mandel as a function so we can call it if need be.
function mandel(xmin=-2.0, xmax=2.0, ymin=-2.0, ymax=2.0, max_iter=256, xres=2048, yres=2048, filename="mandel.pgm")

  rmax = 2.0
  pixels = zeros(Int64, xres, yres)

  for px = 1:xres
    x0 = (Float64(px)/Float64(xres)) * Float64(xmax - xmin) + xmin

    for py = 1:xres
      y0 = (Float64(py)/Float64(yres)) * Float64(ymax - ymin) + ymin

      iter = 0

      x = 0.0
      y = 0.0

      while (((x*x + y*y) < (rmax * rmax)) && (iter < max_iter))
        xt = (x*x) - (y*y) + x0
        y = 2.0 * x * y + y0
        x = xt
        iter = iter + 1
      end

      pixels[px,py] = max_iter - iter
    end

  end

  imageio.writepgm(pixels, max_iter - 1, "mandel.pgm")

end

mandel()

max_stars_repo_name: owainkenwayucl/mandel

max_issues_repo_name: owainkenwayucl/mandel

max_forks_repo_name: owainkenwayucl/mandel


text:
# ---
# title: 495. Teemo Attacking
# id: problem495
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/teemo-attacking/description/>
# hidden: true
# ---
#
# In LOL world, there is a hero called Teemo and his attacking can make his
# enemy Ashe be in poisoned condition. Now, given the Teemo's attacking
# **ascending** time series towards Ashe and the poisoning time duration per
# Teemo's attacking, you need to output the total time that Ashe is in poisoned
# condition.
#
# You may assume that Teemo attacks at the very beginning of a specific time
# point, and makes Ashe be in poisoned condition immediately.
#
# **Example 1:**
#
#
#
#     Input: [1,4], 2
#     Output: 4
#     Explanation: At time point 1, Teemo starts attacking Ashe and makes Ashe be poisoned immediately.
#     This poisoned status will last 2 seconds until the end of time point 2.
#     And at time point 4, Teemo attacks Ashe again, and causes Ashe to be in poisoned status for another 2 seconds.
#     So you finally need to output 4.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: [1,2], 2
#     Output: 3
#     Explanation: At time point 1, Teemo starts attacking Ashe and makes Ashe be poisoned.
#     This poisoned status will last 2 seconds until the end of time point 2.
#     However, at the beginning of time point 2, Teemo attacks Ashe again who is already in poisoned status.
#     Since the poisoned status won't add up together, though the second poisoning attack will still work at time point 2, it will stop at the end of time point 3.
#     So you finally need to output 3.
#
#
#
#
# **Note:**
#
#   1. You may assume the length of given time series array won't exceed 10000.
#   2. You may assume the numbers in the Teemo's attacking time series and his poisoning time duration per attacking are non-negative integers, which won't exceed 10,000,000.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using MultivariatePolynomials
using HomotopyContinuation
using MultivariateMoments
using SemialgebraicSets
using LinearAlgebra
using JuMP

struct TKMPMeasureBuilder
	maxdegree::Int
	on
	vars
end

mutable struct TKMPMeasure
	var_map
	vars
	maxdegree::Int
end

struct NonnegativeOnKConstraintBuilder
	m::TKMPMeasure
	k::AbstractSemialgebraicSet
end

# -----------------------------------------------------------------------------
# utils
# -----------------------------------------------------------------------------

SemialgebraicSets.inequalities(::AbstractSemialgebraicSet) = []

function coeffs(p::AbstractPolynomialLike, vars)
	deg = maxdegree(p)
	ms = monomials(vars, 0:deg) |> reverse
	cs = MultivariatePolynomials.coefficient.(Ref(p), ms, Ref(vars))

	cs, ms
end

function integrate(p::AbstractPolynomialLike, d::TKMPMeasure)
	trimvec(v::AbstractVector, t) = v[begin:min(t, length(v))]

	l = length(d.var_map.data)
	cs, ms = trimvec.(coeffs(p, d.vars), l)
	us = map(m -> d.var_map[m], ms)

	dot(cs, us)
end

integrate(p::Number, d::TKMPMeasure) = integrate(p * only(monomials(d.vars, 0)), d)
integrate(p::AbstractPolynomialLike, d::AbstractVector{TKMPMeasure}) = reduce(integrate, d; init=p)
∫ = integrate

function _moment_matrix(m::TKMPMeasure)
	xd = monomials(m.vars, 0:m.maxdegree ÷ 2) |> reverse
	Md = xd * xd'
	integrate.(Md, Ref(m))
end

function localizing_matrix(m::TKMPMeasure, g::AbstractPolynomialLike)
	round_up_even(x) = cld(x,2)*2

	tg = round_up_even(maxdegree(g))
	xd = monomials(m.vars, 0:(m.maxdegree ÷ 2 - tg ÷ 2)) |> reverse
	Ld = xd * xd' * g
	integrate.(Ld, Ref(m))
end

function JuMP.value(m::TKMPMeasure)
	polyeval(p::AbstractPolynomialLike) = value(p)((variables(p) .=> 1)...)
	b = monomials(m.vars, 0:m.maxdegree ÷ 2) |> reverse

	tol = 1e-3
	res = nothing
	while isnothing(res) && tol <= 1e-1
		solver = SemialgebraicSetsHCSolver(; excess_residual_tol = tol, real_tol = tol, compile = false)
		ma = polyeval.(_moment_matrix(m))

		res = extractatoms(MomentMatrix(ma, b), tol, solver)
		tol *= 1.2
	end
	res
end

# -----------------------------------------------------------------------------
# JuMP Variables
# -----------------------------------------------------------------------------

function JuMP.build_variable(_err::Function, info::VariableInfo, ::Type{TKMPMeasure}; maxdegree, on)
	TKMPMeasureBuilder(maxdegree, on, variables(on))
end

function JuMP.add_variable(model::JuMP.Model, b::TKMPMeasureBuilder, name::String)
	ms = monomials(b.vars, 0:b.maxdegree) |> reverse
	var_map = @variable(model, [ms]; base_name=name)

	μ = TKMPMeasure(var_map, b.vars, b.maxdegree)

	@constraint(model, _moment_matrix(μ) in PSDCone(), base_name=string(name," moment matrix"))
	@constraint(model, [i in inequalities(b.on)], localizing_matrix(μ, i) in PSDCone(), base_name=string(name," localizing matrix inequalities"))
	@constraint(model, [i in equalities(b.on)], localizing_matrix(μ, i) .== 0, base_name=string(name," localizing matrix equalities"))

	μ
end


text:
function string_rotation(str1 :: ASCIIString, str2 :: ASCIIString)
  """
  Check if a string is a sub-string of another string
  """
  function is_substring(str1 :: ASCIIString, str2 :: ASCIIString)
    return contains(str1, str2)
  end

  str :: ASCIIString = string(str1, str1)
  return is_substring(str, str2)
end

# Time  : O(1)
# Space : O(1)
# Note  : assume the function string_rotation takes O(1)

using Base: Test

@test string_rotation("waterbottle", "erbottlewat") == true
@test string_rotation("waterbottle", "erbottlewta") == false

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using HDF5

batch_size = 64

net = Net(batch_size)
data, data_value = MemoryDataLayer(net, :data, (112, 112, 64))
rand!(data_value)

conv1_1 = ConvolutionLayer( :conv1_1, net, data,    128,  3, 1, 1)

init(net)

forward(net; phase=Latte.Test)
backward(net)

forward_task = net.forward_tasks[Latte.Test][end]
forward_args = []
for arg in forward_task.args
    if isa(arg, Symbol)
        push!(forward_args, get_buffer(net,arg))
    else
        push!(forward_args, arg)
    end
end
function forward_bench()
    forward_task.func(forward_args...)
end

backward_task = net.backward_tasks[Latte.Train][end]
backward_args = []
for arg in backward_task.args
    if isa(arg, Symbol)
        push!(backward_args, get_buffer(net,arg))
    else
        push!(backward_args, arg)
    end
end
function backward_bench()
    backward_task.func(backward_args...)
end

for i = 1:3
    forward_bench()
    backward_bench()
end

num_trials = 10

forward_time = 0.0
backward_time = 0.0
for i = 1:num_trials
    tic()
    forward_bench()
    forward_time += toq()
    tic()
    backward_bench()
    backward_time += toq()
end
println("Avg forward time for $num_trials runs: $(forward_time / num_trials * 1000.0)ms")
println("Avg backward time for $num_trials runs: $(backward_time / num_trials * 1000.0)ms")

function convolution_forward(filter::Array, bias::Array, input::Array,
                             output::Array, stride, pad, kernel)
    width, height, channel, num = size(input)
    width_out, height_out, n_filter, _ = size(output)

    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out,
        k = 1:channel, p = 1:kernel, q = 1:kernel

        in_y = (y-1) * stride - pad + p
        in_x = (x-1) * stride - pad + q
        if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
            filter_idx = ((k - 1) * kernel + (p - 1)) * kernel + q
            output[x, y, o, n] += input[in_x, in_y, k, n] *
                                  filter[filter_idx, o]
        end
    end

    # add bias
    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out
        output[x, y, o, n] += bias[o]
    end

    return output
end

input   = data_value
filters = get_buffer(net, :conv1_1weights)
bias    = get_buffer(net, :conv1_1bias)
actual  = get_buffer(net, :conv1_1value)
fill!(actual, 0.0f0)
expected = zeros(actual)
forward_bench()
convolution_forward(filters, bias, input, expected, 1, 1, 3)
ϵ = 1e-4
println("Checking correctness")
if !all(-ϵ .< expected - actual .< ϵ)
    throw("Error: ")
end
println("Done")


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.

"""
Given a flatbuffers metadata type definition (a Field instance from Schema.fbs),
translate to the appropriate Julia storage eltype
"""
function juliaeltype end

finaljuliatype(T) = T
finaljuliatype(::Type{Missing}) = Missing
finaljuliatype(::Type{Union{T, Missing}}) where {T} = Union{Missing, finaljuliatype(T)}

"""
Given a FlatBuffers.Builder and a Julia column or column eltype,
Write the field.type flatbuffer definition of the eltype
"""
function arrowtype end

arrowtype(b, col::AbstractVector{T}) where {T} = arrowtype(b, maybemissing(T))
arrowtype(b, col::DictEncoded) = arrowtype(b, col.encoding.data)
arrowtype(b, col::Compressed) = arrowtype(b, col.data)

function juliaeltype(f::Meta.Field, ::Nothing, convert::Bool)
    T = juliaeltype(f, convert)
    return convert ? finaljuliatype(T) : T
end

function juliaeltype(f::Meta.Field, meta::Dict{String, String}, convert::Bool)
    TT = juliaeltype(f, convert)
    !convert && return TT
    T = finaljuliatype(TT)
    # deprecated begin
    TTT = ArrowTypes.extensiontype(f, meta)
    TTT !== nothing && return TTT
    # end deprecated
    if haskey(meta, "ARROW:extension:name")
        typename = meta["ARROW:extension:name"]
        metadata = get(meta, "ARROW:extension:metadata", "")
        JT = ArrowTypes.JuliaType(Val(Symbol(typename)), maybemissing(TT), metadata)
        if JT !== nothing
            return f.nullable ? Union{JT, Missing} : JT
        else
            @warn "unsupported ARROW:extension:name type: \"$typename\", arrow type = $TT"
        end
    end
    return something(TTT, T)
end

function juliaeltype(f::Meta.Field, convert::Bool)
    T = juliaeltype(f, f.type, convert)
    return f.nullable ? Union{T, Missing} : T
end

juliaeltype(f::Meta.Field, ::Meta.Null, convert) = Missing

function arrowtype(b, ::Type{Missing})
    Meta.nullStart(b)
    return Meta.Null, Meta.nullEnd(b), nothing
end

function juliaeltype(f::Meta.Field, int::Meta.Int, convert)
    if int.is_signed
        if int.bitWidth == 8
            Int8
        elseif int.bitWidth == 16
            Int16
        elseif int.bitWidth == 32
            Int32
        elseif int.bitWidth == 64
            Int64
        elseif int.bitWidth == 128
            Int128
        else
            error("$int is not valid arrow type metadata")
        end
    else
        if int.bitWidth == 8
            UInt8
        elseif int.bitWidth == 16
            UInt16
        elseif int.bitWidth == 32
            UInt32
        elseif int.bitWidth == 64
            UInt64
        elseif int.bitWidth == 128
            UInt128
        else
            error("$int is not valid arrow type metadata")
        end
    end
end

function arrowtype(b, ::Type{T}) where {T <: Integer}
    Meta.intStart(b)
    Meta.intAddBitWidth(b, Int32(8 * sizeof(T)))
    Meta.intAddIsSigned(b, T <: Signed)
    return Meta.Int, Meta.intEnd(b), nothing
end

# primitive types
function juliaeltype(f::Meta.Field, fp::Meta.FloatingPoint, convert)
    if fp.precision == Meta.Precision.HALF
        Float16
    elseif fp.precision == Meta.Precision.SINGLE
        Float32
    elseif fp.precision == Meta.Precision.DOUBLE
        Float64
    end
end

function arrowtype(b, ::Type{T}) where {T <: AbstractFloat}
    Meta.floatingPointStart(b)
    Meta.floatingPointAddPrecision(b, T === Float16 ? Meta.Precision.HALF : T === Float32 ? Meta.Precision.SINGLE : Meta.Precision.DOUBLE)
    return Meta.FloatingPoint, Meta.floatingPointEnd(b), nothing
end

juliaeltype(f::Meta.Field, b::Union{Meta.Utf8, Meta.LargeUtf8}, convert) = String

datasizeof(x) = sizeof(x)
datasizeof(x::AbstractVector) = sum(datasizeof, x)

juliaeltype(f::Meta.Field, b::Union{Meta.Binary, Meta.LargeBinary}, convert) = Vector{UInt8}

juliaeltype(f::Meta.Field, x::Meta.FixedSizeBinary, convert) = NTuple{Int(x.byteWidth), UInt8}

# arggh!
Base.write(io::IO, x::NTuple{N, T}) where {N, T} = sum(y -> Base.write(io, y), x)

juliaeltype(f::Meta.Field, x::Meta.Bool, convert) = Bool

function arrowtype(b, ::Type{Bool})
    Meta.boolStart(b)
    return Meta.Bool, Meta.boolEnd(b), nothing
end

struct Decimal{P, S, T}
    value::T # only Int128 or Int256
end

Base.zero(::Type{Decimal{P, S, T}}) where {P, S, T} = Decimal{P, S, T}(T(0))
==(a::Decimal{P, S, T}, b::Decimal{P, S, T}) where {P, S, T} = ==(a.value, b.value)
Base.isequal(a::Decimal{P, S, T}, b::Decimal{P, S, T}) where {P, S, T} = isequal(a.value, b.value)

function juliaeltype(f::Meta.Field, x::Meta.Decimal, convert)
    return Decimal{x.precision, x.scale, x.bitWidth == 256 ? Int256 : Int128}
end

ArrowTypes.ArrowKind(::Type{<:Decimal}) = PrimitiveKind()

function arrowtype(b, ::Type{Decimal{P, S, T}}) where {P, S, T}
    Meta.decimalStart(b)
    Meta.decimalAddPrecision(b, Int32(P))
    Meta.decimalAddScale(b, Int32(S))
    Meta.decimalAddBitWidth(b, Int32(T == Int256 ? 256 : 128))
    return Meta.Decimal, Meta.decimalEnd(b), nothing
end

Base.write(io::IO, x::Decimal) = Base.write(io, x.value)

abstract type ArrowTimeType end
Base.write(io::IO, x::ArrowTimeType) = Base.write(io, x.x)
ArrowTypes.ArrowKind(::Type{<:ArrowTimeType}) = PrimitiveKind()

struct Date{U, T} <: ArrowTimeType
    x::T
end

const DATE = Date{Meta.DateUnit.DAY, Int32}
Base.zero(::Type{Date{U, T}}) where {U, T} = Date{U, T}(T(0))
storagetype(::Type{Date{U, T}}) where {U, T} = T
bitwidth(x::Meta.DateUnit) = x == Meta.DateUnit.DAY ? Int32 : Int64
Date{Meta.DateUnit.DAY}(days) = DATE(Int32(days))
Date{Meta.DateUnit.MILLISECOND}(ms) = Date{Meta.DateUnit.MILLISECOND, Int64}(Int64(ms))

juliaeltype(f::Meta.Field, x::Meta.Date, convert) = Date{x.unit, bitwidth(x.unit)}
finaljuliatype(::Type{DATE}) = Dates.Date
Base.convert(::Type{Dates.Date}, x::DATE) = Dates.Date(Dates.UTD(Int64(x.x + UNIX_EPOCH_DATE)))
finaljuliatype(::Type{Date{Meta.DateUnit.MILLISECOND, Int64}}) = Dates.DateTime
Base.convert(::Type{Dates.DateTime}, x::Date{Meta.DateUnit.MILLISECOND, Int64}) = Dates.DateTime(Dates.UTM(Int64(x.x + UNIX_EPOCH_DATETIME)))

function arrowtype(b, ::Type{Date{U, T}}) where {U, T}
    Meta.dateStart(b)
    Meta.dateAddUnit(b, U)
    return Meta.Date, Meta.dateEnd(b), nothing
end

const UNIX_EPOCH_DATE = Dates.value(Dates.Date(1970))
Base.convert(::Type{DATE}, x::Dates.Date) = DATE(Int32(Dates.value(x) - UNIX_EPOCH_DATE))

const UNIX_EPOCH_DATETIME = Dates.value(Dates.DateTime(1970))
Base.convert(::Type{Date{Meta.DateUnit.MILLISECOND, Int64}}, x::Dates.DateTime) = Date{Meta.DateUnit.MILLISECOND, Int64}(Int64(Dates.value(x) - UNIX_EPOCH_DATETIME))

ArrowTypes.ArrowType(::Type{Dates.Date}) = DATE
ArrowTypes.toarrow(x::Dates.Date) = convert(DATE, x)
const DATE_SYMBOL = Symbol("JuliaLang.Date")
ArrowTypes.arrowname(::Type{Dates.Date}) = DATE_SYMBOL
ArrowTypes.JuliaType(::Val{DATE_SYMBOL}, S) = Dates.Date
ArrowTypes.fromarrow(::Type{Dates.Date}, x::DATE) = convert(Dates.Date, x)
ArrowTypes.default(::Type{Dates.Date}) = Dates.Date(1,1,1)

struct Time{U, T} <: ArrowTimeType
    x::T
end

Base.zero(::Type{Time{U, T}}) where {U, T} = Time{U, T}(T(0))
const TIME = Time{Meta.TimeUnit.NANOSECOND, Int64}

bitwidth(x::Meta.TimeUnit) = x == Meta.TimeUnit.SECOND || x == Meta.TimeUnit.MILLISECOND ? Int32 : Int64
Time{U}(x) where {U <: Meta.TimeUnit} = Time{U, bitwidth(U)}(bitwidth(U)(x))
storagetype(::Type{Time{U, T}}) where {U, T} = T
juliaeltype(f::Meta.Field, x::Meta.Time, convert) = Time{x.unit, bitwidth(x.unit)}
finaljuliatype(::Type{<:Time}) = Dates.Time
periodtype(U::Meta.TimeUnit) = U === Meta.TimeUnit.SECOND ? Dates.Second :
                               U === Meta.TimeUnit.MILLISECOND ? Dates.Millisecond :
                               U === Meta.TimeUnit.MICROSECOND ? Dates.Microsecond : Dates.Nanosecond
Base.convert(::Type{Dates.Time}, x::Time{U, T}) where {U, T} = Dates.Time(Dates.Nanosecond(Dates.tons(periodtype(U)(x.x))))

function arrowtype(b, ::Type{Time{U, T}}) where {U, T}
    Meta.timeStart(b)
    Meta.timeAddUnit(b, U)
    Meta.timeAddBitWidth(b, Int32(8 * sizeof(T)))
    return Meta.Time, Meta.timeEnd(b), nothing
end

Base.convert(::Type{TIME}, x::Dates.Time) = TIME(Dates.value(x))

ArrowTypes.ArrowType(::Type{Dates.Time}) = TIME
ArrowTypes.toarrow(x::Dates.Time) = convert(TIME, x)
const TIME_SYMBOL = Symbol("JuliaLang.Time")
ArrowTypes.arrowname(::Type{Dates.Time}) = TIME_SYMBOL
ArrowTypes.JuliaType(::Val{TIME_SYMBOL}, S) = Dates.Time
ArrowTypes.fromarrow(::Type{Dates.Time}, x::TIME) = convert(Dates.Time, x)
ArrowTypes.default(::Type{Dates.Time}) = Dates.Time(1,1,1)

struct Timestamp{U, TZ} <: ArrowTimeType
    x::Int64
end

Base.zero(::Type{Timestamp{U, T}}) where {U, T} = Timestamp{U, T}(Int64(0))

function juliaeltype(f::Meta.Field, x::Meta.Timestamp, convert)
    return Timestamp{x.unit, x.timezone === nothing ? nothing : Symbol(x.timezone)}
end

const DATETIME = Timestamp{Meta.TimeUnit.MILLISECOND, nothing}

finaljuliatype(::Type{Timestamp{U, TZ}}) where {U, TZ} = ZonedDateTime
finaljuliatype(::Type{Timestamp{U, nothing}}) where {U} = DateTime
Base.convert(::Type{ZonedDateTime}, x::Timestamp{U, TZ}) where {U, TZ} =
    ZonedDateTime(Dates.DateTime(Dates.UTM(Int64(Dates.toms(periodtype(U)(x.x)) + UNIX_EPOCH_DATETIME))), TimeZone(String(TZ)))
Base.convert(::Type{DateTime}, x::Timestamp{U, nothing}) where {U} =
    Dates.DateTime(Dates.UTM(Int64(Dates.toms(periodtype(U)(x.x)) + UNIX_EPOCH_DATETIME)))
Base.convert(::Type{Timestamp{Meta.TimeUnit.MILLISECOND, TZ}}, x::ZonedDateTime) where {TZ} =
    Timestamp{Meta.TimeUnit.MILLISECOND, TZ}(Int64(Dates.value(DateTime(x)) - UNIX_EPOCH_DATETIME))
Base.convert(::Type{Timestamp{Meta.TimeUnit.MILLISECOND, nothing}}, x::DateTime) =
    Timestamp{Meta.TimeUnit.MILLISECOND, nothing}(Int64(Dates.value(x) - UNIX_EPOCH_DATETIME))

function arrowtype(b, ::Type{Timestamp{U, TZ}}) where {U, TZ}
    tz = TZ !== nothing ? FlatBuffers.createstring!(b, String(TZ)) : FlatBuffers.UOffsetT(0)
    Meta.timestampStart(b)
    Meta.timestampAddUnit(b, U)
    Meta.timestampAddTimezone(b, tz)
    return Meta.Timestamp, Meta.timestampEnd(b), nothing
end

ArrowTypes.ArrowType(::Type{Dates.DateTime}) = DATETIME
ArrowTypes.toarrow(x::Dates.DateTime) = convert(DATETIME, x)
const DATETIME_SYMBOL = Symbol("JuliaLang.DateTime")
ArrowTypes.arrowname(::Type{Dates.DateTime}) = DATETIME_SYMBOL
ArrowTypes.JuliaType(::Val{DATETIME_SYMBOL}, S) = Dates.DateTime
ArrowTypes.fromarrow(::Type{Dates.DateTime}, x::DATETIME) = convert(Dates.DateTime, x)
ArrowTypes.default(::Type{Dates.DateTime}) = Dates.DateTime(1,1,1,1,1,1)

ArrowTypes.ArrowType(::Type{ZonedDateTime}) = Timestamp
ArrowTypes.toarrow(x::ZonedDateTime) = convert(Timestamp{Meta.TimeUnit.MILLISECOND, Symbol(x.timezone)}, x)
const ZONEDDATETIME_SYMBOL = Symbol("JuliaLang.ZonedDateTime")
ArrowTypes.arrowname(::Type{ZonedDateTime}) = ZONEDDATETIME_SYMBOL
ArrowTypes.JuliaType(::Val{ZONEDDATETIME_SYMBOL}, S) = ZonedDateTime
ArrowTypes.fromarrow(::Type{ZonedDateTime}, x::Timestamp) = convert(ZonedDateTime, x)
ArrowTypes.default(::Type{TimeZones.ZonedDateTime}) = TimeZones.ZonedDateTime(1,1,1,1,1,1,TimeZones.tz"UTC")

struct Interval{U, T} <: ArrowTimeType
    x::T
end

Base.zero(::Type{Interval{U, T}}) where {U, T} = Interval{U, T}(T(0))

bitwidth(x::Meta.IntervalUnit) = x == Meta.IntervalUnit.YEAR_MONTH ? Int32 : Int64
Interval{Meta.IntervalUnit.YEAR_MONTH}(x) = Interval{Meta.IntervalUnit.YEAR_MONTH, Int32}(Int32(x))
Interval{Meta.IntervalUnit.DAY_TIME}(x) = Interval{Meta.IntervalUnit.DAY_TIME, Int64}(Int64(x))

function juliaeltype(f::Meta.Field, x::Meta.Interval, convert)
    return Interval{x.unit, bitwidth(x.unit)}
end

function arrowtype(b, ::Type{Interval{U, T}}) where {U, T}
    Meta.intervalStart(b)
    Meta.intervalAddUnit(b, U)
    return Meta.Interval, Meta.intervalEnd(b), nothing
end

struct Duration{U} <: ArrowTimeType
    x::Int64
end

Base.zero(::Type{Duration{U}}) where {U} = Duration{U}(Int64(0))

function juliaeltype(f::Meta.Field, x::Meta.Duration, convert)
    return Duration{x.unit}
end

finaljuliatype(::Type{Duration{U}}) where {U} = periodtype(U)
Base.convert(::Type{P}, x::Duration{U}) where {P <: Dates.Period, U} = P(periodtype(U)(x.x))

function arrowtype(b, ::Type{Duration{U}}) where {U}
    Meta.durationStart(b)
    Meta.durationAddUnit(b, U)
    return Meta.Duration, Meta.durationEnd(b), nothing
end

arrowperiodtype(P) = Meta.TimeUnit.SECOND
arrowperiodtype(::Type{Dates.Millisecond}) = Meta.TimeUnit.MILLISECOND
arrowperiodtype(::Type{Dates.Microsecond}) = Meta.TimeUnit.MICROSECOND
arrowperiodtype(::Type{Dates.Nanosecond}) = Meta.TimeUnit.NANOSECOND

Base.convert(::Type{Duration{U}}, x::Dates.Period) where {U} = Duration{U}(Dates.value(periodtype(U)(x)))

ArrowTypes.ArrowType(::Type{P}) where {P <: Dates.Period} = Duration{arrowperiodtype(P)}
ArrowTypes.toarrow(x::P) where {P <: Dates.Period} = convert(Duration{arrowperiodtype(P)}, x)
const PERIOD_SYMBOL = Symbol("JuliaLang.Dates.Period")
ArrowTypes.arrowname(::Type{P}) where {P <: Dates.Period} = PERIOD_SYMBOL
ArrowTypes.JuliaType(::Val{PERIOD_SYMBOL}, ::Type{Duration{U}}) where {U} = periodtype(U)
ArrowTypes.fromarrow(::Type{P}, x::Duration{U}) where {P <: Dates.Period, U} = convert(P, x)

# nested types; call juliaeltype recursively on nested children
function juliaeltype(f::Meta.Field, list::Union{Meta.List, Meta.LargeList}, convert)
    return Vector{juliaeltype(f.children[1], buildmetadata(f.children[1]), convert)}
end

# arrowtype will call fieldoffset recursively for children
function arrowtype(b, x::List{T, O, A}) where {T, O, A}
    if eltype(A) == UInt8
        if T <: AbstractString || T <: Union{AbstractString, Missing}
            if O == Int32
                Meta.utf8Start(b)
                return Meta.Utf8, Meta.utf8End(b), nothing
            else # if O == Int64
                Meta.largUtf8Start(b)
                return Meta.LargeUtf8, Meta.largUtf8End(b), nothing
            end
        else # if Vector{UInt8}
            if O == Int32
                Meta.binaryStart(b)
                return Meta.Binary, Meta.binaryEnd(b), nothing
            else # if O == Int64
                Meta.largeBinaryStart(b)
                return Meta.LargeBinary, Meta.largeBinaryEnd(b), nothing
            end
        end
    else
        children = [fieldoffset(b, "", x.data)]
        if O == Int32
            Meta.listStart(b)
            return Meta.List, Meta.listEnd(b), children
        else
            Meta.largeListStart(b)
            return Meta.LargeList, Meta.largeListEnd(b), children
        end
    end
end

function juliaeltype(f::Meta.Field, list::Meta.FixedSizeList, convert)
    type = juliaeltype(f.children[1], buildmetadata(f.children[1]), convert)
    return NTuple{Int(list.listSize), type}
end

function arrowtype(b, x::FixedSizeList{T, A}) where {T, A}
    N = ArrowTypes.getsize(ArrowTypes.ArrowKind(ArrowTypes.ArrowType(Base.nonmissingtype(T))))
    if eltype(A) == UInt8
        Meta.fixedSizeBinaryStart(b)
        Meta.fixedSizeBinaryAddByteWidth(b, Int32(N))
        return Meta.FixedSizeBinary, Meta.fixedSizeBinaryEnd(b), nothing
    else
        children = [fieldoffset(b, "", x.data)]
        Meta.fixedSizeListStart(b)
        Meta.fixedSizeListAddListSize(b, Int32(N))
        return Meta.FixedSizeList, Meta.fixedSizeListEnd(b), children
    end
end

function juliaeltype(f::Meta.Field, map::Meta.Map, convert)
    K = juliaeltype(f.children[1].children[1], buildmetadata(f.children[1].children[1]), convert)
    V = juliaeltype(f.children[1].children[2], buildmetadata(f.children[1].children[2]), convert)
    return Dict{K, V}
end

function arrowtype(b, x::Map)
    children = [fieldoffset(b, "entries", x.data)]
    Meta.mapStart(b)
    return Meta.Map, Meta.mapEnd(b), children
end

struct KeyValue{K, V}
    key::K
    value::V
end
keyvalueK(::Type{KeyValue{K, V}}) where {K, V} = K
keyvalueV(::Type{KeyValue{K, V}}) where {K, V} = V
Base.length(kv::KeyValue) = 1
Base.iterate(kv::KeyValue, st=1) = st === nothing ? nothing : (kv, nothing)
ArrowTypes.default(::Type{KeyValue{K, V}}) where {K, V} = KeyValue(default(K), default(V))

function arrowtype(b, ::Type{KeyValue{K, V}}) where {K, V}
    children = [fieldoffset(b, "key", K), fieldoffset(b, "value", V)]
    Meta.structStart(b)
    return Meta.Struct, Meta.structEnd(b), children
end

function juliaeltype(f::Meta.Field, list::Meta.Struct, convert)
    names = Tuple(Symbol(x.name) for x in f.children)
    types = Tuple(juliaeltype(x, buildmetadata(x), convert) for x in f.children)
    return NamedTuple{names, Tuple{types...}}
end

function arrowtype(b, x::Struct{T, S}) where {T, S}
    names = fieldnames(Base.nonmissingtype(T))
    children = [fieldoffset(b, names[i], x.data[i]) for i = 1:length(names)]
    Meta.structStart(b)
    return Meta.Struct, Meta.structEnd(b), children
end

# Unions
function UnionT(f::Meta.Field, convert)
    typeids = f.type.typeIds === nothing ? nothing : Tuple(Int(x) for x in f.type.typeIds)
    UT =  UnionT{f.type.mode, typeids, Tuple{(juliaeltype(x, buildmetadata(x), convert) for x in f.children)...}}
    return UT
end

juliaeltype(f::Meta.Field, u::Meta.Union, convert) = Union{(juliaeltype(x, buildmetadata(x), convert) for x in f.children)...}

function arrowtype(b, x::Union{DenseUnion{S, UnionT{T, typeIds, U}}, SparseUnion{S, UnionT{T, typeIds, U}}}) where {S, T, typeIds, U}
    if typeIds !== nothing
        Meta.unionStartTypeIdsVector(b, length(typeIds))
        for id in Iterators.reverse(typeIds)
            FlatBuffers.prepend!(b, id)
        end
        TI = FlatBuffers.endvector!(b, length(typeIds))
    end
    children = [fieldoffset(b, "", x.data[i]) for i = 1:fieldcount(U)]
    Meta.unionStart(b)
    Meta.unionAddMode(b, T)
    if typeIds !== nothing
        Meta.unionAddTypeIds(b, TI)
    end
    return Meta.Union, Meta.unionEnd(b), children
end


text:
projected_position_ECI2ECEF = PVTSolution(
    ECEF([4.018703954398344e6, 427061.6186769246, 4.918310675640356e6]), -2.1332062204026792e7, 2.346849505593898)

@testset "User position ECI2ECEF" begin
    positions = map(x -> sat_position_ECI_2_ECEF(x), satellite_states)
    pseudo_ranges = PositionVelocityTime.pseudo_ranges(satellite_states)
    out = user_position(positions, pseudo_ranges)
    @test out.pos ≈ projected_position_ECI2ECEF.pos
    @test out.receiver_time_correction ≈ projected_position_ECI2ECEF.receiver_time_correction
    @test out.GDOP ≈ projected_position_ECI2ECEF.GDOP

    out = false
    try
        calc_PVT([sv1_struct])
    catch e
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true
end

projected_position_ECEF = PVTSolution(
    ECEF([4.0187039543983485e6, 427061.6186769258, 4.9183106756403595e6]),
    -2.133206220402679e7,
    2.346849505593898)
@testset "User position ECEF" begin
    positions = map( x -> sat_position_ECEF(x), satellite_states)
    pseudo_ranges = PositionVelocityTime.pseudo_ranges(satellite_states)
    out = user_position(positions, pseudo_ranges)
    @test out.pos ≈ projected_position_ECEF.pos
    @test out.receiver_time_correction ≈ projected_position_ECEF.receiver_time_correction
    @test out.GDOP ≈ projected_position_ECEF.GDOP

    out = false
    try
        calc_PVT([sv1_struct])
    catch e
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true
end


text:
using Documenter, PseudoArcLengthContinuation

makedocs(doctest = false,
	sitename = "Pseudo Arc Length Continuation in Julia",
	format = Documenter.HTML(collapselevel = 1),
	# format = DocumenterLaTeX.LaTeX(),
	authors = "Romain Veltz",
	pages = Any[
		"Home" => "index.md",
		"Tutorials" => [
			"1/ Temperature model" => "tutorials1.md",
			"2/ Temperature model with Spectral Collocation" => "tutorials1b.md",
			"3/ Swift-Hohenberg" => "tutorials2.md",
			"4/ The Swift-Hohenberg on the GPU (non-local)" => "tutorials2b.md",
			"5/ Brusselator 1d" => "tutorials3.md",
			"6/ Period Doubling in BVAM model" => "tutorialsPD.md",
			"7/ Ginzburg-Landau 2d (GPU)" => "tutorialsCGL.md",
			"8/ Ginzburg-Landau 2d (Shooting)" => "tutorialsCGLShoot.md",
		],
		"Advanced Usage" => [
			"Bordered arrays" => "Borderedarrays.md",
			"Linear / Eigen Solvers" => "linearsolver.md",
			"Bordered linear solvers" => "borderedlinearsolver.md",
			"Bifurcations" => "detectionBifurcation.md",
			"Fold / Hopf Continuation" => "codim2Continuation.md",
			"Constrained problem" => "constrainedproblem.md",
			"Periodic Orbits (Finite Differences)" => "periodicOrbitFD.md",
			"Periodic Orbits (Shooting)" => "periodicOrbitShooting.md",
			"DiffEq wrapper" => "diffeq.md",
			"Iterator Interface" => "iterator.md",
		],
		"Frequently Asked Questions" => "faq.md",
		"Library" => "library.md"
	]
	)

deploydocs(
	repo = "github.com/rveltz/PseudoArcLengthContinuation.jl.git",
)


text:
using AbInitioSoftwareBase.Inputs: FormatConfig

@testset "Test `FormatConfig` is mutable" begin
    config = FormatConfig(;
        delimiter = " ",
        newline = "\n",
        indent = ' '^4,
        float = "%14.9f",
        int = "%5d",
        bool = ".%.",
    )
    config.newline = "\r\n"
    config.indent = ' '^5
    config.float = "%15.8f"
    @test config == FormatConfig(;
        delimiter = " ",
        newline = "\r\n",
        indent = ' '^5,
        float = "%15.8f",
        int = "%5d",
        bool = ".%.",
    )
    @test_throws MethodError config.delimiter = ' '
end


text:
#########################################################
#=
Generate figures for the 1-d nonstationary phase example.
Corresponds to Figures 2 and 3 in the paper.
Run this script with the command `include("make_figure_2and3.jl")` entered into the Julia REPL in this directory.
This script saves figure2.pdf and figure3.pdf to the current directory.
To include  figure2.pdf and figure3.pdf in the paper it must be copied to the directory paper/
=#
#########################################################

# --- load modules
using NonstationaryPhase
using PyPlot

# --- set the seed
#seedstart = rand(UInt64)
seedstart = 0xe8bac1438a62066d
srand(seedstart)

# --- save figures to disk or not
savefigures = true

# --- set grid geometry
const d_const = 1
const period_const = 10
const nside_const  = nextprod([2,3,5,7], 10_000)

# --- parameters of Z(x)
ν      = 2.0
tild_ν = 2.1
ρ      = 0.005*period_const
tild_ρ = 0.005*period_const
t_0    = 1.5

# --- noise parameters
σpixl    = 0.0
beamFWHM = 0.0

# --- prior parameters
νΦ = 5.0
ρΦ = 0.15*period_const
σΦ = 0.15*(period_const/2π)^2

#########################################################
#=
Generate the first set of simulations
=#
#########################################################

# --- define an instance of PhaseParms which carries all the parameters
parms = PhaseParms(d_const, period_const, nside_const)
ηkfun, Ckfun, C1kfun, C2kfun = gen_tangentMatern(ν=ν, tild_ν=tild_ν, ρ=ρ, tild_ρ=tild_ρ, t_0=t_0)
NonstationaryPhase.CNNkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})   = CNNkfun(p.k, p.deltx; σpixl=σpixl, beamFWHM=beamFWHM)
Cϕϕkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})  = maternk(p.k; ν=νΦ, ρ=ρΦ, σ=σΦ)
parms.ηk[:]       =  ηkfun(parms.k)
parms.Ck[:]       =  Ckfun(parms.k)
parms.C1k[:]      =  C1kfun(parms.k)
parms.C2k[:]      =  C2kfun(parms.k)
parms.CNNk[:]     =  CNNkfun(parms)
parms.Cϕϕk[:]     = Cϕϕkfun(parms)
parms.CZZmk[:]    =  CZZmkfun(parms.Cϕϕk, parms.Ck, parms.ηk, parms.ξk, parms.x, parms.k, parms.deltx, parms.deltk)
parms.CZZmobsk[:] = parms.CZZmk + parms.CNNk
parms.CZZmobsk[parms.r .< 1]              = Inf
parms.CZZmobsk[parms.r .>= 0.9*parms.nyq] = Inf

# Now a High resolution version of PhaseParms for simulation
parmsHR = PhaseParms(d_const, period_const, nside_const)
parmsHR.ηk[:]   =  ηkfun(parmsHR.k)
parmsHR.Ck[:]   =  Ckfun(parmsHR.k)
parmsHR.CNNk[:] =  CNNkfun(parmsHR)
parmsHR.Cϕϕk[:] =  Cϕϕkfun(parmsHR)

# ---- generate the initial simulation and other plotting quantities
zkobs, zk, zx, zx_noϕ, ϕk, ϕx = simNPhaseGRF(parms, parmsHR)
Dθx          = NonstationaryPhase.ϕx_2_Dθx(ϕx, parms)
@time Cℓvar        = Cℓvarfun(parms)
@time Cℓbias       = Cℓbiasfun(parms)
@time approxCℓbias = approxCℓbiasfun(parms)
estϕk        = estϕkfun(zkobs, parms)
# plot(zx)
# plot(NonstationaryPhase.ϕk_2_Dθx(estϕk.* (0 .< parms.r .< 20*parms.deltk), parms),"k")
# plot(NonstationaryPhase.ϕk_2_Dθx(ϕk, parms),"b")

# --- get multiple quadratic estimates based on the same ϕx and the average estimate
function get_estimates{dm,T1,T2}(sims::Integer, ϕx, parms::PhaseParms{dm,T1,T2}, parmsHR::PhaseParms{dm,T1,T2})
	Eestϕk        = zeros(Complex{Float64}, size(parms.k[1],1))
	estimates_Dθx = zeros(Complex{Float64}, size(parms.k[1],1),sims)
	avetime = 0.0
	for loop = 1:sims
		_zkobs, = simNPhaseGRF(ϕx, parms, parmsHR)
		tic()
		_estϕk  = estϕkfun(_zkobs, parms)
		avetime += toc() / sims
		_estϕx  = ifftdr(_estϕk .* (0 .< parms.r .< 20*parms.deltk), parms.deltk)
		estimates_Dθx[:,loop] = NonstationaryPhase.ϕx_2_Dθx(_estϕx, parms)
		Eestϕk += _estϕk .* (1/sims) .* (0 .< parms.r .< 40*parms.deltk)
	end
	EestDθx = NonstationaryPhase.ϕx_2_Dθx(ifftdr(Eestϕk,parms.deltk), parms)
	return EestDθx, Eestϕk, estimates_Dθx, avetime
end
EestDθx, Eestϕk, estimates_Dθx, avetime1 = get_estimates(100, ϕx, parms, parmsHR)
println("ave quadratic estimate clock speed =  $avetime1")

#########################################################
#=
Generate plots for the first set of simulations
=#
#########################################################

# ------ figure ---- show the quadratic estimate
function makeplots()
	global fig1 = figure(figsize = (10,20))

	# ----- figure: plot Z(x)
	ax1 = subplot2grid((3,2), (0,0), colspan=2)
	plot(fftshift(parms.x[1]), fftshift(zx), "k", linewidth = 1.5)
	axis("tight")
	title("Nonstationary phase simulation")

	# ----- figure: plot the estimates
	ax2 = subplot2grid((3,2), (1,0), colspan=2)
	plot(fftshift(parms.x[1]), fftshift(estimates_Dθx[:,1]),"k", label = L"\hat\theta^\prime(x)", linewidth = 1.5, alpha = 0.4)
	for cntr = 2:min(size(estimates_Dθx,2),5)
		plot(fftshift(parms.x[1]), fftshift(estimates_Dθx[:,cntr]),"k", linewidth = 1.5, alpha = 0.4)
	end
	plot(fftshift(parms.x[1]), fftshift(EestDθx), "k--", label = L"average $\hat\theta^\prime(x)$", linewidth = 3)
	plot(fftshift(parms.x[1]), fftshift(Dθx), "b", label = L"\theta^\prime(x)", linewidth = 2.5,alpha = 0.7)
	axis("tight")
	title(L"Quadratic estimate, average quadratic estimate and the simulation truth $\theta^\prime(x)$")
	legend(loc="best",fontsize = 10)

	# ----- figure ---- Cℓbias and the signal
	ax4 = subplot2grid((3,2), (2,1))
	δ0 = 1 / parms.deltk ^ d_const
	kbins, bias_radpower        = radial_power(abs2(parms.r.*(Eestϕk - ϕk))./δ0, 1, parms)
	kbins, Cbias_radpower       = radial_power(abs2(parms.r).*Cℓbias, 1, parms)
	kbins, approxCbias_radpower = radial_power(abs2(parms.r).*approxCℓbias, 1, parms)
	kbins, CΦΦ_radpower         = radial_power(abs2(parms.r).*parms.Cϕϕk, 1, parms)
	semilogy(kbins[1:25], Cbias_radpower[1:25], "r-", label = L"\ell^2 C^{bias\, \hat\phi}_\ell", linewidth = 2.0)
	semilogy(kbins[1:25], approxCbias_radpower[1:25], "b--", label = L"approx $\ell^2 C^{bias\, \hat\phi}_\ell$", linewidth = 2.0)
	semilogy(kbins[1:25], CΦΦ_radpower[1:25], "k:", label = L"\ell^2C^{\phi\phi}_\ell", linewidth = 3.0)
	semilogy(kbins[1:25], bias_radpower[1:25], "ro", label = L"\ell^2|\phi_\ell - E(\hat\phi_\ell)|^2\delta_0", linewidth = 3.0)
	xlabel(L"frequency $\ell$")
	title("Analytic vrs empirical bias spectra")
	legend(loc="best",fontsize = 10)

	# ----- figure ---- Cℓvar and the signal
	ax3 = subplot2grid((3,2), (2,0), sharex=ax4, sharey=ax4)
	δ0 = 1 / parms.deltk ^ d_const
	kbins, var_radpower  = radial_power(abs2(parms.r.*(Eestϕk - estϕk))./δ0, 1, parms)
	kbins, Cvar_radpower = radial_power(abs2(parms.r).*Cℓvar, 1, parms)
	kbins, CΦΦ_radpower  = radial_power(abs2(parms.r).*parms.Cϕϕk, 1, parms)
	semilogy(kbins[1:25], Cvar_radpower[1:25], "g-", label = L"\ell^2 C^{var\, \hat\phi}_\ell", linewidth = 2.0)
	semilogy(kbins[1:25], CΦΦ_radpower[1:25], "k:", label = L"\ell^2C^{\phi\phi}_\ell", linewidth = 3.0)
	semilogy(kbins[1:25], var_radpower[1:25], "go", label = L"\ell^2|E(\hat\phi_\ell) - \hat\phi_\ell|^2\delta_0", linewidth = 3.0)
	ylabel("spectral density")
	xlabel(L"frequency $\ell$")
	title("Analytic vrs empirical variance spectra")
	legend(loc="best",fontsize = 10)
	axis("tight")

	#ax4[:set_yticklabels]([])
	fig1[:subplots_adjust](wspace=0.2)

end

# --- now evaluate the above `script'
makeplots()

# --- save the figures to disk
if savefigures
    fig1[:savefig]("figure2.pdf", dpi=300, bbox_inches="tight", transparent=true)
    plt[:close](fig1)
end

##############################################################
#=
Compute the minimum cut loci and increase the signal size to bring out the bias
=#
##############################################################

# --- reset the seed to reproduce the same results as above
srand(seedstart)

# --- now decrease t_0 to expose the bias
reduce_t_0_factor = 1/7

# --- define an instance of PhaseParms which carries all the parameters
parms = PhaseParms(d_const, period_const, nside_const)
ηkfun, Ckfun, C1kfun, C2kfun = gen_tangentMatern(ν=ν, tild_ν=tild_ν, ρ=ρ, tild_ρ=tild_ρ, t_0=t_0*reduce_t_0_factor)
NonstationaryPhase.CNNkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})   = CNNkfun(p.k, p.deltx; σpixl=σpixl, beamFWHM=beamFWHM)
Cϕϕkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})  = maternk(p.k; ν=νΦ, ρ=ρΦ, σ=σΦ)
parms.ηk[:]       =  ηkfun(parms.k)
parms.Ck[:]       =  Ckfun(parms.k)
parms.C1k[:]      =  C1kfun(parms.k)
parms.C2k[:]      =  C2kfun(parms.k)
parms.CNNk[:]     =  CNNkfun(parms)
parms.Cϕϕk[:]     = Cϕϕkfun(parms)
parms.CZZmk[:]    =  CZZmkfun(parms.Cϕϕk, parms.Ck, parms.ηk, parms.ξk, parms.x, parms.k, parms.deltx, parms.deltk)
parms.CZZmobsk[:] = parms.CZZmk + parms.CNNk
parms.CZZmobsk[parms.r .< 1]              = Inf
parms.CZZmobsk[parms.r .>= 0.9*parms.nyq] = Inf
# Now a High resolution version of PhaseParms for simulation
parmsHR = PhaseParms(d_const, period_const, nside_const)
parmsHR.ηk[:]   =  ηkfun(parmsHR.k)
parmsHR.Ck[:]   =  Ckfun(parmsHR.k)
parmsHR.CNNk[:] =  CNNkfun(parmsHR)
parmsHR.Cϕϕk[:] =  Cϕϕkfun(parmsHR)

# ---- generate the initial simulation and other plotting quantities
zkobs, zk, zx, zx_noϕ, ϕk, ϕx = simNPhaseGRF(parms, parmsHR)
Dθx          = NonstationaryPhase.ϕx_2_Dθx(ϕx, parms)
Cℓvar        = Cℓvarfun(parms)
Cℓbias       = Cℓbiasfun(parms)
approxCℓbias = approxCℓbiasfun(parms)
estϕk        = estϕkfun(zkobs, parms)

two_sided_cutloci = NonstationaryPhase.cutloci(parms; tmax = 20, deltatmax = 0.01)

# --- get multiple quadratic estimates based on the same ϕx and the average estimate
EestDθx, Eestϕk, estimates_Dθx, avetime2 = get_estimates(100, ϕx, parms, parmsHR)
println("ave quadratic estimate clock speed =  $avetime2")

#########################################################
#=
Generate plots for the first set of simulations
=#
#########################################################

# --- now evaluate the above `script'
makeplots()
# redo the plot of the estimates to show the two sided cut locus
ax2 = subplot2grid((3,2), (1,0), colspan=2)
	plot(fftshift(parms.x[1]), fftshift(estimates_Dθx[:,1]),"k", label = L"\hat\theta^\prime(x)", linewidth = 1.5, alpha = 0.4)
	for cntr = 2:min(size(estimates_Dθx,2),5)
		plot(fftshift(parms.x[1]), fftshift(estimates_Dθx[:,cntr]),"k", linewidth = 1.5, alpha = 0.4)
	end
	plot(fftshift(parms.x[1]), fftshift(EestDθx), "k--", label = L"average $\hat\theta^\prime(x)$", linewidth = 3)
	plot(fftshift(parms.x[1]), fftshift(Dθx), "b", label = L"\theta^\prime(x)", linewidth = 2.5,alpha = 0.7)
	fill_between(fftshift(parms.x[1]), -two_sided_cutloci, two_sided_cutloci, facecolor="k", alpha=0.2)
	title(L"Quadratic estimate, average quadratic estimate and the simulation truth $\theta^\prime(x)$")
	axis("tight")
	legend(loc="best",fontsize = 10)

# --- save the figures to disk
if savefigures
    fig1[:savefig]("figure3.pdf", dpi=300, bbox_inches="tight", transparent=true)
    plt[:close](fig1)
end


text:
module MyModule
function complexfunc(a::Int, b::Int)::Float64
    x = rand(Float64, (a, a))
    for _ in 1:b
        # This is extremely inefficient, if you look at `htop` while
        # `async-tasks` is running you'll notice that three threads are doing
        # a lot of work: two threads are working on the task, the third is the
        # garbage collector.
        x += rand(Float64, (a, a))
    end

    z::Float64 = 0.0
    for j in 1:a
        z += x[j, j]
    end

    z
end
end

max_stars_repo_name: marhkb/jlrs

max_issues_repo_name: marhkb/jlrs

max_forks_repo_name: marhkb/jlrs


text:
using Gadfly

set_default_plot_size(14cm, 10cm)

plot(x=[5100, 5400], y=[0, 1], Geom.blank,
    Guide.yticks(ticks=0:0.1:1),
    Theme(grid_line_style=:solid, grid_line_width=0.5mm,
        grid_color="gray92", panel_stroke="black",
        grid_line_order=-2,  panel_line_width=1mm)
)


text:
"""
Delegates loading of code with functions and parameters specific to different instruments.
Subdirectories of src/instruments include provide functions specialized for each instrument,
typically the file I/O and pre-processing, so data ends up in a common format.
src/instruments/common.jl provides routines that can be shared by instruments.
"""

module InstrumentsCommon

using ..RvSpectMLBase
using DataFrames

include("common.jl")   # Mostly trait functions to be specialized by instruments

#include("io.jl")
#export read_manifest, read_header, read_metadata_from_fits

#include("linelists.jl")
#export read_linelist_espresso, read_linelist_vald

#include("masks.jl")
#export ChunkWidthFixedΔlnλ
#export read_mask_espresso, read_mask_vald

#export predict_intrinsic_stellar_line_width
#export λ_vac_to_air, λ_air_to_vac
#export find_overlapping_chunks
#export merge_chunks

include("param.jl")

export min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default, get_pixel_range
export orders_all, pixels_all, max_pixels_in_spectra       # generic implementations avaliable
export metadata_symbols_default, metadata_strings_default  # need to specialize
export default_ccf_mask_v_width
export get_inst_module

end # module InstrumentsCommon

# Now returning to module EchelleInstrments

#=
For each instrument/data file type, users need to create a sub-type of either AbstractInstrument2D or AbstractInstrument1D and
to implement the functions in neid/traits.jl for each instrument trait.
=#

include("theory/theory.jl")
import .TheoreticalInstrument: TheoreticalInstrument1D, TheoreticalInstrument2D, AnyTheoreticalInstrument
import .TheoreticalInstrument: get_inst_module
#import .TheoreticalInstrument: filter_line_list, find_worst_telluric_in_each_chunk
export TheoreticalInstrument, TheoreticalInstrument1D, TheoreticalInstrument2D, AnyTheoreticalInstrument
#=
import .TheoreticalInstrument: min_order, max_order, min_pixel_in_order, max_pixel_in_order
import .TheoreticalInstrument: orders_to_use_default, min_col_default, max_col_default
import .TheoreticalInstrument: orders_all, pixels_all, max_pixels_in_spectra       # generic implementations avaliable
import .TheoreticalInstrument: metadata_symbols_default, metadata_strings_default  # need to specialize
import .TheoreticalInstrument: default_ccf_mask_v_width
import .TheoreticalInstrument: get_inst_module
=#

#=
include("tellurics.jl")
export find_worst_telluric_in_each_chunk, make_clean_line_list_from_tellurics_expres
export filter_line_list
=#

import .InstrumentsCommon: pixels_all, min_pixel, max_pixel, max_pixels_in_spectra, min_col_default, max_col_default, bad_col_ranges, get_pixel_range
import .InstrumentsCommon: orders_all, min_order, max_order, orders_to_use_default
import .InstrumentsCommon: min_pixel_in_order, max_pixel_in_order, min_pixels_in_chunk
import .InstrumentsCommon: metadata_symbols_default, metadata_strings_default
import .InstrumentsCommon: default_ccf_mask_v_width
import .InstrumentsCommon: make_clean_line_list_from_tellurics, choose_obs_idx_for_init_guess

export min_order, max_order, min_pixel_in_order, max_pixel_in_order, get_pixel_range, bad_col_ranges
export orders_to_use_default, min_col_default, max_col_default
export orders_all, pixels_all, max_pixels_in_spectra       # generic implementations avaliable
export metadata_symbols_default, metadata_strings_default  # need to specialize
export default_ccf_mask_v_width
export make_clean_line_list_from_tellurics
export choose_obs_idx_for_init_guess
export get_inst_module

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
const _GW_TASK_OS_UPFREC = 5.0

## ------------------------------------------------------
function _run_task_os(gwt::GWTask)
    if is_worker_mode(gwt)

        # HANDLE STATUS
        _fatal_err(gwt) do

            _read_task_toml!(gwt)
            # check if pending
            curr_status = _get_task_status(gwt)
            !_is_pending_task(gwt) &&
                error("FATAL ERROR: Task status '$(curr_status)' != '$(_GW_TASK_PENDING_STATUS)'")

            # update
            _up_task_status!(gwt, _GW_TASK_RUNNING_STATUS)

            # PROC REGISTRY

            # SYS MAINTINANCE
            while true

                # PROC REGISTRY

                # FLUSH
                _flush()

                # SLEEP
                sleep(_GW_TASK_OS_UPFREC)
            end

        end

    end
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers