text:
# Constraints
untrackedconstraint(model::Z3Model, constr) = model.ptr[:add](constr)
trackedconstraint(model::Z3Model, constr, trackvar) = model.ptr[:assert_and_track](constr,trackvar)

function constraint(model::Z3Model, constr)
    model.numconstr += 1
    numconstr = model.numconstr
    trackvar = variable(model,:Bool,name="c$numconstr")
    push!(model.constr,trackvar)
    trackedconstraint(model,constr,trackvar)
    return trackvar
end

function objective(m::Z3Model,sense::Symbol,expr)
    objvar = variable(m, name="obj")
    m.objvar = objvar
    if sense == :Min
        constraint(m, objvar >= expr)
    elseif sense == :Max
        constraint(m, objvar <= expr)
    else
        error("Invalid sense, valid senses are :Min, :Max")
    end
end

implies(m::Z3Model, prop1, prop2) = constraint(m, z3[:Implies](prop1,prop2))
atmost(m::Z3Model,arr::Array,k::Int64) = constraint(m, z3[:AtMost](arr..., k))
atleast(m::Z3Model,arr::Array,k::Int64) = constraint(m, z3[:AtLeast](arr..., k))


text:
"""
Sets the data source
data_values(*sym1*=*vec1*, *sym2*=*vec2*, ...)

Adds data vectors (`vec1`, `vec1`,..) and binds each of them to a symbol (`sym1`, `sym2`, ...)
"""
function data_values(;values...)
  length(values)==0 && error("no values given")

  ks = [v[1] for v in values]
  vs = [v[2] for v in values]

  n    = length(values[1][2])
  any(v -> length(v) != n, vs) && error("all values must have the same length")

  vss = [ Dict(zip(ks, [v[i] for v in vs])) for i in  1:n ]
  VegaLiteVis(Dict(:data => Dict(:values => vss)))
end

ext: jl


text:
#=
   https://www.allendowney.com/blog/2018/10/21/the-game-of-ur-problem/
   """
   Here’s a probability puzzle to ruin your week.

   In the Royal Game of Ur, players advance tokens along a track with 14 spaces.
   To determine how many spaces to advance, a player rolls 4 dice with 4 sides. Two corners
   on each die are marked; the other two are not. The total number of marked corners —
   which is 0, 1, 2, 3, or 4 — is the number of spaces to advance.

   For example, if the total on your first roll is 2, you could advance a token to space 2.
   If you roll a 3 on the next roll, you could advance the same token to space 5.

   Suppose you have a token on space 13. How many rolls did it take to get there?
   """

   See:
   https://www.allendowney.com/blog/lions-and-tigers-and-bears/

   Allen Downey's solution:
   http://nbviewer.jupyter.org/github/AllenDowney/ThinkBayes2/blob/master/solutions/game_of_ur_soln.ipynb?flush=true

   cf ~/blog/game_of_ur_problem.blog
      ~/webppl/game_of_ur_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function game_of_ur_problem2()

    function roll1(a)
        ss = sum(a)
        if ss == 13
            return  a
        elseif ss > 13
            return []
        else
            t = rand(DiscreteUniform(0,4))
            return roll1(vcat(a,t))
        end
    end

    # Simpler
    function roll(s,len)
        if s == 13
            return  len
        elseif s > 13
            return 0
        else
            t = rand(DiscreteUniform(0,4))
            return roll(s+t,len+1)
        end
    end

    len ~ DiscreteUniform(0,100)
    # a = roll(Int64[])

    # len = length(a)
    # s = len > 0 ? sum(a) : 0

    # Dirac(len > 0 && s == 13)
    # obs!(_sampler,_varinfo, s == 13) && return
    # len > 0 || begin Turing.@addlogprob! -Inf; return end

    len = roll(0,0)
    return len

end

model = game_of_ur_problem2()

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, MH(), 10_000)
# chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, PG(20), 1_000)

chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 40_000)

# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), 10_000)

# chains = sample(model, NUTS(1000,0.65), 1_000)
# chains = sample(model, HMC(0.1,5), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(10,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(PG(10,:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1_000,0.65,:height)), 1_000)

display(chains)
# display(plot(chains))

println("Distribution of len")
gen = generated_quantities(model, chains)
println("After removing the 0s from gen")
gen2 = [g for g in gen if g != 0]
show_var_dist_pct(gen2, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import .Base: unsafe_convert, lock, trylock, unlock, islocked, wait, notify, AbstractLock

export SpinLock

# Important Note: these low-level primitives defined here
#   are typically not for general usage

##########################################
# Atomic Locks
##########################################

"""
    SpinLock()

Create a non-reentrant, test-and-test-and-set spin lock.
Recursive use will result in a deadlock.
This kind of lock should only be used around code that takes little time
to execute and does not block (e.g. perform I/O).
In general, [`ReentrantLock`](@ref) should be used instead.

Each [`lock`](@ref) must be matched with an [`unlock`](@ref).

Test-and-test-and-set spin locks are quickest up to about 30ish
contending threads. If you have more contention than that, different
synchronization approaches should be considered.
"""
mutable struct SpinLock <: AbstractLock
    handle::Int
    SpinLock() = new(0)
end

import Base.Sys.WORD_SIZE

@eval _xchg!(x::SpinLock, v::Int) =
    llvmcall($"""
             %ptr = inttoptr i$WORD_SIZE %0 to i$WORD_SIZE*
             %rv = atomicrmw xchg i$WORD_SIZE* %ptr, i$WORD_SIZE %1 acq_rel
             ret i$WORD_SIZE %rv
             """, Int, Tuple{Ptr{Int}, Int}, unsafe_convert(Ptr{Int}, pointer_from_objref(x)), v)

@eval _get(x::SpinLock) =
    llvmcall($"""
             %ptr = inttoptr i$WORD_SIZE %0 to i$WORD_SIZE*
             %rv = load atomic i$WORD_SIZE, i$WORD_SIZE* %ptr acquire, align $(gc_alignment(Int))
             ret i$WORD_SIZE %rv
             """, Int, Tuple{Ptr{Int}}, unsafe_convert(Ptr{Int}, pointer_from_objref(x)))

@eval _set!(x::SpinLock, v::Int) =
    llvmcall($"""
             %ptr = inttoptr i$WORD_SIZE %0 to i$WORD_SIZE*
             store atomic i$WORD_SIZE %1, i$WORD_SIZE* %ptr release, align $(gc_alignment(Int))
             ret void
             """, Cvoid, Tuple{Ptr{Int}, Int}, unsafe_convert(Ptr{Int}, pointer_from_objref(x)), v)

# Note: this cannot assert that the lock is held by the correct thread, because we do not
# track which thread locked it. Users beware.
Base.assert_havelock(l::SpinLock) = islocked(l) ? nothing : Base.concurrency_violation()

function lock(l::SpinLock)
    while true
        if _get(l) == 0
            p = _xchg!(l, 1)
            if p == 0
                return
            end
        end
        ccall(:jl_cpu_pause, Cvoid, ())
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Cvoid, ())
    end
end

function trylock(l::SpinLock)
    if _get(l) == 0
        return _xchg!(l, 1) == 0
    end
    return false
end

function unlock(l::SpinLock)
    _set!(l, 0)
    ccall(:jl_cpu_wake, Cvoid, ())
    return
end

function islocked(l::SpinLock)
    return _get(l) != 0
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using QuasiGeostrophy
using LinearAlgebra, Test, FFTW, BenchmarkTools, Plots
include(pwd() * "/sandbox/qg_utils.jl") # timestepping and wrappers here
# Domain
Ω = S¹(0, 2π) × S¹(0, 2π)
Nx = Ny = 2^5;
array = Array # to switch between gpu and cpu
fourier_grid = FourierGrid((Nx, Ny), Ω, arraytype = array)
# Fields
fieldnames = ("q¹", "q²", "ψ¹", "ψ²", "f")
q¹, q², ψ¹, ψ², f = create_fields(names = fieldnames, grid = fourier_grid, arraytype = array)
# Operators
∂x, ∂y = create_operators(fourier_grid, arraytype = array)
Δ = ∂x*∂x + ∂y*∂y  # multiplication for GPU purposes
Δ⁴ = -Δ  #* -Δ * -Δ * -Δ  # multiplication for GPU purposes
# Parameters (λ/L)² = 0.01 => λ = 0.1 L
const λ  = 0.01*2π # Rossby deformation radius, √(gravity * depth) / (coriolis force)
const U  = 1e-4    # background velocity in the top layer, - in bottom
const κ  = 1.0     # bottom drag
const Q  = 0.5     # surface temperature flux

Δtᵐᵃˣ = 1e-1 / Nx     # Largest timestep size
ν = 1e-0/maximum(abs.((Δ⁴.op)))/Δtᵐᵃˣ  # largest stable damping
# initialize stream functions in layers
x, y = fourier_grid.grid
ψ¹( 1.0 * (sin.(x) .* sin.(y) .- 1.0 * sin.(y)) )
ψ²( 1.0 *( (cos.(x) .* sin.(y)) .+ 0.1 * sin.(8*y) .+ 1.0 * sin.(y)  ))
# initialize forcing
f(Q * sin.(y))
# stream function to potential vorticity
q¹.data .= (Δ(ψ¹) + 1 /(2*λ^2) * (ψ² - ψ¹)).data
q².data .= (Δ(ψ²) + 1 /(2*λ^2) * (ψ¹ - ψ²)).data
# potential vorticity to stream function
G = inv( (Δ + -1 /(2*λ^2) ) * (Δ + -1 /(2*λ^2) ) + -(1/(2*λ^2))^2 )
ψ¹.data .= G(Δ(q¹) + 1 /(2*λ^2) * (-q² - q¹)).data
ψ².data .= G(Δ(q²) + 1 /(2*λ^2) * (-q¹ - q²)).data

# Impero
@wrapper q¹=q¹ ψ¹=ψ¹ q²=q² ψ²=ψ² f=f
Δ = Operator(Δ); ∂x = Operator(∂x); ∂y = Operator(∂y)
Δ⁴ = Operator(Δ⁴);
G = Operator(nothing, OperatorMetaData(G, "G"))
∂t = Operator(nothing, OperatorMetaData(nothing, "∂t"))
drag = (-2 * κ) * Δ(ψ²)

# create cutoff filter
fmd = FourierOperatorMetaData("Filter")
dx, dy = create_operators(fourier_grid, arraytype = array)
wavenums = sqrt.(abs.((dx^2 + dy^2).op))
maxk = maximum(wavenums)
op = cutoff.(wavenums, 1/2 * maxk) #max k is about \sqrt(2) larger than max kx
cutoff_filter = FourierOperator(op, fmd)
ℱ = Operator(nothing, OperatorMetaData(cutoff_filter, "ℱ"))

# more dissipative QG
@pde_system qg_system = (
    ψ¹ = G(Δ(q¹) + 1/(2*λ^2) * (-q²-q¹)),
    ψ² = G(Δ(q²) + 1/(2*λ^2) * (-q¹-q²)),
    ∂t(q¹) = -ν*Δ⁴(q¹) + (ℱ(∂y(ψ¹))*ℱ(∂x(q¹)) - ℱ(∂x(ψ¹))*ℱ(∂y(q¹))) - (U/λ^2) * ∂x(ψ¹) - U * ∂x(q¹) + f,
    ∂t(q²) = -ν*Δ⁴(q²) + (ℱ(∂y(ψ²))*ℱ(∂x(q²)) - ℱ(∂x(ψ²))*ℱ(∂y(q²))) + (U/λ^2) * ∂x(ψ²) + U * ∂x(q²) - f + drag
)

## Timestepping in qg_utils
vort1 = []
vort2 = []
for i in 1:10000
    # adaptive timestepping
    Δt = 0.3 * cfl(ψ¹, ψ², ∂x, ∂y)
    Δt = minimum([Δt, Δtᵐᵃˣ])
    evolve_system(qg_system, Δt)

    if i%100==0
        println("hi at "* string(i))

        p1 = plot(q¹.data)
        p2 = plot(q².data)
        display(plot(p1,p2))
        sleep(0.01)
        println("hi at "* string(i))
        println("The norm is ")
        println(norm(compute(∂y(ψ¹)).data)/length(ψ¹.data.data))
        println("The norm is ")
        println(norm(compute(∂y(ψ²)).data)/length(ψ².data.data))
        println("The Δt is ")
        println(Δt)
        println("The CFL is ")
        println(Δt/cfl(ψ¹, ψ², ∂x, ∂y) )

        P = q¹.data.metadata.transform.backward
        push!( vort1 , real.(P * q¹.data.data))
        push!( vort2 , real.(P * q².data.data))
    end
end
@save "qg_stuff.jld2" vort1 vort2

plot(q¹.data)
plot(q².data)

spectrum(q¹.data-q².data)
pyplot()
tmp = q¹.data.metadata.transform.backward * (q¹.data.data - q².data.data)
contourf(real.(tmp'), color = :thermometer, linewidth = 0, contourlevlels = 30)
y = sum(real.(tmp'), dims = 2) ./Nx
plot(x, y, label = "horizontal average")
plot!(x,sin.(x) * maximum(y), label = "sin")

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
#= Extra data for long and short symbols =#

# short symbol: may be ambiguous, take the most commonly-found unicode symbol.
# NB: couldn't find a good free source of this data, so it's very incomplete.
# send a pull request to contribute!
const SHORT_SYMBOL = Dict{Symbol, String}(
    :AED => "د.إ",
    :AFN => "؋",
    :ALL => "L",
    :AMD => "֏",
    :ANG => "ƒ",
    :AOA => "Kz",
    :ARS => "\$",
    :AUD => "\$",
    :AWG => "Afl.",
    :AZN => "₼",
    :BBD => "\$",
    :BSD => "\$",
    :CAD => "\$",
    :CHF => "Fr.",
    :CNY => "¥",
    :EUR => "€",
    :GBP => "£",
    :JPY => "¥",
    :KPW => "₩",
    :KRW => "₩",
    :MGA => "Ar",
    :MXN => "\$",
    :NOK => "kr",
    :NZD => "\$",
    :USD => "\$")

# long symbol: short where possible without being ambiguous
const LONG_SYMBOL = Dict{Symbol, String}(
    :AED => "د.إ",
    :AFN => "؋",
    :ALL => "L",
    :AMD => "֏",
    :ANG => "ƒ",
    :AOA => "Kz",
    :ARS => "AR\$",
    :AUD => "A\$",
    :AWG => "Afl.",
    :AZN => "₼",
    :BBD => "Bds\$",
    :BSD => "B\$",
    :CAD => "C\$",
    :CHF => "Fr.",
    :CNY => "CN¥",
    :EUR => "€",
    :GBP => "GB£",
    :JPY => "JP¥",
    :KPW => "KP₩",
    :KRW => "KR₩",
    :MGA => "Ar",
    :MXN => "MX\$",
    :NZD => "NZ\$",
    :USD => "US\$")


text:
using DataFrames
using MixedModels
using MixedModelsPermutations
using MixedModelsPermutations: permute!
using StableRNGs
using Statistics
using Test

isdefined(@__MODULE__, :io) || const io = IOBuffer()

@testset "LMM" begin
    sleepstudy = MixedModels.dataset(:sleepstudy)
    m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)
    rm1 = permute!(StableRNG(42), deepcopy(m1); residual_method=:signflip)
    rm1 = permute!(StableRNG(42), deepcopy(m1); residual_method=:shuffle)
    # test marking as not fit
    @test_logs (:warn,) show(io, rm1)
    @test_throws ArgumentError permute!(deepcopy(m1); residual_method=:bad)

    H0 = coef(m1)
    H0[2] = 0.0 # slope of days is 0

    perm = permutation(StableRNG(42),1000, m1; β=H0)
    @test perm isa MixedModelPermutation

    df = combine(groupby(DataFrame(perm.allpars), [:type, :group, :names]),
                :value => shortestcovint => :interval)

    # our original slope should lie well outside of the null distribution
    let days = filter(df) do row
            return row.type == "β" && row.names == "days"
        end

        lower, upper = only(days.interval)

        @test coef(m1)[2] > upper
        @test coef(m1)[2] > -lower
    end

    # our original intercept should lie within the range of
    # the intercepts generated via permutation
    let est = filter(df) do row
            return row.type == "β" && row.names == "(Intercept)"
        end

        lower, upper = only(est.interval)
        @test lower < first(H0) < upper
    end

    @testset "permutationtest" begin
        # the total area above and below the observed value should equal 1
        @test all(map(sum,
                  zip(permutationtest(perm, m1, :lesser),
                      permutationtest(perm, m1, :greater))) .== 1)

        # we should have a p-value near 1 since our null distribution
        # was generated from value being tested...
        @test first(permutationtest(perm, m1)) > 0.95
        # we should have a p-value near 0 since this effect is clear
        @test last(permutationtest(perm, m1)) <= 1 / 1000
        @test_throws ArgumentError permutationtest(perm, m1, :bad)
    end

    @testset "olsranef" begin
        permols = permutation(StableRNG(42),1000, m1; β=H0, blup_method=olsranef)
        @test permols isa MixedModelPermutation
        @test last(permutationtest(perm, m1, :greater)) ≈ last(permutationtest(permols, m1, :greater)) atol=0.01
        @test last(permutationtest(perm, m1, :lesser)) ≈ last(permutationtest(permols, m1, :lesser)) atol=0.01
    end

end

@testset "GLMM" begin
    cbpp = MixedModels.dataset(:cbpp)
    gm1 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)),
              cbpp, Binomial(); wts=cbpp.hsz, fast=true)

    @test_throws MethodError permutation(1, gm1)
    @test_throws MethodError permute!(gm1)
end


text:
# ps3000_run_block.jl

export ps3000_run_block

function ps3000_run_block(handle::Int16, no_of_samples, timebase, oversample)
  time_indisposed_ms = Ref{Int32}()
  ps_status = ccall((:ps3000_run_block, ps3000driver), Int16,
    (Int16, Int32, Int16, Int16, Ref{Int32}),
    handle, no_of_samples, timebase, oversample, time_indisposed_ms)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return time_indisposed_ms[]
end


text:
using DelimitedFiles
function GenCurveLoad(filename_GenCurve)
    ## source from: https://www.iso-ne.com/isoexpress/web/reports/operations/-/tree/daily-gen-fuel-type
    temp_data = readdlm(filename_GenCurve, ',',header=true)[1][:,4]
    raw_data = zeros(365,24)
    for dd in 1:365
       for tt in 1:24
          if ~isempty(temp_data[24*(dd-1)+tt])
             raw_data[dd,tt] = temp_data[24*(dd-1)+tt]
          end
       end
    end
    # normalized_data = raw_data/maximum(raw_data)
    # GenCurve = normalized_data
    # return GenCurve
    return raw_data
end

max_stars_repo_name: JipKim/IA

max_issues_repo_name: JipKim/IA

max_forks_repo_name: JipKim/IA


text:
# This is a single line comment

x = 9.4 	# This is an inline comment
println(x)

#=
This is a multiline comment.
Second line of the comments is here.
=#

println("Done")

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
module pwr

using Distributions, Roots, DataStructures

export cohenES, ESh, ESw1, ESw2
export powerAnovaTest, samplesizeAnovaTest, effectsizeAnovaTest, alphaAnovaTest
export powerChisqTest, samplesizeChisqTest, effectsizeChisqTest, alphaChisqTest
export powerF2Test, samplesizeF2Test, effectsizeF2Test, alphaF2Test
export powerNormTest, samplesizeNormTest, effectsizeNormTest, alphaNormTest
export powerPTest, samplesizePTest, effectsizePTest, alphaPTest
export powerTwoPTest, samplesizeTwoPTest, effectsizeTwoPTest, alphaTwoPTest
export powerTwoP2nTest, samplesizeTwoP2nTest, effectsizeTwoP2nTest, alphaTwoP2nTest
export powerRTest, samplesizeRTest, effectsizeRTest, alphaRTest
export powerT2nTest, samplesizeT2nTest, effectsizeT2nTest, alphaT2nTest
export powerTTest, samplesizeTTest, effectsizeTTest, alphaTTest
export plot

include("utils.jl")
include("cohenES.jl")
include("ESh.jl")
include("ESw1.jl")
include("ESw2.jl")
include("AnovaTest.jl")
include("ChisqTest.jl")
include("F2Test.jl")
include("NormTest.jl")
include("PTest.jl")
include("RTest.jl")
include("T2nTest.jl")
include("TTest.jl")
include("Twop2nTest.jl")
include("TwopTest.jl")
include("PowerPlot.jl")

end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
# TODO -- convert to use add_vertex! instead, since edges type must be made also
function addVerticesSubgraph(fgl::FactorGraph,
    fgseg::FactorGraph,
    vertdict::Dict{Int,Graphs.ExVertex})

    for vert in vertdict
      fgseg.g.vertices[vert[1]] = vert[2]
      if haskey(fgl.v,vert[1])
        fgseg.g.vertices[vert[1]] = vert[2]
        fgseg.IDs[Symbol(vert[2].label)] = vert[1]

        # add edges going in opposite direction
        elr = Graphs.out_edges(vert[2], fgl.g)
        len = length(elr)
        keeprm = trues(len)
        j = 0
        for i in 1:len
          if !haskey(vertdict, elr[i].target.index) # a function node in set, so keep ref
            keeprm[i] = false
            j+=1
          end
        end
        if j < len
          elridx = elr[1].source.index
          fgseg.g.inclist[elridx] = elr[keeprm]
        end
      elseif haskey(fgl.f, vert[1])
        fgseg.f[vert[1]] = vert[2] # adding element to subgraph
        fgseg.fIDs[Symbol(vert[2].label)] = vert[1]
        # get edges associated with function nodes and push edges onto incidence list
        el = Graphs.out_edges(vert[2], fgl.g)
        elidx = el[1].source.index
        fgseg.g.inclist[elidx] = el # okay because treating function nodes only
        fgseg.g.nedges += length(el)
      else
        error("Unknown type factor graph vertex type, something is wrong")
      end
    end
    nothing
end

"""
    $SIGNATURES

Construct a new factor graph object as a subgraph of `fgl::FactorGraph` based on the
variable labels `syms::Vector{Symbols}`.

Notes
- Slighly messy internals, but gets the job done -- some room for performance improvement.

Related

getVariableIds
"""
function buildSubgraphFromLabels(dfg::G,
                                 syms::Vector{Symbol},
                                 destType::Type{<:AbstractDFG}=InMemDFGType ) where G <: AbstractDFG

  # data structure for cliq sub graph
  if G <: InMemoryDFGTypes
    #Same type
    cliqSubFg = initfg(G, params=getSolverParams(dfg))
  else
    #Default
    cliqSubFg = initfg(destType, params=getSolverParams(dfg))
  end

  # add a little too many variables (since we need the factors)
  for sym in syms
    DFG.getSubgraphAroundNode(dfg, DFG.getVariable(dfg, sym), 2, false, cliqSubFg)
  end

  # remove excessive variables that were copied by neighbors distance 2
  currVars = DFG.getVariableIds(cliqSubFg)
  toDelVars = setdiff(currVars, syms)
  for dv in toDelVars
    # delete any neighboring factors first
    for fc in DFG.lsf(cliqSubFg, dv)
      DFG.deleteFactor!(cliqSubFg, fc)
    end

    # and the variable itself
    DFG.deleteVariable!(cliqSubFg, dv)
  end

  return cliqSubFg
end
# function buildSubgraphFromLabels(dfg::G, syms::Vector{Symbol}) where G <: AbstractDFG
#   fgseg = initfg() #sessionname=dfg.sessionname, robotname=dfg.robotname)
#
#   for sym in syms
#     vert = DFG.getVariable(dfg, sym) #, api=localapi)
#     st = getSofttype(vert)
#     addVariable!(fgseg, sym, st) #, uid=vert.index)
#     if isInitialized(dfg,sym)
#       manualinit!(fgseg, sym, getKDE(vert))
#     end
#   end
#
#   for sym in syms
#     for fct in DFG.ls(dfg, :x1)
#       if !hasFactor(fgseg, fct)
#         # check all variables are in desired variable set
#         possibleVars = DFG.lsf(dfg, fct)
#         ivars = intersect(possibleVars, syms)
#         @show length(ivars), length(possibleVars)
#         if length(ivars) == length(possibleVars)
#           # fvert = getVert(dfg, fct, api=localapi, nt=:fct)
#           ufc = DFG.getFactor(dfg, fct) # fvert
#
#           addFactor!(fgseg, possibleVars, getData(ufc).fnc.usrfnc!, autoinit=false) #, uid=fvert.index)
#         end
#       end
#     end
#   end
#
#   return fgseg
# end

# NOTICE, nodeIDs and factorIDs are not brough over by this method yet
# must sort out for incremental updates
function genSubgraph(fgl::FactorGraph,
                     vertdict::Dict{Int,Graphs.ExVertex})
  # edgedict::Dict{Int,Graphs.Edge{Graphs.ExVertex}},
  #

  fgseg = initfg(sessionname=fgl.sessionname, robotname=fgl.robotname)

  error("genSubgraph() is obsolete")

  # fgseg = FactorGraph() # new handle for just a segment of the graph
  # fgseg.g = Graphs.inclist(Graphs.ExVertex, is_directed=false)
  #
  # fgseg.v = Dict{Int,Graphs.ExVertex}()
  # fgseg.f = Dict{Int,Graphs.ExVertex}()
  # fgseg.IDs = Dict{AbstractString,Int}()
  # fgseg.fIDs = Dict{AbstractString,Int}()
  #
  # # TODO -- convert to use empty constructor since Graphs.incdict now works
  # fgseg.g.vertices = Array{Graphs.ExVertex,1}(length(fgl.g.vertices))
  # fgseg.g.inclist = Array{Array{Graphs.Edge{Graphs.ExVertex},1},1}(length(fgl.g.inclist))
  #
  # addVerticesSubgraph(fgl, fgseg, vertdict)
  #
  # fgseg.id = fgl.id
  # fgseg.bnid = fgl.bnid
  # fgseg.dimID = fgl.dimID

  return fgseg
end

function getShortestPathNeighbors(fgl::FactorGraph;
    from::Graphs.ExVertex=nothing,
    to::Graphs.ExVertex=nothing,
    neighbors::Int=0 )

  edgelist = shortest_path(fgl.g, ones(num_edges(fgl.g)), from, to)
  vertdict = Dict{Int,Graphs.ExVertex}()
  edgedict = edgelist2edgedict(edgelist)
  expandVertexList!(fgl, edgedict, vertdict) # grow verts
  for i in 1:neighbors
    expandEdgeListNeigh!(fgl, vertdict, edgedict) # grow edges
    expandVertexList!(fgl, edgedict, vertdict) # grow verts
  end
  return vertdict
end

function subgraphShortestPath(fgl::FactorGraph;
                              from::Graphs.ExVertex=nothing,
                              to::Graphs.ExVertex=nothing,
                              neighbors::Int=0  )
  #
  vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
  return genSubgraph(fgl, vertdict)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
function subGraphFromVerts(fgl::FactorGraph,
                           verts::Dict{Int,Graphs.ExVertex};
                           neighbors::Int=0  )
  #
  allverts = Dict{Int,Graphs.ExVertex}()
  allkeys = collect(keys(verts))
  len = length(allkeys)
  # union all shortest path combinations in a vertdict
  for i in 1:len, j in (i+1):len
    from = verts[allkeys[i]]
    to = verts[allkeys[j]]
    vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
    for vert in vertdict
      if !haskey(allverts, vert[1])
        allverts[vert[1]] = vert[2]
      end
    end
  end

  return genSubgraph(fgl, allverts)
end

"""
    $SIGNATURES

Transfer contents of `src::FactorGraph` variables `syms::Vector{Symbol}` to `dest::FactorGraph`.

Notes
- Reads, `dest` := `src`, for all `syms`
"""
function transferUpdateSubGraph!(dest::G1,
                                 src::G2,
                                 syms::Vector{Symbol}=union(ls(src)...),
                                 logger=ConsoleLogger()  ) where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  with_logger(logger) do
    @info "transferUpdateSubGraph! -- syms=$syms"

    # TODO add with DFG v0.4
    # DFG.updateGraphSolverData!(src, dest, syms)
    for sym in syms
      vari = DFG.getVariable(src, sym)
      rc = size(solverData(vari).val)
      pp = getKDE(vari)
      rc2 = size(getPoints(pp))
      @info "sym=$sym, mem size of val=$rc and $(rc2)"
      updateFullVertData!(dest, vari, updateMAPest=true)
    end
  end
  nothing
end


text:
# using LinearAlgebra, Plots, JSON
include("src/SFFM.jl")
include("dev/METools.jl")

# define SFM
T =[-2.0 2.0 0; 1.0 -2.0 1; 4 3 -7]# [-2.0 2.0; 1.0 -1.0]#
C = [1.0; 2.0; 0]
fn(x) = [ones(size(x)) ones(size(x)) ones(size(x))]
model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])
N₋ = sum(C.<=0)
N₊ = sum(C.>=0)
NPhases = length(C)

# construct global approximation

t = 3
τ = SFFM.FixedTime(T=t)
NSim = 400_000
sims = SFFM.SimSFM(model=model,StoppingTime=τ,InitCondition=(φ=2*ones(Int,NSim),X=zeros(NSim)))

p = plot()
vecΔ =  2#[5 2.5 1.25 1.25/2]
let
    vecNBases = [1,3,5,7,11,15,21]#,29]
    errPH = vecNBases
    errME = vecNBases
    errbkwdME = vecNBases
    errπME = vecNBases
    errDG = vecNBases
    c=0
    for Δ in vecΔ
        c+=1
        Nodes = collect(0:Δ:10)
        globalerrPH = []
        globalerrME = []
        globalerrbkwdME = []
        globalerrπME = []
        globalerrDG = []
        for NBases in vecNBases
            mesh = SFFM.MakeMesh(model=model,NBases=NBases,Nodes=Nodes,Basis="lagrange")
            simDist = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")

            # define generator for up approximation
            Erlang = MakeErlang(NBases, mean = Δ)

            NCells = length(Nodes)-1
            Q, B = MakeGlobalApprox(;
                NCells = NCells,
                up = Erlang,
                down = Erlang,
                T = T,
                C = C,
            )

            ME = MakeME(CMEParams[NBases],mean=Δ)

            QME, BME = MakeGlobalApprox(;
                NCells = NCells,
                up = ME,
                down = ME,
                T = T,
                C = C,
            )

            MEBkwd = reversal(ME)
            # αdownMEbkwd = αupME#[αupME[1]; [αupME[3:2:end] αupME[2:2:end]]'[:]]'
            # QdownMEbkwd = QupME#.*sum(-αdownMEbkwd*QupME'^-1)./Δ

            QMEbkwd, BMEbkwd = MakeGlobalApprox(;
                NCells = NCells,
                up = ME,
                down = ME,
                T = T,
                C = C,
                bkwd = false,
                D = [],
                plusI = false,
            )

            QMEπ, BMEπ = MakeGlobalApprox(;
                NCells = NCells,
                up = ME,
                down = ME,
                T = T,
                C = C,
                bkwd = true,
                D = [],
                plusI = false,
            )

            DGMesh = SFFM.MakeMesh(model=model,NBases=1,Nodes=collect(Nodes[1]:Δ/NBases:Nodes[end]),Basis="lagrange")
            All = SFFM.MakeAll(model=model,mesh=DGMesh)

            initDist = zeros(1,size(B,1))
            initDist[1] = 1

            temp = initDist*exp(Matrix(All.B.B)*t)#SFFM.EulerDG(D=All.B.B,y=t,x0=initDist)#
            DGdist_t = SFFM.Coeffs2Dist(model=model,mesh=DGMesh,Coeffs=temp,type="probability")

            dist_t = initDist*exp(B*t)#SFFM.EulerDG(D=B,y=t,x0=initDist)#
            pm_t = dist_t[[1:N₋;(end-N₊+1):end]]
            dist_t = reshape(dist_t[N₋+1:end-N₊],NBases,NPhases,NCells)

            distME_t = initDist*exp(BME*t)#SFFM.EulerDG(D=BME,y=t,x0=initDist)
            pmME_t = distME_t[[1:N₋;(end-N₊+1):end]]
            distME_t = reshape(distME_t[N₋+1:end-N₊],NBases,NPhases,NCells)

            distMEbkwd_t = initDist*exp(BMEbkwd*t)#SFFM.EulerDG(D=BME,y=t,x0=initDist)
            pmMEbkwd_t = distMEbkwd_t[[1:N₋;(end-N₊+1):end]]
            distMEbkwd_t = reshape(distMEbkwd_t[N₋+1:end-N₊],NBases,NPhases,NCells)

            distMEπ_t = initDist*exp(BMEπ*t)#SFFM.EulerDG(D=BME,y=t,x0=initDist)
            pmMEπ_t = distMEπ_t[[1:N₋;(end-N₊+1):end]]
            distMEπ_t = reshape(distMEπ_t[N₋+1:end-N₊],NBases,NPhases,NCells)

            # plot(layout = (NPhases,1))
            localerrPH = sum(abs.(pm_t-simDist.pm))
            localerrME = sum(abs.(pmME_t-simDist.pm))
            localerrbkwdME = sum(abs.(pmMEbkwd_t-simDist.pm))
            localerrπME = sum(abs.(pmMEπ_t-simDist.pm))
            localerrDG = sum(abs.(DGdist_t.pm-simDist.pm))
            for n in 1:NCells
                x = simDist.x[n]
                for i in 1:NPhases
                    if n==1
                        label1 = "PH"
                        label2 = "DG"
                        label3 = "ME"
                        label4 = "bkwdME"
                        label5 = "πME"
                    else
                        label1 = false
                        label2 = false
                        label3 = false
                        label4 = false
                        label5 = false
                    end
                    yvalsPH = [sum(dist_t[:,i,n])]
                    yvalsDG = [sum(DGdist_t.distribution[:,(1:NBases).+NBases*(n-1),i])]
                    yvalsME = [sum(distME_t[:,i,n])]
                    yvalsbkwdME = [sum(distMEbkwd_t[:,i,n])]
                    yvalsπME = [sum(distMEπ_t[:,i,n])]
                    # scatter!([x],yvalsPH,label=label1,subplot=i,color=:red,markershape=:x)
                    # scatter!([x],yvalsDG,label=label2,subplot=i,color=:blue,markershape=:rtriangle)
                    # scatter!([x],yvalsME,label=label3,subplot=i,color=:black,markershape=:ltriangle)
                    # scatter!([x],yvalsbkwdME,label=label4,subplot=i,color=:black,markershape=:utriangle)
                    localerrPH += abs(sum(yvalsPH-simDist.distribution[:,n,i]))
                    localerrDG += abs(sum(yvalsDG-simDist.distribution[:,n,i]))
                    localerrME += abs(sum(yvalsME-simDist.distribution[:,n,i]))
                    localerrbkwdME += abs(sum(yvalsbkwdME-simDist.distribution[:,n,i]))
                    localerrπME += abs(sum(yvalsπME-simDist.distribution[:,n,i]))
                end
            end
            push!(globalerrPH,localerrPH)
            push!(globalerrDG,localerrDG)
            push!(globalerrME,localerrME)
            push!(globalerrbkwdME,localerrbkwdME)
            push!(globalerrπME,localerrπME)
            # for i in 1:NPhases
            #     scatter!(simDist.x,simDist.distribution[:,:,i][:],subplot=i,label="sim")
            # end
            # display(plot!())
        end
        errPH = [errPH globalerrPH]
        errDG = [errDG globalerrDG]
        errME = [errME globalerrME]
        errbkwdME = [errbkwdME globalerrbkwdME]
        errπME = [errπME globalerrπME]
        # plot!(p,vecNBases,log.(globalerrPH),label=false,linestyle=:dot,colour=c,linewidth=2)
        # plot!(p,vecNBases,log.(globalerrME),label=false,linestyle=:dash,colour=c)
        # plot!(p,vecNBases,log.(globalerrbkwdME),label=false,linestyle=:dashdot,colour=c)
        # plot!(p,vecNBases,log.(globalerrπME),label=false,linestyle=:dashdotdot,colour=c)
        plot!(p,vecNBases,log.(globalerrDG),label="DG Δ",color=c)
        plot!(p,legend=:topright,xlabel="Order",ylabel="log(error)")
        display(plot!(p,title="DG"))
        display(globalerrbkwdME[end])
    end
    # plot(log.(vecΔ)[:], log.(errPH[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dot,label=false,linewidth=2)
    # plot!(log.(vecΔ)[:], log.(errME[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dash,label=false)
    # plot!(log.(vecΔ)[:], log.(errbkwdME[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dashdot,label=false)
    # plot!(log.(vecΔ)[:], log.(errDG[:,2:end]'),colour=[1 2 3 4 5],label=vecNBases')
    # plot!(legend=:bottomright,xlabel="log(Δ)",ylabel="log(error)",legendtitle="Order")
    # plot!(title="... PH,  -- ME,  -. pois ME, Solid DG")
end
Δ = vecΔ

#
# let
#     vecNBases = [1,3,5,7,11,15,21]#,29]
#     vecΔ = [1.25/2]# 1.25/2]#[5 2.5 1.25 1.25/2 1.25/4]
#     errPH = vecNBases
#     errCoxian = vecNBases
#     errbkwdCoxian = vecNBases
#     errSomePH = vecNBases
#     errSomePHBkwd = vecNBases
#     errDG = vecNBases
#     plot()
#     c=0
#     for Δ in vecΔ
#         c+=1
#         Nodes = collect(0:Δ:10)
#         globalerrPH = []
#         globalerrCoxian = []
#         globalerrbkwdCoxian = []
#         globalerrSomePH = []
#         globalerrSomePHBkwd = []
#         globalerrDG = []
#         for NBases in vecNBases
#             mesh = SFFM.MakeMesh(model=model,NBases=NBases,Nodes=Nodes,Basis="lagrange")
#             simDist = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")
#
#             # define generator for up approximation
#             Erlang = MakeErlang(NBases, mean = Δ)
#
#             NCells = length(Nodes)-1
#             Q, B = MakeGlobalApprox(;
#                 NCells = NCells,
#                 up = Erlang,
#                 down = Erlang,
#                 T = T,
#                 C = C,
#             )
#
#             Coxian = MakeME(CMEParams[NBases], mean = Δ)#MakeSomeCoxian(NBases, mean = Δ)#
#
#             QCoxian, BCoxian = MakeGlobalApprox(;
#                 NCells = NCells,
#                 up = Coxian,
#                 down = Coxian,
#                 T = T,
#                 C = C,
#             )
#
#             CoxianBkwd = reversal(Coxian)
#
#             QCoxianbkwd, BCoxianbkwd = MakeGlobalApprox(;
#                 NCells = NCells,
#                 up = Coxian,
#                 down = CoxianBkwd,
#                 T = T,
#                 C = C,
#             )
#
#             SomePH = Coxian#MakeSomeCoxian(NBases, mean = Δ)#MakeSomePH(NBases, mean = Δ)
#
#             QSomePH, BSomePH = MakeGlobalApprox(;
#                 NCells = NCells,
#                 up = SomePH,
#                 down = SomePH,
#                 T = T,
#                 C = C,
#             )
#
#             # SomePHBkwd = reversal(SomePH)
#
#             QSomePHBkwd, BSomePHBkwd = MakeGlobalApprox(;
#                 NCells = NCells,
#                 up = SomePH,
#                 down = SomePH,
#                 T = T,
#                 C = C,
#                 bkwd = true,
#                 D = jumpMatrixD(SomePH)
#             )
#
#             DGMesh = SFFM.MakeMesh(model=model,NBases=1,Nodes=collect(Nodes[1]:Δ/NBases:Nodes[end]),Basis="lagrange")
#             All = SFFM.MakeAll(model=model,mesh=DGMesh)
#
#             initDist = zeros(1,size(B,1))
#             initDist[1] = 1
#
#             temp = initDist*exp(Matrix(All.B.B)*t)#SFFM.EulerDG(D=All.B.B,y=t,x0=initDist)#
#             DGdist_t = SFFM.Coeffs2Dist(model=model,mesh=DGMesh,Coeffs=temp,type="probability")
#
#             dist_t = initDist*exp(B*t)#SFFM.EulerDG(D=B,y=t,x0=initDist)#
#             pm_t = dist_t[[1:N₋;(end-N₊+1):end]]
#             dist_t = reshape(dist_t[N₋+1:end-N₊],NBases,NPhases,NCells)
#
#             distCoxian_t = initDist*exp(BCoxian*t)#SFFM.EulerDG(D=BCoxian,y=t,x0=initDist)
#             pmCoxian_t = distCoxian_t[[1:N₋;(end-N₊+1):end]]
#             distCoxian_t = reshape(distCoxian_t[N₋+1:end-N₊],NBases,NPhases,NCells)
#
#             distCoxianbkwd_t = initDist*exp(BCoxianbkwd*t)#SFFM.EulerDG(D=BCoxian,y=t,x0=initDist)
#             pmCoxianbkwd_t = distCoxianbkwd_t[[1:N₋;(end-N₊+1):end]]
#             distCoxianbkwd_t = reshape(distCoxianbkwd_t[N₋+1:end-N₊],NBases,NPhases,NCells)
#
#             distSomePH_t = initDist*exp(BSomePH*t)#SFFM.EulerDG(D=BCoxian,y=t,x0=initDist)
#             pmSomePH_t = distSomePH_t[[1:N₋;(end-N₊+1):end]]
#             distSomePH_t = reshape(distSomePH_t[N₋+1:end-N₊],NBases,NPhases,NCells)
#
#             distSomePHBkwd_t = initDist*exp(BSomePHBkwd*t)#SFFM.EulerDG(D=BCoxian,y=t,x0=initDist)
#             pmSomePHBkwd_t = distSomePHBkwd_t[[1:N₋;(end-N₊+1):end]]
#             distSomePHBkwd_t = reshape(distSomePHBkwd_t[N₋+1:end-N₊],NBases,NPhases,NCells)
#
#             # plot(layout = (NPhases,1))
#             localerrPH = sum(abs.(pm_t-simDist.pm))
#             localerrCoxian = sum(abs.(pmCoxian_t-simDist.pm))
#             localerrbkwdCoxian = sum(abs.(pmCoxianbkwd_t-simDist.pm))
#             localerrSomePH = sum(abs.(pmSomePH_t-simDist.pm))
#             localerrSomePHBkwd = sum(abs.(pmSomePHBkwd_t-simDist.pm))
#             localerrDG = sum(abs.(DGdist_t.pm-simDist.pm))
#             for n in 1:NCells
#                 x = simDist.x[n]
#                 for i in 1:NPhases
#                     if n==1
#                         label1 = "PH"
#                         label2 = "DG"
#                         label3 = "Coxian"
#                         label4 = "bkwdCoxian"
#                         label5 = "SomePH"
#                         label6 = "SomePHBkwd"
#                     else
#                         label1 = false
#                         label2 = false
#                         label3 = false
#                         label4 = false
#                         label5 = false
#                         label6 = false
#                     end
#                     yvalsPH = [sum(dist_t[:,i,n])]
#                     yvalsDG = [sum(DGdist_t.distribution[:,(1:NBases).+NBases*(n-1),i])]
#                     yvalsCoxian = [sum(distCoxian_t[:,i,n])]
#                     yvalsbkwdCoxian = [sum(distCoxianbkwd_t[:,i,n])]
#                     yvalsSomePH = [sum(distSomePH_t[:,i,n])]
#                     yvalsSomePHBkwd = [sum(distSomePHBkwd_t[:,i,n])]
#                     # scatter!([x],yvalsPH,label=label1,subplot=i,color=:red,markershape=:x)
#                     # scatter!([x],yvalsDG,label=label2,subplot=i,color=:blue,markershape=:rtriangle)
#                     # scatter!([x],yvalsCoxian,label=label3,subplot=i,color=:black,markershape=:ltriangle)
#                     # scatter!([x],yvalsbkwdCoxian,label=label4,subplot=i,color=:black,markershape=:utriangle)
#                     localerrPH += abs(sum(yvalsPH-simDist.distribution[:,n,i]))
#                     localerrDG += abs(sum(yvalsDG-simDist.distribution[:,n,i]))
#                     localerrCoxian += abs(sum(yvalsCoxian-simDist.distribution[:,n,i]))
#                     localerrbkwdCoxian += abs(sum(yvalsbkwdCoxian-simDist.distribution[:,n,i]))
#                     localerrSomePH += abs(sum(yvalsSomePH-simDist.distribution[:,n,i]))
#                     localerrSomePHBkwd += abs(sum(yvalsSomePHBkwd -simDist.distribution[:,n,i]))
#                 end
#             end
#             push!(globalerrPH,localerrPH)
#             push!(globalerrDG,localerrDG)
#             push!(globalerrCoxian,localerrCoxian)
#             push!(globalerrbkwdCoxian,localerrbkwdCoxian)
#             push!(globalerrSomePH,localerrSomePH)
#             push!(globalerrSomePHBkwd,localerrSomePHBkwd)
#             # for i in 1:NPhases
#             #     scatter!(simDist.x,simDist.distribution[:,:,i][:],subplot=i,label="sim")
#             # end
#             # display(plot!())
#         end
#         errPH = [errPH globalerrPH]
#         errDG = [errDG globalerrDG]
#         errCoxian = [errCoxian globalerrCoxian]
#         errbkwdCoxian = [errbkwdCoxian globalerrbkwdCoxian]
#         errSomePH = [errSomePH globalerrSomePH]
#         errSomePHBkwd = [errSomePHBkwd globalerrSomePHBkwd]
#         plot!(vecNBases,log.(globalerrPH),label="Erlang "*string(Δ),colour=c)
#         plot!(vecNBases,log.(globalerrDG),label="DG "*string(Δ),color=c,linestyle=:dash)
#         plot!(vecNBases,log.(globalerrCoxian),label="Coxian "*string(Δ),linestyle=:solid,colour=c+1,markershape=:+)
#         plot!(vecNBases,log.(globalerrbkwdCoxian),label="Cox. w rev. "*string(Δ),linestyle=:dash,colour=c+1,markershape=:+)
#         plot!(vecNBases,log.(globalerrSomePH),label="PH "*string(Δ),linestyle=:solid,colour=c+2)
#         plot!(vecNBases,log.(globalerrSomePHBkwd),label="PH. w rev. "*string(Δ),linestyle=:dash,colour=c+2)
#         plot!(legend=:outerright,xlabel="Order",ylabel="log(error)",legendtitle="Δ")
#         display(plot!(title="... PH,  -- Coxian,  -. bkwd Coxian,  Solid DG"))
#     end
#     # plot(log.(vecΔ)[:], log.(errPH[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dot,label=false,linewidth=2)
#     # plot!(log.(vecΔ)[:], log.(errCoxian[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dash,label=false)
#     # plot!(log.(vecΔ)[:], log.(errbkwdCoxian[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dashdot,label=false)
#     # plot!(log.(vecΔ)[:], log.(errSomePH[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dashdot,label=false)
#     # plot!(log.(vecΔ)[:], log.(errSomePHBkwd[:,2:end]'),colour=[1 2 3 4 5],linestyle=:dashdot,label=false)
#     # plot!(log.(vecΔ)[:], log.(errDG[:,2:end]'),colour=[1 2 3 4 5],label=vecNBases')
#     # plot!(legend=:bottomright,xlabel="log(Δ)",ylabel="log(error)",legendtitle="Order")
#     # plot!(title="... PH,  -- Coxian,  -. bkwd Coxian, Solid DG")
# end

#
# ME = MakeME(CMEParams[3],mean=2)#MakeSomeCoxian(21)#
# ~,~,πME = reversal(ME)
# γ = -1
# S = ME.Q-γ*I
# α = ME.α*(-S^-1)
# α = α./sum(α)
# MEpois = (α = α, Q = S)
# MEmix = (α = α, Q = ME.Q)
# αwhat = -ME.α*(ME.Q+γ*I)^-1
# αwhat = αwhat./sum(αwhat)
# MEwhat = (α = αwhat, Q = ME.Q)
# # A = sylvester(ME.Q,-ME.Q,γ*Matrix{Float64}(I(NBases)))
# t = 0:0.1:4
# densityME = density(t,ME)
# densityMEpois = density(t,MEpois)
# # densityMEmix = density(t,MEmix)
# densityMEwhat = density(t,MEwhat)
# plot(t,densityME)
# plot!(t,densityMEpois)
# # plot!(t,densityMEmix)
# plot!(t,densityMEwhat)
1
# a = CMEParams[3]["a"]
# b = CMEParams[3]["b"]
# c = CMEParams[3]["c"]
# omega = CMEParams[3]["omega"]
# kω = 1*omega
# Minv = [2*c*(1+kω^2) 0 0; 0 -b a; 0 a b]./(2*(1+kω^2))
# # Minv = [2*c*(1+kω^2) 0 0; 0 -a b; 0 -b -a]./(2*(1+kω^2))
# α, Q, q = MakeME(CMEParams[3], mean=CMEParams[3]["mu1"])
# println([q[1];q[3];q[2]]'*Minv)
# # println([q[1];-q[2];-q[3]]'*Minv)
# println(α)
#
# M = Minv^-1
#
# Fkkp1 = q*α*M
# display(Fkkp1)
# Fkkplus = -q*q'
# display(Fkkplus)
# display(-Q*M) # = (F+G)
# G = -Fkkplus-Q*M
# display(G) # G = -F + Q*M
# display(M)
#
# # check -(F+G)M⁻¹=Q
# display(Q)
# display(-(Fkkplus+G)*Minv)
#
# # phi at RHS
# display(q')
# # phi at LHS
# display(α*M)
# # inner products
# display(M)
#
# # for negative phases
# Fkkminus = (α*M)'*(α*M)
# display(Fkkminus)
# Q̃ = (G+Fkkminus)*Minv
# display(Q̃)
# display(eigen(Q̃).values)
#
#
# mesh = SFFM.MakeMesh(
#     model = model,
#     Nodes = collect(0:1:10),
#     NBases = 4,
#     Basis = "lagrange",
# )
#
# Mat = SFFM.MakeMatrices(mesh=mesh,model=model)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
# # Meridional diffusion

using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel, PrescribedVelocityFields

using Oceananigans.TurbulenceClosures: Horizontal
using Statistics
using JLD2
using Printf
using GLMakie

Nx = 360

# A spherical domain
grid = LatitudeLongitudeGrid(size = (Nx, 1, 1),
                             radius = 1,
                             latitude = (-60, 60),
                             longitude = (-180, 180),
                             z = (-1, 0))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = :c,
                                    velocities = PrescribedVelocityFields(), # quiescent
                                    closure = HorizontalScalarDiffusivity(κ=1),
                                    buoyancy = nothing)

# Tracer patch for visualization
Gaussian(λ, ϕ, L) = exp(-(λ^2 + ϕ^2) / 2L^2)

# Tracer patch parameters
L = 12 # degree

cᵢ(λ, ϕ, z) = Gaussian(λ, 0, L)

set!(model, c=cᵢ)

c = model.tracers.c

function progress(s)
    c = s.model.tracers.c
    @info "Maximum(c) = $(maximum(c)), time = $(s.model.clock.time) / $(s.stop_time)"
    return nothing
end

ϕᵃᶜᵃ_max = maximum(abs, ynodes(Center, grid))
Δx_min = grid.radius * cosd(ϕᵃᶜᵃ_max) * deg2rad(grid.Δλ)
Δy_min = grid.radius * deg2rad(grid.Δϕ)
Δ_min = min(Δx_min, Δy_min)

# Time-scale for gravity wave propagation across the smallest grid cell
cell_diffusion_time_scale = Δ_min^2

simulation = Simulation(model,
                        Δt = 0.1cell_diffusion_time_scale,
                        stop_time = 100cell_diffusion_time_scale,
                        iteration_interval = 100,
                        progress = progress)

output_fields = model.tracers

output_prefix = "longitudinal_tracer_diffusion_Nx$(grid.Nx)"

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = TimeInterval(cell_diffusion_time_scale),
                                                      filename = output_prefix,
                                                      overwrite_existing = true)

run!(simulation)

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

for iter in iterations
    c = file["timeseries/c/$iter"][:, 1, 1]
    @show maximum(c)
end

λ = xnodes(Center, grid)

iter = Node(0)

title = @lift "Tracer diffusion on a parallel, t = $(file["timeseries/t/" * string($iter)])"

c = @lift file["timeseries/c/" * string($iter)][:, 1, 1]

fig = Figure(resolution = (1080, 540))

ax = fig[1, 1] = Axis(fig, ylabel = "c(λ)", xlabel = "λ")

lines!(ax, λ, c, color=:black)

record(fig, "longitudinal_tracer_diffusion_Nx$(grid.Nx).mp4", iterations, framerate=30) do i
    @info "Animating iteration $i/$(iterations[end])..."
    iter[] = i
end


text:
using CUDArt,Knet

@knet function addtest(x; alpha=1, beta=1, init=nothing)
    b = par(;init=init)
    z = add(b,x; alpha=alpha, beta=beta)
end

@knet function multest(x; alpha=1, beta=1, init=nothing)
    w = par(;init=init)
    z = mul(w,x; alpha=alpha, beta=beta)
end

function rtest(;s1=nothing, s2=nothing, a=nothing, b=nothing, batchsize=8)
    global prog = (rand()<=0 ? addtest : multest)
    global dims1, dims2, ndims1, ndims2
    if s2 != nothing
        dims2 = s2
        ndims2 = length(s2)
    else
        ndims2 = 1
        while rand()<0.5 && ndims2<7
            ndims2 += 1
        end
        dims2 = ntuple(i->rand(1:10), ndims2)
    end
    if s1 != nothing
        dims1 = s1
        ndims1 = length(s1)
    else
        ndims1 = rand()<0.5 ? ndims2 : rand()<0.5 ? 1 : rand(1:ndims2)
        dims1 = ntuple(i->(rand()<0.5 ? dims2[i] : 1), ndims1)
        if prog==addtest && length(dims1)==1 # CUDNN_ADD_SAME_C exception
            dims1 = (length(dims2) == 1 ? dims2 : dims2[end:end])
        end
    end
    @show (dims1,dims2)
    global x = randn(dims2...,batchsize)
    global y = randn(dims2...,batchsize)
    global w = randn(dims1...)
    data = ItemTensor(x,y; batch=batchsize)

    global alpha = (a!=nothing ? a : rand()<0.5 ? 1 : rand()<0.5 ? -1 : rand()<0.5 ? randn() : 0)
    global beta  = (b!=nothing ? b : rand()<0.5 ? 1 : rand()<0.5 ? -1 : rand()<0.5 ? randn() : 0)
    @show (alpha,beta)
    global net = FNN(prog; alpha=alpha, beta=beta, init=w)
    @show net.net.op[end]

    gradcheck(net, data, quadloss)
end


text:
@testset "model exceptions" begin
    @test_throws SVDD.ModelStateException throw(SVDD.ModelStateException(SVDD.model_created, SVDD.model_fitted))
    @test_throws SVDD.ModelInvariantException throw(SVDD.ModelInvariantException("invariant violated"))
end


text:
module RuntimeEnv

# TODO (MAYBE): correct wrong check for path separators: Win ';' Unix: ':'
const PATH_SEPS = r";|:"

#  using LinearAlgebra
#  BLAS.vendor()
#  BLAS.get_num_threads()
#  BLAS.openblas_get_config()  (! if different vendor)
#
#  On Linux: result of `lscpu`

#region TomlDict

# Allows the keep the TOML table in the order they are declared
# TODO: change this DataStructures.OrderedDict
const IDict = Base.ImmutableDict
IDict{K,Any}(KV::Pair{K,V}) where {K,V} = IDict{K,Any}(KV[1], KV[2])
IDict{K,Any}(dict::IDict{K,Any}, KV::Pair{K,V}) where {K,V} = IDict{K,Any}(dict, KV[1], KV[2])
function IDict{K,Any}(ps::AbstractVector{Union{Nothing,Pair{K,V}}}) where {K,V}
    d = IDict{K,Any}()
    for i ∈ reverse(eachindex(ps))
        if !isnothing(ps[i])
            d = IDict{K,Any}(d, ps[i])
        end
    end
    return d
end
function IDict{K,Any}(ps::AbstractVector{Pair{K,V}}) where {K,V}
    d = IDict{K,Any}()
    for i ∈ reverse(eachindex(ps))
        d = IDict{K,Any}(d, ps[i])
    end
    return d
end
const TomlDict = IDict{String, Any}

#endregion

function julia_system()

    # $ export JULIA_LLVM_ARGS=--version
    # $ jr -e 1
    # LLVM (http://llvm.org/):
    #   LLVM version 11.0.1jl
    #   Optimized build.
    #   Default target: x86_64-w64-mingw32
    #   Host CPU: haswell
    #
    # $ export JULIA_LLVM_ARGS=--help
    # $ .... long list ....

    # loadavg                0 bytes typeof(Base.Sys.loadavg)
    # maxrss                 0 bytes typeof(Base.Sys.maxrss)
    # uptime                 0 bytes typeof(Base.Sys.uptime)
    # which                  0 bytes typeof(Base.Sys.which)

    mem_tot = round(Float64(Base.Sys.total_memory())/(1 << 30);digits=1)
    mem_free = round(Float64(Base.Sys.free_memory())/(1 << 30);digits=1)

    nfo = Sys.cpu_info() # Vector{Base.Sys.CPUinfo}, one by HT core

    lsb = "" # code below from stdlib\v1.6\...\InteractiveUtils.jl, versioninfo
    if Sys.islinux()
        try lsb = readchomp(pipeline(`lsb_release -ds`, stderr=devnull)); catch; end
    end
    if Sys.iswindows()
        try lsb = strip(read(`$(ENV["COMSPEC"]) /c ver`, String)); catch; end
    end
    if Sys.isunix()
        lsb *= " uname: " * readchomp(`uname -mprsv`)
    end

    TomlDict([
        "VERSION" => string(Base.VERSION);
        "BINDIR" => Sys.BINDIR;
        "STDLIB" => Sys.STDLIB;
        "JIT" => Sys.JIT;
        "CPU_NAME" => Sys.CPU_NAME;
        "CPU_THREADS" => Sys.CPU_THREADS;
        "WORD_SIZE" => Sys.WORD_SIZE;
        "KERNEL" => string(Sys.KERNEL);
        "ARCH" => string(Sys.ARCH);
        "MACHINE" => Sys.MACHINE;
        "PLATFORM" => lsb;
        "BUILD_STDLIB_PATH" => Sys.BUILD_STDLIB_PATH;
        "WINDOWS_VERSION" => string(Sys.windows_version());
        "PROCESS_TITLE" => Sys.get_process_title();
        "CPU_MODEL" => nfo[1].model;
        "CPU_SPEED" => "$(nfo[1].speed) MHz";
        "CPU_HTCORES" => string(length(nfo));
        "TOTAL_MEM" => "$mem_tot GB";
        "FREE_MEM" => "$mem_free GB";
        "LIBM" =>  Base.libm_name;
        "LLVM" => string("libLLVM-",Base.libllvm_version," (", Sys.JIT, ", ", Sys.CPU_NAME, ")");
    ])
end

function host_tools()
    function wch(name)
        try
            string(Sys.which(name))
        catch e
            string(e)
        end
    end
    TomlDict([
        "julia" => wch("julia");
        "ssh" => wch("ssh");
        "git" => wch("git");
        "cmd" => wch("cmd");
        "bash" => wch("bash");
        "code" => wch("code");
        "wt" => wch("wt");
        "winget" => wch("winget");
        "powershell" => wch("powershell");
    ])
end

function dictof(a::AbstractVector)
    IDict{String,Any}([string(i) => string(a[i]) for i ∈ eachindex(a)])
end

function pathlike(value::String)
    dictof(split(value, PATH_SEPS; keepempty=false))
end

function env_context()
    j_any = Set{String}()
    j_pat = Set{String}()
    paths = Set{String}()
    other = Set{String}()

    for var ∈ keys(ENV)
        if contains(var, "JULIA")
            if contains(ENV[var], PATH_SEPS)
                push!(j_pat, var)
            else
                push!(j_any, var)
            end
        elseif contains(ENV[var], PATH_SEPS)
            push!(paths, var)
        else
            push!(other, var)
        end
    end

    a = ["ENV:JULIA" => TomlDict([ var => ENV[var] for var ∈ sort(collect(j_any))])]
    b = ["ENV:(JULIAPATH):"*var => pathlike(ENV[var])  for var ∈ sort(collect(j_pat))]
    c = ["ENV(PATH):"*var => pathlike(ENV[var])  for var ∈ sort(collect(paths))]
    d = ["ENV:OTHERS" => TomlDict([ var => ENV[var] for var ∈ sort(collect(other))])]
    vcat(a,b,c,d)
end

function julia_loaded_modules()
    # Base.loaded_modules :: Dict{Base.PkgId, Module}
    # Base.PkgId <: Any
    #   uuid::Union{Nothing, Base.UUID}
    #   name::String
    id(x) = isnothing(x) ? "top-level" : string(x)
    ks = sort(collect(keys(Base.loaded_modules)); lt=(x,y)->x.name<y.name)
    TomlDict([k.name => id(k.uuid) for k ∈ ks])
end

function julia_options()
    fs = fieldnames(Base.JLOptions)
    op = Base.JLOptions() # Ignore field commands :: Ptr{Ptr{UInt8}}
    str(x::Ptr{<:Integer})= x == C_NULL ? "NULL" : unsafe_string(x)
    ptr(x) = eltype(typeof(x)) <: Integer ? str(x) : string(typeof(x))
    toml(x) = x isa Ptr ? ptr(x) : x isa Number ? x : "?"
    Dict( string(f) => toml(getproperty(op,f)) for f ∈ fs)
end

function full_context()

    # TODO ?
    # Base.DL_LOAD_PATH ->  -> https://docs.julialang.org/en/v1/stdlib/Libdl/

    a = [
        "PROGRAM_FILE" => PROGRAM_FILE;
        "PROGRAM_ARGS" => ARGS;
        "LOAD_PATH" => LOAD_PATH;
        "ACTIVE_PROJECT" => string(Base.ACTIVE_PROJECT[]);
        "SYSTEM" => julia_system();
        "LOAD_PATH_EXP" => dictof(Base.load_path());
        "DEPOT_PATH" => dictof(DEPOT_PATH);
        "WHICH" => host_tools();
    ]

    b = env_context()

    c = [
        "Base_JLOptions" => julia_options();
        "Base_loaded_modules" => julia_loaded_modules();
    ]

    TomlDict(vcat(a,b,c))
end

# using TOML
# function save(file::String)
#     open(file, "w") do io
#         TOML.print(io, full_context())
#     end
# end
# function load(file::String)
#     TOML.parsefile(file)
# end
# save("env.toml")
# TOML.print(load("env.toml"))

end # module RuntimeEnv


text:
__precompile__()
module NaNMath

using Compat

for f in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
    @eval begin
        ($f)(x::Float64) = ccall(($(string(f)),Base.Math.libm), Float64, (Float64,), x)
        ($f)(x::Float32) = ccall(($(string(f,"f")),Base.Math.libm), Float32, (Float32,), x)
        ($f)(x::Real) = ($f)(float(x))
        Compat.@dep_vectorize_1arg Number $f
    end
end

# Would be more efficient to remove the domain check in Base.sqrt(),
# but this doesn't seem easy to do.
sqrt(x::Real) = x < 0.0 ? NaN : Base.sqrt(x)

# Don't override built-in ^ operator
pow(x::Float64, y::Float64) = ccall((:pow,Base.Math.libm),  Float64, (Float64,Float64), x, y)
pow(x::Float32, y::Float32) = ccall((:powf,Base.Math.libm), Float32, (Float32,Float32), x, y)
pow(x::Number,y::Number) = pow(float(x),float(y))

"""
NaNMath.sum(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
*    Returns the sum of all elements in the array, ignoring NaN's.

##### Examples:
```julia
using NaNMath as nm
nm.sum([1., 2., NaN]) # result: 3.0
```
"""
function sum{T<:AbstractFloat}(x::Vector{T})
    if size(x)[1] == 0
        result = zero(eltype(x))
    else
        result = convert(eltype(x), NaN)
        for i in x
            if !isnan(i)
                if isnan(result)
                    result = i
                else
                    result += i
                end
            end
        end
    end

    if isnan(result)
        Base.warn_once("All elements of the array, passed to \"sum\" are NaN!")
    end
    return result
end

"""
NaNMath.maximum(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
*    Returns the maximum of all elements in the array, ignoring NaN's.

##### Examples:
```julia
using NaNMath as nm
nm.maximum([1., 2., NaN]) # result: 2.0
```
"""
function maximum{T<:AbstractFloat}(x::Vector{T})
    result = convert(eltype(x), NaN)
    for i in x
        if !isnan(i)
            if (isnan(result) || i > result)
                result = i
            end
        end
    end
    return result
end

"""
NaNMath.minimum(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
*    Returns the minimum of all elements in the array, ignoring NaN's.

##### Examples:
```julia
using NaNMath as nm
nm.minimum([1., 2., NaN]) # result: 1.0
```
"""
function minimum{T<:AbstractFloat}(x::Vector{T})
    result = convert(eltype(x), NaN)
    for i in x
        if !isnan(i)
            if (isnan(result) || i < result)
                result = i
            end
        end
    end
    return result
end

"""
NaNMath.mean(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
*    Returns the arithmetic mean of all elements in the array, ignoring NaN's.

##### Examples:
```julia
using NaNMath as nm
nm.mean([1., 2., NaN]) # result: 1.5
```
"""
function mean{T<:AbstractFloat}(x::Vector{T})
    return mean_count(x)[1]
end

"""
Returns a tuple of the arithmetic mean of all elements in the array, ignoring NaN's,
and the number of non-NaN values in the array.
"""
function mean_count{T<:AbstractFloat}(x::Vector{T})
    sum = convert(eltype(x), NaN)
    count = 0
    for i in x
        if !isnan(i)
            if isnan(sum)
                sum = i
                count = 1
            else
                sum += i
                count += 1
            end
        end
    end
    result = sum / count
    return (result, count)
end

"""
NaNMath.var(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
* Returns the sample variance of a vector A. The algorithm will return
  an estimator of the  generative distribution's variance under the
  assumption that each entry of v is an IID drawn from that generative
  distribution. This computation is  equivalent to calculating \\
  sum((v - mean(v)).^2) / (length(v) - 1). NaN values are ignored.

##### Examples:
```julia
using NaNMath as nm
nm.var([1., 2., NaN]) # result: 0.5
```
"""
function var{T<:AbstractFloat}(x::Vector{T})
    mean_val, n = mean_count(x)
    if !isnan(mean_val)
        sum_square = zero(eltype(x))
        for i in x
            if !isnan(i)
                sum_square += (i - mean_val)^2
            end
        end
        return sum_square / (n - one(eltype(x)))
    else
        return mean_val # NaN or NaN32
    end
end

"""
NaNMath.std(A)

##### Args:
* `A`: A one dimensional array of floating point numbers

##### Returns:
* Returns the standard deviation of a vector A. The algorithm will return
  an estimator of the  generative distribution's standard deviation under the
  assumption that each entry of v is an IID drawn from that generative
  distribution. This computation is  equivalent to calculating \\
  sqrt(sum((v - mean(v)).^2) / (length(v) - 1)). NaN values are ignored.

##### Examples:
```julia
using NaNMath as nm
nm.std([1., 2., NaN]) # result: 0.7071067811865476
```
"""
function std{T<:AbstractFloat}(x::Vector{T})
    return sqrt(var(x))
end

end


text:
using Keras
using Pyall
@pyimport numpy

model = Keras.load("structure.json", "weights.h5")


text:
using Test
import LanguageServer: Document, Range, Position,
                       get_text, get_offset, get_line_offsets, get_position_at,
                       get_open_in_editor, set_open_in_editor,
                       is_workspace_file,
                       applytextdocumentchanges, TextDocumentContentChangeEvent

text1 = """
123456
abcde
ABCDEFG
"""
doc1 = Document("untitled", text1, false)

text2 = """
12μ456
abηde
ABCDEFG"""
doc2 = Document("untitled", text2, true)

@testset "Basic operations" begin
    @testset "Document 1" begin

        @test get_text(doc1) == text1
        @test get_offset(doc1, 0, 4) == 4
        @test get_offset(doc1, 1, 2) == 9
        @test get_line_offsets(doc1) == [0, 7, 13, 21]
        @test get_position_at(doc1, 1) == (0, 1)
        @test get_position_at(doc1, 8) == (1, 1)
        @test get_position_at(doc1, 15) == (2, 2)

    end

    @testset "Dcoument 2" begin

        @test get_offset(doc2, 0, 4) == 5
        @test get_offset(doc2, 1, 2) == 10
        @test get_line_offsets(doc2) == [0, 8, 15]
        @test get_position_at(doc2, 1) == (0, 1)
        @test get_position_at(doc2, 9) == (1, 1)
        @test get_position_at(doc2, 17) == (2, 2)

    end

    @testset "Workspace" begin

        @test is_workspace_file(doc1) == false
        @test is_workspace_file(doc2) == true

        @test get_open_in_editor(doc1) == false
        set_open_in_editor(doc1, true)
        @test get_open_in_editor(doc1) == true
        set_open_in_editor(doc1, false)
        @test get_open_in_editor(doc1) == false

    end

    @testset "Document changes" begin

        new_text2 = """
        12μ456
        12abηde
        ABCDEFG"""
        applytextdocumentchanges(doc2, TextDocumentContentChangeEvent(Range(1), 0, "12"))

        @test get_text(doc2) == new_text2
        @test get_line_offsets(doc2) == [0, 8, 17]

    end

    @testset "Breaklines" begin

        text = "\n"
        breakline_doc = Document("untitled", text, false)
        @test get_line_offsets(breakline_doc) == [0,1]

        text = "1234\r\nabcd"
        breakline_doc = Document("untitled", text, false)
        @test_broken get_line_offsets(breakline_doc) == [0, 5]

        text = "1234\nabcd\n"
        breakline_doc = Document("untitled", text, false)
        @test get_line_offsets(breakline_doc) == [0, 5, 10]

    end

    @testset "whitespaces" begin

        text = " \n"
        whitespaces_doc = Document("untitled", text, false)
        @test get_line_offsets(whitespaces_doc) == [0,2]

        text = "\n "
        whitespaces_doc = Document("untitled", text, false)
        @test get_line_offsets(whitespaces_doc) == [0, 1]

        text = " \n \n "
        whitespaces_doc = Document("untitled", text, false)
        @test get_line_offsets(whitespaces_doc) == [0, 2, 4]

    end

end


text:
using BinDeps

@BinDeps.setup

libsteinspanner = library_dependency("libsteinspanner",
                                     aliases=["libstein_spanner.so", "libstein_spanner.1.so", "libstein_spanner.1.0.so", "libstein_spanner.1.0.0.so"])

provides(SimpleBuild,
   (@build_steps begin
        ChangeDirectory(joinpath(dirname(@__FILE__), "../src/discrepancy/spanner"))
        MakeTargets(["libstein_spanner.so"])
        `mkdir -p ../../../deps/usr/lib`
        `cp libstein_spanner.so ../../../deps/usr/lib`
    end), libsteinspanner, os = :Unix)

@BinDeps.install Dict(:libsteinspanner => :libsteinspanner)


text:
#
# CoolingScheduler
#

g = program(:linear, 100.0)
@test g(1) == 50.0
@test g(3) == 25.0
@test g(4) == 20.0
@test g(7) == 12.5

g = program(:quadratic, 100.0)
@test g(1) == 50.0
@test g(3) == 10.0
@test g(7) == 2.0
@test g(1) == 50.0

sch = CoolingScheduler()


text:
"""
Construct Toy Green's functions
Decaying peaks, control number of events, and their positions, depending on bfrac and efrac.
* `afrac` : control amplitude of peaks
"""
function toy_reflec_green!(x; bfrac=0.0, c=1.0, afrac=1.0)
	nt=size(x,1)
	nr=size(x,2)
	c=c*sqrt((nr^2)/(nt-1))

	it0=1+round(Int,bfrac*nt)
	for ir in 1:nr
		it=round(Int, it0+(ir*ir)*inv(c)*inv(c))
		if(it<=nt)
			x[it, ir]=1.0*afrac
		end
	end
	return x
end

function toy_direct_green!(x; bfrac=0.0, c=1.0, afrac=1.0)
	nt=size(x,1)
	nr=size(x,2)
	c=c*((nr)/(nt-1))

	it0=1+round(Int,bfrac*nt)
	for ir in 1:nr
		it=round(Int, it0+(ir)*inv(c))
		if(it<=nt)
			x[it, ir]=1.0*afrac
		end
	end
	return x

end


text:
# Autogenerated wrapper script for SCIP_jll for x86_64-linux-gnu-cxx11
export libscip

using bliss_jll
using GMP_jll
using CompilerSupportLibraries_jll
using Ipopt_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("SCIP")
JLLWrappers.@declare_library_product(libscip, "libscip.so.7.0")
function __init__()
    JLLWrappers.@generate_init_header(bliss_jll, GMP_jll, CompilerSupportLibraries_jll, Ipopt_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libscip,
        "lib/libscip.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##############################################################################
#
# Rayleigh distribution from Distributions Handbook
#
##############################################################################

immutable Rayleigh <: ContinuousUnivariateDistribution
    scale::Float64
    function Rayleigh(s::Real)
        s > zero(s) || error("scale must be positive")
        new(float64(s))
    end
    Rayleigh() = new(1.0)
end

cdf(d::Rayleigh, x::Real) = 1.0 - exp(-x^2 / (2.0 * d.scale^2))

entropy(d::Rayleigh) = 1.0 + log(d.scale) - log(sqrt(2.0)) - digamma(1.0) / 2.0

insupport(::Rayleigh, x::Real) = zero(x) < x < Inf
insupport(::Type{Rayleigh}, x::Real) = zero(x) < x < Inf

kurtosis(d::Rayleigh) = -(6.0 * pi^2 - 24.0 * pi + 16.0) / (4.0 - pi)^2

mean(d::Rayleigh) = d.scale * sqrt(pi / 2.)

median(d::Rayleigh) = d.scale * sqrt(2. * log(2.))

mode(d::Rayleigh) = d.scale
modes(d::Rayleigh) = [d.scale]

function pdf(d::Rayleigh, x::Real)
    insupport(d, x) ? (x / (d.scale^2)) * exp(-(x^2)/(2.0 * (d.scale^2))) : 0.0
end

quantile(d::Rayleigh, p::Real) = sqrt(-2.0 * d.scale^2 * log(1.0 - p))

rand(d::Rayleigh) = d.scale * sqrt(-2.0 * log(rand()))

skewness(d::Rayleigh) = (2.0 * sqrt(pi) * (pi - 3.0)) / (4.0 - pi)^1.5

var(d::Rayleigh) = d.scale^2 * (2.0 - pi / 2.0)


text:
function colmeans(mat::Array{T,2})::Array{Float64,1} where {T<:Number}
    _, p = size(mat)
    means = Array{Float64,1}(undef, p)
    for i = 1:p
        means[i] = Statistics.mean(mat[:, i])
    end
    return means
end

max_stars_repo_name: jbytecode/MCCGA

max_issues_repo_name: jbytecode/MCCGA

max_forks_repo_name: jbytecode/MCCGA


text:
export
    isprime, iscoprime, ismersenneprime, isrieselprime,
    primes, primesmask, primepi,
    nthprime, nextprime, prevprime,
    twinprimes

Nemo.isprime(n::UInt) = is_prime(n)
Nemo.isprime(n::Integer) = typemin(UInt) <= n <= typemax(UInt) ? (isprime ∘ UInt)(n) : (isprime ∘ fmpz)(n)
Nemo.isprime(n::Integer, mask) = n <= length(mask) ? mask[n] : isprime(n)

iscoprime(n::Integer, m::Integer) = n != 0 && m != 0 && gcd(n, m) == 1

Primes.primes(r::UnitRange{<:Integer}) = primes(first(r), last(r))

Primes.primesmask(r::UnitRange{<:Integer}) = primesmask(first(r), last(r))

primepi(a::Integer, b::Integer) = (count ∘ primesmask)(a, b)
primepi(r::UnitRange{<:Integer}) = primepi(first(r), last(r))
primepi(n::Integer) = primepi(1, n)

nthprime(a::Integer, b::Integer) = (n = ceil(typeof(b), b*log(b+2) + b*(log∘log)(b+2)); primes(n)[a:b])
nthprime(r::UnitRange{<:Integer}) = nthprime(first(r), last(r))
nthprime(n::Integer) = nthprime(n, n)[1]

# https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl
nextprime(n::T) where {T<:Integer} = begin
    n <= 1 && return convert(T, 2)
    n == 2 && return convert(T, 3)

    n += iseven(n) ? one(T) : convert(T, 2)
    isprime(n) && return n

    m = mod(n, 3)
    if m == 1
        a = convert(T, 4); b = convert(T, 2)
    elseif m == 2
        a = convert(T, 2); b = convert(T, 4)
    else
        n += 2
        a = convert(T, 2); b = convert(T, 4)
    end

    p = n
    while !isprime(p)
        p += a
        isprime(p) && break
        p += b
    end
    p
end

# https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl
prevprime(n::T) where {T<:Integer} = begin
    n <= 2 && return zero(T)
    n == 3 && return convert(T, 2)

    n -= iseven(n) ? one(T) : convert(T, 2)
    isprime(n) && return n

    m = mod(n, 3)
    if m == 1
        a = convert(T, 2); b = convert(T, 4)
    elseif m == 2
        a = convert(T, 4); b = convert(T, 2)
    else
        n -= 2
        a = convert(T, 2); b = convert(T, 4)
    end

    p = n
    while !isprime(p)
        p -= a
        isprime(p) && break
        p -= b
    end
    p
end

##  Find all twin primes
# https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl
twinprimes(a::Integer, b::Integer) = begin
    P = primes(a, b)
    inds = find(diff(P) .== 2)
    hcat(P[inds], P[inds+1])
end
twinprimes(r::UnitRange{<:Integer}) = twinprimes(first(r), last(r))