text:
using PyPlot, JLD
savefigs = true
savepath = "scripts/saved_run_0x3f77869c073f2346/"

#figure 1
#' Plot
fig, ax = plt[:subplots](1,1, figsize=(6,3.5))
fontsize = 20
x_grid = [x for x in range(-10,21) for y in range(0,21)]
y_grid = [y for x in range(-10,21) for y in range(0,21)]
ax[:scatter](x_grid, y_grid, s=1.3, color="k")
Θ = linspace(0., π, 100)
r0= 4.3
r1= sqrt(2)*r0
r2= sqrt(3)*r0
r3= sqrt(4)*r0
ax[:plot](r0*cos(Θ), r0*sin(Θ))
ax[:plot](r1*cos(Θ), r1*sin(Θ))
ax[:plot](r2*cos(Θ), r2*sin(Θ))
ax[:plot](r3*cos(Θ), r3*sin(Θ))
ax[:text](0, -1, L"$\ell_1$",  horizontalalignment="center", fontsize=fontsize)
ax[:text](-12, 5, L"$\ell_2$",  horizontalalignment="center", rotation=90, fontsize=fontsize)
ax[:text](0, 1, L"$^0B^{\rm del}_{\bf l}$",  horizontalalignment="center", fontsize=fontsize)
ax[:text](0, 4.5, L"$^1B^{\rm del}_{\bf l}$",  horizontalalignment="center", fontsize=fontsize-1)
ax[:text](0, 6.2, L"$^2B^{\rm del}_{\bf l}$",  horizontalalignment="center", fontsize=fontsize-2)
ax[:text](0, 7.6, L"$^3B^{\rm del}_{\bf l}$",  horizontalalignment="center", fontsize=fontsize-3)
ax[:set_xlim]([-11, 11])
ax[:set_ylim]([0, 11])
ax[:set_aspect]("equal")
ax[:axis]("off")

fig[:tight_layout]()

#' Now save
savename  = "f1.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

include("src/PrimordialB.jl")

kbins, N0ℓ_a, cϕϕ = let
    nside     = load(savepath*"out_figure2a.jld", "nside")
    period    = load(savepath*"out_figure2a.jld", "period")
    μKarcminT = load(savepath*"out_figure2a.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2a.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2a.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end
kbins, N0ℓ_b, cϕϕ = let
    nside     = load(savepath*"out_figure2b.jld", "nside")
    period    = load(savepath*"out_figure2b.jld", "period")
    μKarcminT = load(savepath*"out_figure2b.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2b.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2b.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end
kbins, N0ℓ_c, cϕϕ = let
    nside     = load(savepath*"out_figure2c.jld", "nside")
    period    = load(savepath*"out_figure2c.jld", "period")
    μKarcminT = load(savepath*"out_figure2c.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2c.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2c.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end

#figure 2

fig, ax = plt[:subplots](1,1, figsize=(4.5,3.5))
fontsize = 12
ax[:plot](kbins, cϕϕ, lw=2, label = L"$\ell^4 C_\ell^{\phi\phi}$")
ax[:plot](kbins, N0ℓ_a, lw=2, label = L"$\Delta_{\rm T} = 10 \ \mu{\rm K-arcmin}$")
ax[:plot](kbins, N0ℓ_b, lw=2, label = L"$\Delta_{\rm T} = 1 \ \mu{\rm K-arcmin}$")
ax[:plot](kbins, N0ℓ_c, lw=2, label = L"$\Delta_{\rm T} = 0.5 \ \mu{\rm K-arcmin}$")
ax[:set_xscale]("log")
ax[:set_yscale]("log")
ax[:set_xlim]([10.,1.2e3])
ax[:set_ylim]([1.e-8,4e-6])
ax[:set_title](L"$\phi$ uncertainty", fontsize=fontsize)
ax[:set_xlabel](L"multiple $\ell$", fontsize=fontsize)
ax[:legend](loc="best", fontsize=fontsize-2)
#ax[:set_ylim]([1e-8, 1.5e-6])
ax[:grid](true)

#' Now save
savename  = "f2.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#figure 3

rng_true = load(savepath*"out_figure2a.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure2a.jld", "rest_B")
rest_1   = load(savepath*"out_figure2a.jld", "rest_modes")
rest_B_2 = load(savepath*"out_figure2b.jld", "rest_B")
rest_2   = load(savepath*"out_figure2b.jld", "rest_modes")
rest_B_3 = load(savepath*"out_figure2c.jld", "rest_B")
rest_3   = load(savepath*"out_figure2c.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_B_2  = mean((rest_B_2 .- rng_true),2)
bias_B_3  = mean((rest_B_3 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)
bias_L_2  = mean((rest_2 .- rng_true),2)
bias_L_3  = mean((rest_3 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_B_2 = sqrt.(mean((rest_B_2 .- rng_true).^2, 2))
standard_error_B_3 = sqrt.(mean((rest_B_3 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))
standard_error_L_2 = sqrt.(mean((rest_2 .- rng_true).^2, 2))
standard_error_L_3 = sqrt.(mean((rest_3 .- rng_true).^2, 2))

maxylab = maximum(vcat(
    rng_true./standard_error_B_1,
    rng_true./standard_error_B_2,
    rng_true./standard_error_B_3,
    rng_true./standard_error_L_1,
    rng_true./standard_error_L_2,
    rng_true./standard_error_L_3
    ))

minbias, maxbias = vcat(
         bias_B_1./standard_error_B_1,
         bias_B_2./standard_error_B_2,
         bias_B_3./standard_error_B_3,
         bias_L_1./standard_error_L_1,
         bias_L_2./standard_error_L_2,
         bias_L_3./standard_error_L_3,
         ) |> x-> (minimum(x), maximum(x))

#' f3.pdf : plot r / σ(r) and Bias(r)/r
#' ================================================

fig, ax = plt[:subplots](2,3, figsize=(10,6), sharex=true)
fontsize = 13

#  CMB-S3 for ϕ est
ax[1,1][:set_title](L"$\Delta_{\rm T} = 10 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[1,1][:plot](rng_true, rng_true./standard_error_L_1, "b", lw = 2,label = L"${\rm likelihood\ approx}$")
ax[1,1][:plot](rng_true, rng_true./standard_error_B_1, "r.", lw = 2,label = L"${\rm quadratic\ delenser}$")
ax[1,1][:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4, rotation="horizontal")
#ax[1,1][:set_xlabel](L"$r$", fontsize=fontsize)
ax[1,1][:set_ylim]([0,17])
ax[1,1][:set_xticks]([.01, .03, .05, .07, .09])
ax[1,1][:tick_params](axis="y", labelsize = fontsize-2)
ax[1,1][:tick_params](axis="x", labelsize = fontsize-2)
ax[1,1][:legend](loc="best", fontsize=fontsize)
ax[1,1][:grid](true)

#  CMB-S4 for ϕ est
ax[1,2][:set_title](L"$\Delta_{\rm T} = 1 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[1,2][:plot](rng_true, rng_true./standard_error_B_2, "r.",lw = 2,label = L"${\rm quadratic\ delenser}$")
ax[1,2][:plot](rng_true, rng_true./standard_error_L_2, "b-",lw = 2,label = "likelihood approx")
#ax[1,2][:set_xlabel](L"$r$", fontsize=fontsize)
ax[1,2][:set_ylim]([0,17])
ax[1,2][:tick_params](axis="y", left="off", right="off")#labelsize = fontsize-2, pad = 5, labelright="on")
ax[1,2][:tick_params](axis="x", labelsize = fontsize-2)
ax[1,2][:grid](true)

#  CMB-S4 for ϕ est
ax[1,3][:set_title](L"$\Delta_{\rm T} = 0.5 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[1,3][:plot](rng_true, rng_true./standard_error_B_3, "r.",lw = 2,label = "quadratic delenser")
ax[1,3][:plot](rng_true, rng_true./standard_error_L_3, "b-",lw = 2,label = "likelihood approx")
#ax[1,3][:set_xlabel](L"$r$", fontsize=fontsize)
ax[1,3][:set_ylim]([0,17])
ax[1,3][:tick_params](axis="y", labelsize = fontsize-2, labelright="on")
ax[1,3][:tick_params](axis="x", labelsize = fontsize-2, pad = 5)
ax[1,3][:axis]("tight")
ax[1,3][:grid](true)

#  CMB-S3 for ϕ est
#ax[2,1][:set_title](L"$\Delta_{\rm T} = 10 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[2,1][:plot](rng_true, bias_L_1./rng_true, "b-", lw =2, label = L"${\rm likelihood\ approx}$")
ax[2,1][:plot](rng_true, bias_B_1./rng_true, "r.", lw = 2, label = L"${\rm quadratic\ delenser}$")
ax[2,1][:set_ylabel](L"$\frac{{\rm Bias}(r)}{r}$        ", fontsize=fontsize+4, rotation="horizontal")
ax[2,1][:set_xlabel](L"$r$", fontsize=fontsize)
ax[2,1][:set_ylim]([-0.02,0.2])
ax[2,1][:set_xticks]([.01, .03, .05, .07, .09])
ax[2,1][:tick_params](axis="y", labelsize = fontsize-2)
ax[2,1][:tick_params](axis="x", labelsize = fontsize-2)
ax[2,1][:legend](loc="best", fontsize=fontsize)
ax[2,1][:grid](true)
ax[2,1][:axhline](y=.0, "k.")

#  CMB-S4 for ϕ est
#ax[2,2][:set_title](L"$\Delta_{\rm T} = 1 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[2,2][:plot](rng_true, bias_B_2./rng_true, "r.", lw = 2,label = "quadratic delenser")
ax[2,2][:plot](rng_true, bias_L_2./rng_true, "b-", lw = 2,label = "likelihood approx")
ax[2,2][:set_xlabel](L"$r$", fontsize=fontsize)
ax[2,2][:set_ylim]([-0.02,0.2])
ax[2,2][:tick_params](axis="y", left="off", right="off")
ax[2,2][:tick_params](axis="x", labelsize = fontsize-2)
ax[2,2][:grid](true)

#  CMB-S4 for ϕ est
#ax[2,3][:set_title](L"$\Delta_{\rm T} = 0.5 \ \mu{\rm K-arcmin}$", fontsize=fontsize)
ax[2,3][:plot](rng_true, bias_B_3./rng_true, "r.", lw = 2,label = "quadratic delenser")
ax[2,3][:plot](rng_true, bias_L_3./rng_true, "b-", lw = 2,label = "likelihood approx")
ax[2,3][:set_xlabel](L"$r$", fontsize=fontsize)
ax[2,3][:set_ylim]([-0.02,0.2])
ax[2,3][:tick_params](axis="y", labelsize = fontsize-2, labelright="on")
ax[2,3][:tick_params](axis="x", labelsize = fontsize-2)
#ax[2,3][:axis]("tight")
ax[2,3][:grid](true)

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f3.pdf"
if savefigs
      fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#for figure 4 use
dtec_level_full = rng_true./standard_error_B_2
bias_level_full = bias_B_2./rng_true

# figure 4: Explore higher order bias impact on the quadraitc delenser

rng_true = load(savepath*"out_figure3a.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure3a.jld", "rest_B")
rest_1   = load(savepath*"out_figure3a.jld", "rest_modes")
rest_B_2 = load(savepath*"out_figure3b.jld", "rest_B")
rest_2   = load(savepath*"out_figure3b.jld", "rest_modes")

rest_B_3 = load(savepath*"out_figure3c.jld", "rest_B")
rest_3   = load(savepath*"out_figure3c.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_B_2  = mean((rest_B_2 .- rng_true),2)
bias_B_3  = mean((rest_B_3 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)
bias_L_2  = mean((rest_2 .- rng_true),2)
bias_L_3  = mean((rest_3 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_B_2 = sqrt.(mean((rest_B_2 .- rng_true).^2, 2))
standard_error_B_3 = sqrt.(mean((rest_B_3 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))
standard_error_L_2 = sqrt.(mean((rest_2 .- rng_true).^2, 2))
standard_error_L_3 = sqrt.(mean((rest_3 .- rng_true).^2, 2))

maxylab = maximum(vcat(
        rng_true./standard_error_B_1,
        rng_true./standard_error_B_2,
        rng_true./standard_error_B_3,
        rng_true./standard_error_L_1,
        rng_true./standard_error_L_2,
        rng_true./standard_error_L_3
         ))

minbias, maxbias = vcat(
        bias_B_1./standard_error_B_1,
        bias_B_2./standard_error_B_2,
        bias_B_3./standard_error_B_3,
        bias_L_1./standard_error_L_1,
        bias_L_2./standard_error_L_2,
        bias_L_3./standard_error_L_3,
         ) |> x-> (minimum(x), maximum(x))

#' f4a.pdf : plot r / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 13

#  CMB-S3 for ϕ est
ax[:set_title](L"$\Delta_{\rm T} = 1 \ \mu{\rm K-arcmin}$", fontsize=fontsize-1)
ax[:plot](rng_true, rng_true./standard_error_B_1, "k-", lw=2, label = L"${\rm Linear\ order\ (null\ test)}$")
ax[:plot](rng_true, rng_true./standard_error_B_2, "g:", lw=6, label = L"${\rm All\ order\ E\ from\ E}$")
ax[:plot](rng_true, rng_true./standard_error_B_3, "b--", lw=2, label = L"${\rm All\ order\ B\ from\ E}$")
ax[:plot](rng_true, dtec_level_full, "r.", lw=2, label = L"${\rm All\ order\ E/B\ from\ E}$")
ax[:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
#ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:set_xlim]([0.,0.1])

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f4a.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

fig, ax = plt[:subplots](1,1, figsize=(4.5,3.7), sharex=true,  sharey=true)
fontsize = 13

#  CMB-S3 for ϕ est
ax[:set_title](L"$\Delta_{\rm T} = 1 \ \mu{\rm K-arcmin}$", fontsize=fontsize-1)
ax[:plot](rng_true, bias_B_1./rng_true, "k-", lw=2, label = L"${\rm Linear\ order\ (null\ test)}$")
ax[:plot](rng_true, bias_B_2./rng_true, "g:", lw=6, label = L"${\rm All\ order\ E\ from\ E}$")
ax[:plot](rng_true, bias_B_3./rng_true, "b--", lw=2, label = L"${\rm All\ order\ B\ from\ E}$")
ax[:plot](rng_true, bias_level_full, "r.", lw=2, label = L"${\rm All\ order\ E/B\ from\ E}$")
ax[:set_ylabel](L"$\frac{{\rm Bias}(r)}{r}$        ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize-2)
ax[:grid](true)
ax[:set_xlim]([0.,0.1])
#ax[:axhline](y=.0, color = "k", linestyle="dashed")

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f4b.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#figure 5

rng_true = load(savepath*"out_figure4.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure4.jld", "rest_B")
rest_1   = load(savepath*"out_figure4.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))

maxylab = maximum(vcat(
    rng_true./standard_error_B_1,
    rng_true./standard_error_L_1,
    ))

minbias, maxbias = vcat(
         bias_B_1./standard_error_B_1,
         bias_L_1./standard_error_L_1,
         ) |> x-> (minimum(x), maximum(x))

#' f5a.pdf :
#' ================================================

fig = figure(figsize=(4.5,3.7))
fontsize = 13
σmod       = load(savepath*"out_figure4.jld", "σmod")
μKarcminT  = load(savepath*"out_figure4.jld", "μKarcminT")
#imshow(fftshift(σmod .* √2 .* μKarcminT))
imshow(fftshift(σmod))
colorbar()
#title("Non-stationary noise modulation", fontsize=fontsize)
text(0., 220., "Non-stationary noise modulation", fontsize=fontsize)
axis("tight")
axis("off")
#fig[:suptitle]("    Bias", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f5a.pdf"

if savefigs
  fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f5b.pdf : plot r / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 13

#  CMB-S3 for ϕ est
#ax[:set_title](L"Nonstationary noise (baseline $\mu K arcmin = 1$)", fontsize=fontsize-1)
ax[:plot](rng_true, rng_true./standard_error_L_1, "b-", lw=2, label = L"${\rm likelihood\ approx}$")
ax[:plot](rng_true, rng_true./standard_error_B_1, "r.", lw=2, label = L"${\rm quadratic\ delenser}$")
ax[:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:set_xlim]([0.,0.1])

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f5b.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f5c.pdf : plot  Bias(r) / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.5,3.7), sharex=true,  sharey=true)
fontsize = 13

#  CMB-S3 for ϕ est
#ax[:set_title](L"Nonstationary noise (baseline $\mu K arcmin = 1$)", fontsize=fontsize-1)
ax[:plot](rng_true, bias_L_1./rng_true, "b-", lw=2, label = L"${\rm likelihood\ approx}$")
ax[:plot](rng_true, bias_B_1./rng_true, "r.", lw=2, label =L"${\rm quadratic\ delenser}$")
ax[:set_ylabel](L"$\frac{{\rm Bias}(r)}{r}$        ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:set_xlim]([0.,0.1])
#ax[:axhline](y=.0, color = "k", linestyle="dashed")

#fig[:suptitle]("    Bias", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f5c.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end


text:
module GetData

import Base.close
# package code goes here
export gd_open,
    gd_nfields,
    nvectors,
    field_list,
    vector_list,
    getdata

type dirfile
end

function gd_open(dirfilename::AbstractString)
    flags = 0
    f = ccall( (:gd_open, "libgetdata"), Ptr{dirfile},
        (Cstring, Culong), dirfilename, flags)
end

function close(d::Ptr{dirfile})
    ccall( (:gd_close, "libgetdata"), Cint,
        (Ptr{dirfile}, ), d)
end

function gd_nfields(d::Ptr{dirfile})
    ccall( (:gd_nfields, "libgetdata"), Cuint,
        (Ptr{dirfile}, ), d)
end

function nvectors(d::Ptr{dirfile})
    ccall( (:gd_nvectors, "libgetdata"), Cuint,
        (Ptr{dirfile}, ), d)
end

function decode_stringarray(fptr::Ptr{Cstring}, nstr::Integer)
    const ptrsize = sizeof(fptr)
    [unsafe_string(unsafe_load(fptr+i*ptrsize)) for i=0:nstr-1]
end

function field_list(d::Ptr{dirfile})
    fptr = ccall( (:gd_field_list, "libgetdata"), Ptr{Cstring},
        (Ptr{dirfile}, ), d)
    decode_stringarray(fptr, gd_nfields(d))
end

function vector_list(d::Ptr{dirfile})
    fptr = ccall( (:gd_vector_list, "libgetdata"), Ptr{Cstring},
        (Ptr{dirfile}, ), d)
    decode_stringarray(fptr, gd_nfields(d))
end

#Data may be fetched from a vector field in the dirfile (including metafields) with

# size_t gd_getdata(DIRFILE *dirfile, const char *field_code, off_t first_frame, off_t first_sample, size_t num_frames, size_t num_samples, gd_type_t return_type, void *data_out);

#Here return_type identifies the desired type of the returned data, which need not be the same as the type of the data in the dirfile. For SINDIR fields, return_type must be GD_STRING or GD_NULL. For numeric vector fields, any other value may be used for return_type. Type conversion will be performed as necessary. This function returns the number of samples (not frames or bytes) successfully read, or zero on error. If return_type is GD_NULL, no data is returned, and data_out is ignored; however, the number of samples read is still returned.

function getdata(d::Ptr{dirfile}, field_code::AbstractString, first_frame::Csize_t,
    first_sample::Csize_t, num_frames::Csize_t, num_samples::Csize_t)
    return_type::Cint = 0x84 # float 32
    data = Array{Cfloat}(num_frames*num_samples)
    num_fetched = ccall( (:gd_getdata, "libgetdata"), Csize_t,
        (Ptr{dirfile}, Cstring, Csize_t, Csize_t, Csize_t, Csize_t, Cint, Ref{Cfloat}),
        d, field_code, first_frame, first_sample, num_frames, num_samples, return_type, data)
    (num_fetched, data)
end

end # module

ext: jl


text:
using Plots

# Model
include_model("planar_push_block")

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[model.idx_u[1:8]] .= 5.0
_uu[model.idx_λ] .= 1.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[model.idx_u[1:8]] .= -5.0
_ul[model.idx_λ] .= 1.0

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 0.0]
x1 = [q1; q1]
qT = [1.0; 0.0; 0.0]
xT = [qT; qT]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
include_objective("velocity")
obj_velocity = velocity_objective(
    [t > T / 2 ? Diagonal(10.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3]))

obj_tracking = quadratic_tracking_objective(
    [Diagonal(1.0 * ones(model.n)) for t = 1:T],
    # [Diagonal(0.1 * ones(model.m)) for t = 1:T-1],
	[Diagonal([0.1 * ones(model.nu);
		zeros(model.nc);
		ones(model.nb);
		zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_penalty = PenaltyObjective(1.0e4, model.m)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity])

# Constraints
include_constraints(["contact", "stage"])
t_idx = vcat([t for t = 1:T-1])

function control_comp_con!(c, x, u, t)
	q = x[nq .+ (1:nq)]
	s = u[model.idx_s]
    k = control_kinematics_func(model, q)
	k_input = u[model.idx_u][9:10]

	e1 = k[1:2] - k_input
	e2 = k[3:4] - k_input
	e3 = k[5:6] - k_input
	e4 = k[7:8] - k_input
	# e5 = k[9:10] - k_input
	# e6 = k[11:12] - k_input
	# e7 = k[13:14] - k_input
	# e8 = k[15:16] - k_input

	d1 = e1' * e1
	d2 = e2' * e2
	d3 = e3' * e3
	d4 = e4' * e4
	# d5 = e5' * e5
	# d6 = e6' * e6
	# d7 = e7' * e7
	# d8 = e8' * e8

	u_ctrl = u[model.idx_u]

	c[1] = s[1] - u_ctrl[1] * d1
	c[2] = s[1] + u_ctrl[1] * d1

	c[3] = s[1] - u_ctrl[2] * d1
	c[4] = s[1] + u_ctrl[2] * d1

	c[5] = s[1] - u_ctrl[3] * d2
	c[6] = s[1] + u_ctrl[3] * d2

	c[7] = s[1] - u_ctrl[4] * d2
	c[8] = s[1] + u_ctrl[4] * d2

	c[9] = s[1] - u_ctrl[5] * d3
	c[10] = s[1] + u_ctrl[5] * d3

	c[11] = s[1] - u_ctrl[6] * d3
	c[12] = s[1] + u_ctrl[6] * d3

	c[13] = s[1] - u_ctrl[7] * d4
	c[14] = s[1] + u_ctrl[7] * d4

	c[15] = s[1] - u_ctrl[8] * d4
	c[16] = s[1] + u_ctrl[8] * d4

	# c[9] = s[1] - u[9] * d5
	# c[10] = s[1] - u[10] * d5
	#
	# c[11] = s[1] - u[11] * d6
	# c[12] = s[1] - u[12] * d6
	#
	# c[13] = s[1] - u[13] * d7
	# c[14] = s[1] - u[14] * d7
	#
	# c[15] = s[1] - u[15] * d8
	# c[16] = s[1] - u[16] * d8

    nothing
end

n_ctrl_comp = 16
con_ctrl_comp = stage_constraints(control_comp_con!, n_ctrl_comp, (1:16), t_idx)

function control_limits_con!(c, x, u, t)
	q = x[nq .+ (1:nq)]
	θ = q[3]
	R = rotation_matrix(θ)
	u_ctrl = u[model.idx_u]

	c[1] = -1.0 * (R' * u_ctrl[1:2])[1]
	c[2] = (model.μ[end] * u_ctrl[1])^2.0 - u_ctrl[2]^2.0

	c[3] = -1.0 * (R' * u_ctrl[3:4])[2]
	c[4] = (model.μ[end] * u_ctrl[4])^2.0 - u_ctrl[3]^2.0

	c[5] = (R' * u_ctrl[5:6])[1]
	c[6] = (model.μ[end] * u_ctrl[5])^2.0 - u_ctrl[6]^2.0

	c[7] = (R' * u_ctrl[7:8])[2]
	c[8] = (model.μ[end] * u_ctrl[8])^2.0 - u_ctrl[7]^2.0
end

n_ctrl_lim = 8
con_ctrl_lim = stage_constraints(control_limits_con!, n_ctrl_lim, (1:8), t_idx)

con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact, con_ctrl_comp, con_ctrl_lim])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
#open(vis)
visualize!(vis, model,
    q̄, ū,
    Δt = h)

# settransform!(vis["/Cameras/default"],
# 	compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
# setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)

plot(hcat(q̄...)', color = :red, width = 1.0, labels = "")
plot(hcat([ū..., ū[end]]...)[model.idx_u[1:8], :]', linetype = :steppost, color = :black, width = 1.0, labels = "")

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Model A
#
using ClimateMachine.DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder, NumericalFluxSecondOrder, NumericalFluxGradient

using ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient, CentralNumericalFluxSecondOrder

using LinearAlgebra: I, dot, Diagonal

using ClimateMachine.BalanceLaws:
   BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux

using ClimateMachine.VariableTemplates

import ClimateMachine.DGMethods:
     init_state_auxiliary!, update_auxiliary_state!, update_auxiliary_state_gradient!, vars_state, VerticalDirection, boundary_state!, compute_gradient_flux!, init_state_prognostic!, flux_first_order!, flux_second_order!, source!, wavespeed, compute_gradient_argument!

using StaticArrays

"""
 AModel{M} <: BalanceLaw
"""

# Create a new linear model instance
abstract type AbstractAModel <: BalanceLaw end
struct AModel{FT} <: AbstractAModel
 function AModel{FT}(
  ;
 ) where {FT<: AbstractFloat}
    return new{FT}(
    )
 end
end

function ADGModel(
   bl::AModel,
   grid,
   nfnondiff,
   nfdiff,
   gnf;
   kwargs...,
   )

   modeldata=()

   return DGModel(bl,grid,nfnondiff,nfdiff,gnf;kwargs...,modeldata=modeldata,)
end

"""
 Set model state variables and operators
"""

# State variable and initial value, just one for now, θ
##
vars_state(m::AModel, ::Prognostic, FT) = @vars(θ::FT)

init_state_prognostic!( args...) = ( init_state_cp!( args... ) )
function init_state_cp!( m::IVDCModel, Q::Vars, A::Vars, coords, t,)
  @inbounds begin
    Q.θ = 10.
  end
  return nothing
end
##

##
vars_state(m::IVDCModel, ::Auxiliary, FT) = @vars(θ_init::FT)

function init_state_auxiliary!(m::IVDCModel,A::Vars, _...)
  @inbounds begin
    A.θ_init = -0
  end
  return nothing
end
##
#

# Variables and operations used in differentiating first derivatives
##
vars_state(m::IVDCModel, ::Gradient, FT) = @vars(∇θ::FT, ∇θ_init::FT,)
@inline function compute_gradient_argument!(
    m::IVDCModel,
    G::Vars,
    Q::Vars,
    A,
    t,
)
    G.∇θ = Q.θ
    G.∇θ_init = A.θ_init

    return nothing
end
##
#

# Variables and operations used in differentiating second derivatives
##
vars_state(m::IVDCModel, ::GradientFlux, FT) = @vars(κ∇θ::SVector{3, FT})
@inline function compute_gradient_flux!(
    m::IVDCModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    κ = diffusivity_tensor(m, G.∇θ_init[3])
    D.κ∇θ = -κ * G.∇θ
    return nothing
end
##

##
## Set vertical diffusivity profile based on vertical hydrography profile
@inline function diffusivity_tensor(m::IVDCModel, ∂θ∂z)
    κᶻ = m.κᶻ
    κᶜ = m.κᶜ
    ∂θ∂z < 0 ? κ = (@SVector [0, 0, κᶜ]) : κ = (@SVector [0, 0, κᶻ])
    return Diagonal(-κ)
end
##
#

# Function to apply I to state variable

##
@inline function source!(
    m::IVDCModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t,
    direction,
)
    ivdc_dt = m.dt
    @inbounds begin
     S.θ = Q.θ/ivdc_dt
    end

    return nothing
end
##
#

# Numerical fluxes and boundaries

##
function flux_first_order!(::IVDCModel, _...) end

function flux_second_order!(
    ::IVDCModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    H::Vars,
    A::Vars,
    t,
)
    F.θ += D.κ∇θ
end

function wavespeed(m::IVDCModel, n⁻, _...)
    C = abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)
    return C
end
##

function boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient, CentralNumericalFluxGradient},
    m::IVDCModel,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    bctype,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    return nothing
end

###    From -  function numerical_boundary_flux_gradient! , DGMethods/NumericalFluxes.jl
###    boundary_state!(
###        numerical_flux,
###        balance_law,
###        state_conservative⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        aux1⁻,
###    )

function boundary_state!(
    nf::Union{NumericalFluxSecondOrder,CentralNumericalFluxSecondOrder},
    m::IVDCModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    bctype,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    D⁺.κ∇θ = n⁻ * -0
    # D⁺.κ∇θ = Q/(rho*cp)   ( at top boundary for heat flux )
    return nothing
end

###    boundary_state!(
###        numerical_flux,
###        balance_law,
###        state_conservative⁺,
###        state_gradient_flux⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_gradient_flux⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        diff1⁻,
###        aux1⁻,
###    )

###    boundary_flux_second_order!(
###        numerical_flux,
###        balance_law,
###        Grad{S}(flux),
###        state_conservative⁺,
###        state_gradient_flux⁺,
###        state_hyperdiffusive⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_gradient_flux⁻,
###        state_hyperdiffusive⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        diff1⁻,
###        aux1⁻,
###    )

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
module TestExamplesMPI

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive=true)

@testset "MPI" begin
  # TreeMesh tests
  include("test_mpi_tree.jl")

  # P4estMesh tests
  include("test_mpi_p4est.jl")
end # MPI 2D

# Clean up afterwards: delete Trixi output directory
Trixi.mpi_isroot() && @test_nowarn rm(outdir, recursive=true)

end # module


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

for (func,funcs,funch,typ) in ((:test_csrmv,:test_csrmv_symm,:test_csrmv_herm,:CudaSparseMatrixCSR),
                               (:test_cscmv,:test_cscmv_symm,:test_cscmv_herm,:CudaSparseMatrixCSC))
    @eval begin
        function $funcs(elty)
            A = sparse(rand(elty,m,m))
            A = A + A.'
            x = rand(elty,m)
            y = rand(elty,m)
            alpha = rand(elty)
            beta = rand(elty)
            d_x = CudaArray(x)
            d_y = CudaArray(y)
            d_A = Symmetric($typ(A))
            d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
            h_y = to_host(d_y)
            y = alpha * A * x + beta * y
            @test_approx_eq(y,h_y)
            x   = rand(elty,n)
            d_x = CudaArray(x)
            @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
        end

        function $funch(elty)
            A = sparse(rand(elty,m,m))
            A = A + A'
            x = rand(elty,m)
            y = rand(elty,m)
            alpha = rand(elty)
            beta = rand(elty)
            d_x = CudaArray(x)
            d_y = CudaArray(y)
            d_A = Hermitian($typ(A))
            d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
            h_y = to_host(d_y)
            y = alpha * A * x + beta * y
            @test_approx_eq(y,h_y)
            x   = rand(elty,n)
            d_x = CudaArray(x)
            @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
        end

        function $func(elty)
            A = sparse(rand(elty,m,n))
            x = rand(elty,n)
            y = rand(elty,m)
            alpha = rand(elty)
            beta = rand(elty)
            d_x = CudaArray(x)
            d_y = CudaArray(y)
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
            d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
            h_z = to_host(d_z)
            z = alpha * A * x + beta * y
            @test_approx_eq(z,h_z)
            d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
            h_z = to_host(d_z)
            z = A * x + beta * y
            @test_approx_eq(z,h_z)
            d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
            h_z = to_host(d_z)
            z = A * x + y
            @test_approx_eq(z,h_z)
            d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
            h_z = to_host(d_z)
            z = alpha * A * x + y
            @test_approx_eq(z,h_z)
            d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
            h_z = to_host(d_z)
            z = alpha * A * x
            @test_approx_eq(z,h_z)
            d_z = CUSPARSE.mv('N',d_A,d_x,'O')
            h_z = to_host(d_z)
            z = A * x
            @test_approx_eq(z,h_z)
            d_z = d_A*d_x
            h_z = to_host(d_z)
            z = A * x
            @test_approx_eq(z,h_z)
        end
    end
end

function test_bsrmv(elty)
    A = sparse(rand(elty,m,n))
    x = rand(elty,n)
    y = rand(elty,m)
    alpha = rand(elty)
    beta = rand(elty)
    d_x = CudaArray(x)
    d_y = CudaArray(y)
    d_A = CudaSparseMatrixCSR(A)
    d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
    @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
    h_z = to_host(d_z)
    z = alpha * A * x + beta * y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
    h_z = to_host(d_z)
    z = A * x + beta * y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
    h_z = to_host(d_z)
    z = A * x + y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
    h_z = to_host(d_z)
    z = alpha * A * x + y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
    h_z = to_host(d_z)
    z = alpha * A * x
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,'O')
    h_z = to_host(d_z)
    z = A * x
    @test_approx_eq(z,h_z)
    d_z = d_A*d_x
    h_z = to_host(d_z)
    z = A * x
    @test_approx_eq(z,h_z)
    @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
end

function test_hybmv(elty)
    A = sparse(rand(elty,m,n))
    x = rand(elty,n)
    y = rand(elty,m)
    alpha = rand(elty)
    beta = rand(elty)
    d_x = CudaArray(x)
    d_y = CudaArray(y)
    d_A = CudaSparseMatrixCSR(A)
    d_A = CUSPARSE.switch2hyb(d_A)
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
    h_z = to_host(d_z)
    z = alpha * A * x + beta * y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
    h_z = to_host(d_z)
    z = A * x + beta * y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
    h_z = to_host(d_z)
    z = A * x + y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
    h_z = to_host(d_z)
    z = alpha * A * x + y
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
    h_z = to_host(d_z)
    z = alpha * A * x
    @test_approx_eq(z,h_z)
    d_z = CUSPARSE.mv('N',d_A,d_x,'O')
    h_z = to_host(d_z)
    z = A * x
    @test_approx_eq(z,h_z)
    d_z = d_A*d_x
    h_z = to_host(d_z)
    z = A * x
    @test_approx_eq(z,h_z)
    @test_throws(DimensionMismatch, CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O'))
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_csrmv_symm(elty)
    test_cscmv_symm(elty)
    println("mv_symm took ", toq(), " for ", elty)
    tic()
    test_csrmv_herm(elty)
    test_cscmv_herm(elty)
    println("mv_herm took ", toq(), " for ", elty)
    tic()
    test_csrmv(elty)
    test_cscmv(elty)
    test_bsrmv(elty)
    test_hybmv(elty)
    println("mv took ", toq(), " for ", elty)
end


text:
function namedtp_to_vec( x :: NamedTuple )
    l = length(x)
    t = [ ( typeof(xx) <: AbstractFloat) ? asℝ :
           as( (typeof(xx) <: AbstractArray) ? Array : typeof(xx), size(xx) ) for xx in x]
    trans = as(NamedTuple{keys(x)}(t))
    return trans
end


text:
"""
    gridcv(X, Y; segm, score, fun, pars, verbose = false)
Cross-validation (CV) over a grid of parameters.
* `X` : X-data.
* `Y` : Y-data.
* `segm` : Segments of the CV (output of functions
     [`segmts`](@ref), [`segmkf`](@ref) etc.).
* `score` : Function computing the prediction score (= error rate; e.g. MSEP).
* `fun` : Function computing the prediction model.
* `pars` : tuple of named vectors (arguments of `fun`)
    defining the grid of parameters.
* `verbose` : If true, fitting information are printed.

The score is computed over X and Y for each combination of the
grid defined in `pars`.

The vectors in `pars` must have same length.
"""
function gridcv(X, Y; segm, score, fun, pars, verbose = false)
    q = size(Y, 2)
    nrep = length(segm)
    res_rep = list(nrep)
    nco = length(pars[1]) # nb. combinations in pars
    @inbounds for i in 1:nrep
        verbose ? print("/ rept=", i, " ") : nothing
        listsegm = segm[i]       # segments in the repetition
        nsegm = length(listsegm) # segmts: 1; segmkf: K
        zres = list(nsegm)       # results for the repetition
        @inbounds for j = 1:nsegm
            verbose ? print("segm=", j, " ") : nothing
            s = listsegm[j]
            zres[j] = gridscore(
                rmrow(X, s), rmrow(Y, s),
                X[s, :], Y[s, :];
                score = score, fun = fun, pars = pars)
        end
        zres = reduce(vcat, zres)
        dat = DataFrame(rept = fill(i, nsegm * nco),
            segm = repeat(1:nsegm, inner = nco))
        zres = hcat(dat, zres)
        res_rep[i] = zres
    end
    verbose ? println("/ End.") : nothing
    res_rep = reduce(vcat, res_rep)
    gdf = groupby(res_rep, collect(keys(pars)))
    namy = map(string, repeat(["y"], q), 1:q)
    res = combine(gdf, namy .=> mean, renamecols = false)
    (res = res, res_rep = res_rep, )
end

"""
    gridcvlv(X, Y; segm, score, fun, nlv, pars, verbose = false)
* `nlv` : Nb., or collection of nb., of latent variables (LVs).

Same as [`gridcv`](@ref) but specific to (and much faster for) models
using latent variables (e.g. PLSR).

Argument `pars` must not contain `nlv`.
"""
function gridcvlv(X, Y; segm, score, fun, nlv, pars = nothing,
        verbose = false)
    q = size(Y, 2)
    nrep = length(segm)
    res_rep = list(nrep)
    nlv = max(minimum(nlv), 0):maximum(nlv)
    le_nlv = length(nlv)
    @inbounds for i in 1:nrep
        verbose ? print("/ rept=", i, " ") : nothing
        listsegm = segm[i]       # segments in the repetition
        nsegm = length(listsegm) # segmts: 1; segmkf: K
        zres = list(nsegm)       # results for the repetition
        @inbounds for j = 1:nsegm
            verbose ? print("segm=", j, " ") : nothing
            s = listsegm[j]
            zres[j] = gridscorelv(
                rmrow(X, s), rmrow(Y, s),
                X[s, :], Y[s, :];
                score = score, fun = fun, nlv = nlv, pars = pars)
        end
        zres = reduce(vcat, zres)
        ## Case where pars is empty
        if isnothing(pars)
            dat = DataFrame(rept = fill(i, nsegm * le_nlv),
                segm = repeat(1:nsegm, inner = le_nlv))
        else
            nco = length(pars[1]) # nb. combinations in pars
            dat = DataFrame(rept = fill(i, nsegm * le_nlv * nco),
                segm = repeat(1:nsegm, inner = le_nlv * nco))
        end
        zres = hcat(dat, zres)
        res_rep[i] = zres
    end
    verbose ? println("/ End.") : nothing
    res_rep = reduce(vcat, res_rep)
    isnothing(pars) ? namgroup = [:nlv] : namgroup =  [:nlv ; collect(keys(pars))]
    gdf = groupby(res_rep, namgroup)
    namy = map(string, repeat(["y"], q), 1:q)
    res = combine(gdf, namy .=> mean, renamecols = false)
    (res = res, res_rep = res_rep, )
end

"""
    gridcvlb(X, Y; segm, score, fun, lb, pars, verbose = false)
* `nlv` : Nb., or collection of nb., of latent variables (LVs).

Same as [`gridcv`](@ref) but specific to (and much faster for) models
using ridge regularization (e.g. RR).

Argument `pars` must not contain `lb`.
"""
function gridcvlb(X, Y; segm, score, fun, lb, pars = nothing,
        verbose = false)
    q = size(Y, 2)
    nrep = length(segm)
    res_rep = list(nrep)
    lb = sort(unique(lb))
    le_lb = length(lb)
    @inbounds for i in 1:nrep
        verbose ? print("/ rept=", i, " ") : nothing
        listsegm = segm[i]       # segments in the repetition
        nsegm = length(listsegm) # segmts: 1; segmkf: K
        zres = list(nsegm)       # results for the repetition
        @inbounds for j = 1:nsegm
            verbose ? print("segm=", j, " ") : nothing
            s = listsegm[j]
            zres[j] = gridscorelb(
                rmrow(X, s), rmrow(Y, s),
                X[s, :], Y[s, :];
                score = score, fun = fun, lb = lb, pars = pars)
        end
        zres = reduce(vcat, zres)
        ## Case where pars is empty
        if isnothing(pars)
            dat = DataFrame(rept = fill(i, nsegm * le_lb),
                segm = repeat(1:nsegm, inner = le_lb))
        else
            nco = length(pars[1]) # nb. combinations in pars
            dat = DataFrame(rept = fill(i, nsegm * le_lb * nco),
                segm = repeat(1:nsegm, inner = le_lb * nco))
        end
        zres = hcat(dat, zres)
        res_rep[i] = zres
    end
    verbose ? println("/ End.") : nothing
    res_rep = reduce(vcat, res_rep)
    isnothing(pars) ? namgroup = [:lb] : namgroup =  [:lb ; collect(keys(pars))]
    gdf = groupby(res_rep, namgroup)
    namy = map(string, repeat(["y"], q), 1:q)
    res = combine(gdf, namy .=> mean, renamecols = false)
    (res = res, res_rep = res_rep, )
end

####################### Multiblock

function gridcv_mb(X, Y; segm, score, fun, pars, verbose = false)
    q = size(Y, 2)
    nrep = length(segm)
    res_rep = list(nrep)
    nco = length(pars[1]) # nb. combinations in pars
    nbl = length(X)
    @inbounds for i in 1:nrep
        verbose ? print("/ rept=", i, " ") : nothing
        listsegm = segm[i]       # segments in the repetition
        nsegm = length(listsegm) # segmts: 1; segmkf: K
        zres = list(nsegm)       # results for the repetition
        @inbounds for j = 1:nsegm
            verbose ? print("segm=", j, " ") : nothing
            s = listsegm[j]
            zX1 = list(nbl, Matrix{Float64})
            zX2 = list(nbl, Matrix{Float64})
            for k = 1:nbl
                zX1[k] = rmrow(X[k], s)
                zX2[k] = X[k][s, :]
            end
            zres[j] = gridscore(zX1, rmrow(Y, s), zX2, Y[s, :];
                score = score, fun = fun, pars = pars)
        end
        zres = reduce(vcat, zres)
        dat = DataFrame(rept = fill(i, nsegm * nco),
            segm = repeat(1:nsegm, inner = nco))
        zres = hcat(dat, zres)
        res_rep[i] = zres
    end
    verbose ? println("/ End.") : nothing
    res_rep = reduce(vcat, res_rep)
    gdf = groupby(res_rep, collect(keys(pars)))
    namy = map(string, repeat(["y"], q), 1:q)
    res = combine(gdf, namy .=> mean, renamecols = false)
    (res = res, res_rep = res_rep, )
end

function gridcvlv_mb(X, Y; segm, score, fun, nlv, pars = nothing,
        verbose = false)
    q = size(Y, 2)
    nrep = length(segm)
    res_rep = list(nrep)
    nlv = max(minimum(nlv), 0):maximum(nlv)
    le_nlv = length(nlv)
    nbl = length(X)
    @inbounds for i in 1:nrep
        verbose ? print("/ rept=", i, " ") : nothing
        listsegm = segm[i]       # segments in the repetition
        nsegm = length(listsegm) # segmts: 1; segmkf: K
        zres = list(nsegm)       # results for the repetition
        @inbounds for j = 1:nsegm
            verbose ? print("segm=", j, " ") : nothing
            s = listsegm[j]
            zX1 = list(nbl, Matrix{Float64})
            zX2 = list(nbl, Matrix{Float64})
            for k = 1:nbl
                zX1[k] = rmrow(X[k], s)
                zX2[k] = X[k][s, :]
            end
            zres[j] = gridscorelv(zX1, rmrow(Y, s), zX2, Y[s, :];
                score = score, fun = fun, nlv = nlv, pars = pars)
        end
        zres = reduce(vcat, zres)
        ## Case where pars is empty
        if isnothing(pars)
            dat = DataFrame(rept = fill(i, nsegm * le_nlv),
                segm = repeat(1:nsegm, inner = le_nlv))
        else
            nco = length(pars[1]) # nb. combinations in pars
            dat = DataFrame(rept = fill(i, nsegm * le_nlv * nco),
                segm = repeat(1:nsegm, inner = le_nlv * nco))
        end
        zres = hcat(dat, zres)
        res_rep[i] = zres
    end
    verbose ? println("/ End.") : nothing
    res_rep = reduce(vcat, res_rep)
    isnothing(pars) ? namgroup = [:nlv] : namgroup =  [:nlv ; collect(keys(pars))]
    gdf = groupby(res_rep, namgroup)
    namy = map(string, repeat(["y"], q), 1:q)
    res = combine(gdf, namy .=> mean, renamecols = false)
    (res = res, res_rep = res_rep, )
end


text:
module Sequencing

using ..Qubits # Includes Clifford and Waveforms

using InstrumentControl
using InstrumentControl.AWG5014C

export sendSequence
export sendNormalSequence

# From qubits
export Qubit
export QubitWithZ
export QubitNoZ
export Readout
export cosInit
export gaussInit
export generalInit

# From waveforms
export FloatWaveform
export ExactWaveform

# From clifford
export benchmark1Qubit
export benchmark2Qubit
export Pulse

#250 for Hardware Sequencer
const NORMSEQLENGTH = 400
lastSeqNormal = false

# gateNames has its last use here in naming the AWG pulse - no parents need it

#============================ Basic seqeuncing ================================#
# The sequencing also has to be different depending on which type
# of waveform we are using.  In the case of FloatWaveforms, we need to insert
# idle pulses regularly, which are all the same point in the DAC's memory.
# In the case of ExactWaveforms, we need to access distinct points in the DAC's
# memory for each combination of (line, pulse).  The code allows
# for the possibility that the qubit's lines are not all on the same board.

# In preparing for 1-qubit sequencing, the qubit ensures that all pulses in its
# dictionary are up to date, and updates the 10x3 matrix of indices.
# These map the 10 basic single qubit pulses to their corresponding on-board
# waveform index.  That means that by calling q.pulseConvert[seq, n], with some
# sequence generated elsewhere in the code, we get the actual index list for
# sequencing on the nth line.
function prepForSeq(q::Qubit)
  # Communicate with DAC only if necessary
  for p in q.waveforms
    idx = p[1][1]
    if p[2].undefined
      write(q.lineXYI[1], "WLIST:WAV:NEW "*quoted(gateNames[idx]*"I")*", "*string(AWGLENGTH)*", INTEGER")
      write(q.lineXYQ[1], "WLIST:WAV:NEW "*quoted(gateNames[idx]*"Q")*", "*string(AWGLENGTH)*", INTEGER")
      if isa(q, QubitWithZ)
        write(q.lineZ[1], "WLIST:WAV:NEW "*quoted(gateNames[idx])*", "*string(AWGLENGTH)*", INTEGER")
      end
      p[2].undefined = false
    end
    if p[2].dirty
      q.pulseConvert[idx,1] = pushWaveform(p[2].XYI, q.lineXYI[1], gateNames[idx]*"I")
      q.pulseConvert[idx,2] = pushWaveform(p[2].XYQ, q.lineXYQ[1], gateNames[idx]*"Q")
      if isa(q, QubitWithZ)
        q.pulseConvert[idx,3] = pushWaveform(p[2].Z, q.lineZ[1], gateNames[idx])
      end
      p[2].dirty = false
    end
  end
end

function prepForRO(r::Readout)
  if r.ROIdle.undefined
    write(r.lineXYI[1], "WLIST:WAV:NEW \"ROIdle\", "*string(length(r.ROIdle.XYI))*", INTEGER")
    r.ROIdle.undefined = false
  end
  if r.ROWaveform.undefined
    write(r.lineXYI[1], "WLIST:WAV:NEW \"ROI\", "*string(length(r.ROIdle.XYI))*", INTEGER")
    write(r.lineXYQ[1], "WLIST:WAV:NEW \"ROQ\", "*string(length(r.ROIdle.XYI))*", INTEGER")
    r.ROWaveform.undefined = false
  end
  if r.ROIdle.dirty
    pushWaveform(r.ROIdle.XYI, r.lineXYI[1], "ROIdle")
    r.ROIdle.dirty = false
  end
  if r.ROWaveform.dirty
    pushWaveform(r.ROWaveform.XYI, r.lineXYI[1], "ROI")
    pushWaveform(r.ROWaveform.XYQ, r.lineXYQ[1], "ROQ")
    r.ROWaveform.dirty = false
  end
end

# Return the label of the waveform in memory
function pushWaveform(wavedata::Vector{UInt16}, board::Instrument, name="")
  #if length(wavedata) == 0
    # This is one of the off channels of a FloatWaveform.  Just return the
    # stored index of idle pulse
    #return idles[board]
  #end
  if isa(board, InsAWG5014C)
    return awgPush(wavedata, board, name)
  else
    error("Instrument not supported yet, whine at a programmer")
  end
end

# We will be using the VISA protocol and binary writing.  We will also force
# every waveform to 250 points, with front-padded zeros.  A 30-point waveform
# will then have 220 points of offsetValue and then 30 points with signal.
import InstrumentControl.AWG5014C: offsetValue
function awgPush(wavedata::Vector{UInt16}, ins::InsAWG5014C, name::String)
  header = *("WLIST:WAV:DATA ", quoted(name), ",")
  if length(wavedata) < AWGLENGTH
    pad = div(AWGLENGTH - length(wavedata), 2)
    binblockwrite(ins, header, reinterpret(UInt8,
        htol.([fill(offsetValue, AWGLENGTH - length(wavedata) - pad); wavedata; fill(offsetValue, pad)])))
  else
    binblockwrite(ins, header, reinterpret(UInt8,
        htol.(wavedata)))
  end
  name # Return the name, how the AWG will refer to it later.
end

# The following functions take pulse sequences in the form produced by the code,
# and send them to the DACs based on their indices in the DAC memory.  This
# means that sendSequence(q, [Xpi2, I, Xpi2]) for example, will send XYI, XYQ,
# and Z pulses according to the prior definitions of these pulses.
function sendSequence(q::Qubit, sequence::Pulse)
  if size(sequence, 2) == 1
    sendSequence(q, sequence[:])
  else
    error("Please specify only one sequence\n")
  end
end

function sendSequence(q::Qubit, sequence::Vector{Int8})
  # Make sure we aren't going to send obsolete pulses
  # (no writes done if everything is current)
  prepForSeq(q)

  # TODO: More generalized version of this function.  For now, we assume the
  # qubit is controlled by the AWG, that it is a QubitNoZ object, and that
  # the other AWG lines need to idle until the sequence terminates, and then
  # perform a readout pulse.  This allows the sequencing to be done in one
  # (albeit large) command to the AWG.  If this isn't the goal, bug Brett.
  # Someday readout will be a single pulse for multiple qubits.
  ins = q.lineXYI[1]
  mess = "AWGC:RMOD SEQ\nSEQ:LENG 0\nSEQ:LENG "*string(length(sequence))*"\n"

  mess *= map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYI[2])*
                                                    " \""*string(y)*"\"\n",
                            1:length(sequence), q.pulseConvert[sequence, 1])[:]...
  mess *= map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYQ[2])*
                                                    " \""*string(y)*"\"\n",
                            1:length(sequence), q.pulseConvert[sequence, 2])[:]...
  # Get the other channels and instruct them to be idle at this time
  a = 1:4
  a = a[(a.!=q.lineXYI[2])&(a.!=q.lineXYQ[2])] # Select channels not used by XY
  mess *= map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(a[1])*" \"IdleI\"\n", 1:length(sequence))...
  mess *= map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(a[2])*" \"IdleQ\"\n", 1:length(sequence))...
  # Put it all together and send to the AWG, and turn on output
  write(ins, mess)
  println(errors(ins))
  ins[SequenceGOTOTarget,1] = 1
  ins[SequenceGOTOState, length(sequence)] = true
  ins[SequenceWaitTrigger,1] = true
  @allch ins[ChannelOutput] = true
  ins[Output] = true
end

function sendSequence(q::Qubit, sequence::Pulse, r::Readout)
  if size(sequence, 2) == 1
    sendSequence(q, sequence[:], r)
  else
    error("Please specify only one sequence\n")
  end
end

# The flag segnifies whether to make the assumptions in the function below this.
function sendSequence(q::Qubit, sequence::Vector{Int8}, r::Readout, flag::Bool = false)
  # Make sure we aren't going to send obsolete pulses
  # (no writes done if everything is current)
  ins = q.lineXYI[1]
  ins[Output] = false
  prepForSeq(q)
  prepForRO(r)

  # TODO: More generalized version of this function.  For now, we assume the
  # qubit is controlled by the AWG, that it is a QubitNoZ object, and that
  # the other AWG lines need to idle until the sequence terminates, and then
  # perform a readout pulse.  This allows the sequencing to be done in one
  # (albeit large) command to the AWG.  If this isn't the goal, bug Brett.
  # Someday readout will be a single pulse for multiple qubits.
  mess = flag ? "" : "AWGC:RMOD SEQ\nSEQ:LENG 0\nSEQ:LENG "*string(length(sequence) + 1)*"\n"

  mess *= map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYI[2])*
                                                    " \""*string(y)*"\"\n",
                            1:length(sequence), q.pulseConvert[sequence, 1])[:]...
  mess *= map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYQ[2])*
                                                    " \""*string(y)*"\"\n",
                            1:length(sequence), q.pulseConvert[sequence, 2])[:]...
  # Instruct readout channels to idle
  if !flag
      mess *= map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(r.lineXYI[2])*" \"IdleI\"\n", 1:length(sequence))...
      mess *= map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(r.lineXYQ[2])*" \"IdleQ\"\n", 1:length(sequence))...

      # Readout via a defined readout pulse
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":WAV"*string(q.lineXYI[2])*" \"ROIdle\"\n"
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":WAV"*string(q.lineXYQ[2])*" \"ROIdle\"\n"
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":WAV"*string(r.lineXYI[2])*" \"ROI\"\n"
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":WAV"*string(r.lineXYQ[2])*" \"ROQ\"\n"

      # Ensure gotos and triggers are correct
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":GOTO:STAT 1\n"
      mess *= "SEQ:ELEM"*string(length(sequence)+1)*":GOTO:IND 1\n"
      mess *= "SEQ:ELEM1:TWA 1\n"
      global lastSeqNormal
      lastSeqNormal = false # overwritten later if the sequence is normal
  end
  write(ins, mess)
  @allch ins[ChannelOutput] = true
  ins[Output] = true
end

# AWG-friendly sequencer that makes sequences all have a standard form.
# Sequences are assumed to have NORMSEQLENGTH + 1 elements, NORMSEQLENGTH of
# which are length-AWGLENGTH pulses, and the final a ~1e5 point readout.  The
# plus side is that the readout never needs to be changed between sequence
# definitions, nor do the gotos.  Sequences shorter than NORMSEQLENGTH gates
# (excluding readout) are front padded with idles.
import ..Clifford: Idle
function sendNormalSequence(q::Qubit, sequence::Pulse, r::Readout)
    assert(length(sequence) < 1 + NORMSEQLENGTH)
    nseq = [Idle*(NORMSEQLENGTH-length(sequence)); sequence]
    global lastSeqNormal # cache the state of the sequencer
    sendSequence(q, nseq[:], r, lastSeqNormal)
    lastSeqNormal = true
end

function sendMultiROSequence(q::Qubit, sequences::Vector{Pulse}, r::Readout)
    ins = q.lineXYI[1]
    ins[Output] = false
    prepForSeq(q)
    prepForRO(r)

#=
function sendSequence(q::Qubit, sequence::Vector{Int8})
  # Make sure we aren't going to send obsolete pulses
  # (no writes done if everything is current)
  prepForSeq(q)
  # Call a lower level function for each line, using indices on the board
  # that correspond to the sequence requested
  seqLowLevel(q.pulseConvert[sequence, 1], q.lineXYI)
  seqLowLevel(q.pulseConvert[sequence, 2], q.lineXYQ)
  if isa(q, QubitWithZ)
    seqLowLevel(q.pulseConvert[sequence, 3], q.lineZ)
  end
  # Turn on output where along this path?

  # How to manage readout?  Do we send a trigger pulse along with the last
  # pulse in the seqeunce, or do we calculate the length of the pulse sequence
  # and just readout based on timing?  Also i expect not all sequences to be
  # sent are necessarily to be followed by readout, so I guess readout-
  # terminated sequences get their own method.
end
=#
# Nothing here covers readout or two-qubit gates.  Both are heirarchically
# lateral from the qubits themselves - pairs of coupled qubits need to know
# waveforms for 2-qubit gates, but not the qubits themselves.  The readout
# process covers all qubits connected to the readout line.  Whatever controls
# the readout will need to be able to tell which qubit is which, but the qubits
# themselves may remain unaware of the readout process.

end # End module sequencing

max_stars_repo_name: PainterQubits/PulseController

max_issues_repo_name: PainterQubits/PulseController

max_forks_repo_name: PainterQubits/PulseController


text:
num_johnson_subsets(simplex_length::Integer) = 2^simplex_length - 1

"""
Compute all subsets of the points in the simplex in a reliable order. The
order is arbitrary, but was chosen to match the implemenation in
S. Cameron, “Enhancing GJK: computing minimum and penetration distances between
convex polyhedra,”. Specifically, subset i contains point j iff the binary
representation of i has a one at the jth bit.
"""
function johnson_subsets(simplex_length::Integer)
    num_subsets = num_johnson_subsets(simplex_length)
    subsets = falses(simplex_length, num_subsets)
    for i in 1:num_subsets
        for j in 1:simplex_length
            subsets[j, i] = (i >> (j - 1)) & 1
        end
    end
    subsets
end

"""
    weights = projection_weights(simplex)

This function implements Johnson's distance subalgorithm, as described in
E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for
computing the distance between complex objects in three-dimensional space,”
1988. Given a simplex (a length N+1 vector of points of dimension N), it
returns weights such that dot(weights, simplex) yields the point in the convex
hull of the simplex which is closest to the origin.

This is the critical loop of the GJK algorithm, so it has been heavily optimized
to precompute, inline, and unroll as much of the algorithm as possible. For a
more readable (and much slower) implementation, see
projection_weights_reference()
"""
@generated function projection_weights(simplex::SVector{M, SVector{N, T}}; atol=sqrt(eps(T))) where {M,N,T}
    simplex_length = M
    num_subsets = num_johnson_subsets(M)
    subsets = johnson_subsets(M)
    complements = .!subsets

    expr = quote
        deltas = SVector(tuple($([zero(SVector{simplex_length, T}) for i = 1 : num_subsets]...)))
    end

    # Set the weight of every singleton subset to 1.
    for i in 1:simplex_length
        elements = [:(zero($T)) for j in 1:simplex_length]
        elements[i] = :(one($T))
        arg = :(SVector{$simplex_length, $T}(tuple($(elements...))))

        push!(expr.args, quote
            deltas = setindex(deltas, $(arg), $(2^(i - 1)))
        end)
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        push!(expr.args, quote
            viable = true
        end)
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            push!(expr.args, quote
                d = deltas[$s2][$j]
            end)
            for i in (1:M)[subsets[:,s]]
                push!(expr.args, quote
                    d += deltas[$s][$i] *
                    (dot(simplex[$i], simplex[$k]) - dot(simplex[$i], simplex[$j]))
                end)
            end
            push!(expr.args, quote
                deltas = setindex(deltas, setindex(deltas[$s2], d, $j), $s2)
            end)
            push!(expr.args, quote
                if d > atol # See #17
                    viable = false
                end
            end)
        end
        push!(expr.args, quote
            viable = viable && all($(Expr(:tuple,
            [:(deltas[$s][$i] >= 0) for i in (1:M)[subsets[:,s]]]...)))
            if viable
                return deltas[$s] ./ sum(deltas[$s])
            end
        end)
    end
    push!(expr.args, quote
        return deltas[$num_subsets] ./ sum(deltas[$num_subsets])
    end)
    return expr
end

function projection_weights_reference(simplex::SVector{M, SVector{N, T}}) where {M,N,T}
    num_subsets = 2^M - 1
    complements = falses(M, num_subsets)
    subsets = falses(M, num_subsets)
    for i in 1:num_subsets
        for j in 1:M
            subsets[j, i] = (i >> (j - 1)) & 1
            complements[j, i] = !subsets[j, i]
        end
    end

    deltas = zeros(MMatrix{M, num_subsets, T})
    viable = ones(MVector{num_subsets, Bool})

    for i in 1:M
        deltas[i, 2^(i - 1)] = 1
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            delta_j_s2 = 0
            for i in (1:M)[subsets[:,s]]
                delta_j_s2 += deltas[i, s] * (dot(simplex[i], simplex[k]) - dot(simplex[i], simplex[j]))
            end
            if delta_j_s2 > 0
                viable[s] = false
            elseif delta_j_s2 < 0
                viable[s2] = false
            end
            deltas[j, s2] = delta_j_s2
        end
        if viable[s]
            return deltas[:,s] ./ sum(deltas[:,s])
        end
    end
    return deltas[:,end] ./ sum(deltas[:,end])
end


text:
# Integrations with AtomsBase.jl

abstract type AbstractAtom end

"""
    DynamicAtom{D,L<:Unitful.Length,V<:Unitful.Velocity} <: AbstractAtom

An atom representation based on the StaticAtom but with a velocity.

**Type parameters**
- `D`: the dimension of the coordinate space
- `L`: the type for the position components
- `V`: the type for the velocity components

Fields should not be accessed directly. Use the provided accessors instead.
"""
struct DynamicAtom{D,L<:Unitful.Length,V<:Unitful.Velocity} <: AbstractAtom
    position::SVector{D,L}
    velocity::SVector{D,V}
    element::Element
end
function DynamicAtom(position, velocity, element)
    DynamicAtom{length(position),eltype(position),eltype(velocity)}(position, velocity, element)
end
function DynamicAtom(position, velocity, symbol::Union{Integer,AbstractString,Symbol})
    DynamicAtom(position, velocity, elements(symbol))
end

AtomsBase.position(atom::DynamicAtom) = atom.position
AtomsBase.velocity(atom::DynamicAtom) = atom.velocity
AtomsBase.species(atom::DynamicAtom) = atom.element

AtomsBase.atomic_symbol(a::DynamicAtom) = a.element.symbol
AtomsBase.atomic_mass(a::DynamicAtom) = a.element.atomic_mass
AtomsBase.atomic_number(a::DynamicAtom) = a.element.number
AtomsBase.atomic_property(a::DynamicAtom, property::Symbol) = getproperty(a.element, property)

"""
    DynamicSystem{D,AT<:AbstractAtom,L<:Unitful.Length,TT<:Unitful.Time} <: AbstractAtomicSystem{D}

A representation of a system of dynamic atoms which is similar to the FlexibleSystem but with a time field.

**Type parameters**
- `D`: the dimension of the coordinate space
- `A`: the type for the atoms that make up the system
- `L`: the type for the bounding box components
- `T`: the type for time field

Fields should not be accessed directly. Use the provided accessors instead.
"""
struct DynamicSystem{D,A<:AbstractAtom,L<:Unitful.Length,T<:Unitful.Time} <: AbstractAtomicSystem{D}
    box::SVector{D,SVector{D,L}}
    boundary_conditions::SVector{D,BoundaryCondition}
    particles::Vector{A}
    time::T
end
function DynamicSystem(box, boundary_conditions, particles, time)
    D = length(box)
    A = eltype(particles)
    L = eltype(eltype(box))
    T = typeof(time)

    DynamicSystem{D,A,L,T}(box, boundary_conditions, particles, time)
end
function DynamicSystem(particles::Vector{<:AbstractAtom}, box_size::Unitful.Length, time::Unitful.Time = 0.0u"s")
    z = zero(typeof(box_size))
    box = SVector(SVector(box_size, z, z), SVector(z, box_size, z), SVector(z, z, box_size))
    boundary_conditions = SVector(Periodic(), Periodic(), Periodic())
    DynamicSystem(box, boundary_conditions, particles, time)
end

AtomsBase.bounding_box(sys::DynamicSystem) = sys.box
AtomsBase.boundary_conditions(sys::DynamicSystem) = sys.boundary_conditions

Base.size(sys::DynamicSystem) = size(sys.particles)
Base.length(sys::DynamicSystem) = length(sys.particles)
Base.getindex(sys::DynamicSystem, i::Int) = getindex(sys.particles, i)

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
immutable PersistentArrayMap{K, V} <: PersistentMap{K, V}
    kvs::Vector{KVPair{K, V}}

    PersistentArrayMap(kvs::Vector{KVPair{K, V}}) = new(kvs)
    PersistentArrayMap() = new(KVPair{K, V}[])
end
PersistentArrayMap{K, V}(kvs::(@compat Tuple{K, V})...) =
    PersistentArrayMap{K, V}(KVPair{K, V}[KVPair(k, v) for (k, v) in kvs])
PersistentArrayMap(; kwargs...) = PersistentArrayMap(kwargs...)

Base.isequal(m1::PersistentArrayMap, m2::PersistentArrayMap) =
    isequal(Set(m1.kvs), Set(m2.kvs))
==(m1::PersistentArrayMap, m2::PersistentArrayMap) =
    Set(m1.kvs) == Set(m2.kvs)

Base.length(m::PersistentArrayMap)  = length(m.kvs)
Base.isempty(m::PersistentArrayMap) = length(m) == 0

findkeyidx(m::PersistentArrayMap, k) = findfirst((kv) -> kv.key == k, m.kvs)

function _get(m::PersistentArrayMap, k, default, hasdefault::Bool)
    for kv in m.kvs
        kv.key == k && return kv.value
    end
    hasdefault ? default : default()
end

Base.get(m::PersistentArrayMap, k) =
    _get(m, k, ()->error("key not found: $k"), false)
Base.get(m::PersistentArrayMap, k, default) =
    _get(m, k, default, true)
Base.getindex(m::PersistentArrayMap, k) = get(m, k)

Base.haskey(m::PersistentArrayMap, k) = get(m, k, NotFound()) != NotFound()

function assoc{K, V}(m::PersistentArrayMap{K, V}, k, v)
    idx = findkeyidx(m, k)
    idx == 0 && return PersistentArrayMap{K, V}(push!(m.kvs[1:end], KVPair(k, v)))

    kvs = m.kvs[1:end]
    kvs[idx] = KVPair(k, v)
    PersistentArrayMap{K, V}(kvs)
end

function dissoc{K, V}(m::PersistentArrayMap{K, V}, k)
    idx = findkeyidx(m, k)
    idx == 0 && return m

    kvs = m.kvs[1:end]
    splice!(kvs, idx)
    PersistentArrayMap{K, V}(kvs)
end

Base.start(m::PersistentArrayMap)   = 1
Base.done(m::PersistentArrayMap, i) = i > length(m)
Base.next(m::PersistentArrayMap, i) = (convert(Tuple, m.kvs[i]), i+1)

Base.map(f::(@compat Union{DataType, Function}), m::PersistentArrayMap) =
    PersistentArrayMap([f(kv) for kv in m]...)

@compat Base.show{K, V}(io::IO, ::MIME"text/plain", m::PersistentArrayMap{K, V}) =
    print(io, "Persistent{$K, $V}$(m.kvs)")


text:
# From https://github.com/unkcpz/LibSymspg.jl/blob/f342e72/test/runtests.jl#L83-L85
@testset "Niggli reduce" begin
    lattice = [
        4.0 0.0 0.0
        20.0 2.0 0.0
        0.0 0.0 12.0
    ]
    @test niggli_reduce(lattice, 1e-3) ≈ [
        0.0 4.0 0.0
        -2.0 0.0 0.0
        0.0 0.0 12.0
    ]
end
# From https://github.com/unkcpz/LibSymspg.jl/blob/f342e72/test/runtests.jl#L87-89
@testset "Delaunay reduce" begin
    lattice = [
        4.0 0.0 0.0
        20.0 2.0 0.0
        0.0 0.0 12.0
    ]
    @test delaunay_reduce(lattice, 1e-3) ≈ [
        0.0 -4.0 0.0
        2.0 0.0 0.0
        0.0 0.0 12.0
    ]
end


text:
using QuantumCumulants
using Test

@testset "cluster" begin

order = 4
N_c = 3 #number of clusters
N = [Parameter(Symbol(:N_, i)) for i=1:N_c]

hf = FockSpace(:cavity)
ha = [NLevelSpace(Symbol(:atoms, j),3) for j=1:N_c]
ha_c = [ClusterSpace(ha[j],N[j],order) for j=1:N_c]
h = ⊗(hf, ha_c...)
# Define the fundamental operators
a = Destroy(h,:a,1)
S(i,j,c) = Transition(h,Symbol(:σ, c),i, j, 1+c) #c=cluster

@test QuantumCumulants.has_cluster(h)
@test !(QuantumCumulants.has_cluster(hf))

@test S(2,2,1) ≠ S(2,2,2) ≠ S(2,2,3)
@test isequal(S(2,2,2)[1]*S(2,2,1)[2], S(2,2,1)[2]*S(2,2,2)[1])
@test length(S(2,2,1)) == length(S(2,2,2)) == length(S(2,2,3)) == order
@test acts_on(S(2,2,1)[1]) == QuantumCumulants.ClusterAon(2,1)
@test acts_on(S(2,2,2)[2]) == QuantumCumulants.ClusterAon(3,2)
@test acts_on(S(2,2,1)[1]) < acts_on(S(2,2,1)[2]) < acts_on(S(2,2,2)[1])

end #testset


text:
export has

using IntervalSets

function Base.write(db::SQLite.DB, ::Type{Material}, mat::Material; atol=1.0e-4)::Int
    res = find(db, Material, mat.name)
    if res ≠ -1
        dbmat = read(db, Material, res)
        if !isapprox(mat, dbmat, atol=atol)
            error("$(mat.name) has already been defined as $(mat).")
        end
    else
        stmt1 = SQLite.Stmt(db, "INSERT INTO MATERIAL (MATNAME, MATDESCRIPTION, MATDENSITY, MATPEDIGREE) VALUES ( ?, ?, ?, ? );")
        r = DBInterface.execute(stmt1, (name(mat), get(mat, :Description, ""), get(mat, :Density, missing), get(mat, :Pedigree, missing)))
        res = DBInterface.lastrowid(r)
        stmt2 = SQLite.Stmt(db, "SELECT PKEY FROM MATERIAL WHERE MATNAME=?;")
        pkey = (DBInterface.execute(stmt2, (name(mat), )) |> DataFrame)[end,:PKEY]
        stmt3 = SQLite.Stmt(db, "INSERT INTO MASSFRACTION ( MATKEY, MFZ, MFC, MFUC, MFA ) VALUES ( ?, ?, ?, ?, ? )")
        foreach(elm->DBInterface.execute(stmt3, (pkey, z(elm), value(mat[elm]), σ(mat[elm]), get(mat.a, elm, 0.0))), keys(mat))
    end
    return res
end

Base.write(db::SQLite.DB, mat::Material; atol=1.0e-4) = write(db, Material, mat, atol=atol)

function Base.read(db::SQLite.DB, ::Type{Material}, pkey::Int)::Material
    stmt1 = SQLite.Stmt(db, "SELECT * FROM MATERIAL WHERE PKEY=?;")
    q1 = DBInterface.execute(stmt1, (pkey, ))
    if SQLite.done(q1)
        error("No known material with pkey = '$(pkey)'.")
    end
    r1 = SQLite.Row(q1)
    row, den, desc = r1[:PKEY], r1[:MATDENSITY], r1[:MATDESCRIPTION]
    stmt2 = SQLite.Stmt(db, "SELECT * FROM MASSFRACTION WHERE MATKEY=?;")
    q2 = DBInterface.execute(stmt2, (row, ))
    massfrac, aa = Dict{Int,UncertainValue}(), Dict{Int,Float64}()
    for r2 in q2
        @assert r2[:MATKEY]==row
        z, c, uc, a = r2[:MFZ], r2[:MFC], r2[:MFUC], r2[:MFA]
        massfrac[z] = uv(c,uc)
        a>0.0 && (aa[z]=a)
    end
    props=Dict{Symbol,Any}()
    if !ismissing(den)
        props[:Density]=den
    end
    if (!ismissing(desc)) && (lastindex(desc)>=1)
        props[:Description]=desc
    end
    mat = Material(r1[:MATNAME], massfrac, aa, props)
    mat[:Database] = db
    return mat
end

function find(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Int
    stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM MATERIAL WHERE MATNAME=?;")
    q1 = DBInterface.execute(stmt1, (matname, ))
    return SQLite.done(q1) ? -1 : SQLite.Row(q1)[:PKEY]
end

NeXLCore.has(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Bool =
    find(db, Material, matname) ≠ -1

Base.read(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Material =
    read(db, Material, find(db, Material, matname))

function NeXLCore.material(db::SQLite.DB, matname::AbstractString)
    idx = find(db, Material, matname)
    if idx ≠ -1
        return read(db, Material, idx)
    else
        error("Unable to find $(matname) in the material database.")
    end
end

function Base.delete!(db::SQLite.DB, ::Type{Material}, matname::AbstractString)
    stmt1 = SQLite.Stmt(db, "SELECT PKEY, * FROM MATERIAL WHERE MATNAME=?;")
    q1 = DBInterface.execute(stmt1, (matname, ))
    for r1 in q1
        SQLite.transaction(db) do
            stmt1 = SQLite.Stmt(db, "DELETE FROM MASSFRACTION where MATKEY=?;")
            DBInterface.execute(stmt1, (r1[:PKEY], ))
            stmt2 = SQLite.Stmt(db, "DELETE FROM MATERIAL where PKEY=?;")
            DBInterface.execute(stmt2, (r1[:PKEY], ))
        end
    end
end

Base.filter(db::SQLite.DB, ::Type{Material}, prs::Pair{Element, ClosedInterval{Float64}}...)::Vector{<:Material} =
    filter(db, Material, Dict(prs...))

function Base.filter(db::SQLite.DB, ::Type{Material}, filt::Dict{Element, ClosedInterval{Float64}})::Vector{<:Material}
    cmds, args = String[], Any[]
    for (elm, ci) in filt
        push!(cmds, "SELECT MATKEY FROM MASSFRACTION WHERE MFZ=? AND MFC>=? and MFC<=?")
        append!(args, [ z(elm), minimum(ci), maximum(ci) ])
    end
    stmt = SQLite.Stmt(db, join(cmdsol," INTERSECT ")*";")
    q = DBInterface.execute(stmt, args)
    return [ read(db, Material, r[:MATKEY]) for r in q ]
end

function Base.similar(db::SQLite.DB, mat::Material, tol::Float64)::Vector{<:Material}
    interval(el) = max(0.0,mat[el]-tol)..(mat[el]+tol)
    return filter(db, Material, Dict(elm=>interval(elm) for elm in keys(mat)))
end


text:
using MathOptInterface
using JuMP
using SCIP

const MOI = MathOptInterface

export ɛ, settimelimit!, stilltime, @cb_constraint

# tolerance for numerical comparison
const ɛ = 1e-6

"numerically (strictly) nonzero"
nonzero(x::Real; atol=ɛ) = abs(x) > atol

"numerically less than"
isapproxlt(x, y; atol=ɛ) = isapprox(x, y, atol=atol) || x < y

"numerically more than"
isapproxgt(x, y; atol=ɛ) = isapprox(x, y, atol=atol) || x > y

"check that value is within bounds"
function isapproxin(x, lb, ub; atol=ɛ)
    isapproxgt(x, lb, atol=atol) && isapproxlt(x, ub, atol=atol)
end

"check the special ordered set (type 1) property"
function isSOS1(xs; atol=ɛ)
    all([isapproxgt(x, 0.0, atol=atol) for x in xs]) || return false
    sum(xs .>= atol) <= 1
end

"check the special ordered set (type 2) property"
function isSOS2(xs; atol=ɛ)
    all([isapproxgt(x, 0.0, atol=atol) for x in xs]) || return false
    nonzeros = xs .>= atol
    sum(nonzeros) <= 1 || sum(nonzeros) == 2 &&
        (sum(nonzeros[1:end-1] .& nonzeros[2:end]) == 1)
end

const timebuffer = 10.0 # seconds

"update timelimit (in seconds) for internal solver"
function settimelimit!(model::JuMP.Model, solver, limit)
    limit = max(limit, timebuffer) # at least buffer
    if limit < Inf
        MOI.set(JuMP.backend(model), MOI.TimeLimitSec, limit)
    end
end

"is there still enough until final limit?"
function stilltime(finaltime; buffer=timebuffer)
    time() + buffer < finaltime
end

"""Extract solution values for given JuMP variables, for use in callbacks."""
function SCIP.sol_values(o::SCIP.Optimizer,
                         vars::AbstractArray{JuMP.VariableRef},
                         sol::Ptr{SCIP.SCIP_SOL}=C_NULL)
    return SCIP.sol_values(o, [JuMP.index(v) for v in vars], sol)
end

"""
Add JuMP constraint expression to MOI optimizer, converting the types.

This is useful for solver-specific callbacks, e.g, with SCIP.
"""
macro cb_constraint(optimizer, expr)
    code = quote
        jump_cons = JuMP.@build_constraint $expr
        MOI.add_constraint(
            $optimizer,
            JuMP.moi_function(jump_cons.func),
            jump_cons.set
        )
    end
    quote
        let
            $(esc(code))
        end
    end
end


text:
using HDF5
using Test
using Pkg
filter_path = joinpath(dirname(pathof(HDF5)), "..", "filters")
Pkg.develop(PackageSpec(path=joinpath(filter_path, "H5Zblosc")))
Pkg.develop(PackageSpec(path=joinpath(filter_path, "H5Zbzip2")))
Pkg.develop(PackageSpec(path=joinpath(filter_path, "H5Zlz4")))
Pkg.develop(PackageSpec(path=joinpath(filter_path, "H5Zzstd")))

@info "libhdf5 v$(HDF5.API.h5_get_libversion())"

# To debug HDF5.jl tests, uncomment the next line
# ENV["JULIA_DEBUG"] = "Main"

@testset "HDF5.jl" begin

@debug "plain"
include("plain.jl")
@debug "compound"
include("compound.jl")
@debug "custom"
include("custom.jl")
@debug "reference"
include("reference.jl")
@debug "dataspace"
include("dataspace.jl")
@debug "hyperslab"
include("hyperslab.jl")
@debug "readremote"
include("readremote.jl")
@debug "extend_test"
include("extend_test.jl")
@debug "gc"
include("gc.jl")
@debug "external"
include("external.jl")
@debug "swmr"
include("swmr.jl")
@debug "mmap"
include("mmap.jl")
@debug "properties"
include("properties.jl")
@debug "table"
include("table.jl")
@debug "filter"
include("filter.jl")
@debug "chunkstorage"
include("chunkstorage.jl")
@debug "fileio"
include("fileio.jl")
@debug "filter test utils"
include("filters/FilterTestUtils.jl")

using MPI
if HDF5.has_parallel()
    # basic MPI tests, for actual parallel tests we need to run in MPI mode
    include("mpio.jl")
end

# Clean up after all resources
HDF5.API.h5_close()

end


text:
"""
# AccelerometerCalibrationPlots.jl

Debug plots for
[AccelerometerCalibration.jl](https://github.com/notinaboat/AccelerometerCalibration.jl)
"""
module AccelerometerCalibrationPlots

using AccelerometerCalibration
using Plots

offset_series = []
scale_series = []
rotation_series = []

function reset()
    empty!(offset_series)
    empty!(scale_series)
    empty!(rotation_series)
end

macro circle_plot(x, y)
    quote
        plot( x->sqrt(1 - x^2), color=:black, legend=false);
        plot!(x->-sqrt(1 - x^2), color=:black, legend=false);
        scatter!($(esc(x)),
                 $(esc(y)),
                 xaxis=(label=$(string(x))),
                 yaxis=(label=$(string(y))),
                 xlims=(-1.5,1.5),
                 ylims=(-1.5,1.5),
                 legend=false)
    end
end

macro series_plot(x)
    quote
        plot($(esc(x)), legend=false, xaxis=(label=$(string(x))))
    end
end

function calplot(x, y, z, xc, yc, zc,
                 offset, scale, rotation)
     plot(@circle_plot(x, y), @circle_plot(xc, yc), @series_plot(offset),
          @circle_plot(x, z), @circle_plot(xc, zc), @series_plot(scale),
          @circle_plot(y, z), @circle_plot(yc, zc), @series_plot(rotation),
          layout = (3,3), size=(1200,1200))
end

function calplot(c::AbstractArray{AccelerometerCalibration.Calibration})

    c_count = length(c)
    l = length(c[1].points)

    xyz = [[sc.points[j][i] for i in 1:3, j in 1:l] for sc in c]
    xyzc =[sc.rotation * ((xyz[i] .* sc.scale) .- sc.offset)
           for (i, sc) in enumerate(c)]

    x, y, z = (hcat((xyz[i][j,:] for i in 1:c_count)...) for j in 1:3)
    xc, yc, zc = (hcat((xyzc[i][j,:] for i in 1:c_count)...) for j in 1:3)

    push!(offset_series, vec([c[j].offset[i] for i in 1:3, j in 1:c_count]))
    push!(scale_series, vec([c[j].scale[i] for i in 1:3, j in 1:c_count]))
    push!(rotation_series, vcat([[r.theta1, r.theta2, r.theta3]
                               for r in (c[j].rotation for j in 1:c_count)]...))
    o = permutedims(hcat(offset_series...))
    s = permutedims(hcat(scale_series...))
    r = permutedims(hcat(rotation_series...))

    calplot(x, y, z, xc, yc, zc, o, s, r)
end

end # module


text:
module MianalyzerBackend

function allJuliaFiles(dir)
  contents = String[]
  for (root, dirs, files) in walkdir(dir)
    push!.(Ref(contents), joinpath.(root, filter(x -> splitext(x)[2] == ".jl", files)))
  end
  return contents
end

srcdir = @__DIR__
model_dependency_orders = map(["Imports.jl", "Stats.jl"]) do f
  srcdir * "/models/" * f
end
for file in model_dependency_orders
  println(file)
  include(file)
end

for file in allJuliaFiles(srcdir * "/models")
  if !(file in model_dependency_orders)
    println(file)
    include(file)
  end
end

for file in allJuliaFiles(srcdir * "/controllers")
  println(file)
  include(file)
end
end


text:
"""
    bump_stdlibs(julia_repo; kwargs...)
"""
function bump_stdlibs(julia_repo::AbstractString; kwargs...)
    config = Config(; kwargs...)
    bump_stdlibs(julia_repo, config)
    return nothing
end

function bump_stdlibs(julia_repo::AbstractString, config::Config)
    repo_regex = r"^([\w\-\_]*?)\/([\w]*?)$"
    if !occursin(repo_regex, julia_repo)
        msg = "Invalid value for julia_repo: \"$(julia_repo)\". julia_repo must be of the form owner/repo. For example: julia_repo=\"JuliaLang/julia\""
        throw(ArgumentError(msg))
    end
    m = match(repo_regex, julia_repo)
    upstream_julia_repo_owner = m[1]
    upstream_julia_repo_name = m[2]
    upstream_julia_repo = "$(upstream_julia_repo_owner)/$(upstream_julia_repo_name)"
    @debug "Attempting to get upstream `julia` repo" upstream_julia_repo
    upstream_julia_repo_gh = GitHub.repo(upstream_julia_repo; config.auth)
    whoami = whoami_username(; config.auth)
    @debug("My GitHub username is: $(whoami)")
    fork_julia_repo_owner = whoami
    fork_julia_repo_name = upstream_julia_repo_name
    fork_julia_repo = "$(fork_julia_repo_owner)/$(fork_julia_repo_name)"
    fork_julia_repo_gh = create_or_get_fork(fork_julia_repo, upstream_julia_repo_gh; auth = config.auth)
    if config.julia_repo_default_branch isa Nothing
        upstream_julia_repo_default_branch = convert(String, upstream_julia_repo_gh.default_branch)::String
    else
        upstream_julia_repo_default_branch = convert(String, config.julia_repo_default_branch)::String
    end
    update_fork_branch(fork_julia_repo_gh, upstream_julia_repo_gh, upstream_julia_repo_gh.default_branch; auth = config.auth)
    stdlib_list = get_stdlib_list(upstream_julia_repo_gh, upstream_julia_repo_gh.default_branch; auth = config.auth)
    @info "Identified $(length(stdlib_list)) stdlibs that live in external repositories"
    for (i, stdlib) in enumerate(stdlib_list)
        @info "" i stdlib
    end
    filtered_stdlib_list = filter_stdlib_list(stdlib_list, config.stdlibs_to_include)
    @info "Based on the input criteria \"$(config.stdlibs_to_include)\", we will try to update $(length(filtered_stdlib_list)) stdlibs"
    for (i, stdlib) in enumerate(filtered_stdlib_list)
        @info "" i stdlib
    end
    errors = Any[]
    all_pr_branches = String[]
    state = State(;
        fork_julia_repo_gh = fork_julia_repo_gh,
        upstream_julia_repo_default_branch = upstream_julia_repo_default_branch,
        upstream_julia_repo_gh = upstream_julia_repo_gh,
    )
    for (i, stdlib) in enumerate(filtered_stdlib_list)
        @info "Starting to work on" i stdlib length(filtered_stdlib_list)
        try
            _bump_single_stdlib(stdlib, config, state)
            @info "Successfully finished working on" i stdlib length(filtered_stdlib_list)
        catch ex
            push!(errors, ex)
            @error "Encountered an error while working on" exception=(ex, catch_backtrace()) i stdlib length(filtered_stdlib_list)
        end
    end
    if config.close_old_pull_requests
        mktempdir() do temp_dir
            cd(temp_dir) do
                fork_clone_url = "https://x-access-token:$(config.auth.token)@github.com/$(state.fork_julia_repo_gh.full_name).git"
                run(`git clone $(fork_clone_url) FORK`)
                cd("FORK") do
                    run(`git fetch --all --prune`)
                    for (i, stdlib) in enumerate(filtered_stdlib_list)
                        prefix = "BumpStdlibs/$(stdlib.name)-"
                        delete_branches_with_prefix_on_origin_older_than(
                            prefix,
                            config.close_old_pull_requests_older_than;
                            exclude = state.all_pr_branches,
                        )
                    end
                end # cd("FORK")
            end # cd(temp_dir)
        end # end mktempdir()
    end
    if !isempty(errors)
        @error "Encountered at least one error" errors
        throw(ErrorException("Encountered at least one error"))
    end
    return nothing
end

function _bump_single_stdlib(stdlib::StdlibInfo, config::Config, state::State)
    auth = config.auth
    pr_branch_suffix = config.pr_branch_suffix
    fork_julia_repo_gh = state.fork_julia_repo_gh
    upstream_julia_repo_gh = state.upstream_julia_repo_gh
    upstream_julia_repo_default_branch = state.upstream_julia_repo_default_branch
    mktempdir() do temp_dir
        cd(temp_dir) do
            fork_clone_url = "https://x-access-token:$(auth.token)@github.com/$(fork_julia_repo_gh.full_name).git"
            run(`git clone $(fork_clone_url) FORK`)
            cd("FORK") do # we do this just to double-check that the `FORK` directory was created
            end
            run(`git clone $(stdlib.git_url) STDLIB`)
            cd("STDLIB") do
                run(`git checkout $(stdlib.branch)`)
                assert_current_branch_is(stdlib.branch)
                stdlib_latest_commit = strip(read(`git rev-parse HEAD`, String))
                stdlib_latest_commit_short = strip(read(`git rev-parse --short $(stdlib_latest_commit)`, String))
                assert_string_startswith(stdlib_latest_commit, stdlib_latest_commit_short)
                stdlib_current_commit_in_upstream = stdlib.current_shas["sha1"]
                stdlib_current_commit_in_upstream_short = strip(read(`git rev-parse --short $(stdlib_current_commit_in_upstream)`, String))
                assert_string_startswith(stdlib_current_commit_in_upstream, stdlib_current_commit_in_upstream_short)
                run(`git fetch --all --prune`)
                changelog_cmd = `git log --oneline $(stdlib_current_commit_in_upstream_short)..$(stdlib_latest_commit_short)`
                changelog = read(changelog_cmd, String)
                if stdlib_latest_commit == stdlib_current_commit_in_upstream
                    @info "stdlib is already up to date" stdlib stdlib_latest_commit stdlib_current_commit_in_upstream
                else
                    @info "stdlib is not up to date, so I will update it now" stdlib stdlib_latest_commit stdlib_current_commit_in_upstream
                    cd(joinpath(temp_dir, "FORK")) do
                        run(`git checkout $(upstream_julia_repo_default_branch)`)
                        assert_current_branch_is(upstream_julia_repo_default_branch)
                        pr_title_without_emoji = "Bump the $(stdlib.name) stdlib from $(stdlib_current_commit_in_upstream_short) to $(stdlib_latest_commit_short)"
                        pr_title = "🤖 $(pr_title_without_emoji)"
                        commit_message = pr_title
                        pr_branch_suffix_stripped = strip(pr_branch_suffix)
                        if isempty(pr_branch_suffix_stripped)
                            pr_branch_suffix_with_hyphen = ""
                        else
                            pr_branch_suffix_with_hyphen = "-$(pr_branch_suffix_stripped)"
                        end
                        pr_branch = "BumpStdlibs/$(stdlib.name)-$(stdlib_latest_commit_short)$(pr_branch_suffix_with_hyphen)"
                        push!(state.all_pr_branches, pr_branch)
                        git_url_markdown = _git_url_to_formatted_markdown(stdlib.git_url)
                        bumpstdlibs_sender = strip(get(ENV, "BUMPSTDLIBS_SENDER", ""))
                        bumpstdlibs_sender_ping = isempty(bumpstdlibs_sender) ? "unknown user" : "@$(bumpstdlibs_sender)"
                        pr_body_lines = String[
                            "Stdlib: $(stdlib.name)",
                            "URL: $(git_url_markdown)",
                            "Branch: $(stdlib.branch)",
                            "Old commit: $(stdlib_current_commit_in_upstream_short)",
                            "New commit: $(stdlib_latest_commit_short)",
                            "Bump invoked by: $(bumpstdlibs_sender_ping)",
                            "Powered by: [BumpStdlibs.jl](https://github.com/JuliaLang/BumpStdlibs.jl)",
                            "",
                            "```",
                            "\$ $(strip(string(changelog_cmd), '`'))",
                            changelog,
                            "```",
                        ]
                        pr_body = strip(join(strip.(pr_body_lines), "\n"))
                        run(`git checkout -B $(pr_branch)`)
                        assert_current_branch_is(pr_branch)
                        version_filename = joinpath("stdlib", "$(stdlib.name).version")
                        old_version_contents = strip(read(version_filename, String))
                        new_version_contents = replace(old_version_contents, stdlib_current_commit_in_upstream => stdlib_latest_commit)
                        rm(version_filename; force = true, recursive = true)
                        open(version_filename, "w") do io
                            println(io, new_version_contents)
                        end
                        assert_file_contains(version_filename, stdlib_latest_commit)
                        for (root, dirs, files) in walkdir(joinpath(pwd(), "deps", "checksums"))
                            for file in files
                                full_path = joinpath(root, file)
                                startswith(strip(file), "$(stdlib.name)-") && ispath(full_path) && rm(full_path; force = true, recursive = true)
                            end
                        end
                        directories_to_delete = String[]
                        for (root, dirs, files) in walkdir(joinpath(pwd(), "deps", "checksums"))
                            for dir in dirs
                                full_path = joinpath(root, dir)
                                if startswith(strip(dir), "$(stdlib.name)-")
                                    @debug "" root dir full_path
                                    if isdir(full_path)
                                        push!(directories_to_delete, full_path)
                                    end
                                end
                            end
                        end
                        for full_path in directories_to_delete
                            @debug "" full_path
                            if isdir(full_path)
                                rm(full_path; force = true, recursive = true)
                            end
                        end
                        checksum_files_that_need_refreshing = String[]
                        for (root, dirs, files) in walkdir(joinpath(pwd(), "deps", "checksums"))
                            for file in files
                                full_path = joinpath(root, file)
                                contains_matching_lines = delete_checksum_lines(;
                                    checksum_file = full_path,
                                    stdlib_name = stdlib.name,
                                )
                                if contains_matching_lines
                                    push!(checksum_files_that_need_refreshing, file)
                                end
                            end
                        end
                        cd("stdlib") do
                            run(`make`)
                        end
                        file_to_target = Dict(
                            "suitesparse" => "libsuitesparse",
                        )
                        for file in checksum_files_that_need_refreshing
                            if haskey(file_to_target, file)
                                target = file_to_target[file]
                            else
                                target = file
                            end
                            run(`make -f contrib/refresh_checksums.mk $(target)`)
                        end
                        run(`git add -A`)
                        run(`git commit -m "$(commit_message)"`)
                        do_push = true
                        if !(config.push_if_no_changes)
                            if pr_branch in get_origin_branches()
                                if git_diff_is_empty("HEAD", "origin/$(pr_branch)")
                                    do_push = false
                                end
                            end
                        end
                        @info "" do_push
                        if do_push
                            run(`git push --force origin $(pr_branch)`)
                        end
                        whoami = whoami_username(; config.auth)
                        pr_head_with_fork_owner = "$(whoami):$(pr_branch)"
                        pr_state = Dict{String, Any}(
                            "base" => upstream_julia_repo_default_branch,
                            "body" => strip(pr_body),
                            "head" => pr_head_with_fork_owner,
                            "maintainer_can_modify" => true,
                            "title" => pr_title,
                        )
                        @debug "" upstream_julia_repo_gh.full_name
                        @debug "" fork_julia_repo_gh.full_name
                        @debug "" pr_state
                        pr = create_or_update_pull_request(upstream_julia_repo_gh, pr_state; auth = config.auth)
                    end # cd(joinpath(temp_dir, "FORK"))
                end # if
            end # cd("STDLIB"
        end # cd(temp_dir)
    end # mktempdir()
    return nothing
end


text:
# TODO: replace this by a linear solver
import Tulip

import LinearAlgebra.LAPACK

"""
    UnitBlockAngularFactor{Tv}

"""
struct UnitBlockAngularFactor{Tv<:Real} <: Tulip.KKT.AbstractKKTSolver{Tv}

    A::UnitBlockAngularMatrix{Tv}

    θ::Vector{Tv}     # Diagonal scaling
    regP::Vector{Tv}  # Primal regularization
    regD::Vector{Tv}  # Dual   regularization

    D::Vector{Tv}  # Diagonal elements of Cholesky factor (n x 1)
    L::Matrix{Tv}  # Lower block of Cholesky factor (m0 x n)
    C::Matrix{Tv}  # Dense Schur complement (m0 x m0)

    # Local copies, needed to compute rank-k updates of the form C += M*Θ*M'
    _B0::Matrix{Tv}
    _B::Matrix{Tv}
    _L::Matrix{Tv}

    function UnitBlockAngularFactor(A::UnitBlockAngularMatrix{Tv}) where{Tv<:Real}
        return new{Tv}(
            A, Vector{Tv}(undef, A.N, ), Vector{Tv}(undef, A.N), Vector{Tv}(undef, A.M),
            Vector{Tv}(undef, A.R), Matrix{Tv}(undef, A.m0, A.R), Matrix{Tv}(undef, A.m0, A.m0),
            Matrix{Tv}(undef, A.m0, A.n0), Matrix{Tv}(undef, A.m0, A.n), Matrix{Tv}(undef, A.m0, A.R)
        )
    end
end

Tulip.KKT.backend(::UnitBlockAngularFactor) = "UnitBlockAngular"
Tulip.KKT.linear_system(::UnitBlockAngularFactor) = "Normal equations"

function Tulip.KKT.update!(
    kkt::UnitBlockAngularFactor{Tv},
    θ::AbstractVector{Tv},
    regP::AbstractVector{Tv},
    regD::AbstractVector{Tv}
) where{Tv<:Real}

    m0 = kkt.A.m0
    n0 = kkt.A.n0
    n  = kkt.A.n
    R  = kkt.A.R

    # Copy scaling
    kkt.θ .= θ
    kkt.regP .= regP
    kkt.regD .= regD

    # I. Compute normal equations
    # I.1 - Schur complement `Φ = Σ Bi * Θi * Bi'`
    θ_ = one(Tv) ./ (kkt.θ .+ kkt.regP)

    # Linking block
    @views θ0 = θ_[1:(kkt.A.n0)]
    copyto!(kkt._B0, kkt.A.B0)
    rmul!(kkt._B0, Diagonal(sqrt.(θ0)))
    BLAS.syrk!('U', 'N', one(Tv), kkt._B0, zero(Tv), kkt.C)  # C = B0 * Θ0 * B0'

    # Columns block
    @views θ1 = θ_[(1+kkt.A.n0):end]
    copyto!(kkt._B, kkt.A.B)
    rmul!(kkt._B , Diagonal(sqrt.(θ1)))
    BLAS.syrk!('U', 'N', one(Tv), kkt._B, one(Tv), kkt.C)  # C += Σ Bi * Θi * Bi'

    # I.2 Lower block
    kkt.L .= zero(Tv)
    @inbounds for j in 1:n
        z = θ1[j]
        k = kkt.A.blockidx[j]
        @inbounds for i in 1:m0
            kkt.L[i, k] += z * kkt.A.B[i, j]
        end
    end

    # I.2 Compute diagonal
    # Di = Σ { θj | block[j] == i }
    kkt.D .= 0
    @inbounds for (i, θj) in zip(kkt.A.blockidx, θ1)
        kkt.D[i] += θj
    end

    # I.3 Apply dual regularizations
    kkt.D .+= kkt.regD[(kkt.A.m0+1):end]

    # II. Compute the factorization
    kkt._L .= kkt.L
    rmul!(kkt._L, Diagonal(sqrt.(inv.(kkt.D))))
    BLAS.syrk!('U', 'N', -one(Tv), kkt._L, one(Tv), kkt.C)  # C -= (Di)⁻¹ (Biθi) (Biθi)'

    # 6. Factorize C
    # Dual regularizations on C
    @inbounds for i in 1:m0
        kkt.C[i, i] += kkt.regD[i]
    end
    _, ic = LAPACK.potrf!('U', kkt.C)
    ic == 0 || throw(PosDefException(ic))

    # Done
    return nothing
end

function Tulip.KKT.solve!(
    dx, dy,
    kkt::UnitBlockAngularFactor{Tv},
    ξp, ξd
) where{Tv<:Real}
    m0 = kkt.A.m0
    R = kkt.A.R
    # TODO: dimension checks

    # 1. Setup right-hand side for normal equations
    dy .= ξp
    mul!(dy, kkt.A, ξd ./ (kkt.θ .+ kkt.regP), one(Tv), one(Tv))
    @views dy0 = dy[1:m0]
    @views dy1 = dy[(m0+1):end]

    dy1 ./= kkt.D
    mul!(dy0, kkt.L, dy1, -one(Tv), one(Tv))

    # 2. Solve Schur complement
    LAPACK.potrs!('U', kkt.C, dy0)

    # 3. Recover solution to dy1
    dy1 .*= kkt.D
    mul!(dy1, kkt.L', dy0, -one(Tv), one(Tv))
    dy1 ./= kkt.D

    # 4. Recover solution to augmented system
    mul!(dx, kkt.A', dy)
    dx .-= ξd
    dx ./= kkt.θ .+ kkt.regP

    return nothing
end