text:
struct DBFGS{T1, T2} <: QuasiNewton{T1}
   approx::T1
   theta::T2
end

DBFGS(approx) = DBFGS(;approx=approx)
DBFGS(;approx=Inverse(), theta=0.2) = DBFGS(approx, theta)

function update!(scheme::DBFGS{<:Direct, <:Any}, B, s, y)
   # We could write this as
   #     B .+= (y*y')/dot(s, y) - (B*s)*(s'*B)/(s'*B*s)
   #     B .+= (y*y')/dot(s, y) - b*b'/dot(s, b)
   # where b = B*s
   # But instead, we split up the calculations. First calculate the denominator
   # in the first term
   σ = dot(s, y)
   ρ = inv(σ) # scalar
   # Then calculate the vector b
   b = B*s # vector temporary
   sb = dot(s, b)
   if σ ≥ scheme.theta*sb
      θ = 1.0
      # Calculate one vector divided by dot(s, b)
      ρbb = inv(sb)*b
      # And calculate
      B .+= (inv(σ)*y)*y' .- ρbb*b'
   else
      θ = 0.8*sb/(sb-σ)
      r = y*θ + (1-θ)*b
      # Calculate one vector divided by dot(s, b)
      ρbb = inv(dot(s, b))*b
      # And calculate
      B .+= (inv(dot(s, r))*r)*r' .- ρbb*b'
   end
end
function update(scheme::DBFGS{<:Direct, <:Any}, B, s, y)
   # As above, but out of place

   σ = dot(s, y)
   b = B*s
   sb = dot(s, b)
   if σ ≥ scheme.theta*sb
      θ = 1.0
      # Calculate one vector divided by dot(s, b)
      ρbb = inv(sb)*b
      # And calculate
      return B .+= (inv(σ)*y)*y' .- ρbb*b'
   else
      θ = 0.8*sb/(sb-σ)
      r = y*θ + (1-θ)*b
      # Calculate one vector divided by dot(s, b)
      ρbb = inv(dot(s, b))*b
      # And calculate
      return B .+= (inv(dot(s, r))*r)*r' .- ρbb*b'
   end
end
function update(scheme::DBFGS{<:Inverse, <:Any}, H, s, y)
   σ = dot(s, y)
   ρ = inv(σ)
   #   if isfinite(ρ)
   C = (I - ρ*s*y')
   H = C*H*C' + ρ*s*s'
   #   end
   H
end
function update!(scheme::DBFGS{<:Inverse, <:Any}, H, s, y)
   σ = dot(s, y)
   ρ = inv(σ)

   if isfinite(ρ)
      Hy = H*y
      H .= H .+ ((σ+y'*Hy).*ρ^2)*(s*s')
      Hys = Hy*s'
      Hys .= Hys .+ Hys'
      H .= H .- Hys.*ρ
   end
   H
end

function update!(scheme::DBFGS{<:Inverse, <:Any}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end
function update!(scheme::DBFGS{<:Direct, <:Any}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end


text:
## This is an example of MSM that illustrates
#  the problem of "chatter" of the objective
#  function, which makes extremum estimation
#  difficult. It goes on to explore Bayesian
#  MSM, which is pretty easy to implement.

#  We generate data that follows a Poisson
#  distribution, and estimate using some
#  simple simulated moments

# the standard parameterization of Poisson density
function λ(x,θ)
    exp.(x*θ)
end

## generate the "true" data
using Distributions
n = 100
x = randn(n,2)
θ₀ = [1.,1.] # true params
y = rand.(Poisson.(λ(x,θ₀)))

## the statistics, averaged over observations
function k(x,y)
    vec(mean([y sqrt.(y) x.*y x.*sqrt.(y)], dims=1))
end

## the simulated moments: S replications
function simulated_moments(θ, S, x)
    [k(x, rand.(Poisson.(λ(x,θ)))) for _ = 1:S]
end

## Define the MSM criterion
S = 20 # number of simulation replications
m = θ -> k(x,y) - mean(simulated_moments(θ, S, x))
 # sums of squares of moments: corresponds to GMM with identity weight
function obj(θ)
    mean(abs2, m(θ))
end
function obj(a,b) # a method for plots
    obj([a,b])
end

## profile the MSM objective function
using Plots
gr()
θ₁ = range(0.6, length=100, stop=1.2)
θ₂ = θ₁
p2 = contour(θ₁, θ₂, (θ₁,θ₂)->obj(θ₁,θ₂),c=:viridis)
xlabel!("x")
ylabel!("y")
## The objective function looks horrible!
#  Nevertheless, let's attempt to do gradient-based estimation
using Econometrics
fminunc(obj, zeros(2))
# note that theta hat is same as start values, it didn't work

## Now, let's do Bayesian MSM, as suggested by
#  Chernozhukov and Hong (2003)
using Turing, AdvancedMH, Econometrics
k_data = k(x,y) # stats from the real data
#  Define the prior and asymptotic Gaussian likelihood of the moments
@model function MSM(k_data, S, x)
    θ ~ arraydist([Normal() for _=1:2]) # the prior (note: it's biased)
    # sample from the model, at the trial parameter value, and compute statistics
    ks = simulated_moments(θ, S, x)
    kbar = mean(ks) # mean of simulated statistics, over S replications
    Σ = cov(ks)  # estimated covariance of statistics
    k_data ~ MvNormal(kbar, Σ) # this corresponds to CUE GMM, but adding the determinant from the MVN
end

# do MCMC sampling
# look at acceptance rate, below, and adjust
# the tuning to make acceptance around 0.25
S = 100 # make this as large as possible, given computation time, when doing real research
tuning = 0.01
length = 10000
burnin = 1000
chain = sample(MSM(k_data, S, x),
    MH(:θ => AdvancedMH.RandomWalkProposal(MvNormal(zeros(2), tuning*eye(2)))),
    length; discard_initial=burnin)
display(chain)
plot(chain)
corner(chain)
c = Array(chain)
acceptance = size(unique(c[:,1]),1)[1] / size(c,1)
println("Acceptance rate: $acceptance")

ext: jl

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "integration" begin
    include("test_forwarddiff_utils.jl")
end


text:
__precompile__()

module mbPred
export predCount, predEps, ∧

# general iterface:
#		pred are 'until'-predicates, becoming false when a loop should stop
#		functions here are constructors for predicates
#		pred must have init::Bool=false keyword argument
#			on init==true initialize and return true, afterwards behaviour should be same as when newly created
#		? must accept arb number of inputs

# means of putting preds together...
# somehow pipe through v...

# and-connect predicates
∧(p...) = function pred(v...; init::Bool=false)
	return foldl( (A,a) -> A&&a, map( f->f(v..., init=init), p ) )
end

function predCount(max::Integer)
	local cnt = 0
	return function pred(v...; init::Bool=false)
		if init
			cnt = 0
			return true
		end
		cnt += 1
		return cnt <= max
	end
end

function predEps(ϵ,i=1)
	return function pred(v...; init::Bool=false)
		if init return true end
		return norm(v[i]) ≥ ϵ
	end
end

end


text:
"""Abstract cell type"""
abstract type AbstractCell end
"""Dimension of the cell"""
dim(c::AbstractCell) = throw(MethodError(dim,(typeof(c),)))
"""Get cell properties: `:index` & `:values`"""
Base.getindex(c::AbstractCell, k::Symbol) = throw(MethodError(getindex,(typeof(c),)))
"""Set cell properties: `:index` & `:values`"""
Base.setindex!(c::AbstractCell, v, k::Symbol) = throw(MethodError(setindex!,(typeof(c),)))
"""Cell comparison"""
==(a::AbstractCell, b::AbstractCell) = throw(MethodError(==,(typeof(a), typeof(b))))
"""Get cell faces"""
faces(c::AbstractCell) = throw(MethodError(faces,(typeof(c),)))

#=== Simples ===#
mutable struct Simplex{P} <: AbstractCell
    idx::Int
    vs::Vector{P}
    hash::UInt64
    Simplex{P}(idx, vals) where {P} = new(idx, vals, hash(vals))
    Simplex{P}() where {P} = new(0, P[])
end
Simplex(splx::Vector{P}; ordered=true) where {P} = Simplex{P}(0, ordered ? sort!(splx) : splx)
Simplex(splx::P...; ordered=true) where {P} = Simplex(P[splx...], ordered=ordered)

# Private methods

Base.convert(::Type{Simplex{P}}, v::Vector{P}) where {P} = Simplex{P}(0, sort!(v))
Base.hash(splx::Simplex) = splx.hash
Base.show(io::IO, splx::Simplex) = show(io, "Σ($(splx.vs))[$(splx.idx)]")
Base.eltype(::Type{Simplex{P}}) where {P} = P

# Public methods

dim(splx::Simplex) = length(splx.vs)-1

function Base.getindex(splx::Simplex, k::Symbol)
    if k == :index
        return splx.idx
    elseif k == :values
        return splx.vs
    else
        throw(KeyError(k))
    end
end

function Base.setindex!(splx::Simplex, v, k::Symbol)
    if k == :index
        splx.idx = v
    elseif k == :values
        splx.vs = v
    else
        throw(KeyError(k))
    end
end

==(a::Simplex, b::Simplex) = a[:values] == b[:values]

function faces(splx::Simplex{P}) where {P}
    faces = typeof(splx)[]
    for i in 1:dim(splx)+1
        face = copy(splx[:values])
        deleteat!(face,i)
        push!(faces, Simplex(face))
    end
    return faces
end

# Misc. methods

function volume(S::AbstractMatrix)
    d, vc = size(S)
    @assert d == vc-1 "Number of vertexes in simplex should be dim+1"
    v0 = S[:,1]
    return abs(det(S[:,2:end] .- v0))/prod(1:d)
end
volume(s::Simplex{Int}, X::AbstractMatrix) = volume(X[:,s[:values]])

#=== Cube ===#

mutable struct Cube{P<:Vector} <: AbstractCell
    origin::P
    extent::P
    function Cube{P}(o::P, x::P) where {P}
        @assert length(x) <= length(o) "Too many extent elements"
        new(o, x)
    end
end
Base.show(io::IO, c::Cube) = show(io, "Cube[$(c.origin) + $(c.extent)]")

# Public methods

dim(c::Cube) = length(findall(!iszero, c.extent))

function Base.getindex(c::Cube, k::Symbol)
    if k == :index || k == :values
        return (c.origin, c.extent)
    else
        throw(KeyError(k))
    end
end

function Base.setindex!(c::Cube, v, k::Symbol)
    if k == :index || k == :values
        c.origin = v[1]
        c.extent = v[2]
    else
        throw(KeyError(k))
    end
end

==(a::Cube, b::Cube) = a[:values] == b[:values]

# Misc. methods

function volume(c::Cube)
    sizes = similar(c.origin,0)
    for (i,ex) in enumerate(c.extent)
        if ex != 0.
            expoint = copy(c.origin)
            expoint[i] += ex
            l = norm(expoint - c.origin)
            push!(sizes, l)
        end
    end
    return prod(sizes)
end

# Simplex Iterators

struct Simplices{T}
    itr::T
    dim::Int
end
simplices(itr::T, dim::Int=-1) where T = Simplices{T}(itr, dim)
Base.show(io::IO, splxs::Simplices{T}) where T = print(io, "Simplex Iterator", splxs.dim < 0 ? "" : " (d=$(splxs.dim))", " for $T")


text:
using OmegaGrav
using MatterPower
using Dierckx
using DelimitedFiles, Printf

# %% Cosmological parameters
Ωm = 0.315

# %% Specify a redshift
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
ired = 1
z = redshift[ired]

# %% Read in linear and non-linear power spectra computed by CLASS
filename = @sprintf("data/matterpower05_z%1d_pk.dat", ired)
d = readdlm(filename, comments = true)
pk_spl = Spline1D(d[:, 1], d[:, 2])
filename = @sprintf("data/matterpower05_z%1d_pk_nl.dat", ired)
d = readdlm(filename, comments = true)
pknl_spl = Spline1D(d[:, 1], d[:, 2])

# %% Compute `Ωgrav = Ωm * W / 2`
println("redshift = ", z)
# Linear matter power spectrum
Ωgrav_pklin = ograv_pk(x -> pk_spl(x), z, Ωm)
println("Ωgrav(linear pk) = ", Ωgrav_pklin)
# Non-linear matter power spectrum
Ωgrav_pknl = ograv_pk(x -> pknl_spl(x), z, Ωm)
println("Ωgrav(non-linear pk) = ", Ωgrav_pknl)


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    fig = Figure(resolution = (450, 400))
    ax = Axis(fig,  aspect = 1, xlabel = "x", ylabel = "y")
    hmap = heatmap!(rand(20,20), colormap = :Spectral_11)
    cbar = Colorbar(fig, hmap, label = "values", width = 15, flipaxis = false,
                ticksize=15, tickalign = 1, height = Relative(3.55/4))
    fig[1, 1] = cbar
    fig[1, 2] = ax
    colgap!(fig.layout, 7)
    fig
    save(joinpath(@__DIR__, "output", "heatmapCbarLeft.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
export read, write, existsfile, mkdir
export filenames, filepaths, dirnames, dirpaths
export readmat, writemat
export writejls, readjls
# using MAT

existsfile(filename::AbstractString) = (s = stat(filename); s.inode!=0)

import Base.read
read(filename::AbstractString) = open(read, filename)

import Base.write
write(data::AbstractString, filename::AbstractString) = write(data, filename, "w")
function write(data::AbstractString, filename::AbstractString, mode)
    io = open(filename, mode)
    finalizer(close, io)
    write(io,data)
    close(io)
end

function filedirnames(path = pwd(); selector = isdir, hidden = false, withpath = false,
    recursive = false, suffix = "", regex = Nothing)
    # @show path selector hidden withpath
    files = readdir(path)
    r = filter(x->selector(joinpath(path,x)) && (hidden || x[1]!='.'), files)
    # @show r
    r = sort(r)
    r = withpath ? mapvec(r, x->joinpath(path,x)) : r
    f = String
    r = mapvec(r, f)
    if recursive
        f(x) = filedirnames(x;
            selector = selector, hidden = hidden, withpath = withpath, recursive = true)
        r = @p dirpaths path | mapvec f | flatten | concat r _
    end
    if !isempty(suffix)
        suffix = isa(suffix, AbstractString) ? [suffix] : suffix
        r = @p mapvec suffix (suf->@p filter r x->@p takelast x len(suf) | isequal suf) | flatten
    end
    if regex != Nothing
        r = @p filter r x->@p ismatch regex x
    end
    r
end
dirnames(path = pwd(); kargs...) = filedirnames(path; selector = isdir, kargs...)
dirpaths(path = pwd(); kargs...) = dirnames(path; withpath = true, kargs...)
filenames(path = pwd(); kargs...) = filedirnames(path; selector = x->!isdir(x), kargs...)
filepaths(path = pwd(); kargs...) = filenames(path; withpath = true, kargs...)

function readmat(filename)
    return matread(filename)
end

readmat(filename,variables) = readmat(filename, variables...)

function readmat(filename,variables...)
    r = Any[]
    mat = matopen(filename)
    try
        for v in variables
            r[v] = read(mat, v)
        end
    finally
        close(mat)
    end
end

function writejls(a,filename)
    tempfilename = @p concat filename "." randstring(10) ".tmp"
    open(fd->serialize(fd,a),tempfilename,"w")
    mv(tempfilename, filename, force = true)
    filename
end
readjls(filename) = open(Serialization.deserialize, filename, "r")


text:
module TestDimensions

using Mimi
using Test

import Mimi:
    compdef, AbstractDimension, RangeDimension, Dimension, key_type, first_period, last_period,
    ComponentReference, ComponentPath

dim_varargs = Dimension(:foo, :bar, :baz)   # varargs
dim_vec = Dimension([:foo, :bar, :baz]) # Vector
dim_range = Dimension(2010:2100)     # AbstractRange
rangedim = RangeDimension(2010:2100) # RangeDimension type
dim_vals = Dimension(4) # Same as 1:4

# RJP: unclear what was intended here, but you cannot instantiate an abstract type...
# dim_vals_abstract = AbstractDimension(dim_vals) # Abstract

@test key_type(dim_varargs) == Symbol
@test key_type(dim_vec) == Symbol
@test key_type(dim_range) == Int
@test key_type(dim_vals) == Int

@test length(dim_varargs) == 3
@test length(dim_vec) == 3
@test length(dim_vals) == 4
@test length(dim_range) == 91
@test length(rangedim) == 91

# Test iteration
@test [x.first for x in dim_varargs] == collect(keys(dim_varargs))
@test [x.first for x in dim_range]   == collect(keys(dim_range))

@test lastindex(dim_varargs) == :baz
@test lastindex(dim_vec) == :baz
@test lastindex(dim_range) == 2100
@test lastindex(dim_vals) == 4

@test Base.keys(rangedim) == [2010:2100...]
@test Base.values(rangedim) == [1:91...]

for i = 1:length(dim_range)
    @test getindex(dim_range, 2010 + i - 1) == i
end
@test dim_range[:] == [1:91...]
@test dim_varargs[:bar] == 2
@test dim_varargs[:] == [1,2,3]

# @test dim_vals_abstract[1:4...]== [1:4...]

# @test rangedim[2011] == 2 # TODO: this errors..

@test get(dim_varargs, :bar, 999) == 2
@test get(dim_varargs, :new, 4) == 4 #adds a key/value pair
@test get(rangedim, 2010, 1) == 1
# @test get(rangedim, 2101, 92) == 92 # TODO:  this errors ...

#test iteratable
dim_vals2 = Dimension(2:2:8)
intkeys = [2,4,6,8]
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in dim_vals2
        @test dim_vals2[intkeys[index]] == index
        index += 1
    end
end

rangedim2 = RangeDimension(2:2:8)
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in rangedim2   # uses iterate()
        @test get(rangedim2, Base.keys(rangedim2)[index]) == index
        index += 1
    end
end

@test getindex(dim_varargs, :bar) == 2
@test getindex(dim_varargs, :) == [1,2,3]

# Test resetting the time dimension

@defcomp foo2 begin
    x = Parameter(index=[time])
    y = Variable(index=[4])
end

m = Model()
set_dimension!(m, :time, 2000:2100)

# First and last have been disabled...
#@test_throws ErrorException add_comp!(m, foo2; first = 2005, last = 2105)   # Can't add a component longer than a model

foo2_ref = add_comp!(m, foo2)

foo2_ref = ComponentReference(m, :foo2)
my_foo2 = compdef(foo2_ref)

# First and last have been disabled...
# Can't set time more narrowly than components are defined as
# @test_throws ErrorException set_dimension!(m, :time, 1990:2200)
# @test first_period(my_foo2) == 2005
# @test last_period(my_foo2)  == 2095

# Test parameter connections
# @test_throws ErrorException set_param!(m, :foo2, :x, 1990:2200) # too long
# set_param!(m, :foo2, :x, 2005:2095) # Shouldn't throw an error

set_dimension!(m, :time, 2010:2050)

@test first_period(m.md) == 2010
@test last_period(m.md)  == 2050

# Test that d.time returns AbstracTimesteps that can be used as indexes

@defcomp bar begin
    v1 = Variable(index = [time])

    function run_timestep(p, v, d, t)
        if is_first(t)
            for i in d.time
                v.v1[i] = gettime(i)
            end
        end
    end
end

fixed_years = 2000:2010
variable_years = [2000, 2005, 2020, 2050, 2100]

m = Model()
set_dimension!(m, :time, fixed_years)
add_comp!(m, bar)
run(m)
@test m[:bar, :v1] == fixed_years

set_dimension!(m, :time, variable_years)
run(m)
@test m[:bar, :v1] == variable_years

# test variable_dimensions function
comp_path = ComponentPath((:bar,))
dims = [:time]
@test variable_dimensions(m, comp_path, :v1) == dims
@test variable_dimensions(m, :bar, :v1) == dims
@test variable_dimensions(m, (:bar,), :v1) == dims

end #module


text:
using Graphs
using Base.Test

function setofsets(array_of_arrays)
    Set(map(Set, array_of_arrays))
end

function test_cliques(graph, expected)
    # Make test results insensitive to ordering
    setofsets(maximal_cliques(graph)) == setofsets(expected)
end

g = simple_adjlist(3, is_directed=false)
add_edge!(g, 1, 2)
@test test_cliques(g, Array[[1,2], [3]])
add_edge!(g, 2, 3)
@test test_cliques(g, Array[[1,2], [2,3]])

# Test for "pivotdonenbrs not defined" bug
h = simple_inclist(6, is_directed=false)
add_edge!(h, 1, 2)
add_edge!(h, 1, 3)
add_edge!(h, 1, 4)
add_edge!(h, 2, 5)
add_edge!(h, 2, 6)
add_edge!(h, 3, 4)
add_edge!(h, 3, 6)
add_edge!(h, 5, 6)

@test maximal_cliques(h) != []

# test for extra cliques bug

g2 = simple_inclist(7, is_directed=false)
add_edge!(g2,1,3)
add_edge!(g2,2,6)
add_edge!(g2,3,5)
add_edge!(g2,3,6)
add_edge!(g2,4,5)
add_edge!(g2,4,7)
add_edge!(g2,5,7)
@test test_cliques(g2, Array[[7,4,5], [2,6], [3,5], [3,6], [3,1]])


text:
import Pkg

# Julia Packages
using PGFPlots
using Profile
using StatProfilerHTML

# Allow for easy updating
using Revise

# Satellite Dynamics Packages
using SatelliteDynamics

# Load SatelliteTasking
using SatelliteTasking

# Add Map Plots
using PGFMapPlots

##
# Simulation Definition
##

# Set simulation time window
t_start = Epoch(2019, 1, 1, 0, 0, 0, tsys="UTC"); # Start of time span
t_end = Epoch(2019, 1, 2, 0, 0, 0, tsys="UTC"); # End of simulation time span

# Initialize Planning Problem
problem = SatPlanningProblem(t_start=t_start, t_end=t_end);

# Add constraints to problem
push!(problem.constraints, constraint_agility_single_axis)

# Create Spacecraft
line1 = "1     1U          19  1.00000000  .00000000  00000-0  00000-0 0    05"
line2 = "2     1  90.0000   0.0000 0010000   0.0000   0.0000 15.24162312    00"
sc1 = Spacecraft(id=1, tle=TLE(line1, line2), slew_rate=1.0);
push!(problem.spacecraft, sc1);

# Compute Orital Period
T = orbit_period(state(problem.spacecraft[1].tle, problem.spacecraft[1].tle.epoch)[1]);

println("Loadded $(length(problem.spacecraft)) spacecraft")

##
# Load Groundstation Data
##

stations = load_stations("./data/groundstations.json");
num_stations = length(stations);

# Load Image Data - [100, 200, 500, 1000, 2000, 5000, 10000, 12000, 15000]
num_requests = 1000
requests = load_requests("./data/landsat_test_$num_requests.json", id_offset=num_stations);
@assert num_requests == length(requests);

# Set Request Parameters
collect_duration = 30.0
set_collect_duration(requests, collect_duration)

# Add Locations To Problem
# add_locations(problem, stations)
add_locations(problem, requests)

##
# Compute Opportunities for collects and contacts
##

# Remove all existing opportunities
clear_opportunities(problem)

# Compute problem accesses
@time compute_access(problem, orbit_fraction=0.75)

println("Found $(length(problem.contacts)) passes in planning horizon")
println("Found $(length(problem.collects)) collects in planning horizon")
num_accessible_requests = sum([1 for req in problem.requests if length(problem.lt_loc_opps[req.id]) > 0])
pct_accessible_requests = floor(num_accessible_requests/length(problem.requests) * 10000)/100
println("Found $num_accessible_requests out of $(length(problem.requests)) requests have accesses. $pct_accessible_requests%")

##
# Baseline Solve (Single-Threaded)
##

# @time baseline_plan, baseline_reward = satellite_plan_baseline(problem, allow_repeats=false)

# # Check feasibility, and compute actual reward
# baseline_feasible, baseline_realized_reward, baseline_nur, baseline_ndr, baseline_nr, baseline_nc = analyze_plan(problem, baseline_plan);

# println("Computed plan is feasible: $(uppercase(string(baseline_feasible)))")
# println("Computed optimal schedule. Computed reward: $baseline_reward - Realized reward: $baseline_realized_reward.")
# println("Requests Satisfied: $baseline_nur/$(length(problem.requests)), Duplicate Collects: $baseline_ndr")
# println("Ground Contacts Taken: $baseline_nc")

##
# Graph Solve (Single-Threaded)
##

# @time graph_plan, graph_reward = satellite_plan_graph(problem, allow_repeats=false)

# # Check feasibility, and compute actual reward
# graph_feasible, graph_realized_reward, graph_nur, graph_ndr, graph_nr, graph_nc = analyze_plan(problem, graph_plan);

# println("Computed plan is feasible: $(uppercase(string(graph_feasible)))")
# println("Computed optimal schedule. Computed reward: $graph_reward - Realized reward: $graph_realized_reward.")
# println("Requests Satisfied: $graph_nur/$(length(problem.requests)), Duplicate Collects: $graph_ndr")
# println("Ground Contacts Taken: $graph_nc")

##
# Mixed-Integer Linear Programing Solve (Single-Threaded)
##

# @time milp_plan, milp_reward = satellite_plan_milp(problem, allow_repeats=false)

# # Check feasibility, and compute actual reward
# milp_feasible, milp_realized_reward, milp_nur, milp_ndr, milp_nr, milp_nc = analyze_plan(problem, milp_plan);

# println("")
# println("Computed plan is feasible: $(uppercase(string(milp_feasible)))")
# println("Computed optimal schedule. Computed reward: $milp_reward - Realized reward: $milp_realized_reward.")
# println("Requests Satisfied: $milp_nur/$(length(problem.requests)), Duplicate Collects: $milp_ndr")
# println("Ground Contacts Taken: $milp_nc")

##
# Markov Decision Process (No Resources)
##

# # Adjust solve parameters
# problem.solve_discount = 1.0 # Typical values 0.999 - 0.9999
# problem.solve_depth = 5
# problem.solve_breadth = 3
# problem.solve_horizon = T

# @time mdp_fs_plan, mdp_fs_reward = satellite_plan_mdp_fs(problem)

# # Check feasibility, and compute actual reward
# mdp_fs_feasible, mdp_fs_realized_reward, mdp_fs_nur, mdp_fs_ndr, mdp_fs_nr, mdp_fs_nc = analyze_plan(problem, mdp_fs_plan);

# println("Computed plan is feasible: $(uppercase(string(mdp_fs_feasible)))")
# println("Computed optimal schedule. Computed reward: $mdp_fs_reward - Realized reward: $mdp_fs_realized_reward.")
# println("Requests Satisfied: $mdp_fs_nur/$(length(problem.requests)), Duplicate Collects: $mdp_fs_ndr")
# println("Ground Contacts Taken: $mdp_fs_nc")

# # Profile MDP Planning code
# @profilehtml satellite_plan_mdp_fs(problem);

##
# Markov Decision Process - Monte Carlo Tree Search (No Resources)
##

# Adjust solve parameters
problem.solve_allow_repeats = false
problem.mdp_reward_scarcity = false
problem.solve_discount = 1.0 # Typical values 0.99 - 0.9999
problem.solve_depth = 10
problem.solve_breadth = 0
problem.solve_horizon = T
problem.mcts_sim_iterations = 30
problem.mcts_c = 3.0

mdp_mcts_plan, mdp_mcts_reward = satellite_plan_mdp_mcts(problem, parallel=false)

# Check feasibility, and compute actual reward
mdp_mcts_feasible, mdp_mcts_realized_reward, mdp_mcts_nur, mdp_mcts_ndr, mdp_mcts_nr, mdp_mcts_nc = analyze_plan(problem, mdp_mcts_plan);

println("Computed plan is feasible: $(uppercase(string(mdp_mcts_feasible)))")
println("Computed optimal schedule. Computed reward: $mdp_mcts_reward - Realized reward: $mdp_mcts_realized_reward.")
println("Requests Satisfied: $mdp_mcts_nur/$(length(problem.requests)), Duplicate Collects: $mdp_mcts_ndr")
println("Ground Contacts Taken: $mdp_mcts_nc")

# Profile MCTS
@profilehtml satellite_plan_mdp_mcts(problem);


text:
using Plots
using PGFPlotsX
using Printf
using CUDA
using CUDA.CUSPARSE
using OrdinaryDiffEq

include("../physical_params.jl")
include("mms_funcs.jl")
include("../domain.jl")
include("../numerical.jl")
include("../solvers.jl")

function plot_convergence(I_error, D_error, ns)

    @pgf push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{amsmath}")

    @pgf con_plot = LogLogAxis(
        {
            width = "7cm",
            height = "5cm",
            xlabel = raw"$N$",
            ylabel = raw"$||\Delta\boldsymbol{u} ||_{\boldsymbol{H}}$",
            legend_pos = "outer north east",
            legend_entries = [raw"Interseismic $p=2$ ",
                              raw"Interseismic $p=4$",
                              raw"Interseismic $p=6$",
                              raw"Coseismic $p=2$",
                              raw"Coseismic $p=4$",
                              raw"Coseismic $p=6$",],
        },
    )

    @pgf push!(con_plot, Plot(
        {
            color = "red",
            mark = "oplus",
        },
        Table(ns,I_error[1,:])
    )
               )
    @pgf push!(con_plot, Plot(
        {
            color = "blue",
            mark = "oplus",
        },
        Table(ns,I_error[2,:])
    )
               )
    @pgf push!(con_plot, Plot(
        {
            color = "green",
            mark = "oplus",
        },
        Table(ns,I_error[3,:])
    )
               )
        @pgf push!(con_plot, Plot(
        {
            color = "red",
            mark = "asterisk",
        },
        Table(ns,D_error[1,:])
    )
               )
    @pgf push!(con_plot, Plot(
        {
            color = "blue",
            mark = "asterisk",
        },
        Table(ns,D_error[2,:])
    )
               )
    @pgf push!(con_plot, Plot(
        {
            color = "green",
            mark = "asterisk",
        },
        Table(ns,D_error[3,:])
    )
               )

    pgfsave("converge_plot.tex", con_plot)
    #pgfsave("../../Basin_paper/figures/I_converge.tex", I_plot)

 end

function refine(ps, ns, Lw, D, B_p, RS, R, MMS)

    year_seconds = 31556952
    #xt, yt = transforms_e(Lw, .1, .05)
    # expand to (0,Lw) × (0, Lw)
    (x1, x2, x3, x4) = (0, Lw, 0, Lw)
    (y1, y2, y3, y4) = (0, 0, Lw, Lw)
    #xt, yt = transfinite(x1, x2, x3, x4, y1, y2, y3, y4)
    xt, yt = transforms_e(Lw, .2, .2)

    errI = Array{Float64, 2}(undef, (length(ps), length(ns)))
    errD = Array{Float64, 2}(undef, (length(ps), length(ns)))

    for γ in .2:.2#.3:-.05:.1
        @printf "gamma: %f\n\n\n" γ
    for (np, p) in enumerate(ps)

        @printf "Operator order: %d\n\n" p

        for (iter, N) in enumerate(ns)

            nn = N + 1

            @printf "\tNodes per Dimension: %d\n\n" nn
            flush(stdout)

            Nn = nn^2
            mt = @elapsed begin
                metrics = create_metrics(N, N, B_p, μ, ρ, xt, yt)
            end

            #@printf "Got metrics: %s s\n" mt
            facecoord = metrics.facecoord
            x = metrics.coord[1]
            y = metrics.coord[2]

            d_ops = operators(p, N, N, μ, ρ, R, B_p, metrics)

            b = repeat([.02], nn)
            τ̃f = Array{Float64, 1}(undef, nn)
            vf = Array{Float64, 1}(undef, nn)
            v̂_fric = Array{Float64, 1}(undef, nn)

            cpu_operators = (d_ops = d_ops,
                             nn = nn,
                             R = R,
                             fc = metrics.facecoord,
                             coord = metrics.coord,
                             B_p = B_p,
                             MMS = MMS,
                             RS = RS,
                             b = b,
                             sJ = metrics.sJ,
                             τ̃f = τ̃f,
                             v̂ = zeros(nn))

            threads = 512
            GS = 0.0
            GS += (length(d_ops.Λ.nzval) * 8)/1e9
            GS += length(metrics.sJ[1]) * 8/1e9
            GS += length(d_ops.Z̃f[1]) * 8/1e9
            GS += length(b) * 8/1e9
            GS += length(d_ops.L[1].nzval) * 8/1e9
            GS += length(d_ops.H[1].nzval) * 8/1e9
            GS += length(d_ops.JIHP.nzval) * 8/1e9
            GS += length(d_ops.nCnΓ1.nzval) * 8/1e9
            GS += length(d_ops.HIGΓL1.nzval) * 8/1e9

            #@printf "Estimated Gigabytes allocating to the GPU %f\n" GS

            #quit()

            ot = @elapsed begin
                #=
                GPU_operators = (nn = nn,
                threads = threads,
                blocks = cld(nn, threads),
                Λ = CuSparseMatrixCSC(d_ops.Λ),
                sJ = CuArray(metrics.sJ[1]),
                Z̃f = CuArray(d_ops.Z̃f[1]),
                L = CuSparseMatrixCSC(d_ops.L[1]),
                H = CuArray(diag(d_ops.H[1])),
                JIHP = CuSparseMatrixCSC(d_ops.JIHP),
                nCnΓ1 = CuSparseMatrixCSC(d_ops.nCnΓ1),
                nBBCΓL1 = CuSparseMatrixCSC(d_ops.nBBCΓL1),
                RS = CuArray([RS.a, RS.σn, RS.V0, RS.Dc, RS.f0, nn]),
                b = CuArray(b),
                τ̃f = CuArray(zeros(nn)))

                =#
            end
            #@printf "Got Operators: %s s\n" ot

            xf1 = metrics.facecoord[1][1]
            yf1 = metrics.facecoord[2][1]

            t_begin = 35 * year_seconds - .5
            t_final = 35 * year_seconds + .5

            u0 = he(x[:], y[:], t_begin, MMS)
            v0 = he_t(x[:], y[:], t_begin, MMS)
            q1 = [u0;v0]
            for i in 1:4
                q1 = vcat(q1, d_ops.L[i]*u0)
            end
            q1 = vcat(q1, ψe_hd(metrics.facecoord[1][1],
                               metrics.facecoord[2][1],
                               t_begin, B_p, RS, MMS))

            # q3 = CuArray(deepcopy(q1))
            q4 = deepcopy(q1)
            q5 = deepcopy(q1)
            @assert length(q1) == 2nn^2 + 5nn

            for dt_scale = .5 .^ (3:3)
                dt = dt_scale * d_ops.hmin / (sqrt(B_p.μ_out/B_p.ρ_out))
                if dt < (t_final - t_begin)
                    #@show dt_scale, dt
                    t_span = (t_begin, t_final)

                    #=
                    #@printf "Got initial conditions: %s s\n" it
                    st3 = @elapsed begin
                    timestep!(q3, FAULT_GPU!, GPU_operators, dt, t_span)
                    end
                    =#

                    #@printf "Ran GPU to time %s in: %s s \n\n" t_span[2] st3

                    st4 = @elapsed begin
                        timestep!(q4, MMS_FAULT_CPU!, cpu_operators, dt, t_span)
                    end

                    #@printf "Ran CPU MMS to time %s in: %s s \n\n" t_span[2] st4

                    #=
                    st5 = @elapsed begin
                    timestep!(q5, FAULT_CPU!, cpu_operators, dt, t_span)
                    end
                    =#

                    x = metrics.coord[1]
                    y = metrics.coord[2]

                    u_end4 = @view q4[1:Nn]
                    diff_u4 = u_end4 - he(x[:], y[:], t_span[2], MMS)
                    errD[np, iter] = sqrt(diff_u4' * d_ops.JH * diff_u4)

                    @printf "\n\t\tdynamic error with MS: %e\n" errD[np, iter]
                    if iter > 1
                        @printf "\t\tdynamic rate: %f\n" log(2, errD[np, iter - 1]/errD[np, iter])
                    end
                    flush(stdout)
                end
            end

            #@printf "L2 error displacements between CPU and GPU: %e\n\n" norm(u_end5 - u_end3)

            u = zeros(nn^2)
            ge = zeros(nn^2)
            vf = zeros(nn)

            t_final = 70 * year_seconds
            t_begin = 0.0 * year_seconds
            params = (t_final = t_final,
                      year_seconds = year_seconds,
                      reject_step = [false],
                      nn = nn,
                      Δτ = zeros(nn),
                      u = u,
                      ge = ge,
                      vf = vf,
                      RS = RS,
                      MMS = MMS,
                      B_p = B_p,
                      ops = d_ops,
                      b = b,
                      metrics = metrics,
                      counter = [0],
                      f1_source = ue_t)

            δ = 2 * he(xf1,
                       yf1,
                       t_begin,
                       MMS)

            ψ = ψe_hd(xf1,
                     yf1,
                     t_begin,
                     B_p,
                     RS,
                     MMS)

            ψδ = [ψ ; δ]

            t_span = (t_begin, t_final)

            #Q_STATIC_MMS!(params)
            #=
            dt_scale = .5 ^ 14
                dt = dt_scale * d_ops.hmin
            timestep!(ψδ, Q_DYNAMIC_MMS!, params, dt, t_span)
            =#

            prob = ODEProblem(Q_DYNAMIC_MMS_NOROOT!, ψδ, t_span, params)
            plotter = DiscreteCallback(PLOTFACE, terminate!)
            sol = solve(prob, Tsit5();
                        isoutofdomain=stepcheck,
                        atol = 1e-12,
                        rtol = 1e-12,
                        gamma = γ,
                        #dtmax = 1e3,
                        internalnorm=(x,_)->norm(x, Inf))
                        #callback=plotter)

            diff = params.u[:] .- he(x[:], y[:], t_span[2], MMS)

            errI[np, iter] = sqrt(diff' * d_ops.JH * diff)

            #=
            plt1 = contour(x[:, 1], y[1, :],
                           (reshape(params.u, (nn, nn)) .- he(x, y, t_span[2], MMS))',
                            title = "error", fill=true, yflip=true)

            plt2 = contour(x[:, 1], y[1, :],
                           he(x, y, t_span[2], MMS)',
                           fill = true, yflip=true, title = "exact")

            plt3 = contour(x[:, 1], y[1, :],
                           Forcing_h(x, y, t_span[2], B_p, MMS)',
                           fill=true, yflip=true, title = "forcing")
            plt4 = contour(x[:, 1], y[1, :],
                           reshape(params.u, (nn,nn))',
                           fill=true, yflip=true, title = "numerical")
            plot(plt1, plt2, plt3, plt4, layout=4)
            gui()
            =#

            #=
            plt5 = plot((d_ops.L[1] * u - he(xf1, yf1, sol.t[end], MMS)), yf1,
                        yflip = true, title = "face 1 error", legend=false)

            plt6 = plot(he(xf1, yf1, sol.t[end], MMS), yf1,
                        yflip = true, title = "face 1 exact", legend=false, xlims=(0, 1.1))

            plt7 = plot(d_ops.L[1] * u, yf1,
                        yflip = true, title = "face 1 numerical", legend=false, xlims=(0, 1.1))
            plot(plt5, plt6, plt7, layout=3)
            gui()
            =#

            @printf "\n\t\tquasi-dynamic error with MS: %e\n" errI[np, iter]
            if iter > 1
                @printf "\t\tquasi-dynamic rate: %f\n" log(2, errI[np, iter - 1]/errI[np, iter])
            end

            @printf "\t___________________________________\n\n"
            flush(stdout)
        end
    end
    end
    plot_convergence(errD, errI, ns)

end


text:
# Transform a range expression (e.g., `1:10`) into a (more-or-less) equivalent
# construct in JavaScript.
# NOTE: We could also do this be converting a range into an array and including
# that array verbatim in the JS code, but that's not legible if there are a
# significant number of elements in the array (e.g., it's not desirable to
# include an array literal with 100 elements in it in the generated JS code).
function crawl_call(::Val{:(:)}, start::Integer, stop::Integer)
    # Because JavaScript arrays are terrible, we have to fill the array with
    # undefined (otherwise it can't be iterated over) after we construct it with
    # the length specified in the constructor.
    return :(JSTerminal(jsstring($((
        "(",
        "new Array($(stop - start + 1))",
        ".fill(undefined)",
        ".map((_, i) => i + $start)",
        ")",
    )...))))
end

function crawl_call(::Val{:(:)}, args...)
    error("Range syntax cannot be parsed by JSExpr.")
end


text:
module Agents

using Requires
using Distributed
using DataStructures
using LightGraphs
using DataFrames
using Random
import ProgressMeter

import Base.length
import LinearAlgebra

# Core structures of Agents.jl
include("core/agents.jl")
include("core/model.jl")
include("core/space_interaction_API.jl")

# Existing spaces
include("spaces/nothing.jl")
include("spaces/graph.jl")
include("spaces/grid.jl")
include("spaces/discrete.jl")
include("spaces/continuous.jl")
include("spaces/openstreetmap.jl")
include("spaces/utilities.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")
include("simulations/ensemblerun.jl")

# Other features that exist in submodules
include("submodules/pathfinding/all_pathfinders.jl")
include("submodules/schedulers.jl")
include("submodules/io/AgentsIO.jl")
include("deprecations.jl")

# Predefined models
include("models/Models.jl")
export Models

# Update messages:
using Scratch
display_update = true
version_number = "4.5"
update_name = "update_v$(version_number)"

if display_update
    # Get scratch space for this package
    versions_dir = @get_scratch!("versions")
    if !isfile(joinpath(versions_dir, update_name))
        printstyled(
            stdout,
            """
            \nUpdate message: Agents v$(version_number)
            Please see the changelog online. Some key features:

            * `get_spatial_property` and `get_spatial_index` functions have been added for easier usage of spatially distributed properties in `ContinuousSpace`.
            * The old pathfinding API has been deprecated in favour of a complete rework:
              Pathfinding structs are no longer stored by the space. Instead, `AStar` structs should be created by passing in the space and other parameters. All pathfinding functions now require the `AStar` struct to be passed in.
            * Pathfinding is now supported for `ContinuousSpace` using the A* algorithm.
            * Additional utility functions `nearby_walkable` and `random_walkable` for use with pathfinding.
            """;
            color = :light_magenta,
        )
        touch(joinpath(versions_dir, update_name))
    end
end

end # module


text:
module ErrorBudget
using Cuas, ExternalBallistics
using Distributions
using DataFrames
using Interpolations

# Write your package code here.
include("types.jl")
include("paralax.jl")
include("variableBias.jl")
include("totalError.jl")

export μ, AimingError, AlignmentError, LofError, FlightError, σ, UnitEffect, windCrossError,windRangeError,interTable,rangeError,
        elevationError, azimuthError, muzzleVelError, temperatureError, cantError, Error, variableBias

end


text:
# CHA Monitoring
mutable struct CHAMonitor{T,N} <: AbstractMonitor
    # Collection of CHAUncorePMU
    pmus::Record{:CHA,Vector{T}}
    cha_numbers::Vector{IndexZero{Int}}

    # The CPU that we're
    # using for the MSR reading.
    cpu::IndexZero{Int}
    events::NTuple{N, UncoreSelectRegister}

    # The filters for the CHA events
    filter0::CHAFilter0
    filter1::CHAFilter1
end

"""
    CHAMonitor(events, socket, cpu; [program = true], [filter0], [filter1])

Monitor the Caching Home Agent counters for `events` on a **single** selected CPU `socket`.
This can gather information such as number of L3 hits and misses.
Argument `event` should be a `Tuple` of [`CounterTools.UncoreSelectRegister`](@ref) and
`socket` should be either an `Integer` or `IndexZero`.
Further, `cpu` is the CPU that will be actually reading the counters.
For best performance, `cpu` should be located on `socket`.

Filters
=======

The CHA Performance Monitoring Units allow counters to be filtered in various ways such
as issuing Core or Thread ID, request opcode etc.

These can be passed via the `filter0` and `filter1` keyword arguments and correspond to the
CHA filters 0 and 1 repectively.

Note: `filter0` should be a [`CounterTools.CHAFilter0`](@ref) and `filter1` should be a
[`CounterTools.CHAFilter1`](@ref).
"""
function CHAMonitor(
        events,
        socket,
        cpu;
        program = true,
        filter0::CHAFilter0 = CHAFilter0(),
        filter1::CHAFilter1 = CHAFilter1(),
    )

    # Convert the cpu into an IndexZero
    cpu = indexzero(cpu)

    # I used to look at just a subset of the available CHAs, but that seems to misconfigure
    # the last counters.
    #
    # Instead, just running them all seems to do the trick as the kernel/CPU simply
    # returns "0" for reads to MSRS that are out of range.
    cha_numbers = IndexZero.(0:27)

    # Now, we open up a MSR file to the given CPU
    handle = Handle(msrpath(cpu))

    # Construct the CHA PMUs
    #
    # Make a single buffer for all CHA PMUs
    buffer = zeros(UInt8, numbytes(CHA()))
    pmus = map(cha_numbers) do cha
        pmu = CHAUncorePMU(handle, cha, buffer)
        reset!(pmu)
        return pmu
    end |> Record{:CHA}

    monitor = CHAMonitor(
        pmus,
        cha_numbers,
        cpu,
        events,
        filter0,
        filter1
    )

    finalizer(reset!, monitor)
    program && program!(monitor)
    return monitor
end

mapleaves(f, monitor::CHAMonitor) = mapleaves(f, monitor.pmus)

function program!(monitor::CHAMonitor)
    f = x -> program(x, monitor.events, monitor.filter0, monitor.filter1)
    mapleaves(f, monitor)
    return nothing
end

function program(x::CHAUncorePMU, events, filter0, filter1)
    # Apply the filters
    setextra!(x, IndexZero(0), filter0)
    setextra!(x, IndexZero(1), filter1)

    # Program the counters
    for (i, evt) in enumerate(events)
        setcontrol!(x, i, evt)
    end

    # Enable the counters
    enable!(x)
end

reset!(monitor::CHAMonitor) = mapleaves(reset!, monitor)
Base.read(monitor::CHAMonitor) = mapleaves(getallcounters, monitor)

# TODO: Hack for now - find busses by socket for the CAPID 6 register,
# which provides a bitmask for what CHAs are available.
#
# NOTE: The documentation implies that there is a correlation between bit locations
# and the CHA numbers that are available.
#
# However, direct testing and the note found:
# https://lore.kernel.org/patchwork/patch/893681/
#
# Implies that only the number of set bits counts.
function cha_masks(socket = nothing)
    bus_numbers = 0:255
    device = 30
    fn = 3

    socket_to_number = UInt32[]
    for bus in bus_numbers
        path = pcipath(bus, device, fn)
        ispath(path) || continue

        pci = Handle(path)
        #  check that the vendor is Intel
        #seek(pci, 0)
        (read(pci, UInt32, IndexZero(0)) & VENDOR_ID_MASK) == DRV_IS_PCI_VENDOR_ID_INTEL || continue

        #seek(pci, 0x9C)
        value = read(pci, UInt32, IndexZero(0x9C))
        push!(socket_to_number, value)
    end
    return isnothing(socket) ? socket_to_number : socket_to_number[socket]
end


text:
# This file contains the code that allows things to be added to the registry

const ext2sym    = Dict{String, Union{Symbol,Vector{Symbol}}}()
const magic_list = Vector{Pair}()     # sorted, see magic_cmp below
const sym2info   = Dict{Symbol,Any}() # Symbol=>(magic, extension)
const magic_func = Vector{Pair}()     # for formats with complex magic #s

## OS:
abstract type OS end
abstract type Unix <: OS end
struct Windows <: OS end
struct OSX <: Unix end
struct Linux <: Unix end

split_predicates(list) = filter(x-> x <: OS, list), filter(x-> !(x <: OS), list)
applies_to_os(os::Vector) = isempty(os) || any(applies_to_os, os)
applies_to_os(os::Type{<:OS}) = false

applies_to_os(os::Type{<:Unix}) = Sys.isunix()
applies_to_os(os::Type{Windows}) = Sys.iswindows()
applies_to_os(os::Type{OSX}) = Sys.isapple()
applies_to_os(os::Type{Linux}) = Sys.islinux()

## Magic bytes:

# magic_cmp results in magic_list being sorted in order of increasing
# length(magic), then (among tuples with the same length) in
# dictionary order. This ordering has the advantage that you can
# incrementally read bytes from the stream without worrying that
# you'll encounter an EOF yet still have potential matches later in
# the list.
function magic_cmp(p::Pair, t::Tuple)
    pt = first(p)
    lp, lt = length(pt), length(t)
    lp < lt && return true
    lp > lt && return false
    pt < t
end
function magic_cmp(t::Tuple, p::Pair)
    pt = first(p)
    lp, lt = length(pt), length(t)
    lt < lp && return true
    lt > lp && return false
    t < pt
end

canonicalize_magic(m::NTuple{N,UInt8}) where {N} = m
canonicalize_magic(m::AbstractVector{UInt8}) = tuple(m...)
canonicalize_magic(m::String) = canonicalize_magic(codeunits(m))

## Load/Save

struct LOAD end
struct SAVE end

function add_loadsave(format, predicates)
    library = popfirst!(predicates)
    os, loadsave = split_predicates(predicates)
    if applies_to_os(os)
        if isempty(loadsave) || (LOAD in loadsave)
            add_loader(format, library)
        end
        if isempty(loadsave) || (SAVE in loadsave)
            add_saver(format, library)
        end
    end
end

## Add Format:

function add_format(fmt, magic, extension, load_save_libraries...)
    add_format(fmt, magic, extension)
    for library in load_save_libraries
        add_loadsave(fmt, library)
    end
    fmt
end

"""
`add_format(fmt, magic, extention)` registers a new `DataFormat`.
For example:

    add_format(format"TIFF", (UInt8[0x4d,0x4d,0x00,0x2b], UInt8[0x49,0x49,0x2a,0x00]), [".tiff", ".tif"])
    add_format(format"PNG", [0x89,0x50,0x4e,0x47,0x0d,0x0a,0x1a,0x0a], ".png")
    add_format(format"NRRD", "NRRD", [".nrrd",".nhdr"])

Note that extensions, magic numbers, and format-identifiers are case-sensitive.
"""
function add_format(fmt::Type{DataFormat{sym}}, magic::Union{Tuple,AbstractVector,String}, extension) where sym
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    m = canonicalize_magic(magic)
    rng = searchsorted(magic_list, m, lt=magic_cmp)
    if !isempty(m) && !isempty(rng)
        error("magic bytes ", m, " are already registered")
    end
    insert!(magic_list, first(rng), Pair(m, sym))  # m=>sym in 0.4
    sym2info[sym] = (m, extension)
    add_extension(extension, sym)
    fmt
end

# for multiple magic bytes
function add_format(fmt::Type{DataFormat{sym}},
                    magics::Tuple{T,Vararg{T}}, extension) where {sym, T <: Vector{UInt8}}
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    magics = map(canonicalize_magic, magics)
    for magic in magics
        rng = searchsorted(magic_list, magic, lt=magic_cmp)
        if !isempty(magic) && !isempty(rng)
            error("magic bytes ", magic, " are already registered")
        end
        insert!(magic_list, first(rng), Pair(magic, sym))  # m=>sym in 0.4
    end
    sym2info[sym] = (magics, extension)
    add_extension(extension, sym)
    fmt
end

# For when "magic" is supplied as a function (see the HDF5 example in
# registry.jl)
function add_format(fmt::Type{DataFormat{sym}}, magic, extension) where sym
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    push!(magic_func, Pair(magic,sym))  # magic=>sym in 0.4
    sym2info[sym] = (magic, extension)
    add_extension(extension, sym)
    fmt
end

"""
`del_format(fmt::DataFormat)` deletes `fmt` from the format registry.
"""
function del_format(fmt::Type{DataFormat{sym}}) where sym
    magic, extension = sym2info[sym]
    del_magic(magic, sym)
    delete!(sym2info, sym)
    del_extension(extension)
    nothing
end

# Deletes multiple magic bytes
del_magic(magic::Tuple, sym) = for m in magic
    del_magic(m, sym)
end
# Deletes single magic bytes
function del_magic(magic::NTuple{N, UInt8}, sym) where N
    rng = searchsorted(magic_list, magic, lt=magic_cmp)
    if length(magic) == 0
        fullrng = rng
        found = false
        for idx in fullrng
            if last(magic_list[idx]) == sym
                rng = idx:idx
                found = true
                break
            end
        end
        found || error("format ", sym, " not found")
    end
    @assert length(rng) == 1
    deleteat!(magic_list, first(rng))
    nothing
end

function del_magic(magic::Function, sym)
    deleteat!(magic_func, something(findfirst(isequal(Pair(magic,sym)), magic_func), 0))
    nothing
end

## File Extensions:

function add_extension(ext::String, sym)
    if haskey(ext2sym, ext)
        v = ext2sym[ext]
        if isa(v, Symbol)
            ext2sym[ext] = Symbol[v, sym]
        else
            push!(ext2sym[ext], sym)
        end
        return
    end
    ext2sym[ext] = sym
end
function add_extension(ext::Union{Array,Tuple}, sym)
    for e in ext
        add_extension(e, sym)
    end
end

del_extension(ext::String) = delete!(ext2sym, ext)
function del_extension(ext::Union{Array,Tuple})
    for e in ext
        del_extension(e)
    end
end


text:
# simple_LP.jl
# Use our LP formulation with the constraint that all edges have the same velocity
# Authored by Arthur J Delarue on 7/31/15

function simple_LP(
	manhattan::Manhattan,
	travelTimes::Array{Float64,2},
	numRides::Array{Int,2},
	testingData::Array{Float64,2},
	directory::AbstractString;
	model_type::AbstractString="relaxed",
	max_rounds::Int=5,
	turnCost::Float64=10.0,
	turnCostAsVariable::Bool=true,
	delta_bound::Array{Float64}=[0.06,0.05,0.04],
	maxNumPathsPerOD::Int=3,
	computeFinalSP::Bool=false)

	graph = manhattan.network
	roadTimes = deepcopy(manhattan.roadTime)
	distances = manhattan.distances
	positions = manhattan.positions
	roads = edges(graph)
	nodes = vertices(graph)
	out = [copy(out_neighbors(graph,i)) for i in nodes]
	inn = [copy(in_neighbors(graph,i)) for i in nodes]

	# Find nonzero travel times
	pairs = [find(travel_times[i,:]) for i=nodes]
	numDataPoints = sum([length(pairs[i]) for i=nodes])

	TESTDIR = directory

	# Create directory if necessary:
	if !isdir("Outputs/$TESTDIR")
		mkdir("Outputs/$TESTDIR")
	end

	# Tell us where output is going
	println("-- Saving outputs to directory Outputs/$(TESTDIR)/")
	turnCostFile = open("Outputs/$(TESTDIR)/tc-simple.csv", "w")

	# Create path storing array and initialize. By convention totalPaths[i,1] is the equality constraint
	numPaths = zeros(Int, numDataPoints)
	totalPaths = Array(Any, numDataPoints)
	totalNumExpensiveTurns = Array(Any, numDataPoints)
	for i = 1:numDataPoints
		totalPaths[i] = Array{Int}[]
		totalNumExpensiveTurns[i] = Int[]
		sizehint(totalPaths[i], maxNumPathsPerOD)
		sizehint(totalNumExpensiveTurns[i], maxNumPathsPerOD)
	end

	# Run over all pairs of nodes that have data
	srcs = Int[]
	dsts = Int[]
	for i in nodes, j in nodes
		if travelTimes[i,j] > 0
			# Load sources and destinations: will be useful when adding constraints
			push!(srcs, i)
			push!(dsts, j)
		end
	end

	status = 0
	newTimes = roadTimes
	old_objective = 0

	# Compute shortest paths (with turn cost)
	println("**** Computing shortest paths ****")
	@time new_graph, new_edge_dists, new_nodes, new_edge_isExpensive = modifyGraphForDijkstra(graph, newTimes, positions, turn_cost=turnCost)
	old_nodes = getInverseMapping(new_nodes, nv(new_graph))
	@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)

	# l is round number
	l = 1
	while l <= max_rounds
		println("###### ROUND $l ######")
		actual_l = l
		# Update bound on relaxation variables
		if model_type == "relaxed"
			delta = (1/2)^(div(actual_l-1,length(delta_bound)))*delta_bound[((actual_l-1)%length(delta_bound)) + 1]
			# Avoid numerical issues
			if delta < 1e-4 || l == max_rounds
				delta = 0
			end
		end

		# Create JuMP model for LP and QP
		println("**** Creating LP instance ****")
		m = Model(solver=GurobiSolver(TimeLimit=10000, OutputFlag=1, Method=3))

		# Add one variable for each road, for each model
		@defVar(m, t[i=nodes,j=out[i]] >= 0.038 * manhattan.distances[i,j])

		# Add decision variable for turn cost and inverse velocity
		if turnCostAsVariable
			@defVar(m, tc >= 0)
		else
			tc = turnCost
		end

		@defVar(m, w >= 0)
		@addConstraint(m, oneVelocity[i=nodes, j=out[i]], t[i,j] == w * distances[i,j])

		# Epsilon is the variable used for the absolute value difference between T and \hat{T}
		@defVar(m, epsilon[i=nodes,j=pairs[i]] >= 0)

		# Create handles for inequality constraints
		if maxNumPathsPerOD > 1
			@defConstrRef inequalityPath[1:numDataPoints, 1:(max_rounds-1)]
		end

		# Set first LP objective
		@setObjective(m, Min, sum{ sqrt(numRides[i,j]/travelTimes[i,j]) * epsilon[i,j], i=nodes, j=pairs[i]})

		# Add path constraints
		println("**** Adding constraints ****")
		tic()
		paths, numExpensiveTurns = reconstructMultiplePathsWithExpensiveTurnsParallel(new_sp.previous, srcs, dsts, old_nodes, new_sp.real_destinations, new_edge_isExpensive)
		for i=1:numDataPoints
			index = findfirst(totalPaths[i], paths[i])
			# If path already in paths
			if index != 0
				totalPaths[i][1], totalPaths[i][index] = totalPaths[i][index], totalPaths[i][1]
				totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][index] = totalNumExpensiveTurns[i][index], totalNumExpensiveTurns[i][1]			# New path is equality constraint
			# If we still have space to add the path
			elseif numPaths[i] < maxNumPathsPerOD
				numPaths[i] += 1
				if numPaths[i] == 1
					push!(totalPaths[i], paths[i])
					push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
				else
					push!(totalPaths[i], paths[i])
					push!(totalNumExpensiveTurns[i], numExpensiveTurns[i])
					assert(numPaths[i] == length(totalPaths[i]))
					totalPaths[i][numPaths[i]], totalPaths[i][1] = totalPaths[i][1], totalPaths[i][numPaths[i]]
					totalNumExpensiveTurns[i][numPaths[i]], totalNumExpensiveTurns[i][1] = totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][numPaths[i]]
				end
			# If we need to remove a path
			else
				worstIndex = findWorstPathIndex(totalPaths[i], totalNumExpensiveTurns[i], turnCost, newTimes)
				if worstIndex == 1
					totalPaths[i][1] = paths[i]
					totalNumExpensiveTurns[i][1] = numExpensiveTurns[i]
				else
					totalPaths[i][1], totalPaths[i][worstIndex] = paths[i], totalPaths[i][1]
					totalNumExpensiveTurns[i][1], totalNumExpensiveTurns[i][worstIndex] = numExpensiveTurns[i], totalNumExpensiveTurns[i][1]
				end
			end
			# Add inequality constraints
			if numPaths[i] > 1
				for j = 1:(numPaths[i]-1)
					if model_type == "strict"
						inequalityPath[i,j] = @addConstraint(m, sum{t[totalPaths[i][j+1][a],totalPaths[i][j+1][a+1]], a=1:(length(totalPaths[i][j+1])-1)} - sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} >= - tc * totalNumExpensiveTurns[i][j+1] + tc * totalNumExpensiveTurns[i][1])
					elseif model_type == "relaxed"
						inequalityPath[i,j] = @addConstraint(m, sum{t[totalPaths[i][j+1][a],totalPaths[i][j+1][a+1]], a=1:(length(totalPaths[i][j+1])-1)} - sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} >= - tc * totalNumExpensiveTurns[i][j+1] + tc * totalNumExpensiveTurns[i][1] - delta * travelTimes[srcs[i],dsts[i]])
					end
				end
			end
		end
		# Equality constraints (shortest path close to travel time data)
		@addConstraint(m, TLowerBound[i=1:numDataPoints], sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} + tc * totalNumExpensiveTurns[i][1] - travelTimes[srcs[i],dsts[i]] >= - epsilon[srcs[i],dsts[i]])
		@addConstraint(m, TUpperBound[i=1:numDataPoints], sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} + tc * totalNumExpensiveTurns[i][1] - travelTimes[srcs[i],dsts[i]] <= epsilon[srcs[i],dsts[i]])
		toc()

		# Solve LP
		println("**** Solving LP ****")
		status = solve(m)

		# Debug if infeasible
		if status == :Infeasible
			println("!!!! Computing IIS !!!!")
			print_iis_gurobi(m)
			break
		# Prepare output
		elseif status == :Optimal || status == :Suboptimal
			if turnCostAsVariable
				turnCost = getValue(tc)
			end
			println(string("Left turn cost (s): ", turnCost))
			println(string("Average speed (kph): ", 3.6/getValue(w)))
			write(turnCostFile, string(l,",",turnCost, "\n"))
			result = getValue(t)
			newTimes = spzeros(length(nodes), length(nodes))
			for element in result
				newTimes[element[1], element[2]] = element[3]
			end
			# Save updated Manhattan road times to file
			saveRoadTimes(newTimes, "$TESTDIR/manhattan-simple-times-$l")
			objective = getObjectiveValue(m)
			if abs(objective - old_objective)/old_objective < 1e-3
				save("Outputs/$TESTDIR/simple-end.jld", "num_iter", l)
				l = max_rounds
			else
				old_objective = objective
			end
		elseif status == :UserLimit
			println("!!!! User time limit exceeded !!!!")
			break
		end
		if l < max_rounds || computeFinalSP
			println("**** Computing shortest paths ****")
			@time new_graph, new_edge_dists, new_nodes = modifyGraphForDijkstra(graph, newTimes, manhattan.positions, turn_cost=turnCost)
			@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
			compute_error(testingData, new_sp.traveltime)
		end
		l += 1
	end
	println("-- Saved outputs to directory Outputs/$(TESTDIR)/")
	close(turnCostFile)
	return status, newTimes, turnCost, totalPaths, totalNumExpensiveTurns
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
module TestHappyEyeballs

nolatencies(_) = nothing
const INJECTOR = Ref{Any}(nolatencies)

inject_latencies(islast) = INJECTOR[](islast)

include("../../../examples/happy_eyeballs.jl")

function with_latencies(f, injector)
    old = INJECTOR[]
    INJECTOR[] = injector
    try
        f()
    finally
        INJECTOR[] = old
    end
end

function sleep1_if_not_last(islast)
    if !islast
        sleep(1)
    end
end

function random_latencies(_)
    sleep(rand(0.1:0.1:1.5))
end

function test_happy_eyeballs_with_sleep1_if_not_last()
    with_latencies(sleep1_if_not_last) do
        @testset for trial in 1:3
            test_happy_eyeballs()
        end
    end
end

function test_happy_eyeballs_with_random_latencies()
    with_latencies(random_latencies) do
        @testset for trial in 1:3
            test_happy_eyeballs()
        end
    end
end

end  # module


text:
using Dierckx: Spline1D, evaluate

function splinefun(λ, n, name::Symbol)
    # Some of the tabulated data is not sorted
    p = sortperm(λ)
    λ = λ[p]
    n = n[p]
    range = [minimum(λ), maximum(λ)]
    if length(n) == 2
        k = 1
    elseif length(n) == 1
        # If there's only one data point return that value
        fun = quote
            function (λ)
                rangecheck(λ, $range)
                n[1]
            end
        end
        return eval(fun)
    else
        k = 1
    end
    spl = Spline1D(λ, n; k=k)
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            evaluate($spl, λ)
        end
    end
    eval(fun)
end

function f1(range, coeffs)
    # Sellmeier
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n2 = coefficients[1]
            for i = 2:2:length(coefficients)
                n2 += coefficients[i] * λ.^2 ./ (λ.^2 - coefficients[i+1].^2)
            end
            return sqrt.(n2 + 1)
        end
    end
end

function f2(range, coeffs)
    # Sellmeier-2
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n2 = coefficients[1]
            for i = 2:2:length(coefficients)
                n2 += coefficients[i] * λ.^2 ./ (λ.^2 - coefficients[i+1])
            end
            return sqrt.(n2 + 1)
        end
    end
end

function f3(range, coeffs)
    # Polynominal
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n2 = coefficients[1]
            for i = 2:2:length(coefficients)
                n2 += coefficients[i] * λ.^coefficients[i+1]
            end
            return sqrt.(n2)
        end
    end
end

function f4(range, coeffs)
    # RefractiveIndex.INFO
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n2 = coefficients[1]
            for i = 2:4:9
                n2 += coefficients[i] * λ.^(coefficients[i+1]) ./ (λ.^2 - coefficients[i+2]^coefficients[i+3])
            end
            for i = 10:2:length(coefficients)
                n2 += coefficients[i] * λ.^coefficients[i+1]
            end
            return sqrt.(n2)
        end
    end
end

function f5(range, coeffs)
    # Cauchy
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n2 = coefficients[1]
            for i = 2:2:length(coefficients)
                n2 += coefficients[i] * λ.^coefficients[i+1]
            end
            return n2
        end
    end
end

function f6(range, coeffs)
    # Gases
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n = coefficients[1]
            for i = 2:2:length(coefficients)
                n += coefficients[i] ./ (coefficients[i+1] - λ.^-2)
            end
            return n + 1.0
        end
    end
end

function f7(range, coeffs)
    # Herzberger
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n = coefficients[1] + coefficients[2] ./ (λ.^2 - 0.028) + coefficients[3] * (1 ./ (λ.^2 - 0.028)).^2
            for i = 4:length(coefficients)
                n += coefficients[i] * λ.^((i-3)*2)
            end
            return n
        end
    end
end

function f8(range, coeffs)
    # Retro
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            x = coefficients[1] + coefficients[2] * λ.^2 ./ (λ.^2 - coefficients[3]) + coefficients[4] * λ.^2
            return sqrt.((2x + 1) ./ (1 - x))
        end
    end
end

function f9(range, coeffs)
    # Exotic
    fun = quote
        function (λ)
            rangecheck(λ, $range)
            coefficients = $coeffs
            n = coefficients[1] + coefficients[2] ./ (λ.^2 - coefficients[3]) +
                coefficients[4] * (λ - coefficients[5]) ./ ((λ - coefficients[5]).^2 + coefficients[6])
            n = sqrt.(n)
        end
    end
end

function rangecheck(λ, range)
    lowerrange = range[1]
    upperrange = range[2]
    if any(λ .> upperrange) || any(λ .< lowerrange)
        warn("One or more values are out of range. ($lowerrange - $upperrange µm)")
    end
end


text:
"""
A basic mocking module, inspired by Python's [`unittest.mock`](https://docs.python.org/3/library/unittest.mock.html) and implemented with [Cassette](https://github.com/jrevels/Cassette.jl).
"""
module SimpleMock

using Base: Callable, invokelatest, unwrap_unionall
using Base.Iterators: Pairs
using Core: Builtin, kwftype

using Cassette: Cassette, Context, overdub, posthook, prehook, recurse, @context

export
    Call,
    Mock,
    mock,
    calls,
    ncalls,
    called,
    called_once,
    called_with,
    called_once_with,
    has_call,
    has_calls,
    reset!,
    max_depth,
    min_depth,
    excluding,
    including

include("metadata.jl")
include("filters.jl")
include("mock_type.jl")
include("mock_fun.jl")

end


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
    NCvar

Data structure containing information needed to write a NetCDF file. This
includes a list of filenames (see `Bindata`) if the data is not loaded into memory.
"""
struct NCvar
    name::String
    units::String
    dims
    values
    atts::Union{Dict,Nothing}
    backend::Module
end

function NCvar(name::String, units::String, dims::NCvar, values, atts::Union{Dict,Nothing}, backend::Module)
    return NCvar(name, units, [dims], values, atts, backend)
end

"""
    replacevalues(vals,ncvar::NCvar)

Helper function: Replaces the "values" attribute of an ncvar without changing
anything else. Used to apply land mask.
"""
replacevalues(vals,ncvar::NCvar) = NCvar(ncvar.name,ncvar.units,ncvar.dims,vals,ncvar.atts,ncvar.backend)

"""
    BinData

Data structure containing a string or an array of strings (NetCDF
    file names) as well as metadata needed to read a file.
"""
struct BinData # Pointer to data stored in binary files- contains info needed to read in
    fnames::Union{Array{String},String}
    precision::Type
    iosize::Tuple
    fldidx::Int
end

"""
    BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple)

Construct a BinData struct for files that contain one field.
"""
function BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple)
    return BinData(fnames,precision,iosize,1)
end

function BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple,meta::Dict,fldname::String)
    return BinData(fnames,precision,iosize,findfirst(meta["fldList"] .== fldname)[2])
end

"""
    NCData

Data structure containing a string or an array of strings (file names) of
    NetCDF files as well as information needed to read a file.
"""
struct NCData
    fname::AbstractString
    varname::AbstractString
    backend::Module
    precision::Type
end

"""
    TileData{T}

Data structure containing either a `MeshArray` struct or `BinData` struct (see `vals`),
    `MeshArray` structs describing the tile layout (`tileinfo`), and other information for
    reading/writing tile data.
"""
struct TileData{T}
    vals::T
    tileinfo::Dict
    tilesize::Tuple
    precision::Type
    numtiles::Int
end

"""
    TileData(vals,tilesize::Tuple)

Construct a TileData struct. First generate the tileinfo, precision, and numtiles attributes.
"""
function TileData(vals,tilesize::Tuple,grid::gcmgrid)
    ni,nj = tilesize
    if isa(vals,BinData)
        prec = vals.precision
    else
        prec = Float32
    end
    τ=Tiles(grid,ni,nj); TileMap=MeshArray(grid,prec);
    [TileMap[τ[ii]["face"]][τ[ii]["i"],τ[ii]["j"]].=ii for ii in 1:length(τ)]
    tileinfo = Dict("tileNo" => TileMap)
    gridvars=GridLoad(grid)
    tileinfo["XC"]=gridvars["XC"]
    tileinfo["YC"]=gridvars["YC"]
    return TileData(vals,tileinfo,tilesize,prec,Int(maximum(tileinfo["tileNo"])))
end
TileData(vals,tilesize::Tuple,grid::String="LatLonCap") = TileData(vals,tilesize::Tuple,GridSpec(grid))
"""
    replacevalues(vals,td::TileData)

Helper function: Replaces the "values" attribute of a TilData struct without changing
anything else. Used to apply land mask.
"""
replacevalues(vals,td::TileData) = TileData(vals,td.tileinfo,td.tilesize,td.precision,td.numtiles)

import Base: getindex

"""
    getindex(v::NCvar,i::Integer)

Retrive an NCvar at time step i.

Currently only works when values is a BinData, NCdata, or array of time steps.
"""
function getindex(v::NCvar,i::Integer)
    NCvar(v.name,v.units,v.dims,v.values[i],v.atts,v.backend)
end

"""
    getindex(d::BinData,i::Integer)

 Return a new BinData object with the ith file, this being time index i.
"""
function getindex(d::BinData,i::Integer) # Gets file at time index i
    if isa(d.fnames,Array)
        newfnames = [d.fnames[i]]
    elseif i == 1
        newfnames = d.fnames
    end

    BinData(newfnames,d.precision,d.iosize)
end

"""
    getindex(d::NCData,i)

Retrieve the ith time step from the NetCDF file.
"""
function getindex(d::NCData,i) # Gets data at time index i
    readncdata(d,i)
end

"""
    addDim(ds::NCDatasets.Dataset,dimvar::NCvar) # NCDatasets

Add a dimension to a `NCDatasets.Dataset`
"""
function addDim(ds::NCDatasets.Dataset,dimvar::NCvar)
    defDim(ds,dimvar.name,dimvar.dims[1])
end

"""
    addDim(dimvar::NCvar)

Add a dimension to a NetCDF file using `NetCDF.jl`
"""
function addDim(dimvar::NCvar) #NetCDF
    NcDim(dimvar.name,collect(dimvar.values),
    atts = merge(Dict(("units" =>dimvar.units)),dimvar.atts),
    unlimited = dimvar.dims==Inf)
end

"""
    addVar(ds::NCDatasets.Dataset,field::NCvar)

Add a variable to a NetCDF file using `NCDatasets.jl`
"""
function addVar(ds::NCDatasets.Dataset,field::NCvar)
    if ~isempty(field.units)
        atts = merge(Dict(("units" =>field.units)),field.atts)
    else
        atts = field.atts
    end
    if isa(field.values,Array)
        if isa(field.values[1],Array)
            prec = typeof(field.values[1][1])
        else
            prec = typeof(field.values[1])
        end
        fieldvar = defVar(ds,field.name,prec,
                            tuple([f.name for f in field.dims]...),attrib=atts)
    else
        fieldvar = defVar(ds,field.name,field.values.precision,
                            tuple([f.name for f in field.dims]...),attrib=atts)
    end
end

"""
    addVar(field::NCvar,dimlist::Array{NetCDF.NcDim})

Add a variable with dimensions dimlist to a NetCDF file using `NetCDF.jl`
"""
function addVar(field::NCvar,dimlist::Array{NetCDF.NcDim})
    if ~isempty(field.units)
        atts = merge(Dict(("units" =>field.units)),field.atts)
    else
        atts = field.atts
    end
    fieldvar = NcVar(field.name,dimlist,
                        atts = atts,
                        t = field.values.precision)
end

"""
    addVar(field::NCvar})

Add a variable and its dimensions to a NetCDF file using `NetCDF.jl`
"""
function addVar(field::NCvar)
    dimlist = addDim.(field.dims)
    if ~isempty(field.units)
        attributes = merge(Dict(("units" =>field.units)),field.atts)
    else
        attributes = field.atts
    end
    fieldvar = NcVar(field.name,dimlist,
        atts = merge(Dict(("units" =>field.units)),field.atts))
end

"""
    addData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar},var::NCvar)

Fill variable with data in netcdf file using either `NCDatasets.jl`
or `NetCDF.jl`
"""
function addData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar,Array},var::NCvar;startidx=1,land_mask=nothing)
    isBinData = isa(var.values,BinData)
    isNCData = isa(var.values,NCData)
    isTileData = isa(var.values,TileData)
    if isTileData
        isBinData = isa(var.values.vals,BinData)
    end
    nsteps = getnsteps(var.values)
    ndims = length(var.dims)
    if hastimedim(var)
        ndims = ndims-1
    end

    if isBinData || isNCData || isTileData || isa(var.values[1],Array) # Binary files or array of timesteps

        if isBinData && ~ isTileData
            if isa(var.values.fnames,Array)
                fnames = var.values.fnames
            else
                fnames = [var.values.fnames]
            end
        end

        if isTileData
            if isnothing(land_mask)
                gridvars = GridLoad(var.values.tileinfo["XC"].grid)
                if any(occursin.("_s",getfield.(var.dims,:name)))
                    land_mask = gridvars["hFacS"]
                elseif any(occursin.("_w",getfield.(var.dims,:name)))
                    land_mask = gridvars["hFacW"]
                else
                    land_mask = gridvars["hFacC"]
                end
                for f in land_mask.fIndex
                    for d in 1:size(land_mask,2)
                        land_mask[f,d][land_mask[f,d].==0] .= NaN
                        land_mask[f,d][land_mask[f,d].>0] .= 1
                    end
                end
            end
        end

        for i = startidx:nsteps
            if isBinData || isNCData || isTileData
                v0 = readdata(var.values,i)
            else
                v0 = var.values[i]
            end

            if ~isnothing(land_mask)
                if length(size(v0)) == 1 && length(size(land_mask)) == 2
                    v0 = v0 .* land_mask[:,1]
                elseif length(size(v0)) == 2 && length(size(land_mask)) == 3
                    v0 = v0 .* land_mask[:,:,1]
                else
                    v0 = v0 .* land_mask
                end
            end
            if isTileData
                tmpvar = replacevalues(replacevalues(v0,var.values),var)
                writetiles.(v,Ref(tmpvar),1:var.values.numtiles,Ref(i),Ref(land_mask))
            else
                checkdims(v0,var)
                if ndims == 0
                    v[i] = v0
                elseif ndims == 1
                    v[:,i] = v0
                elseif ndims == 2
                    v[:,:,i] = v0
                elseif ndims == 3
                    v[:,:,:,i] = v0
                end
            end
        end

    elseif isa(var.values[1],Number) # Single array of data- just insert it
        ndims = length(size(var.values))
        if ~isnothing(land_mask)
            v0 = var.values .* land_mask
        else
            v0 = var.values
        end
        checkdims(v0,var::NCvar)
        if ndims == 1
            v[:] = v0
        elseif ndims == 2
            v[:,:] = v0
        elseif ndims == 3
            v[:,:,:] = v0
        elseif ndims == 4
            v[:,:,:,:] = v0
        end
    else
        print("Unrecognized values")
    end
    return nothing
end

"""
    writetiles(v,var,tilenum,timeidx=1)

Helper function for writing a tile to a NetCDF file.
"""
function writetiles(v,var,tilenum,timeidx=1,land_mask=nothing)
    if isa(v,Array)
        v = v[findfirst(isequal(var.name),name.(v))]
    end
    tileinfo = var.values.tileinfo; tilesize = var.values.tilesize
    if isa(var.values.vals,MeshArray) || isa(var.values.vals,MeshArrays.gcmfaces)
        v0 = var.values.vals
    else
        v0 = readdata(var.values,timeidx)
    end

    v0 = gettile(v0,tileinfo,tilesize,tilenum)
    checkdims(v0,var::NCvar)
    numdims = ndims(v0)
    if numdims == 1
        v[:,timeidx] = v0
    elseif numdims == 2
        v[:,:,timeidx] = v0
    elseif numdims == 3
        v[:,:,:,timeidx] = v0
    end

end

"""
    addDimData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar,Array},var::NCvar)

Add dimension data to predefined dimensions in a NetCDF file.
"""
function addDimData(ds,dimvar::NCvar)
    atts = merge(Dict(("units" =>dimvar.units)),dimvar.atts)
    defVar(ds,dimvar.name,dimvar.values[:],(dimvar.name,),attrib=atts)
end

"""
    createfile(filename, field::Union{NCvar,Dict{String,NCvar}}, README;
                fillval=NaN, missval=NaN, itile=1, ntile=1)

Create NetCDF file and add variable + dimension definitions
using either `NCDatasets.jl` or `NetCDF.jl`
"""
function createfile(filename, field::Union{NCvar,Dict}, rdm="";
                    fillval=NaN, missval=NaN, itile=1, ntile=1, attribs=nothing)

    if isa(field,Dict)

        dims = getdims(field::Dict)
        #dims = unique(vcat([field[v].dims for v in keys(field)]...))
        #dims = filter( d -> isa(d,NCvar),dims)
        field = collect(values(field))

        fieldnames = getfield.(field[findall(hastimedim.(field))],:name)
        backend = field[1].backend
    else
        dims = field.dims
        fieldnames = field.name
        backend = field.backend
        if backend == NCDatasets
            field = [field]
        end
    end

    if isa(fieldnames,Array)
        fieldnamestring = join(filter(x -> x !== "climatology_bounds", fieldnames),",")
    else
        fieldnamestring = fieldnames
    end

    README=deepcopy(rdm)
    if ~isempty(README) && ~occursin(fieldnamestring*" -- ",README[1])
        README[1] = fieldnamestring*" -- "*README[1]
    end

    if ~isnothing(attribs)
        fillval = pop!(attribs,"_FillValue",NaN)
        missingval = pop!(attribs,"missing_value",NaN)
        haskey(attribs,"itile") ? itile = pop!(attribs,"itile") : nothing
        haskey(attribs,"ntile") ? ntile = pop!(attribs,"ntile") : nothing
        description = [pop!(attribs,k,"") for k in ["description","A","B","C","D","E","F","G","H","I","J"]]
        description = description[.~isempty.(description)]
        if isempty(README)
            README = description
        elseif !(description == README)
            README = vcat(README,description)
        end
        pop!(attribs,"date",nothing)
        if isempty(attribs); attribs=nothing; end
    end

    if ~isempty(README)
        if isa(README,Array) && length(README) > 1
            description = vcat("description" => README[1],
            [string(Char(65+(i-2))) => README[i] for i in 2:length(README)])
        else
            description = ["description" => isa(README,Array) ? README[1] : README]
        end
    else
        description = ""
    end

    file_atts = vcat(["date" => Dates.format(today(),"dd-u-yyyy"),
    "Conventions" => "CF-1.6"], "NCTiles" => NCTiles.version())

    if ~isempty(description)
        file_atts = vcat(file_atts,description)
    end
    file_atts = vcat(file_atts,
    [#string(Char(65+length(README)-1)) => "file created using NCTiles.jl",
    "_FillValue" => fillval,
    "missing_value" => missval,
    "itile" => itile,
    "ntile" => ntile])

    if ~isnothing(attribs)
        file_atts = vcat(file_atts,[k => attribs[k] for k in keys(attribs)])
    end

    if backend == NCDatasets
        ds =  Dataset(filename,"c",attrib=file_atts)
        dimlist = addDim.(Ref(ds),dims)
        fieldvar = addVar.(Ref(ds), field)
        if length(fieldvar) == 1
            fieldvar = fieldvar[1]
        end
        return ds,fieldvar,dimlist
    elseif backend == NetCDF # Not supporting multiple variables right now
        #= for i = 1:length(field)
            f = field[1];
            if isa(f.values,BinData) || isa(f.values,NCData)
                prec = f.values.precision
            else
                prec = eltype(f.values)
            end
            nccreate(filename,f.name,Tuple(vcat([[d.name; [collect(d.values)]; d.atts] for d in f.dims]...))...,atts = f.atts,t = prec)
        end
        #ncclose(filename)
        ds = NetCDF.open(filename)
        if length(field) == 1
            fieldvar = ds[field[1].name]
        else
            fieldvar = getkey.(Ref(ds),[f.name for f in field])
        end

        return ds,fieldvar,collect(values(ds.dim)) =#
        #fieldvar = Array{NcVar,1}
        dimlist = addDim.(dims)
        #for f in field
        #    fieldvar = vcat(fieldvar,addVar(f))
        #end
        field = field
        fieldvar = addVar(field)#[addVar(f) for f in field]
        file_atts = Dict(file_atts)
        #nccreate(filename,)
        return NetCDF.create(filename,fieldvar, gatts = file_atts),fieldvar,dimlist
    end

end

import Base: write
using Pkg
"""
    write(myfld::NCvar,savename::String;README="",globalattribs=Dict())

Creates NetCDF file and writes myfld and all its dimensions to the file.
"""
function write(myfld::NCvar,savename::String;README="",globalattribs=Dict())
    if hastiledata(myfld) # Create one file for each tile

        numtiles = myfld.values.numtiles
        savenames = savename*".".*lpad.(string.(1:numtiles),4,"0").*".nc"

        datasets = [createfile(savenames[tidx],myfld,README, itile = tidx, ntile = length(savenames), attribs = globalattribs) for tidx in 1:length(savenames)]

        ds = [x[1] for x in datasets]
        fldvars = [x[2] for x in datasets]

        addData(fldvars,myfld)

        if myfld.backend == NCDatasets
            dims = myfld.dims

            for dim in dims
                addDimData.(ds,Ref(dim))
            end

            close.(ds)
        else
            if Pkg.installed()["NetCDF"] < v"0.10"
                NetCDF.close.(ds)
            else
                finalize.(ds)
            end
        end
    else # Create one file
        ds,fldvar,dimlist = createfile(savename,myfld,README,attribs=globalattribs)

        # Add field and dimension data
        addData(fldvar,myfld)

        if myfld.backend == NCDatasets
            addDimData.(Ref(ds),myfld.dims)

            # Close the file
            close(ds)
        else
            if Pkg.installed()["NetCDF"] < v"0.10"
                NetCDF.close(ds)
            else
                finalize(ds)
            end
        end
    end
end

"""
    write(myflds::Dict,savename::String;README="",globalattribs=Dict())

Creates NetCDF file and writes myflds and all their dimensions to the file.
"""
function write(myflds::Dict{AbstractString,NCvar},savename::String;README="",globalattribs=Dict())

    if hastiledata(myflds)
        fldnames = collect(keys(myflds))
        tilefld = myflds[fldnames[findfirst([hastiledata(myflds[f]) for f in fldnames])]]
        numtiles = tilefld.values.numtiles

        landidx = findfirst(get.([myflds[f].atts for f in fldnames],"standard_name","none").=="land_binary_mask")
        if ~isnothing(landidx); land_mask = myflds[fldnames[landidx]].values; else land_mask = nothing; end
        if isa(land_mask,TileData); land_mask = land_mask.vals; end
        savenames = savename*".".*lpad.(string.(1:numtiles),4,"0").*".nc"

        datasets = [createfile(savenames[tidx],myflds,README, itile = tidx, ntile = length(savenames), attribs = globalattribs) for tidx in 1:length(savenames)]

        ds = [x[1] for x in datasets]
        fldvars = [x[2] for x in datasets]

        for k in keys(myflds)
            if isa(myflds[k].values,TileData)
                addData(fldvars,myflds[k],land_mask = land_mask)
            else
                tmpfldvars = [fv[findfirst(isequal(k),name.(fv))] for fv in fldvars]
                addData.(tmpfldvars,Ref(myflds[k]))
            end
        end

        dims = unique(vcat([myflds[v].dims for v in keys(myflds)]...))
        dims = filter( d -> isa(d,NCvar),dims)

        for dim in dims
            addDimData.(ds,Ref(dim))
        end

        close.(ds)

    else
        dims = getdims(myflds)

        ## Here down to move to other branch- high level API for writing
        ds,fldvar,dimlist = createfile(savename,myflds,README,attribs=globalattribs)

        # Add field data
        for k in keys(myflds)
            addData(ds[k],myflds[k])
        end

        # Only insert data for dims with data
        dims = filter( d -> ~isempty(d.values[:]),dims)

        # Add dimension data
        addDimData.(Ref(ds),dims)

        # Close the file
        close(ds)
    end
end

write(myflds::Dict,savename::String;README="",globalattribs=Dict()) = write(Dict{AbstractString,NCvar}(myflds),savename;README="",globalattribs=Dict())


text:
using URIs
using Test

include("uri.jl")
include("url.jl")


text:
using Test
using ParameterizedFunctions
using DiffEqBiological

# align
f = @ode_def TestAlignFeedback2 begin
    dx = y*c_1 - x
    dy = x^c_2 - y
end c_1 c_2

@test latexify(f; env=:align, cdot=false) ==
raw"\begin{align}
\frac{dx}{dt} =& y c_{1} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

@test latexify(f; env=:align, cdot=true) ==
raw"\begin{align}
\frac{dx}{dt} =& y \cdot c_{1} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

rn = @reaction_network begin
    hillr(D₂,α,K,n), ∅ --> m₁
    hillr(D₁,α,K,n), ∅ --> m₂
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    (k₊,k₋), 2P₁ ↔ D₁
    (k₊,k₋), 2P₂ ↔ D₂
    (k₊,k₋), P₁+P₂ ↔ T
end α K n δ γ β μ k₊ k₋;

@test latexify(rn; cdot=false) ==
raw"\begin{align}
\frac{dm_1}{dt} =& \frac{\alpha K^{n}}{K^{n} + D_2^{n}} - \delta m_1 + \gamma \\
\frac{dm_2}{dt} =& \frac{\alpha K^{n}}{K^{n} + D_1^{n}} - \delta m_2 + \gamma \\
\frac{dP_1}{dt} =& \beta m_1 - \mu P_1 -2 k_+ \frac{P_1^{2}}{2} + 2 k_- D_1 - k_+ P_1 P_2 + k_- T \\
\frac{dP_2}{dt} =& \beta m_2 - \mu P_2 -2 k_+ \frac{P_2^{2}}{2} + 2 k_- D_2 - k_+ P_1 P_2 + k_- T \\
\frac{dD_1}{dt} =& k_+ \frac{P_1^{2}}{2} - k_- D_1 \\
\frac{dD_2}{dt} =& k_+ \frac{P_2^{2}}{2} - k_- D_2 \\
\frac{dT}{dt} =& k_+ P_1 P_2 - k_- T
\end{align}
"

@test latexify(rn; cdot=true) ==
raw"\begin{align}
\frac{dm_1}{dt} =& \frac{\alpha \cdot K^{n}}{K^{n} + D_2^{n}} - \delta \cdot m_1 + \gamma \\
\frac{dm_2}{dt} =& \frac{\alpha \cdot K^{n}}{K^{n} + D_1^{n}} - \delta \cdot m_2 + \gamma \\
\frac{dP_1}{dt} =& \beta \cdot m_1 - \mu \cdot P_1 -2 \cdot k_+ \cdot \frac{P_1^{2}}{2} + 2 \cdot k_- \cdot D_1 - k_+ \cdot P_1 \cdot P_2 + k_- \cdot T \\
\frac{dP_2}{dt} =& \beta \cdot m_2 - \mu \cdot P_2 -2 \cdot k_+ \cdot \frac{P_2^{2}}{2} + 2 \cdot k_- \cdot D_2 - k_+ \cdot P_1 \cdot P_2 + k_- \cdot T \\
\frac{dD_1}{dt} =& k_+ \cdot \frac{P_1^{2}}{2} - k_- \cdot D_1 \\
\frac{dD_2}{dt} =& k_+ \cdot \frac{P_2^{2}}{2} - k_- \cdot D_2 \\
\frac{dT}{dt} =& k_+ \cdot P_1 \cdot P_2 - k_- \cdot T
\end{align}
"

@test latexify(rn; env=:chemical, cdot=false) ==
raw"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[\frac{\alpha K^{n}}{K^{n} + D_2^{n}}] m_{1}}\\
\ce{ \varnothing &->[\frac{\alpha K^{n}}{K^{n} + D_1^{n}}] m_{2}}\\
\ce{ m_{1} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{2} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{1} &->[\beta] m_{1} + P_{1}}\\
\ce{ m_{2} &->[\beta] m_{2} + P_{2}}\\
\ce{ P_{1} &->[\mu] \varnothing}\\
\ce{ P_{2} &->[\mu] \varnothing}\\
\ce{ 2 P_1 &<=>[k_{+}][k_{-}] D_{1}}\\
\ce{ 2 P_2 &<=>[k_{+}][k_{-}] D_{2}}\\
\ce{ P_{1} + P_{2} &<=>[k_{+}][k_{-}] T}
\end{align}
"

@test latexify(rn; env=:chemical, cdot=true) ==
raw"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + D_2^{n}}] m_{1}}\\
\ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + D_1^{n}}] m_{2}}\\
\ce{ m_{1} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{2} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{1} &->[\beta] m_{1} + P_{1}}\\
\ce{ m_{2} &->[\beta] m_{2} + P_{2}}\\
\ce{ P_{1} &->[\mu] \varnothing}\\
\ce{ P_{2} &->[\mu] \varnothing}\\
\ce{ 2 \cdot P_1 &<=>[k_{+}][k_{-}] D_{1}}\\
\ce{ 2 \cdot P_2 &<=>[k_{+}][k_{-}] D_{2}}\\
\ce{ P_{1} + P_{2} &<=>[k_{+}][k_{-}] T}
\end{align}
"

# tabular
@test latexify(rn.symjac; env=:tabular, cdot=false) ==
raw"\begin{tabular}{ccccccc}
$ - \delta$ & $0$ & $0$ & $0$ & $0$ & $\frac{ - K^{n} n \alpha D_2^{-1 + n}}{\left( K^{n} + D_2^{n} \right)^{2}}$ & $0$\\
$0$ & $ - \delta$ & $0$ & $0$ & $\frac{ - K^{n} n \alpha D_1^{-1 + n}}{\left( K^{n} + D_1^{n} \right)^{2}}$ & $0$ & $0$\\
$\beta$ & $0$ & $ - \mu - 2 k_+ P_1 - k_+ P_2$ & $ - k_+ P_1$ & $2 k_-$ & $0$ & $k_{-}$\\
$0$ & $\beta$ & $ - k_+ P_2$ & $ - \mu - k_+ P_1 - 2 k_+ P_2$ & $0$ & $2 k_-$ & $k_{-}$\\
$0$ & $0$ & $k_+ P_1$ & $0$ & $ - k_-$ & $0$ & $0$\\
$0$ & $0$ & $0$ & $k_+ P_2$ & $0$ & $ - k_-$ & $0$\\
$0$ & $0$ & $k_+ P_2$ & $k_+ P_1$ & $0$ & $0$ & $ - k_-$\\
\end{tabular}
"

@test latexify(rn.symjac; env=:tabular, cdot=true) ==
raw"\begin{tabular}{ccccccc}
$ - \delta$ & $0$ & $0$ & $0$ & $0$ & $\frac{ - K^{n} \cdot n \cdot \alpha \cdot D_2^{-1 + n}}{\left( K^{n} + D_2^{n} \right)^{2}}$ & $0$\\
$0$ & $ - \delta$ & $0$ & $0$ & $\frac{ - K^{n} \cdot n \cdot \alpha \cdot D_1^{-1 + n}}{\left( K^{n} + D_1^{n} \right)^{2}}$ & $0$ & $0$\\
$\beta$ & $0$ & $ - \mu - 2 \cdot k_+ \cdot P_1 - k_+ \cdot P_2$ & $ - k_+ \cdot P_1$ & $2 \cdot k_-$ & $0$ & $k_{-}$\\
$0$ & $\beta$ & $ - k_+ \cdot P_2$ & $ - \mu - k_+ \cdot P_1 - 2 \cdot k_+ \cdot P_2$ & $0$ & $2 \cdot k_-$ & $k_{-}$\\
$0$ & $0$ & $k_+ \cdot P_1$ & $0$ & $ - k_-$ & $0$ & $0$\\
$0$ & $0$ & $0$ & $k_+ \cdot P_2$ & $0$ & $ - k_-$ & $0$\\
$0$ & $0$ & $k_+ \cdot P_2$ & $k_+ \cdot P_1$ & $0$ & $0$ & $ - k_-$\\
\end{tabular}
"

f = @ode_def feedback begin
    dx = y/c_1 - x
    dy = x^c_2 - y
    end c_1 c_2

@test latexify(f; starred=true) ==
raw"\begin{align*}
\frac{dx}{dt} =& \frac{y}{c_{1}} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align*}
"

f = @ode_def TestRaw begin
    dx = y/c_1 - x
    dy = x^c_2 - y
end c_1 c_2
@test latexraw(f) == ["\\frac{dx}{dt} = \\frac{y}{c_{1}} - x", "\\frac{dy}{dt} = x^{c_{2}} - y"]

#==============================================================================#
#================================  Test align  ================================#
#==============================================================================#

f = @ode_def TestAlignFeedback begin
    dx = y/c_1 - x
    dy = x^c_2 - y
end c_1 c_2

g = @ode_def TestAlign2 begin
    dx = p_x - d_x*x - d_x_y*y*x
    dy = p_y - d_y*y
    dz = p_z
end p_x d_x d_x_y p_y d_y p_z

@test latexify(f) ==
raw"\begin{align}
\frac{dx}{dt} =& \frac{y}{c_{1}} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

@test latexify(f; rows=1) ==
raw"\begin{align}
\frac{dx}{dt} =& \frac{y}{c_{1}} - x
\end{align}
"

@test latexify(f; env=:align, field=:symjac) ==
raw"\begin{align}
\frac{dx}{dt} =& -1 \cdot 1 =& c_{1}^{-1} \\
\frac{dy}{dt} =& x^{-1 + c_{2}} \cdot c_{2} =& -1 \cdot 1
\end{align}
"

@test latexalign(["d$(x)/dt" for x = f.syms], f.funcs) ==
raw"\begin{align}
\frac{dx}{dt} =& \frac{y}{c_{1}} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

@test latexalign([f,g]) == "\\begin{align}\n\\frac{dx}{dt}  &=  \\frac{y}{c_{1}} - x  &  \\frac{dx}{dt}  &=  p_{x} - d_{x} \\cdot x - d_{x\\_y} \\cdot y \\cdot x  &  \\\\\n\\frac{dy}{dt}  &=  x^{c_{2}} - y  &  \\frac{dy}{dt}  &=  p_{y} - d_{y} \\cdot y  &  \\\\\n  &    &  \\frac{dz}{dt}  &=  p_{z} \\cdot 1  & \n\\end{align}\n"

@reaction_func hill2(x, v, k) = v*x^2/(k^2 + x^2)

rn = @reaction_network TestAlignMyRnType begin
  hill2(y, v_x, k_x), 0 --> x
  p_y, 0 --> y
  (d_x, d_y), (x, y) --> 0
  (r_b, r_u), x ↔ y
end v_x k_x p_y d_x d_y r_b r_u

@test latexalign(rn; clean=true).s ==
raw"\begin{align}
\frac{dx}{dt} =& \frac{v_{x} \cdot y^{2}}{k_{x}^{2} + y^{2}} - d_{x} \cdot x - r_{b} \cdot x + r_{u} \cdot y \\
\frac{dy}{dt} =& p_{y} - d_{y} \cdot y + r_{b} \cdot x - r_{u} \cdot y
\end{align}
"


text:
using Random
using GLPK
using MathOptInterface

const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities;

using Test

D = 10  # variable dimension
N = 20  # no of inequality constraints

s = rand(N)
s = 2*s.-1
λ = max.(-s, 0)
s = max.(s, 0)
x̂ = rand(D)
A = rand(N, D)
b = A*x̂ + s
c = -A' * λ

# can feed dual problem to optimizer like this:
# model = MOI.instantiate(dual_optimizer(GLPK.Optimizer), with_bridge_type=Float64)

model = GLPK.Optimizer()
x = MOI.add_variables(model, D)

# define objective
objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0)
MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

# will be useful later
constraint_indices = []

# set constraints
for i in 1:N
    push!(constraint_indices, MOI.add_constraint(model,MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:], x), 0.),MOI.LessThan(b[i])))
end

for i in 1:D
    push!(constraint_indices, MOI.add_constraint(model,MOI.SingleVariable(x[i]),MOI.GreaterThan(0.)))
end

MOI.optimize!(model)

@assert MOI.get(model, MOI.TerminationStatus()) in [MOI.LOCALLY_SOLVED, MOI.OPTIMAL]

x̄ = MOI.get(model, MOI.VariablePrimal(), x);  # solution

@assert abs(c'x̄ - c'x̂) <= 1e-8   # sanity check

##
## Checking KKT conditions
##

is_less_than(set::S) where {S<:MOI.AbstractSet} = false
is_less_than(set::MOI.LessThan{T}) where T = true

for con_index in constraint_indices
    con_value = MOI.get(model, MOI.ConstraintPrimal(), con_index)
    set = MOI.get(model, MOI.ConstraintSet(), con_index)
    μ = MOI.get(model, MOI.ConstraintDual(), con_index)

    if is_less_than(set)
        # μ[i]*(Ax - b)[i] = 0
        @assert μ*(con_value - set.upper) < 1e-10
    else
        # μ[j]*x[j] = 0
        @assert μ*(con_value - set.lower) < 1e-10
    end
end


text:
# Exercise 10-1
function nestedsum(t)
   return sum(sum.(t))
end

t = [[1, 2], [3], [4, 5, 6]]
# println(nestedsum(t))

# Exercise 10-2
function cumulsum(t)
    res = Array{Int64, 1}()
    for ele in t
        if ele == t[begin]
            push!(res, ele)
        else
            push!(res, ele + res[end])
        end
    end
    return res
end

function cumulsum!(t)
    for i in eachindex(t)
        if i == firstindex(t)
            continue
        else
            t[i] += t[prevind(t, i)]
        end
    end
    return t
end

t = [1, 2, 3, 4, 5, 6]
# println(cumulsum(t))
cumulsum!(t)
# println(t)

# Exercise 10-3
function interior(t)
    return t[begin+1:end-1]
end

t = [1, 2, 3, 4, 5, 6]
# println(interior(t))

# Exercise 10-4
function interior!(t)
    pop!(t)
    popfirst!(t)
end

interior!(t)
# println(t)

# Exercise 10-5
function issort(t)
    for i in eachindex(t)
        if i == firstindex(t)
            continue
        else
            if t[i] < t[prevind(t, i)]
                return false
            end
        end
    end
    return true
end

# println(issort([1, 2, 2, 3, 4, 5, 6]))
# println(issort(["abcdef", "bacdef", "cbadef"]))
# println(issort(["d", "bacdef", "cbadef"]))

# Exercise 10-6
function isanagram(a, b)
    return sort(collect(a)) == sort(collect(b))
end

# println(isanagram("abcdef", "acbdef"))
# println(isanagram("aaaa", "bbbb"))

# Exercise 10-7
function hasduplicates(t)
    return unique(t) != t
end

dup = [3, 4, 5, 1, 2, 3, 4, 5]
unq = [1, 2, 3, 4, 5]
# println(hasduplicates(dup))
# println(hasduplicates(unq))

# Exercise 10-8
function gensample()
    days = Array{Int64, 1}()
    for i = 1:23
        push!(days, rand(1:365))
    end
    return hasduplicates(days)
end

function ex10_8(n)
    count = 0
    for i = 1:n
        if gensample()
            count += 1
        end
    end
    return count / n
end

# println(ex10_8(10000))

# Exercise 10-9
# reader2 takes longer to run due to repeated creation and allocation of new arrays.
function reader1()
    fin = open("words.txt")
    words = Array{String, 1}()
    for line in eachline("words.txt")
        push!(words, line)
    end
    return words
end

function reader2()
    fin = open("words.txt")
    words = Array{String, 1}()
    for line in eachline("words.txt")
        words = [words..., line]
    end
    return words
end

function ex10_9()
    @time begin
        reader1()
    end
    @time begin
        reader2()
    end
end

# ex10_9()

# Exercise 10-10
function inbisect(arr, val)
    lo, hi = firstindex(arr), lastindex(arr)
    while lo <= hi
        mid = lo + (hi - lo) ÷ 2
        if arr[mid] < val
            lo = mid + 1
        elseif arr[mid] > val
            hi = mid - 1
        else
            return mid
        end
    end
    return -1
end

function ex10_10()
    arr = [1, 3, 5, 7, 9, 10, 12, 14, 17, 20]
    println(inbisect(arr, 17))
    println(inbisect(arr, 25))
    println(inbisect(arr, 11))
    println(inbisect(arr, 5))
    println(inbisect(arr, -5))
end

# ex10_10()

# Exercise 10-11
function inbisect(arr, val, lo, hi)
    while lo <= hi
        mid = lo + (hi - lo) ÷ 2
        if arr[mid] < val
            lo = mid + 1
        elseif arr[mid] > val
            hi = mid - 1
        else
            return mid
        end
    end
    return -1
end

function reversepairs()
    words = reader1()
    n = length(words)
    revpairs = Array{Tuple, 1}()
    for i = 1:n
        curr = words[i]
        index = inbisect(words, reverse(curr), i + 1, n)
        if index != -1
            push!(revpairs, (curr, words[index]))
        end
    end
    return revpairs
end

# println(reversepairs())

# Exercise 10-12
function interlock()
    words = reader1()
    n = length(words)
    interlocks = Array{Tuple, 1}()
    for i = 1:n
        curr = words[i]
        odd, even = curr[begin:2:end], curr[begin+1:2:end]
        if inbisect(words, odd) != -1 && inbisect(words, even) != -1
            push!(interlocks, (odd, even))
        end
    end
    return interlocks
end

# println(interlock())

function interlock3()
    words = reader1()
    n = length(words)
    interlocks = Array{Tuple, 1}()
    for i = 1:n
        curr = words[i]
        a, b, c = curr[begin:3:end], curr[begin+1:3:end], curr[begin+2:3:end]
        if inbisect(words, a) != -1 && inbisect(words, b) != -1 && inbisect(words, c) != -1
            push!(interlocks, (a, b, c))
        end
    end
    return interlocks
end

# println(interlock3())