text:
#=

# Debugging applications

## Debug API usage with validation layers

When a Vulkan application crashes, it can be due to a wrong API usage. Validation layers (provided by Khronos) can be enabled which check all API calls to make sure their invocation and execution is in accordance with the Vulkan specification.
The validation layers are typically used with a debug callback, which prints out messages reported by the validation layers.

First, let's create an instance with the validations layers and the `VK_EXT_debug_utils` extension enabled:

```julia
using Vulkan

instance = Instance(["VK_LAYER_KHRONOS_validation"], ["VK_EXT_debug_utils"])
```

Then, we will define a C function to be called when messages are received. We use the [`default_debug_callback`](@ref) provided by the Vulkan.jl:

```julia
const debug_callback_c = @cfunction(default_debug_callback, UInt32, (DebugUtilsMessageSeverityFlagEXT, DebugUtilsMessageTypeFlagEXT, Ptr{VkCore.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))
```

You need then need to define which message types and logging levels you want to include. Note that only setting e.g. `DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT` will not enable you the other levels, you need to set them all explicitly.

```julia
# Include all severity messages. You can usually leave out
# the verbose and info severities for less noise in the output.
message_severity = |(
    DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
)

# Include all message types.
message_type = |(
    DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT,
    DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT,
    DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT,
)
```

You can now create your debug callback:

```julia
messenger = DebugUtilsMessengerEXT(instance, message_severity, message_type, debug_callback_c)
```

However, make sure that the messenger stays alive as long as it is used. You should typically put it in a config next to the `instance` (for scripts, a global variable works as well).

!!! tip
    To debug the initialization of the instance itself, you can create the DebugUtilsMessengerEXTCreateInfo manually and include it in the `next` parameter of the `Instance`:
    ```julia
    create_info = DebugUtilsMessengerCreateInfoEXT(message_severity, message_type, debug_callback_c)
    instance = Instance(["VK_LAYER_KHRONOS_validation"], ["VK_EXT_debug_utils"]; next = create_info)
    messenger = DebugUtilsMessengerEXT(instance, create_info)
    ```

=#

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
println("**********Testing Distribution**********")

facts("Creating a Distribution") do
  fv = FeatureVector()
  c = Cluster()
  ds = DataSet()
  d1 = Distribution(fv)
  d2 = Distribution(c)
  d3 = Distribution(ds)

  @fact typeof(d1) => Distribution{FeatureVector}
  @fact typeof(d2) => Distribution{Cluster}
  @fact typeof(d3) => Distribution{DataSet}
end

facts("If sent in with a FeatureVector, d.total is set") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact d1.total => 7
end

facts("Get probability of seeing feature in Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  value = probability(d1,"word")

  @fact value => 4/7
end

facts("Get features of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact features(d1) => keys(d1.space)
end

facts("Check isempty on a Distribution") do
  d1 = Distribution{FeatureVector}()
  d2 = Distribution{Cluster}()
  d3 = Distribution{DataSet}()

  @fact isempty(d1) => isempty(d1.space)
  @fact isempty(d2) => isempty(d2.space)
  @fact isempty(d3) => isempty(d3.space)
end

facts("entropy(Distribution) returns 0 if empty Distribution") do
  d1 = Distribution{FeatureVector}()

  @fact entropy(d1) => 0
end

facts("entropy(Distribution) returns entropy of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact entropy(d1) => 0.9852281360342515
end

facts("entropy of different Distribution types") do
  fv1 = FeatureVector(["c"=>1,"b"=>1,"a"=>1,"d"=>1])
  fv2 = FeatureVector(["f"=>1,"g"=>1,"e"=>1,"h"=>1])
  fva = FeatureVector(["1"=>1,"2"=>1,"3"=>1,"4"=>1])
  fvb = FeatureVector(["8"=>1,"7"=>1,"6"=>1,"5"=>1])

  c1 = Cluster()
  c2 = Cluster()
  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c2["fva"] = fva
  c2["fvb"] = fvb

  ds = DataSet()
  ds["c1"] = c1
  ds["c2"] = c2

  fv1_dist = Distribution(fv1)
  fv2_dist = Distribution(fv2)
  fva_dist = Distribution(fva)
  fvb_dist = Distribution(fvb)
  c1_dist = Distribution(c1)
  c2_dist = Distribution(c2)
  ds_dist = Distribution(ds)

  @fact entropy(fv1_dist) => 2
  @fact entropy(fv2_dist) => 2
  @fact entropy(fva_dist) => 2
  @fact entropy(fvb_dist) => 2
  @fact entropy(c1_dist) => 3
  @fact entropy(c2_dist) => 3
  @fact entropy(ds_dist) => 4
end

facts("entropy of certain fv or clust") do
  fv1 = FeatureVector(["c"=>1,"b"=>1,"a"=>1,"d"=>1])
  fv2 = FeatureVector(["f"=>1,"g"=>1,"e"=>1,"h"=>1])
  fva = FeatureVector(["1"=>1,"2"=>1,"3"=>1,"4"=>1])
  fvb = FeatureVector(["8"=>1,"7"=>1,"6"=>1,"5"=>1])

  c1 = Cluster()
  c2 = Cluster()
  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c2["fva"] = fva
  c2["fvb"] = fvb

  ds = DataSet()
  ds["c1"] = c1
  ds["c2"] = c2

  c1_dist = Distribution(c1)
  c2_dist = Distribution(c2)
  ds_dist = Distribution(ds)

  @fact fv_entropy(c1_dist) => 1
  @fact fv_entropy(ds_dist) => 2
  @fact fv_entropy(ds_dist, "c1") => 1
  @fact clust_entropy(ds_dist) => 1
end

facts("Testing prob_* functions") do
  fv1 = FeatureVector(["c"=>1,"b"=>1,"a"=>1,"d"=>1])
  fv2 = FeatureVector(["f"=>1,"g"=>1,"e"=>1,"h"=>1])
  fva = FeatureVector(["1"=>1,"2"=>1,"3"=>1,"4"=>1])
  fvb = FeatureVector(["8"=>1,"7"=>1,"6"=>1,"5"=>1])

  c1 = Cluster()
  c2 = Cluster()
  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c2["fva"] = fva
  c2["fvb"] = fvb

  ds = DataSet()
  ds["c1"] = c1
  ds["c2"] = c2

  c1_dist = Distribution(c1)
  ds_dist = Distribution(ds)

  @fact prob_fv(c1_dist, "fv1") => 4/8
  @fact prob_fv(ds_dist, "c1", "fv1") => 4/16
  @fact prob_cl(ds_dist, "c1") => 8/16
end

facts("Testing cond_prob_* functions") do
  fv1 = FeatureVector(["c"=>1,"b"=>1,"a"=>1,"d"=>1])
  fv2 = FeatureVector(["f"=>1,"g"=>1,"e"=>1,"h"=>1])
  fva = FeatureVector(["1"=>1,"2"=>1,"3"=>1,"4"=>1])
  fvb = FeatureVector(["8"=>1,"7"=>1,"6"=>1,"5"=>1])

  c1 = Cluster()
  c2 = Cluster()
  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c2["fva"] = fva
  c2["fvb"] = fvb

  ds = DataSet()
  ds["c1"] = c1
  ds["c2"] = c2

  c1_dist = Distribution(c1)
  ds_dist = Distribution(ds)

  @fact cond_prob_f_given_fv(c1_dist, "fv1", "d") => 1/4
  @fact cond_prob_f_given_fv(ds_dist, "c2", "fva", "3") => 1/4
  @fact cond_prob_f_given_clust(ds_dist, "c1", "f") => 1/8
  @fact cond_prob_fv_given_clust(ds_dist, "c2", "fva") => 4/8
end

facts("info_gain(Distribution,Distribution) returns 0 if empty Distributions") do
  d1 = Distribution{FeatureVector}()
  d2 = Distribution{FeatureVector}()

  @fact info_gain(d1,d2) => 0
end

facts("ld_info_gain rules out similar features over all decades and outliers that wouldn't tell us much") do
  fv1 = FeatureVector(["the"=>5, "c1"=>5])
  fv2 = FeatureVector(["the"=>5, "c1"=>5, "romeo"=>100])
  fv3 = FeatureVector(["the"=>5, "c1"=>5])
  fv4 = FeatureVector(["the"=>5, "c1"=>5])
  fv5 = FeatureVector(["the"=>5, "c2"=>5])
  fv6 = FeatureVector(["the"=>5, "c2"=>5])
  fv7 = FeatureVector(["the"=>5, "c2"=>5])
  fv8 = FeatureVector(["the"=>5, "c2"=>5])
  fv9 = FeatureVector(["the"=>5, "c3"=>5])
  fv10 = FeatureVector(["the"=>5, "c3"=>5])
  fv11 = FeatureVector(["the"=>5, "c3"=>5])
  fv12 = FeatureVector(["the"=>5, "c3"=>5])
  fv13 = FeatureVector(["the"=>5, "c4"=>5])
  fv14 = FeatureVector(["the"=>5, "c4"=>5])
  fv15 = FeatureVector(["the"=>5, "c4"=>5])
  fv16 = FeatureVector(["the"=>5, "c4"=>5])

  c1 = Cluster()
  c2 = Cluster()
  c3 = Cluster()
  c4 = Cluster()

  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c1["fv3"] = fv3
  c1["fv4"] = fv4
  c2["fv5"] = fv5
  c2["fv6"] = fv6
  c2["fv7"] = fv7
  c2["fv8"] = fv8
  c3["fv9"] = fv9
  c3["fv10"] = fv10
  c3["fv11"] = fv11
  c3["fv12"] = fv12
  c4["fv13"] = fv13
  c4["fv14"] = fv14
  c4["fv15"] = fv15
  c4["fv16"] = fv16

  ds = DataSet()

  ds["c1"] = c1
  ds["c2"] = c2
  ds["c3"] = c3
  ds["c4"] = c4

  d = Distribution(ds)

  for feature in keys(d.space.vector_sum)
    if feature == "the" || feature == "romeo"
      @fact ld_info_gain(d,feature) => 0.0
    else
      @fact ld_info_gain(d,feature) => 1.0
    end
  end
end

facts("info_gain(Distribution,Distribution) returns info_gain of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  dict2 = ["∂" => .1, "happy" => .3]
  fv2 = FeatureVector(dict2)
  d2 = Distribution(fv2)

  @fact info_gain(d1,d2) => 0.17395001157511858
end

facts("Smoothing FeatureVector has correct probabilties before and after") do
  dict1 = ["word" => 7, "another" => 13]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  #no smoothing
  @fact probability(d1,"word") => (7/20)
  @fact probability(d1,"another") => (13/20)
  @fact probability(d1,"unk") => (0/20)

  #delta smoothing
  delta_smoothing!(d1)
  @fact probability(d1,"word") => (8/(3+20))
  @fact probability(d1,"another") => (14/(3+20))
  @fact probability(d1,"unk") => (1/(3+20))

  #good-turing smoothing
  dict3 = ["the" => 15, "of" => 16, "what" => 1, "a" => 1, "unique" => 1, "word" => 1, "twice" => 2, "two" => 2, "party" => 2, "three" => 3]
  fv3 = FeatureVector(dict3)
  d3 = Distribution(fv3)
  goodturing_smoothing!(d3)
  @fact probability(d3,"word") => (((1+1)*(3/4))/44)
  @fact probability(d3,"twice") => (((2+1)*(1/3))/44)
  @fact probability(d3,"three") => (((3+1)*(0/1))/44)
  @fact probability(d3,"unk") => (4/44)

  #removing smoothing
  remove_smoothing!(d1)
  @fact probability(d1,"word") => (7/20)
  @fact probability(d1,"another") => (13/20)
  @fact probability(d1,"unk") => (0/20)
end

ext: jl


text:
# const USE_HWLOC = parse(Bool, get(ENV, "CPUSUMMARY_HWLOC", "true"))
# if !USE_HWLOC
#   run(`$(Base.julia_cmd()) --project=$(Base.active_project()) -e'using CPUSummary; CPUSummary.use_hwloc(false)'`)
# end

using CPUSummary
using Test

@testset "CPUSummary.jl" begin
  @test @inferred(CPUSummary.sys_threads()) == Sys.CPU_THREADS::Int
  @test @inferred(CPUSummary.num_threads()) == Threads.nthreads()
  # if USE_HWLOC
  #   for i ∈ 1:4
  #     @test CPUSummary.redefine_cache(1) === nothing
  #   end
  # else
  #   @test !isdefined(CPUSummary, :redefine_cache)
  # end
end


text:
using LinearAlgebra
using StaticArrays

@inline function get_λ_η(α, β,θ, a)
    λ = -α * sin(θ)
    η = (α^2 - a^2)*cos(θ)^2+β^2
    return λ, η
end

@inline function get_up_um(λ, η, a)
    Δ_θ  = 1/2*(1-(η+λ^2)/a^2)
    up = Δ_θ + √(Δ_θ^2 +η/a^2)
    um = Δ_θ - √(Δ_θ^2 +η/a^2)
    return up, um
end

function radialroots(λ::Real, η::Real, a)
    T = promote_type(typeof(λ), typeof(η))
    return radialroots(Complex{T}(λ), Complex{T}(η), a)
end

"""
    $(SIGNATURES)
Finds the radial roots using the of the geodesic with
energy-scaled angular momentum `λ` and carter constant `η`.
"""
@fastmath function radialroots(λ::Complex, η::Complex, a)
    A = a^2 - η - λ^2
    B = 2*(η+(λ-a)^2)
    C = -a^2 * η
    P = -A^2 / 12 - C
    Q = -A/3 * ((A/6)^2-C)-B^2/8
    H = -9*Q + √(12*P^3 + 81*Q^2)
    z = √((-2*(3^(1/3)*P)+2^(1/3)*H^(2/3))/(2*6^(2/3)*H^(1/3)) - A/6)
    discp = -A/2 - z^2 + B/(4*z)
    dp = sqrt(discp)
    discm = discp - B/(2*z)
    dm = sqrt(discm)
    r1 = -z - dp
    r2 = -z + dp
    r3 =  z - dm
    r4 =  z + dm
    #roots = SVector(r1,r2,r3,r4)
    #ind = sortperm(real.(roots))
    #println(real(λ), " ", real(η))
    return r1,r2,r3,r4#Tuple(roots[ind])
end

function get_Λ(β, χ)
    s,c = sincos(χ)
    γ = 1/sqrt(1-β^2)

    Λ00 = γ
    Λ01 = -γ*β*c
    Λ03 = -γ*β*s
    Λ11 = (γ-1)*c^2 + 1
    Λ13 = (γ-1)*s*c
    Λ33 = (γ-1)*s^2 + 1
    z = zero(typeof(γ))
    o =  one(typeof(γ))
    return @SMatrix[Λ00 Λ01  z  Λ03;
                    Λ01 Λ11  z  Λ13;
                    z   z    o  z
                    Λ03 Λ13  z  Λ33
                   ]
end


text:
@enum EventKind EVENT_INTERVAL EVENT_POINT

struct EventID
    value::UInt64
end

# TODO: wrap?
const LocationID = UInt64

struct Event
    eventkind::EventKind
    # RawEventID
    eventid::EventID  # mixing localid + threadid
    threadid::typeof(Threads.threadid())
    # LocationProxy
    locationid::UInt64
    # IntervalRecord
    approx_taskid::UInt
    start::TimeNS
    stop::Union{TimeNS,Nothing,Missing}
end

struct Location
    # LocationProxy
    locationid::UInt64
    tag::Union{Symbol,Nothing}
    line::Int
    file::Union{Symbol,Nothing}
    _module::Module
end

event_table() = append_events!(Event[])

location_table() = append_locations!(Location[])

function foreach_event(f, record_stores = INTERVAL_RECORD_STORES)
    for records in record_stores
        node = records.head
        block = node.block
        while true
            bound = node.i - 1
            node = node.next
            if node === nothing
                for i in 1:bound
                    x = @inbounds block[i]
                    f(x)
                end
                break
            end
            for i in eachindex(block)
                x = @inbounds block[i]
                f(x)
            end
            block = node.block
        end
    end
end

function eventkind_from_stop(stop::TimeNS)
    if stop === TIME_NOT_AVAILABLE
        return EVENT_POINT
    else
        return EVENT_INTERVAL
    end
end

EventKind(record::IntervalRecord) = eventkind_from_stop(record.stop)

# TODO: is it better to do this inside RawEventID (and get rid of it)?
function EventID(event::RawEventID)
    threadid = event.threadid
    localid = event.localid
    upper = bswap(UInt64(threadid - 1))
    lower = localid % UInt64
    return EventID(upper | lower)
end

function Event(record::IntervalRecord)
    Event(
        EventKind(record),
        EventID(record.event),
        record.event.threadid,
        record.location,
        record.approx_taskid,
        record.start,
        interpret_time(record.stop),
    )
end

function append_events!(table)
    foreach_event() do x
        push!(table, Event(x))
    end
    return table
end

function Location(record::IntervalRecord)
    local location
    GC.@preserve LOCATIONS_ROOT let uints
        uints = unsafe_wrap(Array, Ptr{UInt}(pointer(LOCATIONS_ROOT)), size(LOCATIONS_ROOT))
        # @assert record.location in uints
        location = unsafe_pointer_to_objref(Ptr{Cvoid}(record.location))::LocationProxy
    end
    return Location(
        record.location,
        location.tag === NOTAG ? nothing : location.tag,
        location.line.line,
        location.line.file,
        location._module,
    )
end

"""
    EventTracker.append_locations!(table)
"""
function append_locations!(table)
    seen = Set{UInt}()
    foreach_event() do record::IntervalRecord
        n0 = length(seen)
        push!(seen, record.location)
        n0 == length(seen) && return  # a bit faster than `in`
        push!(table, Location(record))
    end
    return table
end


text:
using Test, NiLangCore

@testset "basic" begin
    @test ~(~sin) === sin
    @test ~(~typeof(sin)) === typeof(sin)
    @test isreflexive(XorEq(NiLangCore.logical_or))
    println(XorEq(*))
    println(PlusEq(+))
    println(MinusEq(-))
    println(MulEq(*))
    println(DivEq(/))
end

@static if VERSION > v"1.5.100"
@testset "composite function" begin
    @i function f1(x)
        x.:1 += x.:2
    end
    @i function f2(x)
        x.:2 += cos(x.:1)
    end
    @i function f3(x)
        x.:1 ↔ x.:2
    end
    x = (2.0, 3.0)
    y = (f3∘f2∘f1)(x)
    z = (~(f3∘f2∘f1))(y)
    @show x, z
    @test all(x .≈ z)
end
end


text:
## Confidence intervals.

"""
Compute confidence intervals for the input realisations of the time series.

  * `matrix`: one scenario per year. Hypothesis: scenarios are only *one*-year-long.
  * `level`: the confidence level at which the confidence interval should be computed.
  * `ci_period`: on which the confidence intervals must be computed. Only recognised periods:
    * `week` (52 elements per line in `matrix`),
    * `day` (365 elements per line in `matrix`).
  * `method`: the method to use to compute the confidence intervals:
    * `:TStudent` (default): a standard statistical method is used (t-Student hypothesis testing).
      Those confidence intervals are computed for each step independently, with no correlation considered
      between the time steps,
    * `Adam2015`: a physically-based statistical method, based on a sequence of probability distributions.
      This implementation is based on http://hdl.handle.net/2268/142572 without the flooding part.

It is impossible to compute intervals on time scales smaller than the input data (i.e., if `ci_period` is `day`,
but the data is only available weekly). The output will always be at the same discretisation step as the input
scenarios, which may involve copying of the confidence interval values.

TODO:

  * implement other methods, including with correlation.
  * for Adam2015: implement flooding. (Should impact the mean and high part of the CI, not the low?)
  * replace the hard-coded part with a hash table. Keep hard-coded doc for built-in generators, though they are added by "standard" ways (such as addgenerator!() function).
"""
function computeCI_sub(r::River, level::Float64=.95; ci_period::Period=Week(1), method::Symbol=:TStudent)
  TS = r.scenarios
  period = getPeriod(r)
  if period == Week(1) && ci_period == Day(1)
    error("Cannot compute confidence intervals on inexistent data.")
  end

  if period != Week(1) && period != Day(1)
    error("Unrecognised time series period. Only allowable: `week` and `day`.")
  end
  if ci_period != Week(1) && ci_period != Day(1)
    error("Unrecognised confidence interval period. Only allowable: `week` and `day`.")
  end

  # Merge time series as required and make a matrix of it.
  if period != ci_period
    period_f = (period == Day(1)) ? day : (period == Week(1)) ? week : error("Unrecognised period.")
    TS = [collapse(ts, values -> sum(values), period=period_f) for ts in TS] # TODO: test this line hard. Doesn't work properly if period == ci_period (last value much higher than the others).
  end

  # Actual computation.
  # TODO: replace this with a hash table symbol <-> function to make it easier to extend? (Even outside this library!)
  if method == :TStudent
    return computeCI_TStudent(r, level, ci_period=ci_period)
  elseif method == :Adam2015
    return computeCI_Adam2015(r, level)
  end
end

function computeCI(r::River, level::Float64; kwargs...)
  # Sizes? means: (T,). confidence_intervals: (T, 2).
  means, ci_low, ci_high = computeCI_sub(r, level; kwargs...)

  r_mean = copy(r, scenarios=[means])
  r_low  = copy(r, scenarios=[ci_low])
  r_high = copy(r, scenarios=[ci_high])
  return r_mean, r_low, r_high
end

function computeCI(reservoir::Reservoir, level::Float64; ci_period::Period=Week, method::Symbol=:TStudent)
  lnr = [computeCI(r, level, ci_period=ci_period, method=method) for r in naturalRivers(reservoir)]
  ldr = [computeCI(r, level, ci_period=ci_period, method=method) for r in divertedRivers(reservoir)]

  rivers_mean = NaturalRiver[r[1] for r in lnr]
  rivers_low  = NaturalRiver[r[2] for r in lnr]
  rivers_high = NaturalRiver[r[3] for r in lnr]
  divertedRivers_mean = DivertedRiver[r[1] for r in ldr]
  divertedRivers_low  = DivertedRiver[r[2] for r in ldr]
  divertedRivers_high = DivertedRiver[r[3] for r in ldr]

  variantBase = ""
  if reservoir.variant != "vanilla"
    variantBase = reservoir.variant * " "
  end
  variantBase *= "CI@" * string(level * 100) * "% "

  # TODO: add empirical method? (95% of scenarios within bounds, in other words: no real statistics.)
  if method == :TStudent
    techniqueSuffix = "t-Student"
  elseif method == :Adam2015
    techniqueSuffix = "Adam 2015"
  else
    error("Unknown CI technique.")
  end

  reservoir_mean = copy(reservoir, variant=variantBase * "(mean)" * ", " * techniqueSuffix, rivers_in=rivers_mean, rivers_diverted=divertedRivers_mean)
  reservoir_low  = copy(reservoir, variant=variantBase *  "(low)" * ", " * techniqueSuffix, rivers_in=rivers_low,  rivers_diverted=divertedRivers_low)
  reservoir_high = copy(reservoir, variant=variantBase * "(high)" * ", " * techniqueSuffix, rivers_in=rivers_high, rivers_diverted=divertedRivers_high)
  return reservoir_mean, reservoir_low, reservoir_high
end

## Scenario generation from statistics.
"""
Generate a new `River` object based on the given `river`, with its scenarios replaced by newly generated ones.
The generation is parameterised by a number of scenarios to generate `n`, a periodicity for the output `period`,
a generation `method`, and a scenario duration `scenarioDuration` (which must be in years). The parameter `iterator`
can be set to `true` to have an iterator over scenarios (a full reservoir or river being computed at a time) when
multiple ones are requested.

The following methods are implemented:

  * `Adam2015`: a physically-based statistical method, based on a sequence of probability distributions.

Implement new scenario generation functions
===========================================

Mandatory interface to make a new scenario generation function:
    `f(r::River, period::Period, currentYear::Int; kwargs...)`
This function then generates a new scenario for the given `river` for one hydrological year (`currentYear` is only used
to produce good-looking time series objects). The generated data should have a periodicity of `period`, i.e. a value
is output each `period` from the beginning of the hydrological year. The generator might use more parameters than those
ones: they are given as keyword arguments (`kwargs`) through this function.

Steps:

  * create a function with the previous prototype
  * decide for a new keyword (such as :Adam2015); if the same statistical principles are used elsewhere,
    the same symbol should be used
  * modify the function `generateScenarios` to recognise the new keyword and associate the new function
  * update the documentation above

TODO: replace the hard-coded part with a hash table (store the function and the suffix in the same place). Keep hard-coded doc for built-in generators, though they are added by "standard" ways (such as addgenerator!() function).
"""
immutable RiverScenarioIterator
  r::River
  n::Int
  period::Period
  method::Symbol
  scenarioDuration::Period
  kwargs
end
Base.start(::RiverScenarioIterator) = 1
Base.next(i::RiverScenarioIterator, s::Int) = (generateScenarios(i.r, 1, i.period, i.method, i.scenarioDuration; i.kwargs...), s + 1)
Base.done(i::RiverScenarioIterator, s::Int) = s >= i.n
Base.eltype(i::RiverScenarioIterator) = River
Base.length(i::RiverScenarioIterator) = i.n
eachscenario(r::RiverScenarioIterator) = r

function generateScenarios(r::River, n::Int=1, period::Period=Day(1), method::Symbol=:Adam2015, scenarioDuration::Period=Year(1); kwargs...)
  @assert ! Dates.periodisless(scenarioDuration, Year(1)) "Scenario generation only allows generation for time periods that are multiple of years."

  if method == :Adam2015
    f = generateScenarios_Adam2015
  else
    error("Unrecognised method: " * string(method))
  end

  hasIterator = length(filter(e -> e[1] == :iterator, kwargs)) > 0 && filter(e -> e[1] == :iterator, kwargs)[1][2]
  hasIterator && deleteat!(kwargs, find(e -> e[1] == :iterator, kwargs))
  if n == 1 || hasIterator
    # No iterator here.
    years_per_scenario = Int(floor(Dates.days(scenarioDuration) / Dates.days(Year(1))))
    total_n = n * years_per_scenario
    generated = TimeSeries.TimeArray{Float64,1,DateTime,Array{Float64,1}}[let pair = f(r, period, y; kwargs...)
                    TimeArray(pair[1], pair[2], colnames(r.scenarios[1]), meta(r.scenarios[1]))
                  end for y in collect(NonLeapYearIterator(1, total_n))]
    return scenarioGeneration(copy(r, scenarios=generated), :Concatenate, years_per_scenario)
  else
    # Provide an iterator.
    return RiverScenarioIterator(r, n, period, method, scenarioDuration, kwargs)
  end
end

immutable ReservoirScenarioIterator
  reservoir::Reservoir
  n::Int
  period::Period
  method::Symbol
  scenarioDuration::Period
  kwargs
end
Base.start(::ReservoirScenarioIterator) = 1
Base.next(i::ReservoirScenarioIterator, s::Int) = (generateScenarios(i.reservoir, 1, i.period, i.method, i.scenarioDuration; i.kwargs...), s + 1)
Base.done(i::ReservoirScenarioIterator, s::Int) = s >= i.n
Base.eltype(i::ReservoirScenarioIterator) = Reservoir
Base.length(i::ReservoirScenarioIterator) = i.n
eachscenario(r::ReservoirScenarioIterator) = r

function generateScenarios(reservoir::Reservoir, n::Int=1, period::Period=Day(1), method::Symbol=:Adam2015, scenarioDuration::Period=Year(1); kwargs...)
  variant = ""
  if reservoir.variant != "vanilla"
    variant = reservoir.variant * " "
  end
  variant *= "generated scenarios"

  if method == :Adam2015
    techniqueSuffix = "Adam 2015"
  else
    error("Unknown scenario generation technique.")
  end
  variant *= ", " * techniqueSuffix

  hasIterator = length(filter(e -> e[1] == :iterator, kwargs)) > 0 && filter(e -> e[1] == :iterator, kwargs)[1][2]
  if n == 1 || hasIterator
    # No iterator here.
    lnr = NaturalRiver[generateScenarios(r, n, period, method, scenarioDuration; kwargs...) for r in naturalRivers(reservoir)]
    ldr = DivertedRiver[generateScenarios(r, n, period, method, scenarioDuration; kwargs...) for r in divertedRivers(reservoir)]
    return copy(reservoir, variant=variant, rivers_in=lnr, rivers_diverted=ldr)
  else
    # Provide an iterator.
    return ReservoirScenarioIterator(reservoir, n, period, method, scenarioDuration, kwargs)
  end
end


text:
using Atom.Markdown

#add a variable and playing with it
messere = "pippo cicuta minatore saldatore"
list = []
#making a for cycle
for i in 1:2
    print("Hello how are you?? \n")
    print("Quante pila ave lu poccu? blabla??"); a = readline()
    append!(list, parse(Int, a))
    print("Lu poccu ave $list pilas")
    end

#Justify the text adding 70 spaces
rightjustify = (" "^65 * "monty")
print(rightjustify)

#printing a nice reticula
reticula = function()
    print("+ - - - - - +" ^2)
    print("|         |         |" ^5 )
    print("+ - - - - - +" ^2)
end

reticula()

using ThinkJulia
🐢 = Turtle()
@svg begin
    forward(🐢, 100)
end

🐢 = Turtle()
@svg begin
    for i in 1:4
        forward(🐢, 100)
        turn(🐢, -90)
    end
end

Encapsulation

The first exercise asks you to put your square-drawing code into a
function definition and then call the function, passing the turtle as a parameter.
Here is a solution:

function square(t)
    for i in 1:4
        forward(t, 100)
        turn(t, -90)
    end
end
🐢 = Turtle()
@svg begin
    square(🐢)
end

x = 4
if 0 < x < 10
    println("x is a positive single-digit number.")
end

function countdown(n)
    if n ≤ 0
        println("Blastoff!")
    else
        println(n, " ")
        countdown(n-1)
    end
end

countdown(3)

#reading keyboard input
text = readline()

#If you expect the user to type an integer, you can try to convert the return value to Int64:
println("What...is the airspeed velocity of an unladen swallow?"); speed = readline()
parse(Int64, speed)

x != y               # x is not equal to y
  x ≠ y                # (\ne TAB)
  x > y                # x is greater than y
  x < y                # x is less than y
  x >= y               # x is greater than or equal to y
  x ≥ y                # (\ge TAB)
  x <= y               # x is less than or equal to y
  x ≤ y                # (\le TAB)

  prefixes = "JKLMNOPQ"
  suffix = "ack"

  for letter in prefixes
      println(letter * suffix)
  end

 a = "sjhdsjfh"
 b = "pippolo"

 for lettera in a
     println(lettera * b)
 end

 word = "bananaaaaaaaaaaaaaaaaahhhhhhhhhhhhdfgsdsjhjahsjahajajajhsjahjs"
 counts = 0
 for letter in word
     if letter == 'a'
         global counts += 1
     end
 end
 println(counts)


text:
#takes a crsp daily file and makes it monthly
function makemonthly(selectionop::Function = (noops(x)=x),
    monthlyop::Function = (noopm(x)=x);
    crspnameout=PARAM[:crspmonthly],
    mainpanel::Bool=true,
    binprefix = "m_",
    computeivol::Bool = PARAM[:computeivol])

  local crsp::DataFrame

  #WARNING- crspd is NOT valid, below is only to create crspm
  crsp = prepcrsp(crsptype=PARAM[:crspdaily],
    csvextension = CSV_EXTENSION,
    incsvstream = IN_CSV_STREAM,
    binprefix = binprefix,
    refreshcrsp=PARAM[:refreshmonthly],
    crspcolumns=[PARAM[:monthlycrspcols]; keys(PARAM[:momentumstrategies]) |> collect],
    crsprequiredcolumns=PARAM[:dailycrsprequired]) do crspd

    @assert (PARAM[:crspfrequencysuffix] == "m") || (!mainpanel)
    sort!(crspd, [:permno, :date]) #cache the sort

    #computerangemeasures!(crspd, sorted=true) #not used, mainly for kyle's lamda
    crspd.dvol = computedollarvol(crspd)
    rename!(crspd, :vol=>:svol)
    #crspd.shares = computeshares(crspd)

    #NOTE: testing note- change to trailingsigmatest! to test vol routine
    #println(" checksum: $(sum(skipmissing(crspd.sigma)))")

    crspd = mergeff(crspd, Fgroup=:permno, datefrequency=:day)
    #crspd = dailyop(crspd)

    if !issorted(crspd, [:permno, :date])
      sort!(crspd, [:permno, :date])
    end

    crspm::DataFrame = makemonthly(selectionop, crspd; computeivol)
    validatedates(crspm.date |> unique |> sort!; frequency=:month, validateexactfrequency=true)
    sort!(crspm, [:permno, :date])
    validatesummaryintegrity(crspm, Fid=:monthid)

    if PARAM[:testmonthlycsv]

      #time series test output
      testpermno::Int = PARAM[:testpermnomult]
      #testcrspm::DataFrame = crspm[crspm.permno .% testpermno .== 0, :]
      testcrspm::DataFrame = crspm[crspm.permno .== 14541, :]
      #if size(testcrspm,1)>1_000 #prevents scenarios where the file is larger than needed
      #  testcrspm=testcrspm[1:1_000]
      #end

      #=permnos::Vector{Int} = unique(crspm.permno)[1:400:min(7000, end)]
      crspm.totest = falses(size(crspm,1))
      for scrspm ∈ groupby(crspm, :permno)
        (first(scrspm.permno) ∈ permnos) && (scrspm.totest .= true)
      end
      testcrspm = crspm[crspm.totest,:]
      select!(crspm, Not(:totest))=#
      rename!(testcrspm, :date=>:datem, :dvol=>:dvolm, :ivol=>:ivolm, :ret=>:retm, :price=>:pricem,
        :tridx=>:tridxm, :shares=>:sharesm, :sigma=>:sigmam, :mc=>:mcm, :plret=>:plretm,
        :exchcd=>:exchcdm, :siccd=>:siccdm,)
      rename!(testcrspm, (k->k=>Symbol(k, :m)).(keys(PARAM[:momentumstrategies])))
      #make sure we merge using the same df as used to create the monthly df
      testcrspout = innerjoin(testcrspm, crspd, on=[:permno, :monthid])
      sort!(testcrspout, [:permno, :monthid])
      @assert size(testcrspout,1) ≤ 30_000

      testcrspout |> CSV.write(
        "$(PARAM[:testpath])\\$(binprefix)dailyformonthly_$(PARAM[:crspdaily]).csv")
    end

    crspm = monthlyop(crspm)

    return crspm
  end #end do block

  #println(describe(crsp))
  return crsp
end

function debugmsg(hdr, df)
  #=meanret = mean(skipmissing(df.ret)) |> (x) -> round(x,sigdigits=12)
  meanprice = mean(skipmissing(df.price)) |> (x) -> round(x,sigdigits=12)
  meanvol = mean(skipmissing(df.dvol)) |> (x) -> round(x,sigdigits=12)
  meanpermno = mean(df.permno) |> (x) -> round(x,sigdigits=12)=#
  print("$hdr hdr: size(df): $(size(df))")
  for s ∈ propertynames(df)
    (!(eltype(df[!,s]) <: Union{Real, Missing})) && continue
    all(df[!,s] .=== missing) && continue
    occursin("placebo", string(s)) && continue
    print(" $s $(mean(skipmissing(df[!, s]))),")
  end
  print("\n")
end

#aggregates a daily dataframe into a monthly dataframe
function makemonthly(selectionop::Function, crspd::AbstractDataFrame;
  computeivol::Bool)

  #Assign each date to a monthid
  crspd.monthid = assignmonthid(crspd.date)
  crspd.retP1 = 1.0 .+ crspd.ret #convert to gross returns for geometric compounding
  crspd.absret = (abs).(crspd.ret)

  issorted(crspd, [:permno, :date]) || error("crspd must be sorted in makemonthly")
  computemonthend!(crspd)
  trailingsigma!(crspd, Fendpoint=:monthend)
  select!(crspd, Not(:monthend))

  #println(describe(crspd))

  computemonthend!(crspd)
  if computeivol
    regressatperiodend!(crspd,
      minpoints=PARAM[:sigmaminpoints],
      maxpoints=PARAM[:sigmamaxpoints], #this is mainly for technical reasons
      calendarinterval = PARAM[:sigmacalendarinterval],
      FXs = PARAM[:ivolfields],
      Fagg = :ivol,
      aggfunc = std,
      Fendpoint = :monthend,
      #=aggfunc = (lm::FMLM)->std(lm.ε)=#)

  else
    @info "ivol not computed"
    crspd.ivol = missings(Float64, nrow(crspd))
  end

  #sanity checks
  #debugmsg("dbg1", crspd)
  #drop rows with too much missing data
  @assert issorted(crspd, [:permno, :date])
  oldnrow = nrow(crspd)
  crspd::DataFrame = crspd[completecases(crspd, PARAM[:monthlycrsprequired]), :]
  @assert (size(crspd,1) ≥ oldnrow * 0.5) "monthlycrsprequired dropped more than 50% of rows!!!
    (oldnrow=$oldnrow, size(crspd)=$(size(crspd)))"
  #debugmsg("dbg2", crspd)
  valid = selectionop(crspd, Ftimegroup=:monthid, datefrequency=:month)
  #debugmsg("dbg2.5", valid)
  #valid = combine(groupby(dvalid[!, [:permno, :monthid, ]], [:permno, :monthid]), nrow => :ndays)
  #valid = ensurecompleteness(valid, Ftimegroup=:monthid, datefrequency=:month)
  crspd = innerjoin(crspd, valid, on=[:permno, :monthid], validate=(false,true))
  #debugmsg("dbg3", crspd)
  if !issorted(crspd, [:permno, :date])
    sort!(crspd, [:permno, :date])
  end

  #need to adjust for industry for certain versions of momentum
  crspd.plret = adjustindustry(crspd)
  #crspd.plret = Vector{Float64}(crspd.plret)
  crspd.lplret = (log).(crspd.plret .+ 1)

  print("Time to compute monthly momentum: ")
  @time momentumstrategies!(crspd, Fconditional=:monthend)
#  throw("Stop-inspect momentum")
  select!(crspd, Not(:monthend))

  if !issorted(crspd, [:permno, :date])
    sort!(crspd, [:permno, :date])
  end

  #scrspd::SubDataFrame = view(crspd, completecases(crspd, PARAM[:monthlycrsprequired]), :)
  #@assert size(crspd,1) .* 0.5 .≤ size(scrspd,1) "
  #  monthlycrsprequired dropped more than 50% of rows!!!"
  #debugmsg("dbg4", crspd)
  #combine by month
  Fmomentums = keys(PARAM[:momentumstrategies]) |> collect
  s2crspds::GroupedDataFrame = groupby(crspd, [:permno, :monthid])
  crspm::DataFrame = combine(s2crspds,
    #:permno => length=> :ndays,
    :ndays => last => :ndays,
    :date => last=> :date,
    :sigma=>last=> :sigma,
    :shares=>last=> :shares,

    :dvol => sum=> :dvol,
    :price => last=> :price,
    :tridx => last=> :tridx,

    :mc => last => :mc,
    :retP1 => prod=> :ret, #WARNING- to be adjusted later

    :ivol => last=> :ivol,
    :exchcd => last => :exchcd,
    :siccd => last => :siccd,

    :lplret => sum => :plret, #WARNING- to be adjusted later
    Fmomentums .=> last .=> Fmomentums,
    #  :ret => std=> :std,
    #yzstd = [:hi, :lo, :open, :price, :Ldistadjclose]=>yzstd
    )

  crspm.ret .-= 1.0
  crspm.plret .= exp.(crspm.plret) .- 1.0
  #debugmsg("dbg5", crspm)
  #sanity check
  testcrspms = groupby(crspm, [:permno, :monthid])
  Threads.@threads for k ∈ keys(s2crspds)
    testscrspm = testcrspms[(permno=k.permno,monthid=k.monthid)]
    @assert sum(log.(1 .+ s2crspds[k].ret)) |> exp ≈ testscrspm.ret[1] + 1.0
  end
  #conform the dates
  scrspms::GroupedDataFrame = groupby(crspm, :monthid)
  Threads.@threads for i ∈ 1:length(scrspms)
    scrspm::SubDataFrame = scrspms[i]
    @assert (maximum(scrspm.date) - minimum(scrspm.date)).value < 31
    scrspm.date .= maximum(scrspm.date)
  end

  #debugmsg("dbg6", crspm)
  return crspm
end

#recommend calling this function whenever we want only the ends of the month, and deleting afterwar
function computemonthend!(crspd::DataFrame; createnew::Bool = true)
  #compute change in month
  (:monthend ∉ names(crspd)) || (error(":monthend already in crspd.
    You probably should have deleted it the last time you used it."))
  crspd.monthend = falses(size(crspd,1))
  @assert issorted(crspd, [:permno, :date])

  scrspds = groupby(crspd, [:permno, :monthid])
  Threads.@threads for i::Int ∈ 1:length(scrspds)
    scrspd::SubDataFrame = scrspds[i]
    scrspd[end, :monthend] = true
  end

  return nothing
end

#assigns a monthid for each date
function assignmonthid(dates::AbstractVector{Date}, mindate::Date=Date(1900,1,1))::Vector{Int}
  day(mindate) == 1 || (
    @warn "day of $mindate ≠ 1- mindate must be first day of month")

  #first compute all possible first days of the month
  maxdate = maximum(dates)
  allmonths = mindate:Month(1):firstdayofmonth(maxdate) |> collect

  #create an index to look up the monthid
  monthidindex = Dict(allmonths[mid] => mid for mid ∈ 1:length(allmonths))

  monthid::Vector{Int} = (d::Date->monthidindex[d |> firstdayofmonth]).(dates)

  return monthid
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function load_parameters(yaml_file::String)
  d = YAML.load_file(yaml_file)

  check_yml_params(d)
  check_param_shapes(d)
  reformat_params!(d)

  outd = Dict(:species => d["species"], :model => d["model"])

  return outd
end


text:
struct Twist{T} <: StaticVector{6,T}
    angular::AngularVelocity{T}
    linear::SVector{3,T}
end
function (::Type{T})(xs::NTuple{6}) where {T<:Twist}
    (x1, x2, x3, x4, x5, x6) = Base.promote(xs...)
    ω = AngularVelocity{typeof(x1)}(x1, x2, x3)
    v = SVector{3,typeof(x4)}(x4, x5, x6)
    return T(ω, v)
end

Base.@propagate_inbounds function Base.getindex(x::Twist, i::Int)
    return i < 4 ? getindex(x.angular, i) : getindex(x.linear, i - 3)
end

@inline Base.Tuple(t::Twist) = (Tuple(t.angular)..., Tuple(t.linear)...)

function groupexp(t::Twist)
    ω, b = t.angular, t.linear
    θ = rotation_angle(ω)
    rotation = groupexp(t.angular)
    translation = if iszero(θ)
        b
    else
        sinθ, cosθ = sincos(θ)
        θ² = θ^2
        p₁, p₂ = (1 - cosθ) / θ², (1 - sinθ / θ) / θ²
        ωxb = cross(ω, b)
        b + p₁ * ωxb + p₂ * cross(ω, ωxb)
    end
    return AffineMap(rotation, translation)
end

function grouplog(x::RigidMotion{3})
    r, t = x.linear, x.translation
    angular = grouplog(r)
    θ = rotation_angle(angular)
    linear = if iszero(θ)
        t
    else
        sinθ, cosθ = sincos(θ)
        θ² = θ^2
        p₂ = 1 / θ² - (1 + cosθ) / 2 / θ / sinθ
        ωxt = cross(angular, t)
        t - ωxt / 2 + p₂ * cross(angular, ωxt)
    end
    return Twist(angular, linear)
end


text:
"""
    SimulationLogs

Log internal variables in a DifferentialEquations.jl simulation with the `@log` macro and
use the `get_log` function to retrieve the logged variables.
"""
module SimulationLogs

using DiffEqBase: CallbackSet, DiscreteCallback, ContinuousCallback, isinplace, solve
using RecipesBase: RecipesBase, @recipe, @series, @userplot
using SciMLBase: AbstractTimeseriesSolution, remake

include("types.jl")
include("macros.jl")
include("logging.jl")
include("plotting.jl")

const GLOBAL_LOG = SimulationLog()

export SimulationLog, is_active, value_dict
export Logged
export @log
export get_log, logged_solve, scope, scope!

end # module

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
#####
##### Fields computed from abstract operations
#####

using KernelAbstractions: @kernel, @index
using Oceananigans.Grids: default_indices
using Oceananigans.Fields: FieldStatus, reduced_dimensions, validate_indices
using Oceananigans.Utils: launch!

import Oceananigans.Fields: Field, compute!

const ComputedField = Field{<:Any, <:Any, <:Any, <:AbstractOperation}

"""
    Field(operand::AbstractOperation; kwargs...)

Return `f::Field` where `f.data` is computed from `f.operand` by
calling compute!(f).

Keyword arguments
=================

data (AbstractArray): An offset Array or CuArray for storing the result of a computation.
                      Must have `total_size(location(operand), grid)`.

boundary_conditions (FieldBoundaryConditions): Boundary conditions for `f`.

recompute_safely (Bool): whether or not to _always_ "recompute" `f` if `f` is
                         nested within another computation via an `AbstractOperation`.
                         If `data` is not provided then `recompute_safely=false` and
                         recomputation is _avoided_. If `data` is provided, then
                         `recompute_safely=true` by default.
"""
function Field(operand::AbstractOperation;
               data = nothing,
               indices = default_indices(3),
               boundary_conditions = FieldBoundaryConditions(operand.grid, location(operand)),
               recompute_safely = true)

    grid = operand.grid
    loc = location(operand)
    indices = validate_indices(indices, loc, grid)

    if isnothing(data)
        data = new_data(grid, loc, indices)
        recompute_safely = false
    end

    status = recompute_safely ? nothing : FieldStatus()

    return Field(loc, grid, data, boundary_conditions, indices, operand, status)
end

"""
    compute!(comp::ComputedField)

Compute `comp.operand` and store the result in `comp.data`.
"""
function compute!(comp::ComputedField, time=nothing)
    # First compute `dependencies`:
    compute_at!(comp.operand, time)

    arch = architecture(comp)
    event = launch!(arch, comp.grid, size(comp), _compute!, comp.data, comp.operand, comp.indices)
    wait(device(arch), event)

    fill_halo_regions!(comp)

    return comp
end

@inline offset_compute_index(::Colon, i) = i
@inline offset_compute_index(range::UnitRange, i) = range[1] + i - 1

"""Compute an `operand` and store in `data`."""
@kernel function _compute!(data, operand, index_ranges)
    i, j, k = @index(Global, NTuple)

    i′ = offset_compute_index(index_ranges[1], i)
    j′ = offset_compute_index(index_ranges[2], j)
    k′ = offset_compute_index(index_ranges[3], k)

    @inbounds data[i′, j′, k′] = operand[i′, j′, k′]
end


text:
module LightTable

using JSON, Lazy, Jewel, LNR, Requires

import Jewel: lines

export server, ltprint, popup, notify

exit_on_sigint(on) = ccall(:jl_exit_on_sigint, Void, (Cint,), on)

const connect_notify = Condition()

function flushpipe(io)
  @schedule begin
    wait(connect_notify)
    while true
      ltprint(readline(io))
    end
  end
end

# -------------------
# Basic Communication
# -------------------

function server(port, id, headless = false)
  global isheadless = headless
  exit_on_sigint(false)
  ltconnect(port, id)
  headless && pushdisplay(LTConsole())
  while isopen(conn)
    try
      handlenext()
    catch e
      warn("LightTable.jl: "sprint(showerror, e, catch_backtrace()))
    end
  end
end

function ltconnect(port, id)
  global conn = connect(port)
  ltwrite(@d("type" => "julia",
             "name" => "Julia",
             "commands" => ["editor.eval.julia", "editor.julia.hints", "editor.julia.doc"],
             "client-id" => id))
end

function ltclose()
  close(conn)
end

function ltwrite(data)
  @assert isopen(conn)
  println(conn, json(data))
end

function ltread()
  @assert isopen(conn)
  data = JSON.parse(conn)
  eof(conn) || read(conn, Char) # newline
  return data
end

raise(object::Integer, event, data = nothing) = ltwrite([object, event, data])

# ----------------
# Command Handling
# ----------------

const cmdhandlers = Dict{AbstractString, Function}()

function handlecmd(data)
  data == nothing && return
  cmd = data[2]
  if haskey(cmdhandlers, cmd)
    cmdhandlers[cmd](data[1], data[3])
  else
    warn("Can't handle command $cmd")
  end
end

handle(f, cmd) = (cmdhandlers[cmd] = f)

const cmdqueue = c()

function queuecmds()
  while nb_available(conn) > 0
    push!(cmdqueue, ltread())
  end
end

handlenext() = handlecmd(!isempty(cmdqueue) ? shift!(cmdqueue) : ltread())

handle("notify-connected") do client, data
  raise(client, "connected")
end

handle("client.close") do req, data
  close(conn)
end

# ------
# Others
# ------

include("commands.jl")
include("interaction/interaction.jl")
include("interaction/collapsible.jl")
include("interaction/scales.jl")
include("interaction/lazy.jl")
include("eval.jl")
include("completions.jl")
include("doc.jl")
include("errorshow.jl")
include("display/display.jl")
include("misc.jl")

end


text:
module LoRDeckCodes

export Deck, CardCodeAndCount, DeckEncoder

include("types.jl")
include("VarintTranslator.jl")
include("Base32.jl")
include("DeckEncoder.jl")

end # module LoRDeckCodes


text:
using Printf
using TimerOutputs
using Oceananigans
using Oceananigans.Solvers
using Oceananigans.Utils

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Ni = 2   # Number of iterations before benchmarking starts.
Nt = 10  # Number of iterations to use for benchmarking time stepping.

# Model resolutions to benchmarks. Focusing on 3D models for GPU benchmarking.
            Ns = [(256, 256, 256)]
   float_types = [Float32, Float64]  # Float types to benchmark.
         archs = [CPU()]             # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]      # Benchmark GPU on systems with CUDA-enabled GPUs.

#####
##### Run benchmarks
#####

for arch in archs, float_type in float_types, N in Ns

    #####
    ##### Horizontally periodic pressure solve on CPU
    #####

    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), extent=(1.0, 2.5, 3.6))
    solver = PoissonSolver(arch, PPN(), grid)
    fbcs = HorizontallyPeriodicBCs()

    RHS = CellField(FT, arch, grid)
    interior(RHS) .= rand(Nx, Ny, Nz)
    interior(RHS) .= interior(RHS) .- mean(interior(RHS))

    RHS_orig = deepcopy(RHS)
    solver.storage .= interior(RHS)
    solve_poisson_3d!(solver, grid)

    Nx, Ny, Nz = N
    Lx, Ly, Lz = 1, 1, 1

    model = Model(architecture=arch, float_type=float_type, grid=RegularCartesianGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz)))
    time_step!(model, Ni, 1)

    bname =  benchmark_name(N, "", arch, float_type)
    @printf("Running static ocean benchmark: %s...\n", bname)
    for i in 1:Nt
        @timeit timer bname time_step!(model, 1, 1)
    end
end


text:
# specialized methods
"""
```julia
dotdot(::Vec, ::SymmetricFourthOrderTensor, ::Vec)
```
Computes a special dot product between two vectors and a symmetric fourth order tensor
such that ``a_k C_{ikjl} b_l``.
"""
@generated function dotdot{dim}(v1::Vec{dim}, S::SymmetricTensor{4, dim}, v2::Vec{dim})
    idx(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    exps = Expr(:tuple)
    for j in 1:dim, i in 1:dim
        ex1, ex2 = Expr[],Expr[]
        for l in 1:dim, k in 1:dim
            push!(ex1, :(get_data(v1)[$k] * get_data(S)[$(idx(i,k,j,l))]))
            push!(ex2, :(get_data(v2)[$l]))
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    return quote
        $(Expr(:meta, :inline))
        @inbounds r = $exps
        Tensor{2, dim}(r)
    end
end


text:
const plot_dict = OrderedDict(
    "lineplot" => plot,
    "scatter" => scatter,
    "bar" => groupedbar,
    "boxplot" => boxplot,
    "violin" => violin,
    "histogram2d" => histogram2d,
    "marginalhist" => marginalhist,
)

function dropdownrow(t::IndexedTable)
    ns = collect(colnames(t))
    x = dropdown(ns, label = "x")
    y = dropdown(vcat([:density, :cumulative, :hazard], ns), label = "y")
    plot_type = dropdown(plot_dict, label = "plot")
    axis_type = dropdown(Symbol.(["auto", "pointbypoint", "discrete", "binned", "continuous"]), label = "axis_type")
    across_strs = vcat(["none", "bootstrap", "across"], "across " .* string.(ns))
    across_vals = map(across_map, across_strs)
    across = dropdown(OrderedDict(zip(across_strs, across_vals)), label = "compute_error")
    package = dropdown(
        OrderedDict("auto" => nothing, "StatPlots" => StatPlotsRecipe, "GroupedErrors" => GroupedError),
        label = "module"
    )
    [x, y, plot_type, axis_type, across, package]
end

function across_map(s)
    if s=="none"
        nothing
    elseif s == "bootstrap"
        (:bootstrap, 1000)
    elseif s == "across"
        (:across, :all)
    else
        (:across, Symbol(match(r"^across (.*)$", s)[1]))
    end
end

selecteditems(c::AbstractObservable) = c[]
selecteditems(c::AbstractArray) = selecteditems.(c)


text:
abstract type AbstractFileType end

struct GMSHfile <: AbstractFileType
  file_name::String
end

# Define constructor for Mesh that takes in GMSH file
function Mesh(file_obj::GMSHfile)
  xy,cm,bndNodes,bndNames,order = GMSH_QUAD_READER(file_obj.file_name)

  bndNodes2 = GMSH_EXTRACTBOUNDARY_ALL(bndNodes,length(bndNames))
  bndNames2 = Symbol.(bndNames)
  bndNames3 = Index(bndNames2)
  bdry = DictList{Int}(bndNodes2,bndNames3)

  # generate element list
  cmEl = [Element(cm[i,:]) for i in 1:size(cm,1)]

  # generate node data structure from xy,xyp
  xyNodes = [Node(xy[i,1],xy[i,2]) for i=1:size(xy,1)]

  # turn order into a symbol (fast comparisons)
  order==1 ? o = Symbol("Linear") : o = Symbol("Quadratic")

  return ScalarMesh(xyNodes,cmEl,bdry,o)
end

# Define constructor for Fluid Mesh that takes in GMSH file
function FluidMesh(file_obj::GMSHfile)
  # generate normal mesh routine -- assumes that order of elements is 2
  xy,cm,bndNodes,bndNames,order = GMSH_QUAD_READER(file_obj.file_name)
  xyp,cmp = pMesh(xy,cm)

  bndNodes2 = GMSH_EXTRACTBOUNDARY_ALL(bndNodes,length(bndNames))
  bndNames2 = Symbol.(bndNames)
  bndNames3 = Index(bndNames2)
  bdry = DictList{Int}(bndNodes2,bndNames3)

  # generate element list
  cmEl  = [Element(cm[i,:])  for i in 1:size(cm,1)]
  cmpEl = [Element(cmp[i,:]) for i in 1:size(cmp,1)]

  # generate node data structure from xy,xyp
  xyNodes  = [Node(xy[i,1],xy[i,2])   for i=1:size(xy,1)]
  xypNodes = [Node(xyp[i,1],xyp[i,2]) for i=1:size(xyp,1)]

  # turn order into a symbol (fast comparisons)
  order==1 ? o = Symbol("Linear") : o = Symbol("Quadratic")

  return FluidMesh(xyNodes,cmEl,bdry,o,xypNodes,cmpEl)
end

"""
  printBlankLine(f)

prints blank line to IOstream handle f
"""
function printBlankLine(f)
  println(f,"")
end

"""
  GMSH_QUAD_READER(file_name) -> xy,cm,bndNodes

reading in GMSH file that assumes 2D quadrilateral elements.

To determine boundaries, you need to look at *.msh file header to see identifying numbers of boundaries, and afterwards group them how you want.
"""
function GMSH_QUAD_READER(file_name::String)
  # ---
  # open file
  # ---

  f          = open(file_name,"r")
	lines      = readlines(f)           # can be large ... might want to find
                                      # alternative method extract file info
	lineLength = length(lines)

	# ---
  # initialize tracking values
  # ---

  surfaceObjID = -1
  lineObjID    = -1
  npe          = -1

  numNodes = 0
  nodesLineStart = 0
  numElements = 0
  elementsLineStart = 0
  numBoundaries = 0
  numLineObjs = 0

  order = 0
  bndNames::Vector{String} = []

  # determine order
  for lineNumber = 1:lineLength
    ln = lines[lineNumber]
    if isElements(ln)
      elmStr = split(lines[lineNumber+2]," ")
      elType = str2int(elmStr[2])
      elType == 1 ? order=1 : order=2
    end
  end

  # ---
  # correct numbering expected for order
  # ---

  if order==1
    surfaceObjID = 3
    lineObjID = 1
    npe = 4
  elseif order==2
    surfaceObjID = 10
    lineObjID = 8
    npe = 9
  else
    #error("invalid order entered"::String)
  end

  # ---
  # First pass through file to obtain parameters
	# ---
	for lineNumber = 1:lineLength
	  ln = lines[lineNumber]
	  # boundary data
    if isBoundary(ln)
      numBoundaries = parse(Int,lines[lineNumber+1])

      for i=(lineNumber+2):(lineNumber+2+numBoundaries-2)
        str = lines[i]
        name = str[6:end-1]
        push!(bndNames,name)
      end
    end
    # node data
    if isNodes(ln)
      nodesLineStart = lineNumber+2
      numNodes = parse(Int,lines[lineNumber+1])
    end
    # element data
    if isElements(ln)
      elementsLineStart = lineNumber+2
      numElements = parse(Int,lines[lineNumber+1])
      for elm = 1:numElements
        elmStr = split(lines[lineNumber+elm+1]," ")
        elType = str2int(elmStr[2])
        if elType == lineObjID
          numLineObjs += 1
        end
      end
    end
  end

  # check that 'interior' surface exists
  ln = lines[5 + numBoundaries]
  name = ln[6:end-1]
  if !(isInterior(name))
    error("invalid mesh: requires PhysicalBoundary surface 'interior'")
  end

  # adjust element indexing to remove line objects
  numElements       -= numLineObjs
  lineObjsStart      = elementsLineStart
  elementsLineStart += numLineObjs

  # ---
  # generate arrays using above parameters
  # ---

  xy = zeros(Float64,(numNodes,2))
  cm = zeros(Int,(numElements,npe))
  # needs check in case numNodes < numElements, which happens for small meshes
  if numNodes < (numElements + numLineObjs)
    bndNodes = zeros(Int,(numNodes+numElements+numLineObjs,numBoundaries))
  else
    bndNodes = zeros(Int,(numNodes,numBoundaries))
  end

  # ---
  # obtain node data
  # ---
  nodeCounter = 1
  for lineNumber = nodesLineStart:(nodesLineStart+numNodes-1)
    ln = lines[lineNumber]
    POINTS = split(ln," ")
    xy[nodeCounter,1] = str2float(POINTS[2])
    xy[nodeCounter,2] = str2float(POINTS[3])
    nodeCounter += 1
  end

  # ---
  # obtain boundary data
  # ---
  bndCounter = 1
  for lineNumber = lineObjsStart:(lineObjsStart+numLineObjs-1)
    ln = lines[lineNumber]
    elmStr = split(ln," ")

    skip = parse(Int,elmStr[3]) + 3

    # boundary array
    physicalEntity = parse(Int,elmStr[4])
    bndNodes[bndCounter  ,physicalEntity] = str2int(elmStr[skip+1])
    bndNodes[bndCounter+1,physicalEntity] = str2int(elmStr[skip+2])
    if lineObjID==8   # quadratic elements have 3 line nodes
      bndNodes[bndCounter+2,physicalEntity] = str2int(elmStr[skip+3])
      bndCounter += 1
    end
    bndCounter += 2
  end

  # ---
  # obtain element connectivity data
  # ---
  elmCounter = 1

  for lineNumber = elementsLineStart:(elementsLineStart+numElements-1)
    ln = lines[lineNumber]
    elmStr = split(ln," ")

    skip = parse(Int,elmStr[3]) + 3

    for k=1:npe
      cm[elmCounter,k] = str2float(elmStr[skip+k])
    end
    elmCounter += 1
  end

  return xy,cm,bndNodes,bndNames,order
end

"""
  isBoundary(str)

Checks whether string is '\$PhysicalNames' function to be used when reading in *.msh files from gmsh. Used in GSMH_QUAD_READER
"""
function isBoundary(str)
    if VERSION==v"0.5.0"
      str = str[1:end-1]
    end
    if str == "\$PhysicalNames"
        return true
    else
        return false;
    end
end

"""
  isElements(str)

Checks whether string is '\$Elements' function to be used when reading in *.msh files from gmsh. Used in GSMH_QUAD_READER
"""
function isElements(str)
    if VERSION==v"0.5.0"
      str = str[1:end-1]
    end
    if str == "\$Elements"
        return true
    else
        return false;
    end
end

"""
  isNodes(str)

Checks whether string is '\$Nodes' function to be used when reading in
*.msh files from gmsh. Used in GSMH_QUAD_READER
"""
function isNodes(str)
    if VERSION==v"0.5.0"
      str = str[1:end-1]
    end
    if str == "\$Nodes"
        return true
    else
        return false;
    end
end

function isInterior(str)
  if VERSION==v"0.5.0"
    str = str[1:end-1]
  end
  if str == "interior"
      return true
  else
      return false;
  end
end

"""
  str2int(str)

parses String into Int
"""
function str2int(str)
  return parse(Int,str)
end

"""
  str2float(str)

parses String into Float64
"""
function str2float(str)
  return parse(Float64,str)
end

"""
  GMSH_BOUNDARYPARSE(N)

Takes in gmsh boundary matrix and returns first N boundaries (sorted and remove 0 element).

Really crappy programming, just a stopgap until a more elegant solution is found
"""
function GMSH_EXTRACTBOUNDARY(bndNodes,i)
  # weird aspect of GSMH reader implementation includes 0 at beginning of sorted array. Need to include check in case implementation is changed later on
  return sort(unique(bndNodes[:,i]))[2:end]
end

"""
  GMSH_EXTRACTBOUNDARY_ALL(bndNodes,total)

completely stupid way to extract boundary node lists, simply used to shorten code length in driver files
"""
function GMSH_EXTRACTBOUNDARY_ALL(bndNodes,total)
  # weird aspect of GSMH reader implementation includes 0 at beginning of sorted array. Need to include check in case implementation is changed later on
  if total==1
    A = sort(unique(bndNodes[:,1]))[2:end]
    return [A]
  elseif total==2
    A = sort(unique(bndNodes[:,1]))[2:end]
    B = sort(unique(bndNodes[:,2]))[2:end]
    return [A,B]
  elseif total==3
    A = sort(unique(bndNodes[:,1]))[2:end]
    B = sort(unique(bndNodes[:,2]))[2:end]
    C = sort(unique(bndNodes[:,3]))[2:end]
    return [A,B,C]
  elseif total==4
    A = sort(unique(bndNodes[:,1]))[2:end]
    B = sort(unique(bndNodes[:,2]))[2:end]
    C = sort(unique(bndNodes[:,3]))[2:end]
    D = sort(unique(bndNodes[:,4]))[2:end]
    return [A,B,C,D]
  elseif total==5
    A = sort(unique(bndNodes[:,1]))[2:end]
    B = sort(unique(bndNodes[:,2]))[2:end]
    C = sort(unique(bndNodes[:,3]))[2:end]
    D = sort(unique(bndNodes[:,4]))[2:end]
    E = sort(unique(bndNodes[:,5]))[2:end]
    return [A,B,C,D,E]
  elseif total==6
    A = sort(unique(bndNodes[:,1]))[2:end]
    B = sort(unique(bndNodes[:,2]))[2:end]
    C = sort(unique(bndNodes[:,3]))[2:end]
    D = sort(unique(bndNodes[:,4]))[2:end]
    E = sort(unique(bndNodes[:,5]))[2:end]
    F = sort(unique(bndNodes[:,6]))[2:end]
    return [A,B,C,D,E,F]
  else
    error("total=$(total) is too large for this file")
  end
end


text:
import Pkg
Pkg.activate(".")
Pkg.instantiate()

using Flutes
using Dates

# input env vars
tuning = parse(Float64, readvariable("FLUTE_TUNING"))
scale = tones("FLUTE_SCALE", tuning)
maxd = floats("FLUTE_MAX_DIAMETERS")
maxp = floats("FLUTE_MAX_PADDING")
rots = floats("FLUTE_HOLE_ROTATIONS")
brk = parse(Int, readvariable("FLUTE_BREAK"))
head_length = parse(Float64, readvariable("FLUTE_HEAD_LENGTH"))
tenon_length = parse(Float64, readvariable("FLUTE_TENON_LENGTH"))
trace = parse(Bool, readvariable("TRACE"))

function rround(value)
  round(value; digits=3)
end

# minimum diameters = 1mm
mind = fill(1, length(maxd))
# calculate minimum padding between holes
minp = []
for h in 1:length(maxd)
  mp = 0
  if h == 1
    # first hole must be on body, not the headjoint
    mp = head_length+maxd[h]
  elseif h == brk+1
    # make room for the tenon/mortise at the break point
    mp = maxd[brk]+tenon_length+maxd[h]
  else
    # minimum hole spacing (sum of max diameters)
    mp = maxd[h-1]+maxd[h]
  end
  append!(minp, mp)
end

# find best fit
# all the magic happens here
println("Optimizing flute parameters: ", ARGS[1])
diameters, error = optimal(scale, mind, maxd, minp, maxp; trace=trace)
lengths = mapflute(scale, diameters)
# end magic

# display hole separations
with_origin = [0.0; lengths]
print("|")
for h in 1:length(lengths)
  print(" ∘ ", round(with_origin[h+1]-with_origin[h]; digits=1))
end
println(" |")

# break holes by foot/body
body_diameters = diameters[1:brk]
body_positions = lengths[1:brk]
body_rotations = rots[1:brk]

foot_diameters = diameters[brk+1:end]
foot_positions = lengths[brk+1:end-1]
foot_rotations = rots[brk+1:end]

# place body/foot joint
spare = max((foot_positions[1] - body_positions[end] - tenon_length)/2, 0)
nofoot = body_positions[end] + spare + tenon_length

flute_length = lengths[end]
body_length = nofoot - head_length
foot_length = flute_length - nofoot

# export parameters to opencad props
params = createscadparameters()
headset = "head.data"
setscadparameter!(params, headset, "HeadLength", head_length)
setscadparameter!(params, headset, "TenonLength", tenon_length)
bodyset = "body.data"
setscadparameter!(params, bodyset, "BodyLength", rround(body_length))
setscadparameter!(params, bodyset, "HoleDiameters", rround.(body_diameters))
setscadparameter!(params, bodyset, "HolePositions", map(bp->rround(bp-head_length), body_positions))
setscadparameter!(params, bodyset, "HoleRotations", body_rotations)
setscadparameter!(params, bodyset, "TenonLength", tenon_length)
setscadparameter!(params, bodyset, "MortiseLength", tenon_length)
footset = "foot.data"
setscadparameter!(params, footset, "FootLength", rround(foot_length))
setscadparameter!(params, footset, "HoleDiameters", rround.(foot_diameters))
setscadparameter!(params, footset, "HolePositions", map(fp->rround(fp-nofoot), foot_positions))
setscadparameter!(params, footset, "HoleRotations", foot_rotations)
setscadparameter!(params, footset, "MortiseLength", tenon_length)
extraset = "extra.data"
setscadparameter!(params, extraset, "CreationDate", now())
setscadparameter!(params, extraset, "Tuning", tuning)
setscadparameter!(params, extraset, "Scale", readvariable("FLUTE_SCALE"))
setscadparameter!(params, extraset, "Score", rround(error))
writescadparameters(params, ARGS[1])

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
module TestExamples3DStructured

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "structured_3d_dgsem")

@testset "Structured mesh" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [0.00013446460962856976],
      linf = [0.0012577781391462928])
  end

  @trixi_testset "elixir_advection_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_free_stream.jl"),
      l2   = [1.830875777528287e-14],
      linf = [7.491784970170556e-13],
      atol = 8e-13, # required to make tests pass on Windows
      )
  end

  @trixi_testset "elixir_advection_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"),
      l2   = [6.522004549411137e-5],
      linf = [0.005554857853361295])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.0281388160824776],
      linf = [0.08740635193023694])
  end

  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [0.01032310150257373, 0.009728768969448439, 0.009728768969448494, 0.009728768969448388, 0.015080412597559597],
      linf = [0.034894790428615874, 0.033835365548322116, 0.033835365548322116, 0.03383536554832034, 0.06785765131417065])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [2.8815700334367128e-15, 9.361915278236651e-15, 9.95614203619935e-15, 1.6809941842374106e-14, 1.4815037041566735e-14],
      linf = [4.1300296516055823e-14, 2.0444756998472258e-13, 1.0133560657266116e-13, 2.0627943797535409e-13, 2.8954616482224083e-13])
  end

  @trixi_testset "elixir_euler_free_stream.jl with FluxRotated(flux_lax_friedrichs)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      surface_flux=FluxRotated(flux_lax_friedrichs),
      l2   = [2.8815700334367128e-15, 9.361915278236651e-15, 9.95614203619935e-15, 1.6809941842374106e-14, 1.4815037041566735e-14],
      linf = [4.1300296516055823e-14, 2.0444756998472258e-13, 1.0133560657266116e-13, 2.0627943797535409e-13, 2.8954616482224083e-13])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [0.0018268326813744103, 0.001745601029521995, 0.001745601029521962, 0.0017456010295218891, 0.003239834454817457],
      linf = [0.014660503198892005, 0.01506958815284798, 0.01506958815283821, 0.015069588152864632, 0.02700205515651044])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.011367083018614027, 0.007022020327490176, 0.006759580335962235, 0.006820337637760632, 0.02912659127566544],
      linf = [0.2761764220925329, 0.20286331858055706, 0.18763944865434593, 0.19313636558790004, 0.707563913727584],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.009082353036644902, 0.007128360240528109, 0.006970330025996491, 0.006898850266874514,
              0.03302008823756457, 0.003203389099143526, 0.003077498677885352, 0.0030740006760477624,
              4.192129696970217e-5],
      linf = [0.2883946030582689, 0.25956437344015054, 0.2614364943543665, 0.24617277938134657,
              1.1370443512475847, 0.1278041831463388, 0.13347391885068594, 0.1457563463643099,
              0.0021174246048172563],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003015476175153681, 0.00145499403283373, 0.0009125744757935803, 0.0017703080480578979,
              0.0013046447673965966, 0.0014564863387645508, 0.0013332311430907598, 0.001647832598455728,
              0.0013647609788548722],
      linf = [0.027510637768610846, 0.02797062834945721, 0.01274249949295704, 0.038940694415543736,
              0.02200825678588325, 0.03167600959583505, 0.021420957993862344, 0.03386589835999665,
              0.01888303191983353])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_lax_friedrichs" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003096205134732006, 0.0014592894240925266, 0.000907526173067474, 0.001821970949880564,
              0.001289067735865978, 0.0014595578872467106, 0.0013269375784248857, 0.0016998139907659089,
              0.0013685168042016343],
      linf = [0.027618475633440998, 0.027093787212065318, 0.012584560784257667, 0.039456640084648914,
              0.020759073985165077, 0.031771018340953416, 0.02059036404759229, 0.03456102393654076,
              0.019663511833857894],
      surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell))
  end

  # TODO: nonconservative terms, remove
  @trixi_testset "elixir_mhd_alfven_wave.jl with old nonconservative stuff" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003015476175153681, 0.00145499403283373, 0.0009125744757935803, 0.0017703080480578979,
              0.0013046447673965966, 0.0014564863387645508, 0.0013332311430907598, 0.001647832598455728,
              0.0013647609788548722],
      linf = [0.027510637768610846, 0.02797062834945721, 0.01274249949295704, 0.038940694415543736,
              0.02200825678588325, 0.03167600959583505, 0.021420957993862344, 0.03386589835999665,
              0.01888303191983353],
      volume_flux = flux_central,
      surface_flux = flux_hll)
  end
end

end # module


text:
using CartesianGeneticProgramming
using ArgParse
using ArcadeLearningEnvironment
include("../graphing/graph_utils.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/atari_ram.yaml"
    "--game"
    help = "game rom name"
    default = "pong"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)

ale = ALE_new()
loadROM(ale, args["game"])
n_in = length(getRAM(ale))
n_out = length(getMinimalActionSet(ale))
ALE_del(ale)

cfg = get_config(args["cfg"]; game=args["game"], n_in=n_in, n_out=n_out)

ind = CGPInd(cfg, read(args["ind"], String))
dot_array = walk_nodes(ind)
dot_file  = "graph.dot"
open(dot_file, "w") do f
   for i in dot_array
      println(f, i)
   end
end
# to view from commandline:
# $ dot -Tpng graph.dot -o graph.png
# $ display graph.png


text:
@resumable function recruitProcess( sim::Simulation, recruitment::Recruitment,
    mpSim::MPsim )

    # Timing of the process.
    processTime = Dates.Millisecond( 0 )
    tStart = now()
    recName = string( "Recruitment process '", recruitment.name, "' to node '",
        recruitment.targetNode, "' " )

    # Preparatory steps.
    priority = recruitment.priority -
        ( recruitment.isAdaptive ? 1 : 0 ) * mpSim.nPriorities
    timeToWait = recruitment.offset

    # If an initial population snapshot is uploaded, and it contains zero-time
    #   events, don't execute zero-time recruitment.
    if ( timeToWait == 0 ) && !mpSim.isVirgin
        timeToWait += recruitment.freq
    end  # if ( timeToWait == 0 ) && ...

    # Process loop.
    while now( sim ) + timeToWait <= mpSim.simLength
        processTime += now() - tStart
        @yield timeout( sim, timeToWait, priority = priority )
        tStart = now()
        recruitmentCycle( recruitment, mpSim )
        timeToWait = recruitment.freq
    end  # while now( sim ) + timeToWait <= mpSim.simLength

    # Timing of the process.
    processTime += now() - tStart

    if mpSim.showInfo
        println( recName, "took ", processTime.value / 1000, " seconds." )
    end  # if mpSim.showInfo

end  # @resumable recruitProcess( sim, recruitment, mpSim )

function recruitmentCycle( recruitment::Recruitment, mpSim::MPsim )
    nToRecruit = generatePoolSize( recruitment,
        mpSim.baseNodeList[recruitment.targetNode], mpSim )
    generatePersons( recruitment, nToRecruit, mpSim )
end  # recruitmentCycle( recruitment, mpSim )

function generatePoolSize( recruitment::Recruitment, targetNode::BaseNode,
    mpSim::MPsim )
    # If the recruitment scheme isn't adaptive, draw from the distribution.
    !recruitment.isAdaptive && return recruitment.recruitmentDist()

    nToOrgTarget = mpSim.personnelTarget > 0 ? mpSim.personnelTarget -
        mpSim.orgSize : typemax( Int )
    nToNodeTarget = targetNode.target >= 0 ? targetNode.target -
        length( targetNode.inNodeSince ) : typemax( Int )
    return max( min( recruitment.maxRecruitment, nToOrgTarget, nToNodeTarget ),
        recruitment.minRecruitment )
end  # generatePoolSize( recruitment, targetNode, mpSim )

function generatePersons( recruitment::Recruitment, nToRecruit::Int,
    mpSim::MPsim )
    nToRecruit == 0 && return

    # Generate attribute values for the new personnel members.
    attributeNames = collect(keys(mpSim.attributeList))
    newPersons = DataFrame( fill( "", nToRecruit,
        length(mpSim.attributeList) ), Symbol.(attributeNames) )
    targetNode = mpSim.baseNodeList[recruitment.targetNode]

    for name in attributeNames
        attribute = mpSim.attributeList[name]

        if haskey( targetNode.requirements, name )
            newPersons[:, name] .= targetNode.requirements[name]
        elseif !isempty(attribute.initValues)
            newPersons[:, name] = generateValues( attribute, nToRecruit )
        end  # if !isempty(attribute.initValues)
    end  # for name in attributeNames

    # Add ID key and ages.
    ids = string.( "Sim", mpSim.dbSize .+ (1:nToRecruit) )
    insertcols!( newPersons, 1, Symbol(mpSim.idKey) => ids )
    insertcols!( newPersons, 2, :recAge => recruitment.ageDist(nToRecruit) )

    # Update target node populaton.
    setindex!.( Ref(targetNode.inNodeSince ), now(mpSim), ids )

    # Entries into databases.
    attrition = mpSim.attritionSchemes[targetNode.attrition]
    timeOfAttrition = now(mpSim) .+ generateTimeToAttrition( attrition,
        nToRecruit )
    insertcols!( newPersons, 3, :attrTime => timeOfAttrition )

    # Personnel database.
    attrCmd = ""

    if !isempty(attributeNames)
        attrCmd = map( attributeNames ) do name
            return string.( ", '", newPersons[:, Symbol( name )], "'" )
        end  # map( attributeNames ) do name

        attrCmd = hcat( attrCmd... )
        attrCmd = map( ii -> join(attrCmd[ii, :]), 1:nToRecruit )
    end  # if !isempty(attributeName)

    sqliteCmd = map( toa -> toa == +Inf ? "NULL" : toa, timeOfAttrition )
    # sqliteCmd = string.( "\n    ('", newPersons[:, Symbol(mpSim.idKey)],
    sqliteCmd = string.( "\n    ('", newPersons[:, mpSim.idKey],
        "', 'active', ", now(mpSim), ", ", newPersons[:, :recAge], ", ",
        sqliteCmd, ", '", targetNode.name, "', ", now(mpSim), attrCmd,
         ")" )
    sqliteCmd = string( "INSERT INTO `", mpSim.persDBname, "` (`", mpSim.idKey,
        "`, status, timeEntered, ageAtRecruitment, ",
        "expectedAttritionTime, currentNode, inNodeSince",
        join( string.( ", `", attributeNames, "`" ) ), ") VALUES",
        join( sqliteCmd, "," ) )
    # println(sqliteCmd)
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    # Transition database.
    sqliteCmd = string.( "\n    ('", ids, "', ", now(mpSim), ", '",
        recruitment.name, "', '", recruitment.targetNode, "')" )
    sqliteCmd = string( "INSERT INTO `", mpSim.transDBname, "` (`", mpSim.idKey,
        "`, timeIndex, transition, targetNode) VALUES", join( sqliteCmd, "," ) )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    # History database.
    if !isempty(attributeNames)
        sqliteCmd = map(attributeNames) do name
            return string.( "\n    ('", ids, "', ", now(mpSim), ", '", name,
                "', '", newPersons[:, Symbol(name)], "')" )
        end  # map(attributeNames) do name

        sqliteCmd = string( "INSERT INTO `", mpSim.histDBname, "` (`",
            mpSim.idKey, "`, timeIndex, attribute, value) VALUES",
            join( vcat(sqliteCmd...), "," ) )
        DBInterface.execute( mpSim.simDB, sqliteCmd )
    end  # if !isempty(attributeNames)

    mpSim.orgSize += nToRecruit
    mpSim.dbSize += nToRecruit
end  # generatePersons( recruitment, nToRecruit, mpSim )


text:
export WitnessSet, witness_set, linear_subspace, system, dim, codim, trace_test

"""
    WitnessSet(F, L, S)

Store solutions `S` of the polynomial system `F(x) = L(x) = 0` into a witness set.
"""
struct WitnessSet{
    S<:AbstractSystem,
    Sub<:AbstractSubspace,
    R<:Union{Vector{ComplexF64},PathResult},
}
    F::S
    L::Sub
    # only non-singular results or solutions
    R::Vector{R}
    projective::Bool
end

function WitnessSet(
    F::AbstractSystem,
    L::LinearSubspace,
    R;
    projective::Bool = is_linear(A) && is_homogeneous(System(F)),
)
    WitnessSet(F, L, R, projective)
end

"""
    system(W::WitnessSet)

Get the system stored in `W`.
"""
system(W::WitnessSet) = W.F

"""
    linear_subspace(W::WitnessSet)

Get the linear subspace stored in `W`.
"""
linear_subspace(W::WitnessSet) = W.L

"""
    solutions(W::WitnessSet)

Get the solutions stored in `W`.
"""
solutions(W::WitnessSet{A,B,PathResult}) where {A,B} = solutions(W.R)
solutions(W::WitnessSet{A,B,Vector{ComplexF64}}) where {A,B} = W.R

"""
    results(W::WitnessSet)

Get the results stored in `W`.
"""
results(W::WitnessSet{<:Any,<:Any,PathResult}) = W.R

"""
    degree(W::WitnessSet)

Returns the degree of the witness set `W`. This equals the number of solutions stored.
"""
ModelKit.degree(W::WitnessSet) = length(W.R)

"""
    dim(W::WitnessSet)

The dimension of the algebraic set encoded by the witness set.
"""
dim(W::WitnessSet) = codim(W.L)

"""
    codim(W::WitnessSet)

The dimension of the algebraic set encoded by the witness set.
"""
codim(W::WitnessSet) = dim(W.L)

function Base.show(io::IO, W::WitnessSet)
    print(io, "Witness set for dimension $(dim(W)) of degree $(degree(W))")
end

### Construct witness sets
"""
    witness_set(F; codim = nvariables(F) - length(F), dim = nothing, options...)

Compute a [`WitnessSet`](@ref) for `F` in the given dimension (resp. codimension)
by sampling a random (affine) linear subspace.
After constructing the system this calls [`solve`](@ref) with the provided `options`.

    witness_set(F, L; options...)

Compute [`WitnessSet`](@ref) for `F` and the (affine) linear subspace `L`.

    witness_set(W::WitnessSet, L; options...)

Compute a new [`WitnessSet`](@ref) with the (affine) linear subspace `L` by moving
the linear subspace stored in `W` to `L`.

### Example
```julia-repl
julia> @var x y;
julia> F = System([x^2 + y^2 - 5], [x, y])
System of length 1
 2 variables: x, y

 -5 + x^2 + y^2

julia> W = witness_set(F)
Witness set for dimension 1 of degree 2
```
"""
witness_set(F::System, args...; compile = COMPILE_DEFAULT[], kwargs...) =
    witness_set(fixed(F; compile = compile), args...; kwargs...)
function witness_set(
    F::AbstractSystem;
    target_parameters = nothing,
    dim = nothing,
    codim = nothing,
    options...,
)
    f = System(F)
    projective = is_homogeneous(f)
    if isnothing(dim) && isnothing(codim)
        dim = corank(F) - projective
    elseif !isnothing(codim) && projective
        codim += 1
    end
    L = rand_subspace(size(F, 2); dim = codim, codim = dim, affine = !projective)
    witness_set(
        F,
        L;
        target_parameters = target_parameters,
        projective = projective,
        options...,
    )
end

function witness_set(F::AbstractSystem, L::LinearSubspace; projective = nothing, options...)
    res = solve(F; target_subspace = L, options...)
    if isnothing(projective)
        WitnessSet(F, L, results(res; only_nonsingular = true))
    else
        WitnessSet(F, L, results(res; only_nonsingular = true); projective = projective)
    end
end

corank(F::AbstractSystem; kwargs...) = nvariables(F) - LA.rank(F)
function LinearAlgebra.rank(F::AbstractSystem; target_parameters = nothing)
    m, n = size(F)
    u = zeros(ComplexF64, m)
    U = zeros(ComplexF64, m, n)
    x = randn(ComplexF64, n)
    evaluate_and_jacobian!(u, U, F, x, target_parameters)
    LA.rank(U)
end

### Move witness sets around
function witness_set(W::WitnessSet, L::LinearSubspace; options...)
    if W.projective && !is_linear(L)
        error(
            "The given space is an affine linear subspace (``b ≠ 0``). " *
            " Expected a linear subspace since the given witness set is projective.",
        )
    end
    res = solve(W.F, W.R; start_subspace = W.L, target_subspace = L, options...)
    WitnessSet(W.F, L, results(res; only_nonsingular = true), is_linear(L) && W.projective)
end

"""
    trace_test(W::WitnessSet; options...)

Performs a trace test [^LRS18] to verify whether the given witness set `W` is complete.
Returns the trace of the witness set which should be theoretically be 0 if `W` is complete.
Due to floating point arithmetic this is not the case, thus is has to be manually checked
that the trace is sufficiently small.
Returns `nothing` if the trace test failed due to path tracking failures.
The `options` are the same as for calls to [`witness_set`](@ref).

```julia-repl
julia> @var x y;
julia> F = System([x^2 + y^2 - 5], [x, y])
System of length 1
 2 variables: x, y

 -5 + x^2 + y^2

julia> W = witness_set(F)
Witness set for dimension 1 of degree 2

julia> trace = trace_test(W)
9.981960497718987e-16
```

[^LRS18]: Leykin, Anton, Jose Israel Rodriguez, and Frank Sottile. "Trace test." Arnold Mathematical Journal 4.1 (2018): 113-125.
APA

"""
function trace_test(W₀::WitnessSet; options...)
    L₀ = linear_subspace(W₀)
    F = system(W₀)
    S₀ = solutions(W₀)
    # if we are in the projective setting, we need to make sure that
    # all solutions are on the same affine chart
    # Therefore make the affine chart now
    if W₀.projective
        F = on_affine_chart(F)
        s₀ = sum(s -> set_solution!(s, F, s), S₀)
    else
        s₀ = sum(S₀)
    end

    v = randn(ComplexF64, codim(L₀))
    L₁ = translate(L₀, v)
    L₋₁ = translate(L₀, -v)

    R₁ = solve(F, S₀; start_subspace = L₀, target_subspace = L₁, options...)
    nsolutions(R₁) == degree(W₀) || return nothing

    R₋₁ = solve(F, S₀; start_subspace = L₀, target_subspace = L₋₁, options...)
    nsolutions(R₋₁) == degree(W₀) || return nothing

    s₁ = sum(solutions(R₁))
    s₋₁ = sum(solutions(R₋₁))

    M = [s₋₁ s₀ s₁; 1 1 1]
    singvals = LA.svdvals(M)
    trace = singvals[3] / singvals[1]

    trace
end


text:
# One-factor Designs
# Y = A + e
#
# author: yizhan miao
# email: yzmiao@protonmail.com
# last update: Oct 25 2018

""" One-way ANOVA

Syntax:
- `RESULT, _F, _p = anova(Y::Array{T, 2}; preview::Bool) where {T <: Number}`
- `RESULT, _F, _p = anova(Y::Array{T, 1}, F::Array{D, 1}; preview::Bool) where {T <: Number, D<:Integer}`

Arguments:
- `Y`: data set in 1d vector, or 2d matrix, columns as treatments and rows as observations
- `F`: corresponding factor markers, as `Array{Integer, 1}`

Keywords:
- `preview`: print anova table [default: `true`]

Returns:
- ANOVA Result
- F statistic
- p value
"""
function anova(Y::Array{T, 2}; preview::Bool=true) where {T <: Number}
    _row, _col = size(Y)
    _mean_total = mean(Y)
    _mean_treat = mean(Y, dims = 1)

    _SS_total = (Y .- _mean_total) .^ 2 |> sum
    _SS_treat = _row * sum((_mean_treat .- _mean_total).^2)
    _SS_error = sum(broadcast(-, Y, _mean_treat).^2)

    _df_total = _row * _col - 1
    _df_treat = _col - 1
    _df_error = _col * (_row - 1)

    _MS_treat = _SS_treat / _df_treat
    _MS_error = _SS_error / _df_error

    _F = _MS_treat / _MS_error
    _F_dist = FDist(_df_treat, _df_error)

    _p = ccdf(_F_dist, _F)

    TOTAL = anova_entry(:total, _SS_total, _df_total, NaN, NaN, NaN)
    ERROR = anova_entry(:error, _SS_error, _df_error, _MS_error, NaN, NaN)
    TREAT = anova_entry(:columns, _SS_treat, _df_treat, _MS_treat, _F, _p)
    RESULT = anova_result("One-way Anova", [TREAT, ERROR, TOTAL])

    preview&&preview_anova_result(RESULT)

    RESULT, _F, _p
end

function anova(Y::Array{T, 1}, F::Array{D, 1}; preview::Bool=true) where {T <: Number, D<:Integer}
    # check balance
    _factor = unique(F)
    _n_factor = zeros(Int64, length(_factor))
    for (idx, each) in enumerate(_factor)
        _n_factor[idx] = F[F .== each] |> length
    end

    @assert(length(unique(_n_factor)) == 1, "input data is not balanced across factors.")
    # TODO: support unbalanced ANOVA

    _mean_total = mean(Y)
    _mean_treat = [mean(Y[F .== _F]) for _F in _factor]

    _SS_total = (Y .- _mean_total) |> x->x.^2 |> sum
    _SS_treat = (_mean_treat .- _mean_total).^2 .* _n_factor[1] |> sum
    _SS_error = _SS_total - _SS_treat

    _df_total = length(Y) - 1
    _df_treat = length(_factor) - 1
    _df_error = _df_total - _df_treat

    ##### one way ANOVA #####
    _MS_treat = _SS_treat / _df_treat
    _MS_error = _SS_error / _df_error

    _F = _MS_treat / _MS_error
    _F_dist = FDist(_df_treat, _df_error)

    _p = ccdf(_F_dist, _F)

    TOTAL = anova_entry(:total, _SS_total, _df_total, NaN, NaN, NaN)
    ERROR = anova_entry(:error, _SS_error, _df_error, _MS_error, NaN, NaN)
    TREAT = anova_entry(:columns, _SS_treat, _df_treat, _MS_treat, _F, _p)
    RESULT = anova_result("One-way Anova", [TREAT, ERROR, TOTAL])

    preview&&preview_anova_result(RESULT)

    RESULT, _F, _p
end


text:
module Portinari

export Area, Context, D3Attr, Line, Restyle, Scale, Shape

include("./js_base.jl")

include("./d3_abstractions.jl")

include("./context.jl")

include("./area.jl")
include("./line.jl")
include("./shape.jl")

end # module


text:
using DataFrames, Resampling

##########
#Fitting Funcs
##########
#build a little train function that just takes indices
function train(indices, lam::Float64, deg::Int, c::Float64,
                demand_data, flow_data, arcs; fcoeffs=nothing)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))
    m, coeffs, reg_term = setUpFitting(deg, c, [k for k=keys(demand_data)], 1:numNodes)

    for a in values(arcs)
        a.obsflow = a.trueflow
    end
    addIncreasingCnsts(m, coeffs, arcs, TOL=1e-8)  #uses the original obs flows

    avgCost = mean( [bpacost(a.obsflow, a.capacity, 1.0) for a in values(arcs)] )
    normalize(m, coeffs, [a.obsflow / a.capacity for a in values(arcs)], avgCost)

    # normalize(m, coeffs, 1., bpacost(1., 1., 1.))

    resids = Variable[]
    for i = indices
        #copy the flow data over to the arcs, demand data to demands (slow)
        for (ix, a) in enumerate(vArcs)
            a.obsflow = flow_data[ix, i]
        end
        for odpair in keys(demands)
            demands[odpair] = demand_data[odpair][i]
        end

        #Dual Feasibility
        ys = addNetworkCnsts(m, coeffs, demands, arcs, numNodes)

        #add the residual for this data point
        push!(resids, addResid(m, coeffs, ys, demands, arcs, 1e6))
    end

    if fcoeffs != nothing
        fixCoeffs(m, fcoeffs, coeffs)
    end
    @setObjective(m, Min, sum{resids[i], i = 1:length(resids)} + lam*reg_term)
    solve(m)
    println(getObjectiveValue(m) - lam * getValue(reg_term) )
    return [getValue(coeffs[i]) for i =1:length(coeffs)]
end

# a function for testing.
function test(fcoeffs, indices, demand_data, flow_data, arcs, g, vArcs)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))

    l1_eps = 0.
    for i = indices
        costs = [polyEval(fcoeffs, flow_data[ix, i]/a.capacity)*a.freeflowtime
                    for (ix, a) in enumerate(vArcs)]
        primal = dot(costs, flow_data[:, i])

        #solve a bunch of shortest path problems to get dual cost
        dual = 0.
        for odpair = keys(demand_data)
            #solve the shortest path problems, and update the total flow
            r = dijkstra_shortest_paths(g, costs, odpair[1] )
            dual += r.dists[odpair[2]] * demand_data[odpair][i]
        end
        l1_eps += abs(primal - dual)
    end
    return l1_eps / length(indices)
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BisectFractionPartition(normal, fraction=0.5, maxiter=10)

A method for partitioning spatial data into two half spaces
defined by a `normal` direction and a `fraction` of points.
The partition is returned within `maxiter` bisection iterations.
"""
struct BisectFractionPartition{T,N} <: PartitionMethod
  normal::SVector{N,T}
  fraction::Float64
  maxiter::Int

  function BisectFractionPartition{T,N}(normal, fraction, maxiter) where {N,T}
    new(normalize(normal), fraction, maxiter)
  end
end

BisectFractionPartition(normal::SVector{N,T}, fraction=0.5, maxiter=10) where {T,N} =
  BisectFractionPartition{T,N}(normal, fraction, maxiter)

BisectFractionPartition(normal::NTuple{N,T}, fraction=0.5, maxiter=10) where {T,N} =
  BisectFractionPartition(SVector(normal), fraction, maxiter)

function partition(object, method::BisectFractionPartition)
  bbox = boundbox(object)
  n = method.normal
  f = method.fraction
  c = center(bbox)
  d = diagonal(bbox)

  # maximum number of bisections
  maxiter = method.maxiter

  iter = 0; p = 0
  a = c - d/2 * n
  b = c + d/2 * n
  while iter < maxiter
    m = (a + b) / 2

    p = partition(object, BisectPointPartition(n, m))
    g = nelms(p[1]) / nelms(object)

    g ≈ f && break
    g > f && (b = m)
    g < f && (a = m)

    iter += 1
  end

  p
end


text:
"""
    GAMBAR

Gyromagnetic ratio for ¹H with units Hz/G.
"""
const GAMBAR = 4258

"""
    GAMMA

Gyromagnetic ratio for ¹H with units rad/s/G.
"""
const GAMMA  = 2π * GAMBAR

"""
    Position(x, y, z)

Create a mutable `Position` object representing a 3D location. Units are cm.

# Properties
- `x::Real`: x position
- `y::Real`: y position
- `z::Real`: z position
"""
mutable struct Position{T<:Real}
    x::T
    y::T
    z::T
end

Position(x, y, z) = Position(promote(x, y, z)...)

Base.show(io::IO, pos::Position) = print(io, "Position(", pos.x, ", ", pos.y, ", ", pos.z, ")")
Base.show(io::IO, ::MIME"text/plain", pos::Position{T}) where {T} =
    print(io, "Position{$T}:\n x = ", pos.x, " cm\n y = ", pos.y, " cm\n z = ", pos.z, " cm")

Base.convert(::Type{Position{T}}, p::Position) where {T} = Position(T(p.x), T(p.y), T(p.z))

Base.:(==)(p1::Position, p2::Position) = p1.x == p2.x && p1.y == p2.y && p1.z == p2.z

const DualFloat64 = Union{Float64,<:ForwardDiff.Dual}

"""
    AbstractSpin

Abstract type for representing individual spins.
"""
abstract type AbstractSpin end

"""
    Spin([M], M0, T1, T2, Δf, [pos]) <: AbstractSpin

Create an object that represents a single spin.

# Properties
- `M::Magnetization = Magnetization(0, 0, M0)`: Magnetization vector
- `M0::Real`: Equilibrium magnetization
- `T1::Real`: Spin-lattice recovery time constant (ms)
- `T2::Real`: Spin-spin recovery time constant (ms)
- `Δf::Real`: Off-resonance frequency (Hz)
- `pos::Position = Position(0, 0, 0)`: Spatial location (cm)
- `N::Int = 1`: Number of compartments

# Examples
```jldoctest
julia> Spin(1, 1000, 100, 0, Position(1, 2, 3))
Spin{Float64}:
 M = Magnetization(0.0, 0.0, 1.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 0.0 Hz
 pos = Position(1.0, 2.0, 3.0) cm

julia> Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0)
Spin{Float64}:
 M = Magnetization(1.0, 2.0, 3.0)
 M0 = 1.0
 T1 = 1000.0 ms
 T2 = 100.0 ms
 Δf = 0.0 Hz
 pos = Position(0.0, 0.0, 0.0) cm
```
"""
struct Spin{T<:DualFloat64} <: AbstractSpin
    M::Magnetization{T}
    M0::T
    T1::T
    T2::T
    Δf::T
    pos::Position{Float64}

    # Default constructor with optional argument pos
    function Spin{T}(
        M::Magnetization,
        M0,
        T1,
        T2,
        Δf,
        pos::Position = Position(0.0, 0.0, 0.0)
    ) where {T<:DualFloat64}

        new{T}(M, M0, T1, T2, Δf, pos)

    end
end

function Spin(M, M0, T1, T2, Δf, pos...)

    args = (M0, T1, T2, Δf)
    T = promote_type(Float64, eltype(M), typeof.(args)...)
    Spin{T}(M, args..., pos...)

end

# If magnetization vector is not specified then use equilibrium
Spin{T}(M0, T1, T2, Δf, pos::Position...) where {T} =
    Spin{T}(Magnetization(0.0, 0.0, M0), M0, T1, T2, Δf, pos...)
Spin(M0, T1, T2, Δf, pos::Position...) =
    Spin(Magnetization(0.0, 0.0, M0), M0, T1, T2, Δf, pos...)

function Base.show(io::IO, ::MIME"text/plain", spin::Spin{T}) where {T}

    print(io, "Spin{$T}:")
    print(io, "\n M = ", spin.M)
    print(io, "\n M0 = ", spin.M0)
    print(io, "\n T1 = ", spin.T1, " ms")
    print(io, "\n T2 = ", spin.T2, " ms")
    print(io, "\n Δf = ", spin.Δf, " Hz")
    print(io, "\n pos = ", spin.pos, " cm")

end

Base.eltype(::Spin{T}) where {T} = T

"""
    SpinMC([M], M0, frac, T1, T2, Δf, τ, [pos]) <: AbstractSpin

Create an object that represents a single spin with multiple compartments.

# Properties
- `M::MagnetizationMC = Meq`: Magnetization vector
- `Meq::MagnetizationMC = MagnetizationMC((0, 0, frac[1] * M0), ...)`:
  Equilibrium magnetization vector
- `M0::Real`: Equilibrium magnetization
- `frac::Tuple{<:Real}`: Water fraction of each compartment
- `T1::Tuple{<:Real}`: Spin-lattice recovery time constants (ms)
- `T2::Tuple{<:Real}`: Spin-spin recovery time constants (ms)
- `Δf::Tuple{<:Real}`: Off-resonance frequencies (Hz)
- `r::Tuple{Tuple{<:Real}}`: Exchange rates (1/ms); `r[i][j]` is the exchange
  rate from compartment `i` to compartment `j`
- `pos::Position = Position(0, 0, 0)`: Spatial location (cm)
- `N::Int = length(frac)`: Number of compartments

## Note
The `SpinMC` constructor takes `τ` (inverse exchange rate, or residence time) as
input, *not* `r`. Furthermore, `τ` is given as a `Tuple` with `N^2 - N`
elements, arranged like (τ12, τ13, ..., τ1N, τ21, τ23, ..., τ2N, ...).

# Examples
```jldoctest
julia> SpinMC(1, (0.2, 0.8), (400, 1000), (20, 100), (15, 0), (100, 25))
SpinMC{Float64,2}:
 M = MagnetizationMC((0.0, 0.0, 0.2), (0.0, 0.0, 0.8))
 M0 = 1.0
 frac = (0.2, 0.8)
 T1 = (400.0, 1000.0) ms
 T2 = (20.0, 100.0) ms
 Δf = (15.0, 0.0) Hz
 r = ((0.0, 0.01), (0.04, 0.0)) 1/ms
 pos = Position(0.0, 0.0, 0.0) cm
```
"""
struct SpinMC{T<:DualFloat64,N} <: AbstractSpin
    M::MagnetizationMC{T,N}
    Meq::MagnetizationMC{T,N}
    M0::T
    frac::NTuple{N,T}
    T1::NTuple{N,T}
    T2::NTuple{N,T}
    Δf::NTuple{N,T}
    r::NTuple{N,NTuple{N,T}}
    pos::Position{Float64}

    function SpinMC{T}(
        M::MagnetizationMC{S,N},
        M0,
        frac,
        T1,
        T2,
        Δf,
        τ,
        pos::Position = Position(0.0, 0.0, 0.0)
    ) where {T<:DualFloat64,S,N}

        N > 1 || error(sprint(print, "SpinMC expects 2 or more compartments, got ", N))
        Meq = MagnetizationMC{T}(N)
        for i = 1:N
            Meq[i].z = frac[i] * M0
        end
        itmp = 0
        r = ntuple(N) do i
            ntuple(N) do j
                if i == j
                    out = zero(T)
                else
                    itmp += 1
                    out = 1 / τ[itmp]
                end
                out
            end
        end
        new{T,N}(M, Meq, M0, frac, T1, T2, Δf, r, pos)

    end
end

function SpinMC(M, M0, frac, T1, T2, Δf, τ, pos...)

    frac = promote(frac...)
    T1 = promote(T1...)
    T2 = promote(T2...)
    Δf = promote(Δf...)
    rtmp = promote((1 ./ τ)...)
    T = promote_type(Float64, eltype(M), typeof(M0),
        eltype(frac), eltype(T1), eltype(T2), eltype(Δf), eltype(rtmp))
    SpinMC{T}(M, M0, frac, T1, T2, Δf, τ, pos...)

end

SpinMC{T}(M::NTuple{N,Magnetization}, args...) where {T,N} = SpinMC{T}(MagnetizationMC(M...), args...)
SpinMC{T}(M::NTuple{N,NTuple{3,Real}}, args...) where {T,N} = SpinMC{T}(MagnetizationMC(M...), args...)
SpinMC{T}(M0::Real, frac, T1, T2, Δf, τ, pos::Position...) where {T} =
    SpinMC{T}(MagnetizationMC((Magnetization(0.0, 0.0, frac[i] * M0) for i = 1:length(frac))...), M0, frac, T1, T2, Δf, τ, pos...)
SpinMC(M::NTuple{N,Magnetization}, args...) where {N} = SpinMC(MagnetizationMC(M...), args...)
SpinMC(M::NTuple{N,NTuple{3,Real}}, args...) where {N} = SpinMC(MagnetizationMC(M...), args...)
SpinMC(M0::Real, frac, T1, T2, Δf, τ, pos::Position...) =
    SpinMC(MagnetizationMC((Magnetization(0.0, 0.0, frac[i] * M0) for i = 1:length(frac))...), M0, frac, T1, T2, Δf, τ, pos...)

function Base.show(io::IO, spin::SpinMC{T,N}) where {T,N}

    print(io, "SpinMC{$T,$N}(", spin.M, ", ")
    print(io, spin.M0, ", ")
    print(io, spin.frac, ", ")
    print(io, spin.T1, ", ")
    print(io, spin.T2, ", ")
    print(io, spin.Δf, ", ")
    print(io, spin.r, ", ")
    print(io, spin.pos, ")")

end

function Base.show(io::IO, ::MIME"text/plain", spin::SpinMC{T,N}) where {T,N}

    print(io, "SpinMC{$T,$N}:")
    print(io, "\n M = ", spin.M)
    print(io, "\n M0 = ", spin.M0)
    print(io, "\n frac = ", spin.frac)
    print(io, "\n T1 = ", spin.T1, " ms")
    print(io, "\n T2 = ", spin.T2, " ms")
    print(io, "\n Δf = ", spin.Δf, " Hz")
    print(io, "\n r = ", spin.r, " 1/ms")
    print(io, "\n pos = ", spin.pos, " cm")

end

Base.eltype(::SpinMC{T,N}) where {T,N} = T

signal(spin::AbstractSpin) = signal(spin.M)

function Base.getproperty(spin::Spin, s::Symbol)

    if s === :N
        return 1
    else
        return getfield(spin, s)
    end

end

function Base.getproperty(spin::SpinMC{T,N}, s::Symbol) where {T,N}

    if s === :N
        return N
    else
        return getfield(spin, s)
    end

end


text:
module DarcyCubedSphereTests

# This program solves Problem in Section 5.1.1 in the following paper
# Rognes, M. E., Ham, D. A., Cotter, C. J., and McRae, A. T. T.
# Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2,
# Geosci. Model Dev., 6, 2099–2119, https://doi.org/10.5194/gmd-6-2099-2013, 2013.

using Test
using Gridap
using GridapGeosciences
using Plots

include("ConvergenceAnalysisTools.jl")

function solve_darcy(model,order,degree)
  g(x) = sin(0.5*π*x[2]) # x[1]*x[2]*x[3]
  RT=ReferenceFE(raviart_thomas,Float64,order)
  DG=ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,RT; conformity=:Hdiv)
  Q = FESpace(model,DG; conformity=:L2)
  T = FESpace(model,DG; conformity=:L2)

  U = TrialFESpace(V)
  P = TrialFESpace(Q)
  R = TrialFESpace(T)

  Y = MultiFieldFESpace([V, Q, T])
  X = MultiFieldFESpace([U, P, R])

  trian = Triangulation(model)
  dΩ = Measure(trian,degree)
  dω = Measure(trian,degree,ReferenceDomain())

  aref((u,p,r),(v,q,t)) = ∫( v⋅u + q*r + t*p )*dΩ +  ∫(q*DIV(u)+DIV(v)*p)*dω
  aphy((u,p,r),(v,q,t)) = ∫( v⋅u + q*r + t*p + q*(∇⋅u) + (∇⋅v)*p)*dΩ
  b((v,q,t)) = ∫(q*g)*dΩ

  op = AffineFEOperator(aref,b,X,Y)
  xh = solve(op)
  uh, ph, rh = xh
  # writevtk(trian, "u", cellfields=["uh"=>uh])
  e = g-rh
  sqrt(sum(∫(e*e)dΩ))
end

# Testing cubed sphere mesh with analytical geometric mapping
@time hs0,k0errors,s0=convergence_study(solve_darcy,generate_n_values(2),0,2)
@test round(s0,digits=3) ≈ 1.0

@time hs1,k1errors,s1=convergence_study(solve_darcy,generate_n_values(2,n_max=50),1,4)
@test round(s1,digits=3) ≈ 1.932

# Testing cubed sphere mesh with polynomial geometric mapping
@time hs0,k0errors,s0=convergence_study(solve_darcy,generate_n_values(2),2,0,2)
@test round(s0,digits=3) ≈ 1.0

@time hs1,k1errors,s1=convergence_study(solve_darcy,generate_n_values(2,n_max=50),2,1,4)
@test round(s1,digits=3) ≈ 1.936

plot([hs0,hs1],[k0errors,k1errors],
     xaxis=:log, yaxis=:log,
     label=["L2 k=0","L2 k=1"],
     shape=:auto,
     xlabel="h",ylabel="L2 error norm")

end # module


text:
using ArgParse2

function main()
    parser = ArgumentParser()
    add_argument(parser, "arg1", help = "a positional argument")
    add_argument(parser, "--opt1", help = "an option with an argument")
    add_argument(parser, "-o", "--opt2", help = "another option with an argument")
    add_argument(parser, "--flag1", action="store_true", default=false, help = "an option without argument, i.e. a flag")

    println("Parsed args:")
    for (arg,val) in pairs(parse_args(parser))
        println("  $arg  =>  $val")
    end
end

main()


text:
export Eulerian
"""
`Eulerian(n,k)` returns the number of permutations of `{1,2,...,n}`
with `k` ascents.
"""
function Eulerian(n::Integer, k::Integer)::BigInt
    @assert (n>=0 && k>=0) "$n,$k must both be nonnegative"

    if n==0
        return big(0)
    end

    if k>n || k==0
        return big(0)
    end

    if n < 2
        return big(1)
    end

    if k==n      # includes the case n=k=0
        return big(1)
    end

    if k==1
        return big(1)
    end

    if _has(Eulerian,(n,k))
        return _get(Eulerian,(n,k))
    end

    val = (n-k+1)*Eulerian(n-1,k-1) + k*Eulerian(n-1,k)
    _save(Eulerian,(n,k),val)
    return val
end

_make(Eulerian,Tuple{Integer,Integer})