text:
"""
Main module for Starlight.jl - a greedy framework for greedy developers.

  # Reexports

  - Base: ReentrantLock, lock, unlock
  - FileIO: load
  - Pkg.Artifacts
  - LazyArtifacts
  - DataStructures: Queue, enqueue!, dequeue!
  - DataFrames
  - Colors
  - SimpleDirectMediaLayer
  - SimpleDirectMediaLayer.LibSDL2
  - Telescope
"""
module Starlight

using Reexport
@reexport using Base: ReentrantLock, lock, unlock
@reexport using FileIO: load
@reexport using Pkg.Artifacts
@reexport using LazyArtifacts
@reexport using DataStructures: Queue, enqueue!, dequeue!
@reexport using DataFrames
@reexport using Colors
@reexport using SimpleDirectMediaLayer
@reexport using SimpleDirectMediaLayer.LibSDL2
@reexport using Telescope

include("Logging.jl")
using Starlight.Log

export handleMessage!, sendMessage, listenFor, unlistenFrom, handleException, dispatchMessage
export App, awake!, shutdown!, system!, run!, on, off
export clk, ecs, inp, ts, phys, scn
export listeners, messages, listener_lock
export systemAwakeOrder, systemShutdownOrder

"""
    const listeners = Dict{DataType, Set{Any}}()

Internal dictionary that pairs event types with their listeners.
"""
const listeners = Dict{DataType, Set{Any}}()

"""
    const messages = Channel(Inf)

Internal Channel used to hold messages.
"""
const messages = Channel(Inf)

"""
    const listener_lock = ReentrantLock()

Internal lock used to synchronize access to the listeners dictionary.
"""
const listener_lock = ReentrantLock()

"""
```julia
handleMessage!(l, m) = nothing
```

Invoked on listener l when message m is received.

Currently there are three ways this function can be invoked for a given listener and message:
  - automatically by the message dispatcher, when the listener has been registered with listenFor and the message is of the corresponding type
  - automatically by the Physics subsystem on two colliding objects with message type TS_CollisionEvent
  - manually by user code (discouraged)

By default does (literally) nothing, but methods can be added for any type combination.
"""
handleMessage!(l, m) = nothing

"""
```julia
function sendMessage(m)
  if haskey(listeners, typeof(m))
    put!(messages, m)
  end
end
```
Add a message to the message queue if the message's type has any listeners. Otherwise drop.

The preferred way to send events/messages in Starlight.
"""
function sendMessage(m)
  # drop if no one's listening
  if haskey(listeners, typeof(m))
    put!(messages, m)
  end
end

"""
```julia
function listenFor(e::Any, d::DataType)
  lock(listener_lock)
  if !haskey(listeners, d) listeners[d] = Set{Any}() end
  push!(listeners[d], e)
  unlock(listener_lock)
end
```

Add a listener e for messages of type d. l's handleMessage! method will be invoked when messages of type d are received.
"""
function listenFor(e::Any, d::DataType)
  lock(listener_lock)
  if !haskey(listeners, d) listeners[d] = Set{Any}() end
  push!(listeners[d], e)
  unlock(listener_lock)
end

"""
```julia
function unlistenFrom(e::Any, d::DataType)
  lock(listener_lock)
  if haskey(listeners, d) delete!(listeners[d], e) end
  unlock(listener_lock)
end
```
Stop listener e from receiving messages of type d.
"""
function unlistenFrom(e::Any, d::DataType)
  lock(listener_lock)
  if haskey(listeners, d) delete!(listeners[d], e) end
  unlock(listener_lock)
end

"""
`handleException() -> Nothing`

Print a stacktrace and exit safely.
"""
function handleException() ::Nothing

    exception_stack = current_exceptions(backtrace=true)
    if isempty(exception_stack)
        return
    end

    exception = first(first(exception_stack)) # sic

    # non-fatal exceptions go here, currently all exceptions are fatal
    non_fatal_exception_types = Type[]
    if typeof(exception) in non_fatal_exception_types
        println(eval(exception))
        Log.@warning "a non-fatal exception occurred: " * string(exception)
        return
    end

    Log.@error "an exception occurred: " * string(exception)
    Log.@error "aborting..."
    shutdown!(App())
    throw(exception_stack)
    return nothing
end

"""
```julia
function dispatchMessage(arg)
  try
    m = take!(messages) # NOTE messages are fully processed in the order they are received
    d = typeof(m)
    if haskey(listeners, d)
      # Threads.@threads doesn't work on raw sets
      # because it needs to use indexing to split
      # up memory, i work around it this way
      Threads.@threads for l in Vector([listeners[d]...])
        handleMessage!(l, m)
      end
    end
  catch
    handleException()
  end
end
```

Take a single message from the event queue and invoke handleMessage! for all listeners in parallel.

If there are no listeners for the message type (i.e. they were removed after the message was sent), drop.

Normally runs in an infinite loop in a background Task started by awake!(::App).

!!! danger

    This function is called internally by Starlight and documented here
    for completeness. Do not invoke it yourself unless you know what you
    are doing.
"""
function dispatchMessage(arg)
  Log.@debug "dispatchMessage"
  try
    m = take!(messages) # NOTE messages are fully processed in the order they are received
    d = typeof(m)
    Log.@debug "dequeued message $(m) of type $(d)"
    lock(listener_lock)
    if haskey(listeners, d)
      # Threads.@threads doesn't work on raw sets
      # because it needs to use indexing to split
      # up memory, i work around it this way
      Threads.@threads for l in Vector([listeners[d]...])
        handleMessage!(l, m)
      end
    end
    unlock(listener_lock)
  catch
    handleException()
  end
end

"""
```julia
awake!(a) = nothing
```

Arbitrary startup function.

By default does literally nothing, but methods can be added for any type.

The preferred way to call listenFor is from inside an awake! method.
"""
awake!(a) = nothing

"""
```julia
shutdown!(a) = nothing
```

Arbitrary shutdown function.

By default does literally nothing, but methods can be added for any type.

The preferred way to call unlistenFrom is from inside a shutdown! method.
"""
shutdown!(a) = nothing

include("Clock.jl")
include("ECS.jl")
include("TS.jl")
include("Entities.jl")
include("Input.jl")
include("Physics.jl")

"""
Struct for the "master App" singleton.

Constructor accepts keyword arguments for all fields except systems, and no positional arguments.

# Fields
  - systems::Dict{DataType, Any} Dictionary of subsystems by type
  - running::Bool Whether awake! has been called
  - wdth::Int Window width
  - hght::Int Window height
  - bgrd::Colorant Window default background color
"""
mutable struct App
  systems::Dict{DataType, Any}
  running::Bool
  wdth::Int
  hght::Int
  bgrd::Colorant
  function App(; wdth::Int=400, hght::Int=400, bgrd=colorant"grey")
    # singleton pattern from Tom Kwong's book
    global app
    global app_lock
    lock(app_lock)
    try
      if !isassigned(app)
        a = new(Dict{DataType, Any}(), false, wdth, hght, bgrd)

        system!(a, Clock())
        system!(a, ECS())
        system!(a, Scene())
        system!(a, TS())
        system!(a, Input())
        system!(a, Physics())

        app[] = finalizer(shutdown!, a)

        # root pid is 0 (default) indicating "here and no further",
        # always needed, also it will technically parent of itself.
        # mutate at your own peril, but remember that a  user can
        # define update!(r::Root, Δ::AbstractFloat) if they want
        instantiate!(Root())
      end

    catch e
      rethrow()
    finally
      unlock(app_lock)
      return app[]
    end
  end
end

const app = Ref{App}()
const app_lock = ReentrantLock()

# TODO: fix this, it feels hacky
"""
    clk() = App().systems[Clock]

Get the current Clock.
"""
clk() = App().systems[Clock]

"""
    scn() = App().systems[Scene]

Get the current Scene graph.
"""
scn() = App().systems[Scene]

"""
    ts() = App().systems[TS]

Get the current Telescope backend subsystem.
"""
ts() = App().systems[TS]

"""
    inp() = App().systems[Input]

Get the current Input manager.
"""
inp() = App().systems[Input]

"""
    ecs() = App().systems[ECS]

Get the current Entity Component System.
"""
ecs() = App().systems[ECS]

"""
    phys() = App().systems[Physics]

Get the current Physics manager.
"""
phys() = App().systems[Physics]

"""
    on(a::App) = a.running

Check whether awake! has been called.
"""
on(a::App) = a.running

"""
    off(a::App) = !a.running

Opposite of on.
"""
off(a::App) = !a.running

# TODO: fix this with an ordered dict or something

"""
```julia
systemAwakeOrder = () -> [clk(), ts(), inp(), phys(), ecs(), scn()]
```

The order in which App subsystems should be awoken in order to not produce bugs.
"""
systemAwakeOrder = () -> [clk(), ts(), inp(), phys(), ecs(), scn()]

"""
```julia
systemShutdownOrder = () -> [clk(), inp(), phys(), scn(), ecs(), ts()]
```

The order in which App subsystems should be shut down in order to not produce bugs.
"""
systemShutdownOrder = () -> [clk(), inp(), phys(), scn(), ecs(), ts()]

# TODO: need an elegant way to remove systems

"""
```julia
system!(a::App, s) = a.systems[typeof(s)] = s
```

Add something to an App's systems dictionary.
"""
system!(a::App, s) = a.systems[typeof(s)] = s

"""
```julia
function awake!(a::App)
  if !on(a)
    job!(clk(), dispatchMessage)
    map(awake!, systemAwakeOrder())
    a.running = true
  end
end
```

If not on, start the message dispatcher call awake! on all subsystems.

Note that if running from a script the app will still exit when Julia exits, it will never block. Figuring out whether/how to keep it alive is on the user. One method is to use run!, see below.
"""
function awake!(a::App)
  if !on(a)

    Log.init() #Base.Filesystem.pwd() * "starlight.log")
    # Log.set_debug_enabled(true)

    job!(clk(), dispatchMessage) # this could be parallelized if not for mqueue_lock
    map(awake!, systemAwakeOrder())
    a.running = true
  end
end

"""
```julia
function shutdown!(a::App)
  if !off(a)
    map(shutdown!, systemShutdownOrder())
    a.running = false
  end
end
```

If not off, call shutdown! on all subsystems.
"""
function shutdown!(a::App)
  if !off(a)
    map(shutdown!, systemShutdownOrder())
    a.running = false
    Log.quit()
  end
end

"""
```julia
function run!(a::App)
  awake!(a)
  if !isinteractive()
    while on(a)
      yield()
    end
  end
end
```

Call awake! and keep alive until switched off.
"""
function run!(a::App)
  awake!(a)
  if !isinteractive()
    while on(a)
      yield()
    end
  end
end

end


text:
"""
AK135 evaluates the seismic velocities, density and attenuation for the
1D Earth model AK135.
"""
module AK135

export
	rho,
	vp,
	vs

"Number of radial points in AK135 model (number of layers + 1)"
const N = 136
"Radius of Earth in AK135 model"
const R = 6371.

const rad = [
       0.,   50.,  101.,  152.,  202.,  253.,  304.,  354.,
     405.,  456.,  507.,  557.,  659.,  710.,  760.,  811.,
     862.,  912.,  963., 1014., 1064., 1115., 1166., 1217.,
    1217., 1267., 1317., 1368., 1418., 1468., 1519., 1569.,
    1670., 1720., 1770., 1821., 1871., 1921., 1972., 2022.,
    2072., 2123., 2173., 2223., 2274., 2324., 2374., 2425.,
    2475., 2525., 2576., 2626., 2676., 2727., 2777., 2827.,
    2878., 2928., 2978., 3029., 3079., 3129., 3180., 3230.,
    3280., 3331., 3381., 3431., 3479., 3479., 3531., 3581.,
    3631., 3631., 3681., 3731., 3779., 3829., 3878., 3928.,
    3977., 4027., 4076., 4126., 4175., 4225., 4274., 4324.,
    4373., 4423., 4472., 4522., 4571., 4621., 4670., 4720.,
    4769., 4819., 4868., 4918., 4967., 5017., 5066., 5116.,
    5165., 5215., 5264., 5314., 5363., 5413., 5462., 5512.,
    5561., 5611., 5661., 5711., 5711., 5761., 5811., 5861.,
    5911., 5961., 5961., 6011., 6061., 6111., 6161., 6161.,
    6206., 6251., 6293., 6336., 6336., 6351., 6351., 6371.] # / km

# Indices of nodes in each major layer
const irad_inner_core = 1:24
const irad_outer_core = 25:69
const irad_mantle = 70:136

const rho_model = [
    13.0122, 13.0117, 13.0100, 13.0074, 13.0036, 12.9988, 12.9929, 12.9859,
    12.9779, 12.9688, 12.9586, 12.9474, 12.9217, 12.9070, 12.8917, 12.8751,
    12.8574, 12.8387, 12.8188, 12.7980, 12.7760, 12.7530, 12.7289, 12.7037,
    12.1391, 12.1133, 12.0867, 12.0593, 12.0311, 12.0001, 11.9722, 11.9414,
    11.8772, 11.8437, 11.8092, 11.7737, 11.7373, 11.6998, 11.6612, 11.6216,
    11.5809, 11.5391, 11.4962, 11.4521, 11.4069, 11.3604, 11.3127, 11.2639,
    11.2137, 11.1623, 11.1095, 11.0555, 11.0001, 10.9434, 10.8852, 10.8257,
    10.7647, 10.7023, 10.6385, 10.5731, 10.5062, 10.4378, 10.3679, 10.2964,
    10.2233, 10.1485, 10.0722,  9.9942,  9.9145,  5.7721,  5.7458,  5.7196,
     5.6934,  5.4387,  5.4176,  5.3962,  5.3748,  5.3531,  5.3313,  5.3092,
     5.2870,  5.2646,  5.2420,  5.2192,  5.1963,  5.1732,  5.1499,  5.1264,
     5.1027,  5.0789,  5.0548,  5.0306,  5.0062,  4.9817,  4.9570,  4.9321,
     4.9069,  4.8817,  4.8562,  4.8307,  4.8050,  4.7790,  4.7528,  4.7266,
     4.7001,  4.6735,  4.6467,  4.6198,  4.5926,  4.5654,  4.5162,  4.4650,
     4.4118,  4.3565,  4.2986,  4.2387,  3.9201,  3.9206,  3.9218,  3.9233,
     3.9273,  3.9317,  3.5068,  3.4577,  3.4110,  3.3663,  3.3243,  3.3243,
     3.3711,  3.4268,  3.3450,  3.3200,  2.9200,  2.9200,  2.7200,  2.7200] # g/cm^3

const vp_model = [
    11.2622, 11.2618, 11.2606, 11.2586, 11.2557, 11.2521, 11.2477, 11.2424,
    11.2364, 11.2295, 11.2219, 11.2134, 11.1941, 11.1830, 11.1715, 11.1590,
    11.1457, 11.1316, 11.1166, 11.0983, 11.0850, 11.0718, 11.0585, 11.0427,
    10.2890, 10.2854, 10.2745, 10.2565, 10.2329, 10.2049, 10.1739, 10.1415,
    10.0768, 10.0439, 10.0103,  9.9761,  9.9410,  9.9051,  9.8682,  9.8304,
     9.7914,  9.7513,  9.7100,  9.6673,  9.6232,  9.5777,  9.5306,  9.4814,
     9.4297,  9.3760,  9.3205,  9.2634,  9.2042,  9.1426,  9.0792,  9.0138,
     8.9461,  8.8761,  8.8036,  8.7283,  8.6496,  8.5692,  8.4861,  8.4001,
     8.3122,  8.2213,  8.1283,  8.0382,  8.0000, 13.6601, 13.6570, 13.6533,
    13.6498, 13.6498, 13.5899, 13.5311, 13.4741, 13.4156, 13.3584, 13.3017,
    13.2465, 13.1895, 13.1337, 13.0786, 13.0226, 12.9663, 12.9093, 12.8524,
    12.7956, 12.7384, 12.6807, 12.6226, 12.5638, 12.5030, 12.4427, 12.3813,
    12.3181, 12.2558, 12.1912, 12.1247, 12.0571, 11.9891, 11.9208, 11.8491,
    11.7768, 11.7020, 11.6265, 11.5493, 11.4704, 11.3897, 11.3068, 11.2228,
    11.1355, 11.0553, 10.9222, 10.7909, 10.2000, 10.0320,  9.8640,  9.6962,
     9.5280,  9.3601,  9.0302,  8.8476,  8.6650,  8.4822,  8.3007,  8.3007,
     8.1750,  8.0505,  8.0450,  8.0400,  6.5000,  6.5000,  5.8000,  5.8000] # km/s

const vs_model = [
    3.6678, 3.6675, 3.6667, 3.6653, 3.6633, 3.6608, 3.6577, 3.6540,
    3.6498, 3.6450, 3.6396, 3.6337, 3.6202, 3.6130, 3.6044, 3.5957,
    3.5864, 3.5765, 3.5661, 3.5551, 3.5435, 3.5314, 3.5187, 3.5043,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 7.2817, 7.2700, 7.2593,
    7.2485, 7.2485, 7.2253, 7.2031, 7.1804, 7.1584, 7.1368, 7.1144,
    7.0932, 7.0722, 7.0504, 7.0286, 7.0069, 6.9852, 6.9625, 6.9416,
    6.9194, 6.8972, 6.8743, 6.8517, 6.8289, 6.8056, 6.7820, 6.7579,
    6.7323, 6.7070, 6.6813, 6.6554, 6.6285, 6.6009, 6.5728, 6.5431,
    6.5131, 6.4822, 6.4514, 6.4182, 6.3860, 6.3519, 6.3164, 6.2799,
    6.2424, 6.2100, 6.0898, 5.9607, 5.6104, 5.5047, 5.3989, 5.2922,
    5.1864, 5.0806, 4.8702, 4.7832, 4.6964, 4.6094, 4.5184, 4.5184,
    4.5090, 4.5000, 4.4900, 4.4800, 3.8500, 3.8500, 3.4600, 3.4600] # km/s

const Qkappa_model = [
      601.27,   601.32,   601.46,   601.70,   602.05,   602.49,   603.04,   603.69,
      604.44,   605.28,   606.26,   607.31,   609.74,   611.18,   612.62,   614.21,
      615.93,   617.78,   619.71,   621.50,   624.08,   626.87,   629.89,   633.26,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00, 57822.00,
    57822.00, 57822.00, 57822.00, 57822.00, 57822.00,   723.12,   725.11,   726.87,
      722.73,   933.21,   940.88,   952.00,   960.36,   968.46,   976.81,   985.63,
      990.77,   999.44,  1008.79,  1018.38,  1032.14,  1042.07,  1048.09,  1058.03,
     1064.23,  1070.38,  1085.97,  1097.16,  1108.58,  1120.09,  1127.02,  1134.01,
     1141.32,  1148.76,  1156.04,  1163.16,  1170.53,  1178.19,  1186.06,  1193.99,
     1202.04,  1210.02,  1217.91,  1226.52,  1234.54,  1243.02,  1251.69,  1260.68,
     1269.44,  1277.93,  1311.17,  1350.54,   428.69,   425.51,   422.55,   419.94,
      417.32,   413.66,   377.93,   366.34,   355.85,   346.37,   338.47,   200.97,
      188.72,   182.57,   182.03,   182.03,   972.77,   972.77,  1368.02,  1368.02]

const Qmu_model = [
     85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,
     85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,
     85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,  85.03,
      0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,
      0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,
      0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,
      0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,
      0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.00,
      0.00,   0.00,   0.00,   0.00,   0.00, 273.97, 273.97, 273.97,
    271.74, 350.88, 354.61, 359.71, 363.64, 367.65, 371.75, 375.94,
    378.79, 383.14, 387.60, 392.16, 398.41, 403.23, 406.50, 411.52,
    414.94, 418.41, 425.53, 431.03, 436.68, 442.48, 446.43, 450.45,
    454.55, 458.72, 462.96, 467.29, 471.70, 476.19, 480.77, 485.44,
    490.20, 495.05, 500.00, 505.05, 510.20, 515.46, 520.83, 526.32,
    531.91, 537.63, 543.48, 549.45, 172.93, 170.82, 168.78, 166.80,
    164.87, 162.50, 146.57, 142.76, 139.38, 136.38, 133.72,  79.40,
     76.55,  76.06,  75.60,  75.60, 403.93, 403.93, 599.99, 599.99]

# Generate functions to evaluate properties at given radius
for (v, doc_name, doc_unit) in zip((:rho, :vp, :vs), ("rho", "Vp", "Vs"), ("g/cm^3", "km/s", "km/s"))
	func_name, var_name = symbol("$v"), symbol("$(v)_model")
	@eval begin
		@doc """
		`$($func_name)(r) -> $($doc_name)`

		Return the value of $($doc_name) in the AK135 model at radius `r` km ($($doc_unit)).
		""" ->
		($func_name)(r) = interp(r, $var_name)
	end
end
for (v, doc_name) in zip((:Qkappa, :Qmu), ("Qkappa", "Qmu"))
	func_name, var_name = symbol("$v"), symbol("$(v)_model")
	@eval begin
		@doc """
		`$($func_name)(r) -> $($doc_name)`

		Return the value of $($doc_name) in the AK135 model at radius `r` km.
		""" ->
		($func_name)(r) = $var_name[findlayer(r)]
	end
end

"""
`interp(r, v) -> V`

Return the interpolated value `V` from the array `v` at radius `r` km,
providing the array `v` is the same length as the number of layers in
AK135 (136, confusingly).
"""
function interp(r, v)
	length(v) == N || error("AK135.interp: Size of `v` is not $(N)")
	# i is the index of the node *below* this radius
	i = findlayer(r)
	# Exactly at a node or at the surface
	if rad[i] == r || i == N
		return v[i]
	# Discontinuity: choose the lower layer
	elseif rad[i] == rad[i+1]
		return v[i]
	else
		return (v[i]*(rad[i+1] - r) + v[i+1]*(r - rad[i]))/(rad[i+1] - rad[i])
	end
end

"""
`findlayer(r) -> i`

Return the node index `i` at or above which the radius `r` occurs.
"""
function findlayer(r)
	rad[1] <= r <= rad[end] ||
		error("AK135.findlayer: Radius must be in range $(rad[1]) to $(rad[end]) km")
	@inbounds for i = 1:N
		r < rad[i] && return i
	end
	N
end

end # module

ext: jl

max_stars_repo_name: anowacki/assorted-julia-modules

max_issues_repo_name: anowacki/assorted-julia-modules

max_forks_repo_name: anowacki/assorted-julia-modules


text:
using Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Levels" begin

@testset "floor pivots" begin
    # values verified by various website calculators
    @test isapprox(values(floorpivots(ohlc)[:r3])[1]   , 123.310, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:r2])[1]   , 119.52, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:r1])[1]   , 115.73, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:pivot])[1], 108.71, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:s1])[1]   , 104.92, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:s2])[1]   , 97.900, atol=.01)
    @test isapprox(values(floorpivots(ohlc)[:s3])[1]   , 94.110, atol=.01)
    @test timestamp(floorpivots(ohlc))[end] == Date(2001,12,31)
end

@testset "woodiespivots" begin
    # @test_approx_eq  97.37500000000001 value(wr4)[2]   # values NEED to be verified with online calculators
    # @test_approx_eq   88.62500000000001 value(ws4)[2]
    @test isapprox(values(woodiespivots(ohlc)[:r3])[1]   , 124.465, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:r2])[1]   , 118.480, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:r1])[1]   , 113.655, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:pivot])[1], 107.670, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:s1])[1]   , 102.845, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:s2])[1]   , 96.8625, atol=.01)
    @test isapprox(values(woodiespivots(ohlc)[:s3])[1]   , 92.035, atol=.01)
    @test timestamp(woodiespivots(ohlc))[end] == Date(2001,12,31)
end

end  # @testset "Levels"


text:
# Randomized ICA

function rica(X):

end


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_CG_smearing.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("calc_Focc.jl")
include("smear_FD.jl")
include("calc_entropy.jl")
include("sum_upto_E_fermi.jl")
include("calc_grad_Haux.jl")
include("andersonmix.jl")

include("printMatrix.jl")

function test_main( Ns )

    Random.seed!(1234)

    a = 5.66/0.52917721 # Lattice constant (converted from angstroms to bohrs)
    LatVecs = a*diagm(0 => ones(3))

    pw = PWGrid( Ns, LatVecs )

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("Ns = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # diamond lattice in cubic cell
    Xpos = zeros(3,1)
    Xpos[:,1] = [0.0, 0.0, 0.0]

    Nspecies = 1
    atmsymb = ["Ge"] # unique list of atomic symbols
    atm2species = [1]    # mapping from atom to species
    Zv = [4.0]    # only valence ?

    Sf = calc_strfact( Xpos, Nspecies, atm2species, pw.gvec.G )

    E_nn = calc_ewald( pw, Sf, Xpos, Nspecies, atm2species, Zv )
    @printf("E_nn = %18.10f\n", E_nn)

    Vps = zeros(ComplexF64, Npoints)

    # Ge pseudopotential
    Z = Zv[1]
    λ = 18.5
    rc = 1.052

    Vps = zeros(ComplexF64,Npoints)
    for i = 2:Npoints
        Gm = sqrt(G2[i])
        Vps[i] = -2*pi*exp(-pi*Gm/λ)*cos(Gm*rc)*(Gm/λ)/(1 - exp(-2*pi*Gm/λ))
        for n = 0:4
            Vps[i] = Vps[i] + (-1.0)^n*exp(-λ*rc*n)/(1 + (n*λ/Gm)^2)
        end
        Vps[i] = Vps[i]*4*pi*Z/Gm^2*(1+exp(-λ*rc)) - 4*pi*Z/Gm^2
    end

    n = collect(1:4)
    Vps[1] = 4*pi*Z*(1+exp(-λ*rc))*(rc^2/2+1/λ^2 *
                     (pi^2/6 + sum((-1).^n.*exp.(-λ*rc*n)./n.^2)))

    Vps[:] = Vps[:]/Ω

    V_ionic = real( G_to_R(Ns, Vps .* Sf) ) * Npoints

    # needed to sum up over Nspecies for more than one species
    V_ionic = reshape( V_ionic, (Npoints) )

    Nstates = 7
    # Initial Focc
    Focc = zeros(Nstates)
    Focc[1] = 2.0
    Focc[2] = 1.0
    Focc[3] = 1.0
    Nelectrons = 4

    Energies, Potentials, psi, evals = KS_solve_Emin_CG_smearing( pw, V_ionic, Focc, Nstates, E_NN=E_nn )

end

test_main( [24,24,24] )


text:
mutable struct EnemySets
    var::Vector{Symbol}
    names::Vector{String}
    sets::Vector{Vector{Int64}}
    EnemySets(var, names, sets) = new(var, names, sets)
    EnemySets() = new(Symbol[], String[], Vector{Vector{Int64}}(undef, 0))
end


text:
import PowerModelsACDC
const _PMACDC = PowerModelsACDC
import PowerModels
const _PM = PowerModels
import Ipopt
import Memento
import JuMP
import Gurobi  # needs startvalues for all variables!
import Juniper
import FlexPlan
const _FP = FlexPlan

file = normpath(@__DIR__,"../../test/data/case67/case67.m")
file_2030 = normpath(@__DIR__,"../../test/data/case67/case67_tnep_2030.m")
file_2040 = normpath(@__DIR__,"../../test/data/case67/case67_tnep_2040.m")
file_2050 = normpath(@__DIR__,"../../test/data/case67/case67_tnep_2050.m")

ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0)
gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer)
juniper = JuMP.optimizer_with_attributes(Juniper.Optimizer)

s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)

resultAC = _PMACDC.run_acdcopf(file, _PM.ACPPowerModel, ipopt; setting = s)
resultDC = _PMACDC.run_acdcopf(file, _PM.DCPPowerModel, gurobi; setting = s)

result2030= _PMACDC.run_tnepopf(file_2030, _PM.DCPPowerModel, gurobi, setting = s)
result2040= _PMACDC.run_tnepopf(file_2040, _PM.DCPPowerModel, gurobi, setting = s)
result2050= _PMACDC.run_tnepopf(file_2050, _PM.DCPPowerModel, gurobi, setting = s)


text:
"""
 Holds the tableau of a weak explicit Runge-Kutta method.

   According to Andreas Rossler, "Second order Runge-Kutta methods for Stratonovich stochastic differential equations",
   BIT Numerical Mathematics (2007) 47, equation (5.1)
"""
struct TableauWERK{T} <: AbstractTableauERK{T}
    name::Symbol
    s::Int

    qdrift0::CoefficientsRK{T}
    qdrift1::CoefficientsRK{T}
    qdrift2::CoefficientsRK{T}

    qdiff0::CoefficientsRK{T}
    qdiff1::CoefficientsRK{T}
    qdiff2::CoefficientsRK{T}
    qdiff3::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift, qdiff are irrelevant and set to 0
    #
    # qdrift0, qdrift1, qdrift2 correspond to A0, A1, A2 in the paper
    # qdiff0, qdiff1, qdiff2, qdiff3 correspond to B0, B1, B2, B3
    # qdrift0.b = alpha
    # qdiff0.b  = beta1
    # qdiff3.b  = beta2
    # qdrift0.c = c0
    # qdrift1.c = c1
    # qdrift2.c = c2

    function TableauWERK{T}(name, qdrift0, qdrift1, qdrift2, qdiff0, qdiff1, qdiff2, qdiff3) where {T}
        @assert qdrift0.s == qdrift1.s == qdrift2.s == qdiff0.s == qdiff1.s == qdiff2.s == qdiff3.s
        @assert qdrift0.c[1] == qdrift1.c[1] == qdrift2.c[1] == qdiff0.c[1] == qdiff1.c[1] == qdiff2.c[1] == qdiff3.c[1] == 0.
        @assert istrilstrict(qdrift0.a)
        @assert istrilstrict(qdrift1.a)
        @assert istrilstrict(qdrift2.a)
        @assert istrilstrict(qdiff0.a)
        @assert istrilstrict(qdiff1.a)
        @assert istrilstrict(qdiff2.a)
        @assert istrilstrict(qdiff3.a)
        @assert !(qdrift0.s==1. && qdrift0.a[1,1] ≠ 0.)
        @assert !(qdrift1.s==1. && qdrift1.a[1,1] ≠ 0.)
        @assert !(qdrift2.s==1. && qdrift2.a[1,1] ≠ 0.)
        @assert !(qdiff0.s==1. && qdiff0.a[1,1] ≠ 0.)
        @assert !(qdiff1.s==1. && qdiff1.a[1,1] ≠ 0.)
        @assert !(qdiff2.s==1. && qdiff2.a[1,1] ≠ 0.)
        @assert !(qdiff3.s==1. && qdiff3.a[1,1] ≠ 0.)
        new(name, qdrift0.s, qdrift0, qdrift1, qdrift2, qdiff0, qdiff1, qdiff2, qdiff3)
    end
end

function TableauWERK(name, qdrift0::CoefficientsRK{T}, qdrift1::CoefficientsRK{T}, qdrift2::CoefficientsRK{T},
                            qdiff0::CoefficientsRK{T}, qdiff1::CoefficientsRK{T}, qdiff2::CoefficientsRK{T}, qdiff3::CoefficientsRK{T}) where {T}
    TableauWERK{T}(name, qdrift0, qdrift1, qdrift2, qdiff0, qdiff1, qdiff2, qdiff3)
end

function TableauWERK(name::Symbol, A0::Matrix{T}, A1::Matrix{T}, A2::Matrix{T},
                                   B0::Matrix{T}, B1::Matrix{T}, B2::Matrix{T}, B3::Matrix{T},
                                    α::Vector{T}, β1::Vector{T}, β2::Vector{T},
                                   c0::Vector{T}, c1::Vector{T}, c2::Vector{T} ) where {T}
    TableauWERK{T}(name, CoefficientsRK(name, 0, A0, α, c0), CoefficientsRK(name, 0, A1, α, c1), CoefficientsRK(name, 0, A2, α, c2),
                         CoefficientsRK(name, 0, B0, β1, c0), CoefficientsRK(name, 0, B1, β1, c1), CoefficientsRK(name, 0, B2, β1, c2), CoefficientsRK(name, 0, B3, β2, c1))
end

"Stochastic Explicit Runge-Kutta integrator."
struct IntegratorWERK{DT,TT,FT} <: StochasticIntegrator{DT,TT}
    equation::SDE{DT,TT,FT}
    tableau::TableauWERK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}

    q ::Matrix{Vector{DT}}
    Q0::Vector{DT}           # Q0[1..n]      - the internal stage H^(0)_i for a given i
    Q1::Vector{Vector{DT}}   # Q1[1..m][1..n] - the internal stages H^(l)_i for a given i
    Q2::Vector{Vector{DT}}   # Q2[1..m][1..n] - the internal stages \hat H^(l)_i for a given i
    V ::Vector{Vector{DT}}
    B1::Vector{Matrix{DT}}   # B1[i] holds the values of the diffusion matrix such that B1[i][:,l] is evaluated at H^(l)_i
    B2::Vector{Matrix{DT}}   # B2[i] holds the values of the diffusion matrix such that B2[i][:,l] is evaluated at \hat H^(l)_i
    tB::Vector{DT}           # the value of the l-th column of the diffusion term evaluated at the internal stage H^(l)_i for given i and l

    function IntegratorWERK{DT,TT,FT}(equation, tableau, Δt) where {DT,TT,FT}
        D = equation.d
        M = equation.m
        NS= equation.ns
        NI= equation.n
        S = tableau.s

        # create solution vectors
        q = create_solution_vector(DT, D, NS, NI)

        # create internal stage vectors
        Q0 = zeros(DT, D)
        Q1 = create_internal_stage_vector(DT, D, M)
        Q2 = create_internal_stage_vector(DT, D, M)
        V  = create_internal_stage_vector(DT, D, S)
        B1 = create_internal_stage_vector(DT, D, M, S)
        B2 = create_internal_stage_vector(DT, D, M, S)

        new(equation, tableau, Δt, zeros(DT,M), zeros(DT,M), q, Q0, Q1, Q2, V, B1, B2, zeros(DT,D))
    end
end

function IntegratorWERK(equation::SDE{DT,TT,FT}, tableau::TableauWERK{TT}, Δt::TT) where {DT,TT,FT}
    IntegratorWERK{DT,TT,FT}(equation, tableau, Δt)
end

function initialize!(int::IntegratorWERK, sol::SolutionSDE, k::Int, m::Int)
    check_solution_dimension_asserts(sol, k, m)

    # copy initial conditions from solution
    get_initial_conditions!(sol, int.q[k,m], k, m)
end

"""
Integrate SDE with explicit Runge-Kutta integrator.
 Calculating the n-th time step of the explicit integrator for the sample path r and the initial condition m
"""
function integrate_step!(int::IntegratorWERK{DT,TT,FT}, sol::SolutionSDE{DT,TT,NQ,NW}, r::Int, m::Int, n::Int) where {DT,TT,FT,NQ,NW}
    local tᵢ::TT
    local ydrift::DT
    local ydiff  = zeros(DT,sol.nm)
    local ydiff2::DT

    # copy the random variables \hat I and \tilde I representing the Wiener process
    if NW==1
        #1D Brownian motion, 1 sample path
        int.ΔW[1] = sol.W.ΔW[n-1]
        int.ΔZ[1] = sol.W.ΔZ[n-1]
    elseif NW==2
        #Multidimensional Brownian motion, 1 sample path
        int.ΔW .= sol.W.ΔW[n-1]
        int.ΔZ .= sol.W.ΔZ[n-1]
    elseif NW==3
        #1D or Multidimensional Brownian motion, r-th sample path
        int.ΔW .= sol.W.ΔW[n-1,r]
        int.ΔZ .= sol.W.ΔZ[n-1,r]
    end

    # THE FIRST INTERNAL STAGES ARE ALL EQUAL TO THE PREVIOUS STEP SOLUTION
    # calculates v(t,tQ0) and assigns to the 1st column of V
    int.equation.v(sol.t[0] + (n-1)*int.Δt, int.q[r,m], int.V[1])
    # calculates B(t,Q) and assigns to the matrix B1[1]
    # no need to calculate the columns of B separately, because all internal stages
    # for i=1 are equal to int.q[r,m]
    int.equation.B(sol.t[0] + (n-1)*int.Δt, int.q[r,m], int.B1[1])
    # copy B1[i] to B2[i] (they're equal for i=1)
    int.B2[1] .= int.B1[1]

    for i in 2:int.tableau.s

        # Calculating the internal stage H^(0)_i
        @inbounds for k in eachindex(int.Q0)
            # contribution from the drift part
            ydrift = 0.
            for j = 1:i-1
                ydrift += int.tableau.qdrift0.a[i,j] * int.V[j][k]
            end

            # ΔW contribution from the diffusion part
            ydiff .= zeros(DT, sol.nm)
            for j = 1:i-1
                ydiff += int.tableau.qdiff0.a[i,j] * int.B1[j][k,:]
            end

            int.Q0[k] = int.q[r,m][k] + int.Δt * ydrift + dot(ydiff,int.ΔW)

        end

        # Calculating the internal stages H^(l)_i for l=1..sol.nm
        @inbounds for k in eachindex(int.Q1[1])
            # contribution from the drift part (same for all noises)
            ydrift = 0.
            for j = 1:i-1
                ydrift += int.tableau.qdrift1.a[i,j] * int.V[j][k]
            end

            # contribution from the diffusion part is different for each noise
            @inbounds for noise_idx in eachindex(int.Q1)

                # ΔW contribution from the diffusion part
                ydiff .= zeros(DT, sol.nm)
                for j = 1:i-1
                    for l = 1:sol.nm
                        if l==noise_idx
                            ydiff[l] += int.tableau.qdiff1.a[i,j] * int.B1[j][k,l]
                        else
                            ydiff[l] += int.tableau.qdiff3.a[i,j] * int.B1[j][k,l]
                        end
                    end
                end

                int.Q1[noise_idx][k] = int.q[r,m][k] + int.Δt * ydrift + dot(ydiff,int.ΔW)
            end
        end

        # Calculating the internal stages \hat H^(l)_i for l=1..sol.nm
        @inbounds for k in eachindex(int.Q2[1])
            # contribution from the drift part (same for all noises)
            ydrift = 0.
            for j = 1:i-1
                ydrift += int.tableau.qdrift2.a[i,j] * int.V[j][k]
            end

            # contribution from the diffusion part is different for each noise
            @inbounds for noise_idx in eachindex(int.Q2)

                # ΔW contribution from the diffusion part
                ydiff2 = 0.0

                for j = 1:i-1
                    for l = 1:sol.nm
                        # calculating the terms with the random variables \hat I_{noise_idx,l}
                        if l<noise_idx
                            ydiff2 += int.tableau.qdiff2.a[i,j] * int.B1[j][k,l] * int.ΔW[noise_idx] * int.ΔZ[l]
                        elseif l>noise_idx
                            ydiff2 += -int.tableau.qdiff2.a[i,j] * int.B1[j][k,l] * int.ΔW[l] * int.ΔZ[noise_idx]
                        end
                    end
                end

                int.Q2[noise_idx][k] = int.q[r,m][k] + int.Δt * ydrift + ydiff2/sqrt(int.Δt)
            end
        end

        # CALCULATING THE NEW VALUES OF V
        tᵢ = sol.t[0] + (n-1)*int.Δt + int.Δt * int.tableau.qdrift0.c[i]
        int.equation.v(tᵢ, int.Q0, int.V[i])

        # CALCULATING THE NEW VALUES OF B1
        # each column of B evaluated at a different internal stage
        tᵢ = sol.t[0] + (n-1)*int.Δt + int.Δt * int.tableau.qdrift1.c[i]

        for l = 1:sol.nm
            int.equation.B(tᵢ, int.Q1[l], int.tB, col=l)
            simd_copy_yx_first!(int.tB, int.B1[i], l)
        end

        #CALCULATING THE NEW VALUES OF B2
        # each column of B evaluated at a different internal stage
        tᵢ = sol.t[0] + (n-1)*int.Δt + int.Δt * int.tableau.qdrift2.c[i]

        for l = 1:sol.nm
            int.equation.B(tᵢ, int.Q2[l], int.tB, col=l)
            simd_copy_yx_first!(int.tB, int.B2[i], l)
        end

    end

    # compute final update
    update_solution!(int.q[r,m], int.V, int.B1, int.B2, int.tableau.qdrift0.b, int.tableau.qdiff0.b, int.tableau.qdiff3.b, int.Δt, int.ΔW)

    # take care of periodic solutions
    cut_periodic_solution!(int.q[r,m], int.equation.periodicity)

    # copy to solution
    set_solution!(sol, int.q[r,m], n, r, m)
end


text:
using Test
using SparseArrays
using LinearAlgebra
import sbp
using sbp.Sparse: diag_block_matrix_2x2
using sbp.Grid: grid_xp, grid_xm, grid_2d_x, grid_2d_y

st = sbp.StaggeredAcoustic
sop = sbp.OP2019

nx=10
ny=20

function mapping(nx, ny)
        x1, h = grid_xp(nx)
        y1, h = grid_xp(ny)

        p = grid_2d_x(x1, ny)
        q = grid_2d_y(y1, nx)
        x = p
        y = (1 .+ 0.1 * sin.(p)).*q
        return x, y
end

fx, fy = mapping(nx, ny)
desc = sop.load_description()
builder = n -> sop.build_operators(desc, n)
ops = st.init_operators(nx, ny, builder, fx, fy)

HJp = st.pressure_norm(ops)
HJ = st.velocity_norm(ops)
Gp = st.contravariant_metric_tensor(ops)
Ek = HJ*Gp
Ap = st.divergence(ops)
Av = st.gradient(ops, Gp)
S = st.pressure_sat(ops, Gp)
H = st.energy_norm(HJp, HJ, Gp)
A = st.spatial_discretization(Ap, Av - S)
d = st.pressure_source(ops, 0.5, 0.5, 4, 4)

@testset "Grids" begin
        xp, yp = st.grids("p", nx, ny)
        @test size(xp,1) == (nx + 1) * (ny + 1)
        @test size(yp,1) == (nx + 1) * (ny + 1)
        x1, y1 = st.grids("v1", nx, ny)
        @test size(x1,1) == (nx + 0) * (ny + 1)
        @test size(y1,1) == (nx + 0) * (ny + 1)
        x2, y2 = st.grids("v2", nx, ny)
        @test size(x2,1) == (nx + 1) * (ny + 0)
        @test size(y2,1) == (nx + 1) * (ny + 0)
        xn, yn = st.grids("node", nx, ny)
        @test size(xn,1) == (nx + 0) * (ny + 0)
        @test size(yn,1) == (nx + 0) * (ny + 0)
end

@testset "Kinetic energy tensor" begin
        # Check symmetry
        @test isapprox(Ek, (Ek)')
        # Check positivity
        lam = real(eigvals(Matrix(Ek)))
        @test minimum(lam) > 0

end

@testset "Stability" begin
        # Check that the discretization is energy conservative to rounding error
        lam = real(eigvals(Matrix(A)))
        @test minimum(lam) > -1e-13
        @test maximum(lam) < 1e-13

        E = H*A + (H*A)'
        lam = real(eigvals(Matrix(E)))
        @test minimum(lam) > -1e-13
        @test maximum(lam) < 1e-13

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
## load in the data
using Test
using PowerSpectra
using Healpix
using JLD2
using DelimitedFiles

@testset "Planck 100 GHz MCM" begin
    nside = 256
    lmax = nside2lmax(nside)

    maskT₁ = PowerSpectra.planck256_mask("100", "hm1", "T")
    maskP₁ = PowerSpectra.planck256_mask("100", "hm1", "P")
    maskT₂ = PowerSpectra.planck256_mask("100", "hm2", "T")
    maskP₂ = PowerSpectra.planck256_mask("100", "hm2", "P")
    flat_map = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )

    M₀₀ = mcm(:TT, map2alm(maskT₁), map2alm(maskT₂))
    @load "data/mcm00.jld2" mcm00
    @test mcm00 ≈ parent(M₀₀)
    ##
    @load "data/mcm20.jld2" mcm20
    A = mcm20[1:2:2*lmax+1, 1:2:2*lmax+1]
    M₂₀ = mcm(:ET, map2alm(maskP₁), map2alm(maskT₂))
    @test A[3:end, 3:end] ≈ parent(M₂₀)[3:end, 3:end]
    M₂₀ = mcm(:BT, map2alm(maskP₁), map2alm(maskT₂))
    @test A[3:end, 3:end] ≈ parent(M₂₀)[3:end, 3:end]
    @load "data/mcm02.jld2" mcm02
    A = mcm02[1:2:2*lmax+1, 1:2:2*lmax+1]
    M₀₂ = mcm(:TE, map2alm(maskT₁), map2alm(maskP₂))
    @test A[3:end, 3:end] ≈ parent(M₀₂)[3:end, 3:end]
    M₀₂ = mcm(:TB, map2alm(maskT₁), map2alm(maskP₂))
    @test A[3:end, 3:end] ≈ parent(M₀₂)[3:end, 3:end]

    ##
    @load "data/mcm22.jld2" mcm22
    M⁺⁺ = mcm(:M⁺⁺, map2alm(maskP₁), map2alm(maskP₂))
    M⁻⁻ = mcm(:M⁻⁻, map2alm(maskP₁), map2alm(maskP₂))

    ##
    A = mcm22[1:4:4*(lmax+1), 1:4:4*(lmax+1)]
    @test A[3:end, 3:end] ≈ parent(M⁺⁺)[3:(lmax+1), 3:(lmax+1)]

    ##
    Δ₁, Δ₂ = 1, 1
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(M⁺⁺)[3:(lmax+1), 3:(lmax+1)]

    Δ₁, Δ₂ = 2, 2
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(M⁺⁺)[3:(lmax+1), 3:(lmax+1)]

    Δ₁, Δ₂ = 3, 3
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(M⁺⁺)[3:(lmax+1), 3:(lmax+1)]

    Δ₁, Δ₂ = 1, 2
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(-M⁻⁻)[3:(lmax+1), 3:(lmax+1)]
    ##
    Δ₁, Δ₂ = 2, 1
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(-M⁻⁻)[3:(lmax+1), 3:(lmax+1)]

    Δ₁, Δ₂ = 0, 3
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(M⁻⁻)[3:(lmax+1), 3:(lmax+1)]

    Δ₁, Δ₂ = 3, 0
    A = mcm22[1+Δ₁:4:4*(lmax+1)+Δ₁, 1+Δ₂:4:4*(lmax+1)+Δ₂]
    @test A[3:end, 3:end] ≈ parent(M⁻⁻)[3:(lmax+1), 3:(lmax+1)]
end

##
@testset "Planck 100 GHz MCM decoupling" begin
    nside = 256
    m₁ = PowerSpectra.planck256_polmap("100", "hm1")
    m₂ = PowerSpectra.planck256_polmap("100", "hm2")
    maskT₁ = PowerSpectra.planck256_mask("100", "hm1", :T)
    maskP₁ = PowerSpectra.planck256_mask("100", "hm1", :P)
    maskT₂ = PowerSpectra.planck256_mask("100", "hm2", :T)
    maskP₂ = PowerSpectra.planck256_mask("100", "hm2", :P)

    # convert to μK
    scale!(m₁, 1e6)
    scale!(m₂, 1e6)
    Cl = master(m₁, maskT₁, maskP₁,
                m₂, maskT₂, maskP₂; lmin=2)

    @load "data/planck_spec.jld2" cl_00 cl_02 cl_20 cl_22

    @test cl_00[1,:] ≈ Cl[:TT][2:end]

    @test cl_02[1,:] ≈ Cl[:TE][2:end]
    @test cl_02[2,:] ≈ Cl[:TB][2:end]
    @test cl_20[1,:] ≈ Cl[:ET][2:end]
    @test cl_20[2,:] ≈ Cl[:BT][2:end]

    @test cl_22[1,:] ≈ Cl[:EE][2:end]
    @test cl_22[2,:] ≈ Cl[:EB][2:end]
    @test cl_22[3,:] ≈ Cl[:BE][2:end]
    @test cl_22[4,:] ≈ Cl[:BB][2:end]
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using SparseArrays

#--------------------------------------
# Testing transformations between bases
#--------------------------------------

@testset "transformations.jl" begin
    n = 5;
    import GalerkinSparseGrids.transform_1D

    @info "Testing 1D transformation between point and nodal bases... "
    @testset "Point to Nodal 1D" begin
        for k in [2,3,5]
            p2n = transform_1D(k, n, "points", "nodal")
            n2p = transform_1D(k, n, "nodal", "points")
            @test norm(p2n*n2p-I) < 1e-15*(10^k)
            @test norm(n2p*p2n-I) < 1e-15*(10^k)
        end
    end

    @info "Testing 1D transformation between point and modal bases... "
    @testset "Point to Modal 1D" begin
        for k in [2,3,5]
            p2m = transform_1D(k, n, "points", "modal")
            m2p = transform_1D(k, n, "modal", "points")
            @test norm(p2m*m2p-I) < 1e-15*(10^k)
            @test norm(m2p*p2m-I) < 1e-15*(10^k)
        end
    end

    @info "Testing 1D transformation between nodal and modal bases... "
    @testset "Nodal to Modal 1D" begin
        for k in [2,3,5]
            n2m = transform_1D(k, n, "nodal", "modal")
            m2n = transform_1D(k, n, "modal", "nodal")
            @test norm(n2m*m2n-I) < 1e-15*(10^k)
            @test norm(m2n*n2m-I) < 1e-15*(10^k)
        end
    end

    @info "Testing 1D transformation from position to modal bases... "
    @testset "Position to Modal 1D" begin
        for k in [2,3,5]
            pos2m = transform_1D(k, n, "pos", "modal")
            m2pos = transform_1D(k, n, "modal", "pos")
            @test norm(pos2m*m2pos-I) < 1e-15*(10^k)
            @test norm(m2pos*pos2m-I) < 1e-15*(10^k)
        end
    end

    D = 2; k = 3; n = 5;
    @info "Testing 2D transformation between nodal and points bases... "
    @testset "Points to Nodal 2D" begin
        p2n = transform(D, k, n, "points", "nodal")
        n2p = transform(D, k, n, "nodal", "points")
        @test norm(p2n*n2p - I) < 1e-10
        @test norm(n2p*p2n - I) < 1e-10
    end

    @info "Testing 2D transformation between modal and points bases... "
    @testset "Points to Modal 2D" begin
        p2m = transform(D, k, n, "points", "modal")
        m2p = transform(D, k, n, "modal", "points")
        @test norm(p2m*m2p - I) < 1e-10
        @test norm(m2p*p2m - I) < 1e-10
    end

    @info "Testing 2D transformation between modal and nodal bases... "
    @testset "Nodal to Modal 2D" begin
        n2m = transform(D, k, n, "nodal", "modal")
        m2n = transform(D, k, n, "modal", "nodal")
        @test norm(m2n*n2m - I) < 1e-10
        @test norm(n2m*m2n - I) < 1e-10
    end
end


text:
struct topic_modeling_hyper <: distribution_hyper_params
    α::AbstractArray{Float64,1}
end

mutable struct topic_modeling_sufficient_statistics <: sufficient_statistics
    N::Float64
    points_sum::AbstractArray{Int64,1}
end

function calc_posterior(prior:: topic_modeling_hyper, suff_statistics::topic_modeling_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    return topic_modeling_hyper(prior.α + suff_statistics.points_sum)
end

function sample_distribution(hyperparams::topic_modeling_hyper)
    cat_dist = rand(Dirichlet(hyperparams.α))
    return topic_modeling_dist(log.(cat_dist))
end

function create_sufficient_statistics(hyper::topic_modeling_hyper,posterior::topic_modeling_hyper,points::AbstractArray{Float64,2}, pts_to_group = 0)
    if length(points) == 0
        return topic_modeling_sufficient_statistics(size(points,2),zeros(Int64,length(hyper.α)))
    end
    pt_count = counts(Int.(points),length(hyper.α))
    return topic_modeling_sufficient_statistics(size(points,2),pt_count)
end

function log_marginal_likelihood(hyper::topic_modeling_hyper, posterior_hyper::topic_modeling_hyper, suff_stats::topic_modeling_sufficient_statistics)
    D = length(suff_stats.points_sum)
    N = Int(ceil(sum(posterior_hyper.α - hyper.α)))
    logpi = log(pi)
    val = sum((x-> logabsgamma(x)[1]).(posterior_hyper.α) - suff_stats.points_sum .* (x-> logabsgamma(x)[1]).(hyper.α)) + logabsgamma(sum(hyper.α))[1] -logabsgamma(N + sum(hyper.α))[1]
    return val
end

function aggregate_suff_stats(suff_l::topic_modeling_sufficient_statistics, suff_r::topic_modeling_sufficient_statistics)
    return topic_modeling_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum)
end


text:
# ---------------------------------------------------------------------------- #
#
#   divido.jl
#
#   Mesh data structures for FEM and lattice structures
#
#   divido
#   Fall 2017 / Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

__precompile__()

"""
    divido

Julia package for mesh data structures for FEM and lattice structures.

Max Opgenoord

Fall 2017 / Spring 2018
"""

module divido

using GroupSlices

# mesh
export Mesh2D, Mesh3D
include("mesh/mesh.jl")

# links to meshers
export writeMsh, readBAMG, runMesher, MesherFEFLOA, MesherBAMG, writeMsh
export readFEFLOA_2D, readFEFLOA_3D, writeMetric
include("meshers/mesher.jl")

# io
include("io/readSU2.jl")

end


text:
# ---
# title: 104. Maximum Depth of Binary Tree
# id: problem104
# author: zhwang
# date: 2022-02-20
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/maximum-depth-of-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _its maximum depth_.
#
# A binary tree's **maximum depth**  is the number of nodes along the longest
# path from the root node down to the farthest leaf node.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/tmp-tree.jpg)
#
#
#
#     Input: root = [3,9,20,null,null,15,7]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: root = [1,null,2]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: root = []
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: root = [0]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-100 <= Node.val <= 100`
#
#
## @lc code=start
using LeetCode

maximum_depth_of_btree(::Nothing) = 0
function maximum_depth_of_btree(root::TreeNode)
    return maximum(maximum_depth_of_btree, [root.left, root.right]) + 1
end
## @lc code=end


text:
# # The shared write problem on GPU

# We will write a GPU version of `axpy!` function.

# ## The main program

using NiLang, NiLang.AD
using CUDA
using KernelAbstractions
CUDA.allowscalar(true)

# so far, this example requires patch: https://github.com/JuliaGPU/KernelAbstractions.jl/pull/52

@i @kernel function axpy_kernel(y!, α, x)
    ## invcheckoff to turn of `reversibility checker`
    ## GPU can not handle errors!
    @invcheckoff begin
        i ← @index(Global)
        y![i] += x[i] * α
        i → @index(Global)
    end
end

@i function cu_axpy!(y!::AbstractVector, α, x::AbstractVector)
    @launchkernel CUDADevice() 256 length(y!) axpy_kernel(y!, α, x)
end

@i function loss(out, y!, α, x)
    cu_axpy!(y!, α, x)
    ## Note: the following code is stupid scalar operations on CuArray,
    ## They are only for testing.
    for i=1:length(y!)
        out += y![i]
    end
end

y! = rand(100)
x = rand(100)
cuy! = y! |> CuArray
cux = x |> CuArray
α = 0.4

# ## Check the correctness of results

using Test
cu_axpy!(cuy!, α, cux)
@test Array(cuy!) ≈ y! .+ α .* x
(~cu_axpy!)(cuy!, α, cux)
@test Array(cuy!) ≈ y!

# Let's check the gradients
lsout = 0.0
@instr Grad(loss)(Val(1), lsout, cuy!, α, cux)

# you will see a correct vector `[0.4, 0.4, 0.4 ...]`
grad.(cux)

# you will see `0.0`.
grad(α)

# ## Why some gradients not correct?
# In the above example, `α` is a scalar, whereas a scalar is not allowed to change in a CUDA kernel.
# What if we change `α` to a CuArray?

# ## This one works: using a vector of `α`
@i @kernel function axpy_kernel(y!, α, x)
    @invcheckoff begin
        i ← @index(Global)
        y![i] += x[i] * α[i]
        i → @index(Global)
    end
end

cuy! = y! |> CuArray
cux = x |> CuArray
cuβ = repeat([0.4], 100) |> CuArray
lsout = 0.0
@instr Grad(loss)(Val(1), lsout, cuy!, cuβ, cux)

# You will see correct answer
grad.(cuβ)

# ## This one has the shared write problem: using a vector of `α`, but shared read.
@i @kernel function axpy_kernel(y!, α, x)
    @invcheckoff begin
        i ← @index(Global)
        y![i] += x[i] * α[i]
        i → @index(Global)
    end
end

cuy! = y! |> CuArray
cux = x |> CuArray
cuβ = repeat([0.4], 100) |> CuArray
lsout = 0.0
cuβ = [0.4] |> CuArray

# Run the following will give you a happy error
#
# > ERROR: a exception was thrown during kernel execution.
# >        Run Julia on debug level 2 for device stack traces.

# ```julia
# @instr Grad(loss)(Val(1), lsout, cuy!, cuβ, cux)
# ```

# Because, shared write is not allowed. We need someone clever enough to solve this problem for us.

# ## Conclusion
# * Shared scalar: the gradient of a scalar will not be updated.
# * Expanded vector: works properly, but costs more memory.
# * Shared 1-element vector: error on shared write.


text:
basedir = realpath(joinpath(@__DIR__, ".."))
floattype = Float64
# Set the input data used in human dispersal models
humanactivity = :gdp
# humanactivity = :pop

# Set simulation replicates and grid aggregation
# Fast but innacurate
nreps = 2
scale = 10
# Slow but accurate: for the paper
# nreps = 100
# scale = 2

include(joinpath(basedir, "src", "spreadsetup.jl"))

# locname =                 :US
# locname =                 :EU
locname =                 :Aus

### Load everything required ###

if locname == :US
    tspan =                   us_tspan
    ruleset =                 us_ruleset
    boolmask =                us_boolmask
    missingmask_ =            us_missingmask
    growthrates_zeromissing = us_growthrates_zeromissing
    populationgrid =          us_populationgrid
    incursionpoints =         us_incursionpoints
    incursion =               us_incursionpoints[:San_Francisco]
elseif locname == :EU
    tspan =                   eu_tspan
    ruleset =                 eu_ruleset
    boolmask =                eu_boolmask
    missingmask_ =            eu_missingmask
    growthrates_zeromissing = eu_growthrates_zeromissing
    populationgrid =          eu_populationgrid
    incursionpoints =         eu_incursionpoints
    #incursion = eu_incursionpoints[:Rome]
elseif locname == :Aus
    tspan =                   DateTime(2021, 1):timestep:DateTime(2026, 1)
    ruleset =                 aus_ruleset
    boolmask =                aus_boolmask
    missingmask_ =            aus_missingmask
    growthrates_zeromissing = aus_growthrates_zeromissing
    populationgrid =          aus_populationgrid
    incursionpoints =         aus_incursionpoints
    incursion = aus_incursionpoints[:Brisbane]
end
aux = (growthrates=growthrates_zeromissing,)
init = (population=populationgrid,)
# Reconstruct with fitted parameters if you have them allready
using JLD2, FileIO, LabelledArrays
stochparams = load(joinpath(basedir, "output/stochparams_$humanactivity.jld2"), "stochparams_$humanactivity")
ruleset = reconstruct(ruleset, stochparams[:full])

# Run the simulation in a live interface :
# unccomment to run in a live interface
# using DynamicGridsInteract
# output = ElectronOutput((population=parent(populationgrid),);
#     ruleset=ruleset,
#     mask=BitArray(collect(boolmask)),
#     aux=aux,
#     tspan=tspan,
#     store=true,
#     processor=inferno,
#     minval=zero(carrycap),
#     maxval=carrycap,
#     fps=50,
# )
# display(output)
# And save the simulation as a gif:
# savegif(joinpath(basedir, "output/dispersal.gif"), output, ruleset; fps=10)

## Mapping incursion-point sensitivity

function incursionreps!(output, ruleset, missingmask, nreps)
    steps_established = similar(output[1][:population], Int)
    steps_established .= 0
    @set! steps_established.refdims = ()
    for rep in 1:nreps
        sim!(output, ruleset)
        println("rep: ", rep)
        map(output) do step
            steps_established .+= step[:population] .>= 1
        end
    end
    GeoArray(steps_established ./ nreps .* missingmask; missingval=missing, name=:established)
end
function incursionplot(established, key, nreps)
    plot(established;
        color=:inferno,
        xlabel="",
        ylabel="",
        title="$key: num timesteps established with $nreps reps"
    )
end
function set_incursion(output, (lat, lon))
    DynamicGrids.init(output)[:population] .= 0
    DynamicGrids.init(output)[:population][Lat(Contains(lat)), Lon(Contains(lon))] = carrycap
end

# Plot the initial incursion:

output = ArrayOutput(init; tspan=tspan, aux=aux, mask=boolmask)
established = incursionreps!(output, ruleset, missingmask_, nreps)
incursionplot(established, locname, nreps)
filename = joinpath(basedir, "output/$(locname)_steps_established_$(nreps)reps")
savefig("$filename.png")
write("$filename.tif", GDALarray, replace_missing(established, typemin(Float64)))
# Check the save worked
GDALarray("$filename.tif"; usercrs=EPSG(4326)) |> plot

# Plot all the incursions in incursionpoints:
plots = []

for (key, loc) in zip(keys(incursionpoints), incursionpoints)
    set_incursion(output, loc)
    steps_established = incursionreps!(output, ruleset, missingmask_, nreps)
    filename = joinpath(basedir, "output/$(locname)_steps_established_$(nreps)reps")
    write("$filename.tif", GDALarray, replace_missing(steps_established, typemin(Float64)))
    push!(plots, incursionplot(steps_established, key, nreps))
end

# Plot all togeth
plot(plots...; size=(1000, 600))
savefig(joinpath(basedir, "output/months_established.png"))

#' We can also save the plots as individual figures.
for (i, key) in enumerate(keys(incursionpoints))
    plot(plots[i])
    savefig(joinpath(basedir, "output/months_established_$key.png"))
end

### Incursion point grid ###

function incursion_point_grid!(
    cellsinvaded6, cellsinvaded12, ruleset, extent, nreps, scale, season;
    writedir="$basedir/output/tiffs", write=false
)
    if write
        rm(writedir; recursive=true)
        mkdir(writedir)
    end
    init = DynamicGrids.init(extent)
    londim, latdim = dims(init[:population], (Lon, Lat))
    latindex = reverse(val(latdim))
    outputs = [ArrayOutput(init; extent=extent) for t in 1:Threads.nthreads()]
    simdata = [DynamicGrids.SimData(extent, ruleset) for t in 1:Threads.nthreads()]
    inits = [deepcopy(init) for t in 1:Threads.nthreads()]
    Threads.@threads for i = 1:scale:size(cellsinvaded12, 1) * scale
        acc = similar(init[:population], UInt8)
        acc .= 0
        @set! acc.missingval = 0x00
        threadinit = inits[Threads.threadid()]
        threaddata = simdata[Threads.threadid()]
        threadoutput = outputs[Threads.threadid()]
        for j = 1:scale:size(cellsinvaded12, 2) * scale
            println("i, j: ", (i, j))
            DynamicGrids.mask(extent)[i, j] || continue
            checkbounds(Bool, acc, i, j) || continue
            lat, lon = ArchGDAL.reproject([[londim[j], latindex[i]]], crs(init[:population]), EPSG(4326))[1]
            # Zero the accumulator array
            acc .= 0
            # Zero init array
            inits[Threads.threadid()][:population] .= 0
            # Set the current cell to maximum population
            inits[Threads.threadid()][:population][i, j] = carrycap
            invaded6 = 0
            invaded12 = 0
            for k in 1:nreps
                sim!(threadoutput, ruleset; init=inits[Threads.threadid()], simdata=threaddata)
                invaded6 += count(x -> x > one(x), threadoutput[7][:population])
                invaded12 += count(x -> x > one(x), threadoutput[13][:population])
                acc .+= threadoutput[13][:population] .> 0
            end
            cellsinvaded6[(i - 1) ÷ scale + 1, (j - 1) ÷ scale + 1] = invaded6 / nreps
            cellsinvaded12[(i - 1) ÷ scale + 1, (j - 1) ÷ scale + 1] = invaded12 / nreps
            if write && invaded12 > 0
                filepath = "$writedir/cells_invaded_from_$(lon)_$(lat)_$season.tif"
                # Remove underscores
                filepath = replace(filepath, r"(.*)\.(.*)\.(.*)\.tif" => s"\1_\2_\3.tif")
                GeoData.write(filepath, GDALarray, acc)
            end
        end
    end
end

function spread_probability!(output, I...;
    gridname,
    initval,
    ruleset=ruleset(output),
    acc=zero(output.extent.init[gridname]),
    init=output.extent.init,
    nreps=1,
    simdata=nothing
)
    invaded = 0
    acc .= 0
    initgrids = deepcopy(init)
    initgrids[gridname] .= 0
    # Set the current cell to maximum population
    initgrids[gridname][I...] = carrycap
    @show nreps
    for k in 1:nreps
        sim!(output, ruleset; init=initgrids, simdata=simdata)
        invaded += count(x -> x > one(x), output[end][gridname])
        acc .+= output[end][gridname] .> 0
    end
    acc ./= nreps
    return acc
end

pyplot()
oneyear = DateTime(2021):Month(1):DateTime(2022)
output = ArrayOutput(init; tspan=oneyear, aux=aux, mask=boolmask)
spreaddir = joinpath(basedir, "output/spread_probability")
mkpath(spreaddir)

for key in keys(incursionpoints)
    println("==== Generating spread probability for $key ======================")
    lat, lon = incursionpoints[key]
    spreadprob = spread_probability!(output, Lon(Contains(lon)), Lat(Contains(lat));
        initval=carrycap,
        gridname=:population,
        ruleset=ruleset,
        simdata=nothing,
        nreps=nreps,
    )
    filename = joinpath(spreaddir, "$(key)_spread_probability_$(nreps)reps")
    plot(spreadprob; title=key, colorbar_title="Spread probability")
    savefig(string(filename, ".png"))
    write(string(filename, ".tif"), GDALarray, spreadprob)
end

# Incursion locations from a spreadsheet

cellsinvaded = GeoData.aggregate(Center(), populationgrid, scale) .* 0
@set! cellsinvaded.name = "Cells invaded"
cellsinvaded6_summer = deepcopy(cellsinvaded)
cellsinvaded12_summer = deepcopy(cellsinvaded)
cellsinvaded6_winter = deepcopy(cellsinvaded)
cellsinvaded12_winter = deepcopy(cellsinvaded)
tspan_summer = DateTime(2021, 1):timestep:DateTime(2022, 1)
tspan_winter = DateTime(2021, 7):timestep:DateTime(2022, 7)
init = (population=zero(populationgrid),)

# Summer sims
extent = DynamicGrids.Extent(init, boolmask, aux, tspan)
@time incursion_point_grid!(
    cellsinvaded6_summer, cellsinvaded12_summer, ruleset, extent, nreps, scale, "summer";
    write=false
)

# Winter sims
@set! extent.tspan = tspan_winter
@time incursion_point_grid!(
    cellsinvaded6_winter, cellsinvaded12_winter, ruleset, extent, nreps, scale, "winter";
    write=false
)

# Now plot and save the combined cells invaded arrays:

function plot_cells_invaded(ci, label, months, nreps, mmask; kwargs...)
    toarea = 9*9*1e-6
    area = ci .* toarea .* GeoData.aggregate(Center(), mmask, scale)
    p = plot(area;
        colorbar_title="Area invaded from incursion pt. (million km²)",
        title="Area invaded after $months months",
        kwargs...
    )
    path = "output/cellsinvaded_$(label)_$(months)months_$(nreps)reps"
    savefig(path * ".png")
    write(path * ".tif", GDALarray, replace_missing(ci, -1.00e+10))
    p
end
plot_cells_invaded(cellsinvaded6_summer, "summer_$humanactivity", 6, nreps, missingmask_)
plot_cells_invaded(cellsinvaded12_summer, "summer_$humanactivity", 12, nreps, missingmask_)
plot_cells_invaded(cellsinvaded6_winter, "winter_$humanactivity", 6, nreps, missingmask_)
plot_cells_invaded(cellsinvaded12_winter, "winter_$humanactivity", 12, nreps, missingmask_)

max_stars_repo_name: rafaqz/DSuzukiiInvasionPaper

max_issues_repo_name: rafaqz/DSuzukiiInvasionPaper


text:
path = "C:\\Users\\vsger\\Downloads\\rosalind_lexf (1).txt"
a = readlines(path)
alphabet = a[1]
l = parse(Int64,a[2])

a = sort(join.(vcat(Iterators.product(ntuple(_ -> split(alphabet," "),l)...)...)))
for i in a
	println(i)
end

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
function rhosp_mixed(probs::Array{Float64,1}, states::Union{Array{State,1},Array{State_fixed,1}})
    if round(sum(probs), digits=15) != 1
        throw(ArgumentError("Your probabilities must add up to 1!"))
    end

    ls = length(states)
    if length(probs) != ls
        throw(ArgumentError("There must be the same ammount of probabilities and states!"))
    end

    #for i in 1:length(states)
    #    stat = st(states[i])
    #    if round(stat'*stat, digits=15) != 1
    #        throw(ArgumentError("All states must be normalized to 1!"))
    #    end
    #end

    rhos = [probs[i]*rhosp(states[i]) for i in 1:ls]
    rho = sum(rhos)
end

function eigensp_mixed(probs::Array{Float64,1}, states::Union{Array{State,1},Array{State_fixed,1}})
    if round(sum(probs), digits=15) != 1
        throw(ArgumentError("Your probabilities must add up to 1!"))
    end

    ls = length(states)
    if length(probs) != ls
        throw(ArgumentError("There must be the same ammount of probabilities and states!"))
    end

    #for i in 1:length(states)
    #    stat = st(states[i])
    #    if round(stat'*stat, digits=15) != 1
    #        throw(ArgumentError("All states must be normalized to 1!"))
    #    end
    #end

    eig = sort(eigvals(Matrix(rhosp_mixed(probs, states))), rev=true)
    return eig
end


text:
# TODO Change Vector{Num} to either Symbolics.Arr{Num, 1} or something more general. May require some finessing

export  MPCCDimSpec,
        MPCCFunctions,
        MPCCPointEvalReq,
        MPCCPointEval,
        MPCCModelTestVector,
        MPCCModelNZMask,
        MPCCModelConfig,
        MPCCModel,
        MPCCJumpFixedFunctions,
        MPCCNewtonPenaltyFunctions

# Trick from https://discourse.julialang.org/t/aliases-for-union-t-nothing-and-union-t-missing/15402/15
const Opt{T} = Union{Missing,T}

"""
MPCCDimSpec

Determines the dimensions of a program.

* n: number of spatial dimensions
* q: columns in F (number of complementarity contraints)
* l: row in F (number of expressions in each complementarity contraint)
* me: numebr of equality constraints
* me: numebr of inequality constraints
* r: number of continuous parameters
* s: number of integer parameters

"""
struct MPCCDimSpec
    n::Int64		# Dimension of spatial variables
    q::Int64		# Columns in F
    l::Int64        # Number of F functions per complementarity variable (usually 2)
    me::Int64		# Number of equality constraints
    mi::Int64		# Number of inequality constraints
	r::Int64		# Number of continuous parameters for PMPCC
	s::Int64		# Number of discrete parameters for PMPCC
end
function MPCCDimSpec(dict::Dict)
    return MPCCDimSpec(
        dict["n"],
        dict["q"],
        dict["l"],
        dict["me"],
        dict["mi"],
        dict["r"],
        dict["s"]
    )
end

"""
struct MPCCDefinition

Definition of the MPCC/NLP containing Symbolics Num expressions for f, ce, ci,
and F.

* f: objective function
* ce: vector of equality constraints ( == 0 )
* ci: vector of inequality constraints ( >= 0 )
* F: matrix l-by-q of complementarity constraints ( >= 0 )

"""
struct MPCCDefinition
    f::Num
    ce::Vector{Num}
    ci::Vector{Num}
    F::Matrix{Num}
end

"""
MPCCFunctions

A MPCCDefinition compiled into functions. This contains both standard and
mutating functions, and individual elements encoded as a function within a
vector or matrix.
"""
struct MPCCFunctions
    f::Function
    f!::Function
    ce::Function                   # Returns vector of length me
    ce!::Function
    ce_i::AbstractVector           # Vector so we can write scalar values in to mutating function argument
    ce_i!::AbstractVector
    ci::Function                   # Returns vector of length mi
    ci!::Function
    ci_i::AbstractVector
    ci_i!::AbstractVector
    F::Function                    # Returns vatrix of dim l x q
    F!::Function
    F_i::AbstractMatrix
    F_i!::AbstractMatrix
    Fq::Function                    # Returns vector of length q
    Fq!::Function
    Fq_i::AbstractVector
    Fq_i!::AbstractVector
end

"""
MPCCPointEvalReq

When making a request to evaluate several things in one request, this bitmask
determines what will be evaluated.
"""
struct MPCCPointEvalReq
    f::Bool
    ce::Bool
    ci::Bool
    F::Bool
    gradf::Bool
    gradce::Bool
    gradci::Bool
    gradF::Bool
    hessf::Bool
    hessce::Bool
    hessci::Bool
    hessF::Bool
    fdp::Bool
    cedp::Bool
    cidp::Bool
    Fdp::Bool
    gradfdp::Bool
    gradcedp::Bool
    gradcidp::Bool
    gradFdp::Bool
end
function MPCCPointEvalReq()
    return MPCCPointEvalReq(    true, true, true, true,
                                true, true, true, true,
                                true, true, true, true,
                                true, true, true, true,
                                true, true, true, true
                            )
end

"""
MPCCPointEval

Contains results of a point evaluation.
"""
struct MPCCPointEval{T <: AbstractFloat}
    f::Opt{T}
    ce::Opt{Vector{T}}                          # Vector of length me
    ci::Opt{Vector{T}}                          # Vector of length mi
    F::Opt{Matrix{T}}                           # Matrix of dim l x q
    gradf::Opt{Vector{T}}                       # Vector of length n
    gradce::Opt{Matrix{T}}                      # Changed 20211029: now matrix of me by n
    gradci::Opt{Matrix{T}}                      # Changed 20211029: now matrix of mi by n
    gradF::Opt{Matrix{Vector{T}}}               # Matrix of dim l x q of n length vectors (the gradients wrt x)
    hessf::Opt{Matrix{T}}                       # Matrix dim n x n
    hessce::Opt{Vector{Matrix{T}}}              # Vector of length me of 2d Hessian matrices
    hessci::Opt{Vector{Matrix{T}}}              # Vector of length mi of 2d Hessian matrices
    hessF::Opt{Matrix{Matrix{T}}}               # Matrix of dim l x q of 2d Hessian matrices
    fdp::Opt{Vector{T}}                         # Vector of length r
    cedp::Opt{Vector{Vector{T}}}                # Vector of length r of vectors of length me
    cidp::Opt{Vector{Vector{T}}}                # Vector of length r of vectors of length mi
    Fdp::Opt{Vector{Matrix{T}}}                 # Vector of length r of matrix of dim l x q
    gradfdp::Opt{Vector{Vector{T}}}             # Vector of length r of gradients of f/dp (vectors of length n)
    gradcedp::Opt{Vector{Matrix{T}}}            # Vector of length r of gradce, each /dp
    gradcidp::Opt{Vector{Matrix{T}}}            # Vector of length r of gradci, each /dp
    gradFdp::Opt{Vector{Matrix{Vector{T}}}}     # Vector of length r of gradF, each /dp
end
function MPCCPointEval()  # Quick way to create empty/zero struct
	return MPCCPointEval(
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing,
        missing
    )
end
Base.:(==)(x::MPCCPointEval, y::MPCCPointEval) = (
	x.f == y.f &&
    x.ce == y.ce &&
    x.ci == y.ci &&
	x.F == y.F &&
	x.gradf == y.gradf &&
    x.gradce == y.gradce &&
    x.gradci == y.gradci &&
    x.gradF == y.gradF &&
    x.hessf == y.hessf &&
    x.hessce == y.hessce &&
    x.hessci == y.hessci &&
    x.hessF == y.hessF &&
    x.fdp == y.fdp &&
    x.cedp == y.cedp &&
    x.cidp == y.cidp &&
    x.Fdp == y.Fdp &&
    x.gradfdp == y.gradfdp &&
    x.gradcedp == y.gradcedp &&
    x.gradcidp == y.gradcidp &&
    x.gradFdp == y.gradFdp )
Base.:(≈)(x::MPCCPointEval, y::MPCCPointEval) = (
    x.f ≈ y.f &&
    x.ce ≈ y.ce &&
    x.ci ≈ y.ci &&
    x.F ≈ y.F &&
    x.gradf ≈ y.gradf &&
    x.gradce ≈ y.gradce &&
    x.gradci ≈ y.gradci &&
    x.gradF ≈ y.gradF &&
    x.hessf ≈ y.hessf &&
    x.hessce ≈ y.hessce &&
    x.hessci ≈ y.hessci &&
    x.hessF ≈ y.hessF &&
    x.fdp ≈ y.fdp &&
    x.cedp ≈ y.cedp &&
    x.cidp ≈ y.cidp &&
    x.Fdp ≈ y.Fdp &&
    x.gradfdp ≈ y.gradfdp &&
    x.gradcedp ≈ y.gradcedp &&
    x.gradcidp ≈ y.gradcidp &&
    x.gradFdp ≈ y.gradFdp )

# An aid for plotting only, determines which variables are used in which constraints
struct MPCCModelNZMask
    ce::Vector{Set{Int64}}
    ci::Vector{Set{Int64}}
    F::Vector{Vector{Set{Int64}}}
end

struct MPCCModelTestVector{R <: Real, S <: Real, T <: Real}
    x_val::Vector{S}
    pr_val::Vector{T}
    ps_val::Vector{Int64}
    eval_test_pt::MPCCPointEval{R}
end

"""
MPCCParameterisationDefn

The definition for a parameterisation of a model, i.e. how the `pr` depends on `t`.
"""
struct MPCCParameterisationDefn{R <: Real}
    pr::Vector{Num}
    # prdt::Num
    tspan::Tuple{R, R}
    descr::String
end

"""
MPCCParameterisationFunctions

Compilation of `MPCCParameterisationFunctions` into Julia functions, including derivative.
"""
struct MPCCParameterisationFunctions
    pr::Function            # t -> Vector pr
    pr!::Function
    prdt::Function          # t -> Vector of d(pr) / dt
    prdt!::Function
end

"""
MPCCParameterisations

Collection of Num `t` along with definitions of parameterisations and their
compilation.
"""
struct MPCCParameterisations
    t::Num  # Symbolics variable
    defns::Vector{MPCCParameterisationDefn}
    fns::Vector{MPCCParameterisationFunctions}
end

"""
MPCCModelConfig

Model config contains the Num variables for `x`, `pr`, `ps` along with the
`MPCCDimSpec`, definition, compiled definition, any test vectors, any known
solutions, and parameterisations. Basically everything that specifies a model.
"""
struct MPCCModelConfig
    x::Vector{Num}      # Symbolics variables
    pr::Vector{Num}     # Symbolics variables
    ps::Vector{Num}     # Symbolics variables
    dimspec::MPCCDimSpec
    defn::MPCCDefinition
    fns::MPCCFunctions
    nzmask::Opt{MPCCModelNZMask}
    testvectors::Vector{MPCCModelTestVector}
    knownsols::Vector{Function}
    parameterisations::MPCCParameterisations
end

"""
MPCCModel

A `MPCCModelConfig` compiled to produce the derivatives we want.
"""
struct MPCCModel
    config::MPCCModelConfig
    f::Function
    f!::Opt{Function}
    ce::Function
    ce!::Opt{Function}
    ci::Function
    ci!::Opt{Function}
    F::Function
    F!::Opt{Function}
    Fq::Function
    Fq!::Opt{Function}
    gradf::Function
    gradf!::Opt{Function}
    gradce::Function
    gradce!::Opt{Function}
    gradci::Function
    gradci!::Opt{Function}
    gradF::Function
    gradF!::Opt{Function}
    gradFq::Function
    gradFq!::Opt{Function}
    hessf::Function
    hessf!::Opt{Function}
    hessce::Function
    hessce!::Opt{Function}
    hessci::Function
    hessci!::Opt{Function}
    hessF::Function
    hessF!::Opt{Function}
    hessFq::Function
    hessFq!::Opt{Function}
    fdp::Function
    fdp!::Opt{Function}
    cedp::Function
    cedp!::Opt{Function}
    cidp::Function
    cidp!::Opt{Function}
    Fdp::Function
    Fdp!::Opt{Function}
    gradfdp::Function
    gradfdp!::Opt{Function}
    gradcedp::Function
    gradcedp!::Opt{Function}
    gradcidp::Function
    gradcidp!::Opt{Function}
    gradFdp::Function
    gradFdp!::Opt{Function}
end

struct MPCCJumpFixedFunctions
    jump_f::Function
    jump_f_pen::Function
    jump_ce::Vector{Function}
    jump_ci::Vector{Function}
    jump_F::Matrix{Function}
end

# Only std return functions for the moment
struct MPCCNewtonPenaltyFunctions
    ϕ::Function
    gradϕ::Function
    hessϕ::Function
end


text:
"""
    struct DijkstraState{T, U}

An [`AbstractPathState`](@ref) designed for Dijkstra shortest-paths calculations.
"""
struct DijkstraState{T <: Real,U <: Integer} <: AbstractPathState
    parents::Vector{U}
    dists::Vector{T}
    predecessors::Vector{Vector{U}}
    pathcounts::Vector{UInt64}
    closest_vertices::Vector{U}
end

"""
    dijkstra_shortest_paths(g, srcs, distmx=weights(g));

Perform [Dijkstra's algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)
on a graph, computing shortest distances between `srcs` and all other vertices.
Return a [`LightGraphs.DijkstraState`](@ref) that contains various traversal information.

### Optional Arguments
- `allpaths=false`: If true, returns a [`LightGraphs.DijkstraState`](@ref) that keeps track of all
predecessors of a given vertex.

### Performance
Use a matrix type for `distmx` that is implemented in [row-major matrix format](https://en.wikipedia.org/wiki/Row-_and_column-major_order)
for better run-time.
Eg. Set the type of `distmx` to `Transpose{Int64, SparseMatrixCSC{Int64,Int64}}`
instead of `SparseMatrixCSC{Int64,Int64}`.
"""
function dijkstra_shortest_paths(g::AbstractGraph,
    srcs::Vector{U},
    distmx::AbstractMatrix{T}=weights(g);
    allpaths=false,
    trackvertices=false
    ) where T <: Real where U <: Integer

    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    visited = zeros(Bool, nvg)

    pathcounts = zeros(UInt64, nvg)
    preds = fill(Vector{U}(), nvg)
    H = PriorityQueue{U,T}()
    # fill creates only one array.

    for src in srcs
        dists[src] = zero(T)
        visited[src] = true
        pathcounts[src] = 1
        H[src] = zero(T)
    end

    closest_vertices = Vector{U}()  # Maintains vertices in order of distances from source
    sizehint!(closest_vertices, nvg)

    while !isempty(H)
        u = dequeue!(H)

        if trackvertices
            push!(closest_vertices, u)
        end

        d = dists[u] # Cannot be typemax if `u` is in the queue
        for v in outneighbors(g, u)
            alt = d + distmx[u, v]

            if !visited[v]
                visited[v] = true
                dists[v] = alt
                parents[v] = u

                pathcounts[v] += pathcounts[u]
                if allpaths
                    preds[v] = [u;]
                end
                H[v] = alt
            elseif alt < dists[v]
                dists[v] = alt
                parents[v] = u
                #615
                pathcounts[v] = pathcounts[u]
                if allpaths
                    resize!(preds[v], 1)
                    preds[v][1] = u
                end
                H[v] = alt
            elseif alt == dists[v]
                pathcounts[v] += pathcounts[u]
                if allpaths
                    push!(preds[v], u)
                end
            end
        end
    end

    if trackvertices
        for s in vertices(g)
            if !visited[s]
                push!(closest_vertices, s)
            end
        end
    end

    for src in srcs
        pathcounts[src] = 1
        parents[src] = 0
        empty!(preds[src])
    end

    return DijkstraState{T,U}(parents, dists, preds, pathcounts, closest_vertices)
end

dijkstra_shortest_paths(g::AbstractGraph, src::Integer, distmx::AbstractMatrix=weights(g); allpaths=false, trackvertices=false) =
dijkstra_shortest_paths(g, [src;], distmx; allpaths=allpaths, trackvertices=trackvertices)


text:
function lbfgs(fcn::Function,
               x0::Vector,
               args...;
               xtol::Real           = 1e-32, # default from Optim.jl
               ftol::Float64        = 1e-14, # Default from csminwel
               grtol::Real          = 1e-8,  # default from Optim.jl
               iterations::Int      = 1000,
               store_trace::Bool    = false,
               show_trace::Bool     = false,
               extended_trace::Bool = false,
               verbose::Symbol      = :none,
               rng::AbstractRNG     = MersenneTwister(),
               autodiff::Bool       = false,
               kwargs...)
    if autodiff
        Optim.optimize(fcn, x0, LBFGS(),
                       Optim.Options(g_tol = grtol, f_tol = ftol, x_tol = xtol,
                                     iterations = iterations, store_trace = store_trace,
                                     show_trace = show_trace,
                                     extended_trace = extended_trace))
    else
        Optim.optimize(fcn, x0, LBFGS(), autodiff=:forward,
                       Optim.Options(g_tol = grtol, f_tol = ftol, x_tol = xtol,
                                     iterations = iterations, store_trace = store_trace,
                                     show_trace = show_trace,
                                     extended_trace = extended_trace))
    end
end


text:
using DataDepsGenerators
using Test

using ReferenceTests

@testset "KNB test" begin
    # TODO: Findout why KBN doesn't work, I keep getting SSL errors.

#    @test_reference "../references/KNB/KnowledgeNetworkforBiocomplexity Forage Fish.txt" generate(KnowledgeNetworkforBiocomplexity(), "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/doi:10.5063/F1T43R7N")
end


text:
using ADCME
using PyPlot

db = Database("simulation.db")
c = execute(db, "select size, forward, backward from acoustic where level = 'cpu_tensor' and size <= 600 order by size")
c = collect(c)
s = [x[1] for x in c]
t1 = [x[2] for x in c]
t2 = [x[3] for x in c]
close("all")
loglog(s, t1, "*-", label = "Forward Computation")
loglog(s, t2, "*-", label = "Gradient Backpropagation")
savefig("result.png")


text:
# Create mock objects used in testing to check that the correct calls are made.
#
# This isn't used by the runtime for WebSockets. This is only used in testing. However, it is also
# used by the DandelionSlack project, for which this WebSocket package was built. All of this should
# ideally be in a separate package, and only used from tests. However, I'd prefer to keep the
# prerequisites for this package to a minimum, at least for now.

export @mock, @mockfunction, @expect, Throws, MockExpectationException,
       MockCall, ValueMatcher, mock_match, TypeMatcher, AbstractMatcher,
       mock_match

"Tells a mock function that it should throw an exception."
struct Throws
    ex::Exception
end

mock_action(action::Throws) = throw(action.ex)
mock_action(x::Any) = x

"An expected function call, along with expected arguments, and an action it should perform."
mutable struct MockCall
    sym::Symbol
    args::Vector{Any}
    action::Any
end

"An exception thrown when an unexpected argument was found or a function was called."
mutable struct MockExpectationException <: Exception
    s::AbstractString
end

MockExpectationException(expected::Any, actual::Any) =
    MockExpectationException("Expected $expected, but got $actual")
show(m::MockExpectationException) = "MockExpectationException: $(m.s)"

"""
Abstract type for objects that match actual arguments against expected arguments. For instance, we
sometimes don't care what value is provided a mock function, only what type it is. Also, this can
be used to parse an actual argument as JSON and compare the resulting object, rather than relying on
a string comparison.

All `AbstractMatcher` types `T` must define a function `mock_match(::T, v::Any)`.
"""
abstract type AbstractMatcher end

"Simply match a value by equality. This is the default."
struct ValueMatcher <: AbstractMatcher
    value::Any
end

mock_match(m::ValueMatcher, v::Any) = m.value == v || throw(MockExpectationException(m.value, v))
show(m::ValueMatcher) = show(m.value)

"Match a value by checking its type only."
struct TypeMatcher <: AbstractMatcher
    typ::DataType
end

mock_match(m::TypeMatcher, v::Any) =
    isa(v, m.typ) || throw(MockExpectationException(m.typ, typeof(v)))
show(m::TypeMatcher) = "TypeMatcher($(m.typ))"

"Define a mock type, given a symbol for its name and a type from which it should inherit."
macro mock(mock_type::Symbol, abstract_type::Symbol)
    esc(
        quote
            # Define the mock type. It should contain a list of the calls it expects.
            struct $mock_type <: $abstract_type
                calls::Vector{MockCall}

                $mock_type() = new([])
            end

            # Define a function that checks that all functions were in fact called.
            function check(m::$mock_type)
                try
                    if !isempty(m.calls)
                        throw(MockExpectationException("Missing calls: $(m.calls)"))
                    end
                finally
                    empty!(m.calls)
                end
            end
        end
    )
end

"Define a name for an argument, if no name is provided."
function name_args_(fdef::Expr)
    function_sym = fdef.args[1]
    for a in fdef.args[2:end]
        if length(a.args) == 1
            unshift!(a.args, gensym(function_sym))
        end
    end
end

"""
Define one or more functions that mock the behaviour the object should have.

The functions defined here are the functions called by the code under test. The functions check that
their arguments match those that are expected, and that the mock action is performed.
"""
macro mockfunction(t::Symbol, fdefs...)
    funs = Vector{Expr}()
    for fdef in fdefs
        # Name all arguments in the function definition.
        name_args_(fdef)
        function_sym = fdef.args[1]
        params = [x.args[1] for x in fdef.args[2:end]]

        # A block of code that each function should have.
        # Basically, we collect all arguments. We get the next expected call, and then we match each
        # arguments against the expected value.
        check_call_block = quote
            args = vcat($(params...))
            called_sym = $(Expr(:quote, function_sym))

            if isempty(($t).calls)
                throw(MockExpectationException("Unexpected call $called_sym, with args $args"))
            end

            mock_call = shift!(($t).calls)

            # Check that the expected function is called.
            if called_sym != mock_call.sym
                throw(MockExpectationException(
                    "Expected call $(mock_call.sym)($(mock_call.args)) but got $(called_sym)($args)"))
            end

            # Check that we have the right number of args.
            if length(mock_call.args) != length(args)
                throw(MockExpectationException(
                    "Expected $(length(mock_call.args)) arguments, but got $(length(args)):" *
                    "Expected $(mock_call.args), but got $(args)"))
            end

            # Match each actual argument against the expected.
            for i in range(1, length(mock_call.args))
                matcher = mock_call.args[i]
                arg = args[i]
                DandelionWebSockets.mock_match(matcher, arg)
            end

            # Perform the action set with the @expect macro.
            return DandelionWebSockets.mock_action(mock_call.action)
        end

        fun = Expr(:function)
        push!(fun.args, fdef)
        push!(fun.args, check_call_block)

        push!(funs, fun)
    end

    esc(quote
        $(funs...)
    end)
end

# For consistency, each expected value is wrapped with a matcher, even when that matcher merely
# checks for equality.
to_matcher(x::AbstractMatcher) = x
to_matcher(x::Any) = ValueMatcher(x)

"Expect a call with some arguments, and perform an action when that call has been made.."
macro expect(mock_object::Symbol, fcall::Expr, rvs...)
    sym = fcall.args[1]
    args = Vector{Any}(fcall.args[2:end])

    rv = nothing
    if !isempty(rvs)
        rv = rvs[1]
    end

    quote
        qargs = vcat($([esc(x) for x in args]...))
        match_args = [to_matcher(x) for x in qargs]
        push!(($(esc(mock_object))).calls, MockCall($(Expr(:quote, sym)), match_args, $(esc(rv))))
    end
end


text:
const CPC = Ptr{Cvoid}
const CPCType = Cstring

abstract type AbstractPC{T} end

mutable struct PC{T} <: AbstractPC{T}
    ptr::Ptr{Cvoid}
end

scalartype(::AbstractPC{T}) where {T} = T

@for_libpetsc begin

    function PC{$PetscScalar}(comm::MPI.Comm)
        pc = PC{$PetscScalar}(C_NULL)
        @chk ccall((:PCCreate, $libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CPC}), comm, pc)
        finalizer(destroy, pc)
        return pc
    end

    function PC(ksp::KSP{$PetscScalar})
        pc = PC{$PetscScalar}(C_NULL)
        @chk ccall((:KSPGetPC, $libpetsc), PetscErrorCode, (CKSP, Ptr{CPC}), ksp, pc)
        incref(pc) # need to manually increment the reference counter
        finalizer(destroy, pc)
        return pc
    end

    function destroy(pc::AbstractPC{$PetscScalar})
        finalized($petsclib) ||
        @chk ccall((:PCDestroy, $libpetsc), PetscErrorCode, (Ptr{CPC},), pc)
        return nothing
    end

    function settype!(pc::AbstractPC{$PetscScalar}, pctype::String)
        @chk ccall((:PCSetType, $libpetsc), PetscErrorCode, (CPC, Cstring), pc, pctype)
        return nothing
    end

    function setpc!(ksp::KSP{$PetscScalar}, pc::AbstractPC{$PetscScalar})
        @chk ccall((:KSPSetPC, $libpetsc), PetscErrorCode, (CKSP, CPC), ksp, pc)
        return nothing
    end

    function gettype(pc::AbstractPC{$PetscScalar})
        t_r = Ref{CPCType}()
        @chk ccall((:PCGetType, $libpetsc), PetscErrorCode, (CPC, Ptr{CPCType}),  pc, t_r)
        return unsafe_string(t_r[])
    end

    function view(pc::AbstractPC{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, getcomm(pc)))
        @chk ccall((:PCView, $libpetsc), PetscErrorCode,
                    (CPC, CPetscViewer),
                pc, viewer);
        return nothing
    end

end

Base.show(io::IO, pc::AbstractPC) = _show(io, pc)