text:
using BenchmarkTools, SNOPT7, Plots
using FileIO, MeshIO, GeometryTypes, CoordinateTransformations, MeshCat

T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    iterations=300)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=40,
    cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-4,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    R_inf=1.0e-8,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    opts=Dict(:max_iter=>10000),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol,
    opts=Dict(:Iterations_limit=>500000,
        :Major_iterations_limit=>1000))

# ALTRO w/ Newton
prob_altro = copy(Problems.quadrotor)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
max_violation_direct(p1)
X1 = to_array(p1.X)
plot(X1[1:3,:]',title="Quadrotor position (ALTRO)")
plot(p1.U,title="Quadrotor control (ALTRO)")

# DIRCOL w/ Ipopt
prob_ipopt = copy(Problems.quadrotor)
rollout!(prob_ipopt)
prob_ipopt = update_problem(prob_ipopt,model=Dynamics.quadrotor_euler) # get continuous time model
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
max_violation_direct(p2)
X2 = to_array(p2.X)
plot(X2[1:3,:]',title="Quadrotor position (Ipopt)")
plot(p2.U,title="Quadrotor control (Ipopt)")

# DIRCOL w/ SNOPT
prob_snopt = copy(Problems.quadrotor)
rollout!(prob_snopt)
prob_snopt = update_problem(prob_snopt,model=Dynamics.quadrotor_euler) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
max_violation_direct(p3)
X3 = to_array(p3.X)
plot(X3[1:3,:]',title="Quadrotor position (SNOPT)")
plot(p3.U,title="Quadrotor control (SNOPT)")

# Visualization
function animate_quadrotor_line(prob)
    vis = Visualizer()
    open(vis)

    traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
    urdf_folder = joinpath(traj_folder, "dynamics","urdf")
    obj = joinpath(urdf_folder, "quadrotor.obj")

    quad_scaling = 0.7
    robot_obj = FileIO.load(obj)
    robot_obj.vertices .= robot_obj.vertices .* quad_scaling

    sphere_small = HyperSphere(Point3f0(0), convert(Float32,0.25)) # trajectory points
    sphere_med = HyperSphere(Point3f0(0), convert(Float32,0.5));
    sphere_quad = HyperSphere(Point3f0(0), convert(Float32,2.0));

    obstacles = vis["obs"]
    traj = vis["traj"]
    robot = vis["robot"]
    setobject!(vis["robot"]["quad"],robot_obj,MeshPhongMaterial(color=RGBA(0, 0, 0, 1.0)));
    # setobject!(vis["robot"]["ball"],sphere_medium,MeshPhongMaterial(color=RGBA(0, 0, 0, 0.5)));

    settransform!(vis["/Cameras/default"], compose(Translation(0., 72., 60.),LinearMap(RotX(pi/7.5)*RotZ(pi/2))))
    traj = vis["traj"]

    for i = 1:prob.N
        setobject!(vis["traj"]["t$i"],sphere_small,MeshPhongMaterial(color=RGBA(0, 1, 0, 1.0)))
        settransform!(vis["traj"]["t$i"], Translation(prob.X[i][1], prob.X[i][2], prob.X[i][3]))
    end

    anim = MeshCat.Animation(24)
    for i = 1:prob.N
        MeshCat.atframe(anim,vis,i) do frame
            settransform!(frame["robot"], compose(Translation(prob.X[i][1:3]...),LinearMap(Quat(prob.X[i][4:7]...))))
        end
    end
    MeshCat.setanimation!(vis,anim)
end

animate_quadrotor_line(p1)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
import DataStructures: DefaultDict

###############################
# @prp
###############################

struct __IsaPrp end
struct __IsnotaPrp end
@inline __isprp(x) = __IsnotaPrp()
@inline __asprp(f) = __asprp(__isprp(f), f)
@inline __asprp(::__IsaPrp, f) = o -> f(o)
@inline __asprp(::__IsnotaPrp, f) = o -> ((a...;ka...) -> f(o, a...; ka...))

export @prp

macro prp(ex)
    # f
    if isa(ex, Symbol)
        f = ex
        esc(:(@inline Wild.__isprp(::typeof($f)) = Wild.__IsaPrp()))
    # function f(x) 2x end
    elseif ex.head == :function
        f = ex.args[1].args[1]
        ff = :(@inline Wild.__isprp(::typeof($f)) = Wild.__IsaPrp())
        esc(Expr(:toplevel, ex, ff))
    elseif ex.head == :(=)
        # f(x) = 2x
        if isa(ex.args[1], Expr) && ex.args[1].head == :call
            f = ex.args[1].args[1]
            ff = :(@inline Wild.__isprp(::typeof($f)) = Wild.__IsaPrp())
            esc(Expr(:toplevel, ex, ff))
        # f = x -> 2x
        else
            f = ex.args[1]
            ff = :(@inline Wild.__isprp(::typeof($f)) = Wild.__IsaPrp())
            esc(Expr(:toplevel, ex, ff))
        end
    else
        f = ex
        esc(:(@inline Wild.__isprp(::typeof($f)) = Wild.__IsaPrp()))
    end
end

@prp Base.methods
@prp Base.inv
@prp Base.real
@prp Base.imag
@prp Base.reim
@prp Base.abs
@prp Base.abs2
@prp Base.conj
@prp Base.angle
@prp Base.cis
@prp Base.extrema
@prp Base.minimum
@prp Base.maximum
@prp Base.one
@prp Base.zero
@prp Base.keys
@prp Base.values
@prp Base.typeof
@prp Base.println
@prp Base.length
@prp Base.adjoint
@prp Base.transpose

###############################
# prpdict
###############################
#=
mutable struct __Tree
    t::Type
    p::Union{Nothing, __Tree}
    c::Vector{__Tree}
end
__Tree(t::Type) = __Tree(t, nothing, __Tree[])
__Tree(t::Type, p) = __Tree(t, p, __Tree[])
__Tree(t::Type, p, c::Type) = __Tree(t, p, __Tree[__Tree(c, t)])

@inline __Tree_is_child(t1::__Tree, t2::__Tree) = t1.t <: t2.t
@inline __Tree_is_child(t1::__Tree, t2::Type) = t1.t <: t2
@inline __Tree_is_child(t1::Type, t2::__Tree) = t1 <: t2.t
@inline __Tree_is_child(t1::Type, t2::Type) = t1 <: t2

@inline __Tree_is_parent(t1::__Tree, t2::__Tree) = t2.t <: t1.t
@inline __Tree_is_parent(t1::__Tree, t2::Type) = t2 <: t1.t
@inline __Tree_is_parent(t1::Type, t2::__Tree) = t2.t <: t1
@inline __Tree_is_parent(t1::Type, t2::Type) = t2 <: t1

__Tree_insert(tr::__Tree, t::Type) =
begin
    __Tree_is_parent(t, tr) && (return __Tree(t, nothing, tr))
    !__Tree_is_child(t, tr) && error("something wrong")

    for (i,x) in enumerate(tr.c)
        if __Tree_is_child(t, x)
            __Tree_insert(x, t)
            return tr
        elseif __Tree_is_parent(t, x)
            oldtree = tr.c[i]
            newtree = __Tree(t, tr, [oldtree])
            oldtree.p = newtree
            tr.c[i] = newtree
            return tr
        end
    end
    push!(tr.c, __Tree(t, tr))
    tr
end

__Tree_path(tr::__Tree, t::Type, path=Type[]) =
begin
    push!(path, tr.t)
    tr.t == t && (return path)
    for tt in tr.c
        t <: tt.t && (return __Tree_path(tt, t, path))
    end
    path
end

getprp_dct_type_tree = __Tree(Any)
getprp_dct_path_dct = Dict{Type, Vector{Type}}()

__prpdct_type_list = [Any,
                      Number,
                      Integer,
                      Complex,
                      Tuple,
                      NamedTuple,
                      AbstractArray,
                      AbstractVector,
                      AbstractRange,
                      AbstractString,
                      Function,
                      Base.Generator,
                      Iterators.ProductIterator,
                      ]

for T in __prpdct_type_list[2:end]
    __Tree_insert(getprp_dct_type_tree, T)
end
for T in __prpdct_type_list
    getprp_dct_path_dct[T] = reverse(__Tree_path(getprp_dct_type_tree, T))
end

__PrpDct = Dict{Symbol, Function}
struct __getprp_dct <: Function _dct::Dict{Type, __PrpDct} end
Base.getindex(d::__getprp_dct, T::Type) = d._dct[T]
Base.setindex!(d::__getprp_dct, x::__PrpDct, T::Type) = Base.setindex!(d._dct, x, T)

getprp_dct = __getprp_dct(Dict{Type, __PrpDct}())

for T in __prpdct_type_list

    getprp_dct[T] = __PrpDct()

    if T == Any continue end

    Base.getproperty(o::T, atr::Symbol) =
        begin
            hasfield(typeof(o), atr) && (return Base.getfield(o, atr))
            for t in getprp_dct_path_dct[T]
                haskey(getprp_dct[t], atr) && (return getprp_dct[t][atr](o))
            end
            __asprp(Base.eval(Base.Main, atr))(o)
        end

    Base.hasproperty(o::T, atr::Symbol) =
        begin
            hasfield(typeof(o), atr) && (return true)
            for t in getprp_dct_path_dct[T]
                haskey(getprp_dct[t], atr) && (return true)
            end
            false
        end
    Base.propertynames(o::T, private::Bool=false) =
        tuple(fieldnames(typeof(o))...,
              vcat([collect(keys(getprp_dct[t]))
                    for t in getprp_dct_path_dct[T]]...)...)
end

#####################################
# for adding a new type to getprp_dct
#####################################

insert_type_to_getprp_dct(T::Type) =
begin
    __Tree_insert(getprp_dct_type_tree, T)
    getprp_dct_path_dct[T] = reverse(__Tree_path(getprp_dct_type_tree, T))
    getprp_dct[T] = __PrpDct()
end
=#

#=
Wild.insert_type_to_getprp_dct(T)

Base.getproperty(o::T, atr::Symbol) =
    begin
        hasfield(typeof(o), atr) && (return Base.getfield(o, atr))
        for t in Wild.getprp_dct_path_dct[T]
            haskey(Wild.getprp_dct[t], atr) && (return Wild.getprp_dct[t][atr](o))
        end
        __asprp(Base.eval(Base.Main, atr))(o)
    end

Base.hasproperty(o::T, atr::Symbol) =
    begin
        hasfield(typeof(o), atr) && (return true)
        for t in Wild.getprp_dct_path_dct[T]
            haskey(Wild.getprp_dct[t], atr) && (return true)
        end
        false
    end
Base.propertynames(o::T, private=false) =
    tuple(fieldnames(typeof(o))...,
          vcat([collect(keys(Wild.getprp_dct[t]))
                for t in Wild.getprp_dct_path_dct[T]]...)...)
=#

ext: jl


text:
#GracePlot: I/O facilities
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const MIMEpng = MIME"image/png"
#const MIMEsvg = MIME"image/svg+xml"
#const MIMEeps = MIME"image/eps"
#const MIMEpdf = MIME"application/pdf"

#==Helper functions
===============================================================================#
#Opens a file for read, after it is done being written:
#(Writing files with grace is a non-blocking operation.)
#timeout in seconds
function openreadafterexport(p::Plot, filepath::String; timeout=15)
	poll_interval = .1 #WANTCONST sec
	timeout = timeout*1_000_000_000
	flushpipe(p) #Make sure last write operation is registered with Grace.
	tstart = time_ns()
	local io

	while true
		try
#			filepath = "nope"
			io = open(filepath, "r")
			break
		catch e
			sleep(poll_interval)
			if time_ns()-tstart > timeout
				rethrow(e)
			end
		end
	end

	szlast = 0
	while true
		sleep(poll_interval)
		sz = filesize(io)
		if sz == 0
			if time_ns()-tstart > timeout
				close(io)
				throw("$filepath access timed out.")
			end
		elseif sz == szlast
			break
		end
		szlast = sz
	end
#	@show filesize(io)/1e6
	return io
end
#-------------------------------------------------------------------------------
function exportplot(p::Plot, filefmt::String, filepath::String)
	sendcmd(p, "HARDCOPY DEVICE \"$filefmt\"")
	sendcmd(p, "PRINT TO \"$filepath\"")
	sendcmd(p, "PRINT")
	flushpipe(p)
	#Sadly, Grace does not wait unil file is saved to return from "PRINT"
end

#==_write interface
===============================================================================#

#Write plot to Grace .agr file:
#-------------------------------------------------------------------------------
function _write(path::String, p::Plot)
	_ensure(!occursin(path, "\""), ArgumentError("File path contains '\"'."))
	sendcmd(p, "SAVEALL \"$path\"")
	flushpipe(p)
end

#==write_FILEFMT interface
===============================================================================#

#Write to PNG:
#-------------------------------------------------------------------------------
function _write_png(path::String, p::Plot, dpi::Int)
	w = round(Int, val(TPoint(p.canvas.width)));
	h = round(Int, val(TPoint(p.canvas.height)));
	sendcmd(p, "DEVICE \"PNG\" DPI $dpi") #Must set before PAGE SIZE
	sendcmd(p, "DEVICE \"PNG\" PAGE SIZE $w, $h")
	exportplot(p, "PNG", path)
end
_write_png(path::String, p::Plot, ::Nothing) =
	_write_png(path, p, p.dpi) #Use dpi setting in plot

#User-level wrapper function:
write_png(path::String, p::Plot; dpi::Union{Int,Nothing}=nothing) =
	_write_png(path, p, dpi)

#Write to EPS:
#-------------------------------------------------------------------------------
function write_eps(path::String, p::Plot)
	sendcmd(p, "DEVICE \"EPS\" OP \"bbox:page\"")
	exportplot(p, "EPS", path)
end

#Write to SVG.  (Fix Grace output according W3C 1999 format):
#TODO: Make more robust... use more try/catch.
#NOTE: Replace xml (svg) statements using Julia v3 compatibility.
#-------------------------------------------------------------------------------
function write_svg(path::String, p::Plot)
	tmpfilepath = "$(tempname())_export.svg"
	#Export to svg, using the native Grace format:
	exportplot(p, "SVG", tmpfilepath)
	#NOTE: Is this the 2001 format?  Most programs do not appear to read it.

	src = openreadafterexport(p, tmpfilepath)
	filedat = read(src, String)
	close(src)

	#Remove <!DOCTYPE svg ...> statement:
	pat = r"^[ \t]*<!DOCTYPE svG.*$"mi
	filedat = replace(filedat, pat => "")

	#Modify <svg ...> statement:
	pat = r"(^[ \t]*<svg) xml:space=\"preserve\" (.*$)"mi
	captures = match(pat, filedat).captures
	cap1 = captures[1]; cap2 = captures[2]
	filedat = replace(filedat, pat => "$cap1 xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" $cap2")

	dest = open(path, "w")
	write(dest, filedat)
	close(dest)

	rm(tmpfilepath)
end

#==MIME support
===============================================================================#
function Base.show(io::IO, ::MIMEpng, p::Plot; dpi::Union{Int,Nothing}=nothing)
	tmpfile = "$(tempname())_export.png"
	_write_png(tmpfile, p, dpi)
	flushpipe(p)
	src = openreadafterexport(p, tmpfile)
	data = read(src, String)
	write(io, data)
	close(src)
	rm(tmpfile)
end

#=
#==FileIO2 support:
===============================================================================#
#	-Deactivated to avoid un-necessary dependency (at this point).
#	-TODO: Re-activate once FileIO2 is registered?

using FileIO2

#Declare file formats (Not exported):
abstract type GraceFmt <: FileIO2.DataFormat end #Grace plot file
abstract type ParamFmt <: FileIO2.DataFormat end #Grace "parameter" (.par template) format

#Could also define shorthand:
#FileIO2.File(::FileIO2.Shorthand{:grace}, path::String) = File{GraceFmt}(path)
#FileIO2.File(::FileIO2.Shorthand{:graceparam}, path::String) = File{ParamFmt}(path)

#Typed _write interface:
#_write(file::File{ParamFmt}, p::Plot) = #TODO: Is there a way to write .par files?
_write(file::File{GraceFmt}, p::Plot) = _write(path::String, p::Plot)
_write(file::File{PNGFmt}, p::Plot) = write_png(path::String, p::Plot)
_write(file::File{SVGFmt}, p::Plot) = write_svg(path::String, p::Plot)
_write(file::File{EPSFmt}, p::Plot) = write_eps(path::String, p::Plot)
=#

#Last line


text:
# TODO: column generation for all these functions.

function _master_mmf_problem(m::Model, rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction)
    rm = basic_routing_model_unitary(m, rd)
    capacity_constraints(rm, rd.traffic_matrix)

    @variable(m, τ)
    rm.tau = τ

    return rm
end

function master_mmf_problem(m::Model, rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::MinMaxFair)
    rm = _master_mmf_problem(m, rd, edge_obj)

    @objective(m, Min, rm.tau)
    @constraint(m, mmf[e in edges(rd)], objective_edge_expression(rm, edge_obj, e) <= rm.tau)
    rm.constraints_mmf = mmf

    return rm
end

function master_mmf_problem(m::Model, rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::MaxMinFair)
    rm = _master_mmf_problem(m, rd, edge_obj)

    @objective(m, Max, rm.tau)
    @constraint(m, mmf[e in edges(rd)], objective_edge_expression(rm, edge_obj, e) >= rm.tau)
    rm.constraints_mmf = mmf

    return rm
end

function compute_routing(rd::RoutingData, edge_obj::EdgeWiseObjectiveFunction, agg_obj::Union{MinMaxFair, MaxMinFair}, type::FormulationType, cg::Val, algo::Automatic, unc::NoUncertaintyHandling, uncparams::NoUncertainty)
    # Implementation of the water-filling algorithm for MMF routings (aka MFMF).
    # https://ica1www.epfl.ch/PS_files/LEB3132.pdf
    # https://www.utc.fr/~dnace/recherche/Publication/Inoc03Nace.pdf

    start = time_ns()

    # Create the master problem.
    m = _create_model(rd)
    rm = master_mmf_problem(m, rd, edge_obj, agg_obj)
    time_create_master_model_ms = (time_ns() - start) / 1_000_000.

    # Memorise the evolution of the algorithm.
    n_iter = 0
    n_new_paths = 0
    routings = Routing[]
    objectives = Float64[]
    times_ms = Float64[]
    end_status = MOI.OPTIMAL

    # Iteratively solve it by removing τ constraints and fixing the values of
    # the corresponding objectives.
    edges_to_do = Set(edges(rd))
    fixed_objectives = Dict{Edge{Int}, Float64}()

    while length(edges_to_do) > 0
        start = time_ns()

        # Optimise for this iteration, possibly using column generation.
        cg_it = 0
        while true
            # Actual optimisation
            optimize!(m)
            status = termination_status(m)
            _export_lp_if_allowed(rd, m, "lp_$(n_iter)_$(cg_it)")

            # Handle abnormal cases.
            if status != MOI.OPTIMAL
                # OPTIMAL status is shown with the objective value.
                rd.logmessage("[$(n_iter)] Status: $(status).")

                # Debug infeasibility if need be.
                _debug_infeasibility_mmf(rd, status, n_iter, edge_obj, agg_obj, type, cg, algo, unc, uncparams)

                # Report other problematic codes.
                if status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED] && n_iter > 0
                    end_status = MOI.ALMOST_OPTIMAL
                    break
                end
            end

            # Quit if no column generation process allowed.
            if cg == Val(false)
                break
            end

            # Perform some column generation if need be.
            @assert typeof(type) <: PathFormulation

            # First, the pricing.
            dual_value_convexity = dual.(rm.constraints_convexity) # Demand -> dual value
            dual_values_capacity = Dict(e => dual(rm.constraints_capacity[e]) for e in edges(rd)) # Edge -> dual value
            dual_values_mmf = Dict(e => dual(rm.constraints_mmf[e]) for e in edges(rd)) # Edge -> dual value

            new_paths = _mmf_solve_pricing_problem(rd, dual_values_capacity, dual_values_mmf, dual_value_convexity)
            n_new_paths += length(new_paths)

            # Add the new columns to the formulation.
            for (d, path) in new_paths
                # Add the new columns to the data object and retrieve a path ID.
                path_id = add_path!(rm, d, path)

                # Get your hands dirty in the formulation.
                # - Create a new variable for this path, add it into the
                #   convexity and capacity constraints.
                add_routing_var!(rm, demand, path_id, constraint_capacity=true)

                # TODO: normalised way of having traffic matrices in rd? Are graph attributes forbidden/unused everywhere? What about passing them as arguments, as helpers in knowncapacities?
                dd = rd.traffic_matrix[d]
                for e in path
                    # - MMF constraint.
                    ce = capacity(rd, e)
                    set_normalized_coefficient(rm.constraints_mmf[e], rm.routing[d, path_id], dd / ce)
                end
            end

            rd.logmessage("[$(n_iter) | $(cg_it)] Added $(length(new_paths)) new paths.")
            cg_it += 1

            if length(new_paths) == 0
                break
            end
        end

        # Report the values for this iteration.
        master_τ = objective_value(m)
        rd.logmessage("[$(n_iter)] Status: $(status). Value: $(master_τ)")
        push!(routings, Routing(rd, value.(rm.routing)))
        push!(objectives, master_τ)

        # Find edges to fix.
        new_edges = _mmf_find_edges_to_fix(rm, master_τ, edges_to_do, agg_obj)
        rd.logmessage("  Number of variables to fix: $(length(new_edges)).")

        # Debugging when no edge met the conditions.
        if length(new_edges) == 0
            if n_iter == 0
                # No variable can be fixed right now, there must be a problem!
                @warn "No variable can be fixed at the first iteration, " *
                      "the MMF process stalls."
            else
                # Every link seems to be saturated, end now.
                rd.logmessage("[$(n_iter)] All links are saturated.")
            end
        end

        # Add the new constraints for the next iteration.
        for (e, value) in new_edges
            # Maintain the data structures.
            fixed_objectives[e] = value
            pop!(edges_to_do, e)

            # Modify the master to indicate that the value is now known.
            _mmf_fix_edge(rm, e, value, agg_obj)
        end

        rd.logmessage("  Edges not yet saturated: $(length(edges_to_do)). Total edges saturated: $(length(fixed_objectives))")

        # Log a few things at the end of the iteration. Stop if no new edge was
        # found.
        n_iter += 1
        push!(times_ms, (time_ns() - start) / 1_000_000.)

        if length(new_edges) == 0
            # The user is already warned.
            break
        end
    end

    # Not necessary to optimise one last time: the remaining edges were fixed
    # to the current value in the master (i.e. the last iteration just added
    # redundant constraints).

    return RoutingSolution(rd, result=end_status,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_create_master_model_ms=time_create_master_model_ms,
                           time_solve_ms=times_ms,
                           n_columns=n_new_paths,
                           objectives=objectives,
                           matrices=rd.traffic_matrix,
                           routings=routings,
                           master_model=rm)
end

function _mmf_solve_pricing_problem(rd::RoutingData, dual_values_capacity, dual_values_mmf, dual_value_convexity)
    paths = Dict{Edge{Int}, Vector{Edge}}()

    # Determine the weight matrix to use for the computations:
    #     \sum_{e\in p} dualcapa_{e} + dualmmf_{e} / capa_e
    weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
    for e in edges(rd)
        value = dual_values_capacity[e] + dual_values_mmf[e] / capacity(rd, e)
        if value <= - CPLEX_REDUCED_COST_TOLERANCE
            weight_matrix[src(e), dst(e)] = value
        end
    end
    weight_matrix .+= maximum(abs.(weight_matrix))

    # Check that at least some weights are nonzero.
    if nnz(weight_matrix) == 0
        # No need for the potentially expensive _check_weight_matrix,
        # as no value in weight_matrix can be below the tolerance by
        # construction (no addition, unlike iterative version of oblivious).
        return paths
    end

    # Perform the pricing for each demand.
    for demand in demands(rd)
        rd.logmessage("Pricing for $demand")

        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
            continue
        end
        reduced_cost = - state.dists[dst(demand)] * rd.traffic_matrix[demand] - dual_value_convexity[demand]
        if reduced_cost >= - CPLEX_REDUCED_COST_TOLERANCE
            continue
        end
        rd.logmessage(reduced_cost)

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
            rd.logmessage("New path $(path) is already in the formulation!")
        else
            paths[demand] = path
        end
    end

    return paths
end

function _mmf_find_edges_to_fix(rm::RoutingModel, master_τ::Float64, edges_to_do, agg_obj)
    # Check which constraints are satisfied at equality (i.e. nonzero dual
    # value, as the constraint is either ≥ or ≤).
    new_edges = Dict{Edge{Int}, Float64}()

    if has_duals(rm.model)
        for e in edges_to_do
            if - agg_obj.ε <= dual(rm.constraints_mmf[e]) <= agg_obj.ε
                new_edges[e] = master_τ
            end
        end
    end

    if length(new_edges) == 0
        # If no edge goes through the dual variable test, check the actual
        # values (this may be due to the way the objective function is
        # modelled, like Fortz-Thorup, or if the solver doesn't give duals).
        for e in edges_to_do
            if - agg_obj.ε <= value(rm.constraints_mmf[e]) <= agg_obj.ε
                new_edges[e] = master_τ
            end
        end
    end

    return new_edges
end

function _mmf_fix_edge(rm::RoutingModel, e::Edge{Int}, value::Float64, rhs::Float64)
    rm.data.logmessage("  Fixing $e to $value ≈ $rhs (relaxed value)")
    set_normalized_coefficient(rm.constraints_mmf[e], rm.tau, 0)
    set_normalized_rhs(rm.constraints_mmf[e], rhs)
end

_mmf_fix_edge(rm::RoutingModel, e::Edge{Int}, value::Float64, agg_obj::MinMaxFair) =
    _mmf_fix_edge(rm, e, value, value * (1 + agg_obj.ε))

_mmf_fix_edge(rm::RoutingModel, e::Edge{Int}, value::Float64, agg_obj::MaxMinFair) =
    _mmf_fix_edge(rm, e, value, value * (1 - agg_obj.ε))

# Debug infeasiblity. Implement the main logic here. Throw an error in case of
# infeasibility.
function _debug_infeasibility_mmf(rd::RoutingData, status::MOI.TerminationStatusCode, n_iter::Int,
                                  edge_obj::EdgeWiseObjectiveFunction, agg_obj::Union{MinMaxFair, MaxMinFair},
                                  type::FormulationType, cg::Val, algo::Automatic,
                                  unc::NoUncertaintyHandling, uncparams::NoUncertainty)
    if status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED] && n_iter == 0
        rd.logmessage("[$(n_iter)] This iteration of MMF yielded an infeasible optimisation program. ")
        d = _debug_infeasibility_mmf(rd, edge_obj, agg_obj, type, cg, algo, unc, uncparams)
        if d !== nothing
            rd.logmessage("[$(n_iter)] The source of infeasibility could be automatically detected: ")
            rd.logmessage("[$(n_iter)] $d")
        else
            rd.logmessage("[$(n_iter)] The source of infeasibility could not be automatically detected.")
        end

        # Exit with an error.
        _export_lp_if_failed(rd, status, m, "error_$(n_iter)", "Current problem could not be solved!")
    end
end

_debug_infeasibility_mmf(::RoutingData, ::EdgeWiseObjectiveFunction,
                         ::MinMaxFair, ::FormulationType, ::Val, ::Automatic,
                         ::NoUncertaintyHandling, ::NoUncertainty) = nothing

function _debug_infeasibility_mmf(rd::RoutingData, ::KleinrockLoad,
                                  ::MinMaxFair, type::FormulationType,
                                  ::Val{false}, ::Automatic,
                                  ::NoUncertaintyHandling, ::NoUncertainty)
    # Main source of infeasibility: impossible to have links with a load less
    # than 100%.
    # Technique: minimise the maximum load (i.e. not Kleinrock), check if it
    # is 100%.
    r = compute_routing(rd, Load(), MinimumMaximum(), type, Val(false),
                        Automatic(), NoUncertaintyHandling(), NoUncertainty())
    if compute_max_load(r) >= 0.99999
        return "This network cannot withstand this traffic matrix: " *
               "its maximum load is 100% or more (i.e. the most loaded link is " *
               "loaded at 100%, when minimising this maximum load). " *
               "Kleinrock function associates an infinite penalty to links " *
               "with such a high load."
    end
    return nothing
end


text:
"Load a `NULL`-terminated list of strings"
function loadstringlist(pstringlist::Ptr{Ptr{UInt8}})
    stringlist = Vector{String}()
    (pstringlist == C_NULL) && return stringlist
    i = 1
    item = unsafe_load(pstringlist, i)
    while item != C_NULL
        push!(stringlist, unsafe_string(item))
        i += 1
        item = unsafe_load(pstringlist, i)
    end
    stringlist
end

# --- Raster Metadata ---

"""
Fetch list of (non-empty) metadata domains. (Since: GDAL 1.11)

### Returns
`NULL` or a string list. Must be freed with `CSLDestroy()`
"""
_getmetadatadomainlist(obj::GDALMajorObjectH) =
    GDALGetMetadataDomainList(obj)::Ptr{Ptr{UInt8}}

getmetadatadomainlist(obj::GDALMajorObjectH) =
    loadstringlist(_getmetadatadomainlist(obj))

"""
Fetch metadata.

The returned string list is owned by the object, and may change at any time.
It is formated as a "Name=value" list with the last pointer value being `NULL`.
Use the CPLStringList functions such as `CSLFetchNameValue()` to manipulate it.

Note that relatively few formats return any metadata at this time.

### Parameters
* `domain`   the domain of interest. Use `""` or `NULL` for the default domain.

### Returns
`NULL` or a string list.
"""
_getmetadata(obj::GDALMajorObjectH, domain::Ptr{UInt8}) =
    GDALGetMetadata(obj, domain)::Ptr{Ptr{UInt8}}

getmetadata(obj::GDALMajorObjectH, domain::String = "") =
    loadstringlist(_getmetadata(obj, pointer(domain)))

"""
Set metadata.

### Parameters
* `metadata`  the metadata in name=value string list format to apply.
* `domain`    the domain of interest. Use "" or NULL for the default domain.

### Returns
`CE_None` on success, `CE_Failure` on failure and `CE_Warning` if the metadata
has been accepted, but is likely not maintained persistently by the underlying
object between sessions.
"""
_setmetadata(obj::GDALMajorObjectH,
             metadata::Ptr{Ptr{UInt8}},
             domain::Ptr{UInt8}) =
    GDALSetMetadata(obj, metadata, domain)::CPLErr

function setmetadata!(obj::GDALMajorObjectH,
                      metadata::Vector{String},
                      domain::String = "")
    result = _setmetadata(obj,
                          Ptr{Ptr{UInt8}}(pointer(metadata)),
                          pointer(domain))
    (result == CE_Failure) && error("Failed to set metadata")
    if result == CE_Warning
        warn("Metadata might not persist between sessions")
    end
    @assert result == CE_None
end

"""
Fetch single metadata item.

### Parameters
* `name`     the key for the metadata item to fetch.
* `domain`   the domain to fetch for, use `NULL` for the default domain.

### Returns
`NULL` on failure to find the key, or a pointer to an internal copy of the
value string on success.
"""
_getmetadataitem(obj::GDALMajorObjectH,
                 name::Ptr{UInt8},
                 domain::Ptr{UInt8}) =
    GDALGetMetadataItem(obj, name, domain)::Ptr{UInt8}

function getmetadataitem(obj::GDALMajorObjectH,
                         name::String,
                         domain::String)
    result = _getmetadataitem(obj, pointer(name), pointer(domain))
    result == C_NULL && error("Unable to find key $name in domain $domain")
    unsafe_string(result)
end

function getmetadataitem(obj::GDALMajorObjectH,
                         name::String)
    result = _getmetadataitem(obj, pointer(name), Ptr{UInt8}(C_NULL))
    result == C_NULL && error("Unable to find key $name")
    unsafe_string(result)
end

"""
Set single metadata item.

### Parameters
* `name`     the key for the metadata item to fetch.
* `value`    the value to assign to the key.
* `domain`   the domain to set within, use `NULL` for the default domain.

### Returns
`CE_None` on success, or an error code on failure.
"""
_setmetadataitem(obj::GDALMajorObjectH,
                 name::Ptr{UInt8},
                 value::Ptr{UInt8},
                 domain::Ptr{UInt8} = Ptr{UInt8}(C_NULL)) =
    GDALSetMetadataItem(obj, name, value, domain)::CPLErr

function setmetadataitem!(obj::GDALMajorObjectH,
                          name::String,
                          value::String,
                          domain::String)
    result = _setmetadataitem(obj, pointer(name), pointer(value),
                              pointer(domain))
    result == CE_None || error("Unable to set metadata item")
end

function setmetadataitem!(obj::GDALMajorObjectH,
                          name::String,
                          value::String)
    result = _setmetadataitem(obj, pointer(name), pointer(value))
    result == CE_None || error("Unable to set metadata item")
end

"""
Fetch object description.

The semantics of the returned description are specific to the derived type.

* For `GDALDatasets` it is the dataset name.
* For `GDALRasterBands` it is actually a description (if supported) or `""`.

### Returns
non-null pointer to internal description string.
"""
_getdescription(obj::GDALMajorObjectH) = GDALGetDescription(obj)::Ptr{UInt8}

getdescription(obj::GDALMajorObjectH) = unsafe_string(_getdescription(obj))

"""
Set object description.

The semantics of the description are specific to the derived type.

For `GDALDatasets` it is the dataset name.
For `GDALRasterBands` it is actually a description (if supported) or `""`.

Normally application code should not set the "description" for `GDALDatasets`.
It is handled internally.
"""
_setdescription(obj::GDALMajorObjectH, desc::Ptr{UInt8}) =
    GDALSetDescription(obj, desc)

setdescription!(obj::GDALMajorObjectH, desc::String) =
    _setdescription(obj, pointer(desc))

"""
Get runtime version information.

Available pszRequest values:

* `"VERSION_NUM"`:  Returns GDAL_VERSION_NUM as a string. (ie. "1170")
Starting with GDAL 1.10, this string will be >4 characters.
* `"RELEASE_DATE"`: Returns GDAL_RELEASE_DATE as a string. (ie. "20020416")
* `"RELEASE_NAME"`: Returns the GDAL_RELEASE_NAME. (ie. "1.1.7")
* `"--version"`:    Returns one line version message suitable for use in
response to --version requests. ie. "GDAL 1.1.7, released 2002/04/16"
* `"LICENSE"`:      Returns the content of the LICENSE.TXT file from the
GDAL_DATA directory.
* `"BUILD_INFO"`:   List of NAME=VALUE pairs separated by newlines with
information on build time options.

### Parameters
`request`:  the type of version info desired, as listed above.

### Returns
an internal string containing the requested information.
"""
_versioninfo(request::Ptr{UInt8}) = GDALVersionInfo(request)::Ptr{UInt8}

versioninfo(request::String) = unsafe_string(_versioninfo(pointer(request)))

"""
Return `TRUE` if GDAL version at runtime matches nVersionMajor.nVersionMinor.

The purpose of this method is to ensure that calling code will run with the
GDAL version it is compiled for. It is primarly intented for external plugins.

### Parameters
* `major`   Major version to be tested against
* `minor`   Minor version to be tested against
* `name`    If not NULL, in case of version mismatch, the method will issue a
failure mentionning the name of the calling component.

### Returns
`TRUE` if GDAL library version at runtime matches nVersionMajor.nVersionMinor,
`FALSE` otherwise.
"""
_checkversion(major::Integer, minor::Integer, name::Ptr{UInt8}=Ptr{UInt8}(C_NULL)) =
    GDALCheckVersion(major, minor, name)::Cint

checkversion(major::Integer, minor::Integer) = Bool(_checkversion(major, minor))

checkversion(major::Integer, minor::Integer, name::String) =
    Bool(_checkversion(major, minor, pointer(name)))

"""
Set maximum cache memory.

This function sets the maximum amount of memory that GDAL is permitted to use
for GDALRasterBlock caching. The unit of the value is bytes.

The maximum value is 2GB, due to the use of a signed 32 bit integer. Use
`GDALSetCacheMax64()` to be able to set a higher value.

### Parameters
`nbytes`     the maximum number of bytes for caching.
"""
_setcachemax(nbytes::Integer) = GDALSetCacheMax(nbytes)

"""
Get maximum cache memory.

Gets the maximum amount of memory available to the `GDALRasterBlock` caching
system for caching GDAL read/write imagery.

The first type this function is called, it will read the `GDAL_CACHEMAX`
configuation option to initialize the maximum cache memory.

This function cannot return a value higher than 2 GB. Use `GDALGetCacheMax64()`
to get a non-truncated value.

### Returns
maximum in bytes.
"""
_getcachemax() = GDALGetCacheMax()::Cint

"""
Get cache memory used.

### Returns
number of bytes of memory currently in use by `GDALRasterBlock` memory caching.
"""
_getcacheused() = GDALGetCacheUsed()::Cint

"""
Try to flush one cached raster block.

This function will search the first unlocked raster block and will flush it to
release the associated memory.

### Returns
`TRUE` if one block was flushed, `FALSE` if there are no cached blocks or if
they are currently locked.
"""
_flushcacheblock() = GDALFlushCacheBlock()::Cint


text:
using Dates
using ArgParse
using Printf
using Random
using GaussianProcesses

# before_all = Dates.now()

include("../../tools/split.jl")
include("../../btg.jl")
include("../../covariatefun.jl")

###
###  Manually enter parsed args entries so we can run this file as a script
###
##
parsed_args = Dict()
push!(parsed_args, "randseed" => 1234)
push!(parsed_args, "lmax" => 0.5)
push!(parsed_args, "lmin" => 0.2)
#push!(parsed_args, "noiselevel" => 1)
push!(parsed_args, "N" => 100)
push!(parsed_args, "p_train" => 0.5)
push!(parsed_args, "p" => 1) #dimension of covariates
push!(parsed_args, "N_plot"=>200) #number pts plotted on graph
push!(parsed_args, "upper_bound_estimate" => 10) #crude estimate for upper bound of function for box-constrained opt
randseed = parsed_args["randseed"]; rng = MersenneTwister(randseed)
push!(parsed_args, "log_scale" => false)

Random.seed!(randseed);

#fun(x) = .5sin(10x) + 1+ 1/((x-5)/4)^2
fun(x) = sin(10x) .+ 2 .+ 0.5*rand(1)[1] #noisy settin
# training set
N = parsed_args["N"]
xs = range(0.001, stop = 4.9999, length = N)
#xs = vcat(range(.001, stop=4.8, length=N), range(4.81, stop = 4.99, length = Int64(round(N/4))))
#xs = vcat(xs, range(4.991, stop = 4.9999, length = Int64(round(N/8))))
inds = collect(1:1:length(xs))
ind_train, x_train, ind_test, x_test  = sample(inds, parsed_args["p_train"])
x_test = reshape(x_test, length(x_test), 1) #reshape
x_train = reshape(x_train, length(x_train), 1) #reshape

n_train = length(ind_train)
n_test = length(ind_test)
ys = fun.(xs)

#data = sort((rand(Float64, (n_train, )) .- 0.5)*2*pi)
y_train = ys[ind_train]
Fx = covariate_fun(x_train, parsed_args["p"])
trainingData0 = trainingData(reshape(x_train, length(x_train), 1), Fx, y_train)
d = getDimension(trainingData0); n = getNumPts(trainingData0); p = getCovDimension(trainingData0)

Fx_test = covariate_fun(reshape(x_test, length(x_test), 1), parsed_args["p"])
y_test = ys[ind_test]
elapsedmin = 0

#parameter setting
# myquadtype = parsed_args["sparse"] ? ["SparseCarlo", "SparseCarlo"] : ["QuasiMonteCarlo", "QuasiMonteCarlo"]
myquadtype = ["Gaussian", "Gaussian"]
rangeλ = [0.5 1.5]
#rangeλ = reshape([1.0], 1, 1)
lmin = parsed_args["lmin"]
lmax = parsed_args["lmax"]
#rangeθ = [1/lmax^2 1/lmin^2]
#rangeθ = reshape([140.0], 1, 1)
#rangeθ = reshape([145.0], 1, 1)
rangeθ = reshape([14.0], 1, 1)
@info "rangeθ, rangel:" rangeθ, [lmin lmax]
# rangeθ = [0.111 25]
# build btg model

function load_synthetic_btg()
    return btg(trainingData0, rangeθ, rangeλ; priorθ = inverseUniform(rangeθ), quadtype = myquadtype)
end

function getseed()
    return randseed
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
abstract type IRISRequest <: SeisRequest end

protocol_string(::IRISRequest) = "irisws"

"Available processing filters and allowable types for them"
const IRISTimeSeries_PROCESSING_FIELDS = Dict(
    :taper => Real,
    :taper_type => AbstractString,
    :envelope => Bool,
    :lpfilter => Real,
    :hpfilter => Real,
    :bpfilter => Tuple{Number,Number},
    :demean => Bool,
    :scale => Union{Real,AbstractString},
    :divscale => Real,
    :correct => Bool,
    :freqlimits => Tuple{Number,Number,Number,Number},
    :autolimits => Union{AbstractString,Tuple{Number,Number}},
    :units => AbstractString,
    :diff => Bool,
    :int => Bool,
    :decimate => Real)

const IRISTimeSeries_OUTPUT_TYPES = ("ascii", "ascii2", "ascii", "geocsv", "geocsv.tspair",
    "geocsv.slist", "audio", "miniseed", "plot", "saca", "sacbb", "sacbl")

"""
    IRISTimeSeries

Create an data query which can be sent to a datacentre which implements the IRIS
Web Services timeseries specification.

## Available options

### Channel Options (required)
The four SCNL parameters (Station – Channel – Network – Location) are used to determine
the channel of interest, and are all required. Wildcards are not accepted.

- `network`: Seismic network name
- `station`: Station name
- `location`: Location code. Use `loc=--` for empty location codes
- `channel`: Channel Code

### Date-Range Options (required)
Time ranges must be defined by specifying a start time with either an end time or a duration
in seconds.

For example, the following two time specifications are valid, and are equivalent.

#### Examples

    IRISTimeSeries([SCNL]..., starttime=Dates.DateTime(2005, 02, 23, 12), endtime=Dates.DateTime(2005, 02, 23, 12, 05))
    IRISTimeSeries([SCNL]..., starttime=Dates.DateTime("2005-02-23T12:00:00"), duration=300)

- `starttime`: Start time
- `endtime`: End time
- `duration`: Duration of requested data, in seconds

### Time Series Processing Options
The following parameters act as filters upon the timeseries. Parameter order matters
because each filter operation is performed in the order given.

#### Examples

    IRISTimeSeries(..., demean=true, lpfilter=2.0) # will demean and then apply a low-pass filter
    IRISTimeSeries(..., lpfilter=2.0, demean=true) # will apply a low-pass filter, and then demean

- `taper`: Apply a time domain symmetric tapering function to the timeseries data. The
  width is specified as a fraction of the trace length from 0 to 0.5.
- `taper_type`: Form of taper to apply.  Supported types: HANNING (default), HAMMING, COSINE
- `envelope`: Calculate the envelope of the time series. This calculation uses a Hilbert
  transform approximated by a time domain filter
- `lpfilter`: Low-pass filter the time-series using an IIR 4th order filter, using this value
  (in Hertz) as the cutoff
- `hpfilter`: High-pass filter the time-series using an IIR 4th order filter, using this value
  (in Hertz) as the cutoff
- `bpfilter`: Band pass frequencies, in Hz. (See NOTE on variable separators1)
- `demean`: Remove mean value from data. true or false (default false).
- `scale`: Scale data samples by specified factor. When scale=AUTO scales by the stage-zero gain
- `divscale`: Scale data samples by the inverse of the specified factor
- `correct`: Apply instrument correction to convert to earth units. true or false. Uses
  either deconvolution3,4 or polynomial response correction. true or false (default false)
- `freqlimits`: Specify an envelope for a spectrum taper for deconvolution. Frequencies are
  specified in Hertz. This cosine taper scales the spectrum from 0 to 1 between f1 and f2 and
  from 1 to 0 between f3 and f4. Can only be used with the correct option. Cannot be used in
  combination with the autolimits option.
- `autolimits`: Automatically determine frequency limits for deconvolution. A pass band is
  determined for all frequencies with the lower and upper corner cutoffs defined in terms of
  dB down from the maximum amplitude. This algorithm is designed to work with flat responses,
  i.e. a response in velocity for an instrument which is flat to velocity. Other combinations
  will likely result in unsatisfactory results. Cannot be used in combination with the
  freqlimits option.
- `units`: Specify output units. Can be DIS, VEL, ACC or DEF, where DEF results in no unit
  conversion.  `units` can only be used with `correct`.
- `diff`: Differentiate using 2 point (uncentered) method. true or false (default false)
- `int`: Integrate using trapezoidal (midpoint) method. true or false (default false)
- `decimate`: Sample-rate to decimate to. See Help for more details. A linear-phase,
  anti-alias filter is applied during decimation.

### Plot Options
- `antialias`: If true, the created image will have anti-aliasing used. Can only be specified
  with the `output=plot` option
- `width`: Width of the output plot (pixels). Can only be specified with the `output=plot` option
- `height`: Height of the output plot (pixels). Can only be specified with the `output=plot` option

### Audio Options
- `audiosamplerate`: The sample rate of the output wav file in Hz. Defaults to 16000.
  Can only be specified with the output=audio option
- `audiocompress`: Apply dynamic range compression to waveform data. This makes more of the
  signal audible. Can only be specified with the output=audio option

## Output options

|Value|Description|
|:----|:----------|
|`ascii1`|ASCII data format, 1 column (values)|
|`ascii2`|ASCII data format, 2 columns (time, value)|
|`ascii`|Same as `ascii2`|
|`geocsv`, `geocsv.tspair`|ASCII GeoCSV data format, 2 column (time, values)|
|`geocsv.slist`|ASCII geocsv data format, 1 columns (value)|
|`audio`|Audio WAV file|
|`miniseed`|FDSN miniSEED format|
|`plot`|A simple plot of the timeseries in PNG format|
|`saca`|SAC – ASCII format|
|`sacbb`|SAC – binary big-endian format|
|`sacbl`|SAC – binary little-endian format|

"""
struct IRISTimeSeries <: IRISRequest
    network::String
    station::String
    location::String
    channel::String
    starttime::DateTime
    endtime::DateTime
    duration::MFloat
    antialias::MBool
    width::MInt
    height::MInt
    audiosamplerate::MInt
    audiocompress::MBool
    output::String
    process::OrderedDict{Symbol,Any} # Because order matters for processing

    function IRISTimeSeries(network, station, location, channel, starttime, endtime, duration,
                            antialias, width, height, audiosamplerate, audiocompress, output,
                            process)
        location == "  " && (location = "--")
        !ismissing(starttime) && !ismissing(endtime) && starttime > endtime &&
            throw(ArgumentError("`starttime` must be before endtime"))
        coalesce(duration, 1) > 0 || throw(ArgumentError("`duration` must be positive"))
        0 <= get(process, :taper, 0.1) <= 1 || throw(ArgumentError("`taper` length must be between 0 and 1"))
        # Check types for processing
        for (f, T) in IRISTimeSeries_PROCESSING_FIELDS
            if haskey(process, f)
                try
                    process[f] = convert(T, process[f])
                catch err
                    err isa MethodError || err isa InexactError &&
                        throw(ArgumentError("field `$f` must be of type $T"))
                    rethrow(err)
                end
            end
        end
        if haskey(process, :taper_type)
            haskey(process, :taper) || throw(ArgumentError("`taper_type` must be specficied with taper"))
            taper_type = uppercase(process[:taper_type])
            taper_type in ("HANNING", "HAMMING", "COSINE") ||
            throw(ArgumentError("taper_type must be one of \"HANNING\", \"HAMMING\" or \"COSINE\""))
        end
        haskey(process, :scale) && haskey(process, :divscale) &&
            throw(ArgumentError("Cannot specify both `scale` and `divscale`"))
        if haskey(process, :scale)
            if process[:scale] isa AbstractString
                process[:scale] = uppercase(process[:scale])
                process[:scale] == "AUTO" || throw(ArgumentError("`scale` must be a float or \"AUTO\""))
            end
        end
        if haskey(process, :units)
            !haskey(process, :correct) || !process[:correct] &&
                throw(ArgumentError("`units` can only be used when correct is true"))
            process[:units] in ("DIS", "VEL", "ACC", "DEF") ||
                throw(ArgumentError("`units` must be one of \"DIS\", \"VEL\", \"ACC\" or \"DEF\""))
        end
        haskey(process, :autolimits) && haskey(process, :freqlimits) &&
            throw(ArgumentError("Cannot specify both `autolimits` and `freqlimits`"))

        # Disallowed combinations
        coalesce(audiocompress, false) || !ismissing(audiosamplerate) && output != "audio" &&
            throw(ArgumentError("Cannot specify `audiocompress` if `output` is not \"audio\""))

        # Disallowed values
        output in IRISTimeSeries_OUTPUT_TYPES ||
            throw(ArgumentError("`output` must be one of $(IRISTimeSeries_OUTPUT_TYPES)"))

        new(network, station, location, channel, starttime, endtime, duration,
            antialias, width, height, audiosamplerate, audiocompress, output, process)
    end
end

function IRISTimeSeries(;
                        network=missing, station=missing, location=missing, channel=missing,
                        starttime=missing, endtime=missing, duration=missing, antialias=missing,
                        width=missing, height=missing, audiosamplerate=missing,
                        audiocompress=missing, output=missing, kwargs...)
    any(ismissing, (network, station, location, channel, starttime, output)) &&
        throw(ArgumentError("network, station, location, channel, starttime " *
                            "and output must all be specified"))
    all(ismissing, (endtime, duration)) &&
        throw(ArgumentError("One of `endtime` or `duration` must be specified"))
    process = OrderedDict{Symbol,Any}()
    for (k, v) in kwargs
        k in keys(IRISTimeSeries_PROCESSING_FIELDS) ||
            throw(ArgumentError("Field `$k` is not a valid IRISTimeSeries processing field"))
        process[k] = v
    end
    IRISTimeSeries(network, station, location, channel, starttime, endtime, duration,
                   antialias, width, height, audiosamplerate, audiocompress, lowercase(output),
                   process)
end

# To be the same as the structs which use Parameters.@with_kw
Base.show(io::IO, p::SeisRequests.IRISTimeSeries) = dump(IOContext(io, :limit => true), p, maxdepth=1)

service_string(::IRISTimeSeries) = "timeseries"

function request_uri(r::IRISTimeSeries; server=DEFAULT_SERVER)
    server = server in keys(SERVERS) ? SERVERS[server] : server
    protocol = protocol_string(r)
    service = service_string(r)
    version = version_string(r)
    uri = join((server, protocol, service, version, "query?"), "/")
    firstfield = true
    for f in fieldnames(typeof(r))
        if f == :process
            for (field, v) in r.process
                field_val = if field == :taper
                    value = "$v"
                    if haskey(r.process, :taper_type)
                        value = join((value, uppercase(r.process[:taper_type])), ",")
                    end
                    "$field=$value"
                elseif field == :taper_type
                    continue
                elseif field in (:bpfilter, :freqlimits, :autolimits)
                    "$field=$(join(v, "-"))"
                else
                    "$field=$v"
                end
                @assert !firstfield # process is the last field
                uri = join((uri, field_val), "&")
            end
        else
            v = getfield(r, f)
            if !ismissing(v)
                if firstfield
                    uri *= "$f=$v"
                    firstfield = false
                else
                    uri = join((uri, "$f=$v"), "&")
                end
            end
        end
    end
    uri
end


text:
function ee(s0, K, vol, iter; dates = 1)
   vals = Array(Float64, dates)
   for ti in 1:dates
       x = rand(Float32, iter) - rand(Float32, iter)
       vol_ti = vol #- (ti * 0.0001)
       S = s0 * exp(x*vol_ti*sqrt(ti) .+ (- (vol_ti^2)/2*ti))
       ee_vals = max(S.-K, 0)
       vals[ti] = sum(ee_vals) / iter
   end
   vals
end

const kernel_source = "
float max(float x, float y){
   if(x >= y)
      return x;
   return y;
}

__kernel void ee(
        const float s0,
        const float K,
        const float vol,
        const int t,
        __global float* data)
{
        int i = get_global_id(0);
        float x = data[i];
        S = s0 * exp(x*vol*sqrt(t) + (- (vol^2)/2*t))
        ee_val = max(S - K, 0)
        vals[ti] = sum(ee_vals) / iter
}
"


text:
(r::Reagent)() = r(nothing)
function (r::Reagent)(a::A) where {A}
    @trace(label = :begin_react, taskid = objectid(current_task()))
    # backoff = Backoff()
    actr = then(r, Commit())
    ntries = 0

    # without offer
    while true
        ans = tryreact!(actr, a, Reaction(), nothing)
        @trace(
            label = :tryreact_without_offer,
            taskid = objectid(current_task()),
            offerid = offerid(nothing),
            ans
        )
        anchor(:tryreact_without_offer, (; ans))
        if ans isa SomehowBlocked
            if ans isa NeedNack
                runpostcommithooks(ans, nothing)
            end
            break
        elseif ans isa Retry
            GC.safepoint()  # once(backoff)?
            if maysync(r)
                break
            end
        else
            return ans
        end
        ntries += 1
        should_limit_retries() && ntries > 1000 && error("too many retries")
    end

    # with offer
    while true
        offer = Waiter{Any}()  # TODO: narrow type
        ans = tryreact!(actr, a, Reaction(), offer)
        @trace(
            label = :tryreact_with_offer,
            taskid = objectid(current_task()),
            offerid = offerid(offer),
            ans
        )
        anchor(:tryreact_with_offer, (; ans, offer))
        if ans isa SomehowBlocked
            wait(offer)
        elseif ans isa Retry
            yield()  # once(backoff)?
        else
            return ans
        end
        y = rescind!(offer)
        if ans isa NeedNack
            runpostcommithooks(ans, nothing)
        end
        y === nothing || return something(y)
        ntries += 1
        should_limit_retries() && ntries > 1000 && error("too many retries")
    end
end

Reagents.try!(r::Reagent, a = nothing) = ((r ⨟ Map(Some)) | Return(nothing))(a)
Reagents.trysync!(r::Reagent, a = nothing) = Reagents.trysyncing(r)(a)
Reagents.trysyncing(r::Reagent) = ReturnIfBlocked(nothing) ⨟ r ⨟ Map(Some)

hascas(actr::Reactor) = hascas(actr.reagent) || hascas(actr.continuation)

function Base.wait(offer::Waiter)
    if cas_weak!(offer.state, Pending(), Waiting()).success
        @trace(
            label = :start_wait,
            taskid = objectid(current_task()),
            offerid = offerid(offer),
        )
        wait()
        @trace(label = :stop_wait, taskid = objectid(current_task()))
    end
    return
end

function tryput!(offer::Waiter{T}, value::T) where {T}
    old, success = cas_weak!(offer.state, Pending(), value)
    success && return true
    old isa Rescinded && return false
    if old isa Waiting
        (_, success) = cas_weak!(offer.state, Waiting(), value)
        if success
            @trace(
                label = :tryput_wake_dual,
                taskid = objectid(current_task()),
                dual_taskid = objectid(offer.task),
            )
            schedule(offer.task)
        end
        return success
    else
        return false
    end
end

rescind!(::Nothing) = nothing
function rescind!(offer::Waiter)
    (old, success) = cas!(offer.state, Pending(), Rescinded())
    # success && return nothing
    old isa WaiterFlags && return nothing
    return Some(old)
end

function tryreact!(::Commit, a, rx::Reaction, offer::Union{Offer,Nothing})
    if offer isa Waiter
        ans = offer.state[]
        if ans isa Pending
            # Rescinding `offer` as part of commit, to keep the `offer` alive
            # until the last moment. This is required when CASes and channel
            # swaps are in multiple branches of `Choice`.
            # TODO: Check if this is OK. It's different from the paper.
            rx = withcas(rx, CAS(offer.state, Pending(), Rescinded()))
        elseif ans isa WaiterFlags
            return Retry()
        else
            @trace(
                label = :took_offer,
                taskid = objectid(current_task()),
                offerid = offerid(offer),
                ans,
            )
            return ans
        end
    end
    anchor(:commit, (; offer, rx))
    issuccess = commit!(rx)
    @trace(
        label = :commit,
        taskid = objectid(current_task()),
        offerid = offerid(offer),
        offer,
        rx,
        ans = a,
        issuccess,
    )
    if issuccess
        runpostcommithooks(rx, a)
        return a
    else
        return Retry()
    end
end

function runpostcommithooks(rx, @nospecialize(a))
    hooks = rx.postcommithooks
    if hooks !== nothing
        for f in hooks
            f(a)
        end
    end
end

hascas(::PostCommit) = false

function tryreact!(
    actr::Reactor{<:PostCommit},
    a,
    rx::Reaction,
    offer::Union{Offer,Nothing},
)
    (; f) = actr.reagent
    after_commit(_) = f(a)
    return tryreact!(actr.continuation, a, withpostcommit(rx, after_commit), offer)
end


text:
using Optimisers, Test
using Zygote, Random
using Statistics

@testset "Optimisers" begin
  Random.seed!(84)
  w′ = (α = rand(3, 3), β = rand(3, 3))
  @testset for o in (Descent(), ADAM(), Momentum(), Nesterov(), RMSProp(),
                     ADAGrad(), AdaMax(), ADADelta(), AMSGrad(), NADAM(),
                     ADAMW(), RADAM(), OADAM(), AdaBelief())
    w = (α = rand(3, 3), β = rand(3, 3))
    st = Optimisers.state(o, w)
    loss(x, y) = mean((x.α .* x.β .- y.α .* y.β) .^ 2)
    l = loss(w, w′)
    for i = 1:10^4
      gs = gradient(x -> loss(x, w′), w)
      st, w = o(st, w, gs...)
    end
    @test loss(w, w′) < 0.01
  end
end

@testset "OptimiserChain" begin
  Random.seed!(84)
  w = randn(10, 10)
  w′ = randn(10, 10)
  loss(x, w, w′) = mean((w*x .- w′*x) .^ 2)
  opt = OptimiserChain(WeightDecay(), ADAM(0.001))
  st = Optimisers.state(opt, w)
  for t = 1:10^5
    x = rand(10)
    gs = gradient(w -> loss(x, w, w′), w)
    st, w = Optimisers.update(opt, st, w, gs...)
  end
  @test loss(rand(10, 10), w, w′) < 0.01
end

@testset "Optimiser Updates" begin
  opt = ADAM()
  new_opt = ADAM(opt, eta = 9.f0)
  @test new_opt.eta == 9.f0
end


text:
function validate_libsize(libsize, source, dim, τ, η, replace)
    n_available_pts = length(source) - dim*τ - abs(η)
    if libsize > n_available_pts
        if replace
            #@warn "libsize = $libsize > n_available_pts = $n_available_pts. Sampling with replacement."
        else
            throw(DomainError(libsize, "libsize = $libsize > n_available_pts = $n_available_pts. Reduce `libsize` (preferably to something much smaller than `n_available_pts = $n_available_pts`) or enable sampling with replacement."))
        end
    end
end

function validate_embedding!(embedding, jitter)
    maxval = abs.(maximum(embedding))
    unique_pts = unique(embedding, dims = 2)
    if size(unique_pts, 2) < size(embedding, 2)
        #@warn "Not all embedding points are unique. Jittering coordinates by `rand(Uniform(-$jitter*maximum(embedding), $jitter*maximum(embedding)))`"

        for i in 1:length(embedding)
            embedding[i] += rand(Uniform(-jitter*maxval, jitter*maxval))
        end
    end
end

function validate_theiler_window!(theiler_window, points_available)
    if theiler_window < 0
        throw(DomainError(theiler_window, "`theiler_window=$theiler_window`. Must be ≧ 0."))
    end
    if theiler_window >= ceil(Int, 0.5*points_available)
        throw(DomainError(theiler_window, "`theiler_window=$theiler_window >= ceil(Int, 0.5*(length(target) - dim*τ))=$points_available`. Please reduce `theiler_window`."))
    end
end

function validate_embedding_params(dim, τ, points_available, theiler_window)
    if dim == 1
        @warn "`dim=$dim`, but must be at least 2 to construct an embedding."
    end
    if dim <= 0
        throw(DomainError(dim, "`dim=$dim` must be at least 2 to construct an embedding."))
    end
    if τ == 0
        @warn "`τ=$τ` does not produce a delay embedding!"
    end
    if τ < 0
        @warn "`τ=$τ` is negative. The cross mapping algorithm is implemented assuming τ is positive."
    end
    if dim*τ >= ceil(Int, 0.5*points_available)
        throw(DomainError(theiler_window, "`theiler_window=$theiler_window >= ceil(Int, 0.5*(length(target) - dim*τ))=$points_available`. Please reduce `theiler_window`."))
    end
end

function validate_uncertainty_measure(uncertainty_measure)
    if uncertainty_measure ∉ [:quantile, :std, :none]
        throw(DomainError(uncertainty_measure, "uncertainty_measure = $uncertainty_measure is invalid. Must be either :quantile, :std or :none."))
    end
end

function validate_average_measure(average_measure)
    if average_measure ∉ [:mean, :median, :none]
        throw(DomainError(average_measure, "average_measure = $average_measure is invalid. Must be either :median, :mean or :none."))
    end
end

function validate_output_selection(average_measure, uncertainty_measure, summarise)
    if ((average_measure, uncertainty_measure) == (:none, :none)) && summarise
        throw(ErrorException("Setting both average_measure and uncertainty_measure to :none while summarise = true will not produce output."))
    end
end

export
validate_libsize,
validate_embedding!,
validate_embedding_params,
validate_theiler_window!,
validate_uncertainty_measure,
validate_average_measure,
validate_output_selection


text:
# License for this file: MIT (expat)
# Copyright 2016-2018, DLR Institute of System Dynamics and Control

"""
    module test_withoutMacros_withoutVariables2

Test function ModiaMath.simulate(..) with models that have events
"""
module  test_withoutMacros_withoutVariables2

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
using ModiaMath.Test

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "SimpleStateEvents.jl"))
import .SimpleStateEvents

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "BouncingBall.jl"))
import .BouncingBall

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "models", "IdealClutch.jl"))
import .IdealClutch

@testset "\nTest ModiaMath: withoutMacros_withoutVariables/*.jl with events" begin

    @testset "Simulate SimpleStateEvents" begin
        m1 = SimpleStateEvents.Model()
        result = ModiaMath.simulate!(m1, stopTime=10.0, interval=0.1)
        time = result["time"]
        f    = result["f"]
        nStateEvents = m1.simulationState.statistics.nStateEvents
        tEvents = [1.62283; 3.36771; 6.51321; 8.04188]
        iEvents = [16 + 2, 33 + 2 + 2, 65 + 2 + 4, 80 + 2 + 6]

        @test length(time) == 101 + 4 * 2
        @test nStateEvents == 4
        @test isapprox(tEvents, time[iEvents]; atol=0.001 )
    end

    @testset "Simulate BouncingBall" begin
        m2 = BouncingBall.Model()
        result = ModiaMath.simulate!(m2, stopTime=3.0, interval=0.01)
        nStateEvents = m2.simulationState.statistics.nStateEvents

        @test nStateEvents == 18
    end

    @testset "Simulate IdealClutch" begin
        m2 = IdealClutch.Model()
        result = ModiaMath.simulate!(m2, stopTime=500.0)

        w1_end = result["inertia1.w"][end]
		w2_end = result["inertia2.w"][end]
		w_end_required = 38.9277466565

        @test isapprox(w1_end, w_end_required; atol=0.001 )
        @test isapprox(w2_end, w_end_required; atol=0.001 )
    end

end

include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_SimpleStateEvents.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_BouncingBall.jl"))
include(joinpath(ModiaMath.path, "examples", "withoutMacros_withoutVariables", "Simulate_IdealClutch.jl"))

end


text:
module PCHIPInterpolation
#
# Simple PCHIP implementation following Fritsch and Carlson, SIAM J. NUMER. ANAL. 17 (1980) 238-246.
# Derivative calculated in a fashion similar to SciPy's PchipInterpolate
#

export Interpolator, integrate

using ArgCheck: @argcheck

struct Interpolator{Xs,Ys,Ds}
    xs::Xs
    ys::Ys
    ds::Ds

    function Interpolator(xs::AbstractVector, ys::AbstractVector)
        @argcheck length(xs) ≥ 2
        @argcheck length(xs) == length(ys) DimensionMismatch
        foldl(xs) do a,b
            @argcheck a < b "xs must be stricly increasing"
            return b
        end

        ds = _initial_ds_scipy(xs, ys)
        new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds)
    end
end

ϕ(t) = 3t^2 - 2t^3
ψ(t) = t^3 - t^2

function _value_with_index(pchip::Interpolator, x::Number, i)
    x1, x2 = pchip.xs[i:i+1]
    @assert x1 ≤ x ≤ x2
    y1, y2 = pchip.ys[i:i+1]
    d1, d2 = pchip.ds[i:i+1]
    h = x2 - x1

    return (y1 * ϕ((x2-x)/h)
           + y2 * ϕ((x-x1)/h)
           - d1*h * ψ((x2-x)/h)
           + d2*h * ψ((x-x1)/h))
end

(pchip::Interpolator)(x::Number) = _value_with_index(pchip, x, _pchip_index(pchip, x))

function _integrate_segment(pchip::Interpolator, i, x1=nothing, x2=nothing)
    if isnothing(x1)
        x1 = pchip.xs[i]
        y1 = pchip.ys[i]
    else
        y1 = _value_with_index(pchip, x1, i)
    end

    if isnothing(x2)
        x2 = pchip.xs[i+1]
        y2 = pchip.ys[i+1]
    else
        y2 = _value_with_index(pchip, x2, i)
    end

    return (x2 - x1)/6 * (y1 + 4*_value_with_index(pchip, (x1 + x2)/2, i) + y2)  # Simpson's rule
end

function integrate(pchip::Interpolator, a::Number, b::Number)
    if b < a
        return -integrate(pchip, b, a)
    end

    i = _pchip_index(pchip, a)
    j = _pchip_index(pchip, b)

    if i == j
        return _integrate_segment(pchip, i, a, b)
    end

    integral = _integrate_segment(pchip, i, a, nothing) + _integrate_segment(pchip, j, nothing, b)

    for k in range(i+1, stop=j-1)
        integral += _integrate_segment(pchip, k)
    end

    return integral
end

function _pchip_index(pchip :: Interpolator, x)
    @argcheck (pchip.xs[1] <= x <= pchip.xs[end]) DomainError
    N = length(pchip.xs)
    if N < 200  # Approximate performance cross-over on my old intel i7-3517U
        i = _pchip_index_linear_search(pchip, x)
    else
        i = _pchip_index_bisectional_search(pchip, x)
    end
    if i == N
        # Treat right endpoint as part of rightmost interval
        i = N-1
    end
    i
end

function _pchip_index_linear_search(pchip :: Interpolator, x)
    i = 1
    N = length(pchip.xs)
    while i < N  &&  x >= pchip.xs[i+1]
        i = i + 1
    end
    i
end

function _pchip_index_bisectional_search(pchip :: Interpolator, x)
    N = length(pchip.xs)
    imin, imax = 1, N

    i = imin + div(imax - imin + 1, 2)
    while imin < imax
        if x < pchip.xs[i]
            imax = i - 1
        elseif x >= pchip.xs[i+1]
            imin = i + 1
        else
            break
        end
        i = imin + div(imax - imin + 1, 2)
    end
    i
end

"Similar to how SciPy's PCHIP does it"
function _initial_ds_scipy(xs::AbstractVector, ys::AbstractVector)
    h(i) = xs[i+1]-xs[i]
    Δ(i) = (ys[i+1]-ys[i]) / h(i)

    N = length(xs)
    ds = similar(ys./xs)
    if N == 2
        ds[:] .= Δ(1)
    else
        Δl = Δ(1)
        hl = h(1)
        for i ∈ 2:N-1
            Δr = Δ(i)
            hr = h(i)
            if sign(Δl) != sign(Δr) || Δl ≈ 0.0 || Δr ≈ 0.0
                ds[i] = 0.0
            else
                wl = 2hl + hr
                wr = hl + 2hr
                ds[i] = (wl + wr) / (wl/Δl + wr/Δr)
            end
            Δl = Δr
            hl = hr
        end
        ds[1] = _edge_derivative(h(1), h(2), Δ(1), Δ(2))
        ds[N] = _edge_derivative(h(N-1), h(N-2), Δ(N-1), Δ(N-2))
    end
    ds
end

function _edge_derivative(h1, h2, Δ1, Δ2)
    d = ((2h1 + h2)*Δ1 - h2*Δ2) / (h1 + h2)
    if sign(d) != sign(Δ1)
        d = 0.0
    elseif sign(Δ1) != sign(Δ2)  &&  abs(d) > abs(3Δ1)
        d = 3Δ1
    end
    d
end

end  # module


text:
using Test
using Libdl
using LinearAlgebra
using SparseArrays
using InplaceOps
using FFTW
using Wavelets
using PyCall
using JOLI

jo_type_mismatch_error_set(false)
#JOLI.jo_type_mismatch_warn_set(false)
verbose=false

# write your own tests here
stime=time()
include("test_joMatrixSingle.jl")
include("test_joMatrixProduct.jl")
include("test_joMatrixSum.jl")
include("test_joMatrixDifference.jl")
include("test_joLinearFunctionSingle.jl")
include("test_joLinearFunctionProduct.jl")
include("test_joLinearFunctionSum.jl")
include("test_joLinearFunctionDifference.jl")
include("test_joLinearFunctionAndMatrixProduct.jl")
include("test_joLinearFunctionAndMatrixSum.jl")
include("test_joLinearFunctionAndMatrixDifference.jl")
include("test_joInplaceOps.jl")
include("test_joMatrixCons.jl")
include("test_joKron.jl")
include("test_joCoreBlock.jl")
include("test_joBlock.jl")
include("test_joBlockDiag.jl")
include("test_joDict.jl")
include("test_joStack.jl")
include("test_joDFT.jl")
include("test_joNFFT.jl")
include("test_joDCT.jl")
include("test_joDWT.jl")
include("test_joRomberg.jl")
include("test_joSWT.jl")
include("test_joSincInterp.jl")
include("test_joCurvelet2D.jl")
include("test_joCurvelet2DnoFFT.jl")
include("test_joExtend.jl")
include("test_joGaussian.jl")
include("test_joOuterProd.jl")
etime=time()
dtime=etime-stime
println("\nTest Total elapsed time: ",round(dtime,digits=1),"s")


text:
#Write a program that prints all prime numbers. (Note: if your programming language does not support arbitrary size numbers, printing all primes up to the largest number you can easily represent is fine too.)

function prime_nums(n)
    i = 0
    while i < n
        i += 1
        if i % 2 == true
            println(i)
        end
    end
end

prime_nums(100)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using ControlSystems
using MAT

@testset "stuffABCD" begin
    # MATLAB output
    test = matread("resources/stuffABCD.mat")

    ntf = synthesizeNTF(5, 42, 1)
    forms = ["CRFB", "CRFF", "CRFBD", "CRFFD", "Stratos", "DSFB"]
    for i = 1:length(forms)
        form = forms[i]
        a, g, b, c = realizeNTF(ntf, form)
        ABCD = stuffABCD(a, g, b, c, form)
        @test isapprox(test["ABCD_$form"], ABCD)
    end

    # "CIFB" and "CIFF" expect the real part of zeros to be 1
    ntf = zpk([1, 1+2im, 1-2im, 1+4im, 1-4im], [2, 2+1im, 2-1im, 2+3im, 2-3im], 1, 1)
    forms = ["CIFB", "CIFF"]
    for i = 1:length(forms)
        form = forms[i]
        a, g, b, c = realizeNTF(ntf, form)
        ABCD = stuffABCD(a, g, b, c, form)
        @test isapprox(test["ABCD_$form"], ABCD)
    end

    # "PFF" expects zeros to be split into two groups separated by the 0.5
    # radian line in the complex plane
    ntf = zpk([0.5, 0.7+0.2im, 0.7-0.2im, 0.3+0.4im, 0.3-0.4im], [2, 2+1im, 2-1im, 2+3im, 2-3im], 1, 1)
    a, g, b, c = realizeNTF(ntf, "PFF")
    ABCD = stuffABCD(a, g, b, c, "PFF")
    @test isapprox(test["ABCD_PFF"], ABCD)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
@testset "losses" begin
    ŷ = collect(reshape(0.1:0.1:1.2, 3, 4))
    c = [1, 1, 3, 2]
    @test nllloss(ŷ, c) == -0.625
    test_rrule(nllloss, ŷ, c)
end


text:
module ISL
using ISL_jll
include("api.jl")

function version()
    strip(Base.unsafe_string(ISL.API.isl_version()))
end

end # module


text:
@testset "Transforms" begin
  for c in [Composition(1,2,3),
            Composition(0,2,3),
            Composition(1,2,0)]
    @test alrinv(alr(c)) == c
    @test clrinv(clr(c)) == c
    @test ilrinv(ilr(c)) == c
  end
end


text:
module MinimalRLCore

using Random

greet() = println("Hello Reinforcement Learning!")

export
    AbstractEnvironment,
    start!,
    step!,
    get_actions,
    get_reward,
    is_terminal
include("environment.jl")

export
    AbstractAgent,
    start!,
    step!
include("agent.jl")

export
    Episode,
    run_episode!
include("episode.jl")

export
    AbstractFeatureConstructors,
    create_features,
    feature_size
include("feature_constructors.jl")

# include("features/HashlessTileCoder.jl")

end # module


text:
using Optim
using Statistics

"""
    k_on, k_off, k_bleach, p_miss, nemitters =
        estimateparams(smd::SMD, clusterdata::Vector{Matrix{Float32}})

Estimate rate parameters from the clusters in `smd` and `clusterdata`.

# Description
The rate parameters `k_on`, `k_off`, and `k_bleach`, the probability of missing
a localization for a visible emitter `p_miss`, and the underlying number of
emitters `nemitters` are estimated from the pre-clusters in
`smd`/`clusterdata`.  This is done by assuming that pre-clusters are, on
average, representative of a single blinking event of a single emitter.
"""
function estimateparams(smd::SMD, clusterdata::Vector{Matrix{Float32}})
    # Compute some quantities needed from each of the preclusters.
    nclusters = size(clusterdata)[1]
    clusterdurations, nobservations = computeclusterinfo(clusterdata)

    # Estimate the sum of the off rate and the bleaching rate, as well as the
    # probability of missing a localization.
    k_offpbleach = -log(1 - 1/Statistics.mean(clusterdurations))
    k_offpbleach = isinf(k_offpbleach) ?  1.0 :  k_offpbleach
    p_miss = 1 - Statistics.mean(nobservations./clusterdurations)

    # Compute other parameters by fitting the observed number of localizations
    # over time.
    nlocperframe = counts(smd.framenum)
    nlocperframe = nlocperframe[nlocperframe .!= 0]
    frames = convert(Array{Float32}, unique(smd.framenum))
    nloccumulative = cumsum(nlocperframe)
    k_off = nclusters / nloccumulative[end]
    k_bleach = maximum([10^-5  k_offpbleach - k_off])
    k(k_on) = k_on + k_offpbleach
    l1(k_on) = k_on * k_bleach / k(k_on)
    l2(k_on) = k(k_on) - l1(k_on)
    model(x) = ceil(x[1]) * (1.0.-p_miss) * (x[2]/k(x[2])) *
        ((1/l1(x[2]))*(1.0.-exp.(-l1(x[2])*(frames.-1.0))) .-
        (1/l2(x[2]))*(1.0.-exp.(-l2(x[2])*(frames.-1.0))))
    costfunction(x) = Statistics.mean((nloccumulative .- model(x)).^2)
    lowerbound = [maximum(nlocperframe) 10^-5]
    upperbound = [nclusters nloccumulative[end]/frames[end]]
    x1_guess = ceil(nclusters * k_bleach)
    x1_guess = (x1_guess>lowerbound[1]) & (x1_guess<upperbound[1]) ?
        x1_guess : (upperbound[1]-lowerbound[1])/2.0
    initialguess = [x1_guess 1.0/frames[end]]
    inner_optimizer = NelderMead()
    optimresults = optimize(costfunction, lowerbound, upperbound, initialguess,
        Fminbox(inner_optimizer))
    xhat = optimresults.minimizer
    nemitters = xhat[1]
    k_on = xhat[2]

    return k_on, k_off, k_bleach, p_miss, nemitters
end


text:
using NiLangCore, BenchmarkTools

bg = BenchmarkGroup()
# pop!/push!
bg["NiLang"] = @benchmarkable begin
    @instr PUSH!(x)
    @instr POP!(x)
end seconds=1 setup=(x=3.0)

# @invcheckoff pop!/push!
bg["NiLang-@invcheckoff"] = @benchmarkable begin
    @instr @invcheckoff PUSH!(x)
    @instr @invcheckoff POP!(x)
end seconds=1 setup=(x=3.0)

# @invcheckoff pop!/push!
bg["NiLang-@invcheckoff-@inbounds"] = @benchmarkable begin
    @instr @invcheckoff @inbounds PUSH!(x)
    @instr @invcheckoff @inbounds POP!(x)
end seconds=1 setup=(x=3.0)

# Julia pop!/push!
bg["Julia"] = @benchmarkable begin
    push!(stack, x)
    x = pop!(stack)
end seconds=1 setup=(x=3.0; stack=Float64[])

# FastStack-inbounds-Any
bg["FastStack-inbounds-Any"] = @benchmarkable begin
    @inbounds push!(stack, x)
    @inbounds pop!(stack)
end seconds=1 setup=(x=3.0; stack=FastStack(10))

# Julia pop!/push!
bg["Julia-Any"] = @benchmarkable begin
    push!(stack, x)
    x = pop!(stack)
end seconds=1 setup=(x=3.0; stack=Any[])

# setindex
bg["setindex"] = @benchmarkable begin
    stack[2] = x
    x = 0.0
    x = stack[2]
end seconds=1 setup=(x=3.0; stack=Float64[1.0, 2.0])

# setindex-inbounds
bg["setindex-inbounds"] = @benchmarkable begin
    stack[2] = x
    x = 0.0
    x = stack[2]
end seconds=1 setup=(x=3.0; stack=Float64[1.0, 2.0])

# FastStack
bg["FastStack"] = @benchmarkable begin
    push!(stack, x)
    x = 0.0
    x = pop!(stack)
end seconds=1 setup=(x=3.0; stack=FastStack{Float64}(10))

# FastStack-inbounds
bg["FastStack-inbounds"] = @benchmarkable begin
    @inbounds push!(stack, x)
    x = 0.0
    @inbounds x = pop!(stack)
end seconds=1 setup=(x=3.0; stack=FastStack{Float64}(10))

tune!(bg)
run(bg)


text:
__precompile__()

module PolyChaosODE

using BenchmarkTools,
    DifferentialEquations,
    Distributions,
    FastGaussQuadrature,
    LaTeXStrings,
    LinearAlgebra,
    Parameters,
    PolyChaos,
    Plots,
    Statistics

include("grids.jl")
include("orthogonal_polynomial_basis.jl")
include("quadrature.jl")
include("stochastic_galerkin_ode.jl")
include("monte_carlo_gsa.jl")
include("utils.jl")

export TensorIntegrationGrid
export OrthogonalPolyBasis, LegendrePolyBasis, HermitePolyBasis
export StochGalerkinODE
export compute_expectation_and_diag_variance,
    mc_gsa_estimate_total_indices,
    plot_with_plus_minus_std,
    compute_total_order_sobol_indices

end

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
using FileIO
using Random
using StatsBase
using FFTViews
using ColorTypes

export generate

# We could just use Images.jl, but Images.jl significantly increases start-up time.
import Base.+
import Base./
+(a::RGB{Float32}, b::RGB{Float32}) = RGB{Float32}(a.r+b.r, a.g+b.g, a.b+b.b)
/(a::RGB{Float32}, b::Int64) = RGB{Float32}(a.r/b, a.g/b, a.b/b)

# Overwrite FFTViews method to use 1-based indexing.
FFTViews.indrange(i) = FFTViews.URange(first(i), last(i))
FFTViews._reindex(::Type{FFTView}, ind, i) = FFTViews.modrange(i, ind)
function Base.similar(A::AbstractArray, T::Type, shape::Tuple{FFTViews.FFTVRange,Vararg{FFTViews.FFTVRange}})
    all(x->first(x)==1, shape) || throw(BoundsError("cannot allocate FFTView with the first element of the range non-zero"))
    FFTViews.FFTView(similar(A, T, map(length, shape)))
end
function Base.similar(f::Union{Function,Type}, shape::Tuple{FFTViews.FFTVRange,Vararg{FFTViews.FFTVRange}})
    all(x->first(x)==1, shape) || throw(BoundsError("cannot allocate FFTView with the first element of the range non-zero"))
    FFTViews.FFTView(similar(f, map(length, shape)))
end
PeriodicArray = FFTView

function getImage(field_patterns, id_to_pattern; average_superpositions=false)
    output_patterns = map(field_patterns) do pattern_ids
        if average_superpositions
            pattern = sum(getindex.(Ref(id_to_pattern), collect(pattern_ids)))[end, 1] / length(pattern_ids)
        else
            pattern = length(pattern_ids) == 1 ? id_to_pattern[first(pattern_ids)] : zero(id_to_pattern[first(pattern_ids)])
            fill(pattern[end, 1], 1, 1)
        end
    end
    vcat([hcat(output_patterns[row, :]...) for row in 1:size(output_patterns)[1]]...)
end

const directions = [(1, 0), (0, 1), (-1, 0), (0, -1)]
const opposite = Dict((1,0) => (-1,0), (0,1) => (0, -1), (-1, 0) => (1, 0), (0, -1) => (0, 1))

function getNeighbors(idx, field_patterns)
    ((idx+CartesianIndex(direction), direction) for direction in directions
        if checkbounds(Bool, field_patterns, idx+CartesianIndex(direction)))
end
function disallowPattern(field_idx, pattern_id, field_patterns, pattern_adjacency, patterns_allowed)
    for (neighbor, direction) in getNeighbors(field_idx, field_patterns),
            adj_pattern_id in pattern_adjacency[direction][pattern_id]
        patterns_allowed[opposite[direction]][Tuple(neighbor)..., adj_pattern_id] -= 1
    end
end

function expandFrontier(field_idx, field_patterns, frontier)
    offset = 3
    min_row = field_idx[1]-offset
    max_row = field_idx[1]+offset
    min_col = field_idx[2]-offset
    max_col = field_idx[2]+offset
    surrounding_fields = [CartesianIndex(row, col) for (row, col) in Iterators.product(min_row:max_row, min_col:max_col)]
    surrounding_fields = filter(f -> checkbounds(Bool, field_patterns, f), surrounding_fields)
    for field in surrounding_fields
        if length(field_patterns[field]) > 1
            push!(frontier, field)
        end
    end
    delete!(frontier, field_idx)
end

function collapseField(collapse_field_idx, chosen_pattern_id, field_patterns, pattern_adjacency, patterns_allowed, fast, frontier)
    #  Constrain patterns_allowed with the new information.
    if fast
        expandFrontier(collapse_field_idx, field_patterns, frontier)
    end
    field = collect(field_patterns[collapse_field_idx])
    for pattern_id in field
        if pattern_id == chosen_pattern_id continue end
        disallowPattern(collapse_field_idx, pattern_id, field_patterns, pattern_adjacency, patterns_allowed)
    end
    field_patterns[collapse_field_idx] = Set([chosen_pattern_id])
    # Enforce new constraint in all other fields.
    field_idxs = [neighbor for (neighbor, _) in getNeighbors(collapse_field_idx, field_patterns)]
    while !isempty(field_idxs)
        field_idx = pop!(field_idxs)
        constrained_valid_patterns = Set([])
        for pattern_id in field_patterns[field_idx]
            if all(getNeighbors(field_idx, field_patterns)) do (_, direction)
                    patterns_allowed[direction][field_idx[1], field_idx[2], pattern_id] > 0
                end
                push!(constrained_valid_patterns, pattern_id)
            else
                disallowPattern(field_idx, pattern_id, field_patterns, pattern_adjacency, patterns_allowed)
            end
        end
        if length(constrained_valid_patterns) != length(field_patterns[field_idx])
            for (neighbor, _) in getNeighbors(field_idx, field_patterns)
                if (!fast) || (neighbor in frontier)
                    push!(field_idxs, neighbor)
                end
            end
        end
        field_patterns[field_idx] = constrained_valid_patterns
        if fast && length(field_patterns[field_idx]) == 1
            expandFrontier(field_idx, field_patterns, frontier)
        end
    end
end

"""
    generate(filename::String
        [; patternsize=2,
           width=16,
           height=16,
           periodic_input=false,
           periodic_output=false,
           mirror_input_horizontally=false,
           mirror_input_vertically=false,
           rotate_input_clockwise=false,
           rotate_input_anticlockwise=false,
           ground=nothing,
           seed=0,
           save_to_gif=false,
           fast=false]
    )

Apply the WaveFunctionCollapse algorithm to generate an image of dimension `width`*`height`
based on the input image at location `filename`. The algorithm splits the input image into
many tiles and then constructs the output image by places similar tiles next to each other.
`patternsize` is the width and height of a pattern. If `periodic_input` is true, the algorithm
will consider patterns in the input image that wrap around vertically and horizontally. If
`periodic_output` is true, the output image can be seamlessly repeated vertically and horizontally.
If the `mirror_input_` and `rotate_input_` parameters are true, the algorithm generates patterns
also from the mirrored and rotated input image. If `ground` is a tuple (y, x), the left bottom field
of the output is set to the input pattern with the upper left corner at (y, x). In practice, this
can enforce regularity on images with ground rows, like `samples/Flowers.png`. If `seed` != 0,
calling the algorithm with the same seed will generate the same output image, provided that the other
parameters are identical. If `save_to_gif` is true, the generation process gets animated and saved
as "output.gif".

If `fast` is true, an experimental version of the algorithm is used. In practice, it is 2-4 times as fast,
but may lead to more failed generation attempts when generating big outputs. The experimental algorithm is
almost like the normal algorithm, but it doesn't propagate pattern constraints through all fields. Instead,
it limits the propagation to a frontier around already collapsed fields. That frontier is currently 3 fields
thick, but that value is subject to experimentation.

    generate(input::Array{<:ColorTypes.RGB,2}
    [; kwargs...]
    )

Like `generate(filename [; kwargs...])`, but you can directly pass an array for the input image.

# Examples
```
generate(
    filename="samples/Flowers.png",
    width=20,
    height=20,
    patternsize=3,
    periodic_output=true,
    mirror_input_horizontally=true,
    seed=726259)
```
"""
function generate(input::Array{<:Any,2};
    patternsize=3,
    width=20,
    height=20,
    periodic_input=false,
    periodic_output=false,
    mirror_input_horizontally=false,
    mirror_input_vertically=false,
    rotate_input_clockwise=false,
    rotate_input_anticlockwise=false,
    ground=nothing,
    seed=0,
    save_to_gif=false,
    fast=false,
)
    generate(input, patternsize, width, height, periodic_input,
        periodic_output, mirror_input_horizontally, mirror_input_vertically,
        rotate_input_clockwise, rotate_input_anticlockwise, ground, seed,
        save_to_gif, fast
    )
end

generate(filename::String; kwargs...) = generate(FileIO.load(filename); kwargs...)

function generate(
    input::Array{<:Any,2}, patternsize, width, height, periodic_input, periodic_output,
    mirror_input_horizontally, mirror_input_vertically, rotate_input_clockwise,
    rotate_input_anticlockwise, ground, seed, save_to_gif, fast,
)
    input = convert(Array{RGB{Float32}, 2}, input)
    pattern_to_id = Dict()
    id_to_pattern = Dict()
    pattern_count = Dict{Int64, Int64}()
    images = []
    frontier = Set{CartesianIndex{2}}()

     # direction => pattern_id => Set([pattern_ids...])
    pattern_adjacency = Dict{Tuple{Int64,Int64},Dict{Int64,Set{Int64}}}()

    # direction => [width x height x patterns]. Pattern i at (x, y) overlaps with patterns
    # from the field at (x, y)+direction patterns_allowed[direction][x, y, i] times.
    patterns_allowed = Dict{Tuple{Int64,Int64},AbstractArray{Int64,3}}()

    if seed != 0
        Random.seed!(seed)
    else
        newseed = rand(1:1_000_000)
        Random.seed!(newseed)
        println("Seed: $newseed")
    end
    bound = patternsize-1
    if periodic_input
        input = PeriodicArray(input)
        bound = 0
    end
    ground_pattern_id = nothing
    for col in 1:size(input)[2]-bound, row in 1:size(input)[1]-bound
        pattern_variations = [input[row:row+patternsize-1, col:col+patternsize-1]]
        if mirror_input_horizontally
            append!(pattern_variations, map(p -> p[1:end, end:-1:1], pattern_variations))
        end
        if mirror_input_vertically
            append!(pattern_variations, map(p -> p[end:-1:1, 1:end], pattern_variations))
        end
        if rotate_input_clockwise
            append!(pattern_variations, map(rotr90, pattern_variations))
        end
        if rotate_input_anticlockwise
            append!(pattern_variations, map(rotl90, pattern_variations))
        end
        for pattern in pattern_variations
            if pattern in keys(pattern_to_id)
                pattern_count[pattern_to_id[pattern]] += 1
            else
                newid = length(id_to_pattern) + 1
                id_to_pattern[newid] = pattern
                pattern_to_id[pattern] = newid
                pattern_count[newid] = 1
            end
        end
        if (row, col) == ground
            ground_pattern_id = pattern_to_id[pattern_variations[1]]
            println("Ground:"); flush(stdout); display(id_to_pattern[ground_pattern_id])
        end
    end

    for direction in directions, p1 in keys(id_to_pattern), p2 in keys(id_to_pattern)
        (dy, dx) = direction
        lap1 = id_to_pattern[p1][(dy==1 ? 2 : 1):(dy==-1 ? end-1 : end), (dx==1 ? 2 : 1):(dx==-1 ? end-1 : end)]
        lap2 = id_to_pattern[p2][(dy==-1 ? 2 : 1):(dy==1 ? end-1 : end), (dx==-1 ? 2 : 1):(dx==1 ? end-1 : end)]
        overlap = lap1 == lap2
        if overlap
            # TODO: Properly initialize pattern_adjacency
            push!(get!(get!(pattern_adjacency, direction, Dict()), p1, Set()), p2)
        end
    end
    for direction in directions
        # WARNING: Requires that pattern IDs go from 1 to n.
        patterns = map(1:length(id_to_pattern)) do id
            length(get!(pattern_adjacency[direction], id, Set()))
        end
        patterns_allowed[direction] = [p for _ in 1:height, _ in 1:width, p in patterns]
        if periodic_output
            patterns_allowed[direction] = PeriodicArray(patterns_allowed[direction])
        end
    end

    field_patterns = [Set(keys(id_to_pattern)) for i in 1:height, j in 1:width]
    if periodic_output
        field_patterns = PeriodicArray(field_patterns)
    end
    if ground_pattern_id != nothing
        collapseField(CartesianIndex(height, 1), ground_pattern_id, field_patterns, pattern_adjacency, patterns_allowed, fast, frontier)
        if any(x -> length(x) == 0, field_patterns)
            error("Impossible to initialize bottom row with ground. Try another ground.")
        end
    end
    while !all(x -> length(x) == 1, field_patterns)
        # Find field with lowest entropy in field_patterns.
        weight_sum = sum(field_patterns) do valid_patterns
            sum(id -> pattern_count[id], valid_patterns)
        end
        field_entropies = map(field_patterns) do valid_patterns
            # Ignore already collapsed fields.
            if length(valid_patterns) <= 1
                return 10e9
            end
            randn()*1e-6 - sum(valid_patterns) do id
                weight = pattern_count[id]
                weight * log(weight/weight_sum)
            end
        end
        (_, idx_min) = findmin(field_entropies)
        field = collect(field_patterns[idx_min])
        choice = sample(field, Weights(getindex.(Ref(pattern_count), field)))
        collapseField(idx_min, choice, field_patterns, pattern_adjacency, patterns_allowed, fast, frontier)
        if any(x -> length(x) == 0, field_patterns)
            error("Impossible to complete generation. Restart the algorithm.")
        end
        if save_to_gif
            push!(images, getImage(field_patterns, id_to_pattern, average_superpositions=true))
        end
    end
    if save_to_gif
        save("output.gif", cat(images..., dims=[3]), fps=10)
    end
    getImage(field_patterns, id_to_pattern)
end


text:
# # Nonlinear bouncing ball

# ## Model

module bouncing_ball_nonlinear #jl

using ReachabilityAnalysis, ModelingToolkit

@taylorize function flow_down!(du, u, params, t)
    du[1] = u[2]
    du[2] = -9.8 + 0.1*(u[2])^2
    return du
end

@taylorize function flow_up!(du, u, params, t)
    du[1] = u[2]
    du[2] = -9.8 - 0.1*(u[2])^2
    return du
end

function bouncingBallNonlinear_model()
    var = @variables x, v

    ## hybrid automaton with state variables x, v
    HA = LightAutomaton(2)

    ## mode 1 ("down")
    X = HPolyhedron([x ≥ 0, v ≤ 0], vars)
    m1 = @system(x' = flow_down!(x), dim: 2, x ∈ X)

    ## mode 2 ("up")
    X = HPolyhedron([x ≥ 0, v ≥ 0], vars)
    m2 = @system(x' = flow_up!(x), dim: 2, x ∈ X)

    ## α transition down → up
    add_transition!(HA, 1, 2, 1)
    G = HalfSpace(x ≤ 0, vars)
    A = [1.0 0.0; 0.0 -0.8]
    Rα = ConstrainedLinearMap(A, G)  # v := -0.8v

    ## β transition up → down
    add_transition!(HA, 2, 1, 2)
    G = HalfSpace(v ≤ 0, vars)
    Rβ = ConstrainedIdentityMap(2, G)

    ## hybrid system
    S = HybridSystem(HA, [m1, m2], [Rα, Rβ], fill(AutonomousSwitching(), 2))

    return S
end

function model(X0)
    H = bouncingBallNonlinear_model()
    return IVP(H, X0)
end

end #jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
struct FieldMatrix{T,N,Lt<:AbstractArray{T},Dt<:AbstractArray{T,N}} <: AbstractMatrix{T}
    L :: Lt # Lower diagonal coefficients
    D :: Dt # Diagonal coefficients
    R :: CartesianIndices{N,NTuple{N,UnitRange{Int}}}
end
Base.size(p::FieldMatrix) = (s = length(p.R); (s,s))
Base.IndexStyle(::Type{<:FieldMatrix}) = IndexCartesian()
function Base.getindex(p::FieldMatrix{T}, i::Int, j::Int) where T
    if i == j
        return p.D[p.R[i]]
    else
        N = length(size(p.R))
        for d ∈ 1:N
            p.R[i]==p.R[j]+δ(d,N) && return p.L[p.R[i],d]
            p.R[i]==p.R[j]-δ(d,N) && return p.L[p.R[j],d]
        end
        return zero(T)
    end
end

@inline δ(i,N::Int) = CartesianIndex(ntuple(j -> j==i ? 1 : 0, N))
@fastmath @inline multL(I::CartesianIndex{N},L,x) where {N} =
    sum(@inbounds(x[I-δ(i,N)]*L[I,i]) for i ∈ 1:N)
@fastmath @inline multU(I::CartesianIndex{N},L,x) where {N} =
    sum(@inbounds(x[I+δ(i,N)]*L[I+δ(i,N),i]) for i ∈ 1:N)
@fastmath @inline mult(I::CartesianIndex{N},L,D,x) where {N} =
    @inbounds(x[I]*D[I])+multL(I,L,x)+multU(I,L,x)

import LinearAlgebra: mul!,dot,diag
mul!(b::FieldVector,p::FieldMatrix,x::FieldVector) = (@loop b[I]=mult(I,p.L,p.D,x); b)
@fastmath function dot(b::FieldVector,p::FieldMatrix,x::FieldVector)
    s = zero(eltype(x))
    @inbounds @simd for I ∈ b.R
        s+= b[I]*mult(I,p.L,p.D,x)
    end
    s
end
diag(p::FieldMatrix) = FieldVector(p.D)

import Base: *
*(p::FieldMatrix,x::FieldVector) = mul!(zero(x),p,x)

"""
    Poisson(L)

Construct a symmetric `FieldMatrix` from lower diagonal coefficients `L` with zero-sum rows.
This is a requirement for conservative Poisson equations: ∫ ∇⋅β∇ϕ dv = ∮ β ∂ϕ/∂n da.
"""
function Poisson(L::AbstractArray{T}) where T
    D = zeros(T,Base.front(size(L)))
    R = inside(D)
    for I ∈ R; D[I] = calcdiag(I,L); end
    FieldMatrix(L,D,R)
end
@fastmath @inline calcdiag(I::CartesianIndex{N},L) where {N} =
    -sum(@inbounds(L[I,i]+L[I+δ(i,N),i]) for i ∈ 1:N)