text:
using AstroTime: TDBEpoch, TTEpoch, centuries, j2000, value, julian_twopart
using ReferenceFrameRotations: angleaxis_to_dcm
using StaticArrays: SVector

using ..Bodies: CelestialBody,
    Mercury,
    Venus,
    EarthBarycenter,
    Mars,
    Jupiter,
    Saturn,
    Uranus,
    Neptune,
    Pluto,
    Sun,
    sun,
    grav_param
using ..Constants: astronomical_unit
using ..EarthAttitude: obliquity_of_ecliptic_06
using ..TwoBody: cartesian, transform, mean_anomaly, true_anomaly
using ..Util: sec2rad, normalize_angle

import ..Interfaces: AbstractEphemeris, position!, state!, velocity!

export SimonBretagnon, simon_bretagnon

struct SimonBretagnon <: AbstractEphemeris end
const simon_bretagnon = SimonBretagnon()

const ϵ = 0.40909280422232897

function semi_major(body::CelestialBody, t)
    throw(ArgumentError("Body '$body' is not supported by the Simon-Bretagnon ephemeris."))
end

semi_major(::Mercury, t) = 0.3870983098
eccentricity(::Mercury, t) = @evalpoly(t,
	0.2056317526,
	0.0002040653,
	-28349e-10,
)
inclination(::Mercury, t) = sec2rad(@evalpoly(t,
	7.00498625 * 3600,
	-214.25629,
	0.28977,
))
ascending_node(::Mercury, t) = sec2rad(@evalpoly(t,
	48.33089304 * 3600,
	-4515.21727,
	-31.79892,
))
longitude_of_perihelion(::Mercury, t) = sec2rad(@evalpoly(t,
	77.45611904 * 3600,
	5719.11590,
	-4.83016,
))
mean_longitude(::Mercury, t) = sec2rad(@evalpoly(t,
	252.25090552 * 3600,
	5381016286.88982,
	-1.92789,
))

semi_major(::Venus, t) = 0.7233298200
eccentricity(::Venus, t) = @evalpoly(t,
	0.0067719164,
	-0.0004776521,
	98127e-10,
)
inclination(::Venus, t) = sec2rad(@evalpoly(t,
	3.39466189 * 3600,
	-30.84437,
	-11.67836,
))
ascending_node(::Venus, t) = sec2rad(@evalpoly(t,
	76.67992019 * 3600,
	-10008.48154,
	-51.32614,
))
longitude_of_perihelion(::Venus, t) = sec2rad(@evalpoly(t,
	131.56370300 * 3600,
	175.48640,
	-498.48184,
))
mean_longitude(::Venus, t) = sec2rad(@evalpoly(t,
	181.97980085 * 3600,
	2106641364.33548,
	0.59381,
))

semi_major(::EarthBarycenter, t) = 1.0000010178
eccentricity(::EarthBarycenter, t) = @evalpoly(t,
	0.0167086342,
	-0.0004203654,
	-0.0000126734,
)
inclination(::EarthBarycenter, t) = sec2rad(@evalpoly(t,
	0.0,
	469.97289,
	-3.35053 ,
))
ascending_node(::EarthBarycenter, t) = sec2rad(@evalpoly(t,
	174.87317577 * 3600,
	-8679.27034,
	15.34191,
))
longitude_of_perihelion(::EarthBarycenter, t) = sec2rad(@evalpoly(t,
	102.93734808 * 3600,
	11612.35290,
	53.27577,
))
mean_longitude(::EarthBarycenter, t) = sec2rad(@evalpoly(t,
	100.46645683 * 3600,
	1295977422.83429,
	-2.04411,
))

semi_major(::Mars, t) = @evalpoly(t,
    1.5236793419,
    3e-10,
)
eccentricity(::Mars, t) = @evalpoly(t,
	0.0934006477,
	0.0009048438,
	-80641e-10,
)
inclination(::Mars, t) = sec2rad(@evalpoly(t,
	1.84972648 * 3600,
	-293.31722,
	-8.11830,
))
ascending_node(::Mars, t) = sec2rad(@evalpoly(t,
	49.55809321 * 3600,
	-10620.90088,
	-230.57416,
))
longitude_of_perihelion(::Mars, t) = sec2rad(@evalpoly(t,
	336.06023395 * 3600,
	15980.45908,
	-62.32800,
))
mean_longitude(::Mars, t) = sec2rad(@evalpoly(t,
	355.43299958 * 3600,
	689050774.93988,
	0.94264,
))

semi_major(::Jupiter, t) = @evalpoly(t,
    5.2026032092,
    19132e-10,
    -39e-10,
)
eccentricity(::Jupiter, t) = @evalpoly(t,
	0.0484979255,
	0.0016322542,
	-0.0000471366,
)
inclination(::Jupiter, t) = sec2rad(@evalpoly(t,
	1.30326698 * 3600,
	-71.55890,
	11.95297,
))
ascending_node(::Jupiter, t) = sec2rad(@evalpoly(t,
	100.46440702 * 3600,
	6362.03561,
	326.52178,
))
longitude_of_perihelion(::Jupiter, t) = sec2rad(@evalpoly(t,
	14.33120687 * 3600,
	7758.75163,
	259.95938,
))
mean_longitude(::Jupiter, t) = sec2rad(@evalpoly(t,
	34.35151874 * 3600,
	109256603.77991,
	-30.60378,
))

semi_major(::Saturn, t) = @evalpoly(t,
    9.5549091915,
    -0.0000213896,
    444e-10,
)
eccentricity(::Saturn, t) = @evalpoly(t,
	0.0555481426,
	-0.0034664062,
	-0.0000643639,
)
inclination(::Saturn, t) = sec2rad(@evalpoly(t,
	2.48887878 * 3600,
	91.85195,
	-17.66225,
))
ascending_node(::Saturn, t) = sec2rad(@evalpoly(t,
	113.66550252 * 3600,
	-9240.19942,
	-66.23743,
))
longitude_of_perihelion(::Saturn, t) = sec2rad(@evalpoly(t,
	93.05723748 * 3600,
	20395.49439,
	190.25952,
))
mean_longitude(::Saturn, t) = sec2rad(@evalpoly(t,
	50.07744430 * 3600,
	43996098.55732,
	75.61614,
))

semi_major(::Uranus, t) = @evalpoly(t,
    19.2184460618,
    -3716e-10,
    979e-10,
)
eccentricity(::Uranus, t) = @evalpoly(t,
	0.0463812221,
	-0.0002729293,
	0.0000078913,
)
inclination(::Uranus, t) = sec2rad(@evalpoly(t,
	0.77319689 * 3600,
	-60.72723,
	1.25759,
))
ascending_node(::Uranus, t) = sec2rad(@evalpoly(t,
	74.00595701 * 3600,
	2669.15033,
	145.93964,
))
longitude_of_perihelion(::Uranus, t) = sec2rad(@evalpoly(t,
	173.00529106 * 3600,
	3215.56238,
	-34.09288,
))
mean_longitude(::Uranus, t) = sec2rad(@evalpoly(t,
    314.05500511 * 3600,
    15424811.93933,
    -1.75083,
))

semi_major(::Neptune, t) = @evalpoly(t,
    30.1103868694,
    -16635e-10,
    686e-10,
)
eccentricity(::Neptune, t) = @evalpoly(t,
	0.0094557470,
	0.0000603263,
)
inclination(::Neptune, t) = sec2rad(@evalpoly(t,
	1.76995259 * 3600,
	8.12333,
	0.08135,
))
ascending_node(::Neptune, t) = sec2rad(@evalpoly(t,
	131.78405702 * 3600,
	-221.94322,
	-0.78728,
))
longitude_of_perihelion(::Neptune, t) = sec2rad(@evalpoly(t,
	48.12027554 * 3600,
	1050.71912,
	27.39717,
))
mean_longitude(::Neptune, t) = sec2rad(@evalpoly(t,
	304.34866548 * 3600,
	7865503.20744,
	0.21103,
))

# Tables for trigonometric terms to be added to the mean elements of
# the semi-major axes

kp(::Mercury) = SVector(69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0)
kp(::Venus) = SVector(21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0)
kp(::EarthBarycenter) = SVector(16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0)
kp(::Mars) = SVector(6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0)
kp(::Jupiter) = SVector(1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454)
kp(::Saturn) = SVector(574, 0, 880, 287, 19, 1760, 1167, 306, 574)
kp(::Uranus) = SVector(204, 0, 177, 1265, 4, 385, 200, 208, 204)
kp(::Neptune) = SVector(0, 102, 106, 4, 98, 1367, 487, 204, 0)

ca(::Mercury) = SVector(4, -13, 11, -9, -9, -3, -1, 4, 0)
ca(::Venus) = SVector(-156, 59, -42, 6, 19, -20, -10, -12, 0)
ca(::EarthBarycenter) = SVector(64, -152, 62, -8, 32, -41, 19, -11, 0)
ca(::Mars) = SVector(124, 621, -145, 208, 54, -57, 30, 15, 0)
ca(::Jupiter) = SVector(-23437, -2634, 6601, 6259, -1507, -1821, 2620, -2115, -1489)
ca(::Saturn) = SVector(62911, -119919, 79336, 17814, -24241, 12068, 8306, -4893, 8902)
ca(::Uranus) = SVector(389061, -262125, -44088, 8387, -22976, -2093, -615, -9720, 6633)
ca(::Neptune) = SVector(-412235, -157046, -31430, 37817, -9740, -13, -7449, 9644, 0)

sa(::Mercury) = SVector(-29, -1, 9, 6, -6, 5, 4, 0, 0)
sa(::Venus) = SVector(-48, -125, -26, -37, 18, -13, -20, -2, 0)
sa(::EarthBarycenter) = SVector(-150, -46, 68, 54, 14, 24, -28, 22, 0)
sa(::Mars) = SVector(-621, 532, -694, -20, 192, -94, 71, -73, 0)
sa(::Jupiter) = SVector(-14614, -19828, -5869, 1881, -4372, -2255, 782, 930, 913)
sa(::Saturn) = SVector(139737, 0, 24667, 51123, -5102, 7429, -4095, -1976, -9566)
sa(::Uranus) = SVector(-138081, 0, 37205, -49039, -41901, -33872, -27037, -12474, 18797)
sa(::Neptune) = SVector(0, 28492, 133236, 69654, 52322, -49577, -26430, -3593, 0)

# Tables giving the trigonometric terms to be added to the mean
# elements of the mean longitudes

kq(::Mercury) = SVector(3086, 15746, 69613, 59899, 75645, 88306, 12661, 2658, 0, 0)
kq(::Venus) = SVector(21863, 32794, 10931, 73, 4387, 26934, 1473, 2157, 0, 0)
kq(::EarthBarycenter) = SVector(10, 16002, 21863, 10931, 1473, 32004, 4387, 73, 0, 0)
kq(::Mars) = SVector(10, 6345, 7818, 1107, 15636, 7077, 8184, 532, 10, 0)
kq(::Jupiter) = SVector(19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454)
kq(::Saturn) = SVector(19, 574, 287, 306, 1760, 12, 31, 38, 19, 574)
kq(::Uranus) = SVector(4, 204, 177, 8, 31, 200, 1265, 102, 4, 204)
kq(::Neptune) = SVector(4, 102, 106, 8, 98, 1367, 487, 204, 4, 102)

cl(::Mercury) = SVector(21, -95, -157, 41, -5, 42, 23, 30, 0, 0)
cl(::Venus) = SVector(-160, -313, -235, 60, -74, -76, -27, 34, 0, 0)
cl(::EarthBarycenter) = SVector(-325, -322, -79, 232, -52, 97, 55, -41, 0, 0)
cl(::Mars) = SVector(2268, -979, 802, 602, -668, -33, 345, 201, -55, 0)
cl(::Jupiter) = SVector(7610, -4997, -7689, -5841, -2617, 1115, -748, -607, 6074, 354)
cl(::Saturn) = SVector(-18549, 30125, 20012, -730, 824, 23, 1289, -352, -14767, -2062)
cl(::Uranus) = SVector(-135245, -14594, 4197, -4030, -5630, -2898, 2540, -306, 2939, 1986)
cl(::Neptune) = SVector(89948, 2103, 8963, 2695, 3682, 1648, 866, -154, -1963, -283)

sl(::Mercury) = SVector(-342, 136, -23, 62, 66, -52, -33, 17, 0, 0)
sl(::Venus) = SVector(524, -149, -35, 117, 151, 122, -71, -62, 0, 0)
sl(::EarthBarycenter) = SVector(-105, -137, 258, 35, -116, -88, -112, -80, 0, 0)
sl(::Mars) = SVector(854, -205, -936, -240, 140, -341, -97, -232, 536, 0)
sl(::Jupiter) = SVector(-56980, 8016, 1012, 1448, -3024, -3710, 318, 503, 3767, 577)
sl(::Saturn) = SVector(138606, -13478, -4964, 1441, -1319, -1482, 427, 1236, -9167, -1918)
sl(::Uranus) = SVector(71234, -41116, 5334, -4935, -1848, 66, 434, -1748, 3780, -701)
sl(::Neptune) = SVector(-47645, 11647, 2166, 3194, 679, 0, -244, -419, -2531, 48)

function state!(pos, vel, ::SimonBretagnon, ep::Epoch, ::Sun, body::CelestialBody)
    t = value(centuries(j2000(TDBEpoch(ep)))) / 10.0

    a = semi_major(body, t)
    e = eccentricity(body, t)
    i = inclination(body, t)
    Ω = mod2pi(ascending_node(body, t))
    ω′ = mod2pi(longitude_of_perihelion(body, t))
    λm = mean_longitude(body, t)

    dμ = 0.35953620 * t;

    kp0 = kp(body)[1:8]
    kq0 = kq(body)[1:8]
    ca0 = ca(body)[1:8]
    sa0 = sa(body)[1:8]
    cl0 = cl(body)[1:8]
    sl0 = sl(body)[1:8]

    arga = kp0 .* dμ
    argl = kq0 .* dμ
    a += sum((ca0 .* cos.(arga) .+ sa0 .* sin.(arga)) * 1e-7)
    λm += sum((cl0 .* cos.(argl) .+ sl0 .* sin.(argl)) * 1e-7)

    kp1 = kp(body)[9:end]
    kq1 = kq(body)[9:end]
    ca1 = ca(body)[9:end]
    sa1 = sa(body)[9:end]
    cl1 = cl(body)[9:end]
    sl1 = sl(body)[9:end]

    arga = kp1 .* dμ
    argl = kq1 .* dμ
    a += sum(t * (ca1 .* cos.(arga) .+ sa1 .* sin.(arga)) * 1e-7)
    λm += sum(t * (cl1 .* cos.(argl) .+ sl1 .* sin.(argl)) * 1e-7)
    λm = normalize_angle(λm, 0.0)

    M = λm - ω′
    ω = ω′ - Ω
    ν = transform(mean_anomaly, true_anomaly, M, e)

    μ = grav_param(sun)

    r, v = cartesian(a * astronomical_unit(), e, i, Ω, ω, ν, μ)

    rot = angleaxis_to_dcm(-ϵ, [1, 0, 0])
    pos .+= rot * r
    vel .+= rot * v

	return pos, vel
end

function position!(pos, ::SimonBretagnon, ep::Epoch, ::Sun, body::CelestialBody)
    state!(pos, zeros(3), simon_bretagnon, ep, sun, body)
    pos
end

function velocity!(vel, ::SimonBretagnon, ep::Epoch, ::Sun, body::CelestialBody)
    state!(zeros(3), vel, simon_bretagnon, ep, sun, body)
    vel
end

function state!(pos, vel, ::SimonBretagnon, ep::Epoch, from::CelestialBody, to::CelestialBody)
    rv1 = state!(zeros(3), zeros(3), simon_bretagnon, ep, sun, from)
    rv2 = state!(zeros(3), zeros(3), simon_bretagnon, ep, sun, to)
    pos .+= rv2[1] .- rv1[1]
    vel .+= rv2[2] .- rv1[2]
    pos, vel
end

function position!(pos, ::SimonBretagnon, ep::Epoch, from::CelestialBody, to::CelestialBody)
    state!(pos, zeros(3), simon_bretagnon, ep, from, to)
    pos
end

function velocity!(vel, ::SimonBretagnon, ep::Epoch, from::CelestialBody, to::CelestialBody)
    state!(zeros(3), vel, simon_bretagnon, ep, from, to)
    vel
end


text:
"""
    $(FUNCTIONNAME)(leaf::LeafProspectProProperties{FT},
                    optis) where {FT<:AbstractFloat}

Computes leaf optical properties (reflectance and transittance) based on PROSPECT-PRO
# Arguments
- `leaf`  [`LeafProspectProProperties`](@ref) type struct which provides leaf composition
- `optis` [`LeafOpticalProperties`](@ref) type struct, which provides absorption cross sections and spectral grid
# Examples
```julia-repl
julia> opti = createLeafOpticalStruct((400.0:5:2400)*u"nm");
julia> leaf = LeafProspectProProperties{Float64}(Ccab=30.0);
julia> T,R = prospect(leaf,opti);
````
"""
function prospect(
            leaf::LeafProspectProProperties{FT},
            optis
) where {FT<:AbstractFloat}
    # ***********************************************************************
    # Jacquemoud S., Baret F. (1990), PROSPECT: a model of leaf optical
    # properties spectra; Remote Sens. Environ.; 34:75-91.
    # Reference:
    # Féret, Gitelson, Noble & Jacquemoud [2017]. PROSPECT-D: Towards modeling
    # leaf optical properties through a complete lifecycle
    # Remote Sensing of Environment; 193:204215
    # DOI: http://doi.org/10.1016/j.rse.2017.03.004
    # The specific absorption coefficient corresponding to brown pigment is()
    # provided by Frederic Baret [EMMAH, INRA Avignon, baret@avignon.inra.fr]
    # & used with his autorization.
    # ***********************************************************************

    @unpack N, Ccab, Ccar, Cbrown, Canth, Cw, Cm, Cprot, Ccbc = leaf;
    @unpack Kcab, Kant, Kb, Kcar, Km, Kw, nᵣ, Kp, Kcbc = optis;

    # This can go into a separate multiple dispatch function as the rest remains constant across versions!
    Kall=(Ccab*Kcab + Ccar*Kcar + Canth*Kant + Cbrown*Kb + Cw*Kw + Cm*Km + Cprot*Kp +Ccbc * Kcbc) / N;

    # Adding eps() here to keep it stable and NOT set to 1 manually when Kall=0 (ForwardDiff won't work otherwise)
    tau = (FT(1) .-Kall).*exp.(-Kall) .+ Kall.^2 .*real.(expint.(Kall.+eps(FT)))

    # ***********************************************************************
    # reflectance & transmittance of one layer
    # ***********************************************************************
    # Allen W.A., Gausman H.W., Richardson A.J., Thomas J.R. (1969)
    # Interaction of isotropic ligth with a compact plant leaf; J. Opt.
    # Soc. Am., 59[10]:1376-1379.
    # ***********************************************************************
    # reflectivity & transmissivity at the interface
    #-------------------------------------------------
    # From Prospect-D, uses 40 here instead of 59 from CVT)

    #talf    = calctav.(59.,nr)
    talf    = calctav.(FT(40),nᵣ)
    ralf    = FT(1) .-talf

    t12     = calctav.(FT(90), nᵣ)
    r12     = FT(1) .-t12
    t21     = t12./(nᵣ.^2)
    r21     = FT(1) .-t21

    # top surface side
    denom   = FT(1) .-r21.*r21.*tau.^2
    Ta      = talf.*tau.*t21./denom
    Ra      = ralf.+r21.*tau.*Ta
    # bottom surface side
    t       = t12.*tau.*t21./denom
    r       = r12+r21.*tau.*t

    # ***********************************************************************
    # reflectance & transmittance of N layers
    # Stokes equations to compute properties of next N-1 layers [N real]
    # Normal case()
    # ***********************************************************************
    # Stokes G.G. (1862), On the intensity of the light reflected from
    # | transmitted through a pile of plates; Proc. Roy. Soc. Lond.
    # 11:545-556.
    # ***********************************************************************
    D       = sqrt.(((FT(1) .+r.+t).*(FT(1) .+r.-t).*(FT(1) .-r.+t).*(FT(1) .-r.-t)).+5eps(FT))
    #println(typeof(D), typeof(r), typeof(t))
    rq      = r.^2
    tq      = t.^2
    a       = (FT(1) .+rq.-tq.+D)./(2r)
    b       = (FT(1) .-rq.+tq.+D)./(2t)

    bNm1    = b.^(N-1);                  #
    bN2     = bNm1.^2
    a2      = a.^2
    denom   = a2.*bN2.-1
    Rsub    = a.*(bN2.-1)./denom
    Tsub    = bNm1.*(a2.-1)./denom

    # Case of zero absorption
    j       = findall(r.+t .>= 1)
    Tsub[j] = t[j]./(t[j]+(1 .-t[j])*(leaf.N-1))
    Rsub[j] = 1 .-Tsub[j]

    # Reflectance & transmittance of the leaf: combine top layer with next N-1 layers
    denom   = FT(1) .-Rsub.*r
    # lambertian Tranmsission
    T    = Ta.*Tsub./denom
    # lambertian Reflectance
    R    = Ra.+Ta.*Rsub.*t./denom

    return T,R
end

"""
    calctav(α::FT, nr::FT) where {FT<:AbstractFloat}

Computes transmission of isotropic radiation across a dielectric surface
(Stern F., 1964; Allen W.A.,Appl. Opt., 3(1):111-113 1973)).
From calctav.m in PROSPECT-D
# Arguments
- `α` angle of incidence [degrees]
- `nr` Index of refraction
"""
function calctav(α::FT, nr::FT) where {FT<:AbstractFloat}
    a   = ((nr+1) ^ 2) / 2;
    a3  = a  ^ 3;
    n2  = nr ^ 2;
    n4  = nr ^ 4;
    n6  = nr ^ 6;
    np  = n2 + 1;
    np2 = np ^ 2;
    np3 = np ^ 3;
    nm2 = (n2 - 1) ^2;
    k   = ((n2-1) ^ 2) / -4;
    k2  = k  ^ 2;
    sa2 = sind(α) ^ 2;

    _b1 = (α==90 ? 0 : sqrt( (sa2 - np/2)^2 + k ));
    _b2 = sa2 - np/2;
    b   = _b1 - _b2;
    b3  = b ^ 3;
    ts  = ( k2 / (6*b3) + k/b - b/2 ) - ( k2 / (6*a3) + k/a - a/2 );
    tp1 = -2 * n2 * (b-a) / (np2);
    tp2 = -2 * n2 * np * log(b/a) / (nm2);
    tp3 = n2 * (1/b - 1/a) / 2;
    tp4 = 16 * n4 * (n4+1) * log((2*np*b - nm2) / (2*np*a - nm2)) / (np3*nm2);
    tp5 = 16 * n6 * (1 / (2*np*b-nm2) - 1 / (2*np*a-nm2)) / (np3);
    tp  = tp1 + tp2 + tp3 + tp4 + tp5;
    tav = (ts + tp) / (2 * sa2);

    return tav
end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Host introspection

print("This host's word size is ", WORD_SIZE, ".")
if ENDIAN_BOM == 0x04030201
    println("And it is a little-endian machine.")
elseif ENDIAN_BOM == 0x01020304
    println("And it is a big-endian machine.")
else
    println("ENDIAN_BOM = ", ENDIAN_BOM, ", which is confusing")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Documenter
using EvolutionaryModelingTools

makedocs(;
    sitename="EvolutionaryModelingTools.jl",
    pages=["index.md", "example.md", "reference.md"],
)

deploydocs(;
    repo="github.com/wangl-cc/EvolutionaryModelingTools.jl.git",
    push_preview=true,
)


text:
struct Context{Ns, Nb}
    β :: Float64

    system_space :: StateSpace{Ns}
    system_state :: Vector{Float64}

    bath_space :: StateSpace{Nb}
    bath_state :: Vector{Float64}

    degeneracies :: Vector{NTuple{2, NTuple{2, Int}}}
    P :: Vector{BigFloat}
    Pid :: BigFloat
    Pcum :: Vector{BigFloat}

    function Context{Ns, Nb}(β::Float64, system_space::StateSpace{Ns}, system_state::AbstractVector,
                             bath_space::StateSpace{Nb}) where{Ns, Nb}
        bath_state = boltzmann(β, bath_space)
        deg, P, Pid = degeneracies(system_space, bath_space)
        Pcum = cumsum(P) + Pid
        new(β, system_space, system_state, bath_space, bath_state, deg, P, Pid, Pcum)
    end
end

function Context(β::Float64, system_space::StateSpace{N}, system_state::AbstractVector) where {N}
    Context{N,N}(β, system_space, system_state, system_space)
end

function timestep(rng::AbstractRNG, ctx::Context)
    p, i = BigFloat(rand(rng)), 0
    if p > ctx.Pid
        i = 1
        while i < length(ctx.Pcum) && p > ctx.Pcum[i]
            i += 1
        end
        if i <= length(ctx.P)
            ((a,b), (u,v)) = ctx.degeneracies[i]
            Pa, Pb = ctx.system_state[a], ctx.bath_state[b]
            Pu, Pv = ctx.system_state[u], ctx.bath_state[v]
            degab = ctx.system_space.deg[a] * ctx.bath_space.deg[b]
            deguv = ctx.system_space.deg[u] * ctx.bath_space.deg[v]
            dP = ((Pu*Pv)/deguv) - ((Pa*Pb)/degab)
            ctx.system_state[a] += dP
            ctx.system_state[u] -= dP
        end
    end
    ctx.system_state, i
end

timestep(ctx::Context) = timestep(Base.GLOBAL_RNG, ctx)


text:
"""
    checkerboard(shape)

Generates a checkerboard level set function.
According to Pascal Getreuer, such a level set function has fast
convergence.

# Arguments
- shape: 2D shape of the output array
"""
function checkerboard(shape)
    x₀ = reshape(collect(0:(shape[1] - 1)), shape[1], 1)
    y₀ = reshape(collect(0:(shape[2] - 1)), 1, shape[2])
    return 𝚽₀ = sin.(pi / 5 .* x₀) .* sin.(pi / 5 .* y₀)
end

"""
    disk(shape)

Generates a disk level set function.
The disk covers the whole image along its smallest dimension
or `factor` can be given to decrease the size of the disk.

# Arguments
- shape: 2D shape of the output array
- factor: size to decrease the disk, default is
    1 which implies no size change
"""
function disk(shape; factor=1)
    res = ones(shape)
    centerY = Int(round((shape[1] - 1) / 2))
    centerX = Int(round((shape[2] - 1) / 2))
    res[centerY, centerX] = 0.0
    radius = Float32(min(centerX, centerY)) / factor
    return disk = (radius .- DistanceTransforms.euclidean(res)) / radius
end


text:
function K_eular_shear(beam::Beam)::Matrix{Float64}
    E,ν=beam.material.E,beam.material.ν
    A,I₂,I₃,J,l=beam.section.A,beam.section.I₂,beam.section.I₃,beam.section.J,beam.l
    As₂,As₃=beam.section.As₂,beam.section.As₃
    G=E/2/(1+ν)
    ϕ₂,ϕ₃=12E*I₃/(G*As₂*l^2),12E*I₂/(G*As₃*l^2)
    K=zeros(12,12)
    K[1,1]=E*A/l
    K[2,2]=12E*I₃/l^3/(1+ϕ₂)
    K[3,3]=12E*I₂/l^3/(1+ϕ₃)
    K[4,4]=G*J/l
    K[5,5]=(4+ϕ₃)*E*I₂/l/(1+ϕ₃)
    K[6,6]=(4+ϕ₂)*E*I₃/l/(1+ϕ₂)
    K[7,7]=E*A/l
    K[8,8]=12E*I₃/l^3/(1+ϕ₂)
    K[9,9]=12E*I₂/l^3/(1+ϕ₃)
    K[10,10]=G*J/l
    K[11,11]=(4+ϕ₃)*E*I₂/l/(1+ϕ₃)
    K[12,12]=(4+ϕ₂)*E*I₃/l/(1+ϕ₂)

    K[3,5]=K[5,3]=-6E*I₂/l^2/(1+ϕ₃)
    K[6,8]=K[8,6]=-6E*I₃/l^2/(1+ϕ₂)
    K[9,11]=K[11,9]=6E*I₂/l^2/(1+ϕ₃)

    K[2,6]=K[6,2]=6E*I₃/l^2/(1+ϕ₂)
    K[5,9]=K[9,5]=6E*I₂/l^2/(1+ϕ₃)
    K[8,12]=K[12,8]=-6E*I₃/l^2/(1+ϕ₂)

    K[7,1]=K[1,7]=-E*A/l
    K[8,2]=K[2,8]=-12E*I₃/l^3/(1+ϕ₂)
    K[9,3]=K[3,9]=-12E*I₂/l^3/(1+ϕ₃)
    K[10,4]=K[4,10]=-G*J/l
    K[11,5]=K[5,11]=(2-ϕ₃)*E*I₂/l/(1+ϕ₃)
    K[12,6]=K[6,12]=(2-ϕ₂)*E*I₃/l/(1+ϕ₂)

    K[3,11]=K[11,3]=-6E*I₂/l^2/(1+ϕ₃)

    K[2,12]=K[12,2]=6E*I₃/l^2/(1+ϕ₂)

    return K
end

function K2_eular_shear(beam::Beam)::Matrix{Float64}
    E,ν=beam.material.E,beam.material.ν
    A,I₂,I₃,J,l=beam.section.A,beam.section.I₂,beam.section.I₃,beam.section.J,beam.l
    As₂,As₃=beam.section.As₂,beam.section.As₃
    G=E/2/(1+ν)
    T=σ*A
    ϕ₂,ϕ₃=12E*I₃/(G*As₂*l^2),12E*I₂/(G*As₃*l^2)
    K=zeros(12,12)

    K[2,2]=(6/5+2ϕ₂+ϕ₂^2)/(1+ϕ₂)^2
    K[3,3]=(6/5+2ϕ₃+ϕ₃^2)/(1+ϕ₃)^2
    K[4,4]=J/A
    K[5,5]=(2*l^2/15+l^2*ϕ₃/6+l^2*ϕ₃^2/12)/(1+ϕ₃)^2
    K[6,6]=(2*l^2/15+l^2*ϕ₂/6+l^2*ϕ₂^2/12)/(1+ϕ₂)^2
    K[8,8]=(6/5+2ϕ₂+ϕ₂^2)/(1+ϕ₂)^2
    K[9,9]=(6/5+2ϕ₃+ϕ₃^2)/(1+ϕ₃)^2
    K[10,10]=J/A
    K[11,11]=(2*l^2/15+l^2*ϕ₃/6+l^2*ϕ₃^2/12)/(1+ϕ₃)^2
    K[12,12]=(2*l^2/15+l^2*ϕ₂/6+l^2*ϕ₂^2/12)/(1+ϕ₂)^2

    K[3,5]=K[5,3]=-(l/10)/(1+ϕ₃)^2
    K[6,8]=K[8,6]=-(l/10)/(1+ϕ₂)^2
    K[9,11]=K[11,9]=(l/10)/(1+ϕ₃)^2

    K[2,6]=K[6,2]=(l/10)/(1+ϕ₂)^2
    K[5,9]=K[9,5]=(l/10)/(1+ϕ₃)^2

    K[2,8]=K[8,2]=-(6/5+2ϕ₂+ϕ₂^2)/(1+ϕ₂)^2
    K[3,9]=K[9,3]=-(6/5+2ϕ₃+ϕ₃^2)/(1+ϕ₃)^2
    K[4,10]=K[10,4]=-J/A
    K[5,11]=K[11,5]=-(l^2/30+l^2*ϕ₃/6+l^2*ϕ₃^2/12)/(1+ϕ₃)^2
    K[6,12]=K[12,6]=-(l^2/30+l^2*ϕ₂/6+l^2*ϕ₂^2/12)/(1+ϕ₂)^2

    K[3,11]=K[11,3]=-(l/10)/(1+ϕ₃)^2

    K[2,12]=K[12,2]=(l/10)/(1+ϕ₂)^2

    K*=T/l
    return K
end


text:
"""
	prep_data_granger(data, seglen, verbose)

Checks and prepares data shape.

### Arguments

  - `data::AbstractArray`: Nx2 array for N data points in 2 channels.
  - `seglen::Integer`: segment length.
  - `verbose::Bool`: if `true`, warnings and info logs would be echoed.

### Returns

  - `data::AbstractArray`: NxM array for N data points in M channels.
"""
function prep_data_granger(data::AbstractArray, seglen::Integer, verbose::Bool)

    if ndims(data) != 2  # data dimension
        data = squeeze(data)
        ndims(data) != 2 && throw(DimensionMismatch("Data must be a 2D-array!"))
        verbose && @info "Data is squeezed to a 2D-array)"
    end
    if size(data, 1) < size(data, 2)  # should be NxM array for N data points in M channels
        verbose && @info "Data is transposed to (#samples, #channels)"
        data = reshape(data, size(data, 2), size(data, 1))
    end
    if size(data, 2) != 2  # should be exactly 2 channels
        throw("Only 2 channels are supported!")
    end
    if size(data, 1) < seglen  # seglen must be smaller than number of samples
        throw(DimensionMismatch("seglen must be smaller than number of samples!"))
    end
    return data
end

"""
	granger_est(data, seglen; order, method, verbose)

Granger's method of causal relation approximation

### Arguments

  - `data::AbstractArray`: Nx2 array for N data points in 2 channels.
  - `seglen::Integer`: segment length.

*optional arguments*

  - `method::String`: standard deviation estimation method (default is `"jackknife"`)
  - `order::Int`: Model order. Assumed time delay order of interest (default is `n_samples / 2`)
  - `verbose::Bool`: if `true`, warnings and info logs would be echoed.

### Returns

  - `Grind::Float64`: Granger causality index
  - `Grind_std::Float64`: estimated standard deviation of error

### Internal variables:

  - Covs: Concatenated covariance mats of different orders
  - Acoef: Design (A) matrix coefficients
  - Perr: Prediction Error
  - Σ: Measures the accuracy of the autoregressive prediction

### References:
  - C. W. J. Granger (1969). Investigating causal relations by econometric models
    and cross-spectral methods, Econometrica, 37, 424-438.
  - M. Ding et al. (2006). Granger Causality: Basic Theory and Application to Neuroscience
    available online: https://arxiv.org/abs/q-bio/0608035

"""
function granger_est(
    data::Array{Float64,2},
    seglen::Int;
    order::Int = 0,
    method::String = "none",
    verbose::Bool = true,
)

    data = prep_data_granger(data, seglen, verbose)

    if order == 0
        order = int(seglen / 2) - 1
    end

    Covs12 = est_sig_covs(data, order)
    Covs1 = est_sig_covs(review(data, 1), order)
    Covs2 = est_sig_covs(review(data, 2), order)

    Acoef12, Perr12 = mvar_est(Covs12)
    Acoef1, _ = mvar_est(Covs1)
    Acoef2, _ = mvar_est(Covs2)

    Σ12 = Perr12[:, end-1:end]
    Σ12 = est_noise_covs(data, Acoef12, seglen)
    Σ1 = est_noise_covs(review(data, 1), Acoef1, seglen)
    Γ1 = est_noise_covs(review(data, 2), Acoef2, seglen)
    Σ2 = Σ12[1, 1]
    Γ2 = Σ12[2, 2]

    X_to_Y = log(Γ1 / Γ2)
    Y_to_X = log(Σ1 / Σ2)

    Grind = (X_to_Y-Y_to_X)[1]

    if lowercase(method) == "jackknife"
        Grind_std = granger_jackknife(data, seglen, order, Covs12, Covs1, Covs2)
        return Grind, Grind_std
    else
        return Grind
    end

end

"""
	est_sig_covs(data, order)

Concatenated covariances of signal to the given order

### Arguments

  - `data::AbstractArray`: Nx2 array for N data points in 2 channels.
  - `order::Int`: Model order (default is nsamples/2 )

### Returns

  - `Covs::Array{Float64, 2}`: Concatenated covariances of 0 to the given order

"""
function est_sig_covs(data::AbstractArray{Float64,2}, order::Int)::Array{Float64,2}

    nsamples, nchan = size(data) # number of rows (samples) and columns (channels)

    Covs = Array{Float64}(undef, nchan, nchan * (order + 1))

    Covs[:, 1:nchan] = covariance(data)
    for k = 1:order
        @views Covs[:, (1:nchan).+k*nchan] =
            covariance(data[k+1:nsamples, :], data[1:nsamples-k, :])
    end
    return Covs
end

"""
	mvar_est(Covs)

MultiVariate AutoRegressive (MVAR) model parameter estimation using
"Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function"

### Arguments

  - `Covs::Array{Float64, 2}`: Concatenated covariances of 0 to the given order

### Returns

  - `ARF::Array{Float64, 2}`: Design (A) matrix coefficients
  - `PE::Array{Float64}` : model Prediction Error

### References:
  - Alois Schlögl et al. <a.schloegl@ieee.org>, part of TSA-toolbox and FieldTrip-toolbox:
    https://github.com/fieldtrip/fieldtrip/blob/master/external/biosig/mvar.m

"""
function mvar_est(Covs::Array{Float64,2})

    n_i, m_j = size(Covs)  # number of rows and columns
    order = int(m_j / n_i) - 1

    ARF = Array{Float64}(undef, n_i, n_i * order)
    ARB = Array{Float64}(undef, n_i, n_i * order)
    RCF = Array{Float64}(undef, n_i, n_i * order)
    RCB = Array{Float64}(undef, n_i, n_i * order)
    PE = Array{Float64}(undef, n_i, n_i * (order + 1))

    PE[:, 1:n_i] = Covs[:, 1:n_i]
    PEF = Covs[:, 1:n_i]  # it's a different normalization in BioSig
    PEB = Covs[:, 1:n_i]  # it's a different normalization in BioSig

    for k = 1:order
        S1 = k * n_i .+ (1:n_i)  # slice
        S2 = k * n_i .+ (1-n_i:0)  # slice

        D = Covs[:, S1]
        for L = 1:k-1
            S3 = L * n_i .+ (1-n_i:0)  # slice
            D -= ARF[:, S3] * Covs[:, (k-L)*n_i.+(1:n_i)]
        end

        ARF[:, S2] = D / PEB
        ARB[:, S2] = D' / PEF
        for L = 1:k-1
            S3 = L * n_i .+ (1-n_i:0)  # slice
            S4 = (k - L) * n_i .+ (1-n_i:0)  # slice
            ARB[:, S4], ARF[:, S3] =
                (ARB[:, S4] - ARB[:, S2] * ARF[:, S3], ARF[:, S3] - ARF[:, S2] * ARB[:, S4])
        end

        RCF[:, S2] = ARF[:, S2]
        RCB[:, S2] = ARB[:, S2]
        PEF = (identity(n_i) - ARF[:, S2] * ARB[:, S2]) * PEF
        PEB = (identity(n_i) - ARB[:, S2] * ARF[:, S2]) * PEB
        PE[:, S1] = PEF
    end

    return ARF, PE
end

"""
	est_noise_covs(data, A, seglen)

Estimation of covariance of noise

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels.
  - `A::Array{Float64, 2}`: Design (A) matrix coefficients
  - `seglen::Integer`: segment length.

### Returns

  - `C::Array{Float64}` : covariance matrix of noise

"""
function est_noise_covs(
    data::AbstractArray{Float64,2},
    A::AbstractArray{Float64,2},
    seglen::Int,
)

    nsamples, nchan = size(data)  # number of rows and columns
    order = int(size(A, 2) / nchan)  # number of columns = order * number of channels
    nsegments = int(nsamples / seglen)  # number of segments

    A = reshape(A, nchan, nchan, :)
    C = zeros(nchan, nchan)

    @inbounds for e = 1:nsegments
        S1 = (e-1)*seglen+1:e*seglen
        S2 = order+1:seglen
        @views X = data[S1, :][S2, :]'
        for i = 1:order
            @views X -= A[:, :, i] * data[S1, :][S2.-i, :]'
        end
        C += X * X'
    end

    C /= nsegments * (seglen - order)

    return C
end

"""
	granger_jackknife(data, seglen, order, Covs12, Covs1, Covs2)

Jackknife sampling method for estimating the error standard deviation

### Arguments

  - `data::AbstractArray`: Nx2 array for N data points in 2 channels.
  - `order::Int`: Model order. Assumed time delay order of interest.
  - `seglen::Integer`: segment length.
  - `Covs12::Array{Float64, 2}`: Concatenated COvariances of 0 to the given order
  - `Covs1::Array{Float64, 2}`: Channel 1 Concatenated variances of 0 to the given order
  - `Covs2::Array{Float64, 2}`: Channel 2 Concatenated variances of 0 to the given order

### Returns

  - `::Float64`: estimated standard deviation of error
"""
function granger_jackknife(
    data::Array{Float64,2},
    seglen::Int,
    order::Int,
    Covs12::AbstractArray{Float64,2},
    Covs1::AbstractArray{Float64,2},
    Covs2::AbstractArray{Float64,2},
)

    nsamples = size(data, 1)  # number of samples
    nsegments = int(nsamples / seglen)  # number of segments
    trials = zeros(nsegments, 1)

    Threads.@threads for i = 1:nsegments
        Segment = view(data, (i-1)*seglen+1:i*seglen, :)

        Covs12_ = est_sig_covs(Segment, order)
        Covs1_ = est_sig_covs(review(Segment, 1), order)
        Covs2_ = est_sig_covs(review(Segment, 2), order)

        Covs12_s = (nsegments * Covs12 - Covs12_) / (nsegments - 1)
        Covs1_s = (nsegments * Covs1 - Covs1_) / (nsegments - 1)
        Covs2_s = (nsegments * Covs2 - Covs2_) / (nsegments - 1)

        Acoef12_s, _ = mvar_est(Covs12_s)
        Acoef1_s, _ = mvar_est(Covs1_s)
        Acoef2_s, _ = mvar_est(Covs2_s)

        SegmentC =
            vcat(view(data, 1:(i-1)*seglen, :), view(data, (i+1)*seglen.+1:nsamples, :))

        Σ12_s = est_noise_covs(SegmentC, Acoef12_s, seglen)
        Σ1_s = est_noise_covs(review(SegmentC, 1), Acoef1_s, seglen)
        Γ1_s = est_noise_covs(review(SegmentC, 2), Acoef2_s, seglen)
        Σ2_s = Σ12_s[1, 1]
        Γ2_s = Σ12_s[2, 2]

        X_to_Y_s = log(Γ1_s / Γ2_s)
        Y_to_X_s = log(Σ1_s / Σ2_s)

        trials[i, :] = X_to_Y_s - Y_to_X_s

    end
    return sqrt(nsegments) * std(trials)
end

"""
	granger_aic(data, seglen, order_range)

	Akaike Information Criterion
"""
function granger_aic(data::Array{Float64,2}, seglen::Int, order_range::UnitRange{Int64})

    nchan = size(data, 2)
    aic_range = Array{Float64}(undef, size(order_range))

    for (i, order) in enumerate(order_range)
        Covs = est_sig_covs(data, order)
        Acoef, _ = mvar_est(Covs)
        Σ = est_noise_covs(data, Acoef, seglen)
        aic_range[i] = 2 * log(det(Σ)) + (2 * order * nchan^2 / seglen)
    end
    return aic_range
end

"""
	granger_bic(data, seglen, order_range)

Bayesian Information Criterion
"""
function granger_bic(data::Array{Float64,2}, seglen::Int, order_range::UnitRange{Int64})

    nchan = size(data, 2)
    bic_range = Array{Float64}(undef, size(order_range))

    for (i, order) in enumerate(order_range)
        Covs = est_sig_covs(data, order)
        Acoef, _ = mvar_est(Covs)
        Σ = est_noise_covs(data, Acoef, seglen)
        bic_range[i] = 2 * log(det(Σ)) + (2 * order * nchan^2 * log(seglen) / seglen)
    end
    return bic_range
end


text:
struct Replicate_Report
    ids::Vector{Chain_ID} #replicates
    state_vecs::Dict{Chain_ID,Matrix{Float64}} #vectors of autotransition probabilities evolving by iterate, concatanated to matrices, indexed by Chain_ID
    emission_arrays::Dict{Chain_ID,Array{Float64}} #iterate*symbol*state_vecs
    sorted_id1_states::Vector{Integer}
    sort_dicts::Dict{Chain_ID,Dict{Integer,Integer}}#emission channels sorted on the basis of euclidean closeness to ids[1]
    sorted_symbols::Dict{Integer,Vector{Integer}}

    Replicate_Report(ids,state_vecs,emission_arrays, sorted_id1_states,sort_dicts,sorted_symbols) = assert_repreport(ids) && new(ids,state_vecs,emission_arrays,sorted_id1_states,sort_dicts,sorted_symbols)
end

function assert_repreport(ids::Vector{Chain_ID})
    Ks=Vector{Int64}()
    orders=Vector{Int64}()
    replicates=Vector{Int64}()
    obsids=Vector{String}()

    for id in ids
        push!(Ks,id.K)
        push!(orders,id.order)
        push!(obsids,id.obs_id)
        push!(replicates, id.replicate)
    end

    length(unique(Ks)) > 1 && throw(ArgumentError("Replicate set includes chains with different K state numbers! All chains must have HMMs with the same number of states."))
    length(unique(orders)) > 1 && throw(ArgumentError("Replicate set includes chains with different order coding numbers! All chains must have HMMs with the same order coding."))
    !allunique(replicates) && throw(ArgumentError("Replicate set includes chains with the same replicate #. Replicates should be unique!"))

    return true
end

function report_replicates(repset::Vector{Chain_ID},chains::Dict{Chain_ID,Vector{EM_step}})
    assert_repreport(repset) #check that the repset will pass its constructor before doing the work
    emission_arrays=Dict{Chain_ID,AbstractArray{AbstractFloat}}()
    state_arrays=Dict{Chain_ID,AbstractArray{AbstractFloat}}()

    for id in repset
        emission_array=zeros(length(chains[id]),4^(id.order+1),id.K)
        state_array=zeros(length(chains[id]),id.K)
        for (it, step) in enumerate(chains[id])
            for k in 1:id.K
                emission_array[it,:,k]=step.hmm.B[k].p
                state_array[it,k]=step.hmm.A[k,k]
            end
        end
        emission_arrays[id]=emission_array
        state_arrays[id]=state_array
    end

    println(repset)
    println(emission_arrays)
    println(state_arrays)
    sort_dicts,sorted_id1_states,sorted_symbols=sort_emitters_by_distance!(repset,emission_arrays)

    return Replicate_Report(repset,state_arrays,emission_arrays,sorted_id1_states, sort_dicts,sorted_symbols)
end

function sort_emitters_by_distance!(repset,emission_arrays)
    id1=repset[1]
    length(repset)==1 && (return Dict{Chain_ID,Dict{Integer,Integer}}(),[1:id1.K...],[1:4^id1.order+1])

    final_emissions = zeros(length(repset),size(emission_arrays[id1])[2:3]...)
    for (n,id) in enumerate(repset)
        final_emissions[n,:,:] = emission_arrays[id][end,:,:]
    end

    distances=zeros(length(repset)-1,id1.K,id1.K)
    for rep in 1:length(repset)-1
        for k in 1:id1.K, j in 1:id1.K
            distances[rep,k,j]=euclidean(final_emissions[1,:,k],final_emissions[1+rep,:,j])
        end
    end

    sortdicts=Dict{Chain_ID,Dict{Integer,Integer}}()
    sorted_id1_states=Vector{Integer}()
    sorted_symbols=Dict{Integer,Vector{Integer}}()
    for rep in 1:length(repset)-1
        sortdicts[repset[rep+1]]=Dict{Integer,Integer}()
    end

    for i in 1:id1.K
        id1_mindist_K=findmin([sum([minimum(distances[rep,k,:]) for rep in 1:length(repset)-1]) for k in 1:id1.K])[2]#find the state from the first replicate that has minimum cumulative distance to the closest state in the other replicates, excluding already-chosen states by spiking their values
        push!(sorted_id1_states,id1_mindist_K)
        for rep in 1:length(repset)-1
            rep_k=findmin(distances[rep,id1_mindist_K,:])[2]
            sortdicts[repset[rep+1]][id1_mindist_K]=rep_k
            distances[rep,:,rep_k].=1.0 #mask
            distances[1,id1_mindist_K,:].=1.0 #mask
        end

        symbol_distances=sum(hcat([euclidean.(final_emissions[1,:,id1_mindist_K],final_emissions[n+1,:,sortdicts[id][id1_mindist_K]]) for (n,id) in enumerate(repset[2:end])]...),dims=2)
        sorted_symbols[id1_mindist_K]=sortperm(symbol_distances[:,1])
    end

    return sortdicts,sorted_id1_states,sorted_symbols
end

function Base.show(io::IO, report::Replicate_Report; top_states::Integer=report.ids[1].K, top_symbols::Integer=2)
    printstyled(io, "REPLICATE CONVERGENCE REPORT\n", color=:green)
    println(" ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶ ̶̶ ̶ ̶ ̶")
    for i in 1:top_states
        id=report.ids[1]
        id1_state=report.sorted_id1_states[i]
        itmax=maximum([size(report.state_vecs[id],1) for id in report.ids])

        autoplt=lineplot(report.state_vecs[id][:,id1_state],
                        title="Autotransition prob. convergence",
                        name="$(report.ids[1]) K$id1_state",
                        xlim=(0,ceil(itmax, sigdigits=2)),
                        ylim=(0,1),
                        xlabel="iterate",
                        ylabel="p")

        symplot=lineplot(report.emission_arrays[id][:,report.sorted_symbols[id1_state][1],id1_state],
                        report.emission_arrays[id][:,report.sorted_symbols[id1_state][2],id1_state],
                        title="Symbol prob. convergence",
                        name="$(report.ids[1]) K$id1_state",
                        xlim=(0,1),
                        ylim=(0,1),
                        xlabel="Symbol 1 p",
                        ylabel="S2 p")

        for (n,id) in enumerate(report.ids)
            if n > 1
                lineplot!(autoplt, report.state_vecs[id][:,report.sort_dicts[id][id1_state]],
                            name="$(report.ids[n]) K$(report.sort_dicts[id][id1_state])")
                lineplot!(symplot,
                            report.emission_arrays[id][:,report.sorted_symbols[id1_state][1],report.sort_dicts[id][id1_state]],
                            report.emission_arrays[id][:,report.sorted_symbols[id1_state][2],report.sort_dicts[id][id1_state]],
                            name="$(report.ids[n]) K$(report.sort_dicts[id][id1_state])")
            end
        end

        show(autoplt)
        println("\n")
        show(symplot)
        println("\n")
    end
end


text:
using Documenter, StochasticPrograms

makedocs(sitename = "StochasticPrograms.jl",
         clean = false,
         pages = [
             "Home" => "index.md",
             "Manual" => Any[
                 "Quick start" => "manual/quickstart.md",
                 "Stochastic data" => "manual/data.md",
                 "Stochastic models" => "manual/model.md",
                 "Distributed stochastic programs" => "manual/distributed.md",
                 "Structured solvers" => "manual/structuredsolvers.md",
                 "Examples" => "manual/examples.md",
             ],
             "Library" =>  Any[
                 "Public interface" => "library/public.md",
                 "Solver interface" => "library/solverinterface.md"
             ]
         ])

deploydocs(
    repo = "github.com/martinbiel/StochasticPrograms.jl.git",
)


text:
export iauFk54z
"""
Convert a J2000.0 FK5 star position to B1950.0 FK4, assuming zero
proper motion in FK5 and parallax.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   - r2000,d2000    double   J2000.0 FK5 RA,Dec (rad)
   - bepoch         double   Besselian epoch (e.g. 1950.0)

Returned:
   - r1950,d1950    double   B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH
   - dr1950,dd1950  double   B1950.0 FK4 proper motions (rad/trop.yr)

Notes:

   1. In contrast to the iauFk524  routine, here the FK5 proper
   motions, the parallax and the radial velocity are presumed zero.

   2. This function converts a star position from the IAU 1976 FK5
   (Fricke) system to the former FK4 (Bessel-Newcomb) system, for
   cases such as distant radio sources where it is presumed there is
   zero parallax and no proper motion.  Because of the E-terms of
   aberration, such objects have (in general) non-zero proper motion
   in FK4, and the present routine returns those fictitious proper
   motions.

   3. Conversion from B1950.0 FK4 to J2000.0 FK5 only is provided for.
   Conversions involving other equinoxes would require additional
   treatment for precession.

   4. The position returned by this routine is in the B1950.0 FK4
   reference system but at Besselian epoch BEPOCH.  For comparison
   with catalogs the BEPOCH argument will frequently be 1950.0. (In
   this context the distinction between Besselian and Julian epoch
   is insignificant.)

   5. The RA component of the returned (fictitious) proper motion is
   dRA/dt rather than cos(Dec)*dRA/dt.

Called:
   - iauAnp       normalize angle into range 0 to 2pi
   - iauC2s       p-vector to spherical
   - iauFk524     FK4 to FK5
   - iauS2c       spherical to p-vector

This revision:   2018 December 5

SOFA release 2019-07-22

Copyright (C) 2019 IAU SOFA Board.  See notes at end.
"""

# void iauFk54z(double r2000, double d2000, double bepoch,
#               double *r1950, double *d1950,
#               double *dr1950, double *dd1950)

function iauFk54z(r2000::Real, d2000::Real, bepoch::Real)
   ref_r1950  = Ref{Float64}(0.0)
   ref_d1950  = Ref{Float64}(0.0)
   ref_dr1950 = Ref{Float64}(0.0)
   ref_dd1950 = Ref{Float64}(0.0)

   ccall((:iauFk54z, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, ),
         convert(Float64, r2000), convert(Float64, d2000),
         convert(Float64, bepoch),
         ref_r1950, ref_d1950, ref_dr1950, ref_dd1950)

   return ref_r1950[], ref_d1950[], ref_dr1950[], ref_dd1950[]
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module FreeBodyRotation

DAE model of a body with unconstrained rotational motion in 3D (using quaternions;)

DAE equations:
  u := A.*sin(2*pi*freqHz*time + phase);
  0 = w - 2*([Q[4],  Q[3], -Q[2], -Q[1];
             -Q[3],  Q[4],  Q[1], -Q[2];
              Q[2], -Q[1],  Q[4], -Q[3]])*der(Q);
  0 = I*der(w) + cross(w, I*w) - u;
  0 = dot(Q,Q)-1.0;

Initialization (inconsistent values)
  Q[1] = 0.1
  Q[2] = 0.5
  Q[3] = 0
  Q[4] = 1.0
  w = zeros(3)

Arguments of getModelResidues! function:
   x    = [Q     ; w]
   derx = [der(Q); der(w)]

For testing:
   Modelica model ModelicaReferenceModels.ODAEs.FreeBody
"""
module FreeBodyRotation

import ..ModiaMath
using  ..LinearAlgebra
using  ..StaticArrays

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    m::Float64
    I::SMatrix{3,3,Float64,9}
    A::SVector{3,Float64}
    freqHz::SVector{3,Float64}
    phase::SVector{3,Float64}

    function Model(;m=1.0, I=Diagonal([1.0,2.0,3.0]), A=[3.0,4.0,5.0], freqHz=[0.3,0.2,0.1],
                  phase=[0,0.5235987755983,1.0471975511966],
                  Q0=[0.1, 0.5, 0.0, 1.0],
                  w0=zeros(3), linearDerivatives=true)
        @assert(m > 0.0)
        @assert(size(I) == (3, 3))
        @assert(isposdef(I))
        @assert(length(A) == 3)
        @assert(length(freqHz) == 3)
        @assert(length(phase) == 3)
        @assert(minimum(A) > 0.0)
        @assert(length(Q0) == 4)
        @assert(length(w0) == 3)
        is_constraint    = fill(false,7)
        is_constraint[7] = true
        simulationState = ModiaMath.SimulationState("FreeBodyRotation", getModelResidues!, Vector{Float64}([Q0;w0]), getVariableName;
                                                    is_constraint = is_constraint,
                                                    structureOfDAE = linearDerivatives ? ModiaMath.DAE_LinearDerivativesAndConstraints :
                                                                                ModiaMath.DAE_NoSpecialStructure)
        new(simulationState, m, SMatrix{3,3,Float64,9}(I), SVector{3,Float64}(A),
                                                           SVector{3,Float64}(freqHz),
                                                           SVector{3,Float64}(phase))
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                               xNames=["Q[1]", "Q[2]", "Q[3]", "Q[4]", "w[1]", "w[2]", "w[3]"])

function getModelResidues!(m::Model, _t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    u = m.A .* sin.(2.0 * pi * m.freqHz * _t + m.phase)

    Q    = SVector{4,Float64}(_x[1:4])
    w    = SVector{3,Float64}(_x[5:7])
    derQ = SVector{4,Float64}(_derx[1:4])
    derw = SVector{3,Float64}(_derx[5:7])

    _r[1:3] = w - 2.0 * ([ Q[4]  Q[3] -Q[2] -Q[1];
                       -Q[3]  Q[4]  Q[1] -Q[2];
                        Q[2] -Q[1]  Q[4] -Q[3]] * derQ)
    _r[4:6] = m.I * derw + cross(w, m.I * w) - u
    _r[7]   = dot(Q, Q) - 1.0

    return nothing
end

end


text:
module PlutoDataTable

using JSON2
using Random
using Tables

export data_table

"""
    data_table(table; items_per_page = 10)

Display `table` as an interactive data table. `table` is any tabular data structure supporting the Tables.jl interface. The
`items_per_page` keyword argument determines the number of rows to display at one time (default is 10).
"""

function data_table(table; items_per_page = 10)
	app_id = randstring('a':'z')
	d = Dict(
        "headers" => [Dict("text" => string(name), "value" => string(name)) for name in Tables.columnnames(table)],
        "data" => [Dict(string(name) => string(row[name]) for name in Tables.columnnames(table)) for row in Tables.rows(table)]
    )
	djson = JSON2.write(d)

	return HTML("""
		<link href="https://cdn.jsdelivr.net/npm/@mdi/font@5.4.55/css/materialdesignicons.min.css" rel="stylesheet">
		<link href="https://cdn.jsdelivr.net/npm/vuetify@2.3.14/dist/vuetify.min.css" rel="stylesheet">

	  <div id="$app_id">
		<v-app>
		  <v-data-table
		  :headers="headers"
          :items="data"
          :items-per-page = $items_per_page
		></v-data-table>
		</v-app>
	  </div>

	  <script src="https://cdn.jsdelivr.net/npm/vue@2.6.12/dist/vue.js"></script>
	  <script src="https://cdn.jsdelivr.net/npm/vuetify@2.3.14/dist/vuetify.js"></script>

	<script>
		new Vue({
		  el: '#$app_id',
		  vuetify: new Vuetify(),
		  data () {
				return $djson
			}
		})
	</script>
	<style>
		.v-application--wrap {
			min-height: 10vh;
		}
		.v-data-footer__select {
			display: none;
		}
	</style>
	""")
end

end


text:
export STANDARD_TUNING, Annotation, Fretboard, get_tuning, get_frets, annotate_matches!

"""
Standard tuning for 6-string guitar.
"""
const STANDARD_TUNING = [pitch"E4", pitch"B3", pitch"G3", pitch"D3", pitch"A2", pitch"E2"]

####
#### character-based fretboard drawing
####

struct Annotation
    "The string which appears when the fretboard is rendered. Preferably something short."
    mark::String
    "Style for the mark (when supported by the terminal)."
    crayon::Crayon
end

"""
    Annotation(mark::AbstractString, [crayon])

A string annotation, with an optional color or style using [Crayons.Crayon](@ref).
"""
Annotation(mark::AbstractString) = Annotation(String(mark), Crayon())

struct Fretboard <: AbstractMatrix{Vector{Annotation}}
    "Tuning of the strings (from highest pitch)"
    tuning::Vector{Pitch}
    "Number of frets"
    frets::Int
    "Canvas for the annotations"
    canvas::OffsetMatrix{Vector{Annotation}}
    @doc """
    $(SIGNATURES)

    Create a fretboard which can be annotated with strings, either by as an
    `AbstractMatrix`, or using [`annotate_matches!`](@ref) (recommended).
    """
    function Fretboard(tuning::Vector{Pitch}, frets::Int)
        @argcheck frets ≥ 1
        @argcheck !isempty(tuning)
        canvas = OffsetMatrix{Vector{Annotation}}(undef, axes(tuning, 1), 0:frets)
        for i in eachindex(canvas)
            canvas[i] = Annotation[]
        end
        new(tuning, frets, canvas)
    end
end

Base.axes(fretboard::Fretboard) = axes(fretboard.canvas)

function Base.getindex(fretboard::Fretboard, string::Int, fret::Int)
    getindex(fretboard.canvas, string, fret)
end

"""
$(SIGNATURES)

Return the number of frets.
"""
get_frets(fretboard::Fretboard) = fretboard.frets

"""
$(SIGNATURES)

Return the tuning.
"""
get_tuning(fretboard::Fretboard) = get_frets(fretboard.fretboard)

"""
$(SIGNATURES)

For each fret on each string, call `f` with a value that has the following properties:

- `string`: the string number, from highest pitch (1), eg `6` for `pitch"E2"`
- `open_pitch`: the pitch of the open string (eg `pitch"E2"`)
- `pitch`: pitch of the fret (eg `pitch"G2")
- `fret`: fret number, eg `3`.

`f` should return an [`Annotation`](@ref) or `nothing`.
"""
function annotate_matches!(f, fretboard::Fretboard; skip_duplicates = true)
    @unpack tuning, frets, canvas = fretboard
    for (string, open_pitch) in pairs(tuning)
        for fret in 0:frets
            annotation = f((string = string, open_pitch = open_pitch,
                            pitch = open_pitch + fret, fret = fret))
            if !isnothing(annotation)
                annotations = canvas[string, fret]
                if !skip_duplicates || !(annotation ∈ annotations)
                    push!(annotations, annotation)
                end
            end
        end
    end
    fretboard
end

"""
$(SIGNATURES)

Shorthand for inserting the same annotation when `predicate` returns `true`. See the other
method of `annotate_matches!` for details.

# Example

```julia
julia> using Fretboards, Fretboards.Intervals

julia> fb = Fretboard(STANDARD_TUNING, 7);

julia> let root = pitch"C"
           for (mark, semitones) in [("1", 0), ("3", M3), ("5", P5)]
               annotate_matches!(x -> x.pitch ≂ root + semitones, Annotation(mark), fb)
           end
       end

julia> fb
          0      1      2      3      4      5      6      7
E4       3|------|------|-----5|------|------|------|------|
B3        |-----1|------|------|------|-----3|------|------|
G3       5|------|------|------|------|-----1|------|------|
D3        |------|-----3|------|------|-----5|------|------|
A2        |------|------|-----1|------|------|------|-----3|
E2       3|------|------|-----5|------|------|------|------|
```
"""
function annotate_matches!(predicate, annotation::Annotation, fretboard::Fretboard; kwargs...)
    annotate_matches!(x -> predicate(x) ? annotation : nothing, fretboard; kwargs...)
end

"""
$(TYPEDEF)

Characters for drawing a [`Fretboard`](@ref).

# Fields

$(FIELDS)
"""
struct DrawingChars
    "Strings between frets."
    string_char::Char
    "Fret-string crossings."
    fret_char::Char
    "Nut (fret 0)."
    nut_char::Char
end

const UNICODE_CHARS = DrawingChars('─','┼','├')

const ASCII_CHARS = DrawingChars('-','|','|')

"""
$(SIGNATURES)

Text width of all the annotations for a fret, with `string_char`.
"""
function _textwidth_ascii_fret_annotations(string_char, annotations)
    isempty(annotations) && return 0
    annotations_width = sum(annotation -> textwidth(annotation.mark), annotations)
    length(annotations) * textwidth(string_char) + annotations_width
end

"""
$(SIGNATURES)

Print ASCII for a single fret annotation, for the space between frets.

# Arguments

- `io`: output
- `string_char`: character for depicting strings, eg `-`,
- `fret_spacing`: the number of characters between frets
- `annotations`: a vector of `Annotation`s
"""
function _print_ascii_fret_annotations(io::IO, string_char, fret_spacing, annotations)
    padding = fret_spacing - _textwidth_ascii_fret_annotations(string_char, annotations)
    @argcheck padding ≥ 0
    print(io, string_char^padding)
    for annotation in annotations
        print(io, string_char, annotation.crayon(annotation.mark))
    end
end

"""
$(SIGNATURES)

Print ASCII for all annotations for a string.

# Arguments

- `io`: output
- `drawing_chars`: character set for drawing
- `fret_spacing`: the number of characters between frets
- `open_pitch`: pitch of the open string (fret 0)
- `all_annotations`: a vector of `Vector{Annotation}`s
"""
function _print_ascii_all_annotations(io::IO, drawing_chars::DrawingChars, fret_spacing,
                                      open_pitch, all_annotations)
    @unpack string_char, fret_char, nut_char = drawing_chars
    print(io, rpad(string(open_pitch), 4, ' '))
    for (fret, annotations) in pairs(all_annotations)
        _print_ascii_fret_annotations(io, fret == 0 ? ' ' : string_char, fret_spacing,
                                      annotations)
        print(io, fret == 0 ? nut_char : fret_char)
    end
    println(io)
end

"""
$(SIGNATURES)

Print fret numbers.

# Arguments

- `io`: output
- `fret_spacing`: the number of characters between frets
- `frets`: the number of frets
"""
function _print_ascii_fret_numbers(io::IO, fret_spacing, frets)
    print(io, ' '^4)
    for fret in 0:frets
        print(io, lpad(string(fret), fret_spacing + 1, ' '))
    end
    println(io)
end

"""
$(SIGNATURES)

Print `fretboard` with the specified `drawing_chars` character set.
"""
function print_ascii(io::IO, fretboard::Fretboard,
                     drawing_chars::DrawingChars = ASCII_CHARS,
                     min_fret_spacing = 4)
    @unpack frets, tuning, canvas = fretboard
    @unpack string_char = drawing_chars
    fret_spacing = max(maximum(a -> _textwidth_ascii_fret_annotations(string_char, a), canvas),
                       min_fret_spacing)
    _print_ascii_fret_numbers(io, fret_spacing, frets)
    for (string, open_pitch) in pairs(tuning)
        _print_ascii_all_annotations(io, drawing_chars, fret_spacing, open_pitch,
                                     @view canvas[string, :])
    end
end

function Base.show(io::IO, fretboard::Fretboard)
    print_ascii(io, fretboard)
end

function Base.show(io::IO, ::MIME"text/plain", fretboard::Fretboard)
    print_ascii(io, fretboard)
end


text:
# Author: Lukasz Bratos

include("methods.jl")

f1(x) = exp(1 - x) - 1
pf1(x) = -exp(1 - x)
f2(x) = x * exp(-x)
pf2(x) = -exp(-x) * (x - 1)
delta = 10^-5
epsilon = 10^-5
maxit = 100

# f1
println("Metoda bisekcji dla f1 na przedziale [0.0, 1.5]: ", Methods.mbisekcji(f1, 0.0, 1.5, delta, epsilon))
println("Metoda bisekcji dla f1 na przedziale [0.0, 2.0]: ", Methods.mbisekcji(f1, 0.0, 2.0, delta, epsilon))
println("Metoda Newtona dla f1 z x0 = 0.0: ", Methods.mstycznych(f1, pf1, 0.0, delta, epsilon, maxit))
println("Metoda siecznych dla f1, gdzie x0 = 0.5, a x1 = 1.5: ", Methods.msiecznych(f1, 0.0, 2.0, delta, epsilon, maxit))

# f2
println("Metoda bisekcji dla f2 na przedziale [-1.0, 1.5]: ", Methods.mbisekcji(f2, -1.0, 1.5, delta, epsilon))
println("Metoda bisekcji dla f2 na przedziale [-1.0, 1.0]: ", Methods.mbisekcji(f2, -1.0, 1.0, delta, epsilon))
println("Metoda Newtona dla f2 z x0 = -1.0: ", Methods.mstycznych(f2, pf2, -1.0, delta, epsilon, maxit))
println("Metoda siecznych dla f2, gdzie x0 = -1.0, a x1 = 1.0: ", Methods.msiecznych(f2, -1.0, 1.0, delta, epsilon, maxit))

# testy Newton
println("Metoda Newtona dla f1 z x0 = 3.0: ", Methods.mstycznych(f1, pf1, 3.0, delta, epsilon, maxit))
println("Metoda Newtona dla f1 z x0 = 15.0: ", Methods.mstycznych(f1, pf1, 15.0, delta, epsilon, maxit))
println("Metoda Newtona dla f2 z x0 = 1.0: ", Methods.mstycznych(f2, pf2, 1.0, delta, epsilon, maxit))
println("Metoda Newtona dla f2 z x0 = 50.0: ", Methods.mstycznych(f2, pf2, 50.0, delta, epsilon, maxit))


text:
module BioParticleSimulation

using Random, StatsBase, Plots

include("genome.jl")

export FieldSet, fillfield, plotfield, addresource, round
    abstract type AbstractField end

    const UP    = ( 1,  0)
    const DOWN  = (-1,  0)
    const LEFT  = ( 0, -1)
    const RIGHT = ( 0,  1)

    mutable struct FieldSet
        field::Matrix{Union{AbstractField, Bool}}
        list::Vector{AbstractField}
        rfield::Matrix{Float64}
        makecell::Function
        function FieldSet(x::Int, y::Int)::FieldSet
            new(Matrix{Union{AbstractField, Bool}}(undef, y, x), Vector{AbstractField}(undef, 0), zeros(Int, y, x), makecell_simple)::FieldSet
        end
    end

    function Base.length(f::FieldSet)
        length(f.list)
    end

struct CellType
    id::Int
    movep::Float32
    killp::Float32
    replp::Float32
    kille::Int
    maxe::Int
    maxl::Int
    replcost::Int
    replinite::Int
    movecost::Int
    idlecost::Int
    function CellType(
        id,
        movep,
        killp,
        replp,
        kille,
        maxe,
        maxl,
        replcost,
        replinite,
        movecost,
        idlecost,
    )::CellType
        new(
            id,
            movep,
            killp,
            replp,
            kille,
            maxe,
            maxl,
            replcost,
            replinite,
            movecost,
            idlecost,
        )::CellType
    end
    function CellType(;
        id = 0,
        movep = 0.3,
        killp = 0.000001,
        replp = 0.0007,
        kille = 100,
        maxe = 2000,
        maxl = 4000,
        replcost = 100,
        replinite = 100,
        movecost = 10,
        idlecost = 1,
    )::CellType
        new(
            id,
            movep,
            killp,
            replp,
            kille,
            maxe,
            maxl,
            replcost,
            replinite,
            movecost,
            idlecost,
        )::CellType
    end
end

    const CELL_A = CellType(0, 0.3, 0.000001, 0.0007, 100, 2000, 4000, 100, 100, 10, 1)

    mutable struct Cell <: AbstractField
        enrg::Float64
        x::Int
        y::Int
        livetime::Int
        field::FieldSet
        type::CellType
        function Cell(enrg, x, y, l, field)::Cell
            new(enrg, x, y, l, field, CELL_A)::Cell
        end
        function Cell(enrg, x, y, l, field, type)::Cell
            new(enrg, x, y, l, field, type)::Cell
        end
    end

    function fillfield(f::FieldSet, p, type)
        x    = size(f.field, 2)
        y    = size(f.field, 1)
        mask = sample([true, false], Weights([p, 1-p]), (x,y))
        for r = 1:y
            for c = 1:x
                if mask[r, c]
                    newcell = Cell(1000, c, r, type.maxl, f, type)
                    push!(f.list, newcell)
                    f.field[r, c] = newcell
                else
                    f.field[r, c] = false
                end
            end
        end
    end
    function fillfield(f::FieldSet, p)
        fillfield(f::FieldSet, p, CELL_A)
    end
"""
    fillfield(f::FieldSet)

p - cell contamination probability;
type -
"""
    function fillfield(f::FieldSet; p, type)
        fillfield(f::FieldSet, p, type)
    end

    function addresource(f::FieldSet, n, val)
        x     = size(f.rfield, 2)
        y     = size(f.rfield, 1)
        for i = 1:n
            rx = rand(2:x-1)
            ry = rand(2:y-1)
            if f.rfield[ry, rx] < 2000 f.rfield[ry, rx] += val end
        end
    end
    function growresource(f::FieldSet, m)
        for r = 2:size(f.rfield, 1)-1
            for c = 2:size(f.rfield, 2)-1
                f.rfield[r, c] += f.rfield[r, c] * m
                if f.rfield[r, c] > 2000 f.rfield[r, c] = 2000 end
                if f.rfield[r, c] > 10
                    if f.rfield[r, c] > f.rfield[r, c + 1] f.rfield[r, c] -=1;  f.rfield[r, c + 1] +=1 end
                    if f.rfield[r, c] > f.rfield[r + 1, c] f.rfield[r, c] -=1;  f.rfield[r + 1, c] +=1 end
                    if f.rfield[r, c] > f.rfield[r, c - 1] f.rfield[r, c] -=1;  f.rfield[r, c - 1] +=1 end
                    if f.rfield[r, c] > f.rfield[r - 1, c] f.rfield[r, c] -=1;  f.rfield[r - 1, c] +=1 end
                end
            end
        end
    end

    function cellonfield(f, x, y)
        if isa(f.field[y, x], Cell)
            if f.field[y, x].enrg > 0 && f.field[y, x].livetime > 0 return true end
        end
        false
    end
    #=
    function resonfield(f, x, y)
        isa(f.field[y, x], Resource)
    end
    =#
    function emptycellonfield(f, x, y)
        if isa(f.field[y, x], Bool)
            return true
        elseif isa(f.field[y, x], Cell)
            if f.field[y, x].enrg <= 0 || f.field[y, x].livetime <= 0 return true end
        end
        false
    end

    function smovecell(c, x, y)
        c.enrg -= c.type.movecost
        setcellfield(c, false)
        c.x = x
        c.y = y
        setcellfield(c, c)
    end
    function movecell(c)
        vec = sample([UP, DOWN, LEFT, RIGHT])
        newx = c.x + vec[1]
        newy = c.y + vec[2]
        if c.enrg > c.type.movecost && c.livetime > 0 && newx > 0 && newy > 0 && newx <= size(c.field.field, 2) && newy <= size(c.field.field, 1)
            if cellonfield(c.field, newx, newy)
                killcell(c, newx, newy)
                return
            end
            #c.enrg -= c.type.movecost
            smovecell(c, newx, newy)
            if c.enrg + getrfield(c, newx, newy) >= c.type.maxe
                setrfield(c, newx, newy, getrfield(c, newx, newy) - c.type.maxe + c.enrg)
                c.enrg = c.type.maxe
            else
                c.enrg += getrfield(c, newx, newy)
                setrfield(c, newx, newy, 0)
            end

        else
            idlecell(c)
        end
    end
    function replcell(c)
        if c.enrg > c.type.replcost && c.livetime > 0
            c.field.makecell(c)
        else
            idlecell(c)
        end
    end
    function makecell_simple(cp::Cell)
        vec = sample([UP, DOWN, LEFT, RIGHT])
        newx = cp.x + vec[1]
        newy = cp.y + vec[2]
        if newx > 0 && newy > 0 && newx <= size(cp.field.field, 2) && newy <= size(cp.field.field, 1) && emptycellonfield(cp.field, newx, newy)
            cp.enrg -= cp.type.replcost
            newcell = Cell(cp.type.replinite, newx, newy, cp.type.maxl, cp.field, cp.type)
            push!(cp.field.list, newcell)
            cp.field.field[newy, newx] = newcell
        else
            idlecell(cp)
        end
    end

    function killcell(c, x, y)
        if sample([true, false], Weights([c.type.killp, 1 - c.type.killp]))
            kc = getfield(c, x, y)
            if c.enrg + kc.type.kille >= c.type.maxe
                c.enrg = c.type.maxe
            else
                c.enrg += kc.type.kille
            end
            #println("KILL!")
            setcellfield(kc, false)
            kc.enrg     = 0.0
            kc.livetime = 0.0

            smovecell(c, x, y)
        else
            idlecell(c)
        end
    end
    function idlecell(c)
        c.enrg -= c.type.idlecost
        if c.enrg <= 0 || c.livetime <= 0
            setcellfield(c, false)
        end
    end

    function setcellfield(c::Cell, val)
        c.field.field[c.y, c.x] = val
    end
    function setfield(c::Cell, x, y, val)
        c.field.field[y, x] = val
    end
    function getfield(c::Cell, x, y)
        c.field.field[y, x]
    end
    function getrfield(c::Cell, x, y)
        c.field.rfield[y, x]
    end
    function setrfield(c::Cell, x, y, val)
        c.field.rfield[y, x] = val
    end

    function cellaction(c::Cell)
        if c.enrg > 0 && c.livetime > 0
            action =  sample([movecell, replcell, idlecell], Weights([c.type.movep, c.type.replp, 1 - c.type.movep - c.type.replp]))
            action(c)
            c.livetime -= 1
            #if c.livetime <= 0
            #    setcellfield(c, false)
            #end
        end
        if c.enrg <= 0 || c.livetime <= 0
            setcellfield(c, false)
        end
    end

    #function cellaction(b::Bool)
    #end

    function round(f::FieldSet)
        for c in f.list
            cellaction(c)
        end
        deleteat!(f.list, findall(x -> x.enrg <= 0.0 || x.livetime <= 0.0, f.list))
        nothing
    end

    function round(f::FieldSet, n::Int)
        for i = 1:n
            round(f)
        end
    end

    function plotfield(f)
        cnt  = 0
        mx   = Matrix{Int}(undef, length(f.list), 2)
        ccnt  = 1

        for r = 1:size(f.field, 1)
            for c = 1:size(f.field, 2)
                if isa(f.field[r, c], Cell)
                    mx[ccnt, 1] = r
                    mx[ccnt, 2] = c
                    ccnt += 1
                end
            end
        end
        mygrad = cgrad([:white, :green])
        p = heatmap(f.rfield, c = mygrad, clim = (0, 500), legend = false)
        #p = heatmap(f.rfield, c = :greens, clim = (0, 2000), legend = false)
        p = scatter!(mx[:,2], mx[:,1], title = "Cells", color = :red, legend = false)
        #p = scatter!(p, mxr[:,1], mxr[:,2], title = "Cells", color = :green, legend = false)
        p
    end
    #=
    function Base.show(io::IO, f::FieldSet)
        x    = size(f.field, 2)
        y    = size(f.field, 1)
        for r = 1:y
            for c = 1:x
                if f.field[r, c] == false
                    print(io, " * ")
                elseif isa(f.field[r, c], Cell)
                    print(io, " C ")
                else
                    print(io, "   ")
                end
            end
            println(io, "")
        end
    end
    =#
    function printall(f)
        x    = size(f.field, 2)
        y    = size(f.field, 1)
        cnt  = 0
        for r = 1:y
            for c = 1:x
                if isa(f.field[r, c], Cell)
                    println("X: $c ($(f.field[r, c].x)) Y: $r ($(f.field[r, c].y)) Energy: $(f.field[r, c].enrg) Live: $(f.field[r, c].livetime)")
                    cnt += 1
                end
            end
        end
        println("Total: $cnt")
    end

end # module


text:
"""

```
grid3D(solute,solute_atoms,mddf_result,output_file; dmin=1.5, ddax=5.0, step=0.5)
```

This function builds the grid of the 3D density function and fills an array of
mutable structures of type Atom, containing the position of the atoms of
grid, the closest atom to that position, and distance.

`solute` is a `ComplexMixtuers.Selection`, defining the solute. `solute_atoms` is the corresponding
vector of `PDBTools.Atom`s, and `mddf_result` is the result of a `mddf_result` calculation with
the correspondign solute.

`dmin` and `dmax` define the range of distance where the density grid will be built, and `step`
defines how fine the grid must be. Be aware that fine grids involve usually a very large (hundreds
of thousands points).

All parameters can be provides as keyword parameters.

### Example

```julia-repl
julia> using ComplexMixtures, PDBTools

julia> pdb = readPDB("./system.pdb");

julia> R = ComplexMixtures.load("./results.json");

julia> protein = select(pdb,"protein");

julia> solute = ComplexMixtures.Selection(protein,nmols=1);

julia> grid = ComplexMixtures.grid3D(solute=solute, solute_atoms=protein, mddf_result=R, output_file="grid.pdb");

```

`grid` will contain a vector of `Atom`s with the information of the MDDF at each grid point, and the
same data will be written in the `grid.pdb` file. This PDB file can be opened in VMD, for example, and contain
in the `beta` field the contribution of each protein residue to the MDDF at each point in space relative to the
protein, and in the `occupancy` field the distance to the protein. Examples of how this information can be
visualized are provided in the user guide of `ComplexMixtures`.

"""
grid3D(;solute=nothing, solute_atoms=nothing, mddf_result=nothing, output_file=nothing,
        dmin=1.5, dmax=5.0, step=0.5) =
  grid3D(solute, solute_atoms, mddf_result, output_file, dmin=dmin, dmax=dmax, step=step)

function grid3D(solute,solute_atoms,mddf_result,output_file;
                dmin=1.5, dmax=5.0, step=0.5)

  if nothing in (solute, solute_atoms, mddf_result)
    error("grid3D requires `solute`, `solute_atoms` and `mddf_result` definitions.")
  end

  # Simple function to interpolate data
  interpolate(x₁,x₂,y₁,y₂,xₙ) = y₁ + (y₂-y₁)/(x₂-x₁)*(xₙ-x₁)

  # Maximum and minimum coordinates of the solute
  lims = maxmin(solute_atoms)
  n = @. ceil(Int,(lims.xlength + 2*dmax)/step + 1)

  # Building the grid with the nearest solute atom information
  igrid = 0
  grid = PDBTools.Atom[]
  for ix in 1:n[1]; x = lims.xmin[1]-dmax + step*(ix-1)
  for iy in 1:n[2]; y = lims.xmin[2]-dmax + step*(iy-1)
  for iz in 1:n[3]; z = lims.xmin[3]-dmax + step*(iz-1)
    rgrid = -1
    iat, r = PDBTools.closest(x,y,z,solute_atoms)
    if (dmin < r < dmax)
      if rgrid < 0 || r < rgrid
        at = solute_atoms[iat]
        # Get contribution of this atom to the MDDF
        c = ComplexMixtures.contrib(solute,mddf_result.solute_atom,[at.index])
        # Interpolate c at the current distance
        iright = findfirst(d -> d > r, mddf_result.d)
        ileft = iright - 1
        cᵣ = interpolate(mddf_result.d[ileft],mddf_result.d[iright],c[ileft],c[iright],r)
        if cᵣ > 0
          gridpoint = Atom(index=at.index,index_pdb=at.index_pdb,
                           name=at.name,chain=at.chain,
                           resname=at.resname,resnum=at.resnum,
                           x=x,y=y,z=z,
                           occup=r, beta=cᵣ,
                           model=at.model,segname=at.segname,)
          if rgrid < 0
            igrid += 1
            push!(grid,gridpoint)
          elseif r < rgrid
            grid[igrid] = gridpoint
          end
          rgrid = r
        end # cᵣ>0
      end # rgrid
    end # dmin/dmax
  end #iz
  end #iy
  end #ix

  # Now will scale the density to be between 0 and 99.9 in the temperature
  # factor column, such that visualization is good enough
  bmin, bmax = +Inf, -Inf
  for gridpoint in grid
    bmin = min(bmin,gridpoint.beta)
    bmax = max(bmax,gridpoint.beta)
  end
  for gridpoint in grid
    gridpoint.beta = (gridpoint.beta - bmin)/(bmax-bmin)
  end

  if output_file != nothing
    writePDB(grid,output_file)
  end
  return grid
end


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

"""
    BlockCrossValidation(sides; loss=Dict())

Cross-validation with blocks of given `sides`. Optionally,
specify `loss` function from `LossFunctions.jl` for some
of the variables. If only one side is provided, then blocks
become cubes.

## References

* Roberts et al. 2017. [Cross-validation strategies for data with
  temporal, spatial, hierarchical, or phylogenetic structure]
  (https://onlinelibrary.wiley.com/doi/10.1111/ecog.02881)
* Pohjankukka et al. 2017. [Estimating the prediction performance
  of spatial models via spatial k-fold cross-validation]
  (https://www.tandfonline.com/doi/full/10.1080/13658816.2017.1346255)
"""
struct BlockCrossValidation{S} <: ErrorEstimationMethod
  sides::S
  loss::Dict{Symbol,SupervisedLoss}
end

BlockCrossValidation(sides; loss=Dict()) =
  BlockCrossValidation{typeof(sides)}(sides, loss)

function error(solver, problem, method::BlockCrossValidation)
  s = method.sides
  n = ncoords(_foldable(problem))
  l = length(s) > 1 ? s : ntuple(i->s, n)

  # uniform weights
  weighting = UniformWeighting()

  # block folds
  folding = BlockFolding(l)

  wcv = WeightedCrossValidation(weighting, folding,
                                lambda=1, loss=method.loss)

  error(solver, problem, wcv)
end


text:
# main.jl

include("./infrastructure.jl")

# Test co
struct AppliPerson
   id::String
   name::String
   email::String
    #constructors
   AppliPerson(id, name) = new(id, name, "")
   AppliPerson(id, name, email) = new(id, name, email)
end

daisy = AppliPerson("A", "Daisy", "")

mickey = AppliPerson("B", "Mickey", "")

#db = connect("test.sqlite")
db = DBInterface.connect(SQLite.DB, "test.sqlite")

archive(db, "subscribers", [daisy, mickey])

#import SQLite.DBInterface

r1 = retrieve( db, "subscribers")
#r1 = SQLite.DBInterface.execute(db, "select * from subscribers") |> DataFrame

println(r1)

r2 = retrieve( db, "subscribers", "key = '" * daisy.id * "'" )
#r2 = SQLite.DBInterface.execute( db, "select * from subscribers where key = '$(daisy.id)'") |> DataFrame

println(r2)

store(db, "users", [daisy, mickey])

#r3 = retrieve( db, "users")
r3 = SQLite.DBInterface.execute( db, "select * from users") |> DataFrame

println(r3)

#r4 = retrieve( db, "users", "id = '" * daisy.id * "'" )
r4 = SQLite.DBInterface.execute(db, "select * from users where id = '$(daisy.id)'") |> DataFrame

println(r4)

cmd = `rm test.sqlite`
run(cmd)


text:
function indegrees(receivers, bins)
    ones = (similar(receivers) .= 1)
    scatter(ones, receivers, bins, +)
end

function indegrees(g)
    indegrees(receivers(g), nnodes(g))
end

function symmetricnorm(g, indegreesv)
    d = max.(indegreesv, 1)

    dᵢ = gather(d, senders(g))
    dⱼ = gather(d, receivers(g))

    # TODO: actual solution for type stability!
    Float32.(sqrt.(1 ./ dᵢ) .* sqrt.(1 ./ dⱼ))
end

function symmetricnorm(g)
    indegreesv = indegrees(receivers(g), nnodes(g))
    symmetricnorm(g, indegreesv)
end

function nodespergraph(g::BatchGraphTuple)
    ones = similar(g.nodes, Int32, size(g.nodes, 2))
    ones .= 1
    scatter(ones, g.node2graph, g.numgraphs, +)
end

Flux.@nograd indegrees
Flux.@nograd symmetricnorm
Flux.@nograd nodespergraph

struct GCN
    l
    σ
end

function GCN(in::Integer, out::Integer, σ=identity)
    return GCN(Dense(in, out, bias=false), σ)
end

function (l::GCN)(g::AbstractGraphTuple; symnorm::Bool=true)
    gathered = gather(g.nodes, g.senders)
    if symnorm gathered = symmetricnorm(g) .* gathered end
    @set g.nodes = l.σ.(l.l(scatter(gathered, g.receivers, size(g.nodes, 2), +)))
end

Flux.@functor GCN

function graphmeanpool(g::BatchGraphTuple)
    nodesperg = nodespergraph(g)
    scatter(g.nodes, g.node2graph, g.numgraphs, +) ./ Flux.unsqueeze(nodesperg, 1)
end

struct NodeLayer
    layer
end

function (l::NodeLayer)(g)
	@set g.nodes = l.layer(g.nodes)
end

Flux.@functor NodeLayer


text:
import JSON

abstract type Filter{T,TENC} end
function getfilters(d::Dict)
    if !haskey(d,"filters")
        return nothing
    else
        if d["filters"] === nothing || isempty(d["filters"])
            return nothing
        end
        f = map(d["filters"]) do f
            getfilter(filterdict[f["id"]], f)
        end
        return (f...,)
    end
end
sourcetype(::Filter{T}) where T = T
desttype(::Filter{<:Any,T}) where T = T

zencode(ain,::Nothing) = ain

"""
    VLenArrayFilter(T)

Encodes and decodes variable-length arrays of arbitrary data type
"""
struct VLenArrayFilter{T} <: Filter{T,UInt8} end

function zdecode(ain, ::VLenArrayFilter{T}) where T
    f = IOBuffer(ain)
    nitems = read(f, UInt32)
    out = Array{Vector{T}}(undef,nitems)
    for i=1:nitems
        len1 = read(f,UInt32)
        out[i] = read!(f,Array{T}(undef,len1 ÷ sizeof(T)))
    end
    close(f)
    out
end

#Encodes Array of Vectors a into bytes
function zencode(ain,::VLenArrayFilter)
    b = IOBuffer()
    nitems = length(ain)
    write(b,UInt32(nitems))
    for a in ain
        write(b, UInt32(length(a) * sizeof(eltype(a))))
        write(b, a)
    end
    take!(b)
end

JSON.lower(::VLenArrayFilter{T}) where T = Dict("id"=>"vlen-array","dtype"=> typestr(T) )

getfilter(::Type{<:VLenArrayFilter}, f) = VLenArrayFilter{typestr(f["dtype"])}()

filterdict = Dict("vlen-array"=>VLenArrayFilter)


text:
using PDBTools
using ComplexMixtures
const CM = ComplexMixtures

# Load PDB file of the system
atoms = readPDB("../Data/system.pdb")

# Select the protein and the GLYC molecules
protein = select(atoms,"protein")
glyc = select(atoms,"resname GLYC")

# Setup solute and solvent structures
solute = CM.Selection(protein,nmols=1)
solvent = CM.Selection(glyc,natomspermol=14)

# Setup the Trajectory structure
trajectory = CM.Trajectory("../Data/glyc50.dcd",solute,solvent)

# Run the calculation and get results
results = CM.mddf(trajectory)

# Save the reults to recover them later if required
CM.save(results,"./glyc50.json")

# We will now load the "true" results file, computed from the
# full trajectory.
results = CM.load("../Data/results_glyc50.json")

# Produce plots
using Plots, Plots.PlotMeasures, LaTeXStrings
default(fontfamily="Computer Modern",
        linewidth=2, framestyle=:box, label=nothing, grid=false)

# The complete MDDF and the Kirkwood-Buff Integral
plot(layout=(1,2))
plot!(results.d,results.mddf,
      xlabel=L"r/\AA",ylabel="mddf",subplot=1)
hline!([1],linestyle=:dash,linecolor=:gray,subplot=1)
plot!(results.d,results.kb/1000, #to L/mol
      xlabel=L"r/\AA",ylabel=L"G_{us}/\mathrm{L~mol^{-1}}",subplot=2)
plot!(size=(800,300),margin=4mm)
savefig("./mddf.pdf")

# Atomic contributions to the MDDF
hydroxyls = ["O1","O2","O3","H1","H2","H3"]
aliphatic = ["C1","C2","HA","HB","HC","HD"]
hydr_contrib = CM.contrib(solvent,results.solvent_atom,hydroxyls)
aliph_contrib = CM.contrib(solvent,results.solvent_atom,aliphatic)

plot(results.d,results.mddf,xlabel=L"r/\AA",ylabel="mddf",size=(600,400))
plot!(results.d,hydr_contrib,label="Hydroxyls")
plot!(results.d,aliph_contrib,label="Aliphatic chain")
hline!([1],linestyle=:dash,linecolor=:gray)
savefig("./mddf_atom_contrib.pdf")


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    int xTESTOUTPUT = 7;
    int wTESTOUTPUT = 7;
    void m() {
      final Condition c = new Condition();
      int y = ({L /c H}) 42;
      int z = 0;
      if (y > 0) {
	  //c = true;
	  z = 4; // level of z should be {L c/ H}, since the PC should hasn't been raised since c hasn't been set.
      }
      this.xTESTOUTPUT = y; // At this point, the level of y should be {L /c H}, since c hasn't been set.
      this.wTESTOUTPUT = z;
    }
    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
#
# Coherence (as the square root of the correlation matrix)
#

immutable Correlation <: PairwiseStatistic; end
Base.eltype{T<:Real}(::Correlation, X::AbstractArray{T}) = T

# Single input matrix
allocwork{T<:Real}(::Correlation, X::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), Array(T, 1, size(X, 2)))
function computestat!{T<:Real}(::Correlation, out::AbstractMatrix{T}, work::Tuple{Matrix{T},Matrix{T}},
                               X::AbstractVecOrMat{T})
    Xmμ, μ = work
    mean!(μ, X)
    broadcast!(-, Xmμ, X, μ)
    cov2coh!(out, Ac_mul_A!(out, Xmμ))
end

# Two input matrices
allocwork{T<:Real}(::Correlation, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), Array(T, size(Y, 1), size(Y, 2)), cov2coh_work(X), cov2coh_work(Y))
function computestat!{T<:Real}(::Correlation, out::AbstractMatrix{T},
                               work::Tuple{Matrix{T}, Matrix{T}, Array{T}, Array{T}},
                               X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    Xmμ, Ymμ, Xμ, Yμ = work
    mean!(Xμ, X)
    broadcast!(-, Xmμ, X, Xμ)
    mean!(Yμ, Y)
    broadcast!(-, Ymμ, Y, Yμ)
    cov2coh!(out, Xmμ, Ymμ, Xμ, Yμ, Ac_mul_B!(out, Xmμ, Ymμ))
end

#
# Jackknifing for Correlation
#
surrogateval(::Correlation, v) = v

# Single input matrix
allocwork{T<:Real}(t::JackknifeSurrogates{Correlation},
                   X::AbstractVecOrMat{T}) = (Array(T, size(X, 1), size(X, 2)),
                                              Array(T, 1, size(X, 2)),
                                              Array(T, div(size(X, 1), jnn(t)), size(X, 2)))
function computestat!{T<:Real}(t::JackknifeSurrogates{Correlation},
                               out::JackknifeSurrogatesOutput,
                               work::NTuple{3,Matrix{T}},
                               X::AbstractVecOrMat{T})
    stat = t.transform
    trueval = out.trueval
    surrogates = out.surrogates
    Xmμ, μ, jninvmoment = work
    chkinput(trueval, X)
    ntrials(X) % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    size(out.surrogates, 1) == div(ntrials(X), jnn(t)) || throw(DimensionMismatch("invalid output size"))

    # Compute means and comoments
    mean!(μ, X)
    broadcast!(-, Xmμ, X, μ)
    Ac_mul_A!(trueval, Xmμ)

    # Compute jackknifed moments
    multiplier = ntrials(X)/(ntrials(X) - 1)
    @inbounds for j = 1:size(X, 2)
        ssq = trueval[j, j]
        @simd for i = 1:size(surrogates, 1)
            v = ssq - abs2(Xmμ[i, j])*multiplier
            jninvmoment[i, j] = sqrt(1/ifelse(v < 0, zero(v), v))
        end
    end

    # Surrogates
    @inbounds for k = 1:size(X, 2)
        for j = 1:k-1
            comoment = trueval[j, k]
            @simd for i = 1:size(surrogates, 1)
                surrogates[i, j, k] = (comoment - Xmμ[i, j]*Xmμ[i, k]*multiplier)*(jninvmoment[i, j]*jninvmoment[i, k])
            end
        end
        for i = 1:size(surrogates, 1)
            surrogates[i, k, k] = 1
        end
    end

    # Finish true value
    cov2coh!(trueval, trueval)

    out
end

# Two input matrices
allocwork{T<:Real}(t::JackknifeSurrogates{Correlation}, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), Array(T, size(Y, 1), size(Y, 2)),
     Array(T, 1, size(X, 2)), Array(T, 1, size(Y, 2)),
     Array(T, div(size(X, 1), jnn(t)), size(X, 2)),
     Array(T, div(size(Y, 1), jnn(t)), size(Y, 2)))
function computestat!{T<:Real}(t::JackknifeSurrogates{Correlation},
                               out::JackknifeSurrogatesOutput,
                               work::NTuple{6,Matrix{T}},
                               X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    stat = t.transform
    trueval = out.trueval
    surrogates = out.surrogates
    Xmμ, Ymμ, Xμ, Yμ, Xjninvmoment, Yjninvmoment = work
    chkinput(trueval, X, Y)
    ntrials(X) % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    size(out.surrogates, 1) == div(ntrials(X), jnn(t)) || throw(DimensionMismatch("invalid output size"))

    # Compute means and comoments
    mean!(Xμ, X)
    broadcast!(-, Xmμ, X, Xμ)
    mean!(Yμ, Y)
    broadcast!(-, Ymμ, Y, Yμ)
    Ac_mul_B!(trueval, Xmμ, Ymμ)

    # Compute jackknifed moments
    multiplier = ntrials(X)/(ntrials(X) - 1)
    for (arrmμ, jninvmoment) in ((Xmμ, Xjninvmoment), (Ymμ, Yjninvmoment))
        @inbounds for j = 1:size(arrmμ, 2)
            ssq = zero(T)
            @simd for i = 1:size(arrmμ, 1)
                ssq += abs2(arrmμ[i, j])
            end
            @simd for i = 1:size(surrogates, 1)
                jninvmoment[i, j] = sqrt(1/(ssq - abs2(arrmμ[i, j])*multiplier))
            end
        end
    end

    # Surrogates
    @inbounds for k = 1:size(Y, 2), j = 1:size(X, 2)
        comoment = trueval[j, k]
        @simd for i = 1:size(surrogates, 1)
            surrogates[i, j, k] = (comoment - Ymμ[i, k]*Xmμ[i, j]*multiplier)*(Yjninvmoment[i, k]*Xjninvmoment[i, j])
        end
    end

    # Finish true value
    cov2coh!(trueval, Xmμ, Ymμ, Xμ, Yμ, trueval)

    out
end


text:
using DSP

"""
    rawphase(intensity0, intensity1, intensity2, intensity3)

Compute the raw phase from four intensity of interference pattern. The reference
in the pattern of intensity `intensityi` is assumed to be shifted of i×π/2.

See also: [`unwrapphase`](@ref), [`retrievephase`](@ref)
"""
function rawphase(
    intensity1::Array{<:Number,2},intensity2::Array{<:Number,2},
    intensity3::Array{<:Number,2},intensity4::Array{<:Number,2}
    )
    atan.(intensity4.-intensity2, intensity1.-intensity3)
end

"""
    rawphase(intensity0, intensity1, intensity2)

Compute the raw phase from four intensity of interference pattern. The reference
in the pattern of intensity `intensityi` is assumed to be shifted of i×π/3.

See also: [`unwrapphase`](@ref), [`retrievephase`](@ref)
"""
function rawphase(
    intensity1::Array{<:Number,2},intensity2::Array{<:Number,2},
    intensity3::Array{<:Number,2}
    )
    atan.(√(3) .* (intensity2.-intensity3), 2*intensity1.-intensity3.-intensity2)
end

"""
    unwrapphase(wrappedphase)

Unwrap the phase. The method used is inspired by the Standard Experiment Laboratory
script from the MSc in Optics and Photonics at Imperial College London. The input
is typically the output of [`rawphase`](@ref).

See also: [`rawphase`](@ref), [`retrievephase`](@ref)
"""
function unwrapphase(wrappedphase::Array{<:Number,2})
    unwrapped = copy(wrappedphase)
    unwrapped[end÷2+1:end,:] = unwrap(unwrapped[end÷2+1:end, :], dims=1)
    unwrapped = reverse(unwrapped, dims=1)
    unwrapped[end÷2+1:end,:] = unwrap(unwrapped[end÷2+1:end, :], dims=1)
    unwrapped[:, end÷2+1:end] = unwrap(unwrapped[:, end÷2+1:end], dims=2)
    unwrapped = reverse(unwrapped, dims=2)
    unwrapped[:, end÷2+1:end] = unwrap(unwrapped[:, end÷2+1:end], dims=2)
    reverse(reverse(unwrapped, dims=2), dims=1)
end

"""
    retrievephase(intensity0, intensity1, intensity2[, intensity3]; correct_aberrations, mask)

Basically compose [`rawphase`](@ref) and [`unwrapphase`](@ref). If `correct_aberrations`
is given, tries to [`correct`](@ref) the given aberrations with the given mask.

See also: [`rawphase`](@ref), [`unwrapphase`](@ref), [`correct`](@ref)
"""
function retrievephase(
    intensity1::Array{<:Number,2},intensity2::Array{<:Number,2},
    intensity3::Array{<:Number,2},intensity4::Array{<:Number,2};
    correct_aberrations=[],
    mask=(x,y)->1
    )
    unwrapped = (unwrapphase∘rawphase)(intensity1, intensity2, intensity3, intensity4)
    for aberration in correct_aberrations
        unwrapped = correct(aberration, unwrapped, mask=mask)
    end
    x,y = range(-1,1,length=size(unwrapped)[1]), range(-1,1,length=size(unwrapped)[2])
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    unwrapped .* mask.(X,Y)
end
function retrievephase(
    intensity1::Array{<:Number,2},intensity2::Array{<:Number,2},
    intensity3::Array{<:Number,2};
    correct_aberrations=[],
    mask=(x,y)->1
    )
    unwrapped = (unwrapphase∘rawphase)(intensity1, intensity2, intensity3)
    for aberration in correct_aberrations
        unwrapped = correct(aberration, unwrapped, mask=mask)
    end
    x,y = range(-1,1,length=size(unwrapped)[1]), range(-1,1,length=size(unwrapped)[2])
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    unwrapped .* mask.(X,Y)
end


text:
mutable struct NER
    worddict::Dict
    chardict::Dict
    tagset
    model
end

function NER()
    words = h5read(wordembeds_file, "key")
    worddict = Dict(words[i] => i for i=1:length(words))
    chardict = Dict{String,Int}()
    NER(worddict, chardict, BIO(), nothing)
end

function encode(ner::NER, words::Vector{String})
    worddict = ner.worddict
    chardict = ner.chardict
    unkword = worddict["UNKNOWN"]
    w = map(w -> get!(worddict,lowercase(w),unkword), words)
    cs = map(words) do w
        chars = Vector{Char}(w)
        map(c -> get!(chardict,string(c),length(chardict)+1), chars)
    end
    w, cs
end

function readdata!(ner::NER, path::String)
    data_w, data_c, data_t = Var[], Var[], Var[]
    words, tags = String[], String[]
    lines = open(readlines, path)
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line) || i == length(lines)
            isempty(words) && continue
            w, cs = encode(ner, words)
            t = encode(ner.tagset, tags)
            push!(data_w, Var(w))
            batchdims = map(length, cs)
            c = cat(1, cs...)
            push!(data_c, Var(c,batchdims))
            push!(data_t, Var(t))
            empty!(words)
            empty!(tags)
        else
            items = split(line, "\t")
            push!(words, String(items[1]))
            push!(tags, String(items[2]))
            #word = replace(word, r"[0-9]", '0')
        end
    end
    data_w, data_c, data_t
end

function train(ner::NER, trainfile::String, testfile::String)
    train_w, train_c, train_t = readdata!(ner, trainfile)
    test_w, test_c, test_t = readdata!(ner, testfile)
    info("# Training sentences:\t$(length(train_w))")
    info("# Testing sentences:\t$(length(test_w))")
    info("# Words:\t$(length(ner.worddict))")
    info("# Chars:\t$(length(ner.chardict))")
    info("# Tags:\t$(length(ner.tagset))")

    wordembeds = h5read(wordembeds_file, "value")
    charembeds = randn(Float32, 20, length(ner.chardict)) * 0.01f0
    tagembeds = randn(Float32, 50, length(ner.tagset)) * 0.01f0
    ner.model = Model(wordembeds, charembeds, tagembeds)
    opt = SGD()
    for epoch = 1:50
        println("Epoch:\t$epoch")
        opt.rate = 0.001 / (1 + 0.05*(epoch-1))
        # opt.rate = 0.00075

        train_data = makebatch(16, train_w, train_c, train_t)
        function train_f(data::Tuple)
            w, c, t = data
            y = ner.model(w, c, t)
            #softmax_crossentropy(t, y)
        end
        #train_data = collect(zip(train_w, train_c, train_t))
        loss = minimize!(train_f, opt, collect(zip(train_data...)))
        println("Loss:\t$loss")

        # test
        println("Testing...")
        function test_f(data::Tuple)
            w, c = data
            y = ner.model(w, c)
            vec(argmax(y.data,1))
        end
        test_data = collect(zip(test_w, test_c))
        pred = cat(1, map(test_f, test_data)...)
        gold = cat(1, map(t -> t.data, test_t)...)
        length(pred) == length(gold) || throw("Length mismatch.")

        ranges_p = decode(ner.tagset, pred)
        ranges_g = decode(ner.tagset, gold)
        fscore(ranges_g, ranges_p)
        println()
    end
end

function fscore(golds::Vector{T}, preds::Vector{T}) where T
    set = intersect(Set(golds), Set(preds))
    count = length(set)
    prec = round(count/length(preds), 5)
    recall = round(count/length(golds), 5)
    fval = round(2*recall*prec/(recall+prec), 5)
    println("Prec:\t$prec")
    println("Recall:\t$recall")
    println("Fscore:\t$fval")
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using UniqueVectors: UniqueVector

abstract type Geo end

struct GeoPoint <: Geo
    x::Float64
    y::Float64
    z::Float64
    mesh_size::Float64
end

struct GeoLine <: Geo
    p1::GeoPoint
    p2::GeoPoint
end

struct GeoCurve <: Geo
    lines::Vector{GeoLine}
end

struct GeoPlane <: Geo
    curve::GeoCurve
end

struct GeoSurfaceDifference <: Geo
    p1::GeoPlane
    p2::GeoPlane
    keep_tool::Bool
end

struct GeoSurfaceExtrusion <: Geo
    x::Float64
    y::Float64
    z::Float64
    plane::GeoPlane
end

struct GeoVolumeDifference <: Geo
    v1::GeoSurfaceExtrusion
    v2::GeoSurfaceExtrusion
    keep_tool::Bool
end

struct GeoSurfaceGroup <: Geo
    name::String
    planes::Vector{GeoPlane}
end

struct GeoVolumeGroup <: Geo
    name::String
    extrusions::Vector{GeoSurfaceExtrusion}
end

export GeoCode
"""
    GeoCode(points::UniqueVector{GeoPoint}, lines::UniqueVector{GeoLine}
        curves::UniqueVector{GeoCurve}, planes::UniqueVector{GeoPlane},
        surface_differences::UniqueVector{GeoSurfaceDifference},
        extrusions::UniqueVector{GeoSurfaceExtrusion},
        volume_differences::UniqueVector{GeoVolumeDifference},
        surface_groups::UniqueVector{GeoSurfaceGroup},
        volume_groups::UniqueVector{GeoVolumeGroup})
    GeoCode()

A representation of the .geo format for constructive solid geometry accepted by Gmsh.
"""
struct GeoCode <: Geo
    points::UniqueVector{GeoPoint}
    lines::UniqueVector{GeoLine}
    curves::UniqueVector{GeoCurve}
    planes::UniqueVector{GeoPlane}
    surface_differences::UniqueVector{GeoSurfaceDifference}
    extrusions::UniqueVector{GeoSurfaceExtrusion}
    volume_differences::UniqueVector{GeoVolumeDifference}
    surface_groups::UniqueVector{GeoSurfaceGroup}
    volume_groups::UniqueVector{GeoVolumeGroup}
end

import Base: ==
function ==(g1::Geo, g2::Geo)
    if typeof(g1) != typeof(g2)
        return false
    end
    for n in fieldnames(typeof(g1))
        if getfield(g1, n) != getfield(g2, n)
            return false
        end
    end
    return true
end

GeoCode() = GeoCode([UniqueVector(eltype(t)[]) for t in fieldtypes(GeoCode)]...)

import Base: append!
function append!(g1::GeoCode, g2::GeoCode)
    for n in fieldnames(GeoCode)
        arr = getfield(g1, n)
        for e in getfield(g2, n)
            if !(e in arr)
                push!(arr, e)
            end
        end
    end
end

function compile(gc::GeoCode, p::GeoPoint)
    t = findfirst(isequal(p), gc.points)
    return "Point($t) = {$(p.x), $(p.y), $(p.z), $(p.mesh_size)};"
end

function compile(gc::GeoCode, l::GeoLine)
    lt = findfirst(isequal(l), gc.lines)
    pt1, pt2 = findfirst(isequal(l.p1), gc.points), findfirst(isequal(l.p2), gc.points)
    return "Line($lt) = {$pt1, $pt2};"
end

function vector_to_string(l::Vector)
    s = "$(tuple(l...))"[2:end-1]
    if length(l) == 1
        return s[1:end-1]
    else
        return s
    end
end

function compile(gc::GeoCode, c::GeoCurve)
    ct = findfirst(isequal(c), gc.curves)
    lts = [findfirst(isequal(l), gc.lines) for l in c.lines]
    return "Curve Loop($ct) = {" * vector_to_string(lts) * "};"
end

function compile(gc::GeoCode, p::GeoPlane)
    pt = findfirst(isequal(p), gc.planes)
    ct = findfirst(isequal(p.curve), gc.curves)
    return "Plane Surface($pt) = {$ct};"
end

function compile(gc::GeoCode, sd::GeoSurfaceDifference)
    pt1 = findfirst(isequal(sd.p1), gc.planes)
    pt2 = findfirst(isequal(sd.p2), gc.planes)
    ending = sd.keep_tool ? "" : " Delete;"
    return "BooleanDifference{Surface{$pt1}; Delete;}{Surface{$pt2};$ending}"
end

function compile(gc::GeoCode, se::GeoSurfaceExtrusion)
    pt = findfirst(isequal(se.plane), gc.planes)
    return "Extrude {$(se.x), $(se.y), $(se.z)} {Surface{$pt};}"
end

function compile(gc::GeoCode, vd::GeoVolumeDifference)
    vt1 = findfirst(isequal(vd.v1), gc.extrusions)
    vt2 = findfirst(isequal(vd.v2), gc.extrusions)
    ending = vd.keep_tool ? "" : " Delete;"
    return "BooleanDifference{Volume{$vt1}; Delete;}{Volume{$vt2};$ending}"
end

function compile(gc::GeoCode, group::GeoSurfaceGroup)
    ts = [findfirst(isequal(p), gc.planes) for p in group.planes]
    return "Physical Surface(\"$(group.name)\") = {" * vector_to_string(ts) * "};"
end

function compile(gc::GeoCode, group::GeoVolumeGroup)
    ts = [findfirst(isequal(p), gc.extrusions) for p in group.extrusions]
    return "Physical Volume(\"$(group.name)\") = {" * vector_to_string(ts) * "};"
end

function compile(gc::GeoCode)
    code = String["//points\n"]
    for p in gc.points
        push!(code, compile(gc, p) * "\n")
    end
    push!(code, "\n//lines\n")
    for l in gc.lines
        push!(code, compile(gc, l) * "\n")
    end
    push!(code, "\n//curves\n")
    for c in gc.curves
        push!(code, compile(gc, c) * "\n")
    end
    push!(code, "\n//planes\n")
    for p in gc.planes
        push!(code, compile(gc, p) * "\n")
    end
    push!(code, "\n//surface differences\n")
    for sd in gc.surface_differences
        push!(code, compile(gc, sd) * "\n")
    end
    push!(code, "\n//surface extrusions\n")
    for se in gc.extrusions
        push!(code, compile(gc, se) * "\n")
    end
    push!(code, "\n//volume differences\n")
    for vd in gc.volume_differences
        push!(code, compile(gc, vd) * "\n")
    end
    push!(code, "\n//surface groups\n")
    for group in gc.surface_groups
        push!(code, compile(gc, group) * "\n")
    end
    push!(code, "\n//volume groups\n")
    for group in gc.volume_groups
        push!(code, compile(gc, group) * "\n")
    end
    return code
end

export geo_write!
"""
    geo_write!(file_name::String, gc::GeoCode=GeoCode(); random_factor=1e-9,
        random_factor_3D=1e-9, characteristic_length_factor::Real=1, footer::String="")

Write the contents of a GeoCode to a file along with additional geo code in `footer`.
The `random_factor` and `random_factor_3D` arguments adjust how Gmsh uses randomness in its
2D and 3D meshing algorithms, respectively. `characteristic_length_factor` scales the
characteristic length of the mesh; smaller values yield denser meshes.
"""
function geo_write!(file_name::String, gc::GeoCode=GeoCode(); random_factor=1e-9,
    random_factor_3D=1e-9, characteristic_length_factor::Real=1, footer::String="")
    code_strings = ["Mesh.RandomFactor=$random_factor;\n";
                    "Mesh.RandomFactor3D=$random_factor_3D;\n";
                    "Mesh.CharacteristicLengthFactor=$characteristic_length_factor;\n";
                    "SetFactory(\"OpenCASCADE\");\n\n";
                    compile(gc);
                    "\n//additional\n";
                    footer]
    file = Base.open(file_name, write=true)
    for s in code_strings
        write(file, s)
    end
    close(file)
end

"""
    GeoCode(points::AbstractVector{<:NTuple{4, <:Real}})

Generate the GeoCode for a closed 3D curve. The 4th element of each point is a mesh size
that Gmsh uses to spatially modulate the mesh density.
"""
function GeoCode(points::AbstractVector{<:NTuple{4, <:Real}})
    gc = GeoCode()
    for p in points
        push!(gc.points, GeoPoint(p...))
    end
    n = length(points)
    for i in 1:n
        p1, p2 = gc.points[i], gc.points[mod1(i+1, n)]
        push!(gc.lines, GeoLine(p1, p2))
    end
    if n > 0
        geo_curve = GeoCurve(gc.lines)
        push!(gc.curves, geo_curve)
        push!(gc.planes, GeoPlane(geo_curve))
    end
    return gc
end


text:
@testset "ConstantMean" begin
    N = 20
    x = rand()
    X = rand(N)
    c = rand()
    μ₀ = ConstantMean(c; opt=Descent(1.0))
    st = AGP.init_priormean_state((;), μ₀)
    @test μ₀ isa ConstantMean{Float64,<:Descent}
    @test repr("text/plain", μ₀) == "Constant Mean Prior (c = $c)"
    @test μ₀(X) == c .* ones(N)
    @test μ₀(x) == c
    g = Zygote.gradient(μ₀) do m
        sum(m(X))
    end
    AGP.update!(μ₀, st, Vector(first(g)[].C))
    @test μ₀.C[1] == (c + first(g)[].C[1])
end