text:
defaultOutputFormat(plt::Plot) = "png"

function png(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "png")
  io = open(fn, "w")
  writemime(io, MIME("image/png"), plt)
  close(io)
end
png(fn::AbstractString) = png(current(), fn)

function svg(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "svg")
  io = open(fn, "w")
  writemime(io, MIME("image/svg+xml"), plt)
  close(io)
end
svg(fn::AbstractString) = svg(current(), fn)

function pdf(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "pdf")
  io = open(fn, "w")
  writemime(io, MIME("application/pdf"), plt)
  close(io)
end
pdf(fn::AbstractString) = pdf(current(), fn)

function ps(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "ps")
  io = open(fn, "w")
  writemime(io, MIME("application/postscript"), plt)
  close(io)
end
ps(fn::AbstractString) = ps(current(), fn)

function tex(plt::Plot, fn::AbstractString)
  fn = addExtension(fn, "tex")
  io = open(fn, "w")
  writemime(io, MIME("application/x-tex"), plt)
  close(io)
end
tex(fn::AbstractString) = tex(current(), fn)

# ----------------------------------------------------------------

@compat const _savemap = Dict(
    "png" => png,
    "svg" => svg,
    "pdf" => pdf,
    "ps"  => ps,
    "tex" => tex,
  )

function getExtension(fn::AbstractString)
  pieces = split(fn, ".")
  length(pieces) > 1 || error("Can't extract file extension: ", fn)
  ext = pieces[end]
  haskey(_savemap, ext) || error("Invalid file extension: ", fn)
  ext
end

function addExtension(fn::AbstractString, ext::AbstractString)
  try
    oldext = getExtension(fn)
    if oldext == ext
      return fn
    else
      return "$fn.$ext"
    end
  catch
    return "$fn.$ext"
  end
end

function savefig(plt::Plot, fn::AbstractString)

  # get the extension
  local ext
  try
    ext = getExtension(fn)
  catch
    # if we couldn't extract the extension, add the default
    ext = defaultOutputFormat(plt)
    fn = addExtension(fn, ext)
  end

  # save it
  func = get(_savemap, ext) do
    error("Unsupported extension $ext with filename ", fn)
  end
  func(plt, fn)
end
savefig(fn::AbstractString) = savefig(current(), fn)

# ---------------------------------------------------------

gui(plt::Plot = current()) = display(PlotsDisplay(), plt)

function Base.display(::PlotsDisplay, plt::Plot)
    prepare_output(plt)
    _display(plt)
end

# override the REPL display to open a gui window
Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", plt::Plot) = gui(plt)

# ---------------------------------------------------------

const _mimeformats = Dict(
    "application/eps"         => "eps",
    "image/eps"               => "eps",
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
    "text/plain"              => "txt",
)

const _best_html_output_type = KW(
    :pyplot => :png,
    :unicodeplots => :txt,
)

# a backup for html... passes to svg or png depending on the html_output_format arg
function Base.writemime(io::IO, ::MIME"text/html", plt::Plot)
    output_type = Symbol(plt.attr[:html_output_format])
    if output_type == :auto
        output_type = get(_best_html_output_type, backend_name(plt.backend), :svg)
    end
    if output_type == :png
        # info("writing png to html output")
        print(io, "")
    elseif output_type == :svg
        # info("writing svg to html output")
        writemime(io, MIME("image/svg+xml"), plt)
    elseif output_type == :txt
        writemime(io, MIME("text/plain"), plt)
    else
        error("only png or svg allowed. got: $output_type")
    end
end

# for writing to io streams... first prepare, then callback
for mime in keys(_mimeformats)
    @eval function Base.writemime(io::IO, m::MIME{Symbol($mime)}, plt::Plot)
        prepare_output(plt)
        _writemime(io, m, plt)
    end
end

# ---------------------------------------------------------
# A backup, if no PNG generation is defined, is to try to make a PDF and use FileIO to convert

if is_installed("FileIO")
    @eval begin
        import FileIO
        function _writemime(io::IO, ::MIME"image/png", plt::Plot)
            fn = tempname()

            # first save a pdf file
            pdf(plt, fn)

            # load that pdf into a FileIO Stream
            s = FileIO.load(fn * ".pdf")

            # save a png
            pngfn = fn * ".png"
            FileIO.save(pngfn, s)

            # now write from the file
            write(io, readall(open(pngfn)))
        end
    end
end

# function html_output_format(fmt)
#     if fmt == "png"
#         @eval function Base.writemime(io::IO, ::MIME"text/html", plt::Plot)
#             print(io, "")
#         end
#     elseif fmt == "svg"
#         @eval function Base.writemime(io::IO, ::MIME"text/html", plt::Plot)
#             writemime(io, MIME("image/svg+xml"), plt)
#         end
#     else
#         error("only png or svg allowed. got: $fmt")
#     end
# end
#
# html_output_format("svg")

# ---------------------------------------------------------
# IJulia
# ---------------------------------------------------------

const _ijulia_output = Compat.ASCIIString["text/html"]

function setup_ijulia()
    # override IJulia inline display
    if isijulia()
        @eval begin
            import IJulia
            export set_ijulia_output
            function set_ijulia_output(mimestr::AbstractString)
                # info("Setting IJulia output format to $mimestr")
                global _ijulia_output
                _ijulia_output[1] = mimestr
            end
            function IJulia.display_dict(plt::Plot)
                global _ijulia_output
                Dict{Compat.ASCIIString, ByteString}(_ijulia_output[1] => sprint(writemime, _ijulia_output[1], plt))
            end
        end
        set_ijulia_output("text/html")
    end
end

# ---------------------------------------------------------
# Atom PlotPane
# ---------------------------------------------------------

function setup_atom()
    # @require Atom begin
    if isatom() && get(ENV, "PLOTS_USE_ATOM_PLOTPANE", false) in (true, 1, "1", "true", "yes")
        # @eval import Atom, Media
        @eval import Atom

        # connects the render function
        for T in (GadflyBackend,ImmerseBackend,PyPlotBackend,GRBackend)
            Atom.Media.media(Plot{T}, Atom.Media.Plot)
        end
        # Atom.Media.media{T <: Union{GadflyBackend,ImmerseBackend,PyPlotBackend,GRBackend}}(Plot{T}, Atom.Media.Plot)

        # Atom.displaysize(::Plot) = (535, 379)
        # Atom.displaytitle(plt::Plot) = "Plots.jl (backend: $(backend(plt)))"

        # this is like "display"... sends an html div with the plot to the PlotPane
        function Atom.Media.render(pane::Atom.PlotPane, plt::Plot)
            Atom.Media.render(pane, Atom.div(Atom.d(), Atom.HTML(stringmime(MIME("text/html"), plt))))
        end

        # function Atom.Media.render(pane::Atom.PlotPane, plt::Plot{PlotlyBackend})
        #     html = Media.render(pane, Atom.div(Atom.d(), Atom.HTML(stringmime(MIME("text/html"), plt))))
        # end
    end
end


text:
macro checked_lib(libname, path)
    const libptr = if VERSION >= v"0.4.0-dev+3844"
        Base.Libdl.dlopen_e(path)
    else
        Base.dlopen_e(path)
    end

    found = libptr != C_NULL
    if !found
        error("Unable to load $libname ($path)\n")
    end

    quote
        const $(esc(symbol(libname,"_found"))) = $found
        const $(esc(libname)) = $libptr
    end
end

macro load_symbol(symname, sym, libname)
    const symptr = if VERSION >= v"0.4.0-dev+3844"
        Base.Libdl.dlsym_e(eval(libname), eval(sym))
    else
        Base.dlsym_e(eval(libname), eval(sym))
    end

    found = symptr != C_NULL
    if !found
        error("Unable to load $sym ($libname)\n")
    end

    quote
        const $(esc(symbol(symname,"_found"))) = $found
        const $(esc(symname)) = $symptr
    end
end

@checked_lib _libinform "libinform"

ext: jl


text:
export JDFFile, @jdf_str, path

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you can apply `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
struct JDFFile{T<:AbstractString}
    path::T
end

"""
    jdf"path/to/JDFfile.jdf"

    JDFFile("path/to/JDFfile.jdf")

Define a JDF file, which you can apply `names` and `size`.

## Example
using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]
ncol(jdf"plsdel.jdf") # 2
size(jdf"plsdel.jdf") # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

size(jdf"plsdel.jdf", 1) # (2, 3)

# clean up
rm("plsdel.jdf", force = true, recursive = true)
"""
macro jdf_str(path)
    return :(JDFFile($path))
end

"""
    path(jdf::JDFFile)

Return the path of the JDF
"""
path(jdf) = getfield(jdf, :path)


text:
# https://github.com/JuliaLabs/brutus/blob/9d87d7aa543090ab6e8b98dc7413476c107a7075/Brutus/src/Brutus.jl#L98-L146

function get_methodinstance(@nospecialize(sig);
                            world=Base.get_world_counter(),
                            interp=Core.Compiler.NativeInterpreter(world))
    ms = Base._methods_by_ftype(sig, 1, Base.get_world_counter())
    @assert length(ms) == 1
    m = ms[1]
    mi = ccall(:jl_specializations_get_linfo,
               Ref{Core.MethodInstance}, (Any, Any, Any),
               m[3], m[1], m[2])
    return mi
end

function code_ircode_by_signature(@nospecialize(sig);
                                  world=Base.get_world_counter(),
                                  interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(ccall(:jl_specializations_get_linfo,
                              Ref{Core.MethodInstance},
                              (Any, Any, Any),
                              data[3], data[1], data[2]);
                        world=world, interp=interp)
            for data in Base._methods_by_ftype(sig, -1, world)]
end

function code_ircode(@nospecialize(f), @nospecialize(types=Tuple);
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(mi; world=world, interp=interp)
            for mi in Base.method_instances(f, types, world)]
end

function code_ircode(mi::Core.Compiler.MethodInstance;
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, interp)
    frame === nothing && return nothing
    if Core.Compiler.typeinf(interp, frame)
        opt_params = Core.Compiler.OptimizationParams(interp)
        opt = Core.Compiler.OptimizationState(frame, opt_params, interp)
        ir = Core.Compiler.run_passes(opt.src, opt.nargs - 1, opt)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return nothing
    # TODO(yhls): Fix this upstream
    resize!(ir.argtypes, opt.nargs)
    return ir => Core.Compiler.widenconst(result.result)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7927ff00-1806-11ec-2787-4db3c2bd8448
# imports
begin
	using Unitful, UnitfulAstro
	using UnitfulRecipes
	using Revise
	using LinearAlgebra
	using Pkg
	Pkg.develop(path="/Users/miles/dev/julia/Transits")
	using Transits
	using Plots
	using Markdown
end

# ╔═╡ 87591ebf-f36b-4643-bbf1-a961ed1825ff
md"""
# Equillibrium temperature in Sirius

When considering a planet orbiting Sirius B, it is interesting to consider the irradiation from Sirius A. This is especially significant due to Sirius B's eccentric orbit (``e=0.6``), leading to an overall range of separations from 7.9-32 AU.

The calculation for equilibrium temperature from a single star is
```math
T_\mathrm{eq} = \left[\frac{L(1 - A_B)}{16\sigma\pi a^2} \right]^{1/4}
```

we assume that adding an additional star will linearly add to the incident flux, so we can extend this to two stars with
```math
T_\mathrm{eq} = \sum_i{\left[\frac{L_i(1 - A_B)}{16\sigma\pi a_i^2} \right]^{1/4}}
```

This does not take into effect any reheating or geometric phase effects, which means this is likely an *upper limit*.
"""

# ╔═╡ b4375a6e-c669-4f64-aa9c-72375bf96868
begin
	plx = 0.3789 # arcsec
	dist = inv(plx)u"pc" # parsec
end

# ╔═╡ 0bdbeded-0d6f-426b-80bd-0d20d058baf8
# don't include i, Ω, or ω because we don't care about the
# line-of-sight rotation, only the relative position
# orbital elements from Bond et al. 2017
orbit = KeplerianOrbit(
	period=50.1284u"yr",
	a=7.4957u"AU" / plx, # convert arcsec to AU
	#incl=136.336u"°",
	#Omega=45.400u"°",
	t_0=1994.5715u"yr",
	ecc=0.59142,
	omega=0u"°"
	#omega=149.161u"°"
)

# ╔═╡ 0dd1df5c-7faa-4795-b4cd-0ebb9ba9e0c0
function relative_position(orbit, t)
	# orbit is in xzy because it wasn't rotated
	xzy = Transits.Orbits.relative_position(orbit, t)
	return [xzy[1], xzy[3]]
end

# ╔═╡ 181ea5bd-18f0-47ae-8911-9c823b8b5531
# plot orbit
begin
	t = range(0u"d", orbit.period, length=1000)
	t_today = 2020u"yr" + 260u"d" |> u"d"
	# calculate separations
	pos = mapreduce(ti -> relative_position(orbit, ti), hcat, t)
	plot(eachrow(pos)...,
		c=:black,
		leg=false,
		xlabel="Δx [AU]", ylabel="Δy [AU]")
	scatter!([0], [0], m=:star, ms=10, c=:black)
	pos_today = relative_position(orbit, t_today)
	scatter!(pos_today[1:1], pos_today[2:2], m=:star, ms=5, c=:black)
end

# ╔═╡ 1235c53f-e01f-457d-9057-7587af56a9f8
function eq_temp(orbit, t, dx, dy)
	pos = relative_position(orbit, t)u"AU"
	d = [dx, dy]
	a1 = norm(d)        # dist to Sirius A
	a2 = norm(d .- pos) # dist to Sirius B
	# assume albedo is negligible
	L1 = 25.4u"Lsun"
	L2 = 0.056u"Lsun"
	Ts = (L1 / (16π * Unitful.σ * a1^2))^(1/4) +
	     (L2 / (16π * Unitful.σ * a2^2))^(1/4)
	return Ts |> u"K"
end

# ╔═╡ 7738d09a-b7c5-4a38-b0fe-aafcba7e092a
begin
	dxs = range(-10, 35, length=500)u"AU"
	dys = range(-17, 17, length=500)u"AU"
end

# ╔═╡ 0dfa9c8b-0dd1-4967-99ce-6e729a3506a9
function eq_temp(orbit, t, dist)
	pos = relative_position(orbit, t)u"AU"
	d = norm(pos) # total orbital distance
	a1 = d - dist # distance to Sirius A
	a2 = dist     # distance to Sirius B
	# assume albedo is negligible
	L1 = 25.4u"Lsun"
	L2 = 0.056u"Lsun"
	T1 = (L1 / (16π * Unitful.σ * a1^2))^(1/4)
	T2 = (L2 / (16π * Unitful.σ * a2^2))^(1/4)
	return (T1 + T2, T1, T2) .|> u"K"
end

# ╔═╡ 57a6fbf8-0b5b-4a0e-b60a-e7d367692446
let temps
	temps = eq_temp.(orbit, t_today, dxs', dys)
	heatmap(dxs, dys, log10.(ustrip.(temps)),
		# levels=100,
		xlim=extrema(dxs),
		ylim=extrema(dys),
		c=:acton,
		cbar_title="log10(Teq)",
		title="today"
	)
	plot!(eachrow(pos)..., c=:white, leg=false, xlabel="Δx [AU]", ylabel="Δy [AU]")
	scatter!([0], [0], m=:star, ms=10, c=:white)
	scatter!(pos_today[1:1], m=:star, ms=5, pos_today[2:2], c=:white)
end

# ╔═╡ 56cd4737-3bdc-43da-a1fc-530baca135a3
let temps
	t_per = orbit.t_0 + orbit.t_ref
	temps = eq_temp.(orbit, t_per, dxs', dys)
	heatmap(dxs, dys, log10.(ustrip.(temps)),
		# levels=100,
		xlim=extrema(dxs),
		ylim=extrema(dys),
		c=:acton,
		cbar_title="log10(Teq)",
		title="periapsis"
	)
	plot!(eachrow(pos)..., c=:white, leg=false, xlabel="Δx [AU]", ylabel="Δy [AU]")
	scatter!([0], [0], m=:star, ms=10, c=:white)
	pos_per = relative_position(orbit, t_per)
	scatter!(pos_per[1:1], pos_per[2:2], m=:star, ms=5, c=:white)
end

# ╔═╡ 0523a4b3-7ecc-444e-853f-3adb46b8c04d
let temps
	t_apo = orbit.t_0 + orbit.t_ref + orbit.period/2
	temps = eq_temp.(orbit, t_apo, dxs', dys)
	heatmap(dxs, dys, log10.(ustrip.(temps)),
		# levels=100,
		xlim=extrema(dxs),
		ylim=extrema(dys),
		c=:acton,
		cbar_title="log10(Teq)",
		title="apoapsis"
	)
	plot!(eachrow(pos)..., c=:white, leg=false, xlabel="Δx [AU]", ylabel="Δy [AU]")
	scatter!([0], [0], m=:star, ms=10, c=:white)
	pos_apo = relative_position(orbit, t_apo)
	scatter!(pos_apo[1:1], pos_apo[2:2], m=:star, ms=5, c=:white)
end

# ╔═╡ cb972d70-1652-464a-9b35-215921674735
begin
	dists = range(0, 2, length=100)u"AU"
	plot(dists, map(d -> eq_temp(orbit, t_today, d)[1], dists),
		lab="",
		ylabel="Teq",
		xlabel="separation from Sirius B",
		title="today",
		xlim=extrema(dists)
	)
	plot!(dists, map(d -> eq_temp(orbit, t_today, d)[2], dists),
		ls=:dash,
		lab="Sirius A"
	)
	plot!(dists, map(d -> eq_temp(orbit, t_today, d)[3], dists),
		ls=:dash,
		lab="Sirius B"
	)
end

# ╔═╡ 6c63961a-f6b6-47e7-8dad-854ff1746aeb
begin
	t_per = orbit.t_0 + orbit.t_ref
	plot(dists, map(d -> eq_temp(orbit, t_per, d)[1], dists),
		lab="",
		ylabel="Teq",
		xlabel="separation from Sirius B",
		title="periapsis",
		xlim=extrema(dists)
	)
	plot!(dists, map(d -> eq_temp(orbit, t_per, d)[2], dists),
		ls=:dash,
		lab="Sirius A"
	)
	plot!(dists, map(d -> eq_temp(orbit, t_per, d)[3], dists),
		ls=:dash,
		lab="Sirius B"
	)
end

# ╔═╡ 860fc3e3-08d5-4cd5-8a4a-3b9db25a2b0e
begin
	t_apo = orbit.t_0 + orbit.t_ref + orbit.period/2
	plot(dists, map(d -> eq_temp(orbit, t_apo, d)[1], dists),
		lab="",
		ylabel="Teq",
		xlabel="separation from Sirius B",
		title="apoapsis",
		xlim=extrema(dists)
	)
	plot!(dists, map(d -> eq_temp(orbit, t_apo, d)[2], dists),
		ls=:dash,
		lab="Sirius A"
	)
	plot!(dists, map(d -> eq_temp(orbit, t_apo, d)[3], dists),
		ls=:dash,
		lab="Sirius B"
	)
end

# ╔═╡ d97e3642-f701-46e2-8c59-13148d2bebae
begin
	curves = hcat(
		map(d -> eq_temp(orbit, t_today, d)[1], dists),
		map(d -> eq_temp(orbit, t_per, d)[1], dists),
		map(d -> eq_temp(orbit, t_apo, d)[1], dists),
	)
	plot(dists, curves,
		lab=["today" "periapsis" "apoapsis"],
		ylabel="Teq",
		xlabel="separation from Sirius B",
		xlim=extrema(dists)
	)
end

# ╔═╡ Cell order:
# ╟─87591ebf-f36b-4643-bbf1-a961ed1825ff
# ╠═b4375a6e-c669-4f64-aa9c-72375bf96868
# ╠═0bdbeded-0d6f-426b-80bd-0d20d058baf8
# ╠═0dd1df5c-7faa-4795-b4cd-0ebb9ba9e0c0
# ╟─181ea5bd-18f0-47ae-8911-9c823b8b5531
# ╠═1235c53f-e01f-457d-9057-7587af56a9f8
# ╠═7738d09a-b7c5-4a38-b0fe-aafcba7e092a
# ╟─57a6fbf8-0b5b-4a0e-b60a-e7d367692446
# ╟─56cd4737-3bdc-43da-a1fc-530baca135a3
# ╟─0523a4b3-7ecc-444e-853f-3adb46b8c04d
# ╠═0dfa9c8b-0dd1-4967-99ce-6e729a3506a9
# ╟─cb972d70-1652-464a-9b35-215921674735
# ╟─6c63961a-f6b6-47e7-8dad-854ff1746aeb
# ╟─860fc3e3-08d5-4cd5-8a4a-3b9db25a2b0e
# ╟─d97e3642-f701-46e2-8c59-13148d2bebae
# ╠═7927ff00-1806-11ec-2787-4db3c2bd8448

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using DelimitedFiles

function solve(frequencyStart,input)
    changes = readdlm(input)
    # part 1
    finalFrequency = frequencyStart + sum(changes)
    println("final frequency: $finalFrequency")
    # part 2
    frequency = frequencyStart
    reached = Set()
    iteration = 0
    num_changes = length(changes)
    while ! in(frequency,reached)
        push!(reached,frequency)
        frequency+=changes[iteration%num_changes+1]
        iteration+=1
    end
    println("first duplicated frequency: $frequency at iteration $iteration")
end

solve(0,"input.txt")


text:
@testset "integration" begin
    prec = 64
    a = Acb(0, prec = prec)
    b = Acb(1, prec = prec)

    # Test with just a plain method
    f1 = sin
    f1! = Arblib.sin!
    res1 = "[0.459697694131860282 +/- 7.22e-19]"
    @test string(Arblib.integrate(f1, a, b, prec = prec)) == res1
    @test string(Arblib.integrate!(f1!, Acb(prec = prec), a, b)) == res1

    # Test with a method that accepts precision as a keyword argument
    f2 = (x; prec) -> Arblib.sin!(Acb(), x, prec = prec)
    f2! = Arblib.sin!
    res2 = "[0.459697694131860282 +/- 7.22e-19]"
    @test string(Arblib.integrate(f2, a, b, take_prec = true, prec = prec)) == res2
    @test string(Arblib.integrate!(f2!, Acb(prec = prec), a, b, take_prec = true)) == res2

    # Test with a method that accepts analytic as a keyword argument
    f3 = (x; analytic) -> Arblib.real_abs!(Acb(prec = prec), x, analytic)
    f3! = (res, x; analytic) -> Arblib.real_abs!(res, x, analytic)
    # FIXME: These are supposed to be identical but due to a bug in
    # Arb the produce slightly different results, see
    # https://github.com/kalmarek/Arblib.jl/issues/70. Once Arb is
    # updated so that they start to produce identical results this can
    # be updated and the issue closed.
    res3 = "[0.50000000000000000 +/- 2.73e-18]"
    @test string(Arblib.integrate(f3, a, b, check_analytic = true, prec = prec)) == res3
    @test string(Arblib.integrate!(f3!, Acb(prec = prec), a, b, check_analytic = true)) ==
          res3

    # Test with a method that accepts both precision and analytic as
    # a keyword arguments
    f4 = (x; analytic, prec) -> Arblib.real_abs!(Acb(), x, analytic, prec = prec)
    f4! = (res, x; analytic, prec) -> Arblib.real_abs!(res, x, analytic, prec = prec)
    # FIXME: See above
    res4 = "[0.50000000000000000 +/- 2.73e-18]"
    @test string(
        Arblib.integrate(f4, a, b, check_analytic = true, take_prec = true, prec = prec),
    ) == res4
    @test string(
        Arblib.integrate!(
            f4!,
            Acb(prec = prec),
            a,
            b,
            check_analytic = true,
            take_prec = true,
        ),
    ) == res4

    # Test with set tolerance
    f5 = x -> sin(exp(x))
    f5! = (res, x) -> Arblib.sin!(res, Arblib.exp!(res, x))
    res5 = "[0.8750 +/- 5.14e-5]"
    res5! = "[0.87495720 +/- 2.77e-9]"
    @test string(Arblib.integrate(f5, a, b, prec = prec, rtol = 1e-4)) == res5
    @test string(Arblib.integrate!(f5!, Acb(prec = prec), a, b, atol = 1e-8)) == res5!
end


text:
using Test
using Pluto.Configuration: CompilerOptions
using Pluto.WorkspaceManager: _merge_notebook_compiler_options
import Pluto: update_save_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell, project_relative_path

@testset "Workspace manager" begin
# basic functionality is already tested by the reactivity tests

    @testset "Multiple notebooks" begin

        fakeclientA = ClientSession(:fakeA, nothing)
        fakeclientB = ClientSession(:fakeB, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclientA.id] = fakeclientA
        🍭.connected_clients[fakeclientB.id] = fakeclientB

        notebookA = Notebook([
            Cell("x = 3")
        ])
        fakeclientA.connected_notebook = notebookA

        notebookB = Notebook([
            Cell("x")
        ])
        fakeclientB.connected_notebook = notebookB

        @test notebookA.path != notebookB.path

        Sys.iswindows() && sleep(.5) # workaround for https://github.com/JuliaLang/julia/issues/39270
        update_save_run!(🍭, notebookA, notebookA.cells[1])
        Sys.iswindows() && sleep(.5) # workaround for https://github.com/JuliaLang/julia/issues/39270
        update_save_run!(🍭, notebookB, notebookB.cells[1])

        @test notebookB.cells[1].errored == true

        Sys.iswindows() && sleep(.5) # workaround for https://github.com/JuliaLang/julia/issues/39270
        WorkspaceManager.unmake_workspace((🍭, notebookA))
        Sys.iswindows() && sleep(.5) # workaround for https://github.com/JuliaLang/julia/issues/39270
        WorkspaceManager.unmake_workspace((🍭, notebookB))
    end
    @testset "Variables with secret names" begin
        fakeclient = ClientSession(:fake, nothing)
        🍭 = ServerSession()
        🍭.options.evaluation.workspace_use_distributed = false
        🍭.connected_clients[fakeclient.id] = fakeclient

        notebook = Notebook([
            Cell("result = 1"),
            Cell("result"),
            Cell("elapsed_ns = 3"),
            Cell("elapsed_ns"),
        ])
        fakeclient.connected_notebook = notebook

        update_save_run!(🍭, notebook, notebook.cells[1:4])
        @test notebook.cells[1].output_repr == "1"
        @test notebook.cells[2].output_repr == "1"
        @test notebook.cells[3].output_repr == "3"
        @test notebook.cells[4].output_repr == "3"

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    @testset "notebook environment" begin
        session_options = CompilerOptions()
        notebook = Notebook([Cell("x")])
        notebook.compiler_options = CompilerOptions(;project="test")
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            joinpath(dirname(notebook.path), "test")

        notebook.compiler_options = CompilerOptions(;project=project_relative_path("test"))
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            project_relative_path("test")

        session_options = CompilerOptions(;project=project_relative_path("test"))
        notebook.compiler_options = CompilerOptions(;project=project_relative_path("Project.toml"))
        @test _merge_notebook_compiler_options(notebook, session_options).project ==
            project_relative_path("Project.toml")
    end
end


text:
@symbol_func function g_p(cur_reactor::AbstractReactor)
  cur_kappa = cur_reactor.kappa_95

  cur_delta = cur_reactor.delta_95

  cur_factor = 1.0

  cur_factor += 2 * cur_delta ^ 2

  cur_factor -= 1.2 * cur_delta ^ 3

  cur_factor *= cur_kappa ^ 2

  cur_factor += 1

  cur_factor /= 2

  cur_factor
end


text:
include("viaCheb.jl")
using .bycheb
using LinearAlgebra, Plots

method="Explicit"
#method="Implicit"
#method="Mixed"; μ = 0.5

α= 1.0
ν = 30.0
s = 30.0
Lrel = 10

#N=40
N = 71

M = 1

TINY = 0.01

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

function eye(n::Integer)
    diagm(ones(n))
end

nn = N÷2 + N%2
r = bycheb.ring(N,(α,ν,s,Lrel))
ws = sqrt(3*ν/s) # wave speed
v=[ones(nn);zeros(nn)]
#Σ = view(v,    1:nn )
#Φ = view(v, nn+1:e)
minAd1 = minimum( abs.(diff(r.y)./(abs.(r.u[1:end-1]).+ws)))
minAd2 = minimum( abs.(diff(r.y)./(abs.(r.u[2:end]).+ws)))
minDiff = minimum( diff(r.y).^2 / (3*ν) )
dt = TINY * 0.5 * min(minAd1, minAd2, minDiff)

Nm = trunc(Int,1.0/dt)
for i in 1:M*Nm
    global v
    local w
    if method=="Explicit"
        w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
        #Base.Threads.@threads
        for j in 1:2nn
            v[j] = w[j] + dt*sum((-r.AA[j,:]+r.B[j,:]).*w)
        end
    elseif method=="Implicit"
        w = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
        v = (eye(2*nn) - dt*(-r.AA+r.B))\w
    elseif method=="Mixed"
        w  = [abs.(v[1:nn]); -abs.(v[nn+1:end])]; w[1]=0.0; w[nn+1]=0.0; w./=sum(w[1:nn]);
        vv = w + μ*dt*(-r.AA+r.B)*w
        v = (eye(2*nn) - μ*dt*(-r.AA+r.B))\vv
    end
end
ynn=r.y[1:nn]
Σ = v[1:nn]
Φ = v[nn+1:end];
L0 = r.Ls["L0"]
Lcrit = r.Ls["Lcrit"]
Lvisc = r.Ls["Lvisc"]

Plots.plot(ynn,Σ)
Plots.vline!([L0-Lcrit],label="crit",legend=:bottomleft)
Plots.vline!([L0-Lvisc],label="visc")
Plots.vline!([L0],label="moon")

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
"""
    GeneralizedLinearMixedModel

Generalized linear mixed-effects model representation

# Fields
- `LMM`: a [`LinearMixedModel`](@ref) - the local approximation to the GLMM.
- `β`: the pivoted and possibly truncated fixed-effects vector
- `β₀`: similar to `β`. Used in the PIRLS algorithm if step-halving is needed.
- `θ`: covariance parameter vector
- `b`: similar to `u`, equivalent to `broadcast!(*, b, LMM.Λ, u)`
- `u`: a vector of matrices of random effects
- `u₀`: similar to `u`.  Used in the PIRLS algorithm if step-halving is needed.
- `resp`: a `GlmResp` object
- `η`: the linear predictor
- `wt`: vector of prior case weights, a value of `T[]` indicates equal weights.
The following fields are used in adaptive Gauss-Hermite quadrature, which applies
only to models with a single random-effects term, in which case their lengths are
the number of levels in the grouping factor for that term.  Otherwise they are
zero-length vectors.
- `devc`: vector of deviance components
- `devc0`: vector of deviance components at offset of zero
- `sd`: approximate standard deviation of the conditional density
- `mult`: multiplier

# Properties

In addition to the fieldnames, the following names are also accessible through the `.` extractor

- `theta`: synonym for `θ`
- `beta`: synonym for `β`
- `σ` or `sigma`: common scale parameter (value is `NaN` for distributions without a scale parameter)
- `lowerbd`: vector of lower bounds on the combined elements of `β` and `θ`
- `formula`, `trms`, `A`, `L`, and `optsum`: fields of the `LMM` field
- `X`: fixed-effects model matrix
- `y`: response vector

"""
struct GeneralizedLinearMixedModel{T <: AbstractFloat} <: MixedModel{T}
    LMM::LinearMixedModel{T}
    β::Vector{T}
    β₀::Vector{T}
    θ::Vector{T}
    b::Vector{Matrix{T}}
    u::Vector{Matrix{T}}
    u₀::Vector{Matrix{T}}
    resp::GLM.GlmResp
    η::Vector{T}
    wt::Vector{T}
    devc::Vector{T}
    devc0::Vector{T}
    sd::Vector{T}
    mult::Vector{T}
end

"""
    deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1)::T where {T}

Return the deviance of `m` evaluated by the Laplace approximation (`nAGQ=1`)
or `nAGQ`-point adaptive Gauss-Hermite quadrature.

If the distribution `D` does not have a scale parameter the Laplace approximation
is the squared length of the conditional modes, `u`, plus the determinant
of `Λ'Z'WZΛ + I`, plus the sum of the squared deviance residuals.
"""
function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1) where {T}
    nAGQ == 1 && return T(sum(m.resp.devresid) + logdet(m) + sum(u -> sum(abs2, u), m.u))
    u = vec(first(m.u))
    u₀ = vec(first(m.u₀))
    copyto!(u₀, u)
    ra = RaggedArray(m.resp.devresid, first(m.LMM.reterms).refs)
    devc0 = sum!(map!(abs2, m.devc0, u), ra)  # the deviance components at z = 0
    sd = map!(inv, m.sd, m.LMM.L[Block(1,1)].diag)
    mult = fill!(m.mult, 0)
    devc = m.devc
    for (z, w) in GHnorm(nAGQ)
        if !iszero(w)
            if iszero(z)  # devc == devc0 in this case
                mult .+= w
            else
                @. u = u₀ + z * sd
                updateη!(m)
                sum!(map!(abs2, devc, u), ra)
                @. mult += exp((abs2(z) + devc0 - devc)/2)*w
            end
        end
    end
    copyto!(u, u₀)
    updateη!(m)
    sum(devc0) - 2 * (sum(log, mult) + sum(log, sd))
end

StatsBase.deviance(m::GeneralizedLinearMixedModel) = deviance(m, m.optsum.nAGQ)

"""
deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)

Update `m.η`, `m.μ`, etc., install the working response and working weights in
`m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate the [`deviance`](@ref).
"""
function deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
    updateη!(m)
    GLM.wrkresp!(m.LMM.y, m.resp)
    reweight!(m.LMM, m.resp.wrkwt)
    deviance(m, nAGQ)
end

GLM.dispersion(m::GeneralizedLinearMixedModel, sqr::Bool=false) =
    dispersion(m.resp, dof_residual(m), sqr)

GLM.dispersion_parameter(m::GeneralizedLinearMixedModel) = dispersion_parameter(m.resp.d)

function StatsBase.dof(m::GeneralizedLinearMixedModel)
    length(m.β) + length(m.θ) + GLM.dispersion_parameter(m.resp.d)
end

fit(::Type{GeneralizedLinearMixedModel}, f::FormulaTerm, tbl,
    d::Distribution = Normal(), l::Link = canonicallink(d);
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    offset = [],
    verbose::Bool = false,
    fast::Bool = false,
    nAGQ::Integer = 1) =
    fit(GeneralizedLinearMixedModel, f, columntable(tbl), d, l,
        wts = wts,
        offset = offset,
        contrasts = contrasts,
        verbose = verbose,
        fast = fast,
        nAGQ = nAGQ)
fit(::Type{GeneralizedLinearMixedModel}, f::FormulaTerm, tbl::Tables.ColumnTable,
    d::Distribution, l::Link = canonicallink(d);
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    offset = [],
    verbose::Bool = false,
    fast::Bool = false,
    nAGQ::Integer = 1) =
    fit!(GeneralizedLinearMixedModel(f, tbl, d, l,
                                     wts = wts,
                                     offset = offset,
                                     contrasts = contrasts),
                                     verbose = verbose,
                                     fast = fast,
                                     nAGQ = nAGQ)
"""
    fit!(m::GeneralizedLinearMixedModel[, verbose = false, fast = false, nAGQ = 1])

Optimize the objective function for `m`.

When `fast` is `true` a potentially much faster but slightly less accurate algorithm, in
which `pirls!` optimizes both the random effects and the fixed-effects parameters,
is used.
"""
function fit!(m::GeneralizedLinearMixedModel{T};
              verbose::Bool = false,
              fast::Bool = false,
              nAGQ::Integer = 1) where {T}
    β = m.β
    lm = m.LMM
    optsum = lm.optsum
    if !fast
        fit!(m, verbose=verbose, fast=true, nAGQ=nAGQ)
        optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
        optsum.initial = vcat(β, m.θ)
        optsum.final = copy(optsum.initial)
        optsum.initial_step = vcat(stderror(m) ./ 3, min.(T(0.05), m.θ ./ 4))
    end
    setpar! = fast ? setθ! : setβθ!
    feval = 0
    function obj(x, g)
        isempty(g) || error("gradient not defined for this model")
        feval += 1
        val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
        feval == 1 && (optsum.finitial = val)
        if verbose
            println("f_", feval, ": ", val, " ", x)
        end
        val
    end
    opt = Opt(optsum)
    NLopt.min_objective!(opt, obj)
    fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial))
    ## check if very small parameter values bounded below by zero can be set to zero
    xmin_ = copy(xmin)
    for i in eachindex(xmin_)
        if iszero(optsum.lowerbd[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    if xmin ≠ xmin_
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    pirls!(setpar!(m, xmin), fast, verbose)
    optsum.nAGQ = nAGQ
    optsum.feval = feval
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXEVAL_REACHED]
        @warn("NLopt optimization failure: $ret")
    end
    m
end

StatsBase.fitted(m::GeneralizedLinearMixedModel) = m.resp.mu

function fixef(m::GeneralizedLinearMixedModel{T}, permuted=true) where {T}
    permuted && return m.β
    Xtrm = first(m.LMM.feterms)
    piv = Xtrm.piv
    v = fill(-zero(T), size(piv))
    copyto!(view(v, 1:Xtrm.rank), m.β)
    invpermute!(v, piv)
end
GeneralizedLinearMixedModel(f::FormulaTerm, tbl,
                            d::Distribution,
                            l::Link = canonicallink(d);
                            wts = [],
                            offset = [],
                            contrasts = Dict{Symbol,Any}()) =
    GeneralizedLinearMixedModel(f, Tables.columntable(tbl),
                                d, l;
                                wts = [],
                                offset = [],
                                contrasts = Dict{Symbol,Any}())
GeneralizedLinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable,
                            d::Normal,
                            l::IdentityLink;
                            wts = [],
                            offset = [],
                            contrasts = Dict{Symbol,Any}()) =
    throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
function GeneralizedLinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable,
                                     d::Distribution,
                                     l::Link = canonicallink(d);
                                     wts = [],
                                     offset = [],
                                     contrasts = Dict{Symbol,Any}())
    if isa(d, Binomial) && isempty(wts)
        d = Bernoulli()
    end
    (isa(d, Normal) && isa(link, IdentityLink)) &&
        throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
    LMM = LinearMixedModel(f, tbl, contrasts = contrasts; wts = wts)
    y = copy(LMM.y)
        # the sqrtwts field must be the correct length and type but we don't know those
        # until after the model is constructed if wt is empty.  Because a LinearMixedModel
        # type is immutable, another one must be created.
    if isempty(wts)
        LMM = LinearMixedModel(LMM.formula, LMM.reterms, LMM.feterms, fill!(similar(y), 1),
            LMM.A, LMM.L, LMM.optsum)
    end
    updateL!(LMM)
        # fit a glm to the fixed-effects only - awkward syntax is to by-pass a test
    gl = isempty(wts) ? glm(LMM.X, y, d, l) : glm(LMM.X, y, d, l, wts = wts)
    β = coef(gl)
    u = [fill(zero(eltype(y)), vsize(t), nlevs(t)) for t in LMM.reterms]
        # vv is a template vector used to initialize fields for AGQ
        # it is empty unless there is a single random-effects term
    vv = length(u) == 1 ? vec(first(u)) : similar(y, 0)

    res = GeneralizedLinearMixedModel(LMM, β, copy(β), LMM.θ, copy.(u), u,
        zero.(u), gl.rr, similar(y), oftype(y, wts), similar(vv),
        similar(vv), similar(vv), similar(vv))
    deviance!(res, 1)
    res
end

function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
    if s == :theta
        m.θ
    elseif s == :coef
        coef(m)
    elseif s == :beta
        m.β
    elseif s ∈ (:σ, :sigma)
        sdest(m)
    elseif s ∈ (:A, :L, :λ, :lowerbd, :optsum, :X, :reterms, :feterms, :formula, :σs, :σρs)
        getproperty(m.LMM, s)
    elseif s == :y
        m.resp.y
    else
        getfield(m, s)
    end
end

function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
    accum = zero(T)
    D = Distribution(m.resp)
    if D <: Binomial
        for (μ, y, n) in zip(m.resp.mu, m.resp.y, m.wt)
            accum += logpdf(D(round(Int, n), μ), round(Int, y * n))
        end
    else
        for (μ, y) in zip(m.resp.mu, m.resp.y)
            accum += logpdf(D(μ), y)
        end
    end
    accum - (mapreduce(u -> sum(abs2, u), + , m.u) + logdet(m)) / 2
end

StatsBase.nobs(m::GeneralizedLinearMixedModel) = length(m.η)

StatsBase.predict(m::GeneralizedLinearMixedModel) = fitted(m)

Base.propertynames(m::GeneralizedLinearMixedModel, private=false) =
    (:A, :L, :theta, :beta, :coef, :λ, :lambda, :σ, :sigma, :X, :y, :lowerbd, :σρs, :σs,
    fieldnames(typeof(m))...)

"""
    pirls!(m::GeneralizedLinearMixedModel)

Use Penalized Iteratively Reweighted Least Squares (PIRLS) to determine the conditional
modes of the random effects.

When `varyβ` is true both `u` and `β` are optimized with PIRLS.  Otherwise only `u` is
optimized and `β` is held fixed.

Passing `verbose = true` provides verbose output of the iterations.
"""
function pirls!(m::GeneralizedLinearMixedModel{T}, varyβ=false, verbose=false; maxiter::Integer=10) where {T}
    u₀ = m.u₀
    u = m.u
    β = m.β
    β₀ = m.β₀
    lm = m.LMM
    for j in eachindex(u)         # start from u all zeros
        copyto!(u₀[j], fill!(u[j], 0))
    end
    varyβ && copyto!(β₀, β)
    obj₀ = deviance!(m) * 1.0001
    if verbose
        print("varyβ = ", varyβ, ", obj₀ = ", obj₀)
        if varyβ
            print(", β = ")
            show(β)
        end
        println()
    end
    for iter in 1:maxiter
        varyβ && ldiv!(adjoint(feL(m)), copyto!(β, lm.L.blocks[end, end - 1]))
        ranef!(u, m.LMM, β, true) # solve for new values of u
        obj = deviance!(m)        # update GLM vecs and evaluate Laplace approx
        verbose && println(lpad(iter,4), ": ", obj)
        nhalf = 0
        while obj > obj₀
            nhalf += 1
            if nhalf > 10
                if iter < 2
                    throw(ErrorException("number of averaging steps > 10"))
                end
                break
            end
            for i in eachindex(u)
                map!(average, u[i], u[i], u₀[i])
            end
            varyβ && map!(average, β, β, β₀)
            obj = deviance!(m)
            verbose && println(lpad(nhalf, 8), ", ", obj)
        end
        if isapprox(obj, obj₀; atol = 0.00001)
            break
        end
        copyto!.(u₀, u)
        copyto!(β₀, β)
        obj₀ = obj
    end
    m
end

"""
    setβθ!(m::GeneralizedLinearMixedModel, v)

Set the parameter vector, `:βθ`, of `m` to `v`.

`βθ` is the concatenation of the fixed-effects, `β`, and the covariance parameter, `θ`.
"""
function setβθ!(m::GeneralizedLinearMixedModel, v)
    setβ!(m, v)
    setθ!(m, view(v, (length(m.β) + 1) : length(v)))
end

function setβ!(m::GeneralizedLinearMixedModel, v)
    β = m.β
    copyto!(β, view(v, 1:length(β)))
    m
end

function setθ!(m::GeneralizedLinearMixedModel, v)
    setθ!(m.LMM, copyto!(m.θ, v))
    m
end

function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y)
    if s == :β
        setβ!(m, y)
    elseif s == :θ
        setθ!(m, y)
    elseif s == :βθ
        setβθ!(m, y)
    else
        setfield!(m, s, y)
    end
end

sdest(m::GeneralizedLinearMixedModel{T}) where {T} = convert(T, NaN)

function Base.show(io::IO, m::GeneralizedLinearMixedModel)
    nAGQ = m.LMM.optsum.nAGQ
    println(io, "Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = $nAGQ)")
    println(io, "  ", m.LMM.formula)
    println(io, "  Distribution: ", Distribution(m.resp))
    println(io, "  Link: ", GLM.Link(m.resp), "\n")
    println(io, string("  Deviance: ", @sprintf("%.4f", deviance(m, nAGQ))), "\n")

    show(io,VarCorr(m))

    print(io," Number of obs: $(length(m.y)); levels of grouping factors: ")
    join(io, nlevs.(m.reterms), ", ")
    println(io)

    println(io, "\nFixed-effects parameters:")
    show(io, coeftable(m))
end

"""
    updateη!(m::GeneralizedLinearMixedModel)

Update the linear predictor, `m.η`, from the offset and the `B`-scale random effects.
"""
function updateη!(m::GeneralizedLinearMixedModel)
    η = m.η
    b = m.b
    u = m.u
    reterms = m.LMM.reterms
    mul!(η, modelmatrix(m), m.β)
    for i in eachindex(b)
        unscaledre!(η, reterms[i], mul!(b[i], reterms[i].λ, u[i]))
    end
    GLM.updateμ!(m.resp, η)
    m
end

varest(m::GeneralizedLinearMixedModel{T}) where {T} = one(T)

            # delegate GLMM method to LMM field
for f in (
    :describeblocks,
    :feL,
    :(LinearAlgebra.logdet),
    :lowerbd,
    :(StatsModels.modelmatrix),
    :(StatsBase.vcov),
    :σs,
    :σρs,
    )
    @eval begin
        $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
    end
end


text:
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
using JLD2
using NNFEM
using ADCME

global nntype="ae_scaled"
stress_scale = 100.0

include("../nnutil.jl")

for tid in [1,3]
    @load "../Data/domain$tid.jld2" domain
    @load "../Data/domain$(tid)_te.jld2" domain_te

    close("all")
    u1 = hcat(domain.history["state"]...)
    u2 = hcat(domain_te.history["state"]...)
    err = zeros(8, 101)
    for i = 1:8
        u1_ = [u1[i,:] u1[i+8,:]]
        u2_ = [u2[i,:] u2[i+8,:]]
        err[i,:] = sqrt.(sum((u1_-u2_).^2,dims=2)[:])
    end

    T = 0.5
    NT = 100
    t = LinRange(0.0,T, NT+1)
    for i = 1:8
        plot(t, err[i,:], label="$i")
    end
    legend()
    xlabel("\$t\$")
    ylabel("\$u_{ref}-u_{est}\$")
    mpl.save("truss2d_loc_diff$tid.tex")
    savefig("truss2d_loc_diff$tid.pdf")

    sess = Session(); init(sess)

    close("all")
    X, Y = prepare_strain_stress_data1D(domain)
    x = constant(X)
    y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
    sess = Session(); init(sess)
    ADCME.load(sess, "../Data/trained_nn_fem.mat")
    out = run(sess, y)
    scatter(X[:,1], Y[:], marker=".",s=5, label="Reference")
    scatter(X[:,1], out[:], marker=".",s=5, label="Estimated")
    xlabel("Strain")
    ylabel("Stress")
    legend()
    grid("on")
    mpl.save("truss2d_stress_diff$tid.tex")
    savefig("truss2d_stress_diff$tid.pdf")

end


text:
# comparing with ROCR
data = CSV.read(joinpath(@__DIR__, "data", "data.csv"))
# comparing with sklearn
roc_sklearn = CSV.read(joinpath(@__DIR__, "data", "roc_sklearn.csv"))

scores = data.predictions
labels = data.labels

fpr_sklearn = roc_sklearn.fpr
tpr_sklearn = roc_sklearn.tpr
op_sklearn = roc_sklearn.op

# testing constructors
curve = ROC(labels, scores)
B = BinaryMetrics(labels, scores)
curve = ROC(B)
println(curve)
@test isapprox( auc(curve),       0.834187    , atol = 1e-4) # AUCROCR
@test isapprox( auc(curve),       0.8341875188423276 ) # sklearn AUCROC
@test norm(fpr_sklearn[2:end]-fpr(curve)) < 1e-5
@test norm(tpr_sklearn[2:end]-tpr(curve)) < 1e-5
@test norm(op_sklearn[2:end]-op(curve)) < 1e-5

OP = op(B)
@test op(curve) == OP
@test auc(curve) == curve.AUC
@test fpr(curve) == fpr(B,OP)
@test tpr(curve) == tpr(B,OP)

# "ROC analysis: web-based calculator for ROC curves' example"
scores = [1 , 2 , 3 , 4 , 6 , 5 , 7 , 8 , 9 , 10]
labels = [0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 1]
res = ROC(labels, scores)

@test auc(res) ≈ 0.96
@test tpr(res) == [0.2, 0.4, 0.6, 0.8, 0.8, 1.0, 1.0, 1.0, 1.0, 1.0]
@test fpr(res) == [0.0, 0.0, 0.0, 0.0, 0.2, 0.2, 0.4, 0.6, 0.8, 1.0]

# Example from https://www.epeter-stats.de/roc-curves-and-ties/
# This package uses the first strategy, i.e.  AUC is equivalent to the
# Mann-Whitney U statistic

labels = [ones(Bool,1000); zeros(Bool,1000)]
scores =  [rand(7:14, 1000); rand(1:8, 1000)]
roc_data = ROC(labels, scores)
@test abs( round(auc(roc_data), digits=2) - 0.97 ) ≤ 0.015

# Example from https://github.com/brian-lau/MatlabAUC/issues/1
data = [-1 1; 1 2; -1 3; -1 4; 1 5; -1 6; 1 7; -1 8; 1 9; 1 10;
        1 11; -1 13; 1 13; 1 14; 1 14]
labels = data[:,1] .> 0
scores = data[:,2]
roc_data = ROC(labels, scores)
spss = [1.000 1.000
        1.000 .833
        .889 .833
        .889 .667
        .889 .500
        .778 .500
        .778 .333
        .667 .333
        .667 .167
        .556 .167
        .444 .167
        .333 .167
        .222 .000
        .000 .000]

@test norm(round.(reverse(spss[:,1])[2:end] - tpr(roc_data), digits=2)) == 0
@test norm(round.(reverse(spss[:,2])[2:end] - fpr(roc_data), digits=2)) == 0

# Are AUC and ROC consistent after permutation of scores and labels?
perm = randperm(100)
scores = rand(100)
labels = [zeros(50); ones(50)]

@test auc(ROC(labels, scores)) == auc(ROC(labels[perm], scores[perm]))


text:
using UUIDs

using TeaFiles.Header: AbstractSection, ContentDescriptionSection, Field,
    ItemSection, NameValue, NameValueSection, TeaFileMetadata, TimeSection,
    MAGIC_VALUE, _tea_size, _write_section, _write_tea, section_id

to_tea_byte_array(x::Real) = reinterpret(UInt8, [x])

function to_tea_byte_array(x)::Vector{UInt8}
    out = IOBuffer()
    _write_tea(out, x)
    return out.data[1:out.size]
end

function _test_name_value(
    name::AbstractString,
    value::Union{Int32,Float64,AbstractString,Base.UUID}
)
    name_value = NameValue(name, value)
    out = IOBuffer()
    written_bytes = _write_tea(out, name_value)
    @test written_bytes == out.size
    @test written_bytes == ((4 + 3) + 4 + _tea_size(value))

    kind = if isa(value, Int32)
        Int32(1)
    elseif isa(value, Float64)
        Int32(2)
    elseif isa(value, AbstractString)
        Int32(3)
    elseif isa(value, Base.UUID)
        Int32(4)
    end

    @test out.data[1:written_bytes] == vcat(
        to_tea_byte_array(name),
        to_tea_byte_array(kind),
        to_tea_byte_array(value)
    )
end

function _test_section(section::AbstractSection)
    out = IOBuffer()
    written_bytes = _write_section(out, section)
    @test written_bytes == out.size
    @test written_bytes == 4 + 4 + _tea_size(section)

    @test out.data[1:written_bytes] == vcat(
        to_tea_byte_array(section_id(section)),
        to_tea_byte_array(_tea_size(section)),
        to_tea_byte_array(section)
    )
end

function _test_metadata(metadata::TeaFileMetadata)
    out = IOBuffer()
    written_bytes = write(out, metadata)
    @test written_bytes == out.size
    @test written_bytes == metadata.item_start

    # Work out the byte array for all sections present.
    sections_bytes = vcat(
        map(
            section -> vcat(
                to_tea_byte_array(section_id(section)),
                to_tea_byte_array(_tea_size(section)),
                to_tea_byte_array(section)
            ),
            metadata.sections
        )...
    )

    # Work out the amount of padding that should have been added.
    unpadded_length = (
        8 + 8 + 8 + 8 +
        if isempty(metadata.sections)
            0
        else
            sum(
                section -> (4 + 4 + _tea_size(section)),
                metadata.sections
            )
        end
    )
    padding_bytes = repeat([0x00], metadata.item_start - unpadded_length)

    @test out.data[1:written_bytes] == vcat(
        to_tea_byte_array(MAGIC_VALUE),
        to_tea_byte_array(metadata.item_start),
        to_tea_byte_array(metadata.item_end),
        to_tea_byte_array(length(metadata.sections)),
        sections_bytes,
        padding_bytes
    )
end

@testset "write" begin
    @testset "write_field" begin
        field = Field(1, 2, "hi!")
        out = IOBuffer()
        written_bytes = _write_tea(out, field)
        @test written_bytes == out.size
        @test written_bytes == (4 + 4 + (4 + 3))

        @test out.data[1:written_bytes] == UInt8[
            1, 0, 0, 0,  # 1::Int32
            2, 0, 0, 0,  # 2::Int32
            3, 0, 0, 0,  # 3::Int32
            0x68,  # "h"
            0x69,  # "i"
            0x21,  # "!"
            ]
    end

    @testset "write_name_value" begin
        @testset "int32" begin _test_name_value("moo", Int32(32)) end
        @testset "float64" begin _test_name_value("moo", 42.0) end
        @testset "string" begin _test_name_value("moo", "lowing") end
        @testset "uuid" begin _test_name_value("moo", UUIDs.uuid1()) end
    end

    @testset "write_item_section" begin
        _test_section(ItemSection("", Field[]))
        _test_section(ItemSection("Plop", _make_fields([Float64, Float64, Int32])))
        _test_section(ItemSection(
            "Mooo",
            _make_fields(["a" => Float64, "b" => Float64, "c" => Int32])
        ))
    end

    @testset "write_time_section" begin
        _test_section(TimeSection(0, 0, Int32[]))
        _test_section(TimeSection(
            719162,  # 1970-01-01
            86400,  # Seconds
            Int32[0]  # First field
        ))
    end

    @testset "write_content_description_section" begin
        _test_section(ContentDescriptionSection(""))
        _test_section(ContentDescriptionSection("Some non-empty description"))
    end

    @testset "write_name_value_section" begin
        _test_section(NameValueSection([]))
        _test_section(NameValueSection([NameValue("a", Int32(1))]))
        _test_section(NameValueSection([
            NameValue("a", Int32(1)),
            NameValue("b", Float64(42.0)),
            NameValue("c", "coo"),
            NameValue("d", UUIDs.uuid1()),
        ]))
    end

    @testset "write_tea_file_metadata" begin
        _test_metadata(TeaFileMetadata(200, 0, []))

        # This is the example in the specification.
        example_metadata = _get_example_metadata()
        _test_metadata(example_metadata)

        # We should pad up to the item_start
        @test write(IOBuffer(), example_metadata) == 200
        @test write(IOBuffer(), TeaFileMetadata(1000, 0, [])) == 1000
    end
end


text:
if normpath(@__DIR__) ∉ LOAD_PATH
    pushfirst!(LOAD_PATH, normpath(@__DIR__, "../.."))
    pushfirst!(LOAD_PATH, normpath(@__DIR__))
end

using DECAES
using BenchmarkTools
using PrettyTables
using LibGit2

function commit_hash()
    repo = LibGit2.GitRepo(normpath(@__DIR__, "../.."))
    hash = LibGit2.GitShortHash(LibGit2.peel(LibGit2.GitCommit, LibGit2.head(repo)))
    return string(hash)
end

function build_suite()
    # Build suite of EPG algorithm benchmark groups
    suite = BenchmarkGroup()
    for ETL in [8,16,32,48,64]
        suite_ETL = suite["ETL=$ETL"] = BenchmarkGroup()
        for T in [Float32, Float64]
            suite_T = suite_ETL["T=$T"] = BenchmarkGroup()
            for alg in DECAES.EPG_Algorithms
                alg == DECAES.EPGWork_ReIm_Generated && continue
                α, TE, T2, T1, β = T(163.0), T(11e-3), T(39e-3), T(1.1), T(151.0)
                θ    = DECAES.EPGOptions((; α, TE, T2, T1, β), Val(ETL), T)
                work = alg(T, ETL)
                suite_T["alg=$(nameof(alg))"] = @benchmarkable DECAES.EPGdecaycurve!($work, $θ)
            end
        end
    end
    return suite
end

# Display results
function default_pretty_table(io, data, header, row_names; backend = :text, kwargs...)
    is_minimum = (data,i,j) -> !isnan(data[i,j]) && data[i,j] ≈ minimum(filter(!isnan, data[i,:]))
    hl = PrettyTables.Highlighter(is_minimum, foreground = :blue, bold = true)
    PrettyTables.pretty_table(io, data, header; backend = :text, row_names, highlighters = (hl,), formatters = (v,i,j) -> round(v, sigdigits = 3), kwargs...)
end

function print_results(io, results)
    names, times = map(leaves(results)) do (name, results)
        name, time(minimum(results)) / 1000
    end |> xs -> ((x->x[1]).(xs), (x->x[2]).(xs))
    alg_names = (x->x[end]).(names) |> unique |> permutedims
    row_names = reshape(names, length(alg_names), :)[1,:] .|> x -> join(x[1:end-1], ", ")
    tbl_data  = reshape(times, length(alg_names), :) |> permutedims
    tbl_hdr   = vcat(alg_names, fill("Time [us]", size(alg_names)...))
    default_pretty_table(io, tbl_data, tbl_hdr, row_names; backend = :text)
end

function main()
    @info "building suite"
    suite = build_suite()

    @info "tuning suite"
    tune!(suite)

    @info "running benchmarks"
    results = run(suite; verbose = true, seconds = 1)

    @info "saving results"
    print_results(stdout, results)
    open(joinpath(@__DIR__, "epg_bench_$(commit_hash()).txt"); write = true) do io
        print_results(io, results)
    end
    return results
end

# v1.5.3
# ┌───────────────────┬──────────────────┬───────────────────────────────┬────────────────────────────┬──────────────────────────────────┬────────────────────────┬─────────────────┬──────────────────┬─────────────────────────────────────┐
# │                   │ alg=EPGWork_ReIm │ alg=EPGWork_Cplx_Vec_Unrolled │ alg=EPGWork_ReIm_DualCache │ alg=EPGWork_ReIm_DualCache_Split │ alg=EPGWork_Basic_Cplx │ alg=EPGWork_Vec │ alg=EPGWork_Cplx │ alg=EPGWork_ReIm_DualCache_Unrolled │
# │                   │        Time [us] │                     Time [us] │                  Time [us] │                        Time [us] │              Time [us] │       Time [us] │        Time [us] │                           Time [us] │
# ├───────────────────┼──────────────────┼───────────────────────────────┼────────────────────────────┼──────────────────────────────────┼────────────────────────┼─────────────────┼──────────────────┼─────────────────────────────────────┤
# │ ETL=48, T=Float64 │            0.833 │                          2.56 │                      0.844 │                            0.831 │                   6.28 │            1.97 │             3.86 │                               0.962 │
# │ ETL=48, T=Float32 │            0.822 │                          3.13 │                      0.825 │                            0.816 │                   6.02 │            1.85 │              3.7 │                               0.954 │
# │ ETL=64, T=Float64 │             1.45 │                          4.31 │                       1.54 │                             1.44 │                   10.5 │            3.33 │             6.62 │                                 1.6 │
# │ ETL=64, T=Float32 │             1.44 │                          4.95 │                       1.44 │                             1.46 │                   10.1 │             3.2 │             6.43 │                                1.61 │
# │ ETL=32, T=Float64 │            0.395 │                          1.34 │                      0.396 │                             0.39 │                   3.13 │           0.957 │             1.93 │                               0.487 │
# │ ETL=32, T=Float32 │            0.391 │                          1.51 │                      0.387 │                            0.381 │                   2.94 │           0.891 │             1.82 │                               0.484 │
# └───────────────────┴──────────────────┴───────────────────────────────┴────────────────────────────┴──────────────────────────────────┴────────────────────────┴─────────────────┴──────────────────┴─────────────────────────────────────┘

# v1.6.0-beta1
# ┌───────────────────┬────────────────────────────────────┬────────────────────────────┬──────────────────────────────────┬────────────────────────┬──────────────────┬──────────────────┬─────────────────────────────────────┐
# │                   │ alg=EPGWork_ReIm_DualMVector_Split │ alg=EPGWork_ReIm_DualCache │ alg=EPGWork_ReIm_DualCache_Split │ alg=EPGWork_Basic_Cplx │ alg=EPGWork_Cplx │ alg=EPGWork_ReIm │ alg=EPGWork_ReIm_DualCache_Unrolled │
# │                   │                          Time [us] │                  Time [us] │                        Time [us] │              Time [us] │        Time [us] │        Time [us] │                           Time [us] │
# ├───────────────────┼────────────────────────────────────┼────────────────────────────┼──────────────────────────────────┼────────────────────────┼──────────────────┼──────────────────┼─────────────────────────────────────┤
# │ ETL=48, T=Float64 │                               0.79 │                      0.808 │                            0.824 │                   10.7 │             2.99 │             0.84 │                               0.919 │
# │ ETL=48, T=Float32 │                              0.814 │                      0.832 │                            0.816 │                   13.0 │             3.47 │             0.83 │                               0.875 │
# │ ETL=64, T=Float64 │                               1.43 │                       1.46 │                             1.44 │                   19.0 │             5.05 │             1.46 │                                1.54 │
# │ ETL=64, T=Float32 │                               1.43 │                       1.45 │                             1.43 │                   23.2 │              7.4 │             1.45 │                                1.49 │
# │ ETL=32, T=Float64 │                              0.389 │                      0.397 │                            0.388 │                   4.94 │             1.54 │            0.399 │                               0.465 │
# │ ETL=32, T=Float32 │                              0.381 │                      0.392 │                            0.381 │                   5.92 │             1.65 │            0.389 │                               0.431 │
# └───────────────────┴────────────────────────────────────┴────────────────────────────┴──────────────────────────────────┴────────────────────────┴──────────────────┴──────────────────┴─────────────────────────────────────┘

# v1.6.1
# ┌───────────────────┬─────────────────────────────────────────┬────────────────────────────────────┬────────────────────────────┬──────────────────────────────────┬────────────────────────────────────────┬────────────────────────┬─────────────────┬──────────────────┐
# │                   │ alg=EPGWork_ReIm_DualPaddedVector_Split │ alg=EPGWork_ReIm_DualMVector_Split │ alg=EPGWork_ReIm_DualCache │ alg=EPGWork_ReIm_DualCache_Split │ alg=EPGWork_ReIm_DualMVector_Vec_Split │ alg=EPGWork_Basic_Cplx │ alg=EPGWork_Vec │ alg=EPGWork_ReIm │
# │                   │                               Time [us] │                          Time [us] │                  Time [us] │                        Time [us] │                              Time [us] │              Time [us] │       Time [us] │        Time [us] │
# ├───────────────────┼─────────────────────────────────────────┼────────────────────────────────────┼────────────────────────────┼──────────────────────────────────┼────────────────────────────────────────┼────────────────────────┼─────────────────┼──────────────────┤
# │ ETL=48, T=Float64 │                                   0.986 │                              0.821 │                      0.818 │                            0.802 │                                   1.17 │                   11.6 │            1.79 │             0.82 │
# │ ETL=48, T=Float32 │                                     1.0 │                              0.793 │                       0.81 │                            0.794 │                                  0.954 │                   13.1 │            1.85 │              0.8 │
# │ ETL=64, T=Float64 │                                    1.71 │                               1.39 │                       1.41 │                             1.39 │                                   2.01 │                   20.2 │            3.02 │             1.41 │
# │ ETL=64, T=Float32 │                                    1.74 │                               1.37 │                       1.41 │                             1.39 │                                   1.68 │                   23.1 │             3.4 │             1.41 │
# │ ETL=32, T=Float64 │                                   0.459 │                              0.375 │                      0.384 │                            0.376 │                                  0.542 │                   5.25 │           0.847 │            0.387 │
# │ ETL=32, T=Float32 │                                   0.473 │                              0.369 │                      0.379 │                            0.368 │                                  0.444 │                    5.8 │           0.794 │            0.377 │
# └───────────────────┴─────────────────────────────────────────┴────────────────────────────────────┴────────────────────────────┴──────────────────────────────────┴────────────────────────────────────────┴────────────────────────┴─────────────────┴──────────────────┘


text:
#==========================================================
This routine initializes an optimization problem to recover
a |0⟩ to |2⟩ swap gate on a single qudit with 3 energy
levels (and 1 guard state). The  drift Hamiltonian in the
rotating frame is
        H0 = 2π*diagm([0 0 -2.2538e-1 -7.0425e-1]).
Here the control Hamiltonian includes the usual symmetric
and anti-symmetric terms
     H_{sym} = p(t)(a + a^†),    H_{asym} = q(t)(a - a^†),
where a is the annihilation operator for the qudit.
The problem parameters for this example are:
                ω_a =  2π × 4.09947    Grad/s,
                ξ_a =  2π × 2.2538e-01 Grad/s.
We use Bsplines with carrier waves with frequencies
0, ξ_a Grad/s.
==========================================================#
using LinearAlgebra
using Plots
using FFTW
using DelimitedFiles
using Printf
using Ipopt
using Random
using JLD2
using FastGaussQuadrature

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

using Juqbox # quantum control module

# Default values for uniform distribution on [-ϵ/2,ϵ/2]
if(!@isdefined(ep_max))
    ep_max = 2*pi*3e-2
end

# Default number of quadrature nodes to evaluate the
# expected value of the objective function via
# Gaussian quadrature, i.e.
#       E[J] = ∑ w[k] J[x[k]]
# where w,x are the weights and nodes on [-ϵ/2,ϵ/2]
if(!@isdefined(nquad))
    nquad = 10
    nquad = 1
end

if nquad == 1
    nodes, weights = gausslegendre(nquad)
    weights = weights./2.0
else
    nodes_tmp, weights_tmp = gausshermite(nquad)

    # Make a larger node list for full distribution
    n_modes = length(mean_vec)
    nodes = zeros(n_modes*nquad)
    weights = copy(nodes)
    inv_n = 1.0/(n_modes*sqrt(pi))

    # Mapping to reference weighting function
    for i = 1:n_modes
        μ = mean_vec[i]
        σ = sig_vec[i]

        offset = (i-1)*nquad
        for j = 1:length(nodes_tmp)
            nodes[j + offset] = sqrt(2)*σ*nodes_tmp[j] + μ
            weights[j+offset] = weights_tmp[j]*inv_n
        end
    end
end

N = 3 # Number of essential energy levels
Nguard = 1 # Number of guard/forbidden energy levels
Ntot = N + Nguard # Total number of energy levels

samplerate = 32 # for output files

T = 275.0 # Duration of gate

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10336
xa = 0.2198
rot_freq = [fa] # Used to calculate the lab frame ctrl function

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number) # xa is in GHz

utarget = zeros(ComplexF64,Ntot,N)

if N == 2
    utarget[2,1] = 1.0
    utarget[1,2] = 1.0
end

# 0>  to  |2>  swap gate
if N >= 3
    utarget[1,1] = 0 #1/sqrt(2)
    utarget[2,1] = 0
    utarget[3,1] = 1 #1/sqrt(2)
#
    utarget[1,2] = 0
    utarget[2,2] = 1
    utarget[3,2] = 0
    #
    utarget[1,3] = 1 #1/sqrt(2)
    utarget[2,3] = 0
    utarget[3,3] = 0 #-1/sqrt(2)
#
end

if N==4
    utarget[4,4] = 1
end

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], [fa])

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the lab frame
vtarget = utarget

startFromScratch = true
# startFile is used when startFromScratch = false
startFile="swap02-baseline-pcof-opt.jld2" # "swap02-pert-4em1-pcof-opt.jld2"

usePrior = false #  true
priorFileName = startFile # Usually makes sense to also use the start file as the prior, but not required

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard lowering matrix
adag = Array(transpose(amat));
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

Ncoupled = length(Hsym_ops) # Number of paired Hamiltonians
Nfreq= 2 # number of carrier frequencies 3 gives a cleaner sol than 2

# setup carrier frequencies
use_bcarrier = true
om = zeros(Ncoupled,Nfreq)
if use_bcarrier
    om[1:Ncoupled,2] .= -2.0*pi *xa       # Note negative sign
end
println("Carrier frequencies [GHz]: ", om[1,:]./(2*pi))
println("H0: ", H0)

maxctrl = 2*pi*1.2e-2

#max amplitude (in angular frequency) 2*pi*GHz
maxamp = zeros(Nfreq)
if Nfreq >= 5
    const_fact = 1.0/(Nfreq-2)
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # max B-spline coefficient amplitude, factor 3.0 is ad hoc
elseif Nfreq >= 4
    const_fact = 0.4
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # max B-spline coefficient amplitude, factor 3.0 is ad hoc
elseif Nfreq >= 3
    const_fact = 0.45
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # max B-spline coefficient amplitude, factor 3.0 is ad hoc
else
    maxamp .= maxctrl/Nfreq
end

maxpar = maximum(maxamp)
nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, [maxctrl])
println("# time steps: ", nsteps)

# Initial conditions for basis
Ident = Matrix{Float64}(UniformScaling(1), Ntot, Ntot)
U0 = Ident[1:Ntot,1:N]

# params = Juqbox.parameters([N], [Nguard], T, nsteps, U0, vtarget, om, H0, Hsym_ops, Hanti_ops)
params = Juqbox.objparams([N], [Nguard], T, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops, wmatScale=1.0)

if usePrior
    Juqbox.setup_prior!(params, priorFileName)
    params.tik0 = 100.0 # increase Tikhonov regularization coefficient
end

Random.seed!(2456)

# initial parameter guess
if startFromScratch
  D1 = 15 # Number of B-spline coefficients per frequency, sin/cos and real/imag
  nCoeff = 2*Ncoupled*Nfreq*D1 # factor '2' is for sin/cos
  pcof0 = (rand(nCoeff) .- 0.5).*maxpar*0.1

  if(nquad == 1)
    D1 = 15 # Number of B-spline coefficients per frequency, sin/cos and real/imag
    nCoeff = 2*Ncoupled*Nfreq*D1 # factor '2' is for sin/cos
    pcof0 = (rand(nCoeff) .- 0.5).*maxpar*0.1
  end
  println("*** Starting from random pcof with amplitude ", maxpar*0.1)
else
    # use if you want to have initial coefficients read from file
    @load startFile pcof
    pcof0 = pcof
    println("*** Starting from B-spline coefficients in file: ", startFile)
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Ncoupled*Nfreq)  # number of B-spline coeff per control function
end

# min and max coefficient values (set first and last two to zero)
useBarrier = true
minCoeff, maxCoeff = Juqbox.assign_thresholds_freq(maxamp, Ncoupled, Nfreq, D1)
zero_start_end!(params, D1, minCoeff, maxCoeff)

println("*** Settings ***")
println("System Hamiltonian coefficients [GHz]: (fa, xa) =  ", fa, xa)
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
println("Max parameter amplitudes: maxpar = ", maxpar)
println("Target rotating frame control amplitude, maxctrl = ", maxctrl, " [rad/ns], ", maxctrl*0.5/pi, " [GHz]")
# params.tik0 = 0
println("Tikhonov coefficients: tik0 = ", params.tik0)

# optional arguments to setup_ipopt_problem()
maxIter = 100
lbfgsMax = 5
ipTol = 1e-5
acceptTol = 1e-5
acceptIter = 15

# Estimate number of terms in Neumann series for time stepping (Default 3)
tol = eps(1.0); # machine precision
Juqbox.estimate_Neumann!(tol, params, [maxpar])

# Allocate all working arrays
wa = Juqbox.Working_Arrays(params,nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter, lbfgsMax, startFromScratch, ipTol,acceptTol, acceptIter, nodes, weights)

println("Initial coefficient vector stored in 'pcof0'")


text:
module TrajPlanning

export TrajPlanningProblem, plan_trajectory, x_path_from_u

using OneVision
using OneVision: ℝ, ℕ, 𝕋, @kwdef
using StaticArrays
using Optim

struct TrajPlanningProblem{H,ΔT,X,U,Dy <: SysDynamics}
    horizon::Val{H}
    control_inverval::Val{ΔT}
    dy::Dy
    x_weights::X
    u_weights::U
    optim_options::Optim.Options
    cache::Ref{Any}
end

TrajPlanningProblem(H::ℕ, ΔT::𝕋, dy, x_weights, u_weights;
    optim_options = Optim.Options(iterations = 100)) =
    TrajPlanningProblem(Val(H), Val(ΔT), dy, x_weights, u_weights,
        optim_options, Ref{Any}(missing))

const check_optimizer_converge = Ref(true)

"""
Returns `(u⋆, objective_value)` where `u⋆[1] = u⋆(τ)`, `length(u⋆) = H`.

# Arguments
    - `x0[τ]`: state at time τ
    - `x_path[t in τ+1: τ+HΔT]`: state trajectory ∈ [τ+1, τ+HΔT]
    - `u_path[t in τ: τ+HΔT-1]`: action trajectory ∈ [τ, τ+HΔT-1]
"""
function plan_trajectory(
    p::TrajPlanningProblem{H,ΔT,X,U,Dy},
    x0::X,
    x_path::SVector{H_ΔT,X},
    u_path::SVector{H_ΔT,U},
    τ::𝕋
)::Tuple where {H,ΔT,X,U,Dy,H_ΔT}
    @assert H_ΔT == H * ΔT

    function loss(uvec::AbstractVector{R})::R where R
        XR, UR = similar_type(X, R), similar_type(U, R)
        us = reinterpret(UR, uvec)
        x::XR = x0
        l_x::XR = zero(XR)
        l_u::UR = zero(UR)

        @inbounds for j = 0:H - 1
            local u0::UR = us[j+1]
            local u::UR = limit_control(p.dy, u0, x, τ + j * ΔT)
            for k = 0:ΔT - 1
                t = j * ΔT + k
                x = sys_forward(p.dy, x, u, τ + t)
                l_x += (x - x_path[t + 1]).^2
                l_u += (u0 - u_path[t + 1]).^2
            end
        end
        (sum(l_x .* p.x_weights) + sum(l_u .* p.u_weights)) / ΔT
    end

    use_warmstart = true

    u0::Vector{ℝ} = (use_warmstart && !ismissing(p.cache[])) ?
                    p.cache[] : zeros(length(U) * H)

    res = Optim.optimize(
        loss,
        u0, Optim.LBFGS(),
        p.optim_options,
        autodiff = :forward,
    )
    if check_optimizer_converge[] && !Optim.converged(res)
        @warn "Optim not converged: $res"
    end

    uvec = Optim.minimizer(res)
    n_u = length(U)
    # û = @SVector[U(uvec[1+n_u*i:n_u*(i+1)]) for i in 1:H]
    û = SVector{H}([U(uvec[i:i+n_u-1]...) for i in 1:n_u:H*n_u])

    if use_warmstart
        uvec[1:end - n_u] = uvec[1 + n_u:end]
        p.cache[] = uvec
    end

    @assert (obj = Optim.minimum(res)) ≥ -sqrt(eps()) "negative objective: $(obj)"
    û, obj
end

function x_path_from_u(x0::X, t0::𝕋, û::SVector{H,U}, dy, ::Val{ΔT}) where {X,U,H,ΔT}
    H_ΔT = H * ΔT
    x̂ = MVector{H_ΔT,X}(undef)
    x::X = x0
    for i = 0:H - 1
        for j = 0:ΔT - 1
            t = t0 + i * ΔT + j
            x = sys_forward(dy, x, û[i + 1], t)
            x̂[1 + i * ΔT + j] = x
        end
    end
    x̂
end

end # module TrajPlanning

@reexport using .TrajPlanning


text:
int main()
{
	printDouble(3.0);
	return 10;
};

void f()
{
	return;
}

void g()
{
	int x = 0;
	printInt(x);
	return;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
struct Interval
    min::Float64
    max::Float64
end

Base.minimum(iv::Interval) = iv.min
Base.maximum(iv::Interval) = iv.max

Base.in(x::Real, iv::Interval) = (x <= iv.max) & (x >= iv.min)

Base.issubset(iv1::Interval, iv2::Interval) = (iv1.max <= iv2.max) & (iv1.min >= iv2.min)

Base.isempty(iv::Interval) = iv.max < iv.min

const ∅ = Interval(1, -1)

Base.:(==)(iv1::Interval, iv2::Interval) = ((iv1.max == iv2.max) & (iv1.min == iv2.min)) | ( isempty(iv1) & isempty(iv2) )

Base.intersect(iv1::Interval, iv2::Interval) = Interval(max(iv1.min,iv2.min), min(iv1.max,iv2.max))

function Base.show(io::IO, mime::MIME"text/plain", iv::Interval)
    if isempty(iv)
        print(io,"∅")
    else
        print(io,"\u301a$(iv.min), $(iv.max)\u301b")
    end
end


text:
using Pkg

# Install required packages
for P in ["HTTP", "PyPlot"]
  !haskey(Pkg.installed(), P) && Pkg.add(P)
end

useCompressedMatrices = true
suffixSM = useCompressedMatrices ? "Small" : "Large"
useFastData = true
suffixMeas = useFastData ? "Fast" : "Slow"

# Download data
include("downloadData.jl")

using MPIReco, PyPlot

################
## Parameters ##
################
alpha = 1.2

# Choose the peak number according to the expected approximate value. Be careful, if there are peaks not corresponding to object motion, division by choosePeak can be false.
# For experiments with different periodic motions it is used to pick one motion frequency (then don't divide by choosePeak)

choosePeak = 4
windowType = 1  # 1: Hann, 2: FT1A05, 3: Rectangle

# Measurement data
datadirMeas = "./data/"
bMeas = MPIFile(datadirMeas*"meas$(suffixMeas).mdf") # high frequency
bBG = MPIFile(datadirMeas*"measBG.mdf") # background measurement

# System matrices
datadirSF = "./data/"
SFall = ["SF$(d)$(suffixSM).mdf" for d=1:4]
bSF = MultiMPIFile(datadirSF.*SFall)

# Background frames
bgFrames = [1:200, 201:400, 401:600, 601:800]

# Selected frames
frames = 1:2

####################
## Reconstruction ##
####################

# Reconstruction parameters
lambda = 0.01
iterations = 1

c = reconstruction(bSF, bMeas, periodicMotionCorrection = true,
                              bEmpty = bBG, bgFrames = bgFrames, #background measurement
                              alpha = alpha, choosePeak = choosePeak, frames = frames,
                              windowType=windowType, higherHarmonic=choosePeak,
                              lambda=lambda, iterations=iterations, # reconstruction parameter
                              SNRThresh=2,minFreq=80e3,recChannels=[1,2,3]) #frequency selection

figure(2)
imshow( maximum(arraydata(c[1,:,:,:,1]),dims=3)[:,:,1] )

##################
### Static Data ##
##################

bStatic = MPIFile(datadirMeas*"measStatic.mdf") # static measurement

lambda = 0.01
iterations = 5

cStatic = reconstructionMultiPatch(bSF, bStatic, bgCorrection=false,
			      SNRThresh=2, minFreq=80e3,recChannels=[1,2,3],
                              λ=lambda, iterations=iterations)

figure(3)
imshow( maximum(arraydata(cStatic[1,:,:,:,1]),dims=3)[:,:,1] )


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

@doc raw"""
    ObjsTeam(;
        archived=nothing,
        avatar_base_url=nothing,
        created=nothing,
        date_create=nothing,
        deleted=nothing,
        discoverable=nothing,
        domain=nothing,
        email_domain=nothing,
        enterprise_id=nothing,
        enterprise_name=nothing,
        external_org_migrations=nothing,
        has_compliance_export=nothing,
        icon=nothing,
        id=nothing,
        is_assigned=nothing,
        is_enterprise=nothing,
        is_over_storage_limit=nothing,
        limit_ts=nothing,
        locale=nothing,
        messages_count=nothing,
        msg_edit_window_mins=nothing,
        name=nothing,
        over_integrations_limit=nothing,
        over_storage_limit=nothing,
        pay_prod_cur=nothing,
        plan=nothing,
        primary_owner=nothing,
        sso_provider=nothing,
    )

    - archived::Bool
    - avatar_base_url::String
    - created::Int32
    - date_create::Int32
    - deleted::Bool
    - discoverable::Any
    - domain::String
    - email_domain::String
    - enterprise_id::DefsEnterpriseId
    - enterprise_name::DefsEnterpriseName
    - external_org_migrations::ObjsExternalOrgMigrations
    - has_compliance_export::Bool
    - icon::ObjsIcon
    - id::DefsWorkspaceId
    - is_assigned::Bool
    - is_enterprise::Int32
    - is_over_storage_limit::Bool
    - limit_ts::Int32
    - locale::String
    - messages_count::Int32
    - msg_edit_window_mins::Int32
    - name::String
    - over_integrations_limit::Bool
    - over_storage_limit::Bool
    - pay_prod_cur::String
    - plan::String
    - primary_owner::ObjsPrimaryOwner
    - sso_provider::ObjsTeamSsoProvider
"""
mutable struct ObjsTeam <: SwaggerModel
    archived::Any # spec type: Union{ Nothing, Bool } # spec name: archived
    avatar_base_url::Any # spec type: Union{ Nothing, String } # spec name: avatar_base_url
    created::Any # spec type: Union{ Nothing, Int32 } # spec name: created
    date_create::Any # spec type: Union{ Nothing, Int32 } # spec name: date_create
    deleted::Any # spec type: Union{ Nothing, Bool } # spec name: deleted
    discoverable::Any # spec type: Union{ Nothing, Any } # spec name: discoverable
    domain::Any # spec type: Union{ Nothing, String } # spec name: domain
    email_domain::Any # spec type: Union{ Nothing, String } # spec name: email_domain
    enterprise_id::Any # spec type: Union{ Nothing, DefsEnterpriseId } # spec name: enterprise_id
    enterprise_name::Any # spec type: Union{ Nothing, DefsEnterpriseName } # spec name: enterprise_name
    external_org_migrations::Any # spec type: Union{ Nothing, ObjsExternalOrgMigrations } # spec name: external_org_migrations
    has_compliance_export::Any # spec type: Union{ Nothing, Bool } # spec name: has_compliance_export
    icon::Any # spec type: Union{ Nothing, ObjsIcon } # spec name: icon
    id::Any # spec type: Union{ Nothing, DefsWorkspaceId } # spec name: id
    is_assigned::Any # spec type: Union{ Nothing, Bool } # spec name: is_assigned
    is_enterprise::Any # spec type: Union{ Nothing, Int32 } # spec name: is_enterprise
    is_over_storage_limit::Any # spec type: Union{ Nothing, Bool } # spec name: is_over_storage_limit
    limit_ts::Any # spec type: Union{ Nothing, Int32 } # spec name: limit_ts
    locale::Any # spec type: Union{ Nothing, String } # spec name: locale
    messages_count::Any # spec type: Union{ Nothing, Int32 } # spec name: messages_count
    msg_edit_window_mins::Any # spec type: Union{ Nothing, Int32 } # spec name: msg_edit_window_mins
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    over_integrations_limit::Any # spec type: Union{ Nothing, Bool } # spec name: over_integrations_limit
    over_storage_limit::Any # spec type: Union{ Nothing, Bool } # spec name: over_storage_limit
    pay_prod_cur::Any # spec type: Union{ Nothing, String } # spec name: pay_prod_cur
    plan::Any # spec type: Union{ Nothing, String } # spec name: plan
    primary_owner::Any # spec type: Union{ Nothing, ObjsPrimaryOwner } # spec name: primary_owner
    sso_provider::Any # spec type: Union{ Nothing, ObjsTeamSsoProvider } # spec name: sso_provider

    function ObjsTeam(;archived=nothing, avatar_base_url=nothing, created=nothing, date_create=nothing, deleted=nothing, discoverable=nothing, domain=nothing, email_domain=nothing, enterprise_id=nothing, enterprise_name=nothing, external_org_migrations=nothing, has_compliance_export=nothing, icon=nothing, id=nothing, is_assigned=nothing, is_enterprise=nothing, is_over_storage_limit=nothing, limit_ts=nothing, locale=nothing, messages_count=nothing, msg_edit_window_mins=nothing, name=nothing, over_integrations_limit=nothing, over_storage_limit=nothing, pay_prod_cur=nothing, plan=nothing, primary_owner=nothing, sso_provider=nothing)
        o = new()
        validate_property(ObjsTeam, Symbol("archived"), archived)
        setfield!(o, Symbol("archived"), archived)
        validate_property(ObjsTeam, Symbol("avatar_base_url"), avatar_base_url)
        setfield!(o, Symbol("avatar_base_url"), avatar_base_url)
        validate_property(ObjsTeam, Symbol("created"), created)
        setfield!(o, Symbol("created"), created)
        validate_property(ObjsTeam, Symbol("date_create"), date_create)
        setfield!(o, Symbol("date_create"), date_create)
        validate_property(ObjsTeam, Symbol("deleted"), deleted)
        setfield!(o, Symbol("deleted"), deleted)
        validate_property(ObjsTeam, Symbol("discoverable"), discoverable)
        setfield!(o, Symbol("discoverable"), discoverable)
        validate_property(ObjsTeam, Symbol("domain"), domain)
        setfield!(o, Symbol("domain"), domain)
        validate_property(ObjsTeam, Symbol("email_domain"), email_domain)
        setfield!(o, Symbol("email_domain"), email_domain)
        validate_property(ObjsTeam, Symbol("enterprise_id"), enterprise_id)
        setfield!(o, Symbol("enterprise_id"), enterprise_id)
        validate_property(ObjsTeam, Symbol("enterprise_name"), enterprise_name)
        setfield!(o, Symbol("enterprise_name"), enterprise_name)
        validate_property(ObjsTeam, Symbol("external_org_migrations"), external_org_migrations)
        setfield!(o, Symbol("external_org_migrations"), external_org_migrations)
        validate_property(ObjsTeam, Symbol("has_compliance_export"), has_compliance_export)
        setfield!(o, Symbol("has_compliance_export"), has_compliance_export)
        validate_property(ObjsTeam, Symbol("icon"), icon)
        setfield!(o, Symbol("icon"), icon)
        validate_property(ObjsTeam, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsTeam, Symbol("is_assigned"), is_assigned)
        setfield!(o, Symbol("is_assigned"), is_assigned)
        validate_property(ObjsTeam, Symbol("is_enterprise"), is_enterprise)
        setfield!(o, Symbol("is_enterprise"), is_enterprise)
        validate_property(ObjsTeam, Symbol("is_over_storage_limit"), is_over_storage_limit)
        setfield!(o, Symbol("is_over_storage_limit"), is_over_storage_limit)
        validate_property(ObjsTeam, Symbol("limit_ts"), limit_ts)
        setfield!(o, Symbol("limit_ts"), limit_ts)
        validate_property(ObjsTeam, Symbol("locale"), locale)
        setfield!(o, Symbol("locale"), locale)
        validate_property(ObjsTeam, Symbol("messages_count"), messages_count)
        setfield!(o, Symbol("messages_count"), messages_count)
        validate_property(ObjsTeam, Symbol("msg_edit_window_mins"), msg_edit_window_mins)
        setfield!(o, Symbol("msg_edit_window_mins"), msg_edit_window_mins)
        validate_property(ObjsTeam, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ObjsTeam, Symbol("over_integrations_limit"), over_integrations_limit)
        setfield!(o, Symbol("over_integrations_limit"), over_integrations_limit)
        validate_property(ObjsTeam, Symbol("over_storage_limit"), over_storage_limit)
        setfield!(o, Symbol("over_storage_limit"), over_storage_limit)
        validate_property(ObjsTeam, Symbol("pay_prod_cur"), pay_prod_cur)
        setfield!(o, Symbol("pay_prod_cur"), pay_prod_cur)
        validate_property(ObjsTeam, Symbol("plan"), plan)
        setfield!(o, Symbol("plan"), plan)
        validate_property(ObjsTeam, Symbol("primary_owner"), primary_owner)
        setfield!(o, Symbol("primary_owner"), primary_owner)
        validate_property(ObjsTeam, Symbol("sso_provider"), sso_provider)
        setfield!(o, Symbol("sso_provider"), sso_provider)
        o
    end
end # type ObjsTeam

const _property_map_ObjsTeam = Dict{Symbol,Symbol}(Symbol("archived")=>Symbol("archived"), Symbol("avatar_base_url")=>Symbol("avatar_base_url"), Symbol("created")=>Symbol("created"), Symbol("date_create")=>Symbol("date_create"), Symbol("deleted")=>Symbol("deleted"), Symbol("discoverable")=>Symbol("discoverable"), Symbol("domain")=>Symbol("domain"), Symbol("email_domain")=>Symbol("email_domain"), Symbol("enterprise_id")=>Symbol("enterprise_id"), Symbol("enterprise_name")=>Symbol("enterprise_name"), Symbol("external_org_migrations")=>Symbol("external_org_migrations"), Symbol("has_compliance_export")=>Symbol("has_compliance_export"), Symbol("icon")=>Symbol("icon"), Symbol("id")=>Symbol("id"), Symbol("is_assigned")=>Symbol("is_assigned"), Symbol("is_enterprise")=>Symbol("is_enterprise"), Symbol("is_over_storage_limit")=>Symbol("is_over_storage_limit"), Symbol("limit_ts")=>Symbol("limit_ts"), Symbol("locale")=>Symbol("locale"), Symbol("messages_count")=>Symbol("messages_count"), Symbol("msg_edit_window_mins")=>Symbol("msg_edit_window_mins"), Symbol("name")=>Symbol("name"), Symbol("over_integrations_limit")=>Symbol("over_integrations_limit"), Symbol("over_storage_limit")=>Symbol("over_storage_limit"), Symbol("pay_prod_cur")=>Symbol("pay_prod_cur"), Symbol("plan")=>Symbol("plan"), Symbol("primary_owner")=>Symbol("primary_owner"), Symbol("sso_provider")=>Symbol("sso_provider"))
const _property_types_ObjsTeam = Dict{Symbol,String}(Symbol("archived")=>"Bool", Symbol("avatar_base_url")=>"String", Symbol("created")=>"Int32", Symbol("date_create")=>"Int32", Symbol("deleted")=>"Bool", Symbol("discoverable")=>"Any", Symbol("domain")=>"String", Symbol("email_domain")=>"String", Symbol("enterprise_id")=>"DefsEnterpriseId", Symbol("enterprise_name")=>"DefsEnterpriseName", Symbol("external_org_migrations")=>"ObjsExternalOrgMigrations", Symbol("has_compliance_export")=>"Bool", Symbol("icon")=>"ObjsIcon", Symbol("id")=>"DefsWorkspaceId", Symbol("is_assigned")=>"Bool", Symbol("is_enterprise")=>"Int32", Symbol("is_over_storage_limit")=>"Bool", Symbol("limit_ts")=>"Int32", Symbol("locale")=>"String", Symbol("messages_count")=>"Int32", Symbol("msg_edit_window_mins")=>"Int32", Symbol("name")=>"String", Symbol("over_integrations_limit")=>"Bool", Symbol("over_storage_limit")=>"Bool", Symbol("pay_prod_cur")=>"String", Symbol("plan")=>"String", Symbol("primary_owner")=>"ObjsPrimaryOwner", Symbol("sso_provider")=>"ObjsTeamSsoProvider")
Base.propertynames(::Type{ ObjsTeam }) = collect(keys(_property_map_ObjsTeam))
Swagger.property_type(::Type{ ObjsTeam }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsTeam[name]))}
Swagger.field_name(::Type{ ObjsTeam }, property_name::Symbol) =  _property_map_ObjsTeam[property_name]

const _allowed_ObjsTeam_plan = ["", "std", "plus", "compliance", "enterprise"]

function check_required(o::ObjsTeam)
    (getproperty(o, Symbol("domain")) === nothing) && (return false)
    (getproperty(o, Symbol("email_domain")) === nothing) && (return false)
    (getproperty(o, Symbol("icon")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsTeam }, name::Symbol, val)
    if name === Symbol("plan")
        Swagger.validate_param(name, "ObjsTeam", :enum, val, _allowed_ObjsTeam_plan)
    end
end


text:
using DynamicHMCModels

Random.seed!(1)

cd(@__DIR__)
include("sdt_functions.jl")
include("simulateSDT.jl")

struct SDTProblem
    hits::Int64
    fas::Int64
    Nd::Int64
end

function (problem::SDTProblem)(θ)
    @unpack hits,fas,Nd=problem   # extract the data
    @unpack d,c=θ
    logpdf(SDT(d,c),[hits,fas,Nd])+logpdf(Normal(0,1/sqrt(2)),d) +
    logpdf(Normal(0,1/sqrt(2)),c)
end

# Define problem with data and inits.

data = simulateSDT(;Nd=100)

p = SDTProblem(data...)
p((d=2.0,c=.0))

# Write a function to return properly dimensioned transformation.
problem_transformation(p::SDTProblem) =
     as((d=asℝ,c=asℝ))

# Use Flux for the gradient.
P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P);

#import Zygote
#∇P = ADgradient(:Zygote, P);

# FSample from the posterior.
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 4000);

# Undo the transformation to obtain the posterior from the chain.
posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
  get_position.(chain));

chns = nptochain(posterior,NUTS_tuned)


text:
# Autogenerated wrapper script for CUTENSOR_jll for x86_64-linux-gnu-cuda+10.2
export libcutensor

using CUDA_loader_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CUTENSOR")
JLLWrappers.@declare_library_product(libcutensor, "libcutensor.so.1")
function __init__()
    JLLWrappers.@generate_init_header(CUDA_loader_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcutensor,
        "lib/libcutensor.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
type ModPC <: PairCop
    cop::PairCop
    mod::PCModification
end

##########################
## function evaluations ##
##########################

function pdf(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return pdf(cop.cop, cop.mod, u1, u2)
end

function cdf(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return cdf(cop.cop, cop.mod, u1, u2)
end

function hfun(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return hfun(cop.cop, cop.mod, u1, u2)
end

function vfun(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return vfun(cop.cop, cop.mod, u1, u2)
end

function hinv(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return hinv(cop.cop, cop.mod, u1, u2)
end

function vinv(cop::ModPC,
             u1::FloatVec, u2::FloatVec)
    return vinv(cop.cop, cop.mod, u1, u2)
end


text:
module Deploy

module Docker

using Revise
using Genie, Genie.FileTemplates

function dockerfile(path::String = "."; user::String = "genie", env::String = "dev",
                    filename::String = "Dockerfile", port::Int = 8000, dockerport::Int = 80, force::Bool = false)
  filename = normpath(joinpath(path, filename))
  isfile(filename) && force && rm(filename)
  isfile(filename) && error("File $(filename) already exists. Use the `force = true` option to overwrite the existing file.")

  open(filename, "w") do io
    write(io, FileTemplates.dockerfile(user = user, env = env, filename = filename,
                                        port = port, dockerport = dockerport))
  end

  "Docker file successfully written at $(abspath(filename))" |> println
end

function build(path::String = "."; appname = "genie")
  `docker build -t "$appname" $path` |> Base.run

  "Docker container successfully built" |> println
end

function run(; containername::String = "genieapp", hostport::Int = 80, containerport::Int = 8000, appdir::String = "/home/genie/app",
                mountapp::Bool = false, image::String = "genie", command::String = "bin/server", rm::Bool = true, it::Bool = true)
  options = []
  it && push!(options, "-it")
  rm && push!(options, "--rm")
  push!(options, "-p")
  push!(options, "$hostport:$containerport")
  push!(options, "--name")
  push!(options, "$containername")
  if mountapp
    push!(options, "-v")
    push!(options,  "$(pwd()):$appdir")
  end
  push!(options, image)
  isempty(command) || push!(options, command)

  "Starting docker container with `docker run $options`" |> println

  `docker run $options` |> Base.run
end

end # end module Docker

########

module Heroku

function createapp(appname::String; region::String = "us")
  `heroku create $(lowercase(appname)) --region $region` |> Base.run
end

function push(appname::String; apptype::String = "web")
  `heroku container:push $apptype -a $appname` |> Base.run
end

function release(appname::String; apptype::String = "web")
  `heroku container:release $apptype -a $appname` |> Base.run
end

function open(appname::String)
  `heroku open -a $appname` |> Base.run
end

function login()
  `heroku container:login` |> Base.run
end

function logs(appname::String; lines::Int = 1_000)
  `heroku logs --tail -a $appname -n $lines` |> Base.run
end

end # end module Heroku

end # end module Deploy


text:
# Script to replicate table 3 (regression from equation 12)

using DataFrames, FixedEffectModels, RegressionTables, StatsBase, Statistics, LinearAlgebra, XLSX, PrettyTables

# -------------- Import data --------------------------------------------------------------------------------------------------------------------------------

# import data from local drive
t_tab3 = DataFrame(XLSX.readtable("clean/t_tab2.xlsx", "Sheet1")...) # still a lot of NaN
transform!(t_tab3, :country .=> ByRow(string) .=> :country, renamecols=false) # change type to String
transform!(t_tab3, [:year, :sector] .=> ByRow(Int) .=> [:year, :sector], renamecols=false) # change type to Int64
transform!(t_tab3, [:FD_GO, :VA_GO, :U, :D] .=> ByRow(Float64) .=> [:FD_GO, :VA_GO, :U, :D], renamecols=false) # change type to Float64

# need to take out rows with NaN, considerably improves match with original numbers (too much correlation otherwise)
data_reg = subset(t_tab3, :VA_GO => ByRow(x -> !isnan(x)))
data_reg[:, :country_sector] .= data_reg[:, :country] .* "_" .* string.(data_reg[:, :sector]) # using seperate FE does not change anything?

function reg_tab3(df::DataFrame, year::Int64, y, x)
    df = subset(df, :year => ByRow(x -> x == year))

    rr_1 = FixedEffectModels.reg(df, term(y) ~ term(x), Vcov.cluster(:country_sector), save=true)
    rr_2 = FixedEffectModels.reg(df, term(y) ~ term(x) + fe(term(:country)), Vcov.cluster(:country_sector), save=true)
    rr_3 = FixedEffectModels.reg(df, term(y) ~ term(x) + fe(term(:country)) + fe(term(:sector)), Vcov.cluster(:country_sector), save=true)

    reg = RegressionTables.regtable(rr_1, rr_2, rr_3;
    renderSettings = asciiOutput(), regression_statistics=[:nobs, :r2], print_fe_section=true, estimformat="%0.4f") # output on the console

    return reg
end

# regression results in ascii format
reg_tab3(data_reg, 1995, :FD_GO, :VA_GO)
reg_tab3(data_reg, 2011, :FD_GO, :VA_GO)
reg_tab3(data_reg, 1995, :U, :D)
reg_tab3(data_reg, 2011, :U, :D)


text:
# @copyright (c) 2019 RWTH Aachen. All rights reserved.
#
# @file docs/make.jl
# @version 1.0.0
# @author Aleksandr Mikhalev
# @date 2019-12-10

using Documenter
using Maxvol

makedocs(
    sitename = "Maxvol",
    format = Documenter.HTML(),
    modules = [Maxvol]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
# ---
# title: 875. Koko Eating Bananas
# id: problem875
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/koko-eating-bananas/description/>
# hidden: true
# ---
#
# Koko loves to eat bananas.  There are `N` piles of bananas, the `i`-th pile
# has `piles[i]` bananas.  The guards have gone and will come back in `H` hours.
#
# Koko can decide her bananas-per-hour eating speed of `K`.  Each hour, she
# chooses some pile of bananas, and eats K bananas from that pile.  If the pile
# has less than `K` bananas, she eats all of them instead, and won't eat any
# more bananas during this hour.
#
# Koko likes to eat slowly, but still wants to finish eating all the bananas
# before the guards come back.
#
# Return the minimum integer `K` such that she can eat all the bananas within
# `H` hours.
#
#
#
# **Example 1:**
#
#
#
#     Input: piles = [3,6,7,11], H = 8
#     Output: 4
#
#
# **Example 2:**
#
#
#
#     Input: piles = [30,11,23,4,20], H = 5
#     Output: 30
#
#
# **Example 3:**
#
#
#
#     Input: piles = [30,11,23,4,20], H = 6
#     Output: 23
#
#
#
#
# **Constraints:**
#
#   * `1 <= piles.length <= 10^4`
#   * `piles.length <= H <= 10^9`
#   * `1 <= piles[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Types
    import ...SSGrB: Abstract_GrB_Struct, AbstractUnaryOp, AbstractMonoid, AbstractSelectOp,
    AbstractSemiring, AbstractBinaryOp, AbstractDescriptor
    using ..libgb
end

mutable struct GBScalar{T} <: Abstract_GrB_Struct
    p::libgb.GxB_Scalar
    function GBScalar{T}(p::libgb.GxB_Scalar) where {T}
        s = new(p)
        function f(scalar)
            libgb.GxB_Scalar_free(Ref(scalar.p))
        end
        return finalizer(f, s)
    end
end

mutable struct GBVector{T} <: AbstractSparseArray{T, UInt64, 1}
    p::libgb.GrB_Vector
    function GBVector{T}(p::libgb.GrB_Vector) where {T}
        v = new(p)
        function f(vector)
            libgb.GrB_Vector_free(Ref(vector.p))
        end
        return finalizer(f, v)
    end
end

mutable struct GBMatrix{T} <: AbstractSparseArray{T, UInt64, 2}
    p::libgb.GrB_Matrix
    function GBMatrix{T}(p::libgb.GrB_Matrix) where {T}
        A = new(p)
        function f(matrix)
            libgb.GrB_Matrix_free(Ref(matrix.p))
        end
        return finalizer(f, A)
    end
end

const GBVecOrMat = Union{GBVector{T}, GBMatrix{T}} where T
const OpUnion = Union{
    libgb.GrB_Monoid,
    libgb.GrB_UnaryOp,
    libgb.GrB_Semiring,
    libgb.GrB_BinaryOp,
    libgb.GxB_SelectOp
}

const MonoidSemiringOrBinary = Union{
    libgb.GrB_Monoid,
    libgb.GrB_Semiring,
    libgb.GrB_BinaryOp,
    AbstractSemiring,
    AbstractBinaryOp,
    AbstractMonoid
}
const GBMatOrTrans = Union{<:GBMatrix, Transpose{<:Any, <:GBMatrix}}
const GBVecMatOrTrans = Union{<:GBVector, GBMatOrTrans}

function getoperator(op, t)
    #Default Semiring should be LOR_LAND for boolean
    if op == Semirings.PLUS_TIMES_SEMIRING
        if t == Bool
            op = Semirings.LOR_LAND_SEMIRING
        end
    end
    #Default BinaryOp should be LAND or LOR for boolean
    if op == BinaryOps.TIMES
        if t == Bool
            op = BinaryOps.LAND
        end
    elseif op == BinaryOps.PLUS
        if t == Bool
            op = BinaryOps.LOR
        end
    end
    #Default monoid should be LAND/LOR
    if op == Monoids.TIMES_MONOID
        if t == Bool
            op = Monoids.LAND_MONOID
        end
    elseif op == Monoids.PLUS_MONOID
        if t == Bool
            op = Monoids.LOR_MONOID
        end
    end

    if op isa AbstractOp
        return op[t]
    elseif op isa OpUnion
        return op
    else
        error("Not a valid GrB op/semiring")
    end
end