text:
#=
Handling of time in data as a physical quantity, and time-related data processing
=#
using Statistics, StatsBase
export monthday_indices, maxyearspan, daymonth, realtime_days, realtime_milliseconds
export temporal_sampling
export timemean, timeagg
export monthlyagg, yearlyagg, temporalrange, seasonalyagg, season
export DAYS_IN_ORBIT, HOURS_IN_ORBIT
export seasonality
#########################################################################
# Datetime related
#########################################################################
using Dates

const DAYS_IN_ORBIT = 365.26
const HOURS_IN_ORBIT = 365.26*24

"daymonth(t) = day(t), month(t)"
daymonth(t) = day(t), month(t)

maxyearspan(A::AbstractDimArray, tsamp = temporal_sampling(A)) =
maxyearspan(dims(A, Time).val, tsamp)

"""
    temporal_sampling(x) → symbol
Return the temporal sampling type of `x`, which is either an array of `Date`s or
a dimensional array (with `Time` dimension).

Possible return values are:
- `:hourly`, where the temporal difference between successive entries is exactly 1 hour.
- `:daily`, where the temporal difference between successive entries is exactly 1 day.
- `:monthly`, where all dates have the same day, but different month.
- `:yearly`, where all dates have the same month and day, but different year.
- `:other`, which means that `x` doesn't fall to any of the above categories.
"""
temporal_sampling(A::AbstractDimArray) = temporal_sampling(dims(A, Time).val)
temporal_sampling(t::Dimension) = temporal_sampling(t.val)

function temporal_sampling(t::AbstractVector{<:TimeType})
    function issame(f)
        x0 = f(t[1])
        return all(i -> f(t[i]) == x0, 2:length(t))
    end
    samehour, sameday, samemonth, sameyear = map(issame, (hour, day, month, year))
    constdailydiff = all(i -> day(t[i])-day(t[i-1]) ∈ (1, -30, -27, -28, -29), 2:length(t))

    if !samehour
        # check if they have exactly 1 hour difference
        dh = getproperty.(diff(t), :value)
        return all(isequal(3600000), dh) ? :hourly : :other
    elseif !sameday && samemonth
        return :daily
    elseif samehour && !sameday && !samemonth && constdailydiff
        # this clause checks daily data that span more than 1 month
        return :daily
    elseif sameday && !samemonth
        return :monthly
    elseif sameday && samemonth && !sameyear
        return :yearly
    else
        return :other
    end
end
temporal_sampling(t::AbstractVector) = :other
temporal_sampling(t::StepRange{<:Any,Month}) = :monthly
temporal_sampling(t::StepRange{<:Any,Year}) = :yearly
temporal_sampling(t::StepRange{<:Any,Day}) = :daily
temporal_sampling(t::StepRange{<:Any,Hour}) = :hourly
temporal_sampling(t::StepRange{<:Any,<:Any}) = :other

"return the appropriate subtype of `Dates.Period` or `nothing`."
function tsamp2period(tsamp)
    tsamp == :monthly && return Month(1)
    tsamp == :yearly && return Year(1)
    tsamp == :daily && return Day(1)
    tsamp == :hourly && return Hour(1)
    tsamp == :other && return nothing
end

"""
    maxyearspan(A::ClimArray) = maxyearspan(dims(A, Time))
    maxyearspan(t::Vector{<:DateTime}) → i
Find the maximum index `i` of `t` so that `t[1:i]` covers exact(*) multiples of years.

(*) For monthly spaced data `i` is a multiple of `12` while for daily data we find
the largest possible multiple of `DAYS_IN_ORBIT`.
"""
function maxyearspan(times, tsamp = temporal_sampling(times))
    l = length(times)
    if tsamp == :monthly
        l % 12 == 0 && return l
        m = month(times[1])
        x = findlast(i -> month(times[i]) == m, 1:l) - 1
        if x != 0
            return x
        elseif x == 0
            @warn "Caution: data does not cover a full year."
            return l
        end
    elseif tsamp == :yearly
        return length(times)
    elseif tsamp == :daily
        n_max = Int(l÷DAYS_IN_ORBIT)
        nb_years = findlast(n -> round(Int, n * DAYS_IN_ORBIT) ≤ l, 1:n_max)
        if nb_years != nothing
            return round(Int,nb_years * DAYS_IN_ORBIT)-1
        elseif nb_years == nothing
            @warn "Caution: data does not cover a full year."
            return l
        end
    elseif tsamp == :hourly
        n_max = Int(l÷HOURS_IN_ORBIT)
        nb_years = findlast(n -> round(Int, n * HOURS_IN_ORBIT) ≤ l, 1:n_max)
        if nb_years != nothing
            return round(Int,nb_years * HOURS_IN_ORBIT)-1
        elseif nb_years == nothing
            @warn "Caution: data does not cover a full year."
            return l
        end
    else
        error("maxyearspan: not implemented yet for $(tsamp)-sampled data")
    end
end

"""
    monthday_indices(times, date = times[1])
Find the indices in `times` (which is a `Vector{Date}`) at which
the date in `times` gives the same day and month as `date`.
"""
function monthday_indices(times, date = times[1])
    d1, m1 = daymonth(date)
    a = findall(i -> daymonth(times[i]) == (d1, m1), 1:length(times))
end

"""
    monthamount(x::AbstractArray{<:AbstractDateTime})
Convert `x` to the amount of months (starting from 0).
"""
function monthamount(t)
    L = length(t)
    yd = [year(t[i]) - year(t[1]) for i in 1:L]
    md = [month(t[i]) - month(t[1]) for i in 1:L]
    [yd[i]*12 + md[i] for i in 1:L]
end

"""
    monthspan(t::TimeType)
Return a vector of daily spaced `Date`s that span the entire month that `t` belongs in.
"""
function monthspan(t::TimeType)
    m = mod1(month(t), 12)
    n = mod1(m+1, 12)
    y = year(t)
    u = m == 12 ? y+1 : y
    d = collect(Date(y, m, 1):Day(1):Date(u, n, 1))[1:end-1]
end

"""
    realtime_days(t::AbstractVector{<:TimeType}, T = Float32)
Convert the given _sequential_ date time vector `t` in a vector in a format of "real time",
where time is represented by real numbers, increasing cumulatively, as is the case when
representing a timeseries `x(t)`.
As only differences matter in this form, the returned vector always starts from 0.
The measurement unit of time here is days.

For temporal sampling less than daily return `realtime_milliseconds(t) ./ (24*60*60*1000)`.

Example:
```juliarepl
julia> t = Date(2004):Month(1):Date(2004, 6)
Date("2004-01-01"):Month(1):Date("2004-06-01")

julia> realtime_days(t)
6-element Vector{Float32}:
   0.0
  29.0
  60.0
  90.0
 121.0
 151.0
```
"""
function realtime_days(t::AbstractArray{<:TimeType}, T = Float32)
    @assert issorted(t)
    ts = temporal_sampling(t)
    if ts == :monthly
        truetime = cumsum(daysinmonth.(t))
        return T.(truetime .- truetime[1])
    elseif ts == :yearly
        truetime = cumsum(daysinyear.(t))
        return T.(truetime .- truetime[1])
    elseif ts == :daily
        return T.(0:length(t)-1)
    else
        return T.(realtime_milliseconds(t) ./ 86400000)
    end
end
realtime_days(A) = realtime_days(dims(A, Ti).val)

"""
    realtime_milliseconds(t::AbstractArray{<:TimeType}, T = Float64)
Similar with [`realtime_days`](@ref), but now the measurement unit is millisecond.
For extra accuracy, direct differences in `t` are used.
"""
function realtime_milliseconds(t::AbstractArray{<:TimeType}, T = Float64)
    @assert issorted(t)
    r = cumsum([T(x.value) for x in diff(t)])
    pushfirst!(r, 0)
    return r
end
realtime_milliseconds(A) = realtime_milliseconds(dims(A, Ti).val)

#########################################################################
# temporal statistics
#########################################################################
"""
    timemean(A::ClimArray [, w]) = timeagg(mean, A, w)
Temporal average of `A`, see [`timeagg`](@ref).
"""
timemean(A::ClimArray, w = nothing) = timeagg(mean, A, w)

"""
    timeagg(f, A::ClimArray, W = nothing)
Perform a proper temporal aggregation of the function `f` (e.g. `mean, std`) on `A` where:
* Only full year spans of `A` are included, see [`maxyearspan`](@ref)
  (because most processes are affected by yearly cycle,
  and averaging over an uneven number of cycles typically results in artifacts)
* Each month in `A` is weighted with its length in days (for monthly sampled data)

If you don't want these features, just do [`dropagg`](@ref)`(f, A, Time, W)`.
This is also done in the case where the time sampling is unknown.

`W` are possible statistical weights that are used in conjuction to the temporal weighting,
to weight each time point differently.
If they are not a vector (a weight for each time point), then they have to be a dimensional
array of same dimensional layout as `A` (a weight for each data point).

See also [`monthlyagg`](@ref), [`yearlyagg`](@ref), [`seasonalyagg`](@ref).

    timeagg(f, t::Vector, x::Vector, w = nothing)
Same as above, but for arbitrary vector `x` accompanied by time vector `t`.
"""
function timeagg(f, A::AbDimArray, w = nothing)
    !hasdim(A, Time) && error("Array does not have `Time` dimension!")
    w isa AbDimArray && @assert dims(w) == dims(A)
    w isa Vector && @assert length(w) == size(A, Time)
    tsamp = temporal_sampling(A)
    if tsamp == :other
        return dropagg(f, A, Time) # TODO: Extend dropagg to arbitrary weights
    end
    r = if tsamp == :daily
        timeagg_daily(f, A, w)
    elseif tsamp == :monthly
        timeagg_monthly(f, A, w)
    elseif tsamp == :yearly
        timeagg_yearly(f, A, w)
    end
    R = ClimArray(r, otherdims(A, Time()), A.name)
end

function timeagg_yearly(f, A, w)
     if isnothing(w)
        dropagg(f, A, Time)
    elseif w isa AbDimArray
        map(i -> f(view(A, i), weights(view(W, i))), otheridxs(A, Time()))
    elseif w isa Vector
        fw = weights(w)
        map(i -> f(view(A, i), fw), otheridxs(A, Time()))
    end
end

function timeagg_monthly(f, A::AbDimArray, w)
    t = dims(A, Time).val
    mys = maxyearspan(t, :monthly)
    tw = daysinmonth.(t)
    W = if isnothing(w)
        tw
    elseif w isa Vector
        tw .* w
    elseif w isa AbDimArray
        _w = dimwise(*, w, ClimArray(tw, (Time(t),)))
        @view _w[Time(1:mys)]
    end
    other = otherdims(A, Time)
    _A = @view A[Time(1:mys)]
    if w isa AbDimArray
        r = map(i -> f(view(_A, i), weights(view(W, i))), otheridxs(A, Time()))
    else
        fw = weights(view(W, 1:mys))
        r = map(i -> f(view(_A, i), fw), otheridxs(A, Time()))
    end
end

function timeagg_daily(f, A::AbDimArray, w)
    t = dims(A, Time).val
    mys = maxyearspan(t)
    _A = view(A, Time(1:mys))
    if w isa AbDimArray
        W = view(w, Time(1:mys))
        r = map(i -> f(view(_A, i), weights(view(W, i))), otheridxs(A, Time()))
    elseif w isa Vector
        fw = weights(view(w, 1:mys))
        r = map(i -> f(view(_A, i), fw), otheridxs(A, Time()))
    elseif isnothing(w)
        r = dropagg(f, _A, Time)
    end
end

function timeagg(f, a::AbDimArray{T, 1}, w = nothing) where {T} # version with only time dimension
    t = dims(a, Time).val
    return timeagg(f, t, Array(a), w)
end

function timeagg(f, T::AbstractVector{<:TimeType}, a::Vector, w = nothing) # version with just vectors
    tsamp = temporal_sampling(T)
    mys = maxyearspan(T, tsamp)
    t = view(T, 1:mys)
    if tsamp == :monthly
        dimw = float.(daysinmonth.(t))
        !isnothing(w) && (dimw .*= view(w, 1:mys))
        return f(view(a, 1:mys), weights(dimw))
    else
        if isnothing(w)
            return f(view(a, 1:mys))
        else
            return f(view(a, 1:mys), weights(view(w, 1:mys)))
        end
    end
end

#########################################################################
# Monthly/yearly/daily/seasonal means
#########################################################################
"""
    monthlyagg(A::ClimArray, f = mean; mday = 15) -> B
Create a new array where the temporal information has been aggregated into months
using the function `f`.
The dates of the new array always have day number of `mday`.
"""
function monthlyagg(A::ClimArray, f = mean; mday = 15)
    t0 = dims(A, Time).val
    startdate = Date(year(t0[1]), month(t0[1]), mday)
    finaldate = Date(year(t0[end]), month(t0[end]), mday+1)
    t = startdate:Month(1):finaldate
    tranges = temporalrange(t0, Dates.month)
    return timegroup(A, f, t, tranges, "monthly")
end

"""
    yearlyagg(A::ClimArray, f = mean) -> B
Create a new array where the temporal information has been aggregated into years
using the function `f`.
By convention, the dates of the new array always have month and day number of `1`.
"""
function yearlyagg(A::ClimArray, f = mean)
    t0 = dims(A, Time).val
    startdate = Date(year(t0[1]), 1, 1)
    finaldate = Date(year(t0[end]), 2, 1)
    t = startdate:Year(1):finaldate
    tranges = temporalrange(t0, Dates.year)
    return timegroup(A, f, t, tranges, "yearly")
end

function timegroup(A, f, t, tranges, name)
    other = otherdims(A, Time)
    n = A.name == Symbol("") ? A.name : Symbol(A.name, ", $(name)")
    B = ClimArray(zeros(eltype(A), length.(other)..., length(t)), (other..., Time(t)), n)
    for i in 1:length(tranges)
        B[Time(i)] .= dropagg(f, view(A, Time(tranges[i])), Time)
    end
    return B
end

"""
    temporalrange(A::ClimArray, f = Dates.month) → r
    temporalrange(t::AbstractVector{<:TimeType}}, f = Dates.month) → r
Return a vector of ranges so that each range of indices are values of `t` that
belong in either the same month, year, day, or season, depending on `f`.
`f` can take the values: `Dates.year, Dates.month, Dates.day` or `season`
(all are functions).

Used in e.g. [`monthlyagg`](@ref), [`yearlyagg`](@ref) or [`seasonalyagg`](@ref).
"""
function temporalrange(t::AbstractArray{<:TimeType}, f = Dates.month)
    @assert issorted(t) "Sorted time required."
    L = length(t)
    r = Vector{UnitRange{Int}}()
    i, x = 1, f(t[1]) # previous entries
    for j in 2:L
        y = f(t[j])
        x == y && continue
        push!(r, i:(j-1))
        i, x = j, y
    end
    push!(r, i:L) # final range not included in for loop
    return r
end
temporalrange(A::AbstractDimArray, f = Dates.month) = temporalrange(dims(A, Time).val, f)

"""
    seasonalyagg(A::ClimArray, f = mean) -> B
Create a new array where the temporal information has been aggregated into seasons
using the function `f`.
By convention, seasons are represented as Dates spaced 3-months apart, where only the
months December, March, June and September are used to specify the date, with day 1.
"""
function seasonalyagg(A::ClimArray, f = mean)
    t0 = dims(A, Time).val
    startdate = to_seasonal_date(t0[1])
    finaldate = to_seasonal_date(t0[end])
    t = startdate:Month(3):finaldate
    tranges = temporalrange(t0, season)
    return timegroup(A, f, t, tranges, "seasonaly")
end

function to_seasonal_date(t)
    y, m = year(t), month(t)
    if m ∈ 3:5
        return Date(y, 3, 1)
    elseif m ∈ 6:8
        return Date(y, 6, 1)
    elseif m ∈ 9:11
        return Date(y, 9, 1)
    elseif m == 12
        return Date(y, 12, 1)
    elseif m ∈ 1:2
        return Date(y-1, 12, 1)
    end
end

"""
    season(date) → s::Int
Return the season of the given date, 1 for DJF, 2 for MAM, 3 for JJA, 4 for SON.
Complements functions like `Dates.year, Dates.month, Dates.day`.
"""
function season(t::Dates.AbstractTime)
    m = month(t)
    if m ∈ 3:5
        2
    elseif m ∈ 6:8
        3
    elseif m ∈ 9:11
        4
    else
        1
    end
end

#########################################################################
# Other advanced functions
#########################################################################
"""
    seasonality(t, x; y0 = year(t[1])) → dates, vals
Calculate the "seasonality" of a vector `x` defined with respect to a datetime
vector `t` and return `dates, vals`.
`dates` are all unique dates present in `t` *disregarding the year* (so only the
month and day are compared). The `dates` have as year entry `y0`.
`vals` is a vector of vectors, where `vals[i]` are all the values of `x` that have
day and month same as `dates[i]`. The elements of `vals` are sorted as encountered in `x`.

Typically one is interested in `mean.(vals)`, which actually is the seasonality,
and `std.(vals)` which is the interannual variability at each date.

    seasonality(A::ClimArray) → dates, vals
If given a `ClimArray`, then the array must have only one dimension (time).
"""
function seasonality(t, xs; y0 = year(t[1]))
    dates = unique(daymonth.(t))
    vals = [eltype(xs)[] for u in dates]
    for (τ, x) in zip(t, xs)
        i = findfirst(isequal(daymonth(τ)), dates)
        push!(vals[i], x)
    end
    outdates = [Date(y0, mon, day) for (day, mon) in dates]
    si = sortperm(outdates)
    return outdates[si], vals[si]
end

function seasonality(A::ClimArray; kwargs...)
    @assert length(dims(A)) == 1
    return seasonality(dims(A, Time).val, A.data; kwargs...)
end


text:
export iso, iso_matrix, is_iso, info_map, uhash, degdeg, fast_iso_test

const iso_err_msg = "The graphs are not isomorphic"

"""
`iso_matrix(G,H)` returns a permutation matrix `P` such that
`A*P==P*B` where `A` is the adjacency matrix of `G` and `B` is the
adjacency matrix of `H`.  If the graphs are not isomorphic, an
error is raised.
"""
function iso_matrix(G::SimpleGraph, H::SimpleGraph)
    n = NV(G)

    if !fast_iso_test_basic(G,H)
        error(iso_err_msg)
    end

    m = Model(with_optimizer(_SOLVER.Optimizer;_OPTS...))
    @variable(m,P[1:n,1:n],Bin)
    A = adjacency(G)
    B = adjacency(H)

    for i=1:n
        for k=1:n
            @constraint(m,
                           sum(A[i,j]*P[j,k] for j=1:n) ==
                           sum(P[i,j]*B[j,k] for j=1:n)
                           )
        end
        @constraint(m, sum(P[i,j] for j=1:n) == 1)
        @constraint(m, sum(P[j,i] for j=1:n) == 1)
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error(iso_err_msg)
    end

    return Int.(value.(P))
end

"""
`iso(G,H)` finds an isomorphism from `G` to `H` (or throws an error if
the graphs are not isomorphic). Returns a dictionary mapping the
vertices of `G` to `H`.
"""

"""
`is_iso(G,H,d)` checks if `d` is an isomorphism from `G` to `H`.
"""
function is_iso(G::SimpleGraph, H::SimpleGraph, d::Dict)::Bool
    n = NV(G)

    # standard quick check
    if n != NV(H) || NE(G) != NE(H) || deg(G) != deg(H) || length(d) != n
        return false
    end

    # check keys and values in d match VG and VH, respectively
    for k in keys(d)
        if !has(G,k)
            return false
        end
    end

    for v in values(d)
        if !has(H,v)
            return false
        end
    end

    # Check edges match. Only need to go one way since we determined
    # that the two graphs have the same number of edges.
    EE = elist(G)
    for e in EE
        v = e[1]
        w = e[2]
        if !has(H,d[v],d[w])
            return false
        end
    end

    return true
end

"""
`is_iso(G,H)` returns `true` if the two graphs are isomorphic and `false`
otherwise.
"""
function is_iso(G::SimpleGraph, H::SimpleGraph)::Bool
    try
        f = iso(G,H)
    catch
        return false
    end
    return true
end

"""
`fast_iso_test_basic(G,H)` is a very quick test to rule out graphs being
isomorphic by considering the number of vertices, number of edges, and
degree sequences. Returns `false` if the graphs fail this very basic
test of isomorphism. A `true` result does *not* imply the graphs are
isomorphic.
"""
function fast_iso_test_basic(G::SimpleGraph, H::SimpleGraph)::Bool
    if NV(G)!=NV(H) || NE(G) != NE(H) || deg(G) != deg(H)
        return false
    end
    return true
end

"""
`fast_iso_test(G,H)` applies various heuristics to see if
the graphs *might* be isomorphic. A `false` return certifies the
graphs are **not** isomorphic; a `true` result indicates the
graphs might be (indeed, likely are) isomorphic.
"""
function fast_iso_test(G::SimpleGraph, H::SimpleGraph)::Bool
    return fast_iso_test_basic(G,H) && uhash(G) == uhash(H)
end

"""
`degdeg(G)` creates an n-by-n matrix where the nonzero entries in a
row are the degrees of the neighbors of that vertex. The rows are
lexicographically sorted.
"""
function degdeg(G::SimpleGraph)
    if cache_check(G,:degdeg)
      return cache_recall(G,:degdeg)
    end

    n = NV(G)
    result = zeros(Int,n,n)
    VV = vlist(G)

    for k=1:n
        v = VV[k]
        dv = [ deg(G,w) for w in G[v] ]
        dv = sort(dv)
        dv = [ dv; zeros(Int,n-deg(G,v)) ]
        result[k,:] = dv
    end

    result = sortslices(result,dims=1)
    cache_save(G,:degdeg,result)
    return result
end

"""
`info_map(G)`:
We create a dictionary mapping the vertices of `G` to 128-bit integer
values in such a way that twin vertices *will* have the same value
but, we hope, nontwin vertices will have different values. (By *twin*
we mean a pair of vertices such that there is an automorphism of the
graph mapping one to the other.)
"""
function info_map(G::SimpleGraph)
    if cache_check(G,:info_map)
      return cache_recall(G,:info_map)
    end
    n = NV(G)
    VV = vlist(G)
    H = G'

    # Section 1: Neighborhood degrees
    M1 = zeros(Int,n,n-1)
    for k=1:n
        v = VV[k]
        rv = [ deg(G,w) for w in G[v] ]
        rv = [ sort(rv) ; zeros(Int,n-1-deg(G,v)) ]
        M1[k,:] = rv
    end

    # Section 2: Complement degrees
    M2 = zeros(Int,n,n-1)
    for k=1:n
        v = VV[k]
        rv = [ deg(H,w) for w in H[v] ]
        rv = [ sort(rv) ; zeros(Int,n-1-deg(H,v)) ]
        M2[k,:] = rv
    end

    # Section 3: Distances
    M3 = zeros(Int,n,n-1)
    for k=1:n
        v = VV[k]
        dv = dist(G,v)
        M3[k,:] = sort(collect(values(dv)))[1:n-1]
    end

    # Section 4: Distances in complement
    M4 = zeros(Int,n,n-1)
    for k=1:n
        v = VV[k]
        dv = dist(H,v)
        M4[k,:] = sort(collect(values(dv)))[1:n-1]
    end

    M = [M1 M2 M3 M4]

    # println(sortrows(M))  # debug

    T = vertex_type(G)
    d = Dict{T,Int128}()
    for k=1:n
        v = VV[k]
        d[v] = Int128(hash(M[k,:]))
    end
    cache_save(G,:info_map,d)
    return d
end

"""
`iso(G,H)` finds an isomorphism from `G` to `H` (or throws an error if
the graphs are not isomorphic). Returns a dictionary mapping the
vertices of `G` to `H`.
"""
function iso(G::SimpleGraph, H::SimpleGraph)

    # quickly rule out some easily detected nonisomorphic graphs
    if !fast_iso_test(G,H) || uhash(G) != uhash(H)
        error(iso_err_msg)
    end

    VG = vlist(G)
    VH = vlist(H)
    n = NV(G)

    TG = vertex_type(G)
    TH = vertex_type(H)

    dG = info_map(G)
    dH = info_map(H)

    # Create mappings from vertex key values back to the vertices themselves

    xG = Dict{Int128,Set{TG}}()
    xH = Dict{Int128,Set{TH}}()

    valsG = sort(collect(values(dG)))
    valsH = sort(collect(values(dH)))

    for x in valsG
        xG[x] = Set{TG}()
        xH[x] = Set{TH}()
    end

    for v in VG
        x = dG[v]
        push!(xG[x],v)
    end

    for v in VH
        x = dH[v]
        push!(xH[x],v)
    end

    MOD = Model(with_optimizer(_SOLVER.Optimizer;_OPTS...))

    @variable(MOD, x[VG,VH],Bin)

    for v in VG
        @constraint(MOD, sum(x[v,VH[k]] for k=1:n)==1)
    end

    for w in VH
        @constraint(MOD, sum(x[VG[k],w] for k=1:n)==1)
    end

    for v in VG
        for w in VH
            @constraint(MOD,
                           sum( has(G,v,VG[k])*x[VG[k],w] for k=1:n ) ==
                           sum( x[v,VH[k]]*has(H,VH[k],w) for k=1:n )
                           )
        end
    end

    # Add contraints based on vertex numbers

    for val in unique(valsG)
        SG = collect(xG[val])
        SH = collect(xH[val])
        s  = length(SG)
        @constraint(MOD,
                       sum(x[u,v] for u in SG for v in SH)==s
                       )
    end

    optimize!(MOD)
    status = Int(termination_status(MOD))

    if status != 1
        error(iso_err_msg)
    end

    X = value.(x)

    result = Dict{vertex_type(G), vertex_type(H)}()

    for v in VG
        for w in VH
            if X[v,w] > 0
                result[v] = w
            end
        end
    end
    return result
end

using LinearAlgebra

function matrix_moments(A, max_exp::Int = 10)
    result = zeros(Int, max_exp)
    B = copy(A)

    for t=1:max_exp
        B = A*B
        result[t] = tr(B)
    end
    return result
end

"""
`uhash(G)` creates an `UInt64` hash of the graph such that isomorphic
graphs will produce the same value. We hope that nonisomorphic graphs
will create different values, but, alas, that need not be the case.
"""
function uhash(G::SimpleGraph)
    if cache_check(G,:uhash)
      return cache_recall(G,:uhash)
    end

    v1 = sort(collect(values(info_map(G))))

    A = adjacency(G)
    B = deepcopy(A)
    depth = min(10,NV(G))

    v2 = matrix_moments(adjacency(G))
    v3 = matrix_moments(laplace(G))

    vv = [v1; v2; v3]

    x = hash(vv)
    cache_save(G,:uhash,x)
    return x
end

ext: jl


text:
export IRBasis, FiniteTempBasis, MatsubaraSampling, TauSampling
export FiniteTempBasisSet
export fermion, boson, statistics_sign
export evaluate, fit, overlap, finite_temp_bases
export LogisticKernel, RegularizedBoseKernel
export SparsePoleRepresentation, from_IR, to_IR
export LegendreBasis
export CompositeBasis


text:
# This file contains a function to reset
# an instance to default values

"""
    reset!(Gen::GenStruct)

Reset an instance of [`GenStruct`](@ref) to default values.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Gen.example(file)
s.Gen.read!(file)
s.Gen.reset!()

# output

```
"""
function reset!(Gen::GenStruct)

    # Reset all elements
    Gen.N = Gen.r = 0
    Gen.Δt = Gen.q = Gen.α = Gen.β = Gen.γ = 0
    Gen.A = Gen.ν = Gen.ϕ = []

    return nothing

end


text:
function foundNear(n)
    near = true
    for j in 1:n
        psi,vort = randVortexField(1)
        vortfound = findvortices(psi)
        vfdata = rawData(vortfound)
        vdata = rawData(vort)
        near *= isapprox(vdata,vfdata,rtol = 0.2)
    end
    return near
end
function findwhere(A)
    I = findall(!iszero,A)
    v = A[I]
    ix = [I[i][1] for i in eachindex(I)]
    iy = [I[i][2] for i in eachindex(I)]
    return ix,iy,v
end

function findvortices_jumps(psi::Field;shift=true)
    @unpack x,y,ψ = psi
    phase = angle.(ψ)

    # Compute all nearest neighbour jumps
    diffx = phasejumps(phase,1); diffy = phasejumps(phase,2)

    # Compute all plaquette loop integrals with in-place memory recycling
    circshift!(phase,diffx,(0,1))
    diffx .-= phase; diffx .-= diffy
    circshift!(phase,diffy,(1,0))
    diffx .+= phase

    # Find all windings
    ixp,iyp,vp = findwhere(diffx .> 0.0)
    xp = x[ixp]; yp = y[iyp]; np = length(vp)
    ixn,iyn,vn = findwhere(diffx .< 0.0)
    xn = x[ixn]; yn = y[iyn]; nn = length(vn)

    if shift
        dx = x[2]-x[1]; dy = y[2] - y[1]
        xp .-= dx/2; yp .-= dy/2; xn .-= dx/2; yn .-= dy/2
    end

    vortices = [xn yn -vn; xp yp vp]

    return PointVortex(vortices)
end

function findvortices_grid(psi::Torus;shift=true)
    vort = findvortices_jumps(psi,shift=shift)
    rawvort = rawData(vort)
    vort = sortslices(rawvort,dims=1)
    return PointVortex(vort)
end

function findvortices_grid(psi::Sphere;shift=true)
    @unpack ψ = psi
    windvals = phasejumps(angle.(ψ),2)
    vort = findvortices_jumps(psi,shift=shift)
    rawvort = rawData(vort)

    # North pole winding. - sign means polar vortex co-rotating with w > 0
    w1 = sum(windvals[1,:])
    (sign(w1) != 0) && (rawvort = [rawvort; 0.0 0.0 -w1])
    # South pole winding. + sign means co-rotating with w > 0
    w2 = sum(windvals[end,:])
    (sign(w2) != 0) && (rawvort = [rawvort; pi 0.0 w2])

    vort = sortslices(vort,dims=1)
    return PointVortex(vort)
end

function findvortices_interp(psi::Field)
    vort = findvortices_grid(psi,shift=true)
    vort = removeedgevortices(vort,psi)
    #TODO: allow for interp with periodic data (here edges are stripped)

    for (j,vortex) in enumerate(vort)
        try
        vortz,psiz = corezoom(vortex,psi)
        vortz,psiz = corezoom(vortz,psiz)
        vortz,psiz = corezoom(vortz,psiz)
        vort[j] = vortz
        catch nothing
        end
        return vort
    end
end

"""
    vortices = findvortices(ψ<:Field,interp=true)

Locate vortices as `2π` phase windings around plaquettes on a cartesian spatial grid.

Requires a 2D wavefunction `ψ(x,y)` on a cartesian grid specified by vectors `x`, `y`.

`vortices` - array of vortex coordinates `xv,yv` and charges `qv`.

Each row is of the form `[xv, yv, cv]`, and the array is sorted into lexical order
according to the `xv` coordinates.
"""
function findvortices(psi::Field,interp=true)
    if interp
        return findvortices_interp(psi)
    else
        return findvortices_grid(psi)
    end
end

"""
    vortices = removeedgevortices(vort::Array{PointVortex,1},x,y,edge=1)

Strip artifact edgevortices arising from periodic phase differencing.
"""
function removeedgevortices(vort::Array{PointVortex,1},psi::Field,edge=1)
    @unpack x,y = psi; dx=x[2]-x[1]; dy=y[2]-y[1]
    keep = []
    for j = 1:length(vort)
        xi,yi,qi = rawData(vort[j])
        xedge = isapprox(xi,x[1],atol=edge*dx) || isapprox(xi,x[end],atol=edge*dx)
        yedge = isapprox(yi,y[1],atol=edge*dy) || isapprox(yi,y[end],atol=edge*dy)
        not_edge = !(xedge || yedge)
        not_edge && push!(keep,j)
    end
    return vort[keep]
end

"""
    vortz,psiz = corezoom(vortex,ψ<:Field,winhalf=2,Nz=30)

Uses local interpolation to resolve core location.
"""
function corezoom(vortex::PointVortex,psi::T,winhalf=2,Nz=30) where T<:Field
    @unpack ψ,x,y = psi
    xv,yv,qv = rawData(vortex)
    dx=x[2]-x[1]; dy=y[2]-y[1]
    ixv = isapprox.(x,xv,atol=dx) |> findlast
    iyv = isapprox.(y,yv,atol=dy) |> findlast
    ixwin = (ixv-winhalf):(ixv+winhalf-1)
    iywin = (iyv-winhalf):(iyv+winhalf-1)
    xw = x[ixwin]; yw = y[iywin]; psiw = ψ[ixwin,iywin]
    xz = LinRange(xw[1],xw[end],Nz)
    yz = LinRange(yw[1],yw[end],Nz)
    knots = (xw,yw)
    itp = interpolate(knots, psiw, Gridded(Linear()))
    psiz = itp(xz,yz)
    ψv = T(psiz,xz |> Vector,yz |> Vector)
    vortz = findvortices_grid(ψv,shift=false)
    vortz = removeedgevortices(vortz,ψv)[1]
    return vortz,ψv
end

corezoom(vortex::Array{PointVortex,1},psi::Field,winhalf=2,Nz=30) = corezoom(vortex[1],psi,2,30)

"""
    jumps = phasejumps(ϕ,dim)

Count phase jumps greater than `π` in phase `ϕ` along dimension `dim`.

See also: [`unwrap`](@ref), [`unwrap`](@ref), [`unwrap!`](@ref)
"""
function phasejumps(phase,dim=1)
    @assert (dim==1 || dim==2)
    Nx,Ny = size(phase)
    pdiff = zero(phase)

    if dim == 1
    @inbounds for j in 1:Ny
        @inbounds for i in 2:Nx
            Δ = phase[i,j] - phase[i-1,j]
            abs(Δ) > π && (pdiff[i,j] += sign(Δ))
        end
            Δ = phase[1,j] - phase[Nx,j]
            abs(Δ) > π && (pdiff[1,j] += sign(Δ))
    end

    elseif dim == 2
    @inbounds for j in 2:Ny
        @inbounds for i in 1:Nx
            Δ = phase[i,j] - phase[i,j-1]
            abs(Δ) > π && (pdiff[i,j] += sign(Δ))
        end
    end
        @inbounds for i in 1:Nx
            Δ = phase[i,1] - phase[i,Ny]
            abs(Δ) > π && (pdiff[i,1] += sign(Δ))
        end
    end
  return pdiff
end

function phasejumps!(pdiff,phase,dim=1)
    @assert (dim==1 || dim==2)
    Nx,Ny = size(phase)

    if dim == 1
    @inbounds for j in 1:Ny
        @inbounds for i in 2:Nx
            Δ = phase[i,j] - phase[i-1,j]
            abs(Δ) > π && (pdiff[i,j] += sign(Δ))
        end
            Δ = phase[1,j] - phase[Nx,j]
            abs(Δ) > π && (pdiff[1,j] += sign(Δ))
    end

    elseif dim == 2
    @inbounds for j in 2:Ny
        @inbounds for i in 1:Nx
            Δ = phase[i,j] - phase[i,j-1]
            abs(Δ) > π && (pdiff[i,j] += sign(Δ))
        end

    end
        @inbounds for i in 1:Nx
            Δ = phase[i,1] - phase[i,Ny]
            abs(Δ) > π && (pdiff[i,1] += sign(Δ))
        end
    end
end

"""
    unwrapped = unwrap(ϕ,dim=1)

Unwraps 2d array `ϕ` along dimension `dim`, acting periodically to give back array of same size as `ϕ`.

See also: [`unwrap!`](@ref)
"""
function unwrap(phase::Array{Float64,2},dim=1)
    @assert (dim==1 || dim==2)
    Nx,Ny = size(phase)
    uphase = copy(phase)

    if dim == 1
    @inbounds for j in 1:Ny
        @inbounds for i in 2:Nx
        (uphase[i,j] - uphase[i-1,j] >= π) && (uphase[i,j] -= 2π)
        (uphase[i,j] - uphase[i-1,j] <= -π) && (uphase[i,j] += 2π)
        end
        (uphase[1,j] - uphase[Nx,j] >= π) && (uphase[1,j] -= 2π)
        (uphase[1,j] - uphase[Nx,j] <= -π) && (uphase[1,j] += 2π)
    end

    elseif dim == 2
    @inbounds for j in 2:Ny
        @inbounds for i in 1:Nx
        (uphase[i,j] - uphase[i,j-1] >= π) && (uphase[i,j] -= 2π)
        (uphase[i,j] - uphase[i,j-1] <= -π) && (uphase[i,j] += 2π)
        end

    end
        @inbounds for i in 1:Nx
        (uphase[i,1] - uphase[i,Ny] >= π) && (uphase[i,1] -= 2π)
        (uphase[i,1] - uphase[i,Ny] <= -π) && (uphase[i,1] += 2π)
        end
    end

    return uphase
end

"""
    unwrap!(unwrapped,ϕ,dim)

Unwraps 2d phase array `ϕ` along dimension `dim`,
acting periodically and writing the unwrapped array `ϕ`
in place.

See also: [`unwrap`](@ref)
"""
function unwrap!(uphase::Array{Float64,2},phase::Array{Float64,2},dim=1)
    @assert (dim==1 || dim==2)
    Nx,Ny = size(phase)
    uphase .= phase

    if dim == 1
    @inbounds for j in 1:Ny
        @inbounds for i in 2:Nx
        (uphase[i,j] - uphase[i-1,j] >= π) && (uphase[i,j] -= 2π)
        (uphase[i,j] - uphase[i-1,j] <= -π) && (uphase[i,j] += 2π)
        end
        (uphase[1,j] - uphase[Nx,j] >= π) && (uphase[1,j] -= 2π)
        (uphase[1,j] - uphase[Nx,j] <= -π) && (uphase[1,j] += 2π)
    end

    elseif dim == 2
    @inbounds for j in 2:Ny
        @inbounds for i in 1:Nx
        (uphase[i,j] - uphase[i,j-1] >= π) && (uphase[i,j] -= 2π)
        (uphase[i,j] - uphase[i,j-1] <= -π) && (uphase[i,j] += 2π)
        end

    end
        @inbounds for i in 1:Nx
        (uphase[i,1] - uphase[i,Ny] >= π) && (uphase[i,1] -= 2π)
        (uphase[i,1] - uphase[i,Ny] <= -π) && (uphase[i,1] += 2π)
        end
    end
end

function unwrap(phase::Array{Float64,1})
    uphase = copy(phase)
    Nx = length(phase)
        @inbounds for i in 2:Nx
        (uphase[i] - uphase[i-1] >= π) && (uphase[i] -= 2π)
        (uphase[i] - uphase[i-1] <= -π) && (uphase[i] += 2π)
        end
        (uphase[1] - uphase[Nx] >= π) && (uphase[1] -= 2π)
        (uphase[1] - uphase[Nx] <= -π) && (uphase[1] += 2π)
        return uphase
end

function unwrap!(uphase::Array{Float64,1},phase::Array{Float64,1})
    Nx = length(phase)
        @inbounds for i in 2:Nx
        (uphase[i] - uphase[i-1] >= π) && (uphase[i] -= 2π)
        (uphase[i] - uphase[i-1] <= -π) && (uphase[i] += 2π)
        end
        (uphase[1] - uphase[Nx] >= π) && (uphase[1] -= 2π)
        (uphase[1] - uphase[Nx] <= -π) && (uphase[1] += 2π)
end


text:
function define_indices_d(::Type{NoInterp}, d, pad)
    symix, symx = Symbol("ix_",d), Symbol("x_",d)
    :($symix = convert(Int, $symx))
end

function coefficients(::Type{NoInterp}, N, d)
    :()
end

function index_gen{IT<:DimSpec}(::Type{NoInterp}, ::Type{IT}, N::Integer, offsets...)
    if (length(offsets) < N)
        return :($(index_gen(IT, N, offsets..., 0)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end

padding(::Type{NoInterp}) = Val{0}()

# How many non-NoInterp dimensions are there?
count_interp_dims(::Type{NoInterp}, N) = 0
count_interp_dims{IT<:InterpolationType}(::Type{IT}, N) = N
function count_interp_dims{IT<:Tuple{Vararg{InterpolationType}}}(it::Type{IT}, N)
    n = 0
    for p in it.parameters
        n += count_interp_dims(p, 1)
    end
    n
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using JLD2
using BenchmarkTools
using Benchmarks
using Plots
pyplot()
# Benchmark parameters

decomposition = Slab()

ranks = (1, 2, 4, 8, 16)

grid_size(R, decomposition) = (256, 256, 256)

rank_size(R, ::Slab) = (1, R, 1)
rank_size(R, ::Pencil) = Int.((1, √R, √R))

# Run benchmarks

print_system_info()

for R in ranks
    Nx, Ny, Nz = grid_size(R, decomposition)
    Rx, Ry, Rz = rank_size(R, decomposition)
    @info "Benchmarking distributed incompressible model strong scaling with $(typeof(decomposition)) decomposition [N=($Nx, $Ny, $Nz), ranks=($Rx, $Ry, $Rz)]..."
    julia = Base.julia_cmd()
    run(`mpiexec -np $R $julia --project strong_scaling_incompressible_model_single.jl $(typeof(decomposition)) $Nx $Ny $Nz $Rx $Ry $Rz`)
end

# Collect and merge benchmarks from all ranks

suite = BenchmarkGroup(["size", "ranks"])

for R in ranks
    Nx, Ny, Nz = grid_size(R, decomposition)
    Rx, Ry, Rz = rank_size(R, decomposition)
    case = ((Nx, Ny, Nz), (Rx, Ry, Rz))

    for local_rank in 0:R-1
        filename = string("strong_scaling_incompressible_model_$(R)ranks_$(typeof(decomposition))_$local_rank.jld2")
        jldopen(filename, "r") do file
            if local_rank == 0
                suite[case] = file["trial"]
            else
                merged_trial = suite[case]
                local_trial = file["trial"]
                append!(merged_trial.times, local_trial.times)
                append!(merged_trial.gctimes, local_trial.gctimes)
            end
        end
    end
end

# Summarize benchmarks

plot_keys = collect(keys(suite))
sort!(plot_keys, by = v -> v[2][2])
plot_num = length(plot_keys)
rank_num = zeros(Int64, plot_num)
run_times = zeros(Float64, plot_num)
eff_ratio = zeros(Float64, plot_num)
for i in 1:plot_num
    rank_num[i] = plot_keys[i][2][2]
    run_times[i] = mean(suite[plot_keys[i]].times) / 1.0e6
    eff_ratio[i] = median(suite[plot_keys[1]].times) /(rank_num[i] * median(suite[plot_keys[i]].times))
end

plt = plot(rank_num, run_times, lw=4, xaxis=:log2, legend=:none,
          xlabel="Cores", ylabel="Times (ms)", title="Strong Scaling Incompressible Times")
display(plt)
savefig(plt, "ss_incompressible_times.png")

plt2 = plot(rank_num, eff_ratio, lw=4, xaxis=:log2, legend=:none, ylims=(0,1.1),
            xlabel="Cores", ylabel="Efficiency", title="Strong Scaling Incompressible Efficiency")
display(plt2)
savefig(plt2, "ss_incompressible_efficiency.png")

df = benchmarks_dataframe(suite)
sort!(df, :ranks)
benchmarks_pretty_table(df, title="Incompressible model strong scaling benchmark")

base_case = (grid_size(1, decomposition), rank_size(1, decomposition))
suite_Δ = speedups_suite(suite, base_case=base_case)
df_Δ = speedups_dataframe(suite_Δ, efficiency=:strong, base_case=base_case, key2rank=k->prod(k[2]))
sort!(df_Δ, :ranks)
benchmarks_pretty_table(df_Δ, title="Incompressible model strong scaling speedup")


text:
using Gurobi: getcoeff

@testset "File Tests" begin
    @testset for file in filter(
            s -> endswith(s, ".jl"), readdir(joinpath(@__DIR__, "C")))
        evalfile(joinpath(@__DIR__, "C", file))
    end
end

@testset "Empty constraints (Issue #142)" begin
    @testset "No variables, no constraints" begin
        model = Gurobi.Model(Gurobi.Env(), "model")
        A = Gurobi.get_constrmatrix(model)
        @test size(A) == (0, 0)
    end
    @testset "One variable, no constraints" begin
        model = Gurobi.Model(Gurobi.Env(), "model")
        Gurobi.add_cvar!(model, 0.0)
        Gurobi.update_model!(model)
        A = Gurobi.get_constrmatrix(model)
        @test size(A) == (0, 1)
    end
end

@testset "changing coefficients" begin
    A = [1.0 2.0; 3.0 0.0]
    b = [0.1, 0.2]
    f = [0.0, 0.0]
    model = gurobi_model(Gurobi.Env(); f=f, A=A, b=b)

    @test getcoeff(model, 1, 1) == 1.0

    # Verify that we can pass any kind of `::Integer` to `getcoeff`
    @test getcoeff(model, 1, Int8(1)) == 1.0
    @test getcoeff(model, UInt32(1), UInt32(1)) == 1.0

    @test getcoeff(model, 1, 2) == 2.0
    @test getcoeff(model, 2, 1) == 3.0
    @test getcoeff(model, 2, 2) == 0.0

    chg_coeffs!(model, 2, 1, 1.5)
    # before updating the model, we still get the old coefficient value
    @test getcoeff(model, 2, 1) == 3.0
    # after updating the model, we see the new coefficient value
    update_model!(model)
    @test getcoeff(model, 2, 1) == 1.5

    # Verify that we are automatically converting the scalars to
    # Cint and Float64 as necessary
    chg_coeffs!(model, Int8(2), Int128(1), Float32(2.0))
    update_model!(model)
    @test getcoeff(model, 2, 1) == 2.0

    chg_coeffs!(model, [1, 2], [1, 1], [5.0, 6.5])
    update_model!(model)
    @test getcoeff(model, 1, 1) == 5.0
    @test getcoeff(model, 2, 1) == 6.5
    @test getcoeff(model, 1, 2) == 2.0
    @test getcoeff(model, 2, 2) == 0.0

    chg_coeffs!(model, Int8[1, 2], Int8[1, 1], Float32[1.0, 2.0])
    update_model!(model)
    @test getcoeff(model, 1, 1) == 1.0
    @test getcoeff(model, 2, 1) == 2.0
    @test getcoeff(model, 1, 2) == 2.0
    @test getcoeff(model, 2, 2) == 0.0
end

@testset "Unicode support for names" begin
    model = Gurobi.Model(Gurobi.Env(), "model")
    Gurobi.add_cvar!(model, 0.0)
    for str in [
        "αβ", "Aα", "αA", "x̂", "xx̂", "x̂x", "¹", "x¹", "¹x", "₂", "x₂", "₂x",
        "xß", "̂", "̄"
    ]
        Gurobi.set_strattr!(model, "ModelName", str)
        Gurobi.update_model!(model)
        @test str == Gurobi.get_strattr(model, "ModelName")

        Gurobi.set_strattrelement!(model, "VarName", 1, str)
        Gurobi.update_model!(model)
        @test str == Gurobi.get_strattrelement(model, "VarName", 1)
    end
end


text:
import Base.@kwdef

struct BruteForce end

@kwdef struct TreecodeStretch{T}
    n::Int
    N0::Int
    eta::T
end

@kwdef struct TreecodeAvgRestFrame{T}
    n::Int
    N0::Int
    eta::T
end

@kwdef struct TreecodeUniform{T}
    n::Int
    N0::Int
    eta::T
end


text:
# Bogumił Kamiński, 2019-03-25

using DataFrames

function parselog(lines)
    author = r"^Author: ([^<]*) <"
    insc = r"^.+changed, ([0-9]+) insertion"
    delc = r"^.+changed.*, ([0-9]+) deletion"
    authordata = DataFrame(author=String[], action=String[], count=Int[])
    curauthor = ""
    for line in lines
        m = match(author, line)
        m === nothing || (curauthor = m[1])
        m = match(insc, line)
        m === nothing || push!(authordata,
                               [curauthor, "insertion", parse(Int, m[1])])
        m = match(delc, line)
        m === nothing || push!(authordata,
                               [curauthor, "deletion", parse(Int, m[1])])
    end
    authorstats = by(authordata, [:author, :action], count = :count=>sum)
    unstack(authorstats, :action, :count)
end

function gitstats(dir)
    if isdir(dir)
        println("\nAnalyzing")
        cd(dir) do
            try
                res = read(`git log --stat`, String)
                lines = split(res, ['\r', '\n'], keepempty=false)
                df = parselog(lines)
                df.all = coalesce.(df.deletion, 0) .+
                         coalesce.(df.insertion, 0)
                display(sort!(df, :all, rev=true))
            catch
                error("Running git log failed")
            end
        end
    else
        error("$dir is not a directory")
    end
end

max_stars_repo_path: 5. Logparser/logparser_2.jl

max_issues_repo_path: 5. Logparser/logparser_2.jl

max_forks_repo_path: 5. Logparser/logparser_2.jl


text:
@testset "213.house-robber-ii.jl" begin
    @test rob_213([200, 3, 140, 20, 10]) == 340
    @test rob_213([2, 3, 2]) == 3
    @test rob_213([1, 2, 3, 1]) == 4
    @test rob_213([1, 2, 3]) == 3
    @test rob_213([3]) == 3
    @test rob_213([1, 2]) == 2
end


text:
module Defs

using Comm

# Global definitions, aka parameter file

const run_type = :debug     # :debug, :benchmark, :production

const ncells_per_grid = 10

const xmin = 0.0
const xmax = 1.0
const cfl = 0.5
const tmin = 0.0
const tmax = 1.0

const wait_every = 0

if run_type == :debug
    const rho = 1               # resolution scale
    const nsteps = -1
    const info_every = 10
    const file_every = 0
elseif run_type == :benchmark
    const rho = 100             # resolution scale
    const nsteps = 10
    const info_every = 0
    const file_every = -1
elseif run_type == :production
    @assert false
else
    @assert false
end

ncells() = (rho * ncells_per_grid * Comm.nprocs())::Int
dx() = ((xmax - xmin) / ncells())::Float64
dt() = (cfl * dx())::Float64

function do_this_time(iteration::Integer, do_every::Integer;
                      is_first::Bool=false, was_last::Bool=false)
    if do_every < 0 return false end
    if is_first return true end
    should_output = do_every > 0 && iteration % do_every == 0
    return should_output != was_last
end

end


text:
using RedPitayaDAQServer
using PyPlot

#rp = RedPitaya("rp-f04972.local")
rp = RedPitaya("192.168.178.46")

dec = 32
modulus = 4800
base_frequency = 125000000
samples_per_period = div(modulus, dec)
periods_per_frame = 10

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)

modeDAC(rp, "STANDARD")

frequencyDAC(rp, 1, 1, base_frequency / modulus)
println(" frequency = $(frequencyDAC(rp,1,1))")

signalTypeDAC(rp, 1 , "SINE")

amplitudeDAC(rp, 1, 1, 2000)
offsetDAC(rp, 1, 0)
phaseDAC(rp, 1, 1, 0.0 )

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "INTERNAL")
masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)
uFirstPeriod = readData(rp, 0, 1)
sleep(0.1)
fr = currentFrame(rp)
@show fr

uCurrentPeriod = readData(rp, fr, 1)

figure(1)
clf()
plot(vec(uFirstPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,1,:,:]))
legend(("first period", "current period"))


text:
#!/usr/bin/env julia

using PhaseSpaceTools, Test

function randuv()
    u,v = randnc(2)
    nrm = abs2(u)-abs2(v)
    if nrm < 0
        new = [0 1;1 0]*[u; v]
        u,v = new
    end
    u /=sqrt(abs(nrm))
    v /=sqrt(abs(nrm))
    return u,v
end

@testset "Arg tests" begin

    # Coherent
    state = Coherent(10.0)
    @test typeof(state) <: State
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(glauberP(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Fock
    n = 100
    state = Fock(n)
    @test typeof(state) == typeof(Fock(n))
    @test typeof(state) <: State
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Bogoliubov
    state = Bogoliubov(.1,.1,1)
    @test typeof(state) <: State
    @test typeof(state) == typeof(Bogoliubov(.1,.1,1))
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(glauberP(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Squeezed
    state = Squeezed(randnc(),randnc())
    @test typeof(state) <: State
    @test typeof(state) == typeof(Squeezed(randnc(),randnc()))
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Thermal
    state = Thermal(.1,1)
    @test typeof(state) <: State
    @test typeof(state) == typeof(Thermal(.1,1))
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(glauberP(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Crescent
    state = Crescent(1.0,0.1,0.1)
    @test typeof(state) <: State
    @test typeof(state) == typeof(Crescent(1.0,0.1,0.1))
    @test typeof(wigner(state,1)[1][1]) == Complex{Float64}
    @test typeof(husimiQ(state,1)[1][1]) == Complex{Float64}
    @test typeof(positiveP(state,1)[1][1]) == Complex{Float64}

    # Helpers
    f(x) = exp(-x^2)
    @test typeof(reject(f,[-1,1],1,1.1)[1]) == Float64
    @test typeof(plaguerre(.1,1)) == Float64
    @test typeof(plaguerre_asymptotic(.1,1)) == Float64
    @test typeof(laguerren(.1,1)) == Float64
    @test typeof(plaguerre_asymptotic(.1,1)) == Float64

end

@testset "Noises" begin

    N = 100_000
    a = randnc(N)
    mean(abs2.(a))

    @test isapprox(mean(abs2.(a)),1.0,rtol=5e-2)
    @test isapprox(abs(mean(a)),0.0,atol=1e-2)
end

@testset "Coherent W" begin

    # Coherent
    N = 100000
    α = 100
    state = Coherent(α)
    a,ā = wigner(state,N)

    meana = mean(a)
    n̄ = real(mean(@. ā*a)-.5)
    Vn = mean(@. ā^2*a^2)-mean(@. a*ā)-n̄^2 |> abs

    #test
    @test norm(meana - α)/abs(α) < 0.01
    @test abs(n̄ - abs(α).^2)/abs(α)^2 < 0.01
    @test abs(Vn - abs(α)^2)/abs(α)^2 < 0.05

end

@testset "Coherent +P" begin

    N = 1000
    α = 100*randnc()
    state = Coherent(α)
    a,ā = positiveP(state,N)

    meana = mean(a)
    n̄ = real(mean(@. ā*a))
    Vn = mean(@. ā^2*a^2)-mean(@. a*ā)-n̄^2 |> abs

    #test
    @test n̄ ≈ abs2(α)
    @test norm(meana - α)/abs(α) < 0.01
    @test abs(n̄ - abs(α).^2)/abs(α)^2 < 0.01
    @test abs(Vn - abs(α)^2)/abs(α)^2 < 0.05

end

@testset "Thermal P,Q,W" begin

    N = 100001
    n = 50
    β = 0.
    s = Thermal(β,n)

    a,ā = glauberP(s,N)
    na = real(mean(a.*ā))
    @test isapprox(na,n,atol=.9)

    a,ā = wigner(s,N)
    na = real(mean(a.*ā))-0.5
    @test isapprox(na,n,atol=.9)

    a,ā = husimiQ(s,N)
    na = real(mean(a.*ā))-1.0
    @test isapprox(na,n,atol=.9)

end

@testset "Fock W" begin

    # Fock
    n = 100
    N = 100000
    state = Fock(n)
    a,ā = wigner(state,N)

    meana = mean(a)
    absa = abs(meana)
    n̄ = mean(abs2.(a))-.5
    Vn= abs(mean(abs2.(a).^2)-mean(abs2.(a))-n̄.^2)
    rel_num_var = sqrt(abs(Vn))/abs(n̄);

    #test
    @test absa < 0.1
    @test n̄ - n < 0.01
    @test Vn < 0.01
    @test rel_num_var < 0.001

end

@testset "Fock +W" begin

#test +W

    n = 99
    N = 100_000
    state = Fock(n)
    a,ā = positiveW(state,N)

    meana = mean(a)
    absa = abs(meana)
    n̄ = real(mean(a.*ā)-.5)
    Vn = abs(mean(a.*a.*ā.*ā)-abs(mean(a.*ā))-n̄.^2)
    rel_num_var = sqrt(abs(Vn))/abs(n̄);

    #test
    @test absa < 0.1
    @test n̄ - n < 1
    @test Vn < 50
    @test rel_num_var < 0.1

end

@testset "Fock +W (n≫1)" begin

#test +W for large n

    n = 101
    state = Fock(n)
    N = 100000

    a,ā = positiveW(state,N)

    av_a = mean(a)
    absa = abs(av_a)
    n̄ = real(mean(a.*ā)-.5)
    Vn= abs(mean(a.*a.*ā.*ā)-mean(a.*ā)-n̄.^2)
    rel_num_var = sqrt(abs(Vn))/abs(n̄)

    #test
    @test absa < 0.1
    @test (n̄ - n)/n < 0.01
    @test Vn < 30
    @test rel_num_var < 0.1

end

@testset "Fock +P" begin

    n = 99
    N = 100_000
    state = Fock(n)
    a,ā = positiveP(state,N)

    meana = mean(a)
    absa = abs(meana)
    n̄ = mean(a.*ā) |> real
    Vn = abs(mean(a.*a.*ā.*ā))-abs(mean(a.*ā)).^2
    rel_num_var = sqrt(abs(Vn))/abs(n̄);

    #test
    @test absa < 0.1
    @test n̄ - n < 1
    @test Vn < 50
    @test rel_num_var < 0.15

end

@testset "Squeezed +P " begin

    #Squeezed
    β = 10
    ϕ = π/16
    r = 2
    ϵ = r*exp(2*im*ϕ)
    state = Squeezed(β,ϵ)
    N = 10000

    a,ā = positiveP(state,N)

    av_a = mean(a);absa = abs(av_a)
    n̄ = real(mean(a.*ā))
    nbar = sinh(abs(ϵ)).^2 + abs2(β)

    #test
    @test norm(av_a - β)/abs(β) < 0.02
    @test abs(n̄ - nbar)/abs(β)^2 < 0.02

end

@testset "Squeezed W" begin

    # Squeezed
    β = 10
    ϕ = π/16
    r = 1.5
    ϵ = r*exp(2*im*ϕ)
    state = Squeezed(β,ϵ)
    N = 10000
    a,ā = wigner(state,N)

    av_a = mean(a);absa = abs(av_a)
    n̄ = real(mean(a.*ā)-.5)
    nbar = sinh(abs(ϵ)).^2+abs2(β)

    #test
    @test norm(av_a - β)/abs(β) < 0.01
    @test abs(n̄ - nbar)/abs(β)^2 < 0.01

end

@testset "Bogoliubov W" begin

    # Bogoliubov
    N = 100000
    n̄ = 10

    u,v = randuv()
    abs2(u)-abs2(v) ≈ 1.0

    #thermal state
    state = Bogoliubov(u,v,n̄)
    a,ā = wigner(state,N)

    # particle mode population
    na = real(mean(a.*ā)) - 0.5
    nth = (abs2(u) + abs2(v))*(n̄ + 0.5) - 0.5
    @test isapprox(na,nth,rtol=1e-1)

    #test vacuum limit
    state = Bogoliubov(u,v,0)
    a,ā = wigner(state,N)

    na = real(mean(a.*ā)) - 0.5
    nvac = abs2(v)
    @test isapprox(na,nvac,atol=0.3)

end

@testset "Bogoliubov Q" begin

    # Bogoliubov
    N = 100000
    n̄ = 10

    u,v = randuv()
    abs2(u)-abs2(v) ≈ 1.0

    #thermal state
    state = Bogoliubov(u,v,n̄)
    a,ā = husimiQ(state,N)

    # particle mode population
    na = real(mean(a.*ā)) - 1.0
    nth = (abs2(u) + abs2(v))*(n̄ + 1.0) - 1.0
    @test isapprox(na,nth,rtol=1e-1)

    #test vacuum limit
    state = Bogoliubov(u,v,0)
    a,ā = husimiQ(state,N)

    na = real(mean(a.*ā)) - 1.0
    nvac = abs2(v)
    @test isapprox(na,nvac,atol=0.3)

end

@testset "Bogoliubov +P" begin

    # Bogoliubov
    N = 100000
    n̄ = 10

    u,v = randuv()
    @test abs2(u)-abs2(v) ≈ 1.0

    #thermal state
    state = Bogoliubov(u,v,n̄)
    a,ā = positiveP(state,N)

    # particle mode population
    na = mean(a.*ā) |> real
    nth = (abs2(u) + abs2(v))*n̄ + abs(v)
    @test isapprox(na,nth,rtol=1e-1)

    #test vacuum limit
    state = Bogoliubov(u,v,0)
    a,ā = positiveP(state,N)

    na = mean(a.*ā) |> real
    nvac = abs2(v)
    @test isapprox(na,nvac,atol=0.3)

end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e8e2c9d5-15a6-4a3a-81b7-a10c32290f79
# environment setup
begin
	using DrWatson
	quickactivate(findproject())

	using Pkg
	Pkg.instantiate()

	using Colors, ColorVectorSpace, ImageShow, FileIO
	using ImageShow.ImageCore
	using ColorSchemes
	using InteractiveUtils, PlutoUI
	using Plots
	using Distributions
	using LinearAlgebra

end

# ╔═╡ c797f938-ab2a-4c0b-99a1-dc1d8a975193
include(srcdir("IFS.jl"));

# ╔═╡ 03386e5e-937b-4e32-b034-e1689834e0dc
TableOfContents()

# ╔═╡ fe69d7ee-aa66-11eb-0e4e-e16fc28be586
md"# Fractal Drawing"

# ╔═╡ 1f0068b8-4df2-44dd-a3a3-072d3c43d66a
md"## Probabilistic method"

# ╔═╡ df738342-cb99-4a8e-9565-2b9422caaee2
md"### Basic usage"

# ╔═╡ cffe70b9-17bc-467a-b2c1-7a01e46ce494
begin
	# set the parameters
	linearIFS = [[0  0; 0 0.16],
				 [0.85  0.04; -0.04 0.85 ],
				 [0.2  -0.26; 0.23 0.22],
				 [-0.15  0.28; 0.26 0.24]]
	transIFS = [[0, 0],
				[0, 1.6],
				[0, 1.6],
				[0, 0.44]]
	weights = [0.01, 0.84, 0.08, 0.07]
	# create the weighted IFS
	myWIFS = WIFS(linearIFS, transIFS, weights)
end

# ╔═╡ 87b32b70-a405-42db-87bd-16a44644e334
md"### Gallery"

# ╔═╡ 24f874d6-93f6-4f53-941c-4cc3cf56963b
md"""
numPts: $(@bind maxNumPts Slider(1000:1000:300000; default=100000, show_value=true)) $ \quad $ Example: $(@bind selectIFS MultiSelect([ "SierpinskiTriangle" => "Sierpinski Triangle",
							"BedMcCarpet" => "Bedford-MacMullen Carpet",
							"BaranskiCarpet" => "Baranski Carpet",
							"HeighwayDragon" => "Heighway Dragon",
							"Twindragon" => "Twindradon",
							"Terdragon" => "Terdragon",
 							"BarnsleyFern" => "Barnsley Fern",
							"myIFSNonlinear" => "nonlinear"];
							default=["SierpinskiTriangle"])) $ \quad $ Color: $(@bind myColor ColorStringPicker(default="#009AFA"))
"""

# ╔═╡ 5f8abed9-0be7-43c7-a6a5-8f93a94eb5ce
md"Bed-MacMullen carpet setup:"

# ╔═╡ e63a2bee-5b49-4b73-8598-fcf8c17b478a
posBM = [1 0 1 1;
		1 0 1 0;
		0 1 1 0;]

# ╔═╡ 106cbdcb-8041-40e4-924a-4693c8d2ea10
md"Baranski carpet setup:"

# ╔═╡ e46ee341-23bb-45a0-8e00-58afdd825634
 vv = [0.6, 0.1, 0.3]

# ╔═╡ 0ba7d829-0a5d-4891-9b50-3203891faa22
hh = [0.5, 0.1, 0.1 , 0.3]

# ╔═╡ 2ca1aa8e-4450-43a2-a15a-ba5d5e2a652f
posB = [1 1 0 1;
	   1 0 0 1;
	   0 1 1 1]

# ╔═╡ d3e449f3-6318-4d9d-87ee-12d0733b65c9
begin
	f(x) = [0.5*x[1]^2 + x[2], x[2] ^ 3]
	g(x) = [0.2* x[1] * x[2] + 1, -0.3 * x[1] ]
	myIFSNonlinear = IFSNonlinear([f, g], 2)
end;

# ╔═╡ dc9d671f-5279-4c8f-a4ef-9ee31ce49e76
begin
	if selectIFS[1] == "myIFSNonlinear"
		ptsSet = itrPtsProb(myIFSNonlinear, maxNumPts)
	else
		itrIFS = eval(Meta.parse("PredefinedIFS." * selectIFS[1]))
		if selectIFS[1] in ["SierpinskiTriangle", "BarnsleyFern", "HeighwayDragon", "Twindragon", "Terdragon"]
			ptsSet = itrPtsProb(itrIFS, maxNumPts)
		elseif selectIFS[1] == "BedMcCarpet"
			ptsSet = itrPtsProb(itrIFS(posBM), maxNumPts)
		else
			ptsSet = itrPtsProb(itrIFS(vv, hh, posB), maxNumPts)
		end
	end

	println("")
end

# ╔═╡ aad51fe8-0ba2-4a68-9860-3d6f55dce7d6
md"## Deterministic method"

# ╔═╡ 99984492-63dd-423f-9db7-1c65370bfcfb
md"""
!!! note
	The method has been implemented via matlab script [`PlanarAffineIFS.m`](https://github.com/zfengg/FractalDraw/blob/master/scripts/matlab/2DIFSPlot/PlanarAffineIFS.m).
"""

# ╔═╡ 68185989-796e-4182-a779-3581b3ce42d1
md"## Lindenmayer System"

# ╔═╡ e0606efb-fdab-4812-8e8f-b64142dc3a56
md"""
!!! note
	The first implementation is based upon [`Lindenmayer.jl`](https://github.com/cormullion/Lindenmayer.jl).
"""

# ╔═╡ 3a84c3cd-e519-4398-89f8-a49c608fdb98
md"### Basic usage"

# ╔═╡ b63e0adb-7125-4fdd-a654-87963bf13282

# ╔═╡ f96a33a5-e70f-4807-8974-d4362ef424db
md"### Gallery"

# ╔═╡ 3ba4d982-b95b-40e0-9844-b4801b7a1e89
md"## Appendix"

# ╔═╡ 7afa357b-4483-4cff-8b12-9401d42ea53b
begin
function plot_ptsSet(ptsSet::Vector{Vector{Float64}}, mc::String="#009AFA")
	xPtsSet = [x[1] for x in ptsSet]
	yPtsSet = [x[2] for x in ptsSet]
	scatter(xPtsSet, yPtsSet,
			leg=false,
			markershape=:circle,
			markeralpha=0.7;
			markersize=1,
			markercolor=mc,
			# markerstroke=false,
    		markerstrokewidth = 0,
    		# markerstrokealpha = 0.2,
			# grid=false,
			showaxis=false,
			ticks=false,
			xlims=extrema(xPtsSet),
			ylims=extrema(yPtsSet),
			size=(680, 500))
end
# function plot_ptsSet(ptsSet::Vector{Vector{Float64}}, makercolor::String)
# 	xPtsSet = [x[1] for x in ptsSet]
# 	yPtsSet = [x[2] for x in ptsSet]
# 	scatter(xPtsSet, yPtsSet,
# 			leg=false,
# 			markershape=:circle,
# 			markeralpha=0.7;
# 			markersize=1,
# 			# markercolor = :purple,
# 			# markerstroke=false,
#     		markerstrokewidth = 0,
#     		# markerstrokealpha = 0.2,
# 			# grid=false,
# 			showaxis=false,
# 			ticks=false,
# 			xlims=extrema(xPtsSet),
# 			ylims=extrema(yPtsSet),
# 			size=(680, 500))
# end

end

# ╔═╡ 5c10f0b0-2275-459d-a6a2-357c8dfb1af1
# generate and plot the points
plot_ptsSet(itrPtsProb(myWIFS, 300000))

# ╔═╡ 5a59bd1c-71e2-44dd-a143-2b3594ddc738
plot_ptsSet(ptsSet, myColor)

# ╔═╡ ac80253a-3422-4457-a51b-6baa5b48ce4f
begin
	show_image(M) = get.([ColorSchemes.rainbow], M ./ maximum(M))
	show_image(x::AbstractVector) = show_image(x')
end

# ╔═╡ fbfb70b3-95a8-45cd-9111-648de10e683c
show_image(posBM)

# ╔═╡ 466de16d-7eeb-427b-b81e-eddb9b93efb3
show_image( (vv * reshape(hh, 1, :)) .* posB)

# ╔═╡ 210c8db0-5083-470c-811b-1aca62ee24b2
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ d726d487-35bc-4233-b613-93e4c4af2984
# ingredients(srcdir("IFS.jl"))

# ╔═╡ 177eb41d-4c00-492d-8657-c8c9d65b9bd8
# @bind selectIFS MultiSelect([ "SierpinskiTriangle" => "Sierpinski Gasket",
# 							"BedMcCarpet" => "Bedford-MacMullen Carpet",
# 							"BaranskiCarpet" => "Baranski Carpet",
# 							"HeighwayDragon" => "Heighway Dragon",
# 							"Twindragon" => "Twindradon",
# 							"Terdragon" => "Terdragon",
#  							"BarnsleyFern" => "Barnsley Fern",
# 							"myIFSNonlinear" => "nonlinear"];
# 							default=["SierpinskiTriangle"])

# ╔═╡ Cell order:
# ╟─e8e2c9d5-15a6-4a3a-81b7-a10c32290f79
# ╟─c797f938-ab2a-4c0b-99a1-dc1d8a975193
# ╟─03386e5e-937b-4e32-b034-e1689834e0dc
# ╟─fe69d7ee-aa66-11eb-0e4e-e16fc28be586
# ╟─1f0068b8-4df2-44dd-a3a3-072d3c43d66a
# ╟─df738342-cb99-4a8e-9565-2b9422caaee2
# ╠═cffe70b9-17bc-467a-b2c1-7a01e46ce494
# ╟─5c10f0b0-2275-459d-a6a2-357c8dfb1af1
# ╟─87b32b70-a405-42db-87bd-16a44644e334
# ╟─24f874d6-93f6-4f53-941c-4cc3cf56963b
# ╟─5a59bd1c-71e2-44dd-a143-2b3594ddc738
# ╟─dc9d671f-5279-4c8f-a4ef-9ee31ce49e76
# ╟─5f8abed9-0be7-43c7-a6a5-8f93a94eb5ce
# ╠═fbfb70b3-95a8-45cd-9111-648de10e683c
# ╠═e63a2bee-5b49-4b73-8598-fcf8c17b478a
# ╟─106cbdcb-8041-40e4-924a-4693c8d2ea10
# ╟─466de16d-7eeb-427b-b81e-eddb9b93efb3
# ╟─e46ee341-23bb-45a0-8e00-58afdd825634
# ╟─0ba7d829-0a5d-4891-9b50-3203891faa22
# ╟─2ca1aa8e-4450-43a2-a15a-ba5d5e2a652f
# ╟─d3e449f3-6318-4d9d-87ee-12d0733b65c9
# ╟─aad51fe8-0ba2-4a68-9860-3d6f55dce7d6
# ╟─99984492-63dd-423f-9db7-1c65370bfcfb
# ╟─68185989-796e-4182-a779-3581b3ce42d1
# ╟─e0606efb-fdab-4812-8e8f-b64142dc3a56
# ╟─3a84c3cd-e519-4398-89f8-a49c608fdb98
# ╠═b63e0adb-7125-4fdd-a654-87963bf13282
# ╟─f96a33a5-e70f-4807-8974-d4362ef424db
# ╟─3ba4d982-b95b-40e0-9844-b4801b7a1e89
# ╟─7afa357b-4483-4cff-8b12-9401d42ea53b
# ╟─ac80253a-3422-4457-a51b-6baa5b48ce4f
# ╟─210c8db0-5083-470c-811b-1aca62ee24b2
# ╟─d726d487-35bc-4233-b613-93e4c4af2984
# ╟─177eb41d-4c00-492d-8657-c8c9d65b9bd8


text:
function solve_069()
    phi = zeros(1000000)
    for i = 1:1000000
        phi[i] = totient(i)
    end
    range = 1:1000000
    fun = range ./ phi
    _, id = findmax(fun)
    return id
end


text:
module TextGrids

using StringEncodings

export
    # types
    Interval, Tier, TextGrid,

    # readwrite
    read_TextGrid, write_TextGrid,

    # tiers
    rename!, relabel!,
    insert_tier!, insert_interval_tier!, insert_point_tier!,
    extract_tier,
    remove_tier!, remove_empty_tiers!,
    is_empty,

    # query
    find_interval,
    find_low_interval,
    find_high_interval,
    find_interval_by_edge,

    find_low_point,
    find_high_point,
    find_nearest_point,

    # boundaries
    insert_boundary!,
    insert_interval!,
    insert_intervals!,

    remove_boundary!,
    remove_left_boundary!,
    remove_right_boundary!,
    remove_interval!,

    move_left_boundary!,
    move_right_boundary!,
    move_boundary!,
    move_interval!,

    # points
    insert_point!,
    insert_points!,

    remove_point!,
    remove_points!

include("types.jl")
include("readwrite.jl")
include("tiers.jl")
include("query.jl")
include("boundaries.jl")
include("points.jl")
include("utils.jl")

end


text:
## Related to constraint handlers

"""
Recover Julia object which is attached to user constraint (via constraint data).
"""
function user_constraint(cons_::Ptr{SCIP_CONS})
    @assert cons_ != C_NULL
    consdata_::Ptr{SCIP.SCIP_CONSDATA} = SCIPconsGetData(cons_)
    return unsafe_pointer_to_objref(consdata_)
end

"""
Extract solution values for given variables.
"""
function sol_values(o::Optimizer, vars::AbstractArray{VI}, sol::Ptr{SCIP_SOL}=C_NULL)
    return [SCIPgetSolVal(o, sol, var(o, vi)) for vi in vars]
end


text:
module TestLocalProperties

using Econometrics
using Base.Test

## include("/home/chris/research/julia/Econometrics/src/localProperties.jl")

## test getInd2RangeFunc
##----------------------

edges = [5:1:12.]
ind2RangeFunc = Econometrics.getInd2RangeFunc(edges)

@test ind2RangeFunc(0) == edges[1]
@test ind2RangeFunc(6) == edges[end-1]

end


text:
function _ermak_loop!(
    pos::AbstractArray,
    forces::AbstractArray,
    rnd_matrix::AbstractArray,
    ttime::AbstractFloat,
    boxl::AbstractFloat;
    pbc::Bool = true,
)
    for j  = axes(pos, 2)
        @inbounds @fastmath for i = axes(pos, 1)
            pos[i, j] += (forces[i, j] * ttime) + rnd_matrix[i, j]
            if pbc
                pos[i, j] -= boxl * round(pos[i, j] / boxl)
            end
        end
    end
end

function ermak!(
    positions::AbstractArray,
    forces::AbstractArray,
    ttime::AbstractFloat,
    boxl::AbstractFloat,
    rnd_matrix::AbstractArray;
    pbc::Bool = true,
)
    𝝈 = √(2.0 * ttime)
    𝝈 = oftype(ttime, 𝝈)
    rnd_matrix .*= 𝝈
    _ermak_loop!(positions, forces, rnd_matrix, ttime, boxl; pbc = pbc)
end


text:
using Distributed, Test, CUDA

println("Envronment variables for testing")
println(ENV)

@testset "AdvancedHMC" begin
    tests = [
        "metric",
        "hamiltonian",
        "integrator",
        "trajectory",
        "adaptation",
        "sampler",
        "sampler-vec",
        "demo",
        "models",
    ]

    if CUDA.functional()
        @eval module TestCUDA
            include("cuda.jl")
        end
    else
        @warn "Skipping GPU tests because no GPU available."
    end

    res = map(tests) do t
        @eval module $(Symbol("Test_", t))
            include($t * ".jl")
        end
        return
    end
end


text:
using Mass2
using ReferenceMicrocalFiles

# type GetPulsesStep{T} <: AbstractStep
# 	pulse_source::T
# 	outputs::Vector{Symbol}
# 	previous_pulse_index::Int # used to keep track of location in an ljh file for example
# 	max_pulses_per_step::Int
# end
# inputs(s::GetPulsesStep) = Symbol[]
# function dostep(s::GetPulsesStep{LJHGroup},c::Channel)
# 	r = s.previous_pulse_index+1:min(s.previous_pulse_index+s.max_pulses_per_step, length(s.pulse_source))
# 	length(r)==0 && (return r)
# 	pulses, rowstamps = collect(s.pulse_source[r])
# 	pulses_out, rowstamps_out = perpulse_outputs(s,c)
# 	append!(pulses_out, pulses)
# 	assert(length(pulses_out)==last(r))
# 	append!(rowstamps_out, rowstamps)
# 	assert(length(rowstamps_out)==last(r))
# 	s.previous_pulse_index=last(r)
# 	r
# end
# graphlabel(s::GetPulsesStep) = repr(typeof(s))

# type LJHPulseEmitter
# 	ljh_group::LJHGroup
# 	max_pulses_per_step::Int
# 	donethru::Int
# end
# type LJHPulseEmitterSlice
# 	ljhpulsemitter::LJHPulseEmitter
# 	r::UnitRange{Int64}
# end
# LJHPulseEmitter(filename,max_pulses_per_step=1000) = LJHPulseEmitter(LJHGroup(filename), max_pulses_per_step, max_pulses_per_step)
# donethru(x::LJHPulseEmitter) = x.donethru
# Base.getindex(x::LJHPulseEmitter, r::UnitRange{Int64}) = LJHPulseEmitterSlice(x, r)
# function getpulses!(ljhpulsemitterslight::LJHPulseEmitterSlice, max_pulses_per_step::Int)
# 	pulses = Array(Vector{UInt16},length(ljh_group_slice))
# 	rowstamps = Array(UInt64, length(ljh_group_slice))
# 	ljh_group_slice = ljhpulsemitterslice.ljhpulsemitter.ljh_group[r]
# 	ljhpulsemitterslice.ljhpulsemitter.donethru = last(ljhpulsemitter.r)
# 	for (i,(pulse, rowstamp)) in enumerate(ljh_group_slice)
# 		pulses[i]=pulse
# 		rowstamps[i]=rowstamp
# 	end
# 	pulses, rowstamps
# end

ljhpaths = ReferenceMicrocalFiles.getljhpaths()
ljh = LJHGroup(ljhpaths[1])

steps = AbstractStep[]
# generate a spectra with multiple peaks for calibration testing
dist_energies = [4000,4700,5500,6000,6800,7500,8200,9600,10000]
dist_filt_values = 5*dist_energies.^0.8
dist_use_inds = [1,3,6,7,8,9]
dist_use_filt_values = dist_filt_values[dist_use_inds]
dist_use_energies = dist_energies[dist_use_inds]
filt_value_distribution = MixtureModel([Normal(mu,1) for mu in dist_filt_values])
push!(steps, GetPulsesStep(ljh, [:pulse, :rowstamp], 0,500))
s=steps[1]

push!(steps, PerPulseStep(compute_summary, [:pulse, :pre_samples, :frame_time],
	[:pretrig_mean, :pretrig_rms, :pulse_average, :pulse_rms, :rise_time, :postpeak_deriv, :peak_index, :peak_value, :min_value]))
push!(steps, PerPulseStep(selectfromcriteria, [:pretrig_rms, :pretrig_rms_criteria, :peak_index, :peak_index_criteria, :postpeak_deriv, :postpeak_deriv_criteria], [:selection_good]))
push!(steps, PerPulseStep(filter1lag, [:pulse, :whitenoise_filter], [:filt_value]))
push!(steps, HistogramStep(update_histogram!, [:filt_value_hist, :selection_good, :filt_value]))
push!(steps, ThresholdStep(calibrate_nofit, [:filt_value_dc_hist,:known_energies],[:calibration],:filt_value_dc_hist, counted, 5000, true))
push!(steps, ThresholdStep(compute_whitenoise_filter, [:pulse, :selection_good], [:whitenoise_filter], :selection_good, sum, 100, true))
push!(steps, ThresholdStep(calc_dc, [:pretrig_mean, :filt_value, :selection_good], [:ptm_dc],:filt_value_hist, counted, 3000, true))
push!(steps, PerPulseStep(applycorrection, [:ptm_dc, :pretrig_mean, :filt_value], [:filt_value_dc]))
push!(steps, HistogramStep(update_histogram!, [:filt_value_dc_hist, :selection_good, :filt_value_dc]))
push!(steps, PerPulseStep(apply_calibration, [:calibration, :filt_value_dc], [:energy]) )
push!(steps, HistogramStep(update_histogram!, [:energy_hist, :selection_good, :energy]))

# push!(steps, ToJLDStep([:filt_value,:pretrig_rms, :energy]))
# push!(steps, FreeMemoryStep())

push!(perpulse_symbols, :filt_value, :selection_good, :energy, :pulse, :rowstamp,
	:pretrig_mean, :pretrig_rms, :pulse_average, :pulse_rms, :rise_time, :postpeak_deriv,
	:peak_index, :peak_value, :min_value, :selection_good, :filt_value_dc,:ljh_group)

c=MassChannel()
# c[:ljh_group] = ljh
c[:pretrig_mean] = RunningVector(Float64)
c[:pretrig_rms] = RunningVector(Float64)
c[:pulse_average] = RunningVector(Float64)
c[:pulse_rms] = RunningVector(Float64)
c[:rise_time] = RunningVector(Float64)
c[:postpeak_deriv] = RunningVector(Float64)
c[:peak_index] = RunningVector(UInt16)
c[:peak_value] = RunningVector(UInt16)
c[:min_value] = RunningVector(UInt16)
c[:filt_value] = RunningVector(Float64)
c[:filt_value_dc] = RunningVector(Float64)
c[:pretrig_rms] = RunningVector(Float64)
c[:selection_good] = RunningSumBitVector()
c[:energy] = RunningVector(Float64)
c[:filt_value_hist] = Histogram(0:1:20000)
c[:filt_value_dc_hist] = Histogram(0:1:20000)
c[:energy_hist] = Histogram(0:1:20000)
c[:pulse] = RunningVector(Vector{UInt16})
c[:rowstamp] = RunningVector(Int)
c[:pre_samples] = 100
c[:frame_time] = 1/100000
c[:peak_index_criteria] = (2280,2360)
c[:pretrig_rms_criteria] = (0.0,10.)
c[:postpeak_deriv_criteria] = (0.0,100.0)
c[:known_energies] = [4952, 5898, 6929, 8040]

workstat(n, s::MockPulsesStep, t) = "MockPulsesStep "*@sprintf("%0.2f pulses/s",n/t)
workstat(n, s::PerPulseStep, t) = "PerPulse:$(graphlabel(s)) "*@sprintf("%0.2f pulses/s",n/t)
workstat(n, s::GetPulsesStep, t) = "GetPulsesStep:$(graphlabel(s)) "*@sprintf("%0.2f pulses/s",n/t)
workstat(n, s::ToJLDStep, t) = "ToJLDStep $n executions at "*@sprintf("%0.2f executions/s",n/t)
workstat(n, s::HistogramStep, t) = "HistogramStep:$(inputs(s)[1]) "*@sprintf("%0.2f pulses/s",n/t)
workstat(n, s::ThresholdStep, t) = "ThresholdStep:$(graphlabel(s)) $n executions at "*@sprintf("%0.2f executions/s",n/t)
workstat(n, s::FreeMemoryStep, t) = "FreeMemoryStep $n executions at "*@sprintf("%0.2f executions/s",n/t)

workdone_cumulative = zeros(Int, length(steps))
stepelapsed_cumulative = Array(Float64, length(steps))
errors = Any[]
close(jldopen(filename(c),"w")) # wipe the test file
for i = 1:30
	println("** loop iteration $i")
	#do steps
	for (j,s) in enumerate(steps)
		tstart = time()
		try
			wu = workunits(dostep(s,c))
			stepelapsed = time()-tstart
			workdone_cumulative[j] += wu
			stepelapsed_cumulative[j] += stepelapsed
		catch ex
			println((i,j))
			println(sprint(io->showerror(io, ex)))
			println(sprint(io->Base.show_backtrace(io, catch_backtrace())))
			debug(s,c)
			break
		end

	end

end
[println("step $i: "*workstat(workdone_cumulative[i], steps[i], stepelapsed_cumulative[i])) for i = 1:length(steps)];
[println("$n donethru $(donethru(c[n]))") for n in keys(c)];
jld = jldopen(filename(c), "r+")


text:
ALL_TASK_OFFSETS[("windows-x86_64", v"1.6.1")] = Dict(
    :END => 384,
    :_isexception => 58,
    :_state => 56,
    :bufsz => 360,
    :copy_stack => 368,
    :ctx => 80,
    :donenotify => 24,
    :eh => 72,
    :excstack => 64,
    :gcstack => 376,
    :logstate => 40,
    :next => 0,
    :prio => 62,
    :queue => 8,
    :result => 32,
    :sizeof_ctx => 272,
    :start => 48,
    :sticky => 57,
    :stkbuf => 352,
    :tid => 60,
    :tls => 16,
    :tls_base_context => 6624,
)


text:
# no closure / lifted

Hash = Dict{String, String}

abstract type AbstractScope end

struct Scope <: AbstractScope
    count    :: Ref{Int}
    freevars :: Hash
    bounds   :: Hash
    parent   :: Union{Nothing, AbstractScope}
end

function lookup(d::Hash, k::String)
    get(d, k, nothing)
end

function lookup(d::Scope, k::String) where T
    v = lookup(d.bounds, k)
    v !== nothing && return v
    v = lookup(d.freevars, k)
    v !== nothing && return v
    nothing
end

function new!(scope::Scope)
    Scope(scope.count, Hash(), Hash(), scope)
end

function enter!(scope::Scope, k)
    c = scope.count.x
    scope.count.x = c + 1
    name = "$k._$(string(c, base=32))"
    scope.bounds[k] = name
    name
end

function enter!(scope::Scope, k, v)
    scope.bounds[k] = v
    v
end

function require!(scope::Scope, k::String)
    v = lookup(scope, k)
    v !== nothing && return v
    required_scopes = [scope]
    scope = scope.parent

    while scope !== nothing
        v = lookup(scope, k)
        v !== nothing && begin
            for each in required_scopes
                each.freevars[k] = v
            end
            return v
        end
        push!(required_scopes, scope)
        scope = scope.parent
    end
    throw(NameUnresolved(k))
end

@data DExp begin
    DAssign(reg::String, val::DExp)
    DClosure(freevars::Vector{String}, argnames::Vector{String}, body::DExp)
    DIf(cond::DExp, br1::DExp, br2::DExp)
    DConst{T} :: (v::T) => DExp
    DVar(sym::String)
    DBlock(elts::Vector{DExp})
    DCall(f::DExp, args::Vector{DExp})
    DAttr(subject::DExp, attr::String)
    DList(elts::Vector{DExp})
    DLoc(lineno::Int, colno::Int, val::DExp)
    DStaged(val :: Any)
    DImport(paths::Vector{String}, name::String, actual::String)
    DModule(modname::String, exports::Vector{Tuple{String, String}}, stmts::Vector{DExp})
end

global_scope() = Scope(Ref(0), Dict{String, String}(), Hash(), nothing)

mutable struct ModuleSpec
    name           :: String # for inspect
    path           :: String
    exports        :: Dict{String, String}
    op_prec        :: Dict{String, Int}
    op_asoc        :: Dict{String, Bool}
end

new_module_spec(name::String, path :: String) = ModuleSpec(name, path, Dict{String, String}(), Dict{String, Int}(), Dict{String, Bool}())

macro position_monitor(l, exp)
    quote
        try $exp
        catch e
            if e isa JMLCompilerError
            else
                e = SimpleMessage(string(e))
            end
            throw(Positioned($l.lineno, $l.colno, e))
        end
    end |> esc
end

function scoping_analysis(
    scope::Scope,
    lexp::LExp,
    modules::OrderedDict{String, ModuleSpec},
    curPath::String = "."
)
    JMLPATH::String =
        get(ENV, "JMLPATH") do
            ""
        end
    JMLPATHS = map(string ∘ strip, split(JMLPATH, ";"))
    pushfirst!(JMLPATHS, curPath)

    main = modules["main"]
    if !isa(lexp, LModule)
        lexp = LModule("main", String[], LExp[lexp])
    end

    sa_mod(scope, lexp, main, modules, JMLPATHS)
end

"""
scoping analysis for module
"""
function sa_mod(
    scope::Scope,
    lexp::LExp,
    modulespec::ModuleSpec,
    modules::OrderedDict{String, ModuleSpec},
    JMLPATHS :: Vector{String},
)
    sa(scope, lexp) = sa_mod(scope, lexp, modulespec, modules, JMLPATHS)

    @match lexp begin
        LImport(_) && limport  => import!(limport, false, modulespec, scope, modules, JMLPATHS)
        LStaged(v)             => DStaged(v)
        LAttr(value, attr)     => DAttr(sa(scope, value), attr)
        LList(elts)            => DList([sa(scope, elt) for elt in elts])
        LCall(f, args)         => DCall(sa(scope, f), [sa(scope, arg) for arg in args])
        LBlock(elts)           => DBlock([sa(scope, elt) for elt in elts])
        LVar(s)                => let s = require!(scope, s); DVar(s) end
        LConst(v)              => DConst(v)
        LIf(a, b, c)           => DIf(sa(scope, a), sa(scope, b), sa(scope, c))
        LBin(seq)              => sa(scope, binop_reduce(seq, modulespec.op_prec, modulespec.op_asoc))
        LLoc(l, v)             => DLoc(l.lineno, l.colno, @position_monitor l sa(scope, v))
        LLet(rec, binds, body) =>
            if rec
                ns = new!(scope)
                block = DExp[]
                binds2 :: Vector{Tuple{String, LExp}} = [(enter!(ns, k), v) for (k, v) in binds]
                for (unmangled_name, v) in binds2
                    push!(block, DAssign(unmangled_name, sa(ns, v)))
                end
                push!(block, sa(ns, body))
                DBlock(block)
            else
                ns = new!(scope)
                block :: Vector{DExp} = [
                    DAssign(enter!(ns, k), sa(new!(scope), v)) for (k, v) in binds
                ]
                push!(block, sa(ns, body))
                DBlock(block)
            end
        LFun(args, body) =>
            let ns = new!(scope),
                args = [enter!(ns, arg) for arg in args],
                body = sa(ns, body)
                DClosure(collect(values(ns.freevars)), args, body)
            end

        LModule(modname, args, seq) =>
            let paths = split(modulespec.name, '.'),
                is_current = modname == paths[end],
                modulespec = is_current ? modulespec : begin
                    qualified_name = join([paths..., modname], ".")
                    get!(modules, qualified_name) do
                        new_module_spec(qualified_name, modulespec.path)
                    end
                end,
                scope = is_current ? scope : new!(scope),
                block = [],
                pairs = Dict{String, String}()
                sort!(seq, by=@λ begin
                    LLoc(_, ::LOps) => 1
                    _               => 0
                end) # put LOps to tail

                function analysis_top(each :: LExp)
                    @match each begin
                        LLoc(l, a) =>
                                let f = @position_monitor l analysis_top(a)
                                    () -> DLoc(l.lineno, l.colno, @position_monitor l f())
                                end
                        LDefine(s, v) =>
                                let ss = enter!(scope, s),
                                    _ = pairs[s] = ss,
                                    v = v
                                    () ->
                                    let ns = new!(scope)
                                        DAssign(ss, sa_mod(ns, v, modulespec, modules, JMLPATHS))
                                    end
                                end
                        LForeign(paths, name) =>
                                let ss = enter!(scope, name)
                                    () -> DImport(paths, name, ss)
                                end
                        LImport(_) && limport =>
                                let a = import!(limport, true, modulespec, scope, modules, JMLPATHS)
                                    () -> a
                                end
                        LInfix(opname, prec, is_right) =>
                                let _ = modulespec.op_prec[opname] = prec,
                                    _ = modulespec.op_asoc[opname] = is_right
                                    () -> DConst(nothing)
                                end
                        LOps(names) =>
                                let _ = for each in names
                                            modulespec.exports[each] = require!(scope, each)
                                        end
                                    () -> DConst(nothing)
                                end
                        a => let a = a; () -> sa_mod(scope, a, modulespec, modules, JMLPATHS) end
                    end
                end

                let fs = map(analysis_top, seq)
                    for each in fs
                        push!(block, each())
                    end
                    DModule(modname, [(a, b) for (a, b) in pairs] , block)
                end
            end
    end
end

function import!(
    limport::LImport,
    is_top_level::Bool,
    modulespec :: ModuleSpec,
    scope::Scope,
    modules::OrderedDict{String, ModuleSpec},
    JMLPATHS::Vector{String}
)
    @when LImport(is_qual, paths, name) = limport begin
        qualified_name = join([paths..., name], '.')
        for JMLPATH in JMLPATHS # start loop
            path = joinpath(JMLPATH, paths..., name * ".jml")
            m = get(modules, qualified_name, nothing)
            if m === nothing
                (m, src_code) =
                    try open(path) do f
                        src_code = read(f, String)
                        new_module_spec(qualified_name, path), src_code
                        end
                    catch e
                        continue
                    end
                try
                    rexp = runparser(src_code, :rexp)
                catch e
                    base = SimpleMessage("parser failed for $path")
                    if e isa ParserFailed
                        base = ComposeExc(
                            SimpleMessage(e.msg),
                            base
                        )
                    end
                    throw(
                        ComposeExc(
                            ModulePathNotFound(qualified_name),
                            base
                        )
                    )
                end
                modules[qualified_name] = m
                sa_mod(global_scope(), to_lexp(rexp), m, modules, JMLPATHS)
            end
            if is_qual
                enter!(scope, name)
            else
                for (k, v) in m.exports
                    enter!(scope, k, v)
                end
            end
            if is_top_level
                for (each, _) in m.exports
                    if haskey(m.op_prec, each)
                        modulespec.op_prec[each] = m.op_prec[each]
                        modulespec.op_asoc[each] = m.op_asoc[each]
                    end
                end
            end
            return DConst(nothing)
        end # end loop
        throw(ModulePathNotFound(qualified_name))
    @otherwise
        throw(SimpleMessage("hmmm, you encountered an internal error."))
    end
end

max_forks_repo_name: thautwarm/Rupy


text:
using MyWorkflow
using Documenter
using Weave

# Generate markdown with GitHub flavor
weave("docs/src/weavesample.jmd", "github")

DocMeta.setdocmeta!(MyWorkflow, :DocTestSetup, :(using MyWorkflow); recursive=true)

makedocs(;
    modules=[MyWorkflow],
    authors="Satoshi Terasaki <terasakisatoshi.math@gmail.com>",
    repo="https://github.com/terasakisatoshi/MyWorkflow.jl/blob/{commit}{path}#L{line}",
    sitename="MyWorkflow.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://terasakisatoshi.github.io/MyWorkflow.jl",
        assets=[asset("assets/theorem.css",islocal=true)],
    ),
    pages=[
        "Home" => "index.md",
        "Example(日本語)" => "example.md",
        "MyWorkflow.jl" => "myworkflow.md",
        "weavesample.md",
    ],
)

deploydocs(;
    repo="github.com/terasakisatoshi/MyWorkflow.jl",
)


text:
import Oceananigans.Grids: required_halo_size

"""
    ScalarBiharmonicDiffusivity{F, N, K}

Holds viscosity and diffusivities for models with prescribed isotropic diffusivities.
"""
struct ScalarBiharmonicDiffusivity{F, N, K} <: AbstractScalarBiharmonicDiffusivity{F}
    ν :: N
    κ :: K

    function ScalarBiharmonicDiffusivity{F}(ν::N, κ::K) where {F, N, K}
        return new{F, N, K}(ν, κ)
    end
end

# Aliases that allow specify the floating type, assuming that the discretization is Explicit in time
          ScalarBiharmonicDiffusivity(FT::DataType;         kwargs...) = ScalarBiharmonicDiffusivity(ThreeDimensionalFormulation(), FT; kwargs...)
  VerticalScalarBiharmonicDiffusivity(FT::DataType=Float64; kwargs...) = ScalarBiharmonicDiffusivity(VerticalFormulation(), FT; kwargs...)
HorizontalScalarBiharmonicDiffusivity(FT::DataType=Float64; kwargs...) = ScalarBiharmonicDiffusivity(HorizontalFormulation(), FT; kwargs...)

required_halo_size(::ScalarBiharmonicDiffusivity) = 2

"""
    ScalarBiharmonicDiffusivity(FT=Float64; νh=0, κh=0, νz=nothing, κz=nothing)

Returns parameters for a scalar biharmonic diffusivity model.

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

  - `ν`: Viscosity. `Number`, `AbstractArray`, or `Function(x, y, z, t)`.

  - `κ`: Diffusivity. `Number`, `AbstractArray`, or `Function(x, y, z, t)`, or
         `NamedTuple` of diffusivities with entries for each tracer.

  - `discrete_form`: `Boolean`.

  - `formulation`: formulation used for the discretization of the diffusivity operator.
                   Options are `VerticalFormulation`, `HorizontalFormulation` and
                   `ThreeDimensionalFormulation`.

"""
function ScalarBiharmonicDiffusivity(formulation=ThreeDimensionalFormulation(), FT=Float64;
                                     ν=0, κ=0,
                                     discrete_form = false)

    ν = convert_diffusivity(FT, ν, Val(discrete_form))
    κ = convert_diffusivity(FT, κ, Val(discrete_form))
    return ScalarBiharmonicDiffusivity{typeof(formulation)}(ν, κ)
end

function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F}) where {F}
    κ = tracer_diffusivities(tracers, closure.κ)
    return ScalarBiharmonicDiffusivity{F}(closure.ν, κ)
end

@inline viscosity(closure::ScalarBiharmonicDiffusivity, args...) = closure.ν
@inline diffusivity(closure::ScalarBiharmonicDiffusivity, ::Val{tracer_index}, args...) where tracer_index = closure.κ[tracer_index]

calculate_diffusivities!(diffusivities, closure::ScalarBiharmonicDiffusivity, args...) = nothing

Base.show(io::IO, closure::ScalarBiharmonicDiffusivity{F}) where {F} =
    print(io, "ScalarBiharmonicDiffusivity: " *
              "(ν=$(closure.ν), κ=$(closure.κ)), " *
              "formulation: $F")


text:
struct CompileError
  transform
  args
  err
end

struct Self end
const self = Self()

function Base.showerror(io::IO, err::CompileError)
  println(io, "Error compiling @dynamo $(err.transform) on $(err.args):")
  showerror(io, err.err)
end

function transform end
function refresh end

function dynamo(f, args...)
  try
    ir = transform(f, args...)::Union{IR,Expr,Nothing}
  catch e
    rethrow(CompileError(f, args, e))
  end
  ir isa Expr && return ir
  ir == nothing && return :(args[1](args[2:end]...))
  m = ir.meta::Meta
  ir = varargs!(m, ir)
  argnames!(m, :args)
  _self = splicearg!(m, ir, Symbol("#self#"))
  prewalk!(x -> x === self ? _self : x, ir)
  return update!(m.code, ir)
end

unesc(x) = prewalk(x -> isexpr(x, :escape) ? x.args[1] : x, x)

macro dynamo(ex)
  @capture(shortdef(ex), (name_(args__) = body_) |
                         (name_(args__) where {Ts__} = body_)) ||
    error("@dynamo needs a function definition.")
  Ts = Ts == nothing ? [] : esc.(Ts)
  f, T = isexpr(name, :(::)) ?
    (length(name.args) == 1 ? (esc(gensym()), esc(name.args[1])) : esc.(name.args)) :
    (esc(gensym()), :(Core.Typeof($(esc(name)))))
  gendef = :(@generated ($f::$T)($(esc(:args))...) where $(Ts...) = return IRTools.dynamo($f, args...))
  quote
    $(isexpr(name, :(::)) || esc(:(function $name end)))
    function IRTools.transform(::Type{<:$T}, $(esc.(args)...)) where $(Ts...)
      $(esc(body))
    end
    $gendef
    IRTools.refresh(::$T) where $(Ts...) = (Core.eval($__module__, $(QuoteNode(unesc(gendef)))); return)
  end
end

macro code_ir(dy, ex)
  @capture(ex, f_(args__)) || error("@code_dynamo f(x...)")
  :(transform(typeof($(esc(dy))), meta(typesof($(esc(f)), $(esc.(args)...)))))
end

function recurse!(ir)
  for (x, st) in ir
    isexpr(st.expr, :call) || continue
    ir[x] = Expr(:call, self, st.expr.args...)
  end
end


text:
using Genie, Genie.Router, Genie.Renderer.Html, Stipple
using StippleUI.Button, StippleUI.Drawer

import StippleUI.Layout

Base.@kwdef mutable struct Model <: ReactiveModel
  process::R{Bool} = false
  output::R{String} = ""
  input::R{String} = ""
end

model = Stipple.init(Model())

on(model.process) do _
  if (model.process[])
    model.output[] = model.input[] |> reverse
    model.process[] = false
  end
end

function ui()
  page(
    root(model), class="container",
      [
         Layout.layout(view="lhh LpR lff",style="height: 500px", class="shadow-2 rounded-borders",[
              Drawer.drawer(showifabove="", side="left", bordered="", dashedwidth="200"  , [
                  p("Here comes the mouse :-)")
              ]),
              Layout.page_container([
                Layout.page(padding="",[
                    p([
                      "Input "
                      input("", @bind(:input), @on("keyup.enter", "process = true"))
                    ]),
                    p([
                      Button.button("Action!", click=("process = true"), color="primary")
                    ]),
                    p([
                      "Output: "
                      span("", @text(:output))
                    ])
                ])
              ])
          ])
     ]
  ) |> html
end

route("/", ui)

up()


text:
module TDGL
    using Plots

    struct parameters
        L::Array{Float64,1}
        N::Array{Int64,1}
        dr::Array{Float64,1}
        Hz::Float64
        δt :: Float64
        Maxstep ::Int64
        dirname ::String
        κ::Float64
        calc_A::Bool
    end

    function init_system(L,N,Hz;δt = 0.01,Maxstep=1000,name="output",κ=2.0,calc_A=false)
        dr = L ./ (N .-1)
        return parameters(L,N,dr,Hz,δt,Maxstep,name,κ,calc_A)
    end

    function tdGLsimulation(Δ0,param,T,U)
        Nx = param.N[1]
        Ny = param.N[2]
        δt = param.δt
        name = param.dirname
        Δ = copy(Δ0)
        dΔdt = zeros(ComplexF64,Nx,Ny)
        dUdt = zeros(Float64,2,Nx,Ny)
        count = 0
        if param.calc_A
            U[:,2:Nx+1,2:Ny+1] .= 1
        end
        println("TDGL simulation starts!")
        println("Num | ","\t"," Δ at the center | ","\t","|Δ| at the center | ","\t","dΔdt")
        maps = @animate for i=1:param.Maxstep
            calc_dΔdt!(dΔdt,U,Δ,param,T)
            if param.calc_A
                calc_dUdt!(dUdt,U,Δ,param,T)
            end
            Δ += δt*dΔdt
            if param.calc_A
                U[:,2:Nx+1,2:Ny+1] = U[:,2:Nx+1,2:Ny+1].*exp.(-im*δt.*dUdt[:,:,:]) #-δt*dUdt
            end

            p = heatmap(1:Nx, 1:Ny, abs.(Δ),aspect_ratio=:equal);

            if i % 40 ==1
                count += 1
                savefig(p,"./"*name*"/amp_"*string(count))
                println(i,"\t",Δ[div(Nx,2),div(Ny,2)],"\t",abs(Δ[div(Nx,2),div(Ny,2)]),"\t",dΔdt[div(Nx,2),div(Ny,2)])
            end

        end every 40
        gif(maps, "./"*name*"/GL.gif", fps = 15)
        return Δ,U
    end

    function init_U(param)
        Nx = param.N[1]
        Ny = param.N[2]
        U = zeros(ComplexF64,2,Nx+2,Ny+2) #which includes boundaries
        dr = param.dr
        for ix=0:Nx+1
            for iy=0:Ny+1
                dx = 1
                jx = ix + dx/2
                jy = iy
                r = ixiy2r(jx,jy,param)
                  U[1,ix+1,iy+1] = exp(im*dr[1]*Ax(r,param))
                dy = 1
                jx = ix
                jy = iy +dy/2
                r = ixiy2r(jx,jy,param)
                U[2,ix+1,iy+1] = exp(im*dr[2]*Ay(r,param))
            end
        end
        return U
    end

    Ax(r,param) = -param.Hz*r[2]/2
    Ay(r,param) = param.Hz*r[1]/2

    function ixiy2r(ix,iy,param)
        rx = (ix-1)*param.dr[1]-param.L[1]/2
        ry = (iy-1)*param.dr[2]-param.L[2]/2
        return rx,ry
    end

    function get_Ux(ix,iy,U)
        return get_U(1,ix,iy,U)
    end

    function get_Uy(ix,iy,U)
        return get_U(2,ix,iy,U)
    end

    function get_U(xy,ix,iy,U)
        return U[xy,ix+1,iy+1]
    end

    function plaq(ix,iy,U)
        U1 = get_Ux(ix,iy,U)
        U4 = get_Uy(ix,iy,U)

        jx = ix + 1
        jy = iy
        U2 = get_Uy(jx,jy,U)

        jx = ix
        jy = iy+1
        U3 = get_Ux(jx,jy,U)

        Lzij = U1*U2*conj(U3)*conj(U4)
        return Lzij
    end

    function calc_dΔdt!(dΔdt,U,Δ,param,T)
        Nx = param.N[1]
        Ny = param.N[2]
        dr = param.dr
        for ix=1:Nx
            for iy=1:Ny
                dΔdt[ix,iy] = (1-T)*(1-abs(Δ[ix,iy])^2)*Δ[ix,iy]
                for dx =-1:1
                    jx = ix + dx
                    jy = iy
                    if jx > Nx || jx < 1
                        t = 0
                    else
                        if dx == 0
                            t = -2*Δ[ix,iy]/dr[1]^2
                        else
                            r = ixiy2r(jx,jy,param)
                            U1 = get_Ux(jx,jy,U)
                            t = Δ[jx,jy]*U1^(dx)/dr[1]^2
                        end
                    end
                    dΔdt[ix,iy] += t
                end

                for dy =-1:1
                    jx = ix
                    jy = iy + dy
                    if jy > Ny || jy < 1
                        t = 0
                    else
                        if dy == 0
                            t = -2*Δ[ix,iy]/dr[2]^2
                        else
                            r = ixiy2r(jx,jy,param)
                            U1 = get_Uy(jx,jy,U)
                            t = Δ[jx,jy]*U1^(dy)/dr[2]^2
                        end
                    end
                    dΔdt[ix,iy] += t
                end
            end
        end
    end

    function calc_dUdt!(dUdt,U,Δ,param,T)
        Nx = param.N[1]
        Ny = param.N[2]
        dr = param.dr
        for ix=1:Nx
            for iy=1:Ny
                Lzij = plaq(ix,iy,U)
                dx = 1
                jx = ix + 1
                jy = iy
                if jx <= Nx

                    U1 = get_Ux(ix,iy,U)
                    Jx = imag(conj(Δ[ix,iy])*U1*Δ[jx,jy])
                    #Jx = 0.0
                else
                    Jx = 0.0
                end
                Lzijm = plaq(ix,iy-1,U)

                dUdt[1,ix,iy] = real((1-T)*Jx+((param.κ^2/dr[2]^2)*(Lzij*conj(Lzijm)-1)/im))

                dy = 1
                jx = ix
                jy = iy  + dy
                if jy <= Ny

                    U1 = get_Uy(ix,iy,U)
                    Jy = imag(conj(Δ[ix,iy])*U1*Δ[jx,jy])
                    #Jy = 0.0
                else
                    Jy = 0.0
                end
                Lzimj = plaq(ix-1,iy,U)
                dUdt[2,ix,iy] = real((1-T)*Jy+((param.κ^2/dr[1]^2)*(conj(Lzij)*Lzimj-1)/im))
            end
        end
    end

    function test2()
        L = [25.0,25.0]
        N = [51,51]
        Hz = 0.1
        param = init_system(L,N,Hz,Maxstep = 40000)
        Δ0 = 0.9*ones(ComplexF64,N[1],N[2])
#        Δ0[5,5] += im
        T = 0.2
        U = init_U(param)
        Δ,U = tdGLsimulation(Δ0,param,T,U)
    end
end # module


text:
# =================================================================
# Coupled Van der Pol model
# See https://easychair.org/publications/paper/nrdD
# =================================================================

using ReachabilityAnalysis, Plots

const ey₁ = [0.0, 1.0, 0.0, 0.0]
const ey₂ = [0.0, 0.0, 0.0, 1.0]

@taylorize function vdp_N2_mu1!(dx, x, p, t)
    local μ = 1.0
    x₁, y₁, x₂, y₂ = x

    aux0 = x₂ - x₁
    #
    aux11 = 1 - x₁^2
    aux12 = μ * y₁
    aux13 = aux11 * aux12
    aux15 = aux13 + aux0
    #
    aux21 = 1 - x₂^2
    aux22 = μ * y₂
    aux23 = aux21 * aux22
    aux25 = aux23 - aux0

    dx[1] = y₁
    dx[2] = aux15 - x₁ # (1 - x₁^2)*(μ * y₁) + x₂ - x₁ - x₁
    dx[3] = y₂
    dx[4] = aux25 - x₂ # (1 - x₂^2) * (μ * y₂) - x₂ + x₁ - x₂

    return dx
end

@taylorize function vdp_N2_mu2!(dx, x, p, t)
    local μ = 2.0
    x₁, y₁, x₂, y₂ = x

    aux0 = x₂ - x₁
    #
    aux11 = 1 - x₁^2
    aux12 = μ * y₁
    aux13 = aux11 * aux12
    aux15 = aux13 + aux0
    #
    aux21 = 1 - x₂^2
    aux22 = μ * y₂
    aux23 = aux21 * aux22
    aux25 = aux23 - aux0

    dx[1] = y₁
    dx[2] = aux15 - x₁ # (1 - x₁^2)*(μ * y₁) + x₂ - x₁ - x₁
    dx[3] = y₂
    dx[4] = aux25 - x₂ # (1 - x₂^2) * (μ * y₂) - x₂ + x₁ - x₂

    return dx
end

function vanderpolN2(; μ=1.0)
    if μ == 1.0
        X0 = (1.25 .. 1.55) × (2.35 .. 2.45) × (1.25 .. 1.55) × (2.35 .. 2.45)
        prob = @ivp(x' = vdp_N2_mu1!(x), dim: 4, x(0) ∈ X0)
    elseif μ == 2.0
        X0 = (1.55 .. 1.85) × (2.35 .. 2.45) × (1.55 .. 1.85) × (2.35 .. 2.45)
        X0 = split(X0, [7, 1, 7, 1])
        prob = @ivp(x' = vdp_N2_mu2!(x), dim: 4, x(0) ∈ X0)
    else
        error("the value of μ = $μ is not implemented")
    end

    return prob
end


text:
using DrWatson
@quickactivate("HyperspectraWithNeXL")
# Load the necessary libraries
using NeXLSpectrum
using Gadfly, DataFrames
using PrettyTables
# using Cairo, Fontconfig

# Load the unknown spectra from disk
specpath = joinpath(datadir(),"exp_raw","K2496")
k2496 = loadspectrum.(joinpath(specpath,"K2496_$(i).msa") for i in 1:3)
# Construct and apply the detector to the unknown spectra
det = matching(k2496[1], 132.0, 110, Dict(MShell=>n"Cs"))
k2496 = map(s->apply(s, det), k2496)
# Load the reference spectra and construct the filtered references
refs = references([
  reference(n"Si", joinpath(specpath, "Sanbornite std.msa"), mat"BaSi2O5"),
  reference(n"Ba", joinpath(specpath, "Sanbornite std.msa"), mat"BaSi2O5"),
  reference(n"Ba", joinpath(specpath, "BaCl2 std.msa"), mat"BaCl2"),
  reference(n"O", joinpath(specpath, "MgO std.msa"), mat"MgO"),
  reference(n"Ti", joinpath(specpath, "Ti std.msa"), mat"Ti") ], det)
# Fit the unknown spectra with the filtered references
fs=fit_spectrum(k2496, refs)
# Plot the raw spectrum and the residual
p = plot(fs[1])
#mkpath(joinpath(papersdir(), "Figures"))
#p |> PDF(joinpath(papersdir(),"Figures", "k2496residual.pdf"), 8inch, 3inch)
set_default_plot_size(8inch, 3inch)
display(p)

# Matrix correct the k-ratios and tabulate the results
q=quantify.(fs)
nom = parse(Material, "0.323*O+0.2291*Si+0.018*Ti+0.4299*Ba",name="Nominal")
df = asa(DataFrame, q, nominal=nom) # or `vcat( [ asa(DataFrame, qq) for qq in q ]...)`
df