text:
using Base.Threads

export maptosites!, maptosites_d!

function mt_split(niter::TI, maxthreads=MAX_THREADS[1]) where TI
   nt = minimum([maxthreads, nthreads(), niter])
   # nn = ceil.(TI, linspace(1, niter+1, nt+1))
   nn = ceil.(TI, range(1, stop=niter+1, length=nt+1))
   rgs = [nn[i]:(nn[i+1]-1) for i = 1:nt]
   return nt, rgs
end

function mt_split_interlaced(niter::TI, maxthreads=MAX_THREADS[1]) where TI
   nt = minimum([maxthreads, nthreads(), niter])
   rgs = [ j:nt:niter for j = 1:nt ]
   return nt, rgs
end

function _mt_map_!(f::FT, out, it, inner_loop) where FT
   nt, rg = mt_split(length(it))
   if nt == 1
      inner_loop(f, out, it, 1:length(it))
   else
      OUT = [[out]; [zeros(out) for i = 2:nt]]
      @threads for i = 1:nt
         inner_loop(f, OUT[i], it, rg[i])
      end
      for it = 2:nt
         out .+= OUT[it]
      end
   end
   return out
end

maptosites!(f, out::AbstractVector, it::AbstractIterator) =
   _mt_map_!(f, out, it, maptosites_inner!)

maptosites_d!(f, out::AbstractVector, it::AbstractIterator) =
   _mt_map_!(f, out, it, maptosites_d_inner!)

# ============ assembly over pairs

"""
`mapreduce_sym!{S, T}(out::AbstractVector{S}, f, it::PairIterator{T})`

symmetric variant of `mapreduce!{S, T}(out::AbstractVector{S}, ...)`, summing only
over bonds (i,j) with i < j and adding f(R_ij) to both sites i, j.
"""
function maptosites_inner!(f::FT, out, it::PairIterator, rg) where FT
   nlist = it.nlist
   for n in rg
      if  nlist.i[n] < nlist.j[n]
         f_ = f(nlist.r[n], nlist.R[n]) / 2
         out[nlist.i[n]] += f_
         out[nlist.j[n]] += f_
      end
   end
   return out
end

"""
`mapreduce_antisym!{T}(out::AbstractVector{SVec{T}}, df, it::PairIterator{T})`

anti-symmetric variant of `mapreduce!{S, T}(out::AbstractVector{S}, ...)`, summing only
over bonds (i,j) with i < j and adding f(R_ij) to site j and
-f(R_ij) to site i.
"""
function maptosites_d_inner!(f::FT, out, it::PairIterator, rg) where FT
   nlist = it.nlist
   for n in rg
      if nlist.i[n] < nlist.j[n]
         f_ = f(nlist.r[n], nlist.R[n])
         out[nlist.j[n]] += f_
         out[nlist.i[n]] -= f_
      end
   end
   return out
end

# ============ assembly over sites

function maptosites!(f::FT, out::AbstractVector, it::SiteIterator) where FT
   @threads for i = 1:nsites(it.nlist)
      _, R = neigs(it.nlist, i)
      out[i] = f(R)
   end
   return out
end

function maptosites_d!(df::FT, out::AbstractVector, it::SiteIterator) where FT
   nt = nthreads()
   OUT = [out; [zeros(out) for n = 2:nt]]
   @threads for i = 1:nsites(it.nlist)
      j, R = neigs(it.nlist, i)
      df_ = df(R)
      OUT[threadid()][j] += df_
      OUT[threadid()][i] -= sum(df_)
   end
   for it = 2:n
      out .+= OUT[it]
   end
   return out
end


text:
include("src/OSMPBF.jl")

using .OSMPBF

function main()
    n_nodes = 0
    n_ways = 0
    n_rels = 0

    w::Float64 = Inf64
    e::Float64 = -Inf64
    s::Float64 = Inf64
    n::Float64 = -Inf64

    function node_handler(node)
        n_nodes += 1
        w = min(w, node.lon)
        e = max(e, node.lon)
        n = max(n, node.lat)
        s = min(n, node.lat)
    end

    @time scan_pbf(ARGS[1], nodes=node_handler, ways=w -> n_ways += 1, relations=r -> n_rels += 1)

    @info "PBF file has $n_nodes nodes, $n_ways ways, and $n_rels relations, bbox $n $e $s $w"
end

main()

ext: jl


text:
######### Stan program example  ###########

using Compat, Stan, Test

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoulli = "
  data {
    int<lower=0> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
      y ~ bernoulli(theta);
  }
  "

  bernoullidata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  global stanmodel, rc, optim
  stanmodel = Stanmodel(Optimize(), name="bernoulli",
    model=bernoulli, useMamba=false);

  rc, optim = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    optim[1] |> display
    println()

    println()
    println("Test round(optim[1][\"optimize\"][\"theta\"][1], digits=1) ≈ 0.3")
    @test round(optim[1]["optimize"]["theta"][1], digits=1) ≈ 0.3
  end
end # cd


text:
type Mmul <: Op; w; x; ybuf; dx; Mmul(p::KUparam)=new(p); end

Mmul(d...; o...)=Mmul(KUparam(d...; o...))
Mmul(n::Integer; o...)=Mmul(n, 0; o...)

params(l::Mmul)=Any[l.w]
ninputs(l::Mmul)=1
overwrites(l::Mmul)=false
back_reads_x(l::Mmul)=true
back_reads_y(l::Mmul)=false

# TODO: consolidate part of this code in common with conv.jl
# TODO: upgrade conv.jl, pool.jl, CUDNN.jl to the new CUDNN library

forw(l::Mmul, ::Void; o...)=nothing

function forw(l::Mmul, x; y=nothing, o...)
    l.x = x
    (y, w, x) = initforw(l, x, y; o...)
    A_mul_B!(y, w, x) # y = w * x
end

function initforw(l::Mmul, x, y; train=true, o...)
    (xrows, xcols) = size2(x)
    (wrows, wcols) = size(l.w)
    if isempty(l.w)
        nz(l.w,:init,nothing) || (l.w.init = randn!; l.w.initp = (0,0.01))
        wcols=xrows
        init(l.w, eltype(x), (wrows, wcols))
    end
    ndims(l.w) == 2 || error("ndims(w)!=2")
    eltype(l.w) == eltype(x) || error("eltype mismatch")
    xrows == wcols || error("xrows!=wcols")

    y == nothing && (y = dsimilar!(l, :ybuf, x, (wrows, xcols)))
    # atype(x) == atype(y) || error("atype mismatch")
    eltype(x) == eltype(y) || error("eltype mismatch")
    size(y) == (wrows, xcols) || error("ysize mismatch: $(size(y)) $((wrows,xcols))")
    return ((!train && nz(l.w, :average, false)) ?
            (y, l.w.avg, x) :
            (y, l.w.arr, x))
end

# dw = dy * x'
# dx = w' * dy
function back(l::Mmul, dy; dx=nothing, x=l.x, incr=false, returndx=true, o...)
    if !isdefined(l.w,:diff)
        if x != nothing && issparse(x)
            l.w.diff = CudaSparseMatrixCSR(spzeros(eltype(x), size(l.w)...))
        else
            l.w.diff = fill!(similar(l.w.arr),0)
        end
    end
    if dy == nothing || x == nothing
        !incr && fill!(l.w.diff,0)
    elseif incr
        !isdefined(l.w,:inc) && (l.w.inc = similar(l.w.diff))
        A_mul_Bt!(l.w.inc, dy, x)
        axpy!(1, l.w.inc, l.w.diff)
    else
        A_mul_Bt!(l.w.diff, dy, x)
    end
    if returndx
        if dy == nothing
            error()
        else
            dx == nothing && (dx = similar!(l, :dx, dy, size(x))) # x can be sparse, dy/dx always dense
            At_mul_B!(dx, l.w.arr, dy)
        end
        return dx
    end
end

function ysize(l::Mmul,x)
    (xrows,xcols) = size2(x)
    (wrows,wcols) = size(l.w)
    wcols==0 || wcols==xrows || error("Bad input size")
    return (wrows,xcols)
end

ysize(l::Mmul, ::Void)=nothing

### DEAD CODE:

# ;; This buffer is for notes you don't want to save, and for Lisp evaluation.
# ;; If you want to create a file, visit that file with C-x C-f,
# ;; then enter the text in that file's own buffer.

# x(ndims,ninst)  x(ninst,ndims)
# y  =  w * x	y  =  x * w
# dw = dy * x'    dw = x' * dy
# dx = w' * dy    dx = dy * w'

# function forwT(l::Mmul, xT::CudaSparseMatrixCSR; y=nothing, o...)
#     # xT has been transposed!
#     @show map(summary, (l.w, xT, y))
#     l.x = xT
#     (y, yT, w, xT) = initforwT(l, xT, y; o...)
#     CUSPARSE.csrmm2!('N','T',1f0,xT,w,0f0,yT,'O') # yT = xT * w'
#     CUBLAS.geam!('T','T',1f0,yT,0f0,yT,y) # y = yT' need temp variable :(
# end

# function initforwT(l::Mmul, xT::CudaSparseMatrixCSR, y; train=true, o...)
#     (xcols, xrows) = size2(xT)  # xT has been transposed!
#     (wrows, wcols) = size(l.w)
#     if isempty(l.w)
#         nz(l.w,:init,nothing) || (l.w.init = randn!; l.w.initp = (0,0.01))
#         wcols=xrows
#         init(l.w, eltype(xT), (wrows, wcols))
#     end
#     ndims(l.w) == 2 || error("ndims(w)!=2")
#     eltype(l.w) == eltype(xT) || error("eltype mismatch")
#     xrows == wcols || error("xrows!=wcols")

#     y == nothing && (y = dsimilar!(l, :ybuf, xT, (wrows, xcols)))
#     yT = dsimilar!(l, :ytbuf, xT, (xcols, wrows))
#     # atype(x) == atype(y) || error("atype mismatch")
#     eltype(xT) == eltype(y) == eltype(yT) || error("eltype mismatch")
#     size(y) == (wrows, xcols) || error("ysize mismatch: $(size(y)) $((wrows,xcols))")
#     size(yT) == (xcols, wrows) || error("yTsize mismatch: $(size(yT)) $((xcols,wrows))")
#     return ((!train && nz(l.w, :average, false)) ?
#             (y, yT, l.w.avg, xT) :
#             (y, yT, l.w.arr, xT))
# end

#     if isa(x, CudaSparseMatrixCSR)
#         returndx && error("Cannot return dx for sparse x") # TODO: are we sure?
#         return backT(l,dy,x,incr;o...)
#     end
# # sparse back with xT
# function backT(l::Mmul, dy, xT::CudaSparseMatrixCSR, incr; o...)
#     # dw = dy * x'
#     dyS = sparse(dy)            # TODO: get rid of alloc
#     isdefined(l.w, :diff) || (l.w.diff = CudaSparseMatrixCSR(spzeros(eltype(xT), size(l.w)...)))
#     if incr
#         isdefined(l.w, :inc) || (l.w.inc = similar(l.w.diff))
#         CUSPARSE.gemm!('N','N',dy,xT,l.w.inc,'O')
#         axpy!(1, l.w.inc, l.w.diff) # TODO: check impl, sparsity patterns may be the same
#     else
#         CUSPARSE.gemm!('N','N',dy,xT,l.w.diff,'O')
#     end
# end

# function ysize(l::Mmul,x::CudaSparseMatrixCSR)
#     # x is transposed
#     (xcols,xrows) = size2(x)
#     (wrows,wcols) = size(l.w)
#     wcols==0 || wcols==xrows || error("Bad input size")
#     return (wrows,xcols)
# end

# function initbackd(l::Mmul, dy, x)
#     atype(dy) == atype(x) || error("atype mismatch")
#     eltype(dy) == eltype(x) || error("eltype mismatch")
#     size(dy) == ysize(l,x) || error("ysize mismatch")
#     similar!(l.w, :diff, l.w.arr)
#     incr && similar!(l.w, :inc, l.w.arr)
# end

# function initbacki(l::Mmul, dy, x, incr)
#     atype(dy) == atype(x) || error("atype mismatch")
#     eltype(dy) == eltype(x) || error("eltype mismatch")
#     size(dy) == ysize(l,x) || error("ysize mismatch")
#     similar!(l.w, :diff, l.w.arr)
#     incr && similar!(l.w, :inc, l.w.arr)
# end

# function initbackx(l::Mmul, x, dx)
#     dx == nothing && (dx = similar!(l, :dx, x))
#     issimilar(dx,x) || error("Gradient mismatch")
#     return dx
# end


text:
using Test
using QuantumOptics, CollectiveSpins
using LinearAlgebra

@testset "interaction" begin

cs = CollectiveSpins

r = [1.0,1.0,0.0]
dips = [1.0,0.0,0.0]

pos = cs.geometry.chain(0.4, 3)
dips = [[0.0,0.0,1.0],[0.0,0.0,1.0],[0.0,1.0,0.0]]

S = cs.SpinCollection(pos, dips)
Ω = cs.interaction.OmegaMatrix(S)
Γ = cs.interaction.GammaMatrix(S)

# Collective effects
function F(ri, rj, µi_, µj_)
    rij = ri - rj
    μi = normalize(µi_)
    μj = normalize(µj_)
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    if rij_norm == 0
        2/3.
    else
        β = 2π*rij_norm
        µi'*(µj*(sin(β)/β + cos(β)/β^2 - sin(β)/β^3) + (rijn'*µj)*rijn*(-sin(β)/β - 3*cos(β)/β^2 + 3*sin(β)/β^3))
    end
end
function G(ri, rj, µi_, µj_)
    rij = ri - rj
    μi = normalize(µi_)
    μj = normalize(µj_)
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    if rij_norm == 0
        0.0
    else
        β = 2π*rij_norm
        µi'*(µj*(-cos(β)/β + sin(β)/β^2 + cos(β)/β^3) + (rijn'*µj)*rijn*(cos(β)/β - 3*sin(β)/β^2 - 3*cos(β)/β^3))
    end
end

# Test circular polarizations (complex μ)
N = 10
pos = geometry.ring(0.1,N;distance=true)

dips = normalize.([[1, im^i, 0] for i=1:N])

Γmat1 = [CollectiveSpins.interaction.Gamma(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]
Ωmat1 = [CollectiveSpins.interaction.Omega(pos[i], pos[j], dips[i], dips[j]) for i=1:N, j=1:N]

end #testset


text:
"""

```
merge(r::Vector{ComplexMixutures.Result})
```

This function merges the results of MDDF calculations obtained by running the same
analysis on multiple trajectories, or multiple parts of the same trajectory. It returns
a Result structure of the same type, with all the functions and counters representing averages
of the set provided weighted by the number of frames read in each Result set.

"""
function Base.merge(r::Vector{Result})

    nr = length(r)
    nframes_read = r[1].nframes_read
    error = false
    for ir = 2:nr
        nframes_read += r[ir].nframes_read
        if r[ir].nbins != r[1].nbins
            println(
                "ERROR: To merge Results, the number of bins of the histograms of both sets must be the same.",
            )
        end
        if (r[ir].cutoff - r[1].cutoff) > 1.e-8
            println(
                "ERROR: To merge Results, cutoff distance of the of the histograms of both sets must be the same.",
            )
        end
    end
    if error
        error(" Incompatible set of results to merge. ")
    end

    # List of files and weights
    nfiles = 0
    for ir = 1:nr
        nfiles += length(r[ir].files)
    end
    files = Vector{String}(undef, nfiles)
    weights = Vector{Float64}(undef, nfiles)

    # Final resuls
    R = Result(
        options = r[1].options,
        nbins = r[1].nbins,
        dbulk = r[1].dbulk,
        cutoff = r[1].cutoff,
        irefatom = r[1].irefatom,
        lastframe_read = r[nr].lastframe_read,
        nframes_read = nframes_read,
        autocorrelation = r[1].autocorrelation,
        solute = r[1].solute,
        solvent = r[1].solvent,
        files = files,
        weights = weights,
    )

    # Average results weighting the data considering the number of frames of each data set
    @. R.d = r[1].d
    ifile = 0
    for ir = 1:nr

        w = r[ir].nframes_read / nframes_read

        @. R.mddf += w * r[ir].mddf
        @. R.kb += w * r[ir].kb

        @. R.rdf += w * r[ir].rdf
        @. R.kb_rdf += w * r[ir].kb_rdf

        @. R.md_count += w * r[ir].md_count
        @. R.md_count_random += w * r[ir].md_count_random

        @. R.sum_md_count += w * r[ir].sum_md_count
        @. R.sum_md_count_random += w * r[ir].sum_md_count_random

        @. R.solute_atom += w * r[ir].solute_atom
        @. R.solvent_atom += w * r[ir].solvent_atom

        @. R.rdf_count += w * r[ir].rdf_count
        @. R.rdf_count_random += w * r[ir].rdf_count_random

        @. R.sum_rdf_count += w * r[ir].sum_rdf_count
        @. R.sum_rdf_count_random += w * r[ir].sum_rdf_count_random

        R.density.solute += w * r[ir].density.solute
        R.density.solvent += w * r[ir].density.solvent
        R.density.solvent_bulk += w * r[ir].density.solvent_bulk

        R.volume.total += w * r[ir].volume.total
        R.volume.bulk += w * r[ir].volume.bulk
        R.volume.domain += w * r[ir].volume.domain
        R.volume.shell += w * r[ir].volume.shell

        for j = 1:length(r[ir].files)
            ifile += 1
            R.files[ifile] = r[ir].files[j]
            R.weights[ifile] = w * r[ir].weights[j]
        end

    end

    return R

end


text:
"""
    effects!(reference_grid::DataFrame, model::RegressionModel;
             eff_col=nothing, err_col=:err, typical=mean)

Compute the `effects` as specified in `formula`.

Effects are the model predictions made using values given via the reference
grid. For terms present in the model, but not in the reference grid, then
the typical value of those predictors is used. (In other words, effects are
conditional on the typical value.) The function for computing typical values
is specified via `typical`. Note that this is also applied to categorical contrasts,
thus yielding an average of the contrast, weighted by the balance of levels in the data
set used to fit the model.

By default, the column corresponding to the response variable in the formula
is overwritten with the effects, but an alternative column for the effects can
be specified by `eff_col`. If this column does not exist, it is created.
Pointwise standard errors are written into the column specified by `err_col`.

The reference grid must contain columns for all predictors in the formula.
(Interactions are computed automatically.) Contrasts must match the contrasts
used to fit the model; using other contrasts will lead to undefined behavior.

Interaction terms are computed in the same way for any regression model: as the
product of the lower-order terms. Typical values of lower terms thus propagate up
into the interaction term in the same way that any value would.

The use of typical values for excluded effects differs from other approaches
such as "partial effects" used in R packages like [`remef`](https://github.com/hohenstein/remef/).
The R package [`effects`](https://cran.r-project.org/web/packages/effects/)
is similar in approach, but due to differing languages and licenses,
no source code was inspected and there is no attempt at API compatibility or
even similarity.

The approach for computing effect is based on the effects plots described here:

Fox, John (2003). Effect Displays in R for Generalised Linear Models.
Journal of Statistical Software. Vol. 8, No. 15
"""
function effects!(reference_grid::DataFrame, model::RegressionModel;
                  eff_col=nothing, err_col=:err, typical=mean)
    # right now this is written for a RegressionModel and implicitly assumes
    # no link function and the existence of an appropriate formula method
    form = formula(model)
    form_typical = typify(reference_grid, form, modelmatrix(model); typical=typical)
    X = modelcols(form_typical, reference_grid)
    eff = X * coef(model)
    err = sqrt.(diag(X * vcov(model) * X'))
    reference_grid[!, something(eff_col, form.lhs.sym)] = eff
    reference_grid[!, err_col] = err
    return reference_grid
    # XXX remove DataFrames dependency
    # this doesn't work for a DataFrame and isn't mutating
    # return (; reference_grid..., depvar => eff, err_col => err)
end

function _reference_grid(design)
    colnames = tuple(keys(design)...)
    rowtab = NamedTuple{colnames}.(product(values(design)...))

    return DataFrame(vec(rowtab))
end

"""
    effects(design::AbstractDict, model::RegressionModel;
            eff_col=nothing, err_col=:err, typical=mean,
            lower_col=:lower, upper_col=:upper)

Compute the `effects` as specified by the `design`.

This is a convenience wrapper for [`effects!`](@ref). Instead of specifying a
reference grid, a dictionary containing the levels/values of each predictor
is specified. This is then expanded into a reference grid representing a
fully-crossed design. Additionally, two extra columns are created representing
the lower and upper edge of the error band (i.e. [resp-err, resp+err]).
"""
function effects(design::AbstractDict, model::RegressionModel;
                 eff_col=nothing, err_col=:err, typical=mean,
                 lower_col=:lower, upper_col=:upper)
    grid = _reference_grid(design)
    dv = something(eff_col, formula(model).lhs.sym)
    effects!(grid, model; eff_col=dv, err_col=err_col, typical=typical)
    # XXX DataFrames dependency
    grid[!, lower_col] = grid[!, dv] - grid[!, err_col]
    grid[!, upper_col] = grid[!, dv] + grid[!, err_col]
    return grid
    # up_low = let dv = getproperty(reference_grid, dv), err = getproperty(reference_grid, err_col)
    #     (; lower_col => dv .- err, upper_col => dv .+ err)
    # end
    # return (; reference_grid..., up_low...)
end


text:
using SDFResultViewer
using Unitful
using Test

@testset "SDFResultViewer.jl" begin
    λ = 800u"nm"
    dir = pwd()
    file = "0002.sdf"
    times = get_times(dir, λ)
    @test length(times) == 1

    m_Lx = mean_Lx(dir, λ)
    @test length(m_Lx) == 1

    Lx, y, z = Lx_slice(joinpath(dir, file), λ, 0)
    @test length(Lx) > 0
    @test length(Lx) == length(y) == length(z)
end


text:
module Memoize
import ..OmegaCore: Var
using ..Traits: traits, trait
using ..Tagging: Mem, tag
using ..Basis:AbstractΩ
export mem, defcache
# # Memoization

"Default cache"
defcache() = Dict{Any, Any}()

@inline tagmem(ω, cache = defcache()) = tag(ω, (mem = cache,))

# cache traits
cache(f, ω) = true

@inline function memapl(traits::T, f::F, ω::Ω) where {T, F, Ω}
  if cache(f, ω)
    result = get(ω.tags.mem, f, 0)
    if result === 0
      res = Var.prepostapply(traits, f, ω)
      ω.tags.mem[f] = res
      return res
    else
      # ::(Core.Compiler).return_type(f, Tuple{Ω}) # This seems to work too but i don't know why
      return ω.tags.mem[f]::(Core.Compiler).return_type(Var.prepostapply, Tuple{T, F, Ω})
    end
  else
    return Var.prepostapply(traits, f, ω)
  end
end

@inline Var.dispatch(traits::trait(Mem), f, ω::AbstractΩ) = memapl(traits, f, ω)

"""
`mem(x)`

A 'memoized' version of `x` such that the result of calls to `x(ω)` are cached.

```julia
using OmegaCore, LinearAlgebra, BenchmarkTools
x = dimsnth(StdNormal(), (1000, 1000))
h(x) = (println("call!"); svd(x).S)
y(ω) = h(x(ω))
const y_ = Variable(y)
vars(ω) = (y_(ω), y_(ω)*10, y_(ω)*20)
@benchmark vars(ω) setup = (ω = randω())
@benchmark mem(vars)(ω) setup = (ω = randω())
```
"""
@inline mem(f, cache = defcache()) = ω -> f(tagmem(ω, cache))

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using SymPy
using Test

include("tests.jl")
include("test-math.jl")
include("test-matrix.jl")
include("test-ode.jl")
include("test-logical.jl")
#include("test-specialfuncs.jl")  ## XXX NEEDS WORK
#include("test-permutations.jl")
include("test-physics.jl")
include("test-external-module.jl")


text:
mutable struct EquilibrationData{T <: Real, S} <: PreconditionerDataK2{T, S}
  P::Diagonal{T, S}
  C_equi::Diagonal{T, S}
end

function Equilibration(
  id::QM_IntData,
  fd::QM_FloatData{T},
  regu::Regularization{T},
  D::AbstractVector{T},
  K::Union{LinearOperator{T}, AbstractMatrix{T}},
) where {T <: Real}
  P = Diagonal(similar(D, id.nvar + id.ncon))
  P.diag .= one(T)
  C_equi = Diagonal(similar(D, id.nvar + id.ncon))
  return EquilibrationData(P, C_equi)
end

function update_preconditioner!(
  pdat::EquilibrationData{T},
  pad::PreallocatedData{T},
  itd::IterData{T},
  pt::Point{T},
  id::QM_IntData,
  fd::QM_FloatData{T},
  cnts::Counters,
) where {T <: Real}
  equilibrate_K2!(
    fd.Q.data,
    fd.A,
    pad.D,
    pad.regu.δ,
    id.nvar,
    id.ncon,
    pad.pdat.P,
    pad.pdat.C_equi,
    fd.uplo;
    ϵ = T(1.0e-4),
    max_iter = 100,
  )
  pdat.P.diag .= pdat.P.diag .^ 2
end


text:
@testset "LockFreeRingBuffer Tests" begin
    @testset "basic read/write" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        result = zeros(Int32, 10)
        @test write(buf, data) == 10
        @test read!(buf, result) == 10
        @test result == data
        close(buf)
    end

    @testset "read/write that wraps the ringbuf" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        result = zeros(Int32, 10)
        write(buf, data)
        read!(buf, result)
        data = rand(Int32, 10)
        @test write(buf, data) == 10
        @test read!(buf, result) == 10
        @test result == data
        close(buf)
    end

    @testset "overflow" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 20)
        result = zeros(Int32, 16)
        @test write(buf, data) == 16
        @test read!(buf, result) == 16
        @test result == data[1:16]
        close(buf)
    end

    @testset "underflow" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        result = zeros(Int32, 16)
        @test write(buf, data) == 10
        @test read!(buf, result) == 10
        @test result[1:10] == data
        close(buf)
    end

    @testset "writing raw array doesn't allocate" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        dataptr = pointer(data)
        # call once to trigger compilation
        write(buf, dataptr, length(data))
        alloc = @allocated write(buf, dataptr, length(data))
        @test alloc == 0
        close(buf)
    end

    @testset "reading raw array doesn't allocate" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        result = zeros(Int32, 10)
        resultptr = pointer(result)
        # call once to trigger compilation
        write(buf, data)
        read!(buf, resultptr, length(result))
        write(buf, data)
        alloc = @allocated read!(buf, resultptr, length(result))
        @test alloc == 0
        close(buf)
    end

    @testset "wake up waiters when written to" begin
        buf = LockFreeRingBuffer(Int32, 16)
        startcond = Condition()
        endcond = Condition()
        data = rand(Int32, 10)
        woke = false
        @schedule begin
            notify(startcond)
            wait(buf)
            woke = true
            notify(endcond)
        end
        wait(startcond)
        @test !woke
        write(buf, data)
        wait(endcond)
        # if we get here then the waiter was successfully woken
        @test woke
        close(buf)
    end

    @testset "wake up waiters when read from" begin
        buf = LockFreeRingBuffer(Int32, 16)
        startcond = Condition()
        endcond = Condition()
        data = rand(Int32, 10)
        result = zeros(Int32, 10)
        write(buf, data)
        wait(buf) # we need to make sure that the internal AsyncConditino gets triggered
        woke = false
        @schedule begin
            notify(startcond)
            wait(buf)
            woke = true
            notify(endcond)
        end
        wait(startcond)
        @test !woke
        read!(buf, result)
        wait(endcond)
        # if we get here then the waiter was successfully woken
        @test woke
        close(buf)
    end

    @testset "supports nreadable and nwritable" begin
        buf = LockFreeRingBuffer(Int32, 16)
        data = rand(Int32, 10)
        result = rand(Int32, 5)
        @test nreadable(buf) == 0
        @test nwritable(buf) == 16
        write(buf, data)
        @test nreadable(buf) == 10
        @test nwritable(buf) == 6
        read!(buf, result)
        @test nreadable(buf) == 5
        @test nwritable(buf) == 11
        close(buf)
    end

    @testset "supports length" begin
        buf = LockFreeRingBuffer(Int32, 16)
        @test length(buf) == 16
    end
end


text:
@testset "Krylov K1" begin
  stats4 = ripqp(
    QuadraticModel(qps4),
    display = false,
    iconf = InputConfig(sp = K1KrylovParams(), solve_method = :IPF, history = true),
    itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
  )
  @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
  @test stats4.status == :acceptable

  stats4 = ripqp(
    QuadraticModel(qps4),
    display = false,
    iconf = InputConfig(sp = K1KrylovParams(uplo = :U), solve_method = :PC, history = true),
    itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
  )
  @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
  @test stats4.status == :acceptable
end


text:
module DelftBlue

using Gridap
using GridapDistributed
using GridapGmsh
using PartitionedArrays
using Gridap.FESpaces: zero_free_values, interpolate!
using Gridap.Fields: meas
using LineSearches: BackTracking
using GridapPETSc
using GridapPETSc: PETSC

include("TGV.jl")
include("Cylinder.jl")

function main(n::Int,order::Int,np::Int,dt::Real,tf::Real;test="TGV",mesh_file="mesh.msh")
  prun(mpi,(np,np,np)) do parts
    options = "-snes_type newtonls -snes_linesearch_type basic  -snes_linesearch_damping 1.0 -snes_rtol 1.0e-12 -snes_atol 0.0 -snes_monitor -ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
    if test=="TGV"
      GridapPETSc.with(args=split(options)) do
        run_TGV(parts,n,order,dt,tf)
      end
    elseif test=="Cylinder"
      GridapPETSc.with(args=split(options)) do
        run_Cylinder(parts,order,dt,tf,mesh_file)
      end
    end
  end
end

end


text:
# readers.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
export read_tweets, read_usenet, filelines, zopen

# -------------------------------------------------------------------------------------------------------------------------
# Basic utilities
# -------------------------------------------------------------------------------------------------------------------------
function zopen(fn :: AbstractString)
  return ismatch(r"^.*\.gz$", fn) ? gzopen(fn) : open(fn)
end

type FileLines
  name :: AbstractString
end
Base.iteratorsize(::Type{FileLines}) = Base.SizeUnknown()

start(itr :: FileLines) = zopen(itr.name)
function done(itr :: FileLines, stream)
  if !eof(stream)
    return false
  end
  close(stream)
  return true
end
function next(itr :: FileLines, stream)
  x = readline(stream, chomp=false)
  return (x, stream)
end
eltype(itr :: FileLines) = ByteString

# get a file line iterator from a file name, open with gzip as needed
filelines(fn :: AbstractString) = FileLines(fn)
streamlines(f) = eachline(f) # convenience

# -------------------------------------------------------------------------------------------------------------------------
# Text format readers
# -------------------------------------------------------------------------------------------------------------------------

# read collection of tweets from a file
function read_tweets(fn :: AbstractString; stopList=english_stoplist, header=news_email_header, tokenizer=tenglish_tokenizer, limit=-1, keepFn=x->true, lg=STDERR)
  ret = Dict{AbstractString, Float32}[]
  rlat = Float32[]
  rlong = Float32[]

  tic()
  f = zopen(fn)
  tweets = JSON.parse(f)
  close(f)

  for tweet in tweets #filelines(fn)
    # json parser is wicked slow
    #tweet = parse(l)
    text = tweet["text"]
    if keepFn(text) # was vec (see below)

      # extract text and location
      # m = match(r"^.*?\"?text\"?\s*:\s*\"(.*?)\".*$", l)
      # if (m == nothing) continue; end
      # text = m.captures[1]
      # text = unescape_string(text)
      valid = true

      # validate text
      for c in text
        if 0xd800 <= c <= 0xdfff || '\U10ffff' < c
          valid = false
        end
      end

      # check for valid geo
      #m = match(r"^.*?\"?geo\"?\s*:\s*(\{.*?\}|null).*$", l)
      #geo = m.captures[1]
      #lat, long = -1000, -1000
      geo = tweet["geo"]
      lat, long = -1000, -1000
      if (geo != nothing && !haskey(geo, "coordinates"))
        #mm = match(r"^{.*\"?coordinates\"?\s*:\s*\[(.*?)\s*,\s*(.*?)\s*\]\s*.*\}$", geo)
        lat  = geo["coordinates"][1]
        long = geo["coordinates"][2]
        #if (m != nothing)
        #  lat, long = mm.captures
        #end
      end

      if (valid)
        text = replace(text, r"\s+"s, " ")
        vec = Dict{AbstractString, Float32}()
        for w in filter(x -> !stopList[x], tokenizer(lowercase(strip(text))))
          if haskey(vec, w)
            vec[w] += 1
          else
            vec[w] = 1
          end
        end

        push!(ret, vec)
        push!(rlat, float32(lat))
        push!(rlong, float32(long))
      end
      if (limit > 0 && size(ret, 1) > limit)
        break
      end
    end
  end
  t = toq()
  @info lg @sprintf("finished reading %60s [n = %6d, time = %10.3f]", fn, length(ret), t)
  return ret, rlat, rlong
end

# usenet/email single document reader -- 20ng
function read_usenet(fn :: AbstractString; stopList=english_stoplist, header=news_email_header, tokenizer=english_tokenizer, lg=STDERR)
  ignore = false
  @info lg @sprintf("reading: %s", fn)
  vec = Dict{AbstractString, Float32}()
  for l in eachline(`iconv -f latin1 -t utf8 $fn`)
    cl = lowercase(strip(l))
    if (ismatch(r"--+", cl)) ignore = true end
    if (!ignore && cl != "" && !ismatch(r"^((jj)?>|:|\$).*$", cl) && !ismatch(header, cl))
      for w in tokenizer(cl) #split(cl, r"\s+")
        if !stopList[w]
          if (has(vec, w))
            vec[w] += 1
          else
            vec[w] = 1
          end
        end
      end
    end
  end
  return [ vec ]
end


text:
module LongMultiplication

using Compat

function addwithcarry!(r, addend, addendpos)
    while true
        pad = max(0, addendpos - lastindex(r))
        append!(r, fill(0, pad))
        addendrst = addend + r[addendpos]
        addend, r[addendpos] = divrem(addendrst, 10)
        iszero(addend) && break
        addendpos += 1
    end
    return r
end

function longmult(mult1::AbstractVector{T}, mult2::AbstractVector{T}) where T <: Integer
    r = T[]
    for (offset1, digit1) in enumerate(mult1), (offset2, digit2) in zip(eachindex(mult2) + offset1 - 1, mult2)
        single_multrst = digits(digit1 * digit2)
        for (addoffset, rstdigit) in zip(eachindex(single_multrst) + offset2 - 1, single_multrst)
            addwithcarry!(r, rstdigit, addoffset)
        end
    end
    return r
end

function longmult(a::T, b::T)::T where T <: Integer
    mult1 = digits(a)
    mult2 = digits(b)
    r = longmult(mult1, mult2)
    return sum(d * T(10) ^ (e - 1) for (e, d) in enumerate(r))
end

function longmult(a::AbstractString, b::AbstractString)
    if !ismatch(r"^\d+", a) || !ismatch(r"^\d+", b)
        throw(ArgumentError("string must contain only digits"))
    end
    mult1 = reverse(collect(Char, a) .- '0')
    mult2 = reverse(collect(Char, b) .- '0')
    r = longmult(mult1, mult2)
    return reverse(join(r))
end

end  # module LongMultiplication

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PhaseRetrieval
using Plots, Images
default(aspect_ratio =:equal)

# repeating my simulations in MMA
dim =520 # total grid size for padding
s=5.2
p=150
# cell=p/s
cell = 64
d=1.2
r=cell*8/dim # aperture of 512 pixels, relative size

ph1 = SHdiversity(dim, cell, (12,-8.5), .1); # offset and rotation are chosen arbitrary
ph1 = 2*pi*ph1./(4* dim)

heatmap(phwrap.(ph1), c=:cyclic_mygbm_30_95_c78_n256, size=(dim,dim)) |> display

ap1 = [ i^2 + j^2 < r^2 for i in range(-1,1,length= dim), j in range(-1,1, length=dim)]
display(heatmap(ap1))
hp1ref = psf(ap1, ph1)
heatmap(hp1ref, c=:grays, size = (540,540), cb=:none)

# using Images instead of heatmap
Gray.(PhaseRetrieval.rescale(phwrap.(ph1)))
Gray.(PhaseRetrieval.rescale(hp1ref))

# for 4 times more lenslets (32 x 32) we need upsample it 4 times
dim=1080
q = 4
dimup = dim*q
phup = SHdiversity(dimup, cell, (12,-8.5), .1); # offset and rotation are chosen arbitrary
phup = 2*pi*phup./(2* dimup)
Gray.(PhaseRetrieval.rescale(phwrap.(phup)))
apup = [ i^2 + j^2 < r^2 for i in range(-1,1,length= dimup), j in range(-1,1, length=dimup)]
hpupref = imresize(psf(apup, phup*2), ratio=1/q)
Gray.(PhaseRetrieval.rescale(hpupref))
hp2=restrict(restrict(psf(apup, -phup*2)))
Gray.(PhaseRetrieval.rescale(hp2))

# in MMA I had dimgrid 5850, q = 5, ap= 4805, number of pixel per s/a is 145
dim=1170
q = 5
dimup = dim*q
r= 961/1170
cellup =145
phup = SHdiversity(dimup, cellup, (12,-8.5), .1); # offset and rotation are chosen arbitrary
phup = 2*pi*phup./(2* dimup)
Gray.(PhaseRetrieval.rescale(phwrap.(phup)))
apup = [ i^2 + j^2 < r^2 for i in range(-1,1,length= dimup), j in range(-1,1, length=dimup)]
hpupref = imresize(psf(apup, phup), ratio=1/q)
Gray.(PhaseRetrieval.rescale(hpupref))
save("p150f10_ref.png",Gray.(PhaseRetrieval.rescale(hpupref)))


text:
__precompile__(true)

module YAML

import Base: isempty, length, show

import Base: iterate

using Base64: base64decode
using Dates
using Printf

include("scanner.jl")
include("parser.jl")
include("composer.jl")
include("constructor.jl")
include("writer.jl") # write Julia dictionaries to YAML files

const _constructor = Union{Nothing, Dict}
const _dicttype = Union{Type,Function}

load(ts::TokenStream, more_constructors::_constructor=nothing; dicttype::_dicttype=Dict{Any, Any}) =
    construct_document(Constructor(_patch_constructors(more_constructors, dicttype)), compose(EventStream(ts)))

# add a dicttype-aware version of construct_mapping to the constructors
function _patch_constructors(more_constructors::_constructor, dicttype::_dicttype)
    if more_constructors == nothing
        more_constructors = Dict{String,Function}()
    else
        more_constructors = copy(more_constructors) # do not change the outside world
    end
    if !haskey(more_constructors, "tag:yaml.org,2002:map")
        more_constructors["tag:yaml.org,2002:map"] = custom_mapping(dicttype) # map to the custom type
    elseif dicttype != Dict{Any,Any} # only warn if another type has explicitly been set
        @warn "dicttype=$dicttype has no effect because more_constructors has the key \"tag:yaml.org,2002:map\""
    end
    return more_constructors
end

load(input::IO, more_constructors::_constructor=nothing; kwargs...) =
    load(TokenStream(input), more_constructors; kwargs...)

mutable struct YAMLDocIterator
    input::IO
    ts::TokenStream
    more_constructors::_constructor
    next_doc

    function YAMLDocIterator(input::IO, more_constructors::_constructor=nothing; dicttype::_dicttype=Dict{Any, Any})
        it = new(input, TokenStream(input), _patch_constructors(more_constructors, dicttype), nothing)
        it.next_doc = eof(it.input) ? nothing : load(it.ts, it.more_constructors)
        return it
    end
end

# Old iteration protocol:
start(it::YAMLDocIterator) = nothing

function next(it::YAMLDocIterator, state)
    doc = it.next_doc
    if eof(it.input)
        it.next_doc = nothing
    else
        reset!(it.ts)
        it.next_doc = load(it.ts, it.more_constructors)
    end
    return doc, nothing
end

done(it::YAMLDocIterator, state) = it.next_doc === nothing

# 0.7 iteration protocol:
iterate(it::YAMLDocIterator) = next(it, start(it))
iterate(it::YAMLDocIterator, s) = done(it, s) ? nothing : next(it, s)

load_all(input::IO, more_constructors::_constructor=nothing; kwargs...) =
    YAMLDocIterator(input, more_constructors; kwargs...)

load(input::AbstractString, more_constructors::_constructor=nothing; kwargs...) =
    load(IOBuffer(input), more_constructors; kwargs...)

load_all(input::AbstractString, more_constructors::_constructor=nothing; kwargs...) =
    load_all(IOBuffer(input), more_constructors; kwargs...)

load_file(filename::AbstractString, more_constructors::_constructor=nothing; kwargs...) =
    open(filename, "r") do input
        load(input, more_constructors; kwargs...)
    end

load_all_file(filename::AbstractString, more_constructors::_constructor=nothing; kwargs...) =
    open(filename, "r") do input
        load_all(input, more_constructors; kwargs...)
    end

end  # module


text:
@testset "1D-Lagrange-Mass" begin
    domain, _ = line(5)
    pfunc = outer(mkbasis(domain, Lagrange, 1))

    func = compile(pfunc)
    mass = integrate(func, domain, 2)
    @test mass ≈ diagm(
        -1 => fill(1/6, (5,)),
        0 => vcat([1/3], fill(2/3, (4,)), [1/3]),
        1 => fill(1/6, (5,)),
    )

    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 2)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 16
    @test mass ≈ sparse(
        [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6],
        [1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6],
        [2, 1, 1, 4, 1, 1, 4, 1, 1, 4, 1, 1, 4, 1, 1, 2] / 6,
    )
end

@testset "1D-Lagrange-Lapl-Ref" begin
    domain, _ = line(5)

    pfunc = outer(grad(mkbasis(domain, Lagrange, 1), 1)[:,1])
    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 1)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 16
    @test mass ≈ sparse(
        [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6],
        [1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6],
        [1, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 1],
    )
end

@testset "1D-Lagrange-Lapl" begin
    domain, geom = line(5)

    pfunc = outer(grad(mkbasis(domain, Lagrange, 1), geom)[:,1])
    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 1)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 16
    @test mass ≈ sparse(
        [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6],
        [1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6],
        [1, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 1],
    )

    pfunc = outer(grad(mkbasis(domain, Lagrange, 1), 2 .* geom)[:,1])
    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 1)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 16
    @test mass ≈ sparse(
        [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6],
        [1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6],
        [1, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 1] / 4,
    )

    pfunc = outer(grad(mkbasis(domain, Lagrange, 1), 0.5 .* geom)[:,1])
    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 1)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 16
    @test mass ≈ sparse(
        [1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6],
        [1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6],
        [1, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 2, -1, -1, 1] * 4,
    )
end

@testset "1D-Load" begin
    domain, _ = line(5)

    func = compile(mkbasis(domain, Lagrange, 1); dense=false)
    load = integrate(func, domain, 1)
    @test isa(load, Vector{Float64})
    @test load ≈ [0.5, 1.0, 1.0, 1.0, 1.0, 0.5]
end

@testset "2D-Lagrange-Mass" begin
    domain, _ = rectilinear(2, 2)
    pfunc = outer(mkbasis(domain, Lagrange, 1))

    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 3)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 49
    @test Matrix(mass) ≈ sparse(
        [1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
         6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9],
        [1, 2, 4, 5, 1, 2, 3, 4, 5, 6, 2, 3, 5, 6, 1, 2, 4, 5, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9,
         2, 3, 5, 6, 8, 9, 4, 5, 7, 8, 4, 5, 6, 7, 8, 9, 5, 6, 8, 9],
        [4, 2, 2, 1, 2, 8, 2, 1, 4, 1, 2, 4, 1, 2, 2, 1, 8, 4, 2, 1, 1, 4, 1, 4, 16, 4, 1, 4, 1,
         1, 2, 4, 8, 1, 2, 2, 1, 4, 2, 1, 4, 1, 2, 8, 2, 1, 2, 2, 4] / 36,
    )
end

@testset "2D-Lagrange-Lapl" begin
    domain, geom = rectilinear(2, 2)
    pfunc = dropdims(sum(outer(grad(mkbasis(domain, Lagrange, 1), geom)), (3,)), (3,))

    func = compile(pfunc; dense=false)
    mass = integrate(func, domain, 3)
    @test isa(mass, SparseMatrixCSC)
    @test nnz(mass) == 49
    @test mass ≈ sparse(
        [1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
         6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9],
        [1, 2, 4, 5, 1, 2, 3, 4, 5, 6, 2, 3, 5, 6, 1, 2, 4, 5, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9,
         2, 3, 5, 6, 8, 9, 4, 5, 7, 8, 4, 5, 6, 7, 8, 9, 5, 6, 8, 9],
        [4, -1, -1, -2, -1, 8, -1, -2, -2, -2, -1, 4, -2, -1, -1, -2, 8, -2, -1, -2, -2, -2, -2, -2, 16, -2, -2,
         -2, -2, -2, -1, -2, 8, -2, -1, -1, -2, 4, -1, -2, -2, -2, -1, 8, -1, -2, -1, -1, 4] / 6,
    )
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
include((@__DIR__)*"/../preamble.jl")

function map_to_t(d0,B,approx_type,t,h=0.005,limit=false,kwargs...)
    if (approx_type==:dg)&&(limit)
        dt = integrate_time(d0,B,t,StableRK4(h); limiter=GeneralisedMUSCL)
    elseif (approx_type==:dg)&&(!limit)
        dt = integrate_time(d0,B,t,StableRK4(h); limiter=NoLimiter)
    else
        dt = integrate_time(d0,B,t,StableRK4(h))
    end
    return dt
end

point_mass_initial_condition(n,rng) = (X = fill(eps(),n), φ = fill(1,n))
d0_map_point_mass(dq) = interior_point_mass(eps(),1,dq)

function exp_initial_condition(n,rng)
    i = rand(rng,1:2,n)
    X = Vector{Float64}(undef,n)
    for k in 1:n
        x = Inf
        while x>=model.b
            x = -log(rand(rng))
        end
        X[k]=x
    end
    return (X=X,φ=i)
end
d0_map_exp(dq) = SFMDistribution((x,i)->(i∈(1:2))*exp(-x)/(1-exp(-model.b))/2,dq)

function re_run_transient(
    orders,approx_types,transient_time,n_sims,n_err_evals,n_boot,qtiles,
    model,m_str,ic,ic_str,d0_map,rng,
)

    # load model and parameters
    model_str = (@__DIR__)*"/"*m_str

    # make folders to write to
    pth = model_str*"/transient_distribution/"*ic_str
    mkpath(pth*"/data")
    mkpath(pth*"/figs")

    # construct all the approximations
    models = make_approximations(orders,approx_types,d0_map,map_to_t)
    sims = simulate_model(model,fixed_time(transient_time),n_sims,ic,pth,rng)

    # evaluate models
    _phases = 1:4
    grid = range(0.0,sims.model.b,length=n_err_evals)
    nodes = 0:1:10#models[orders[1]]["dg1"]["dq"].mesh.nodes
    centres = (nodes[1:end-1] + nodes[2:end])./2
    # @evaluate!(models,(cdf,pdf),grid,_phases)
    evaluate_cell_probs!(models,nodes,_phases)

    # evlauate sims
    truth = @evaluate(sims,(cdf,),grid,_phases)
    truth["cell_probs_evaluated"] = cell_probs(sims,nodes).(centres[:],_phases[:]')
    truth = merge!(truth,bootstrap_cis(sims,n_boot,grid,nodes,_phases,qtiles,rng))

    # parameters for plotting
    _ylims = (-0.01,0.65)#(-0.01,0.4)
    _yticks = 0:0.1:0.65#(0:0.1:0.4,
    _xlims = (-eps(),model.b+eps())
    _xticks = 0:5:10
    make_df_plots!(models,truth,"cdf",grid,_phases,pth*"/figs/cdfs",_xlims,_ylims,_xticks,_yticks)
    _ylims = (-0.01,0.4)
    _yticks = 0:0.1:0.4
    make_df_plots!(models,truth,"pdf",grid,_phases,pth*"/figs/pdfs",_xlims,_ylims,_xticks,_yticks)

    # analysis on cdfs and cell_probs only
    err_funs = ((kolmogorov_smirnov,"cdf"),(L1,"cdf"),(L2,"cdf"),(L1_cell_probs,"cell_probs"))
    titles = [
        "KS error between the true CDF and approximations";
        "L¹ error between the true CDF and approximations";
        "L² error between the true CDF and approximations";
        "L¹ error between the cell masses";
    ]
    analysis(err_funs,models,truth,grid,qtiles,titles,pth;w=true)
    # write all the data for back up
    write_models(models,pth)
    write_truth(truth,pth)
    return nothing
end


text:
## SliceSamplerState holds the internal state ("local variables") of the slice sampler for multivariate parameters

mutable struct SliceSamplerState{F<:VariateForm} <: MCSamplerState{F}
  lstate::ParameterState{Continuous, F}
  rstate::ParameterState{Continuous, F}
  primestate::ParameterState{Continuous, F}
  tune::MCTunerState
  loguprime::Real
  runiform::Real
end

SliceSamplerState(
  lstate::ParameterState{Continuous, F},
  rstate::ParameterState{Continuous, F},
  primestate::ParameterState{Continuous, F},
  tune::MCTunerState=BasicMCTune()
) where {F<:VariateForm} =
  SliceSamplerState(lstate, rstate, primestate, tune, NaN, NaN)

### Slice sampler

struct SliceSampler <: MCSampler
  widths::RealVector # Step sizes for initially expanding the slice
  stepout::Bool # Protects against the case of passing in small widths

  function SliceSampler(widths::RealVector, stepout::Bool)
    @assert all(i -> i > 0, widths) "All widths must be positive"
    new(widths, stepout)
  end
end

SliceSampler(widths::RealVector) = SliceSampler(widths, true)

SliceSampler(widths::Real=1., n::Integer=1, stepout::Bool=true) = SliceSampler(fill(widths, n), stepout)

### Initialize slice sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::SliceSampler,
  outopts::Dict
) where F<:VariateForm
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
end

## Initialize SliceSampler state

sampler_state(
  parameter::Parameter{Continuous, F},
  sampler::SliceSampler,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, F},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
) where {F<:VariateForm} =
  SliceSamplerState(
    generate_empty(pstate),
    generate_empty(pstate),
    generate_empty(pstate),
    tuner_state(parameter, sampler, tuner)
  )

## Reset parameter state

function reset!(pstate::MultivariateParameterState, x::RealVector, parameter::MultivariateParameter, sampler::SliceSampler)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

function reset!(
  sstate::SliceSamplerState{F},
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::MCSampler,
  tuner::MCTuner
) where F<:VariateForm
  tune.proposed, tune.totproposed = (0, tuner.period)
end

show(io::IO, sampler::SliceSampler) = print(io, "Slice sampler")


text:
package java.io;

public class ByteArrayInputStream extends InputStream {
    private byte[] bytes;

    public ByteArrayInputStream(byte[] input) {
	bytes = input;
    }

    public int read() {
	return bytes[0];
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text: include("cameras/projective.jl")

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using Gadfly
using CSV
using DataFrames
using Random
using Dates
import Cairo, Fontconfig

filepath = ARGS[1]
data = CSV.read(filepath, DataFrame)

conns = unique(data[!, :connection])

inds = randperm(length(conns))[1:5]

@show someconns = conns[inds]

someconns = [6474, 2264, 2241, 458, 2372]

filter!(row -> row[:connection] in someconns, data)

p = plot(
	data,
	x=:time,
	y=:flow,
	color=:connection,
	Geom.line,
	Guide.xticks(ticks=[DateTime(2021), DateTime(2026), DateTime(2031)]),
	Guide.xlabel("Year"),
	Guide.ylabel("Power flow (MWh)")
)
img = PNG("conn_flow.png", 6inch, 4inch)
draw(img, p)


text:
# BEGIN 0.1 deprecation

@deprecate restrict(A::AbstractArray, region::Vector{Int}) restrict(A, (region...,))

@deprecate meanfinite(A::AbstractArray, region) meanfinite(A; dims=region)

@deprecate minfinite(A; kwargs...) minimum_finite(A; kwargs...)
@deprecate maxfinite(A; kwargs...) maximum_finite(A; kwargs...)
@deprecate maxabsfinite(A; kwargs...) maximum_finite(abs, A; kwargs...)

# These two symbols are exported by previous ImageBase versions and now organized in the
# FiniteDiff submodule.
@deprecate fdiff(args...; kwargs...) ImageBase.FiniteDiff.fdiff(args...; kwargs...)
@deprecate fdiff!(args...; kwargs...) ImageBase.FiniteDiff.fdiff!(args...; kwargs...)

# END 0.1 deprecation


text:
using DataFrames
import FileIO:save, load

# Top-level runner functions

# Top runner
function runner_compare()
        # ------- Uncomment to replicate the paper's simulations ---------------
        # run_compare_multinomial()
        # run_compare_multinomial(nsteps = 10^6, changeprobabilities = [0.001, 0.0001])
        # run_compare_gaussian()
        # run_compare_gaussian(nsteps = 10^6, changeprobabilities = [0.001, 0.0001])
        # ----------------------------------------------------------------------
        run_compare_multinomial(nsteps = 10, changeprobabilities = [0.01])
        run_compare_gaussian(nsteps = 10, changeprobabilities = [0.01])
end

# Runner for Categorical task (Section 2.3.2)
function run_compare_multinomial(;foldername = "multi", filename = "multi_",
            paramDictpath = "./parameterDicts/",
            # nsteps = 10^6,
            nsteps = 10^5,
            # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
            changeprobabilities = [0.1, 0.05, 0.01, 0.005] # With 10^5 steps
            )

    println("Starting: " * Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS"))
    nseeds = 10

    specs = SpecsMultinomial(
                    # ------- Uncomment to replicate the paper's simulations ---
                    # stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01],
                    stochasticities = [2., 1.],
                    # ----------------------------------------------------------
                    changeprobabilities = changeprobabilities,
                    nseeds = nseeds, seedenvpolicy = collect(1:nseeds),
                    seedlearner = collect(1:nseeds) * 99,
                    ns = 5, nsteps = nsteps)
    # parameterDict = Dict(
    # 		(stochasticityvalue1, pcvalue1) => optimalvalueparameter1
    # 		(stochasticityvalue2, pcvalue2) => optimalvalueparameter2
    # 		(stochasticityvalue3, pcvalue3) => optimalvalueparameter3
    #                 ...
    #                 )
    # # parameterDict = Dict((0.1, 0.9) => 0.3), (0.1, 0.99) => 0.33, (0.1, 0.999) => 0.34)
    paramDict = load(joinpath(paramDictpath, "parameterDict_multi_varsmile.jld2"))
    parameterDict_varsmile = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_multi_leaky.jld2"))
    parameterDict_leaky = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_multi_smileOriginal.jld2"))
    parameterDict_smileOriginal = paramDict["parameterDict"]

    learnerslist = [
                   (name = "leaky", type = TLeakyIntegrator, param = (0.,),
                    parameterDict = parameterDict_leaky),
                    (name = "varsmile", type = TVarSmile, param = (0.,0.),
                    parameterDict = parameterDict_varsmile),
                    # ------ Run VarSmiLe with some other parameter values -----
                    # (name = "varsmile", type = TVarSmile, param = (0.3,0.1),
                    # parameterDict = missing),
                    # ----------------------------------------------------------
                    (name = "pf20", type = TParticleFilter, param = (20, 0., 0.),
                    parameterDict = missing),
                    (name = "pf10", type = TParticleFilter, param = (10, 0., 0.),
                     parameterDict = missing),
                    (name = "pf1", type = TParticleFilter, param = (1, 0., 0.),
                     parameterDict = missing),
                    (name = "bayesfilter", type = TBayesFilter, param = (0.,0.),
                    parameterDict = missing),
                    (name = "tmpn1", type = TMPN, param = (1, 0., 0.),
                    parameterDict = missing),
                    (name = "tmpn10", type = TMPN, param = (10, 0., 0.),
                    parameterDict = missing),
                    (name = "tmpn20", type = TMPN, param = (20, 0., 0.),
                    parameterDict = missing),
                    (name = "smileOriginal", type = TSmileOriginal, param = (0., 0.),
                    parameterDict = parameterDict_smileOriginal),
                    (name = "tsor1", type = TSOR, param = (1, 0., 0.),
                    parameterDict = missing),
                    (name = "tsor10", type = TSOR, param = (10, 0., 0.),
                    parameterDict = missing),
                    (name = "tsor20", type = TSOR, param = (20, 0., 0.),
                    parameterDict = missing)
                    ]
    runner_eachlearner(learnerslist, specs, foldername, filename)
end

# Runner for Gaussian task (Section 2.3.1)
function run_compare_gaussian(;foldername = "gaussian", filename = "gaussian_",
                paramDictpath = "./parameterDicts/",
                # nsteps = 10^6,
                nsteps = 10^5,
                # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps
                # ------- Uncomment to replicate the paper's simulations -------
                # sigmas = [5., 2., 1., 0.5, 0.1],
                sigmas = [0.5, 0.1],
                # --------------------------------------------------------------
                )

    println("Starting: " * Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS"))
    nseeds = 10

    specs = SpecsGaussian(; mu_0 = 0., sigma_0 = 1.,
                    sigmas = sigmas,
                    changeprobabilities = changeprobabilities,
                    nseeds = nseeds, seedenvpolicy = collect(1:nseeds),
                    seedlearner =  collect(1:nseeds) * 99,
                    nsteps = nsteps)
    # parameterDict = Dict(
    # 		(stochasticityvalue1, pcvalue1) => optimalvalueparameter1
    # 		(stochasticityvalue2, pcvalue2) => optimalvalueparameter2
    # 		(stochasticityvalue3, pcvalue3) => optimalvalueparameter3
    #                 ...
    #                 )
    # # parameterDict = Dict((0.1, 0.9) => 0.3)#, (0.1, 0.99) => 0.33, (0.1, 0.999) => 0.34)
    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_varsmile.jld2"))
    parameterDict_varsmile = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_leaky.jld2"))
    parameterDict_leaky = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_smileOriginal.jld2"))
    parameterDict_smileOriginal = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_gnassar10.jld2"))
    parameterDict_gnassar10 = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_gnassar12.jld2"))
    parameterDict_gnassar12 = paramDict["parameterDict"]
    #
    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_gnassar10original.jld2"))
    parameterDict_gnassar10original = paramDict["parameterDict"]

    paramDict = load(joinpath(paramDictpath, "parameterDict_gaussian_gnassar12original.jld2"))
    parameterDict_gnassar12original = paramDict["parameterDict"]

    learnerslist = [
                    (name = "leaky", type = GLeakyIntegrator, param = (0.,),
                    parameterDict = parameterDict_leaky),
                    (name = "varsmile", type = GVarSmile, param = (0.,0.),
                    parameterDict = parameterDict_varsmile),
                    # ------ Run VarSmiLe with some other parameter values -----
                    # (name = "varsmile", type = GVarSmile, param = (0.3,0.1),
                    # parameterDict = missing),
                    # ----------------------------------------------------------
                    (name = "pf20", type = GParticleFilter, param = (20, 0., 0.),
                    parameterDict = missing),
                    (name = "pf10", type = GParticleFilter, param = (10, 0., 0.),
                    parameterDict = missing),
                    (name = "pf1", type = GParticleFilter, param = (1, 0., 0.),
                    parameterDict = missing),
                    (name = "bayesfilter", type = GBayesFilter, param = (0.,0.),
                    parameterDict = missing),
                    (name = "gmpn1", type = GMPN, param = (1, 0., 0.),
                    parameterDict = missing),
                    (name = "gmpn10", type = GMPN, param = (10, 0., 0.),
                    parameterDict = missing),
                    (name = "gmpn20", type = GMPN, param = (20, 0., 0.),
                    parameterDict = missing),
                    (name = "gsor1", type = GSOR, param = (1, 0., 0.),
                    parameterDict = missing),
                    (name = "gsor10", type = GSOR, param = (10, 0., 0.),
                    parameterDict = missing),
                    (name = "gsor20", type = GSOR, param = (20, 0., 0.),
                    parameterDict = missing),
                    (name = "smileOriginal", type = GSmileOriginal, param = (0., 0.),
                    parameterDict = parameterDict_smileOriginal),
                    (name = "gnassar10", type = GNas10, param = (0.,0.),
                    parameterDict = parameterDict_gnassar10),
                    (name = "gnassar12", type = GNas12, param = (0.,0.),
                    parameterDict = parameterDict_gnassar12),
                    # (name = "gnassar10original", type = GNassar10original, param = (0.,0.),
                    # parameterDict = parameterDict_gnassar10original),
                    # (name = "gnassar12original", type = GNassar12original, param = (0.,0.),
                    #  parameterDict = parameterDict_gnassar12original)
                    ]
    runner_eachlearner(learnerslist, specs, foldername, filename)
end


text:
"""
    hyperbola2(x, s1, s2, Δy, xoff=0.0, yoff=0.0)

One branch of a hyperbola as a function of x given by

- s1: asymptote slope for x<0
- s2: asymptote slope for x>0
- Δy: absolute value of abcissa at x=0
- optionally shift the origin by xoff, yoff

Ref https://maurow.bitbucket.io/docs/hyperbola.html
"""
function hyperbola(x, s1, s2, Δy, xoff=0.0, yoff=0.0)
    x = x-xoff
    # Deal with degenerate cases and select the right branch of the hyperbola
    if Δy==0 # return two lines
        i1 = x.<0
        i2 = x.>=0
        return vcat(s1*x[i1], s2*x[i2]) + yoff
    elseif s1==s2 # return one line
        branch = 0
    elseif s1>s2
        branch = -1
        s1,s2 = s2,s1
    else
        branch = +1
    end
    B = (s1+s2)/2
    C = -Δy^2
    D = - (s2-s1)^2/4
    return B*x + branch * sqrt(-D*x.^2-C) + yoff
end

"""
    hyperbola01(x, Δy, xoff=0.0, yoff=0.0)

A hyperbola with asymptotes y=0 for x<0 and y=x for x>0, and
abcissa Δy.  Its origin can be moved with xoff and yoff.
"""
function hyperbola01(x, Δy, xoff=0.0, yoff=0.0)
    x = x-xoff
    B = 1/2
    C = -Δy^2
    D = -1/4
    return B*x + sqrt(-D*x^2-C) + yoff
end

"""
    max_smooth(x0, x, delta)

to replace

    max(x0, x)

for a fixed x0 and varying x.  `delta` gives the difference:
delta = max_smooth(x0,x0,delta)-x0,
which is the maximal difference for any x.

Note, all returned values >= x0 (== for delta==0).
"""
max_smooth(x0, x, delta) =
    delta==0 ? max(x0,x) : hyperbola01(x, delta, x0, x0)
"""
    min_smooth(x0, x, delta)

see max_smooth
"""
min_smooth(x0, x, delta) = -max_smooth(-x0, -x, delta)

####
# Sigmoid smoothing

"""
    sigmoid(x, x0, w)

Smooth transition function from 0 to 1.  At x0+w its value is about 0.99,
at x0-w about 0.01.
"""
sigmoid(x, x0, w) = 1/(1+exp(-(x-x0)/w*5))

"""
    fn_cap(x0, width, fn, x, args...)

Instead of

    max(fn(x0,args...), fn(x, args...))

where at x0+/-width the value is to 1% correct.

There is also the not recommended:

    fn_cap_(f0, width, fn, args...)

where at f0+/-width the value is to 1% correct.
"""
function fn_cap(x0, width, fn, x, args...)
    ff = fn(x, args...)
    f0 = fn(x0, args...)
    si = sigmoid(x, x0, width)
    f0*(1-si) + ff*si
end
function fn_cap_(f0, width, fn, args...)
    ff = fn.(args...)
    si = sigmoid.(ff, f0, width)
    f0.*(1-si) + ff.*si
end

"""
    val_cap(x0,x,width)
instead of
    max(x0, x)

with the return value within 1% for x=x0+/-width.
"""
function val_cap(x0,x,width)
    si = sigmoid(x,x0,width)
    x0*(1-si)+x*si
end


text:
using PNG
using FileIO: DataFormat, @format_str, Stream, File, filename, stream
using Test
using FixedPointNumbers
using ImageCore
using Logging
using Random

#logger = ConsoleLogger(stdout, Logging.Debug)
logger = ConsoleLogger(stdout, Logging.Info)
global_logger(logger)

tmpdir = joinpath(@__DIR__,"temp")
@testset "PNG" begin
    isdir(tmpdir) && rm(tmpdir, recursive = true)
    mkdir(tmpdir)

    img = rand(Bool, 5, 5, 5, 5)
    filepath = joinpath(tmpdir, "5x5x5x5.png")
    @test_throws ErrorException PNG.save(File{DataFormat{:PNG}}(filepath), img)

    @testset "Binary Image" begin
        a = rand(Bool, 11, 10)
        filepath = joinpath(tmpdir, "binary1.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), a)
        b1 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test b1 == convert(Array{Gray{N0f8}}, a)

        a = bitrand(5,5)
        filepath = joinpath(tmpdir, "binary2.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), a)
        b2 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test b2 == convert(Array{Gray{N0f8}}, a)

        a = colorview(Gray, a)
        filepath = joinpath(tmpdir, "binary3.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), a)
        b3 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test b3 == convert(Array{Gray{N0f8}}, a)
    end

    @testset "Gray image" begin
        gray = vcat(fill(Gray(1.0), 10, 10), fill(Gray(0.5), 10, 10), fill(Gray(0.0), 10, 10))
        filepath = joinpath(tmpdir, "gray1.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), gray)
        g1 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test g1 == convert(Array{Gray{N0f8}}, gray)

        gray = rand(Gray{N0f8}, 5, 5)
        filepath = joinpath(tmpdir, "gray2.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), gray)
        g2 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test g2 == gray
    end

    @testset "Color - RGB" begin
        #rgb8 = rand(RGB{N0f8}, 10, 5)
        rgb8 = reshape(range(RGB{N0f8}(1,0,0),RGB{N0f8}(0,1,1), length=10*5), 10, 5)
        filepath = joinpath(tmpdir, "rgb_n0f8.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), rgb8)
        r1 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test r1 == rgb8

        #rgb16 = rand(RGB{N0f16}, 10, 5)
        rgb16 = reshape(range(RGB{N0f16}(1,0,0),RGB{N0f16}(0,1,1), length=10*5), 10, 5)
        filepath = joinpath(tmpdir, "rgb_n0f16.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), rgb16)
        r2 = PNG.load(File{DataFormat{:PNG}}(filepath))
        PNG.save(File{DataFormat{:PNG}}(joinpath(tmpdir, "rgb_n0f16_resave.png")), r2)
        @test r2 == rgb16
    end

    @testset "Alpha" begin
        # RGBA
        r = RGBA(1.0,0.0,0.0, 0.2)
        g = RGBA(0.0,1.0,0.0, 0.8)
        b = RGBA(0.0,0.0,1.0, 1.0)
        rgba1 = vcat(fill(r, 50,100), fill(g, 50,100), fill(b, 50,100))
        filepath = joinpath(tmpdir, "rgba1.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), rgba1)
        r1 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test r1 == rgba1

        # GrayA
        r = GrayA(1.0, 0.25)
        g = GrayA(0.5, 0.5)
        b = GrayA(0.0, 0.75)
        graya = vcat(fill(r, 50,100), fill(g, 50,100), fill(b, 50,100))
        filepath = joinpath(tmpdir, "graya1.png")
        PNG.save(File{DataFormat{:PNG}}(filepath), graya)
        g1 = PNG.load(File{DataFormat{:PNG}}(filepath))
        @test g1 == convert(Array{GrayA{N0f8}}, graya)
    end
    # TODO implement palette
end


text:
abstract type HMCMetric end

struct DiagEuclideanMetric <: HMCMetric end

function ahmc_metric(metric::DiagEuclideanMetric, dim::Integer)
    return AdvancedHMC.DiagEuclideanMetric(dim)
end

struct UnitEuclideanMetric <: HMCMetric end

function ahmc_metric(metric::UnitEuclideanMetric, dim::Integer)
    return AdvancedHMC.UnitEuclideanMetric(dim)
end

struct DenseEuclideanMetric <: HMCMetric end

function ahmc_metric(metric::DenseEuclideanMetric, dim::Integer)
    return AdvancedHMC.DenseEuclideanMetric(dim)
end


text:
# output functions

@doc """
### EulerEquationMod.printSolution

  This function prints the solution vector u to a file named solution.dat.
  Each line contains: element number, local node number, and the 4 conservative
  variable values.

  Inputs:
   * mesh :: PumiMesh2
   * u  : vector of values
"""->
function printSolution(mesh::AbstractMesh, u::AbstractVector)
# print solution to file
# format = for each element, node, print u rho*u rho*v E

println("entered printSolution")
f = open("solution.dat", "w")

for i=1:mesh.numEl
  dofnums_i = getGlobalNodeNumbers(mesh, i)

  for j=1:3
    u_vals = u[dofnums_i[:, j]]
    str = @sprintf("%d %d %16.15e %16.15e %16.15e %16.15e \n", i, j, u_vals[1], u_vals[2], u_vals[3], u_vals[4] )  # print element number

    print(f, str)
  end
#  print(f, "\n")
end

close(f)
return nothing

end

@doc """
### EulerEquationMod.printCoordinates

  This function prints the mesh vertex coordinates in the same format
  as printSolution prints the solution vector.

  Inputs:
    * mesh :: PumiMesh2
"""->
function printCoordinates(mesh::AbstractMesh)
# print solution to file
# format = for each element, node, print u rho*u rho*v E

myrank = mesh.myrank
println("entered printCoordinates")
writedlm("coords_output_$myrank.dat", mesh.coords)

return nothing

end

@doc """
### EulerEquationMod.printSolution

  This function prints the solution vector to the file with the given name.
  The solution vector is printed one value per line.

  Inputs:
    * name : AbstractString naming file, including extension
    * u  : vector to print
"""->
function printSolution(name::AbstractString, u::AbstractVector)

  f = open(name, "a+")

  for i=1:length(u)
    write(f, string(u[i], "\n"))
  end

  close(f)

  return nothing
end

@doc """
### EulerEquationMod.printMatrix

  This function prints a 2D matrix of Float64s to a file with the given name

  Inputs:
    * name  : AbstractString naming file, including extension
    * u  : matrix to be printed
"""->
function printMatrix(name::AbstractString, u::AbstractArray{Float64, 2})
# print a matrix to a file, in a readable format

(m,n) = size(u)

f = open(name, "a+")

for i=1:m
  for j=1:n
    str = @sprintf("%16.15e", u[i,j])
    print(f, str)
    if j < n  # don't put space at end of line
      print(f, " ")
    end
  end

  print(f, "\n")
end

close(f)
return nothing

end

function printMatrix(name::AbstractString, u::AbstractArray{T, 3}) where T
# print a matrix to a file, in a quasireadable format

#println("printing 3d matrix")
(tmp, m,n) = size(u)
#println("numel = ", n)
#println("nnodes = ", m)

f = open(name, "a+")

for i=1:n  # loop over elements
  for j=1:m  # loop over nodes on element
#    str = @sprintf("%16.15e", u[i,j])
    println(f, "el ", i, " node ", j, " flux = ", u[:, j, i] )
#    if j < n  # don't put space at end of line
#      print(f, " ")
#    end
  end

#  print(f, "\n")
end

close(f)
return nothing

end

# print_qvec_coords
#
# This will print out a table of all dof's, with the below columns. Useful for debugging.
#
# Assumes you have eqn and mesh properly initialized.

function print_qvec_coords(mesh, sbp, eqn, opts; to_file=false, filename="null", other_field=0)

  if to_file
    if filename == "null"
      error("filename has not been specified but print_qvec_coords has been passed the to_file=true arg.")
    end
    fh = open(filename, "w")
    write(fh,"x-coord    y-coord    elnum    nodenum  dofnum   val\n")
  else
    println("x-coord    y-coord    elnum    nodenum  dofnum   val")
  end

  for dof_ix = 1:mesh.numDof
    # st:  subscript_tuple
    st = ind2sub((mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl), dof_ix)

    dofnum = st[1]      # dof number in this node
    nodenum = st[2]     # node number in this element
    elnum = st[3]       # element number

    if getindex(eqn.q, dofnum, nodenum, elnum) != eqn.q_vec[dof_ix]
      error("problem with ind2sub")
    end

    coords_this_dof = getindex(mesh.coords, :, nodenum, elnum)

    x_coord = coords_this_dof[1]
    y_coord = coords_this_dof[2]

    if other_field == 0
      val_this_dof = eqn.q_vec[dof_ix]
    else
      val_this_dof = other_field[dof_ix]      # assuming ordering of other_field is the same as eqn.q_vec
    end

    if to_file
      str = @sprintf("%-8.5f   %-8.5f   %-6d   %-6d   %-6d   %s\n", x_coord, y_coord, elnum, nodenum, dofnum, val_this_dof)
      write(fh, str)
    else
      @printf("%-8.5f   %-8.5f   %-6d   %-6d   %-6d   %s\n", x_coord, y_coord, elnum, nodenum, dofnum, val_this_dof)
    end

  end

  if to_file
    close(fh)
  end

end   # end function print_qvec_coords


text:
using HDF5,JLD,KUnet
isdefined(:xtrn) || (@date @load "zn11oparse1.jld")

d0 = 6f0
c0 = 1f0
g0 = .1f0
nc = size(ytrn,1)
niters=1000
nbatch=100
ntest=10000
net=net1=net2=net3=nothing
KUnet.gpu(true)
for y in (:ytrn, :ydev, :ytst)
    @eval $y=full($y)
end # This speeds up the accuracy fn
#xtrn1=xtrn[:,1:100000]
#ytrn1=ytrn[:,1:100000]

net = Layer[KPerceptron(nc, KUnet.kgauss, [0.12f0])]
for epoch=1:100
    @show epoch
    @date train(net, xtrn, ytrn)
    gc()
    @date zdev = predict(net,xdev)
    @date ztst = predict(net,xtst)
    @show (epoch, size(net[1].s), accuracy(ydev,zdev), accuracy(ytst,ztst))
end

if false

# 100x100 best at g=3 g0=0.125 a=0.8464 s=2358
# 250x100 best at g=3 g0=0.125 a=0.8843 s=4548
# 500x100 best at g=3 g0=0.125 a=0.9063 s=7416
# 500x100 best at g=3.5 g0=0.088 a=0.91 s=7719
# 1000x100 best at g=3 g0=0.125 a=0.9221 s=12633
# 1000x100 best at g=3.5 g0=0.088 a=0.9228 s=12895
# 1000x100 best at g0=0.12 a=0.9226 s=12572
for g0=.11:.005:.13
    println("")
    @show g0 # = 2f0^-g
    net = Layer[KPerceptron(nc, KUnet.kgauss, [g0])]
    @date train(net, xtrn, ytrn; iters=niters, batch=nbatch)
    @show size(net[1].s)
    @date y = predict(net,xdev)
    @show accuracy(ydev,y)
end

    for kernel in (
                   (:kgauss, [g0]),
                   # (:klinear, nothing),
                   # (:kpoly, [c0,d0]),
                   )
        for loc in (
                    :gpu,
                    # :cpu,
                    )
            KUnet.gpu(loc==:gpu)
            for i=1:2
                println("")
                @show (loc, kernel)
                @show net = Layer[KPerceptron(nc, KUnet.(kernel[1]), kernel[2])]
                @date train(net, xtrn, ytrn; iters=niters, batch=nbatch)
                @show size(net[1].s)
                @date y = predict(net,xtst[:,1:ntest])
                @show accuracy(ytst[:,1:ntest],y)
            end
        end
    end

    KUnet.gpu(true)
    for i=1:2
        println("")
        @show net1 = Layer[KPerceptron(nc, KUnet.klinear)]
        @date train(net1, xtrn, ytrn; iters=niters, batch=nbatch)
        @show size(net1[1].s)
        @date y = predict(net1,xtst[:,1:ntest])
        @show accuracy(ytst[:,1:ntest],y)
    end

    for i=1:2
        println("")
        @show net1 = Layer[KPerceptron(nc, KUnet.klinear4)]
        @date train(net1, xtrn, ytrn; iters=niters, batch=nbatch)
        @show size(net1[1].s)
        @date y = predict(net1,xtst[:,1:ntest])
        @show accuracy(ytst[:,1:ntest],y)
    end

    # this shows the difference is due to the algorithm, same is observed on the cpu
    # we get cpu/gpu divergence between 16 and 17: actually it is a difference of algorithms: limit bug? 1838/1843 nsv @17.
    # turns out the data was corrupt (unsorted rowval), the algorithms are ok
    KUnet.gpu(false)
    for i=1:1
        @show net3 = Layer[KPerceptron(nc, KUnet.klinear)]
        @date train(net3, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 1.698705077 seconds (541901976 bytes allocated)
        @show (size(net3[1].s),accuracy(ytst[:,1:1000],predict(net3,xtst[:,1:1000]))) # size(net3[1].s) => (283246,5341)
    end

    KUnet.gpu(false)
    for i=1:1
        @show net4 = Layer[KPerceptron(nc, KUnet.klinear3)]
        @date train(net4, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 5.800463573 seconds (762111972 bytes allocated, 1.68% gc time)
        @show (size(net4[1].s),accuracy(ytst[:,1:1000],predict(net4,xtst[:,1:1000]))) # size(net4[1].s) => (283246,5341)
    end

    for i=1:2
        @show net1 = Layer[KPerceptron(nc, KUnet.klinear0)]
        @date train(net1, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 2.249490525 seconds (984203260 bytes allocated, 19.72% gc time)
        @show size(net1[1].s)                                     # size(net1[1].s) => (283246,5341)
    end

    for i=1:2
        @show net1 = Layer[KPerceptron(nc, KUnet.klinear)]
        @date train(net1, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 2.006203976 seconds (481354216 bytes allocated, 11.61% gc time)
        @show size(net1[1].s)                                     # size(net1[1].s) => (283246,5341)
    end

    for i=1:2
        @show net2 = Layer[KPerceptron(nc, KUnet.kpoly0, [c0,d0])]
        @date train(net2, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 3.646227942 seconds (1101788224 bytes allocated, 12.16% gc time)
        @show size(net2[1].s)                                     # size(net2[1].s) => (283246,5248)
    end

    for i=1:2
        @show net2 = Layer[KPerceptron(nc, KUnet.kpoly1, [c0,d0])]
        @date train(net2, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 3.433380693 seconds (606331904 bytes allocated, 8.84% gc time)
        @show size(net2[1].s)                                     # size(net2[1].s) => (283246,5248)
    end

    for i=1:2
        @show net2 = Layer[KPerceptron(nc, KUnet.kpoly, [c0,d0])]
        @date train(net2, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 3.336944494 seconds (474105748 bytes allocated, 6.95% gc time)
        @show size(net2[1].s)                                     # size(net2[1].s) => (283246,5248)
    end

    for i=1:2
        @show net3 = Layer[KPerceptron(nc, KUnet.kgauss0, [g0])]
        @date train(net3, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 3.335641623 seconds (1467981384 bytes allocated, 21.92% gc time)
        @show size(net3[1].s)                                     # (283246,5039)
    end

    for i=1:2
        @show net3 = Layer[KPerceptron(nc, KUnet.kgauss, [g0])]
        @date train(net3, xtrn, ytrn; iters=niters, batch=nbatch) # elapsed time: 2.678601132 seconds (988978008 bytes allocated, 19.58% gc time)
        @show size(net3[1].s)                                     # (283246,5039)
    end

end # if false


text:
# This script takes an AdaGram model
# and generates a file with all word senses and their nearest neighbors
# ASuMa 2018

#push!(LOAD_PATH, "./src/")

using ArgParse
using AdaGram

s = ArgParseSettings()

@add_arg_table s begin
  "AdaGramFile"
    help = "File with AdaGram model"
    arg_type = AbstractString
    required = true
  "outputFile"
    help = "File to output word senses"
    arg_type = AbstractString
    required = true
  "--neighbors"
    help = "Number or nearest neighbors to print"
    arg_type = Int64
    default = 5
  "--min-prob"
    help = "minimum probability value to print a word sense"
    arg_type = Float64
    default = 0.05
  "--voc-fraction"
    help = "fraction of the vocabulary to print out"
    arg_type = Float64
    default = 1.0
end

args = parse_args(ARGS, s)

vm, dict = load_model(args["AdaGramFile"]);

fo = open(args["outputFile"], "w")

numWords = round(Int,V(vm) * args["voc-fraction"])
letters = "abcdefghjklmnopqrstuvwxyz" # "i" skipped on purpose bc LG dict

for v in 1:numWords
    probs = expected_pi(vm, v)
    for s in 1:T(vm)
        if probs[s] > args["min-prob"] && vm.counts[s, v] >= 1
            nn = nearest_neighbors(vm, dict, vec(vm, v, s), args["neighbors"]; exclude=[(Int32(v), s)])
            tag_main = string(letters[s % length(letters)])
            @printf(fo, "%s@%s\t%.3f\t", dict.id2word[v], tag_main, probs[s])
    	    for neighbor in nn
        		#@printf(fo, "%s@%s@%.3f ", neighbor[1], neighbor[2], neighbor[3])
                tag_neighbor = string(letters[neighbor[2] % length(letters)])
                @printf(fo, "%s@%s\t", neighbor[1], tag_neighbor)
    	    end
            for neighbor in nn
                @printf(fo, "%.3f\t", neighbor[3])
            end
            println(fo)
            for iDim in 1:M(vm)
                @printf(fo, "%s ", vm.In[:,s,v][iDim])
            end
            println(fo)
        end
    end
end

close(fo)

max_stars_repo_name: akolonin/language-learning

max_issues_repo_name: akolonin/language-learning

max_forks_repo_name: akolonin/language-learning


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

"""
    BallSampler(radius, [maxsize])

A method for sampling isolated points from spatial objects using
a ball neighborhood of given `radius`. The maximum size `maxsize`
of the sample can be specified, but is not required.
"""
struct BallSampler{B<:BallNeighborhood} <: AbstractSampler
  ball::B
  maxsize::Union{Int,Nothing}
end

BallSampler(ball::BallNeighborhood) = BallSampler(ball, nothing)

BallSampler(radius::Real, maxsize=nothing) =
  BallSampler(BallNeighborhood(radius), maxsize)

function sample(object::AbstractSpatialObject{T,N}, sampler::BallSampler) where {T,N}
  npts = npoints(object)
  ball = sampler.ball
  size = sampler.maxsize ≠ nothing ? sampler.maxsize : Inf

  # neighborhood search method with ball
  searcher = NeighborhoodSearcher(object, ball)

  # pre-allocate memory for coordinates
  coords = MVector{N,T}(undef)

  locations = Vector{Int}()
  notviewed = trues(npts)
  while length(locations) < size && any(notviewed)
    location = rand(findall(notviewed))
    coordinates!(coords, object, location)

    # neighbors (including the location)
    neighbors = search(coords, searcher)

    push!(locations, location)
    notviewed[neighbors] .= false
  end

  view(object, locations)
end


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

@userplot ExcessPlot

@recipe function f(p::ExcessPlot)
  # get user input
  xs = p.args[1]

  ks = 2:length(xs)
  ξs = meanexcess(xs, ks)

  seriestype --> :path
  xguide --> "number of maxima"
  yguide --> "extreme value index"

  ks, ξs
end


text:
"""
Based on
https://hyperopt.github.io/hyperopt/
"""
module TestQuantLogNormal

using Test
using TreeParzen

# Simple quantlognormal
# Minimize the objective over the space
space = Dict(:x => HP.QuantLogNormal(:x, 0.0, 1.0, 2.0))
best = fmin(params -> params[:x]^2, space, 100)
@test 0 <= best[:x] <= 1

# What if properties are the same?
space = Dict(:x => HP.QuantLogNormal(:x, 0.0, 1.0, 1.0))
fmin(params -> params[:x]^2, space, 100)

end # module TestQuantLogNormal
true


text:
using IcanProj
using NetCDF
using DataFrames
using ProgressMeter
using PyCall
using StatsBase
using CSV
using Statistics

import IcanProj.project_results

function project_results(df::DataFrame, on::Symbol)

    map_values = fill(NaN, 1550, 1195)

    map_values[df[:ind_julia]] = df[on]

    return map_values

end

function plot_map(df, variable, cb_label, cb_unit, title_label, figpath)

    resmap = project_results(df, variable)

    file = joinpath(figpath, String(variable) * ".png")

    cb_text =  cb_label * " " * cb_unit

    py"""

    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()

    #fig.figsize = (8, 12)

    ax.set_title($(title_label))

    img_plot = ax.imshow($(resmap), cmap = 'jet')

    cbar = fig.colorbar(img_plot)

    cbar.set_label($(cb_text))

    ax.grid(linestyle='dotted')

    plt.tick_params(axis='both', left=False, top=False, right=False, bottom=False,
                    labelleft=False, labeltop=False, labelright=False, labelbottom=False)

    fig.savefig($(file), dpi = 200, bbox_inches='tight')

    plt.close('all')

    """

end

function add_statistics(df, variables, spaceres, func, postfix)

    colnames = names(df)

    for v in variables, s in spaceres

        ikeep = occursin.(Ref(Regex("$(v)_cfg.*_$(s)_mean")), String.(colnames))

        tmp1 = convert(Array{Float64}, df[colnames[ikeep]])

        tmp2 = [func(tmp1[i,:]) for i in 1:size(tmp1,1)]

        df[Symbol(v * "_" * s * "_" * postfix)] = tmp2

    end

    return df

end

function add_turbulent_fluxes(df)

    for cfg in 1:32, spaceres in ["1km", "50km"]

        hatmo_name = Symbol("hatmo_cfg$(cfg)_$(spaceres)_mean")

        latmo_name = Symbol("latmo_cfg$(cfg)_$(spaceres)_mean")

        exchng_name = Symbol("exchng_cfg$(cfg)_$(spaceres)_mean")

        df[exchng_name] = df[hatmo_name] + df[latmo_name]

    end

    return df

end

# Load table with averaged model results

df_all = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "table_results.txt"), delim = ",") |> DataFrame

# Add turbulent heat exchange components

df_all = add_turbulent_fluxes(df_all)

# Compute cov between configurations

variables = ["swe", "snowdepth", "rnet"]

spaceres = ["1km", "50km"]

func = variation

postfix = "cov"

df_all = add_statistics(df_all, variables, spaceres, func, postfix)

# Compute std between configurations

variables = ["hatmo", "latmo", "exchng"]

spaceres = ["1km", "50km"]

func = std

postfix = "std"

df_all = add_statistics(df_all, variables, spaceres, func, postfix)

# Compute average between configurations

variables = ["swe"]

spaceres = ["1km", "50km"]

func = mean

postfix = "mean"

df_all = add_statistics(df_all, variables, spaceres, func, postfix)

# Plot maps

figpath = joinpath(dirname(pathof(IcanProj)), "..", "plots", "variability_cfgs")

variables = [:swe_1km_cov,
             :swe_50km_cov,
             :snowdepth_1km_cov,
             :snowdepth_50km_cov,
             :rnet_1km_cov,
             :rnet_50km_cov,
             :hatmo_1km_std,
             :hatmo_50km_std,
             :latmo_1km_std,
             :latmo_50km_std,
             :exchng_1km_std,
             :exchng_50km_std,
             :swe_1km_mean,
             :swe_50km_mean]

cb_labels = ["Coefficient of variation",
             "Coefficient of variation",
             "Coefficient of variation",
             "Coefficient of variation",
             "Coefficient of variation",
             "Coefficient of variation",
             "Standard deviation",
             "Standard deviation",
             "Standard deviation",
             "Standard deviation",
             "Standard deviation",
             "Standard deviation",
             "Average",
             "Average"]

cb_units = ["(-)",
            "(-)",
            "(-)",
            "(-)",
            "(-)",
            "(-)",
            "(W/m2)",
            "(W/m2)",
            "(W/m2)",
            "(W/m2)",
            "(W/m2)",
            "(W/m2)",
            "(mm)",
            "(mm)"]

title_labels = ["SWE - 1km",
                "SWE - 50km",
                "Snowdepth - 1km",
                "Snowdepth - 50km",
                "Net radiation - 1km",
                "Net radiation - 50km",
                "Sensible heat fluxes - 1km",
                "Sensible heat fluxes - 50km",
                "Latent heat fluxes - 1km",
                "Latent heat fluxes - 50km",
                "Turbulent heat fluxes - 1km",
                "Turbulent heat fluxes - 50km",
                "SWE - 1km",
                "SWE - 50km"]

for (v, cl, cu, t) in zip(variables, cb_labels, cb_units, title_labels)

    plot_map(df_all, v, cl, cu, t, figpath)

end


text:
using Documenter, Joseki, HTTP, JSON

makedocs(
    format = :html,
    sitename = "Joseki.jl"
)


text:
using Distributed
addprocs(2)

@everywhere begin
    using MLJ
    using MLJBase
    using Test
    using Random
end

@testset "ensembles" begin
    @test include("ensembles.jl")
end

@testset "matching models to data" begin
    @test include("model_matching.jl")
end

@testset "scitypes" begin
    @test include("scitypes.jl")
end


text:
using KUnet

net = [Conv(5,5,1,20), Bias(20), Relu(), Pool(2),
       Conv(5,5,20,50), Bias(50), Relu(), Pool(2),
       Mmul(500,800), Bias(500), Relu(),
       Mmul(10,500), Bias(10), XentLoss()]

require("mnist.jl")
using MNIST: xtrn, ytrn, xtst, ytst
xtrn2 = reshape(xtrn, 28, 28, 1, size(xtrn, 2))
xtst2 = reshape(xtst, 28, 28, 1, size(xtst, 2))

using KUnet: accuracy
setparam!(net; lr=0.1)
for i=1:100
    train(net, xtrn2, ytrn)
    println((i, accuracy(ytst, predict(net, xtst2)),
                accuracy(ytrn, predict(net, xtrn2))))
end


text:
using WaveletInformation
using StatsBase

r = MersenneTwister(5)
testmat = rand(r, 100, 1024)
testids = sample(r, [1,2], 100)

@time WaveletInformation.mutual_info(testmat,testids)
@time WaveletInformation.mutual_info(testmat,testids)

@show  WaveletInformation.mutual_info(testmat,testids)[1:5]

L = [1; [2^n for n in 1:10]];
@time WaveletInformation.mutual_info_thresh(testmat,testids,5, 95, L)


text:
export NonlinearSystem

struct NLEq
    rhs::Expression
end
function Base.convert(::Type{NLEq}, eq::Equation)
    isequal(eq.lhs, Constant(0)) || return NLEq(eq.rhs - eq.lhs)
    return NLEq(eq.rhs)
end
Base.:(==)(a::NLEq, b::NLEq) = a.rhs == b.rhs
get_args(eq::NLEq) = Expression[eq.rhs]

struct NonlinearSystem <: AbstractSystem
    eqs::Vector{NLEq}
    vs::Vector{Variable}
    ps::Vector{Variable}
end

function NonlinearSystem(eqs)
    vs, ps = extract_elements(eqs, [_is_unknown, _is_known])
    NonlinearSystem(eqs, vs, ps)
end

function calculate_jacobian(sys::NonlinearSystem)
    rhs = [eq.rhs for eq in sys.eqs]
    jac = expand_derivatives.(calculate_jacobian(rhs, sys.vs))
    return jac
end

function generate_jacobian(sys::NonlinearSystem; version::FunctionVersion = ArrayFunction)
    jac = calculate_jacobian(sys)
    return build_function(jac, sys.vs, sys.ps; version = version)
end

function generate_function(sys::NonlinearSystem; version::FunctionVersion = ArrayFunction)
    rhss = [eq.rhs for eq ∈ sys.eqs]
    return build_function(rhss, sys.vs, sys.ps; version = version)
end


text:
const DEFAULT_INTERRUPT_TEXT = """
                       `InterruptException` recieved.

                       Probably you know about this already.
                       """

default_exception_text(exception, backtrace) = """
            :alert: Error occured! :alert:

            ```
            $(sprint(Base.display_error, exception, backtrace))
            ```
            """

"""
    slack_log_exception(f, thread::AbstractSlackThread; interrupt_text=DEFAULT_INTERRUPT_TEXT,
                        exception_text=default_exception_text)
    slack_log_exception(thread::AbstractSlackThread, exception, backtrace;
                        interrupt_text=DEFAULT_INTERRUPT_TEXT,
                        exception_text=default_exception_text)

Log an exception (and stacktrace) to a `AbstractSlackThread` by constructing
a textual message and passing it to the `thread` with
`SlackThreads.send_exception_message`. Can be used with a zero-argument
function `f`, as in

```julia
slack_log_exception(thread) do
    sqrt(-1)
end
```

or given an `exception` and `backtrace` explictly, as in
```julia
try
    sqrt(-1)
catch exception
    slack_log_exception(thread, exception, catch_backtrace())
    # do other error handling...
end
```

## Optional arguments

* Pass `interrupt_text::AbstractString` to customize the message sent when an `InterruptException` is caught. See `SlackThreads.DEFAULT_INTERRUPT_TEXT` for the default text.
* Pass a two-argument function `exception_text` to customize the message shown for a general exception, where the arguments are `exception` and `backtrace`.
"""
function slack_log_exception(f, thread::AbstractSlackThread;
                             interrupt_text=DEFAULT_INTERRUPT_TEXT,
                             exception_text=default_exception_text)
    try
        f()
    catch exception
        slack_log_exception(thread, exception, catch_backtrace(); interrupt_text,
                            exception_text)
        rethrow()
    end
end

function slack_log_exception(thread::AbstractSlackThread, exception, backtrace;
                             interrupt_text=DEFAULT_INTERRUPT_TEXT,
                             exception_text=default_exception_text)
    @maybecatch begin
        msg = exception isa InterruptException ? interrupt_text :
              exception_text(exception, backtrace)
        send_exception_message(thread, msg)
    end "Error when logging exception to Slack."
    return nothing
end

send_exception_message(thread::SlackThread, msg) = send_message(thread, msg)


text:
domains = getalldomains!(client)

@testset "List all Domains" begin
    @test length(domains) == 1
    @test domains[1].name == "example.com"
    @test domains[1].ttl == 1800
    @test occursin("IN SOA ns1.digitalocean.com.", domains[1].zonefile)
end;

domain = getdomain!(client, "example.com")

@testset "Retrieve an existing Domain" begin
    @test domain.name == "example.com"
    @test domain.ttl == 1800
    @test occursin("IN SOA ns1.digitalocean.com.", domain.zonefile)
end;

domain = createdomain!(client; name="example.com", ip_address="1.2.3.4")

@testset "Create a new Domain" begin
    @test domain.name == "example.com"
    @test domain.ttl == 1800
    @test domain.zonefile == nothing
end;

@testset "Delete a Domain" begin
    @test deletedomain!(client, domain)
end;


text:
@testset "SynthGaussian" begin
    img, p, q, Z = generate_surface_data(SynthGaussian(75), 0.5, [0,0,1])
    @test maximum(img) ≈ 0.5
    @test size(img)[1] ≈ 151
    @test img[10,10] ≈ 0.5
    @test size(p)[1] ≈ 151
    @test size(q)[2] ≈ 151
    @test p[32,54] ≈ 0.03755573897731636
    @test p[79,100] ≈ -0.5570279158609365
    @test q[50,77] ≈ -0.5294367047595221
    @test q[42,42] ≈ -0.06767685688560976
    @test Z[10,10] ≈ 7.796845410874991e-5
    @test Z[75,75] ≈ 18.696742446650838
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "finufft"
version = v"2.0.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/flatironinstitute/finufft/archive/v$(version).zip", "fe83dbc22a1d13799e9db48568d8417d3028e7b831800f2a626ea19b4f0461fe")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/finufft*/
# Overwrite LIBSFFT such that we do not require fftw3_threads or fftw3_omp for OMP support. Since the libraries in FFTW_jll already provide for threading, we do not loose anything.
# Make use of the -DFFTW_PLAN_SAFE flag to allow for multiple threads using finufft at the same time.
make lib CFLAGS="-fopenmp -fPIC -O3 -funroll-loops -fcx-limited-range -Iinclude" CXXFLAGS="-fopenmp -fPIC -O3 -funroll-loops -fcx-limited-range -Iinclude -std=c++14 -DFFTW_PLAN_SAFE" LIBSFFT="-lfftw3 -lfftw3f -lm"
mv lib/libfinufft.so "${libdir}/libfinufft.${dlext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libfinufft", :libfinufft)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil