text:
export ContinousMeasurementProtocol, ContinousMeasurementProtocolParams
"""
Parameters for the MPIMeasurementProtocol
"""
Base.@kwdef mutable struct ContinousMeasurementProtocolParams <: ProtocolParams
  "Pause between measurements"
  pause::typeof(1.0u"s") = 0.2u"s"
  "If set the tx amplitude and phase will be set with control steps"
  controlTx::Bool = false
  "Sequence to measure"
  sequence::Union{Sequence, Nothing} = nothing
end
function ContinousMeasurementProtocolParams(dict::Dict, scanner::MPIScanner)
  sequence = nothing
  if haskey(dict, "sequence")
    sequence = Sequence(scanner, dict["sequence"])
    dict["sequence"] = sequence
    delete!(dict, "sequence")
  end
  params = params_from_dict(ContinousMeasurementProtocolParams, dict)
  params.sequence = sequence
  return params
end
ContinousMeasurementProtocolParams(dict::Dict) = params_from_dict(ContinousMeasurementProtocolParams, dict)

Base.@kwdef mutable struct ContinousMeasurementProtocol <: Protocol
  name::AbstractString
  description::AbstractString
  scanner::MPIScanner
  params::ContinousMeasurementProtocolParams
  biChannel::Union{BidirectionalChannel{ProtocolEvent}, Nothing} = nothing
  executeTask::Union{Task, Nothing} = nothing

  latestMeas::Array{Float32, 4} = zeros(Float32, 0, 0, 0, 0)
  stopped::Bool = false
  cancelled::Bool = false
  finishAcknowledged::Bool = false
  measuring::Bool = false
  txCont::Union{TxDAQController, Nothing} = nothing
  counter::Int64 = 0
end

function requiredDevices(protocol::ContinousMeasurementProtocol)
  result = [AbstractDAQ]
  if protocol.params.controlTx
    push!(result, TxDAQController)
  end
  return result
end

function _init(protocol::ContinousMeasurementProtocol)
  if isnothing(protocol.params.sequence)
    throw(IllegalStateException("Protocol requires a sequence"))
  end
  protocol.stopped = false
  protocol.cancelled = false
  protocol.finishAcknowledged = false
  if protocol.params.controlTx
    controllers = getDevices(protocol.scanner, TxDAQController)
    if length(controllers) > 1
      throw(IllegalStateException("Cannot unambiguously find a TxDAQController as the scanner has $(length(controllers)) of them"))
    end
    protocol.txCont = controllers[1]
    protocol.txCont.currTx = nothing
  else
    protocol.txCont = nothing
  end
  protocol.counter = 0
end

function timeEstimate(protocol::ContinousMeasurementProtocol)
  est = "∞"
  return est
end

function enterExecute(protocol::ContinousMeasurementProtocol)
  protocol.stopped = false
  protocol.cancelled = false
  protocol.finishAcknowledged = false
end

function _execute(protocol::ContinousMeasurementProtocol)
  @info "Measurement protocol started"

  while true
    @info "Taking new measurement"
    measurement(protocol)
    protocol.counter+=1

    measPauseOver = false
    waitTimer = Timer(ustrip(u"s", protocol.params.pause))
    @async begin
      wait(waitTimer)
      measPauseOver = true
    end

    notifiedStop = false
    while !measPauseOver || protocol.stopped
      handleEvents(protocol)
      if !notifiedStop && protocol.stopped
        put!(protocol.biChannel, OperationSuccessfulEvent(StopEvent()))
        notifiedStop = true
      end
      if notifiedStop && !protocol.stopped
        put!(protocol.biChannel, OperationSuccessfulEvent(ResumeEvent()))
        notifiedStop = false
      end
      protocol.cancelled && throw(CancelException())
      sleep(0.05)
    end
    close(waitTimer)
  end

  @info "Protocol finished."
  close(protocol.biChannel)
end

function measurement(protocol::ContinousMeasurementProtocol)
  # Start async measurement
  protocol.measuring = true
  measState = asyncMeasurement(protocol)
  producer = measState.producer
  consumer = measState.consumer

  # Handle events
  while !istaskdone(consumer)
    handleEvents(protocol)
    protocol.cancelled && throw(CancelException())
    sleep(0.05)
  end
  protocol.measuring = false

  # Check tasks
  ex = nothing
  if Base.istaskfailed(producer)
    @error "Producer failed"
    stack = Base.catch_stack(producer)[1]
    @error stack[1]
    @error stacktrace(stack[2])
    ex = stack[1]
  end
  if  Base.istaskfailed(consumer)
    @error "Consumer failed"
    stack = Base.catch_stack(consumer)[1]
    @error stack[1]
    @error stacktrace(stack[2])
    if isnothing(ex)
      ex = stack[1]
    end
  end
  if !isnothing(ex)
    throw(ErrorException("Measurement failed, see logged exceptions and stacktraces"))
  end

  protocol.latestMeas = copy(measState.buffer)
end

function asyncMeasurement(protocol::ContinousMeasurementProtocol)
  scanner = protocol.scanner
  sequence = protocol.params.sequence
  prepareAsyncMeasurement(scanner, sequence)
  if protocol.params.controlTx
    controlTx(protocol.txCont, sequence, protocol.txCont.currTx)
  end
  scanner.seqMeasState.producer = @tspawnat scanner.generalParams.producerThreadID asyncProducer(scanner.seqMeasState.channel, scanner, sequence, prepTx = !protocol.params.controlTx)
  bind(scanner.seqMeasState.channel, scanner.seqMeasState.producer)
  scanner.seqMeasState.consumer = @tspawnat scanner.generalParams.consumerThreadID asyncConsumer(scanner.seqMeasState.channel, scanner)
  return scanner.seqMeasState
end

function cleanup(protocol::ContinousMeasurementProtocol)
  # NOP
end

function stop(protocol::ContinousMeasurementProtocol)
  protocol.stopped = true
end

function resume(protocol::ContinousMeasurementProtocol)
  protocol.stopped = false
end

function cancel(protocol::ContinousMeasurementProtocol)
  protocol.cancelled = true
  #put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent()))
  # TODO stopTx and reconnect for pipeline and so on
end

function handleEvent(protocol::ContinousMeasurementProtocol, event::DataQueryEvent)
  data = protocol.latestMeas
  put!(protocol.biChannel, DataAnswerEvent(data, event))
end

function handleEvent(protocol::ContinousMeasurementProtocol, event::ProgressQueryEvent)
  reply = ProgressEvent(protocol.counter, 0, "Measurements", event)
  put!(protocol.biChannel, reply)
end

handleEvent(protocol::ContinousMeasurementProtocol, event::FinishedAckEvent) = protocol.finishAcknowledged = true


text:
"""

The chemostat steady state impose a contraint over the exchanges:

``intake <= c/ξ``

where ``c`` is the concentration of the metabolite in the medium and ``ξ`` is the cell-specific dilution rate.
Together with the thermodinamic bound we have

stst_bound(c, ξ, b) = sign(b) * min(c/ ξ, abs(b))

"""
stst_bound(c, ξ, b) = sign(b) * min(c/ ξ, abs(b))

ext: jl

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat


text:
tpf = DataFrame(
  incl = 30:5:60,
  tp = 5.2:0.1:5.8
)

actual = DataFrame(
  incls = [40.0, 42.3, 44.0, 44.65, 45.7, 47.25, 49.5],
  atp = [2.55, 2.65, 0.4, 1.15, 1.18, 2.15, 2.35],
  rpm = [120, 120, 60, 70, 80, 80, 140]
)

exp_02 = @pgf TikzPicture(
    Axis({
	    	height = "6cm", width = "8cm",
	    	xmin = 0, xmax = 60,
	    	ymin = 0, ymax = 6,
	        axis_y_line = "left",
	        xlabel = "Inclination [°]",
      		ylabel = "TP [°/100ft]",
      		title = "Actual vs. Isotropic TP",
      		legend_pos = "south west",
		},
	    Plot({ color="blue", mark = "." },
	      Table(tpf.incl[1:7], tpf.tp[1:7])),
	    LegendEntry("TP, wob=30")),
   Axis({
	    	height = "6cm", width = "8cm",
	    	xmin = 0, xmax = 60,
	    	ymin = 0, ymax = 150,
	    	hide_x_axis,
	        axis_y_line = "right",
      		ylabel = "RPM [revs/min]",
	    },
	    Plot({ color="red", mark = "." },
	      Table(actual.incls, actual.rpm)))
)

print_tex(exp_02)

savefig("exp_02", exp_02, @__DIR__)


text:
using Test
using Dorothy

@DorothyAll()

@testset "PDB" begin

	@testset "Read" begin
		m1 = readf("$(Dorothy.datapath)/1BTL.pdb")
		@test length(m1) == 2236
		(a,b,c), (α,β,γ) = pbcgeometry(m1.header.cell)
		@test [a,b,c] ≈ [43.1, 64.4, 91.2]
		@test m1.header.title ==
				"CRYSTAL STRUCTURE OF ESCHERICHIA COLI TEM1 BETA-LACTAMASE" *
				" AT 1.8 ANGSTROMS RESOLUTION"
		@test m1.resids[1] == 26
		@test m1.resnames[1] == "HIS"
		@test m1.names[1] == "N"
		@test m1.ids[1] == 1
		@test m1.chainids[1] == "A"
		@test m1.occupancies[1] == 1.0
		@test m1.bfactors[1] == 14.51
		@test m1.R[end] ≈ [-12.176, -11.865, 48.715]
		@test ! haskey(m1, :V)
	end

	@testset "Write" begin
		m1 = readf("$(Dorothy.datapath)/1BTL.pdb")
		io = IOBuffer()
		write(specify(io, "structure/x-pdb"), m1)
		seekstart(io)
		m2 = read(specify(io, "structure/x-pdb"))
		@test length(m2) == 2236
		(a,b,c), (α,β,γ) = pbcgeometry(m2.header.cell)
		@test [a,b,c] ≈ [43.1, 64.4, 91.2]
		@test m2.header.title ==
				"CRYSTAL STRUCTURE OF ESCHERICHIA COLI TEM1 BETA-LACTAMASE" *
				" AT 1.8 ANGSTROMS RESOLUTION"
		@test m2.resids[1] == 26
		@test m2.resnames[1] == "HIS"
		@test m2.names[1] == "N"
		@test m2.ids[1] == 1
		@test m2.chainids[1] == "A"
		@test m2.occupancies[1] == 1.0
		@test m2.bfactors[1] == 14.51
		@test m2.R[end] ≈ [-12.176, -11.865, 48.715]
		@test ! haskey(m2, :V)
	end

	@testset "Convert" begin
		m1 = readf("$(Dorothy.datapath)/1BTL.pdb")
		io = IOBuffer()
		write(specify(io, "structure/x-gro"), m1)
		seekstart(io)
		m2 = read(specify(io, "structure/x-gro"))
		@test length(m2) == 2236
		(a,b,c), (α,β,γ) = pbcgeometry(m2.header.cell)
		@test [a,b,c] ≈ [43.1, 64.4, 91.2]
		@test m2.header.title ==
				"CRYSTAL STRUCTURE OF ESCHERICHIA COLI TEM1 BETA-LACTAMASE" *
				" AT 1.8 ANGSTROMS RESOLUTION"
		@test m2.resids[1] == 26
		@test m2.resnames[1] == "HIS"
		@test m2.names[1] == "N"
		@test m2.ids[1] == 1
		@test m2.R[end] ≈ [-12.18, -11.87, 48.72]
		@test ! haskey(m2, :V)
	end

end # @testset


text:
using JOLI
using WAVEFORM
using PyPlot
using LinearAlgebra
using Flux, CUDA
using NeuralOperators

n = 101*[1;1]
d = 10.0*[1;1]
o = 0.0*[1;1]
v0 = 2000
freq = 20.0
λ = v0/freq

t0 = 0.0
f0 = 20.0
unit = "m/s"
xsrc = [(n.*d)[1]/2]
ysrc = [0.0]
zsrc = [100.0]
xrec = range(0.0,stop=(n.*d)[1]::Float64,length=n[1])
yrec = [0.0]
zrec = [100.0]
freqs = [freq]
nsrc = length(xsrc)*length(ysrc)*length(zsrc)
nfreq = length(freqs)
model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec)

comp_n = n
comp_d = d
comp_o = o
npml = round(Int,λ/minimum(comp_d))*[1 1; 1 1];
scheme = WAVEFORM.helm2d_chen9p
cut_pml = true
implicit_matrix = true
srcfreqmask = trues(nsrc,nfreq)
misfit = WAVEFORM.least_squares
lsopts = LinSolveOpts(solver=:lufact);
opts = PDEopts{Int64,Float64}(scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

v = v0*ones(Float64,n...)
v[51:end, :] .= 3000
v = vec(v)

(H,comp_grid,T,DT_adj,P) = helmholtz_system(v,model,freq,opts)
nt = comp_grid.comp_n
q = zeros(eltype(H),tuple(nt...))
q[div.(nt,2)...] = 1.0
q = vec(q)
u = H\q
u = reshape(u,nt...)

clim = 10
imshow(real(u),vmin=-clim,vmax=clim)
savefig("waveform")

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

model = Chain(
    Dense(2, 64),
    FourierOperator(64=>64, (16, 16), relu),
    FourierOperator(64=>64, (16, 16), relu),
    FourierOperator(64=>64, (16, 16), relu),
    FourierOperator(64=>64, (16, 16)),
    Dense(64, 128, relu),
    Dense(128, 2)
)

# Grid
n = [121, 121]
#x = range(0.0, stop=(n[1] - 1).*d[1], length=n[1]) .* ones(1, n[2])
#y = range(0.0, stop=(n[2] - 1).*d[2], length=n[2])' .* ones(n[1], 1)

#x = reshape(x, 1, n[1], n[2], 1)
#y = reshape(y, 1, n[1], n[2], 1)
q = reshape(q, 121, 121)
x = reshape(real(q), 1, n[1], n[2], 1)
y = reshape(real(q), 1, n[1], n[2], 1)

X = cat(x, y, dims=1)
ur = reshape(real(u), 1, n[1], n[2], 1)
ui = reshape(imag(u), 1, n[1], n[2], 1)
Y = cat(ur, ui, dims=1)
Y_ = model(X)

loss(x, y) = sum(abs2, y .- model(x)) / size(x)[end]
opt = Flux.Optimiser(WeightDecay(1f-4), Flux.ADAM(1f-3))

X = X |> gpu
Y = Y |> gpu
model = model |> gpu

epochs = 500
ps = params(model)
loss_history = []
for i=1:epochs
    gs = Flux.gradient(() -> loss(X, Y), ps)
    Flux.update!(opt, ps, gs)
    global loss_history = push!(loss_history, loss(X, Y))
    print("Iter ", i, "; ", loss_history[end], "\n")
end

Y_ = model(X)
imshow(Y_[1,:,:,1], vmin=-10, vmax=10); savefig("_Y_REAL")
imshow(Y[1,:,:,1], vmin=-10, vmax=10); savefig("Y_REAL")
imshow(Y_[1,:,:,1] - Y[1,:,:,1], vmin=-1, vmax=1); savefig("Diff")

u0 = zeros(ComplexF64, length(q))
u1, res1 = FGMRES(H, q, u0; maxiter=20, precond=identity, tol=1e-2)

u0 = pre(q)#zeros(ComplexF64, length(q))
u2, res2 = FGMRES(H, q, u0; maxiter=20, precond=identity, tol=1e-2)

u0 = zeros(ComplexF64, length(q))
u3, res3 = FGMRES(H, q, u0; maxiter=20, precond=pre, tol=1e-2)

u4 = pre(u3)
subplot(2,2,1)
imshow(reshape(real(uT), 121, 121), vmin=-10, vmax=10)
subplot(2,2,2)
imshow(reshape(real(u1), 121, 121), vmin=-10, vmax=10)
subplot(2,2,3)
imshow(reshape(real(u2), 121, 121), vmin=-10, vmax=10)
subplot(2,2,4)
imshow(reshape(real(u4), 121, 121), vmin=-10, vmax=10)
savefig("comparison")

plot(res1)
plot(res2)
plot(res3)
legend(["Orig", "Apriori", "Precon"])
savefig("convergence")

u = reshape(u, 121, 121)
imshow(real.(u)); savefig("u")

function pre(x)
    xr = reshape(real(x), 1, 121, 121, 1)
    xi = reshape(imag(x), 1, 121, 121, 1)
    x = cat(xr, xi, dims=1)
    y = model(x)
    return vec(y[1,:,:,1] + im*y[2,:,:,1])
end


text:
#EasyPlot: docstring definitions
#-------------------------------------------------------------------------------

#==attributes()
===============================================================================#

@doc """`EasyPlot.attributes(...)`

Create a `AttributeChangeSpec[]` vector describing which plot attributes are
to be overwritten.

## Construct using function call:
```julia-repl
lin_log = EasyPlot.attributes(
    xyaxes = set(xscale=:lin, yscale=:log)
)
```

## Construct using function `cons()` interface:
```julia-repl
lin_log = cons(:a,
    xyaxes = set(xscale=:lin, yscale=:log)
)
```

## Alternative accepted symbols for `cons()` interface:
```julia-repl
cons(:attribute_list, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
cons(:attr, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
cons(:a, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
```

# See also:
[`cons`](@ref),
[`EasyPlot.PlotCollection`](@ref), [`EasyPlot.Plot`](@ref), [`EasyPlot.Waveform`](@ref).
[`EasyPlot.TextAnnotation`](@ref), [`EasyPlot.vmarker`](@ref), [`EasyPlot.hmarker`](@ref),
""" attributes

@doc """`EasyPlot.Plot([attr]; title="TITLE", [attr_kwargs])`

Create an object that describes a plot.

Alternative construction using `cons()` interface:
```julia-repl
plot = cons(:plot, title="Plot Title",
	xyaxes=set(xscale=:lin, yscale=:lin)
	labels=set(xaxis="X-Axis Label", yaxis="Y-Axis Label")
)
```

The following are supported attribute keyword arguments (`attr_kwargs`):
 - `xyaxes=...`, `xaxis=...`, `nstrips=...`, `ystrip=...`, `xfolded=...`, `labels...`

Alternatively, `attr_kwargs` can be bundled into `AttributeChangeSpec` vectors with the [`EasyPlot.attributes`](@ref)
function call or its `cons(:a, ...)` alias. `attr=AttributeChangeSpec` vectors can then be passed to the constructor.

```julia-repl
linlin = cons(:a, xyaxes=set(xscale=:lin, yscale=:lin))
alabels = cons(:a, labels=set(xaxis="X-Axis Label", yaxis="Y-Axis Label"))

plot = cons(:plot, linlin, alabels, title="Plot Title")
```

# Supported attributes

## X/Y Axes (`xyaxes=`): Produce single y-strip.

The following example lists which attributes of x/y axes values can be overwritten:

TODO: List valid scales
```julia-repl
lin_log = cons(:a,
    xyaxes = set(xscale=:lin, yscale=:log)
)

axesattr = cons(:a,
    xyaxes = set(xscale=:lin, yscale=:log, xmin=3, xmax=5, ymin=0, ymax=10)
)
```

## X-Axis (`xaxis=`)

The following example shows how x-axis attributes can be overwritten:

TODO: List valid scales
```julia-repl
logscale = cons(:a, xaxis = set(scale=:log))

xattr = cons(:a,
	xaxis = set(scale=:log, min=-1, max=4, label="Frequency (Hz)")
)
```

## Multi-Y-Strips (`nstrips=`, `ystrip=`)
The following example lists which y-strip attributes can be overwritten:

TODO: List valid scales
```julia-repl
ystripattributes = cons(:a,
    nstrips = 3,
    ystrip1 = set(min=0, max=.5, scale=:lin, axislabel="Voltage", striplabel="Strip 1 Label"),
    #Leave defaults for strip #2
    ystrip3 = set(min=1, max=1000, scale=:log, striplabel="Strip 3 Label"),
)

#=Note:
kwargs `ystripi` are only supporte from i=1...9. To set properties for i>9, you
can set the base `ystrip=` parameter (only 1 can be set per call to attributes()):
=#
ystripattributes2 = cons(:a,
    nstrips = 15,
    ystrip = set(12, min=0, max=.5, scale=:lin, axislabel="Voltage", striplabel="Strip 12 Label"),
)
```

## Folded x-axis (`xfolded=`)
Folded x=axis overlays (folds) all x-data onto same x-values.
**Used to create eye diagrams.**

The following example lists which folded x-axis attributes can be overwritten:
```julia-repl
foldedxaxis = cons(:a,
    xfolded = set(1e-9, xstart=5e-9, xmin=-1e9, xmax=1e-9)
)
```

## Labels (`labels=`)
The following example lists which labels can be overwritten:
```julia-repl
axislabels = cons(:a,
    labels = set(title="Plot Title", xaxis="X-Axis Label", yaxis="Y-Axis Label")
)
```

# See also:
[`cons`](@ref), [`EasyPlot.attributes`](@ref),
[`EasyPlot.PlotCollection`](@ref), [`EasyPlot.Waveform`](@ref).
""" Plot

@doc """`EasyPlot.Waveform(data::DataMD, [attr]; [attr_kwargs], label="idstr")`

Creates a waveform from `data` from a multi-dimensional (`DataMD`) dataset.

Alternative construction using `cons()` interface:

```julia-repl
wfrm = cons(:wfrm, dataset1, line=set(style=:solid, color=:red), label="Dataset 1")
```

The following are supported attribute keyword arguments (`attr_kwargs`):
 - `line=...`, `glyph=...`

Alternatively, `attr_kwargs` can be bundled into `AttributeChangeSpec` vectors with the [`EasyPlot.attributes`](@ref)
function call or its `cons(:a, ...)` alias. `attr=AttributeChangeSpec` vectors can then be passed to the constructor.

```julia-repl
dfltline = cons(:a, line=set(style=:solid, color=:red))
dfltglyph = cons(:a, glyph=set(shape=:square, size=3))

wfrm = cons(:wfrm, dataset1, dfltline, dfltglyph, label="Dataset 1")
```

# Supported attributes

## Line Attributes (`line=`)
The following example lists which line attributes can be overwritten:

TODO: List valid styles
```julia-repl
lineattributes = cons(:a,
    line = set(style=:dash, width=3, color=:red)
)
```

## Glyph/Marker Attributes (`glyph=`)
The following example lists which glyph attributes can be overwritten:\r
TODO: List valid shapes
```julia-repl
glyphattributes = cons(:a,
    line = set(shape=:circle, size=3, color=:red, fillcolor=:white)
)
```

# See also:
[`cons`](@ref), [`EasyPlot.attributes`](@ref),
[`EasyPlot.PlotCollection`](@ref), [`EasyPlot.Plot`](@ref).
""" Waveform

@doc """`TextAnnotation("Some text", [attr]; [attr_kwargs])`

Text annotation object to overlay on plot.

Alternative construction using `cons()` interface:
```julia-repl
annot = cons(:atext, "Some text", x=1, y=2, angle=0, align=:bl, strip=1)
```

The following are supported attribute keyword arguments (`attr_kwargs`):
 - `prop=...`, `offset=...`, `reloffset=...`

Alternatively, `attr_kwargs` can be bundled into `AttributeChangeSpec` vectors with the [`EasyPlot.attributes`](@ref)
function call or its `cons(:a, ...)` alias. `attr=AttributeChangeSpec` vectors can then be passed to the constructor.

```julia-repl
atextprop = cons(:a, #Construct `attributes` object
	prop=set(angle=0, align=:bl, strip=1),
	offset=set(x=-1, y=-1), #Device units
	reloffset=set(x=.5, y=.5), #Normalized to [0,1] graph bounds; depends on zoom level
)

annot = cons(:atext, "Some text", atextprop, x=1, y=2)
```

# Supported attributes

#TODO: Describe in more detail

# See also:
[`cons`](@ref) [`attributes`](@ref)
""" TextAnnotation

@doc """`cons(obj_symbol, args...; kwargs...)`

"Construct" interface created to minimize number of exported symbols.

## Construct `AttributeChangeSpec` vectors (call `EasyPlot.attributes()`):
```julia-repl
cons(:attribute_list, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
cons(:attr, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
cons(:a, args...; kwargs...) = EasyPlot.attributes(args...; kwargs...)
```

## Construct plot objects:
```julia-repl
cons(:plot_collection, args...; kwargs...) = EasyPlot.PlotCollection(args...; kwargs...)
cons(:plotcoll, args...; kwargs...) = EasyPlot.PlotCollection(args...; kwargs...)
cons(:plot, args...; kwargs...) = EasyPlot.Plot(args...; kwargs...)
cons(:wfrm, args...; kwargs...) = EasyPlot.Waveform(args...; kwargs...)
cons(:fldaxis, args...; kwargs...) = EasyPlot.FoldedAxis(args...; kwargs...)
```

## Construct annotation objects:
```julia-repl
cons(:atext, args...; kwargs...) = EasyPlot.TextAnnotation(args...; kwargs...)
cons(:vmarker, args...; kwargs...) = EasyPlot.vmarker(args...; kwargs...)
cons(:hmarker, args...; kwargs...) = EasyPlot.hmarker(args...; kwargs...)
```

# See also:
[`EasyPlot.attributes`](@ref),
[`EasyPlot.PlotCollection`](@ref), [`EasyPlot.Plot`](@ref), [`EasyPlot.Waveform`](@ref).
[`EasyPlot.TextAnnotation`](@ref), [`EasyPlot.vmarker`](@ref), [`EasyPlot.hmarker`](@ref),
""" cons
#Last line


text:
## Dyestuff data from lme4
const ds = DataFrame(Yield = [1545.,1440.,1440.,1520.,1580.,1540.,1555.,1490.,1560.,1495.,
                              1595.,1550.,1605.,1510.,1560.,1445.,1440.,1595.,1465.,1545.,
                              1595.,1630.,1515.,1635.,1625.,1520.,1455.,1450.,1480.,1445.],
                     Batch = pool(rep('A':'F',1,5)))

## Dyestuff2 data from lme4
const ds2 = DataFrame(Yield = [7.298,3.846,2.434,9.566,7.99,5.22,6.556,0.608,11.788,-0.892,
                               0.11, 10.386,13.434,5.51,8.166,2.212,4.852,7.092,9.288,4.98,
                               0.282,9.014,4.458,9.446,7.198,1.722,4.782,8.106,0.758,3.758],
                      Batch = pool(rep('A':'F',1,5)))

## sleepstudy data from lme4

const slp = DataFrame(Reaction =
                      [249.56,258.7047,250.8006,321.4398,356.8519,414.6901,382.2038,
                       290.1486,430.5853,466.3535,222.7339,205.2658,202.9778,204.707,
                       207.7161,215.9618,213.6303,217.7272,224.2957,237.3142,199.0539,
                       194.3322,234.32,232.8416,229.3074,220.4579,235.4208,255.7511,
                       261.0125,247.5153,321.5426,300.4002,283.8565,285.133,285.7973,
                       297.5855,280.2396,318.2613,305.3495,354.0487,287.6079,285.0,
                       301.8206,320.1153,316.2773,293.3187,290.075,334.8177,293.7469,
                       371.5811,234.8606,242.8118,272.9613,309.7688,317.4629,309.9976,
                       454.1619,346.8311,330.3003,253.8644,283.8424,289.555,276.7693,
                       299.8097,297.171,338.1665,332.0265,348.8399,333.36,362.0428,
                       265.4731,276.2012,243.3647,254.6723,279.0244,284.1912,305.5248,
                       331.5229,335.7469,377.299,241.6083,273.9472,254.4907,270.8021,
                       251.4519,254.6362,245.4523,235.311,235.7541,237.2466,312.3666,
                       313.8058,291.6112,346.1222,365.7324,391.8385,404.2601,416.6923,
                       455.8643,458.9167,236.1032,230.3167,238.9256,254.922,250.7103,
                       269.7744,281.5648,308.102,336.2806,351.6451,256.2968,243.4543,
                       256.2046,255.5271,268.9165,329.7247,379.4445,362.9184,394.4872,
                       389.0527,250.5265,300.0576,269.8939,280.5891,271.8274,304.6336,
                       287.7466,266.5955,321.5418,347.5655,221.6771,298.1939,326.8785,
                       346.8555,348.7402,352.8287,354.4266,360.4326,375.6406,388.5417,
                       271.9235,268.4369,257.2424,277.6566,314.8222,317.2135,298.1353,
                       348.1229,340.28,366.5131,225.264,234.5235,238.9008,240.473,
                       267.5373,344.1937,281.1481,347.5855,365.163,372.2288,269.8804,
                       272.4428,277.8989,281.7895,279.1705,284.512,259.2658,304.6306,
                       350.7807,369.4692,269.4117,273.474,297.5968,310.6316,287.1726,
                       329.6076,334.4818,343.2199,369.1417,364.1236],
                      Days = rep(0:9,18),
                      Subject = pool(rep(1:18,1,10)))

const bb = rep('A':'J',1,6)
const cc = rep('a':'c',10,2)

## Pastes data from the lme4 package
const psts = DataFrame(Strength = [62.8,62.6,60.1,62.3,62.7,63.1,60.0,61.4,57.5,56.9,61.1,58.9,
                                   58.7,57.5,63.9,63.1,65.4,63.7,57.1,56.4,56.9,58.6,64.7,64.5,
                                   55.1,55.1,54.7,54.2,58.8,57.5,63.4,64.9,59.3,58.1,60.5,60.0,
                                   62.5,62.6,61.0,58.7,56.9,57.7,59.2,59.4,65.2,66.0,64.8,64.1,
                                   54.8,54.8,64.0,64.0,57.7,56.8,58.3,59.3,59.2,59.2,58.9,56.6],
                       Batch = pool(bb),
                       Cask = pool(cc),
                       Sample = pool(ASCIIString[string(b)*string(c) for (b,c) in zip(bb,cc)]))

## Penicillin data from the lme4 package
const pen = DataFrame(Diameter = [27,23,26,23,23,21,27,23,26,23,23,21,25,21,25,24,24,20,26,23,25,
                                  23,23,20,25,22,26,22,23,20,24,22,25,23,22,19,24,20,23,21,22,19,
                                  26,22,26,24,24,21,24,21,24,22,22,20,24,21,24,23,22,19,26,23,26,
                                  24,24,21,25,22,26,24,24,20,26,24,26,24,25,22,26,23,26,23,23,20,
                                  26,23,25,24,24,22,25,22,25,23,23,20,25,21,24,23,23,20,25,22,24,
                                  23,23,19,24,21,23,21,21,19,26,23,26,24,24,21,25,21,24,22,22,18,
                                  25,22,25,22,22,20,24,21,24,22,24,19,24,21,24,22,21,18],
                      Plate = pool(rep('a':'x',1,6)),
                      Sample = pool(rep('A':'F',24,1)))

## InstEval data from the lme4 package


text:
module PlasmaParameters

using Unitful
using ..Units
using ..Constants: q, me, ϵ0

function electron_cyclotron_frequency(B::Unitful.BField)
    return uconvert(rad/s, q*B/me)
end

function electron_cyclotron_frequency(B::U) where U
    B_units = isa(B, Unitful.Quantity) ? B : B*T
    return ustrip(electron_cyclotron_frequency(B_units))
end

function ion_cyclotron_frequency(B::Unitful.BField, Zi::Int, mi::Unitful.Mass)
    return uconvert(rad/s, Zi*q*B/mi)
end

function ion_cyclotron_frequency(B, Zi, mi)
    B_units = isa(B, Unitful.Quantity) ? B : B*T
    mi_units = isa(mi, Unitful.Quantity) ? mi : mi*u
    return ustrip(ion_cyclotron_frequency(B_units, Zi, mi_units))
end

export electron_cyclotron_frequency, ion_cyclotron_frequency

function electron_plasma_frequency(ne::NumberDensity)
    return uconvert(rad/s,sqrt((ne*q^2)/(me*ϵ0)))
end

function electron_plasma_frequency(ne)
    ne_units = isa(ne, Unitful.Quantity) ? ne : ne*cm^-3
    return ustrip(electron_plasma_frequency(ne_units))
end

function ion_plasma_frequency(ni::NumberDensity, Zi::Int, mi::Unitful.Mass)
    return uconvert(rad/s,sqrt((ni*(Zi*q)^2)/(mi*ϵ0)))
end

function ion_plasma_frequency(ni, Zi, mi)
    ni_units = isa(ni, Unitful.Quantity) ? ni : ni*cm^-3
    mi_units = isa(mi, Unitful.Quantity) ? mi : mi*u
    return ustrip(ion_plasma_frequency(ni_units, Zi, mi_units))
end

export electron_plasma_frequency, ion_plasma_frequency

function electron_thermal_speed(Te::Unitful.Energy)
    return uconvert(m/s,sqrt(Te/me))
end

function electron_thermal_speed(Te)
    Te_units = isa(Te, Unitful.Quantity) ? Te : Te*keV
    return ustrip(electron_thermal_speed(Te_units))
end

function ion_thermal_speed(Ti::Unitful.Energy, mi::Unitful.Mass)
    return uconvert(m/s,sqrt(Ti/mi))
end

function ion_thermal_speed(Ti, mi)
    Ti_units = isa(Ti, Unitful.Quantity) ? Ti : Ti*keV
    mi_units = isa(mi, Unitful.Quantity) ? mi : mi*u
    return ustrip(ion_thermal_speed(Ti_units, mi_units))
end

export electron_thermal_speed, ion_thermal_speed

function electron_gyro_radius(Te::Unitful.Energy, B::Unitful.BField)
    return uconvert(m, sqrt(2*me*Te)/(q*B))
end

function electron_gyro_radius(Te, B)
    Te_units = isa(Te, Unitful.Quantity) ? Te : Te*keV
    B_units = isa(B, Unitful.Quantity) ? B : B*T
    return ustrip(electron_gyro_radius(Te_units,B_units))
end

function ion_gyro_radius(Ti::Unitful.Energy, B::Unitful.BField, Zi::Int, mi::Unitful.Mass)
    return uconvert(m, sqrt(2*mi*Ti)/(Zi*q*B))
end

function ion_gyro_radius(Ti, B, Zi, mi)
    Ti_units = isa(Ti, Unitful.Quantity) ? Ti : Ti*keV
    B_units = isa(B, Unitful.Quantity) ? B : B*T
    mi_units = isa(mi, Unitful.Quantity) ? mi : mi*u
    return ustrip(ion_gyro_radius(Ti_units,B_units, Zi, mi_units))
end

export electron_gyro_radius, ion_gyro_radius

end #module


text:
"""
data_harvest(BaseDirectory; ptable=false)
"""
function data_harvest(BaseDirectory; ptable = false)
    cvals = [5, 99, 1]
    itvals = [10, 10, 30]
    half = "no"
    levels = 5
    for ic = 1:3
        c = cvals[ic]
        maxit = itvals[ic]
        workingdir = string(BaseDirectory, "c=", string(c))
        if ptable == false
            figtitle = string("Figure", " ", string(ic))
            figure(ic)
        end
        cd(workingdir)
        plotnsold(half, c, maxit, levels; table = ptable)
        #   plotnsold(half, c, maxit, levels; bigtitle=figtitle)
        cd("../..")
    end
    half = "yes"
    for ic = 1:3
        c = cvals[ic]
        maxit = itvals[ic]
        workingdir = string(BaseDirectory, "c=", string(c), "/Float16")
        cd(workingdir)
        if ptable == false
            figure(ic + 3)
            figtitle = string("Figure", " ", string(ic + 3))
        end
        plotnsold(half, c, maxit, levels; table = ptable)
        #   plotnsold(half, c, maxit, levels; bigtitle=figtitle)
        cd("../../..")
    end

end

max_stars_repo_name: ctkelley/MPResults


text:
"""
    ms_mle_step(hmm::AbstractHMM{F}, observations) where F

Perform one step of the EM (Baum-Welch) algorithm.

# Example
```julia
hmm, log_likelihood = ms_mle_step(hmm, observations)
```
"""

function bw_step(hmm::BHMM, observations, obs_lengths)
    lls = bw_llhs(hmm, observations)
    log_α = messages_forwards_log(hmm.a, hmm.A, lls, obs_lengths)
    log_β = messages_backwards_log(hmm.A, lls, obs_lengths)
    log_A = log.(hmm.A)

    K,Tmaxplus1,O = size(lls) #the last T value is the 0 end marker of the longest T

    #transforms to cut down log_ξ, log_γ assignment times
    lls = permutedims(lls, [2,3,1]) # from (K,T,O) to (T,O,K)
    log_α = permutedims(log_α, [2,3,1])
    log_β = permutedims(log_β, [2,3,1])

    # E-step
    log_ξ = fill(-Inf, Tmaxplus1,O,K,K)
    log_γ = fill(-Inf, Tmaxplus1,O,K)
    log_pobs = zeros(O)

   @inbounds for o = 1:O
        log_pobs[o] = logsumexp(lps.(log_α[1,o,:], log_β[1,o,:]))
    end

    Threads.@threads for idx in [(i,j,o) for i=1:K, j=1:K, o=1:O]
        i,j,o = idx[1],idx[2],idx[3]
        obsl = obs_lengths[o]
        @inbounds for t = 1:obsl-1 #log_ξ & log_γ calculated to T-1 for each o
           log_ξ[t,o,i,j] = lps(log_α[t,o,i], log_A[i,j], log_β[t+1,o,j], lls[t+1,o,j], -log_pobs[o])
           log_γ[t,o,i] = lps(log_α[t,o,i], log_β[t,o,i], -log_pobs[o])
        end
        t=obsl #log_ξ @ T = 0
        log_ξ[t,o,i,j] = 0
        log_γ[t,o,i] = lps(log_α[t,o,i], log_β[t,o,i], -log_pobs[o])
    end

    ξ = similar(log_ξ)
    ξ .= exp.(log_ξ)
    ∑k_ξ = sum(ξ, dims=[3,4])
    nan_mask = ∑k_ξ .== 0
    ∑k_ξ[nan_mask] .= Inf #prevent NaNs in dummy renorm
    ξ  ./=  ∑k_ξ #dummy renorm across K to keep numerical creep from causing isprobvec to fail on new new_A during hmm creation

    γ = similar(log_γ)
    γ .= exp.(log_γ)
    ∑k_γ = sum(γ, dims=3)
    ∑k_γ[nan_mask[:,:,:]] .= Inf #prevent NaNs in dummy renorm
    γ ./= ∑k_γ #dummy renorm

    # M-step
    new_A = zeros(K,K)
    for i=1:K, j=1:K
        ∑otξ_vec = zeros(O)
        ∑otγ_vec = zeros(O)
        Threads.@threads for o in 1:O
            ∑otξ_vec[o] = sum(ξ[1:obs_lengths[o]-1,o,i,j])
            ∑otγ_vec[o] = sum(γ[1:obs_lengths[o]-1,o,i])
        end
        new_A[i,j] = sum(∑otξ_vec) / sum(∑otγ_vec)
    end
    new_A ./= sum(new_A, dims=2) #dummy renorm
    new_a = (sum(γ[1,:,:], dims=1)./sum(sum(γ[1,:,:], dims=1)))[1,:]
    new_a ./= sum(new_a) #dummy renorm

    obs_mask = .!nan_mask
    obs_collection = observations[obs_mask[:,:]]

    B = Categorical[]
    @inbounds for (i, d) in enumerate(hmm.B)
        γ_d = γ[:,:,i]
        push!(B, t_categorical_mle(Categorical, d.support[end], obs_collection, γ_d[obs_mask[:,:]]))
    end

    return typeof(hmm)(new_a, new_A, B), lps(log_pobs)
end
                function bw_llhs(hmm, observations)
                    lls = zeros(length(hmm.B), size(observations)...)
                    Threads.@threads for d in 1:length(hmm.B)
                        lls[d,:,:] = logpdf.(hmm.B[d], observations)
                    end
                    return lls
                end
                #Multisequence competent log implementations of forward and backwards algos
                function messages_forwards_log(init_distn, trans_matrix, log_likelihoods, obs_lengths)
                    log_alphas = zeros(size(log_likelihoods))
                    log_trans_matrix = log.(trans_matrix)
                    log_alphas[:,1,:] = log.(init_distn) .+ log_likelihoods[:,1,:]
                    Threads.@threads for o in 1:size(log_likelihoods)[3]
                        @inbounds for t in 2:obs_lengths[o], j in 1:size(log_likelihoods)[1]
                            log_alphas[j,t,o]=logsumexp(log_alphas[:,t-1,o] .+ log_trans_matrix[:,j]) + log_likelihoods[j,t,o]
                        end
                    end
                    return log_alphas
                end

                function messages_backwards_log(trans_matrix, log_likelihoods, obs_lengths)
                    log_betas = zeros(size(log_likelihoods))
                    log_trans_matrix = log.(trans_matrix)
                    Threads.@threads for o in 1:size(log_likelihoods)[3]
                        @inbounds for t in obs_lengths[o]-1:-1:1
                            tmp = view(log_betas, :, t+1, o) .+ view(log_likelihoods, :, t+1, o)
                            @inbounds for i in 1:size(log_likelihoods)[1]
                                log_betas[i,t,o] = logsumexp(view(log_trans_matrix, i,:) .+ tmp)
                            end
                        end
                    end
                    return log_betas
                end
                #subfunc derived from Distributions.jl categorical.jl fit_mle, threaded
                function t_categorical_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}, w::AbstractArray{F}) where T<:Integer where F<:AbstractFloat
                    Categorical(t_pnormalize!(t_add_categorical_counts!(zeros(k), x, w)))
                end

                t_pnormalize!(v::Vector) = (v ./= sum(v); v)

                function t_add_categorical_counts!(h::Vector{F}, x::AbstractArray{T}, w::AbstractArray{F}) where T<:Integer where F<:AbstractFloat
                    n = length(x)
                    if n != length(w)
                        throw(DimensionMismatch("Inconsistent array lengths."))
                    end
                    hlock = ReentrantLock()
                    Threads.@threads for i=1:n
                        @inbounds xi = x[i]
                        @inbounds wi = w[i]
                        lock(hlock)
                        h[xi] += wi   # cannot use @inbounds, as no guarantee that x[i] is in bound
                        unlock(hlock)
                    end
                    return h
                end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2, 3, 4, 2, true, false, false, true, true, false),),
    ((2, 3, 4, 2, true, false, false, true, true, true),),
    ((2, 3, 4, 2, false, true, false, true, true, true),),
    ((2, 3, 4, 2, false, false, true, true, true, true),),
    ((2, 3, 4, 2, true, false, false, false, false, true),),
    ((2, 3, 4, 2, false, true, false, false, false, true),),
    ]
insts["fast"] = [
    ((3, 5, 7, 2, true, false, false, true, true, false),),
    ((3, 5, 7, 2, true, false, false, true, true, true),),
    ((3, 5, 7, 2, false, true, false, true, true, true),),
    ((3, 5, 7, 2, false, false, true, true, true, true),),
    ((3, 5, 7, 2, true, false, false, false, false, true),),
    ((3, 5, 7, 2, false, true, false, false, false, true),),
    ((8, 15, 25, 5, true, false, false, true, true, false),),
    ((8, 15, 25, 5, true, false, false, true, true, true),),
    ((8, 15, 25, 5, false, true, false, true, true, true),),
    ((8, 15, 25, 5, false, false, true, true, true, true),),
    ((8, 15, 25, 5, true, false, false, false, false, true),),
    ((8, 15, 25, 5, false, true, false, false, false, true),),
    ((25, 75, 125, 5, true, false, false, true, true, false),),
    ((25, 75, 125, 5, true, false, false, true, true, true),),
    ((25, 75, 125, 5, false, true, false, true, true, true),),
    ((25, 75, 125, 5, false, false, true, true, true, true),),
    ((25, 75, 125, 5, true, false, false, false, false, true),),
    ((25, 75, 125, 5, false, true, false, false, false, true),),
    ]
return (DOptimalDesignNative, insts)


text:
module VSCodeDebugger

import Sockets

include("../../../error_handler.jl")

include("../../CodeTracking/src/CodeTracking.jl")
include("../../JSON/src/JSON.jl")

module JuliaInterpreter
    using ..CodeTracking

    include("../../JuliaInterpreter/src/packagedef.jl")
end

module JSONRPC
    import ..JSON
    import UUIDs

    include("../../JSONRPC/src/packagedef.jl")
end

module DebugAdapter
    import ..JuliaInterpreter
    import ..JSON
    import ..JSONRPC
    import ..JSONRPC: @dict_readable, Outbound

    include("../../DebugAdapter/src/packagedef.jl")
end

function startdebugger()
    pipenames = DebugAdapter.clean_up_ARGS_in_launch_mode()
    try
        @debug "Trying to connect to debug adapter."
        socket = Sockets.connect(pipenames[1])
        try
            DebugAdapter.startdebug(socket, (err, bt)->global_err_handler(err, bt, pipenames[2], "Debugger"))
        finally
            close(socket)
        end
    catch err
        global_err_handler(err, catch_backtrace(), pipenames[2], "Debugger")
    end
end

end


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

@doc raw"""TokenRequestStatus is the result of a token request.

    IoK8sApiAuthenticationV1TokenRequestStatus(;
        expirationTimestamp=nothing,
        token=nothing,
    )

    - expirationTimestamp::IoK8sApimachineryPkgApisMetaV1Time : ExpirationTimestamp is the time of expiration of the returned token.
    - token::String : Token is the opaque bearer token.
"""
mutable struct IoK8sApiAuthenticationV1TokenRequestStatus <: SwaggerModel
    expirationTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: expirationTimestamp
    token::Any # spec type: Union{ Nothing, String } # spec name: token

    function IoK8sApiAuthenticationV1TokenRequestStatus(;expirationTimestamp=nothing, token=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1TokenRequestStatus, Symbol("expirationTimestamp"), expirationTimestamp)
        setfield!(o, Symbol("expirationTimestamp"), expirationTimestamp)
        validate_property(IoK8sApiAuthenticationV1TokenRequestStatus, Symbol("token"), token)
        setfield!(o, Symbol("token"), token)
        o
    end
end # type IoK8sApiAuthenticationV1TokenRequestStatus

const _property_map_IoK8sApiAuthenticationV1TokenRequestStatus = Dict{Symbol,Symbol}(Symbol("expirationTimestamp")=>Symbol("expirationTimestamp"), Symbol("token")=>Symbol("token"))
const _property_types_IoK8sApiAuthenticationV1TokenRequestStatus = Dict{Symbol,String}(Symbol("expirationTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("token")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1TokenRequestStatus }) = collect(keys(_property_map_IoK8sApiAuthenticationV1TokenRequestStatus))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1TokenRequestStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1TokenRequestStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1TokenRequestStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1TokenRequestStatus[property_name]

function check_required(o::IoK8sApiAuthenticationV1TokenRequestStatus)
    (getproperty(o, Symbol("expirationTimestamp")) === nothing) && (return false)
    (getproperty(o, Symbol("token")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1TokenRequestStatus }, name::Symbol, val)
end


text:
# ## Using Unfold.jl from Python
# it is straight forward to call Unfold from Python using `JuliaCall`.
# Please find a notebook [here](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/src/HowTo/pyjulia_unfold.ipynb)

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
@testset "polynomial" begin
    rng = MersenneTwister(123456)
    x = rand(rng)*2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    c = randn(rng)
    @testset "LinearKernel" begin
        k = LinearKernel()
        @test kappa(k,x) ≈ x
        @test k(v1,v2) ≈ dot(v1,v2)
        @test kappa(LinearKernel(),x) == kappa(k,x)
        @test metric(LinearKernel()) == KernelFunctions.DotProduct()
        @test metric(LinearKernel(c=2.0)) == KernelFunctions.DotProduct()
        @test repr(k) == "Linear Kernel (c = 0.0)"

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x->LinearKernel(c=x[1]), [c])
        test_params(LinearKernel(; c=c), ([c],))
    end
    @testset "PolynomialKernel" begin
        k = PolynomialKernel()
        @test kappa(k,x) ≈ x^2
        @test k(v1,v2) ≈ dot(v1,v2)^2
        @test kappa(PolynomialKernel(),x) == kappa(k,x)
        @test repr(k) == "Polynomial Kernel (c = 0.0, d = 2.0)"

        # Coherence tests.
        @test kappa(PolynomialKernel(d=1.0,c=c),x) ≈ kappa(LinearKernel(c=c),x)
        @test metric(PolynomialKernel()) == KernelFunctions.DotProduct()
        @test metric(PolynomialKernel(d=3.0)) == KernelFunctions.DotProduct()
        @test metric(PolynomialKernel(d=3.0,c=2.0)) == KernelFunctions.DotProduct()

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        # test_ADs(x->PolynomialKernel(d=x[1], c=x[2]),[2.0,  c])
        @test_broken "All, because of the power"
        test_params(PolynomialKernel(; d=x, c=c), ([x], [c]))
    end
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

function dijkstra_cut(A1::AbstractArray, A2::AbstractArray, dim::Integer)
  # permute dimensions so that the algorithm is
  # the same for cuts in x, y and z directions
  B = abs.(A1 - A2)
  if dim == 1
    B = permutedims(B, [1,2,3])
  elseif dim == 2
    B = permutedims(B, [2,1,3])
  elseif dim == 3
    B = permutedims(B, [3,2,1])
  end

  mx, my, mz = size(B)

  # accumulation cube and overlap mask
  E = zeros(size(B)); M = falses(size(B))

  # pad accumulation cube with +inf
  Epad(i,j,k) = all(0 .< [i,j,k] .≤ [mx,my,mz]) ? E[i,j,k] : Inf

  # forward accumulation along 3D cube
  E[:,:,1] = B[:,:,1]
  for k=2:mz, i=1:mx, j=1:my
    square = [Epad(i-1 , j-1 , k-1) , Epad(i-1 , j , k-1) , Epad(i-1 , j+1 , k-1),
              Epad(i   , j-1 , k-1) , Epad(i   , j , k-1) , Epad(i   , j+1 , k-1),
              Epad(i+1 , j-1 , k-1) , Epad(i+1 , j , k-1) , Epad(i+1 , j+1 , k-1)]
    E[i,j,k] = B[i,j,k] + minimum(square)
  end

  # forward accumulation along 2D square (last slice of the cube)
  zslice = view(E,:,:,mz)
  for j=2:my
    zslice[1,j] += minimum(zslice[1:2,j-1])
    for i=2:mx-1
      zslice[i,j] += minimum(zslice[i-1:i+1,j-1])
    end
    zslice[mx,j] += minimum(zslice[mx-1:mx,j-1])
  end

  # backward search along last slice
  idx = argmin(zslice[:,my])
  mslice = view(M,:,:,mz)
  mslice[1:idx,my] = trues(idx)
  idxvec = zeros(Int, my); idxvec[my] = idx # keep track of indexes
  width = isodd(mx) ? mx+1 : mx # avoid zig-zag artifact
  for j=my-1:-1:1
    for i=1:width
      if idx < mx && minimum(zslice[max(idx-1,1):idx+1,j]) == zslice[idx+1,j]
        idx += 1
      elseif idx > 1 && zslice[idx-1,j] ≤ zslice[idx,j]
        idx -= 1
      end
    end
    mslice[1:idx,j] .= true
    idxvec[j] = idx
  end

  # backward search along cube
  for j=1:my
    yslice = view(E,:,j,:)
    idx = idxvec[j]
    for k=mz-1:-1:1
      for i=1:width
        if idx < mx && minimum(yslice[max(idx-1,1):idx+1,k]) == yslice[idx+1,k]
          idx += 1
        elseif idx > 1 && yslice[idx-1,k] ≤ yslice[idx,k]
          idx -= 1
        end
      end
      M[1:idx,j,k] .= true
    end
  end

  # permute back to original shape
  if dim == 1
    M = permutedims(M, [1,2,3])
  elseif dim == 2
    M = permutedims(M, [2,1,3])
  elseif dim == 3
    M = permutedims(M, [3,2,1])
  end

  M
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(StateSpace.call, (Type{StateSpace.LinearGaussianSSM}, Function, Function, Function, Function, Function,))
    precompile(StateSpace.confirm_matrix_sizes, (Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2},))
    precompile(StateSpace.ispossemidef, (Array{Float64, 2},))
end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey


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

using Markdown
using InteractiveUtils

# ╔═╡ 92be606a-0998-11eb-053e-cb42246e0598
using Pkg, DrWatson

# ╔═╡ 92be94e0-0998-11eb-03fc-4b2293b5e32a
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 8cc034b8-0998-11eb-2c71-578b14536d5d
md"## Clip-05-28-44.jl"

# ╔═╡ 92c90f38-0998-11eb-086f-4b4900891b2f
md"### snippets 5.28 - 5.32"

# ╔═╡ 94042790-099b-11eb-1a63-6dc1af913def
md"##### Notice that `missing` values are dropped!"

# ╔═╡ 92d87388-0998-11eb-1881-11abc4118c34
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; missingstring = "NA")
	dropmissing!(df, [:kcal_per_g, :neocortex_perc, :mass])
	df.K = zscore(df.kcal_per_g)
	df.N = zscore(df.neocortex_perc)
	df.M = zscore(log.(df.mass))
end;

# ╔═╡ 55142898-0999-11eb-1cff-6d82e287b64a
md"### snippet 5.33"

# ╔═╡ fdde27d4-099a-11eb-2156-9b15bba277f2
@model function m5_5_draft(N, K)
    a ~ Normal(0, 1)
    bN ~ Normal(0, 1)
    σ ~ Exponential(1)
    μ = lin(a, bN, N)
    K ~ MvNormal(μ, σ)
end

# ╔═╡ fdde6780-099a-11eb-0fd8-b597b7b34d15
begin
	m5_5_draftt = m5_5_draft(df.N, df.K)
	quap5_5_draftt = quap(m5_5_draftt)
end

# ╔═╡ e17b5820-099d-11eb-242b-c1f8c9e2fc14
md"### snippet 5.34"

# ╔═╡ d03a085c-09a2-11eb-3ba3-d18810ca8ab5
xseq = -2:0.1:2

# ╔═╡ e4d38704-099d-11eb-0c30-9d7f0556f5b1
begin
	prior = sample(m5_5_draftt, Prior(), 50) |> DataFrame
	μ = lin(prior.a', prior.bN', xseq)

	fig1 = plot(
		xlab="neocortex percent (std)", ylab="kcal per g (std)",
		title="a ~ Normal(0, 1)\nbN ~ Normal(0, 1)",
		legend = false, ylims = extrema(xseq))
	for c in eachcol(μ)
		plot!(xseq, c, color = :black, alpha = 0.3)
	end
	plot!()
end

# ╔═╡ e4d3b814-099d-11eb-2d2d-bbf106ef33b6
md"### snippet 5.35"

# ╔═╡ e4d468b6-099d-11eb-2b83-3b2d265a4fad
@model function m5_5(N, K)
    a ~ Normal(0, 0.2)
    bN ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bN, N)
    K ~ MvNormal(μ, σ)
end

# ╔═╡ e4ea98f4-099d-11eb-0ca5-bd416c74dbab
begin
	m5_5t = m5_5(df.N, df.K)
	chns5_5t = sample(m5_5t, NUTS(0.65), 2000)
end

# ╔═╡ 82cf498e-09a3-11eb-1edf-4f32d126c20e
begin
	post = sample(m5_5_draftt, NUTS(0.65), 200) |> DataFrame
	μ_post = lin(post.a', post.bN', xseq)

	fig2 = plot(
		xlab="neocortex percent (std)", ylab="kcal per g (std)",
		title="a ~ Normal(0, 0.2)\nbN ~ Normal(0, 0.5)",
		legend = false, ylims = extrema(xseq))
	for c in eachcol(μ_post)
		plot!(xseq, c, color = :black, alpha = 0.3)
	end
	plot(fig1, fig2, layout=(1, 2))
end

# ╔═╡ e4ebf03e-099d-11eb-16e3-93fea8ed8f8a
md"### snippet 5.36"

# ╔═╡ e4f60d74-099d-11eb-324e-31bef9d5181f
begin
	post5_5t = DataFrame(Array(chns5_5t), names(chns5_5t, [:parameters]))
	Text(precis(post5_5t; io=String))
end

# ╔═╡ 26000e2c-099f-11eb-05cc-bb0bede678d1
begin
	quap5_5t = quap(m5_5t)
	Text(precis(quap5_5t; io=String))
end

# ╔═╡ 4010b8b4-09a1-11eb-173a-2bb20afc882d
md"### snippet 5.37"

# ╔═╡ 44a91bfa-09a1-11eb-3bd9-af4110eab915
begin
	post1 = DataFrame(rand(quap5_5t.distr, 1000)', quap5_5t.params)
	μ1 = lin(post1.a', post1.bN', xseq) |> meanlowerupper

	fig3 = plot(xlab="neocortex percent (std)", ylab="kcal per gram (std)")
	scatter!(df.N, df.K, legend = false)
	plot!(xseq, μ1.mean, ribbon = (μ1.mean .- μ1.lower, μ1.upper .- μ1.mean))
end

# ╔═╡ 44a9500c-09a1-11eb-36e3-750c47541694
md"### snippet 5.38"

# ╔═╡ 44aa2266-09a1-11eb-3d17-5d3f2da14256
@model function m5_6(M, K)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bM, M)
    K ~ MvNormal(μ, σ)
end

# ╔═╡ 44bf965a-09a1-11eb-2d3b-95b5e58e1edd
begin
	m5_6t = m5_6(df.M, df.K)
	Text(precis(quap(m5_6t); io=String))
end

# ╔═╡ 753a9688-09a4-11eb-1bce-654aeb39c454
begin
	quap5_6t = quap(m5_6t)
	post2 = DataFrame(rand(quap5_6t.distr, 1000)', quap5_6t.params)
	μ2 = lin(post2.a', post2.bM', xseq) |> meanlowerupper

	fig4 = plot(xlab="log body mass (std)", ylab="kcal per gram (std)")
	scatter!(df.N, df.K, legend = false)
	plot!(xseq, μ2.mean, ribbon = (μ2.mean .- μ2.lower, μ2.upper .- μ2.mean))
end

# ╔═╡ 44c06a94-09a1-11eb-243e-e94c4ea7b657
md"### snippet 5.39"

# ╔═╡ 44d48ee8-09a1-11eb-0dd6-c5aa78fc84ec
@model function m5_7(M, N, K)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bN ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, bM, M, bN, N)
    K ~ MvNormal(μ, σ)
end

# ╔═╡ 44d79f14-09a1-11eb-2122-916596c22454
begin
	m5_7t = m5_7(df.M, df.N, df.K)
	quap5_7t = quap(m5_7t)
	post5_7t = DataFrame(rand(quap5_7t.distr, 1000)', quap5_7t.params)
	Text(precis(post5_7t; io=String))
end

# ╔═╡ 2aa59598-09a2-11eb-0854-5b7b1a9aa1b2
md"### snippet 5.40"

# ╔═╡ 99bfbc84-09a1-11eb-01a6-35c26e589d7b
begin
	quap_res, fig = plotcoef([m5_5t, m5_6t, m5_7t], [:a, :bN, :bM, :σ])
	plot(fig)
end

# ╔═╡ 1982c47a-09a2-11eb-357c-1b16ba7aacb1
quap_res

# ╔═╡ 59c2eef4-09a5-11eb-2a39-b75a6de102ad
md"### snippet 5.41"

# ╔═╡ 5e7f6ba2-09a5-11eb-1168-7bbd0f34e93b
begin
	mu2 = lin(post5_7t.a', post5_7t.bN', xseq, post5_7t.bM', zeros(length(xseq))) |> meanlowerupper
	fig5 = plot(xlab="neocortex percent (std)", ylab="kcal per gram (std)",
		title="Counterfactual (M=0)", leg=false)
	plot!(xseq, mu2.mean, ribbon = (mu2.upper .- mu2.mean, mu2.mean .- mu2.lower))

	mu3 = lin(post5_7t.a', post5_7t.bM', xseq, post5_7t.bN', zeros(length(xseq))) |> meanlowerupper
	fig6 = plot(xlab="neocortex percent (std)", ylab="log body mass (std)",
		title="Counterfactual (N=0)", leg=false)
	plot!(xseq, mu3.mean, ribbon = (mu3.upper .- mu3.mean, mu3.mean .- mu3.lower))

	plot(fig3, fig4, fig5, fig6, layout=(2,2))
end

# ╔═╡ 405a080e-0995-11eb-01c0-3ff35fe8c925
md"### snippet 5.42"

# ╔═╡ 99aa57b4-09a7-11eb-2080-01a888bd0999
n = 100

# ╔═╡ 99aa96fa-09a7-11eb-17c6-bb516e055371
begin
	# M -> K <- N
	# M -> N
	M1 = randn(n)
	N1 = rand.(Normal.(M1))
	K1 = rand.(Normal.(N1 .- M1))
	d_sim1 = DataFrame((; K1, N1, M1))
	Text(precis(d_sim1; io=String))
end

# ╔═╡ 99b7614e-09a7-11eb-21b4-716994e5df9d
md"### snippet 5.43"

# ╔═╡ 99ba1a76-09a7-11eb-1729-339be9f12454
begin
	# M -> K <- N
	# N -> M
	N2 = randn(n)
	M2 = rand.(Normal.(N2))
	K2 = rand.(Normal.(N2 .- M2))
	d_sim2 = DataFrame((; K2, N2, M2))
	Text(precis(d_sim2; io=String))
end

# ╔═╡ 99d1457a-09a7-11eb-18c8-bb41b89014f3
begin
	# M -> K <- N
	# M -> U <- N
	U3 = randn(n)
	N3 = rand.(Normal.(U3))
	M3 = rand.(Normal.(U3))
	K3 = rand.(Normal.(N3 .- M3))
	d_sim3 = DataFrame((; K3, N3, M3))
	Text(precis(d_sim3; io=String))
end

# ╔═╡ 99d510ba-09a7-11eb-0fcf-97d4242bb925
md"### snippet 5.44"

# ╔═╡ 99e1a37a-09a7-11eb-3839-87157105995e
md"##### TODO working with DAGs."

# ╔═╡ 99eecdfe-09a7-11eb-12c8-fd7d822403c1
md"## End of clip-05-28-44.jl"

# ╔═╡ Cell order:
# ╟─8cc034b8-0998-11eb-2c71-578b14536d5d
# ╠═92be606a-0998-11eb-053e-cb42246e0598
# ╠═92be94e0-0998-11eb-03fc-4b2293b5e32a
# ╟─92c90f38-0998-11eb-086f-4b4900891b2f
# ╟─94042790-099b-11eb-1a63-6dc1af913def
# ╠═92d87388-0998-11eb-1881-11abc4118c34
# ╟─55142898-0999-11eb-1cff-6d82e287b64a
# ╠═fdde27d4-099a-11eb-2156-9b15bba277f2
# ╠═fdde6780-099a-11eb-0fd8-b597b7b34d15
# ╟─e17b5820-099d-11eb-242b-c1f8c9e2fc14
# ╠═d03a085c-09a2-11eb-3ba3-d18810ca8ab5
# ╠═e4d38704-099d-11eb-0c30-9d7f0556f5b1
# ╟─e4d3b814-099d-11eb-2d2d-bbf106ef33b6
# ╠═e4d468b6-099d-11eb-2b83-3b2d265a4fad
# ╠═e4ea98f4-099d-11eb-0ca5-bd416c74dbab
# ╠═82cf498e-09a3-11eb-1edf-4f32d126c20e
# ╟─e4ebf03e-099d-11eb-16e3-93fea8ed8f8a
# ╠═e4f60d74-099d-11eb-324e-31bef9d5181f
# ╠═26000e2c-099f-11eb-05cc-bb0bede678d1
# ╟─4010b8b4-09a1-11eb-173a-2bb20afc882d
# ╠═44a91bfa-09a1-11eb-3bd9-af4110eab915
# ╟─44a9500c-09a1-11eb-36e3-750c47541694
# ╠═44aa2266-09a1-11eb-3d17-5d3f2da14256
# ╠═44bf965a-09a1-11eb-2d3b-95b5e58e1edd
# ╠═753a9688-09a4-11eb-1bce-654aeb39c454
# ╟─44c06a94-09a1-11eb-243e-e94c4ea7b657
# ╠═44d48ee8-09a1-11eb-0dd6-c5aa78fc84ec
# ╠═44d79f14-09a1-11eb-2122-916596c22454
# ╟─2aa59598-09a2-11eb-0854-5b7b1a9aa1b2
# ╠═99bfbc84-09a1-11eb-01a6-35c26e589d7b
# ╠═1982c47a-09a2-11eb-357c-1b16ba7aacb1
# ╟─59c2eef4-09a5-11eb-2a39-b75a6de102ad
# ╠═5e7f6ba2-09a5-11eb-1168-7bbd0f34e93b
# ╟─405a080e-0995-11eb-01c0-3ff35fe8c925
# ╠═99aa57b4-09a7-11eb-2080-01a888bd0999
# ╠═99aa96fa-09a7-11eb-17c6-bb516e055371
# ╟─99b7614e-09a7-11eb-21b4-716994e5df9d
# ╠═99ba1a76-09a7-11eb-1729-339be9f12454
# ╠═99d1457a-09a7-11eb-18c8-bb41b89014f3
# ╟─99d510ba-09a7-11eb-0fcf-97d4242bb925
# ╟─99e1a37a-09a7-11eb-3839-87157105995e
# ╟─99eecdfe-09a7-11eb-12c8-fd7d822403c1


text:
function solve()
    s = readline()
    if s[1] == s[2] && s[2] == s[3]
        "Won"
    else
        "Lost"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#=
A fifo with a read and write buffer where data can be pre-loaded and recorded.

Useful for testing individual cores without needing to invoke an entire
simulation.
=#
@with_kw mutable struct TestFifo{T} <: AbstractFifo
    read_data       :: Vector{T} = T[]
    read_pointer    :: Int       = 1
    write_data      :: Vector{T} = T[]

    # Indicate a write is pending for the next "clock cycle"
    write_pending :: Bool = false
    next_write_data :: T  = zero(T)

    # Indicate need to increment read pointer next clock cycle.
    read_increment :: Bool = false
end

isreadready(t::TestFifo) = t.read_pointer <= length(t.read_data)
peek(t::TestFifo) = t.read_data[t.read_pointer]

function read(t::TestFifo)
    t.read_increment = true
    return peek(t)
end

function readupdate!(t::TestFifo)
    if t.read_increment
        t.read_pointer += 1
        t.read_increment = false
    end
    return nothing
end

iswriteready(t::TestFifo) = true

function write(t::TestFifo, data)
    t.write_pending = true
    t.next_write_data = data
    return nothing
end

function writeupdate!(t::TestFifo)
    if t.write_pending
        push!(t.write_data, t.next_write_data)
        t.write_pending = false
    end
    return nothing
end


text:
"""
    resample(rng, x[, log_weights], ndraws)

Draw `ndraws` samples from `x`, with replacement.

If `log_weights` is provided, perform Pareto smoothed importance resampling.
"""
function resample(rng, x, log_ratios, ndraws)
    result = PSIS.psis(log_ratios)
    weights = result.weights
    pweights = StatsBase.ProbabilityWeights(weights, one(eltype(weights)))
    return StatsBase.sample(rng, x, pweights, ndraws; replace=true)
end
resample(rng, x, ndraws) = StatsBase.sample(rng, x, ndraws; replace=true)


text:
using Test
using LinearAlgebra
using FlightMechanics
using FlightMechanics.Aircrafts
using FlightMechanics.Models

@testset "aerodynamics utils" begin
    ac = F16()

    @test typeof(ac)<:Aircraft
    @test get_name(ac) == "F16"
    @test typeof(ac.aerodynamics) == F16Aerodynamics

    import FlightMechanics.Aircrafts: get_interp_idx
    k, l, da = get_interp_idx(2.5, 0.2, -1, 8)
    @test (k, l) == (3, 4)
    @test isapprox(da, 0.5)

    k, l, da = get_interp_idx(1, 0.2, -1, 8)
    @test (k, l) == (3, 4)
    @test isapprox(da, 0.2)

    k, l, da = get_interp_idx(-1, 0.2, -1, 8)
    @test (k, l) == (2, 3)
    @test isapprox(da, 0.8)

    k, l, da = get_interp_idx(0.0, 0.2, -1, 8)
    @test (k, l) == (3, 3)
    @test isapprox(da, 0.0)

    k, l, da = get_interp_idx(-10.0, 0.2, -1, 8)
    @test (k, l) == (2, 1)
    @test isapprox(da, -1.0)

    k, l, da = get_interp_idx(45.0, 0.2, -1, 8)
    @test (k, l) == (11, 12)
    @test isapprox(da, 1.0)

    import FlightMechanics.Aircrafts: interp1d, CXq_data
    r = interp1d(0.0, 0.2, -1, 8, CXq_data)
    @test isapprox(r, 0.308)

    r = interp1d(45, 0.2, -1, 8, CXq_data)
    @test isapprox(r, 1.210)

    r = interp1d(-10, 0.2, -1, 8, CXq_data)
    @test isapprox(r, -0.267)

    import FlightMechanics.Aircrafts: interp2d, CX_data
    r = interp2d(0.0, 0.0, 0.2, 1/12., -1, -1, 8, 1, CX_data)
    @test isapprox(r, -0.021)

    r = interp2d(45.0, 24.0, 0.2, 1/12., -1, -1, 8, 1, CX_data)
    @test isapprox(r, 0.040)

    r = interp2d(-10.0, -24.0, 0.2, 1/12., -1, -1, 8, 1, CX_data)
    @test isapprox(r, -0.099)

    import FlightMechanics.Aircrafts: interp2d2, Cl_data, Cn_data
    r = interp2d2(0.0, 0.0, .2, .2, -1, 1, 8, 5, Cl_data)
    @test isapprox(r, 0.0)

    r = interp2d2(-10.0, 0.0, .2, .2, -1, 1, 8, 5, Cl_data)
    @test isapprox(r, 0.0)

    r = interp2d2(45.0, 30.0, .2, .2, -1, 1, 8, 5, Cl_data)
    @test isapprox(r, -0.076)

    r = interp2d2(0.0, 0.0, .2, .2, -1, 1, 8, 5, Cn_data)
    @test isapprox(r, 0.0)

    r = interp2d2(-10.0, 0.0, .2, .2, -1, 1, 8, 5, Cn_data)
    @test isapprox(r, 0.0)

    r = interp2d2(45.0, 30.0, .2, .2, -1, 1, 8, 5, Cn_data)
    @test isapprox(r, -0.001)
end

@testset "aerodynamics" begin
    ac = F16()

    # TEST against Morelli
    input_data = [
        # tas    α    β    p    q    r    de    da    dr    xcg
        500    0.5  -0.2  0.7  -0.8  0.9  20    -15  -20    0.35;
        500    0.5  -0.2  0.0   0.0  0.0  20    -15  -20    0.35;
        500    0.5  -0.2  0.0   0.0  0.0   0      0    0    0.35;
        500    0.0   0.0  0.0   0.0  0.0   0      0    0    0.35;
        500    0.5  -0.2 -0.8   0.0  0.0  20    -15  -20    0.35;
    ]

    # Data obatined from Morelli
    out_data = [
        -0.0666571605333333 -0.0485321333333333  -0.0196159964939023  0.0213716333333333 -0.116409952000000  0.0272818333333333;
        -0.0629333333333333 -0.0690833333333333  -0.283598396493902   0.0287903333333333 -0.163800000000000  0.0365518333333333;
        -0.0193000000000000  0.00400000000000000 -0.131598396493902   3.35999999999999e-04 -0.00860000000000000   -7.23999999999999e-04;
        -0.0210000000000000  0.0                 -0.100000000000000   0                    -0.00900000000000000   0;
        -0.0629333333333333 -0.0652865333333333 -0.283598396493902  0.0393671333333333 -0.163800000000000  0.0354574333333333;
    ]

    att = Attitude(-1, 1, -1)
    pos = EarthPosition(1000*FT2M, 900*FT2M, -10000*FT2M)
    env = Environment(pos, atmos="ISA1978", wind="NoWind", grav="const")
    fcs = get_fcs(ac)

    @testset "Case $ii" for ii=1:size(input_data, 1)
        tas, α, β, p, q, r, de, da, dr, xcg = input_data[ii, :]
        cx, cy, cz, cl, cm, cn = out_data[ii, :]

        u, v, w = wind2body(tas*FT2M, 0, 0, α*DEG2RAD, β*DEG2RAD)
        state = State(pos,
                      att,
                      [u, v, w],
                      [p, q, r],
                      [0., 0., 0.],
                      [0., 0., 0.]
                     )
        # Act directly on FCS surfaces
        set_value!(fcs.de, de*DEG2RAD)
        set_value!(fcs.da, da*DEG2RAD)
        set_value!(fcs.dr, dr*DEG2RAD)

        aerostate = AeroState(state, env)

        # Calculate aerodynamics
        aero = get_aerodynamics(ac)
        aero = calculate_aerodynamics(ac, aero, aerostate, state)

        pfm = get_body_adim_pfm(aero)

        @test isapprox(pfm.forces[1], cx)
        @test isapprox(pfm.forces[2], cy)
        @test isapprox(pfm.forces[3], cz)
        @test isapprox(pfm.moments[1], cl)
        @test isapprox(pfm.moments[2], cm)
        @test isapprox(pfm.moments[3], cn)
    end
end

@testset "aerodynamics 2" begin
    # TODO: Fix breaking tests if Morelli ac is implemented
    ac = F16()

    # TEST against Morelli
    input_data = [
        # tas    α    β    p    q    r    de    da    dr    xcg
        500    28.64788975654116  -11.459155902616466  0.7  -0.8  0.9  20    -15   20    0.35;
        500    28.64788975654116  -11.459155902616466  0.7  -0.8  0.9  20    -15  -20    0.35;
    ]

    # Results from Morelli
    out_data = [
         0.04247190703  0.297332165194  -1.661302153   0.05809649017 -0.05637730554  -0.0432532826901365;
         0.04247190703  0.182665498527  -1.661302153   0.05809649017 -0.05637730554  -0.0045883093820057;
    ]

    att = Attitude(-1, 1, -1)
    pos = EarthPosition(1000*FT2M, 900*FT2M, -10000*FT2M)
    env = Environment(pos, atmos="ISA1978", wind="NoWind", grav="const")
    fcs = get_fcs(ac)

    @testset "Case $ii" for ii=1:size(input_data, 1)
        tas, α, β, p, q, r, de, da, dr, xcg = input_data[ii, :]
        cx, cy, cz, cl, cm, cn = out_data[ii, :]

        u, v, w = wind2body(tas*FT2M, 0, 0, α*DEG2RAD, β*DEG2RAD)
        state = State(pos,
                      att,
                      [u, v, w],
                      [p, q, r],
                      [0., 0., 0.],
                      [0., 0., 0.]
                     )
        set_value!(fcs.de, de*DEG2RAD)
        set_value!(fcs.da, da*DEG2RAD)
        set_value!(fcs.dr, dr*DEG2RAD)

        aerostate = AeroState(state, env)

        # Calculate aerodynamics
        aero = get_aerodynamics(ac)
        aero = calculate_aerodynamics(ac, aero, aerostate, state)

        pfm = get_body_adim_pfm(aero)

        @test isapprox(pfm.forces[1], cx)
        @test isapprox(pfm.forces[2], cy)
        @test isapprox(pfm.forces[3], cz)
        @test_broken isapprox(pfm.moments[1], cl)  # Differences in Morelli
        @test isapprox(pfm.moments[2], cm)
        @test isapprox(pfm.moments[3], cn)
    end
end

@testset "engine" begin

    test_data = [
        # thtl (-)   alt(ft)   tas (ft/s)    thrust(N)
           0.0        10000      500       -940.891895587440;
           0.1        10000      500       4659.15413340486;
           0.2        10000      500      10259.2001623972;
           0.3        10000      500      15859.2461913895;
           0.4        10000      500      21459.2922203818;
           0.5        10000      500      27059.3382493741;
           0.6        10000      500      32659.3842783664;
           0.7        10000      500      38259.4303073587;
           0.8        10000      500      46841.5273717983;
           0.9        10000      500      62386.4689396095;
           1.0        10000      500      77931.4105074207
    ]

    ac = F16()
    eng = F16Engine()
    fcs = get_fcs(ac)
    att = Attitude(0., 0., 0.)
    pos = EarthPosition(0, 0, 0)
    env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")

    @testset "Case $ii" for ii=1:size(test_data, 1)
        pos = EarthPosition(0, 0, -test_data[ii, 2]*FT2M)
        env = calculate_environment(env, pos)
        u, v, w = wind2body(test_data[ii, 3]*FT2M, 0., 0., 0., 0.)
        set_thtl!(fcs, test_data[ii, 1])
        state = State(pos, att, [u, v, w], [0., 0., 0.], [0., 0., 0.], [0., 0., 0.])
        aerostate = AeroState(state, env)

        eng = calculate_engine(eng, fcs, aerostate, state)

        thrust = get_pfm(eng).forces[1]

        @test isapprox(thrust, test_data[ii, 4], atol=1e-3)
    end
end

@testset "trim test - level flight" begin
    # Integration test for aircraft just making sure everything runs
    ac = F16()

    fcs = get_fcs(ac)
    controls = StickPedalsLeverControls(0.0, 0.5, 0.5, 0.8)
    set_controls!(ac, controls)

    h = 0.0 * M2FT
    psi = 0.0  # rad
    gamma = 0.0
    turn_rate = 0.0

    pos = EarthPosition(0, 0, -h)
    env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")

    #  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
    #  and simulation: dynamics, controls design, and autonomous systems. John Wiley
    #  & Sons. (page 193 table 3.6-2)
    trim_test_data = [
    #   TAS   thtl    AOA     DE
    #   ft/s  unit    deg     deg
        130	  0.816   45.6    20.1 ;
        140   0.736   40.3   -1.36 ;
        150   0.619   34.6    0.173;
        170   0.464   27.2    0.621;
        640   0.23    0.742  -0.871;
        800   0.378  -0.045  -0.943;
        200   0.287   19.7    0.723;
        260   0.148   11.6   -0.09 ;
        300   0.122   8.49   -0.591;
        350   0.107   5.87   -0.539;
        400   0.108   4.16   -0.591;
        440   0.113   3.19   -0.671;
        500   0.137   2.14   -0.756;
        540   0.16    1.63   -0.798;
        600   0.2     1.04   -0.846;
        700   0.282   0.382  -0.9  ;
        ]

    # Cases with small deviations in DE
    @testset "trim tas=$(trim_test_data[ii, 1]) ft/s" for ii=1:4
        tas =  trim_test_data[ii, 1] * FT2M

        exp_thtl = trim_test_data[ii, 2]
        exp_α    = trim_test_data[ii, 3]
        exp_de   = trim_test_data[ii, 4]

        # Initial conditions for trimmer
        α0 = 0.5*exp_α*DEG2RAD
        β0 = 0.
        stick_lon0 = 0.
        thtl0 = 0.5

        controls = StickPedalsLeverControls(
            stick_lon0,
            0.5,
            0.5,
            thtl0
        )

        set_controls!(ac, controls)

        ac, aerostate, state, controls = steady_state_trim(
            ac, controls, env, tas, pos, psi, gamma, turn_rate, α0, 0.0, show_trace=false
            )

        fcs = get_fcs(ac)
        @test isapprox(ac.pfm.forces, zeros(3), atol=1e-7)
        @test isapprox(ac.pfm.moments, zeros(3), atol=1e-7)
        @test isapprox(get_tas(aerostate), tas)
        @test isapprox(get_beta(aerostate), 0.0, atol=1e-13)

        @test isapprox(fcs.thtl.value, exp_thtl, atol=1e-3)
        @test isapprox(aerostate.alpha*RAD2DEG, exp_α, atol=10.0^(-length(split(string(exp_α), ".")[2])))
        @test_broken isapprox(fcs.de.value*RAD2DEG, exp_de, atol=10.0^(-length(split(string(exp_de), ".")[2])))
    end

    # Cases with small deviations in α
    @testset "trim tas=$(trim_test_data[ii, 1]) ft/s" for ii=5:5
        tas =  trim_test_data[ii, 1] * FT2M

        exp_thtl = trim_test_data[ii, 2]
        exp_α    = trim_test_data[ii, 3]
        exp_de   = trim_test_data[ii, 4]

        # Initial conditions for trimmer
        α0 = 0.5*exp_α*DEG2RAD
        β0 = 0.
        stick_lon0 = 0.
        thtl0 = 0.5

        controls = StickPedalsLeverControls(
            stick_lon0,
            0.5,
            0.5,
            thtl0
        )

        set_controls!(ac, controls)

        ac, aerostate, state, controls = steady_state_trim(
            ac, controls, env, tas, pos, psi, gamma, turn_rate, α0, 0.0, show_trace=false
            )

        fcs = get_fcs(ac)

        @test isapprox(ac.pfm.forces, zeros(3), atol=1e-7)
        @test isapprox(ac.pfm.moments, zeros(3), atol=1e-7)
        @test isapprox(get_tas(aerostate), tas)
        @test isapprox(get_beta(aerostate), 0.0, atol=1e-13)

        @test isapprox(fcs.thtl.value, exp_thtl, atol=1e-3)
        @test_broken isapprox(aerostate.alpha*RAD2DEG, exp_α, atol=10.0^(-length(split(string(exp_α), ".")[2])))
        @test isapprox(fcs.de.value*RAD2DEG, exp_de, atol=10.0^(-length(split(string(exp_de), ".")[2])))
    end

    # Cases matching the test values
    @testset "trim tas=$(trim_test_data[ii, 1]) ft/s" for ii=6:size(trim_test_data, 1)
        tas =  trim_test_data[ii, 1] * FT2M

        exp_thtl = trim_test_data[ii, 2]
        exp_α    = trim_test_data[ii, 3]
        exp_de   = trim_test_data[ii, 4]

        # Initial conditions for trimmer
        α0 = 0.5*exp_α*DEG2RAD
        β0 = 0.
        stick_lon0 = 0.
        thtl0 = 0.5

        controls = StickPedalsLeverControls(
            stick_lon0,
            0.5,
            0.5,
            thtl0
        )

        set_controls!(ac, controls)

        ac, aerostate, state, controls = steady_state_trim(
            ac, controls, env, tas, pos, psi, gamma, turn_rate, α0, 0.0, show_trace=false
            )

        fcs = get_fcs(ac)

        @test isapprox(ac.pfm.forces, zeros(3), atol=1e-7)
        @test isapprox(ac.pfm.moments, zeros(3), atol=1e-7)
        @test isapprox(get_tas(aerostate), tas)
        @test isapprox(get_beta(aerostate), 0.0, atol=1e-13)

        @test isapprox(fcs.thtl.value, exp_thtl, atol=1e-3)
        @test isapprox(aerostate.alpha*RAD2DEG, exp_α, atol=10.0^(-length(split(string(exp_α), ".")[2])))
        @test isapprox(fcs.de.value*RAD2DEG, exp_de, atol=10.0^(-length(split(string(exp_de), ".")[2])))
    end
end

@testset "trim test - coord turn" begin
    #  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
    #  and simulation: dynamics, controls design, and autonomous systems. John Wiley
    #  & Sons. (page 192)
    pos = EarthPosition(0, 0, 0)
    env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")

    ac = F16()

    h = 0.0 * M2FT
    psi = 0.234077 # rad
    gamma = 0.0
    turn_rate = 0.3  # rad/s
    tas_fts = 502.0  # ft/s

    tas =  tas_fts * FT2M

    exp_α    =  2.392628e-1  # rad
    exp_β    =  5.061803e-4  # rad
    exp_ϕ    =  1.366289     # rad
    exp_θ    =  5.000808e-2  # rad
    exp_ψ    =  2.340769e-1  # rad
    exp_p    = -1.499617e-2  # rad/s
    exp_q    =  2.933811e-1  # rad/s
    exp_r    =  6.084932e-2  # rad/s
    exp_x    =  0.0  # ft
    exp_y    =  0.0  # ft
    exp_z    =  0.0  # ft

    exp_thtl =  0.8349601
    exp_de   = -1.481766     # deg
    exp_da   =  9.553108e-2  # deg
    exp_dr   = -4.118124e-1  # deg

    # Initial conditions for trimmer
    α0 = 0.5*exp_α*DEG2RAD
    β0 = 0.
    stick_lon0 = 0.0
    stick_lat0 = 0.0
    pedals0 = 0.0
    thtl0 = 0.5

    controls = StickPedalsLeverControls(
            stick_lon0,
            stick_lat0,
            pedals0,
            thtl0
        )

    set_controls!(ac, controls)

    ac, aerostate, state, controls = steady_state_trim(
        ac, controls, env, tas, pos, psi, gamma, turn_rate, α0, β0, show_trace=false
        )

    fcs = get_fcs(ac)

    @test isapprox(get_tas(aerostate), tas)
    @test isapprox(get_alpha(aerostate), exp_α, atol=5e-4)
    @test isapprox(get_beta(aerostate), exp_β, atol=5e-5)

    @test isapprox(get_body_ang_velocity(state), [exp_p, exp_q, exp_r], atol=5e-5)
    @test isapprox(get_euler_angles(state), [exp_ψ, exp_θ, exp_ϕ], atol=5e-5)

    # TODO: maybe these tolerances are too broad. Take into account that trim in
    # Stevens is only based in the number of iterations.
    @test isapprox(fcs.thtl.value, exp_thtl, atol=5e-4)
    @test isapprox(fcs.de.value*RAD2DEG, exp_de, rtol=1e-2)
    @test isapprox(fcs.da.value*RAD2DEG, exp_da, rtol=1e-2)
    @test isapprox(fcs.dr.value*RAD2DEG, exp_dr, rtol=1e-1)
end

@testset "model test case" begin
    #  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
    #  and simulation: dynamics, controls design, and autonomous systems. John Wiley
    #  & Sons. (page 185 table 3.5-2)

    # xcg = 0.4 * CMA
    # Note: empty CG cannot be modified as a parameter so the only option is to
    # modify the getter which has its value hardcoded
    # TODO: modify when ac structure is refactored to do this
    FlightMechanics.Aircrafts.get_empty_cg(ac::F16) = [-0.40 * get_chord(ac), 0, 0]

    # INPUTS
    # U(1) = THTL =   0.9  [-]
    # U(2) = ELEV =  20    [DEG]
    # U(3) = DAIL = -15    [DEG]
    # U(4) = RDR  = -20    [DEG]

    # STATE
    # IDX     name      X        UNITS        XDOT           XDOT MORELLI
    # 1       VT       500       [ft/s]      -75.23724       -77.57521
    # 2       ALPHA      0.5     [RAD]       -0.8813491      -0.88123
    # 3       BETA      -0.2     [RAD]       -0.4759990      -0.45276
    # 4       PHI       -1       [RAD]        2.505734        0.70000
    # 5       THETA      1       [RAD]        0.3250820       0.32508
    # 6       PSI       -1       [RAD]        2.145926        2.14593
    # 7       P          0.7     [RAD/S]     12.62679        12.91108
    # 8       Q         -0.8     [RAD/S]      0.9649671      0.97006
    # 9       R          0.9     [RAD/S]      0.5809759      -0.55450
    # 10      X       1000       [FT]       342.4439         342.44390
    # 11      Y        900       [FT]      -266.7707         -266.77068
    # 12      ALT    10000       [FT]       248.1241          248.12412
    # 13      POW       90       [%]        -58.6899         -58.6900

    # XXX: From morelli implementation
    # udot, vdot, wdot = 100.84. -206.45, -436.98  ft/s
    # thrust =    1.5912e+04  lbf
    # qbar =  219.72  psf
    # mach =  0.46436

    # Integration test for aircraft just making sure everything runs
    ac = F16()

    att = Attitude(-1, 1, -1)
    pos = EarthPosition(1000*FT2M, 900*FT2M, -10000*FT2M)
    u, v, w = wind2body(500*FT2M, 0, 0, 0.5, -0.2)
    state = State(pos,
                  att,
                  [u, v, w],
                  [0.7, -0.8, 0.9],
                  [0., 0., 0.],
                  [0., 0., 0.]
                 )

    # Set controls to obtain the following deflexions
    # fcs = get_fcs(ac)
    # set_value!(fcs.de, 20.0*DEG2RAD)
    # set_value!(fcs.da, -15.0*DEG2RAD)
    # set_value!(fcs.dr, -20*DEG2RAD)
    # set_value!(fcs.cpl, 90.0)
    controls = StickPedalsLeverControls(
        0.5 * (1 + 20 / 25),  # According to DE_MAX
        0.5 * (1 - 15 / 20),  # According to DA_MAX
        0.5 * (1 - 20 / 30),  # According to DR_MAX
        (90.0 + 117.38) / 217.38  # According to tgear
    )

    env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")
    aerostate = AeroState(state, env)
    grav = get_gravity(env)

    ac = calculate_aircraft(ac, controls, aerostate, state, grav; consume_fuel=false)
    pfm = ac.pfm
    mass_props = get_mass_props(ac)
    mass = mass_props.mass
    inertia = mass_props.inertia

    six_dof_euler_fixed_mass_ds = convert(SixDOFEulerFixedMass, state)
    x = get_x(six_dof_euler_fixed_mass_ds)
    f = get_state_equation(six_dof_euler_fixed_mass_ds)
    # Evaluate x_dot given x, u, parameters
    r = f(x, mass, inertia, pfm.forces, pfm.moments, [160.0*SLUGFT2_2_KGM2, 0.0, 0.0])
    u_dot, v_dot, w_dot, p_dot, q_dot, r_dot, ψ_dot, θ_dot, ϕ_dot, x_dot, y_dot, z_dot = r

    # u, v, w are used in this test, so they need to be transformed to VT, α, β
    VT = 500 * FT2M
    VT_dot_est, α_dot_est, β_dot_est = tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

    @testset "Stevens values" begin
        # XDOT Stevens
        # [RAD/S]
        α_dot_exp = -0.8813491
        β_dot_exp = -0.4759990
        ψ_dot_exp = 2.145926
        θ_dot_exp = 0.3250820
        ϕ_dot_exp = 2.505734
        p_dot_exp = 12.62679
        q_dot_exp = 0.9649671
        r_dot_exp = 0.5809759
        # [FT/S]
        VT_dot_exp = -75.23724
        x_dot_exp = 342.4439
        y_dot_exp = -266.7707
        z_dot_exp = 248.1241

        @test isapprox(VT_dot_est*M2FT, VT_dot_exp, atol=0.05)
        @test isapprox(α_dot_est, α_dot_exp, atol=0.0005)
        @test isapprox(β_dot_est, β_dot_exp, atol=0.06)

        @test isapprox(ψ_dot, ψ_dot_exp, atol=1e-6)  # ψ_dot [rad/s]
        @test isapprox(θ_dot, θ_dot_exp, atol=1e-7)  # θ_dot [rad/s]
        @test isapprox(ϕ_dot, ϕ_dot_exp, atol=1e-6)  # ϕ_dot [rad/s]

        @test isapprox(p_dot, p_dot_exp, atol=0.06)   # p_dot [rad/s]
        @test isapprox(q_dot, q_dot_exp, atol=0.005)  # q_dot [rad/s]
        @test isapprox(r_dot, r_dot_exp, atol=0.005)  # r_dot [rad/s]

        @test isapprox(x_dot*M2FT, x_dot_exp, atol=1e-4)   # x_dot (ft/s)
        @test isapprox(y_dot*M2FT, y_dot_exp, atol=1e-4)  # y_dot (ft/s)
        @test isapprox(-z_dot*M2FT, z_dot_exp, atol=1e-4)   # z_dot (ft/s)
    end

    @testset "Morelli" begin
        # Note that conditions are less strict for Morelli because there are
        # differences in Morelli aerodynamic data wrt Stevens.

        # XDOT MORELLI
        # [RAD/S]
        α_dot_exp = -0.88123
        β_dot_exp = -0.45276
        ψ_dot_exp = 2.145926
        θ_dot_exp =  0.32508
        ϕ_dot_exp = 2.14593
        p_dot_exp = 12.91108
        q_dot_exp = 0.97006
        r_dot_exp = -0.55450
        # [FT/S]
        VT_dot_exp = -77.57521
        x_dot_exp = 342.44390
        y_dot_exp = -266.77068
        z_dot_exp = 248.12412
        # From morelli implementation
        u_dot_exp = 100.84
        v_dot_exp = -206.45
        w_dot_exp = -436.98

        thrust_exp = 1.5912e+04
        qbar_exp =  219.72  # psf
        mach_exp =  0.46436

        thrust = get_pfm(get_propulsion(ac)).forces[1]
        @test isapprox(thrust, thrust_exp*LBF2N, atol=0.5)
        mach = get_mach(aerostate)
        @test isapprox(mach, mach_exp, atol=0.000001)
        qbar = get_qinf(aerostate)
        @test isapprox(qbar, qbar_exp*PSF2PA, atol=0.5)

        @test isapprox(u_dot, u_dot_exp*FT2M, atol=0.01)
        @test isapprox(v_dot, v_dot_exp*FT2M, atol=5)
        @test isapprox(w_dot, w_dot_exp*FT2M, atol=0.0005)

        @test_broken isapprox(VT_dot_est*M2FT, VT_dot_exp, atol=0.05)
        @test isapprox(VT_dot_est*M2FT, VT_dot_exp, atol=5)
        @test isapprox(α_dot_est, α_dot_exp, atol=0.00005)
        @test isapprox(β_dot_est, β_dot_exp, atol=0.05)

        @test isapprox(ψ_dot, ψ_dot_exp, atol=1e-6)  # ψ_dot [rad/s]
        @test_broken isapprox(r[8], θ_dot_exp, atol=1e-7)  # θ_dot [rad/s]
        @test isapprox(θ_dot, θ_dot_exp, atol=1e-5)  # θ_dot [rad/s]
        @test_broken isapprox(r[9], ϕ_dot_exp, atol=1e-6)  # ϕ_dot [rad/s]
        @test isapprox(ϕ_dot, ϕ_dot_exp, atol=0.5)  # ϕ_dot [rad/s]

        @test_broken isapprox(p_dot, p_dot_exp, atol=0.06)   # p_dot [rad/s]
        @test isapprox(p_dot, p_dot_exp, atol=0.5)   # p_dot [rad/s]
        @test isapprox(q_dot, q_dot_exp, atol=0.005)  # q_dot [rad/s]
        @test isapprox(r_dot, r_dot_exp, atol=1.5)   # r_dot [rad/s]

        @test isapprox(x_dot*M2FT, x_dot_exp, atol=1e-4)   # x_dot (ft/s)
        @test isapprox(y_dot*M2FT, y_dot_exp, atol=1e-4)  # y_dot (ft/s)
        @test isapprox(-z_dot*M2FT, z_dot_exp, atol=1e-4)   # z_dot (ft/s)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libcap"
version = v"2.27"

# Collection of sources required to build libcap-ng
sources = [
    "https://mirrors.kernel.org/pub/linux/libs/security/linux-privs/libcap2/libcap-$(version.major).$(version.minor).tar.gz" =>
    "260b549c154b07c3cdc16b9ccc93c04633c39f4fb6a4a3b8d1fa5b8a9c3f5fe8",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libcap-*/

# GNU (not musl) 64-bit archs put their loader in `/lib64`
if [[ ${target} == *gnu* ]] && [[ ${nbits} == 64 ]]; then
    LOADER_DIR=lib64
else
    LOADER_DIR=lib
fi

make -j${nproc} BUILD_CC=${BUILD_CC}
make install DESTDIR=${prefix} prefix=/ lib=${LOADER_DIR} RAISE_SETFCAP=no
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = [p for p in supported_platforms() if p isa Linux]

# The products that we will ensure are always built
products = [
    LibraryProduct("libcap", :libcap),
    FileProduct("include/sys/capability.h", :capability_h),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Metaheuristics
using Test

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


text:
module SizedStrings

import Base:
    iterate, length, getindex, lastindex, sizeof,
    ncodeunits, codeunit, isvalid, read, write

using Base: @propagate_inbounds

export SizedString, @sized_str

include("sized.jl")
include("padded.jl")

end


text:
using CSV
using DataFrames
using Dates
using Distributions
using Random

function read_in_cases(
	date,
	rng;
	apply_inc = false,
	omicron_dominant_date = nothing
)

	# read in the reff file
	case_file_name = "data/interim_linelist_"*date*".csv"
	# drop the first column
	df = CSV.read(case_file_name, DataFrame)
	# indicator for the NA dates
	is_confirmation = df.date_onset .== "NA"
	# get the confirmation dates
	confirm_dates = convert(Vector, df.date_confirmation)
	confirm_dates = Date.(confirm_dates)
	# shift them appropriately
	shape_rd = 1.28
	scale_rd = 2.31
	# sample from delay distribtuion
	rd = rand(rng, Gamma(shape_rd, scale_rd), length(confirm_dates))
	confirm_dates = confirm_dates - round.(rd) * Day(1)
	# adjust confirmation dates to get to onset
	# rd = ceil(Int, mean(Gamma(shape_rd, scale_rd)))
	# confirm_dates = confirm_dates .- rd * Day(1)

	# initialise array for complete_onset_dates
	complete_dates = deepcopy(confirm_dates)
	# fill the array with the most informative date
	complete_dates[.!is_confirmation] = Date.(convert(Vector, df.date_onset[.!is_confirmation]))
	complete_dates[is_confirmation] = confirm_dates[is_confirmation]

	# if want to apply delay, subtract an incubation period per individual
	if apply_inc
		shape_inc = 5.807  # taken from Lauer et al. 2020 #1.62/0.418
		scale_inc = 0.948  # 0.418
		inc = rand(rng, Gamma(shape_inc, scale_inc), length(complete_dates))

		# now apply different incubation period for Omicron
		apply_omicron = !isnothing(omicron_dominant_date) ? true : false
		if apply_omicron
			is_omicron = complete_dates .>= Dates.Date(omicron_dominant_date)
			shape_inc_omicron = 3.33
			scale_inc_omicron = 1.34
			inc_omicron = rand(rng, Gamma(shape_inc_omicron, scale_inc_omicron), length(complete_dates))
			# add the incubation for omicron dates in
			inc = (1-is_omicron)*inc + is_omicron*inc_omicron
		end

		complete_dates = complete_dates - round.(inc) * Day(1)
	end

	# other useful fields
	state = df.state
	import_status = convert.(Int, df.import_status .== "imported")

	# construct the full timeseries of counts
	dates_since_start = range(
		Date("2020-03-01"),
		maximum(complete_dates),
		step = Day(1)
	)

	complete_dates_local = complete_dates[import_status .== 0]
	complete_dates_import = complete_dates[import_status .== 1]
	state_local = state[import_status .== 0]
	state_import = state[import_status .== 1]

	# vectors to hold the number of cases each day
	local_cases = zeros(Int, length(dates_since_start))
	import_cases = zeros(Int, length(dates_since_start))

	# construct df to hold all the linelists
	new_df = DataFrame(
		"date" => dates_since_start
	)

	# construct df to hold all the linelists
	local_case_dict = Dict()
	import_case_dict = Dict()
	local_case_dict["date"] = dates_since_start
	import_case_dict["date"] = dates_since_start

	# loop over states and then sum the number of cases on that day
	for s in unique(state)
		# filter by state
		complete_dates_local_tmp = complete_dates_local[state_local .== s]
		complete_dates_import_tmp = complete_dates_import[state_import .== s]
		# get cases on each day
		local_cases = sum.(@views complete_dates_local_tmp .== dss for dss in dates_since_start)
		import_cases = sum.(@views complete_dates_import_tmp .== dss for dss in dates_since_start)
		# append to the df with a deepcopy to avoid a 0 d
		local_case_dict[s] = deepcopy(local_cases)
		import_case_dict[s] = deepcopy(import_cases)
		# reset the case vectors
		local_cases .= 0
		import_cases .= 0
	end

    return (local_case_dict, import_case_dict)

end


text:
########
# Set up parameters
# - ICTH = parameter for Thatcher & Jablonowski initial conditions
#   (see http://www-personal.umich.edu/~cjablono/DCMIP-2016_TestCaseDocument_10June2016.pdf)
########
parameters = (
    a = get_planet_parameter(:planet_radius),
    Ω = get_planet_parameter(:Omega),
    g = get_planet_parameter(:grav),
    κ = get_planet_parameter(:kappa_d),
    R_d = get_planet_parameter(:R_d),
    R_v = get_planet_parameter(:R_v),
    cv_d = get_planet_parameter(:cv_d),
    cv_v = get_planet_parameter(:cv_v),
    cv_l = get_planet_parameter(:cv_l),
    cv_i = get_planet_parameter(:cv_i),
    cp_d = get_planet_parameter(:cp_d),
    cp_v = get_planet_parameter(:cp_v),
    cp_l = get_planet_parameter(:cp_l),
    cp_i = get_planet_parameter(:cp_i),
    molmass_ratio = get_planet_parameter(:molmass_dryair) / get_planet_parameter(:molmass_water),
    γ = get_planet_parameter(:cp_d) / get_planet_parameter(:cv_d),
    pₒ = get_planet_parameter(:MSLP),
    pₜᵣ = get_planet_parameter(:press_triple),
    Tₜᵣ = get_planet_parameter(:T_triple),
    T_0 = get_planet_parameter(:T_0),
    LH_v0 = get_planet_parameter(:LH_v0),
    e_int_v0 = get_planet_parameter(:e_int_v0),
    e_int_i0 = get_planet_parameter(:e_int_i0),
    H = 30e3,   # domain height [m]
    k = 3.0,    # ICTH: power used for temperture field
    Γ = 0.005,  # ICTH: lapse rate [K / m]
    T_E = 300,    # ICTH: surface temperature horizontally averaged [K]
    T_P = 271.0,  # ICTH: surface temperature at the pole [K]
    b = 2.0,    # ICTH: half-width parameter
    z_t = 15e3,   # ICTH: max height of the zonal wind perturbation [m]
    λ_c = π / 9,  # ICTH: longitude of the zonal wind perturbation centerpoint
    ϕ_c = 2π / 9, # ICTH: latitude of the zonal wind perturbation centerpoint
    V_p = 1.0,    # ICTH: max amplitude of the zonal wind perturbation [m / s]
    ϕ_w = 2π / 9,   # ICTH: specific humidity latitudinal width parameter
    p_w = 3.4e4,  # ICTH: specific humidity vertical pressure width parameter
    q₀ = 0.018,  # ICTH: max specific humidity amplitude
    qₜ = 1e-12,  # ICTH: specific humidity above artificial tropopause
    Mᵥ = 0.608,  # ICTH: constant for virtual temperature conversion
    ΔT = 29.0,   # equator-pole SST difference: fixed SST runs only
    Tₘᵢₙ = 271.0,  # polar SST [K]: fixed SST runs only
    Δϕ = 26π / 180.0,  # latitudinal width (standard deviation) of the SST kernel: fixed SST runs only
    day = 86400,  # length of day [s]
    T_ref = 255,    # reference temperature [K]
    τ_precip = 100.0,  # precipitation timescale [s]
    p0 = 1e5,    # surface pressure [Pa]
    Cₑ = 0.0015, # bulk transfer coefficient for sensible heat
    Cₗ = 0.0015, # bulk transfer coefficient for latent heat
    c_o = 3.93e3,   # specific heat for ocean  [J / K / kg]
    T_h = 280.0,      # initial temperature of surface ocean layer [K]
    h_o = 100.0,      # depth of the modelled ocean mixed layer [m]
    F_sol = 1361.0,   # incoming solar TOA radiation [W / m^2]
    τ = 0.9,        # atmospheric transmissivity
    α = 0.5,        # surface albedo
    σ = 5.67e-8,    # Steffan Boltzmann constant [kg / s^3 / K^4]
    ϵ = 0.98,       # broadband emissivity / absorptivity
    F_a = 0.0,      # LW flux from the atmosphere [W / m^2]
    ρ_o = 1026.0,     # density for ocean [kg / m^3]
    Q_0 = 10.0,       # Q-flux amplitude [W / m^2]
    L_w = 16.0,       # width of the Q-flux region (radians)
    diurnal_period = 24 * 60 * 60, # idealized daily cycle period [s]
)

# Mask to pick out the coupled boundary in the MPIStateArrays (here at altitude = 0 m)
epss = sqrt(eps(Float64))
boundary_mask(𝒫, xc, yc, zc) = @. abs((xc^2 + yc^2 + zc^2)^0.5 - 𝒫.a) < epss

########
# Set up inital conditions
########

# 1. Land (ocean) initial condition
T_sfc₀(𝒫, x, y, z) = boundary_mask(𝒫, x, y, z) * 0.5 * (𝒫.T_E + 𝒫.T_P)

# 2. Atmos (single stack) initial conditions
# additional initial condition parameters
T₀(𝒫) = 0.5 * (𝒫.T_E + 𝒫.T_P)
A(𝒫) = 1.0 / 𝒫.Γ
B(𝒫) = (T₀(𝒫) - 𝒫.T_P) / T₀(𝒫) / 𝒫.T_P
C(𝒫) = 0.5 * (𝒫.k + 2) * (𝒫.T_E - 𝒫.T_P) / 𝒫.T_E / 𝒫.T_P
H(𝒫) = 𝒫.R_d * T₀(𝒫) / 𝒫.g
d_0(𝒫) = 𝒫.a / 6

# convenience functions that only depend on height
τ_z_1(𝒫, r) = exp(𝒫.Γ * (r - 𝒫.a) / T₀(𝒫))
τ_z_2(𝒫, r) = 1 - 2 * ((r - 𝒫.a) / 𝒫.b / H(𝒫))^2
τ_z_3(𝒫, r) = exp(-((r - 𝒫.a) / 𝒫.b / H(𝒫))^2)
τ_1(𝒫, r) = 1 / T₀(𝒫) * τ_z_1(𝒫, r) + B(𝒫) * τ_z_2(𝒫, r) * τ_z_3(𝒫, r)
τ_2(𝒫, r) = C(𝒫) * τ_z_2(𝒫, r) * τ_z_3(𝒫, r)
τ_int_1(𝒫, r) = A(𝒫) * (τ_z_1(𝒫, r) - 1) + B(𝒫) * (r - 𝒫.a) * τ_z_3(𝒫, r)
τ_int_2(𝒫, r) = C(𝒫) * (r - 𝒫.a) * τ_z_3(𝒫, r)
F_z(𝒫, r) = (1 - 3 * ((r - 𝒫.a) / 𝒫.z_t)^2 + 2 * ((r - 𝒫.a) / 𝒫.z_t)^3) * ((r - 𝒫.a) ≤ 𝒫.z_t)

# convenience functions that only depend on longitude and latitude
d(𝒫, λ, ϕ) = 𝒫.a * acos(sin(ϕ) * sin(𝒫.ϕ_c) + cos(ϕ) * cos(𝒫.ϕ_c) * cos(λ - 𝒫.λ_c))
c3(𝒫, λ, ϕ) = cos(π * d(𝒫, λ, ϕ) / 2 / d_0(𝒫))^3
s1(𝒫, λ, ϕ) = sin(π * d(𝒫, λ, ϕ) / 2 / d_0(𝒫))
cond(𝒫, λ, ϕ) = (0 < d(𝒫, λ, ϕ) < d_0(𝒫)) * (d(𝒫, λ, ϕ) != 𝒫.a * π)

# base-state thermodynamic variables
I_T(𝒫, ϕ, r) = (cos(ϕ) * r / 𝒫.a)^𝒫.k - 𝒫.k / (𝒫.k + 2) * (cos(ϕ) * r / 𝒫.a)^(𝒫.k + 2)
Tᵥ(𝒫, ϕ, r) = (τ_1(𝒫, r) - τ_2(𝒫, r) * I_T(𝒫, ϕ, r))^(-1) * (𝒫.a / r)^2
p(𝒫, ϕ, r) = 𝒫.pₒ * exp(-𝒫.g / 𝒫.R_d * (τ_int_1(𝒫, r) - τ_int_2(𝒫, r) * I_T(𝒫, ϕ, r)))
q(𝒫, ϕ, r) = (p(𝒫, ϕ, r) > 𝒫.p_w) ? 𝒫.q₀ * exp(-(ϕ / 𝒫.ϕ_w)^4) * exp(-((p(𝒫, ϕ, r) - 𝒫.pₒ) / 𝒫.p_w)^2) : 𝒫.qₜ

# base-state velocity variables
U(𝒫, ϕ, r) =
    𝒫.g * 𝒫.k / 𝒫.a * τ_int_2(𝒫, r) * Tᵥ(𝒫, ϕ, r) * ((cos(ϕ) * r / 𝒫.a)^(𝒫.k - 1) - (cos(ϕ) * r / 𝒫.a)^(𝒫.k + 1))
u(𝒫, ϕ, r) = -𝒫.Ω * r * cos(ϕ) + sqrt((𝒫.Ω * r * cos(ϕ))^2 + r * cos(ϕ) * U(𝒫, ϕ, r))
v(𝒫, ϕ, r) = 0.0
w(𝒫, ϕ, r) = 0.0

# velocity perturbations
δu(𝒫, λ, ϕ, r) =
    -16 * 𝒫.V_p / 3 / sqrt(3) *
    F_z(𝒫, r) *
    c3(𝒫, λ, ϕ) *
    s1(𝒫, λ, ϕ) *
    (-sin(𝒫.ϕ_c) * cos(ϕ) + cos(𝒫.ϕ_c) * sin(ϕ) * cos(λ - 𝒫.λ_c)) / sin(d(𝒫, λ, ϕ) / 𝒫.a) * cond(𝒫, λ, ϕ)
δv(𝒫, λ, ϕ, r) =
    16 * 𝒫.V_p / 3 / sqrt(3) * F_z(𝒫, r) * c3(𝒫, λ, ϕ) * s1(𝒫, λ, ϕ) * cos(𝒫.ϕ_c) * sin(λ - 𝒫.λ_c) /
    sin(d(𝒫, λ, ϕ) / 𝒫.a) * cond(𝒫, λ, ϕ)
δw(𝒫, λ, ϕ, r) = 0.0

# CliMA prognostic variables
# compute the total energy
uˡᵒⁿ(𝒫, λ, ϕ, r) = u(𝒫, ϕ, r) + δu(𝒫, λ, ϕ, r)
uˡᵃᵗ(𝒫, λ, ϕ, r) = v(𝒫, ϕ, r) + δv(𝒫, λ, ϕ, r)
uʳᵃᵈ(𝒫, λ, ϕ, r) = w(𝒫, ϕ, r) + δw(𝒫, λ, ϕ, r)

# cv_m and R_m for moist experiment
cv_m(𝒫, ϕ, r) = 𝒫.cv_d + (𝒫.cv_v - 𝒫.cv_d) * q(𝒫, ϕ, r)
R_m(𝒫, ϕ, r) = 𝒫.R_d * (1 + (𝒫.molmass_ratio - 1) * q(𝒫, ϕ, r))

T(𝒫, ϕ, r) = Tᵥ(𝒫, ϕ, r) / (1 + 𝒫.Mᵥ * q(𝒫, ϕ, r))
e_int(𝒫, λ, ϕ, r) = cv_m(𝒫, ϕ, r) * (T(𝒫, ϕ, r) - 𝒫.T_0) + q(𝒫, ϕ, r) * 𝒫.e_int_v0
e_kin(𝒫, λ, ϕ, r) = 0.5 * (uˡᵒⁿ(𝒫, λ, ϕ, r)^2 + uˡᵃᵗ(𝒫, λ, ϕ, r)^2 + uʳᵃᵈ(𝒫, λ, ϕ, r)^2)
e_pot(𝒫, λ, ϕ, r) = 𝒫.g * r

ρ₀(𝒫, λ, ϕ, r) = p(𝒫, ϕ, r) / R_m(𝒫, ϕ, r) / T(𝒫, ϕ, r)
ρuˡᵒⁿ(𝒫, λ, ϕ, r) = ρ₀(𝒫, λ, ϕ, r) * uˡᵒⁿ(𝒫, λ, ϕ, r)
ρuˡᵃᵗ(𝒫, λ, ϕ, r) = ρ₀(𝒫, λ, ϕ, r) * uˡᵃᵗ(𝒫, λ, ϕ, r)
ρuʳᵃᵈ(𝒫, λ, ϕ, r) = ρ₀(𝒫, λ, ϕ, r) * uʳᵃᵈ(𝒫, λ, ϕ, r)
ρe(𝒫, λ, ϕ, r) = ρ₀(𝒫, λ, ϕ, r) * (e_int(𝒫, λ, ϕ, r) + e_kin(𝒫, λ, ϕ, r) + e_pot(𝒫, λ, ϕ, r))
ρq(𝒫, λ, ϕ, r) = ρ₀(𝒫, λ, ϕ, r) * q(𝒫, ϕ, r)

# Cartesian Representation (boiler plate really)
ρ₀ᶜᵃʳᵗ(𝒫, x...) = ρ₀(𝒫, lon(x...), lat(x...), rad(x...))
ρu⃗₀ᶜᵃʳᵗ(𝒫, x...) = (
    ρuʳᵃᵈ(𝒫, lon(x...), lat(x...), rad(x...)) * r̂(x...) +
    ρuˡᵃᵗ(𝒫, lon(x...), lat(x...), rad(x...)) * ϕ̂(x...) +
    ρuˡᵒⁿ(𝒫, lon(x...), lat(x...), rad(x...)) * λ̂(x...)
)
ρeᶜᵃʳᵗ(𝒫, x...) = ρe(𝒫, lon(x...), lat(x...), rad(x...))
ρqᶜᵃʳᵗ(𝒫, x...) = ρq(𝒫, lon(x...), lat(x...), rad(x...))

#####
# Held-Suarez Forcing
#####
struct HeldSuarezForcing{S} <: AbstractPhysicsComponent
    parameters::S
end

FT = Float64
day = 86400
held_suarez_parameters = (;
    k_a = FT(1 / (40 * day)),
    k_f = FT(1 / day),
    k_s = FT(1 / (4 * day)),
    ΔT_y = FT(65),
    Δθ_z = FT(10),
    T_equator = FT(294),
    T_min = FT(200),
    σ_b = FT(7 / 10),
    R_d = parameters.R_d,
    day = parameters.day,
    grav = parameters.g,
    cp_d = parameters.cp_d,
    cv_d = parameters.cv_d,
    MSLP = parameters.p0,
)

######
# Modified Held-Suarez Forcing
######
function calc_component!(source, hsf::HeldSuarezForcing, state, aux, physics)
    FT = eltype(state)

    _R_d = hsf.parameters.R_d
    _day = hsf.parameters.day
    _grav = hsf.parameters.grav
    _cp_d = hsf.parameters.cp_d
    _cv_d = hsf.parameters.cv_d
    _p0 = hsf.parameters.MSLP

    # Extract the state
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    Φ = aux.Φ

    x = aux.x
    y = aux.y
    z = aux.z
    coord = @SVector[x, y, z]

    p = calc_pressure(physics.eos, state, aux, physics.parameters)
    T = p / (ρ * _R_d)

    # Held-Suarez parameters
    k_a = hsf.parameters.k_a
    k_f = hsf.parameters.k_f
    k_s = hsf.parameters.k_s
    ΔT_y = hsf.parameters.ΔT_y
    Δθ_z = hsf.parameters.Δθ_z
    T_equator = hsf.parameters.T_equator
    T_min = hsf.parameters.T_min
    σ_b = hsf.parameters.σ_b

    # Held-Suarez forcing
    φ = @inbounds asin(coord[3] / norm(coord, 2))

    #TODO: replace _p0 with dynamic surfce pressure in Δσ calculations to account
    #for topography, but leave unchanged for calculations of σ involved in T_equil
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor

    # horizontal projection
    k = coord / norm(coord)
    P = I - k * k'

    # Apply Held-Suarez forcing
    source.ρu -= k_v * P * ρu
    source.ρe -= k_T * ρ * _cv_d * (T - T_equil)
    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module NullArrays

using Base.Broadcast: Broadcasted, DefaultArrayStyle

export NullArray

mutable struct NullArray{N, Axes<:Tuple{Vararg{Any, N}}} <: AbstractArray{Any, N}
    axes::Axes
end

Base.axes(arr::NullArray) = arr.axes
Base.size(arr::NullArray) = map(length(axes(arr)))
Base.setindex!(::NullArray, x, i...) = x
Base.setindex!(::NullArray, x, i::Integer) = x
Base.setindex!(::NullArray, x, i::CartesianIndex) = x

end