text:
export translate_resp!, translate_resp, remove_resp!, remove_resp

# =====================================================================

@doc """
    translate_resp!(S, resp_new[, chans=CC, wl=γ])
    translate_resp(S, resp_new[, chans=CC, wl=γ])

Translate the instrument response of seismic data channels `CC` in `S` to
`resp_new`. Replaces field `:resp` with `resp_new` for all affected channels.

    remove_resp!(S, chans=CC, wl=γ])
    remove_resp(S, chans=CC, wl=γ])

Remove (flatten to DC) the instrument response of seismic data channels `cha`
in `S`. Replaces field `:resp` with the appropriate (all-pass) response.

    translate_resp!(Ch, resp_new[, wl=γ])
    translate_resp(Ch, resp_new[, wl=γ])

Translate the instrument response of seismic data in SeisChannel object `Ch` to
`resp_new`. Replaces field `:resp` with `resp_new`.

    remove_resp!(Ch, wl=γ])
    remove_resp(Ch, wl=γ])

Remove (flatten to DC) the instrument response of seismic data in `Ch`.
Replaces field `:resp` with the appropriate (all-pass) response.

Keywords:
* **chans=CC** restricts response translation to channel(s) `CC`. `chans` can be
an Integer, UnitRange, or Array{Int64,1}. By default, all seismic data channels
have their responses translated to `resp_new`.
* **wl=γ** sets the waterlevel to γ (default: `γ` = eps(Float32) ≈ ~1f-7)

### Interaction with the :resp field

`translate_resp` and `remove_resp` only work on a channel `i` that satisfies
`S.resp[i] <: PZResp, PZResp64, MultiStageResp`. In the last case,
`S.resp[i].stage[1]` must be a PZResp or PZResp64, only the first stage of
the response is changed, and the stage gain is ignored; instead, the
sensitivity `S.resp[i].stage[1].a0` is used.

### Poles and zeros should be rad/s

Always check when loading from an unsupported data format. Responses read
from station XML are corrected to rad/s automatically (most use rad/s);
responses read from a SACPZ or SEED RESP file already use rad/s.
""" translate_resp!
function translate_resp!(S::GphysData,
                         resp_new::Union{PZResp, PZResp64};
                         chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
                         wl::Float32=eps(Float32))

  # first ensure that there is something to do
  chans = mkchans(chans, S.n)
  @inbounds for i in chans
    if S.resp[i] != resp_new
      break
    end
    if i == last(chans)
      @info(string(timestamp(), ": nothing done (no valid responses to translate)."))
      return nothing
    end
  end

  # remove channels with inappropriate response types
  k = Int64[]
  for (n,i) in enumerate(chans)
    if (typeof(S.resp[i]) <: Union{PZResp, PZResp64, MultiStageResp}) == false
      push!(k,n)
    end
  end
  deleteat!(chans, k)

  # initialize complex "work" vectors to the largest size we need
  Nx      = nx_max(S, chans)
  N2      = nextpow(2, Nx)
  Xw      = Array{Complex{Float32},1}(undef, N2)
  ff_old  = Array{Complex{Float32},1}(undef, N2)
  ff_new  = Array{Complex{Float32},1}(undef, N2)
  f       = Array{Float32,1}(undef, N2)

  GRPS = get_unique(S, ["fs", "resp", "units"], chans)
  for grp in GRPS

    # get fs, resp
    j = grp[1]
    RT = typeof(getindex(getfield(S, :resp), j))

    # no translating instrument responses unless we're dealing with seismometers
    codes = inst_codes(S)
    kill = falses(length(grp))
    for i = 1:length(grp)
      if codes[i] in seis_inst_codes
        continue
      else
        kill[i] = true
      end
    end
    deleteat!(grp, kill)
    isempty(grp) && continue
    j = grp[1]
    uu = lowercase(S.units[j])

    # onward
    resp_old = RT == MultiStageResp ? deepcopy(S.resp[j].stage[1]) : deepcopy(S.resp[j])
    fs = Float32(S.fs[j])

    if fs > 0.0f0 && resp_old != resp_new

      # we'll need this for logging
      resp_str = string("translate_resp!, wl = ", wl, ", resp_old = ", typeof(resp_old),
        "(a0=", resp_old.a0, ", f0=", resp_old.f0, ", p=", resp_old.p, ", z=", resp_old.z, ")")

      # for accelerometers, working *in* acceleration units, we check: are there any poles below the nyquist?
      if uu == "m/s2"
        P = resp_old.p
        k = trues(length(P))
        for (i,p) in enumerate(P)
          if abs(p)/pi ≤ fs # equivalent to abs(p)/2pi ≤ fn/2 which is the true condition
            k[i] = false
          end
        end
        deleteat!(P, k)

        # Most accelerometers have two complex zeros at the origin, but many XML resp files don't show it.
        if isempty(resp_old.z)
          resp_old.z = map(eltype(resp_old.z), [0.0 + 0.0im, 0.0 - 0.0im])
        end
      end

      # get views to segments from each target channel
      (L,X) = get_views(S, grp)

      # initialize Nx, N2, xre
      Nx = first(L)
      N2 = nextpow(2, Nx)
      xfl, xre = update_resp_vecs!(Xw, f, ff_old, ff_new, N2)
      update_resp!(f, ff_old, ff_new, N2, fs, resp_old, resp_new, wl)

      j = 0
      while j < length(L)
        j = j + 1
        if L[j] != Nx
          Nx = L[j]
          N2 = nextpow(2, Nx)
          xfl, xre = update_resp_vecs!(Xw, f, ff_old, ff_new, N2)
          update_resp!(f, ff_old, ff_new, N2, fs, resp_old, resp_new, wl)
        end

        # copy X[j] to Xw and compute FFT
        fill!(Xw, zero(Complex{Float32}))
        copyto!(Xw, 1, X[j], 1, Nx)
        fft!(Xw)
        broadcast!(*, Xw, Xw, ff_new)
        ifft!(Xw)
        copyto!(X[j], 1, xre, 1, Nx)
      end

      # post-processing: set resp and log to :notes
      for k in grp
        if RT == MultiStageResp
          S.resp[k].stage[1] = deepcopy(resp_new)
        else
          setindex!(S.resp, deepcopy(resp_new), k)
        end
        note!(S, k, resp_str)
      end
    end
  end
 return nothing
end

@doc (@doc translate_resp!)
function translate_resp( S::GphysData,
                    resp_new::Union{PZResp, PZResp64};
                    chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
                    wl::Float32=eps(Float32))

  U = deepcopy(S)
  translate_resp!(U, resp_new, chans=chans, wl=wl)
  return U
end

function translate_resp!(C::GphysChannel,
                    resp_new::Union{PZResp, PZResp64};
                    wl::Float32=eps(Float32))

  # first ensure that there is something to do
  uu = lowercase(C.units)
  fs = Float32(C.fs)
  if any([C.resp == resp_new,
          uu in ("m/s", "m/s2", "m") == false,
          fs ≤ 0.0f0,
          inst_code(C) in seis_inst_codes == false])
    @info(string(timestamp(), ": nothing done (no valid responses to translate)."))
    return nothing
  end

  # initialize complex "work" vectors to the largest size we need
  Nx      = nx_max(C)
  N2      = nextpow(2, Nx)
  Xw      = zeros(Complex{Float32}, N2)
  ff_old  = Array{Complex{Float32},1}(undef, N2)
  ff_new  = Array{Complex{Float32},1}(undef, N2)
  f       = Array{Float32,1}(undef, N2)

  # for accelerometers, working *in* acceleration units, we check: are there any poles below the nyquist?
  if uu == "m/s2"
    P = C.resp.p
    k = trues(length(P))
    for (i,p) in enumerate(P)
      if abs(p)/pi ≤ fs # equivalent to abs(p)/2pi ≤ fn/2 which is the true condition
        k[i] = false
      end
    end
    deleteat!(P, k)

    # Most accelerometers have two complex zeros at the origin, but many XML resp files don't show it.
    if isempty(C.resp.z)
      C.resp.z = map(eltype(C.resp.z), [0.0 + 0.0im, 0.0 - 0.0im])
    end
  end

  # Get views
  if size(C.t,1) == 2
    xfl, xre = update_resp_vecs!(Xw, f, ff_old, ff_new, N2)
    update_resp!(f, ff_old, ff_new, N2, fs, C.resp, resp_new, wl)

    # copy X[j] to Xw and compute FFT
    copyto!(Xw, 1, C.x, 1, Nx)
    fft!(Xw)
    broadcast!(*, Xw, Xw, ff_new)
    ifft!(Xw)
    copyto!(C.x, 1, xre, 1, Nx)
  else
    (L,X) = get_views(C)
    xfl = reinterpret(Float32, Xw)
    xre = view(xfl, 1:2:2*N2-1)
    update_resp!(f, ff_old, ff_new, N2, fs, C.resp, resp_new, wl)

    j = 0
    while j < length(L)
      j = j + 1

      if L[j] != Nx
        Nx = L[j]
        N2 = nextpow(2, Nx)
        xfl, xre = update_resp_vecs!(Xw, f, ff_old, ff_new, N2)
        update_resp!(f, ff_old, ff_new, N2, fs, C.resp, resp_new, wl)
      end

      if j > 1
        fill!(Xw, zero(Complex{Float32}))
      end

      copyto!(Xw, 1, X[j], 1, Nx)
      fft!(Xw)
      broadcast!(*, Xw, Xw, ff_new)
      ifft!(Xw)
      copyto!(X[j], 1, xre, 1, Nx)
    end
  end
  C.resp = deepcopy(resp_new)
  return nothing
end

function translate_resp(C::GphysChannel,
                        resp_new::Union{PZResp, PZResp64};
                        wl::Float32=eps(Float32))

  U = deepcopy(C)
  translate_resp!(U, resp_new, wl=wl)
  return U
end

@doc (@doc translate_resp!)
remove_resp!(S::GphysData;
             chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
             wl::Float32=eps(Float32)) = translate_resp!(S, flat_resp, chans=chans, wl=wl)

@doc (@doc translate_resp!)
remove_resp(S::GphysData;
            chans::Union{Integer, UnitRange, Array{Int64,1}}=Int64[],
            wl::Float32=eps(Float32)) = translate_resp(S, flat_resp, chans=chans, wl=wl)

remove_resp!(Ch::GphysChannel;
             wl::Float32=eps(Float32)) = translate_resp!(Ch, flat_resp, wl=wl)

remove_resp(Ch::GphysChannel;
            wl::Float32=eps(Float32)) = translate_resp(Ch, flat_resp, wl=wl)


text:
function calculate_curvature(xs, fs, i)
    Δf_iph = fs[i+1] - fs[i]    # i plus a half
    Δf_imh = fs[i] - fs[i-1]    # i minus a half
    Δx_iph = xs[i+1] - xs[i]
    Δx_imh = xs[i] - xs[i-1]
    Δx_ic = xs[i+1] - xs[i-1]   # i central

    dfdx_iph = Δf_iph / Δx_iph
    dfdx_imh = Δf_imh / Δx_imh

    2 * (dfdx_iph - dfdx_imh) * (Δx_ic)
end

function calculate_curvature_integral(xs, fs, p, interval, n_intervals, f_range)
    tot_w = 0.0
    result = 0.0
    # Do a small convolution
    for (q,w) in ((-1, 0.25), (0, 0.5), (1, 0.25))
        interval == 1 && q == -1 && continue
        interval == n_intervals && q == 1 && continue
        tot_w += w
        i = p + q
        # Estimate integral of second derivative over interval, use that as a refinement indicator
        # https://mathformeremortals.wordpress.com/2013/01/12/a-numerical-second-derivative-from-three-points/
        result += abs(calculate_curvature(xs, fs, i)) / f_range * w
    end
    result /= tot_w
end

function calculate_refinement(f, xs, fs, n_intervals, intervals_to_refine, n_tot_refinements)
    @assert length(xs) == length(fs)
    n_points = length(xs)

    n_intervals_to_refine = length(intervals_to_refine)
    n_new_points = 2*n_intervals_to_refine
    # Do division of the intervals
    new_xs = zeros(eltype(xs), n_points + n_new_points)
    new_fs = zeros(eltype(fs), n_points + n_new_points)
    new_tot_refinements = zeros(Int, n_intervals + n_intervals_to_refine)
    k = 0
    kk = 0
    for i in 1:n_points
        if iseven(i) # This is a point in an interval
            interval = i ÷ 2
            if interval in intervals_to_refine
                kk += 1
                new_tot_refinements[interval - 1 + kk] = n_tot_refinements[interval] + 1
                new_tot_refinements[interval + kk] = n_tot_refinements[interval] + 1

                k += 1
                new_xs[i - 1 + k] = (xs[i] + xs[i-1]) / 2
                new_fs[i - 1 + k] = f(new_xs[i-1 + k])

                new_xs[i + k] = xs[i]
                new_fs[i + k] = fs[i]

                new_xs[i + 1 + k] = (xs[i+1] + xs[i]) / 2
                new_fs[i + 1 + k] = f(new_xs[i + 1 + k])
                k += 1
            else
                new_tot_refinements[interval + kk] = n_tot_refinements[interval]
                new_xs[i + k] = xs[i]
                new_fs[i + k] = fs[i]
            end
        else
            new_xs[i + k] = xs[i]
            new_fs[i + k] = fs[i]
        end
    end

    return (new_xs, new_fs, new_tot_refinements)
end

"""
    adapted_grid(f, minmax::Tuple{Number, Number}; max_recursions = 7)

Computes a grid `x` on the interval [minmax[1], minmax[2]] so that
`plot(f, x)` gives a smooth "nice" plot.
The method used is to create an initial uniform grid (21 points) and refine intervals
where the second derivative is approximated to be large. When an interval
becomes "straight enough" it is no longer divided. Functions are never evaluated
exactly at the end points of the intervals.

The parameter `max_recusions` computes how many times each interval is allowed to
be refined.
"""
function adapted_grid(f, minmax::Tuple{Real, Real}; max_recursions = 7)
    if minmax[1] >= minmax[2]
        throw(ArgumentError("interval must be given as (min, max)"))
    end

    # When an interval has curvature smaller than this, stop refining it.
    max_curvature = 0.05

    # Initial number of points
    n_points = 21
    n_intervals = n_points ÷ 2
    @assert isodd(n_points)

    xs = collect(range(minmax[1]; stop=minmax[2], length=n_points))
    # Move the first and last interior points a bit closer to the end points
    xs[2] = xs[1] + (xs[2] - xs[1]) * 0.25
    xs[end-1] = xs[end] - (xs[end] - xs[end-1]) * 0.25

    # Wiggle interior points a bit to prevent aliasing and other degenerate cases
    rng = MersenneTwister(1337)
    rand_factor = 0.05
    for i in 2:length(xs)-1
        xs[i] += rand_factor * 2 * (rand(rng) - 0.5) * (xs[i+1] - xs[i-1])
    end

    n_tot_refinements = zeros(Int, n_intervals)

    # We only evaluate the function on interior points
    fs = [NaN; [f(x) for x in xs[2:end-1]]; NaN]
    while true
        curvatures = zeros(n_intervals)
        active = falses(n_intervals)
        isfinite_f = isfinite.(fs)
        min_f, max_f = any(isfinite_f) ? extrema(fs[isfinite_f]) : (0.0, 0.0)
        f_range = max_f - min_f
        # Guard against division by zero later
        if f_range == 0 || !isfinite(f_range)
            f_range = one(f_range)
        end
        # Skip first and last interval
        for interval in 2:n_intervals-1
            p = 2 * interval
            if n_tot_refinements[interval] >= max_recursions
                # Skip intervals that have been refined too much
                active[interval] = false
            elseif !all(isfinite.(fs[[p-1,p,p+1]]))
                active[interval] = true
            else
                curvatures[interval] = calculate_curvature_integral(
                    xs, fs, p, interval, n_intervals, f_range)
                # Only consider intervals with a high enough curvature
                active[interval] = curvatures[interval] > max_curvature
            end
        end
        # Approximate end intervals as being the same curvature as those next to it.
        # This avoids computing the function in the end points
        curvatures[1] = curvatures[2]
        active[1] = active[2]
        curvatures[end] = curvatures[end-1]
        active[end] = active[end-1]

        if all(x -> x >= max_recursions, n_tot_refinements[active])
            break
        end

        n_target_refinements = n_intervals ÷ 2
        interval_candidates = collect(1:n_intervals)[active]
        n_refinements = min(n_target_refinements, length(interval_candidates))
        perm = sortperm(curvatures[active])
        intervals_to_refine = sort(interval_candidates[perm[length(perm) - n_refinements + 1:end]])

        (new_xs, new_fs, new_tot_refinements) = calculate_refinement(
            f, xs, fs, n_intervals, intervals_to_refine, n_tot_refinements)

        xs = new_xs
        fs = new_fs
        n_tot_refinements = new_tot_refinements
        n_points = length(xs)
        n_intervals = n_points ÷ 2
    end

    return xs[2:end-1]
end

ext: jl


text:
#=
The goal here is to provide a demonstration of the method we use to
extract the "fractal dimension", the slope of some quantity vs ε
=#

using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

N = 20_000
D = 5
dt = 0.1
ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = 24.0)
εs = ℯ .^ range(-1, 4, length = 40)
X = trajectory(ds, N*dt; dt, Ttr = 100.0)

midpoints(x) = [(x[i] + x[i+1])/2 for i in 1:length(x)-1]
localslopes(x, y) = [(y[i+1] - y[i])/(x[i+1] - x[i]) for i in 1:length(x)-1]

midpoints(x, n::Int) = [(x[i] + x[i+n])/2 for i in 1:length(x)-n]
localslopes(x, y, n::Int) = [linreg(x[i:i+n], y[i:i+n])[2] for i in 1:length(x)-n]

q = 2
H = genentropy.(Ref(X), εs; q)
C = boxed_correlationsum(X, εs; w = 10, q)

# i = findfirst(z -> z > 0, diff(-H)) - 1
# Y = -H

# i = findfirst(z -> z > 0, C)
# Y = log.(C)
# %%
fig, axs = subplots(3, 2; sharex = true, figsize = (figx, 3figy/4))
for (k, Y) in enumerate((-H, log.(C)))
    i = k==1 ? (findfirst(z -> z > 0, diff(Y)) - 1) : findfirst(isinf, Y)
    i = isnothing(i) ? 1 : max(1, i)
    x, y = log.(εs)[i:end], (Y[i:end])

    tol = 0.2
    dxi = 2

    lrs, slopes = linear_regions(x, y; tol, dxi)

    axs[1, k].plot(x, y)
    for i in 1:length(lrs)-1
        axs[1, k].plot(
            x[lrs[i]:lrs[i+1]], y[lrs[i]:lrs[i+1]];
            marker = "o", ms = 5
        )
    end

    (i1, i2), s = linear_region(x, y; tol, dxi)

    for j in 1:3
        axs[j, k].axvspan(x[i1], x[i2]; alpha = 0.25, color = "gray")
    end

    em = midpoints(x)
    δ = localslopes(x,y)
    axs[2, k].plot(em, δ, marker = "s", ms = 5)

    em = midpoints(x, 5)
    δ = localslopes(x,y, 5)
    axs[3, k].plot(em, δ, marker = "s", ms = 5)
end

axs[1,1].set_ylabel("\$-H_{$q}\$")
axs[1,1].set_yticks(0:-3:-9)
axs[1,2].set_ylabel("\$\\log(C_{$q})\$")
axs[1,2].set_yticks(0:-5:-15)
for i in 2:3, j in 1:2
    axs[i,j].set_yticks(0:2:6)
    axs[i,j].set_ylabel("\$\\delta_$(i==2 ? 1 : 5)\$")
    axs[3,j].set_xlabel("\$\\log(\\varepsilon)\$")
end

fig.tight_layout(;pad = 0.3)
# wsave(plotsdir("paper", "demonstration"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# fastby multi
tic()
#######################################################################
# setting up
#######################################################################
using Revise
using FastGroupBy, DataFrames, BenchmarkTools, SortingLab, ShortStrings
# import Base: getindex, similar, setindex!, size
N = 100_000_000; K = 100
srand(1);
# val = rand(round.(rand(K)*100,4), N);
val = rand(1:5, N);
pool = "id".*dec.(1:100,3);
df = DataFrame(id1 = rand(pool,N), id2 = rand(pool,N), val = val);
byvec = ([df[:id1], df[:id2]]...);
# valvec = val;
fn = sum;

#######################################################################
# convert to ShortStrings
#######################################################################
# normal sort is twice as slow
x = byvec[1] .|> ShorterString;
@benchmark sort($x, by = x->x.size_content)
# BenchmarkTools.Trial:
#   memory estimate:  114.44 MiB
#   allocs estimate:  11
#   --------------
#   minimum time:     594.169 ms (0.11% GC)
#   median time:      1.026 s (41.83% GC)
#   mean time:        905.181 ms (32.13% GC)
#   maximum time:     1.120 s (39.28% GC)
#   --------------
#   samples:          6
#   evals/sample:     1

using SortingAlgorithms
@benchmark sort($x, by = x->x.size_content, alg=RadixSort)
# BenchmarkTools.Trial:
#   memory estimate:  154.59 MiB
#   allocs estimate:  26
#   --------------
#   minimum time:     187.913 ms (0.90% GC)
#   median time:      695.344 ms (72.81% GC)
#   mean time:        608.867 ms (66.90% GC)
#   maximum time:     768.958 ms (72.30% GC)
#   --------------
#   samples:          9
#   evals/sample:     1

x1 = (x->x.size_content).(x)
@benchmark sort(x1, alg=RadixSort) # roughly the same as the by version above

@time sort!(x1, alg= RadixSort)
@time maximum(sizeof, x)
@time minimum(sizeof, x)

byvecss = (ShorterString.(byvec[1]), ShorterString.(byvec[2]))

SortingLab.uint_mapping(o, s::ShorterString) = s.size_content
SortingLab.uint_mapping(o, s::ShortString) = s.size_content

function SortingLab.uint_hist(bits::Vector{ShorterString}, RADIX_SIZE = 16, RADIX_MASK = 0xffff)
    SortingLab.uint_hist(SortingLab.uint_mapping.(Base.Forward, bits))
end

function shortstring_groupby(byvecss, val)
    x1 = copy(byvecss[1]);
    vi = FastGroupBy.ValIndexVector(val, x1);

    x2 = copy(byvecss[2]);
    @time SortingLab.sorttwo!(x2, vi)
    vi
end

@time shortstring_groupby(byvecss, val)
function shortstring_groupby_msd!(byvecss, val)
    vi = FastGroupBy.ValIndexVector(byvecss[2], val);
    @time SortingLab.sorttwo!(byvecss[1], vi)
    vi
end
shortstring_groupby_msd!(byvecss, val)

@benchmark shortstring_groupby_msd!($byvecss, $val)

#######################################################################
# FastGroupBy
# too slow
#######################################################################
@time y1 = FastGroupBy.fastby1(fn, byvec, valvec); # 4.5s 10m
@time y2 = FastGroupBy.fastby2(fn, byvec, valvec); # 4.5s 10m

index = fcollect(length(byvec[1]))

@benchmark grouptwo!($(byvec[1]), $index)

@benchmark radixsort($(byvec[1]))

using SortingLab
@time fsortperm(byvec[1])
@benchmark radixsort($(byvec[1]))

function grp2(v)
    cv = copy(v)
    i = fcollect(length(cv))
    grouptwo!(cv, i)
end
@benchmark grp2(byvec[2])
@code_warntype grp2(byvec[2])

#######################################################################
# Pkg.clone("https://github.com/JuliaData/SplitApplyCombine.jl.git")
# using SplitApplyCombine
#######################################################################

import Base: ht_keyindex2!

function sac(fn, byvec, valvec::Vector{T}) where T
    cm = Dict{Tuple{eltype.(byvec)...}, T}()
    for i = 1:length(valvec)
        v = Tuple(bv[i] for bv in byvec)
        index = ht_keyindex2!(cm, v)
        if index > 0
            @inbounds cm.vals[index] += valvec[i]
        else
            @inbounds Base._setindex!(cm, valvec[i], v, -index)
        end
    end
    l = length(cm)

    # bvout = Tuple(Vector{eltype(bv)}(l) for bv in byvec)
    # for (i, c) in enumerate(keys(cm))
    #     for j in 1:length(c)
    #         bvout[j][i] = c[j]
    #     end
    # end
    # (bvout, values(cm))
end

@benchmark sac($byvec, $valvec)

# @code_warntype sac(sum, byvec, valvec)
res = sac(sum, byvec, valvec)

# byvec = [:id1,:id2]
# valvec = setdiff(names(df), byvec)
# @time y = fastby(fn, df, byvec, :val); # 0.4
# @code_warntype fastby(fn, df, byvec, :val)
@time x = DataFrames.aggregate(df[1,:], [:id1,:id2], sum); # test
@benchmark x = DataFrames.aggregate($df, [:id1,:id2], sum)

using JuliaDB

@time t = table(df[:id1], df[:id2], df[:val], names = [:id1, :id2, :val]);
@time groupreduce(+, t, (:id1,:id2), select = :val)
@time IndexedTables.groupby(sum, t, (:id1,:id2), select = :val)

@time t1 = reindex(t, (:id1,:id2), :val);
@benchmark groupreduce(+, t1, (:id1,:id2), select = :val);
@benchmark IndexedTables.groupby(sum, t1, (:id1,:id2), select = :val);

# @time IndexedTables.summarize(+, t, (:id1,:id2), select = :val)

ti = table(df[:id1], df[:id2], df[:val], names = [:id1, :id2, :val], pkey=(:id1, :id2))
@benchmark groupreduce(+, ti, (:id1,:id2), select = :val)
@benchmark IndexedTables.groupby(sum, ti, (:id1,:id2), select = :val)

using RCall
@rput df
R"""
library(data.table)
setDT(df)
system.time(df[,sum(val), keyby="id1,id2"])
"""

R"""
library(data.table)
setDT(df)
setkey(df, id1, id2)
system.time(df[,sum(val), keyby="id1,id2"])
"""

@rput df
R"""
library(data.table)
setDT(df)
df[,id1:=as.factor(id1)]
df[,id2:=as.factor(id2)]
system.time(df[,sum(val), keyby="id1,id2"])
""" # 0.45 10m

@rput df
R"""
library(data.table)
setDT(df)
df[,id1:=as.factor(id1)]
df[,id2:=as.factor(id2)]
setkey(df, id1, id2)
system.time(df[,sum(val), keyby="id1,id2"])
""" # 0.45 10m

#######################################################################
# setting up
#######################################################################
using Revise
using FastGroupBy, DataFrames, BenchmarkTools, SortingLab, ShortStrings
using SortingAlgorithms
N = 100_000_000; K = 100
srand(1);
r = rand(1:N÷K, N)

@time sortperm(r, alg=RadixSort)

a = rand(rand(Float64, N÷K), N)
@time sort(a)
@time sort(a, alg = RadixSort)

@time fsort(a)

rr = "id".*dec.(1:N÷K, 10)
r = rand(rr, N)

@benchmark fsort($r) samples = 5 seconds = 120
@benchmark fsort($r, false, (11, 0x007ff)) samples = 5 seconds = 120

@time radixsort(r)
@time fsort(r, (11, 0x007ff))

@benchmark fsort($r) samples = 5 seconds = 120
@benchmark fsort($r, (11, 0x007ff)) samples = 5 seconds = 120
@benchmark fsort($r, (13, 0x00001fff)) samples = 5 seconds = 120
@benchmark fsort($r, (22, 0x003fffff)) samples = 5 seconds = 120

@benchmark sort($r, alg=RadixSort)
@benchmark sort($r)

rr = rand(1:Int(N÷2), N)
rr1 = rand(1:Int(N÷2)-1, N)
@time sort(rr)
@time sort(rr1)

using SortingAlgorithms
@benchmark sortperm($r)
@benchmark sortperm($r, alg = RadixSort)

rs = [randstring(8) for i = 1:1_000_000]
rs1 = rand(rs, 100_000_000)

aa(s) = (s |> pointer |> Ptr{UInt32} |> unsafe_load) >> 16

@time sort(rs1, by = aa);
@time radixsort(rs1);
@time SortingLab.fsort(rs1, radix)

radix_opts = (16, 0xffff)

Base.sort!(rs1)

ordr = Base.ord(isless, identity, false, Base.Forward)

using SortingLab, FastGroupBy

function sortandperm!(v::Vector{Int64})
    a = collect(1:length(v))
    SortingLab.sorttwo!(v, a)
    return (v,a)
end

sortandperm(v) = sortandperm!(copy(v))

@benchmark sortandperm(r)
issorted(r)

function sortandperm2!(v::Vector{Int64})
    a = fcollect(length(v))
    SortingLab.sorttwo!(v, a)
    return (v,a)
end

sortandperm2(v) = sortandperm2!(copy(v))

@benchmark sortandperm2($r)

function fsortperm8(v)
    @time a = collect(UInt32(1):UInt32(length(v)))
    @time v1 .= (v .<< 32) .| a
    @time SortingLab.sort32!(v1)
    @time v1 .= ((v .>> 32) .<< 32) .| ((v1 .<< 32) .>> 32)
    @time SortingLab.sort32!(v1)
    @time v1 .& 0xffffffff
end

@time fsortperm8(r);

@benchmark fsortperm8($r)

a = fsortperm8(r)
issorted(r[a])

cv = copy(r)
r_sorted, r_perm = sortandperm!(r)
issorted(r_sorted)
issorted(cv[r_perm])

v = r
min, max = extrema(v)
(diff, o1) = Base.sub_with_overflow(max, min)
(rangelen, o2) = Base.add_with_overflow(diff, oneunit(diff))
if !o1 && !o2 && rangelen < div(n,2)
    return sort_int_range!(v, rangelen, min)
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_parse_demo_settings" begin
     @test isnothing(demo(
                         "%COMMAND%"=>"onoff",
                         "camera_type"=>"per",
                         "camera_position"=>"[1 , 2 , 3]",
                         "alpha"=>30.,
                         "width"=>40,
                         "height"=>30,

                         "normalization"=>0.18,
                         "gamma"=>1.0,
                         "avg_lum"=>0.15,

                         "world_type"=>"B" ,
                         "set_pfm_name"=>"prova.pfm",
                         "set_png_name"=>"prova.png",
                         "samples_per_pixel"=>16,
                         "ONLY_FOR_TESTS"=>true,
                    ))

      @test isnothing(demo(
                         "alpha"=>"30.",
                         "width"=>"40",
                         "height"=>"30",
                         "ONLY_FOR_TESTS"=>true
                    ))

     @test isnothing(demo("ONLY_FOR_TESTS"=>true))

     @test_throws ArgumentError demo("infile"=>"prova.pfm", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("camera_type"=>"new", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("camera_position"=>"[1 , 2 , 3" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("alpha"=>"pi greco" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("width"=>"13" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("height"=>"14.5" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("world_type"=>"C", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("samples_per_pixel"=>3 , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo("world_type"=>"C", "ONLY_FOR_TESTS"=>true)
end

@testset "test_parse_demo_renderer_settings" begin
     @test isnothing(demo(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"onoff",
                    "onoff"=>Dict(
                         "background_color"=>"<1,2,3>",
                         "color"=>"<4,  5,  6>",
                         )
                    ))

     @test isnothing(demo(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"flat",
                    "flat"=>Dict(
                         "background_color"=>"<1,2,3>",
                         )
                    ))

     @test isnothing(demo(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"pathtracer",
                    "pathtracer"=>Dict(
                         "init_state"=>1,
                         "init_seq"=>1,
                         "background_color"=>"<1,2,3>",
                         "num_of_rays"=>15,
                         "max_depth"=>5,
                         "russian_roulette_limit"=>3,
                         )
                    ))

     @test isnothing(demo(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"pointlight",
                    "pointlight"=>Dict(
                         "background_color"=>"<1,2,3>",
                         "ambient_color"=>"<4,  5,  6>",
                         "dark_parameter"=>"1.0",
                         )
                    ))

     @test isnothing(demo("%COMMAND%"=>"onoff", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo("%COMMAND%"=>"flat", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo("%COMMAND%"=>"pathtracer", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo("%COMMAND%"=>"pointlight", "ONLY_FOR_TESTS"=>true))

     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"onoff", "onoff"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"onoff", "onoff"=>Dict("color"=>"[1,2,3]"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"flat", "flat"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"flat", "flat"=>Dict("background_color"=>"[1,2,3]"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("init_state"=>"3.14"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("init_seq"=>-5))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("num_of_rays"=>-1))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("max_depth"=>π))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("russian_roulette_limit"=>-3))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pointlight", "pointlight"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pointlight", "pointlight"=>Dict("ambient_color"=>"[1,2,3]"))
end

##########################################################################################92

@testset "test_parse_demoanimation_settings" begin
     @test isnothing(demo_animation(
                         "%COMMAND%"=>"onoff",
                         "camera_type"=>"per",
                         "camera_position"=>"[1 , 2 , 3]",
                         "width"=>40,
                         "height"=>30,

                         "normalization"=>0.18,
                         "gamma"=>1.0,
                         "avg_lum"=>0.15,

                         "world_type"=>"B" ,
                         "set_anim_name"=>"prova.mp4",
                         "samples_per_pixel"=>16,
                         "ONLY_FOR_TESTS"=>true,
                    ))

     @test isnothing(demo_animation(
                         "width"=>"40",
                         "height"=>"30",
                         "ONLY_FOR_TESTS"=>true
                    ))

     @test isnothing(demo_animation("ONLY_FOR_TESTS"=>true))

     @test_throws ArgumentError demo_animation("infile"=>"prova.pfm", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("camera_type"=>"new", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("camera_position"=>"[1 , 2 , 3" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("alpha"=>"pi greco" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("width"=>"13" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("height"=>"14.5" , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("world_type"=>"C", "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("samples_per_pixel"=>3 , "ONLY_FOR_TESTS"=>true)
     @test_throws ArgumentError demo_animation("world_type"=>"C", "ONLY_FOR_TESTS"=>true)
end

@testset "test_parse_demoanimation_renderer_settings" begin
     @test isnothing(demo_animation(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"onoff",
                    "onoff"=>Dict(
                         "background_color"=>"<1,2,3>",
                         "color"=>"<4,  5,  6>",
                         )
                    ))

     @test isnothing(demo_animation(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"flat",
                    "flat"=>Dict(
                         "background_color"=>"<1,2,3>",
                         )
                    ))

     @test isnothing(demo_animation(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"pathtracer",
                    "pathtracer"=>Dict(
                         "init_state"=>1,
                         "init_seq"=>1,
                         "background_color"=>"<1,2,3>",
                         "num_of_rays"=>15,
                         "max_depth"=>5,
                         "russian_roulette_limit"=>3,
                         )
                    ))

     @test isnothing(demo_animation(
                    "ONLY_FOR_TESTS"=>true,
                    "%COMMAND%"=>"pointlight",
                    "pointlight"=>Dict(
                         "background_color"=>"<1,2,3>",
                         "ambient_color"=>"<4,  5,  6>",
                         "dark_parameter"=>"1.0",
                         )
                    ))

     @test isnothing(demo_animation("%COMMAND%"=>"onoff", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo_animation("%COMMAND%"=>"flat", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo_animation("%COMMAND%"=>"pathtracer", "ONLY_FOR_TESTS"=>true))
     @test isnothing(demo_animation("%COMMAND%"=>"pointlight", "ONLY_FOR_TESTS"=>true))

     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"onoff", "onoff"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"onoff", "onoff"=>Dict("color"=>"[1,2,3]"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"flat", "flat"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"flat", "flat"=>Dict("background_color"=>"[1,2,3]"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("init_state"=>"3.14"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("init_seq"=>-5))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("num_of_rays"=>-1))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("max_depth"=>π))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pathtracer", "pathtracer"=>Dict("russian_roulette_limit"=>-3))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pointlight", "pointlight"=>Dict("background_color"=>"<1,2,3"))
     @test_throws ArgumentError demo_animation("ONLY_FOR_TESTS"=>true, "%COMMAND%"=>"pointlight", "pointlight"=>Dict("ambient_color"=>"[1,2,3]"))
end


text:
"""
Updates the accumulated cost.
"""
function update_cost(acc_cost::Float64, new_cost::Union{Int64,Float64},
                     reduction_type::Symbol)
    if reduction_type == :sum
        return acc_cost + new_cost
    elseif reduction_type == :max
        return max(acc_cost, new_cost)
    end
    throw("Unsupported cost reduction strategy: $(reduction_type)")
end

"""
Numerically stable softmax.
"""
function softmax(x, θ::AbstractFloat=1.0, dim::Integer=1)
    x_exps = exp.((x .- maximum(x)) * θ)
    return x_exps ./ sum(x_exps, dims=dim)
end

"""
Various strategies for computing action selection probabilities.

TODO(kykim): Consider putting these out into a separate jl.
"""
function ucb_probs(ucb_scores)
    return ucb_scores /= sum(ucb_scores)
end

function ucb_softmax_probs(ucb_scores)
    return softmax(ucb_scores, θ=1.0)
end

function expected_cost_probs(est_α_probs, est_α_costs)
    exp_α_costs = est_α_probs .* est_α_costs
    return exp_α_costs /= sum(exp_α_costs)
end

function adaptive_probs(est_α_probs, est_α_costs, nominal_probs, β, γ)
    max_α_cost = maximum(est_α_costs)
    max_α_prob = maximum(est_α_probs)
    max_nominal_prob = maximum(nominal_probs)

    cvar_strategy = est_α_costs .* nominal_probs .+ (max_α_cost / 20) .+ eps()
    cdf_strategy = est_α_probs .* nominal_probs .+ (max_α_prob * max_nominal_prob / 20) .+ eps()
    cvar_strategy /= sum(cvar_strategy)
    cdf_strategy /= sum(cdf_strategy)

    # Mixture weighting.
    a_probs = β * nominal_probs .+ γ * cdf_strategy .+ (1 - β - γ) * cvar_strategy
    return a_probs
end


text:
# # Electric Field

function E(r, t, laser)
    R = geometry(laser).rotation_matrix
    r′ = rotate_coords(R, r)
    inv_c = immutable_cache(laser, :inv_c)
    if unit(eltype(r)) ≠ NoUnits
        λ = laser.λ
        r′ = uconvert.(unit(λ), r′)
    end

    ElectricField = real(E(r′, laser) * g(r′[3], t, laser; inv_c))

    return inv_rotate(R, ElectricField)
end

function E(x::SVector{4}, laser::AbstractLaser)
    inv_c = immutable_cache(laser, :inv_c)
    r = x[SVector{3}(2,3,4)]
    t = x[begin] * inv_c
    E(r, t, laser)
end

function E(r, laser::AbstractLaser)
    @assert length(r) == 3 "The laser is only defined in 3D"
    fill!(mutable_cache(laser), r)
    coords = required_coords(laser, r)

    E_x = Ex(laser, coords)
    update_cache!(laser, :Ex, E_x)
    E_y = Ey(laser, coords)
    update_cache!(laser, :Ey, E_y)
    E_z = Ez(laser, coords)
    update_cache!(laser, :Ez, E_z)

    return SVector{3}(E_x, E_y, E_z)
end

@inline function Ey(laser, coords)
    @unpack ξx, ξy = polarization(laser)
    @unpack Ex = mutable_cache(laser)

    return ξy / ξx * Ex
end


text:
#=
Binomial migration model from Mallet & Barton 1989
max_migration: the maximum number of demes that get migration from the source (n in Mallet & Barton 1989)
migration_distance: distance (in demes) between the source and the receiving deme (i in MAllet & BArton 1989)
α is a coefficient that control for symmetry in the migration (both sides of demes)
α = 0.5 is the default and ensures symmetric migration
The function produces the probability of migration that a deme receives from
itself (migration_distance = 0, offspring does not migrate), and neighboring demes (-migration_distance to +migration_distance)
The sum of binomial_migration(-migration_distance to +migration_distance) should = 1.
The larger max_migration is, the lowest binomial_migration with the same migration_distance
=#
function binomial_migration(max_migration, migration_distance, α = 0.5)
    # Binomial distribution with n = 2 * max_migration and p = α
    prob = Binomial(2 * max_migration, α)
    # Distance of the binomial migration a deme located at migration_distance demes from the source deme.
    # Values larger than max_migration get evaluated at 0.0
    dist = max_migration - migration_distance
    return(pdf(prob, dist))
end

function migration(migration_distance, cline, alpha = 0.5)
    # Object to catch results
    migration_matrix = zeros(length(cline), length(cline))

    # Produces the matrix with the initial structure for migration
    for i in 1:length(cline)               # row
        for j in 1:length(cline)           # column
            if -migration_distance <= (i - j) <= migration_distance
                migration_matrix[i, j] = binomial_migration(migration_distance, i-j, alpha)
            end
        end
    end

    next_generation = (cline * migration_matrix) + correction(migration_distance, cline, alpha)

    return(next_generation)
end

function correction(migration_distance, cline, alpha = 0.5)
    #=
    Correct for borders (imaginary demes beyond the studied ones that contribute to
    migration)
    =#
    correction = zeros(1, length(cline))
    correction_right = map(sum, correction_vector(migration_distance, cline, alpha) * last(cline))
    for i = 1:migration_distance
        correction[i] = correction_right[i]
    end
    correction = reverse(correction; dims = 2)
    correction_left = map(sum, correction_vector(migration_distance, cline, alpha) * first(cline))
    for i = 1:migration_distance
        correction[i] = correction_left[i]
    end
    return(correction)
end

function correction_vector(migration_distance, cline, alpha = 0.5)
    correction_vec = fill(Float64[], 1, migration_distance)
    n = length(cline)
    deme_vec = []

    for i in reverse(1:migration_distance)
        deme_vec = push!(deme_vec, binomial_migration(migration_distance, i, alpha))
        correction_vec[i] = deme_vec
    end
    return(correction_vec)
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export LeftStringWeight

"""
`LeftStringWeight(x)`

| Set                     |           ``\\oplus``  |  ``\\otimes``  | ``\\bar{0}`` | ``\\bar{1}`` |
|:-----------------------:|:----------------------:|:--------------:|:------------:|:------------:|
|``L^*\\cup\\{\\infty\\}``|  longest common prefix |     ``\\cdot`` |``\\infty``   |``\\epsilon`` |

where ``L^*`` is Kleene closure of the set of characters ``L`` and ``\\epsilon`` the empty string.
"""
struct LeftStringWeight <: Semiring
  x::String
  iszero::Bool
end
LeftStringWeight(x::String) = LeftStringWeight(x,false)
reversetype(::Type{S}) where {S <: LeftStringWeight} = RightStringWeight

zero(::Type{LeftStringWeight}) = LeftStringWeight("",true)
one(::Type{LeftStringWeight})  = LeftStringWeight("",false)

# longest common prefix
function lcp(a::LeftStringWeight,b::LeftStringWeight)
  if a.iszero
    return b
  elseif b.iszero
    return a
  else
    cnt = 0
    for x in zip(a.x,b.x)
      if x[1] == x[2]
        cnt +=1
      else
        break
      end
    end
    return LeftStringWeight(a.x[1:cnt])
  end
end

*(a::T, b::T) where {T<:LeftStringWeight}=
((a.iszero) | (b.iszero)) ? zero(T) : T(a.x * b.x)
+(a::LeftStringWeight, b::LeftStringWeight) = lcp(a,b)
function /(a::LeftStringWeight, b::LeftStringWeight)
  if b == zero(LeftStringWeight)
    throw(ErrorException("Cannot divide by zero"))
  elseif a == zero(LeftStringWeight)
    return zero(LeftStringWeight)
  else
    str = ""
    for i=1:length(a.x)
      if i > length(b.x)
        str *= a.x[i]
      end
    end
    return LeftStringWeight(str)
  end
end
reverse(a::LeftStringWeight) = RightStringWeight(reverse(a.x),a.iszero)

#properties
isidempotent(::Type{W}) where {W <: LeftStringWeight} = true
isleft(::Type{W}) where {W <: LeftStringWeight}= true
isweaklydivisible(::Type{W}) where {W <: LeftStringWeight}= true


text:
dilate_ref(img, dims, r) = dilate_ref(img, strel_box(img, dims; r))
dilate_ref(img, se) = extreme_filter(max, img, se)

@testset "dilate" begin
    test_types = [Bool, Int, Float64, Gray{N0f8}, Gray{Float64}]
    for T in test_types
        for N in (1, 2, 3)
            sz = ntuple(_ -> 32, N)
            img = rand(T, sz...)

            out = dilate(img)
            @test out == dilate(img, strel_box(img, ntuple(identity, N); r=1))
            @test out == dilate_ref(img, ntuple(identity, N), 1)

            @test dilate(img; dims=(1,), r=2) == dilate_ref(img, (1,), 2)

            se = centered(rand(Bool, ntuple(_ -> 3, N)))
            @test dilate(img, se) == dilate_ref(img, se)
        end
    end

    img = rand(1:5, 7, 7)
    out = similar(img)
    dilate!(out, img)
    @test out == dilate(img)

    img = rand(RGB, 7, 7)
    @test_throws ArgumentError("color image is not supported") dilate(img)

    # ImageMeta
    @testset "ImageMeta" begin
        A = rand(1:5, 7, 7)
        @test arraydata(dilate(ImageMeta(A))) == dilate(A)
        @test arraydata(dilate(ImageMeta(A); dims=1:2)) == dilate(A; dims=1:2)
    end
end


text:
using SurfaceGeometry

isdefined(:N) || (N=100)
isdefined(:PasiveStabilisation) || (PasiveStabilisation=true)
isdefined(:TopologyStabilisation) || (TopologyStabilisation=true)
isdefined (:NormalVelocities) || (NormalVelocities=true)
isdefined(:ViewerActive) || (ViewerActive=false)

using JLD
data = load("sphere.jld")
points = data["points"]
faces = data["faces"]

function interpolate(p1,v1,p2)
    v2 = zeros(p2)
    for xkey in 1:size(p2,2)
        d = Inf
        vi = 0
        for ykey in 1:size(p1,2)
            dist = norm(p1[:,ykey] - p2[:,xkey])
            if dist<d
                d = dist
                vi = ykey
            end
        end
        v2[:,xkey] = v1[:,vi]
    end
end

function velocity(t,points,faces)
    v = Array(Float64,size(points)...)
    if NormalVelocities==true
        normals = Array(Float64,size(points)...);
        NormalVectors!(normals,points,faces,i->FaceVRing(i,faces))

        for xkey in 1:size(points,2)
            v[:,xkey] = velocity(t,points[:,xkey])
            nx = normals[:,xkey]
            v[:,xkey] = nx*dot(nx,v[:,xkey])
        end
    else
        for xkey in 1:size(points,2)
            v[:,xkey] = velocity(t,points[:,xkey])
        end
    end
    return v
end

function velocity(t,pos)
    x,y,z = pos

    x = x*0.15 + 0.35
    y = y*0.15 + 0.35
    z = z*0.15 + 0.35

    u = 2*sin(pi*x)^2 * sin(2*pi*y) * sin(2*pi*z) * sin(2/3*pi*t)
    v = - sin(2*pi*x) * sin(pi*y)^2 * sin(2*pi*z) * sin(2/3*pi*t)
    w = - sin(2*pi*x) * sin(2*pi*y) * sin(pi*z)^2 * sin(2/3*pi*t)

    [u,v,w] /0.15
end

memory = []

h = 3/N

t = 0.001
p = points
f = faces
push!(memory,(t,p,f))

vp = zeros(p)
pp = zeros(p)

v1 = zeros(p)
v2 = zeros(p)

zc = SurfaceGeometry.Erdmanis2016(C=0.01,ftol=1e-3)

scale = 0.2
par = Elparameters(
                   m_use_fraction = false,
                   m_min_edge_length = 0.7*scale,
                   m_max_edge_length = 1.5*scale,
                   m_max_volume_change = 0.1*scale^3,
                   m_min_curvature_multiplier = 1.0,
                   m_max_curvature_multiplier = 1.0,
                   m_merge_proximity_epsilon = 0.5*scale,
                   m_proximity_epsilon = 0.00001,
m_perform_improvement = true,
m_collision_safety = false,
m_min_triangle_angle = 15,
m_max_triangle_angle = 120,
m_allow_vertex_movement = true, ###
m_use_curvature_when_collapsing = false,
m_use_curvature_when_splitting = false,
m_dt = h
)

for i in 1:35
    println("step $i just started")

    v1 = velocity(t,p,f)
    if PasiveStabilisation==true
        res = stabilise(p,f,v1;vp=interpolate(pp,vp,p),method=zc)
        println("Finit $(res.Finit) \t Fres $(res.Fres)")
        v1 = res.vres
    end

    v2 = velocity(t + h/2,p + h/2*v1,f)
    if PasiveStabilisation==true
        res = stabilise(p + h/2*v1,f,v2;vp=v1,method=zc)
        println("Finit $(res.Finit) \t Fres $(res.Fres)")
        v2 = res.vres
    end

    vp = copy(v2)
    pp = copy(p)

    ### Adams Bachforth oneliner
    #(t1,p1,v1), (t2,p2) = (t2,p2,v2), (t2+h,p2 + h/2/(t2-t1)*( (2*(t2-t1) + h)*v2 - h*v1))

    if TopologyStabilisation==true
        actualdt,p,f = improvemeshcol(p,f,p + h*v2,par)
        t = t + actualdt
    else
        p = p + h*v2
        t = t + h
    end

    push!(memory,(t,p,f))
end

using Escher
if !ViewerActive
    include(Pkg.dir("Escher", "src", "cli", "serve.jl"))
    @spawn escher_serve(5555,Pkg.dir("SurfaceGeometry","examples","viewers"))
    ViewerActive = true
end


text:
using Test, Documenter, COBREXA

ex = quote
    using COBREXA, Tulip
    include(joinpath(dirname(pathof(COBREXA)), "..", "test", "data", "test_models.jl"))
    model = test_LP()
end

# set module-level metadata
DocMeta.setdocmeta!(COBREXA, :DocTestSetup, ex; recursive = true)

@testset "Documentation tests" begin
    doctest(COBREXA; manual = true)
end


text:
#
#----------------------------------------------------------------------
export Stokes
#----------------------------------------------------------------------
"""
 Steady stokes flow

 -ν∇²u⃗ + ∇p = f⃗
        ∇⋅u⃗ = 0

 variational formulation:

 find u⃗, p such that ∀ v⃗,p

          | trial | test |
 -----------------------------------
 vel space|   u⃗   |  v⃗   | ∈ H¹₀(Ω)ᵈ
 -----------------------------------
 pr  space|   p   |  q   | ∈ H¹₀(Ω)

 H(v⃗,u⃗) + (v⃗,∇p) = (v⃗,∇p)   ---(1)
 (q,∇⋅u⃗)         = (q,0)    ---(2)

 [ HH -DD']*[u̲⃗] = [f̲⃗]
 [-DD     ] [p̲̲̲]   [0̲]

 where HH * u̲⃗ = [H  ] * [ux],
                [  H]   [uy]

 is the Diffusion operator acting diagonally on X, Y velocity
 components.

 DD * u̲⃗ = DD * [ux]
               [uy]

 is the divergence operator on the pressure mesh.

 To sove the coupled system of equations, we diagonalize the linear system
 with Schur complements as follows.

 (2) -> (2) + DD HH^-1 *(1)

 [HH -DD']*[u̲⃗] = [f̲⃗]
 [    EE ] [p̲̲̲]   [g]

 where
 EE = -DD HH^-1 * DD'
 g  = -DD HH^-1 * f̲⃗

 is the pressure system
"""
struct Stokes{T,U} <: Equation
    vx::ConvectionDiffusion{T,U}
    vy::ConvectionDiffusion{T,U}

    pr::Field{T}

    px ::Array{T} # pressure forcing (on mshV)
    py ::Array{T} #
    rhs::Array{T} # Stokes rhs

    JrPV::Array{T,2}
    JsPV::Array{T,2}

    mshV::Mesh{T} # velocity mesh
    mshP::Mesh{T} # pressure mesh
    mshD::Mesh{T} # dealias  mesh

    set0!::Function # initial condition
    set∂!::Function # dirichlet boundary condition
    setF!::Function # forcing function
    setν!::Function # viscosity
end
#--------------------------------------#
function Stokes(
                bcVX::Array{Char,1},bcVY::Array{Char,1}#,bcPS::Array{Char,1}
               ,mshV::Mesh,mshD::Mesh,mshP::Mesh
               ;Ti=0.,Tf=0.,dt=0.,k=3)

    # create set functions for vx, vy

    # create convection diffusion objects
    velX = Field(bcVX,mshV)
    velY = Field(bcVY,mshV)

    tstep = TimeStepper(Ti,Tf,dt,k)

    vx = ConvectionDiffusion("vx",velX,velX.uh[1],velY.uh[1],tstep,mshD
                            ,set0!,set∂!,setF!,setν!)
    vy = ConvectionDiffusion("vy",velY,velX.uh[1],velY.uh[1],tstep,mshD,
                            ,set0!,set∂!,setF!,setν!)

#   ps = [
#         ConvectionDiffusion("ps$i",Field(bcPS,mshV),velX.u,velY.u,mshD
#                            ,Tf=1.0,dt=5e-3)
#        for i in 1:k]

    pr = Field(['N','N','N','N'],mshP;k=k-1)

    JrPV = interpMat(mshV.zr,mshP.zr)
    JsPV = interpMat(mshV.zs,mshP.zs)

    return Stokes(vx,vy,pr
                 ,tstep
                 ,JrPV,JsPV
                 ,mshV,mshP,mshD)
end
#----------------------------------------------------------------------
function opStokesLHS(q::Array,sks::Stokes)
    @unpack mshV, mshP, JrPV, JsPV = sks
    @unpack bdfB = sks.tstep

    Mvx = vx.fld.M
    Mvy = vy.fld.M

    Eq = stokesOp(q,mshV,Mvx,Mvy,JrPV,JsPV)
    Eq = gatherScatter(Eq,mshP)

    return Eu
end

function opPrecond(u::Array,sks::Stokes)
    Mu = u
    return Mu
end

function makeStokesRHS!(sks::Stokes)
    @unpack rhs, mshV = sks
    @unpack vx,vy,JrPV,JsPV = sks
    @unpack bdfB = sks.tstep

    fx = approxHlmzInv(vx.f,mshV,bdfB[1]) .+ vx.fld.ub
    fy = approxHlmzInv(vy.f,mshV,bdfB[1]) .+ vy.fld.ub

    rhs  .= diver(fx,fy,mshV,JrPV,JsPV)

#   rhs .= mask(rhs,pr.M) # all Neumann BC
    rhs .= gatherScatter(rhs,mshV)
    return
end

function solveStokes!(sks::Stokes)
    @unpack rhs, mshV, fld = sks
    @unpack u,ub = fld

    opL(u) = opLHS(u,sks)
    opM(u) = opPrecond(u,sks)

    δp = pcg(rhs,opL;opM=opM,mult=mshV.mult,ifv=false)

    return δp
end
#----------------------------------------------------------------------
"""
 Project velocity field to divergence
 free subspace
"""
function pressureProject!(sts::Stokes)

    makeStokesRHS!(sks)

    δp = solveStokes!(sks)

    px,py = diverᵀ(δp,mshV,JrPV,JsPV)

    # Ainv approximate
    px = approxHlmzInv(px,mshV)
    py = approxHlmzInv(py,mshV)

    # correction
    vx .+= px
    vy .+= py
    pr .+= δp

    return
end
#----------------------------------------------------------------------
function update!(sks::Stokes
                ,setBC! =fixU!
                ,setForcing! =fixU!
                ,setVisc! =fixU!)

#   updateHist!(ps)
    updateHist!(vx)
    updateHist!(vy)
    updateHist!(pr)

    update!(tstep)

    return
end
#----------------------------------------------------------------------
export evolve!
#----------------------------------------------------------------------
function evolve!(sks::Stokes
                ,setBC! =fixU!
                ,setForcing! =fixU!
                ,setVisc! =fixU!)

    @unpack vx,vy,pr, tstep = sks
    @unpack time = tstep

    update!(sks)

    # explicit pressure update
    pr .= zero(pr)
    for i=1:pr.k
        pr .+= pr.tstep[i]*pr.uh[i]
    end

    # pressure forcing
    px,py = diverᵀ(pr,mshV,JrPV,JsPV)

    evolve!(vx,setBC!,setForcing! .+ px,setVisc!)
    evolve!(vy,setBC!,setForcing! .+ py,setVisc!)

    # project velocity to divergence free subspace with pressure
    pressureProject(sks)

    return
end
#----------------------------------------------------------------------
export simulate!
#----------------------------------------------------------------------
function simulate!(sks::Stokes,callback!::Function
                  ,setIC! =fixU!
                  ,setBC! =fixU!
                  ,setForcing! =fixU!
                  ,setVisc! =fixU!)

    @unpack fld, mshV = sks
    @unpack time, istep, dt, Tf = sks.tstep

    setIC!(fld.u,mshV.x,mshV.y,time[1])

    callback!(sks)
    while time[1] <= Tf[1]

        evolve!(sks,setBC!,setForcing!,setVisc!)

        callback!(sks)

        if(time[1] < 1e-12) break end

    end

    return
end
#----------------------------------------------------------------------
#


text:
function mytauchen(nz,ρ,σ,s::Int=3)
    """
    need package Distributions
    Discretize the following AR(1) process.
    z_t = ρz_{t-1} + ϵ_t

    Arguments:
    nz::Int64 :grid size # step 1
    ρ::Float64 :persistence of AR(1) process
    σ::Float64 :standard diviation of random components
    s::Int64 :decide maximum and minimum value of node often used s=3.

    return:
    z::array{float64,1} :discretized grid
    Π::array{float64,2} :transition matrix of z
    """
    d = Normal() # standard normal distribution

    ###########################################
    # step 2. decide max and min value of grid
    ##########################################

    σ_z = sqrt(σ^2/(1.0-ρ^2))
    zmin = -s*σ_z
    zmax = s*σ_z
    z = collect(linspace(zmin,zmax,nz))

    ###########################################
    # step 3. store width of between grid points
    ###########################################

    w = z[2] -z[1] #width of between grid points

    ###########################################
    # step 4. Get transition matrix
    ###########################################

    Π = zeros(nz,nz) # transition matrix

    for row in 1:nz
        # Do end points first
        @inbounds Π[row, 1] = cdf(d,(z[1] - ρ*z[row] + w/2) / σ)

        @inbounds Π[row, nz] = ccdf(d,(z[nz] - ρ*z[row] - w/2) / σ)
                               #1.0- cdf((z[nz] - ρ*z[row] - w/2) / σ)
    end

    # fill in the middle columns
    for row in 1:nz
        for col in 2:nz-1
            @inbounds Π[row, col] = (cdf(d,(z[col] - ρ*z[row] + w/2) / σ)
                                      -cdf(d,(z[col] - ρ*z[row] - w/2) / σ))
        end
    end

    return z,Π
end

max_stars_repo_name: tokupooh/Quant-Macro

max_issues_repo_name: tokupooh/Quant-Macro

max_forks_repo_name: tokupooh/Quant-Macro


text:
import SequentialMonteCarlo.SMCModel
using SMCExamples.FiniteFeynmanKac
import SMCExamples.Particles.Int64Particle

setRNGs(12345)

function _getFreqs(model::SMCModel, lM::F, ccsmcio::CCSMCIO, algorithm::Symbol,
  m::Int64, d::Int64) where F<:Function
  ccsmcio.ref1 .= FiniteFeynmanKac.Int642Path(1, d, ccsmcio.n)
  ccsmcio.ref2 .= FiniteFeynmanKac.Int642Path(d^ccsmcio.n, d, ccsmcio.n)
  counts1 = zeros(Int64, d^ccsmcio.n)
  counts2 = zeros(Int64, d^ccsmcio.n)
  result1 = Vector{Float64}(undef, d^ccsmcio.n)
  result2 = Vector{Float64}(undef, d^ccsmcio.n)
  for i = 1:m
    ccXpf!(model, lM, ccsmcio, algorithm)

    counts1[FiniteFeynmanKac.Path2Int64(ccsmcio.ref1, d)] += 1
    counts2[FiniteFeynmanKac.Path2Int64(ccsmcio.ref2, d)] += 1
  end
  result1 .= counts1 ./ m
  result2 .= counts2 ./ m
  return result1, result2
end

function testccsmc()
  d = 4
  n = 4
  ffk = FiniteFeynmanKac.randomFiniteFK(d, n)

  model = FiniteFeynmanKac.makeSMCModel(ffk)
  lM = FiniteFeynmanKac.makelM(ffk)
  densities = Vector{Float64}(undef, d^n)
  for i = 1:length(densities)
    densities[i] = FiniteFeynmanKac.fullDensity(ffk, FiniteFeynmanKac.Int642Path(i, d, n))
  end
  densities ./= sum(densities)

  nsamples = 2^12

  ccsmcio = CCSMCIO{model.particle, model.pScratch}(4, model.maxn)
  freqs1, freqs2 = _getFreqs(model, lM, ccsmcio, :BS, nsamples, d)

  testapproxequal(freqs1, densities, 0.05, false)
  testapproxequal(freqs2, densities, 0.05, false)

  freqs1, freqs2 = _getFreqs(model, error, ccsmcio, :AT, nsamples, d)

  testapproxequal(freqs1, densities, 0.05, false)
  testapproxequal(freqs2, densities, 0.05, false)

  freqs1, freqs2 = _getFreqs(model, lM, ccsmcio, :AS, nsamples, d)

  testapproxequal(freqs1, densities, 0.05, false)
  testapproxequal(freqs2, densities, 0.05, false)
end

@time @testset "Finite FK: ccsmc" begin
  testccsmc()
end


text:
@testset "Mincut" begin

    g = lg.complete_digraph(5)
    cap1 = [
       0.0 2.0 2.0 0.0 0.0
       0.0 0.0 0.0 0.0 3.0
       0.0 1.0 0.0 3.0 0.0
       0.0 0.0 0.0 0.0 1.0
       0.0 0.0 0.0 0.0 0.0
    ];
    (part1, part2, value) = LightGraphsFlows.mincut(g,1,5,cap1,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 4.0
    @test part1 == [1]
    @test sort(part2) == collect(2:5)
    cap2 = [
       0.0 3.0 2.0 0.0 0.0
       0.0 0.0 0.0 0.0 3.0
       0.0 1.0 0.0 3.0 0.0
       0.0 0.0 0.0 0.0 1.5
       0.0 0.0 0.0 0.0 0.0
    ];
    (part1, part2, value) = LightGraphsFlows.mincut(g,1,5,cap2,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 4.5
    @test sort(part1) == collect(1:4)
    @test part2 == [5]

    g2 = lg.DiGraph(5)
    lg.add_edge!(g2,1,2)
    lg.add_edge!(g2,1,3)
    lg.add_edge!(g2,3,4)
    lg.add_edge!(g2,3,2)
    lg.add_edge!(g2,2,5)

    (part1, part2, value) = LightGraphsFlows.mincut(g2,1,5,cap1,LightGraphsFlows.PushRelabelAlgorithm())
    @test value ≈ 3.0
    @test sort(part1) == [1,3,4]
    @test sort(part2) == [2,5]

    #non regression test
    flow_graph = lg.DiGraph(7)
    capacity_matrix = zeros(7,7)
    flow_edges = [
    (1,2,2),(1,3,2),(2,4,4),(2,5,4),
    (3,5,4),(3,6,4),(4,7,1),(5,7,1),(6,7,1)
    ]
    for e in flow_edges
       u, v, f = e
    lg.add_edge!(flow_graph, u, v)
    capacity_matrix[u,v] = f
    end
    (part1, part2, value) = LightGraphsFlows.mincut(flow_graph, 1, 7, capacity_matrix, EdmondsKarpAlgorithm())
    @test value ≈ 3.0
    @test sort(part1) == [1,2,3,4,5,6]
    @test sort(part2) == [7]

end


text:
# Residual Sum of Squares
function compute_objval_rss(
    sim_data::Vector{Float64},
    exp_data::Vector{Float64})::Float64
    error::Float64 = 0.0
    for i in eachindex(exp_data)
        @inbounds error += (sim_data[i] - exp_data[i])^2
    end
    return error
end

# Cosine similarity
function compute_objval_cos(
    sim_data::Vector{Float64},
    exp_data::Vector{Float64})::Float64
    error::Float64 = 1.0 - dot(sim_data, exp_data) / (norm(sim_data) * norm(exp_data))
    return error
end

function conditions_index(condition_name::String)::Int
    if !(condition_name in Sim.conditions)
        error("$condition_name is not defined in Sim.conditions")
    end
    return findfirst(isequal(condition_name), Sim.conditions)
end

function diff_sim_and_exp(
    sim_matrix::Matrix{Float64},
    exp_dict::Dict{String,Array{Float64,1}},
    exp_timepoint::Vector{Float64},
    conditions::Vector{String};
    sim_norm_max::Float64)::Tuple{Vector{Float64},Vector{Float64}}
    sim_result::Vector{Float64} = []
    exp_result::Vector{Float64} = []

    for (idx, condition) in enumerate(conditions)
        if condition in keys(exp_dict)
            append!(sim_result, sim_matrix[Int.(exp_timepoint .+ 1),idx])
            append!(exp_result, exp_dict[condition])
        end
    end

    return (sim_result ./ sim_norm_max, exp_result)
end

# Define an objective function to be minimized.
function objective(indiv_gene)::Float64
    indiv::Vector{Float64} = decode_gene2val(indiv_gene)

    (p, u0) = update_param(indiv)

    if Sim.simulate!(p, u0) isa Nothing
        error::Vector{Float64} = zeros(length(observables))
        for (i, obs_name) in enumerate(observables)
            if isassigned(Exp.experiments, i)
                if length(Sim.normalization) > 0
                    norm_max::Float64 = (
                    Sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
                        Sim.simulations[
                            i,
                            Sim.normalization[obs_name]["timepoint"],
                            [conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
                        ]
                    ) : maximum(
                        Sim.simulations[
                            i,
                            :,
                            [conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
                        ]
                    )
                )
                end
                error[i] = compute_objval_rss(
                diff_sim_and_exp(
                    Sim.simulations[i,:,:],
                    Exp.experiments[i],
                    Exp.get_timepoint(obs_name),
                    Sim.conditions,
                    sim_norm_max=ifelse(
                        length(Sim.normalization) == 0, 1.0, norm_max
                    )
                )...
            )
            end
        end
        return sum(error) # < 1e12
    else
        return 1e12
    end
end


text:
# Autogenerated wrapper script for glslang_jll for x86_64-w64-mingw32-cxx11
export glslangValidator

JLLWrappers.@generate_wrapper_header("glslang")
JLLWrappers.@declare_executable_product(glslangValidator)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        glslangValidator,
        "bin\\glslangValidator.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("output.jl")
parameters = Dict{Symbol, Any}()
parameters[:N] = 25
parameters[:bp_weights] = [0.774074394803123, -0.201004684236153, -0.135080548288772,-0.050951964876636]
stats = ImpvolOutput.write_results(parameters, ARGS[1], :real_GDP, true)


text:
# LSAC JUNE 2007 Form 8LSN75
# https://www.lsac.org/sites/default/files/legacy/docs/default-source/jd-docs/sampleptjune.pdf

# Questions 6-10
using AbstractLogic

"""
Exactly three films—Greed, Harvest, and Limelight—are
shown during a film club’s festival held on Thursday, Friday,
and Saturday. Each film is shown at least once during the
festival but never more than once on a given day. On each day
at least one film is shown. Films are shown one at a time. The
following conditions apply:
"""

# Three films Greed, Harvest, and Limelight as well as 0 no film
# t1 is the first film, t2 the second, and t3 are the last film on Thursday
# f1, f2, f3 and s1, s2, s3 are similarly coded for friday and saturday
logicset = logicalparse("t.1, t.2, t.3, f.1, f.2, f.3, s.1, s.2, s.3  ∈  _, Greed, Harvest, Limelight")

# On each day at least one film is shown. Since 3 is last slot then .3 can't be empty.
logicset = logicalparse("{{j}}.3 != '_'", logicset)

# Let's impose the non-binding contrainst that the highest number 3 is the
# last film shown. Thus if t2|t1 != 0 => t3 != 0
# logicset = logicalparse("t1 != 0 ==> t2 !=0", logicset)
logicset = logicalparse("{{j}}.1 != '_' ==> {{j}}.2 != '_'", logicset)

# Each film is shown at least once during the festival
logicset = logicalparse("{{i}} == 'Greed' {{1,}}", logicset)
logicset = logicalparse("{{i}} == 'Harvest' {{1,}}", logicset)
logicset = logicalparse("{{i}} == 'Limelight' {{1,}}", logicset)

# Never more than once on a given day. If t1 does not equal '_' then t1 != t2.
# Since we know that t3 can't be '_' then we can assert that t2 != t3
logicset = logicalparse("{{j}}.1 != '_' ==> {{j}}.1 != {{j}}.2 &&& {{j}}.1 != {{j}}.3 &&& {{j}}.2 != {{j}}.3", logicset)

# On Thursday Harvest is shown, and no film is shown after it on that day
logicset = logicalparse("t.3 = 'Harvest'", logicset)

# On Friday either Greed or Limelight, but not both, and no film is
# shown after it on that day.
logicset = logicalparse(["f.3 = 'Greed'|'Limelight'; f.2,f.1 != 'Greed','Limelight'"], logicset);

# On Saturday either Greed or Harvest, but not both, is shown, and no film is
# shown after it on that day.
logicset = logicalparse(["s.3 = 'Greed'|'Harvest'; s.1,s.2 != 'Greed','Harvest'"], logicset)

showfeasible(logicset)

range(logicset)

#6. Which one of the following could be a complete and accurate description of
# the order in which the films are shown at the festival?
#(A) Thursday: Limelight, then Harvest; Friday: Limelight; Saturday: Harvest
checkfeasible("t.2 = 'Limelight' ; t.3 = 'Harvest' ; f.3 = 'Limelight' ; s.3 = 'Harvest' ; t.1,f.2,f.1,s.2,s.1='_'", logicset)
# false 0 out of 72 possible combinations true

#(B) Thursday: Harvest; Friday: Greed, then Limelight; Saturday: Limelight, then Greed
checkfeasible("t.3 = 'Harvest' ; f.2 = 'Greed' ; f.3 = 'Limelight' ; s.2 = 'Limelight' ; s.3 = 'Greed'", logicset)
# false 0 out of 72 possible combinations true

#(C) Thursday: Harvest; Friday: Limelight; Saturday: Limelight, then Greed
checkfeasible("t.3 = 'Harvest' ; f.3 = 'Limelight' ; s.2 = 'Limelight' ; s.3 = 'Greed'; t.1, t.2, f.1, f.2, s.1='_' ", logicset)
# true 1 out of 72 possible combinations true

#(D) Thursday: Greed, then Harvest, then Limelight; Friday: Limelight; Saturday: Greed
checkfeasible("t.1 = 'Greed' ; t.2 = 'Harvest' ; t.3 = 'Limelight' ; f.3 = 'Limelight' ; s.3 = 'Greed'", logicset)
# false, 0 out of 72 possible combinations 'true'.

#(E) Thursday: Greed, then Harvest; Friday: Limelight, then Harvest; Saturday: Harvest
checkfeasible("t.2 = 'Greed' ; t.3 = 'Harvest' ; f.2 = 'Limelight' ; f.3 = 'Harvest' ; s.3 = 'Harvest' ", logicset)
# false, 0 out of 72 possible combinations 'true'.

#7. Which one of the following CANNOT be true?
#(A) Harvest is the last film shown on each day of the festival.
checkfeasible("t.3,f.3,s.3 = 'Harvest'", logicset)
# false, 0 out of 64 possible combinations 'true'

# (B) Limelight is shown on each day of the festival.
checkfeasible("{{i}} == 'Limelight' {{3,}}", logicset)
# possible,  10 out of 72 possible combinations 'true'.

#(C) Greed is shown second on each day of the festival.
checkfeasible("f.1,f.2,f.3 |= 'Greed'", logicset)
# possible,  32 out of 72 possible combinations 'true'

# (D) A different film is shown first on each day of the festival.
## This will be a little trickier. Define .0 as first film
logicset2 = logicalparse("t.0, f.0, s.0 ∈ Greed, Harvest, Limelight", logicset)
logicset2 = logicalparse("{{j}}.0 == {{j}}.1 if {{j}}.1 != '_'", logicset2)
logicset2 = logicalparse("{{j}}.1 == '_' &&& {{j}}.2 != '_' ==> {{j}}.0 == {{j}}.2", logicset2)
logicset2 = logicalparse("{{j}}.0 == {{j}}.3 if {{j}}.2 == '_'", logicset2)
checkfeasible("t.0 != f.0, s.0 ; f.0 != s.0", logicset2)
# possible,  8 out of 18 possible combinations 'true'.

# (E) A different film is shown last on each day of the festival.
checkfeasible("t.3, f.3 != s.3; t.3 != f.3", logicset)
# possible,  24 out of 74 possible combinations 'true'.

# If Limelight is never shown again during the festival once Greed is shown,
logicsetreduced = logicalparse("{{i}} = 'Greed' ==> {{>i}} != 'Limelight'", logicset)

# which one of the following is the maximum number of film showings that could
# occur during the festival?
#(A) three
checkfeasible("{{i}} != '_' {{3}}", logicsetreduced, countany=true)
# true,  1 out of 20 possible combinations 'true'.

#(B) four
checkfeasible("{{i}} != '_' {{4}}", logicsetreduced, countany=true)
# true,  6 out of 20 possible combinations 'true'.

#(C) five
checkfeasible("{{i}} != '_' {{5}}", logicsetreduced, countany=true)
# true,  9 out of 20 possible combinations 'true'.

# (D) six
checkfeasible("{{i}} != '_' {{6}}", logicsetreduced, countany=true)
# true,  4 out of 20 possible combinations 'true'.
# * Answer

# (E) seven
checkfeasible("{{i}} != '_' {{7}}", logicsetreduced, countany=true)
# true, 0 out of 20 possible combinations 'true'.

# If Greed is shown exactly three times, Harvest is shown exactly twice, and
# Limelight is shown exactly once, then which one of the following must be true?
logicsetreduced = logicalparse(
  "{{i}} = 'Greed' {{3}}; {{i}} = 'Harvest' {{2}}; {{i}} = 'Limelight' {{1}}", logicset)

#(A) All three films are shown on Thursday.
checkfeasible("t.1 != '_'", logicsetreduced, force=true)
# false, 2 out of 3 possible combinations 'true'.

#(B) Exactly two films are shown on Saturday.
checkfeasible("s.1 = '_'", logicsetreduced, force=true)
# true, 3 out of 3 possible combinations 'true'.

#(C) Limelight and Harvest are both shown on Thursday.
checkfeasible("t.1,t.2,t.3 |= 'Limelight'; t.1,t.2,t.3 |= 'Harvest'", logicsetreduced, force=true)
# false, 2 out of 3 possible combinations 'true'.

#(D) Greed is the only film shown on Saturday.
checkfeasible("s.1,s.2 != '_'; s.3 = 'Greed'", logicsetreduced, force=true)
# false, 0 out of 3 possible combinations 'true'.

#(E) Harvest and Greed are both shown on Friday
checkfeasible("s.1,s.2,s.3 |= 'Harvest'; s.1,s.2,s.3 |= 'Greed'", logicsetreduced, force=true)
# false, 0 out of 3 possible combinations 'true'.

# If Limelight is shown exactly three times, Harvest is shown exactly twice, and
# Greed is shown exactly once...
logicsetreduced = logicalparse(
  "{{i}} = 'Limelight' {{3}}; {{i}} = 'Harvest' {{2}}; {{i}} = 'Greed' {{1}}", logicset)

# Reintroduce first film
logicsetreduced2 = logicalparse("t.0, f.0, s.0 ∈ Greed, Harvest, Limelight", logicsetreduced)
logicsetreduced2 = logicalparse("{{j}}.0 == {{j}}.1 if {{j}}.1 != '_'", logicsetreduced2)
logicsetreduced2 = logicalparse("{{j}}.2 == '_' ==> {{j}}.0 == {{j}}.3", logicsetreduced2)
# Order of operations helps here as the &&& block is evaluated before the if operator
logicsetreduced2 = logicalparse("{{j}}.2 != '_' &&& {{j}}.1 == '_' ==> {{j}}.0 == {{j}}.2", logicsetreduced2)

# then which one of the following is a complete and accurate list of the films
# that could be the first film shown on Thursday?
range(logicsetreduced2)[Symbol("t.0")]
# Greed, Limelight

# (A) Harvest
# (B) Limelight
# (C) Greed, Harvest
# (D) Greed, Limelight
# (E) Greed, Harvest, Limelight


text:
using LyceumBase
using LyceumDevTools.TestUtil
using LyceumMuJoCo
using Test

const LYCEUM_SUITE = [
    (LyceumMuJoCo.PointMass, (), ()),
    (LyceumMuJoCo.ArmHandPickup, (), ()),
]

const GYM = [
    (LyceumMuJoCo.SwimmerV2, (), ()),
    (LyceumMuJoCo.HopperV2, (), ()),
    (LyceumMuJoCo.Walker2DV2, (), ())
]

const DMC = [
    (LyceumMuJoCo.CartpoleSwingup, (), ()),
]

function test_group(group)
    @testset "Testing $etype\n    Args: $args.\n    Kwargs: $kwargs" for (etype, args, kwargs) in group
        testenv_correctness(etype, args...; kwargs...)
    end
end

@testset "LyceumMuJoCo.jl" begin
    @testset "Environments" begin
        @testset "Lyceum Suite" begin test_group(LYCEUM_SUITE) end
        @testset "Gym" begin test_group(GYM) end
        @testset "DMC" begin test_group(DMC) end
    end
end


text:
using Test
using FluidTensors
using LinearAlgebra

@testset "Vector Operations" begin
    a = Vec(1.,2.,3.)
    b = Vec(4.,5.,6.)

    @test a == [1.,2.,3.]
    @test a+b === Vec(5.,7.,9.)
    @test b-a === Vec(3.,3.,3.)
    @test a⋅b ≈ 32.
    @test norm(a) ≈ 3.7416573867739413
    @test norm(im*a) ≈ 3.7416573867739413
    @test 2*a === Vec(2.,4.,6.)
    @test a/2 === Vec(0.5,1.,1.5)
    @test cross(a,b) === Vec(-3.,6.,-3.)
    @test cross(b,a) === Vec(3.,-6.,3.)
end

@testset "Tensor Operations" begin
    A = SymTen(1.,2.,3.,4.,5.,6.)
    @test A == [1. 2. 3.;
                2. 4. 5.;
                3. 5. 6.]
    B = SymTen(2.,3.,4.,5.,6.,7.)
    @test A+B === SymTen(3.,5.,7.,9.,11.,13.)
    @test B-A === SymTen(1.,1.,1.,1.,1.,1.)
    @test 2*A === SymTen(2.,4.,6.,8.,10.,12.)
    @test A/2 === SymTen(0.5,1.,1.5,2.,2.5,3.)
    @test tr(A) == 11.
    @test A:B == tr(A*B)
    @test norm(A) == sqrt(2*tr(A*A))
    @test traceless(SymTen(1.,2.,3.,1.,5.,1.)) == [0. 2. 3.;
                                                   2. 0. 5.;
                                                   3. 5. 0.]
end

@testset "Tensor x Vector operation" begin
    a = Vec(1.,2.,3.)
    b = Vec(4.,5.,6.)
    A = SymTen(1.,2.,3.,4.,5.,6.)
    @test A⋅a === Vec(14.,25.,31.)
    @test symouter(a,b) === SymTen(4.,6.5,9.,10.,13.5,18.)
    @test antisymouter(a,b) === AntiSymTen(-1.5,-3.0,-1.5)
    @test outer(a,b) === Ten(4.,8.,12.,5.,10.,15.,6.,12.,18.)
end

@testset "UniformScaling" begin
    A = SymTen(1,1,1,1,1,1)
    @test A - I === SymTen(0,1,1,0,1,0)
    @test A + I === SymTen(2,1,1,2,1,2)
    @test I - A === SymTen(0,-1,-1,0,-1,0)
end

@testset "VecArrays" begin
    v = VecArray([1. 2.;3. 4.], [5. 6.; 7. 8.], [9. 10.;-1. -2.])
    @test eltype(v) === Vec{Float64}
    @test size(v) == (2,2)
    @test v[2] === Vec{Float64}(3.,7.,-1.)
    @test v[1,2] === Vec{Float64}(2., 6., 10.)
    @test v[2,2] === Vec{Float64}(4., 8., -2.)
    v[1] = Vec(-1.,-2.,-3.)
    @test v.x[1] == -1 && v.y[1] == -2 && v.z[1] == -3.
end

@testset "Ten Mixed type operations" begin
    v = Vec(rand(3)...)
    for a in (SymTen(rand(6)...), AntiSymTen(rand(3)...), Ten(rand(9)...))
        for b in (SymTen(rand(6)...), AntiSymTen(rand(3)...), Ten(rand(9)...))
            for op in (+,-)
                @test op(a,b) ≈ op(Matrix(a),Matrix(b))
            end
            @test a⋅b ≈ a*b
            @test a:b ≈ tr(Matrix(a)'*Matrix(b)) atol=2e-15
            @test Lie(a,b) ≈ (Matrix(a)*Matrix(b) - Matrix(b)*Matrix(a)) atol=2e-15
        end
        @test a⋅v ≈ Matrix(a)*Vector(v)
        @test v⋅a ≈ Matrix(a)'*Vector(v)
        @test -a ≈ -Matrix(a)
        @test a + I ≈ Matrix(a) + I
        @test a - I ≈ Matrix(a) - I
        @test I + a ≈ I + Matrix(a)
        @test I - a ≈ I - Matrix(a)

        for op in (+,-,tr,det)
            @test op(a) ≈ op(Matrix(a)) atol=2e-15
        end

        @test square(a) ≈ Matrix(a)^2

        @test symmetric(a) ≈ (Matrix(a) + Matrix(a)')/2 atol=2e-15
        @test antisymmetric(a) ≈ (Matrix(a) - Matrix(a)')/2 atol=2e-15

    end
end

@testset "Eigenvalues and vectors of SymTen" begin
    for a in (SymTen(rand(6)...), SymTen(rand(6)...), SymTen(rand(6)...))
        l,e = eigvec(a)
        @test all(isapprox.(l, eig(a))) || all(isapprox.((l[3],l[2],l[1]), eig(a)))
        b = eigen(Matrix(a))
        lv = b.values
        ev = b.vectors
        @test all(isapprox.(l, lv))
        @test (e[1] ≈ ev[:,1] || -e[1] ≈ ev[:,1])
        @test (e[2] ≈ ev[:,2] || -e[2] ≈ ev[:,2])
        @test (e[3] ≈ ev[:,3] || -e[3] ≈ ev[:,3])

        @test cross(e[1],e[2]) ≈ e[3]
    end
end


text:
const PY = begin
    if "PYTHON3" ∈ keys(ENV)
        ENV["PYTHON3"]
    else
        if Sys.iswindows()
            "py -3"
        else
            "python3"
        end
    end
end
const PIP = begin
    if "PIP3" ∈ keys(ENV)
        ENV["PIP3"]
    else
        if Sys.iswindows()
            "py -3 -m pip"
        else
            "pip3"
        end
    end
end
const NODE = begin
    if "NODE" ∈ keys(ENV)
        ENV["NODE"]
    else
        NodeJS.nodejs_cmd()
    end
end

#=
Pre-rendering
- In order to prerender KaTeX, the only thing that is required is `node` and then we can just
require `katex.min.js`
- For highlights, we need `node` and also to have the `highlight.js` installed via `npm`.
=#
const JD_CAN_PRERENDER = try success(`$NODE -v`); catch; false; end
const JD_CAN_HIGHLIGHT = try success(`$NODE -e "require('highlight.js')"`); catch;
                                 false; end

#=
Minification
- We use `css_html_js_minify`. To use it, you need python3 and to install it.
- Here we check there is python3, and pip3, and then if we fail to import, we try to
use pip3 to install it.
=#
const JD_HAS_PY3    = try success(`$([e for e in split(PY)]) -V`); catch; false; end
const JD_HAS_PIP3   = try success(`$([e for e in split(PIP)]) -V`); catch; false; end
const JD_CAN_MINIFY = JD_HAS_PY3 && JD_HAS_PIP3 &&
                      try success(`$([e for e in split(PY)]) -m "import css_html_js_minify"`) ||
                          success(`$([e for e in split(PIP)]) install css_html_js_minify`); catch;
                              false; end

#=
Information to the user
- what JuDoc couldn't find
- that the user should do a `build` step after installing
=#
JD_CAN_HIGHLIGHT || begin
    JD_CAN_PRERENDER || begin
        println("""✘ Couldn't find node.js (`$NODE -v` failed).
                → It is required for pre-rendering KaTeX and highlight.js but is not necessary to run JuDoc (cf docs).""")
    end
    println("""✘ Couldn't find highlight.js (`$NODE -e "require('highlight.js')"` failed).
            → It is required for pre-rendering highlight.js but is not necessary to run JuDoc (cf docs).""")
end

JD_CAN_MINIFY || begin
    if JD_HAS_PY3
        println("✘ Couldn't find css_html_js_minify (`$([e for e in split(PY)]) -m \"import css_html_js_minify\"` failed).\n" *
                """→ It is required for minification but is not necessary to run JuDoc (cf docs).""")
    else
        println("""✘ Couldn't find python3 (`$([e for e in split(PY)]) -V` failed).
                → It is required for minification but not necessary to run JuDoc (cf docs).""")
    end
end

all((JD_CAN_HIGHLIGHT, JD_CAN_PRERENDER, JD_CAN_MINIFY, JD_HAS_PY3, JD_HAS_PIP3)) || begin
    println("→ After installing any missing component, please re-build the package (cf docs).")
end


text:
#  Pseudorandom numbers from a truncated Gaussian distribution.

#  This implements the ziggurat algorithm from
#  N. Chopin, "Fast simulation of truncated Gaussian distributions",
#  Stat Comput (2011) 21:275-288
#
#  The code is based on the implementation by Guillaume Dollé, Vincent Mazet
#  available from http:#miv.u-strasbg.fr/mazet/rtnorm/. In particular, the
#  rtchopin_twosided and rtchopin_onesided are derived from their code.
#  Giada Spaccapanico Proietti adapted these functions to the Julia environment using
#  the random number generators implemented in the `Distributions` package.
#

#------------------------------------------------------------
# Pseudorandom numbers from a truncated Gaussian distribution
# The Gaussian has parameters μ (default 0) and σ (default 1)
# and is truncated on the interval [lower,upper].

include("constants.jl")

function Distributions.rand(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.Normal{T},Distributions.Continuous},
) where {T<:Real}
    lower = d.lower
    upper = d.upper
    mu = d.untruncated.μ
    sigma = d.untruncated.σ

    a = (lower - mu) / sigma
    b = (upper - mu) / sigma

    # return(lower[0])

    if (abs(a) > abs(b))

        r = -Distributions.rand(rng, Distributions.TruncatedNormal(mu, sigma, -b, -a))

    elseif (a == -Inf)
        r = Distributions.rand(rng, Distributions.Normal(0, 1))
        # If a in the right tail (a > xmax), use rejection algorithm with a truncated exponential proposal
    elseif (a > _XMAX)
        if (b == Inf)
            r = c_rtexp_onesided(rng, a)
        else
            r = c_rt_exp_c(rng, a, b)
        end
        # If a in the left tail (a < xmin), use rejection algorithm with a Gaussian proposal
    elseif (a < _XMIN)

        if (b == Inf)
            r = c_rtnaive_onesided(rng, a)
        else
            r = c_rtnaive(rng, a, b)
        end
        # In other cases (xmin < a < xmax), use Chopin's algorithm
    else

        if (b == Inf)
            r = rtchopin_onesided(rng, a)
        else
            r = rtchopin_twosided(rng, a, b)
        end
    end

    return sigma * r + mu
end

#------------------------------------------------------------
# Compute y_l from y_k #IDX_OK
function yl(k::Int64)
    if (k == 1) # 1 was 0
        return 0.053513975472   # y_l of the leftmost rectangle
    elseif (k == _N - 1) #_N was _N - 1
        return 0.000914116389555  # y_l of the rightmost rectangle
    elseif (k <= 1954) # 1954 was 1953
        return _YU[k-1]
    else
        return _YU[k+1]
    end
end

#------------------------------------------------------------
# Rejection algorithm with a truncated exponential proposal
function c_rt_exp_c(rng::Distributions.AbstractRNG, a::Float64, b::Float64)

    stop = false
    twoasq = 2 * a^2
    expab = _exp_c(-a * (b - a)) - 1
    z = 0.0
    while (!stop)
        z = _log_c(1 + Random.rand(rng) * expab)
        e = -_log_c(Random.rand(rng))
        stop = (twoasq * e > z^2)
    end
    return a - z / a
end

#------------------------------------------------------------
# One-sided rejection algorithm with a truncated exponential proposal
function c_rtexp_onesided(rng::Distributions.AbstractRNG, a::Float64)

    stop = false
    lambda = 0.5 * (a + sqrt(a^2 + 4))
    z = 0.0
    while (!stop)
        z = a - _log_c(Random.rand(rng)) / lambda
        stop = (-2 * _log_c(Random.rand(rng)) > (z - lambda)^2)
    end
    return z
end

#------------------------------------------------------------
# Naive accept reject (two-sided)
function c_rtnaive(rng::Distributions.AbstractRNG, a::Float64, b::Float64)
    stop = false
    r = 0.0
    while (!stop)
        r = Distributions.rand(rng, Distributions.Normal(0, 1))
        stop = (r >= a) && (r <= b)
    end
    return r
end

#------------------------------------------------------------
# Naive accept reject  (one-sided)
function c_rtnaive_onesided(rng::Distributions.AbstractRNG, a::Float64)
    stop = false
    r = 0.0
    while (!stop)
        r = Distributions.rand(rng, Distributions.Normal(0, 1))
        stop = (r >= a)
    end

    return r
end

#------------------------------------------------------------
# Rejection algorithm with a truncated Rayleigh proposal
function rtrayleigh(rng::Distributions.AbstractRNG, a::Float64)

    stop = false
    asq = (a^2) * 0.5
    v = 0.0
    x = 0.0
    while (!stop)
        v = Random.rand(rng)
        x = asq - _log_c(Random.rand(rng))
        stop = (v^2 * x < a)
    end

    return sqrt(2 * x)
end

#------------------------------------------------------------
# Chopin's one-sided ziggurat algorithm
function rtchopin_onesided(rng::Distributions.AbstractRNG, a::Float64)

    stop = false

    if (a < _XMIN)  #if a is too small, use simple reject algorithm

        r = c_rtnaive_onesided(rng, a)
        stop = true

    elseif (a > _XMAX) # if a is too large, use Devroye's algorithm

        r = c_rtexp_onesided(rng, a)
        stop = true

    else
        i = Int64(_I0 + floor(a * _INVH)) #was without +1
        ka = _NCELL[i] + 1#was without +1
        while (!stop)
            # sample integer between ka and N
            k = Int64(floor(Random.rand(rng) * (_N - ka + 1)) + ka)  # was without +1
            if (k >= _N)  #right tail (last box on the right) # was ==
                r = c_rtexp_onesided(rng, _XCONST[_N])
                stop = true
            elseif (k <= ka + 1)
                sim = _XCONST[k] + (_XCONST[k+1] - _XCONST[k]) * Random.rand(rng)
                if (sim >= a)
                    simy = _YU[k] * Random.rand(rng)
                    if ((simy < yl(k)) || (sim * sim + 2.0 * _log_c(simy) + _ALPHA) < 0)
                        r = sim
                        stop = true
                    end
                end
            else

                u = Random.rand(rng)
                simy = _YU[k] * u
                d = _XCONST[k+1] - _XCONST[k]
                ylk = yl(k)
                if (simy < ylk)  # That's what happens most of the time
                    r = _XCONST[k] + u * d * _YU[k] / ylk
                    stop = true
                else

                    sim = _XCONST[k] + d * Random.rand(rng)

                    # Otherwise, check you're below the pdf curve
                    if ((sim * sim + 2 * _log_c(simy) + _ALPHA) < 0)

                        r = sim
                        stop = true
                    end
                end
            end
        end # end while loop
    end

    return r
end #IDXOK

#------------------------------------------------------------
# Chopin's two-sided ziggurat algorithm
function rtchopin_twosided(rng::Distributions.AbstractRNG, a::Float64, b::Float64)

    kmin = 5                                 # if kb-ka < kmin then use a rejection algorithm
    stop = false

    # Compute ka
    i = Int64(_I0 + floor(a * _INVH)) # was without +1
    ka = _NCELL[i] + 1 #was without +1

    # Compute kb
    if (b >= _XMAX)
        kb = _N
    else
        i = Int64(_I0 + floor(b * _INVH))  # was without +1
        kb = _NCELL[i] + 1 # was without +1
    end

    # If |b-a| is small, use rejection algorithm with a truncated exponential proposal
    if (abs(kb - ka) < kmin)
        r = c_rt_exp_c(rng, a, b)
        stop = true
    end

    while (!stop)

        # Sample integer between ka and kb
        k = Int64(ceil(Distributions.rand(rng, Distributions.Uniform(ka, kb))))  # was without +1

        if (k >= _N) # was ==

            # Right tail
            lbound = _XCONST[_N]
            z = -_log_c(Random.rand(rng))
            e = -_log_c(Random.rand(rng))
            z = z / lbound

            if ((z^2 <= 2 * e) && (z < b - lbound))

                # Accept this proposition, otherwise reject
                r = lbound + z
                stop = true
            end

        elseif ((k <= ka + 1) || (k >= kb - 1 && b < _XMAX))

            # Two leftmost and rightmost regions
            sim = _XCONST[k] + (_XCONST[k+1] - _XCONST[k]) * Random.rand(rng)

            if ((sim >= a) && (sim <= b))

                # Accept this proposition, otherwise reject
                simy = _YU[k] * Random.rand(rng)
                if ((simy < yl(k)) || (sim * sim + 2 * _log_c(simy) + _ALPHA) < 0)

                    r = sim
                    stop = true
                end
            end
        else # All the other boxes

            u = Random.rand(rng)
            simy = _YU[k] * u
            d = _XCONST[k+1] - _XCONST[k]
            ylk = yl(k)
            if (simy < ylk)  # That's what happens most of the time

                r = _XCONST[k] + u * d * _YU[k] / ylk
                stop = true
            else

                sim = _XCONST[k] + d * Random.rand(rng)

                # Otherwise, check you're below the pdf curve
                if ((sim * sim + 2 * _log_c(simy) + _ALPHA) < 0)

                    r = sim
                    stop = true
                end
            end
        end
    end

    return r
end


text:
# run this from test directory on a system with LTspice is_ltspice_installed
# to regenerate all of the test log files.

using LTspice

filelist = [
  "test1.asc",
  "test2.asc",
  "test3.asc",
  "test4.asc",
  "test5.asc",
  "test6.asc",
  "test7.asc",
  "test8.asc",
  "test9.asc",
  "test10.asc",
  "test11.asc",
  "test12.asc",
  "test13.asc",
  "test14.asc",
  "testinc1.asc",
  ]

for file in filelist
  sim = LTspiceSimulation(file)
  run!(sim)
end

#cp("test3.log","temp/test3.log", force=true)


text:
using NLopt

"""

	shimoptim(HA::Array, f0::Vector, shimlims::Tuple; s0::Vector, loss::Function)

Optimize shims using constrained nonlinear optimization and arbitrary loss function

# Inputs
- `HA`: Spherical harmonic basis H times calibration matrix A
- `f0`: Acquired field map
- `shimlims`
  - shimlims[1]: (vector) max shim current on each linear shim coil
  - shimlims[2]: (vector) max shim current on high-order shim coils
  - shimlims[3]: (scalar) max combined high-order shim current
- `s0`: initial guess
- `loss`: loss function

# Output
- `s`: optimized shim currents (including center frequency offset)
"""
function shimoptim(HA::Array, f0::Vector, shimlims::Tuple;
	s0::Vector = zeros(size(HA,2),),
	loss::Function = (s, HA, f0) -> 1/2 * norm(HA*s + f0)^2,
	ftol_rel = 1e-5,
	cflim = 2000    # empirical observation: needs to be finite for numerical stability
	)

	(lin_max, hos_max, hos_sum_max) = shimlims

	opt = Opt(:LN_COBYLA, length(s0))
	opt.ftol_rel = ftol_rel
	opt.min_objective = (s, grad) -> loss(s, HA, f0)

	# shim current limits (box constraints)
	if length(s0) > 4
		opt.lower_bounds = [-cflim; -lin_max; -hos_max]
		opt.upper_bounds = [ cflim;  lin_max;  hos_max]
	elseif length(s0) > 1
		opt.lower_bounds = [-cflim; -lin_max]
		opt.upper_bounds = [ cflim;  lin_max]
	else
		opt.lower_bounds = [-cflim]
		opt.upper_bounds = [ cflim]
	end

	# limit on total HO shim current
	if ~isempty(hos_sum_max) && length(s0) > 4
		inequality_constraint!(opt, (s, grad) -> sum(abs.(s[5:end])) - hos_sum_max)
	end

	(minf,mins,ret) = optimize(opt, s0)

	#numevals = opt.numevals # the number of function evaluations
	#println("Loss=$minf after $numevals iterations (returned $ret)")

	return mins
end

"""
	Optimize shims based on loss involving fieldmap gradient
"""
function shimoptim(dHxA::Array, dHyA::Array, dHzA::Array, g0x::Vector, g0y::Vector, g0z::Vector, shimlims::Tuple, loss::Function;
	s0::Vector = zeros(size(dHxA,2),),
	ftol_rel = 1e-5,
	cflim = 2000    # empirical observation: needs to be finite for numerical stability
	)

	(lin_max, hos_max, hos_sum_max) = shimlims

	opt = Opt(:LN_COBYLA, length(s0))
	opt.ftol_rel = ftol_rel
	opt.min_objective = (s, grad) -> loss(s, dHxA, dHyA, dHzA, g0x, g0y, g0z)

	# shim current limits (box constraints)
	if length(s0) > 4
		opt.lower_bounds = [-cflim; -lin_max; -hos_max]
		opt.upper_bounds = [ cflim;  lin_max;  hos_max]
	else
		opt.lower_bounds = [-cflim; -lin_max]
		opt.upper_bounds = [ cflim;  lin_max]
	end

	# limit on total HO shim current
	if ~isempty(hos_sum_max) && length(s0) > 4
		inequality_constraint!(opt, (s, grad) -> sum(abs.(s[5:end])) - hos_sum_max)
	end

	(minf,mins,ret) = optimize(opt, s0)

	return mins
end


text:
const _GIT_LINK_KEY = :_git_link_key
function gitlink(gw::GitWorker)
    get(gw, _GIT_LINK_KEY) do
        GitLinks.GitLink(gitlink_dir(gw), remote_url(gw))
    end
end

import Base.lock
lock(f::Function, gw::GitWorker; kwargs...) = lock(f, gitlink(gw); kwargs...)

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using Test

# Test that all models run
outdir = joinpath(@__DIR__, "../output")
samplesize = 10
rm(outdir, recursive=true) # NB !!cleans out the 'output' folder!!
mkdir(outdir)

include("../analysis/allscc/runmodel.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) == 60
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_growth.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) ==  4622
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_annual.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) ==  132
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_variability.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) ==  132
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_ARvariability.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) ==  124
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_annualGrowth.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) ==  313
rm(outdir, recursive=true)
mkdir(outdir)

include("../analysis/allscc/runmodel_annualGrowth_ARvariability.jl")
@test sum([length(files) for (root, dirs, files) in walkdir(outdir)]) == 313
rm(outdir, recursive=true)

mkdir(outdir)


text:
################################################################################
###################### DEFINE MASTER PARAMETERS FOR THIS FILE ##################
################################################################################

using Mimi
using Distributions
using CSVFiles
using DataFrames
using CSV
using Random
using StatsBase
using Statistics

include("main_model_growth.jl")

# define the model regions in the order that Mimi returns stuff
myregions = ["EU", "USA", "Other OECD","Former USSR","China","Southeast Asia","Africa","Latin America"]
myyears = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]

# define the output directory
dir_output = joinpath(@__DIR__, "../../output/")

# define number of Monte Carlo runs
samplesize = 50000

# define the seed
masterseed = 22081994

# define the pulse size
scc_pulse_size = 75000.

# define loop parameters that are held constant across rho and adaptation rate values
jj_scen = "RCP4.5 & SSP2"
jj_page09damages = false
jj_permafr = true
jj_seaice = true
jj_civvalue = 1.

# define the persistence (which should differ across runs) - VARY THIS ACROSS WORKERS
jj_ge = parse(Float64, ARGS[1]) * 0.1 # range to be covered: [0:0.1:1;]

# define the adaptation parameter
jj_geadapt = parse(Float64, ARGS[2]) * 0.005 # range to be covered: [0:0.005:0.04] - VARY THIS ACROSS WORKERS

# throw an error to halt the process if zero persistence is combined with non-zero adaptation (which is redundant)
if jj_geadapt > 0. && jj_ge == 0.
    error("This combination of persistence and adaptation is redundant and the process has been interrupted")
end

################################################################################
################### MC MODEL RUNS    ###########################################
################################################################################

# define the output for the Monte Carlo files
dir_MCoutput = string(dir_output, "montecarlo-singleGE/scen", jj_scen, "_permafr", jj_permafr, "_seaice", jj_seaice,
                                    "_page09", jj_page09damages, "/", "ge", jj_ge,
                                    "_geadapt", jj_geadapt,
                                    "_civ", jj_civvalue, "/")

# fix the seed and calculate the SCC
Random.seed!(masterseed)
scc_mcs_object = get_scc_mcs(samplesize, 2020, dir_MCoutput,
                            scenario=jj_scen,
                            pulse_size=scc_pulse_size,
                            use_permafrost=jj_permafr,
                            use_seaice=jj_seaice,
                            use_page09damages=jj_page09damages,
                            ge_minimum=jj_ge,
                            ge_maximum=jj_ge + 10^(-20),
                            ge_mode=jj_ge,
                            civvalue_multiplier=jj_civvalue,
                            geadrate = jj_geadapt)

# save out all relevant SCC summary stats into a Dataframe
df_sccMC_singleGE = DataFrame(damagePAGE09=jj_page09damages,
                            permafr=jj_permafr,
                            seaice=jj_seaice,
                            ge=jj_ge,
                            scen=jj_scen,
                            pulse_size=scc_pulse_size,
                            civvalue=jj_civvalue,
                            ge_adapt = jj_geadapt,
                            mean=mean(scc_mcs_object[:, 1]),
                            median=median(scc_mcs_object[:, 1]),
                            min=minimum(scc_mcs_object[:, 1]),
                            max=maximum(scc_mcs_object[:, 1]),
                            perc25=StatsBase.percentile(scc_mcs_object[:, 1], 25),
                            perc75=StatsBase.percentile(scc_mcs_object[:, 1], 75),
                            sd=Statistics.std(scc_mcs_object[:, 1]),
                            varcoeff=Statistics.std(scc_mcs_object[:, 1]) / mean(scc_mcs_object[:, 1]),
                            perc05=StatsBase.percentile(scc_mcs_object[:, 1], 5),
                            perc95=StatsBase.percentile(scc_mcs_object[:, 1], 95),
                            perc10=StatsBase.percentile(scc_mcs_object[:, 1], 10),
                            perc90=StatsBase.percentile(scc_mcs_object[:, 1], 90),
                            share_zeroSCC = count(scc_mcs_object .== 0.) / samplesize)

# write out the SCC distribution
writedlm(string(dir_MCoutput, "SCC_MC.csv"),
                        scc_mcs_object, ",")

# export the results DataFrame
CSV.write(string(dir_MCoutput, "MimiPageGrowthEffectsResults_SCC_fixedGE.csv"), df_sccMC_singleGE)


text:
# Calculate the temperature of the earth using the simplest model
# @jake
# https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/121

using Unitful, Plots
p_sun = 386e24u"W"      # power output of the sun
radius_a = 6378u"km"    # semi-major axis of the earth
radius_b = 6357u"km"    # semi-minor axis of the earth
orbit_a = 149.6e6u"km"  # distance from the sun to earth
orbit_e = 0.017         # eccentricity of r = a(1-e^2)/(1+ecos(θ))  & time ≈ 365.25 * θ / 360 where θ is in degrees
a = 0.75                # absorptivity of the sun's radiation
e = 0.6                 # emmissivity of the earth (very dependent on cloud cover)
σ = 5.6703e-8u"W*m^-2*K^-4"  # Stefan-Boltzman constant
temp_sky = 3u"K"        # sky temperature

t = (0:0.25:365.25)u"d"       # day of year in 1/4 day increments
θ = 2*π/365.25u"d" .* t       # approximate angle around the sun
r = orbit_a * (1-orbit_e^2) ./ (1 .+ orbit_e .* cos.(θ))    # distance from sun to earth
area_projected = π * radius_a * radius_b    # area of earth facing the sun
ec = sqrt(1-radius_b^2/radius_a^2)  # eccentricity of earth

area_surface = 2*π*radius_a^2*(1 + radius_b^2/(ec*radius_b^2)*atanh(ec)) # surface area of the earth

q_in = p_sun * a * area_projected ./ (4 * π .* r.^2) # total heat impacting the earth

temp_earth = (q_in ./ (e*σ*area_surface) .+ temp_sky^4).^0.25 # temperature of the earth

plot(t*u"d^-1", temp_earth*u"K^-1" .- 273.15, label = false, title = "Temperature of Earth", xlabel = "Day", ylabel = "Temperature [C]")