text:
using Sockets

"""
    @msg(expression)

Calls `msg` with the name of the top-level function in `expression` as the first
argument and the arguments of that function as the remaining arguments. Note that
the toplevel function call will *not* be evaluated.

Example:
`@msg(sin(3))` will call `msg("sin", 3)`, i.e.
```
macroexpand(:(Atom.@msg(sin(3)))) == :((Atom.msg)("sin",3))
```
"""
macro msg(ex)
  @capture(ex, f_(args__)) || error("@msg requires function call syntax")
  :(msg($(string(f)), $(map(esc, args)...)))
end

"""
    @rpc(expression)

Calls `rpc` with the name of the top-level function in `expression` as the first
argument and the arguments of that function as the remaining arguments. Note that
the toplevel function call will *not* be evaluated.

Example:
`@rpc(sin(3))` will call `rpc("sin", 3)`, i.e.
```
macroexpand(:(Atom.@rpc(sin(3)))) == :((Atom.rpc)("sin",3))
```
"""
macro rpc(ex)
  @capture(ex, f_(args__)) || error("@rpc requires function call syntax")
  :(rpc($(string(f)), $(map(esc, args)...)))
end

"""
    sock

Socket over which the communication between server and client happens. Initialized
by `serve`.
"""
global sock = nothing

isactive(sock::Nothing) = false

"""
    isactive(sock)

Check whether the socked `sock` exists and is open.
"""
isactive(sock) = isopen(sock)

"""
    @ierrs(expression)

Catches errors internal to Atom.jl and renders them in the client.
"""
macro ierrs(ex)
  :(try
      $(esc(ex))
    catch e
      ee = EvalError(e, stacktrace(catch_backtrace()))
      Atom.msg("error", Dict(:msg         => "Julia Client – Internal Error",
                             :detail      => string(ee),
                             :dismissable => true))
      nothing
    end)
end

"""
    initialise(; welcome = false)

Sets up the environment for Atom.jl: Stop `SIGINT`s from killing Julia and send
a welcome message to Atom if `welcome == true`.
"""
function initialise(; welcome = false)
  Juno.isprecompiling() && return
  Juno.setactive!(true)
  exit_on_sigint(false)
  welcome && @msg welcome()
end

"""
    serve(port; kws...)

Start an asynchronous Julia server that listens on `localhost:port` and handle
the messages sent over that socket. The `kws...` are passed through to `initialise`.
"""
function serve(port; kws...)
  server = listen(ip"127.0.0.1", port)
  print(stderr, "juno-msg-ready")
  global sock = accept(server)
  @async while isopen(sock)
    @ierrs let
      msg = JSON.parse(sock)
      @async @ierrs handlemsg(msg...)
    end
  end
  initialise(; kws...)
end

function connect(host, port; kws...)
  global sock = Sockets.connect(host, port)
  @async while isopen(sock)
    @ierrs let
      msg = JSON.parse(sock)
      @async @ierrs handlemsg(msg...)
    end
  end
  initialise(; kws...)
end

connect(port; kws...) = connect(ip"127.0.0.1", port; kws...)

"""
    msg(typ, args...)

Prints a json string to `sock` that contains all function arguments in an array.
"""
function msg(t, args...)
  isactive(sock) || return
  println(sock, json(c(t, args...)))
end

const handlers = Dict{String, Function}()

"""
    handle(func, typ)

Set the handler for messages of type `typ` to `func`.
"""
handle(f, t) = handlers[t] = f

id = 0
const callbacks = Dict{Int,Condition}()

"""
    rpc(typ, args...)

Sends a message to `sock` via `msg`. Blocks until the client returns a message
(which needs to send back the `callback` field).
"""
function rpc(t, args...)
  i, c = (global id += 1), Condition()
  callbacks[i] = c
  msg(d(:type => t, :callback => i), args...)
  return wait(c)
end

"""
    handlemsg(typ, args...)

Tries to call the message handler corresponding to `typ` (which can be set with
`handle`).
"""
function handlemsg(t, args...)
  callback = nothing
  isa(t, AbstractDict) && ((t, callback) = (t["type"], t["callback"]))
  if haskey(handlers, t)
    try
      result = handlers[t](args...)
      isa(callback, Integer) && msg("cb", callback, result)
    catch e
      isa(callback, Integer) && msg("cancelCallback", callback, string(e))
      rethrow()
    end
  else
    @warn("Atom.jl: unrecognised message $t.")
    callback ≠ nothing && msg("cancelCallback", callback)
  end
end

handle("cb") do id, result
  notify(callbacks[id], result)
  delete!(callbacks, id)
end

handle(() -> nothing, "junorc")

isconnected() = sock ≠ nothing && isopen(sock)


text:
include("Track.jl")
include("Gaussian.jl")

@doc """
A class representing a set of Tracks
which define the parameters of gaussians
""" ->
type MogTrack
	# Number of points in each track
	num_points::Int64

	# Tracks which determine gaussian parameters
	tracks::Vector{Track}

	# Shift and scale tracks
	mu::Vector{Float64}
	sig::Vector{Float64}
end

@doc """
Constructor, takes Track parameters as input
""" ->
function MogTrack(num_points::Int64)
	return MogTrack(num_points, fill(Track(num_points), 6), zeros(6), ones(6))
end

@doc """
Generate a MogTrack from the prior.
""" ->
function from_prior!(mogtrack::MogTrack)
	for(i in 1:size(mogtrack.tracks)[1])
		calculate_C!(mogtrack.tracks[i])
		generate_y!(mogtrack.tracks[i])
	end
	return nothing
end

@doc """
Evaluate a MogTrack at position (x, y)
""" ->
function evaluate(mogtrack::MogTrack, x::Float64, y::Float64)
	f = 0.0
	params = Array(Float64, (6, ))

	# Loop over the tracks
	for(i in 1:mogtrack.num_points)
		for(j in 1:size(mogtrack.tracks)[1])
			params[j] = mogtrack.mu[j] + mogtrack.sig[j]*mogtrack.tracks[j].y[i]
		end

		params[1] = exp(params[1])
		params[2] = exp(params[2])/(1.0 + exp(params[2]))
		params[5] = exp(params[5])

		gaussian::Gaussian = Gaussian(params)
		f += evaluate(gaussian, x, y)
	end

	return f
end

@doc """
Do a Metropolis proposal.
""" ->
function perturb!(mogtrack::MogTrack)
	which = rand(1:6)
	logH = perturb!(mogtrack.tracks[6])
	return logH
end

ext: jl

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
@testset "random" begin

    abm = @agent(0); s = SimulationFree(abm)

    @test begin
        m = AgentBasedModels.Program_(abm,s)
        r = :(Normal(2,3) += Gamma(3,4))
        r = AgentBasedModels.randomAdapt_(m, r, "cpu")

        r == :(rand(Normal(2,3)) += rand(Gamma(3,4)))
    end

    @test begin
        m = AgentBasedModels.Program_(abm,s)
        r = :(Normal(2,3) += Uniform(3,4))
        r = AgentBasedModels.randomAdapt_(m, r, "gpu")

        r == :(AgentBasedModels.NormalCUDA(rand(),2,3) += AgentBasedModels.UniformCUDA(rand(),3,4))
    end

    @test_throws ErrorException begin
        m = AgentBasedModels.Program_(abm,s)
        r = :(Normal(2,3) += Gamma(3,4))
        r = AgentBasedModels.randomAdapt_(m, r, "gpu")

        r == :(rand(Normal(2,3)) += rand(Gamma(3,4)))
    end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
"""
    Nucleotide_Essentials.FastqtoFasta

Converts a FastqRecord to a FastaRecord. Can also input and convert a .fastq file to a FastaRecord in the same function.

FastqtoFasta(Fastq::Union{String, FastqRecord})
.fastq file => FastqtoFasta(Fastq) => FastaRecord(ID, sequence, filename)
FastqRecord(ID, sequence, quality, filename) => FastqtoFasta(Fastq) => FastaRecord(ID, sequence, filename)

Supported keyword arguments include:

* ```Fastq::Union{String, FastqRecord}```:
    * The full or relative path to a .fastq file
    * A FastqRecord

# Example:

```julia
# Supply the path to a .fastq file that you would like to convert to a FastRecord
myfasta = FastqtoFasta("myfastq.fastq");

# Alternatively, a FastqRecord can be used as the input
myfasta = FastqtoFasta(myFastqRecord);
```
"""
function FastqtoFasta(Fastq::Union{String, FastqRecord})
    if typeof(Fastq) == String
        input = readFastq(Fastq);
        newfasta = FastaRecord(input.ID, input.sequence, input.filename);
    else
        newfasta = FastaRecord(Fastq.ID, Fastq.sequence, Fastq.filename);
    end
    return(newfasta)
end


text:
struct ProximalGradient{T} <: Optimizer
    backtracking::Bool
    extrapolation::T
end
function ProximalGradient(; backtracking = true, extrapolation = VanillaProxGrad())
    return ProximalGradient(backtracking, extrapolation)
end

Base.summary(o::ProximalGradient) = string(
    "Proximal Gradient ",
    o.backtracking && "(bt)",
    summary(o.extrapolation)
)

mutable struct ProximalGradientState{Tx,Te} <: OptimizerState
    it::Int64
    x::Tx
    x_old::Tx
    M::Manifold
    M_old::Manifold
    f_x::Float64
    g_x::Float64
    ∇f_x::Tx
    temp::Tx
    γ::Float64
    extrapolation_state::Te
    ncalls_f::Int64
    ncalls_g::Int64
    ncalls_∇f::Int64
    ncalls_proxg::Int64
end
function ProximalGradientState(
    o::ProximalGradient,
    x::Tx,
    g::R;
    γ = 1e5,
    extra_state::Te = extrapolation_state(o.extrapolation, x, g),
) where {Tx,Te,R}
    return ProximalGradientState{Tx,Te}(
        0,
        copy(x),
        copy(x),
        wholespace_manifold(g, x),
        wholespace_manifold(g, x),
        0.0,
        0.0,
        zero(x),
        zero(x),
        γ,
        extra_state,
        0, 0, 0, 0
    )
end

function print_header(pg::ProximalGradient)
    println("---------------------")
    println("--- Proximal Gradient")
    println(" - backtracking:       ", pg.backtracking)
    println(" - extrapolation:       ", pg.extrapolation)
    return
end

function update_fg∇f!(state::ProximalGradientState, pb)
    state.f_x = f(pb, state.x)
    state.g_x = g(pb, state.x)
    ∇f!(pb, state.∇f_x, state.x)
    state.it += 1
    state.ncalls_∇f += 1
    state.ncalls_f += 1
    state.ncalls_g += 1
    return
end

initial_state(o::ProximalGradient, x, reg; manifold=nothing) = ProximalGradientState(o, x, reg)

function update_iterate!(state::ProximalGradientState, pb, optimizer::ProximalGradient)

    ## Run backtracking line search to update estimate of gradient Lipschitz constant
    if optimizer.backtracking
        state.γ = backtrack_f_lipschitzgradient!(state, pb, state.γ)
    end

    extrapolation!(optimizer.extrapolation, state, state.extrapolation_state, pb)
    return
end

function build_optimstate(::ProximalGradient, state, iteration, time, normstep, minsubgradient_tan, minsubgradient_norm, optimstate_extensions)
    return OptimizationState(
        it = iteration,
        time = time,
        f_x = state.f_x,
        g_x = state.g_x,
        norm_step = normstep,
        ncalls_f = state.ncalls_f,
        ncalls_g = state.ncalls_g,
        ncalls_∇f = state.ncalls_∇f,
        ncalls_proxg = state.ncalls_proxg,
        # ncalls_gradₘF = 0,
        # ncalls_HessₘF = 0,
        norm_minsubgradient_tangent = minsubgradient_tan,
        norm_minsubgradient_normal = minsubgradient_norm,
        additionalinfo = (;
            zip(
                [osextension.key for osextension in optimstate_extensions],
                [
                    copy(osextension.getvalue(state))
                    for osextension in optimstate_extensions
                ],
            )...
        ),
    )
end

display_logs_header_post(::ProximalGradient) = print("γ")

function display_logs_post(::ProximalGradient, state, pb)
    @printf "%.3e" state.γ
    return
end

function backtrack_f_lipschitzgradient!(state, pb, γ)
    τ = 1.2

    ∇f_norm2 = norm(state.∇f_x)^2

    itmax = 200

    it_ls = 0
    while it_ls <= itmax
        state.temp .= state.x .- γ .* state.∇f_x

        # f(pb, state.temp) ≤ state.f_x - 1/(2*γ) * norm(state.temp-state.x)^2 && break
        state.ncalls_f += 1
        f(pb, state.temp) ≤ state.f_x - γ / 2 * ∇f_norm2 && break

        γ = γ / τ
        it_ls += 1
    end

    if it_ls > itmax
        @warn "Gradient backtracking: reached iterations limits."
    end

    return γ
end


text:
# synchronization

export synchronize_threads

"""
     synchronize_threads(state)

in CUDA terms `__synchronize`
in OpenCL terms: `barrier(CLK_LOCAL_MEM_FENCE)`
"""
function synchronize_threads(state)
    error("Not implemented") # COV_EXCL_LINE
end


text:
using MLJLinearModels
using BenchmarkTools, Random, LinearAlgebra
DO_COMPARISONS = false; include("../testutils.jl")

n, p = 50_000, 500
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)

# =============== #
# RIDGE FUNCTIONS #
#
# Hv!           ✅
# =============== #

# No fit_intercept
R.allocate(n, p)
ridge = RidgeRegression(0.5; fit_intercept=false)
Hv! = R.Hv!(ridge, X, y)
v   = randn(p)
Hv  = similar(v)

# Sept 13, 2019 :: 15.300 ms (0 allocations: 0 bytes)
@btime Hv!($Hv, $v)

# With fit_intercept
R.allocate(n, p+1)
ridge = RidgeRegression(λ)
Hv! = R.Hv!(ridge, X, y)
v  = randn(p+1)
Hv = similar(v)

# Sept 13, 2019 :: 26.246 ms (5 allocations: 4.22 KiB)
# (decrease in perf due to views, and need a sum on an array)
@btime Hv!($Hv, $v)

# =============== #
# LASSO FUNCTIONS #
#
# smooth_fg!
# =============== #

# No fit_intercept
R.allocate(n, p)

lasso = LassoRegression(0.5, fit_intercept=false)
smooth_fg! = R.smooth_fg!(lasso, X, y)
v = randn(p)
g = similar(v)

# Sept 13, 2019 :: 15.508 ms (3 allocations: 390.72 KiB)
@btime smooth_fg!($g, $v);

@btime R.get_residuals!($X, $θ, $y)

r = R.get_residuals!(X, θ, y)

@btime R.apply_Xt!(g, X, r)

@btime (g .+= 0.5 .* θ)

# =======
=======
((X, _, _), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.1)

λ  = 50
lr = LassoRegression(λ)
J  = objective(lr, X, y1)
>>>>>>> master

fista = FISTA()
ista = ISTA()

@btime θ_fista = fit($lr, $X, $y1, solver=$fista)

# <old> initial: ~650ms ; 1725 allocations;  104.24 MiB
# using cache: ; ~650ms ; 1666 allocations;   95.02 MiB

@btime θ_ista  = fit($lr, $X, $y1, solver=$ista)

# using cache ;  ~590ms ; 1584 allocations;   91.36 MiB

rr = RidgeRegression(λ)
θ_ridge = fit(rr, X, y1)
cg = CG()

@btime θ_ridge_cg = fit($rr, $X, $y1, solver=$cg)

# <old> initial: ~ 220ms; 100 allocations; 3.5 MiB
# using cache:   ~ 216ms;  76 allocations; 458.5KiB


text:
struct TrialFancyOutputPre end

function prep_trial_fancy_data(t::Trial)
    length(t) > 1 || return nothing

    perm = sortperm(t.times)
    times = t.times[perm]
    gctimes = t.gctimes[perm]

    med = median(t)
    avg = mean(t)
    std = Statistics.std(t)
    min = minimum(t)
    max = maximum(t)

    medtime, medgc = prettytime(time(med)), prettypercent(gcratio(med))
    avgtime, avggc = prettytime(time(avg)), prettypercent(gcratio(avg))
    stdtime, stdgc = prettytime(time(std)), prettypercent(Statistics.std(gctimes ./ times))
    mintime, mingc = prettytime(time(min)), prettypercent(gcratio(min))
    maxtime, maxgc = prettytime(time(max)), prettypercent(gcratio(max))

    memorystr = string(prettymemory(memory(min)))
    allocsstr = string(iallocs(min))

    return (perm, times, gctimes,
            med, avg, std, min, max,
            medtime, medgc,
            avgtime, avggc,
            stdtime, stdgc,
            mintime, mingc,
            maxtime, maxgc,
            memorystr, allocsstr)
end

function show_output(::TrialFancyOutputPre, io, ::MIME"text/plain", t::Trial)
    pad = get(io, :pad, "")
    print(io, "BenchmarkExt.Trial: ", length(t), " sample", if length(t) > 1 "s" else "" end,
          " with ", t.params.evals, " evaluation", if t.params.evals > 1 "s" else "" end ,".\n")

    if length(t) > 1
       (perm, times, gctimes,
            med, avg, std, min, max,
            medtime, medgc,
            avgtime, avggc,
            stdtime, stdgc,
            mintime, mingc,
            maxtime, maxgc,
            memorystr, allocsstr) = prep_trial_fancy_data(t)
    elseif length(t) == 1
        perm = sortperm(t.times)
        times = t.times[perm]
        gctimes = t.gctimes[perm]

        print(io, pad, " Single result which took ")
        printstyled(io, prettytime(times[1]); color=:blue)
        print(io, " (", prettypercent(gctimes[1]/times[1]), " GC) ")
        print(io, "to evaluate,\n")
        print(io, pad, " with a memory estimate of ")
        printstyled(io, prettymemory(t.memory[1]); color=:yellow)
        print(io, ", over ")
        printstyled(io, t.allocs[1]; color=:yellow)
        print(io, " allocations.")
        return
    else
        print(io, pad, " No results.")
        return
    end

    lmaxtimewidth = maximum(length.((medtime, avgtime, mintime)))
    rmaxtimewidth = maximum(length.((stdtime, maxtime)))
    lmaxgcwidth = maximum(length.((medgc, avggc, mingc)))
    rmaxgcwidth = maximum(length.((stdgc, maxgc)))

    # Main stats

    print(io, pad, " Range ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "min"; color=:cyan, bold=true)
    print(io, " … ")
    printstyled(io, "max"; color=:magenta)
    printstyled(io, "):  "; color=:light_black)
    printstyled(io, lpad(mintime, lmaxtimewidth); color=:cyan, bold=true)
    print(io, " … ")
    printstyled(io, lpad(maxtime, rmaxtimewidth); color=:magenta)
    print(io, "  ")
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "min … max")
    printstyled(io, "): "; color=:light_black)
    print(io, lpad(mingc, lmaxgcwidth), " … ", lpad(maxgc, rmaxgcwidth))

    print(io, "\n", pad, " Time  ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "median"; color=:blue, bold=true)
    printstyled(io, "):     "; color=:light_black)
    printstyled(io, lpad(medtime, lmaxtimewidth), rpad(" ", rmaxtimewidth + 5); color=:blue, bold=true)
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "median")
    printstyled(io, "):    "; color=:light_black)
    print(io, lpad(medgc, lmaxgcwidth))

    print(io, "\n", pad, " Time  ")
    printstyled(io, "("; color=:light_black)
    printstyled(io, "mean"; color=:green, bold=true)
    print(io, " ± ")
    printstyled(io, "σ"; color=:green)
    printstyled(io, "):   "; color=:light_black)
    printstyled(io, lpad(avgtime, lmaxtimewidth); color=:green, bold=true)
    print(io, " ± ")
    printstyled(io, lpad(stdtime, rmaxtimewidth); color=:green)
    print(io, "  ")
    printstyled(io, "┊"; color=:light_black)
    print(io, " GC ")
    printstyled(io, "("; color=:light_black)
    print(io, "mean ± σ")
    printstyled(io, "):  "; color=:light_black)
    print(io, lpad(avggc, lmaxgcwidth), " ± ", lpad(stdgc, rmaxgcwidth))

end

struct TrialFancyHistogram end

function show_output(::TrialFancyHistogram, io, ::MIME"text/plain", t::Trial)
    length(t) <= 1 && return
    pad = get(io, :pad, "")
    (perm, times, gctimes,
     med, avg, std, min, max,
     medtime, medgc,
     avgtime, avggc,
     stdtime, stdgc,
     mintime, mingc,
     maxtime, maxgc,
     memorystr, allocsstr) = prep_trial_fancy_data(t)

    lmaxtimewidth = maximum(length.((medtime, avgtime, mintime)))
    rmaxtimewidth = maximum(length.((stdtime, maxtime)))
    lmaxgcwidth = maximum(length.((medgc, avggc, mingc)))
    rmaxgcwidth = maximum(length.((stdgc, maxgc)))

    histquantile = 0.99
    # The height and width of the printed histogram in characters.
    histheight = 2
    histwidth = 42 + lmaxtimewidth + rmaxtimewidth

    histtimes = times[1:round(Int, histquantile*end)]
    histmin = get(io, :histmin, first(histtimes))
    histmax = get(io, :histmax, last(histtimes))
    logbins = get(io, :logbins, nothing)
    bins = bindata(histtimes, histwidth - 1, histmin, histmax)
    append!(bins, [1, floor((1-histquantile) * length(times))])
    # if median size of (bins with >10% average data/bin) is less than 5% of max bin size, log the bin sizes
    if logbins === true || (logbins === nothing && median(filter(b -> b > 0.1 * length(times) / histwidth, bins)) / maximum(bins) < 0.05)
        bins, logbins = log.(1 .+ bins), true
    else
        logbins = false
    end
    hist = asciihist(bins, histheight)
    hist[:,end-1] .= ' '
    maxbin = maximum(bins)

    delta1 = (histmax - histmin) / (histwidth - 1)
    if delta1 > 0
        medpos = 1 + round(Int, (histtimes[length(times) ÷ 2] - histmin) / delta1)
        avgpos = 1 + round(Int, (mean(times) - histmin) / delta1)
    else
        medpos, avgpos = 1, 1
    end

    print(io, "\n")
    for r in axes(hist, 1)
        print(io, "\n", pad, "  ")
        for (i, bar) in enumerate(view(hist, r, :))
            color = :default
            if i == avgpos color = :green end
            if i == medpos color = :blue end
            printstyled(io, bar; color=color)
        end
    end

    remtrailingzeros(timestr) = replace(timestr, r"\.?0+ " => " ")
    minhisttime, maxhisttime = remtrailingzeros.(prettytime.(round.([histmin, histmax], sigdigits=3)))

    print(io, "\n", pad, "  ", minhisttime)
    caption = "Histogram: " * ( logbins ? "log(frequency)" : "frequency" ) * " by time"
    if logbins
        printstyled(io, " " ^ ((histwidth - length(caption)) ÷ 2 - length(minhisttime)); color=:light_black)
        printstyled(io, "Histogram: "; color=:light_black)
        printstyled(io, "log("; bold=true, color=:light_black)
        printstyled(io, "frequency"; color=:light_black)
        printstyled(io, ")"; bold=true, color=:light_black)
        printstyled(io, " by time"; color=:light_black)
    else
        printstyled(io, " " ^ ((histwidth - length(caption)) ÷ 2 - length(minhisttime)), caption; color=:light_black)
    end
    print(io, lpad(maxhisttime, ceil(Int, (histwidth - length(caption)) / 2) - 1), " ")
    printstyled(io, "<"; bold=true)
end

struct TrialFancyOutputPost end

function show_output(::TrialFancyOutputPost, io, ::MIME"text/plain", t::Trial)
    # Memory info
    length(t) <= 1 && return
    pad = get(io, :pad, "")
    (perm, times, gctimes,
     med, avg, std, min, max,
     medtime, medgc,
     avgtime, avggc,
     stdtime, stdgc,
     mintime, mingc,
     maxtime, maxgc,
     memorystr, allocsstr) = prep_trial_fancy_data(t)

    print(io, "\n\n", pad, " Memory estimate")
    printstyled(io, ": "; color=:light_black)
    printstyled(io, memorystr; color=:yellow)
    print(io, ", allocs estimate")
    printstyled(io, ": "; color=:light_black)
    printstyled(io, allocsstr; color=:yellow)
    print(io, ".")
end


text:
using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, SpecialFunctions, FillArrays, Plots
import ApproxFun: blockbandwidths, Vec, PiecewiseSegment
import MultivariateOrthogonalPolynomials: DirichletTriangle

#######
# Example 1 Poisson
#######

S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
Δ = Laplacian(S)

N = 500
M = sparse(Δ[Block.(1:N), Block.(1:N)])
F = lu(M)
u_1 = F \ pad(coefficients(Fun(1, rangespace(Δ))), size(M,1))
u_2 = F \ pad(coefficients(Fun((x,y) -> x^2 + y^2, rangespace(Δ))), size(M,1))
u_3 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y), rangespace(Δ))), size(M,1))
u_4 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y) *cos(300x*y), rangespace(Δ))), size(M,1))
u_5 = F \ pad(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(Δ))), size(M,1))

gr()

plot(abs.(u_1)[1:10_000]; yscale=:log10, xscale=:log10, label="1")
plot!(abs.(u_2)[1:10_000]; yscale=:log10, xscale=:log10, label="x^2 + y^2")
plot!(abs.(u_3)[1:10_000]; yscale=:log10, xscale=:log10, label="xy(1-x-y)")
plot!(abs.(u_4)[1:10_000]; yscale=:log10, xscale=:log10, label="xy(1-x-y)cos(300xy)")

#######
# Example 2 Variable coefficient Helmholts
######

S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
x,y = Fun(Triangle())

V = 1-(3*(x-1)^4 + 5*y^4)

f = Fun((x,y) -> x*y*exp(x) , JacobiTriangle(1.,1.,1.))

k = 100
L = Laplacian(S) + (k^2 * V)
N = ceil(Int,2k)
@time M̃ = L[Block.(1:N), Block.(1:N)]
@time M = sparse(M̃)
@time F = lu(M)
@time u_c = F \ pad(coefficients(f), size(M,1))
u = Fun(domainspace(L), u_c)

#######
# Example 3 Bi-harmonic
#######

S = TriangleWeight(2,2,2, JacobiTriangle(2,2,2))
Δ = Laplacian(S)
L = Δ^2

N = 1001
M = sparse(L[Block.(1:N), Block.(1:N)])
F = lu(M)

u_1 = F \ pad(coefficients(Fun(1, rangespace(L))), size(M,1))
u_2 = F \ pad(coefficients(Fun((x,y) -> x^2 + y^2, rangespace(L))), size(M,1))
u_3 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y), rangespace(L))), size(M,1))
# u_4 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y) *cos(300x*y), rangespace(L))), size(M,1))
u_4 = F \ pad(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(L))), size(M,1))

bnrms = u -> [norm(PseudoBlockArray(u, 1:N)[Block(K)],Inf) for K=1:N]

plot(bnrms(u_1); xscale=:log10, yscale=:log10, linewidth=3, label="1")
plot!(bnrms(u_2); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")
plot!(bnrms(u_3); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")
plot!(bnrms(u_4); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")

plot(abs.(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(L)))); scale=:log10, yscale=:log10)

#######
# Example 4 Laplace's equation
#######

S = DirichletTriangle{1,1,1}()
B = Dirichlet() : S
Δ = Laplacian() : S
# L = [B; Δ]

fs = ((x,y) -> exp(x)*cos(y),
      (x,y) -> x^2,
      (x,y) ->  x^3*(1-x)^3*(1-y)^3,
      (x,y) -> x^3*(1-x)^3*(1-y)^3 * cos(100x*y))

N = 1000;
@time BM = sparse(B[Block.(1:N), Block.(1:N)]);
@time ΔM = sparse(Δ[Block.(1:N), Block.(1:N)]);
@time M = [BM; ΔM];
M = [sparse(I, size(M,1), 2) M] # add tau term
@time F = factorize(M);

# RHS
k = 1
∂u = Fun(fs[k], rangespace(B))
r = pad(coefficients(∂u, rangespace(B)), size(Fs[N],1))
@time u_cfs = F \ r

#######
# Example 5 Transport equation
#######

S = DirichletTriangle{0,1,0}()
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0)))
B = I : S → ∂S

Dx = Derivative([1,0]) : S
Dt = Derivative([0,1]) : S

x,_ = Fun(∂S)
u₀ = x*(1-x)*exp(x)

# u_t = u_x
L = [B; Dx - Dt]
N = 20; M = sparse(L[Block.(1:N), Block.(1:N)]);
u_c = M \ pad(coefficients([u₀; 0], rangespace(L)), size(M,1))
u = Fun(domainspace(L), u_c)

# 0.5u_t = u_x
S = DirichletTriangle{0,1,1}()
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0))) , Legendre(Segment(Vec(0.0,1.0),Vec(1.0,0)))
B_x = I : S → ∂S[1]
B_z = I : S → ∂S[2]
Dx = Derivative([1,0]) : S
Dt = Derivative([0,1]) : S

L = [B_x; B_z; Dx - 0.5Dt]
N = 100; M = sparse(L[Block.(1:N), Block.(1:N)])
u₀_x = Fun((x,_) -> x*exp(x-1), rangespace(B_x))
u₀_z = Fun((x,y) -> (1-y), rangespace(B_z))
u_c = M \ pad(coefficients([u₀_x; u₀_z; 0], rangespace(L)), size(M,1))
u = Fun(domainspace(L), u_c)

# u_t + u_x = 0
S = DirichletTriangle{1,1,0}()
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0))) , Legendre(Segment(Vec(0.0,0.0),Vec(0.0,1.0)))
B_x = I : S → ∂S[1]
B_y = I : S → ∂S[2]
Dx = Derivative([1,0]) : S
Dt = Derivative([0,1]) : S

L = [B_x; B_y; Dt + Dx]
N = 100; M = sparse(L[Block.(1:N), Block.(1:N)])
u₀_x = Fun((x,_) -> (1-x)*exp(x), rangespace(B_x))
u₀_y = Fun((x,y) -> (1-y), rangespace(B_y))

u_c = M \ pad(coefficients([u₀_x; u₀_y; 0], rangespace(L)), size(M,1))
u = Fun(domainspace(L), u_c)

#######
# Example 7 Helmholtz in a polygon
#######
d = Triangle() ∪ Triangle(Vec(1,1),Vec(1,0),Vec(0,1)) ∪
                 Triangle(Vec(1,1),Vec(0,1),Vec(0,2)) ∪
                 Triangle(Vec(0,0),Vec(0,1),Vec(-1,1.5))

f = Fun((x,y) -> cos(x*sin(y)), d)

∂d = PiecewiseSegment([Vec(0.,0), Vec(1.,0), Vec(1,1), Vec(0,2), Vec(0,1), Vec(-1,1.5), Vec(0,0)])

S = DirichletTriangle{1,1,1}.(components(d))
∂S = Legendre.(components(∂d))

∂S12 = Legendre(Segment(Vec(1.,0),Vec(0,1)))
∂S23 = Legendre(Segment(Vec(1.,1),Vec(0,1)))
∂S14 = Legendre(Segment(Vec(0.,0),Vec(0,1)))

Dx = Derivative([1,0])
Dy = Derivative([0,1])

k = 20

L = [(I : S[1] → ∂S[1]) 0                 0                     0                  0                0                       0                    0                  0                   0                    0                  0;
     0                  0                 0                    (I : S[2] → ∂S[2])  0                0                       0                    0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       (I : S[3] → ∂S[3])   0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       (I : S[3] → ∂S[4])   0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       0                    0                  0                   (I : S[4] → ∂S[5])   0                  0;
     0                  0                 0                     0                  0                0                       0                    0                  0                   (I : S[4] → ∂S[6])   0                  0;
     (I : S[1] → ∂S12)  0                 0                   -(I : S[2] → ∂S12)   0                0                       0                    0                  0                    0                   0                  0;
     0                  (I : S[1] → ∂S12) (I : S[1] → ∂S12)     0                -(I : S[2] → ∂S12) -(I : S[2] → ∂S12)      0                    0                  0                    0                   0                  0;
     0                  0                 0                    (I : S[2] → ∂S23)   0                0                      -(I : S[3] → ∂S23)    0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                (I : S[2] → ∂S23)       0                    0                  -(I : S[3] → ∂S23)   0                   0                  0;
     (I : S[1] → ∂S14)  0                 0                     0                  0                0                       0                    0                  0                   -(I : S[4] → ∂S14)   0                  0;
     0                  (I : S[1] → ∂S14) 0                     0                  0                0                       0                    0                  0                    0                  -(I : S[4] → ∂S14)   0;
     Dx                -I                 0                     0                  0                0                       0                    0                  0                    0                   0                  0;
     Dy                 0                -I                     0                  0                0                       0                    0                  0                    0                   0                  0;
     0                  0                 0                     Dx                -I                0                       0                    0                  0                    0                   0                  0;
     0                  0                 0                     Dy                 0               -I                       0                    0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       (Dx : S[3])         -(I : S[3])         0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       Dy                   0                 -I                    0                   0                  0;
     0                  0                 0                     0                  0                0                       0                    0                  0                    Dx                 -I                  0;
     0                  0                 0                     0                  0                0                       0                    0                  0                    (Dy : S[4])         0                 -(I : S[4])
    k^2*I              Dx                Dy                    0                  0                0                       0                    0                  0                    0                   0                  0;
     0                  0                 0                     k^2*I                  Dx               Dy                      0                    0                  0                    0                   0                  0;
     0                  0                 0                     0                  0                0                       k^2*I                    Dx                 Dy                   0                   0                  0;
     0                  0                 0                     0                  0                0                       0                    0                  0                    k^2*I                  Dx                 Dy
     ]

f = ones(∂d)

N = 101
@time M = sparse(L[Block.(1:N), Block.(1:N)]);
@time F = qr(M);
@time u_c = F \ pad(coefficients(vcat(components(f)..., Zeros(18)...), rangespace(L)), size(M,1));
u = Fun(domainspace(L), u_c)


text:
using Tmatrix
using Zygote
using Plots
function objective_function(r_array_input, θ_array_input)
    # println(r_array)
    # println("I am <<objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")

    # quadrupling_mesh
    r_θ_array = Tmatrix.quadruple_mesh_density(r_array_input, θ_array_input)
    r_array = r_θ_array[:, 1]
    θ_array = r_θ_array[:, 2]
    println(r_array)
    println(θ_array)
    ϕ_array = zeros(size(θ_array))

    T = Tmatrix.T_matrix_SeparateRealImag_arbitrary_mesh(
        n_max,
        wl_or_freq_input,
        input_unit,
        Eps_r_r_1,
        Eps_r_i_1,
        Mu_r_r_1,
        Mu_r_i_1,
        Eps_r_r_2,
        Eps_r_i_2,
        Mu_r_r_2,
        Mu_r_i_2,
        (collect(r_array)),
        θ_array,
        ϕ_array,
        rotationally_symmetric,
        symmetric_about_plane_perpendicular_z,
        BigFloat_precision,
    )
    k1_complex = Tmatrix.get_WaveVector(
        wl_or_freq_input;
        input_unit = input_unit,
        Eps_r = Complex(Eps_r_r_1, Eps_r_i_1),
        Mu_r = Complex(Mu_r_r_1, Mu_r_i_1),
    )
    # return abs(T[1]) # if we need to minimize the first element of T-matrix

    # return sum(abs.(T)) # if we need to minimize the sum of all real and imag

    return Tmatrix.get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, k1_complex) # if we need to minimize scattering cross section

    """
    return -1 * Tmatrix.get_OrentationAv_emissivity_from_Tmatrix(
        T,
        Complex(k1_r, k1_i),
        Tmatrix.calculate_surface_area_of_axisymmetric_particle(r_array, θ_array)
    ) # if we need to minimize emissivity
    """
end

function ∂objective_function(r_array, θ_array)
    # println("I am <<∂objective_function>>, and the type of r_array is $(typeof(r_array))")
    # println("r_array = $r_array")
    return Zygote.gradient(
        r_array -> objective_function(r_array, θ_array),
        (collect(r_array)),
    )[1]
end

wl_or_freq_input = 1e-6;
input_unit = "m";
n_max = 1
Eps_r_r_1 = 1.0;
Eps_r_i_1 = 0.0;
Eps_r_r_2 = 1.5;
Eps_r_i_2 = 0.01;
Mu_r_r_1 = 1.0;
Mu_r_i_1 = 0.0;
Mu_r_r_2 = 1.0;
Mu_r_i_2 = 0.0;

n_θ_points = 10
θ_array = collect(LinRange(1e-6, pi - 1e-6, n_θ_points))
r_array = 0.5 * wl_or_freq_input * ones(size(θ_array))
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
BigFloat_precision = nothing

learning_rate = 0.5e-4
global loss_here = -1e6
global ∂loss = 0
println("Starting optimization for T-matrix ...")
loss_array = []
loss_old = 9999999999999
for n_iteration in 1:50
    global r_array = r_array .- learning_rate .* ∂loss

    global loss_here = objective_function(r_array, θ_array)
    global ∂loss = ∂objective_function(r_array, θ_array)
    append!(loss_array, loss_here)

    println()
    println("iteration #$n_iteration: loss_here = $loss_here, ∂loss = $∂loss")
    println("r_array = $r_array")

    xyz = vcat(
        Tmatrix.convert_coordinates_Sph2Cart.(r_array, θ_array, zeros(size(r_array)))...,
    )
    p1 = plot!(
        xyz[:, 1],
        xyz[:, 3],
        aspect_ratio = :equal,
        label = "particle after $n_iteration iterations",
        title = "scattering cross section = $loss_here m",
    )
    p2 = plot(
        1:length(loss_array),
        loss_array,
        xlabel = "iteration #",
        ylabel = "scattering cross section (m)",
    )
    fig = plot(p1, p2, layout = (1, 2), size = (1200, 800))
    mkpath("cache/iteration_particle_plots/maximizing_emissivity")
    savefig(
        fig,
        "cache/iteration_particle_plots/maximizing_emissivity/particle_geom_iteration_$(n_iteration).png",
    )
    # if (loss_here - loss_old) < 1e-6; break; end
    global loss_old = loss_here
end


text:
# Autogenerated wrapper script for DynarePreprocessor_jll for powerpc64le-linux-gnu-cxx11
export dynare_preprocessor

JLLWrappers.@generate_wrapper_header("DynarePreprocessor")
JLLWrappers.@declare_executable_product(dynare_preprocessor)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        dynare_preprocessor,
        "bin/dynare-preprocessor",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function opposite_relationship(obj::Relationship)
    return jcall(obj, "oppositeRelationship", Relationship, ())
end

function to_string(obj::Relationship)
    return jcall(obj, "toString", JString, ())
end

function value_of(::Type{Relationship}, arg0::JString)
    return jcall(Relationship, "valueOf", Relationship, (JString,), arg0)
end

function values(::Type{Relationship})
    return jcall(Relationship, "values", Vector{Relationship}, ())
end


text:
my_f(x,y) = 2x + 3y
derivative_of_my_f(x,y) = ForwardDiff.derivative(x -> my_f(x,y), x)


text:
"""
    feynman_kac(𝕋 [; t, f, ψ, v])

With direction = :backward
Solve the following PDE:
u(x, t[end]) = ψ(x)
0 = u_t + 𝕋u - v(x, t)u + f(x, t)

With direction = :forward
Solve the following PDE:
u(x, t[1]) = ψ(x)
u_t = 𝕋u - v(x, t)u + f(x, t)
"""
function feynman_kac(𝕋;
    t::AbstractVector = range(0, 100, step = 1/12),
    f::Union{AbstractVector, AbstractMatrix} = zeros(size(𝕋, 1)),
    ψ::AbstractVector = ones(size(𝕋, 1)),
    v::Union{AbstractVector, AbstractMatrix} = zeros(size(𝕋, 1)),
    direction= :backward)
    if direction == :backward
        u = zeros(size(𝕋, 1), length(t))
        u[:, end] = ψ
        if isa(f, AbstractVector) && isa(v, AbstractVector)
            if isa(t, AbstractRange)
                dt = step(t)
                B = factorize(I + (Diagonal(v) - 𝕋) * dt)
                for i in (length(t)-1):(-1):1
                    ψ = ldiv!(B, u[:, i+1] .+ f .* dt)
                    u[:, i] = ψ
                end
            else
                for i in (length(t)-1):(-1):1
                    dt = t[i+1] - t[i]
                    B = I + (Diagonal(v) - 𝕋) * dt
                    u[:, i] = B \ (u[:, i+1] .+ f .* dt)
                end
            end
        elseif isa(f, AbstractMatrix) && isa(v, AbstractMatrix)
            for i in (length(t)-1):(-1):1
                dt = t[i+1] - t[i]
                B = I + (Diagonal(view(v, :, i)) - 𝕋) * dt
                u[:, i] = B \ (u[:, i+1] .+ f[:, i] .* dt)
            end
        else
            error("f and v must be both AbstractVectors or both AbstractMatrices")
        end
        return u
    elseif direction == :forward
        u = feynman_kac(𝕋; t = - reverse(t), ψ = ψ, f = f, v = v, direction = :backward)
        return u[:,end:-1:1]
    else
        error("Direction must be :backward or :forward")
    end
end


text:
function add_constraints!(
   mn_data,
   ref_subnet,
   subnet_idx,
   model,
   va, vm, pg, qg, p, q, p_dc, q_dc, g, b, b_fr, b_to, g_fr, g_to, l_series,
   c_inv_series, f, b_shunt, cost_pg, cost_dcline, area_pg, zone_pg, redispatch_pg,
   redispatch_qg, redispatch_vm, redispatch_va, p_i, q_i,
   constraints,
   obj, gen_areas, area_interface,
   gen_zones, zone_interface,
   direct_pq,
   dc_current_limit,
   master_subnet
   )
   # Add Constraints
   # ---------------
   # println("Adding constraints in subnet $subnet_idx.")
   # check if DC subnet
   dc_subnet = ((!ref_subnet[:variable_f]) && (ref_subnet[:f_base] == 0))
   if (:f_max in keys(ref_subnet))
      println("f_max = $(ref_subnet[:f_max])")
      dc_subnet = dc_subnet || (ref_subnet[:f_max] == 0)
   end
   if dc_subnet
      println("Subnetwork $subnet_idx is DC. Using power flow equations for a DC network.")
   end
   # Fix the voltage angle to zero at the reference bus in each subnetwork
   # If no direct PQ control, fix the angle only for the master subnet
   if !dc_subnet && (direct_pq || (subnet_idx==master_subnet))
      for (i,bus) in ref_subnet[:ref_buses]
         if !(subnet_idx in keys(constraints[:theta_ref]))
            constraints[:theta_ref][subnet_idx] = Dict{Int64, Any}()
         end
         constraints[:theta_ref][subnet_idx][i] = @constraint(model, va[subnet_idx][i] == 0)
         # println("set angle reference bus $i in subnet $subnet_idx")
      end
   end

   # Nodal power balance constraints
   # Defining frequency dependent shunts

   # For each subnet, the bus shunts are indexed by tuple (i,j),
   # denoting the jth shunt at bus i.
   b_shunt[subnet_idx] = Dict{Tuple{Int64,Int64},Any}()
   for (i,bus) in ref_subnet[:bus]
      # Define converter interface power flow variables by bus
      if haskey(bus, "converter_index")
         bus_interfaces = bus["converter_index"]
      else
         bus_interfaces = Array{Int64,1}()
      end
      # Build a list of the loads and shunt elements connected to the bus i
      bus_loads = [ref_subnet[:load][l] for l in ref_subnet[:bus_loads][i]]
      bus_shunts = [ref_subnet[:shunt][s] for s in ref_subnet[:bus_shunts][i]]
      if ref_subnet[:variable_f]
         for (j,shunt) in enumerate(bus_shunts)
            if dc_subnet
               b_shunt[subnet_idx][i,j] = 0
            else
               c_bus_shunt = shunt["bs"]/(2*pi*ref_subnet[:f_base])
               b_shunt[subnet_idx][i,j] = @variable(
                  model,
                  base_name = "b_shunt_$(subnet_idx)_$(i)_$(j)",
               )
               @constraint(
                  model,
                  (2pi)*c_bus_shunt*f[subnet_idx] == b_shunt[subnet_idx][i,j]
               )
               # print("shunt $(j) = variable\n")
            end
         end
      else #if not a variable frequency subnetwork
         if :f_fixed in keys(ref_subnet)
            f_ratio = ref_subnet[:f_fixed]/ref_subnet[:f_base]
         else
            f_ratio = 1
         end
         for (j,shunt) in enumerate(bus_shunts)
            b_shunt[subnet_idx][i,j] = shunt["bs"]*f_ratio
            # print("shunt $(j) = $(shunt["bs"])\n")
         end
      end
      # print("defined shunts\n")
      # Active power balance at node i
      @constraint(
         model,
         sum(p[subnet_idx][a] for a in ref_subnet[:bus_arcs][i]) +                  # sum of active power flow on lines from bus i +
         sum(p_dc[subnet_idx][a_dc] for a_dc in ref_subnet[:bus_arcs_dc][i]) ==     # sum of active power flow on HVDC lines from bus i =
         sum(pg[subnet_idx][g] for g in ref_subnet[:bus_gens][i]) +                 # sum of active power generation at bus i +
         sum(p_i[conv_idx[1][1]][(subnet_idx, i)]
             for conv_idx in bus_interfaces) -              # sum of active power from interface to bus i -
         sum(load["pd"] for load in bus_loads) -                                    # sum of active load consumption at bus i -
         sum(shunt["gs"] for shunt in bus_shunts)*vm[subnet_idx][i]^2               # sum of active shunt element injections at bus i
      )
      # println("busgens: $(ref_subnet[:bus_gens][i])")
      # Reactive power balance at node i
      # println("bus_interfaces: $bus_interfaces")
      if !dc_subnet
         if ref_subnet[:variable_f]
            @NLconstraint(
               model,
               sum(q[subnet_idx][a] for a in ref_subnet[:bus_arcs][i]) +                                    # sum of reactive power flow on lines from bus i +
               sum(q_dc[subnet_idx][a_dc] for a_dc in ref_subnet[:bus_arcs_dc][i]) ==                       # sum of reactive power flow on HVDC lines from bus i =
               sum(qg[subnet_idx][g] for g in ref_subnet[:bus_gens][i]) +                                   # sum of reactive power generation at bus i +
               sum(q_i[conv_idx[1][1]][(subnet_idx, i)]
                   for conv_idx in bus_interfaces) -                                # sum of reactive power from interface to bus i -
               sum(load["qd"] for load in bus_loads) -                                                      # sum of reactive load consumption at bus i -
               sum(-b_shunt[subnet_idx][i,j] for (j,shunt) in enumerate(bus_shunts))*vm[subnet_idx][i]^2    # sum of reactive shunt element injections at bus i
            )
         else
            @constraint(
               model,
               sum(q[subnet_idx][a] for a in ref_subnet[:bus_arcs][i]) +                                    # sum of reactive power flow on lines from bus i +
               sum(q_dc[subnet_idx][a_dc] for a_dc in ref_subnet[:bus_arcs_dc][i]) ==                       # sum of reactive power flow on HVDC lines from bus i =
               sum(qg[subnet_idx][g] for g in ref_subnet[:bus_gens][i]) +                                   # sum of reactive power generation at bus i +
               sum(q_i[conv_idx[1][1]][(subnet_idx, i)]
                   for conv_idx in bus_interfaces) -                                # sum of reactive power from interface to bus i -
               sum(load["qd"] for load in bus_loads) -                                                      # sum of reactive load consumption at bus i -
               sum(-b_shunt[subnet_idx][i,j] for (j,shunt) in enumerate(bus_shunts))*vm[subnet_idx][i]^2    # sum of reactive shunt element injections at bus i
            )
         end
      end
   end

   # Branch power flow physics and limit constraints
   g[subnet_idx] = Dict{Int64,Any}()
   b[subnet_idx] = Dict{Int64,Any}()
   b_fr[subnet_idx] = Dict{Int64,Any}()
   b_to[subnet_idx] = Dict{Int64,Any}()
   g_fr[subnet_idx] = Dict{Int64,Any}()
   g_to[subnet_idx] = Dict{Int64,Any}()
   l_series[subnet_idx] = Dict{Int64,Any}()
   c_inv_series[subnet_idx] = Dict{Int64,Any}()
   for (i,branch) in ref_subnet[:branch]
      # Build the from variable id of the i-th branch, which is a tuple given by (branch id, from bus, to bus)
      f_idx = (i, branch["f_bus"], branch["t_bus"])
      # Build the to variable id of the i-th branch, which is a tuple given by (branch id, to bus, from bus)
      t_idx = (i, branch["t_bus"], branch["f_bus"])
      # note: it is necessary to distinguish between the from and to sides of a branch due to power losses

      p_fr = p[subnet_idx][f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
      p_to = p[subnet_idx][t_idx]                     # p_to is a reference to the optimization variable p[t_idx]
      # note: adding constraints to p_fr is equivalent to adding constraints to p[f_idx], and so on

      vm_fr = vm[subnet_idx][branch["f_bus"]]         # vm_fr is a reference to the optimization variable vm on the from side of the branch
      vm_to = vm[subnet_idx][branch["t_bus"]]         # vm_to is a reference to the optimization variable vm on the to side of the branch

      # Compute the branch parameters and transformer ratios from the data
      g_base, b_base = PowerModels.calc_branch_y(branch)
      tr, ti = PowerModels.calc_branch_t(branch)
      g_fr[subnet_idx][i] = branch["g_fr"]
      b_fr_base = branch["b_fr"]
      g_to[subnet_idx][i] = branch["g_to"]
      b_to_base = branch["b_to"]
      tm = branch["tap"]^2
      # note: tap is assumed to be 1.0 on non-transformer branches
      # print("g_base: $(g_base)\n")
      # print("b_base: $(b_base)\n")
      if dc_subnet && branch["shift"] != 0
         throw(ArgumentError("Subnetwork $subnet_idx is specified as a DC network, but branch $(branch["index"]) has a nonzero phase shift: $(branch["shift"])"))
      end

      if ref_subnet[:variable_f]
         f_base = ref_subnet[:f_base]
         # Series parameters

         # Identify series L and C from line data
         # Assumption: the given base line impedance is either
         # purely inductive (positive) or purely capacitive (negative).
         # This means that series compensation should be represented in the
         # model as a separate branch in series
         x_base = branch["br_x"]
         if x_base >= 0
            l_series[subnet_idx][i] = x_base / (2pi*f_base)
            c_inv_series[subnet_idx][i] = 0
         else
            l_series[subnet_idx][i] = 0
            c_inv_series[subnet_idx][i] = - x_base * (2pi*f_base)
         end
         r = branch["br_r"]
         # Shunt parameters
         c_fr = b_fr_base / (2pi*f_base)
         c_to = b_to_base / (2pi*f_base)
         # print("c_fr: $(c_fr)\n")
         # print("c_to: $(c_to)\n")

         if dc_subnet
            g[subnet_idx][i] = 1/r
            b[subnet_idx][i] = 0
            b_fr[subnet_idx][i] = 0
            b_to[subnet_idx][i] = 0
         else
            g[subnet_idx][i] = @variable(
               model,
               base_name = "g_$(subnet_idx)_$i",
            )
            @NLconstraint(
               model,
               g[subnet_idx][i] == r/(r^2 + (
                  l_series[subnet_idx][i]^2 * (2pi)^2*f[subnet_idx]^2
                  + c_inv_series[subnet_idx][i]^2 / ((2pi)^2*f[subnet_idx]^2)
                  - 2 * l_series[subnet_idx][i] * c_inv_series[subnet_idx][i]
                  ))
            )
            # print("defined G and equality constraint\n")
            b[subnet_idx][i] = @variable(
               model,
               base_name = "b_$(subnet_idx)_$i",
            )
            @NLconstraint(
               model,
               b[subnet_idx][i] == -(
                  (
                     2pi*f[subnet_idx]*l_series[subnet_idx][i]
                     - c_inv_series[subnet_idx][i] / ((2pi)*f[subnet_idx])
                     )
                  /(r^2 + (
                     l_series[subnet_idx][i]^2 * (2pi)^2*f[subnet_idx]^2
                     + c_inv_series[subnet_idx][i]^2 / ((2pi)^2*f[subnet_idx]^2)
                     - 2 * l_series[subnet_idx][i] * c_inv_series[subnet_idx][i]
                     )
                  )
               )
            )
            b_fr[subnet_idx][i] = @variable(
               model,
               base_name = "b_fr_$(subnet_idx)_$i",
            )
            @constraint(
               model,
               b_fr[subnet_idx][i] == 2pi*c_fr*f[subnet_idx]
            )
            b_to[subnet_idx][i] = @variable(
               model,
               base_name = "b_to_$(subnet_idx)_$i",
            )
            @constraint(
               model,
               b_to[subnet_idx][i] == 2pi*c_to*f[subnet_idx]
            )
            # print("defined variable G, B\n")
         end
      else #if !ref_subnet[:variable_f]
         # The frequency is fixed at the base value and all remaining line parameters are defined as fixed based on the base values at that frequency
         # This also applies if the fixed base frequency is DC and the parameters have been defined at DC.
         # If the subnetwork is DC but the line parameters are defined at a different base frequency,
         # the frequency limits should be set to 0 and variable_f true, which is handled in the clause above.
         if :f_fixed in keys(ref_subnet)
            f_ratio = ref_subnet[:f_fixed]/ref_subnet[:f_base]
         else
            f_ratio = 1
         end

         if f_ratio == 1
            g[subnet_idx][i] = g_base
            b[subnet_idx][i] = b_base
            b_fr[subnet_idx][i] = b_fr_base
            b_to[subnet_idx][i] = b_to_base
         elseif f_ratio == 0 # DC subnetwork
            r = branch["br_r"]
            g[subnet_idx][i] = 1/r
            b[subnet_idx][i] = 0
            b_fr[subnet_idx][i] = 0
            b_to[subnet_idx][i] = 0
         else
            f_base = ref_subnet[:f_base]
            # Series parameters

            # Identify series L and C from line data
            # Assumption: the given base line impedance is either
            # purely inductive (positive) or purely capacitive (negative).
            # This means that series compensation should be represented in the
            # model as a separate branch in series
            x_base = branch["br_x"]
            if x_base >= 0
               l_series[subnet_idx][i] = x_base / (2pi*f_base)
               c_inv_series[subnet_idx][i] = 0
            else
               l_series[subnet_idx][i] = 0
               c_inv_series[subnet_idx][i] = - x_base * (2pi*f_base)
            end
            r = branch["br_r"]
            # Shunt parameters
            c_fr = b_fr_base / (2pi*f_base)
            c_to = b_to_base / (2pi*f_base)

            g[subnet_idx][i] = r/(r^2 + (
               l_series[subnet_idx][i]^2 * (2pi)^2*ref_subnet[:f_fixed]^2
               + c_inv_series[subnet_idx][i]^2 / ((2pi)^2*ref_subnet[:f_fixed]^2)
               - 2 * l_series[subnet_idx][i] * c_inv_series[subnet_idx][i]
               ))
            b[subnet_idx][i] = -(
               (
                  2pi*ref_subnet[:f_fixed]*l_series[subnet_idx][i]
                  - c_inv_series[subnet_idx][i] / ((2pi)*ref_subnet[:f_fixed])
                  )
               /(r^2 + (
                  l_series[subnet_idx][i]^2 * (2pi)^2*ref_subnet[:f_fixed]^2
                  + c_inv_series[subnet_idx][i]^2 / ((2pi)^2*ref_subnet[:f_fixed]^2)
                  - 2 * l_series[subnet_idx][i] * c_inv_series[subnet_idx][i]
                  )
               )
            )
            b_fr[subnet_idx][i] = 2pi*c_fr*ref_subnet[:f_fixed]
            b_to[subnet_idx][i] = 2pi*c_to*ref_subnet[:f_fixed]
         end

      end

      g_br = g[subnet_idx][i]          # g_br is a reference to the branch series conductance optimization variable or fixed parameter g[subnet_idx][i]
      b_br = b[subnet_idx][i]          # b_br is a reference to the branch series susceptance optimization variable or fixed parameter b[subnet_idx][i]
      b_fr_br = b_fr[subnet_idx][i]    # b_fr_br is a reference to the branch "from" side shunt susceptance optimization variable or fixed parameter b_fr_br[subnet_idx][i]
      b_to_br = b_to[subnet_idx][i]    # b_to_br is a reference to the branch "to" side shunt susceptance optimization variable or fixed parameter b_to_br[subnet_idx][i]
      g_fr_br = g_fr[subnet_idx][i]    # g_fr_br is a reference to the branch "from" side shunt conductance optimization variable or fixed parameter g_fr_br[subnet_idx][i]
      g_to_br = g_to[subnet_idx][i]    # g_to_br is a reference to the branch "to" side shunt conductance optimization variable or fixed parameter g_to_br[subnet_idx][i]

      # print("adding power flow constraints\n")

      if dc_subnet
         # DC power flow
         if !(:k_cond in keys(ref_subnet))
            ref_subnet[:k_cond] = sqrt(3)
         end
         if !(:k_ins in keys(ref_subnet))
            ref_subnet[:k_ins] = 1
         end

         k_cond = ref_subnet[:k_cond]
         k_ins = ref_subnet[:k_ins]

         println("k_cond: $k_cond")
         println("k_ins: $k_ins")
         println("g_br: $g_br")
         ## From side of the branch flow
         @NLconstraint(
            model,
            p_fr == k_cond*k_ins^2/sqrt(3)*(
               (g_br+(4/3)*g_fr_br)*vm_fr^2
               + (-g_br)*(vm_fr*vm_to)
            )
         )
         # @constraint(
         #    model,
         #    q_fr == 0
         # )

         # To side of the branch flow
         @NLconstraint(
            model,
            p_to == k_cond*k_ins^2/sqrt(3)*(
               (g_br+(4/3)*g_to_br)*vm_to^2
               + (-g_br)*(vm_to*vm_fr)
            )
         )
         # @constraint(
         #    model,
         #    q_to == 0
         # )

         # @constraint(
         #    model,
         #    va_fr == 0
         # )
         # @constraint(
         #    model,
         #    va_to == 0
         # )

      else
         q_fr = q[subnet_idx][f_idx]                     # q_fr is a reference to the optimization variable q[f_idx]
         q_to = q[subnet_idx][t_idx]                     # q_to is a reference to the optimization variable q[t_idx]
         va_fr = va[subnet_idx][branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
         va_to = va[subnet_idx][branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch
         # From side of the branch flow
         @NLconstraint(
            model,
            p_fr == (
               (g_br+g_fr_br)/tm*vm_fr^2
               + (-g_br*tr+b_br*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to))
               + (-b_br*tr-g_br*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to))
            )
         )
         @NLconstraint(
            model,
            q_fr == (
               -(b_br+b_fr_br)/tm*vm_fr^2
               - (-b_br*tr-g_br*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to))
               + (-g_br*tr+b_br*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to))
            )
         )

         # To side of the branch flow
         @NLconstraint(
            model,
            p_to == (
               (g_br+g_to_br)*vm_to^2
               + (-g_br*tr-b_br*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr))
               + (-b_br*tr+g_br*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr))
            )
         )
         @NLconstraint(
            model,
            q_to == (
               -(b_br+b_to_br)*vm_to^2
               - (-b_br*tr+g_br*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to))
               + (-g_br*tr-b_br*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr))
            )
         )

         # Voltage angle difference limit
         if !(subnet_idx in keys(constraints[:theta_ulim]))
            constraints[:theta_ulim][subnet_idx] = Dict{Int64,Any}()
         end
         if !(subnet_idx in keys(constraints[:theta_llim]))
            constraints[:theta_llim][subnet_idx] = Dict{Int64,Any}()
         end
         if (("ignore_anglims" in keys(branch)) && (branch["ignore_anglims"])) #|| branch["transformer"]
            # Angle limits ignored on this branch. Set to +/- 90 for convergence.
            println("Setting +/- 90.0 deg angle limits on branch $(branch["index"])")
            constraints[:theta_ulim][subnet_idx][i] = @constraint(model, va_fr - va_to <= pi/2)
            constraints[:theta_llim][subnet_idx][i] = @constraint(model, va_fr - va_to >= -pi/2)
         else
            # If an alternative angle limit is defined, as in the case of series compensation with an intermediate bus
            if ("alt_ang_lim" in keys(branch)) && (length(branch["alt_ang_lim"])==2)
               (alt_f_bus, alt_t_bus) = branch["alt_ang_lim"]
               # println("ref branch: $branch")
               alt_va_fr = va[subnet_idx][alt_f_bus]
               alt_va_to = va[subnet_idx][alt_t_bus]
               # println("subnet $subnet_idx, branch $(branch["index"])")
               # println("alt_va_fr: $alt_va_fr")
               # println("alt_va_to: $alt_va_to")
               constraints[:theta_ulim][subnet_idx][i] = @constraint(model, alt_va_fr - alt_va_to <= branch["angmax"])
               constraints[:theta_llim][subnet_idx][i] = @constraint(model, alt_va_fr - alt_va_to >= branch["angmin"])

               # Set limits for the branch to +/- 90 for convergence.
               @constraint(model, va_fr - va_to <= pi/2)
               @constraint(model, va_fr - va_to >= -pi/2)
            else
               constraints[:theta_ulim][subnet_idx][i] = @constraint(model, va_fr - va_to <= branch["angmax"])
               constraints[:theta_llim][subnet_idx][i] = @constraint(model, va_fr - va_to >= branch["angmin"])
            end
         end
      end

      # Apparent power limit, from side and to side
      if "rate_a" in keys(branch)
         if !(subnet_idx in keys(constraints[:s_ulim]))
            constraints[:s_ulim][subnet_idx] = Dict{Int64,Any}()
         end
         if !(subnet_idx in keys(constraints[:s_llim]))
            constraints[:s_llim][subnet_idx] = Dict{Int64,Any}()
         end
         if dc_subnet && dc_current_limit
            constraints[:s_ulim][subnet_idx][i] = @constraint(model, p_fr^2 <= k_ins*k_cond*branch["rate_a"]^2)
            constraints[:s_llim][subnet_idx][i] = @constraint(model, p_to^2 <= k_ins*k_cond*branch["rate_a"]^2)
         elseif dc_subnet
            constraints[:s_ulim][subnet_idx][i] = @constraint(model, p_fr^2 <= branch["rate_a"]^2)
            constraints[:s_llim][subnet_idx][i] = @constraint(model, p_to^2 <= branch["rate_a"]^2)
         else
            constraints[:s_ulim][subnet_idx][i] = @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
            constraints[:s_llim][subnet_idx][i] = @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
         end
      end
      # print("added power flow constraints\n")
   end

   # HVDC line constraints
   for (i,dcline) in ref_subnet[:dcline]
      # Build the from variable id of the i-th HVDC line, which is a tuple given by (hvdc line id, from bus, to bus)
      f_idx = (i, dcline["f_bus"], dcline["t_bus"])
      # Build the to variable id of the i-th HVDC line, which is a tuple given by (hvdc line id, to bus, from bus)
      t_idx = (i, dcline["t_bus"], dcline["f_bus"])   # index of the ith HVDC line which is a tuple given by (line number, to bus, from bus)
      # note: it is necessary to distinguish between the from and to sides of a HVDC line due to power losses

      # Constraint defining the power flow and losses over the HVDC line
      @constraint(model, (1-dcline["loss1"])*p_dc[subnet_idx][f_idx] + (p_dc[subnet_idx][t_idx] - dcline["loss0"]) == 0)
   end
   # print("finished adding all constraints for network $(subnet_idx)\n")
   # print("adding cost variables\n")

   # assumes costs are given as quadratic functions
   cost_pg[subnet_idx] = @variable(model, base_name = "cost_pg_$(subnet_idx)")
   @constraint(
      model,
      cost_pg[subnet_idx] == sum(
      gen["cost"][1]*pg[subnet_idx][i]^2
      + gen["cost"][2]*pg[subnet_idx][i]
      + gen["cost"][3]
         for (i,gen) in ref_subnet[:gen]
      )
   )
   # index representing which side the HVDC line is starting
   from_idx = Dict(arc[1] => arc for arc in ref_subnet[:arcs_from_dc])
   cost_dcline[subnet_idx] = @variable(model, base_name = "cost_dcline_$(subnet_idx)")
   @constraint(
      model,
      cost_dcline[subnet_idx] == sum(
      dcline["cost"][1]*p_dc[subnet_idx][from_idx[i]]^2
      + dcline["cost"][2]*p_dc[subnet_idx][from_idx[i]]
      + dcline["cost"][3]
         for (i,dcline) in ref_subnet[:dcline]
      )
   )
   if (obj=="mincost") || (obj=="areagen") || (obj=="zonegen")
      if (obj=="areagen")
         areagens = keys(filter(p->(ref_subnet[:bus][last(p)["gen_bus"]]["area"] in gen_areas), ref_subnet[:gen]))
         # println("subnet $subnet_idx")
         # println("areagens: $(areagens)")
         area_pg[subnet_idx] = @variable(model, base_name = "area_pg_$(subnet_idx)")
         @constraint(
            model,
            area_pg[subnet_idx] == sum(
               pg[subnet_idx][i]
               for i in areagens
            )
         )
      end
      if (obj=="zonegen")
         zonegens = keys(filter(p->(ref_subnet[:bus][last(p)["gen_bus"]]["zone"] in gen_zones), ref_subnet[:gen]))
         # println("subnet $subnet_idx")
         # println("zonegens: $(zonegens)")
         zone_pg[subnet_idx] = @variable(model, base_name = "zone_pg_$(subnet_idx)")
         @constraint(
            model,
            zone_pg[subnet_idx] == sum(
               pg[subnet_idx][i]
               for i in zonegens
            )
         )
      end
   elseif obj=="minredispatch"
      # Get base case dispatch
      redispatch_pg[subnet_idx] = @variable(model, base_name = "dev_pg_$(subnet_idx)")
      @constraint(
         model,
         redispatch_pg[subnet_idx] == sum(
            (gen["pg"] - pg[subnet_idx][i])^2
            for (i,gen) in ref_subnet[:gen]
         )
      )
      redispatch_qg[subnet_idx] = @variable(model, base_name = "dev_qg_$(subnet_idx)")
      @constraint(
         model,
         redispatch_qg[subnet_idx] == sum(
            (gen["qg"] - qg[subnet_idx][i])^2
            for (i,gen) in ref_subnet[:gen]
         )
      )
      # redispatch_vm[subnet_idx] = @variable(model, base_name = "dev_vm_$(subnet_idx)")
      # @constraint(
      #    model,
      #    redispatch_vm[subnet_idx] == sum(
      #       (ref_subnet[:bus][gen["gen_bus"]]["vm"] - vm[subnet_idx][gen["gen_bus"]])^2
      #       for (i,gen) in ref_subnet[:gen]
      #    )
      # )
      redispatch_vm[subnet_idx] = @variable(model, base_name = "dev_vm_$(subnet_idx)")
      @constraint(
         model,
         redispatch_vm[subnet_idx] == sum(
            (bus["vm"] - vm[subnet_idx][i])^2
            for (i,bus) in ref_subnet[:bus] if (bus["bus_type"] in [1,2,3])
         )
      )
      redispatch_va[subnet_idx] = @variable(model, base_name = "dev_va_$(subnet_idx)")
      @constraint(
         model,
         redispatch_va[subnet_idx] == sum(
            (bus["va"] - va[subnet_idx][i])^2
            for (i,bus) in ref_subnet[:bus] if (bus["bus_type"] in [1,2,3])
         )
      )
      # println([gen["vg"] for (i,gen) in ref_subnet[:gen]])
      # println([ref_subnet[:bus][gen["gen_bus"]]["vm"] for (i,gen) in ref_subnet[:gen]])
   end
end

function add_gen_zone_scaling_constraint!(
   ref_data,
   model,
   pg,
   constraints,
   unbounded_pg,
   gen_zones
   )

   alpha_upstream = @variable(
      model,
      base_name="alpha_upstream",
      lower_bound=0,
      start=1
   )
   alpha_downstream = @variable(
      model,
      base_name="alpha_downstream",
      lower_bound=0,
      start=1
   )
   for (subnet_idx, ref_subnet) in ref_data
      for i in keys(ref_subnet[:gen])
         if ref_subnet[:bus][ref_subnet[:gen][i]["gen_bus"]]["zone"] in gen_zones
            # This is a downstream generator (objective is to minimize downstream generation)
            println("adding downstream constraint for gen $i: starting value $(ref_subnet[:gen][i]["pg"])")
            @constraint(
               model,
               pg[subnet_idx][i] == alpha_downstream * ref_subnet[:gen][i]["pg"]
            )
         else
            # This is an upstream generator
            println("adding upstream constraint for gen $i: starting value $(ref_subnet[:gen][i]["pg"])")
            @constraint(
               model,
               pg[subnet_idx][i] == alpha_upstream * ref_subnet[:gen][i]["pg"]
            )
         end
      end
      if unbounded_pg
         delete_upper_bound.(pg[subnet_idx])
         delete!(constraints[:pg_ulim],subnet_idx)
         delete_lower_bound.(pg[subnet_idx])
         delete!(constraints[:pg_llim],subnet_idx)
      end
   end
   return alpha_upstream, alpha_downstream
end


text:
######### StanSample example  ###########

using StanSample

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = [
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 0, 1, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 1, 0]),
  Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 1, 0, 0, 1]),
]

# Keep tmpdir identical across multiple runs to prevent re-compilation
stanmodel = SampleModel(
  "bernoulli", bernoulli_model;
  method =  StanSample.Sample(adapt=StanSample.Adapt(delta=0.85)))

rc = stan_sample(stanmodel, data=bernoulli_data)

if success(rc)

  # Fetch the same output in the `sdf` ChainDataFrame
  sdf = read_summary(stanmodel)

  @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05

end


text:
# clone test_suite repo
sbml_test_suite()

N = 100
println("****SBML TEST SUITE TESTING****")
suite_fns = get_sbml_suite_fns()[1:N]
fn = suite_fns[1]

@test isfile(fn)
@test readSBML(fn, doc -> begin
        set_level_and_version(3, 2)(doc)
        convert_simplify_math(doc)
    end) isa SBML.Model

now_fmtd = Dates.format(now(), dateformat"yyyy-mm-dd\THH-MM-SS")
suite_df = lower_fns(suite_fns; write_fn="test_suite_$(now_fmtd).csv", verbose=true)

num_good = nrow(filter(:retcode => x-> x==5, suite_df))
@info "suite num_good: $num_good / $N"

ds = filter(isdir, readdir(joinpath(datadir, "sbml-test-suite", "semantic"); join=true))[1:N]
df = verify_all(ds)
CSV.write(joinpath(logdir, "suite_verified_$(now_fmtd).csv"), df)


text:
module InversePreSolve

import PartiallyObservedInverseGames.ForwardOptimalControl
using JuMP: @objective

export pre_solve

function pre_solve(
    y_obs,
    u_obs;
    control_system,
    observation_model = (; expected_observation = identity),
    u_regularization = 0,
    inner_solver_kwargs...,
)
    function presolve_objective(; x, u, y_obs, u_obs)
        y_expected = observation_model.expected_observation(x)
        sum(el -> el^2, y_expected - y_obs) +
        (isnothing(u_obs) ? 0 : sum(el -> el^2, u - u_obs)) +
        (iszero(u_regularization) ? 0 : u_regularization * sum(el -> el^2, u))
    end

    reconstruction_cost_model = (;
        # TODO: make configurable perhaps
        add_objective! = function (opt_model, x, u; kwargs...)
            @objective(opt_model, Min, presolve_objective(; x, u, y_obs, u_obs))
        end,
        weights = (),
    )
    T = size(y_obs, 2)

    ForwardOptimalControl.solve_optimal_control(
        control_system,
        reconstruction_cost_model,
        nothing,
        T;
        inner_solver_kwargs...,
    )
end

end


text:
using ControlSystems
using LinearAlgebra

##
# system of the form
# 1 | 2 3
# ———————
# 4 | 5 6
# 7 | 8 9
sys1 = ControlSystems.PartionedStateSpace(ss(1.0, [2.0 3.0], [4.0; 7.0], [5.0 6.0; 8.0 9.0]), 1, 1)

matrix(x::Number) = fill(x, 1, 1)

sys1.A == matrix(1.0)
sys1.B1 == matrix(2.0)
sys1.B2 == matrix(3.0)
sys1.C1 == matrix(4.0)
sys1.C2 == matrix(7.0)
sys1.D11 == matrix(5.0)
sys1.D12 == matrix(6.0)
sys1.D21 == matrix(8.0)
sys1.D22 == matrix(9.0)

##
sys2 = ControlSystems.PartionedStateSpace(ss(fill(1.0, 2, 2), fill(2.0, 2, 5), fill(3.0, 7, 2), fill(4.0, 7, 5)), 2, 3)

@test sys2.A == fill(1.0, 2, 2)
@test sys2.B1 == fill(2.0, 2, 2)
@test sys2.B2 == fill(2.0, 2, 3)
@test sys2.C1 == fill(3.0, 3, 2)
@test sys2.C2 == fill(3.0, 4, 2)
@test sys2.D11 == fill(4.0, 3, 2)
@test sys2.D12 == fill(4.0, 3, 3)
@test sys2.D21 == fill(4.0, 4, 2)
@test sys2.D22 == fill(4.0, 4, 3)

# TODO: Add some tests for interconnections, implicitly tested through delay system implementations though
@test (sys1 + sys1).P[1, 1] == (sys1.P[1,1] + sys1.P[1,1])


text:
using LinearAlgebra

## Test standard basis polynomials with (nearly) the same tests

#   compare upto trailing  zeros
function  upto_tz(as, bs)
    n,m = findlast.(!iszero, (as,bs))
    n == m || return false
    n == nothing &&  return true
    for i in 1:n
        as[i] != bs[i] && return false
    end
    true
end

upto_z(as, bs) = upto_tz(filter(!iszero,as), filter(!iszero,bs))

# compare upto trailing zeros infix operator
==ᵗ⁰(a,b) = upto_tz(a,b)
==ᵗᶻ(a,b) = upto_z(a,b)

Ps = (ImmutablePolynomial, Polynomial, SparsePolynomial, LaurentPolynomial)

@testset "Construction" for coeff in [
                                      Int64[1, 1, 1, 1],
                                      Float32[1, -4, 2],
                                      ComplexF64[1 - 1im, 2 + 3im],
                                      [3 // 4, -2 // 1, 1 // 1]
                                     ]

    for P in Ps
        p = P(coeff)
        @test coeffs(p) ==ᵗ⁰ coeff
        @test degree(p) == length(coeff) - 1
        @test p.var == :x
        P == Polynomial && @test length(p) == length(coeff)
        P == Polynomial && @test size(p) == size(coeff)
        P == Polynomial && @test size(p, 1) == size(coeff, 1)
        P == Polynomial && @test typeof(p).parameters[1] == eltype(coeff)
        @test eltype(p) == eltype(coeff)
        @test all([-200, -0.3, 1, 48.2] .∈ domain(p))
    end
end

@testset "Mapdomain" begin
    for P in Ps
        x = -30:20
        mx = mapdomain(P, x)
        @test mx == x

        x = 0.5:0.01:0.6
        mx = mapdomain(P, x)
        @test mx == x
    end
end

@testset "Other Construction" begin
    for P in (ImmutablePolynomial, Polynomial, SparsePolynomial, LaurentPolynomial)

        # Leading 0s
        p = P([1, 2, 0, 0])
        @test coeffs(p) ==ᵗ⁰ [1, 2]
        P == Polynomial && @test length(p) == 2

        # different type
        p = P{Float64}(ones(Int32, 4))
        @test coeffs(p) ==ᵗ⁰ ones(Float64, 4)

        p = P(30)
        @test coeffs(p) ==ᵗ⁰ [30]

        p = zero(P{Int})
        @test coeffs(p) ==ᵗ⁰ [0]

        p = one(P{Int})
        @test coeffs(p) ==ᵗ⁰ [1]

        pNULL = P(Int[])
        @test iszero(pNULL)
        P != LaurentPolynomial && @test degree(pNULL) == -1

        p0 = P([0])
        @test iszero(p0)
        P != LaurentPolynomial && @test degree(p0) == -1

        # P(2) is  2 (not  2p₀)  connvert(Polynomial, P(s::Number)) = Polynomial(s)
        @test convert(Polynomial, P(2)) ≈ Polynomial(2)
        @test P(2)  ≈ 2*one(P)

        # variable(), P() to generate `x` in given basis
        @test degree(variable(P)) == 1
        @test variable(P)(1) == 1
        @test degree(P()) == 1
        @test P()(1) == 1
        @test variable(P, :y) == P(:y)

        # test degree, isconstant
        P != LaurentPolynomial &&  @test degree(zero(P)) ==  -1
        @test degree(one(P)) == 0
        @test degree(P(1))  ==  0
        @test degree(P([1]))  ==  0
        @test degree(P(:x)) ==  1
        @test degree(variable(P)) == 1
        @test degree(Polynomials.basis(P,5)) == 5
        @test Polynomials.isconstant(P(1))
        @test !Polynomials.isconstant(variable(P))
    end
end

@testset "Arithmetic" begin

    for  P in Ps
        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test p3 == P([1,2,1])
        @test pN * 10 == P([2760, 30, 870, 150, 240])
        @test pN / 10.0 == P([27.6, 0.3, 8.7, 1.5, 2.4])
        @test 10 * pNULL + pN == pN
        @test 10 * p0 + pN == pN
        @test p5 + 2 * p1 == P([3,4,6,4,1])
        @test 10 * pNULL - pN == -pN
        @test p0 - pN == -pN
        @test p5 - 2 * p1 == P([-1,4,6,4,1])
        @test p2 * p2 * p2 == p4
        @test p2^4 == p5
        @test pNULL^3 == pNULL
        @test pNULL * pNULL == pNULL

        @test pNULL + 2 == p0 + 2 == 2 + p0 == P([2])
        @test p2 - 2 == -2 + p2 == P([-1,1])
        @test 2 - p2 == P([1,-1])

    end

    for P in Ps # ensure promotion of scalar +,*,/
        p = P([1,2,3])
        @test p + 0.5 == P([1.5, 2.0, 3.0])
        @test p / 2  == P([1/2, 1.0, 3/2])
        @test p * 0.5 == P([1/2, 1.0, 3/2])
    end

    # ensure  promotion of +,*; issue 215
    for P in Ps
        p,q = P([1,2,3]), P(im, :θ)
        @test p+q == P([1+im, 2, 3])
        @test p*q == P(im*[1,2,3])
    end
end

@testset "Divrem" begin
    for P in  Ps

        p0 = P([0])
        p1 = P([1])
        p2 = P([5, 6, -3, 2 ,4])
        p3 = P([7, -3, 2, 6])
        p4 = p2 * p3
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test divrem(p4, p2) == (p3, zero(p3))
        @test p3 % p2 == p3
        @test all((map(abs, coeffs(p2 ÷ p3 - P([1 / 9,2 / 3])))) .< eps())
        @test divrem(p0, p1) == (p0, p0)
        @test divrem(p1, p1) == (p1, p0)
        @test divrem(p2, p2) == (p1, p0)
        @test divrem(pR, pR) == (one(pR), zero(pR))
        @test_throws DivideError p1 ÷ p0
        @test_throws DivideError divrem(p0, p0)

        # issue #235
        num = P([0.8581454436924945, 0.249671302254737, 0.8048498901050951, 0.1922713965697087]) # degree 3 polynomial
        den = P([0.9261520696359462, 0.07141031902098072, 0.378071465860349]) # degree 2 polynomial
        q, r = divrem(num,den)  # expected degrees: degree(q) = degree(num)-degree(den) = 1, degree(r) = degree(den)-1 = 1
        @test num ≈ den*q+r  # true
        @test degree(q) == 1 # true
        degree(r) < degree(den)
    end
end

@testset "Comparisons" begin
    for P in Ps
        pX = P([1, 2, 3, 4, 5])
        pS1 = P([1, 2, 3, 4, 5], "s")
        pS2 = P([1, 2, 3, 4, 5], 's')
        pS3 = P([1, 2, 3, 4, 5], :s)
        @test pX != pS1
        @test pS1 == pS2
        @test pS1 == pS3
        @test_throws ErrorException pS1 + pX
        @test_throws ErrorException pS1 - pX
        @test_throws ErrorException pS1 * pX
        @test_throws ErrorException pS1 ÷ pX
        @test_throws ErrorException pS1 % pX

        # Testing copying.
        pcpy1 = P([1,2,3,4,5], :y)
        pcpy2 = copy(pcpy1)
        @test pcpy1 == pcpy2

        # Check for isequal
        p1 = P([1.0, -0.0, 5.0, Inf])
        p2 = P([1.0,  0.0, 5.0, Inf])
        p3 = P([0, NaN])

        P != SparsePolynomial && (@test p1 == p2 && !isequal(p1, p2))  # SparsePolynomial doesn't store -0.0,  0.0.
        @test p3 === p3 && p3 ≠ p3 && isequal(p3, p3)

        p = fromroots(P, [1,2,3])
        q = fromroots(P, [1,2,3])
        @test hash(p) == hash(q)

        p1s = P([1,2], :s)
        p1x = P([1,2], :x)
        p2s = P([1], :s)

        @test p1s == p1s
        @test p1s ≠ p1x
        @test p1s ≠ p2s

        @test_throws ErrorException p1s ≈ p1x
        @test p1s ≉ p2s
        @test p1s ≈ P([1,2.], :s)

        @test p2s ≈ 1.0 ≈ p2s
        @test p2s == 1.0 == p2s
        @test p2s ≠ 2.0 ≠ p2s
        @test p1s ≠ 2.0 ≠ p1s

        @test nnz(map(P, sparse(1.0I, 5, 5))) == 5

        @test P([0.5]) + 2 == P([2.5])
        @test 2 - P([0.5]) == P([1.5])

        # check ≈ for P matches usage for Vector{T} (possibly padded with trailing zeros)
        @test (P([NaN]) ≈ P([NaN]))               == ([NaN] ≈ [NaN]) # false
        @test (P([NaN]) ≈ NaN)                    == (false)
        @test (P([Inf]) ≈ P([Inf]))               == ([Inf] ≈ [Inf]) # true
        @test (P([Inf]) ≈ Inf)                    == (true)
        @test (P([1,Inf]) ≈ P([0,Inf]))           == ([1,Inf] ≈ [0,Inf]) # false
        @test (P([1,NaN,Inf]) ≈ P([0,NaN, Inf]))  == ([1,NaN,Inf] ≈ [0,NaN, Inf]) #false
        @test (P([eps(), eps()]) ≈ P([0,0]))      == ([eps(), eps()] ≈ [0,0]) # false
        @test (P([1,eps(), 1]) ≈ P([1,0,1]))      == ([1,eps(), 1] ≈ [1,0,1]) # true
        @test (P([1,2]) ≈ P([1,2,eps()]))         == ([1,2,0] ≈ [1,2,eps()])

        # check how ==, ===, isapprox ignore variable mismatch when constants are involved, issue #217, issue #219
        @test zero(P, :x) == zero(P, :y)
        @test one(P, :x) == one(P, :y)
        @test !(variable(P, :x) == variable(P,:y))

        @test !(zero(P, :x) === zero(P, :y))
        @test !(one(P, :x) === one(P, :y))
        @test !(variable(P, :x) === variable(P,:y))

        @test zero(P, :x) ≈ zero(P, :y)
        @test one(P, :x) ≈ one(P, :y)
        @test (variable(P, :x) ≈ variable(P, :x))
        @test_throws ErrorException variable(P, :x) ≈ variable(P, :y)

    end
end

@testset "Fitting" begin
    for P in Ps
        xs = range(0, stop = π, length = 10)
        ys = sin.(xs)

        p = fit(P, xs, ys)
        y_fit = p.(xs)
        abs_error = abs.(y_fit .- ys)
        @test maximum(abs_error) <= 0.03

        p = fit(P, xs, ys, 2)
        y_fit = p.(xs)
        abs_error = abs.(y_fit .- ys)
        @test maximum(abs_error) <= 0.03

        # Test weighted
        for W in [1, ones(size(xs)), diagm(0 => ones(size(xs)))]
            p = fit(P, xs, ys, 2, weights = W)
            @test p.(xs) ≈ y_fit
        end

        # Getting error on passing Real arrays to polyfit #146
        xx = Real[20.0, 30.0, 40.0]
        yy = Real[15.7696, 21.4851, 28.2463]
        fit(P, xx, yy, 2)

        # issue #214 --  should error
        @test_throws MethodError fit(Polynomial, rand(2,2), rand(2,2))

    end

    # test default   (issue  #228)
    fit(1:3,  rand(3))
end

@testset "Values" begin
    for P in Ps
        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        @test fromroots(P, Int[])(2.) == 1.
        @test pN(-.125) == 276.9609375
        @test pNULL(10) == 0
        @test p0(-10) == 0
        @test fromroots(P, [1 // 2, 3 // 2])(1 // 2) == 0 // 1

        # Check for Inf/NaN operations
        p1 = P([Inf, Inf])
        p2 = P([0, Inf])
        @test p1(Inf) == Inf
        @test isnan(p1(-Inf))
        @test isnan(p1(0))
        @test p2(-Inf) == -Inf

        # issue #189
        p = P([0,1,2,3])
        A = [0 1; 0  0];
        @test  p(A) == A  + 2A^2 + 3A^3

        # issue #209
        ps  = [P([0,1]), P([0,0,1])]
        @test Polynomials.evalpoly.(1/2, ps) ≈ [p(1/2)  for  p  in ps]

    end

    # constant polynomials and type
    Ts = (Int, Float32, Float64, Complex{Int}, Complex{Float64})
    for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
        for T in Ts
            for S in Ts
                c = 2
                p = P{T}(c)
                x = one(S)
                y = p(x)
                @test y === c * one(T)*one(S)
                q = P{T}([c,c])
                @test typeof(q(x)) == typeof(p(x))
            end
        end
    end

    for P in Ps
        p = P(1)
        x = [1 0; 0 1]
        y = p(x)
        @test y == x

        # Issue #208 and  type of output
        p1=P([1//1])
        p2=P([0, 0.9])
        p3=p1(p2)
        @test isa(p3, P)
        @test eltype(p3) == eltype(p2)
    end

    # compensated_horner
    # polynomial evaluation for polynomials with large condition numbers
    for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
        x = variable(P{Float64})
        f(x) = (x - 1)^20
        p = f(x)
        e₁ = abs( (f(4/3) - p(4/3))/ p(4/3) )
        e₂ = abs( (f(4/3) - Polynomials.compensated_horner(p, 4/3))/ p(4/3) )
        λ = cond(p, 4/3)
        u = eps()/2
        @test λ > 1/u
        @test e₁ <= 2 * 20 * u * λ
        @test e₁ > u^(1/4)
        @test e₂ <= u + u^2 * λ * 100
    end
end

@testset "Conversion" begin
    for P in (Polynomial, SparsePolynomial, LaurentPolynomial)
        p = P([0,one(Float64)])
        @test P{Complex{Float64}} == typeof(p + 1im)
        @test P{Complex{Float64}} == typeof(1im - p)
        @test P{Complex{Float64}} == typeof(p * 1im)
    end

    for P in (ImmutablePolynomial,)
        p = P([0,one(Float64)])
        N=2
        @test P{Complex{Float64},N} == typeof(p + 1im)
        @test P{Complex{Float64},N} == typeof(1im - p)
        @test P{Complex{Float64},N} == typeof(p * 1im)
    end

    # unnecessary copy in convert #65
    p1 = Polynomial([1,2])
    p2 = convert(Polynomial{Int64}, p1)
    p2[3] = 3
    @test p1[3] == 3

end

@testset "Roots" begin
    for P in Ps

        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        # From roots
        r = [2, 3]
        @test fromroots(r) == Polynomial([6, -5, 1])
        p = fromroots(P, r)
        @test p == P([6, -5, 1])
        @test sort(roots(p)) ≈ r

        @test roots(p0) == roots(p1) == roots(pNULL) == []
        @test P == LaurentPolynomial ? roots(variable(P)) == [0.0] : roots(P([0,1,0])) == [0.0]

        @test roots(p2) == [-1]
        a_roots = [c for c in coeffs(copy(pN))]
        @test all(map(abs, sort(roots(fromroots(a_roots))) - sort(a_roots)) .< 1e6)
        @test length(roots(p5)) == 4
        @test roots(pNULL) == []
        @test sort(roots(pR)) == [1 // 2, 3 // 2]

        @test sort(roots(LaurentPolynomial([24,10,-15,0,1],-2,:z)))≈[-4.0,-1.0,2.0,3.0]

        A = [1 0; 0 1]
        @test fromroots(A) == Polynomial(Float64[1, -2, 1])
        p = fromroots(P, A)
        @test p == P(Float64[1, -2, 1])
        @test roots(p) ≈ sort!(eigvals(A), rev = true)

        x = variable()
        plarge = 8.362779449448982e41 - 2.510840694154672e57x + 4.2817430781178795e44x^2 - 1.6225927682921337e31x^3 + 1.0x^4  # #120
        @test length(roots(plarge)) == 4

        @test begin
            a = P([1,1,1])*P([1,0.5,1])*P([1,1])    # two complex conjugate pole pairs and one real pole
            r = roots(a)
            b = fromroots(r)
            (b ≈ a) & isreal(coeffs(b))    # the coeff should be real
        end
    end
end

@testset "Integrals and Derivatives" begin
    # Integrals derivatives
    for P in Ps

        pNULL = P(Int[])
        p0 = P([0])
        p1 = P([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        p2 = P([1,1,0,0])
        p3 = P([1,2,1,0,0,0,0])
        p4 = P([1,3,3,1,0,0])
        p5 = P([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
        pN = P([276,3,87,15,24,0])
        pR = P([3 // 4, -2 // 1, 1 // 1])

        c = [1, 2, 3, 4]
        p = P(c)

        der = derivative(p)
        @test coeffs(der) ==ᵗ⁰ [2, 6, 12]
        int = integrate(der, 1)
        @test coeffs(int) ==ᵗ⁰ c

        @test derivative(pR) == P([-2 // 1,2 // 1])
        @test derivative(p3) == P([2,2])
        @test derivative(p1) == derivative(p0) == derivative(pNULL) == pNULL
        @test_throws ErrorException derivative(pR, -1)
        @test integrate(P([1,1,0,0]), 0, 2) == 4.0

        @test derivative(integrate(pN)) == convert(P{Float64}, pN)
        @test integrate(pNULL, 1) == convert(P{Float64}, p1)
        rc = Rational{Int64}[1,2,3]
        @test integrate(P(rc)) == P{eltype(rc)}([0, 1, 1, 1])

        for i in 1:10
            p = P(rand(1:5, 6))
            @test degree(truncate(p - integrate(derivative(p)), atol=1e-13)) <= 0
            @test degree(truncate(p - derivative(integrate(p)), atol=1e-13)) <= 0
        end

        # Handling of `NaN`s
        p     = P([NaN, 1, 5])
        pder  = derivative(p)
        pint  = integrate(p)

        @test isnan(p(1)) # p(1) evaluates to NaN
        @test isequal(pder, P([NaN]))
        @test isequal(pint, P([NaN]))

        pint  = integrate(p, 0.0im)
        @test isequal(pint, P([NaN]))

        # Issue with overflow and polyder Issue #159
        @test derivative(P(BigInt[0, 1])^100, 100) == P(factorial(big(100)))
    end
end

@testset "Elementwise Operations" begin
    for P in Ps
        p1  = P([1, 2])
        p2  = P([3, 1.])
        p   = [p1, p2]
        q   = [3, p1]
        if P !=  ImmutablePolynomial
            @test q isa Vector{typeof(p1)}
            @test p isa Vector{typeof(p2)}
        else
            @test q isa Vector{P{eltype(p1)}} # ImmutablePolynomial{Int64,N} where {N}, different  Ns
            @test p isa Vector{P{eltype(p2)}} # ImmutablePolynomial{Int64,N} where {N}, different  Ns
        end

        psum  = p .+ 3
        pprod = p .* 3
        pmin  = p .- 3
        @test psum  isa Vector{typeof(p2)}
        @test pprod isa Vector{typeof(p2)}
        @test pmin  isa Vector{typeof(p2)}
    end
end

@testset "Chop and Truncate" begin
    # chop and truncate
    for P in Ps
        if P == Polynomial
            p = P([1, 1, 1, 1])
            coeffs(p)[end] = 0
            @assert coeffs(p) == [1, 1, 1, 0]
            p = chop(p)
        else
            p = P([1, 1, 1, 0])
        end

        @test coeffs(p) ==ᵗ⁰ [1, 1, 1]
        ## truncation
        p1 = P([1,1] / 10)
        p2 = P([1,2] / 10)
        p3 = P([1,3] / 10)
        psum = p1 + p2 - p3
        @test degree(psum) == 1         # will have wrong degree
        @test degree(truncate(psum)) == 0 # the degree should be correct after truncation

        @test truncate(P([2,1]), rtol = 1 / 2, atol = 0) == P([2])
        @test truncate(P([2,1]), rtol = 1, atol = 0)   == P([0])
        @test truncate(P([2,1]), rtol = 0, atol = 1)   == P([2])

        pchop = P([1, 2, 3, 0, 0, 0])
        pchopped = chop(pchop)
        @test roots(pchop) == roots(pchopped)

    end
end

@testset "As matrix elements" begin

    for P in Ps
        p = P([1,2,3], :x)
        A = [1 p; p^2 p^3]
        @test !issymmetric(A)
        @test issymmetric(A*transpose(A))
        diagm(0 => [1, p^3], 1=>[p^2], -1=>[p])
    end

    # issue 206 with mixed variable types and promotion
    for P in Ps
        λ = P([0,1],:λ)
        A = [1 λ; λ^2 λ^3]
        @test A ==  diagm(0 => [1, λ^3], 1=>[λ], -1=>[λ^2])
        @test all([1 -λ]*[λ^2 λ; λ 1] .== 0)
        @test [λ 1] + [1 λ] == (λ+1) .* [1 1] # (λ+1) not a number, so we broadcast
    end
end

@testset "Linear Algebra" begin
    for P in Ps
        p = P([3, 4])
        @test norm(p) == 5
        p = P([-1, 3, 5, -2])
        @test norm(p) ≈ 6.244997998398398
        p = P([1 - 1im, 2 - 3im])
        p2 = conj(p)
        @test coeffs(p2) ==ᵗ⁰ [1 + 1im, 2 + 3im]
        @test transpose(p) == p
        P != ImmutablePolynomial && @test transpose!(p) == p
        @test adjoint(Polynomial(im)) == Polynomial(-im) # issue 215
        @test conj(Polynomial(im)) == Polynomial(-im) # issue 215

        @test norm(P([1., 2.])) == norm([1., 2.])
        @test norm(P([1., 2.]), 1) == norm([1., 2.], 1)
    end

    ## Issue #225 and different meanings for "conjugate"
    P = LaurentPolynomial
    p = P(rand(Complex{Float64}, 4), -1:2)
    z = rand(Complex{Float64})
    s = imag(z)*im
    @test conj(p)(z) ≈ (conj ∘ p ∘ conj)(z)
    @test Polynomials.paraconj(p)(z) ≈ (conj ∘ p ∘ conj ∘ inv)(z)
    @test Polynomials.cconj(p)(s) ≈ (conj ∘ p)(s)

end

@testset "Indexing" begin
    # Indexing
    for P in Ps
        # getindex
        p = P([-1, 3, 5, -2])
        @test p[0] == -1
        @test p[[1, 2]] ==ᵗ⁰ [3, 5]
        @test p[1:2] ==ᵗ⁰ [3, 5]
        @test p[:] ==ᵗ⁰ [-1, 3, 5, -2]

        if P != ImmutablePolynomial
            # setindex
            p1  = P([1,2,1])
            p1[5] = 1
            @test p1[5] == 1
            @test p1 == P([1,2,1,0,0,1])

            @test p[end] == coeffs(p)[end]

            if P != SparsePolynomial
                p[1] = 2
                @test coeffs(p)[2] == 2
                p[2:3] = [1, 2]
                @test coeffs(p)[3:4] == [1, 2]
                p[0:1] = 0
                @test coeffs(p)[1:2] ==ᵗ⁰ [0, 0]

                p[:] = 1
                @test coeffs(p) ==ᵗ⁰ ones(4)
            end

            p[:] = 0
            @test chop(p) ≈ zero(p)
        end

        p1 = P([1,2,0,3])
        for term in p1
            @test isa(term, P)
        end

        @test eltype(p1) == Int
        for P in Ps
            p1 = P([1,2,0,3])
            @test eltype(collect(p1)) <: P{Int}
            @test eltype(collect(P{Float64}, p1)) <: P{Float64}
            @test_throws InexactError collect(P{Int}, P([1.2]))
        end

        p1 = P([1,2,0,3])
        @test length(collect(p1)) == degree(p1) + 1

        @test [p1[idx] for idx in eachindex(p1)] ==ᵗᶻ [1,2,0,3]
    end
end

@testset "Copying" begin
    for P in Ps
        pcpy1 = P([1,2,3,4,5], :y)
        pcpy2 = copy(pcpy1)
        @test pcpy1 == pcpy2
    end
end

@testset "GCD" begin
    for P in Ps
        p1 = P([2.,5.,1.])
        p2 = P([1.,2.,3.])

        @test degree(gcd(p1, p2)) == 0          # no common roots
        @test degree(gcd(p1, P(5))) == 0          # ditto
        @test degree(gcd(p1, P(eps(0.)))) == 0          # ditto
        @test degree(gcd(p1, P(0))) == degree(p1) # P(0) has the roots of p1
        @test degree(gcd(p1 + p2 * 170.10734737144486, p2)) == 0          # see, c.f., #122

        p1 = fromroots(P, [1.,2.,3.])
        p2 = fromroots(P, [1.,2.,6.])
        res = roots(gcd(p1, p2))
        @test 1. ∈ res
        @test 2. ∈ res
    end

    # issue 240
    if VERSION >= v"1.2.0" # rank with keywords; require_one_based_indexing

        P = Polynomial

        a = P([0.8457170323029561, 0.47175077674705257,  0.9775441940117577]);
        b = P([0.5410010714904849, 0.533604905984294]);
        d = P([0.5490673726445683, 0.15991109487875477]);
        @test degree(gcd(a*d,b*d)) == 0
        @test degree(gcd(a*d, b*d, atol=sqrt(eps()))) > 0
        @test  degree(gcd(a*d,b*d, method=:noda_sasaki)) == degree(d)
        @test degree(gcd(a*d,b*d, method=:numerical)) == degree(d)
        l,m,n = (5,5,5) # realiable, though for larger l,m,n only **usually** correct
        u,v,w = fromroots.(rand.((l,m,n)))
        @test degree(gcd(u*v, u*w, method=:numerical)) == degree(u)

        # Example of Zeng
        x = variable(P{Float64})
        p = (x+10)*(x^9 + x^8/3 + 1)
        q = (x+10)*(x^9 + x^8/7 - 6//7)

        @test degree(gcd(p,q)) == 0
        (@test degree(gcd(p,q, method=:noda_sasaki)) == 1)
        @test degree(gcd(p,q, method=:numerical)) == 1

        # more bits don't help Euclidean
        x = variable(P{BigFloat})
        p = (x+10)*(x^9 + x^8/3 + 1)
        q = (x+10)*(x^9 + x^8/7 - 6//7)
        @test degree(gcd(p,q)) == 0

        # Test 1 of Zeng
        x =  variable(P{Float64})
        alpha(j,n) = cos(j*pi/n)
        beta(j,n) = sin(j*pi/n)
        r1, r2 = 1/2, 3/2
        U(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in 1:n)
        V(n) = prod( (x-r2*alpha(j,n))^2 + r2^2*beta(j,n)^2 for j in 1:n)
        W(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in (n+1):2n)
        for n in 2:2:20
            p = U(n) * V(n); q = U(n) * W(n)
            @test degree(gcd(p,q, method=:numerical)) == degree(U(n))
        end

        # Test 5 of Zeng
        x =  variable(P{Float64})
        for ms in ((2,1,1,0), (3,2,1,0), (4,3,2,1), (5,3,2,1), (9,6,4,2),
                   (20, 14, 10, 5), (80,60,40,20), (100,60,40,20)
                  )

            p = prod((x-i)^j for (i,j) in enumerate(ms))
            dp = derivative(p)
            @test degree(gcd(p,dp, method=:numerical)) == sum(max.(ms .- 1, 0))
        end

        # fussy pair
        x =  variable(P{Float64})
        for n in (2,5,10,20,50, 100)
            p = (x-1)^n * (x-2)^n * (x-3)
            q = (x-1) * (x-2) * (x-4)
            @test degree(gcd(p,q, method=:numerical)) == 2
        end
    end
end

@testset "Showing" begin

    p = Polynomial{Rational}([1, 4])
    @test sprint(show, p) == "Polynomial(1//1 + 4//1*x)"

    p = Polynomial{Rational{Int}}([1, 4])
    @test sprint(show, p) == "Polynomial(1//1 + 4//1*x)"

    for P in (Polynomial, ImmutablePolynomial)
        p = P([1, 2, 3])
        @test sprint(show, p) == "$P(1 + 2*x + 3*x^2)"

        p = P([1.0, 2.0, 3.0])
        @test sprint(show, p) == "$P(1.0 + 2.0*x + 3.0*x^2)"

        p = P([1 + 1im, -2im])
        @test sprint(show, p) == "$P(1 + im - 2im*x)"

        p = P([1,2,3,1])  # leading coefficient of 1
        @test repr(p) == "$P(1 + 2*x + 3*x^2 + x^3)"
        p = P([1.0, 2.0, 3.0, 1.0])
        @test repr(p) == "$P(1.0 + 2.0*x + 3.0*x^2 + 1.0*x^3)"
        p = P([1, im])
        @test repr(p) == "$P(1 + im*x)"
        p = P([1 + im, 1 - im, -1 + im, -1 - im])# minus signs
        @test repr(p) == "$P(1 + im + (1 - im)x - (1 - im)x^2 - (1 + im)x^3)"
        p = P([1.0, 0 + NaN * im, NaN, Inf, 0 - Inf * im]) # handle NaN or Inf appropriately
        @test repr(p) == "$P(1.0 + NaN*im*x + NaN*x^2 + Inf*x^3 - Inf*im*x^4)"

        p = P([1,2,3])

        @test repr("text/latex", p) == "\$1 + 2\\cdot x + 3\\cdot x^{2}\$"
        p = P([1 // 2, 2 // 3, 1])
        @test repr("text/latex", p) == "\$\\frac{1}{2} + \\frac{2}{3}\\cdot x + x^{2}\$"
        p = P([complex(1,1),complex(0,1),complex(1,0),complex(1,1)])
        @test repr("text/latex", p) == "\$1 + i + i\\cdot x + x^{2} + (1 + i)x^{3}\$"

        # customized printing with printpoly
        function printpoly_to_string(args...; kwargs...)
            buf = IOBuffer()
            printpoly(buf, args...; kwargs...)
            return String(take!(buf))
        end
        @test printpoly_to_string(P([1,2,3], "y")) == "1 + 2*y + 3*y^2"
        @test printpoly_to_string(P([1,2,3], "y"), descending_powers = true) == "3*y^2 + 2*y + 1"
        @test printpoly_to_string(P([2, 3, 1], :z), descending_powers = true, offset = -2) == "1 + 3*z^-1 + 2*z^-2"
        @test printpoly_to_string(P([-1, 0, 1], :z), offset = -1, descending_powers = true) == "z - z^-1"
        @test printpoly_to_string(P([complex(1,1),complex(1,-1)]),MIME"text/latex"()) == "1 + i + (1 - i)x"
    end
end

@testset "Plotting" begin
    p = fromroots([-1, 1]) # x^2 - 1
    r = -1.4:0.055999999999999994:1.4
    rec = apply_recipe(Dict{Symbol,Any}(), p)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "-1 + x^2")
    @test rec[1].args == (r, p.(r))

    r = -1:0.02:1
    rec = apply_recipe(Dict{Symbol,Any}(), p, -1, 1)
    @test getfield(rec[1], 1) == Dict{Symbol,Any}(:label => "-1 + x^2")
    @test rec[1].args == (r, p.(r))

    p = ChebyshevT([1,1,1])
    rec = apply_recipe(Dict{Symbol,Any}(), p)
    r = -1.0:0.02:1.0  # uses domain(p)
    @test rec[1].args == (r, p.(r))

end

@testset "Promotion"  begin

    # Test different types work together
    Ps = (Polynomial,  ImmutablePolynomial, SparsePolynomial, LaurentPolynomial)
    for P₁ in Ps
        for   P₂ in Ps
            p₁, p₂ = P₁(rand(1:5, 4)), P₂(rand(1:5, 5))
            p₁ + p₂
            p₁ * p₂

            p₁, p₂ = P₁(rand(1:5, 4)), P₂(5) # constant
            p₁ + p₂
            p₁ * p₂

            p₁, p₂ = P₁(rand(1:5, 4)), P₂(5, :y) # constant, but wrong variable
            if !(promote_type(P₁, P₂) <: Polynomial || promote_type(P₁, P₂) <: Polynomials.StandardBasisPolynomial)
                p₁ + p₂
                p₁ * p₂
            end
        end
    end

    # P{T}(vector{S}) will promote to P{T}
    for Ts in ((Int32, Int,  BigInt),
               (Int,  Rational{Int}, Float64),
               (Float32, Float64, BigFloat)
              )

        n = length(Ts)
        for i in 1:n-1
            T1,T2 = Ts[i],Ts[i+1]
            for P in (Polynomial, SparsePolynomial, LaurentPolynomial)
                p = P{T2}(T1.(rand(1:3,3)))
                @test typeof(p) == P{T2}
            end
            for P in (ImmutablePolynomial,)
                N = 3
                p = P{T2}(T1.(rand(1:3,N)))
                @test typeof(p) == P{T2,N}
            end

        end
    end

    # test P{T}(...) is P{T}
    for P in (Polynomial,   SparsePolynomial, LaurentPolynomial)
        for  T in (Int32, Int64, BigInt)
            p₁ =  P{T}(Float64.(rand(1:3,5)))
            @test typeof(p₁) == P{T} # conversion works
            @test_throws InexactError  P{T}(rand(5))
        end
    end
    for P in (ImmutablePolynomial,)
        for  T in (Int32, Int64, BigInt)
            N = 5
            p₁ =  P{T}(Float64.(rand(1:3,5)))
            @test typeof(p₁) == P{T,5} # conversion works
            @test_throws InexactError  P{T}(rand(5))
        end
    end

end