text:
# mainly use to prevent display_eval_output stealing rprint output
_output_is_locked = false

const output_buffer = PipeBuffer()
const error_buffer = PipeBuffer()

"Print the value of an Sxp using R's printing mechanism"
function rprint(io::IO, s::Ptr{S}) where S<:Sxp
    global _output_is_locked
    protect(s)
    # Rf_PrintValue can cause segfault if a S3/S4 object has custom
    # print function as it doesn't use R_tryEval
    # below mirrors Rf_PrintValue
    env = protect(newEnvironment(Const.GlobalEnv))
    defineVar(:x, s, env)
    _output_is_locked = true
    status = Ref{Cint}()
    if isObject(s) || isFunction(s)
        if isS4(s)
            tryEval(rlang_p(findNamespace("methods")[:show], :x), env, status)
        else
            tryEval(rlang_p(Const.BaseNamespace[:print], :x), env, status)
        end
    else
        # Rf_PrintValueRec not found on unix!?
        # ccall((:Rf_PrintValueRec,libR),Nothing,(Ptr{S},Ptr{EnvSxp}),s, Const.GlobalEnv)
        tryEval(rlang_p(Const.BaseNamespace[Symbol("print.default")], :x), env, status)
    end
    defineVar(:x, Const.NilValue, env)
    try
        handle_eval_stdout(io=io, force=true)
        handle_eval_stderr(as_warning=(status[] == 0))
    finally
        _output_is_locked = false
        unprotect(2)
    end
    # ggplot2's plot is displayed after `print` function is invoked,
    # so we have to clear any displayed plots.
    isdefined(Main, :IJulia) && Main.IJulia.inited && ijulia_displayplots()
    nothing
end
rprint(io::IO, r::RObject) = rprint(io::IO, r.p)
rprint(r::Ptr{S}) where S<:Sxp = rprint(stdout, r)
rprint(r::RObject) = rprint(stdout, r)

function show(io::IO,r::RObject)
    println(io, typeof(r))
    rprint(io, r)
end

abstract type RException <: Exception end
showerror(io::IO, e::RException, bt; backtrace=false) = showerror(io ,e)

# status = 2
struct RParseIncomplete <: RException
    msg::AbstractString
    RParseIncomplete(msg::AbstractString) = new(rstrip(msg))
end
showerror(io::IO, e::RParseIncomplete) = print(io, "RParseIncomplete: " * e.msg)

# status = 3
struct RParseError <: RException
    msg::AbstractString
    RParseError() = new("unknown parse error")
    RParseError(msg::AbstractString) = new(rstrip(msg))
end
showerror(io::IO, e::RParseError) = print(io, "RParseError: " * e.msg)

# status = 4
struct RParseEOF <: RException
    msg::AbstractString
    RParseEOF() = new("read end of file")
    RParseEOF(msg::AbstractString) = new(rstrip(msg))
end
showerror(io::IO, e::RParseEOF) = print(io, "RParseEOF: " * e.msg)

struct REvalError <: RException
    msg::AbstractString
    REvalError() = new("")
    REvalError(msg::AbstractString) = new(rstrip(msg))
end
showerror(io::IO, e::REvalError) = print(io, "REvalError: " * e.msg)

"""
R API callback to write console output.
"""
function write_console_ex(buf::Ptr{UInt8},buflen::Cint,otype::Cint)
    if otype == 0
        unsafe_write(output_buffer, buf, buflen)
    else
        unsafe_write(error_buffer, buf, buflen)
    end
    return nothing
end

function handle_eval_stdout(;io::IO=stdout, force::Bool=false)
    if (!_output_is_locked || force) && bytesavailable(output_buffer) != 0
        write(io, String(take!(output_buffer)))
    end
end

function handle_eval_stderr(;as_warning::Bool=false)
    if bytesavailable(error_buffer) != 0
        s = String(take!(error_buffer))
        if as_warning
            @warn "RCall.jl: " * s
        else
            throw(REvalError(s))
        end
    end
end


text:
#This example assumes that you ran through the first one.
# It will show how to load a job from a server, and set flags.

using DFControl

server_dir = "Si/"
local_dir  = "/home/ponet/Documents/Si"
job        = DFJob(server_dir, local_dir)

#set flags
setflags!(job, :ecutwfc => 25.0)  #more `Pair{Symbol, Any}`'s can be given as varargs.
#setflags! will only go through already set flags to set them,
#if you want to set new flags simply do (it can also be used to set already set flags):

setflags!(job, :ecut_rho => 80.0) #won't work, it goes through the QE documentation to find the allowed flags for the calculation files in the job, and also tries to convert the given value to what it should be.

setflags!(job, :ecutrho => 80.0, :diagonalization => "david")

#now we might not want to run all the calculations again.
#one can set the job "flow":
setflow!(job, "scf" => false) #again, more `Pair{String, Bool}`'s can be given as varargs
#This looks through all the calculations inside the job, if the `calculation.name` (derived from it's filename) contains a given `String`, it will set `calculation.run` to the corresponding `Bool`.
#If a job has calculations/calculations that don't need to run, their line will be commented out in `job.tt` script.

#now we can resubmit the job with the setd flags
submit(job)

ext: jl


text:
# A forward experiment, where the seismic data are generated
# using some models and acquisition parameters from our gallery.
# Forward modeling consists of a finite-difference simulation of the acoustic wave-equation.
# We specifically aim to save the snapshots, at given time steps in `SeisForwExpt`.

# ### Loading some packages

using GeoPhyInv
using Statistics

# ### Setting up the variables necessary to create the `Expt`

model=GIPh.Gallery.Seismic(:acou_homo1); # load a simple homogeneous acoustic model from the gallery
GIPh.Models.Seismic_addon!(model, randn_perc=0.01); # add some random noise to the model
acqgeom=GIPh.Gallery.Geom(model.mgrid,:xwell); # load a simple acquisition geometry using `mgrid` of the seismic model
tgrid = range(0.0,stop=2.0,length=2000) # generate a time grid
wav = GIPh.Utils.Wavelets.ricker(10.0, tgrid, tpeak=0.25,); # ricker wavelet
acqsrc=GIPh.Acquisition.Src_fixed(acqgeom.nss,1,[:P],wav,tgrid); # distribute the same source wavelet to all the supsersources
@info "We are ready for the modeling."

# ### Final step

# One can plot the model, source and receivers using these commands:
# `using Plots;`
# `p1=JP.seismic(model);`
# `JP.geom!(acqgeom);`
# `plot(p1);`
# Now we have all the required variables to create `SeisForwExpt` object and
# prepare the forward modelling.
# While creating, we switched the `snaps_flag` on, and instructed recording field at
# `tsnaps`.
# Once the `Expt` object is created, do the modelling "without approximately any
# memory allocations" using `mod!`

paE=SeisForwExpt(model=model,
	acqgeom=[acqgeom], acqsrc=[acqsrc],
	snaps_flag=true,
	tsnaps=[0.3, 0.4, 0.5],
	tgridmod=tgrid, verbose=true);

@time mod!(paE);

# ### Extracting snaps from `Expt`
snaps=paE[:snaps,1]; # extracting snaps of the first supersource
snaps=paE[:snaps,2]; # second supersource
@info string("The dimensions of the snaps are (nz,nx,nt)=", size(snaps))

# We can now plot snapshots using these commands:
# `p1=[heatmap(snaps[:,:,ii]) for ii in 1:3];`
# `plot(p1..., layout=(1,3), aspect_ratio=:equal)`


text:
# # R0031/3 Composite plate test

# Source code: [`R0031-3-Composite-benchmark_tut.jl`](R0031-3-Composite-benchmark_tut.jl)

# ## Description

# This is a test recommended by the National Agency for Finite Element Methods and Standards (U.K.): Test R0031/3 from NAFEMS publication R0031, “Composites Benchmarks,” February 1995. It is  a composite  (sandwich) plate of square shape, simply supported along all four edges. Uniform transverse loading is applied to the top skin. The modeled part is one quarter of the full plate here. The serendipity  quadratic hexahedra  are used, with full integration.

# The solution  can be compared with the benchmark results  in the Abaqus manual ["Abaqus Benchmarks Guide"](http://130.149.89.49:2080/v6.7/books/bmk/default.htm?startat=ch04s09anf83.html).

# We begin  by `using` the toolkit `FinEtools`.
using FinEtools

# Further, we use the linear-deformation application package.
using FinEtoolsDeforLinear

# The problem will be solved with a pre-packaged algorithm.
using FinEtoolsDeforLinear.AlgoDeforLinearModule

# For some basic statistics.
import Statistics: mean

# The material parameters are specified for an orthotropic material model.  The units are attached using the `phun` function which can take the specification of the units and spits out the numerical multiplier. Here the benchmark specifies the input parameters in the English Imperial units. The skin  material:

E1s = 1.0e7*phun("psi")
E2s = 0.4e7*phun("psi")
E3s = 0.4e7*phun("psi")
nu12s = 0.3
nu13s = 0.3
nu23s = 0.3
G12s = 0.1875e7*phun("psi")
G13s = 0.1875e7*phun("psi")
G23s = 0.1875e7*phun("psi");

# The core material:
E1c = 10.0*phun("psi")
E2c = 10.0*phun("psi")
E3c = 10e4.*phun("psi")
nu12c = 0.
nu13c = 0.
nu23c = 0.
G12c = 10.0*phun("psi")
G13c = 3.0e4*phun("psi")
G23c = 1.2e4*phun("psi");

# The magnitude  of the distributed uniform transfers loading is
tmag = 100*phun("psi");

# Now we generate the mesh.   The sandwich plate volume is divided  into a regular Cartesian grid in the $X$ and $Y$ direction in the plane of the plate, and  in the thickness direction  it is divided  into three layers, with each layer again subdivided into multiple  elements.

L = 10.0*phun("in") # side of the square plate
nL = 8 # number of elements along the side of the plate
xs = collect(linearspace(0.0, L/2, nL+1))
ys = collect(linearspace(0.0, L/2, nL+1));;

# The thicknesses are specified from the bottom of the plate: skin, core, and then again skin.
ts = [0.028; 0.75; 0.028]*phun("in")
nts = [2; 3;  2]; # number of elements through the thickness for each layer

# The `H8layeredplatex` meshing function generates the mesh and marks the elements  with a label identifying  the layer to which they belong.  We will use the label to create separate regions, with their own separate materials.
fens,fes = H8layeredplatex(xs, ys, ts, nts)

# The linear hexahedra are subsequently converted to serendipity (quadratic) elements.
fens,fes = H8toH20(fens,fes);

# The model reduction  here simply says this is a fully three-dimensional model.  The two orthotropic materials are created.
MR = DeforModelRed3D
skinmaterial = MatDeforElastOrtho(MR,
  0.0, E1s, E2s, E3s,
  nu12s, nu13s, nu23s,
  G12s, G13s, G23s,
  0.0, 0.0, 0.0)
corematerial = MatDeforElastOrtho(MR,
  0.0, E1c, E2c, E3c,
  nu12c, nu13c, nu23c,
  G12c, G13c, G23c,
  0.0, 0.0, 0.0);

# Now we are ready to create three material regions:  one for the bottom skin, one for the core, and one for the top skin. The selection of the finite elements assigned to each of the three regions is based on the label. Full Gauss quadrature  is used.
rl1 = selectelem(fens, fes, label=1)
skinbot = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes, rl1), GaussRule(3, 3)), skinmaterial))

rl3 = selectelem(fens, fes, label=3)
skintop = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes, rl3), GaussRule(3, 3)), skinmaterial))

rl2 = selectelem(fens, fes, label=2)
core = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes, rl2), GaussRule(3, 3)), corematerial));

# Note that since we did not specify the material coordinate system,  the default is assumed  (which is identical to the global Cartesian coordinate system).
@show skinbot["femm"].mcsys

# Next we select the nodes to which  essential boundary conditions  will be applied.  A node is selected  if it is within the specified box  which for the purpose of the test  is inflated in all directions by `tolerance`. The  nodes on the planes of symmetry need to be selected, and also  the nodes  along the edges (faces) to be simply supported  need to be identified.
tolerance = 0.0001*phun("in")
lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
lxL2 = selectnode(fens, box=[L/2 L/2 -Inf Inf -Inf Inf], inflate=tolerance)
ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
lyL2 = selectnode(fens, box=[-Inf Inf L/2 L/2 -Inf Inf], inflate=tolerance);

# We have four sides  of the quarter of the plate, two on each plane of symmetry, and two  along the circumference. Hence we create  four essential boundary condition definitions, one for each of the sides of the plate.
ex0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
exL2 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lxL2 )
ey0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>ly0 )
eyL2 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lyL2 );

# The traction on the top surface of the top skin is applied to the subset  of the surface mesh of the entire domain. First we compute the  boundary mesh, and then from the boundary mesh we select the surface finite elements that  "face" upward (along the positive $Z$ axis).
bfes = meshboundary(fes)
ttopl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 1.0])
Trac = FDataDict("traction_vector"=>[0.0; 0.0; -tmag],
    "femm"=>FEMMBase(IntegDomain(subset(bfes, ttopl), GaussRule(2, 3))));

# The model data  is composed of the  finite element nodes, an array  of the regions, an array of the essential boundary condition definitions, and  an array of  the traction (natural) boundary condition definitions.
modeldata = FDataDict("fens"=>fens,
 "regions"=>[skinbot, core, skintop], "essential_bcs"=>[ex0, exL2, ey0, eyL2],
 "traction_bcs"=> [Trac]
 );

# With the model data assembled,  we can now call the algorithm.
modeldata = AlgoDeforLinearModule.linearstatics(modeldata);

# The  computed solution can now be postprocessed. The displacement is reported at the center of the plate, along the line in the direction of the loading. We select all the nodes along this line.
u = modeldata["u"]
geom = modeldata["geom"]
lcenter = selectnode(fens, box=[L/2 L/2  L/2 L/2 -Inf Inf], inflate=tolerance);

# The variation of the displacement along this line  can be plotted  as (the bottom surface of the shell is at $Z=0$):
ix = sortperm(geom.values[lcenter, 3])

# Plot the data
using Gnuplot

Gnuplot.gpexec("reset session")
@gp  "set terminal windows 0 "  :-
@gp  :- geom.values[lcenter, 3][ix] u.values[lcenter, 3][ix]./phun("in") " lw 2 with lp title 'cold leg' "  :-
@gp  :- "set xlabel 'Z coordinate [in]'" :-
@gp  :- "set ylabel 'Vert displ [in]'"

# A reasonable single number to report for the deflection at the center is the average of the displacements at the nodes at the center of the plate (-0.136348):
cdis = mean(u.values[lcenter, 3])/phun("in");
println("Center node displacements $(cdis) [in]; NAFEMS-R0031-3 reference: –0.123 [in]")

# The reference displacement at the center of -0.123 [in] reported for the benchmark is evaluated from  an analytical formulation that neglects transverse  (pinching) deformation. Due to the soft core, significant pinching is observed. The solution to the benchmark  obtained in Abaqus  with incompatible hexahedral elements (with the same number of elements as in the stacked continuum shell solution) is -0.131 [in], which is close to our own solution. Hence, our own solution is probably more accurate than the reference solution because it includes an effect neglected in the benchmark solution.

# The deformed shape can be investigated  visually in `paraview` (uncomment the line at the bottom if you have `paraview` in your  PATH):
File =  "NAFEMS-R0031-3-plate.vtk"
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.VTK.H20;
    scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
# @async run(`"paraview.exe" $File`);

# Note that the  VTK file will contain element labels (which can help us distinguish between the layers) as scalar field, and the displacements as a vector field.


text:
command = `julia --trace-compile=statements_file.jl call_all_functions.jl`
run(command)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

"""
	flatten(dict::Dict{String,T}; delim = "_") where T <: Any
Given a dict of nested dicts with string keys, return a dict that is flattened (not nested) with keys between root and leaf concatenated by `delim`.

# Example
```julia-repl
julia> flatten(Dict("a" => 1, "b" => Dict("c" => 2)))
Dict{String,Any} with 2 entries:
  "b_c" => 2
  "a"   => 1
```
"""
function flatten(dict::Dict{String,T}; delim = "_") where T <: Any
	dflat = Dict{String,Any}()
	function recurse(d; s = nothing)
		if isa(d, Dict{String,T} where T <: Any)
			for (k,v) in d
				recurse(v; s = s == nothing ? k : string(s, delim, k))
			end
		else
			dflat[s] = d
		end
	end
	recurse(dict)
	return dflat
end


text:
"""
    struct CertificateMonomials <: MOI.AbstractConstraintAttribute end

A constraint attribute for the monomials indexing the
[`GramMatrixAttribute`](@ref) and [`MomentMatrixAttribute`](@ref) certificates.
"""
struct CertificateMonomials <: MOI.AbstractConstraintAttribute end

"""
    struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute end

A constraint attribute for the [`GramMatrix`](@ref) of a constraint, that is,
the positive semidefinte matrix `Q` indexed by the monomials in the vector `X`
such that ``X^\\top Q X`` is the sum-of-squares certificate of the constraint.
The
"""
struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute end

"""
    struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute end

A constraint attribute fot the `MomentMatrix` of a constraint.
"""
struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute end

"""
    struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute end

A constraint attribute fot the `LagrangianMultipliers` assiciated to the
inequalities of the domain of a constraint. There is one multiplier per
inequality and no multiplier for equalities as the equalities are handled by
reducing the polynomials over the ideal they generate instead of explicitely
creating multipliers.
"""
struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute end

function MOI.is_set_by_optimize(::Union{CertificateMonomials,
                                        GramMatrixAttribute,
                                        MomentMatrixAttribute,
                                        LagrangianMultipliers})
    return true
end


text:
# ---
# title: 1348. Tweet Counts Per Frequency
# id: problem1348
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design
# link: <https://leetcode.com/problems/tweet-counts-per-frequency/description/>
# hidden: true
# ---
#
# Implement the class `TweetCounts` that supports two methods:
#
# 1.` recordTweet(string tweetName, int time)`
#
#   * Stores the `tweetName` at the recorded `time` (in **seconds** ).
#
# 2.` getTweetCountsPerFrequency(string freq, string tweetName, int startTime,
# int endTime)`
#
#   * Returns the total number of occurrences for the given `tweetName` per **minute** , **hour** , or **day** (depending on `freq`) starting from the `startTime` (in **seconds** ) and ending at the `endTime` (in **seconds** ).
#   * `freq` is always **minute** _,_ **hour** _  or **day**_ , representing the time interval to get the total number of occurrences for the given `tweetName`.
#   * The first time interval always starts from the `startTime`, so the time intervals are `[startTime, startTime + delta*1>,  [startTime + delta*1, startTime + delta*2>, [startTime + delta*2, startTime + delta*3>, ... , [startTime + delta*i, **min** (startTime + delta*(i+1), endTime + 1)>` for some non-negative number `i` and `delta` (which depends on `freq`).
#
#
#
# **Example:**
#
#
#
#     **Input**
#     ["TweetCounts","recordTweet","recordTweet","recordTweet","getTweetCountsPerFrequency","getTweetCountsPerFrequency","recordTweet","getTweetCountsPerFrequency"]
#     [[],["tweet3",0],["tweet3",60],["tweet3",10],["minute","tweet3",0,59],["minute","tweet3",0,60],["tweet3",120],["hour","tweet3",0,210]]
#
#     **Output**
#     [null,null,null,null,[2],[2,1],null,[4]]
#
#     **Explanation**
#     TweetCounts tweetCounts = new TweetCounts();
#     tweetCounts.recordTweet("tweet3", 0);
#     tweetCounts.recordTweet("tweet3", 60);
#     tweetCounts.recordTweet("tweet3", 10);                             // All tweets correspond to "tweet3" with recorded times at 0, 10 and 60.
#     tweetCounts.getTweetCountsPerFrequency("minute", "tweet3", 0, 59); // return [2]. The frequency is per minute (60 seconds), so there is one interval of time: 1) [0, 60> - > 2 tweets.
#     tweetCounts.getTweetCountsPerFrequency("minute", "tweet3", 0, 60); // return [2, 1]. The frequency is per minute (60 seconds), so there are two intervals of time: 1) [0, 60> - > 2 tweets, and 2) [60,61> - > 1 tweet.
#     tweetCounts.recordTweet("tweet3", 120);                            // All tweets correspond to "tweet3" with recorded times at 0, 10, 60 and 120.
#     tweetCounts.getTweetCountsPerFrequency("hour", "tweet3", 0, 210);  // return [4]. The frequency is per hour (3600 seconds), so there is one interval of time: 1) [0, 211> - > 4 tweets.
#
#
#
#
# **Constraints:**
#
#   * There will be at most `10000` operations considering both `recordTweet` and `getTweetCountsPerFrequency`.
#   * `0 <= time, startTime, endTime <= 10^9`
#   * `0 <= endTime - startTime <= 10^4`
#
#
## @lc code=start
using LeetCode

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


text:
########################################
#              Preamble
########################################
using PyCall
using FFTW
using LinearAlgebra
using JLD2
using YAML
using CosmoMMF

config = YAML.load_file("../config.yaml")

########################################
#            Load in Data
########################################

@load config["input_directory"]*config["run_name"]*".jld2" den
print("\n Loading in Density Cube \n")

# den = den[1:512, 1:512, 1:512] #temporary

########################################
#  Run n' Save NEXUS/NEXUS+ Algorithm(s)
########################################

Rs = (√2) .^ (0:config["num_scales"])
print("\n loading different scales \n")

save_name = config["run_name"]

if ARGS[1] == "NEXUS"
    @time combined_NEXUS = CosmoMMF.maximum_signature(Rs, den; alg=:NEXUS);
    CosmoMMF.save_max_sigs(config["output_directory"], "Sigs_Nexus_"*save_name, combined_NEXUS)
    combined_NEXUS = nothing
elseif ARGS[1] == "NEXUSPLUS"
    @time combined_NEXUSPLUS = CosmoMMF.maximum_signature(Rs, den; alg=:NEXUSPLUS);
    CosmoMMF.save_max_sigs(config["output_directory"], "Sigs_NexusPlus_"*save_name, combined_NEXUSPLUS)
    combined_NEXUSPLUS = nothing
elseif ARGS[1] == "BOTH"
    @time combined_NEXUS = CosmoMMF.maximum_signature(Rs, den; alg=:NEXUS);
    CosmoMMF.save_max_sigs(config["output_directory"], "Sigs_Nexus_"*save_name, combined_NEXUS)
    combined_NEXUS = nothing

    @time combined_NEXUSPLUS = CosmoMMF.maximum_signature(Rs, den; alg=:NEXUSPLUS);
    CosmoMMF.save_max_sigs(config["output_directory"], "Sigs_NexusPlus_"*save_name, combined_NEXUSPLUS)
    combined_NEXUSPLUS = nothing
else
    printstyled("Please give a valid algorithm input argument: \n NEXUS, NEXUSPLUS, or BOTH"; color = :red)
end


text:
## Farfield dynamics
#==========================================================================================#

"""
    trefftz_velocity(r_i, r_j, Γ_j)

Compute the induced velocity at point ``r_i`` of the wake in the Trefftz plane due to the strength ``\\Gamma_j`` at ``r_j``.
"""
trefftz_velocity(r_i, r_j, Γ_j) = let r = r_i - r_j; Γ_j / 2π * SVector(1, 0, 0) × r / norm(r)^2 end

"""
    trefftz_influence_matrix(centers, normals, points)

Compute the aerodynamic influence coefficient matrix of the wake in the Trefftz plane given the center points, normal vectors, and the points of the wake panels.
"""
trefftz_influence_matrix(centers, normals, points) = [ dot(trefftz_velocity(r_i, r_j2, 1.) - trefftz_velocity(r_i, r_j1, 1.), n_i) for (r_i, n_i) in zip(centers, normals), (r_j2, r_j1) in zip(points[2:end], points) ]

"""
    doublet_normal_derivatives(wake_lines :: Vector{<: Line}, Δφs, normals)

Compute the normal derivative strengths of the doublets given the wake `Line`s, net doublet strengths ``\\Delta \\phis``, and associated normal vectors.
"""
function doublet_normal_derivatives(wake_lines :: Vector{<: Line}, Δφs, normals)
    centers = center.(wake_lines)
    pts 	= points(wake_lines)
    AIC   	= trefftz_influence_matrix(centers, normals, pts)
    ∂φ_∂n 	= AIC * Δφs
end

project_vector(vector, U) = vector - dot(U, vector) * U
normal(wake_proj_vector, U) = normalize(U × wake_proj_vector)
dihedral(wake_proj_vec)  = atan(wake_proj_vec[3], wake_proj_vec[2])

"""
    trefftz_preprocessing(horseshoes :: AbstractArray{<: Horseshoe}, α, β)

Project `Horseshoe`s into the Trefftz plane aligned with the wind axes angles ``\\alpha,~\\beta``, and compute normal vectors, projection angles and lengths.
"""
function trefftz_preprocessing(horseshoes :: AbstractArray{<: Horseshoe}, α, β)
    # Reference velocity for broadcasting
    U_ref = (Ref ∘ SVector)(1, 0, 0)

    # Transform to wind axes
    wake_lines     = @. body_to_wind_axes(bound_leg(horseshoes[end,:][:]), α, β)

    # Project trailing edge horseshoes' bound legs into Trefftz plane along wind axes
    wake_vectors   = @. vector(wake_lines);
    wake_proj_vecs = @. project_vector(wake_vectors, U_ref)

    # Normals, dihedral angles, and projection lengths
    wake_normals   = @. normal(wake_proj_vecs, U_ref)
    wake_dihedrals = @. dihedral(wake_proj_vecs)
    wake_lengths   = @. norm(wake_proj_vecs)

    wake_lines, wake_normals, wake_dihedrals, wake_lengths
end

"""
    trefftz_compute(Δφs, Δs, ∂φ_∂n, θs, V, ρ)

Compute the aerodynamic forces in the Trefftz plane given cumulative doublet strengths ``\\Delta \\phi``s, Trefftz panel lengths ``Δs``, doublet-normal directional derivatives ``\\partial \\phi / \\partial n``, Trefftz panel angles ``\\theta``s, the freestream speed and density ``V, \\rho``.
"""
function trefftz_compute(Δφs, Δs, ∂φ_∂n, θs, V, ρ)
    D_i = - 1/2 * ρ * sum(@. Δφs * Δs * ∂φ_∂n)
    Y   = - ρ * V * sum(@. Δφs * Δs * sin(θs))
    L   = ρ * V * sum(@. Δφs * Δs * cos(θs))

    SVector(D_i, Y, L)
end

# ifelse(symmetry, SVector(D_i, 0, 2L), SVector(D_i, Y, L))

"""
    trefftz_forces(Γs, horseshoes, freestream, ρ)

Compute the aerodynamic forces in the Trefftz plane normal to the freestream given horseshoes, their associated strengths Γs, and a density ρ.
"""
function trefftz_forces(Γs, horseshoes :: AbstractArray{<: Horseshoe}, speed, α, β, ρ)
    # Get projections of horseshoes into Trefftz plane with the associated normals, dihedral angles and lengths
    wake_lines, normals, dihedrals, Δs = trefftz_preprocessing(horseshoes, α, β)

    # Compute directional derivatives of doublets in the normal direction
    Δφs   = vec(sum(Γs, dims = 1))
    ∂φ_∂n = doublet_normal_derivatives(wake_lines, Δφs, normals)

    # Compute forces
    trefftz_compute(Δφs, Δs, ∂φ_∂n, dihedrals, speed, ρ)
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
function lift_parent_attribute(scene, attr::Symbol, default_value)
    if haskey(scene.attributes, attr)
        lift(identity, scene[attr])
    else
        lift_parent_attribute(scene.parent, attr, default_value)
    end
end

function lift_parent_attribute(::Nothing, attr::Symbol, default_value)
    default_value
end

function default_attributes(::Type{Axis}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The xlabel string."
        xlabel = ""
        "The ylabel string."
        ylabel = ""
        "The axis title string."
        title = ""
        "The font family of the title."
        titlefont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The title's font size."
        titlesize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The gap between axis and title."
        titlegap = 10f0
        "Controls if the title is visible."
        titlevisible = true
        "The horizontal alignment of the title."
        titlealign = :center
        "The color of the title"
        titlecolor = :black
        "The font family of the xlabel."
        xlabelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The font family of the ylabel."
        ylabelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The color of the xlabel."
        xlabelcolor = RGBf0(0, 0, 0)
        "The color of the ylabel."
        ylabelcolor = RGBf0(0, 0, 0)
        "The font size of the xlabel."
        xlabelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The font size of the ylabel."
        ylabelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "Controls if the xlabel is visible."
        xlabelvisible = true
        "Controls if the ylabel is visible."
        ylabelvisible = true
        "The padding between the xlabel and the ticks or axis."
        xlabelpadding = 15f0
        "The padding between the ylabel and the ticks or axis."
        ylabelpadding = 15f0 # because of boundingbox inaccuracies of ticklabels
        "The font family of the xticklabels."
        xticklabelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The font family of the yticklabels."
        yticklabelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The color of xticklabels."
        xticklabelcolor = RGBf0(0, 0, 0)
        "The color of yticklabels."
        yticklabelcolor = RGBf0(0, 0, 0)
        "The font size of the xticklabels."
        xticklabelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The font size of the yticklabels."
        yticklabelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "Controls if the xticklabels are visible."
        xticklabelsvisible = true
        "Controls if the yticklabels are visible."
        yticklabelsvisible = true
        "The space reserved for the xticklabels."
        xticklabelspace = AbstractPlotting.automatic
        "The space reserved for the yticklabels."
        yticklabelspace = AbstractPlotting.automatic
        "The space between xticks and xticklabels."
        xticklabelpad = 5f0
        "The space between yticks and yticklabels."
        yticklabelpad = 5f0
        "The counterclockwise rotation of the xticklabels in radians."
        xticklabelrotation = 0f0
        "The counterclockwise rotation of the yticklabels in radians."
        yticklabelrotation = 0f0
        "The horizontal and vertical alignment of the xticklabels."
        xticklabelalign = (:center, :top)
        "The horizontal and vertical alignment of the yticklabels."
        yticklabelalign = (:right, :center)
        "The size of the xtick marks."
        xticksize = 10f0
        "The size of the ytick marks."
        yticksize = 10f0
        "Controls if the xtick marks are visible."
        xticksvisible = true
        "Controls if the ytick marks are visible."
        yticksvisible = true
        "The alignment of the xtick marks relative to the axis spine (0 = out, 1 = in)."
        xtickalign = 0f0
        "The alignment of the ytick marks relative to the axis spine (0 = out, 1 = in)."
        ytickalign = 0f0
        "The width of the xtick marks."
        xtickwidth = 1f0
        "The width of the ytick marks."
        ytickwidth = 1f0
        "The color of the xtick marks."
        xtickcolor = RGBf0(0, 0, 0)
        "The color of the ytick marks."
        ytickcolor = RGBf0(0, 0, 0)
        "Locks interactive panning in the x direction."
        xpanlock = false
        "Locks interactive panning in the y direction."
        ypanlock = false
        "Locks interactive zooming in the x direction."
        xzoomlock = false
        "Locks interactive zooming in the y direction."
        yzoomlock = false
        "Controls if rectangle zooming affects the x dimension."
        xrectzoom = true
        "Controls if rectangle zooming affects the y dimension."
        yrectzoom = true
        "The width of the axis spines."
        spinewidth = 1f0
        "Controls if the x grid lines are visible."
        xgridvisible = true
        "Controls if the y grid lines are visible."
        ygridvisible = true
        "The width of the x grid lines."
        xgridwidth = 1f0
        "The width of the y grid lines."
        ygridwidth = 1f0
        "The color of the x grid lines."
        xgridcolor = RGBAf0(0, 0, 0, 0.1)
        "The color of the y grid lines."
        ygridcolor = RGBAf0(0, 0, 0, 0.1)
        "The linestyle of the x grid lines."
        xgridstyle = nothing
        "The linestyle of the y grid lines."
        ygridstyle = nothing
        "Controls if the bottom axis spine is visible."
        bottomspinevisible = true
        "Controls if the left axis spine is visible."
        leftspinevisible = true
        "Controls if the top axis spine is visible."
        topspinevisible = true
        "Controls if the right axis spine is visible."
        rightspinevisible = true
        "The color of the bottom axis spine."
        bottomspinecolor = :black
        "The color of the left axis spine."
        leftspinecolor = :black
        "The color of the top axis spine."
        topspinecolor = :black
        "The color of the right axis spine."
        rightspinecolor = :black
        "The forced aspect ratio of the axis. `nothing` leaves the axis unconstrained, `DataAspect()` forces the same ratio as the ratio in data limits between x and y axis, `AxisAspect(ratio)` sets a manual ratio."
        aspect = nothing
        "The vertical alignment of the axis within its suggested bounding box."
        valign = :center
        "The horizontal alignment of the axis within its suggested bounding box."
        halign = :center
        "The width of the axis."
        width = nothing
        "The height of the axis."
        height = nothing
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The relative margins added to the autolimits in x direction."
        xautolimitmargin = (0.05f0, 0.05f0)
        "The relative margins added to the autolimits in y direction."
        yautolimitmargin = (0.05f0, 0.05f0)
        "The xticks."
        xticks = AbstractPlotting.automatic
        "Format for xticks."
        xtickformat = AbstractPlotting.automatic
        "The yticks."
        yticks = AbstractPlotting.automatic
        "Format for yticks."
        ytickformat = AbstractPlotting.automatic
        "The button for panning."
        panbutton = AbstractPlotting.Mouse.right
        "The key for limiting panning to the x direction."
        xpankey = AbstractPlotting.Keyboard.x
        "The key for limiting panning to the y direction."
        ypankey = AbstractPlotting.Keyboard.y
        "The key for limiting zooming to the x direction."
        xzoomkey = AbstractPlotting.Keyboard.x
        "The key for limiting zooming to the y direction."
        yzoomkey = AbstractPlotting.Keyboard.y
        "The position of the x axis (`:bottom` or `:top`)."
        xaxisposition = :bottom
        "The position of the y axis (`:left` or `:right`)."
        yaxisposition = :left
        "Controls if the x spine is limited to the furthest tick marks or not."
        xtrimspine = false
        "Controls if the y spine is limited to the furthest tick marks or not."
        ytrimspine = false
        "The background color of the axis."
        backgroundcolor = :white
        "Controls if the ylabel's rotation is flipped."
        flip_ylabel = false
        "Constrains the data aspect ratio (`nothing` leaves the ratio unconstrained)."
        autolimitaspect = nothing
        targetlimits = BBox(0, 100, 0, 100)
        "The align mode of the axis in its parent GridLayout."
        alignmode = Inside()
        "Controls if the y axis goes upwards (false) or downwards (true)"
        yreversed = false
        "Controls if the x axis goes rightwards (false) or leftwards (true)"
        xreversed = false
    end

    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Axis has the following attributes:

$(let
    _, docs, defaults = default_attributes(Axis, nothing)
    docvarstring(docs, defaults)
end)
"""
Axis

function default_attributes(::Type{Colorbar}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The color bar label string."
        label = ""
        "The label color."
        labelcolor = RGBf0(0, 0, 0)
        "The label font family."
        labelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The label font size."
        labelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "Controls if the label is visible."
        labelvisible = true
        "The gap between the label and the ticks."
        labelpadding = 15f0
        "The font family of the tick labels."
        ticklabelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The font size of the tick labels."
        ticklabelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "Controls if the tick labels are visible."
        ticklabelsvisible = true
        "The color of the tick labels."
        ticklabelcolor = RGBf0(0, 0, 0)
        "The size of the tick marks."
        ticksize = 10f0
        "Controls if the tick marks are visible."
        ticksvisible = true
        "The ticks."
        ticks = AbstractPlotting.automatic
        "Format for ticks."
        tickformat = AbstractPlotting.automatic
        "The space reserved for the tick labels."
        ticklabelspace = AbstractPlotting.automatic
        "The gap between tick labels and tick marks."
        ticklabelpad = 5f0
        "The alignment of the tick marks relative to the axis spine (0 = out, 1 = in)."
        tickalign = 0f0
        "The line width of the tick marks."
        tickwidth = 1f0
        "The color of the tick marks."
        tickcolor = RGBf0(0, 0, 0)
        "The horizontal and vertical alignment of the tick labels."
        ticklabelalign = (:left, :center)
        "The line width of the spines."
        spinewidth = 1f0
        "Controls if the top spine is visible."
        topspinevisible = true
        "Controls if the right spine is visible."
        rightspinevisible = true
        "Controls if the left spine is visible."
        leftspinevisible = true
        "Controls if the bottom spine is visible."
        bottomspinevisible = true
        "The color of the top spine."
        topspinecolor = RGBf0(0, 0, 0)
        "The color of the left spine."
        leftspinecolor = RGBf0(0, 0, 0)
        "The color of the right spine."
        rightspinecolor = RGBf0(0, 0, 0)
        "The color of the bottom spine."
        bottomspinecolor = RGBf0(0, 0, 0)
        "The vertical alignment of the colorbar in its suggested bounding box."
        valign = :center
        "The horizontal alignment of the colorbar in its suggested bounding box."
        halign = :center
        "Controls if the colorbar is oriented vertically."
        vertical = true
        "Flips the axis to the right if vertical and to the top if horizontal."
        flipaxisposition = true
        "Flips the colorbar label if the axis is vertical."
        flip_vertical_label = false
        "The width setting of the colorbar."
        width = nothing
        "The height setting of the colorbar."
        height = nothing
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The colormap that the colorbar uses."
        colormap = :viridis
        "The range of values depicted in the colorbar."
        limits = (0f0, 1f0)
        "The align mode of the colorbar in its parent GridLayout."
        alignmode = Inside()
        "The number of steps in the heatmap underlying the colorbar gradient."
        nsteps = 100
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Colorbar has the following attributes:

$(let
    _, docs, defaults = default_attributes(Colorbar, nothing)
    docvarstring(docs, defaults)
end)
"""
Colorbar

function default_attributes(::Type{Label}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The displayed text string."
        text = "Text"
        "Controls if the text is visible."
        visible = true
        "The color of the text."
        color = RGBf0(0, 0, 0)
        "The font size of the text."
        textsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The font family of the text."
        font = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The vertical alignment of the text in its suggested boundingbox"
        valign = :center
        "The horizontal alignment of the text in its suggested boundingbox"
        halign = :center
        "The counterclockwise rotation of the text in radians."
        rotation = 0f0
        "The extra space added to the sides of the text boundingbox."
        padding = (0f0, 0f0, 0f0, 0f0)
        "The height setting of the text."
        height = Auto()
        "The width setting of the text."
        width = Auto()
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The align mode of the text in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Label has the following attributes:

$(let
    _, docs, defaults = default_attributes(Label, nothing)
    docvarstring(docs, defaults)
end)
"""
Label

function default_attributes(::Type{Box}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "Controls if the rectangle is visible."
        visible = true
        "The color of the rectangle."
        color = RGBf0(0.9, 0.9, 0.9)
        "The vertical alignment of the rectangle in its suggested boundingbox"
        valign = :center
        "The horizontal alignment of the rectangle in its suggested boundingbox"
        halign = :center
        "The extra space added to the sides of the rectangle boundingbox."
        padding = (0f0, 0f0, 0f0, 0f0)
        "The line width of the rectangle's border."
        strokewidth = 1f0
        "Controls if the border of the rectangle is visible."
        strokevisible = true
        "The color of the border."
        strokecolor = RGBf0(0, 0, 0)
        "The width setting of the rectangle."
        width = nothing
        "The height setting of the rectangle."
        height = nothing
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The align mode of the rectangle in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Box has the following attributes:

$(let
    _, docs, defaults = default_attributes(Box, nothing)
    docvarstring(docs, defaults)
end)
"""
Box

function default_attributes(::Type{Button}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The horizontal alignment of the button in its suggested boundingbox"
        halign = :center
        "The vertical alignment of the button in its suggested boundingbox"
        valign = :center
        "The extra space added to the sides of the button label's boundingbox."
        padding = (10f0, 10f0, 10f0, 10f0)
        "The font size of the button label."
        textsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The text of the button label."
        label = "Button"
        "The font family of the button label."
        font = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The width setting of the button."
        width = Auto()
        "The height setting of the button."
        height = Auto()
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The radius of the rounded corners of the button."
        cornerradius = 4
        "The number of poly segments used for each rounded corner."
        cornersegments = 10
        "The line width of the button border."
        strokewidth = 2f0
        "The color of the button border."
        strokecolor = :transparent
        "The color of the button."
        buttoncolor = RGBf0(0.94, 0.94, 0.94)
        "The color of the label."
        labelcolor = :black
        "The color of the label when the mouse hovers over the button."
        labelcolor_hover = :black
        "The color of the label when the mouse clicks the button."
        labelcolor_active = :white
        "The color of the button when the mouse clicks the button."
        buttoncolor_active = COLOR_ACCENT[]
        "The color of the button when the mouse hovers over the button."
        buttoncolor_hover = COLOR_ACCENT_DIMMED[]
        "The number of clicks that have been registered by the button."
        clicks = 0
        "The align mode of the button in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Button has the following attributes:

$(let
    _, docs, defaults = default_attributes(Button, nothing)
    docvarstring(docs, defaults)
end)
"""
Button

function default_attributes(::Type{LineAxis})
    Attributes(
        endpoints = (Point2f0(0, 0), Point2f0(100, 0)),
        trimspine = false,
        limits = (0f0, 100f0),
        flipped = false,
        flip_vertical_label = false,
        ticksize = 10f0,
        tickwidth = 1f0,
        tickcolor = RGBf0(0, 0, 0),
        tickalign = 0f0,
        ticks = AbstractPlotting.automatic,
        tickformat = AbstractPlotting.automatic,
        ticklabelalign = (:center, :top),
        ticksvisible = true,
        ticklabelrotation = 0f0,
        ticklabelsize = 20f0,
        ticklabelcolor = RGBf0(0, 0, 0),
        ticklabelsvisible = true,
        spinewidth = 1f0,
        label = "label",
        labelsize = 20f0,
        labelcolor = RGBf0(0, 0, 0),
        labelvisible = true,
        ticklabelspace = AbstractPlotting.automatic,
        ticklabelpad = 5f0,
        labelpadding = 15f0,
        reversed = false,
    )
end

function default_attributes(::Type{Slider}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The horizontal alignment of the slider in its suggested bounding box."
        halign = :center
        "The vertical alignment of the slider in its suggested bounding box."
        valign = :center
        "The width setting of the slider."
        width = Auto()
        "The height setting of the slider."
        height = Auto()
        "The range of values that the slider can pick from."
        range = 0:10
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The start value of the slider or the value that is closest in the slider range."
        startvalue = 0
        "The current value of the slider."
        value = 0
        "The width of the slider line"
        linewidth = 15
        "The color of the slider when the mouse hovers over it."
        color_active_dimmed = COLOR_ACCENT_DIMMED[]
        "The color of the slider when the mouse clicks and drags the slider."
        color_active = COLOR_ACCENT[]
        "The color of the slider when it is not interacted with."
        color_inactive = RGBf0(0.94, 0.94, 0.94)
        "Controls if the slider has a horizontal orientation or not."
        horizontal = true
        "The align mode of the slider in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Slider has the following attributes:

$(let
    _, docs, defaults = default_attributes(Slider, nothing)
    docvarstring(docs, defaults)
end)
"""
Slider

function default_attributes(::Type{Toggle}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The horizontal alignment of the toggle in its suggested bounding box."
        halign = :center
        "The vertical alignment of the toggle in its suggested bounding box."
        valign = :center
        "The width of the toggle."
        width = 60
        "The height of the toggle."
        height = 28
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The number of poly segments in each rounded corner."
        cornersegments = 15
        # strokewidth = 2f0
        # strokecolor = :transparent
        "The color of the border when the toggle is inactive."
        framecolor_inactive = RGBf0(0.94, 0.94, 0.94)
        "The color of the border when the toggle is hovered."
        framecolor_active = COLOR_ACCENT_DIMMED[]
        # buttoncolor = RGBf0(0.2, 0.2, 0.2)
        "The color of the toggle button."
        buttoncolor = COLOR_ACCENT[]
        "Indicates if the toggle is active or not."
        active = false
        "The duration of the toggle animation."
        toggleduration = 0.15
        "The border width as a fraction of the toggle height "
        rimfraction = 0.33
        "The align mode of the toggle in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Toggle has the following attributes:

$(let
    _, docs, defaults = default_attributes(Toggle, nothing)
    docvarstring(docs, defaults)
end)
"""
Toggle

function default_attributes(::Type{Legend}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The horizontal alignment of the legend in its suggested bounding box."
        halign = :center
        "The vertical alignment of the legend in its suggested bounding box."
        valign = :center
        "The width setting of the legend."
        width = Auto()
        "The height setting of the legend."
        height = Auto()
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = false
        "The font family of the legend group titles."
        titlefont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The font size of the legend group titles."
        titlesize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The horizontal alignment of the legend group titles."
        titlehalign = :center
        "The vertical alignment of the legend group titles."
        titlevalign = :center
        "Controls if the legend titles are visible."
        titlevisible = true
        "The color of the legend titles"
        titlecolor = :black
        "The group title positions relative to their groups. Can be `:top` or `:left`."
        titleposition = :top
        "The font size of the entry labels."
        labelsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "The font family of the entry labels."
        labelfont = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "The color of the entry labels."
        labelcolor = :black
        "The horizontal alignment of the entry labels."
        labelhalign = :left
        "The vertical alignment of the entry labels."
        labelvalign = :center
        "The additional space between the legend content and the border."
        padding = (10f0, 10f0, 10f0, 10f0)
        "The additional space between the legend and its suggested boundingbox."
        margin = (0f0, 0f0, 0f0, 0f0)
        "The background color of the legend."
        bgcolor = :white
        "The color of the legend border."
        framecolor = :black
        "The line width of the legend border."
        framewidth = 1f0
        "Controls if the legend border is visible."
        framevisible = true
        "The size of the rectangles containing the legend markers."
        patchsize = (20f0, 20f0)
        "The color of the border of the patches containing the legend markers."
        patchstrokecolor = :transparent
        "The line width of the border of the patches containing the legend markers."
        patchstrokewidth = 1f0
        "The color of the patches containing the legend markers."
        patchcolor = :transparent
        "The default entry label."
        label = "undefined"
        "The number of banks in which the legend entries are grouped. Columns if the legend is vertically oriented, otherwise rows."
        nbanks = 1
        "The gap between the label of one legend entry and the patch of the next."
        colgap = 20
        "The gap between the entry rows."
        rowgap = 4
        "The gap between the patch and the label of each legend entry."
        patchlabelgap = 5
        "The default points used for LineElements in normalized coordinates relative to each label patch."
        linepoints = [Point2f0(0, 0.5), Point2f0(1, 0.5)]
        "The default line width used for LineElements."
        linewidth = 3
        "The default marker points used for MarkerElements in normalized coordinates relative to each label patch."
        markerpoints = [Point2f0(0.5, 0.5)]
        "The default marker size used for MarkerElements."
        markersize = 12
        "The default marker stroke width used for MarkerElements."
        markerstrokewidth = 1
        "The default poly points used for PolyElements in normalized coordinates relative to each label patch."
        polypoints = [Point2f0(0, 0), Point2f0(1, 0), Point2f0(1, 1), Point2f0(0, 1)]
        "The default poly stroke width used for PolyElements."
        polystrokewidth = 1
        "The orientation of the legend (:horizontal or :vertical)."
        orientation = :vertical
        "The gap between each group title and its group."
        titlegap = 15
        "The gap between each group and the next."
        groupgap = 30
        "The horizontal alignment of entry groups in their parent GridLayout."
        gridshalign = :center
        "The vertical alignment of entry groups in their parent GridLayout."
        gridsvalign = :center
        "The align mode of the legend in its parent GridLayout."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
Legend has the following attributes:

$(let
    _, docs, defaults = default_attributes(Legend, nothing)
    docvarstring(docs, defaults)
end)
"""
Legend

function attributenames(::Type{LegendEntry})
    (:label, :labelsize, :labelfont, :labelcolor, :labelhalign, :labelvalign,
        :patchsize, :patchstrokecolor, :patchstrokewidth, :patchcolor,
        :linepoints, :markerpoints, :markersize, :markerstrokewidth, :linewidth,
        :polypoints, :polystrokewidth)
end

function extractattributes(attributes::Attributes, typ::Type)
    extracted = Attributes()
    for name in attributenames(typ)
        if haskey(attributes, name)
            extracted[name] = attributes[name]
        end
    end
    extracted
end

function default_attributes(::Type{LScene}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The height setting of the scene."
        height = nothing
        "The width setting of the scene."
        width = nothing
        "Controls if the parent layout can adjust to this element's width"
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height"
        tellheight = true
        "The horizontal alignment of the scene in its suggested bounding box."
        halign = :center
        "The vertical alignment of the scene in its suggested bounding box."
        valign = :center
        "The alignment of the scene in its suggested bounding box."
        alignmode = Inside()
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
LScene has the following attributes:

$(let
    _, docs, defaults = default_attributes(LScene, nothing)
    docvarstring(docs, defaults)
end)
"""
LScene

function default_attributes(::Type{Textbox}, scene)
    attrs, docdict, defaultdict = @documented_attributes begin
        "The height setting of the textbox."
        height = Auto()
        "The width setting of the textbox."
        width = Auto()
        "Controls if the parent layout can adjust to this element's width."
        tellwidth = true
        "Controls if the parent layout can adjust to this element's height."
        tellheight = true
        "The horizontal alignment of the textbox in its suggested bounding box."
        halign = :center
        "The vertical alignment of the textbox in its suggested bounding box."
        valign = :center
        "The alignment of the textbox in its suggested bounding box."
        alignmode = Inside()
        "A placeholder text that is displayed when the saved string is nothing."
        placeholder = "Click to edit..."
        "The currently stored string."
        stored_string = nothing
        "The currently displayed string (for internal use)."
        displayed_string = nothing
        "Controls if the displayed text is reset to the stored text when defocusing the textbox without submitting."
        reset_on_defocus = false
        "Controls if the textbox is defocused when a string is submitted."
        defocus_on_submit = true
        "Text size."
        textsize = lift_parent_attribute(scene, :fontsize, 20f0)
        "Text color."
        textcolor = :black
        "Text color for the placeholder."
        textcolor_placeholder = RGBf0(0.5, 0.5, 0.5)
        "Font family."
        font = lift_parent_attribute(scene, :font, "DejaVu Sans")
        "Color of the box."
        boxcolor = :transparent
        "Color of the box when focused."
        boxcolor_focused = :transparent
        "Color of the box when focused."
        boxcolor_focused_invalid = RGBAf0(1, 0, 0, 0.3)
        "Color of the box when hovered."
        boxcolor_hover = :transparent
        "Color of the box border."
        bordercolor = RGBf0(0.80, 0.80, 0.80)
        "Color of the box border when hovered."
        bordercolor_hover = COLOR_ACCENT_DIMMED[]
        "Color of the box border when focused."
        bordercolor_focused = COLOR_ACCENT[]
        "Color of the box border when focused and invalid."
        bordercolor_focused_invalid = RGBf0(1, 0, 0)
        "Width of the box border."
        borderwidth = 2f0
        "Padding of the text against the box."
        textpadding = (10, 10, 10, 10)
        "If the textbox is focused and receives text input."
        focused = false
        "Corner radius of text box."
        cornerradius = 8
        "Corner segments of one rounded corner."
        cornersegments = 20
        "Validator that is called with validate_textbox(string, validator) to determine if the current string is valid. Can by default be a RegEx that needs to match the complete string, or a function taking a string as input and returning a Bool. If the validator is a type T (for example Float64), validation will be `tryparse(string, T)`."
        validator = str -> true
        "Restricts the allowed unicode input via is_allowed(char, restriction)."
        restriction = nothing
        "The color of the cursor."
        cursorcolor = :transparent
    end
    (attributes = attrs, documentation = docdict, defaults = defaultdict)
end

@doc """
    Textbox(parent::Scene; bbox = nothing, kwargs...)
Textbox has the following attributes:
$(let
    _, docs, defaults = default_attributes(Textbox, nothing)
    docvarstring(docs, defaults)
end)
"""
Textbox


text:
function MakeModel(A,b,C,XLB,XUB,ObjInd;Cuts = [])
    m,n = size(A);
    ObjNum = size(C)[1];
    MyModel = Model();
    @variable(MyModel, XLB[i] <= X[i=1:n] <= XUB[i])
    @objective(MyModel, Min, sum(C[ObjInd,i]*X[i] for i=1:n))
    for i = 1:m
        @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    # for c = 1:length(Cuts)
    #     @constraint(MyModel,sum(Cuts[c].Coeff[i]*(sum(C[i,j]*X[j] for j = 1:n )+d[i])  for i = 1:ObjNum) >= Cuts[c].RHS)
    # end
    return MyModel;
end


text:
function test_sparse_vanilla_same(m)
    vanilla_solver = ValueIterationSolver(max_iterations=500, belres=1e-8)
    sparse_solver = SparseValueIterationSolver(max_iterations=500, belres=1e-8)
    vanilla_policy = solve(vanilla_solver, m)
    sparse_policy = solve(sparse_solver, m)
    @test vanilla_policy.policy == sparse_policy.policy
    @test isapprox(vanilla_policy.util, sparse_policy.util; atol=1e-3)
    for s in states(m)
        si = stateindex(m, s)
        for a in actions(m, s)
            ai = actionindex(m, a)
            @test isapprox(vanilla_policy.qmat[si, ai], sparse_policy.qmat[si, ai]; atol=1e-3)
        end
    end
end

m1 = SimpleGridWorld(size=(10,10))
test_sparse_vanilla_same(m1)

m2 = SpecialGridWorld()
test_sparse_vanilla_same(m2)

## Two state MDP to test for terminal states

struct TwoStatesMDP <: MDP{Int, Int} end

POMDPs.states(mdp::TwoStatesMDP) = 1:2
POMDPs.stateindex(mdp::TwoStatesMDP, s) = s
POMDPs.actions(mdp::TwoStatesMDP) = 1:2
POMDPs.actionindex(mdp::TwoStatesMDP, a) = a
POMDPs.discount(mdp::TwoStatesMDP) = 0.95
POMDPs.transition(mdp::TwoStatesMDP, s, a) = SparseCat([a], [1.0])
POMDPs.reward(mdp::TwoStatesMDP, s, a, sp) = float(sp == 2)
POMDPs.isterminal(mdp::TwoStatesMDP, s) = s == 2

mdp = TwoStatesMDP()
solver = ValueIterationSolver(verbose = true)
policy = solve(solver, mdp)

sparsesolver = SparseValueIterationSolver(verbose=true)
sparsepolicy = solve(sparsesolver, mdp)

@test sparsepolicy.qmat == policy.qmat
@test value(sparsepolicy, 2) ≈ 0.0
@test_throws String solve(SparseValueIterationSolver(), TigerPOMDP())

@inferred solve(SparseValueIterationSolver(verbose=false), mdp)


text:
# Absolute Percentage Error
ape(yhat, y) = abs((y-yhat) ./ y)

# Median Absolute Percentage Error
medape(yhat, y) = 100.0 * median(ape(yhat, y))

# Mean Absolute Percentage Error
mape(yhat, y) = 100.0 * mean(ape(yhat, y))

# Mean Absolute Scaled Error, see https://en.wikipedia.org/wiki/Mean_absolute_scaled_error
# Note that it is not defined if the y vector is constant => all forecast errors are 0.0.
function mase{N <: Number}(yhat::Vector{N}, y::Vector{N})
    T = length(y)
    @assert length(yhat) == T
    sum_abs_errors = sumabs(yhat .- y)
    sum_abs_forecast_errors = 0.0
    for t in 2:T
        sum_abs_forecast_errors += abs(y[t] - y[t-1])
    end
    return sum_abs_errors / (T/(T-1) * sum_abs_forecast_errors)
end

# R2 is R squared which is how much more variation is explained by the
# fitted values (yhat) compared to just using the mean as the model?
r2(yhat, y) = 1.0 - ( sum((y .- yhat).^2) / sum((y .- mean(y)).^2) )

# R2adj is R squared adjusted for the number of estimated variables K.
function r2adjusted(yhat, y, k)
  n = length(y)
  1.0 - ( (sum((y .- yhat).^2)/(n-k)) / (sum((y .- mean(y)).^2)/(n-1)) )
end

# Describe the goodness of fit in terms a user can easily digest.
# If k, the number of estimated variables of the model, is given we use
# R2adjusted instead of R2.
function describe_goodness_of_fit(yhat, y, k = false)
  if k != false
    @sprintf("MAPE = %.2f%%, R2adj = %.2f, MASE = %.2f", mape(yhat, y), r2adjusted(yhat, y, k), mase(yhat, y))
  else
    @sprintf("MAPE = %.2f%%, R2 = %.2f, MASE = %.2f", mape(yhat, y), r2(yhat, y), mase(yhat, y))
  end
end


text:
#=
    Comparison tests of the kernel output (for each data size)
        against the baseline tests
=#

using Test

inputfiles = readdir("../test/test_output", join=true)

outputfiles = readdir("test/test_output", join=true)

for (input, output) in zip(inputfiles, outputfiles)
    datasize = input[21:32]
    @testset "$datasize Comparison " begin

        @test success(`cmp --quiet $input $output`)
    end
end

return true


text:
function loadsky(obj::ASCIIString,mypsf,nu::Array)

    skyst = init_SKY()
    sky0 = lecture(obj)
    if size(sky0)[3] == 1
        sky0 = squeeze(sky0/maximum(sky0),3)
        skyst.sky,skyst.alpha = sky2cube(sky0,nu)
        skyconv = cubefilter(skyst.sky,mypsf)
        skyst.sig = sqrt(mean(skyconv.^2)/100)
        skyst.var = skyst.sig*skyst.sig
        skyst.noise = skyst.sig*randn( size(skyst.sky)[1], size(skyst.sky)[1], size(mypsf)[3])
        skyst.mydata = skyconv + skyst.noise
        for z in 1:length(nu)
            push!(skyst.sumsky2, sum(skyst.sky[:,:,z].^2))
        end
    elseif size(sky0)[3] != 1
        skyst = init_SKY()
        skyst.sky = lecture(obj)
        skyst.mydata = cubefilter(skyst.sky,mypsf)
    end
    return skyst
end

function loadsky_dirty(obj::ASCIIString,mypsf,nu::Array)
    skyst = init_SKY_dirty()
    skyst.sky = lecture(obj)
    skyst.mydata = lecture(obj)
    for z in 1:length(nu)
        push!(skyst.sumsky2, sum(skyst.sky[:,:,z].^2))
    end
    return skyst
end

function loadpsf(psf::ASCIIString,M::Int)
    psfst = init_PSF()
    psfcube = lecture(psf)
    d = round((size(psfcube)[1])/2)
    psfcube = psfcube[d-128:d+127,d-128:d+127,:]
    psfavg = cubeaverage(psfcube,M)
    psfst.nu, psfst.nu0 = cubefreq(psf,psfcube,M)
    psfst.mypsf = cropcubexy(psfavg,size(psfcube)[1])
    psfst.mypsfadj = flipdim(flipdim(psfst.mypsf,1),2)
    return psfst
end

function loadpsf_dirty(psf::ASCIIString)
    psfst = init_PSF_dirty()
    psfcube = lecture(psf)
    psfst.nu, psfst.nu0 = cubefreqchiara(size(psfcube)[3])
    psfst.mypsf = cropcubexy(psfcube,size(psfcube)[1])
    psfst.mypsfadj = flipdim(flipdim(psfst.mypsf,1),2)
    return psfst
end

function loadparam(nspat,nfreq,nspec,nxy,niter,lastiter,nitermax)
    algost = init_Algoparam()
    algost.nspat = nspat
    algost.nfreq = nfreq
    algost.nspec = nspec
    algost.nxy = nxy
    algost.niter = niter
    algost.lastiter = lastiter
    algost.nitermax = nitermax
    return algost
end

function loadarray(rhop,rhot,rhov,rhos,μt,μv,mueps,nspat,nfreq,nxy,mydata,mypsfadj,parallel)
    admmst = init_Admmarray(parallel)
    admmst.s = zeros(Float64,nxy,nxy,nfreq)
    admmst.taus = zeros(Float64,nxy,nxy,nfreq)
    admmst.sh = zeros(Float64,nxy,nxy,nfreq)
    admmst.tauv = zeros(Float64,nxy,nxy,nfreq)
    admmst.v = zeros(Float64,nxy,nxy,nfreq)

    if parallel == "true"
        println("Parallel muffin")
        admmst.taup = SharedArray(Float64,nxy,nxy,nfreq)
        admmst.p = SharedArray(Float64,nxy,nxy,nfreq)
        admmst.t = SharedArray(Float64,nxy,nxy,nfreq,nspat)
        admmst.taut = SharedArray(Float64,nxy,nxy,nfreq,nspat)
        admmst.wlt = SharedArray(Float64,nxy,nxy,nfreq)
    else
        admmst.taup = zeros(Float64,nxy,nxy,nfreq)
        admmst.p = zeros(Float64,nxy,nxy,nfreq)
        admmst.t = zeros(Float64,nxy,nxy,nfreq,nspat)
        admmst.taut = zeros(Float64,nxy,nxy,nfreq,nspat)
        admmst.wlt = zeros(Float64,nxy,nxy,nfreq)
    end

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

    admmst.x = copy(mydata)
    admmst.xmm = zeros(Float64,nxy,nxy,nfreq)
    admmst.spectralwlt = zeros(Float64,nxy,nxy,nfreq)
    admmst.fty = cubefilter(mydata,mypsfadj)
    admmst.rhop = rhop
    admmst.rhot = rhot
    admmst.rhov = rhov
    admmst.rhos = rhos
    admmst.μt = μt
    admmst.μv = μv
    admmst.mueps = mueps
    admmst.tt = rhot*nspat
    admmst.mu = mueps + rhop + admmst.tt + rhos
    return admmst
end

function loadtools(nitermax,nfreq,nxy)
    toolst = init_TOOLS()
    toolst.snr = Float64[]
    toolst.tol1 = Float64[]
    toolst.tol2 = Float64[]
    toolst.tol3 = Float64[]
    toolst.tol4 = Float64[]
    toolst.tol5 = Float64[]
    toolst.err = zeros(Float64,nitermax,nfreq)
    toolst.errorrec = zeros(Float64,nxy,nxy,nfreq)
    toolst.errorest = zeros(Float64,nfreq)
    toolst.errorraw = zeros(Float64,nfreq)
    return toolst
end

function loadarray_ws(rhop,rhot,rhov,rhos,μt,μv,mueps,nspat,nfreq,nxy,mydata,mypsfadj,file)
    println("warm start")
    admmst = load(file,"admmst")
    return admmst
end


text:
# TODO: optimize

# Q(Θ) = Θ'AΘ + b'Θ + c
# A is stochastic average of .25 * X'X
# b is stochastic average of ∇f(Θ) - .25 * X'X * Θ
type LogRegMM{W <: Weight} <: OnlineStat{XYInput}
	β::VecF
    H::Matrix{Float64}  # Storage
    A::Matrix{Float64}
    b::VecF
    η::Float64
    weight::W
end
function LogRegMM(p::Integer, wt::Weight = LearningRate(); η::Float64 = 1.0)
    @assert 0.0 < η <= 1.0
    LogRegMM(zeros(p), zeros(p, p), zeros(p, p), zeros(p), η, wt)
end
function LogRegMM(x::AMat, y::AVec, wt::Weight = LearningRate())
    o = LogRegMM(size(x, 2), wt)
    fit!(o, x, y)
    o
end
function Base.show(io::IO, o::LogRegMM)
    printheader(io, name(o))
    print_item(io, "β", coef(o))
    print_item(io, "nobs", nobs(o))
end

function _fit!(o::LogRegMM, x::AVec, y::Real, γ::Float64)
    γ *= o.η
    BLAS.syrk!('U', 'N', γ, x / 2., 1.0 - γ, o.H)
    H = Symmetric(o.H)
    smooth!(o.A, H, γ)
    smooth!(o.b, (y - predict(o, x)) * x + H * o.β, γ)
    # try
    #     o.β = o.A \ o.b
    # end
    o.β = o.β + γ * inv(o.H) * (y - predict(o, x)) * x
end

value(o::LogRegMM) = coef(o)
coef(o::LogRegMM) = o.β
linpred(o::LogRegMM, x::AVec) = dot(x, o.β)
linpred(o::LogRegMM, x::AMat) = x * o.β
predict(o::LogRegMM, x) = predict(LogisticRegression(), linpred(o, x))
loss(o::LogRegMM, x, y) = loss(LogisticRegression(), y, linpred(o, x))


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# VectorAffineDecisionFunction #
# ========================== #
mutable struct VectorAffineDecisionConstraintBridge{T, S} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.VectorAffineFunction{T}, S}
    decision_function::VectorAffineDecisionFunction{T}
    set::S
end

function MOIB.Constraint.bridge_constraint(::Type{VectorAffineDecisionConstraintBridge{T, S}},
                                           model,
                                           f::VectorAffineDecisionFunction{T},
                                           set::S) where {T, S}
    # All decisions have been mapped to the variable part terms
    # at this point.
    F = MOI.VectorAffineFunction{T}
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    f.variable_part,
                                    set)
    # Save the constraint index, the decision function, and the set, to allow modifications
    return VectorAffineDecisionConstraintBridge{T, S}(constraint, f, set)
end

function MOI.supports_constraint(::Type{<:VectorAffineDecisionConstraintBridge{T}},
                                 ::Type{<:VectorAffineDecisionFunction{T}},
                                 ::Type{<:MOI.AbstractVectorSet}) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:VectorAffineDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:VectorAffineDecisionConstraintBridge{T, S}}) where {T, S}
    return [(MOI.VectorAffineFunction{T}, S)]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:VectorAffineDecisionConstraintBridge{T}},
                                              ::Type{<:VectorAffineDecisionFunction{T}},
                                              S::Type{<:MOI.AbstractVectorSet}) where T
    return VectorAffineDecisionConstraintBridge{T, S}
end

MOI.get(b::VectorAffineDecisionConstraintBridge{T, S}, ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, S}) where {T, S} = 1
MOI.get(b::VectorAffineDecisionConstraintBridge{T, S}, ::MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, S}) where {T, S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return bridge.decision_function
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return bridge.set
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return bridge.constraint
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::VectorAffineDecisionConstraintBridge{T,S},
                 f::VectorAffineDecisionFunction{T}) where {T,S}
    # Update bridge functions and function constant
    bridge.decision_function = f
    # Modify constraint function
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            f.variable_part)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorAffineDecisionConstraintBridge{T,S}, change::S) where {T,S}
    f = bridge.decision_function
    bridge.set = change
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::MOI.VectorConstantChange) where {T,S}
    f = bridge.decision_function
    # Modify variable part of decision function
    f.variable_part.constants .= change.new_constant
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            f.variable_part)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::MOI.MultirowChange) where {T,S}
    f = bridge.decision_function
    # Modify variable part of decision function
    modify_coefficients!(f.variable_part.terms, change.decision, change.new_coefficients)
    # Modify the variable part of the constraint as usual
    MOI.modify(model, bridge.constraint, change)
    return nothing
end


text:
module Awaits

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) Awaits

export @taskgroup, @cancelscope, @in, @go, @await, @check, cancel!

using MacroTools: postwalk, @capture

struct Cancelled <: Exception end

function locking(f, l)
    lock(l)
    try
        return f()
    finally
        unlock(l)
    end
end

"""
    SyncSeq([values])

A minimalistic wrapper of `Vector{Any}` that can be used from multiple
threads.
"""
struct SyncSeq
    values::Vector{Any}
    lock::Threads.SpinLock
end

SyncSeq() = SyncSeq([])
SyncSeq(values::Vector{Any}) = SyncSeq(values, Threads.SpinLock())

Base.push!(v::SyncSeq, x) = locking(v.lock) do
    push!(v.values, x)
end

Base.iterate(v::SyncSeq, i=1) = locking(v.lock) do
    iterate(v.values, i)
end

# For `collect` in test code:
Base.IteratorSize(::Type{<:SyncSeq}) = Base.SizeUnknown()

struct TaskContext
    tasks::SyncSeq
    listening::Vector{Threads.Atomic{Bool}}
    cancelables::Vector{Threads.Atomic{Bool}}
end

function TaskContext()
    c = Threads.Atomic{Bool}(false)
    return TaskContext(SyncSeq(), [c], [c])
end

function newcontext!(ctx::TaskContext)
    c = Threads.Atomic{Bool}(false)
    push!(ctx.cancelables, c)

    return TaskContext(ctx.tasks, push!(copy(ctx.listening), c), [c])
end

function shouldstop(ctx::TaskContext)
    for c in ctx.listening
        c[] && return true
    end
    return false
end

"""
    cancel!(context::TaskContext)

Cancel (stop) execution of task `context`.
"""
function cancel!(ctx::TaskContext)
    for c in ctx.cancelables
        c[] = true
    end
end

_copy(ctx::TaskContext) =
    TaskContext(ctx.tasks, copy(ctx.listening), copy(ctx.cancelables))

_tasks(ctx::TaskContext) =
    (begin
         shouldstop(ctx) && cancel!(ctx)
         task
     end
     for task in ctx.tasks)

const ctx_varname = gensym("taskcontext")

"""
    @taskgroup(body)

Create a new task context.  This is an extension of the `@sync` block.
"""
macro taskgroup(body)
    ctx_var = esc(ctx_varname)
    sync_var = esc(Base.sync_varname)
    quote
        # TODO: "inherit" old task context when a new context is created
        let $ctx_var = TaskContext(),
            $sync_var = $ctx_var.tasks

            v = $(esc(body))
            Base.sync_end(_tasks($ctx_var))
            v
        end
    end
end

"""
    @in context body

Enter into the task `context` so that other macros like `@go` works as
expected.  This should be used inside a function in which the task
`context` is passed by a caller (for example, using [`@go f(_)`](@ref
@go) macro).
"""
macro in(ctx, body)
    ctx_var = esc(ctx_varname)
    sync_var = esc(Base.sync_varname)
    quote
        let $ctx_var = $(esc(ctx)),
            $sync_var = $ctx_var.tasks

            $(esc(body))
        end
    end
end

"""
    @cancelscope(body) :: TaskContext

Create a task context with a new scope of cancellation.  Tasks
launched inside a `@cancelscope` can be cancelled independently from
the tasks outside to this scope.  This scope returns a context object
that can be used to cancel the task by calling [`cancel!`](@ref) on
this object.
"""
macro cancelscope(body)
    ctx_var = esc(ctx_varname)
    quote
        let $ctx_var = newcontext!($ctx_var)
            $(esc(body))
            $ctx_var
        end
    end
end

substitute_context(ctx, body) =
    postwalk(body) do ex
        if @capture(ex, f_(xs__)) && :_ in xs
            :($f($((x == :_ ? ctx : x for x in xs)...)))
        else
            ex
        end
    end

"""
    @await body

Handle an exception returned from `body`.  If `body` evaluates to an
`Exception`, it cancels all the tasks attached to the current context.

The variable `_` in code like `@await f(_, x, y)` is replaced by the
current task context.  It must be invoked inside `@taskgroup` or `@in`
macros.
"""
macro await(body)
    quote
        ans = $(esc(substitute_context(ctx_varname, body)))
        if ans isa Exception
            cancel!($(esc(ctx_varname)))
            return ans
        end
        ans
    end
end

"""
    @go body

A wrapper of `Threads.@spawn` that makes cancellation work.

The variable `_` in code like `@go f(_, x, y)` is replaced by the
current task context.  It must be invoked inside `@taskgroup` or `@in`
macros.
"""
macro go(body)
    quote
        # Copy task context to avoid touching `Vector` without a lock.
        # It means that the `@spawn`ed task will miss future additions
        # of cancellation tokens in this task.  However, since current
        # task context and the copied task context share at least one
        # cancellation token (the root one), the cancellation state
        # will propagate eventually.
        let $ctx_varname = $_copy($ctx_varname)
            $Threads.@spawn $(@__MODULE__).@await begin
                $(substitute_context(ctx_varname, body))
            end
        end
    end |> esc
end

"""
    @check [context]

Check `context` and exit the function (do `return Cancelled()`) if it
is required.
"""
macro check()
    quote
        @check($ctx_varname)
    end |> esc
end

macro check(ctx)
    ctx_var = esc(ctx)
    quote
        shouldstop($ctx_var) && return Cancelled()
    end
end

end # module


text:
using StatsBase
using Flux
using Images

FILTER_SPACE = [16, 32, 64, 128, 256, 512]
SKIP_CONNECTION_TYPES = ["add", "concat"]

abstract type MGene end
abstract type MBlockGene <: MGene end

###################################
########### Convolution ###########
###################################
mutable struct MConv <: MGene
    in_dims::Tuple
    out_dims::Tuple
    kernel_size::Tuple
    padding::Tuple
end

function MConv(in_dims::Tuple)
    out_dims = in_dims
    out_channels = sample(FILTER_SPACE)
    out_dims[3] = out_channels
    kernel_size = 3, 3
    padding = 1, 1
    MConv(in_dims, out_dims, kernel_size, padding)
end

function MConv(in_dims::Tuple, out_dims::Tuple)
    kernel_size = 3, 3
    padding = 1, 1
    MConv(in_dims, out_dims, kernel_size, padding)
end

function get_phenotype(prev_layer::MGene, cur_layer::MConv)
    return Conv(cur_layer.kernel_size, prev_layer.out_dims[3]=>cur_layer.out_dims[3], pad=cur_layer.padding)
end

###################################
# Depthwise Separable Convolution #
###################################
mutable struct MDepthwiseConv <: MGene
    in_dims::Tuple
    out_dims::Tuple
    kernel_size::Tuple
    padding::Tuple
end

function MDepthwiseConv(in_dims::Tuple)
    out_dims = in_dims
    out_channels = sample(FILTER_SPACE)
    out_dims[3] = out_channels
    kernel_size = 3, 3
    padding = 1, 1
    MDepthwiseConv(in_dims, out_dims, kernel_size, padding)
end

function MDepthwiseConv(in_dims::Tuple, out_dims::Tuple)
    kernel_size = 3, 3
    padding = 1, 1
    MDepthwiseConv(in_dims, out_dims, kernel_size, padding)
end

function get_phenotype(prev_layer::MGene, cur_layer::MDepthwiseConv)
    return DepthwiseConv(cur_layer.kernel_size, prev_layer.out_dims[3]=>cur_layer.out_dims[3], pad=cur_layer.padding)
end

###################################
####### Batch Normalization #######
###################################
mutable struct MBatchNorm <: MGene
    in_dims::Tuple
    out_dims::Tuple
end

function MBatchNorm(in_dims::Tuple)
    MBatchNorm(in_dims, in_dims)
end

function get_phenotype(prev_layer::MGene, cur_layer::MBatchNorm)
    return BatchNorm(cur_layer.out_dims[3])
end

###################################
######### ReLU Activation #########
###################################

mutable struct MReLU <: MGene
    in_dims::Tuple
    out_dims::Tuple
end

function MReLU(in_dims::Tuple)
    MReLU(in_dims, in_dims)
end

function get_phenotype(prev_layer::MGene, cur_layer::MReLU)
    return x -> relu.(x)
end

###################################
########### Max Pooling ###########
###################################
mutable struct MMaxPool <: MGene
    in_dims::Tuple
    out_dims::Tuple
    kernel_size::Tuple
end

function MMaxPool(in_dims::Tuple)
    out_dims = in_dims
    out_dims[1] = in_dims[1] / 2
    out_dims[2] = in_dims[2] / 2
    kernel_size = 2, 2
    MMaxPool(in_dims, out_dims, kernel_size)
end

function get_phenotype(prev_layer::MGene, cur_layer::MMaxPool)
    return MaxPool(cur_layer.kernel_size)
end

###################################
######### Skip Connection #########
###################################
mutable struct MSkipConnection <: MGene
    in_dims::Tuple
    out_dims::Tuple
    type::String
end

function MSkipConnection(in_dims::Tuple, out_dims::Tuple)
    # type = sample(SKIP_CONNECTION_TYPES)
    type = "concat"
    MSkipConnection(in_dims, out_dims, type)
end

function get_phenotype(prev_layer::MGene, cur_layer::MSkipConnection, last_pheno)

    if cur_layer.type == "add"
        return Nothing
    elseif cur_layer.type == "concat"
        # input_dims = 40, 40,
        # a = randn(Float32, 2, 2, 4, 1)
        # b = padarray(a, Fill(0, (1, 1, 0, 0)))
        return SkipConnection(last_pheno, (mx, x) -> cat(mx, x, dims=3))
    end
end

max_stars_repo_name: mahmoud-al-najar/neuroevolution

max_issues_repo_name: mahmoud-al-najar/neuroevolution

max_forks_repo_name: mahmoud-al-najar/neuroevolution


text:
using AtomicLevels

function build_args(kwargs)
    map(kwargs) do kv
        "$(kv[1])=$(kv[2])"
    end
end

macro bsr_exe(exe)
    exe_str = string(exe)
    @eval begin
        function ($exe)(args...; kwargs...)
            exe = string(bsr, "/", $exe_str)
            cmd = `$exe $args $(build_args(kwargs))`
            info("Executing $cmd")
            run(cmd)
            info("Finished executing $cmd")
        end
        export $exe
    end
end

macro bsr_mpi_exe(exe)
    exe_str = string(exe)
    @eval begin
        function ($exe)(np, args...; kwargs...)
            exe = string(bsr, "/", $exe_str)
            cmd = `mpiexec -np $np $exe $args $(build_args(kwargs))`
            info("Executing $cmd")
            run(cmd)
            info("Finished executing $cmd")
        end
        export $exe
    end
end

typealias StateDict Dict{Term,Vector{AbstractString}}
state_dict() = StateDict()


text:
using LeafOptics
using Test

@testset verbose = true "LeafOptics Test" begin
    @testset "Spectra" begin
        for FT in [Float32, Float64]
            wls = WaveLengthSet{FT}();
            bio = LeafBiophysics{FT}(wls);
            lha = HyperspectralAbsorption{FT}(wls);

            leaf_spectra!(bio, wls, lha);
            @test true;
            leaf_spectra!(bio, wls, lha; APAR_car=false);
            @test true;
            leaf_spectra!(bio, wls, lha; APAR_car=false, α=FT(59));
            @test true;

            leaf_spectra!(bio, wls, FT(0.1), FT(0.45), FT(0.05), FT(0.25));
            @test true;
        end;
    end;

    @testset "PAR & APAR" begin
        for FT in [Float32, Float64]
            wls = WaveLengthSet{FT}();
            bio = LeafBiophysics{FT}(wls);
            rad = HyperspectralRadiation{FT}();

            par,apar = leaf_PAR(bio, wls, rad);
            @test true;
            par,apar = leaf_PAR(bio, wls, rad; APAR_car=false);
            @test true;
        end;
    end;

    @testset "SIF" begin
        for FT in [Float32, Float64]
            wls = WaveLengthSet{FT}();
            bio = LeafBiophysics{FT}(wls);
            rad = HyperspectralRadiation{FT}();

            sif_b,sif_f = leaf_SIF(bio, wls, rad, FT(0.01));
            @test true;
            sif_b,sif_f = leaf_SIF(bio, wls, rad, FT(0.01); ϕ_photon=false);
            @test true;
        end;
    end;
end;


text:
module SevenDayStreamflowForecasts

import Base.show
using Dates: DateTime, DateFormat
import HTTP
using CSV: File
import JSON
using DataFrames: DataFrame, select!, insertcols!

const SDF_URL = "http://www.bom.gov.au/water/7daystreamflow/"
const SITES_URL = SDF_URL * "content/config/site_config.json"

const SITE_PROPERTIES = ["station_name", "bom_id", "awrc_id",
                         "product_status", "product_id_public",
                         "area", "lat", "lon"]

const HEADER_DELIM = "#"

export SDF, get_forecasts, close!

"""
    sdf = SDF()

Open the connection and download the information of SDF service sites e.g. name, ID, AWRC ID and description.

# Fields
* `sites`: Site information table

# Examples
```julia
julia> sdf = SDF();
julia> sdf.sites
208×8 DataFrame. Omitted printing of 5 columns
│ Row │ station_name                                 │ bom_id │ awrc_id  │
│     │ String                                       │ String │ String   │
├─────┼──────────────────────────────────────────────┼────────┼──────────┤
│ 1   │ Timbarra River D/S Wilkinson Creek           │ 584008 │ 223212   │
│ 2   │ St. Pauls River above Avoca                  │ 592003 │ 18311    │
│ 3   │ Macquarie River D/S Elizabeth River Junction │ 093026 │ 18312    │
│ 4   │ Macquarie River at Fosterville               │ 093059 │ 18313    │
...
```
"""
struct SDF
	sites::DataFrame

	function SDF()
		sites = get_sites()
		new(sites)
	end
end

"""
    sites = get_sites()

Download site information e.g. AWRC ID and names.
"""
function get_sites()::DataFrame
	r = HTTP.get(SITES_URL)

	features = JSON.parse(String(r.body))["stations"]["features"]

	sites = DataFrame()
	for key in SITE_PROPERTIES
		sites[!, key] = [x["properties"][key] for x in features]
	end

    return sites
end

"""
    data, header = get_forecasts(sdf::SDF, site_id::AbstractString)

Return seasonal forecasts of a site.

# Arguments
* `sdf` : SDF object
* `site_id`: AWRC ID of the site. The ID can found in the table from `get_sites()`

# Examples
```julia
julia> using Dates
julia> data, header = get_forecasts(sdf,"410730");
julia> data
264×14 DataFrame. Omitted printing of 11 columns
│ Row │ DateTime            │ Time                  │ Observed Rainfall (mm/hour) │
│     │ DateTime            │ String                │ Union{Missing, Float64}     │
├─────┼─────────────────────┼───────────────────────┼─────────────────────────────┤
│ 1   │ 2020-09-16T10:00:00 │ 2020-09-16 10:00 AEST │ 0.0                         │
│ 2   │ 2020-09-16T11:00:00 │ 2020-09-16 11:00 AEST │ 0.0                         │
│ 3   │ 2020-09-16T12:00:00 │ 2020-09-16 12:00 AEST │ 0.0                         │
...
julia> println(header)
Australian Bureau of Meteorology
Short term streamflow forecasts
hourly_barplot_ens
Cotter River at Gingera (410730)
Catchment area:", 130.0,"km^2
...
```
"""
function get_forecasts(sdf::SDF,
                       site_id::AbstractString)::Tuple{DataFrame,String}
    row = sdf.sites[sdf.sites.awrc_id .== site_id, :]
    isempty(row) && return (DataFrame(), "")

    data_url = get_url(site_id, row.product_id_public[1])

	r = HTTP.get(data_url)

	body_buf = IOBuffer(String(r.body))

	header = extract_header!(body_buf, HEADER_DELIM)
	new_header = prune_header(header, HEADER_DELIM)

	body_buf = seek(body_buf, 0)
	data = DataFrame(File(body_buf, comment=HEADER_DELIM, missingstring="-"))

    df = DateFormat("yyyy-mm-ddTHH:MM")
    function to_datetime(dt)
        return DateTime(join(split(dt)[1:2], "T"), df)
    end
    insertcols!(data, 1, "DateTime"=>to_datetime.(data.Time))

	return data, new_header
end

function show(io::IO, sdf::SDF)
    isempty(sdf.sites) && return
    println("The SDF Sites:")
    show(io, first(sdf.sites, 6))
    println(io)
    println(io, "...")
    return
end

"""
    close!(sdf::SDF)

Close the connection and reset the SDF object.
"""
function close!(sdf::SDF)
    select!(sdf.sites, Int[])
    return
end

"""
    url = get_url(site_id::AbstractString, fc_date::Date,
                  drainage::AbstractString, basin::AbstractString)

Return the URL to download data of a site.
"""
function get_url(site_id::AbstractString, product_id::AbstractString)::AbstractString
    url = "http://www.bom.gov.au/fwo/$(product_id)/data/$(site_id)/$(site_id)_hourly_barplot_ens.csv"
    return url
end

"""
    header = extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}}, delim::AbstractString)

Extract the header document of the data file. Note that it moves the position of body_buf.
"""
function extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}},
                         delim::AbstractString)::Array{String,1}
    header = String[]
    for line in eachline(body_buf)
        startswith(line, delim) ? push!(header, line) : break
    end
    return header
end

"""
    new_header = prune_header(header::AbstractString, delim::AbstractString)

Prune the extracted header document. It drops blank lines and double quotations from the raw document.
"""
function prune_header(header::Array{String,1}, delim::AbstractString)::String
    new_header = String[]
    for line in header
        #TODO: Give more information about the format error.
        startswith(line, delim) || throw(IOError("A header line does not start with $(delim)."))

        the_line = rstrip(line[3:end])
        isempty(the_line) && continue
        istart = 1
        i = findfirst('\"', the_line)
        (i != nothing) && (istart = i+1)
        iend = length(the_line)
        (the_line[end] == '\"') && (iend -= 1)

        push!(new_header, the_line[istart:iend])
    end
    return join(new_header, "\n")
end

end # module


text:
using EasyEvents
using Test

include("utils.jl")

@testset "EasyEvents.jl" begin
    include("FileMTimeEvent_tests.jl")
    include("FileSizeEvent_tests.jl")
    include("FileContentEvent_tests.jl")
    include("PropertyChangedEvent_test.jl")
end


text:
struct EigenStates{T}
    states::Matrix{T}
    energies::Vector{Float64}
end

struct KEigenStates{D, T}
    kmesh::KMesh{D}
    states::Vector{Matrix{T}}
    energies::Vector{Vector{Float64}}
end

struct KEigenAmplitudes{D, T, N}
    kmesh::KMesh{D}
    amplitudes::Vector{Array{N, T}}
    energies::Vector{Vector{Float64}}
end

struct KBandStates{D, T}
    kpath::KPath{D}
    states::Vector{Matrix{T}}
    energies::Vector{Vector{Float64}}
end

struct  KBandAmplitudes{D, T, N}
    kpath::KPath{D}
    amplitudes::Vector{Array{N, T}}
    energies::Vector{Vector{Float64}}
end


text:
@define_scenario SimpleScenario = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end

ξ₁ = SimpleScenario(24.0, 28.0, 500.0, 100.0, probability = 0.4)
ξ₂ = SimpleScenario(28.0, 32.0, 300.0, 300.0, probability = 0.6)

simple_sp = StochasticProgram([ξ₁, ξ₂], Deterministic(), GLPK.Optimizer)

@first_stage simple_sp = begin
    @decision(model, x₁ >= 40)
    @decision(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end

@second_stage simple_sp = begin
    @known x₁ x₂
    @uncertain q₁ q₂ d₁ d₂ from SimpleScenario
    @recourse(model, 0 <= y₁ <= d₁)
    @recourse(model, 0 <= y₂ <= d₂)
    @objective(model, Max, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end