text:
module ScenTreesMakie

using Plots
import Plots.plot
import Plots.plot!

using StatsBase
using Distributions, Random
rng = MersenneTwister(01012019);
using Statistics, LinearAlgebra
using ProgressMeter
using JuMP, Clp
gr()

include("TreeStructure.jl")
include("LatticeStructure.jl")
include("Approx.jl")
include("StochPaths.jl")

include("KernelDensityEstimation.jl")
include("plots.jl")
include("helper.jl")
include("Wasserstein.jl")

export Tree, Lattice,
        tree_approximation!,lattice_approximation,

        stage,height,leaves,nodes,root, children, part_tree,build_probabilities!, checkTree,
        structure, part_lattice, combine_tree,
        gaussian_path, running_maximum, kernel_scenarios,
        sample_path,  plot_path!, plot, plot!,

        path_ident,

        ontoSimplex!, summary,

        distFunction, Wasserstein, Sinkhorn, nestedWasserstein

end


text:
ENV["GKSwstype"]="100" #to enable saving figure in server

using DataFrames,CSV,Random,DelimitedFiles,Statistics
using Plots

Random.seed!(123)
methodALL=["omics_gblup","omics_bayescpi"]
main_path="/group/qtlchenggrp/tianjing/res_andrew_1000snp/missing_pattern_legarra/"
pctAll=["0.1","0.3","0.5","0.7","0.9"]
nRep=20
chain_length=5000
#read bv
bv=vec(readdlm("/group/qtlchenggrp/tianjing/omics_data/andrew_1000snp/data1/bv.txt"));

for method in methodALL
    for pct in pctAll
        cd(main_path*"$pct/$method")
        final_accuracy_all=ones(nRep,2)*999.0

        for i in 1:nRep
            #read results for repi
            if ispath("rep$i/results/MCMC_samples_EBV_NonLinear.txt")
                mcmc=readdlm("rep$i/results/MCMC_samples_EBV_NonLinear.txt",',')[2:end,:]
                if size(mcmc,1)==1000
                    mcmc_mean=zeros(size(mcmc))
                    for i in 1:1000
                        mcmc_mean[i,:]=mean(mcmc[1:i,:],dims=1)
                    end
                    accuracy_datai=zeros(1000)
                    for i in 1:1000
                        accuracy_datai[i]=cor(mcmc_mean[i,:][956:1055],bv[956:1055])
                    end

                    #save results
                    open(main_path*"$pct/accuracy_$method/accuracy.$method.rep$i.txt", "w") do io
                        writedlm(io, accuracy_datai)
                    end
                    final_accuracy=round(accuracy_datai[end],digits=3)
                    final_accuracy_all[i,:]=[i,accuracy_datai[end]]

                    #plot results
                    myfig=plot(collect(chain_length/1000:chain_length/1000:chain_length),accuracy_datai,xlabel="iteration",ylabel="accuracy",title="ind-level missing.pct$pct.rep$i.accu=$final_accuracy")
                    savefig(myfig,main_path*"$pct/accuracy_$method/accuracy.$method.pct$pct.rep$i.png")
                end
            end
        end
        open(main_path*"$pct/accuracy_$method/accuracy_final.$method.pct$pct.txt", "w") do io
            writedlm(io, final_accuracy_all)
        end
    end
end

ext: jl

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
using ImplicitGraphs, Primes

"""
    iPaley(p::Int)

Create an implicit Paley graph on `p` vertices.
Here `p` must be a prime congruent to 1 modulo 4.
The vertices of the graph are integers in the range
`0:p-1` and two vertices are adjacent if their difference
is a quadratic residue mod `p`.
"""
function iPaley(p::Int)::ImplicitGraph
    if !isprime(p) || p % 4 != 1
        error("Argument ($p) must be prime and congruent to 1 mod 4")
    end

    vcheck(v::Int) = (0 <= v < p)

    function outs(v::Int)
        return unique((v + k^2) % p for k = 1:p-1)
    end

    return ImplicitGraph{Int}(vcheck, outs)
end


text:
function harmonic_mean(x...)
    sum = 0
    for i = 1:length(x)
        sum = 1/x[i] + sum
    end
    mean = length(x) / sum
    return mean
end

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ATK_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ATK")
JLLWrappers.@generate_main_file("ATK", UUID("7b86fcea-f67b-53e1-809c-8f1719c154e8"))
end  # module ATK_jll


text:
"""
    function addCheckBounds_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with boundary checking.
"""
function addCheckBounds_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if !isempty([i for i in keys(p.update["Local"])])
        ##Add bound code if bound symbols update
        code = returnBound_(space.box,p)

        if !emptyquote_(code)

            #Construct functions
            f = simpleFirstLoopWrapInFunction_(platform,:boundsCheck_!,code)
            f = vectorize_(abm,f,p)

            push!(p.declareF.args,
                f)

            push!(p.execInloop.args,
                    :(@platformAdapt boundsCheck_!(ARGS_))
                )

        end

    end

    return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using BHAPtfem
using Test

code_tests = [
  #"ex01.jl",
  "materialtable.jl",
  "mediatable.jl",
  #"createmesh.jl",
  #"fisher_a1_01.jl",
  "tp01.jl"
]

println("\n\nRunning BottomHoleAssemblyAnalysis/BHatp.jl tests:\n\n")

@testset "BHAPtfem.jl" begin

  for test in code_tests
      println("\n\n  * $(test) *\n")
      include("test_"*test)
      println()
  end

end


text:
using Revise
using ADCME
using PyPlot
matplotlib.use("macosx")
X,_,_,_ = mnist()
X = X/127.5 .-1.
Y = zeros(size(X,1), 28*28)
for i = 1:size(X,1)
    Y[i,:] = vec(X[i,:,:,1])
end
X = copy(Y)

function discriminator(x, gan)
    y = ae(x, [256,128,64,1])
    y = sigmoid(y)
    return y
end

function generator(x, gan)
    y = ae(x, [128,256,784])
    return y
end

# function loss(t, f, gan)
#     t_logits = discriminator(t, gan)
#     f_logits = discriminator(f, gan)
#     t_loss = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=t_logits, labels=tf.ones_like(t_logits)))
#     f_loss = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=f_logits, labels=tf.zeros_like(f_logits)))
#     d_loss = t_loss + f_loss
#     g_loss = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=f_logits, labels=tf.ones_like(f_logits)))
#     d_loss, g_loss
# end

reset_default_graph()
gan = GAN(X, generator, discriminator, "klgan"; latent_dim=100, batch_size=32)
opt_d = RMSPropOptimizer(5e-5).minimize(gan.d_loss, var_list=gan.d_vars)
opt_g = RMSPropOptimizer(5e-5).minimize(gan.g_loss, var_list=gan.g_vars)
img = sample(gan, 25)

# if gantype=="wgan"
#     op = clip(gan.d_vars[1:2:end], -0.001, 0.001)
#     [push!(gan.update, o) for o in op]
# end
sess = Session(); init(sess)
rm("images/", force=true, recursive=true)
mkdir("images/")
for i = 1:200000
    feed = Dict(gan.noise=>randn(gan.batch_size, gan.latent_dim), gan.ids=>Int32.(rand(1:size(X,1), gan.batch_size)))
    _, dl, _ = run(sess, [opt_d, gan.d_loss,gan.update], feed_dict=feed)
    for i = 1:4
        run(sess, opt_g, feed_dict=feed)
    end
    _, gl = run(sess, [opt_g, gan.g_loss], feed_dict=feed)
    if i%100==1
        println("iter=$i, dloss=$dl, gloss=$gl")
    end
    if i%500==1
        val = run(sess, img)
        fig, axs = subplots(5,5)
        cnt = 1
        for i = 1:5
            for j = 1:5
                axs[i,j].imshow(reshape(val[cnt,:], 28, 28), cmap="gray")
                axs[i,j].axis("off")
                cnt += 1
            end
        end
        savefig("images/$i.png")
        close("all")
    end
end


text:
using Cuba

push!(LOAD_PATH,"src")
using DalitzPlotAnalysis

using amplitudes_compass

function the_big_integral(s, i::Int64)
  function integrand(x,f)
    s1min = 4.0 * mπ2
    s1max = (sqrt(s) - mπ)^2

    s1 =  s1min + x[1] * (s1max - s1min)
    cosθ1 = -1.0 + 2.0*x[2]
    ϕ1 = -π + 2.0*π*x[3]
    cosθ23 = -1.0 + 2.0*x[4]
    ϕ23 = -π + 2.0*π*x[5]
    a = COMPASS_wave(i, s, s1, cosθ1, ϕ1, cosθ23, ϕ23)
    λλ = λ(s,s1,mπ2)*λ(s1,mπ2,mπ2)
    (λλ < 0) && error("λλ < 0: check masses")
    f[1] = abs2(a) * sqrt(λ(s,s1,mπ2)) * sqrt(λ(s1,mπ2,mπ2)) / (s*s1)*(s1max-s1min)
  end
  return cuhre(integrand, 5, 1)
end

# index 2 for ρπ S-wave
the_big_integral(1.1,2)

function the_sml_integral(s, i::Int64)
  function integrand(x,f)
    s1min = 4.0 * mπ2
    s1max = (sqrt(s) - mπ)^2

    s1 =  s1min + x[1] * (s1max - s1min)
    z = -1.0 + 2.0*x[2]
    ampl2 = 0.0
    for lm = -1:1
        ampl2 += abs2(COMPASS_wave_short(i, lm, s, s1, z))
    end
    f[1] = ampl2 * sqrt(λ(s,s1,mπ2)) * sqrt(λ(s1,mπ2,mπ2)) / (s*s1)*(s1max-s1min)
  end
  return cuhre(integrand, 2, 1)
end

the_small_integral(1.1,2)[1][1]

##########################
@time [the_big_integral(1.1,2)[1][1] for i in 1:3]
@time [the_sml_integral(1.1,2)[1][1] for i in 1:3]

@profile [the_big_integral(1.1,2)[1][1] for i=1:2]
Juno.profiler()

#########################
using Plots
theme(:juno)
let ev = linspace(0.5,2.5,10)
    big_cal = [the_big_integral(e^2,2)[1][1] for e in ev]
    sht_cal = [the_sml_integral(e^2,2)[1][1] for e in ev]
     plot(ev, big_cal, lab = "the 5d-integral")
    plot!(ev, sht_cal, lab = "the 2d-integral")
end

#########################
let ev = linspace(0.5,2.5,50)
    sht_cal = [the_sml_integral(e^2,2)[1][1] for e in ev]
    plot(ev, sht_cal, lab = "the 2d-integral")
end


text:
using RNG.Xorshifts

seed = RNG.gen_seed(UInt64, 2)
r = Xoroshiro128(seed)

r1 = WrappedRNG(r, UInt32)
r2 = WrappedRNG(Xoroshiro128, UInt32, seed)

@test rand(r1, UInt64) == rand(r2, UInt64)
@test rand(r1, UInt32) == rand(r2, UInt32)

rand(r2, UInt32)

r3 = WrappedRNG(r, UInt128)
@test rand(r3, UInt128) == rand(r2, UInt128)


text:
using DelayDiffEq
using PGFPlotsX

# create Mackey-Glass equation
f_mackey_glass(u, h, p, t) = 2 * h(p, t - 2) / (1 + h(p, t - 2)^9.65) - u

"""
    plot_mackey_glass(u₀, ts)

Plot numerical solutions of the Mackey-Glass equation
```math
u'(t) = 2 u(t - 2) / (1 + u(t - 2)^9.65) - u(t)
```
for ``t ∈ ts`` with ``a ∈ as`` and ``u(t) = u₀`` for ``t ≤ min ts``.
"""
function plot_mackey_glass(u₀, ts)
    # obtain time points
    t₀, t₁ = extrema(ts)

    # define history function
    h(p, t) = u₀

    # define DDE problem
    prob = DDEProblem(f_mackey_glass, u₀, h, (t₀, t₁); constant_lags = [2])

    # define DDE solver
    alg = MethodOfSteps(Rosenbrock23())

    # compute solution
    sol = solve(prob, alg, saveat = ts)

    # plot solution
    @pgf ax = Axis(
        {
            xlabel = raw"$t$", ylabel = raw"$x(t)$",
            xmin = t₀, xmax = t₁,
            no_marks, grid = "major"
        },
        PlotInc(Table(; x = ts, y = sol.u))
    )

    ax
end

"""
    plot_mackey_glass_embedding(u₀, ts)

Plot time delay embeddings, i.e., plots of ``u(t)`` versus ``u(t - 2)``,
of the Mackey-Glass equation
```math
u'(t) = 2 u(t - 2) / (1 + u(t - 2)^9.65) - u(t)
```
for ``t ∈ ts`` with ``a ∈ as`` and ``u(t) = u₀`` for ``t ≤ min ts``.
"""
function plot_mackey_glass_embedding(u₀, ts)
    # obtain time points
    t₀, t₁ = extrema(ts)

    # compute time points for time delay embedding
    Ts = filter(x -> x ≥ t₀ + 2, ts)

    # define history function
    h(p, t) = u₀

    # define DDE problem
    prob = DDEProblem(f_mackey_glass, u₀, h, (t₀, t₁); constant_lags = [2])

    # define DDE solver
    alg = MethodOfSteps(Tsit5())

    # compute solution
    sol = solve(prob, alg)

    # plot solution
    @pgf ax = Axis(
        {
            xlabel = raw"$x(t)$", ylabel = raw"$x(t - \tau)$",
            no_marks, grid = "major"
        },
        PlotInc(Table(; x = sol(Ts), y = sol(Ts .- 2)))
    )

    ax
end

pgfsave(joinpath(@__DIR__, "..", "figures", "generated", "mackey_glass.tex"),
        plot_mackey_glass(0.5, range(0; stop = 600, length = 10_000));
        include_preamble = false)
pgfsave(joinpath(@__DIR__, "..", "figures", "generated", "mackey_glass_embedding.tex"),
        plot_mackey_glass_embedding(0.5, range(0; stop = 600, length = 10_000));
        include_preamble = false)


text:
# Store an alias table
immutable DiscreteDistributionTable
	table::Vector{Vector{Int64}}
	bounds::Vector{Int64}
end

# TODO: Test if bit operations can speed up Base64 mod's and fld's
function DiscreteDistributionTable{T <: Real}(probs::Vector{T})
	# Cache the cardinality of the outcome set
	n = length(probs)

	# Convert all Float64's into integers
	vals = Array(Int64, n)
	for i in 1:n
		vals[i] = int(probs[i] * 64^9)
	end

	# Allocate digit table and digit sums as table bounds
	table = Array(Vector{Int64}, 9)
	bounds = zeros(Int64, 9)

	# Special case for deterministic distributions
	for i in 1:n
		if vals[i] == 64^9
			table[1] = Array(Int64, 64)
			for j in 1:64
				table[1][j] = i
			end
			bounds[1] = 64^9
			for j in 2:9
				table[j] = Array(Int64, 0)
				bounds[j] = 64^9
			end
			return DiscreteDistributionTable(table, bounds)
		end
	end

	# Fill tables
	multiplier = 1
	for index in 9:-1:1
		counts = Array(Int64, 0)
		for i in 1:n
			digit = mod(vals[i], 64)
			# vals[i] = fld(vals[i], 64)
			vals[i] >>= 6
			bounds[index] += digit
			for itr in 1:digit
				push!(counts, i)
			end
		end
		bounds[index] *= multiplier
		table[index] = counts
		# multiplier *= 64
		multiplier <<= 6
	end

	# Make bounds cumulative
	bounds = cumsum(bounds)

	return DiscreteDistributionTable(table, bounds)
end

function draw(table::DiscreteDistributionTable)
	# 64^9 - 1 == 0x003fffffffffffff
	i = rand(1:(64^9 - 1))
	# if i == 64^9
	# 	return table.table[9][rand(1:length(table.table[9]))]
	# end
	bound = 1
	while i > table.bounds[bound] && bound < 9
		bound += 1
	end
	if bound > 1
		index = fld(i - table.bounds[bound - 1] - 1, 64^(9 - bound)) + 1
	else
		index = fld(i - 1, 64^(9 - bound)) + 1
	end
	return table.table[bound][index]
end


text:
using NeuralAttentionlib
using Test
using Flux
using MacroTools: @forward, @capture

# taken from
# https://github.com/chengchingwen/NeuralAttentionlib.jl/blob/master/test/mha.jl

const Abstract3DTensor{T} = AbstractArray{T, 3}

macro toNd(ex, outref::Int=1)
    if @capture ex f_(xs__; kw__)
        kwe = Expr(:parameters, kw...)
    else
        @capture ex f_(xs__)
    end
    rxs = map(xs) do x
        :(reshape($x, size($x, 1), :))
    end
    newcall = kw === nothing ? Expr(:call, f, rxs...) : Expr(:call, f, kwe, rxs...)
    :(reshape($newcall, :, Base.tail(size($(xs[outref])))...)) |> esc
end

abstract type AbstractAttention end

function atten(mha, x, mask=nothing)
    q = mha.iqproj(x)
    k = mha.ikproj(x)
    v = mha.ivproj(x)
    a = NeuralAttentionlib.multihead_qkv_attention(4, q, k, v, mask)
    return mha.oproj(a)
end

struct MultiheadAttention{Q<:Dense, K<:Dense, V<:Dense, O<:Dense, DP<:Dropout} <: AbstractAttention
    head::Int
    future::Bool
    iqproj::Q
    ikproj::K
    ivproj::V
    oproj::O
    drop::DP
end

Flux.functor(mh::MultiheadAttention) = (mh.iqproj, mh.ikproj, mh.ivproj, mh.oproj), m -> MultiheadAttention(mh.head, mh.future, m..., mh.drop)

"""
    MultiheadAttention(head::Int, is::Int, hs::Int, os::Int;
                       future::Bool=true, pdrop = 0.1)
Multihead dot product Attention Layer, `head` is the number of head,
`is` is the input size, `hs` is the hidden size of input projection layer of each head,
`os` is the output size. When `future` is `false`, the k-th token can't see tokens at > k.
`pdrop` is the dropout rate.
"""
MultiheadAttention(head::Int,
                   is::Int,
                   hs::Int,
                   os::Int;
                   future::Bool=true, pdrop = 0.1) = MultiheadAttention(head,
                                                                        future,
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(hs*head, os),
                                                                        Dropout(pdrop),
                                                                        )

function (mh::MultiheadAttention)(query::A1,
                                  key::A2,
                                  value::A3;
                                  mask=nothing) where {T,
                                                       A1 <: Abstract3DTensor{T},
                                                       A2 <: Abstract3DTensor{T},
                                                       A3 <: Abstract3DTensor{T}}
  qs = size(query)
  ks = size(key)
  vs = size(value)

  #size(ipq) == (h, q_seq_len, batch)
  ipq = @toNd mh.iqproj(query)
  ipk = @toNd mh.ikproj(key)
  ipv = @toNd mh.ivproj(value)

  h = size(ipq, 1)
  hs = div(h, mh.head)

  #size(ipq) == (hs, q_seq_len, head, batch)
  ipq = permutedims(reshape(ipq, hs, mh.head, qs[2], qs[3]), [1, 3, 2, 4])
  ipk = permutedims(reshape(ipk, hs, mh.head, ks[2], ks[3]), [1, 3, 2, 4])
  ipv = permutedims(reshape(ipv, hs, mh.head, vs[2], vs[3]), [1, 3, 2, 4])

  #size(ipq) == (hs, q_seq_len, head * batch)
  ipq = reshape(ipq, hs, qs[2], :)
  ipk = reshape(ipk, hs, ks[2], :)
  ipv = reshape(ipv, hs, vs[2], :)

  atten = attention(ipq,ipk,ipv,
                    mask,
                    mh.future,
                    mh.drop)

  atten = permutedims(reshape(atten, hs, qs[2], mh.head, qs[3]), [1, 3, 2, 4]) #size(atten) == (hs, head, ql, b)
  atten = reshape(atten, h, qs[2], qs[3]) #size(atten) == (h, ql, b)

  out = @toNd mh.oproj(atten)
  out #size(out) == (h, q_seq_len, batch)
end


text:
struct IIF1MConstantCache{vecuType,rhsType,nl_rhsType} <: StochasticDiffEqConstantCache
  uhold::vecuType
  rhs::rhsType
  nl_rhs::nl_rhsType
end

@cache struct IIF1MCache{uType,vecuType,DiffCacheType,rhsType,nl_rhsType,rateType,rateNoiseType,rateNoiseCollectionType,NoiseTmpType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uhold::vecuType
  dual_cache::DiffCacheType
  tmp::uType
  rhs::rhsType
  nl_rhs::nl_rhsType
  rtmp1::rateType
  rtmp2::rateNoiseType
  rtmp3::rateNoiseCollectionType
  noise_tmp::NoiseTmpType
end

function alg_cache(alg::IIF1M,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  uhold = Vector{typeof(u)}(undef, 1)
  tmp = zero(rate_prototype)
  rhs = RHS_IIF1M_Scalar(f,tmp,t,t,p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  IIF1MConstantCache(uhold,rhs,nl_rhs)
end

function alg_cache(alg::IIF1M,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  tmp = similar(u,axes(u)); rtmp1 = zero(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
  uhold = vec(u) # this makes uhold the same values as integrator.u
  rhs = RHS_IIF1(f,tmp,t,t,dual_cache,size(u),p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  noise_tmp = tmp

  rtmp2 = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    rtmp3 = rtmp2
  else
    rtmp3 = zero(rate_prototype)
  end
  IIF1MCache(u,uprev,uhold,dual_cache,tmp,rhs,nl_rhs,rtmp1,rtmp2,rtmp3,noise_tmp)
end

struct IIF2MConstantCache{vecuType,rhsType,nl_rhsType} <: StochasticDiffEqConstantCache
  uhold::vecuType
  rhs::rhsType
  nl_rhs::nl_rhsType
end

@cache struct IIF2MCache{uType,vecuType,DiffCacheType,rhsType,nl_rhsType,rateType,rateNoiseType,rateNoiseCollectionType,NoiseTmpType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uhold::vecuType
  dual_cache::DiffCacheType
  tmp::uType
  rhs::rhsType
  nl_rhs::nl_rhsType
  rtmp1::rateType
  rtmp2::rateNoiseType
  rtmp3::rateNoiseCollectionType
  noise_tmp::NoiseTmpType
end

function alg_cache(alg::IIF2M,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  uhold = Vector{typeof(u)}(undef, 1)
  tmp = zero(rate_prototype)
  rhs = RHS_IIF2M_Scalar(f,tmp,t,t,p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  IIF2MConstantCache(uhold,rhs,nl_rhs)
end

function alg_cache(alg::IIF2M,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  tmp = similar(u,axes(u)); rtmp1 = zero(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
  uhold = vec(u) # this makes uhold the same values as integrator.u
  rhs = RHS_IIF2(f,tmp,t,t,dual_cache,size(u),p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  noise_tmp = tmp

  rtmp2 = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    rtmp3 = rtmp2
  else
    rtmp3 = zero(rate_prototype)
  end
  IIF2MCache(u,uprev,uhold,dual_cache,tmp,rhs,nl_rhs,rtmp1,rtmp2,rtmp3,noise_tmp)
end

struct IIF1MilConstantCache{vecuType,rhsType,nl_rhsType} <: StochasticDiffEqConstantCache
  uhold::vecuType
  rhs::rhsType
  nl_rhs::nl_rhsType
end

@cache struct IIF1MilCache{uType,vecuType,DiffCacheType,rhsType,nl_rhsType,rateType,rateNoiseType,rateNoiseCollectionType,NoiseTmpType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uhold::vecuType
  dual_cache::DiffCacheType
  tmp::uType
  rhs::rhsType
  nl_rhs::nl_rhsType
  rtmp1::rateType
  rtmp2::rateNoiseType
  rtmp3::rateNoiseCollectionType
  noise_tmp::NoiseTmpType
  gtmp::rateNoiseType
  gtmp2::rateNoiseType
end

function alg_cache(alg::IIF1Mil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  uhold = Vector{typeof(u)}(undef, 1)
  C = zero(rate_prototype)
  rhs = RHS_IIF1M_Scalar(f,C,t,t)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  IIF1MilConstantCache(uhold,rhs,nl_rhs)
end

function alg_cache(alg::IIF1Mil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  tmp = similar(u,axes(u)); rtmp1 = zero(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
  uhold = vec(u) # this makes uhold the same values as integrator.u
  rhs = RHS_IIF1(f,tmp,t,t,dual_cache,size(u))
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  noise_tmp = zero(noise_rate_prototype)
  gtmp = zero(noise_rate_prototype); gtmp2 = zero(noise_rate_prototype)
  rtmp2 = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    rtmp3 = rtmp2
  else
    rtmp3 = zero(rate_prototype)
  end
  IIF1MilCache(u,uprev,uhold,dual_cache,tmp,rhs,nl_rhs,rtmp1,rtmp2,rtmp3,noise_tmp,gtmp,gtmp2)
end


text:
function wordcount(sentence::AbstractString)

end


text:
##############################################################################
##
## FixedEffectCoefficients : vector x in A'Ax = A'b
##
## We define these methods used in lsmr! (duck typing):
## copyto!, fill!, rmul!, axpy!, norm
##
## Do not define iteration on each fixedeffect since it would conflict with eltype
##
##############################################################################

struct FixedEffectCoefficients{T}
	x::Vector{<:AbstractVector{T}}
end

Base.eltype(fecoef::FixedEffectCoefficients{T}) where {T} = T
Base.length(fecoef::FixedEffectCoefficients) = sum(length(x) for x in fecoef.x)
LinearAlgebra.norm(fecoef::FixedEffectCoefficients) = sqrt(sum(sum(abs2, x) for x in fecoef.x))

function Base.fill!(fecoef::FixedEffectCoefficients, α::Number)
	for x in fecoef.x
		fill!(x, α)
	end
	return fecoef
end

function LinearAlgebra.rmul!(fecoef::FixedEffectCoefficients, α::Number)
	for x in fecoef.x
		rmul!(x, α)
	end
	return fecoef
end

function Base.copyto!(fecoef1::FixedEffectCoefficients, fecoef2::FixedEffectCoefficients)
	for (x1, x2) in zip(fecoef1.x, fecoef2.x)
		copyto!(x1, x2)
	end
	return fecoef1
end

function LinearAlgebra.axpy!(α::Number, fecoef1::FixedEffectCoefficients, fecoef2::FixedEffectCoefficients)
	for (x1, x2) in zip(fecoef1.x, fecoef2.x)
		axpy!(α, x1, x2)
	end
	return fecoef2
end

##############################################################################
##
## FixedEffectLinearMap
##
## Model matrix of categorical variables
## mutiplied by diag(1/sqrt(∑w * interaction^2, ..., ∑w * interaction^2) (Jacobi preconditoner)
##
## We define these methods used in lsmr! (duck typing):
## eltyp
## size
## mul!
##
##############################################################################

struct FixedEffectLinearMap{T}
	fes::Vector{<:FixedEffect}
	sqrtw::AbstractVector{T}
	colnorm::Vector{<:AbstractVector}
	caches::Vector{<:AbstractVector}
end

function FixedEffectLinearMap{T}(fes::Vector{<:FixedEffect}, weights::AbstractVector, ::Type{Val{:lsmr}}) where {T}
	sqrtw = convert(AbstractVector{T}, sqrt.(values(weights)))
	colnorm = [colnorm!(zeros(T, fe.n), fe.refs, fe.interaction, sqrtw) for fe in fes]
	caches = [cache!(zeros(T, length(sqrtw)), fe.interaction, sqrtw, scale, fe.refs) for (fe, scale) in zip(fes, colnorm)]
	return FixedEffectLinearMap{T}(fes, sqrtw, colnorm, caches)
end

function colnorm!(fecoef::AbstractVector, refs::AbstractVector, interaction::AbstractVector, sqrtw::AbstractVector)
	@inbounds @simd for i in eachindex(refs)
		fecoef[refs[i]] += abs2(interaction[i] * sqrtw[i])
	end
	fecoef .= sqrt.(fecoef)
end

function cache!(y::AbstractVector, interaction::AbstractVector, sqrtw::AbstractVector, fecoef::AbstractVector, refs::AbstractVector)
	@inbounds @simd for i in eachindex(y)
		y[i] = interaction[i] * sqrtw[i] / fecoef[refs[i]]
	end
	return y
end

LinearAlgebra.adjoint(fem::FixedEffectLinearMap) = Adjoint(fem)

function Base.size(fem::FixedEffectLinearMap, dim::Integer)
	(dim == 1) ? length(fem.fes[1].refs) : (dim == 2) ? sum(fe.n for fe in fem.fes) : 1
end
Base.eltype(x::FixedEffectLinearMap{T}) where {T} = T

function LinearAlgebra.mul!(y::AbstractVector, fem::FixedEffectLinearMap,
			  fecoefs::FixedEffectCoefficients, α::Number, β::Number)
	rmul!(y, β)
	for (fecoef, fe, cache) in zip(fecoefs.x, fem.fes, fem.caches)
		demean!(y, α, fecoef, fe.refs, cache)
	end
	return y
end

function demean!(y::AbstractVector, α::Number, fecoef::AbstractVector, refs::AbstractVector, cache::AbstractVector)
	@simd ivdep for i in eachindex(y)
		@inbounds y[i] += α * fecoef[refs[i]] * cache[i]
	end
end

function LinearAlgebra.mul!(fecoefs::FixedEffectCoefficients, Cfem::Adjoint{T, FixedEffectLinearMap{T}},
				y::AbstractVector, α::Number, β::Number) where {T}
	fem = adjoint(Cfem)
	rmul!(fecoefs, β)
	for (fecoef, fe, cache) in zip(fecoefs.x, fem.fes, fem.caches)
		mean!(fecoef, fe.refs, α, y, cache)
	end
	return fecoefs
end

function mean!(fecoef::AbstractVector, refs::AbstractVector, α::Number, y::AbstractVector, cache::AbstractVector)
	@simd ivdep for i in eachindex(y)
		@inbounds fecoef[refs[i]] += α * y[i] * cache[i]
	end
end


text:
module WaveToyDagger

using Dagger
using StaticArrays

################################################################################

struct Just{T}
    value::T
end
just(x::Just) = x.value
Base.wait(::Just) = nothing

struct Future{T}
    thunk::Union{Thunk,Dagger.EagerThunk,Dagger.Chunk{T},Just{T}}
    # Enable this to serialize all calculations
    # Future{T}(x) where {T} = new{T}(Dagger.tochunk(get_result(T, x)))
end
Future{T}(x::T) where {T} = Future{T}(Dagger.tochunk(x))
Base.eltype(::Type{Future{T}}) where {T} = T

get_result(::Type{T}, thunk::Thunk) where {T} = collect(thunk)::T
get_result(::Type{T}, thunk::Dagger.EagerThunk) where {T} = fetch(thunk)::T
get_result(::Type{T}, thunk::Dagger.Chunk{T}) where {T} = collect(thunk)::T
get_result(::Type{T}, thunk::Just{T}) where {T} = just(thunk)::T
Base.fetch(f::Future{T}) where {T} = get_result(T, f.thunk)
Base.wait(f::Future) = wait(f.thunk)

# make_ready_future(value::T) where {T} = Future{T}(Dagger.tochunk(value))
# macro future(T, exprs...)
#     isempty(exprs) && error("Syntax: @future <type> {<option>} <expression>")
#     opts = exprs[1:(end - 1)]
#     expr = exprs[end]
#     quote
#         task() = $(esc(expr))
#         Future{$(esc(T))}(Dagger.spawn(task; $((esc(opt) for opt in opts)...)))
#     end
# end

################################################################################

export linterp
function linterp(x1, y1, x2, y2, x)
    R = typeof(zero(y1) + zero(y2))
    return R(x2 - x) ./ R(x2 - x1) .* y1 + R(x - x1) ./ R(x2 - x1) .* y2
end

const D = 3

IND(i::SVector) = CartesianIndex(Tuple(i))
VEC(i::CartesianIndex) = SVector(Tuple(i))

const vones = SVector(ntuple(d -> 1, D))

const npoints = SVector(ntuple(d -> 256, D))
const ngrids = SVector(ntuple(d -> 2, D))
const nghosts = SVector(ntuple(d -> 1, D))
@assert all(npoints .% ngrids .== 0)

const gsh = npoints + 2 * nghosts
const lsh = npoints .÷ ngrids + 2 * nghosts
lbnd(gpos::SVector{D,Int}) = (gpos .- 1) .* (lsh - 2 * nghosts) .+ 1

# TODO: Introduce `local2global`, `global2local`

const extr = IND(vones):IND(lsh)
const intr = IND(vones + nghosts):IND(lsh - nghosts)
function bnds(dirs::SVector{D,Int})
    @assert all(-1 ≤ dir ≤ 1 for dir in dirs)
    intrlo, intrhi = first(intr), last(intr)
    lbnd = CartesianIndex(ntuple(d -> dirs[d] < 0 ? 1 : dirs[d] > 0 ? intrhi[d] + 1 : intrlo[d], D))
    ubnd = CartesianIndex(ntuple(d -> dirs[d] < 0 ? intrlo[d] - 1 : dirs[d] > 0 ? lsh[d] : intrhi[d], D))
    return lbnd:ubnd
end
function srcs(dirs::SVector{D,Int})
    intrlo, intrhi = first(intr), last(intr)
    lbnd = CartesianIndex(ntuple(d -> dirs[d] < 0 ? intrhi[d] - nghosts[d] + 1 : dirs[d] > 0 ? intrlo[d] : intrlo[d], D))
    ubnd = CartesianIndex(ntuple(d -> dirs[d] < 0 ? intrhi[d] : dirs[d] > 0 ? intrlo[d] + nghosts[d] - 1 : intrhi[d], D))
    return lbnd:ubnd
end

const xmin = SVector{D}(-1.0 for d in 1:D)
const xmax = SVector{D}(+1.0 for d in 1:D)
const dx = (xmax - xmin) ./ npoints
xcoord(ipos::SVector{D,Int}) = linterp(1 .+ nghosts, xmin, gsh - nghosts .+ 1, xmax, ipos)

################################################################################

mutable struct Grid{D,T}
    array::Array{T,D}
    gpos::SVector{D,Int}
end
Base.ndims(::Grid{D}) where {D} = D
Base.eltype(::Type{Grid{D,T}}) where {D,T} = T
Base.size(grid::Grid) = size(grid.array)
Base.similar(grid::Grid{D}, ::Type{T}=eltype(grid)) where {D,T} = Grid{D,T}(similar(grid.array), grid.gpos)

mutable struct Domain{D,S,T}
    grids::Array{Future{Grid{D,T}},D}
    time::S
end
Base.ndims(::Domain{D}) where {D} = D
Base.eltype(::Type{Domain{D,S,T}}) where {D,S,T} = T
function Base.similar(domain::Domain{D,S}, ::Type{T}=eltype(domain)) where {D,S,T}
    ## Dagger.@spawn
    return Domain{D,S,T}([Future{Grid{D,T}}(Just{Grid{D,T}}(similar(fetch(grid)))) for grid in domain.grids], domain.time)
end

Base.wait(dom::Domain) = wait.(dom.grids)

@generated function Base.map!(f, r::Grid{D}, x::Grid{D}, ys::Grid{D}...) where {D}
    quote
        $((:(@assert x.gpos == ys[$i].gpos) for i in 1:length(ys))...)
        map!(f, r.array, x.array, $([:(ys[$i].array) for i in 1:length(ys)]...))
        r.gpos = x.gpos
        return r
    end
end
@generated function Base.map(f, x::Grid{D}, ys::Grid{D}...) where {D}
    quote
        $((:(@assert x.gpos == ys[$i].gpos) for i in 1:length(ys))...)
        array = map(f, x.array, $([:(ys[$i].array) for i in 1:length(ys)]...))
        T = eltype(array)
        return Grid{D,T}(array, x.gpos)
    end
end

@inline function Base.map!(f, r::Domain{D,S}, x::Domain{D,S}, ys::Domain{D,S}...) where {D,S}
    f′(xs′...) = map!(f, (fetch(x′) for x′ in xs′)...)
    foreach(f′, r.grids, x.grids, (y.grids for y in ys)...)
    r.time = f(x.time, (y.time for y in ys)...)
    return r
end
@inline function Base.map(f, x::Domain{D,S}, ys::Domain{D,S}...) where {D,S}
    function f′(xs′...)
        grid = map(f, (fetch(x′) for x′ in xs′)...)
        T = eltype(grid)
        grid::Grid{D,T}
        return Future{Grid{D,T}}(Just{Grid{D,T}}(grid))
    end

    grids = map(f′, x.grids, (y.grids for y in ys)...)
    time = f(x.time, (y.time for y in ys)...)
    T = eltype(eltype(eltype(grids)))
    grids::Array{Future{Grid{D,T}},D}
    return Domain{D,S,T}(grids, time)
end

Base.:+(x::Grid{D}, ys::Grid{D}...) where {D} = map(+, x, ys...)
Base.:-(x::Grid{D}, ys::Grid{D}...) where {D} = map(-, x, ys...)
Base.:*(a::Number, x::Grid) = map(b -> a * b, x)

Base.:+(x::Domain{D,S}, ys::Domain{D,S}...) where {D,S} = map(+, x, ys...)
Base.:-(x::Domain{D,S}, ys::Domain{D,S}...) where {D,S} = map(-, x, ys...)
Base.:*(a::Number, x::Domain) = map(b -> a * b, x)

################################################################################

function solution(t::S, x::SVector{D,S}) where {D,S}
    @assert D == 3
    k = SVector{D,S}(1, 1, 1)
    ω = sqrt(sum(k .^ 2))
    u = cospi(ω * t) * sinpi(k[1] * x[1]) * sinpi(k[2] * x[2]) * sinpi(k[3] * x[3])
    ρ = -S(π) * ω * sinpi(ω * t) * sinpi(k[1] * x[1]) * sinpi(k[2] * x[2]) * sinpi(k[3] * x[3])
    v = SVector(S(π) * k[1] * cospi(ω * t) * cospi(k[1] * x[1]) * sinpi(k[2] * x[2]) * sinpi(k[3] * x[3]),
                S(π) * k[2] * cospi(ω * t) * sinpi(k[1] * x[1]) * cospi(k[2] * x[2]) * sinpi(k[3] * x[3]),
                S(π) * k[3] * cospi(ω * t) * sinpi(k[1] * x[1]) * sinpi(k[2] * x[2]) * cospi(k[3] * x[3]))
    return SVector(u, ρ, v...)
end

################################################################################

function initialize_grid(::Type{Grid{D,T}}, gpos::SVector{D,Int}, time::S) where {D,S,T}
    @assert S ≡ eltype(T)
    array = Array{T,D}(undef, Tuple(lsh))
    for i in extr
        ipos = lbnd(gpos) + VEC(i) .- 1
        x = SVector{D,S}(xcoord(ipos))
        array[i] = solution(time, x)
    end
    return Grid{D,T}(array, gpos)
end

function initialize_domain(::Type{Domain{D,S,T}}, time::S) where {D,S,T}
    grids = Array{Future{Grid{D,T}},D}(undef, Tuple(ngrids))
    for gi in IND(vones):IND(ngrids)
        g = VEC(gi)
        ## grids[IND(g)] = Future{Grid{D,T}}(Dagger.@spawn initialize_grid(Grid{D,T}, g))
        grids[IND(g)] = Future{Grid{D,T}}(Just(initialize_grid(Grid{D,T}, g, time)))
    end
    return Domain{D,S,T}(grids, time)
end

function rhs_intr(res::Grid{D,T}, grid::Grid{D,T}) where {D,T}
    array = res.array
    if size(array) ≠ size(grid.array)
        array = similar(grid.array)
    end
    @assert D == 3
    di = SVector(CartesianIndex(1, 0, 0), CartesianIndex(0, 1, 0), CartesianIndex(0, 0, 1))
    for i in intr
        u, ρ, vx, vy, vz = grid.array[i]
        # TODO: Move pointwise rhs into its own function
        u̇ = ρ
        ρ̇ = (grid.array[i + di[1]][2 + 1] - grid.array[i - di[1]][2 + 1]) / 2dx[1] +
             (grid.array[i + di[2]][2 + 2] - grid.array[i - di[2]][2 + 2]) / 2dx[2] +
             (grid.array[i + di[3]][2 + 3] - grid.array[i - di[3]][2 + 3]) / 2dx[3]
        vẋ = (grid.array[i + di[1]][2] - grid.array[i - di[1]][2]) / 2dx[1]
        vẏ = (grid.array[i + di[2]][2] - grid.array[i - di[2]][2]) / 2dx[2]
        vż = (grid.array[i + di[3]][2] - grid.array[i - di[3]][2]) / 2dx[3]
        array[i] = SVector(u̇, ρ̇, vẋ, vẏ, vż)
    end
    return Grid{D,T}(array, grid.gpos)
end

function rhs_intr(res::Domain{D,S,T}, domain::Domain{D,S,T}, cache::Ref{Any}) where {D,S,T}
    grids = res.grids
    if size(grids) ≠ size(res.grids)
        grids = similar(res.grids)
    end
    if !(typeof(cache[]) isa Array{Future{Grid{D,T}},D} && size(cache[]) == size(domain.grids))
        cache[] = similar(domain).grids
    end
    arrays = cache[]
    for gi in IND(vones):IND(ngrids)
        g = VEC(gi)
        res = arrays[IND(g)]
        ## grids[IND(g)] = Future{Grid{D,T}}(Dagger.@spawn rhs_intr(domain.grids[IND(g)].thunk))
        grids[IND(g)] = Future{Grid{D,T}}(Just(rhs_intr(fetch(res), fetch(domain.grids[IND(g)]))))
    end
    return Domain{D,S,T}(grids, S(1))
end

function calc_error(domain::Domain{D,S,T}) where {D,S,T}
    domain₀ = initialize_domain(Domain{D,S,T}, domain.time)
    return domain - domain₀
end

get_ghosts(grid::Grid{D,T}, dirs::SVector{D,Int}) where {D,T} = grid.array[srcs(dirs)]

function set_ghosts!(grid::Grid{D,T}, ghosts::Array{Future{Array{T,D}},D}) where {D,T}
    @assert D == 3
    for dirsi in CartesianIndex(-1, -1, -1):CartesianIndex(1, 1, 1)
        dirs = VEC(dirsi)
        if !iszero(dirs)
            gh = fetch(ghosts[IND(dirs .+ 2)])
            # Check region sizes
            @assert size(gh) == size(bnds(dirs))
            @assert size(srcs(dirs)) == size(bnds(dirs))
            # Set ghost zvones
            grid.array[bnds(dirs)] = gh
        end
    end
    return grid
end

function exchange_ghosts(domain::Domain{D,S,T}) where {D,S,T}
    grids = Array{Future{Grid{D,T}},D}(undef, Tuple(ngrids))
    for gi in IND(vones):IND(ngrids)
        g = VEC(gi)
        ghosts = Array{Future{Array{T,D}},D}(undef, ntuple(d -> 3, D))
        @assert D == 3
        for dirsi in CartesianIndex(-1, -1, -1):CartesianIndex(1, 1, 1)
            dirs = VEC(dirsi)
            if !iszero(dirs)
                ## ghosts[IND(dirs .+ 2)] = Future{Array{T,D}}(Dagger.@spawn get_ghosts(domain.grids[IND(mod1.(g + dirs, ngrids))].thunk,
                ##                                                                      dirs))
                ghosts[IND(dirs .+ 2)] = Future{Array{T,D}}(Just(get_ghosts(fetch(domain.grids[IND(mod1.(g + dirs, ngrids))]),
                                                                            dirs)))
            else
                ghosts[IND(dirs .+ 2)] = Future{Array{T,D}}(Just(zeros(T, ntuple(d -> 0, D))))
            end
        end
        ## grids[IND(g)] = Future{Grid{D,T}}(Dagger.@spawn set_ghosts!(domain.grids[IND(g)].thunk, ghosts))
        grids[IND(g)] = Future{Grid{D,T}}(Just(set_ghosts!(fetch(domain.grids[IND(g)]), ghosts)))
    end
    return Domain{D,S,T}(grids, domain.time)
end

################################################################################

function maxabs(grid::Grid{D,T}) where {D,T}
    abs′(x) = abs.(x)
    max′(x, y) = max.(x, y)
    return mapreduce(abs′, max′, grid.array; init=zero(T))::T
end

function maxabs(domain::Domain{D,S,T}) where {D,S,T}
    thunks = Future{T}[]
    for gi in IND(vones):IND(ngrids)
        g = VEC(gi)
        ## push!(thunks, Future{T}(Dagger.@spawn maxabs(domain.grids[IND(g)].thunk)))
        push!(thunks, Future{T}(Just(maxabs(fetch(domain.grids[IND(g)])))))
    end
    res = maximum(fetch(th) for th in thunks)
    return res
end

################################################################################

function rhs(res::Domain{D,S,T}, dom::Domain{D,S,T}, cache::Ref{Any}) where {D,S,T}
    res = rhs_intr(res, dom, cache)
    res = exchange_ghosts(res)
    return res
end

# function rk2(f, u0, h)
#     k1 = f(u0)
#     u1 = u0 + (h / 2) * k1
#     k2 = f(u1)
#     u2 = u0 + h * k2
#     return u2
# end

# function rk2(f, u0, h)
#     k1 = f(u0)
#     u1 = map((u0, k1) -> u0 + (h / 2) * k1, u0, k1)
#     k2 = f(u1)
#     u2 = map((u0, k2) -> u0 + h * k2, u0, k2)
#     return u2
# end

struct rk2Cache{V}
    k1::V
    u1::V
    k2::V
    fcache::Ref{Any}
end
function rk2(u2::V, f, u0::V, h, cache::Ref{Any}) where {V}
    if !(cache[] isa rk2Cache{V})
        cache[] = rk2Cache{V}(similar(u0), similar(u0), similar(u0), Ref{Any}(nothing))
    end
    k1, u1, k2, fcache = cache[].k1, cache[].u1, cache[].k2, cache[].fcache

    k1 = f(k1, u0, fcache)
    map!((u0, k1) -> u0 + (h / 2) * k1, u1, u0, k1)
    k2 = f(k2, u1, fcache)
    map!((u0, k2) -> u0 + h * k2, u2, u0, k2)

    return u2
end

function main()
    # ctx = Context()
    # log = Dagger.LocalEventLog()
    # ctx.log_sink = log

    S = Float64
    T = SVector{D + 2,S}
    h = minimum(dx) / 2
    niters = 10

    println("Initial conditions...")
    t = S(0)
    dom = initialize_domain(Domain{D,S,T}, t)
    dom = exchange_ghosts(dom)

    rhscache = Ref{Any}(nothing)
    dom′ = similar(dom)

    for iter in 1:niters
        println("Iteration $iter...")
        dom′ = rk2(dom′, rhs, dom, h, rhscache)

        # Limit parallelism
        wait(dom)
        dom, dom′ = dom′, dom
    end

    # TODO: Use `yield` (?) to show a progress bar or something

    err = calc_error(dom)
    ma = maxabs(err)
    println("maxabs[error]=$ma")

    println("Done.")

    # logs = Dagger.get_logs!(log)
    # Dagger.show_plan(stdout, logs)

    return nothing
end

end


text:
#const survival = 0.5
const survival = 0.9

function evolve(population, board_size, max_moves, nplays, nrounds, path)
    known = Set{UInt64}()
    exists(path) && rm(path)
    dump(known, path, population)
    for i in 1:nrounds
        stats = reset_stats()
        for j in 1:nplays
            a, b = pick_competitors(length(population))
            black, white  = players(a, b, population)
            result = play(population[black], population[white],
                          board_size, max_moves, null_display, null_display)
            display_result(i, nrounds, j, nplays, population, a, b, result)
            apply_result!(population, a, b, result)
            stats = update_stats(stats, result)
        end
        for j in 1:min(3, length(population))
            println("$(name(population[j], j))")
            dump_expression(population[j])
        end
        println("")
        if i != nrounds
            population = update_population(population, stats, board_size)
            dump(known, path, population)
        end
        println("")
    end
    population
end

function pick_competitors(n)
    # surprise() below assumes these are ordered
    b = rand(2:n)
    # biased to top, with one always in surviving fraction
    rand(1:min(b-1, floor(Int, survival * n))), b
end

const name_length = 16
name(data) = sha1(data)[1:name_length]
name(data, rank) = "$(name(data)):$(rank)"

function isless(a::Vector{UInt8}, b::Vector{UInt8})
    for (aa, bb) in zip(a, b)
        if aa < bb
            return true
        elseif aa > bb
            return false
        end
    end
    return false
end

# returns black, white
players(a, b, popn) = popn[a] < popn[b] ? (a, b) : (b, a)

function surprise(a, b, population, result)
    black, white = players(a, b, population)
    if abs(result.stats.score) < 0.1
        return false
    else
        return (black < white) != (result.stats.score > 0)
    end
end

function display_result(i, n, j, m, population, a, b, result)
    # a, b are indices, where a < b
    # a lost to b if surprise(a, b, population, result)
    # we display black on the left, white on the right
    # we indicate the winner with <, > or = which is preceded by !
    # if the results is a surprise
    black, white = players(a, b, population)
    @printf("%3d/%d %3d/%d %s", i, n, j, m, surprise(a, b, population, result) ? "!" : " ")
    @printf(" %s:%-3d ", name(population[black]), black)
    if result.stats.score > 0
        print(">")
    elseif result.stats.score == 0
        print("=")
    else
        print("<")
    end
    @printf(" %s:%-3d ", name(population[white]), white)
    @printf("%6.1f sc  %3d mv  %2d sp  %2d st\n", result.stats.score, result.stats.moves, result.stats.owned, result.stats.stones)
end

function apply_result!(population, a, b, result)
    if surprise(a, b, population, result)
        # this loses good individuals too quickly
#        population[a], population[b] = population[b], population[a]
        a, b = sort([a, b])
        tmp = population[b]
        population[a+1:b] = population[a:b-1]
        population[a] = tmp
    end
end

# track fraction of owned space as an indicator of quality

reset_stats() = 0, 0.0

function update_stats(stats, result)
    n, sum = stats
    # ignore games with one stone since i has an artificially high owned
    if result.stats.stones > 1
        n+1, sum+result.stats.owned
    else
        n, sum
    end
end

function distinct_indices(n; start=1)
    i = rand(start:n)
    j = i
    while i == j
        j = rand(start:n)
    end
    i, j
end

function biased_pair(population)
    a = rand(2:length(population))
    b = rand(1:(a-1))
    population[a], population[b]
end

function distinct_pair(population)
    a, b = distinct_indices(length(population))
    population[a], population[b]
end

biased_single(population) = biased_pair(population)[2]

random_single(population) = rand(population)

function rotate(a::Vector)
    n = rand(lheader+2:length(a))
    c = vcat(a[1:lheader], a[n:end], a[lheader+1:n-1])
    print("rotate $(n-lheader+1) $(name(a)) => $(name(c)) + ")
    c
end

function rotate_chunks(a::Vector)
    k = rand(0:cld(length(a) - lheader, chunk)-1)
    n = lheader + chunk * k + 1
    c = vcat(a[1:lheader], a[n:end], a[lheader+1:n-1])
    print("rotate chunk $k $(name(a)) => $(name(c)) + ")
    c
end

function merge(a::Vector, b::Vector)
    a, b = shuffle(Vector[a, b])
    c, i, j, attempts = a, 0, 0, 0
    while c in (a, b) && attempts < 10
        i, j = sort([distinct_indices(length(a); start=lheader)...])
        c = vcat(a[1:i], b[i+1:j], a[j+1:end])
        attempts += 1
    end
    println("merge $(j-i) $(name(a)), $(name(b)) => $(name(c))")
    c
end

function merge_with_rotate(a::Vector, b::Vector)
    a, b = shuffle(Vector[a, b])
    merge(a, rotate(b))
end

function merge_with_rotate_chunks(a::Vector, b::Vector)
    a, b = shuffle(Vector[a, b])
    merge(a, rotate_chunks(b))
end

function random_bytes(fraction)
    function change(a::Vector{UInt8})
        b, count = copy(a), 0
        while b == a
            for i in lheader+1:length(b)
                if rand(Float64) < fraction
                    b[i] = rand(UInt8)
                    count += 1
                end
            end
        end
        println("$(count) random bytes $(name(a)) => $(name(b))")
        b
    end
end

function random_bits(fraction)
    function change(a::Vector{UInt8})
        b, count = copy(a), 0
        while b == a
            for i in lheader+1:length(b)
                if rand(Float64) < fraction
                    b[i] = b[i] $ (0x01 << rand(0:7))
                    count += 1
                end
            end
        end
        println("$(count) random bits $(name(a)) => $(name(b))")
        b
    end
end

function random_deltas(fraction)
    function change(a::Vector{UInt8})
        b, count = copy(a), 0
        while b == a
            for i in lheader+1:length(b)
                if rand(Float64) < fraction
                    b[i] = UInt8(mod(b[i] + rand([1,-1]), 0x100))
                    count += 1
                end
            end
        end
        println("$(count) random deltas $(name(a)) => $(name(b))")
        b
    end
end

function weighted_rand(ops)
    total = sum([op[1] for op in ops])
    limit = rand(Float64) * total
    acc, i = ops[1][1], 1
    while acc < limit
        acc += ops[i][1]
        i += 1
    end
    ops[i][2]
end

function build_ops(length, temp)
    ops = Tuple{Number, Function}[
              (10, x -> random_bits(temp/100)(random_single(x))),
              (10, x -> random_deltas(temp/100)(random_single(x))),
              (10, x -> merge(biased_pair(x)...)),
              (1, x -> merge_with_rotate_chunks(biased_pair(x)...)),
              (1, x -> random_bytes(temp/100)(random_single(x)))]
    if temp > 0.25
        push!(ops, (10, x -> random_bytes(temp/100)(random_single(x))))
    end
    if temp > 0.5
        push!(ops, (10, x -> merge_with_rotate(biased_pair(x)...)))
    end
    ops
end

function update_population(prev_pop, stats, board_size)
    t, sum = stats
    # temperature varies roughly from 1 (hottest) to 0 (coldest)
    temp = t == 0 ? 1 : max(0, min(1, 1 - (2 * sum / t) / board_size^2))
    @printf("temp %4.3f\n", temp)
    ops = build_ops(length(prev_pop[1]), temp)
    n = length(prev_pop)
    m = Int(survival * n)
    next_pop = Vector{UInt8}[]
    while length(next_pop) < m && length(prev_pop) > 1
        push!(next_pop, prev_pop[1])
        prev_pop = prev_pop[2:end]
    end
    while length(next_pop) < n
        push!(next_pop, weighted_rand(ops)(next_pop))
    end
    next_pop
end

function dump(known, path, population)
    open(path, "a") do io
        for data in population
            n = name(data)
            tag = reinterpret(UInt64, hex2bytes(n))[1]
            if !(tag in known)
                print(io, string(n, ":"))
                print(io, bytes2hex(data))
                print(io, "\n")
                push!(known, tag)
            end
        end
    end
end

function undump(path)
    named = Dict{AbstractString, Vector{UInt8}}()
    ordered = Vector{UInt8}[]
    open(path, "r") do io
        for line in eachline(io)
            name, data = split(line, ":")
            net = hex2bytes(data[1:end-1])  # drop newline
            named[name] = net
            push!(ordered, net)
        end
    end
    (named, ordered)
end

function evolve2(population, board_size, max_moves, nplays, nrounds, path)
    known = Set{UInt64}()
    exists(path) && rm(path)
    dump(known, path, population)
    for i in 1:nrounds
        stats = reset_stats()
        for j in 1:nplays
            a, b = pick_competitors(length(population))
            black, white  = players(a, b, population)
            result = play(population[black], population[white],
                          board_size, max_moves, null_display, null_display)
            display_result(i, nrounds, j, nplays, population, a, b, result)
            apply_result!(population, a, b, result)
            stats = update_stats(stats, result)
        end
        for j in 1:min(3, length(population))
            println("$(name(population[j], j))")
            dump_expression(population[j])
        end
        println("")
        if i != nrounds
            population = update_population2(known, i, max_moves, nplays,
                                            nrounds, path,
                                            population, stats, board_size)
        end
        println("")
    end
    population
end

function update_population2(known, i, max_moves, nplays, nrounds, path,
                            prev_pop, stats, board_size)
    t, sum = stats
    # temperature varies roughly from 1 (hottest) to 0 (coldest)
    temp = t == 0 ? 1 : max(0, min(1, 1 - (2 * sum / t) / board_size^2))
    @printf("temp %4.3f\n", temp)
    ops = build_ops(length(prev_pop[1]), temp)
    n = length(prev_pop)
    m = Int(survival * n)
    next_pop = Vector{UInt8}[]
    while length(next_pop) < m && length(prev_pop) > 1
        push!(next_pop, prev_pop[1])
        prev_pop = prev_pop[2:end]
    end
    m, j = length(next_pop), nplays
    while m < n
        push!(next_pop, weighted_rand(ops)(next_pop))
        a, b, j = length(next_pop), rand(1:m), j+1
        black, white  = players(a, b, next_pop)
        result = play(next_pop[black], next_pop[white],
                      board_size, max_moves, null_display, null_display)
        display_result(i, nrounds, j, nplays, next_pop, a, b, result)
        if surprise(a, b, next_pop, result)
            apply_result!(next_pop, a, b, result)
            m += 1
        end
        dump(known, path, next_pop)
    end
    next_pop[1:n]
end


text:
@testset "Ripple" begin
    img, p, q, Z = generate_surface_data(Ripple(5), 0.5, [0,0,1])
    @test maximum(img) ≈ 0.5
    @test size(img)[1] ≈ 151
    @test img[10,10] ≈ 0.37078406641129436
    @test size(p)[1] ≈ 151
    @test size(q)[2] ≈ 151
    @test p[32,54] ≈ -0.8993519378588697
    @test p[79,100] ≈ 4.922700205722557
    @test q[50,77] ≈ 4.404977131317488
    @test q[42,42] ≈ 0.6742231223764109
    @test Z[10,10] ≈ 24.587377136574045
    @test Z[75,75] ≈ 24.006648914263153
end


text:
const BACKUP_STDIN_TERMIOS = Ref{TERMIOS.termios}()
const BACKUP_STDOUT_TERMIOS = Ref{TERMIOS.termios}()

const WIDTH = Ref{Int}()
const HEIGHT = Ref{Int}()
const MODE = Ref{Symbol}(:default)

function terminal_size(io)
    ws = IOCTL.ioctl(io, IOCTL.TIOCGWINSZ)
    # width, height
    return (Int(ws.ws_col), Int(ws.ws_row))
end

function terminal_size()
    ds = displaysize(stdout)
    return (last(ds), first(ds))
end

terminal_size(io, coord::Int) = terminal_size(io)[coord]
terminal_size(coord::Int) = terminal_size(stdout, coord)

const END = 2 + 1
const CONTROL_SEQUENCE_TIMEOUT = 100e-3

const INDICES = String["$i" for i in 1:1000]

struct CSI
    parameter::String
    immediate::String
    final::String
    _str::String
    CSI(parameter, immediate, final) =
        new(parameter, immediate, final, "$(Terminals.CSI)$(parameter)$(immediate)$(final)")
end

CSI(; parameter = "", immediate = "", final = "") = CSI(parameter, immediate, final)

Base.print(io::IO, csi::CSI) = print(io, csi._str)

const CLEARSCREEN = CSI(final = "2J")
const CLEARBELOWCURSOR = CSI(final = "J")
const CLEARABOVECURSOR = CSI(final = "1J")

const CLEARLINE = CSI(final = "2K")
const CLEARAFTERCURSOR = CSI(final = "K")
const CLEARBEFORECURSOR = CSI(final = "1K")

const HOMEPOSITION = CSI(final = "H")

const BELL = '\x07'

const HIDECURSOR = CSI(final = "?25l")
const SHOWCURSOR = CSI(final = "?25h")

const SAVECURSOR = CSI(final = "s")
const RESTORECURSOR = CSI(final = "u")

const LOCATECURSOR = CSI(final = "6n")

const RESET = "\ec"

const TUIMODE = CSI(final = "?1049h")
const DEFAULTMODE = CSI(final = "?1049l")

const CURSORBLINKINGBLOCK = CSI(final = "1 q")
const CURSORSTEADYBLOCK = CSI(final = "2 q")
const CURSORBLINKINGUNDERLINE = CSI(final = "3 q")
const CURSORSTEADYUNDERLINE = CSI(final = "4 q")
const CURSORBLINKINGIBEAM = CSI(final = "5 q")
const CURSORSTEADYIBEAM = CSI(final = "6 q")

clear_screen() = print(stdout, CLEARSCREEN)
clear_screen_from_cursor_up() = print(stdout, CLEARBELOWCURSOR)
clear_screen_from_cursor_down() = print(stdout, CLEARABOVECURSOR)

clear_line() = print(stdout, CLEARLINE)
clear_line_from_cursor_right() = print(stdout, CLEARAFTERCURSOR)
clear_line_from_cursor_left() = print(stdout, CLEARBEFORECURSOR)

hide_cursor() = print(stdout, HIDECURSOR)
show_cursor() = print(stdout, SHOWCURSOR)

save_cursor() = print(stdout, SAVECURSOR)
restore_cursor() = print(stdout, RESTORECURSOR)

move_cursor(row, col) = print(stdout, Terminals.CSI, INDICES[row], ';', INDICES[col], 'H')
move_cursor_up(row = 1) = print(stdout, Terminals.CSI, INDICES[row], 'A')
move_cursor_down(row = 1) = print(stdout, Terminals.CSI, INDICES[row], 'B')
move_cursor_right(col = 1) = print(stdout, Terminals.CSI, INDICES[col], 'C')
move_cursor_left(col = 1) = print(stdout, Terminals.CSI, INDICES[col], 'D')
move_cursor_home() = print(stdout, HOMEPOSITION)

change_cursor_to_blinking_block() = print(stdout, CURSORBLINKINGBLOCK)
change_cursor_to_steady_block() = print(stdout, CURSORSTEADYBLOCK)
change_cursor_to_blinking_underline() = print(stdout, CURSORBLINKINGUNDERLINE)
change_cursor_to_steady_underline() = print(stdout, CURSORSTEADYUNDERLINE)
change_cursor_to_blinking_ibeam() = print(stdout, CURSORBLINKINGIBEAM)
change_cursor_to_steady_ibeam() = print(stdout, CURSORSTEADYIBEAM)

reset() = print(stdout, RESET)

tui_mode() = print(stdout, TUIMODE)
default_mode() = print(stdout, DEFAULTMODE)

function with_raw_mode(f::Function)
    # if already in in raw mode don't enter raw mode again.
    # entering raw mode
    entered_raw_mode = false
    if MODE[] != :raw
        enable_raw_mode()
        entered_raw_mode = true
    end

    f()

    entered_raw_mode ? disable_raw_mode() : nothing
end

function detect_color(color::Integer)
    with_raw_mode() do
        print(stdout, "\x1B]4;$(color);?$(BELL)")
        Base.flush(stdout)
        total_read = 0
        now = time()
        c = 0x00
        bell = UInt8(BELL)
        channel = Channel(1)
        t = @async begin
            while c != bell
                c = read(stdin, 1)[]
                total_read += 1
            end
        end
        while c != bell && time() - now < CONTROL_SEQUENCE_TIMEOUT
            sleep(1e-3)
        end
        if !istaskdone(t)
            t.exception = InterruptException()
        end
    end
    return total_read > 0
end

function locate_cursor()
    location = UInt8[]
    with_raw_mode() do
        print(stdout, LOCATECURSOR)
        Base.flush(stdout)
        total_read = 0
        now = time()
        c = 0x00
        bell = UInt8(BELL)
        channel = Channel(1)
        t = @async begin
            while c != UInt8('R')
                c = read(stdin, 1)[]
                push!(location, c)
                total_read += 1
            end
        end
        while c != UInt8('R') && time() - now < CONTROL_SEQUENCE_TIMEOUT
            sleep(1e-3)
        end
        if !istaskdone(t)
            t.exception = InterruptException()
        end
    end
    row, col = split(String(Char.(location)[3:end-1]), ";")
    return parse(Int, row), parse(Int, col)
end

function enable_raw_mode()
    termios = TERMIOS.termios()
    TERMIOS.tcgetattr(stdin, termios)
     # Disable ctrl-c, disable CR translation, disable stripping 8th bit (unicode), disable parity
    termios.c_iflag &= ~(TERMIOS.BRKINT | TERMIOS.ICRNL | TERMIOS.INPCK | TERMIOS.ISTRIP | TERMIOS.IXON)
    # Disable output processing
    termios.c_oflag &= ~(TERMIOS.OPOST)
    # Disable parity
    termios.c_cflag &= ~(TERMIOS.CSIZE | TERMIOS.PARENB)
    # Set character size to 8 bits (unicode)
    termios.c_cflag |= (TERMIOS.CS8)
     # Disable echo, disable canonical mode (line mode), disable input processing, disable signals
    termios.c_lflag &= ~(TERMIOS.ECHO | TERMIOS.ICANON | TERMIOS.IEXTEN | TERMIOS.ISIG)
    termios.c_cc[TERMIOS.VMIN] = 0
    termios.c_cc[TERMIOS.VTIME] = 1
    TERMIOS.tcsetattr(stdin, TERMIOS.TCSANOW, termios)
end

function disable_raw_mode()
    TERMIOS.tcsetattr(stdin, TERMIOS.TCSANOW, BACKUP_STDIN_TERMIOS[])
    MODE[] = :default
end

function backup_termios()
    BACKUP_STDIN_TERMIOS[] = TERMIOS.termios()
    BACKUP_STDOUT_TERMIOS[] = TERMIOS.termios()
    TERMIOS.tcgetattr(stdin, BACKUP_STDIN_TERMIOS[])
    TERMIOS.tcgetattr(stdout, BACKUP_STDOUT_TERMIOS[])
end

function initialize()
    backup_termios()
    tui_mode()
    hide_cursor()
    enable_raw_mode()
    clear_screen()
    move_cursor_home()
end

function cleanup()
    move_cursor_home()
    clear_screen()
    disable_raw_mode()
    show_cursor()
    default_mode()
end


text:
# This file includes utiliti functions for Systems module

macro siminfo(msg...)
    quote
        @info "$(now()) $($msg...)"
    end
end

#=
    @def begin name
        code
    end

Copy paste macro
=#
macro def(name, code)
    quote
        macro $(esc(name))()
            esc($(Meta.quot(code)))
        end
    end
end

hasargs(func, n) = n + 1 in [method.nargs for method in methods(func)]

function unwrap(container, etype; depth=10)
    for i in 1 : depth
        container = vcat(container...)
        eltype(container) == etype && break
    end
    container
end

launchport(iport) = @async while true
    all(take!(iport) .=== NaN) && break
end

# Equips `comp` to make it launchable. Equipment is done by constructing and connecting signalling pins (i.e. `trigger`
# and `handshake`), input and output ports (if necessary)
function equip(comp, kickoff::Bool=true)
    oport = typeof(comp) <: AbstractSource ?
        nothing : (typeof(comp.input) === nothing  ? nothing : Outport(length(comp.input)))
    iport = typeof(comp) <: AbstractSink ?
        nothing : (typeof(comp.output) === nothing ? nothing : Inport(length(comp.output)))
    trg = Outpin()
    hnd = Inpin{Bool}()
    oport === nothing || connect!(oport, comp.input)
    iport === nothing || connect!(comp.output, iport)
    connect!(trg, comp.trigger)
    connect!(comp.handshake, hnd)
    if kickoff
        comptask, outputtask = launch(comp), launchport(iport)
    else
        comptask, outputtask = nothing, nothing
    end
    oport, iport, trg, hnd, comptask, outputtask
end


text:
module RocketGenerateObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "GenerateObservable" begin

    println("Testing: generate")

    @testset begin
        source = generate(1, x -> x < 2, x -> x + 1)
        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("GenerateObservable", printed)
        @test occursin(string(eltype(source)), printed)
    end

    run_testset([
        (
            source      = generate(1, x -> x < 3, x -> x + 1),
            values      = @ts([ 1, 2, c ]),
            source_type = Int,
        ),
        (
            source      = generate(1.0, x -> x < 0.5, x -> x + 1.0),
            values      = @ts(c),
            source_type = Float64,
        ),
        (
            source      = generate(1.0, x -> false, x -> x),
            values      = @ts(c),
            source_type = Float64,
        ),
        (
            source      = generate(1, x -> x < 2, x -> x + 1, scheduler = AsyncScheduler(0)),
            values      = @ts([ 1 ] ~ c),
            source_type = Int
        )
    ])

end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ParameterDefinitions <: SwaggerModel

    function ParameterDefinitions(;)
        o = new()
        o
    end
end # type ParameterDefinitions

const _property_map_ParameterDefinitions = Dict{Symbol,Symbol}()
const _property_types_ParameterDefinitions = Dict{Symbol,String}()
Base.propertynames(::Type{ ParameterDefinitions }) = collect(keys(_property_map_ParameterDefinitions))
Swagger.property_type(::Type{ ParameterDefinitions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ParameterDefinitions[name]))}
Swagger.field_name(::Type{ ParameterDefinitions }, property_name::Symbol) =  _property_map_ParameterDefinitions[property_name]

function check_required(o::ParameterDefinitions)
    true
end

function validate_property(::Type{ ParameterDefinitions }, name::Symbol, val)
end


text:
using MLSt
using Test

@testset "MLSt.jl" begin
    # Write your tests here.
end


text:
module MyExample

using ForwardDiff
include("somefile.jl")
end


text:
module GNSS

export readsp3

function readsp3(filename::String)
    # open the file
    f = open(filename)

    data = []
    lineno = 0
    for line in eachline(f)
        lineno += 1

        # check file type
        if lineno == 1
            if line[2] != 'c'
                error("only SP3-c format is supported")
            end
            if line[3] != 'P'
                error("only position-only files are supported")
            end
            nepochs = int(line[33:39])
            nepoch = 0
            data = Array(Real, nepochs, 10, 32)
        end

        # check time system
        if lineno == 13
            if line[4:5] != "G "
                error("only GPS data files are supported")
            end
            if line[10:12] != "GPS"
                error("only GPS time is supported")
            end
        end

        if line[1:2] == "* "
            nepoch += 1
            t = [
                int(line[4:7])     # year
                int(line[9:10])    # month
                int(line[12:13])   # day
                int(line[15:16])   # hour
                int(line[18:19])   # minute
                float(line[21:31]) # second
                ]'
        end

        if line[1] == 'P'
            prn = int(line[3:4])
            p = [
                float(line[5:18])  # x (km)
                float(line[19:32]) # y (km)
                float(line[33:46]) # z (km)
                float(line[47:60]) # clock (μs)
                ]'
            data[nepoch, :, prn] = [t p]
        end
    end

    close(f)
    return data

end

end # module


text:
#Function for calculating two particle matrix elements required for second order perturbation theory evalutation of the blockaed shift
#25/07/2017

function matrixel2p(atom,nni,lli,jji,mji,nn1,ll1,jj1,nn2,ll2,jj2,defect,theta)

  #Predefinitions for efficiency
  sint = sin(theta)
  cost = cos(theta)
  sin2t = sint.^2
  cos2t = cost.^2

  #First calculate angular factors as vectors [minus,zero,plus] for each single particle interaction stae

  angvec1 = angcalc(lli,jji,mji,ll1,jj1)

  angvec2 = angcalc(lli,jji,mji,ll2,jj2)

  #Now the raidla factors

  matelem1i = radiel(atom,nni,lli,jji,nn1,ll1,jj1)

  matelem2i = radiel(atom,nni,lli,jji,nn2,ll2,jj2)

#Calculate terms in two particle matrix element using Wigner-Eckhart theorem (Pritchard, 2012), (Reinhard et al., 2007)
line1 = angvec1[3]*angvec2[3] + angvec1[1]*angvec2[3] + angvec1[2]*angvec2[2]*(1-3*cos2t)
line2 = -(3/2)*sin2t*(angvec1[3]*angvec2[3] + angvec1[3]*angvec2[1] + angvec1[1]*angvec2[3] + angvec1[1]*angvec2[1])
line3 = -(3/sqrt(2))*sint.*cost*(angvec1[3]*angvec2[2] + angvec1[1]*angvec2[2] + angvec1[2]*angvec2[3] + angvec1[2]*angvec2[1])

matrixelement = matelem1i*matelem2i*(line1+line2+line3)

matrixel2P = (matrixelement.^2)./defect

return matrixel2P
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using DashBootstrapComponents

accordion = html_div(
    dbc_accordion(
        [
            dbc_accordionitem("This is the content of the first section", title = "Item 1"),
            dbc_accordionitem(
                "This is the content of the second section",
                title = "Item 2",
            ),
            dbc_accordionitem("This is the content of the third section", title = "Item 3"),
        ],
        always_open = true,
    ),
)