text:
import JSON

# JSON.jl doesn't define a serialization method for MIME and UUID objects, so we these ourselves:
import JSON: lower
JSON.lower(m::MIME) = string(m)
JSON.lower(u::UUID) = string(u)

function serialize_message_to_stream(io::IO, message::UpdateMessage)
    to_send = Dict(:type => message.type, :message => message.message)
    if message.notebook !== nothing
        to_send[:notebook_id] = message.notebook.notebook_id
    end
    if message.cell !== nothing
        to_send[:cell_id] = message.cell.cell_id
    end
    if message.initiator !== missing
        to_send[:initiator_id] = message.initiator.client_id
        to_send[:request_id] = message.initiator.request_id
    end

    JSON.print(io, to_send)
end

function serialize_message(message::UpdateMessage)
    sprint(serialize_message_to_stream, message)
end

function change_remote_cellinput!(notebook, cell, newcode; initiator::Union{Initiator,Missing}=missing)
    # i.e. Ctrl+Enter was pressed on this cell
    # we update our `Notebook` and start execution

    # don't reparse when code is identical (?)
    if cell.code != newcode
        cell.code = newcode
        cell.parsedcode = nothing
    end

    # TODO: feedback to user about File IO
    save_notebook(notebook)

    putnotebookupdates!(notebook, clientupdate_cell_input(notebook, cell, initiator=initiator))
end

responses[:connect] = (body, notebook = nothing; initiator::Union{Initiator,Missing}=missing)->begin
    putclientupdates!(initiator, UpdateMessage(:👋, Dict(
        :notebookExists => (notebook != nothing),
        :ENV => filter(p->startswith(p.first, "PLUTO"), ENV),
    ), nothing, nothing, initiator))
end

responses[:getversion] = (body, notebook = nothing; initiator::Union{Initiator,Missing}=missing)->begin
    putclientupdates!(initiator, UpdateMessage(:versioninfo, Dict(
        :pluto => PLUTO_VERSION_STR,
        :julia => JULIA_VERSION_STR,
    ), nothing, nothing, initiator))
end

# TODO: actions on the notebook are not thread safe
responses[:addcell] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    new_index = body["index"] + 1 # 0-based index (js) to 1-based index (julia)

    new_cell = Cell("")
    new_cell.output_repr = "" # we 'run' the code and get this output

    insert!(notebook.cells, new_index, new_cell)

    putnotebookupdates!(notebook, clientupdate_cell_added(notebook, new_cell, new_index, initiator=initiator))
end

responses[:deletecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    to_delete = cell

    # replace the cell's code with "" and do a reactive run
    change_remote_cellinput!(notebook, to_delete, "", initiator=initiator)
    runtask = run_reactive_async!(notebook, cell)

    # wait for the reactive run to finish, then delete the cells
    # we wait async, to make sure that the web server remains responsive
    @async begin
        wait(runtask)

        filter!(c->c.cell_id ≠ to_delete.cell_id, notebook.cells)
        putnotebookupdates!(notebook, clientupdate_cell_deleted(notebook, to_delete, initiator=initiator))
        save_notebook(notebook) # this might be "too late", but it will save the latest version of `notebook` anyways
    end
end

responses[:movecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    to_move = cell

    # Indexing works as if a new cell is added.
    # e.g. if the third cell (at julia-index 3) of [0, 1, 2, 3, 4]
    # is moved to the end, that would be new julia-index 6

    new_index = body["index"] + 1 # 0-based index (js) to 1-based index (julia)
    old_index = findfirst(isequal(to_move), notebook.cells)

    # Because our cells run in _topological_ order, we don't need to reevaluate anything.
    if new_index < old_index
        deleteat!(notebook.cells, old_index)
        insert!(notebook.cells, new_index, to_move)
    elseif new_index > old_index + 1
        insert!(notebook.cells, new_index, to_move)
        deleteat!(notebook.cells, old_index)
    end

    putnotebookupdates!(notebook, clientupdate_cell_moved(notebook, to_move, new_index, initiator=initiator))
end

responses[:changecell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    newcode = body["code"]

    change_remote_cellinput!(notebook, cell, newcode, initiator=initiator)
    run_reactive_async!(notebook, cell)
end

responses[:foldcell] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    newfolded = body["folded"]
    cell.code_folded = newfolded
    save_notebook(notebook)

    putnotebookupdates!(notebook, clientupdate_cell_folded(notebook, cell, newfolded, initiator=initiator))
end

responses[:run] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    run_reactive_async!(notebook, cell)
end

responses[:runmultiple] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    indices = cellindex_fromID.([notebook], UUID.(body["cells"]))
    cells = [notebook.cells[i] for i in indices if i !== nothing]
    run_reactive_async!(notebook, cells)
end

responses[:getinput] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    putclientupdates!(initiator, clientupdate_cell_input(notebook, cell, initiator=initiator))
end

responses[:setinput] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    change_remote_cellinput!(notebook, cell, body["code"], initiator=initiator)
end

responses[:getoutput] = (body, notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)->begin
    putclientupdates!(initiator, clientupdate_cell_output(notebook, cell, initiator=initiator))
end

responses[:getallcells] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    # TODO: the client's update channel might get full
    update = UpdateMessage(:cell_list,
        Dict(:cells => [Dict(
                :cell_id => string(cell.cell_id),
                ) for cell in notebook.cells]), nothing, nothing, initiator)

    putclientupdates!(initiator, update)
end

responses[:getallnotebooks] = (body, notebook = nothing; initiator::Union{Initiator,Missing}=missing)->begin
    putplutoupdates!(clientupdate_notebook_list(notebooks, initiator=initiator))
end

responses[:movenotebookfile] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    newpath = tamepath(body["path"])
    result = try
        if isfile(newpath)
            (success = false, reason = "File exists already - you need to delete the old file manually.")
        else
            move_notebook(notebook, newpath)
            putplutoupdates!(clientupdate_notebook_list(notebooks))
            (success = true, reason = "")
        end
    catch ex
        showerror(stderr, stacktrace(catch_backtrace()))
        (success = false, reason = sprint(showerror, ex))
    end

    update = UpdateMessage(:move_notebook_result, result, notebook, nothing, initiator)
    putclientupdates!(initiator, update)
end

responses[:interruptall] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    success = WorkspaceManager.interrupt_workspace(notebook)
    # TODO: notify user whether interrupt was successful (i.e. whether they are using a `ProcessWorkspace`)
end

responses[:shutdownworkspace] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    listeners = putnotebookupdates!(notebook) # TODO: shutdown message
    if body["remove_from_list"]
        delete!(notebooks, notebook.notebook_id)
        putplutoupdates!(clientupdate_notebook_list(notebooks))
        for client in listeners
            @async close(client.stream)
        end
    end
    success = WorkspaceManager.unmake_workspace(notebook)
end

responses[:bond_set] = (body, notebook::Notebook; initiator::Union{Initiator,Missing}=missing)->begin
    bound_sym = Symbol(body["sym"])
    new_val = body["val"]

    body["any_dependents"] = any_dependents = !isempty(where_assigned(notebook, Set{Symbol}([bound_sym])))
    putnotebookupdates!(notebook, UpdateMessage(:bond_update, body, notebook, nothing, initiator))

    if any_dependents
        function custom_deletion_hook(notebook::Notebook, to_delete_vars::Set{Symbol}, funcs_to_delete::Set{Vector{Symbol}}, to_reimport::Set{Expr}; to_run::Array{Cell,1})
            push!(to_delete_vars, bound_sym) # also delete the bound symbol
            WorkspaceManager.delete_vars(notebook, to_delete_vars, funcs_to_delete, to_reimport)
            WorkspaceManager.eval_in_workspace(notebook, :($bound_sym = $new_val))
        end
        to_reeval = where_referenced(notebook, Set{Symbol}([bound_sym]))
        run_reactive_async!(notebook, to_reeval; deletion_hook=custom_deletion_hook)
    end

end


text:
## A template for defining a new optimizer model

# Define a mutable struct for storing infinite model to optimizer model mappings
# plus other needed information
mutable struct NewReformData
    # Variable mapping (REPLACE AND REFORMAT BELOW AS NEEDED)
    infvar_to_vars::Dict{InfiniteVariableRef, Vector{JuMP.VariableRef}}
    holdvar_to_var::Dict{HoldVariableRef, JuMP.VariableRef}
    ptvar_to_var::Dict{PointVariableRef, JuMP.VariableRef}

    # Map other variable info (REPLACE AND REFORMAT BELOW AS NEEDED)
    infvar_to_supports::Dict{InfiniteVariableRef, Vector{<:Tuple}}

    # Constraint mapping (REPLACE AND REFORMAT BELOW AS NEEDED)
    infcon_to_constrs::Dict{InfiniteConstraintRef, Vector{JuMP.ConstraintRef}}
    meascon_to_constrs::Dict{MeasureConstraintRef, Vector{JuMP.ConstraintRef}}
    fincon_to_constr::Dict{FiniteConstraintRef, JuMP.ConstraintRef}

    # Map other constraint info (REPLACE AND REFORMAT BELOW AS NEEDED)
    infcon_to_supports::Dict{InfiniteConstraintRef, Vector{<:Tuple}}
    meascon_to_supports::Dict{MeasureConstraintRef, Vector{<:Tuple}}
    infcon_to_params::Dict{InfiniteConstraintRef, Tuple}
    meascon_to_params::Dict{MeasureConstraintRef, Tuple}

    # ADD OTHER MAPPING/MODEL INFORMATION HERE
    # default constructor
    function NewReformData() # REFORMAT BELOW IN ACCORDANCE WITH ABOVE ATTRIBUTES
        return new(Dict{InfiniteVariableRef, Vector{JuMP.VariableRef}}(),
                   Dict{HoldVariableRef, JuMP.VariableRef}(),
                   Dict{PointVariableRef, JuMP.VariableRef}(),
                   Dict{InfiniteVariableRef, Vector{Tuple}}(),
                   Dict{InfiniteConstraintRef, Vector{JuMP.ConstraintRef}}(),
                   Dict{MeasureConstraintRef, Vector{JuMP.ConstraintRef}}(),
                   Dict{FiniteConstraintRef, JuMP.ConstraintRef}(),
                   Dict{InfiniteConstraintRef, Vector{Tuple}}(),
                   Dict{MeasureConstraintRef, Vector{Tuple}}(),
                   Dict{InfiniteConstraintRef, Tuple}(),
                   Dict{MeasureConstraintRef, Tuple}())
    end
end

# Make a constructor for new optimizer model type (extension of JuMP.Model)
function NewReformModel(args...; kwargs...)::JuMP.Model # ADD EXPLICT ARGS AS NEEDED
    # initialize the JuMP Model
    model = JuMP.Model(args...; kwargs...)
    # ADD ADDITIONAL OPERATIONS IF NEEDED
    # add the extension data with a chosen optimizer model key and return
    model.ext[:ReformData] = NewReformData() # REPLACE WITH ACTUAL KEY AND DATA
    return model
end

# Extend build_optimizer_model! (enables build_optimizer_model! and optimize!)
function InfiniteOpt.build_optimizer_model!(model::InfiniteModel,
                                            key::Val{:ReformData}; # INSERT NEW KEY
                                            my_kwarg::Bool = true) # ADD KEYWORD ARGUMENTS AS NEEDED
    # clear the model for a build/rebuild
    reform_model = clear_optimizer_model_build!(model)

    # IT MAY BE USEFUL TO CALL expand_all_measures! TO HANDLE MEASURES FIRST

    # REPLACE LINES 62-69 WITH OPERATIONS TO BUILD A NewReformModel BASED ON `model`
    reform_model.ext[:TransData] = TranscriptionData()
    InfiniteOpt.TranscriptionOpt._build_transcription_model!(reform_model, model)
    reform_data = reform_model.ext[:ReformData]
    for i in eachindex(fieldnames(TranscriptionData))
        value = copy(getfield(transcription_data(reform_model), fieldnames(TranscriptionData)[i]))
        setfield!(reform_data, fieldnames(NewReformData)[i], value)
    end
    delete!(reform_model.ext, :TransData)

    # update the optimizer model status
    set_optimizer_model_ready(model, true)
    return
end

# Make function for extracting the data from the model (optional)
function reform_data(model::JuMP.Model)::NewReformData
    # REPLACE THE KEY BELOW WITH ACTUAL OPTIMIZER MODEL KEY
    haskey(model.ext, :ReformData) || error("Model is not a NewReformModel.")
    return model.ext[:ReformData]
end

# Extend optimizer_model_variable if appropriate to enable variable related queries
function InfiniteOpt.optimizer_model_variable(vref::InfOptVariableRef,
                                              key::Val{:ReformData}; # REPLACE WITH KEY
                                              my_kwarg::Bool = true) # ADD KEY ARGS AS NEEDED
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE OPTIMIZER MODEL VARIABLE(S)
    model = optimizer_model(JuMP.owner_model(vref))
    if vref isa InfiniteVariableRef
        map_dict = reform_data(model).infvar_to_vars
    elseif vref isa PointVariableRef
        map_dict = reform_data(model).ptvar_to_var
    else
        map_dict = reform_data(model).holdvar_to_var
    end
    haskey(map_dict, vref) || error("Variable $vref not used in the optimizer model.")
    return map_dict[vref]
end

# Extend optimizer_model_constraint if appropriate to enable constraint related queries
function InfiniteOpt.optimizer_model_constraint(cref::GeneralConstraintRef,
                                                key::Val{:ReformData}; # REPLACE WITH KEY
                                                my_kwarg::Bool = true) # ADD KEY ARGS AS NEEDED
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE OPTIMIZER MODEL CONSTRAINT(S)
    model = optimizer_model(JuMP.owner_model(cref))
    if cref isa InfiniteConstraintRef
        map_dict = reform_data(model).infcon_to_constrs
    elseif cref isa MeasureConstraintRef
        map_dict = reform_data(model).meascon_to_constrs
    else
        map_dict = reform_data(model).fincon_to_constr
    end
    haskey(map_dict, cref) || error("Constraint $cref not used in the optimizer model.")
    return map_dict[cref]
end

# If appropriate extend variable_supports (enables support queries of infinite variables)
function InfiniteOpt.variable_supports(model::JuMP.Model,
                                       vref::InfiniteVariableRef,
                                       key::Val{:ReformData}; # REPLACE WITH KEY
                                       my_kwarg::Bool = true) # ADD KEY ARGS AS NEEDED)
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    map_dict = reform_data(model).infvar_to_supports
    haskey(map_dict, vref) || error("Variable $vref not used in the optimizer model.")
    return map_dict[vref]
end

# If appropriate extend constraint_supports (enables support queries of infinite constraints)
function InfiniteOpt.constraint_supports(model::JuMP.Model,
                                         cref::GeneralConstraintRef,
                                         key::Val{:ReformData}; # REPLACE WITH KEY
                                         my_kwarg::Bool = true) # ADD KEY ARGS AS NEEDED)
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    if cref isa InfiniteConstraintRef
        map_dict = reform_data(model).infcon_to_supports
    else
        map_dict = reform_data(model).meascon_to_supports
    end
    haskey(map_dict, cref) || error("Constraint $cref does not have supports in the optimizer model.")
    return map_dict[cref]
end

# If appropriate extend constraint_parameter_refs (enables parameter reference queries of infinite constraints)
function InfiniteOpt.constraint_parameter_refs(model::JuMP.Model,
                                              cref::GeneralConstraintRef,
                                              key::Val{:ReformData}; # REPLACE WITH KEY
                                              my_kwarg::Bool = true) # ADD KEY ARGS AS NEEDED)
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    if cref isa InfiniteConstraintRef
        map_dict = reform_data(model).infcon_to_params
    else
        map_dict = reform_data(model).meascon_to_params
    end
    haskey(map_dict, cref) || error("Constraint $cref does not have parameter references in the optimizer model.")
    return map_dict[cref]
end

# If it doesn't make sense to extend optimizer_model_variable and/or optimizer_model_constraint
# then you'll need to extend the following:
# - InfiniteOpt.map_value to enable JuMP.value
# - InfiniteOpt.map_optimizer_index to enable JuMP.optimizer_index
# - InfiniteOpt.map_dual to enable JuMP.dual
# - InfiniteOpt.map_shadow_price to enable JuMP.shadow_price
# - InfiniteOpt.map_lp_rhs_perturbation_range to enable JuMP.lp_rhs_perturbation_range
# - InfiniteOpt.map_lp_objective_perturbation_range to enable JuMP.lp_objective_perturbation_range

ext: jl


text:
#Demo 2: Reflection coefficients
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import Graphics: width, height

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)
μ0 = 4pi*1e-7 #F/m
ϵ0 = 8.854e-12 #H/m

#Input data
#-------------------------------------------------------------------------------
fmax = 10e9
fstep = 50e6
f = collect(0:fstep:fmax)
ℓ = 50e-3 #Meters
ZL = Float64[1e6, 377, 60.0]
_colors = [blue, red, green]

#==Equations
===============================================================================#
Γ(Z; Zref::Real=50.0) = (Z .- Zref) ./ (Z .+ Zref)
#ZC: Characteristic impedance
#ZL: Load impendance (termination)
function Zline(ℓ::Real, f::Vector, ZL::Number; ZC::Number=50.0, α::Real=0, μ::Real=μ0, ϵ::Real=ϵ0)
	j = im
	β = f*(2pi*sqrt(μ*ϵ))
	γ = α .+ j*β
	tanh_γℓ = tanh.(γ*ℓ)
	return ZC*(ZL .+ ZC*tanh_γℓ)./(ZC .+ ZL*tanh_γℓ)
end
function Γline(ℓ::Real, f::Vector, ZL::Number; ZC::Number=50.0, Zref::Number=50.0, α::Real=0, μ::Real=μ0, ϵ::Real=ϵ0)
	return Γ(Zline(ℓ, f, ZL; ZC=ZC, α=α, μ=μ, ϵ=ϵ), Zref=Zref)
end

#Calculations
#-------------------------------------------------------------------------------
Γload = []
for ZLi in ZL
	_Γ = Γline(ℓ, f, ZLi, ZC=40)
	push!(Γload, _Γ)
end

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Transmission Line Example")
mplot.layout[:ncolumns] = 2

smithext = InspectDR.PExtents1D(min=-1.2,max=1.2) #Padded a bit

plot_linf = InspectDR.Plot2D(:lin, :dB20)
	graph_linf = plot_linf.strips[1]
	graph_linf.yext_full = InspectDR.PExtents1D(max=5)
#	plot_linf.layout[:enable_legend] = true

plot_logf = InspectDR.Plot2D(:log10, :dB20)
	graph_logf = plot_logf.strips[1]
	plot_logf.xext_full = InspectDR.PExtents1D(min=10e6) #Avoid issues with log scale
	graph_logf.yext_full = InspectDR.PExtents1D(max=5)
	graph_logf.grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=true)
plot_ysmith = InspectDR.smithchart(:Y, title="Y-Smith Chart")
	graph_ysmith = plot_ysmith.strips[1]
	plot_ysmith.xext_full = smithext
	graph_ysmith.yext_full = smithext
	plot_ysmith.layout[:enable_legend] = true
plot_zsmith = InspectDR.smithchart(:Z, ref=50, title="Z-Smith Chart")
	graph_zsmith = plot_zsmith.strips[1]
	plot_zsmith.xext_full = smithext
	graph_zsmith.yext_full = smithext
	plot_zsmith.layout[:enable_legend] = true

let plot, a #HIDEWARN_0.7
for plot in [plot_linf, plot_logf]
	a = plot.annotation
	a.title = "Reflection Coefficient (Γ)"
	a.xlabel = "Frequency (Hz)"
	a.ylabels = ["Magnitude (dB)"]
end
end

#Select which plots to actually display:
plotlist = [plot_linf, plot_logf, plot_ysmith, plot_zsmith]
#plotlist = [plot_zsmith]

let plot, wfrm #HIDEWARN_0.7
for plot in plotlist
	for i in 1:length(Γload)
		wfrm = add(plot, f, Γload[i], id="ZL=$(ZL[i])")
		wfrm.line = line(color=_colors[i], width=2)
	end

	add(mplot, plot)
end
end

gplot = display(InspectDR.GtkDisplay(), mplot)

#==Save multi-plot to file
===============================================================================#

maximize_square = true
if maximize_square
	#Target plot size to get square Smith plots without gaps:
	#(HACK: Method not suggested)
	lyt = plot_zsmith.layout
		lyt[:halloc_data] = 500
		lyt[:valloc_data] = 500
	bb = InspectDR.plotbounds(lyt.values, InspectDR.grid1(plot_zsmith)) #Required
		mplot.layout[:halloc_plot] = width(bb)
		mplot.layout[:valloc_plot] = height(bb)
end

if false
	#Test several attributes:
	mplot.frame.line.style=:solid
	mplot.frame.line.width=30
	mplot.frame.line.color=red
	mplot.frame.fillcolor=white
end

InspectDR.write_png("export_multiplot.png", mplot)
InspectDR.write_svg("export_multiplot.svg", mplot)
InspectDR.write_eps("export_multiplot.eps", mplot)
InspectDR.write_pdf("export_multiplot.pdf", mplot)

:DONE


text:
module TestSource

using TinyModia
using DifferentialEquations
@usingModiaPlot

SineSource = Model(
    f = 2.0,
    equations = :[y = sin(2*3.14*f*time)]
)

sineSource = @instantiateModel(SineSource, unitless=true)

simulate!(sineSource, Tsit5(), stopTime = 1.0, log=false, logStates=false)

plot(sineSource, "y")

end


text:
module kuramoto_sivashinsky

using SparseArrays
using LinearAlgebra

"""
    B = calcBernoulli(n)

Return the `n`th Bernoulli number.  Follows the convention that `n=1`
Bernoulli number is equal to -1/2.
"""
function calcBernoulli(n)
    A = Vector{Rational{BigInt}}(undef, n + 1)
    for m = 0:n
        A[m + 1] = 1 // (m + 1)
        for j = m:-1:1
            A[j] = j * (A[j] - A[j + 1])
        end
    end
    if n == 1
        return -A[1]
    else
        return A[1]
    end
end

"""
    quad = buildQuadrature(order [, T=Float64])

Return a Gregory-type quadrture rule for uniformly spaced data.
"""
function buildQuadrature(order::Int)
    # Form the linear system to solve
    numbnd = order - 1
    A = zeros(Rational{BigInt}, order-1, numbnd)
    b = zeros(Rational{BigInt}, order-1)
    for j = 1:order-1
        for i = 1:numbnd
            A[j,i] = j*(numbnd - (i-1))^(j-1)
        end
        b[j] = numbnd^j - ((-1)^j)*calcBernoulli(j)
    end
    w = A\b
    return w
    # if order <= 2
    #     quad = [1//2]
    # elseif order <= 4
    #     quad = [17//48; 59//48; 43//48; 49//48]
    # elseif order <= 6
    #     quad = [13649//43200; 12013//8640; 2711//4320; 5359//4320; 7877//8640;
    #             43801//43200]
    # else
    #     quad = [1498139//5080320; 1107307//725760; 20761//80640;
    #             1304999//725760; 299527//725760; 103097//80640;
    #             670091//725760; 5127739//5080320]
    # end
    # return convert(Array{T,1}, quad)
end

"""
Data structure for finite-difference operators on uniformly spaced points
"""
struct FiniteDiff{T}
    """Defines the nodes in the stencil; node 0 is rel. location of operator"""
    stencil::Array{Int,1}
    """Actual finite-difference coefficients"""
    coeffs::Array{T,1}
end

"""
    fd = buildFiniteDifference(deriv, order [, T=Float64, h=1.0])

Returns a finite difference operator struct for a derivative of order `deriv`
and order of accuracy of `order`. The optional argument `T` controls the
precision used, and `h` sets the (uniform) mesh spacing.
"""
function buildFiniteDifference(deriv::Int, order::Int, T::Type=Float64;
                               h=convert(T,1) )
    num_pts = deriv + order - 1
    if iseven(num_pts) num_pts += 1 end
    stencil = Array{Int}(-div(num_pts-1,2):div(num_pts-1,2))
    V = zeros(T, num_pts, num_pts)
    for i = 1:num_pts
        V[i,:] = stencil.^(i-1)
    end
    b = zeros(T, num_pts)
    b[deriv+1] = factorial(deriv)
    c = V\b
    c /= (h^deriv)
    return FiniteDiff{T}(stencil, c)
end

"""
    interp = buildInterp(order [, T=Float64])

Returns a centered interpolation operator of order `order` that does **not**
use data at the central node.
"""
function buildInterp(order::Int, T::Type=Float64)
    num_pts = order
    if isodd(num_pts) num_pts += 1 end
    stencil = Array{Int}([-div(num_pts,2):-1;1:div(num_pts,2)])
    V = zeros(T, num_pts, num_pts)
    for i = 1:num_pts
        V[i,:] = stencil.^(i-1)
    end
    b = zeros(T, num_pts)
    b[1] = 1.0
    c = V\b
    return FiniteDiff{T}(stencil, c)
end

"""
    num_ghost, num_bndry, ghost_op = buildGhostOp(order [, T=Float64])

Constructs a matrix operator `ghost_op` that, when applied to `num_bndry` nodes
at the boundary, produces `order` accurate `num_ghost` values.  **Note**, this
operator is designed to enforce a zero function value and zero derivative value
at the boundary, as per the Blonigan and Wang paper.

The ghost operator is constructed by solving a linearly constrained minimum
norm problem for the ghost node values.
"""
function buildGhostOp(order::Int, T::Type=Float64)
    # first, build the constraint Jacobian that enforces u=0 and du/dx = 0
    # A*u_ghost - B*u_true = c
    fd = buildFiniteDifference(1, order, T)
    #interp = buildInterp(order, T)
    num_nodes = length(fd.stencil)
    num_ghost = div(num_nodes-1,2) + 1 # this includes node at boundary
    num_bndry = num_nodes - num_ghost
    A = zeros(T, 2, num_ghost)
    A[1,end] = 1.0 # last node is the node on the boundary
    #A[1,1:end-1] = interp.coeffs[1:num_ghost-1]
    A[2,:] = fd.coeffs[1:num_ghost]
    B = zeros(T, 2, num_bndry)
    #B[1,:] = interp.coeffs[num_ghost:end]
    B[2,:] = fd.coeffs[num_ghost+1:end]
    return num_ghost, num_bndry, -A'*( (A*A') \ B )

    if false
        # This version uses a TVD-type regularization
        # first, build the constraint Jacobian that enforces u=0 and du/dx = 0
        # A*u_ghost = B*u_true + c
        fd = buildFiniteDifference(1, order, T)
        #interp = buildInterp(order, T)
        num_nodes = length(fd.stencil)
        num_ghost = div(num_nodes-1,2) + 1 # this includes node at boundary
        num_bndry = num_nodes - num_ghost
        A = zeros(T, 2, num_ghost)
        A[1,end] = 1.0 # last node
        #A[1,1:end-1] = interp.coeffs[1:num_ghost-1]
        A[2,:] = fd.coeffs[1:num_ghost]
        B = zeros(T, 2, num_bndry)
        #B[1,:] = interp.coeffs[num_ghost:end]
        B[2,:] = fd.coeffs[num_ghost+1:end]

        D = zeros(T, num_ghost, num_ghost)
        D += diagm(2*ones(T,num_ghost))
        D -= diagm(ones(T, num_ghost-1), 1)
        D -= diagm(ones(T, num_ghost-1), -1)
        K = [D A'; A zeros(T, 2, 2)]
        G = -inv(K)[1:num_ghost,end-1:end]*B
        return num_ghost, num_bndry, G
    end
end

"""
Data for KS function evaluation to avoid reallocation
"""
struct KSData{T}
    """Number of internal nodes (nodes we are solving for)"""
    num_nodes::Int
    """Number of ghost nodes at each end; include x=0 and x=L"""
    num_ghost::Int
    """Number of internal nodes used to determine ghost values"""
    num_bndry::Int
    """Stores the boundary first-derivative operator"""
    bndry_op::Array{T,2}
    """Stores the linear operatrors"""
    linear_op::SparseMatrixCSC{T,Int}
    """finite difference operator for the first-derivative"""
    fd::FiniteDiff{T}
end

"""
    ks = buildKSDataPeriodic(order, num_nodes [, T=Float64])

Returns a KSData structure that can be used to solve the Kuramoto-Sivashinsky
equations on _periodic domains_, while avoiding reallocation of arrays.
"""
function buildKSDataPeriodic(order::Int, num_nodes::Int, T::Type=Float64)
    # domain length used in Blonigan and Wang paper, but here periodic
    Lx = convert(T, 128)
    h = Lx/num_nodes
    # generate the matrix used for the linear operators
    A = zeros(T, num_nodes, num_nodes)
    # second, add the Laplacian contributions
    fd = buildFiniteDifference(2, order, T)
    fac = 1/(h*h)
    for i = 1:num_nodes
        A[i, mod.(fd.stencil .+ i .- 1, num_nodes) + 1] -= fd.coeffs*fac
    end
    # finally, add the BiLaplacian terms
    fd = buildFiniteDifference(4, order, T)
    fac *= fac
    for i = 1:num_nodes
        A[i, mod.(fd.stencil .+ i .- 1, num_nodes) + 1] -= fd.coeffs*fac
    end
    Asparse = sparse(A)
    # finally, get the first-derivative operator for nonlinear terms
    fd = kuramoto_sivashinsky.buildFiniteDifference(1, order, T, h=h)
    return KSData{T}(num_nodes, 0, 0, zeros(0,0), Asparse, fd)
end

"""
    getKSFuncPeriodic!(ks, t, u, du)

A function that returns the perodic spatial discretization of the
Kuramoto-Sivashinsky PDE for use in a time marching method.
"""
function getKSFuncPeriodic!(ks::KSData{Tks}, t, u::AbstractArray{T,1},
                            du::AbstractArray{T,1}) where {Tks, T}
    # apply the linear operators
    du[:] = ks.linear_op*u
    # add the nonlinear terms; this uses a skew-symmetric split-form
    third = convert(T, 1//3)
    for i = 1:ks.num_nodes
        for k = 1:size(ks.fd.stencil,1)
            j = mod(ks.fd.stencil[k]+i-1, ks.num_nodes) + 1
            du[i] -= third*u[i]*ks.fd.coeffs[k]*u[j]
            du[i] -= third*ks.fd.coeffs[k]*u[j]*u[j]
        end
        # This version uses central differencing for the nonlinear term
        #du[i] -= u[i]*dot(ks.fd.coeffs,
        #                  u[mod.(ks.fd.stencil+i-1,ks.num_nodes)+1])
    end
    return nothing
end

"""
    getKSJacPeriodic!(ks, t, u, Jac)

A function that returns the Jacobian of the perodic spatial discretization of
the Kuramoto-Sivashinsky PDE for use in a time marching method.
"""
function getKSJacPeriodic!(ks::KSData{T}, t, u::AbstractArray{T,1},
                           Jac::SparseMatrixCSC{T,Int}) where {T}
    Jac[:,:] = sparse(ks.linear_op)
    third = convert(T, 1//3)
    for i = 1:ks.num_nodes
        for k = 1:size(ks.fd.stencil,1)
            j = mod(ks.fd.stencil[k]+i-1, ks.num_nodes) + 1
            #du[i] -= fac*u[i]*ks.fd.coeffs[k]*u[j]
            Jac[i,i] -= third*ks.fd.coeffs[k]*u[j]
            Jac[i,j] -= third*u[i]*ks.fd.coeffs[k]
            #du[i] -= fac*ks.fd.coeffs[k]*u[j]*u[j]
            Jac[i,j] -= 2*third*ks.fd.coeffs[k]*u[j]
        end
        # This version uses central differencing for the nonlinear term
        #Jac[i,i] -= dot(ks.fd.coeffs,
        #                u[mod.(ks.fd.stencil+i-1,ks.num_nodes)+1])
        #Jac[i,mod.(ks.fd.stencil+i-1,ks.num_nodes)+1] -= u[i]*ks.fd.coeffs
    end
end

"""
    ks = buildKSData(order, num_nodes [, T=Float64])

Returns a KSData structure that can be used to solve the Kuramoto-Sivashinsky
equations while avoiding reallocation of arrays.  This version imposes
homogeneous Dirichlet and Neumann conditions at the ends of the spatial domain.
"""
function buildKSData(order::Int, num_nodes::Int, T::Type=Float64)
    # domain length used in Blonigan and Wang paper
    Lx = convert(T, 128)
    h = Lx/(num_nodes+1)
    # get the ghost node operator, call it Bg
    num_ghost, num_bndry, Bg = buildGhostOp(order, T)
    # println("Bg = ", Bg)

    # generate the matrix used for the linear operators
    A = zeros(T, num_nodes, num_nodes)
    # first, add the Laplacian contributions
    fd = buildFiniteDifference(2, order, T)
    fac = 1/(h*h)
    for i = 1:num_nodes
        idx = [k for k = 1:size(fd.stencil,1)
               if fd.stencil[k] + i >= 1 && fd.stencil[k] + i <= num_nodes]
        A[i, fd.stencil[idx] .+ i] -= fd.coeffs[idx]*fac
    end

    # eliminate ghost node contributions for Laplacian
    Ag = zeros(T, num_bndry, num_ghost)
    for i = 1:num_bndry
        idx = [k for k = 1:size(fd.stencil,1) if fd.stencil[k] + i < 1]
        Ag[i, fd.stencil[idx] .+ num_ghost .+ i] -= fd.coeffs[idx]*fac
    end
    # println("Ag (Laplacian) = ", Ag)
    G = Ag*Bg
    A[1:num_bndry,1:num_bndry] += G
    A[end:-1:end-num_bndry+1,end:-1:end-num_bndry+1] += G

    # next, add the BiLaplacian terms
    fd = buildFiniteDifference(4, order, T)
    fac = 1/(h*h*h*h)
    for i = 1:num_nodes
        idx = [k for k = 1:size(fd.stencil,1)
               if fd.stencil[k] + i >= 1 && fd.stencil[k] + i <= num_nodes]
        A[i, fd.stencil[idx] .+ i] -= fd.coeffs[idx]*fac
    end

    # eliminate ghost node contributions for BiLaplacian
    fill!(Ag, zero(T))
    for i = 1:num_bndry
        idx = [k for k = 1:size(fd.stencil,1) if fd.stencil[k] + i < 1]
        Ag[i, fd.stencil[idx] .+ num_ghost .+ i] -= fd.coeffs[idx]*fac
    end
    # println("Ag (BiLaplacian) = ", Ag)
    G = Ag*Bg
    A[1:num_bndry,1:num_bndry] += G
    A[end:-1:end-num_bndry+1,end:-1:end-num_bndry+1] += G

    # println("A[1,1:5] = ",A[1,1:5])
    # println("A[end,end-4:end] = ",A[end,end-4:end])

    Asparse = sparse(A)

    # finally, get the first-derivative operator for nonlinear terms
    fd = kuramoto_sivashinsky.buildFiniteDifference(1, order, T, h=h)
    # Get the first-derivative operator at the boundary
    fill!(Ag, zero(T))
    G = zeros(T, num_bndry, div(size(fd.stencil,1)+1,2) + num_bndry - 1 )
    for i = 1:num_bndry
        idx = [k for k = 1:size(fd.stencil,1) if fd.stencil[k] + i < 1]
        Ag[i, fd.stencil[idx] .+ num_ghost .+ i] += fd.coeffs[idx]
        idx = [k for k = 1:size(fd.stencil,1) if fd.stencil[k] + i >= 1]
        G[i, fd.stencil[idx] .+ i] += fd.coeffs[idx]
    end
    # println("Ag (first-derivative) = ", Ag)
    # println("G before = ",G)
    G[:,1:num_bndry] += Ag*Bg
    # println("G after = ",G)

    return KSData{T}(num_nodes, num_ghost, num_bndry, G, Asparse, fd)
end

"""
    getKSFunc!(ks, t, u, du)

A function that returns the spatial discretization of the
Kuramoto-Sivashinsky PDE for use in a time marching method.  Homogeneous BCs
for u and dudx are used at both ends of the domain.
"""
function getKSFunc!(ks::KSData{Tks}, t, u::AbstractArray{T,1},
                    du::AbstractArray{T,1}) where {Tks, T}
    # apply the linear operators
    du[:] = ks.linear_op*u
    # add the nonlinear terms to interior nodes
    fac = convert(T, 1//3)
    for i = ks.num_bndry + 1:ks.num_nodes - ks.num_bndry
        for k = 1:size(ks.fd.stencil,1)
            j = mod(ks.fd.stencil[k]+i-1, ks.num_nodes) + 1
            du[i] -= fac*u[i]*ks.fd.coeffs[k]*u[j]
            du[i] -= fac*ks.fd.coeffs[k]*u[j]*u[j]
        end
        #du[i] -= u[i]*dot(ks.fd.coeffs,
        #                  u[mod.(ks.fd.stencil+i-1,ks.num_nodes)+1])
    end
    # add the nonlinear terms to the boundary nodes
    for i = 1:size(ks.bndry_op,1)
        for k = 1:size(ks.bndry_op,2)
            du[i] -= fac*u[i]*ks.bndry_op[i,k]*u[k]
            du[i] -= fac*ks.bndry_op[i,k]*u[k]*u[k]
            du[end-i+1] += fac*u[end-i+1]*ks.bndry_op[i,k]*u[end-k+1]
            du[end-i+1] += fac*ks.bndry_op[i,k]*u[end-k+1]*u[end-k+1]
        end
    end
    return nothing
end

"""
    getKSJac!(ks, t, u, Jac)

A function that returns the Jacobian of the spatial discretization of
the Kuramoto-Sivashinsky PDE for use in a time marching method. Homogeneous BCs
for u and dudx are used at both ends of the domain.
"""
function getKSJac!(ks::KSData{T}, t, u::AbstractArray{T,1},
                   Jac::SparseMatrixCSC{T,Int}) where {T}
    # Jacobian contribution due to linear terms
    Jac[:,:] = sparse(ks.linear_op)
    #Jac[:,:] = deepcopy(ks.linear_op)
    # Jacobian due to nonlinear terms on interior
    fac = convert(T, 1//3)
    for i = ks.num_bndry+1:ks.num_nodes - ks.num_bndry
        for k = 1:size(ks.fd.stencil,1)
            j = mod(ks.fd.stencil[k]+i-1, ks.num_nodes) + 1
            # du[i] -= fac*u[i]*ks.fd.coeffs[k]*u[j]
            Jac[i,i] -= fac*ks.fd.coeffs[k]*u[j]
            Jac[i,j] -= fac*u[i]*ks.fd.coeffs[k]
            # du[i] -= fac*ks.fd.coeffs[k]*u[j]*u[j]
            Jac[i,j] -= 2*fac*ks.fd.coeffs[k]*u[j]
        end
        # This version uses central differencing for the nonlinear term
        #Jac[i,i] -= dot(ks.fd.coeffs,
        #                u[mod.(ks.fd.stencil+i-1,ks.num_nodes)+1])
        #Jac[i,mod.(ks.fd.stencil+i-1,ks.num_nodes)+1] -= u[i]*ks.fd.coeffs
    end
    # Jacobian due to nonlinear terms on boundary
    for i = 1:size(ks.bndry_op,1)
        for k = 1:size(ks.bndry_op,2)
            # du[i] -= fac*u[i]*ks.bndry_op[i,k]*u[k]
            Jac[i,i] -= fac*ks.bndry_op[i,k]*u[k]
            Jac[i,k] -= fac*u[i]*ks.bndry_op[i,k]
            # du[i] -= fac*ks.bndry_op[i,k]*u[k]*u[k]
            Jac[i,k] -= 2*fac*ks.bndry_op[i,k]*u[k]
            # du[end-i+1] += fac*u[end-i+1]*ks.bndry_op[i,k]*u[end-k+1]
            Jac[end-i+1,end-i+1] += fac*ks.bndry_op[i,k]*u[end-k+1]
            Jac[end-i+1,end-k+1] += fac*u[end-i+1]*ks.bndry_op[i,k]
            # du[end-i+1] += fac*ks.bndry_op[i,k]*u[end-k+1]*u[end-k+1]
            Jac[end-i+1,end-k+1] += 2*fac*ks.bndry_op[i,k]*u[end-k+1]

            # This version uses central differencing for the nonlinear term
            #du[i] -= u[i]*ks.bndry_op[i,k]*u[k]
            #  Jac[i,i] -= ks.bndry_op[i,k]*u[k]
            #  Jac[i,k] -= u[i]*ks.bndry_op[i,k]
            #du[end-i+1] -= u[end-i+1]*ks.bndry_op[i,k]*u[end-k+1]
            #  Jac[end-i+1,end-i+1] += ks.bndry_op[i,k]*u[end-k+1]
            #  Jac[end-i+1,end-k+1] += u[end-i+1]*ks.bndry_op[i,k]
        end
    end
end

"""
    Jac = getJacobian(u, ks)

This function uses the complex-step method to get the Jacobian; used for
verification of getKSJac!.  **Do not use** in simulations, due to cost.
"""
function getJacobian(u::AbstractArray{T,1}, ks::KSData{T}) where{T}
    fd_c = FiniteDiff{Complex128}(ks.fd.stencil, ks.fd.coeffs)
    ks_c = KSData{Complex128}(ks.num_nodes, ks.num_ghost, ks.num_bndry,
                              ks.bndry_op, ks.linear_op, fd_c)
    Jac = zeros(size(u,1),size(u,1))
    u_c = complex.(u, 0.0)
    du_c = zeros(u_c)
    ceps = 1e-60
    for i = 1:ks.num_nodes
        u_c[i] += complex(0.0, ceps)
        getKSFunc!(ks_c, 0.0, u_c, du_c)
        Jac[:,i] = imag.(du_c)/ceps
        u_c[i] -= complex(0.0, ceps)
    end
    return Jac
end

"""
    sol = solveUsingMidpoint(ks, Time, num_stpes, u)

Solve the KS problem on the nonperiodic domain using the midpoint rule.  The
problem is run from t=0 to t=`Time` using `num_steps` steps.  The initial
condition is `u`.  Returns `sol`, the solution over all space and temporal
nodes; `sol[i,n]` stores the `i`th space node at time step `n`.  Note that
`sol` stores the boundary nodes, whereas `u` does not.
"""
function solveUsingMidpoint(ks::KSData{T}, Time::T, num_steps::Int,
                            u::Array{T,1}) where {T}
    num_nodes = size(u,1)
    # storage for solution over all time steps
    dt = Time/num_steps
    sol = zeros(num_nodes+2, num_steps+1)
    sol[2:end-1,1] = u

    t = convert(T, 0.0)
    Jac = similar(ks.linear_op)
    F = lu(ks.linear_op) # initialize memory for the lu factorization
    u_mid = zeros(size(u))
    u_old = zeros(size(u))
    r = zeros(size(u))
    newt_tol = 1e-10 # We will need to ensure this is sufficiently small
    max_newt_iter = 50
    for n = 1:num_steps
        #println("Iteration ",n,": time ",t)
        # Newton iterations
        u_old = u
        norm0 = one(T)
        for k = 1:max_newt_iter
            u_mid = 0.5*u + 0.5*u_old
            kuramoto_sivashinsky.getKSFunc!(ks, t, u_mid, r)
            r = u - u_old - dt*r
            #println("\tNewton iter ",k,": res norm = ",norm(r))
            if k == 1
                norm0 = norm(r)
            else
                if norm(r) < newt_tol*norm0
                    break
                end
            end
            kuramoto_sivashinsky.getKSJac!(ks, t, u_mid, Jac)
            Jac *= 0.5*dt
            Jac -= I #speye(T, num_nodes)
            lu!(F, Jac)
            du = F\r
            u += du
        end
        t += dt
        sol[2:end-1,n+1] = u
    end
    return sol
end

"""
    avg = calcSolutionAverage(order, sol, func)

Returns the space and time average of `func(sol[:,:])` using an `order` order
accurate diagonal norm scheme.
"""
function calcSolutionAverage(order::Int, sol::AbstractArray{T,2},
                             func::Function) where {T}
    avg = zero(T)
    quad = convert(Vector{T}, buildQuadrature(order))
    num_bndry = size(quad,1)
    for n = 1:size(sol,2)
        wt = one(T)
        if n <= num_bndry
            wt *= quad[n]
        elseif n >= size(sol,2) - num_bndry + 1
            wt *= quad[size(sol,2)-n+1]
        end
        for i = 1:size(sol,1)
            if i <= num_bndry
                avg += func(sol[i,n])*wt*quad[i]
            elseif i >= size(sol,1) - num_bndry + 1
                avg += func(sol[i,n])*wt*quad[size(sol,1)-i+1]
            else
                avg += func(sol[i,n])*wt
            end
        end
    end
    avg /= (size(sol,1)-1)*(size(sol,2)-1)
    return avg
end

# function factorBandedSparse(A::SparseMatrixCSC{T,Int64}) where {T}

#     fac = one(T)
#     # loop over the elementary matrix operations
#     for k = 1:A.m-1
#         # update column k of the lower triangular matrix
#         # This assumes the rows are sorted, which Julia matrices should be
#         for ptr = A.colptr[k]:A.colptr[k+1]-1
#             if rowval[ptr] > k
#                 continue
#             elseif rowval[ptr] == k
#                 fac = one(T)/A.nzval[ptr]
#             else
#                 A.nzval[ptr] *= fac
#             end
#         end
#         # apply the elementary matrix to the lower block matrix
#         for j = k+1:A.n
#             if A.rowval[A.colptr[j]] > k
#                 break # there are no further nonzeros in row k beyond col j
#             end
#             for ptr = A.colptr[j]:A.colptr[j+1]-1
#                 i = A.rowval[ptr]
#                 if i <= k continue end
#                 A.nzval[ptr] -= A[i,k]*A[k,j] # this is not efficient
#             end
#         end
#     end
# end

end # module

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos


text:
using GLFW

include("GLAbstraction/GLAbstraction.jl")

using .GLAbstraction

const atlas_texture_cache = Dict{Any, Tuple{Texture{Float16, 2}, Function}}()

function get_texture!(atlas)
    # clean up dead context!
    filter!(atlas_texture_cache) do (ctx, tex_func)
        if GLAbstraction.context_alive(ctx)
            true
        else
            AbstractPlotting.remove_font_render_callback!(tex_func[2])
            false
        end
    end
    tex, func = get!(atlas_texture_cache, GLAbstraction.current_context()) do
        tex = Texture(
                atlas.data,
                minfilter = :linear,
                magfilter = :linear,
                # TODO: Consider alternatives to using the builtin anisotropic
                # samplers for signed distance fields; the anisotropic
                # filtering should happen *after* the SDF thresholding, but
                # with the builtin sampler it happens before.
                anisotropic = 16f0,
                mipmap = true
        )
        # update the texture, whenever a new font is added to the atlas
        function callback(distance_field, rectangle)
            ctx = tex.context
            if GLAbstraction.context_alive(ctx)
                prev_ctx = GLAbstraction.current_context()
                ShaderAbstractions.switch_context!(ctx)
                tex[rectangle] = distance_field
                ShaderAbstractions.switch_context!(prev_ctx)
            end
        end
        AbstractPlotting.font_render_callback!(callback)
        return (tex, callback)
    end
    return tex
end

include("GLVisualize/GLVisualize.jl")
using .GLVisualize

include("glwindow.jl")
include("screen.jl")
include("rendering.jl")
include("events.jl")
include("drawing_primitives.jl")

function AbstractPlotting.backend_display(x::GLBackend, scene::Scene)
    screen = global_gl_screen(size(scene), AbstractPlotting.use_display[])
    display_loading_image(screen)
    AbstractPlotting.backend_display(screen, scene)
    return screen
end

"""
    scene2image(scene::Scene)

Buffers the `scene` in an image buffer.
"""
function scene2image(scene::Scene)
    screen = global_gl_screen(size(scene), false)
    AbstractPlotting.backend_display(screen, scene)
    AbstractPlotting.colorbuffer(screen)
end

raw_io(io::IO) = io
raw_io(io::IOContext) = raw_io(io.io)

function AbstractPlotting.backend_show(::GLBackend, io::IO, m::MIME"image/png", scene::Scene)
    img = scene2image(scene)
    FileIO.save(FileIO.Stream(FileIO.format"PNG", raw_io(io)), img)
end

function AbstractPlotting.backend_show(::GLBackend, io::IO, m::MIME"image/jpeg", scene::Scene)
    img = scene2image(scene)
    FileIO.save(FileIO.Stream(FileIO.format"JPEG", raw_io(io)), img)
end


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

JLLWrappers.@generate_main_file_header("LinuxKernelHeaders")
JLLWrappers.@generate_main_file("LinuxKernelHeaders", UUID("d754bdc3-4c87-55c0-a004-f56830624a1e"))
end  # module LinuxKernelHeaders_jll


text:
using JamesPlots
using Test
using CairoMakie

function create_plot()
    fig = Figure()
    ax = MyAxis(fig[1,1], xlabel="numbers", ylabel="more numbers")
    a = rand(100)
    lines!(ax, a)
    scatter!(ax, a)
    b = rand(100)
    lines!(ax, b)
    scatter!(ax, b)
    fig
end

save_figure(joinpath(@__DIR__, "test.pdf"), create_plot)


text:
"""
Enumerations for, well, SAC header enumerations. Names are the same as those
defined in the SAC manual. See
<http://ds.iris.edu/files/sac-manual/manual/file_format.html> for details on
each enum. For a list of values, run `instances(HeaderEnum)`.
"""
@enum(HeaderEnum, undefined=-12345, itime=1, irlim, iamph, ixy, iunkn, idisp, ivel, iacc, ib,
      iday, io, ia, it0, it1, it2, it3, it4, it5, it6, it7, it8, it9, iradnv,
      itannv, iradev, itanev, inorth, ieast, ihorza, idown, iup, illlbb, iwwsn1,
      iwwsn2, ihglp, isro, inucl, ipren, ipostn, iquake, ipreq, ipostq, ichem,
      iother, igood, iglch, idrop, ilowsn, irldta, ivolts, imb=52, ims, iml, imw,
      imd, imx, ineic, ipdeq, ipdew, ipde, iisc, ireb, iusgs, ibrk, icaltech,
      illnl, ievloc, ijsop, iuser, iunknown, iqb, iqb1, iqb2, iqbx, iqmt, ieq,
      ieq1, ieq2, ime, iex, inu, inc, io_, il, ir, it, iu, ieq3, ieq0, iex0,
      iqc, iqb0, igey, ilit, imet, iodor, ios=103)

# Undefined values for different data types, as given in the SAC manual.
const HEADER_UNDEFINED_VAL = Dict{Type,Any}(Float32     => Float32(-12345.0),
                                            Int32       => Int32(-12345),
                                            HeaderEnum  => undefined,
                                            Bool        => false,
                                            String => "-12345  " :: String)

const SAC_HDR_VERSION = 6
const SAC_HDR_NWORDS = 158
# Starting characters for each section of the header (1 indexed). Includes
# newline characters.
const SAC_HDR_ASCII_FLOAT_START = 1
const SAC_HDR_ASCII_INT_START = 1065
const SAC_HDR_ASCII_ENUM_START = 1218
const SAC_HDR_ASCII_BOOL_START = 1422
const SAC_HDR_ASCII_ALPHA_START = 1473
const SAC_HDR_ASCII_END = 1672

"""
Description
===========

`Header` represents the header of a SAC file. Fields are initialised with
the default undefined values for Floats, Ints, Enumerations, Bools and Strings
as described in the SAC manual.

Fields
======

Fields have the same name as their header variables in SAC. See
<http://ds.iris.edu/files/sac-manual/manual/file_format.html> for details on
each one. The different header types become the following Julia types:

- Floating     => Float32
- Integer      => Int32
- Enumerated   => HeaderEnum
- Logical      => Bool
- Alphanumeric => String

Constructors
============

The best way to create a new header is to use:

`Header(npts::Int32, beginning::Float32, ending::Float32, ftype::HeaderEnum, even::Bool, delta::Float32; version::Int32=Int32(6))`

This will create a *valid* SAC header with the required fields initialised and
other fields set to their undefined values. You can create an empty header using
`Header()` but required fields will be initialised to their undefined
values so the header will be invalid.

See Also
========

`HeaderEnum` for information on the possible enumerated values.
"""
type Header
    "Increment between evenly spaced samples."
    delta::Float32
    "Minimum value of the dependent variable."
    depmin::Float32
    "Maximum value of the dependent variable."
    depmax::Float32
    "Scale factor for dependent variable (not used by SAC)."
    scale::Float32
    "Observed increment if different from `delta`."
    odelta::Float32
    "Beginning value of the independent variable."
    b::Float32
    "Ending value of the independent variable."
    e::Float32
    "Event origin time relative to the reference time in seconds."
    o::Float32
    "First arrival time relative to reference time in seconds."
    a::Float32
    internal1::Float32
    "User defined time pick or marker."
    t0::Float32
    "User defined time pick or marker."
    t1::Float32
    "User defined time pick or marker."
    t2::Float32
    "User defined time pick or marker."
    t3::Float32
    "User defined time pick or marker."
    t4::Float32
    "User defined time pick or marker."
    t5::Float32
    "User defined time pick or marker."
    t6::Float32
    "User defined time pick or marker."
    t7::Float32
    "User defined time pick or marker."
    t8::Float32
    "User defined time pick or marker."
    t9::Float32
    "End of event time in seconds relative to reference time."
    f::Float32
    "Instrument response parameter (not used by SAC)."
    resp0::Float32
    "Instrument response parameter (not used by SAC)."
    resp1::Float32
    "Instrument response parameter (not used by SAC)."
    resp2::Float32
    "Instrument response parameter (not used by SAC)."
    resp3::Float32
    "Instrument response parameter (not used by SAC)."
    resp4::Float32
    "Instrument response parameter (not used by SAC)."
    resp5::Float32
    "Instrument response parameter (not used by SAC)."
    resp6::Float32
    "Instrument response parameter (not used by SAC)."
    resp7::Float32
    "Instrument response parameter (not used by SAC)."
    resp8::Float32
    "Instrument response parameter (not used by SAC)."
    resp9::Float32
    "Station latitude with degrees north being positive."
    stla::Float32
    "Station longitude with degrees east being positive."
    stlo::Float32
    "Station elevation above sea level in metres."
    stel::Float32
    "Station depth below the surface in metres."
    stdp::Float32
    "Event latitude with degrees north being positive."
    evla::Float32
    "Event longitude with degrees east being positive."
    evlo::Float32
    "Event elevation in metres."
    evel::Float32
    "Event depth below surface in kilometres."
    evdp::Float32
    "Event magnitude."
    mag::Float32
    "User defined variable storage area."
    user0::Float32
    "User defined variable storage area."
    user1::Float32
    "User defined variable storage area."
    user2::Float32
    "User defined variable storage area."
    user3::Float32
    "User defined variable storage area."
    user4::Float32
    "User defined variable storage area."
    user5::Float32
    "User defined variable storage area."
    user6::Float32
    "User defined variable storage area."
    user7::Float32
    "User defined variable storage area."
    user8::Float32
    "User defined variable storage area."
    user9::Float32
    "Station to event distance in kilometres."
    dist::Float32
    "Event to station azimuth in degrees."
    az::Float32
    "Station to event azimuth in degrees."
    baz::Float32
    "Station to event great circle length (degrees)>"
    gcarc::Float32
    internal2::Float32
    internal3::Float32
    "Mean value of the dependent variable."
    depmen::Float32
    "Component azimuth in degrees measured clockwise from north."
    cmpaz::Float32
    "Component incident angle measured in degrees from the vertical."
    cmpinc::Float32
    "Minimum value of X (spectral files only)."
    xminimum::Float32
    "Maximum value of X (spectral files only)."
    xmaximum::Float32
    "Minimum value of Y (spectral files only)."
    yminimum::Float32
    "Maximum value of Y (spectral files only)."
    ymaximum::Float32
    unused1::Float32
    unused2::Float32
    unused3::Float32
    unused4::Float32
    unused5::Float32
    unused6::Float32
    unused7::Float32
    # Integers
    "GMT year corresponding to reference (zero) time in file."
    nzyear::Int32
    "GMT julian day corresponding to reference (zero) time in file."
    nzjday::Int32
    "GMT hour corresponding to reference (zero) time in file."
    nzhour::Int32
    "GMT minute corresponding to reference (zero) time in file."
    nzmin::Int32
    "GMT second corresponding to reference (zero) time in file."
    nzsec::Int32
    "GMT millisecond corresponding to reference (zero) time in file."
    nzmsec::Int32
    "Header version number"
    nvhdr::Int32
    "Origin ID (CSS 3.0)"
    norid::Int32
    "Event ID (CSS 3.0)"
    nevid::Int32
    "Number of points per data component"
    npts::Int32
    internal4::Int32
    "Waveform ID (CSS 3.0)"
    nwfid::Int32
    "Spectral length (spectral files only)."
    nxsize::Int32
    "Spectral width (spectral files only)."
    nysize::Int32
    unused8::Int32
    # Enumerations
    """Type of file. Must be one of:

     - `itime` (Time series file)
     - `irlim` (Spectral file with real and imaginary components)
     - `iamph` (Spectral file with amplitude and phase components)
     - `ixy` (General x vs. y data)
     - `ixyz` (General xyz data)
    """
    iftype::HeaderEnum
    """Type of dependent variable. Possible values include:

    - `iunkn` (Unknown)
    - `idisp` (Displacement in nm)
    - `ivel` (Velocity in nm/sec)
    - `ivolts` (Velocity in volts)
    - `iacc` (Acceleration in nm/sec/sec)
    """
    idep::HeaderEnum
    """Reference time type. Possible values include:

    - `iunkn` (Unknown)
    - `ib` (Begin time)
    - `iday` (Midnight of reference GMT day)
    - `io` (Event origin time)
    - `ia` (First arrival time)
    - `itn` where n = 0:9. (User defined time pick)
    """
    iztype::HeaderEnum
    unused9::HeaderEnum
    "Type of recording instrument (not used by SAC)."
    iinst::HeaderEnum
    "Station geographic region (not used by SAC)."
    istreg::HeaderEnum
    "Event geographic region (not used by SAC)."
    ievreg::HeaderEnum
    """Type of event. One of:

    - `iunkn` (Unknown)
    - `inucl` (Nuclear event)
    - `ipren` (Nuclear pre-shot event)
    - `ipostn` (Nuclear post-shot event)
    - `iquake` (Earthquake)
    - `ipreq` (Foreshock)
    - `ipostq` (Aftershock)
    - `ichem` (Chemical explosion)
    - `iqb` (Quarry or mine blast confirmed by quarry)
    - `iqb1` (Quarry/mine blast with designed shot info-ripple fired)
    - `iqb2` (Quarry/mine blast with observed shot info-ripple fired)
    - `iqbx` (Quarry or mine blast - single shot)
    - `iqmt` (Quarry/mining-induced events: tremors and rockbursts)
    - `ieq` (Earthquake)
    - `ieq1` (Earthquakes in a swarm or aftershock sequence)
    - `ieq2` (Felt earthquake)
    - `ime` (Marine explosion)
    - `iex` (Other explosion)
    - `inu` (Nuclear explosion)
    - `inc` (Nuclear cavity collapse)
    - `io` (Other source of unknown origin)
    - `il` (Local event of unknown origin)
    - `ir` (Regional event of unknown origin)
    - `it` (Teleseismic event of unknown origin)
    - `iu` (Undetermined or conflicting information)
    - `iother` (Other)
    """
    ievtyp::HeaderEnum
    """Quality of data (not used by SAC). Can be one of:

    - `igood` (Good data)
    - `iglch` (Glitches)
    - `idrop` (Dropouts)
    - `ilowsn` (Low signal to noise ratio)
    - `iother` (Other)
    """
    iqual::HeaderEnum
    "Synthetic data flag (not used by SAC). Can be `irldta` (real data)"
    isynth::HeaderEnum
    """Magnitude type. One of:

    - `imb` (Body wave magnitude)
    - `ims` (Surface wave magnitude)
    - `iml` (Local magnitude)
    - `imw` (Moment magnitude)
    - `imd` (Duration magnitude)
    - `imx` (User defined magnitude)
    """
    imagtyp::HeaderEnum
    """Source of magnitude information. One of:

    - `ineic` (National Earthquake Information Center)
    - `ipde` (Preliminary Determination of Epicentre)
    - `iisc` (International Seismological Centre)
    - `ireb` (Review Event Bulletin)
    - `iusgs` (United States Geological Survey)
    - `ibrk` (UC Berkeley)
    - `icaltech` (California Institute of Technology)
    - `illnl` (Lawrence Livermore National Laboratory)
    - `ievloc` (Event Location (computer program))
    - `ijsop` (Join Seismic Observation Program)
    - `iuser` (The individual using SAC2000)
    - `iunknown` (Unknown)
    """
    imagsrc::HeaderEnum
    unused10::HeaderEnum
    unused11::HeaderEnum
    unused12::HeaderEnum
    unused13::HeaderEnum
    unused14::HeaderEnum
    unused15::HeaderEnum
    unused16::HeaderEnum
    unused17::HeaderEnum
    # Logical Variables
    "`true` if data is evenly spaced."
    leven::Bool
    "`true` if station components have a positive polarity (left-hand rule)."
    lpspol::Bool
    "`true` if it is OK to overwrite this file on disk."
    lovrok::Bool
    "True if `dist`, `az`, `baz` and `gcarc` are to be calculated from station
    and event coordinates."
    lcalda::Bool
    unused18::Bool
    # Alphanumerics
    "Station name."
    kstnm::String
    "Event name."
    kevnm::String
    "If nuclear related then the hole identifier. Otherwise, the location identifier."
    khole::String
    "Event origin time identification."
    ko::String
    "First arrival time identification."
    ka::String
    "User defined time pick identifier."
    kt0::String
    "User defined time pick identifier."
    kt1::String
    "User defined time pick identifier."
    kt2::String
    "User defined time pick identifier."
    kt3::String
    "User defined time pick identifier."
    kt4::String
    "User defined time pick identifier."
    kt5::String
    "User defined time pick identifier."
    kt6::String
    "User defined time pick identifier."
    kt7::String
    "User defined time pick identifier."
    kt8::String
    "User defined time pick identifier."
    kt9::String
    "End of event time identification."
    kf::String
    "User defined variable storage area."
    kuser0::String
    "User defined variable storage area."
    kuser1::String
    "User defined variable storage area."
    kuser2::String
    "Channel name."
    kcmpnm::String
    "Name of the seismic network."
    knetwk::String
    "Date data was read onto computer."
    kdatrd::String
    "Generic name of recording instrument."
    kinst::String

    function Header(npts::Integer, b::AbstractFloat, e::AbstractFloat, ftype::HeaderEnum,
                    even::Bool, delta::AbstractFloat)
        hdr = Header()
        hdr.npts = npts
        hdr.b = b
        hdr.e = e
        hdr.iftype = ftype
        hdr.leven = even
        hdr.delta = delta
        hdr.nvhdr = 6

        return hdr
    end

    function Header(npts::Integer, b::AbstractFloat, e::AbstractFloat, ftype::HeaderEnum, even::Bool)
        # This exists because the SAC documentation says that the `delta` header
        # variable is required even though it doesn't make sense for a lot of
        # data types. With this, we can just initialise it to an undefined
        # variable, which is what SAC seems to do.
        Header(npts, b, e, ftype, even, HEADER_UNDEFINED_VAL[Float32])
    end

    Header() = (hdr = new(); set_undefinedvars!(hdr); cleanhdr!(hdr);)
end

"""
Set all the fields of `hdr::Header` to their undefined values as
specified in the SAC manual.
"""
function set_undefinedvars!(hdr::Header)
    for (field, T) in zip(fieldnames(Header), Header.types)
        setfield!(hdr, field, HEADER_UNDEFINED_VAL[T])
    end
    return hdr
end

"""
    cleanhdr!(hdr::Header)

Clean the string fields of `hdr`. This removes all leading/trailing whitespace
and strips null characters.
"""
function cleanhdr!(hdr::Header)
    for field in fieldnames(Header)[end-22:end]
        setfield!(hdr, field, strip(getfield(hdr, field)))
        setfield!(hdr, field, replace(getfield(hdr, field), '\0', ""))
    end

    return hdr
end


text:
@testset "Deprecations" begin
    @testset "copy/similar" begin
        img = ImageMeta(rand(3,5); prop1 = 1, prop2 = [1,2,3])
        img2 = copy(img)
        @test img2.data == img.data
        img2[2,2] = -1
        @test img2[2,2] < 0
        @test img[2,2] >= 0
        img2["prop2"][2] = -2
        @test img2["prop2"] == [1,-2,3]
        @test img["prop2"] == [1,2,3]
        img2 = similar(img)
        @test img2["prop1"] == 1
        @test img2["prop2"] == [1,2,3]
        @test img2.data != img.data
        img2["prop3"] = 7
        @test !haskey(img, "prop3")
        img2 = similar(img, RGB{Float16}, (Base.OneTo(5),))
        @test img2["prop1"] == 1
        @test img2["prop2"] == [1,2,3]
        @test eltype(img2) == RGB{Float16}
        @test size(img2) == (5,)
        @test img2.data != img.data
        img2["prop3"] = 7
        @test !haskey(img, "prop3")

        A = AxisArray(rand(3,5), :y, :x)
        B = ImageMeta(A, info="blah")
        C = similar(B)
        @test isa(C, ImageMeta)
        @test isa(C.data, AxisArray)
        @test eltype(C) == Float64
        C = similar(B, RGB{Float16})
        @test isa(C, ImageMeta)
        @test isa(C.data, AxisArray)
        @test eltype(C) == RGB{Float16}
    end

    @testset "copy/shareproperties/viewim" begin
        img = ImageMeta(rand(3,5); prop1 = 1, prop2 = [1,2,3])
        @test !isempty(properties(img))
        v = view(img, 1:2, 1:2)
        c = img[1:2, 1:2]
        @test v["prop1"] == 1
        @test c["prop1"] == 1
        img2 = copyproperties(img, reshape(1:15, 5, 3))
        @test size(img2) == (5,3)
        img2["prop1"] = -1
        @test img["prop1"] == 1
        img2 = shareproperties(img, reshape(1:15, 5, 3))
        @test size(img2) == (5,3)
        img2["prop1"] = -1
        @test img["prop1"] == -1
        @test v["prop1"] == -1
        @test c["prop1"] == 1
        imgb = ImageMeta(rand(RGB{N0f8}, 2, 2), propa = "hello", propb = [1,2])
        copy!(img, imgb, "propa", "propb")
        @test img["propa"] == "hello"
        @test img["propb"] == [1,2]
        img["propb"][2] = 10
        @test img["propb"] == [1,10]
        @test imgb["propb"] == [1,2]
        delete!(img, "propb")
        @test  haskey(img, "propa")
        @test !haskey(img, "propb")

        img = ImageMeta(AxisArray(rand(3,5,8),
                                  Axis{:x}(1:3),
                                  Axis{:y}(1:5),
                                  Axis{:time}(0.1:0.1:0.8));
                        prop1 = 1, prop2 = [1,2,3])
        v = view(img, Axis{:time}(0.25..0.5))
        c = img[Axis{:time}(0.25..0.5)]
        @test v["prop1"] == 1
        @test c["prop1"] == 1
    end

    @testset "traits" begin
        img = ImageMeta(AxisArray(rand(3,5,8),
                                  Axis{:x}(1:3),
                                  Axis{:y}(1:5),
                                  Axis{:time}(0.1:0.1:0.8)))
        @test @inferred(timedim(img)) == 3
        @test @inferred(pixelspacing(img)) === (1,1)
        @test spacedirections(img) === ((1,0),(0,1))
        img["spacedirections"] = "broken"
        @test spacedirections(img) == "broken"
        @test sdims(img) == 2
        @test @inferred(coords_spatial(img)) == (1,2)
        @test spatialorder(img) == (:x, :y)
        @test nimages(img) == 8
        @test @inferred(size_spatial(img)) == (3,5)
        @test @inferred(indices_spatial(img)) == (Base.OneTo(3), Base.OneTo(5))
        assert_timedim_last(img)
        @test timeaxis(img) == Axis{:time}(0.1:0.1:0.8)
    end

    @testset "spatialprops" begin
        img = ImageMeta(rand(3,5),
                        spatialproperties=Set(["vector","matrix","tuple"]),
                        vector=[1,2],
                        matrix=[1 3; 2 4],
                        tuple=(1,2))
        for imgp in (img', permutedims(img, (2,1)), permuteddimsview(img, (2,1)))
            @test imgp.data == img.data'
            @test imgp["vector"] == [2,1]
            @test imgp["matrix"] == [4 2; 3 1]
            @test imgp["tuple"]  == (2,1)
            @test img["vector"] == [1,2]
            @test img["matrix"] == [1 3; 2 4]
            @test img["tuple"]  == (1,2)
        end
    end

    @testset "AbstractDicts" begin
        img = ImageMeta(AxisArray(rand(3,5,8),
                                  Axis{:x}(1:3),
                                  Axis{:y}(1:5),
                                  Axis{:time}(0.1:0.1:0.8)),
                        IdDict{String,Any}("pixelspacing" => (1,1)))
        @test @inferred(pixelspacing(img)) === (1,1)

        rand_img = rand(10,10)
        dict = Dict("attr1" => 3, "attr2" => 4)
        @test ImageMeta(rand_img, dict)[1, 1] == rand_img[1, 1]
    end
end

nothing


text:
using Base.Test

tests = [
    "timearray",
    "split",
    "apply",
    "combine",
    "modify",
    "meta",
    "readwrite",
    "timeseriesrc",
    "deprecated",
]

@testset "TimeSeries" begin
    println("Running tests:")

    for test ∈ tests
        println("\t* $test ...")
        include("$test.jl")
    end
end


text:
@testset "$TEST $G" begin

vm = ConstVertexMap(0)
@test typeof(vm) <: AVertexMap
@test get(vm, 1, 1) == 0
@test get(vm, -1, 1) == 0
vm[1] = 1
@test vm[1] == 0
@test haskey(vm, 1)
@test haskey(vm, -1)
@test length(vm) == typemax(Int)

g = G()
vm = VertexMap(g, rand(1:10,10))
@test typeof(vm) <: AVertexMap
@test haskey(vm, 2)
@test !haskey(vm, -1)
@test valtype(vm) == Int
@test get(vm, 1, -100) != -100
@test length(Vector(vm)) == 0

g = G(10,20)
v = rand(10)
vm = VertexMap(g, v)
@test Vector(vm) == v

vm = VertexMap(g, Float64)
@test valtype(vm) == Float64
@test typeof(vm) <: AVertexMap

vm[1] = 2.
@test vm[1] == 2
@test haskey(vm, 1)
@test !haskey(vm, 2)
@test get(vm, -1, -100) == -100

vm2 = deepcopy(vm)
@test vm2 == vm

# @test sprint(show, vm) == "VertexMap: $(vm.data)"

g = G(10,20)
m = VertexMap(g, i -> i^2)
@test typeof(m.data) == Vector{V}
@test length(m.data) == 10
@test length(m) == 10
for i=1:10
    @test m[i] == i^2
end

m = VertexMap(g, i -> rand(2))
@test typeof(m.data) == Vector{Vector{Float64}}
@test length(m.data) == length(m) == nv(g)
@test  eltype(m) == valtype(m) == eltype(m.data) == Vector{Float64}

end # testset


text:
## This script is modified from:
##   `https://github.com/evetion/GeoArrays.jl/blob/master/src/crs.jl`
## Copyright (c) 2018 Maarten Pronk, MIT license

const TypeINDEX = Union{Colon,AbstractRange,Integer}
"""
    getindex(ga::AbstractGeoArray, i::AbstractRange, j::AbstractRange, k::TypeINDEX = :)
    TypeINDEX = Union{Colon,AbstractRange,Integer}

Index a `AbstractGeoArray` with `AbstractRange`s to get a cropped
AbstractGeoArray with the correct `AffineMap` set.

# Examples

```julia-repl
julia> ga[2:3,2:3,1]
2x2x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [1.0, 1.0]) and undefined CRS
```
"""
function Base.getindex(ga::AbstractGeoArray, i::AbstractRange, j::AbstractRange, k::TypeINDEX = :)
    A = @view ga.A[i, j, k]
    x, y = first(i) - 1, first(j) - 1
    t = ga.f(SVector(x, y))

    GeoArray(A, AffineMap(ga.f.linear, t), ga.crs)
end

function Base.getindex(ga::AbstractGeoArray, name::AbstractString)
    k = findall(ga.names .== name)
    rast(@view(ga.A[:, :, k]), ga.f, ga.crs, ga.names[k])
end

function Base.getindex(ga::AbstractGeoArray, k::TypeINDEX=:)
    GeoArray(ga, vals=ga.A[:, :, k])
end

function Base.getindex(ga::AbstractGeoArray, ::Colon, ::Colon, args...)
    rast(ga, vals = @view(ga.A[:, :, args...]))
end

# Base.getindex(ga::AbstractGeoArray, I::Vararg{<:Integer,2}) = begin
#     getindex(ga.A, I..., :)
# end
# Base.getindex(ga::AbstractGeoArray, I::Vararg{<:Integer,3}) = begin
#     println("d3", I)
#     getindex(ga.A, I...)
# end
Base.setindex!(ga::AbstractGeoArray, v, I::Vararg{<:AbstractFloat,2}) = setindex!(ga, v, SVector{2}(I))
Base.setindex!(ga::AbstractGeoArray, v, I::Vararg{TypeINDEX,2}) = setindex!(ga.A, v, I..., :)
Base.setindex!(ga::AbstractGeoArray, v, I::Vararg{TypeINDEX,3}) = setindex!(ga.A, v, I...)


text:
using JuMP
using SCIP

solver = GndStr.CallbackTopo(
    JuMP.optimizer_with_attributes(
        SCIP.Optimizer,
        "display/width" => 139,
        "limits/memory" => 5000.0
    ),
    JuMP.optimizer_with_attributes(
        SCIP.Optimizer,
        "display/width" => 139,
        "limits/memory" => 5000.0
    ),
    timelimit=3600.0)


text:
"""
    mean(nums)

Find mean of a vector of numbers

# Example

```julia
mean([3, 6, 9, 12, 15, 18, 21])      # returns 12.0
mean([5, 10, 15, 20, 25, 30, 35])    # returns 20.0
mean([1, 2, 3, 4, 5, 6, 7, 8])       # returns 4.5
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function mean(nums)
    return sum(nums) / length(nums)
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
#=
The series, 1^1 + 2^2 + 3^3 + ... + 10^10 = 10405071317.

Find the last ten digits of the series, 1^1 + 2^2 + 3^3 + ... + 1000^1000.
=#
function calc()
  a = [i^BigInt(i) for i in 1:1000]
  d = reverse(digits(sum(a))[1:10])
  result = ""
  for digit in d
    result *= "$digit"
  end
  result
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
struct NILSSSensitivityFunction{iip,F,Alg,
     PGPU,PGPP,CONFU,CONGP,DGVAL,DG,jType,RefType} <: DiffEqBase.AbstractODEFunction{iip}
  f::F
  alg::Alg
  numparams::Int
  numindvar::Int
  pgpu::PGPU
  pgpp::PGPP
  pgpu_config::CONFU
  pgpp_config::CONGP
  dg_val::DGVAL
  dg::DG
  jevery::jType # if concrete_solve interface for discrete cost functions is used
  cur_time::RefType
end

function NILSSSensitivityFunction(sensealg,f,u0,p,tspan,g,dg,jevery=nothing,cur_time=nothing)

  numparams = length(p)
  numindvar = length(u0)

  # compute gradients of objective
  if dg !== nothing
    pgpu = nothing
    pgpp = nothing
    pgpu_config = nothing
    pgpp_config = nothing
    if dg isa Tuple && length(dg) == 2
      dg_val = (similar(u0, numindvar),similar(u0, numparams))
      dg_val[1] .= false
      dg_val[2] .= false
    else
      dg_val = similar(u0, numindvar) # number of funcs size
      dg_val .= false
    end
  else
    pgpu = UGradientWrapper(g,tspan[1],p) # ∂g∂u
    pgpp = ParamGradientWrapper(g,tspan[1],u0) #∂g∂p
    pgpu_config = build_grad_config(sensealg,pgpu,u0,tspan[1])
    pgpp_config = build_grad_config(sensealg,pgpp,u0,tspan[1])
    dg_val = (similar(u0, numindvar),similar(u0, numparams))
    dg_val[1] .= false
    dg_val[2] .= false
  end

  NILSSSensitivityFunction{isinplace(f),typeof(f),typeof(sensealg),
                             typeof(pgpu),typeof(pgpp),typeof(pgpu_config),typeof(pgpp_config),typeof(dg_val),typeof(dg),typeof(jevery),typeof(cur_time)}(
                             f,sensealg,numparams,numindvar,pgpu,pgpp,pgpu_config,pgpp_config,dg_val,dg,jevery,cur_time)
end

struct NILSSProblem{A,CacheType,FSprob,probType,u0Type,vstar0Type,w0Type,
    TType,dtType,gType,yType,vstarType,
    wType,RType,bType,weightType,CType,dType,BType,aType,vType,xiType,
    G,DG,resType}
  sensealg::A
  diffcache::CacheType
  forward_prob::FSprob
  prob::probType
  u0::u0Type
  vstar0::vstar0Type
  w0::w0Type
  nus::Int
  T_seg::TType
  dtsave::dtType
  gsave::gType
  y::yType
  dudt::yType
  dgdu::yType
  vstar::vstarType
  vstar_perp::vstarType
  w::wType
  w_perp::wType
  R::RType
  b::bType
  weight::weightType
  Cinv::CType
  d::dType
  B::BType
  a::aType
  v::vType
  v_perp::vType
  ξ::xiType
  g::G
  dg::DG
  res::resType
end

function NILSSProblem(prob, sensealg::NILSS, t=nothing, dg = nothing;
                            kwargs...)

  @unpack f, p, u0, tspan = prob
  @unpack nseg, nstep, nus, rng, g = sensealg  #number of segments on time interval, number of steps saved on each segment

  numindvar = length(u0)
  numparams = length(p)

  # some shadowing sensealgs require knowledge of g
  check_for_g(sensealg,g)

  # integer dimension of the unstable subspace
  if nus === nothing
    nus = numindvar - one(numindvar)
  end
  (nus >= numindvar) && error("`nus` must be smaller than `numindvar`.")

  isinplace = DiffEqBase.isinplace(f)

  p === nothing && error("You must have parameters to use parameter sensitivity calculations!")
  !(u0 isa AbstractVector) && error("`u` has to be an AbstractVector.")

  # segmentation: determine length of segmentation and spacing between saved points
  T_seg = (tspan[2]-tspan[1])/nseg # length of each segment
  dtsave = T_seg/(nstep-1)

  # assert that dtsave is chosen such that all ts are hit if concrete solve interface/discrete costs are used
  if t!==nothing
    @assert t isa StepRangeLen
    dt_ts = step(t)
    @assert dt_ts >= dtsave
    @assert T_seg >= dt_ts
    jevery = Int(dt_ts/dtsave) # will throw an inexact error if dt_ts is not a multiple of dtsave. (could be more sophisticated)
    cur_time = Ref(1)
  else
    jevery = nothing
    cur_time = nothing
  end

  # inhomogenous forward sensitivity problem
  chunk_size = determine_chunksize(numparams,sensealg)
  autodiff = alg_autodiff(sensealg)
  difftype = diff_type(sensealg)
  autojacvec = sensealg.autojacvec
  # homogenous + inhomogenous forward sensitivity problems
  forward_prob = ODEForwardSensitivityProblem(f,u0,tspan,p,ForwardSensitivity(chunk_size=chunk_size,autodiff=autodiff,
                                                diff_type=difftype,autojacvec=autojacvec);nus=nus, kwargs...)

  sense = NILSSSensitivityFunction(sensealg,f,u0,p,tspan,g,dg,jevery,cur_time)

  # pre-allocate variables
  gsave = Matrix{eltype(u0)}(undef, nstep, nseg)
  y = Array{eltype(u0)}(undef, numindvar, nstep, nseg)
  dudt = similar(y)
  dgdu = similar(y)
  vstar = Array{eltype(u0)}(undef, numparams, numindvar, nstep, nseg) # generalization for several parameters numindvar*numparams
  vstar_perp = similar(vstar)
  w = Array{eltype(u0)}(undef, numindvar, nstep, nseg, nus)
  w_perp = similar(w)

  # assign initial values to y, v*, w
  y[:,1,1] .= u0
  for i=1:numparams
    _vstar = @view vstar[i,:,1,1]
    copyto!(_vstar, zero(u0))
  end
  for ius=1:nus
    _w = @view w[:,1,1,ius]
    rand!(rng,_w)
    normalize!(_w)
  end

  vstar0 = zeros(eltype(u0), numindvar*numparams)
  w0 = vec(w[:,1,1,:])

  R = Array{eltype(u0)}(undef, numparams, nseg-1, nus, nus)
  b = Array{eltype(u0)}(undef, numparams, (nseg-1)*nus)

  # a weight matrix for integration, 0.5 at interfaces
  weight = ones(1,nstep)
  weight[1] /= 2
  weight[end] /= 2

  # Construct Schur complement of the Lagrange multiplier method of the NILSS problem.
  # See the paper on FD-NILSS
  # find C^-1
  Cinv = Matrix{eltype(u0)}(undef, nseg*nus, nseg*nus)
  Cinv .*= false
  d = Vector{eltype(u0)}(undef, nseg*nus)
  B = Matrix{eltype(u0)}(undef, (nseg-1)*nus, nseg*nus)
  B .*= false

  a = Vector{eltype(u0)}(undef, nseg*nus)
  v = Array{eltype(u0)}(undef, numindvar, nstep, nseg)
  v_perp = similar(v)

  # only need to use last step in each segment
  ξ = Matrix{eltype(u0)}(undef, nseg, 2)

  res = similar(u0, numparams)

  NILSSProblem{typeof(sensealg),typeof(sense),typeof(forward_prob),typeof(prob),
    typeof(u0), typeof(vstar0), typeof(w0),
    typeof(T_seg),typeof(dtsave),typeof(gsave),typeof(y),typeof(vstar),typeof(w),typeof(R),
    typeof(b),typeof(weight),typeof(Cinv),typeof(d),typeof(B),typeof(a),typeof(v),typeof(ξ),
    typeof(g),typeof(dg),typeof(res)}(sensealg,sense,forward_prob,prob,u0,vstar0,w0,
    nus,T_seg,dtsave,gsave,y,dudt,dgdu,vstar,vstar_perp,w,w_perp,R,b,weight,Cinv,d,
    B,a,v,v_perp,ξ,g,dg,res)
end

function (NS::NILSSForwardSensitivityFunction)(du,u,p,t)
  @unpack S, nus = NS
  y = @view u[1:S.numindvar] # These are the independent variables
  dy = @view du[1:S.numindvar]
  S.f(dy,y,p,t) # Make the first part be the ODE

  # Now do sensitivities
  # Compute the Jacobian

  if !S.isautojacvec
    if has_original_jac(S.f)
      S.original_jac(S.J,y,p,t) # Calculate the Jacobian into J
    else
      S.uf.t = t
      jacobian!(S.J, S.uf, y, S.f_cache, S.alg, S.jac_config)
    end
  end

  if DiffEqBase.has_paramjac(S.f)
    S.paramjac(S.pJ,y,p,t) # Calculate the parameter Jacobian into pJ
  else
    S.pf.t = t
    S.pf.u .= y
    jacobian!(S.pJ, S.pf, p, S.f_cache, S.alg, S.paramjac_config)
  end

  # Compute the parameter derivatives
  for j=1:nus+1
    for i in eachindex(p)
      indx1 = (j-1)*S.numindvar*1 + i*S.numindvar+1
      indx2 = (j-1)*S.numindvar*1 + (i+1)*S.numindvar
      Sj = @view u[indx1:indx2]
      dp = @view du[indx1:indx2]
      if !S.isautojacvec
        mul!(dp,S.J,Sj)
      else
        jacobianvec!(dp, S.uf, y, Sj, S.alg, S.jac_config)
      end
      if j == nus+1
        # inhomogenous (otherwise homogenous tangent solution)
        dp .+= @view S.pJ[:,i]
      end
    end
  end
  return nothing
end

function forward_sense(prob::NILSSProblem,nilss::NILSS,alg)
  #TODO determine a good dtsave (ΔT in paper, see Sec.4.2)
  @unpack nus, T_seg, dtsave, vstar, vstar_perp, w, w_perp, R, b, y, dudt, gsave, dgdu, forward_prob, u0, vstar0, w0 = prob
  @unpack p, f = forward_prob
  @unpack S, sensealg = f
  @unpack nseg, nstep = nilss
  @unpack numindvar, numparams = S

  # push forward
  t1 = forward_prob.tspan[1]
  t2 = forward_prob.tspan[1]+T_seg
  _prob = ODEForwardSensitivityProblem(S.f,u0,(t1,t2),p,sensealg;nus=nus,w0=w0,v0=vstar0)

  for iseg=1:nseg
    # compute y, w, vstar
    # _sol is a numindvar*(1+nus+1) x nstep matrix

    dt = (t2 - t1) / (nstep-1)
    _sol = Array(solve(_prob, alg, saveat=t1:dt:t2))

    store_y_w_vstar!(y, w, vstar, _sol, nus, numindvar, numparams, iseg)

    # store dudt, objective g (gsave), and its derivative wrt. to u (dgdu)
    dudt_g_dgdu!(dudt, gsave, dgdu, prob, y, forward_prob.p, iseg)

    # calculate w_perp, vstar_perp
    perp!(w_perp, vstar_perp, w, vstar, dudt, iseg, numparams, nstep, nus)

    # update sense problem
    if iseg < nseg
      # renormalize at interfaces
      renormalize!(R,b,w_perp,vstar_perp,y,vstar,w,iseg,numparams,nus)
      t1 = forward_prob.tspan[1]+iseg*T_seg
      t2 = forward_prob.tspan[1]+(iseg+1)*T_seg
      _prob = ODEForwardSensitivityProblem(S.f,y[:,1,iseg+1],(t1,t2),p,sensealg; nus=nus,
                                 w0=vec(w[:,1,iseg+1,:]),v0=vec(vstar[:,:,1,iseg+1]))
    end

  end
end

function store_y_w_vstar!(y, w, vstar, sol, nus, numindvar, numparams, iseg)
  # fill y
  _y = @view y[:,:,iseg]
  copyto!(_y, (@view sol[1:numindvar,:]))

  # fill w
  # only calculate w one time, w can be reused for each parameter
  for j=1:nus
   indx1 = (j-1)*numindvar*1 + numindvar+1
   indx2 = (j-1)*numindvar*1 + 2*numindvar

   _w = @view w[:,:,iseg, j]
   copyto!(_w, (@view sol[indx1:indx2,:]))
  end

  # fill vstar
  for i=1:numparams
    indx1 = nus*numindvar*1 + i*numindvar+1
    indx2 = nus*numindvar*1 + (i+1)*numindvar
    _vstar = @view vstar[i,:,:,iseg]
    copyto!(_vstar, (@view sol[indx1:indx2,:]))
  end

  return nothing
end

function dudt_g_dgdu!(dudt, gsave, dgdu, nilssprob::NILSSProblem, y, p, iseg)
  @unpack sensealg, diffcache, dg, g, prob = nilssprob
  @unpack prob = nilssprob
  @unpack jevery, cur_time = diffcache # akin to ``discrete"

  _y = @view y[:,:,iseg]

  for (j,u) in enumerate(eachcol(_y))
    _dgdu = @view dgdu[:,j,iseg]
    _dudt = @view dudt[:,j,iseg]

    # compute dudt
    if isinplace(prob)
      prob.f(_dudt,u,p,nothing)
    else
      copyto!(_dudt,prob.f(u,p,nothing))
    end

    # compute objective
    gsave[j,iseg] = g(u,p,nothing)

    #  compute gradient of objective wrt. state
    if jevery!==nothing
      # only bump on every jevery entry
      # corresponds to (iseg-1)* value of dg
      if (j-1) % jevery == 0
        accumulate_cost!(_dgdu, dg, u, p, nothing, sensealg, diffcache, cur_time[])
        cur_time[] += one(jevery)
      end
    else
      # continuous cost function
      accumulate_cost!(_dgdu, dg, u, p, nothing, sensealg, diffcache, j)
    end
  end
  jevery !== nothing && (cur_time[] -= one(jevery)) # interface between segments gets two bumps
  return nothing
end

function perp!(w_perp, vstar_perp, w, vstar, dudt, iseg, numparams, nsteps, nus)
  for indx_steps=1:nsteps
    _dudt = @view dudt[:,indx_steps,iseg]
    for indx_nus=1:nus
      _w_perp = @view w_perp[:,indx_steps,iseg,indx_nus]
      _w = @view w[:,indx_steps,iseg,indx_nus]
      perp!(_w_perp, _w, _dudt)
    end
    for indx_params=1:numparams
      _vstar_perp = @view vstar_perp[indx_params,:,indx_steps,iseg]
      _vstar = @view vstar[indx_params,:,indx_steps,iseg]
      perp!(_vstar_perp, _vstar, _dudt)
    end
  end

  return nothing
end

function perp!(v1, v2, v3)
  v1 .= v2 - dot(v2, v3)/dot(v3, v3) * v3
end

function renormalize!(R,b,w_perp,vstar_perp,y,vstar,w,iseg,numparams,nus)
  for i=1:numparams
    _b = @view b[i,(iseg-1)*nus+1:iseg*nus]
    _R = @view R[i,iseg,:,:]
    _w_perp = @view w_perp[:,end,iseg,:]
    _vstar_perp = @view vstar_perp[i,:,end,iseg]
    _w = @view w[:,1,iseg+1,:]
    _vstar = @view vstar[i,:,1,iseg+1]

    Q_temp, R_temp = qr(_w_perp)
    b_tmp = @view (Q_temp'*_vstar_perp)[1:nus]

    copyto!(_b, b_tmp)
    copyto!(_R, R_temp)
    # set new initial values
    copyto!(_w, (@view Q_temp[:,1:nus]))
    copyto!(_vstar, _vstar_perp - Q_temp*b_tmp)
  end
  _yend = @view y[:,end,iseg]
  _ystart = @view y[:,1,iseg+1]
  copyto!(_ystart, _yend)

  return nothing
end

function compute_Cinv!(Cinv,w_perp,weight,nseg,nus,indxp)
  # construct Schur complement of Lagrange multiplier
  _weight = @view weight[1,:]
  for iseg=1:nseg
    _C = @view Cinv[(iseg-1)*nus+1:iseg*nus, (iseg-1)*nus+1:iseg*nus]
    for i=1:nus
      wi = @view w_perp[:,:,iseg,i]
      for j =1:nus
        wj = @view w_perp[:,:,iseg,j]
        _C[i,j] = sum(wi .* wj * _weight)
      end
    end
    invC = inv(_C)
    copyto!(_C, invC)
  end
  return nothing
end

function compute_d!(d,w_perp,vstar_perp,weight,nseg,nus,indxp)
  # construct d
  _weight = @view weight[1,:]
  for iseg=1:nseg
    _d = @view d[(iseg-1)*nus+1:iseg*nus]
    vi = @view vstar_perp[indxp,:,:,iseg]
    for i=1:nus
      wi = @view w_perp[:,:,iseg,i]
      _d[i] = sum(wi .* vi * _weight)
    end
  end
  return nothing
end

function compute_B!(B,R,nseg,nus,indxp)
  for iseg=1:nseg-1
    _B = @view B[(iseg-1)*nus+1:iseg*nus, (iseg-1)*nus+1:iseg*nus]
    _R = @view R[indxp,iseg,:,:]
    copyto!(_B, -_R)
    # off diagonal one
    for i=1:nus
      B[(iseg-1)*nus+i, iseg*nus+i] = one(eltype(R))
    end
  end
  return nothing
end

function compute_a!(a,B,Cinv,b,d,indxp)
  _b = @view b[indxp,:]

  lbd = (-B*Cinv*B') \ (B*Cinv*d + _b)
  a .= -Cinv*(B'*lbd + d)
  return nothing
end

function compute_v!(v,v_perp,vstar,vstar_perp,w,w_perp,a,nseg,nus,indxp)
  _vstar = @view vstar[indxp,:,:,:]
  _vstar_perp = @view vstar_perp[indxp,:,:,:]

  copyto!(v, _vstar)
  copyto!(v_perp, _vstar_perp)

  for iseg=1:nseg
    vi = @view v[:,:,iseg]
    vpi = @view v_perp[:,:,iseg]
    for i=1:nus
      wi = @view w[:,:,iseg,i]
      wpi = @view w_perp[:,:,iseg,i]

      vi .+= a[(iseg-1)*nus+i]*wi
      vpi .+= a[(iseg-1)*nus+i]*wpi
    end
  end

  return nothing
end

function compute_xi(ξ,v,dudt,nseg)
  for iseg=1:nseg
    _v = @view v[:,1,iseg]
    _dudt = @view dudt[:,1,iseg]
    ξ[iseg,1] = dot(_v,_dudt)/dot(_dudt,_dudt)

    _v = @view v[:,end,iseg]
    _dudt = @view dudt[:,end,iseg]
    ξ[iseg,2] = dot(_v,_dudt)/dot(_dudt,_dudt)
  end
  # check if segmentation is chosen correctly
  _ξ = ξ[:,1]
  all(_ξ.<1e-4) || @warn "Detected a large value of ξ at the beginning of a segment."
  return nothing
end

function accumulate_cost!(_dgdu, dg, u, p, t, sensealg::NILSS, diffcache::NILSSSensitivityFunction, j)
  @unpack dg_val, pgpu, pgpu_config, pgpp, pgpp_config = diffcache

  if dg===nothing
    if dg_val isa Tuple
      DiffEqSensitivity.gradient!(dg_val[1], pgpu, u, sensealg, pgpu_config)
      copyto!(_dgdu, dg_val[1])
    else
      DiffEqSensitivity.gradient!(dg_val, pgpu, u, sensealg, pgpu_config)
      copyto!(_dgdu, dg_val)
    end
  else
    if dg_val isa Tuple
      dg[1](dg_val[1],u,p,nothing,j)
      @. _dgdu = -dg_val[1]
    else
      dg(dg_val,u,p,nothing,j)
      @. _dgdu = -dg_val
    end
  end

  return nothing
end

function shadow_forward(prob::NILSSProblem,alg; sensealg=prob.sensealg)
  shadow_forward(prob,sensealg,alg)
end

function shadow_forward(prob::NILSSProblem,sensealg::NILSS,alg)
  @unpack nseg, nstep = sensealg
  @unpack res, nus, dtsave, vstar, vstar_perp, w, w_perp, R, b, dudt,
          gsave, dgdu, forward_prob, weight, Cinv, d, B, a, v, v_perp, ξ = prob
  @unpack numindvar, numparams = forward_prob.f.S

  # reset dg pointer
  @unpack jevery, cur_time = prob.diffcache
  jevery !== nothing && (cur_time[] = one(jevery))

  # compute vstar, w
  forward_sense(prob,sensealg,alg)

  # compute avg objective
  gavg = sum(prob.weight*gsave)/((nstep-1)*nseg)

  # reset gradient
  res .*= false

  # loop over parameters
  for i=1:numparams
    compute_Cinv!(Cinv,w_perp,weight,nseg,nus,i)
    compute_d!(d,w_perp,vstar_perp,weight,nseg,nus,i)
    compute_B!(B,R,nseg,nus,i)
    compute_a!(a,B,Cinv,b,d,i)
    compute_v!(v,v_perp,vstar,vstar_perp,w,w_perp,a,nseg,nus,i)
    compute_xi(ξ,v,dudt,nseg)

    _weight = @view weight[1,:]

    for iseg=1:nseg
      _dgdu = @view dgdu[:,:,iseg]
      _v = @view v[:,:,iseg]
      res[i] += sum((_v.*_dgdu)*_weight)/((nstep-1)*nseg)
      res[i] += ξ[iseg,end]*(gavg-gsave[end,iseg])/(dtsave*(nstep-1)*nseg)
    end
  end

  return res
end

check_for_g(sensealg::NILSS,g) = (g===nothing && error("To use NILSS, g must be passed as a kwarg to `NILSS(g=g)`."))


text:
using Pkg
current_path = @__DIR__
Pkg.activate(current_path * "/../.");
using HarmonicBalance
using Test

files = [
    "powers.jl",
    "harmonics.jl",
    "fourier.jl",
    "load.jl",
    "parametron.jl"
    ]

for file in files
    include(file)
    printstyled(file * ":    OK\n"; color = :green)
end

printstyled("\nALL TESTS PASSED!\n"; color = :green)


text:
using Random

include("orders")
include("event")

#Initiate an empty list of portfolio/accounts

accountlist = []

mutable struct Portfolio

   """

   Class Portfolio:

   Create a portfolio object that is being used to pass order and to follow up on the equity when a strategy is being

   traded.

       Parameters

       ----------

       number: {int64}, optional, default None

           The account number that is being used to identify the account/portfolio.

       equity: {float64}, optional, default 10,000

           The amount of equity in the portfolio.

       Attributes

       ----------

       number: int64

           The account number that is being used to identify the account/portfolio.

       equity: float64

           The amount of equity in the portfolio.

       orderslist: array

           The array that will contain all the "Entry" orders created for that portfolio

   Example:

      Portfolio()

      """
   number::Int64
   equity::Float64
   orderslist::Array
   Portfolio(number = nothing, equity = 10000) = begin
      if number == nothing
         number = rand(Int, 1)
      else
         number = number
      end
      equity = equity
      orderslist = []
      newportfolio = new(number, equity, orderslist)
          push!(accountlist, newportfolio)
      end
end

function findentryorder(portfolio, assetname, strategy)

      """ Function for finding an "Entry" order in the portfolio orders list when a fill event is in the event queue. It

      returns the orders so that it can be deleted from the last as soon as an "Exit" order has been created """

   for i in portfolio.orderslist
      if i.assetname == assetname && i.strategy == strategy && i.type == "Entry"
            return i
            break
      else
         continue
      end
   end
end

function eventorder(portfolio, strategy)

      """ Verify the events queue for order or fill events, if one is found the portfolio creates a new Order, either

      an "Entry" one or an "Exit" one """

   for i in event.eventsqueue
      if i.type == "orderevent"
         if i.name == "BUY"
            orders.Order(i.asset, portfolio, "Entry", "Long", i.quantity * portfolio.equity, strategy)
         else:
            orders.Order(i.asset, portfolio, "Entry", "Short", i.quantity * portfolio.equity, strategy)
         end
         remove_event(i)
      elseif i.type == "fillevent"
         if i.name == "BUY"
            orders.Order(i.asset, portfolio, "Exit", "Long", i.quantity * portfolio.equity, strategy)
         else
            orders.Order(i.asset, portfolio, "Exit", "Short", i.quantity * portfolio.equity, strategy)
         end
         remove_event(i)
      end
   end
end

function gainloss(portfolio, amount)
      """ Add the profit or the loss to the portfolio equity """
   portfolio.equity = portfolio.equity + amount
end

max_stars_repo_name: SRelieved/JuliaBacktesting

max_issues_repo_name: MartinGenesse/JuliaBacktesting

max_forks_repo_name: MartinGenesse/JuliaBacktesting


text:
function parseInput(file)
    allNumbers::Vector{Int} = []
    open(file) do file
        for line in eachline(file)
            push!(allNumbers, parse(Int,line))
        end
    end
    return allNumbers
end

function testNumberAtIndex(allNumbers, numberIndex, preambleLength)
    result = 0
    numberToTest = allNumbers[numberIndex]
    lastIx = numberIndex-1
    startIx = numberIndex-preambleLength
    for i in startIx:1:lastIx
        for j in startIx:1:lastIx
            a = allNumbers[i]
            b = allNumbers[j]
            # println("a+b=", a+b)
            if a != b && a + b == numberToTest
                result += 1
            end
        end
    end
    return result > 0
end

function findWeaknessA(file, preambleLength, allNumbers)
    for n in preambleLength+1:length(allNumbers)
        # println(allNumbers[n])
        if !testNumberAtIndex(allNumbers,n,preambleLength)
            println("Found weaknes in ", allNumbers[n], " at index ", n)
            return(n,allNumbers[n])
        end
    end
end

function findWeaknessB(file, preambleLength)
    allNumbers = parseInput(file)
    ixA, numberToTest = findWeaknessA(file, preambleLength, allNumbers)
    contiguousNumbers = []
    for i in 1:ixA-1
        for j in i+1:ixA-1
            view = @view allNumbers[i:j]
            sum = reduce(+, view)
            if sum == numberToTest
                println("SUCCESS! ", view, "=", sum )
                return extremeSum(view)
            end
        end
    end
end

function extremeSum(array)
    min = findmin(array)[1]
    max = findmax(array)[1]
    return min+max
end

println("ExtremeSum=",findWeaknessB("data.txt", 25))


text:
counts = Dict{AbstractString, UInt64}()

for line in eachline(STDIN)
  for word in matchall(r"[a-z\']+", lowercase(line))
    counts[word] = get(counts, word, 0) + 1
  end
end
for (word, count) in sort(collect(counts))
  println("$word $count")
end