text:
# # Introduction

# We first load [DrWatson](https://juliadynamics.github.io/DrWatson.jl/dev/),
# a scientific project management software that helps reproducibility.
using DrWatson

# We then activate this project, which also loads a bundle of packages
# (most notably ClimateBase.jl), and includes all files in the `src` folder.
@quickactivate :AlbedoProperties

# We then load the style file for the figures (text size, colors, etc.)
include(srcdir("plotting", "advanced.jl"))
include(papersdir("paperstyle.jl"))

# Load all necessary CERES files explicitly
EBAF_TOA = datadir("CERES", "EBAF_TOA.nc")
EBAF_SFC = datadir("CERES", "EBAF_SFC.nc")
SEAICE = datadir("CERES", "SYN1deg_AUX.nc")

# And from these files load all spatiotemporal fields necessary
R = ncread(EBAF_TOA, "toa_sw_all_mon")
K = ncread(EBAF_TOA, "toa_sw_clr_c_mon")
I = ncread(EBAF_TOA, "solar_mon")
F = ncread(EBAF_TOA, "cldarea_total_daynight_mon") ./ 100
T = ncread(EBAF_TOA, "cldtau_total_day_mon")
O = ncread(SEAICE, "aux_ocean_mon") ./ 100
E = ncread(SEAICE, "aux_snow_mon") ./ 100
L = clamp.(1 .- O .- E, 0, 1) # ice-free land
R

# The time vector is
t = dims(R, Time).val

# and the month counter is
tn = 1:length(t)

# And the following index is the most amount of data we can use for full-years
mys = maxyearspan(t)

# Notice that not all CERES data have this time span. Surface fluxes and ice fractions
# cover less amount of time, thus for operations involving this data we will infact be
# forced to use a smaller time span.

# Time in cumulative number of days per month is
truetime = daysinmonth.(t)
truetime = Float32.(cumsum(truetime))

# Ιn most figures we will be showing a sub-set of the total time (purely for visual clarity,
# and actual span chosen below is completely arbitrary)
tspan = (Date(2008), Date(2015))

# This is how we decompose
# some fields into northern-southern hemisphere average timeseries.
Ιn, Is = hemispheric_means(I)
Rn, Rs = hemispheric_means(R)
Rn

# This is a good point to demonstrate how critical it is for one to do proper temporal
# averaging. Integrated over a full yearly cycle, both hemispheres receive
# the same amount of insolation. This is a direct consequence of Kepler's second law.
timemean(Ιn) - timemean(Is)

# Using a standard mean instead, while still only including full years, gives
mean(Ιn[1:mys]) - mean(Is[1:mys])

# which is a clear violation of Kepler's second law.

###################################################################### #src
###################################################################### #src
# # Timeseries analysis of reflected solar rad.
###################################################################### #src
###################################################################### #src

# %% Seasonal decomposition of outgoing rad #src
# ## Seasonal decomposition
fig, axs = subplots(2, 1; figsize = (figx, 8))

sinusoidal = Sinusoidal(Float32.([1, 2] ./ DAYS_IN_ORBIT))
Rn_sea, Rn_res = decompose(Rn, sinusoidal)
Rs_sea, Rs_res = decompose(Rs, sinusoidal)

seasonals = (Rn_sea, Rs_sea)
residuals = (Rn_res, Rs_res)

n, s = seasonals
g = (n .+ s) ./ 2
axs[1].plot(truetime[1:24], n[1:24], label = "\$\\mathcal{N}\$", marker = "o", ms = 5)
axs[1].plot(truetime[1:24], s[1:24]; ls = "--", label = "\$\\mathcal{S}\$", marker = "o", ms = 5)
axs[1].set_xlabel("days")

n, s = residuals
g = (n .+ s) ./ 2
axs[2].plot(t, n, label = "\$\\mathcal{N}\$")
axs[2].plot(t, s; ls = "--", label = "\$\\mathcal{S}\$")
axs[2].legend(loc="upper right", ncol = 3)
sn, = plottrend!(axs[2], t, n; color = "C0", addlabel=false)
ss, = plottrend!(axs[2], t, s; color = "C1", addlabel=false, ls = "dashed")

axs[1].set_ylabel("\$Y\$ (W m⁻²)")
axs[2].set_ylabel("\$R'\$ (W m⁻²)")
fig.tight_layout()
fig.align_ylabels(axs)
fig.subplots_adjust(hspace = 0.4, left = 0.13, bottom = 0.1)
wsave(papersdir("figures", "seasonal"), fig) #src

# Slopes of residuals
sn, ss

# Percentage of variance of R that is due to seasonal components:
var(Rn_sea)/var(Rn), var(Rs_sea)/var(Rs)

# Percentage of variance of R that is due to residual components:
var(Rn_res)/var(Rn), var(Rs_res)/var(Rs)
# Percentage of variance of R that is due to covariability of seasonal-residual
cov(Rn_sea, Rn_res)/var(Rn), cov(Rs_sea, Rs_res)/var(Rs)

# ## Residuals are noise
# %% #src
using Random, TimeseriesSurrogates, Distributions
fig = figure(figsize = (figx, 2figy/3))

ax_gaus = fig.add_subplot(1,2,1)
ax_mi = fig.add_subplot(1,2,2)

x = Rn_res # change to Rs_res for SH
# detrend before plotting distribution
_, _, trend = linreg(1:length(x), x)
xdet = x .- trend
# fit normal distribution
gaussian = Distributions.fit(Normal, xdet)
# q-q plot
q1 = quantile.(Ref(xdet), 0:0.01:1)
q2 = quantile.(gaussian, 0:0.01:1)
ax_gaus.plot(q1, q2; marker = "o", ls = "none")
ax_gaus.plot(-3:3, -3:3; lw = 1, color = "C1")
ax_gaus.set_xlim(-3, 3)
ax_gaus.set_ylim(-3, 3)
ax_gaus.set_xlabel("quantile res.")
ax_gaus.set_ylabel("quantile G.")

# ### Calculate mutual information
using Entropies
sg = surrogenerator(x, TFTS(0.1))

x = Array(Rn_res)[1:end-1]
y = Array(Rn_res)[2:end]
function mutualinfo(x, y, ε)
    Hx = genentropy(Dataset(x), ε)
    Hy = genentropy(Dataset(y), ε)
    Hxy = genentropy(Dataset(x,y), ε)
    return Hx + Hy - Hxy
end
bins = std(x)/4 # bins used for MI
mi = mutualinfo(x, y, bins)
M = 10000;
mis = zeros(M)

for i in 1:M
    local s = sg()
    local x = s[1:end-1]
    local y = s[2:end]
    mis[i] = mutualinfo(x, y, bins)
end

μ, σ = meanstd(mis)
ax_mi.axvspan(μ-2σ, μ+2σ, color = "C2", alpha = 0.25)
ax_mi.hist(mis, 50, color = "C3", density = true, label = "surrogate")
ax_mi.set_yticks([])
ax_mi.grid(false)
ax_mi.axvline(mi, color = "C0", label = "real")
ax_mi.legend(loc = "lower right")
ax_mi.set_xlabel("self-MI")
ax_mi.set_ylabel("pdf")

fig.tight_layout()
add_identifiers!(fig)
fig.subplots_adjust(wspace = 0.2, hspace = 0.3, left = 0.14, bottom = 0.2, top =0.92)
wsave(papersdir("figures", "residuals_noise"), fig) #src
fig

# ## Distributions for ΔR values
# %% #src
using TimeseriesSurrogates, Random
fig = figure(;figsize = (figx, 3figy/2))
ax_pps = fig.add_subplot(3,1,1)
ax_ar = fig.add_subplot(3,1,2; sharex = ax_pps)
ax_difpps = fig.add_subplot(3,2,5)
ax_difar = fig.add_subplot(3,2,6)

# ### Pseudoperiodic surrogates
τ = 3 # delay time = when AC becomes zero
ρ = 0.5 # evaluated using `noiseradius` function
d = 3 # embedding dimension
pp = PseudoPeriodic(d, τ, ρ, false)
sgnh = surrogenerator(Array(Rn), pp)
sgsh = surrogenerator(Array(Rs), pp)

Random.seed!(5151525)
ax_pps.plot(t, Rn, label = "\$\\mathcal{N}(R)\$", marker = "o", ms = 4)
ax_pps.plot(t, sgnh(), color = "C3", ls = "dashed", label = "surrogate", marker = "o", ms = 4)
ax_pps.legend(ncol=2)
ax_pps.set_ylabel("W m⁻²")

M = 5000
differences = zeros(M)
for i in 1:M
    rnh = sgnh()
    rsh = sgsh()
    differences[i] = timeagg(mean, t, rnh .- rsh)
end

μ, σ = meanstd(differences)
ax_difpps.axvspan(μ-2σ, μ+2σ, color = "C2", alpha = 0.25)
ax_difpps.hist(differences, 50, color = "C3", density = true, label = "surrogate")
ax_difpps.axvline(timemean(Rn .- Rs), color = "C0", label = "real")
ax_difpps.set_yticks([])
ax_difpps.grid(false)
ax_difpps.set_xticks(-0.25:0.25:0.25)
ax_difpps.set_ylabel("pdf")
ax_difpps.set_xlabel("avg. diff. \$\\Delta R\$ (W m⁻²)")
ax_difpps.legend(loc = "lower left")

# ### Autoregressive approximation
using StatsBase, ARFIMA, DSP, Random

D = Rn_res .- Rs_res
n = length(D)
σ = std(D)

# Now we estimate the parameters for the autoregressive process.
# There are several ways to do it. One can use the Yule-Walker equations, or
# hyper-optimization on a predefined length of parameters `p`. Here we will use
# linear predictive code estimation.

lpcmethod = DSP.LPCLevinson()
φ, e = DSP.lpc(D, 12, lpcmethod)
φ = -SVector(φ...) # definition of `φ` in ARFIMA is with -

rng = Random.MersenneTwister(612354152) # reproducibility
X = arma(rng, n+length(φ), σ, φ)[length(φ)+1:end]

ax_ar.plot(t, D, label = "\$D = R_\\mathcal{N}' - R_\\mathcal{S}'\$")
ax_ar.plot(t, X, label = "autoregressive", color = "C3")
ax_ar.set_xlim(tspan...)
ax_ar.legend(loc = "lower right", ncol = 2)
ax_ar.set_ylabel("W m⁻²")
ax_ar.set_yticks(-4:2:4)

N = 100000
X = arma(N, σ, φ) # simulate ultra long AR process
X .-= mean(X) # assert mean 0 for clarity
possible = Float64[]
for i in n:N-n
    push!(possible, mean(view(X, i:i+n)))
end
q1 = quantile(possible, 0.05)
q2 = quantile(possible, 0.95)
ax_difar.axvspan(q1, q2, color = "C2", alpha = 0.25)
ax_difar.hist(possible, color = "C3", 100, label = "autoregr.")
ax_difar.axvline(0.1, color = "C0", label = "real")
ax_difar.legend(loc="lower left")
ax_difar.set_yticks([])
ax_difar.set_xticks(-0.25:0.25:0.25)
ax_difar.grid(false)
ax_difar.set_xlabel("avg. diff. \$\\Delta R\$ (W m⁻²)")

fig.tight_layout()
fig.subplots_adjust(wspace = 0.1, hspace = 0.3, left = 0.12, bottom = 0.1, top =0.95)
add_identifiers!(fig)
wsave(papersdir("figures", "difference_is_zero"), fig) #src
fig

# ## Mean albedo decomposition + Fourier (only R)
# %% #src
fig, axs = subplots(2,1)

i = 1; X = R
Xn, Xs = hemispheric_means(zonalmean(X))
solarn, internaln = decompose(t, Xn.data, MeanAlbedo(Ιn.data))
solars, internals = decompose(t, Xs.data, MeanAlbedo(Is.data))

αn = timemean(Rn)/timemean(Ιn)
αs = timemean(Rs)/timemean(Is)

axs[i].plot(t, internaln, label = "\$\\mathcal{N}\$")
axs[i].plot(t, internals; ls = "--", label = "\$\\mathcal{S}\$")

in, is = var(internaln)/var(Xn), var(internals)/var(Xs)
println("Percentage of variance of internal for $(X.name)")
@show in, is

cn, cs = 2cov(solarn, internaln)/var(Xn), 2cov(solars, internals)/var(Xs)
println("Percentage of variance of covariability with sun for $(X.name)")
@show cn, cs

sn, ss = var(solarn)/var(Xn), var(solars)/var(Xs)
println("Percentage of variance of covariability with sun for $(X.name)")
@show sn, ss
println()

axs[1].legend(loc="upper right"; ncol=1)
axs[1].set_ylabel("\$R_\\mathrm{internal}\$")
axs[1].set_yticks(-5:5:10)
axs[1].set_xlim(Date(2008), Date(2015))

method = :FFTW
for (j, r) in enumerate((internaln[1:mys], internals[1:mys]))
    ## This can be made into Lomb-Scargle periodogram for more accuracy
    ## since the time axis is non-equi-spaced.
    if method == :FFTW
        P = abs2.(rfft(r .- mean(r)))
        f = FFTW.rfftfreq(length(r))
        ps = 1 ./ Float64.(f)
    elseif method == :LombScargle
        truetime = daysinmonth.(t[1:mys])
        truetime = Float32.(cumsum(truetime))
        periodogram = lombscargle(truetime, Array(r))
        ps, P = periodpower(periodogram)
        ps ./= 30 # divide by number of days in month (approx.)
    end
    spi = sortperm(ps)
    ps = ps[spi]; P = P[spi]
    ## Limit periods shown to `maxmonth`
    maxmonth = 2*12 + 2
    ℓ, k = findnearest.((2, maxmonth), Ref(ps))
    ps = ps[ℓ:k]
    P = P[ℓ:k]
    P = abs.(P)
    axs[2].plot(ps, log10.(P),
    color = "C$(j-1)", ls = j == 1 ? "-" : "--")
end

axs[2].set_ylabel("\$\\log_{10}(|\\mathcal{F}|)\$")
axs[2].set_xlabel("period (months)")

fig.tight_layout()
fig.subplots_adjust(left = 0.12, bottom = 0.12, hspace=0.22, top = 0.96)
wsave(papersdir("figures", "internal_new"), fig) #src
fig

# ## Longitude split (quadrant of hemisphere)
# %% #src
λs = dims(R, Lon)[1:5:size(R, Lon)]
d = zeros(length(λs), 4)
X, Y = hemispheric_functions(R)
for (i, λ) in enumerate(λs)
    X1, X2 = split_at_longitude(X, λ)
    Y1, Y2 = split_at_longitude(Y, λ)
    NE, NW = timemean.(spacemean.((X1, X2)))
    SE, SW = timemean.(spacemean.((Y1, Y2)))

    d[i, 1] = NE - NW
    d[i, 2] = NE - SE
    d[i, 3] = NE - SW
    d[i, 4] = SE - SW
end

fig = figure()
ax = gca()
for (i, n) in zip(1:4, ("NE - NW", "NE - SE", "NE - SW", "SE - SW", "NW - SW", "SE - SW"))
    z = d[:, i]
    a = mean(abs.(z))
    ax.plot(λs, z, label = "$(n), a = $(rdspl(a))")

end

ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc="lower left",
          ncol=2, mode="expand", borderaxespad=0.)
ax.set_xlabel("split longitude \$\\lambda\$ (ᵒ)")
ax.set_ylabel("avg. diff. (W m⁻²)")
ax.set_yticks(-6:3:6)
fig.tight_layout()

wsave(papersdir("figures", "longitude_split"), fig) #src

###################################################################### #src
###################################################################### #src
# # Cloudiness
###################################################################### #src
###################################################################### #src

# %% Compare the two different measures of Cloudiness on a MAP #src
# ## Definition and comparison

Tcorrected = sinusoidal_continuation(T, [1, 2.0]; Tmin = 0)
C1 = cloud_effective_albedo(F, Tcorrected, 0.9)
C1_timeavg = timemean(C1)

# use `\:arrow_down:` and tab for arrows
F_s_⬆ = ncread(EBAF_SFC, "sfc_sw_up_all_mon")
F_s_⬇ = ncread(EBAF_SFC, "sfc_sw_down_all_mon")
F_s_⬆_K = ncread(EBAF_SFC, "sfc_sw_up_clr_c_mon")
F_s_⬇_K = ncread(EBAF_SFC, "sfc_sw_down_clr_c_mon");

# Because surface data do not have the same amount of points, we are forced to use a
# smaller time span. (even if we time-average fields, this is still the accurate thing to do)
l = size(F_s_⬆, Time)
argsall = timemean.((I[Time(1:l)], R[Time(1:l)], F_s_⬆, F_s_⬇))
argsclr = timemean.((I[Time(1:l)], K[Time(1:l)], F_s_⬆_K, F_s_⬇_K))
α_ATM, α_SFC = surface_atmosphere_contributions(argsall...) # all sky
α_ATM_K, α_K_SFC = surface_atmosphere_contributions(argsclr...) # clear sky

C2_timeavg = α_ATM .- α_ATM_K

if false # plot spatial fields
    cmap = "gist_earth"
    as = (α_ATM, α_ATM_K, C2_timeavg, C1_timeavg, C1_timeavg .- C2_timeavg, α_ATM .+ α_SFC)
    nmes = ("α_ATM", "α_ATM_K", "C₂", "C₁", "C₁ - C₂", "α_TOA")
    kwargs = (vmin = 0, levels = 21, vmax = 0.7, cmap = "gist_earth")
    for (a, name) in zip(as, nmes)
        @show name, extrema(a)
        if name != "C₁ - C₂"
            fig, ax, = earthsurface(a; cmap = cmap, vmin = 0, vmax = 0.7)
        else
            fig, ax, = earthsurface(a; cmap = cmap, vmin = -0.2, vmax = 0.2)
        end
        ax.set_title(name)
    end
end

# We now need to estimate `g` so that the two versions match in value.
# Ιn the paper we call this the energetically consistent modification.
a = C2_timeavg; τ = timemean(Tcorrected); f = timemean(F)
_x = -2a ./ (τ .* a - f .* τ)
g = clamp.(1 .- _x ./ √3, 0, 1)

C = cloud_effective_albedo(F, Tcorrected, g)

# Let's produce maps of f, τ, g, C (for appendix)
for (label, X) in zip(("f", "\\tau", "g", "C"), (f, τ, g, timemean(C)))
    m = spacemean(X)
    fig, ax, = earthsurface(X; cmap = "gist_earth")
    ax.set_title("\$$(label)\$, mean = $(rdspl(m))")
    wsave(papersdir("figures", "map_$label"), fig) #src
end

# Average cloud albedo per hemisphere:
Cn, Cs = hemispheric_means(C)
timemean.((Cn, Cs))

# ## Plot of zonally average cloudiness, along with average values
# %% #src
fig = figure(figsize = (figx, 5))
Cz = timemean(zonalmean(C))
ax = gca()
ax.plot(sind.(dims(Cz, Lat)), Cz, label = "\$\\mathcal{T}(C)\$")
Fz = timemean(zonalmean(F))
Fz = (Fz .- mean(Fz)) ./ std(Fz)
Fz = Fz .* std(Cz) .+ mean(Cz)
ax.plot(sind.(dims(Cz, Lat)), Fz, ls = "dashed", color = "C2", label = "\$\\mathcal{T}(f)\$")
ax.legend(ncol = 2)

lats = dims(Cz, Lat) |> Array

latr = [-90, -30, 0, 30, 90]
s = 0.0
for i in 1:4
    l1, l2 = latr[i], latr[i+1]
    mv = latmean(Cz[Lat(Between(l1, l2))])
    x = (sind(l1) + sind(l2))/2
    ax.text(x, 0.025, rdspl(mv, 3), ha = :center)
    global s += mv
    ax.axvspan(sind(l1), sind(l2); alpha = 0.1, color = "C$(i+1)")
end
ax.set_ylim(0, 0.4)
llats = -70:20:70
ax.set_xticks(sind.(llats))
ax.set_xticklabels(llats)
ax.set_xlim(-1,1)
ax.set_xlabel("latitude");

fig.tight_layout()
fig.subplots_adjust(left = 0.1, bottom = 0.2, top = 0.95)
wsave(papersdir("figures", "cloud_lat"), fig) #src

# ## Plot of cloudiness timeseries and correlation with solar cycle
# %% #src

fig, axs = subplots(2,1; figsize = (figx, 7))
ax = axs[1]
ax.plot(t, Cn, label = "\$\\mathcal{N}(C)\$")
ax.plot(t, Cs; ls = "dashed", label = "\$\\mathcal{S}(C)\$")
ax.legend(loc = "upper right")
ax.set_yticks(0.14:0.03:0.20)
ax.set_xlim(tspan...)
ax.set_ylabel("\$C\$")

Cx = (Cn, Cs)
Ix = (Ιn, Is)
ax = axs[2]
lags = -13:13

monthlpv = Sinusoidal([1/12])
for j in 1:2
    cc = crosscor(Cx[j], Ix[j], lags)
    decompose(float.(lags), cc, monthlpv)
    φ = mod2pi(monthlpv.φ[2])
    φ = 12φ/2π # in months
    ax.plot(lags, cc, ls = j==1 ? "-" : "--",
    label="\$\\mathcal{$(j==1 ? "N" : "S")}, d \\approx $(rdspl(φ, 2))\$")
end
ax.set_xlabel("delay (months)"); ax.set_xticks(-12:3:12);
ax.legend(loc = "lower right")
ax.set_ylabel("crosscor.")

fig.tight_layout()
fig.subplots_adjust(left = 0.13, bottom = 0.15, top = 0.95)
add_identifiers!(fig)
fig.align_ylabels()
wsave(papersdir("figures", "cloud_and_sun"), fig) #src

# Weighted integrals of sunlight and cloudiness
D_with_sun = timeagg(mean, Cs, Is) - timeagg(mean, Cn, Ιn)
G_with_sun = 0.5(timeagg(mean, Cs, Is) + timeagg(mean, Cn, Ιn))
P_with_sun = 100 * D_with_sun/G_with_sun # proportion in percentage

D_no_sun = timeagg(mean, Cs) - timeagg(mean, Cn)
G_no_sun = 0.5(timeagg(mean, Cs) + timeagg(mean, Cn))
P_no_sun = 100 * D_no_sun/G_no_sun # proportion in percentage

P_with_sun, P_no_sun

# Alright, we now quantify the extra reflection one gets from the time delay
# by shifting the NH timeseries until to matches the sun.
using Interpolations
cloud_reflection_n = timemean(Cn .* Ιn)
cloud_reflection_s = timemean(Cs .* Is)
t = dims(Ιn, Time) |> Array
interp = CubicSplineInterpolation(1:length(t), Cn.data; extrapolation_bc = Line())

result = []
ds = 0:0.1:1.5
for d in ds
    r = timemean([interp(i) for i in (1:length(t)) .+ d] .* Ιn)
    push!(result, r)
end

fig = figure()
plot(ds, result, label = "Cloud reflection of NH")
axhline(cloud_reflection_n, ls = "dashed", label = "measured NH reflection", color = "C3")
axvline(1.2, color = "C1", label = "SH delay")
axvline(1.5, color = "C2", label = "0 delay")
ylabel("cloud reflection")
xlabel("delay shift (months)")
legend()
tight_layout()

# The effect of shifting the timeseries of NH by 1.2 months (so that they have same
# delay as SH) compensates the asymmetry by:
effect = result[13] - result[1]

# ## Cloudiness over specific land type
# %% #src
# Because `O, L, E` have a shorter timespan than `C`, we need to modify `C` to be
# in the same time span.

function hemispheric_weighted_avgs(X, fields)
    @assert size(X) == size(fields[1])
    Xn, Xs = hemispheric_functions(X)
    Xnt, Xst = timemean.((Xn, Xs))
    northpie = []
    southpie = []
    northcov = []
    southcov = []

    for W in fields
        Wn, Ws = hemispheric_functions(W)
        Gn, Gs = spacemean(Xn, Wn), spacemean(Xs, Ws)
        Gna, Gsa = timemean.((Gn, Gs))
        Wnt, Wst = timemean.((Wn, Ws))
        push!(northpie, (Gna, spacemean(Wnt)))
        push!(southpie, (Gsa, spacemean(Wst)))
        ## calcualte no covariability
        push!(northcov, spacemean(Xnt, Wnt))
        push!(southcov, spacemean(Xst, Wst))
    end
    return northpie, southpie, northcov, southcov
end

ol = size(O, Time)
northpie, southpie, northcov, southcov = hemispheric_weighted_avgs(C[Time(1:ol)], (O, E, L));

# Put the numbers into bar plot
# %% Barplot over surface #src
fig = figure()
ax = fig.add_subplot()
colors = ["C2", "C3", "C0"]
fig.tight_layout()
prefix = 1.2

for (i, data) in enumerate((northpie, southpie))
    pre = i == 1 ? 0.0 : prefix
    widths = [s[2] for s in data]
    x = [0, cumsum(widths)[1:end-1]...] .+ pre
    h = [s[1] for s in data]
    hcov = (northcov, southcov)[i]

    ax.bar(x, h; width = widths,
    color = colors, align = "edge")
    ax.bar(x, hcov; width = widths, linewidth = 3.0,
    color = colors, align = "edge", ls = "dashed", fill = false, edgecolor = lighten.(colors, 0.55))
    pos = x .+ widths./2
    for (j, (c, w)) in enumerate(data)
        str = "$(rdspl(c, 2))\n$(rdspl(w, 2))"
        off = + 0.015
        i == 2 && j == 2 && (off += 0.01)
        ax.text(pos[j], h[j] + off, string(rdspl(c, 2)), ha = :center, color = colors[j])
        botoff = (i == 2 && j == 2) ? +0.01 : -0.02
        ax.text(pos[j], botoff, string(rdspl(w, 2)), ha = :center)
    end
end
ax.spines["right"].set_visible(false)
ax.spines["top"].set_visible(false)
ax.set_xticks([])
ax.spines["bottom"].set_visible(false)

ax.set_ylim(0, 0.25)
ax.set_ylabel("\$\\langle C \\rangle\$ over surface")
ax.text(0.5, 0.25, "northern", ha = :center)
ax.text(0.5 + prefix, 0.25, "southern", ha = :center)
fig.tight_layout()

## Create legend
mpatches = matplotlib.patches
handles = [mpatches.Patch(color=c, label=F) for (c, F) in zip(colors, (:ocean, :ice, :land))]

ax.legend(handles=handles; ncol=3, bbox_to_anchor=(0.0, -0.3, 1, .102), loc="lower left",
            mode="expand", borderaxespad=0.)

fig.tight_layout()
fig.subplots_adjust(bottom = 0.24)

wsave(papersdir("figures", "cloud_over_surface_bar"), fig) #src

# %% total albedo given surface and cloud #src
# ## Total albedo given surface and cloud albedo and transmittance
# Main goal here is to see impact of having cloud albedo difference of `0.02`.
transmitance = 0.53 # some given trasmittance

asfc = 0:0.01:1.0
aatm = 0:0.01:transmitance

achoice = [0.18, 0.2, 0.22]
schoice = [0.1, 0.3, 0.6]

fig = figure()
ax1 = fig.add_subplot(1, 2, 1)
atoa = total_toa_albedo.(aatm, asfc', transmitance)
ax1.pcolormesh(aatm, asfc, atoa')
ax1.set_title("transmittance = $(transmitance)")
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 4)

for i in 1:length(achoice)
    a = achoice[i]
    ax1.axvline(a, color = "C$i", ls = "dashed")
    atotal = total_toa_albedo.(a, asfc, transmitance)
    ax2.plot(asfc, atotal, color = "C$i", label = "$(a)")
    ax2.set_xlabel("a_sfc")
    ax2.legend()
end

for i in 1:length(schoice)
    a = schoice[i]
    ax1.axhline(a, color = "C$i", ls = "dashed")
    atotal = total_toa_albedo.(aatm, a, transmitance)
    ax3.plot(aatm, atotal, color = "C$i", label = "$a")
    ax3.set_xlabel("a_atm")
    ax3.legend()
end

fig.tight_layout()

# So for a broad range of values (and you have to zoom into the figure), we find a
# total change of 0.02. This is also the surface abledo assymetry
Kn, Ks = hemispheric_means(K)
δαK = timemean(Kn ./ Ιn) - timemean(Ks ./ Is)

###################################################################### #src
###################################################################### #src
# # Appendices
###################################################################### #src
###################################################################### #src
# ## Continuation of value of τ
# %% #src
count(ismissing, T)/length(T) # about 6% of the points are missing

fig, axs = subplots(2, 1; figsize = (figx, 6), sharex = true)
for (i, l) in enumerate((-82, 82))
    ax = axs[i]
    example = T[Lat(Near(l)), Lon(Near(5))]
    example = Vector{Any}(Array(example))

    lpv = Sinusoidal(Float32.([1, 2] ./ DAYS_IN_ORBIT))
    nmi = findall(!ismissing, example)
    mi = findall(ismissing, example)
    x0 = Array{Float32}(example[nmi])
    t0 = truetime[nmi]
    sea, res = decompose(t0, x0, lpv)
    c = ClimateBase.cosines(truetime, lpv)
    acosine = copy(example)
    acosine[mi] .= c[mi]
    example[mi] .= nothing

    ax.plot(t, example, marker = "o", label = "lat = $(l)", ms = 4)
    ax.plot(t, c, ls = "dotted", alpha = 0.75, label = "cosine cont.")
    ax.legend(loc = "upper right")
    ax.set_ylabel("\$\\tau\$")
end
axs[2].set_xlim(tspan...)
fig.tight_layout()
fig.subplots_adjust(left = 0.1, right = 0.95, bottom = 0.14, hspace=0.1, top = 0.95)
wsave(papersdir("figures", "tau_continuation"), fig) #src

# Finally, make a Jupyter notebook out of the entire file. #src
using Literate #src
Literate.notebook(papersdir("figures.jl"), projectdir(); execute = false, credit = false) #src

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
using MyMeshIO
using Test

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

ext: jl

max_stars_repo_name: powerzf/MyMeshIO

max_issues_repo_name: powerzf/MyMeshIO

max_forks_repo_name: powerzf/MyMeshIO


text:
"""
    mae(ŷ, y; agg=mean)

Return the loss corresponding to mean absolute error:

    agg(abs.(ŷ .- y))
"""
mae(ŷ, y; agg=mean) = agg(abs.(ŷ .- y))

"""
    mse(ŷ, y; agg=mean)

Return the loss corresponding to mean square error:

    agg((ŷ .- y).^2)
"""
mse(ŷ, y; agg=mean) = agg((ŷ .- y).^2)

"""
    msle(ŷ, y; agg=mean, ϵ=eps(ŷ))

The loss corresponding to mean squared logarithmic errors, calculated as

    agg((log.(ŷ .+ ϵ) .- log.(y .+ ϵ)).^2)

The `ϵ` term provides numerical stability.
Penalizes an under-estimation more than an over-estimatation.
"""
msle(ŷ, y; agg=mean, ϵ=epseltype(ŷ)) = agg((log.((ŷ .+ ϵ) ./ (y .+ ϵ))).^2)

"""
    huber_loss(ŷ, y; δ=1, agg=mean)

Return the mean of the [Huber loss](https://en.wikipedia.org/wiki/Huber_loss)
given the prediction `ŷ` and true values `y`.

                 | 0.5 * |ŷ - y|^2,            for |ŷ - y| <= δ
    Huber loss = |
                 |  δ * (|ŷ - y| - 0.5 * δ), otherwise
"""
function huber_loss(ŷ, y; agg=mean, δ=ofeltype(ŷ, 1))
   abs_error = abs.(ŷ .- y)
   #TODO: remove dropgrad when Zygote can handle this function with CuArrays
   temp = Zygote.dropgrad(abs_error .<  δ)
   x = ofeltype(ŷ, 0.5)
   agg(((abs_error.^2) .* temp) .* x .+ δ*(abs_error .- x*δ) .* (1 .- temp))
end

"""
    crossentropy(ŷ, y; dims=1, ϵ=eps(ŷ), agg=mean)

Return the cross entropy between the given probability distributions;
calculated as

    agg(-sum(y .* log.(ŷ .+ ϵ); dims=dims))

Cross entropy is tipically used as a loss in multi-class classification,
in which case the labels `y` are given in a one-hot format.
`dims` specifies the dimension (or the dimensions) containing the class probabilities.
The prediction `ŷ` is supposed to sum to one across `dims`,
as would be the case with the output of a [`softmax`](@ref) operation.

Use of [`logitcrossentropy`](@ref) is recomended over `crossentropy` for
numerical stability.

See also: [`Flux.logitcrossentropy`](@ref), [`Flux.binarycrossentropy`](@ref), [`Flux.logitbinarycrossentropy`](@ref)
"""
function crossentropy(ŷ, y; dims=1, agg=mean, ϵ=epseltype(ŷ))
    agg(.-sum(xlogy.(y, ŷ .+ ϵ); dims=dims))
end

"""
    logitcrossentropy(ŷ, y; dims=1, ϵ=eps(ŷ), agg=mean)

Return the crossentropy computed after a [`Flux.logsoftmax`](@ref) operation;
calculated as

    agg(.-sum(y .* logsoftmax(ŷ; dims=dims); dims=dims))

`logitcrossentropy(ŷ, y)` is mathematically equivalent to
[`Flux.Losses.crossentropy(softmax(ŷ), y)`](@ref) but it is more numerically stable.

See also: [`Flux.Losses.crossentropy`](@ref), [`Flux.Losses.binarycrossentropy`](@ref), [`Flux.Losses.logitbinarycrossentropy`](@ref)
"""
function logitcrossentropy(ŷ, y; dims=1, agg=mean)
    agg(.-sum(y .* logsoftmax(ŷ; dims=dims); dims=dims))
end

"""
    binarycrossentropy(ŷ, y; agg=mean, ϵ=eps(ŷ))

Return the binary cross-entropy loss, computed as

    agg(@.(-y*log(ŷ + ϵ) - (1-y)*log(1-ŷ + ϵ)))

The `ϵ` term provides numerical stability.

Typically, the prediction `ŷ` is given by the output of a [`sigmoid`](@ref) activation.

Use of `logitbinarycrossentropy` is recomended over `binarycrossentropy` for numerical stability.

See also: [`Flux.Losses.crossentropy`](@ref), [`Flux.Losses.logitcrossentropy`](@ref), [`Flux.Losses.logitbinarycrossentropy`](@ref)
"""
function binarycrossentropy(ŷ, y; agg=mean, ϵ=epseltype(ŷ))
    agg(@.(-xlogy(y, ŷ+ϵ) - xlogy(1-y, 1-ŷ+ϵ)))
end
# Re-definition to fix interaction with CuArrays.
# CUDA.@cufunc binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

"""
    logitbinarycrossentropy(ŷ, y; agg=mean)

Mathematically equivalent to
[`Flux.binarycrossentropy(σ(ŷ), y)`](@ref) but is more numerically stable.

See also: [`Flux.Losses.crossentropy`](@ref), [`Flux.Losses.logitcrossentropy`](@ref), [`Flux.Losses.binarycrossentropy`](@ref)
```
"""
function logitbinarycrossentropy(ŷ, y; agg=mean)
    agg(@.((1-y)*ŷ - logσ(ŷ)))
end
# Re-definition to fix interaction with CuArrays.
# CUDA.@cufunc logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

"""
    kldivergence(ŷ, y; agg=mean)

Return the
[Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence)
between the given probability distributions.

KL divergence is a measure of how much one probability distribution is different
from the other.
It is always non-negative and zero only when both the distributions are equal
everywhere.
"""
function kldivergence(ŷ, y; dims=1, agg=mean, ϵ=epseltype(ŷ))
  entropy = agg(sum(xlogx.(y), dims=dims))
  cross_entropy = crossentropy(ŷ, y; dims=dims, agg=agg, ϵ=ϵ)
  return entropy + cross_entropy
end

"""
    poisson_loss(ŷ, y)

# Return how much the predicted distribution `ŷ` diverges from the expected Poisson
# distribution `y`; calculated as `sum(ŷ .- y .* log.(ŷ)) / size(y, 2)`.
REDO
[More information.](https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/poisson_loss).
"""
poisson_loss(ŷ, y; agg=mean) = agg(ŷ .- xlogy.(y, ŷ))

"""
    hinge_loss(ŷ, y; agg=mean)

Return the [hinge_loss loss](https://en.wikipedia.org/wiki/Hinge_loss) given the
prediction `ŷ` and true labels `y` (containing 1 or -1); calculated as
`sum(max.(0, 1 .- ŷ .* y)) / size(y, 2)`.

See also: [`squared_hinge_loss`](@ref)
"""
hinge_loss(ŷ, y; agg=mean) = agg(max.(0, 1 .-  ŷ .* y))

"""
    squared_hinge_loss(ŷ, y)

Return the squared hinge_loss loss given the prediction `ŷ` and true labels `y`
(containing 1 or -1); calculated as `sum((max.(0, 1 .- ŷ .* y)).^2) / size(y, 2)`.

See also: [`hinge_loss`](@ref)
"""
squared_hinge_loss(ŷ, y; agg=mean) = agg((max.(0, 1 .- ŷ .* y)).^2)

"""
    dice_coeff_loss(ŷ, y; smooth=1)

Return a loss based on the dice coefficient.
Used in the [V-Net](https://arxiv.org/abs/1606.04797) image segmentation
architecture.
Similar to the F1_score. Calculated as:

    1 - 2*sum(|ŷ .* y| + smooth) / (sum(ŷ.^2) + sum(y.^2) + smooth)
"""
dice_coeff_loss(ŷ, y; smooth=ofeltype(ŷ, 1.0)) = 1 - (2*sum(y .* ŷ) + smooth) / (sum(y.^2) + sum(ŷ.^2) + smooth) #TODO agg

"""
    tversky_loss(ŷ, y; β=0.7)

Return the [Tversky loss](https://arxiv.org/abs/1706.05721).
Used with imbalanced data to give more weight to false negatives.
Larger β weigh recall more than precision (by placing more emphasis on false negatives)
Calculated as:
    1 - sum(|y .* ŷ| + 1) / (sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1)
"""
function tversky_loss(ŷ, y; β=ofeltype(ŷ, 0.7))
    #TODO add agg
    num = sum(y .* ŷ) + 1
    den = sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1
    1 - num / den
end


text:
"""
    DateSelector

Determines how to [`partition`](@ref) a date set into disjoint validation and holdout sets.
"""
abstract type DateSelector end

"""
    partition(dates::AbstractInterval{Date}, s::DateSelector)
    partition(dates::StepRange{Date, Day}, selector::DateSelector)

Partition the set of `dates` into disjoint `validation` and `holdout` sets according to the
`selector` and return a `NamedTuple({:validation, :holdout})` of iterators.
"""
function Iterators.partition(dates::AbstractInterval{Date}, s::DateSelector)
    _dates = _interval2daterange(dates)
    return partition(_dates, s)
end

"""
    _getdatesets(st, ed, dates) -> NamedTuple{(:validation, :holdout)}

Construct the NamedTuple of iterators for the validation and holdout date sets.
"""
function _getdatesets(all_dates, holdout_dates)
    return (
        validation=sort(setdiff(all_dates, holdout_dates)),
        holdout=sort(holdout_dates)
    )
end

_getdatesets(all_dates, date::Date) = _getdatesets(all_dates, [date])

"""
    _interval2daterange(dates::AbstractInterval{Day}) -> StepRange{Date, Day}

Helper function to turn an AbstractInterval into a StepRange taking the inclusivity into
account.
"""
function _interval2daterange(dates::AbstractInterval{Date})
    fd = _firstdate(dates)
    ld =_lastdate(dates)
    return fd:Day(1):ld
end

# TODO: remove this once https://github.com/invenia/Intervals.jl/issues/137
# is addressed.
_firstdate(dates::AbstractInterval{Date,Closed}) = first(dates)
_firstdate(dates::AbstractInterval{Date,Open}) = first(dates) + Day(1)
_lastdate(dates::AbstractInterval{Date,<:Bound,Closed}) = last(dates)
_lastdate(dates::AbstractInterval{Date,<:Bound,Open}) = last(dates) - Day(1)

"""
    _initial_date(s::S, dates) where S<: DateSelector

Determines when we start counting from when breaking `dates` up into blocks.
Checks that that initial time is valid for the given `dates`.
"""
function _initial_date(s::S, dates) where S<: DateSelector
    sd, ed = extrema(dates)

    # We would like to start from over 100 years ago
    # 1900, which has the convient feature of starting on a Monday
    beginning_of_time = Date(1900)
    initial_time = beginning_of_time + s.offset
    sd < initial_time && throw(DomainError(
        sd,
        "$S with offset $(s.offset) cannot be used before $(initial_time)",
    ))
    return initial_time
end

function Base.parse(DT::Type{<:DateSelector}, s::AbstractString)
    expr = Meta.parse(s)

    if expr isa Expr && expr.head == :call && isdefined(@__MODULE__, expr.args[1])
        T = getfield(@__MODULE__, expr.args[1])
        T isa Type && T <: DT && return eval(expr)
    end

    throw(ArgumentError("Could not parse \"$s\" as a `$DT`"))
end

if VERSION < v"1.5"
    function Base.show(io::IO, s::S) where S<:DateSelector
        args = _stringarg.(getfield.(Ref(s), fieldnames(S)))

        print(io, "$S($(join(args, ", ")))")
    end

    # Periods print in a format that breaks our parsing
    # So we undo their formatting
    _stringarg(p::Period) = "$(typeof(p))($(p.value))"
    _stringarg(arg) = repr(arg)
end


text:
"""
    watch(f, cli, key; wait_index=-1, recursive=false)

Creates an asynchronous `Task` which watches the `key` on the etcd cluster and runs
function `f` on the response.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#waiting-for-a-change)
for more details.
"""
function watch(f::Function, cli::Client, key::String; wait_index::Int=-1, recursive::Bool=false)
    t = @async begin
        opts = Dict{String, Any}("wait" => true)

        if recursive
            opts["recursive"] = recursive
        end

        if wait_index > 0
            opts["waitIndex"] = wait_index
        end

        resp = get(cli, key, opts; recursive=recursive)
        f(resp)
    end

    return t
end

"""
    watchloop(f, cli, key, [p]; wait_index=-1, recursive=false)

Creates an asynchrous `Task` which continously watches the `key` on the etcd cluster and
runs function `f` on the response.
The predicate function `p` represents a termination condition to exit the loop.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#waiting-for-a-change)
for more details.
"""
function watchloop(
    f::Function, cli::Client, key::String, p::Function=(resp) -> false;
    wait_index::Int=-1, recursive::Bool=false
)
    t = @async begin
        while true
            opts = Dict{String, Any}("wait" => true)

            if recursive
                opts["recursive"] = recursive
            end

            if wait_index > 0
                opts["waitIndex"] = wait_index
            end

            resp = get(cli, key, opts; recursive=recursive)
            f(resp)

            if p(resp)
                break
            end

            # reset wait_index to modified index of the response
            # - to avoid duplicates when the etcd index is incremented by updates elsewhere (outside the watched path)
            # - and missed updates when wait_index was not initially supplied (for ones that happen between calls to `get`)
            if ("node" in keys(resp)) && ("modifiedIndex" in keys(resp["node"]))
                wait_index = resp["node"]["modifiedIndex"]
            end

            if wait_index > 0
                wait_index += 1
            end
        end
    end

    return t
end


text:
#=

utilities for tribe

Ignacio Quintero Mächler

t(-_-t)

May 01 2017

=#

"""
    σ²ϕprop()

Generate proportional proposals for σ²
using random samples from **LogNormal** distributions.
"""
σ²ϕprop() = exp(randn() - 0.1)

"""
    λϕprop()

Generate proportional proposals for λ
using random samples from **LogNormal** distributions.
"""
function λϕprop()

  lg = exp(randn())

  return lg*exp(randn()*0.3 - 0.044),
         lg*exp(randn()*0.3 - 0.044)
end

"""
    coinsamp(p0::Float64)

Generate one sample from a random **Binomial** variable
with probability of failure `p0`.
"""
coinsamp(p0::Float64) = rand() < p0 ? 0 : 1

"""
    rexp(λ::Float64)

Generate one random sample from a **Exponential** distribution
with mean `λ`.
"""
rexp(λ::Float64) = (randexp()/λ)::Float64

"""
    normlize(pt1::Float64, pt2::Float64)

Normalize probabilities to 1.
"""
normlize(p1::Float64, p2::Float64) = p1/(p1 + p2)

"""
    normlize(p::Tuple{Float64,Float64})

Normalize probabilities to 1.
"""
normlize(p::Tuple{Float64,Float64}) = p[1]/(p[1] + p[2])

"""
    Ptrfast(λ1::Float64, λ0::Float64, t::Float64)

Return Markov chain probabilities for
two states through fast analytical solution
after time `t` for rates of gain `λ1` and loss `λ0`.
"""
function Ptrfast(λ1::Float64,
                 λ0::Float64,
                 t::Float64)

    ex  ::Float64 = exp(-(λ1 + λ0)*t)
    sd1 ::Float64 = 1.0/(λ1 + λ0)
    λ1ex::Float64 = λ1*ex
    λ0ex::Float64 = λ0*ex

    ((sd1*(λ0 + λ1ex), sd1*(λ1 - λ1ex)),
     (sd1*(λ0 - λ0ex), sd1*(λ1 + λ0ex)))
end

"""
    Ptrfast_start(λ1::Float64, λ0::Float64, t::Float64, state::Int64)

Return Markov chain probabilities for
two states through fast analytical solution
after time `t` for rates of gain `λ1` and loss `λ0`,
for start value of 0.
"""
function Ptrfast_start(λ1   ::Float64,
                       λ0   ::Float64,
                       t    ::Float64,
                       ::Type{Val{0}})

  ex  ::Float64 = exp(-(λ1 + λ0)*t)
  sd1 ::Float64 = 1.0/(λ1 + λ0)
  λ1ex::Float64 = λ1*ex

  return sd1*(λ0 + λ1ex), sd1*(λ1 - λ1ex)
end

"""
    Ptrfast_start(λ1::Float64, λ0::Float64, t::Float64, state::Int64)

Return Markov chain probabilities for
two states through fast analytical solution
after time `t` for rates of gain `λ1` and loss `λ0`,
for start value of 1.
"""
function Ptrfast_start(λ1   ::Float64,
                       λ0   ::Float64,
                       t    ::Float64,
                       ::Type{Val{1}})

  ex  ::Float64 = exp(-(λ1 + λ0)*t)
  sd1 ::Float64 = 1.0/(λ1 + λ0)
  λ0ex::Float64 = λ0*ex

  return sd1*(λ0 - λ0ex), sd1*(λ1 + λ0ex)
end

"""
    Ptrfast_end(λ1::Float64, λ0::Float64, t::Float64, state::Int64)

Return Markov chain probabilities for
two states through fast analytical solution
after time `t` for rates of gain `λ1` and loss `λ0`,
for end value of 0.
"""
function Ptrfast_end(λ1   ::Float64,
                     λ0   ::Float64,
                     t    ::Float64,
                     ::Type{Val{0}})

  ex ::Float64 = exp(-(λ1 + λ0)*t)
  sd1::Float64 = 1.0/(λ1 + λ0)

  return sd1*(λ0 + λ1*ex), sd1*(λ0 - λ0*ex)
end

"""
    Ptrfast_end(λ1::Float64, λ0::Float64, t::Float64, state::Int64)

Return Markov chain probabilities for
two states through fast analytical solution
after time `t` for rates of gain `λ1` and loss `λ0`,
for end value of 1.
"""
function Ptrfast_end(λ1   ::Float64,
                     λ0   ::Float64,
                     t    ::Float64,
                     ::Type{Val{1}})

  ex ::Float64 = exp(-(λ1 + λ0)*t)
  sd1::Float64 = 1.0/(λ1 + λ0)

  return sd1*(λ1 - λ1*ex), sd1*(λ1 + λ0*ex)
end

"""
    idxlessthan(x::Array{Float64,1}, val::Float64)

Get index in sorted vector `x` corresponding to the value
that is closest to but less than `val` in sorted arrays
using a sort of uniroot algorithm.
"""
function idxlessthan(x::Array{Float64,1}, val::Float64)

  @inbounds begin

    a::Int64 = 1
    b::Int64 = lastindex(x)

    if x[b] < val
      return b
    end

    mid::Int64 = div(b,2)

    while b-a > 1
      val < x[mid] ? b = mid : a = mid
      mid = div(b + a, 2)
    end

  end

  return a
end

"""
    indmindif_sorted(x::Array{Float64,1}, val::Float64)

Get index in sorted vector `x` corresponding to the value
that is closest to `val` in sorted arrays
using a sort of uniroot algorithm.
"""
function indmindif_sorted(x::Array{Float64,1}, val::Float64)
  a::Int64   = 1
  b::Int64   = lastindex(x)
  mid::Int64 = div(b,2)

  while b-a > 1
    val < x[mid] ? b = mid : a = mid
    mid = div(b + a, 2)
  end

  abs(x[a] - val) < abs(x[b] - val) ? a : b
end

"""
    make_edgeind(childs::Array{Int64,1}, B::Array{Float64,2})

Make ragged array with indexes for each edge in `Y`.
"""
function make_edgeind(childs::Array{Int64,1}, B::Array{Float64,2}, ntip::Int64)

  Bli = LinearIndices(B)
  bridx = UnitRange{Int64}[]
  for b in childs
    bindices = Bli[findall(isequal(b), B)]
    if b != (ntip+1)
      pushfirst!(bindices, bindices[1] - 1)
    end
    bidx = bindices[1]:bindices[end]
    push!(bridx, bidx)
  end

  bridx
end

"""
    make_edgeδt(bridx::Array{Array{Int64,1},1}, δt::Array{Float64,1}, m ::Int64)

Make ragged array of the cumulative δtimes for each branch.
"""
function make_edgeδt(bridx::Array{UnitRange{Int64},1},
                     δt   ::Array{Float64,1},
                     m    ::Int64)

  brδt = Array{Float64,1}[]

  for j in 1:(length(bridx)-1)
    bi = collect(bridx[j][1:(end-1)])
    for i in eachindex(bi)
      bi[i] = rowind(bi[i], m)
    end
    push!(brδt, pushfirst!(cumsum(δt[bi]),0))
  end

  brδt
end

"""
    create_wcol(X::Array{Float64,2})

Returns indices for columns along `m` timesteps.
"""
function create_wcol(X::Array{Float64,2})

  wNaN_x = map(!isnan, X)

  # make ragged array for non-NaN columns
  wcol = Array{Int64,1}[]
  for i = Base.OneTo(size(X,1))
    push!(wcol,findall(wNaN_x[i,:]))
  end

  return wcol
end

"""
    Pc(λi::Float64, λj::Float64, δt::Float64)

Estimate probability of collision.
"""
function Pc(λ1::Float64, λ0::Float64, δt::Float64)
    λt = (λ1 + λ0)*δt
    er = exp(-λt)
    return 1.0 - er - λt*er
end

"""
    collision(λ1::Float64, λ0::Float64, nδts  ::Array{Float64,1})

Estimates the probability of at least one collision
along the whole tree.
"""
function collision(λ1   ::Float64,
                   λ0   ::Float64,
                   nδts  ::Array{Float64,1})

  p = 1.0
  for t in nδts
    p *= 1.0 - Pc(λ1, λ0, t)
  end

  return 1.0 - p
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
"""
    range(; start, stop, length, step)

All keyword form of Base.range
"""
Base.range(;start=1, stop=nothing, length=nothing, step=nothing) =
    Base._range(start, step, stop, length)

Base._range(::Nothing, step, stop, length) = Base.range(stop - (length-1)*step; length=length, step=step)

"""
    range(:bels; b, e, l, s)
    range(:[bels]*, args...)

Specify range properties using the following single letter abbreviations
[b]egin: start
[e]nd: stop
[l]ength: length
[s]step: step

If :bels is the first argument, use the abbreviations as keywords.
For symbols using some permutation of b, e, l, or s of one, two, or Three
members, specify the arguments in the order of the symbols used.

Defaults:
[b]egin: 1 or [`oneunit`](@ref) of given parameter
[e]nd: typemax or floatmax of b
[s]tep: 1
"""
Base.range(s::Symbol, args... ; kwargs...) = Base.range(Val(s), args... ; kwargs...)
### THE ABOVE LINE OF CODE IS CRITICAL FOR THE FUNCTIONALITY BELOW ###

Base.range(::Val{:bels}; b=nothing, e=nothing, l=nothing, s=nothing) =
    Base._range(b, s, e, l)

# Single argument

Base.range(::Val{:b}, b) = b:typemax(b)
Base.range(::Val{:e}, e) = Base.OneTo(e)
Base.range(::Val{:l}, l) = Base.OneTo(l)
Base.range(::Val{:s}, s) = oneunit(s):s:typemax(s)

Base.range(::Val{:b}, b::AbstractFloat) = b:typemax(Int(b))
Base.range(::Val{:s}, s::AbstractFloat) = oneunit(s):s:typemax(Int(s))

# Two argument, bels order

Base.range(::Val{:be}, b, e) = b:e
Base.range(::Val{:bl}, b, l) = range(b; length=l)
Base.range(::Val{:bs}, b, s) = b:s:typemax(s)

Base.range(::Val{:el}, e, l=e) = range(e-l+1, e; length=l)
Base.range(::Val{:es}, e, s=nothing) = oneunit(e):s:e

Base.range(::Val{:ls}, l, s=nothing) = range(1; step=s, length=l)

Base.range(::Val{:bs}, b, s::AbstractFloat) = b:s:floatmax(s)

# Three argument, bels order (+4)

Base.range(::Val{:bel}, b, e, l=e) = Base.range(b, e; length=l)
Base.range(::Val{:bes}, b, e, s=nothing) = Base.range(b, e; step=s)
Base.range(::Val{:bls}, b, l, s=nothing) = Base.range(b; length=l, step=s)
Base.range(::Val{:els}, e, l=e, s=nothing) = Base.range(e - (l-1)*s; length=l, step=s)

# Two argument, reversed order

Base.range(::Val{:eb}, e, b=1) = b:e
Base.range(::Val{:lb}, l, b=1) = range(b; length=l)
Base.range(::Val{:sb}, s, b=1) = b:s:typemax(s)

Base.range(::Val{:le}, l, e=l) = Base.range(e-l+1, e; length=l)
Base.range(::Val{:se}, s, e) = oneunit(e):s:e

Base.range(::Val{:sl}, s, l) = Base.range(1; step=s, length=l)

# Three argument, swap last two (+4)

Base.range(::Val{:ble}, b, l, e=nothing) = Base.range(b, e; length=l)
Base.range(::Val{:bse}, b, s, e) = Base.range(b, e; step=s)
Base.range(::Val{:bsl}, b, s, l) = Base.range(b; length=l, step=s)
Base.range(::Val{:esl}, e, s=nothing, l=e) = Base.range(e - (l-1)*s; length=l, step=s)

# Three argument, end first, ebls (+4)

Base.range(::Val{:ebl}, e, b=1, l=e) = Base.range(b, e; length=l)
Base.range(::Val{:ebs}, e, b=1, s=nothing) = Base.range(b, e; step=s)
#Base.range(::Val{:bls}, b, l, s) = Base.range(b; length=l, step=s)
#Base.range(::Val{:els}, e, l, s) = Base.range(e - (l-1)*s; length=l, step=s)

Base.range(::Val{:elb}, e, l=e, b=1) = Base.range(b, e; length=l)
Base.range(::Val{:esb}, e, s=nothing, b=1) = Base.range(b, e; step=s)
#Base.range(::Val{:bsl}, b, s, l) = Base.range(b; length=l, step=s)
#Base.range(::Val{:esl}, e, s, l) = Base.range(e - (l-1)*s; length=l, step=s)

# Three argument, length first (+6)

Base.range(::Val{:lbe}, l, b=1, e=nothing) = Base.range(b, e; length=l)
#Base.range(::Val{:bes}, b, e, s) = Base.range(b, e; step=s)
Base.range(::Val{:lbs}, l, b=1, s=nothing) = Base.range(b; length=l, step=s)
Base.range(::Val{:les}, l, e=l, s=1) = Base.range(e - (l-1)*s; length=l, step=s)

Base.range(::Val{:leb}, l, e=l, b=1) = Base.range(b, e; length=l)
#Base.range(::Val{:bes}, b, s, e) = Base.range(b, e; step=s)
Base.range(::Val{:lsb}, l, s=nothing, b=1) = Base.range(b; length=l, step=s)
Base.range(::Val{:lse}, l, s=1, e=l) = Base.range(e - (l-1)*s; length=l, step=s)

# Three argument, step first (+6)

#Base.range(::Val{:bel}, b, e, l) = Base.range(b, e; length=l)
Base.range(::Val{:sbe}, s, b, e) = Base.range(b, e; step=s)
Base.range(::Val{:sbl}, s, b, l) = Base.range(b; length=l, step=s)
Base.range(::Val{:sel}, s, e, l=e) = Base.range(e - (l-1)*s; length=l, step=s)

#Base.range(::Val{:bel}, b, l, e) = Base.range(b, e; length=l)
Base.range(::Val{:seb}, s, e, b) = Base.range(b, e; step=s)
Base.range(::Val{:slb}, s, l, b) = Base.range(b; length=l, step=s)
Base.range(::Val{:sle}, s, l, e=l) = Base.range(e - (l-1)*s; length=l, step=s)


text:
## ExoplanetsSysSim/examples/syssim_summary_stats.jl
## (c) 2015 Eric B. Ford

# Calculaties summary statistics based on an input "parameter file" and outputs results to stdout or a file
# usage:  syssim_summary_stats.jl param_file [output_file]
# output: a binary (HDF5/JLD) file containing the summary statistics computed based on the parameter file.

verbose = true

if length(ARGS)<2
   warn("usage:  syssim_summary_stats.jl param_file output_file ")
   exit(255)
end

using Compat
import Compat: UTF8String, ASCIIString

param_file = convert(Compat.ASCIIString,ARGS[1])
if !isfile(param_file)
   warn(string("Can't read ",param_file))
   exit(255)
end

if verbose println("# Loading ExoplanetsSysSim module."); end
using ExoplanetsSysSim

if verbose println("# Reading parameter file ", param_file,"."); end
include(param_file)
if verbose println("# Calling setup_sim_param( ", ARGS[3:end], " )."); end
sim_param = setup_sim_param( convert(Array{Compat.ASCIIString,1},ARGS[3:end]) )
if verbose println("# Active parameter values: ", make_vector_of_sim_param(sim_param) ); end

if haskey(sim_param,"rng_seed")
   seed = get_int(sim_param,"rng_seed")
     if verbose println("# Seeding RNG with seed = ", seed); end
   srand(seed)
end

if verbose println("# Generating physical catalog."); end
cat_phys = generate_kepler_physical_catalog(sim_param)

if verbose println("# Generating observational catalog."); end
cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param)

ss = CatalogSummaryStatistics()
if haskey(sim_param,"use obs data only") &&  convert(Bool,get(sim_param,"use obs data only"))
  if verbose println("# Calculating summary statistics based on all data."); end
  ss = calc_summary_stats_obs_demo(cat_obs,sim_param)
else
  if verbose println("# Calculating summary statistics based on observational data only."); end
  ss = calc_summary_stats_sim_pass_one_demo(cat_obs,cat_phys,sim_param)
  ss = calc_summary_stats_sim_pass_two_demo(cat_obs,cat_phys,ss,sim_param)
end

if length(ARGS)>=2
  output_file = convert(Compat.ASCIIString,ARGS[2])
  if verbose println("# Writing summary statistics to ", output_file,"."); end
  try
     save_sim_results(output_file,sim_param,summary_stats=ss)
  catch
     warn(string("Error writing to ",output_file))
  end
else  # Intetionally don't let you get here, since currently the "summary statistics" are so long, it's obnoxious to print.
      # in principle, could specialize how to print a summary statistics object in a useful way.  (I think the relevant function might be "show".)
  if verbose println("# Summary statistics:"); end
  println(ss)
end


text:
#! /usr/bin/julia

# Rosetta Code, Find the missing permutation

function find_missing_permutations{T<:String}(a::Array{T,1})
    std = unique(sort(split(a[1], "")))
    needsperm = trues(factorial(length(std)))
    for s in a
        b = split(s, "")
        p = map(x->findfirst(std, x), b)
        isperm(p) || throw(DomainError())
        needsperm[nthperm(p)] = false
    end
    mperms = T[]
    for i in findn(needsperm)[1]
        push!(mperms, join(nthperm(std, i), ""))
    end
    return mperms
end

test = ["ABCD", "CABD", "ACDB", "DACB", "BCDA", "ACBD",
        "ADCB", "CDAB", "DABC", "BCAD", "CADB", "CDBA",
        "CBAD", "ABDC", "ADBC", "BDCA", "DCBA", "BACD",
        "BADC", "BDAC", "CBDA", "DBCA", "DCAB"]

missperms = find_missing_permutations(test)

print("The test list is:\n    ")
i = 0
for s in test
    print(s, " ")
    i += 1
    i %= 5
    i != 0 || print("\n    ")
end
i == 0 || println()
if length(missperms) > 0
    println("The following permutations are missing:")
    for s in missperms
        println("    ", s)
    end
else
    println("There are no missing permutations.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# MathOptInterface NLP

include("reshop_nlp_moi.jl")

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    # Prevent this function getting called twice
    # This ay happen with EMP
    if (model.len_nl_cons > 0) return end
    MOI.initialize(nlp_data.evaluator, [:ExprGraph])

    # Process constraints
    load_nlp_constraints(model, nlp_data)

    if nlp_data.has_objective
        load_nlp_objective(model, nlp_data)
    end

end

# Keep loading of NLP constraints apart to load all NLP model all in once
# inside ReSHOP.
function load_nlp_constraints(model::Optimizer, nlp_data::MOI.NLPBlockData)
    # This heavily relies on some convention
    model.start_nl_cons = Int(ctx_numequ(model.ctx))
    model.len_nl_cons = length(nlp_data.constraint_bounds)
    for i=1:length(nlp_data.constraint_bounds)
        c = MOI.constraint_expr(nlp_data.evaluator, i)

        ub = +Inf
        lb = -Inf

        # Remove relations and bounds from constraint expressions
        if length(c.args) == 3
            expected_head = :call
            expr_index = 2

            @assert c.head == expected_head
            # Single relation constraint: expr rel bound
            rel = c.args[1]
            if rel in [:<=, :(==)]
                ub = c.args[3]
            end
            if rel in [:>=, :(==)]
                lb = c.args[3]
            end
            c = c.args[expr_index]
            @assert isfinite(lb) || isfinite(ub)
        else
            # Double relation constraint: bound <= expr <= bound
            @assert c.head == :comparison
            error("doubly bounded NLP constraint is not yet supported")
        end

        # Convert non-linear expression to non-linear, linear and constant
        lin_part = Dict{Int32, Float64}()
        c, constant, conlinearities = process_expression!(c, lin_part)

        # Update bounds on constraint
        lb -= constant
        ub -= constant

        eidx = rhp_add_equ(model.ctx)

        tree, node = reshop_get_treedata(model.ctx, eidx)
        reshop_add_nlexpr(model.ctx, tree, node, c)
        lin_part = filter( x -> (abs(last(x)) > eps(0.)), lin_part)
        if length(lin_part) > 0
            avar = _ensure_avar(model)
            rhp_avar_set(avar, Vector{Int32}(keys(lin_part) .- 1))
            rhp_equ_add_linear_chk(model.ctx, eidx, avar, collect(values(lin_part)))
        end
        reshop_set_equtype(model.ctx, eidx, relation_to_reshop[rel])
        isfinite(lb) && reshop_set_rhs(model.ctx, eidx, lb)
        isfinite(ub) && reshop_set_rhs(model.ctx, eidx, ub)
    end

 end

function load_nlp_objective(model::Optimizer, nlp_data::MOI.NLPBlockData)
# Process objective
    obj = MOI.objective_expr(nlp_data.evaluator)

    if length(obj.args) < 2
        obj = 0
    else
        # Convert non-linear expression to non-linear, linear and constant
        lin_part = Dict{Int32, Float64}()
        obj, constant, objlinearity = process_expression!(obj, lin_part)

        # Add constant back into non-linear expression
        if constant != 0
            obj = add_constant(obj, constant)
        end
    end

    eidx = rhp_add_equ(model.ctx)
    tree, node = reshop_get_treedata(model.ctx, eidx)
    reshop_add_nlexpr(model.ctx, tree, node, obj)
    lin_part = filter( x -> (abs(last(x)) > eps(0.)), lin_part)
    if length(lin_part) > 0
       avar = _ensure_avar(model)
       rhp_avar_set(avar, Vector{Int32}(keys(lin_part) .- 1))
       rhp_equ_add_linear_chk(model.ctx, eidx, avar, collect(values(lin_part)))
    end
    rhp_set_objeqn(model.ctx, eidx)
end


text:
export optimize

function optimize(C::SparseRadialMapComponent, X, poff::Union{Nothing, Int64}, pdiag::Union{Int64, Nothing}, maxfamiliesoff::Union{Nothing, Int64, String},
                  ; λ::Float64 = 0.0, δ::Float64 = 1e-8, γ::Float64 = 2.0,
                  maxpatience::Int64 = 10^5, verbose::Bool = false, maxfeatures::Int64 = 10^4)
    Nx = C.Nx
    # By default the diagonal component is selected
    if typeof(maxfamiliesoff) <: Nothing
        C, error = optimize(C, X, nothing, λ, δ)
    elseif typeof(maxfamiliesoff) <: Int64
        C, error =  greedyfit(Nx, poff, pdiag, X, maxfamiliesoff, λ, δ, γ;
                              verbose = verbose)

    elseif maxfamiliesoff ∈ ("kfold", "kfolds", "Kfold", "Kfolds")
        # Define cross-validation splits of data
        n_folds = 5
        folds = kfolds(1:size(X,2), k = n_folds)

        # Run greedy approximation
        if pdiag == 0
            maxfamily = min(Nx-1, ceil(Int64, (sqrt(size(X,2))-(pdiag+1))/(poff+1)), maxfeatures)
        elseif pdiag > 0
            maxfamily = min(Nx-1, ceil(Int64, (sqrt(size(X,2))-(pdiag+3))/(poff+1)), maxfeatures)
        else
            error("Wrong value for pdiag")
        end

        valid_error = zeros(maxfamily+1, n_folds)
        # if typeof(P) <: Serial
        @inbounds for i=1:n_folds
            if verbose == true
                println("Fold "*string(i)*":")
            end
            idx_train, idx_valid = folds[i]

            C, error = greedyfit(Nx, poff, pdiag, X[:,idx_train], X[:,idx_valid], maxfamily,
                                 λ, δ, γ;
                                 maxpatience = maxpatience, verbose  = verbose)

            # error[2] contains the history of the validation error
            valid_error[:,i] .= copy(error[2])
        end
        # elseif typeof(P) <: Thread
          # @inbounds   Threads.@threads for i=1:n_folds
          #     if verbose == true
          #         println("Fold "*string(i)*":")
          #     end
          #     idx_train, idx_valid = folds[i]
          #
          #     C, error = greedyfit(Nx, poff, pdiag, X[:,idx_train], X[:,idx_valid], maxfamily,
          #                          λ, δ, γ;
          #                          maxpatience = maxpatience, verbose  = verbose)
          #
          #     # error[2] contains the history of the validation error
          #     valid_error[:,i] .= copy(error[2])
          # end
        # end

        # Find optimal numbers of terms
        mean_valid_error = mean(valid_error, dims  = 2)[:,1]
        _, opt_families = findmin(mean_valid_error)
        # opt_families corresponds to the number of off-diagonal components
        opt_families -= 1

        if verbose == true
            println("Optimization on the entire data set:")
        end
        # Run greedy fit up to opt_nterms on all the data
        C, error = greedyfit(Nx, poff, pdiag, X, opt_families, λ, δ, γ; verbose = verbose)

    elseif maxfamiliesoff ∈ ("split", "Split")
        nvalid = ceil(Int64, floor(0.2*size(X,2)))
        X_train = X[:,nvalid+1:end]
        X_valid = X[:,1:nvalid]

        # Run greedy approximation
		if pdiag == 0
			maxfamily = min(Nx-1, ceil(Int64, (sqrt(size(X,2))-(pdiag+1))/(poff+1)), maxfeatures)
		elseif pdiag > 0
			maxfamily = min(Nx-1, ceil(Int64, (sqrt(size(X,2))-(pdiag+3))/(poff+1)), maxfeatures)
		else
			error("Wrong value for pdiag")
		end

        C, error = greedyfit(Nx, poff, pdiag, X_train, X_valid, maxfamily, λ, δ, γ;
                             maxpatience = maxpatience, verbose  = verbose)
    else
        println("Argument max_terms is not recognized")
        error()
    end
    return C, error
end

function optimize(S::SparseRadialMap, X::AbstractMatrix{Float64}, poff::Union{Int64, Nothing}, pdiag::Union{Int64, Vector{Int64}, Nothing},
	              maxfamilies::Union{Int64, Nothing, String}; apply_rescaling::Bool=true, start::Int64=1,maxpatience::Int64=10^5,
	              verbose::Bool=false, P::Parallel=serial, maxfeatures::Int64 = 10^4)

	NxX, Ne = size(X)
	@get S (Nx, p, γ, λ, δ, κ)

	@assert NxX == Nx "Wrong dimension of the ensemble matrix X"

	# We can apply the rescaling to all the components once
	if apply_rescaling == true
		transform!(S.L, X)
	end

	# Compute centers and widths
	center_std!(S, X)
	# Optimize coefficients
	# Skip the identity components of the map
	if typeof(P)==Serial
		@inbounds for i=start:Nx
			if typeof(maxfamilies) <: Nothing && S.C[i].p[end] != -1
				S.C[i], _ = optimize(S.C[i], X[1:i,:], nothing, nothing, maxfamilies; λ = λ, δ = δ, γ = γ,
									 maxpatience = maxpatience, verbose = verbose)
				copy!(S.p[i], S.C[i].p)
			else
				if typeof(pdiag) <: Int64 && pdiag != -1
						S.C[i], _ = optimize(S.C[i], X[1:i,:], poff, pdiag, maxfamilies; λ = λ, δ = δ, γ = γ,
						                     maxpatience = maxpatience, verbose = verbose, maxfeatures = maxfeatures)
						copy!(S.p[i], S.C[i].p)
				elseif typeof(pdiag) <: Vector{Int64} && pdiag[i-start+1] != -1
						S.C[i], _ = optimize(S.C[i], X[1:i,:], poff, pdiag[i-start+1], maxfamilies; λ = λ, δ = δ, γ = γ,
											 maxpatience = maxpatience, verbose = verbose, maxfeatures = maxfeatures)
						copy!(S.p[i], S.C[i].p)
				end
			end
		end
	else
		@inbounds for i=start:final
			if typeof(maxfamilies) <: Nothing
				S.C[i], _ = optimize(S.C[i], X[1:i,:], nothing, nothing, maxfamilies; λ = λ, δ = δ, γ = γ,
									 maxpatience = maxpatience, verbose = verbose)
				copy!(S.p[i], S.C[i].p)
			else
				if typeof(pdiag) <: Int64 && pdiag != -1
						S.C[i], _ = optimize(S.C[i], X[1:i,:], poff, pdiag, maxfamilies; λ = λ, δ = δ, γ = γ,
											 maxpatience = maxpatience, verbose = verbose, maxfeatures = maxfeatures)
						copy!(S.p[i], S.C[i].p)
				elseif typeof(pdiag) <: Vector{Int64} && pdiag[i-start+1] != -1
						S.C[i], _ = optimize(S.C[i], X[1:i,:], poff, pdiag[i-start+1], maxfamilies; λ = λ, δ = δ, γ = γ,
											 maxpatience = maxpatience, verbose = verbose, maxfeatures = maxfeatures)
						copy!(S.p[i], S.C[i].p)
				end
			end
		end
	end

	# We can apply the rescaling to all the components once
	if apply_rescaling == true
		itransform!(S.L, X)
	end

	return S
end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fb56c064-b90d-11eb-124f-bd2478c9a5bf
using WAV, SampledSignals, Plots,PlutoUI

# ╔═╡ 56d13e40-49bc-4f45-9033-0a4003070da8
using UnicodePlots

# ╔═╡ a6fcf1e7-e477-41d3-b33a-08f521db740c
@bind a Slider(1:10)

# ╔═╡ 6813d955-d54d-47ac-8cec-5be4ea1847fe
begin
# a =  1   # Amplitude
f = 1# Frequency # in hz # No of oscillations of medium particles per unit time
t = 1/f   # Time Period
l = 100   # WaveLength
w = 2*pi*f # omega
v = f * l # Velocity
k = 1 # 2* pi / l # Wave number
padiff = 10
phdiff = k * padiff # phase difference
padiff = phdiff / k # path difference
phdiff = 10
m = 100 # mass I think
# e = 1/2 * (m * (w * a)^2) # energy of wave
fs = 8000 # sample rate
i = 0.0:1/fs:prevfloat(1.0)
alpha = pi # alpha

end

# ╔═╡ c5413acd-7828-4da0-9009-65fd61264637
begin

# some examples from WaveForm.jl
squarewave(x::Real) = ifelse(mod2pi(x) < π, 1.0, -1.0)

#Compute ``2\pi``-periodic square wave of `x` with a duty cycle `θ`and a peak amplitude ``1``.

function squarewave(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod2pi(x) < 2π * θ, 1.0, -1.0)
end

#Compute ``1``-periodic square wave of `x` with a peak amplitude ``1``.
squarewave1(x::Real) = ifelse(mod(x, 1) < 1/2, 1.0, -1.0)

# Compute ``1``-periodic square wave of `x` with a duty cycle `θ` and a peak amplitude ``1``.
function squarewave1(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave1(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod(x, 1) < θ, 1.0, -1.0)
end

#Compute ``2\pi``-periodic triangle wave of `x` with a peak amplitude ``1``.
function trianglewave(x::Real)
    modx = mod2pi(x + π/2)
    ifelse(modx < π, 2modx/π - 1, -2modx/π + 3)
end

#Compute ``1``-periodic triangle wave of `x` with a peak amplitude ``1``.
function trianglewave1(x::Real)
    modx = mod(x + 1/4, 1.0)
    ifelse(modx < 1/2, 4modx - 1, -4modx + 3)
end
#Compute ``2\pi``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
sawtoothwave(x::Real) = rem2pi(x, RoundNearest) / π
#Compute ``1``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
sawtoothwave1(x::Real) = rem(x, 1.0, RoundNearest) * 2
end

# ╔═╡ df0fbb6c-bad6-435d-a931-2a4cd9ea935c
begin

vpart = a*w*cos.(w*i .- k*phdiff .+alpha) #velocity of a particle
apart = -a*(w^2)*sin.(w*i .- k*phdiff .+alpha) #velocity of a particle
# y = a*sin(w*t-k*phdiff +alpha) # general equation
# y1 = a*sin(w*(i - phdiff/v)) # 1st Equation
# y2 = a*sin((2*pi/l)*(v*i-phdiff)) # 2nd Equation
# y3 = a*sin((2*pi)*(i/t - phdiff/l)) # different forms of same equation

wave = a*sin.(w*i .- k*phdiff .+alpha)
wave1 = a*sin.(w*i .- k*phdiff .+alpha) + a*sin.(w*i .- k*phdiff .+alpha .+pi )
wave2 = a*sin.(w*i .- k*phdiff .+alpha) + a*sin.(w*i .- k*phdiff .+alpha .+pi ) .+ a*squarewave.(w*i .- k*phdiff .+alpha .+pi )
wave3 = a*squarewave.(w*i .- k*phdiff .+alpha .+pi )
# Wave Source # Lossless medium
# fs = 8e3
# t = 0.0:1/fs:prevfloat(1.0)
# f = 1e3
y = sin.(2pi * f * i) * 0.1
wavplay(y, fs)
end

# ╔═╡ 683fa835-7339-4b7c-89fe-6511da652e2d
length(wave)

# ╔═╡ 17af4d9f-687a-48c8-9a68-10d9bb713187
lineplot(i, wave)

# ╔═╡ 263973d5-9d5c-427b-94cb-bd74073a7a29
lineplot(i, wave1)

# ╔═╡ 8c401e62-30b2-492d-a896-b0dceb2c8508
lineplot(i, wave2)

# ╔═╡ aa1d8f3a-09cd-4fee-b6a0-e959650b0faa
lineplot(i, wave3)

# ╔═╡ a5298c48-a277-4c1c-bb35-e129cb205df2
lineplot(i, y)

# ╔═╡ e4b55da5-2cba-46b0-a084-2ebbb93a31f9
# plot(i, y, seriestype = :scatter, title = "My Scatter Plot")

# ╔═╡ 5ffbc4bf-cc16-474e-8985-2ebb81dc821c
# plot(i, wave, seriestype = :scatter, title = "My Scatter Plot")

# ╔═╡ 3ee38394-fce7-456f-86f8-0f889b2c4f1a
# plot(i, wave1, seriestype = :scatter, title = "My Scatter Plot")

# ╔═╡ 37c49652-2500-4236-8ffb-0906ff0248e5
# plot(i, wave2, seriestype = :scatter, title = "My Scatter Plot")

# ╔═╡ 1df63790-0d76-4c41-aaf9-5cf39ec6cb98
@bind fs1 Slider(100:100:10000,default=8e3, show_value=true)

# ╔═╡ 446e0ba3-4c3a-4e15-b1a1-5a8046c02c3c
@bind f1 Slider(0:100:10000, default=1e3, show_value=true)

# ╔═╡ dd194f36-68c0-440b-86f8-9f3ec2a78647
begin
t1 = 0.0:1/fs1:prevfloat(1.0)
#y1 = squarewave.(f1*t1) * 1
#y1 = (sin.(2pi * f1* t1) .*0.1)
#y1 = trianglewave.(f1* t1) * 1
y1 = (sin.(2pi * f1* t1) .*t1) .+ (sin.(2pi * f1*2* t1 .+pi) .*t1/2)
wavplay(y1, fs1)
end

# ╔═╡ b9a462bd-67b0-426a-a912-9422c3524885
lineplot(t1, y1)

# ╔═╡ 545eb713-70f5-4e04-adf5-bbb6100a453f
@bind a1 Slider(0:0.1:2,default=0.1, show_value=true)

# ╔═╡ 8577ca66-db1e-4f9b-a775-adb4b404e4e8
@bind x Slider(0:100:10000,default=100, show_value=true)

# ╔═╡ 2d064c8f-d7f4-48ed-985d-c61172088a62
@bind L Slider(0:100:10000,default=100, show_value=true)

# ╔═╡ ee1eaf9e-1459-4184-85e6-0736053304ef
@bind fw Slider(0:100:10000,default=100, show_value=true)

# ╔═╡ ed916679-4bca-41cc-9fae-c6b8362855a0
@bind n Slider(1:1:10,default=2, show_value=true)

# ╔═╡ 236f96dd-8ed0-4332-9e78-41bd0b57c5ac
@bind vibfs Slider(0:100:10000,default=100, show_value=true)

# ╔═╡ 58f1d291-6cd3-42a2-b48a-44a4fd590dd4
# making vibrations
begin
tn = 0.0:1/vibfs:prevfloat(1.0)
vib =  2 * a1 * sin.(n*pi*x/L).*cos.(2*pi*fw*tn)
end

# ╔═╡ d8396378-e8b5-413f-87b4-d77479e130e7
lineplot(tn, vib)

# ╔═╡ d8d86fdc-8adf-43db-ae77-a04dabb2c63a
# wavplay(vib, vibfs)

# ╔═╡ fd76e076-2b61-4e55-9cd3-e59340168f4d

# ╔═╡ Cell order:
# ╠═fb56c064-b90d-11eb-124f-bd2478c9a5bf
# ╠═a6fcf1e7-e477-41d3-b33a-08f521db740c
# ╠═6813d955-d54d-47ac-8cec-5be4ea1847fe
# ╠═df0fbb6c-bad6-435d-a931-2a4cd9ea935c
# ╠═c5413acd-7828-4da0-9009-65fd61264637
# ╠═56d13e40-49bc-4f45-9033-0a4003070da8
# ╠═683fa835-7339-4b7c-89fe-6511da652e2d
# ╠═17af4d9f-687a-48c8-9a68-10d9bb713187
# ╠═263973d5-9d5c-427b-94cb-bd74073a7a29
# ╠═8c401e62-30b2-492d-a896-b0dceb2c8508
# ╠═aa1d8f3a-09cd-4fee-b6a0-e959650b0faa
# ╠═a5298c48-a277-4c1c-bb35-e129cb205df2
# ╠═e4b55da5-2cba-46b0-a084-2ebbb93a31f9
# ╠═5ffbc4bf-cc16-474e-8985-2ebb81dc821c
# ╠═3ee38394-fce7-456f-86f8-0f889b2c4f1a
# ╠═37c49652-2500-4236-8ffb-0906ff0248e5
# ╠═1df63790-0d76-4c41-aaf9-5cf39ec6cb98
# ╠═446e0ba3-4c3a-4e15-b1a1-5a8046c02c3c
# ╠═dd194f36-68c0-440b-86f8-9f3ec2a78647
# ╠═b9a462bd-67b0-426a-a912-9422c3524885
# ╠═545eb713-70f5-4e04-adf5-bbb6100a453f
# ╠═8577ca66-db1e-4f9b-a775-adb4b404e4e8
# ╠═2d064c8f-d7f4-48ed-985d-c61172088a62
# ╠═ee1eaf9e-1459-4184-85e6-0736053304ef
# ╠═ed916679-4bca-41cc-9fae-c6b8362855a0
# ╠═236f96dd-8ed0-4332-9e78-41bd0b57c5ac
# ╠═58f1d291-6cd3-42a2-b48a-44a4fd590dd4
# ╠═d8396378-e8b5-413f-87b4-d77479e130e7
# ╠═d8d86fdc-8adf-43db-ae77-a04dabb2c63a
# ╠═fd76e076-2b61-4e55-9cd3-e59340168f4d

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
"""
    _abstractify_grouping(f::FormulaTerm)

Remove concrete levels associated with a schematized FormulaTerm.

Returns the formula with the grouping variables made abstract again
and a Dictionary of `Grouping()` contrasts.
"""
function _abstractify_grouping(f::FormulaTerm)
    fe = filter(x -> !isa(x, AbstractReTerm), f.rhs)
    re = filter(x -> isa(x, AbstractReTerm), f.rhs)
    contr = Dict{Symbol,AbstractContrasts}()
    re = map(re) do trm
        if trm.rhs isa InteractionTerm
            rhs = mapreduce(&, trm.rhs.terms) do tt
                # how to define Grouping() for interactions on the RHS?
                # contr[tt.sym] = Grouping()
                return Term(tt.sym)
            end
        else
            contr[trm.rhs.sym] = Grouping()
            rhs = Term(trm.rhs.sym)
        end
        return trm.lhs | rhs
    end
    return (f.lhs ~ sum(fe) + sum(re)), contr
end

"""
    isconstant(x::Array)
    isconstant(x::Tuple)

Are all elements of the iterator the same?  That is, is it constant?
"""
function isconstant(x; comparison=isequal)::Bool
    # the ref is necessary in case the elements of x are themselves arrays
    return isempty(x) ||
           all(ismissing, x) ||
           coalesce(all(comparison.(x, Ref(first(x)))), false)
end

isconstant(x::Vector{Bool})::Bool = !any(x) || all(x)

isconstant(x...; comparison=isequal) = isconstant(x; comparison=comparison)

"""
    average(a::T, b::T) where {T<:AbstractFloat}

Return the average of `a` and `b`
"""
average(a::T, b::T) where {T<:AbstractFloat} = (a + b) / T(2)

"""
    cpad(s::AbstractString, n::Integer)

Return a string of length `n` containing `s` in the center (more-or-less).
"""
cpad(s::String, n::Integer) = rpad(lpad(s, (n + textwidth(s)) >> 1), n)

"""
densify(S::SparseMatrix, threshold=0.1)

Convert sparse `S` to `Diagonal` if `S` is diagonal or to `Array(S)` if
the proportion of nonzeros exceeds `threshold`.
"""
function densify(A::SparseMatrixCSC, threshold::Real=0.1)
    m, n = size(A)
    if m == n && isdiag(A)  # convert diagonal sparse to Diagonal
        # the diagonal is always dense (otherwise rank deficit)
        # so make sure it's stored as such
        Diagonal(Vector(diag(A)))
    elseif nnz(A) / (m * n) ≤ threshold
        A
    else
        Array(A)
    end
end
densify(A::AbstractMatrix, threshold::Real=0.1) = A

densify(A::SparseVector, threshold::Real=0.1) = Vector(A)
function densify(A::Diagonal{T,SparseVector{T,Ti}}, threshold::Real=0.1) where {T,Ti}
    return Diagonal(Vector(A.diag))
end

"""
    RaggedArray{T,I}

A "ragged" array structure consisting of values and indices

# Fields
- `vals`: a `Vector{T}` containing the values
- `inds`: a `Vector{I}` containing the indices

For this application a `RaggedArray` is used only in its `sum!` method.
"""
struct RaggedArray{T,I}
    vals::Vector{T}
    inds::Vector{I}
end

function Base.sum!(s::AbstractVector{T}, a::RaggedArray{T}) where {T}
    for (v, i) in zip(a.vals, a.inds)
        s[i] += v
    end
    return s
end

function rownormalize(A::AbstractMatrix)
    A = copy(A)
    for r in eachrow(A)
        # all zeros arise in zerocorr situations
        if !iszero(r)
            normalize!(r)
        end
    end
    return A
end

function rownormalize(A::LowerTriangular{T,Diagonal{T,Vector{T}}}) where {T}
    return one(T) * I(size(A, 1))
end

# from the ProgressMeter docs
_is_logging(io) = isa(io, Base.TTY) == false || (get(ENV, "CI", nothing) == "true")

"""
    replicate(f::Function, n::Integer; use_threads=false)

Return a vector of the values of `n` calls to `f()` - used in simulations where the value of `f` is stochastic.

`hide_progress` can be used to disable the progress bar. Note that the progress
bar is automatically disabled for non-interactive (i.e. logging) contexts.

!!! warning
    If `f()` is not thread-safe or depends on a non thread-safe RNG,
    then you must set `use_threads=false`. Also note that ordering of replications
    is not guaranteed when `use_threads=true`, although the replications are not
    otherwise affected for thread-safe `f()`.
"""
function replicate(f::Function, n::Integer; use_threads=false, hide_progress=false)
    # no macro version yet: https://github.com/timholy/ProgressMeter.jl/issues/143
    # and we want some advanced options
    p = Progress(n; output=Base.stderr, enabled=!hide_progress && !_is_logging(stderr))
    # get the type
    rr = f()
    next!(p)
    # pre-allocate
    results = [rr for _ in Base.OneTo(n)]
    if use_threads
        Threads.@threads for idx in 2:n
            results[idx] = f()
            next!(p)
        end
    else
        for idx in 2:n
            results[idx] = f()
            next!(p)
        end
    end
    return results
end

"""
    sdcorr(A::AbstractMatrix{T}) where {T}

Transform a square matrix `A` with positive diagonals into an `NTuple{size(A,1), T}` of
standard deviations and a tuple of correlations.

`A` is assumed to be symmetric and only the lower triangle is used.  The order of the
correlations is row-major ordering of the lower triangle (or, equivalently, column-major
in the upper triangle).
"""
function sdcorr(A::AbstractMatrix{T}) where {T}
    m, n = size(A)
    m == n || throw(ArgumentError("matrix A must be square"))
    indpairs = checkindprsk(m)
    rtdiag = sqrt.(NTuple{m,T}(diag(A)))
    return (
        rtdiag,
        ntuple(kchoose2(m)) do k
            i, j = indpairs[k]
            A[i, j] / (rtdiag[i] * rtdiag[j])
        end,
    )
end


text:
# Uses new formulation from limbdark paper.
using Limbdark
import Limbdark: Transit_Struct, s2

#include("s2_stable.jl")
#include("../../../src/s2.jl")
#include("../../../src/define_constants.jl")

# Run some tests:

using PyPlot
nb = 1001; nr = 1001
b=range(0.0,stop=2,length=nr)
r=range(0.0,stop=2,length=nr)
fgrid = zeros(Float64,nr,nb)
fgrid_old = zeros(Float64,nr,nb)
fgrid_big = zeros(Float64,nr,nb)
for ib=1:nb
#  if mod(ib,10) == 0
#    println("Finished: ",ib/nr*100,"%")
#  end
  for ir=1:nr
    fgrid[ir,ib]=s2(r[ir],b[ib])
    fgrid_big[ir,ib]=convert(Float64,s2(big(r[ir]),big(b[ib])))
  end
end

clf()
figure(figsize=(5, 4))
img=imshow(1e15 * (fgrid-fgrid_big),interpolation="nearest",origin="lower",
           extent=[0,2,0,2],cmap="plasma",vmin=-1.5, vmax=1.5)
plot([0,1],[1,2],"k-")
plot([1,2],[0,1],"k-")
plot([0,0.38],[0,0.38],"k-")
plot([0.5,2],[0.5,2],"k-")
plot([0,1],[1,0],"k-")
annotate(L"$r = 1 + b$", xy=(0.5, 1.5), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="k", rotation=45)

annotate(L"$r = b$", xy=(1, 1), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="k", rotation=45)

annotate(L"$r = b - 1$", xy=(1.5, 0.5), xytext=(0, 0), xycoords="data",
         textcoords="offset points", ha="center", va="bottom",
         fontsize=12, color="k", rotation=45)

annotate(L"$r = 1 - b$", xy=(0.5, 0.5), xytext=(-9, -9), xycoords="data",
         textcoords="offset points", ha="center", va="center",
         fontsize=12, color="k", rotation=-45)

xlabel("Impact parameter")
ylabel("Radius ratio")
ticks = [-1.5, -1, -0.5, 0, 0.5, 1, 1.5]
cbar = colorbar(img, ticks=ticks)
cbar[:ax][:set_yticklabels]([L"$-1.5\times 10^{-15}$",
                             L"$-1.0\times 10^{-15}$",
                             L"$-0.5\times 10^{-15}$",
                             L"$0$",
                             L"$0.5\times 10^{-15}$",
                             L"$1.0\times 10^{-15}$",
                             L"$1.5\times 10^{-15}$"])

savefig("s2_residuals.pdf")

max_issues_repo_name: rodluger/limbdark


text:
#Define a very simple TargetModel for all tests below
t1 = linspace(0.0,10.0,100)
p1 = parameters([:a],[1.0],[5.0],[3.0])
d1 = data(:array,t1,sin(3*t1))
n1 = noise(:gaussian,[0.01])
m1 = model(:target,p1,d1,n1,(t,p)->sin(p[1]*t);name="JobSegmentsTestModel")

#Define an MHSampler with normal proposal distribution
s1 = sampler(:mh,:normal,0.1,eye(1))

###################Test a runner with 1 proposal per iteration (normal Metropolis-Hastings)
println("==================")
println("Testing GMHSegment")
println("==================")
nprops1 = 2
nparas1 = length(p1)
policy1 = policy(:mh,nprops1)
indicatorstate1 = samplerstate(s1,1,policy1.sampletype,policy1.calculationtype)
auxiliarystate1 = samplerstate(s1,nprops1,policy1.sampletype,policy1.calculationtype)
copy!(indicatorstate1.proposals.values,[3.0])

segment1 = segment(policy1,m1,s1,nprops1)

@test numproposals(segment1) == 2
@test segment1.samplerstate.from.values == [0.0]
srand(234) ; prepareauxiliary!(indicatorstate1,auxiliarystate1) ; propose!(auxiliarystate1) ; geometry!(m1,proposals(auxiliarystate1)) ; a1 = acceptance!(auxiliarystate1)
srand(234) ; a2 = iterate!(segment1,indicatorstate1)
@test segment1.samplerstate.from.values == [3.0]
@test a1 == a2

@test getsamples(segment1,[2,1,2]) == copy(proposals(auxiliarystate1),[2,1,2])
@test getsamples(segment1,2) == copy(proposals(auxiliarystate1),2)
@test getsamples(segment1,[]) == samples(:base,1,0,Float64,Float64)

@test_approx_eq cov(segment1.samplerstate.density.distribution) cov(auxiliarystate1.density.distribution)
@test_approx_eq cov(segment1.samplerstate.density.distribution) 0.01*eye(1)
tune!(segment1,2.0)
tune!(auxiliarystate1,2.0)
@test_approx_eq cov(segment1.samplerstate.density.distribution) cov(auxiliarystate1.density.distribution)
@test_approx_eq cov(segment1.samplerstate.density.distribution) 0.04*eye(1)

println()
println("====================")
println("Test show() function")
show(segment1)
println("End  show() function")
println("====================")
println()

println("======================")
println("Testing RemoteSegments")
println("======================")

#test of number of job segments
@test GeneralizedMetropolisHastings._numjobsegments(Val{:procs}) == nprocs()
@test GeneralizedMetropolisHastings._numjobsegments(Val{:workers}) == nworkers()
@test GeneralizedMetropolisHastings._numjobsegments(Val{:test}) == 3

#test the process numbers on which job segments will run
@test GeneralizedMetropolisHastings._processnumbers(Val{:procs}) == procs()
@test GeneralizedMetropolisHastings._processnumbers(Val{:workers}) == workers()
@test GeneralizedMetropolisHastings._processnumbers(Val{:test}) == workers()

#test how many proposals per segment
nproposals = [2,3,4,5,6,7]
nsegments = [2,3,3,3,3,3]
npropsperseg = [1,1,2,2,2,3]
for i=1:length(nproposals)
    @test GeneralizedMetropolisHastings._numproposalspersegment(nproposals[i],nsegments[i]) == npropsperseg[i]
end

function gettestprocessnumbers(p::Int)
    if p == 1
        return [1,1,1]
    end
    if p == 2
        w = workers()[1]
        return [w,w,w]
    end
    if p == 3
        w1 = workers()[1]
        w2 = workers()[2]
        return [w1,w2,w1]
    end
    return workers()[1:3] #if there are many workers
end

rprops1 = 4
rpolicy1 = policy(:mh,nprops1,jobsegments=:test)
rsegments1 = remotesegments(rpolicy1,m1,s1,rprops1)
@test numsegments(rsegments1) == 3
@test numproposalspersegment(rsegments1) == 2
@test numtotalproposals(rsegments1) == 6
@test GeneralizedMetropolisHastings._numjobsegments(rpolicy1,rprops1) == 3
@test collect(GeneralizedMetropolisHastings._processnumbers(rpolicy1,3)) == gettestprocessnumbers(nprocs())
@test GeneralizedMetropolisHastings._insegmentindex(rsegments1,[1,4,5,6]) == Array{Int,1}[[1],[2],[1,2]] #create the in-segment index for overall proposal index
@test GeneralizedMetropolisHastings._insegmentindex(rsegments1,[1,4,5,6,1,3,4,5,5]) == Array{Int,1}[[1],[1,2],[1,2]] #repetitions are filtered out
@test rsegments1.prop2collected == Dict{Int,Tuple{Int,Int}}()
GeneralizedMetropolisHastings._insegmentindex2collected!(rsegments1,GeneralizedMetropolisHastings._insegmentindex(rsegments1,[1,4,5,6]))
@test rsegments1.prop2collected == Dict(1=>(1,1),4=>(1,2),5=>(1,3),6=>(2,3))

#prepare for following tests
indicatorstate1 = samplerstate(s1,1,policy1.sampletype,policy1.calculationtype)
initialize!(trait(:initialize,:default),m1,indicatorstate1.proposals)
geometry!(m1,indicatorstate1.proposals)
lsegments1 = [segment(policy1,m1,s1,numproposalspersegment(rsegments1)) for i=1:numsegments(rsegments1)]

#test the iterate function
srand(435) ; a1 = iterate!(rsegments1,indicatorstate1)
srand(435) ; a2 = map((s)->iterate!(s,indicatorstate1),lsegments1)
@test nprocs() > 1 || a1 == a2 #only test this if all segments run on the same process

#test the prepare! function which copies a sample back into the indicator state
prepare!(rsegments1,indicatorstate1,2)
@test  nprocs() > 1 || indicatorstate1.from == copy(lsegments1[1].samplerstate.proposals,2)
prepare!(rsegments1,indicatorstate1,6)
@test  nprocs() > 1 || indicatorstate1.from == copy(lsegments1[3].samplerstate.proposals,2)

#test the retrievesamples! function which accesses the remote segments to retrieve samples from
samples1 = retrievesamples!(rsegments1,[1,4,5,6])
samples2 = map(getsamples,lsegments1,Array{Int,1}[[1],[2],[1,2]])
@test  nprocs() > 1 || samples1 == samples2

#test the getsamples function for an index previously retrieved
@test  nprocs() > 1 || getsamples(rsegments1,4) == getsamples(lsegments1[2],2)

#test the getsamples function for an index not previously retrieved
@test  nprocs() > 1 || getsamples(rsegments1,3) == getsamples(lsegments1[2],1)

#test storing the smaples into the chain
c1 = chain(:standard,1,3,policy1.sampletype,policy1.calculationtype)
c2 = chain(:standard,1,3,policy1.sampletype,policy1.calculationtype)
store!(c1,indicatorstate1.proposals,1)
store!(c2,indicatorstate1.proposals,1)
store!(rsegments1,c1,1)
store!(c2,lsegments1[1].samplerstate.proposals,1)
@test  nprocs() > 1 || c1.values == c2.values

#test tuning the segments
tune!(rsegments1,0.5)
map((s)->tune!(s,0.5),lsegments1)
@test_approx_eq cov(lsegments1[1].samplerstate.density.distribution) cov(fetch(rsegments1.remote[1]).samplerstate.density.distribution)

println()
println("====================")
println("Test show() function")
show(rsegments1)
println("End  show() function")
println("====================")
println()

nothing


text:
module SnakeLadder

include("dices.jl")
include("traps.jl")
include("markovDecision.jl")
include("simulate.jl")
include("perf_plots.jl")
include("experiments.jl")

end # module

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
# Turn a collection of markdown files into a browsable multi-page manual.

# Files to be weaved when generating a packages documentation.
const ext_doc_pat = r"\.(md|txt|rst)$"i

# Pattern for matching github package urls
const pkgurl_pat = r"github.com/(.*)\.git$"

# Generate documentation from the given package.
function collate(package::String; template::String="default")
    declarations = harvest(package)

    pkgver = ""
    try
        pkgver = Pkg.Dir.cd(() -> Pkg.Read.installed()[package][1])
    catch
    end
    pkgver_out = IOBuffer()
    print(pkgver_out, pkgver)
    pkgver = takebuf_string(pkgver_out)

    pkgurl = "#"
    try
        pkgurl_mat = match(pkgurl_pat, Pkg.Dir.cd(() -> Pkg.Read.url(package)))
        if pkgurl_mat != nothing
            pkgurl = string("http://github.com/", pkgurl_mat.captures[1])
        end
    catch
    end

    filenames = {}
    for filename in walkdir(joinpath(Pkg.dir(package), "doc"))
        if match(ext_doc_pat, filename) != nothing
            push!(filenames, filename)
        end
    end

    outdir = joinpath(Pkg.dir(package), "doc", "html")
    if !isdir(outdir)
        mkdir(outdir)
    end

    collate(filenames, template=template, outdir=outdir, pkgname=package,
            pkgver=pkgver, pkgurl=pkgurl, declarations=declarations)
end

# Generate documentation from multiple files.
function collate(filenames::Vector;
                 template::String="default",
                 outdir::String=".",
                 declarations::Dict=Dict(),
                 pkgname=nothing,
                 pkgver=nothing,
                 pkgurl=nothing,
                 pandocargs=nothing)
    toc = Dict()
    titles = Dict()

    declaration_markdown = generate_declaration_markdown(declarations)

    if !isdir(template)
        template = joinpath(Pkg.dir("Judo"), "templates", template)
        if !isdir(template)
            error("Can't find template $(template)")
        end
    end

    # make any expansions necessary in the original document
    docs = Dict()
    for filename in filenames
        name = choose_document_name(filename)
        docs[name] = expand_declaration_docs(readall(filename),
                                             declaration_markdown)
    end

    # dry-run to collect the section names in each document
    for (name, doc) in docs
        metadata, sections = weave(IOBuffer(doc), IOBuffer(), dryrun=true)
        title = get(metadata, "title", name)
        titles[name] = title
        part = get(metadata, "part", nothing)
        if !haskey(toc, part)
            toc[part] = {}
        end

        push!(toc[part], (parseint(get(metadata, "order", "0")), name, title, sections))
    end
    for part_content in values(toc)
        sort!(part_content)
    end

    pandoc_template = joinpath(template, "template.html")

    keyvals = Dict()

    if pkgname != nothing
        keyvals["pkgname"] = pkgname
    end

    if pkgver != nothing
        keyvals["pkgver"] = pkgver
    end

    if pkgurl != nothing
        keyvals["pkgurl"] = pkgurl
    end

    for (name, doc) in docs
        println(STDERR, "weaving ", name)
        fmt = :markdown
        title = titles[name]
        outfilename = joinpath(outdir, string(name, ".html"))
        outfile = open(outfilename, "w")
        keyvals["table-of-contents"] = table_of_contents(toc, title)
        metadata, sections = weave(IOBuffer(doc), outfile,
                                   name=name, template=pandoc_template,
                                   toc=true, outdir=outdir,
                                   keyvals=keyvals,
                                   pandocargs=pandocargs)
        close(outfile)
    end

    # copy template files
    run(`cp -r $(joinpath(template, "js")) $(outdir)`)
    run(`cp -r $(joinpath(template, "css")) $(outdir)`)
end

# pattern for choosing document names
const fileext_pat = r"^(.+)\.([^\.]+)$"

# Choose a documents name from its file name.
function choose_document_name(filename::String)
    filename = basename(filename)
    mat = match(fileext_pat, filename)
    name = filename
    if !is(mat, nothing)
        name = mat.captures[1]
        if mat.captures[2] == "md"
            fmt = :markdown
        elseif mat.captures[2] == "rst"
            fmt = :rst
        elseif mat.captures[2] == "tex"
            fmt = :latex
        elseif mat.captures[2] == "html" || mat.captures[2] == "htm"
            fmt = :html
        else
            name = filename
        end
    end

    name
end

# Generate a table of contents for the given document.
function table_of_contents(toc, selected_title::String)
    parts = collect(keys(toc))
    part_order = [minimum([order for (order, name, title, sections) in toc[part]])
                  for part in parts]
    out = IOBuffer()
    write(out, "
")
    for part in parts[sortperm(part_order)]
        if part != nothing
            write(out,
                """

$(part)

                """)
        end

        for (order, name, title, sections) in toc[part]
            iscurrent = title == selected_title

            classes = iscurrent ?
                "toc-item toc-current-doc" : "toc-item"

            write(out,
                """

                    $(title)

                """)

                write(out, table_of_contents_sections(name, sections, iscurrent,
                    maxlevel=iscurrent ? 2 : 0))
        end
    end
    write(out, "
")
    takebuf_string(out)
end

function table_of_contents_sections(parent, sections, iscurrent; maxlevel=2)
    if isempty(sections)
        return ""
    end

    out = IOBuffer()
    current_level = 0
    for (level, section) in sections
        if level > maxlevel
            continue
        end

        while level > current_level
            current_level += 1
        end

        while level < current_level
            current_level -= 1
        end

        href = iscurrent ?
            "#$(section_id(section))" :
            "$(parent).html#$(section_id(section))"

        write(out,
            """

                $(section)

            """)
    end
    return takebuf_string(out)
end

# Turn a section name into an html id.
function section_id(section::String)
    # Keep only unicode letters, _ and -
    cleaned = replace(section, r"[^\p{L}_\-\s]", "")
    lowercase(replace(cleaned, r"\s+", "-"))
end


text:
################################################################################
# Visuals
################################################################################
function visualize(env::Environment, traj::Vector{Vector{T}}; build::Bool=true) where T
	@assert size(traj[1]) == size(env.state)
    storage = generate_storage(env.mechanism, [env.representation == :minimal ? minimal_to_maximal(env.mechanism, x) : x for x in traj])
    Dojo.visualize(env.mechanism, storage, vis=env.vis, build=build)
end

################################################################################
# Miscellaneous
################################################################################
type2symbol(H) = Symbol(lowercase(String(H.name.name)))

function get_control_mask(n_inputs, indices)
    n_outputs = length(indices)
    m = zeros(n_outputs, n_inputs)
    for (i,ind) in enumerate(indices)
        m[i,ind] = 1.0
    end
    return m
end

function set_springs!(joints, spring::Real)
    if spring > 0
        for joint in joints[2:end]
            joint.spring = true
            joint.translational.spring=spring
            joint.rotational.spring=spring
        end
    end
end

function set_springs!(joints, springs::AbstractArray)
    for (i,spring) in enumerate(springs)
        if spring > 0
            joints[i+1].spring = true
            joints[i+1].translational.spring=spring
            joints[i+1].rotational.spring=spring
        end
    end
end

function set_dampers!(joints, damper::Real)
    if damper > 0
        for joint in joints[2:end]
            joint.damper = true
            joint.translational.damper=damper
            joint.rotational.damper=damper
        end
    end
end

function set_dampers!(joints, dampers::AbstractArray)
    for (i,damper) in enumerate(dampers)
        if damper > 0
            joints[i+1].damper = true
            joints[i+1].translational.damper=damper
            joints[i+1].rotational.damper=damper
        end
    end
end


text:
"""
    DPI

Pi
"""
const DPI = 3.141592653589793

"""
    D2PI

2Pi
"""
const D2PI = 6.283185307179586

"""
    DR2D

Radians to degrees
"""
const DR2D = 57.29577951308232

"""
    DD2R

Degrees to radians
"""
const DD2R = 0.017453292519943295

"""
    DR2AS

Radians to arcseconds
"""
const DR2AS = 206264.80624709636

"""
    DAS2R

Arcseconds to radians
"""
const DAS2R = 4.84813681109536e-6

"""
    DS2R

Seconds of time to radians
"""
const DS2R = 7.27220521664304e-5

"""
    TURNAS

Arcseconds in a full circle
"""
const TURNAS = 1.296e6

"""
    DMAS2R

Milliarcseconds to radians
"""
const DMAS2R = DAS2R / 1000.0

"""
    DTY

Length of tropical year B1900 (days)
"""
const DTY = 365.242198781

"""
    DAYSEC

Seconds per day
"""
const DAYSEC = 86400.0

"""
    DJY

Days per Julian year
"""
const DJY = 365.25

"""
    DJC

Days per Julian century
"""
const DJC = 36525.0

"""
    DJM

Days per Julian millennium
"""
const DJM = 365250.0

"""
    DJ00

Reference epoch (J2000.0), Julian Date
"""
const DJ00 = 2.451545e6

"""
    DJM0

Julian Date of Modified Julian Date zero
"""
const DJM0 = 2.4000005e6

"""
    DJM00

Reference epoch (J2000.0), Modified Julian Date
"""
const DJM00 = 51544.5

"""
    DJM77

1977 Jan 1.0 as MJD
"""
const DJM77 = 43144.0

"""
    TTMTAI

TT minus TAI (s)
"""
const TTMTAI = 32.184

"""
    DAU

Astronomical unit (m, IAU 2012)
"""
const DAU = 1.4959787e11

"""
    CMPS

Speed of light (m/s)
"""
const CMPS = 2.99792458e8

"""
    AULT

Light time for 1 au (s)
"""
const AULT = 499.004782

"""
    DC

Speed of light (au per day)
"""
const DC = DAYSEC / AULT

"""
    ELG

L_G = 1 - d(TT)/d(TCG)
"""
const ELG = 6.969290134e-10

"""
    L_B

L_B = 1 - d(TDB)/d(TCB)
"""
const ELB = 1.550519768e-8

"""
    TDB0

TDB (s) at TAI 1977/1/1.0
"""
const TDB0 = -6.55e-5

"""
    SRS

Schwarzschild radius of the Sun (au) =
2 * 1.32712440041e20 / (2.99792458e8)^2 / 1.49597870700e11
"""
const SRS = 1.97412574336e-8

"""
    Ellipsoid

Reference ellipsoids

- [`WGS84`](@ref)
- [`GRS80`](@ref)
- [`WGS72`](@ref)
"""
@enum Ellipsoid WGS84=1 GRS80=2 WGS72=3

"""
    WGS84

Enum constant for the WGS84 reference ellipsoid
"""
WGS84

"""
    GRS80

Enum constant for the GRS80 reference ellipsoid
"""
GRS80

"""
    WGS72

Enum constant for the WGS72 reference ellipsoid
"""
WGS72

mutable struct ASTROM
    pmt::Cdouble
    eb::NTuple{3,Cdouble}
    eh::NTuple{3,Cdouble}
    em::Cdouble
    v::NTuple{3,Cdouble}
    bm1::Cdouble
    bpn::NTuple{9,Cdouble}
    along::Cdouble
    phi::Cdouble
    xpl::Cdouble
    ypl::Cdouble
    sphi::Cdouble
    cphi::Cdouble
    diurab::Cdouble
    eral::Cdouble
    refa::Cdouble
    refb::Cdouble
end

function ASTROM()
    ASTROM(
        0.0,
        zeros(Cdouble, 3),
        zeros(Cdouble, 3),
        0.0,
        zeros(Cdouble, 3),
        0.0,
        zeros(Cdouble, 3, 3),
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
    )
end

function Base.getproperty(a::ASTROM, field::Symbol)
    val = getfield(a, field)
    if field in (:eb, :eh, :v)
        return collect(val)
    elseif field == :bpn
        return permutedims(reshape(collect(val), 3, 3))
    end
    return val
end

struct LDBODY
    bm::Cdouble
    dl::Cdouble
    pv::NTuple{6,Cdouble}
end

struct ERFAException <: Exception
    msg::String
end

Base.showerror(io::IO, ex::ERFAException) = print(io, ex.msg)

macro checkdims(m::Int, n::Int, arr::Symbol...)
    ex = :()
    for a in arr
        name = string(a)
        expr = quote
            m1, n1 = size($(esc(a)))
            if (m1, n1) != ($m, $n)
                msg = string("`$($name)` must be a $($m)x$($n) matrix but is ", m1, "x", n1, ".")
                throw(ArgumentError(msg))
            end
        end
        push!(ex.args, expr)
    end
    :($ex; nothing)
end

macro checkdims(len::Int, arr::Symbol...)
    ex = :()
    for a in arr
        name = string(a)
        expr = quote
            n = length($(esc(a)))
            if n != $(esc(len))
                throw(ArgumentError("`$($name)` must have $($(esc(len))) elements but has $n."))
            end
        end
        push!(ex.args, expr)
    end
    :($ex; nothing)
end

function array_to_cmatrix(array; n=0)
    any(length.(array) .!= n) && throw(ArgumentError("Expected each vector to have $n elements."))
    return hcat(array...)
end

function cmatrix_to_array(matrix)
    return [matrix[:, i] for i in 1:size(matrix, 2)]
end

function ASTROM(pmt, eb::AbstractArray, eh::AbstractArray, em, v::AbstractArray, bm1, bpn::AbstractArray, along, phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb)
    ASTROM(pmt,
              (eb[1], eb[2], eb[3]),
              (eh[1], eh[2], eh[3]),
              em,
              (v[1], v[2], v[3]),
              bm1,
              (bpn[1], bpn[2], bpn[3], bpn[4], bpn[5], bpn[6], bpn[7], bpn[8], bpn[9]),
              along,
              phi,
              xpl,
              ypl,
              sphi,
              cphi,
              diurab,
              eral,
              refa,
              refb)
end

function LDBODY(bm, dl, pv::AbstractArray)
    LDBODY(bm, dl, (pv...,))
end


text:
using Test #, Random

@testset "Testing solution to Exercise 2" begin

@testset "Running ex2.jl" begin
   include("../ex2.jl")
end;

@testset "Testing that variables exist" begin
   @test @isdefined(response_2a)
   @test @isdefined(response_2b)
   @test @isdefined(response_2c)
   @test @isdefined(response_2d)
   #@test @isdefined(response_2e)
   @test @isdefined(response_2f)
   @test @isdefined(response_2g)
   @test @isdefined(response_2h)
end;

@testset "Testing that variables are not missing" begin
   @test check_type_isa(:response_2a,response_2a,Markdown.MD)
   @test check_type_isa(:response_2b,response_2b,Markdown.MD)
   @test check_type_isa(:response_2c,response_2c,Markdown.MD)
   @test check_type_isa(:response_2d,response_2d,Markdown.MD)
   #@test check_type_isa(:response_2e,response_2e,Markdown.MD)
   @test check_type_isa(:response_2f,response_2f,Markdown.MD)
   @test check_type_isa(:response_2g,response_2g,Markdown.MD)
   @test check_type_isa(:response_2h,response_2h,Markdown.MD)
end;

end; # Exercise 2


text:
@testset "cl5" begin
    data = open(readdlm, joinpath(@__DIR__, "data", "Cl5.txt"))

    for r in 1:size(data, 1)
        row      = data[r, :]
        x        = row[1]
        expected = row[2]

        @test ClausenFunctions.cl5(x) ≈ expected atol=1e-14
        @test ClausenFunctions.cl5(-x) ≈ expected atol=1e-14
        @test ClausenFunctions.cl5(x - 2.0*pi) ≈ expected atol=1e-13
        @test ClausenFunctions.cl5(x + 2.0*pi) ≈ expected atol=1e-13
    end
end