text:
function add_vars!(
   mn_data,
   ref_subnet,
   subnet_idx,
   model,
   va, vm, pg, qg, p, q, p_dc, q_dc, g, b, b_fr, b_to, g_fr, g_to, l_series,
   c_inv_series, f, b_shunt, cost_pg, cost_dcline, area_pg, zone_pg, redispatch_pg,
   redispatch_qg, redispatch_vm, redispatch_va, p_i, q_i,
   constraints
   )
   # Add Optimization and State Variables
   # ------------------------------------

   # println("Adding variables for subnet $subnet_idx.")
   # Check if DC subnet
   # Automatically DC if fixed frequency and the base frequency is 0
   dc_subnet = ((!ref_subnet[:variable_f]) && (ref_subnet[:f_base] == 0))
   if (:f_max in keys(ref_subnet))
      # println("ref_subnet[:f_max]: $(ref_subnet[:f_max])")
      # Also DC if max frequency is 0
      dc_subnet = dc_subnet || (ref_subnet[:f_max] == 0)
   end
   if dc_subnet
      println("Subnetwork $subnet_idx is DC. Using power flow equations for a DC network.")
   end
   # Add voltage angles va for each bus
   # print("subnet $(subnet_idx). Adding variables\n")
   if !dc_subnet
      va[subnet_idx] = @variable(
         model,
         [i in keys(ref_subnet[:bus])],
         base_name="va_$(subnet_idx)",
         start = mn_data["sn"]["$(subnet_idx)"]["bus"]["$i"]["va"]
         # start = 0.0
      )
   end
   # println([mn_data["sn"]["$(subnet_idx)"]["bus"]["$i"]["va"] for i in keys(ref_subnet[:bus])])
   # note: [i in keys(ref[:bus])] adds one `va` variable for each bus in the subnetwork
   # Add voltage magnitudes vm for each bus
   vm[subnet_idx] = @variable(
      model,
      [i in keys(ref_subnet[:bus])],
      base_name = "vm_$(subnet_idx)",
      upper_bound = ref_subnet[:bus][i]["vmax"],
      lower_bound = ref_subnet[:bus][i]["vmin"],
      start = mn_data["sn"]["$(subnet_idx)"]["bus"]["$i"]["vm"]
      # start = 1.0
   )
   for bus_i in keys(ref_subnet[:bus])
      if !(subnet_idx in keys(constraints[:vm_llim]))
         constraints[:vm_llim][subnet_idx] = Dict{Int64,Any}()
      end
      if !(subnet_idx in keys(constraints[:vm_ulim]))
         constraints[:vm_ulim][subnet_idx] = Dict{Int64,Any}()
      end
      constraints[:vm_llim][subnet_idx][bus_i] = LowerBoundRef(vm[subnet_idx][bus_i])
      constraints[:vm_ulim][subnet_idx][bus_i] = UpperBoundRef(vm[subnet_idx][bus_i])
   end
   # note: this variable also includes the voltage magnitude limits and a starting value

   # Add active power generation variable pg for each generator (including limits)
   # NOTE: If any generator limit is at the PSS/E default of +/- 9999.0, it is applied as +/- Inf here
   pg[subnet_idx] = @variable(
      model,
      [i in keys(ref_subnet[:gen])],
      base_name="pg_$(subnet_idx)",
      lower_bound=ref_subnet[:gen][i]["pmin"]*ref_subnet[:baseMVA] != -9999.0 ? ref_subnet[:gen][i]["pmin"] : -Inf,
      upper_bound=ref_subnet[:gen][i]["pmax"]*ref_subnet[:baseMVA] != 9999.0 ? ref_subnet[:gen][i]["pmax"] : Inf,
      start = mn_data["sn"]["$(subnet_idx)"]["gen"]["$i"]["pg"]
   )
   # Add reactive power generation variable qg for each generator (including limits)

   if !dc_subnet
      qg[subnet_idx] =  @variable(
         model,
         [i in keys(ref_subnet[:gen])],
         base_name="qg_$(subnet_idx)",
         lower_bound=ref_subnet[:gen][i]["qmin"],
         upper_bound=ref_subnet[:gen][i]["qmax"],
         start = mn_data["sn"]["$(subnet_idx)"]["gen"]["$i"]["qg"]
      )
   end

   for gen_i in keys(ref_subnet[:gen])
      if !(subnet_idx in keys(constraints[:pg_llim]))
         constraints[:pg_llim][subnet_idx] = Dict{Int64,Any}()
      end
      if !(subnet_idx in keys(constraints[:pg_ulim]))
         constraints[:pg_ulim][subnet_idx] = Dict{Int64,Any}()
      end
      constraints[:pg_llim][subnet_idx][gen_i] = LowerBoundRef(pg[subnet_idx][gen_i])
      constraints[:pg_ulim][subnet_idx][gen_i] = UpperBoundRef(pg[subnet_idx][gen_i])

      if !dc_subnet
         if !(subnet_idx in keys(constraints[:qg_llim]))
            constraints[:qg_llim][subnet_idx] = Dict{Int64,Any}()
         end
         if !(subnet_idx in keys(constraints[:qg_ulim]))
            constraints[:qg_ulim][subnet_idx] = Dict{Int64,Any}()
         end
         constraints[:qg_llim][subnet_idx][gen_i] = LowerBoundRef(qg[subnet_idx][gen_i])
         constraints[:qg_ulim][subnet_idx][gen_i] = UpperBoundRef(qg[subnet_idx][gen_i])
      end
   end

   # if dc_subnet
   #    delete_lower_bound.(qg[subnet_idx])
   #    delete_upper_bound.(qg[subnet_idx])
   #    for i in keys(ref_subnet[:gen])
   #       @constraint(model, qg[subnet_idx][i] == 0)
   #    end
   # end

   # Add power flow variables p to represent the active power flow for each branch
   # Apply thermal limits here to bound the variables. Give very generous thermal limits (5x the angle limit at the base frequency) if no thermal limits are defined
   # rate_a = ["rate_a" in keys(ref_subnet[:branch][l]) ? ref_subnet[:branch][l]["rate_a"] : (5*(1.1)^2)/ref_subnet[:branch][l]["br_x"]*sin(max(ref_subnet[:branch][l]["angmax"],ref_subnet[:branch][l]["angmin"])) for (l, i, j) in ref_subnet[:arcs]]
   p[subnet_idx] = @variable(
      model,
      [(l, i, j) in ref_subnet[:arcs]],
      base_name = "p_$(subnet_idx)",
      # upper_bound = abs("rate_a" in keys(ref_subnet[:branch][l]) ? ref_subnet[:branch][l]["rate_a"] : (5*(1.1)^2)/ref_subnet[:branch][l]["br_x"]*sin(max(ref_subnet[:branch][l]["angmax"],ref_subnet[:branch][l]["angmin"]))),
      # lower_bound = -abs("rate_a" in keys(ref_subnet[:branch][l]) ? ref_subnet[:branch][l]["rate_a"] : (5*(1.1)^2)/ref_subnet[:branch][l]["br_x"]*sin(max(ref_subnet[:branch][l]["angmax"],ref_subnet[:branch][l]["angmin"])))
   )
   # Add power flow variables q to represent the reactive power flow for each branch
   if !dc_subnet
      q[subnet_idx] = @variable(
         model,
         [(l, i, j) in ref_subnet[:arcs]],
         base_name = "q_$(subnet_idx)",
         # upper_bound = abs("rate_a" in keys(ref_subnet[:branch][l]) ? ref_subnet[:branch][l]["rate_a"] : (5*(1.1)^2)/ref_subnet[:branch][l]["br_x"]*sin(max(ref_subnet[:branch][l]["angmax"],ref_subnet[:branch][l]["angmin"]))),
         # lower_bound = -abs("rate_a" in keys(ref_subnet[:branch][l]) ? ref_subnet[:branch][l]["rate_a"] : (5*(1.1)^2)/ref_subnet[:branch][l]["br_x"]*sin(max(ref_subnet[:branch][l]["angmax"],ref_subnet[:branch][l]["angmin"])))
      )
   end
   # note: ref_subnet[:arcs] includes both the from (i,j) and the to (j,i) sides of a branch

   # Add power flow variables p_dc to represent the active power flow for each HVDC line
   p_dc[subnet_idx] = @variable(
      model,
      [a in ref_subnet[:arcs_dc]],
      base_name = "p_dc_$(subnet_idx)",
   )
   # Add power flow variables q_dc to represent the reactive power flow at each HVDC terminal
   if !dc_subnet
      q_dc[subnet_idx] = @variable(
         model,
         [a in ref_subnet[:arcs_dc]],
         base_name = "q_dc_$(subnet_idx)",
      )
   end

   for (l,dcline) in ref_subnet[:dcline]
      f_idx = (l, dcline["f_bus"], dcline["t_bus"])
      t_idx = (l, dcline["t_bus"], dcline["f_bus"])

      JuMP.set_lower_bound(p_dc[subnet_idx][f_idx], dcline["pminf"])
      JuMP.set_upper_bound(p_dc[subnet_idx][f_idx], dcline["pmaxf"])
      JuMP.set_lower_bound(q_dc[subnet_idx][f_idx], dcline["qminf"])
      JuMP.set_upper_bound(q_dc[subnet_idx][f_idx], dcline["qmaxf"])

      JuMP.set_lower_bound(p_dc[subnet_idx][t_idx], dcline["pmint"])
      JuMP.set_upper_bound(p_dc[subnet_idx][t_idx], dcline["pmaxt"])
      JuMP.set_lower_bound(q_dc[subnet_idx][f_idx], dcline["qmint"])
      JuMP.set_upper_bound(q_dc[subnet_idx][f_idx], dcline["qmaxt"])
   end

   if ref_subnet[:variable_f]
      ref_subnet[:f_min] = max(ref_subnet[:f_min],0) # don't allow negative frequencies
      ref_subnet[:f_min] = min(ref_subnet[:f_min],ref_subnet[:f_max]) # don't allow minimum larger than maximum
      if ref_subnet[:f_min] == ref_subnet[:f_max]
         ref_subnet[:variable_f] = false
         ref_subnet[:f_fixed] = ref_subnet[:f_min]
      else
         f_init = (ref_subnet[:f_min]<=ref_subnet[:f_base]<=ref_subnet[:f_max]) ? ref_subnet[:f_base] : (ref_subnet[:f_min] + ref_subnet[:f_max])/2
         f[subnet_idx] = @variable(
            model,
            base_name="f_$(subnet_idx)",
            lower_bound=ref_subnet[:f_min],
            upper_bound=ref_subnet[:f_max],
            start=f_init
         )
      end
   end
end


text:
"""
Main module for `Oceananigans.jl` -- a Julia software for fast, friendly, flexible,
data-driven, ocean-flavored fluid dynamics on CPUs and GPUs.
"""
module Oceananigans

if VERSION < v"1.6"
    error("This version of Oceananigans.jl requires Julia v1.6 or newer.")
end

export
    # Architectures
    CPU, GPU,

    # Logging
    OceananigansLogger,

    # Grids
    Center, Face,
    Periodic, Bounded, Flat,
    FullyConnected, LeftConnected, RightConnected,
    RectilinearGrid,
    LatitudeLongitudeGrid,
    ConformalCubedSphereFaceGrid,
    xnodes, ynodes, znodes, nodes,

    # Immersed boundaries
    ImmersedBoundaryGrid, GridFittedBoundary, GridFittedBottom, ImmersedBoundaryCondition,

    # Advection schemes
    CenteredSecondOrder, CenteredFourthOrder, UpwindBiasedFirstOrder, UpwindBiasedThirdOrder, UpwindBiasedFifthOrder, WENO5,
    VectorInvariant, EnergyConservingScheme, EnstrophyConservingScheme,

    # Boundary conditions
    BoundaryCondition,
    FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition,
    FieldBoundaryConditions,

    # Fields and field manipulation
    Field, CenterField, XFaceField, YFaceField, ZFaceField,
    Average, Integral, Reduction, BackgroundField,
    interior, set!, compute!, regrid!, location,

    # Forcing functions
    Forcing, Relaxation, LinearTarget, GaussianMask, AdvectiveForcing,

    # Coriolis forces
    FPlane, ConstantCartesianCoriolis, BetaPlane, NonTraditionalBetaPlane,

    # BuoyancyModels and equations of state
    Buoyancy, BuoyancyTracer, SeawaterBuoyancy,
    LinearEquationOfState, TEOS10,
    BuoyancyField,

    # Surface wave Stokes drift via Craik-Leibovich equations
    UniformStokesDrift,

    # Turbulence closures
    VerticalScalarDiffusivity,
    HorizontalScalarDiffusivity,
    ScalarDiffusivity,
    VerticalScalarBiharmonicDiffusivity,
    HorizontalScalarBiharmonicDiffusivity,
    ScalarBiharmonicDiffusivity,
    SmagorinskyLilly,
    AnisotropicMinimumDissipation,
    ConvectiveAdjustmentVerticalDiffusivity,
    RiBasedVerticalDiffusivity,
    IsopycnalSkewSymmetricDiffusivity,
    FluxTapering,
    VerticallyImplicitTimeDiscretization,

    # Lagrangian particle tracking
    LagrangianParticles,

    # Models
    NonhydrostaticModel,
    HydrostaticFreeSurfaceModel,
    ShallowWaterModel, ConservativeFormulation, VectorInvariantFormulation,
    PressureField,
    fields,

    # Hydrostatic free surface model stuff
    VectorInvariant, ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface,
    HydrostaticSphericalCoriolis,
    PrescribedVelocityFields,

    # Time stepping
    Clock, TimeStepWizard, time_step!,

    # Simulations
    Simulation, run!, Callback, iteration, stopwatch,
    iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded,
    erroring_NaNChecker!,

    # Diagnostics
    StateChecker, CFL, AdvectiveCFL, DiffusiveCFL,

    # Output writers
    NetCDFOutputWriter, JLD2OutputWriter, Checkpointer,
    TimeInterval, IterationInterval, AveragedTimeInterval, SpecifiedTimes,
    AndSchedule, OrSchedule,

    # Output readers
    FieldTimeSeries, FieldDataset, InMemory, OnDisk,

    # Abstract operations
    ∂x, ∂y, ∂z, @at, KernelFunctionOperation,

    # MultiRegion and Cubed sphere
    MultiRegionGrid, XPartition,
    ConformalCubedSphereGrid,

    # Utils
    prettytime, apply_regionally!, construct_regionally, @apply_regionally, MultiRegionObject

using Printf
using Logging
using Statistics
using LinearAlgebra
using CUDA
using Adapt
using DocStringExtensions
using OffsetArrays
using FFTW
using JLD2
using NCDatasets

using Base: @propagate_inbounds
using Statistics: mean

import Base:
    +, -, *, /,
    size, length, eltype,
    iterate, similar, show,
    getindex, lastindex, setindex!,
    push!

#####
##### Abstract types
#####

"""
    AbstractModel

Abstract supertype for models.
"""
abstract type AbstractModel{TS} end

"""
    AbstractDiagnostic

Abstract supertype for diagnostics that compute information from the current
model state.
"""
abstract type AbstractDiagnostic end

"""
    AbstractOutputWriter

Abstract supertype for output writers that write data to disk.
"""
abstract type AbstractOutputWriter end

#####
##### Place-holder functions
#####

function run_diagnostic! end
function write_output! end
function location end
function instantiated_location end
function tupleit end
function fields end
function prognostic_fields end
function tracer_tendency_kernel_function end

#####
##### Include all the submodules
#####

# Basics
include("Architectures.jl")
include("Units.jl")
include("Grids/Grids.jl")
include("Utils/Utils.jl")
include("Logger.jl")
include("Operators/Operators.jl")
include("BoundaryConditions/BoundaryConditions.jl")
include("Fields/Fields.jl")
include("AbstractOperations/AbstractOperations.jl")
include("Advection/Advection.jl")
include("Solvers/Solvers.jl")
include("Distributed/Distributed.jl")

# Physics, time-stepping, and models
include("Coriolis/Coriolis.jl")
include("BuoyancyModels/BuoyancyModels.jl")
include("StokesDrift.jl")
include("TurbulenceClosures/TurbulenceClosures.jl")
include("LagrangianParticleTracking/LagrangianParticleTracking.jl")
include("Forcings/Forcings.jl")

include("ImmersedBoundaries/ImmersedBoundaries.jl")
include("TimeSteppers/TimeSteppers.jl")
include("Models/Models.jl")

# Output and Physics, time-stepping, and models
include("Diagnostics/Diagnostics.jl")
include("OutputWriters/OutputWriters.jl")
include("OutputReaders/OutputReaders.jl")
include("Simulations/Simulations.jl")

# Abstractions for distributed and multi-region models
include("MultiRegion/MultiRegion.jl")
include("CubedSpheres/CubedSpheres.jl")

#####
##### Needed so we can export names from sub-modules at the top-level
#####

using .Logger
using .Architectures
using .Utils
using .Advection
using .Grids
using .BoundaryConditions
using .Fields
using .Coriolis
using .BuoyancyModels
using .StokesDrift
using .TurbulenceClosures
using .LagrangianParticleTracking
using .Solvers
using .Forcings
using .ImmersedBoundaries
using .Distributed
using .Models
using .TimeSteppers
using .Diagnostics
using .OutputWriters
using .OutputReaders
using .Simulations
using .AbstractOperations
using .MultiRegion
using .CubedSpheres

function __init__()
    threads = Threads.nthreads()
    if threads > 1
        @info "Oceananigans will use $threads threads"

        # See: https://github.com/CliMA/Oceananigans.jl/issues/1113
        FFTW.set_num_threads(4*threads)
    end

    if CUDA.has_cuda()
        @debug "CUDA-enabled GPU(s) detected:"
        for (gpu, dev) in enumerate(CUDA.devices())
            @debug "$dev: $(CUDA.name(dev))"
        end

        CUDA.allowscalar(false)
    end
end

end # module

ext: jl


text:
    word_list = (
        "alias", "consequatur", "aut", "perferendis", "sit", "voluptatem",
        "accusantium", "doloremque", "aperiam", "eaque", "ipsa", "quae", "ab",
        "illo", "inventore", "veritatis", "et", "quasi", "architecto",
        "beatae", "vitae", "dicta", "sunt", "explicabo", "aspernatur", "aut",
        "odit", "aut", "fugit", "sed", "quia", "consequuntur", "magni",
        "dolores", "eos", "qui", "ratione", "voluptatem", "sequi", "nesciunt",
        "neque", "dolorem", "ipsum", "quia", "dolor", "sit", "amet",
        "consectetur", "adipisci", "velit", "sed", "quia", "non", "numquam",
        "eius", "modi", "tempora", "incidunt", "ut", "labore", "et", "dolore",
        "magnam", "aliquam", "quaerat", "voluptatem", "ut", "enim", "ad",
        "minima", "veniam", "quis", "nostrum", "exercitationem", "ullam",
        "corporis", "nemo", "enim", "ipsam", "voluptatem", "quia", "voluptas",
        "sit", "suscipit", "laboriosam", "nisi", "ut", "aliquid", "ex", "ea",
        "commodi", "consequatur", "quis", "autem", "vel", "eum", "iure",
        "reprehenderit", "qui", "in", "ea", "voluptate", "velit", "esse",
        "quam", "nihil", "molestiae", "et", "iusto", "odio", "dignissimos",
        "ducimus", "qui", "blanditiis", "praesentium", "laudantium", "totam",
        "rem", "voluptatum", "deleniti", "atque", "corrupti", "quos",
        "dolores", "et", "quas", "molestias", "excepturi", "sint",
        "occaecati", "cupiditate", "non", "provident", "sed", "ut",
        "perspiciatis", "unde", "omnis", "iste", "natus", "error",
        "similique", "sunt", "in", "culpa", "qui", "officia", "deserunt",
        "mollitia", "animi", "id", "est", "laborum", "et", "dolorum", "fuga",
        "et", "harum", "quidem", "rerum", "facilis", "est", "et", "expedita",
        "distinctio", "nam", "libero", "tempore", "cum", "soluta", "nobis",
        "est", "eligendi", "optio", "cumque", "nihil", "impedit", "quo",
        "porro", "quisquam", "est", "qui", "minus", "id", "quod", "maxime",
        "placeat", "facere", "possimus", "omnis", "voluptas", "assumenda",
        "est", "omnis", "dolor", "repellendus", "temporibus", "autem",
        "quibusdam", "et", "aut", "consequatur", "vel", "illum", "qui",
        "dolorem", "eum", "fugiat", "quo", "voluptas", "nulla", "pariatur",
        "at", "vero", "eos", "et", "accusamus", "officiis", "debitis", "aut",
        "rerum", "necessitatibus", "saepe", "eveniet", "ut", "et",
        "voluptates", "repudiandae", "sint", "et", "molestiae", "non",
        "recusandae", "itaque", "earum", "rerum", "hic", "tenetur", "a",
        "sapiente", "delectus", "ut", "aut", "reiciendis", "voluptatibus",
        "maiores", "doloribus", "asperiores", "repellat"
    )


text:
@testset "Sparse ambiguity methods" begin
    mat = sparse(diagm([1,2,3]))
    v = [2,3,4]
    tv = v.'

    @test (tv/mat)::TransposedVector ≈ [2/1  3/2  4/3]

    @test_throws Exception mat\tv
end


text:
using DataFrames
using YAML

type Stream
  source::AbstractString
  sink::AbstractString
  process::AbstractString
end
function Stream(s::AbstractString)
  source, pt = split(s, "-{")
  process, sink = split(pt, "}->")
  return(Stream(source, sink, process))
end

function load_streams(streamfile::AbstractString)
  return([Stream(s) for s in YAML.load(open(streamfile))["stream"]])
end

function evaluate(S::Stream)
  source = S.source
  sink = S.sink
  process = S.process
  eval(:(parse("$process")(source, sink)))
end

function copyfile(infile, outfile)
  a = readtable(infile)
  writetable(outfile, a)
end

S = load_streams("stream.yml")[1]

max_stars_repo_name: el-uhu/PPSIM

max_issues_repo_name: el-uhu/PPSIM

max_forks_repo_name: el-uhu/PPSIM


text:
function detrend!(x::Vector, order=0, t = 1:length(x))
    x[:] .-= mean(x)
    if order == 1
        k = x\t
        x[:] .-= k*t
    end
    x
end

"""detrend(x, order=0, t = 1:length(x))
Removes the trend of order `order`, i.e, mean and, if order=1, the slope of the signal `x`
If `order` = 1, then the sampling points of `x` can be supplied as `t` (default: `t = 1:length(x)`)
"""
function detrend(x,args...)
    y = copy(x)
    detrend!(y,args...)
end

"""basis_activation_func(V,Nv,normalize,coulomb)

Returns a func v->ϕ(v) ∈ ℜ(Nv) that calculates the activation of `Nv` basis functions spread out to cover V nicely. If coulomb is true, then we get twice the number of basis functions, 2Nv
"""
function basis_activation_func(V,Nv,normalize,coulomb)
    if coulomb # If Coulomb setting is activated, double the number of basis functions and clip the activation at zero velocity (useful for data that exhibits a discontinuity at v=0, like coulomb friction)
        vc      = range(0, stop=maximum(abs.(V)), length=Nv+2)
        vc      = vc[2:end-1]
        vc      = [-vc[end:-1:1]; vc]
        Nv      = 2Nv
        gamma   = Nv/(abs(vc[1]-vc[end]))
        K       = normalize ? V -> _Kcoulomb_norm(V,vc,gamma) : V -> _Kcoulomb(V,vc,gamma) # Use coulomb basis function instead
    else
        vc      = range(minimum(V), stop=maximum(V), length=Nv)
        gamma   = Nv/(abs(vc[1]-vc[end]))
        K       = normalize ? V -> _K_norm(V,vc,gamma) : V -> _K(V,vc,gamma)
    end
end

"""`ridgereg(A,b,λ)`
Accepts `λ` to solve the ridge regression problem using the formulation `[A;λI]\\[b;0]. λ should be given with the same dimension as the columns of A, i.e. if λ represents an inverse standard deviation, then 1/λ = σ, not 1/λ = σ²`"""
function ridgereg(A,b,λ)
    n = size(A,2)
    [A; λ*I]\[b;zeros(n)]
end

"""`real_complex_bs(A,b, λ=0)`
Replaces the backslash operator For complex arguments. Expands the A-matrix into `[real(A) imag(A)]` and performs the computation using real arithmetics. Optionally accepts `λ` to solve the ridge regression problem using the formulation `[A;λI]\\[b;0]. λ should be given with the same dimension as the columns of A, i.e. if λ represents a standard deviation, then λ = σ, not λ = σ²`
"""
function real_complex_bs(A,b, λ=0)
    n  = size(A,2)
    Ar = [real(A) imag(A)]
    xr = λ > 0 ? [Ar; λ*I]\[b;zeros(2n)] : Ar\b
    x  = complex.(xr[1:n], xr[n+1:end])
end

function fourier_solve(A,y,zerofreq,λ=0)
    n  = size(A,2)
    x = λ > 0 ? [A; λ*I]\[y;zeros(n)] : A\y
    fourier2complex(x,zerofreq)
end

function fourier2complex(x,zerofreq)
    n = length(x)÷2
    if zerofreq === nothing
        return complex.(x[1:n], x[n+1:end])
    else
        x0 = x[zerofreq]
        x = deleteat!(copy(x), zerofreq)
        x = complex.(x[1:n], x[n+1:end])
        insert!(x,zerofreq,x0)
        return x
    end
end

""" Returns params as a [nω × N] matrix"""
reshape_params(x,Nf) = reshape(x, Nf,round(Int,length(x)/Nf))

## Complex Normal
import Base.rand

mutable struct ComplexNormal{T<:Complex}
    m::AbstractVector{T}
    Γ::Cholesky
    C::Symmetric{T}
end

function ComplexNormal(X::AbstractVecOrMat,Y::AbstractVecOrMat)
    @assert size(X) == size(Y)
    mc  = complex.(mean(X,dims=1)[:], mean(Y,dims=1)[:])
    V   = Symmetric(cov([X Y]))
    Γ,C = cn_V2ΓC(V)
    ComplexNormal(mc,Γ,C)
end

function ComplexNormal(X::AbstractVecOrMat{T}) where T<:Complex
    ComplexNormal(real.(X),imag.(X))
end

function ComplexNormal(m::AbstractVector{T},V::AbstractMatrix{T}) where T<:Real
    n   = Int(length(m)/2)
    mc  = complex.(m[1:n], m[n+1:end])
    Γ,C = cn_V2ΓC(V)
    ComplexNormal(mc,Γ,C)
end

function ComplexNormal(mc::AbstractVector{Tc},V::AbstractMatrix{Tr}) where {Tr<:Real, Tc<:Complex}
    Γ,C = cn_V2ΓC(V)
    ComplexNormal(mc,Γ,C)
end

function cn_V2ΓC(V::Symmetric{T}) where T<:Real
    n   = size(V,1)÷2
    Vxx = V[1:n,1:n]
    Vyy = V[n+1:end,n+1:end]
    Vxy = V[1:n,n+1:end]
    Vyx = V[n+1:end,1:n]
    Γ   = cholesky(complex.(Vxx + Vyy, Vyx - Vxy))
    C   = Symmetric(complex.(Vxx - Vyy, Vyx + Vxy))
    Γ,C
end

cn_V2ΓC(V::AbstractMatrix{T}) where {T<:Real} = cn_V2ΓC(Symmetric(V))

@inline cn_Vxx(Γ,C) = cholesky(real.(Matrix(Γ)+C)/2)
@inline cn_Vyy(Γ,C) = cholesky(real.(Matrix(Γ)-C)/2)
@inline cn_Vxy(Γ,C) = cholesky(imag.(-Matrix(Γ)+C)/2)
@inline cn_Vyx(Γ,C) = cholesky(imag.(Matrix(Γ)+C)/2)

@inline cn_fVxx(Γ,C) = real.(Matrix(Γ)+C)/2
@inline cn_fVyy(Γ,C) = real.(Matrix(Γ)-C)/2
@inline cn_fVxy(Γ,C) = imag.(-Matrix(Γ)+C)/2
@inline cn_fVyx(Γ,C) = imag.(Matrix(Γ)+C)/2

@inline cn_Vs(Γ,C) = cn_Vxx(Γ,C),cn_Vyy(Γ,C),cn_Vxy(Γ,C),cn_Vyx(Γ,C)
@inline cn_fV(Γ,C) = [cn_fVxx(Γ,C) cn_fVxy(Γ,C); cn_fVyx(Γ,C) cn_fVyy(Γ,C)]
@inline cn_V(Γ,C) = cholesky(cn_fV(Γ,C))
@inline Σ(cn::ComplexNormal) = Matrix(cn.Γ) # TODO: check this

for f in [:cn_Vxx,:cn_Vyy,:cn_Vxy,:cn_Vyx,:cn_fVxx,:cn_fVyy,:cn_fVxy,:cn_fVyx,:cn_Vs,:cn_V,:cn_fV]
    @eval ($f)(cn::ComplexNormal) = ($f)(cn.Γ,cn.C)
end

"""
`f(cn::ComplexNormal, z)`

Probability density Function `f(z)` for a complex normal distribution.
This can probably be more efficiently implemented
"""
function pdf(cn::ComplexNormal, z)
    k = length(cn.m)
    R = conj(cn.C)'*inv(cn.Γ)
    P = Matrix(cn.Γ)-R*cn.C # conj(Γ) = Γ for Γ::Cholesky
    cm = conj(cn.m)
    cz = conj(z)
    zmm = z-cn.m
    czmm = cz-cm
    ld = [czmm' zmm']
    rd = [zmm; czmm]
    S = [Matrix(cn.Γ) cn.C;conj(cn.C) Matrix(cn.Γ)] # conj(Γ) = Γ for Γ::Cholesky
    1/(π^k*sqrt(det(cn.Γ)*det(P))) * exp(-0.5* ld*(S\rd))
end

affine_transform(cn::ComplexNormal, A,b) = ComplexNormal(A*cn.m+b, cholesky(A*Matrix(cn.Γ)*conj(A')), Symmetric(A*cn.C*A'))

function rand(cn::ComplexNormal,s::Integer)
    L = cn_V(cn).U
    m = [real(cn.m); imag(cn.m)]
    n = length(cn.m)
    z = (m' .+ randn(s,2n)*L)
    return complex.(z[:,1:n],z[:,n+1:end])
end


text:
using PIEMetData, TidalFluxExampleData, TidalFluxConfigurations, Base.Dates
using Base.Test

TidalFluxConfigurations.config[:_METDATA_DIR] = Pkg.dir("TidalFluxExampleData","data","met")
@testset "parsemet" begin
    M = parsemet(2016)
    @test year.(M[:DateTime])[1] == 2016
    @test year.(M[:DateTime])[end] == 2017
end

@testset "Error handling" begin
    @test_throws ErrorException parsemet(2010)
end


text:
@testset "Discrete identity system" begin
    for sd in 1:3
        s = DiscreteIdentitySystem(sd)
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Discrete constrained identity system" begin
    for sd in 1:3
        X = Singleton(ones(sd))
        s = ConstrainedDiscreteIdentitySystem(sd, X)
        @test statedim(s) == sd
        @test stateset(s) == X
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Discrete linear system" begin
    for sd in 1:3
        s = LinearDiscreteSystem(zeros(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Discrete linear control system" begin
    for sd in 1:3
        s = LinearControlDiscreteSystem(zeros(sd, sd), ones(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == sd
        @test islinear(s) && isaffine(s)
    end
end

@testset "Discrete constrained linear system" begin
    A = [1. 1; 1 -1]
    X = Line([1., -1], 0.) # line x = y
    s = ConstrainedLinearDiscreteSystem(A, X)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 0
    @test islinear(s) && isaffine(s)
end

@testset "Discrete constrained linear control system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    X = Line([1., -1], 0.)
    U = Interval(0.9, 1.1)
    s = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 1
    @test inputset(s) == U
    @test islinear(s) && isaffine(s)
end

@testset "Discrete linear algebraic system" begin
    for sd in 1:3
        s = LinearAlgebraicDiscreteSystem(zeros(sd, sd), zeros(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Discrete constrained linear algebraic system" begin
    A = [1. 1; 1 -1]
    E = [0. 1; 1 0]
    X = LinearConstraint([0, -1.], 0.) # the set y ≥ 0
    s = ConstrainedLinearAlgebraicDiscreteSystem(A, E, X)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 0
    @test islinear(s) && isaffine(s)
end

@testset "Constant input in a discrete constrained linear control system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    X = Line([1., -1], 0.)
    U = ConstantInput(Hyperrectangle(low=[0.9], high=[1.1]))
    s = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    s = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    @test islinear(s) && isaffine(s)
    if VERSION < v"0.7-"
        @test next(s.U, 1)[1] isa Hyperrectangle
    else
        @test iterate(s.U)[1] isa Hyperrectangle
    end
end

@testset "Varying input in a discrete constrained linear control system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    X = Line([1., -1], 0.)
    U = VaryingInput([Hyperrectangle(low=[0.9], high=[1.1]),
                      Hyperrectangle(low=[0.99], high=[1.0])])
    s = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    @test length(s.U) == 2
    for ui in s.U
        @test ui isa Hyperrectangle
    end
    @test islinear(s) && isaffine(s)
end

@testset "Polynomial system in discrete time" begin
    @polyvar x y
    p = 2x^2 - 3x + y
    s = PolynomialDiscreteSystem(p)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test !islinear(s) && !isaffine(s)
end

@testset "Polynomial system in discrete time with state constraints" begin
    @polyvar x y
    p = 2x^2 - 3x + y
    X = BallInf(zeros(2), 0.1)
    s = ConstrainedPolynomialDiscreteSystem(p, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test dim(stateset(s)) == dim(X)
    @test !islinear(s) && !isaffine(s)
end

@testset "Implicit discrete system" begin
    add_one(x) = x + 1
    s = BlackBoxDiscreteSystem(add_one, 1)
    x = 1.0
    @test s.f(x) ≈ 2.0
end


text:
# # TEST FV32: Cantilevered tapered membrane, free vibration

# Source code: [`FV32_tut.jl`](FV32_tut.jl)

# ## Description

# FV32: Cantilevered tapered membrane is a test recommended by the National
# Agency for Finite Element Methods and Standards (U.K.): Test FV32 from NAFEMS
# publication TNSB, Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.

# Reference solution: 44.623, 130.03, 162.70, 246.05, 379.90, 391.44 Hz for the
# first six modes.

# The benchmark is originally for plane stress conditions. We simulate the
# plane-stress conditions with a three-dimensional mesh that is constrained
# along one plane of nodes to effect the constrained motion only in the plane
# of the trapezoidal membrane.

# ![](FV32-mesh.png)

# ## References

# [1] Test FV32 from NAFEMS publication TNSB, Rev. 3, “The Standard NAFEMS
# Benchmarks,” October 1990.

# ## Goals

# - Show how to generate hexahedral mesh in a rectangular block  and shape it
#   into a trapezoid.
# - Set up model data for the solution algorithms.
# - Use two different finite element model machines to evaluate the stiffness
#   and the mass.
# - Execute the modal algorithm and export the results with another algorithm.

##
# ## Definitions

# Bring in the required support from the basic linear algebra, eigenvalue
# solvers, and the finite element tools.
using LinearAlgebra
using Arpack
using FinEtools
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear: AlgoDeforLinearModule

# The input data is given by the benchmark.
E = 200*phun("GPA");
nu = 0.3;
rho= 8000*phun("KG/M^3");
L = 10*phun("M");
W0 = 5*phun("M");
WL = 1*phun("M");
H = 0.05*phun("M");

# We shall generate a three-dimensional mesh. It should have 1 element through
# the thickness, and 8 and 4 elements in the plane of the membrane.
nL, nW, nH = 8, 4, 1;# How many element edges per side?
# The reference frequencies are obtained from [1].
Reffs = [44.623 130.03  162.70  246.05  379.90  391.44]

# The three-dimensional mesh of 20 node serendipity hexahedral should correspond
# to the plane-stress quadratic serendipity quadrilateral (CPS8R) used in the
# Abaqus benchmark. We simulate the plane-stress conditions with a
# three-dimensional mesh that is constrained along one plane of nodes to effect
# the constrained motion only in the plane of the trapezoidal membrane. No
# bending out of plane!
# First we generate mesh of a rectangular block.
fens,fes = H20block(1.0, 2.0, 1.0, nL, nW, nH)
# Now distort the rectangular block into the tapered plate.
for i in 1:count(fens)
    xi, eta, theta = fens.xyz[i,:];
    eta = eta - 1.0
    fens.xyz[i,:] = [xi*L eta*(1.0 - 0.8*xi)*W0/2 theta*H/2];
end
# We can visualize the mesh with Paraview (for instance).
File =  "FV32-mesh.vtk"
vtkexportmesh(File, fens, fes)
@async run(`"paraview.exe" $File`)

# The simulation will be executed with the help of algorithms defined in the
# package `FinEtoolsDeforLinear`. The algorithms accept a dictionary of model
# data. The model data dictionary will be built up as follows.

# First we make the interior region. The model reduction is for a three-dimensional finite element model.
MR = DeforModelRed3D
material = MatDeforElastIso(MR, rho, E, nu, 0.0)

# We shall create two separate finite element model machines. They are
# distinguished by the quadrature rule. The mass rule, in order to evaluate the
# mass matrix accurately, needs to be of higher order than the one we prefer to
# use for the stiffness.
region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material), "femm_mass"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material))

# Select nodes that will be clamped.
nl1 = selectnode(fens; plane=[1.0 0.0 0.0 0.0], thickness=H/1.0e4)
ebc1 = FDataDict("node_list"=>nl1, "component"=>1, "displacement"=>0.0)
ebc2 = FDataDict("node_list"=>nl1, "component"=>2, "displacement"=>0.0)
ebc3 = FDataDict("node_list"=>nl1, "component"=>3, "displacement"=>0.0)

# Export  a VTK file to visualize the selected points. Choose the
# representation "Points", and select color and size approximately 4. These
# notes should correspond to the clamped base of the membrane.
File =  "FV32-nl1.vtk"
vtkexportmesh(File, fens, FESetP1(reshape(nl1, length(nl1), 1)))

# Select all nodes on the plane Z = 0. This will be prevented from moving in the
# Z direction.
nl4 = selectnode(fens; plane=[0.0 0.0 1.0 0.0], thickness=H/1.0e4)
ebc4 = FDataDict("node_list"=>nl4, "component"=>3, "displacement"=>0.0)

# Export  a VTK file to visualize the selected points. Choose the
# representation "Points", and select color and size approximately 4. These
# points all should be on the bottom face of the three-dimensional domain.
File =  "FV32-nl4.vtk"
vtkexportmesh(File, fens, FESetP1(reshape(nl4, length(nl4), 1)))

# Make model data: the nodes, the regions, the boundary conditions, and the
# number of eigenvalues are set. Note that the number of eigenvalues needs to
# be set to 6+N,  where 6 is the number of rigid body modes, and N is the
# number of deformation frequencies we are interested in.
neigvs = 10                  # how many eigenvalues
modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "essential_bcs"=>[ebc1 ebc2 ebc3 ebc4], "neigvs"=>neigvs)

# Solve using an algorithm: the modal solver. The solver will supplement the
# model data with the geometry and displacement fields, and the solution
# (eigenvalues, eigenvectors), and the data upon return can be extracted from
# the dictionary.
modeldata = AlgoDeforLinearModule.modal(modeldata)

# Here we extract the angular velocities corresponding to the natural frequencies.
fs = modeldata["omega"]/(2*pi)
println("Eigenvalues: $fs [Hz]")
println("Percentage frequency errors: $((vec(fs[1:6]) - vec(Reffs))./vec(Reffs)*100)")

# The problem was solved for instance with Abaqus, using plane stress eight node
# elements. The results were:

# | Element | Frequencies (relative errors) |
# | ------- | ---------------------------- |
# | CPS8R  | 44.629 (0.02)   130.11 (0.06)   162.70 (0.00)   246.42 (0.15)   381.32 (0.37)   391.51 (0.02) |

# Compared these numbers with those computed by our three-dimensional model.

# The mode shapes may be visualized with `paraview`. Here is for instance mode
# 8:

# ![](FV32-mode-8.png)

# The algorithm to export the mode shapes expects some input. We shall specify
# the filename and the numbers of modes to export.
modeldata["postprocessing"] = FDataDict("file"=>"FV32-modes", "mode"=>1:neigvs)
modeldata = AlgoDeforLinearModule.exportmode(modeldata)
# The algorithm attaches a little bit to the name of the exported file. If
# `paraview.exe` is installed, the command below should bring up the
# postprocessing file.
@async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

# To animate the mode shape in `Paraview` do the following:
# - Apply the filter "Warp by vector".
# - Turn on the "Animation view".
# - Add the mode shape data set ("WarpByVector1") by clicking the "+".
# - Double-click the line with the data set. The "Animation Keyframes" dialog
#   will come up. Double-click "Ramp" interpolation, and change it
#   to "Sinusoid". Set the frequency to 1.0. Change the "Value" from 0 to 100.
# - In the animation view, set the mode to "Real-time", and the duration to 4.0
#   seconds.
# - Click on the "Play" button. If you wish, click on the "Loop" button.

true


text:
######################################################################
# GLFWAbstraction central include hub
# -----
# Licensed under MIT License
module GLFWAbstraction

const Optional{T} = Union{T, Nothing}

abstract type GLFWWrapper{T} end
wrapped(w::GLFWWrapper) = w.handle
wrapped_type(::Type{GLFWWrapper{T}}) where T = T

"""`lhs×rhs` creates a 2-tuple (lhs, rhs) for a familiar and convenient notation of 2D measures."""
(×)(lhs::Real, rhs::Real) = (lhs, rhs)

include("./monitors.jl")
include("./windows.jl")
include("./input.jl")

end # module GLFWAbstraction


text:
function isotonic_regression!(y::Vector{Float64}, weights::Vector{Float64})

    n = length(y)
    if n <= 1
        return y
    end
    if n != length(weights)
        throw(DimensionMismatch("Lengths of values and weights mismatch"))
    end

    @inbounds begin
        n -= 1
        while true
            i = 1
            pooled = 0
            while i <= n
                k = i
                while k <= n && y[k] >= y[k+1]
                    k += 1
                end

                # Find a decreasing subsequence, and update
                # all points in the sequence to the weighted average.
                if y[i] != y[k]
                    numerator = 0.0
                    denominator = 0.0
                    for j in i : k
                        numerator += y[j] * weights[j]
                        denominator += weights[j]
                    end

                    for j in i : k
                        y[j] = numerator / denominator
                    end
                    pooled = 1
                end
                i = k + 1
            end
            if pooled == 0
                break
            end
        end
    end
    return y
end

isotonic_regression!(y::Vector{Float64}) = isotonic_regression!(y, ones(size(y, 1)))

# non-mutating versions
isotonic_regression(y::Vector{Float64}, weights::Vector{Float64}) = isotonic_regression!(copy(y), weights)
isotonic_regression(y::Vector{Float64}) = isotonic_regression(y, ones(size(y,1)))


text:
@testset "Heuristic tests" begin
    # Load a graph to test on.
    g = graph_from_gr("../examples/test.gr")

    # Test the min fill heuristic.
    order, tw = min_fill(g)
    @test length(order) == lg.nv(g)

    # Test min fill sampling
    order, tw = min_fill(g, 10)
    @test length(order) == lg.nv(g)

    # # Test the min width heuristic.
    # order, tw = min_width(g)
    # @test length(order) == lg.nv(g)

    # # Test sampling version of min width.
    # order, tw = min_width(g, 10)
    # @test length(order) == lg.nv(g)

    # Test maximum minimum degree heuristic for a lower bound.
    lb = mmd(g)
    @test lb <= tw
end


text:
@inline _F(t, x, P) = B(t, P)*x + β(t, P)
@inline _dHinv(t, K, P) = B(t, P)*K + K*B(t, P)' - a(t, P)
@inline _dK(t, K, P) = B(t, P)*K + K*B(t, P)' + a(t, P)
@inline _dPhi(t, Phi, P) = B(t, P)*Phi

"""
    gpHinv!(K::SamplePath, P, KT=zero(T))

Precompute ``K = H^{-1}`` from ``(d/dt)K = BK + KB' + a`` for a guided proposal.
"""
gpHinv!(K::SamplePath{T}, P, KT=zero(T)) where {T} = _solvebackward!(R3(), _dHinv, K, KT, P)
gpH♢! = gpHinv!

"""
gpV!(K::SamplePath, P, KT=zero(T))

Precompute `V` from ``(d/dt)V = BV + β``, ``V_T = v`` for a guided proposal.
"""
gpV!(V::SamplePath{T}, P, v::T) where {T} = _solvebackward!(R3(), _F, V, v, P)

gpmu(tt, u, P) = solve(R3(), _F, tt, u, P)
gpK(tt, u, P) = solve(R3(), _dK, tt, u, P)

"""
    fundamental_matrix(tt, P)

Compute fundamental solution.
"""
fundamental_matrix(tt, P, Phi0 = one(outertype(P))) = solve(R3(), _dPhi, tt, Phi0, P)


text:
abstract type AbstractModel end

struct PoissonLapseUniform <: AbstractModel
    params::Vector{Float64}
    max::Int
end
function init(::Type{<:PoissonLapseUniform}, cm::Dict)
    m = maximum(key.Omissions_plus_one for key in keys(cm))
    return PoissonLapseUniform([10, 2, 0.1], m)
end

Distributions.params(p::PoissonLapseUniform) = p.params
Base.maximum(m::PoissonLapseUniform) = m.max

function Distributions.cdf(p::PoissonLapseUniform, x::Int)
    T, r, ϵ = params(p)
    d = Gamma(T, 1/r)
    # ϵ = 0
    # 1-ϵ e faccio Gamma, ϵ e faccio uniform, sommare e' uguale a "oppure" in probabilita'
    return (1-ϵ)*cdf(d, x) + ϵ*x/p.max
end

struct RhoComparison <: AbstractModel
    params::Vector{Float64}
    env_rho::Float64
    patch_rho::Float64
end

function init(::Type{<:RhoComparison}, df)
    final_env_rho = sum(df.REWARD/df[end,:TIME])
    return RhoComparison([
        rand(0.0:0.0001:0.1), #=alpha_environment = [0,0.1] learning rate of environment=#
        rand(0.0:0.0001:0.2), #=alpha_patch = [0, 0.2] learning rate of the patch=#
        rand(2:0.0001:10), #=beta = [2, 10] noise parameter or inverse temperature=#
        rand(0.0:0.0001:0.5), #=reset = [0, 0.5] patch env rew rate at the entrance of a trials, because it has decayed in the last trial=#
        0.1 #=bias = [2, 10] bias towards staying in the softmax=#
        ], final_env_rho, final_env_rho)
end

Distributions.params(m::RhoComparison) = m.params
# you need to calculate the likelihood for the entire session, then sum and minimize -log
# function Log_LikeliHood(m::AbstractModel, rew, leaves)
#
# end


text:
# viewdef generation

using Mustache

const GEN_FIELDS = """

{{#fields}}
  {{{.}}}
{{/fields}}

"""

const GEN_TUPLE = """

{{#fields}}
  {{{.}}}
{{/fields}}

"""

const GEN_SHOW = """

"""

get_path(cmd::VarCommand, var::Var) = get_path(VarCommand{:get, ()}(cmd; var))

struct Field{owner, namespace, field, field_type, first}
end

@kwdef mutable struct ListEditor{Namespace}
    id = 1
    name = ""
    list
    selection = nothing
    #add_enabled::Function = (; fields...)-> true
    #check_enabled::Function
    #add_tooltip::AbstractString = ""
end

function ListEditor(con, list::T; args...) where T
    if !haskey(args, :id)
        args = (; args..., id = "list-$(con === nothing ? 1 : con.sequence += 1)")
    end
    ListEditor(; list, args...)
end

# used for creation
list_type(::T) where T = base_list_type(T)
list_type(::Type{T}) where T = base_list_type(T)
list_type(ed::ListEditor) = base_list_type(typeof(ed.list))

function base_list_type(::Type{T}) where T
    et = eltype(T)
    et == T ? nothing : et
end

#add_item(ed::ListEditor) =

add_enabled(ed::ListEditor) = list_type(ed) !== nothing && ed.add_enabled()

edit_enabled(ed::ListEditor) = ed.selection !== nothing

selection_index(ed::ListEditor) = something(findfirst(x-> x == ed.selection, ed.list), 0)
selection_index(ed::ListEditor, index::Number) =
    selection(ed, get(ed.list, index, nothing))

selection(ed::ListEditor, item) = ed.selection = item

show(cmd::VarCommand, ed::ListEditor) = foreach(e-> show(cmd, e), ed.fields)

function generate_viewdef(cmd::VarCommand)
    var = cmd.var
    (var.internal_value === nothing || !haskey(var.metadata, :genview)) && return
    namespace = var.metadata[:genview]
    viewdef = generate_viewdef(cmd, typeof(var.internal_value), Symbol(namespace))
    println("GENERATED VIEWDEF:\n", viewdef)
    set_metadata(cmd, :viewdef, viewdef)
    if cmd.config.output_dir != ""
        filebase = "$(typename(cmd.var.internal_value))"
        namespace != "" && (filebase *= "-$namespace")
       open(joinpath(cmd.config.output_dir, "$filebase.html"); write=true) do io write(io,  viewdef) end
    end
end

"""
    generate_viewdef(item, namespace)

generate a viewdef for item
namespace indicates what type of view to generate (default is editor, could
be list, link, etc.)
"""
function generate_viewdef(cmd::VarCommand, ::Type{T}, ns::Symbol = Symbol("")) where T
    !isstructtype(T) && return GEN_SHOW
    render(GEN_FIELDS, (; fields = [gen_field(cmd, Field{T, ns, field, fieldtype(T, field), i == 1})
                                    for (i, field) in enumerate(fieldnames(T))]))
end

function template(cmd::VarCommand, template, namespace)
    k = (Symbol(template), Symbol(namespace))
    if !haskey(cmd.config.templates, k)
        name = string(namespace) == "" ? template : "$template-$namespace"
        cmd.config.templates[k] = read(joinpath(cmd.config.templates_dir, "$name.html"), String)
    end
    cmd.config.templates[k]
end

function all_presentable(t::Type)
####### HERE
end

#generate_viewdef(cmd::VarCommand, ed::ListEditor; ns) =
#    Mustache.render(template(cmd, :list, ns), tmpldata(ed))

"""
    gen_field(cmd, owner_type, Field{owner, namespace, field, field_type, first})

Generate a field; override this to customize a particular field or type of field.

Generated field types: Vector, Enum, AbstractString, Number
"""
# radio button
function gen_field(cmd::VarCommand, fld::Type{Field{owner, ns, field, field_type, first}}) where {
    owner, ns, field, field_type <: Enum, first
}
    println("FIELD TYPE: $(field_type)")
    render(template(cmd, :radio, ns),
           (;
            introspect(fld)...,
            type = typename(field_type),
            id = cmd.connection.sequence += 1,
            values = [(;name = pretty(e), value = string(e)) for e in instances(field_type)]))
end

# text field
gen_field(cmd::VarCommand, fld::Type{Field{owner, ns, field, field_type, first}}) where {
    owner, ns, field, field_type <: Union{AbstractString, Number}, first
} =
    render(template(cmd, :textfield, ns), introspect(fld))

# checkbox
gen_field(cmd::VarCommand, fld::Type{Field{owner, ns, field, field_type, first}}) where {
    owner, ns, field, field_type <: Bool, first
} =
    render(template(cmd, :checkbox, ns),
           (;
            introspect(fld)...,
            type = typename(field_type)))

#gen_field(cmd::VarCommand, fld::Type{Field{owner, ns, field, field_type, first}}) where {
#    owner, ns, field, field_type <: Tuple, first
#} =
#    render(GEN_FIELDS,
#           field,
#           fields

#list
#gen_field(cmd::VarCommand, fld::Type{Field{owner, ns, field, field_type, first}}) where {
#    ns, owner, field, field_type <: Vector, first
#} =
#    render(template(cmd, :list_field, ns),
#           (; gen_introspect(fld)...,
#            name = pretty(field),
#            first))

"""
    introspect(fld::Type{Field})

computes PATH, PRETTY_NAME, and FIRST for fld
returns (; field = PATH, name = PRETTY_NAME, first = FIRST)
"""
function introspect(::Type{Field{owner, _ns, field, ft, first}}) where {owner, _ns, field, ft, first}
    try
        println("CHECKING METHODS FOR $field")
        func = parentmodule(owner).eval(field)
        println("FUNC: $func$((hasmethod(func, (owner,)) ||
            hasmethod(func, (VarCommand, owner))) ? " [getter]" : "")")
        (hasmethod(func, (owner,)) ||
            hasmethod(func, (VarCommand, owner))) && return (; field = "$field()")
    catch
    end
    (; field, name = pretty(field), first)
end

tmpldata(ed::ListEditor) = to_dict((;
    id = ed.id,
    name = ed.name,
    fields = [(;
               first = i == 1,
               name = pretty(field),
               field = field,
               ) for (i, field) in enumerate(fieldnames(list_type(ed)))]))

render(tmpl, value) = Mustache.render(tmpl, to_dicts(value))

to_dicts(val) = val
to_dicts(val::Array) = map(to_dicts, val)
to_dicts(tup::NamedTuple) = Dict((string(k)=>to_dicts(v) for (k, v) in pairs(tup))...)

# capitalize first letter, use spaces instead of underbars and camelcase transitions
function pretty(val)
    str = string(val)
    str[1] * replace(str[2:end], "_"=>" ", r"(\p{Ll})(\p{Lu}+)"=>(l, u)-> "$l $(lowercase.(split(u, "")))")
end

function edit_list_field()
end


text:
# The purpose of this script is to bootstrap a library of planning subproblems
# that are representative of (what should be) a high performance planner

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics

close()

data_file = "./data/rollout_library.jld2"
reprocess = false

steps = 60
num_robots = [4, 6, 8, 10]

trials = 1:20

horizon = SubmodularMaximization.default_horizon

# Keep samples from only a few (8) relevant time-steps in every trial
sample_steps = 20:5:steps

# At this time, the appropriate number of samples for convergence is unknown
# Using a relatively large number can improve results
num_iterations = 3000

solver = solve_sequential

function get_data()
  data_mutex = Mutex()
  data = Dict()

  global num_robots, trials

  if !isfile(data_file) || reprocess

    println("Running trials")

    tests = collect(product(num_robots, trials))
    @threads for (num_robots, trial) in tests

      configs = MultiRobotTargetTrackingConfigs(num_robots,
                                                horizon=horizon,
                                                solver_iterations=
                                                num_iterations
                                               )

      @time rollout = target_tracking_experiment(steps=steps,
                                           num_robots=num_robots,
                                           configs=configs,
                                           solver=solver
                                          )

      # Samples correspond to joint system states where I will wish to compare
      # solvers
      samples = rollout[sample_steps]

      lock(data_mutex)
      data[num_robots, trial] =
      (rollout=rollout, samples=samples, configs=configs)
      unlock(data_mutex)

      println("Robots ", num_robots, " Trial ", trial, " (done)")
    end

    @save data_file data num_robots trials
  else
    @load data_file data num_robots trials
  end

  data
end

@time data = get_data()

entropies = Dict(key=>map(x->entropy(x.target_filters), value.samples)
                 for (key, value) in data)

concatenated_entropy = map(num_robots) do num_robots
  vcat(map(trials) do trial
    entropies[num_robots, trial]
  end...)
end[:]

objective_values = map(num_robots) do num_robots
  @show num_robots
  vcat(map(trials) do trial
    map(x->x.objective, data[num_robots, trial].samples)
  end...)
end[:]

for (ii, num_robots) in enumerate(num_robots)
  mean_objective = mean(objective_values[ii])
  mean_entropy = mean(concatenated_entropy[ii])

  println("Robots (", num_robots, ")")
  println("  Mean entropy: ", mean_entropy,
          " (", mean_entropy / num_robots, " per robot)")
  println("  Mean objective: ", mean_objective,
          " (", mean_objective / num_robots, " per robot)")
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text: litres(time) = time÷2

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
```
    NonmaximaSuppression <: AbstractEdgeThinningAlgorithm
    NonmaximaSuppression(;) #TODO Decide on keywords

    thin_edges([T,] img, f::NonmaximaSuppression)
    thin_edges!([out,] img, f::NonmaximaSuppression)
```

Returns an image that contains a single response per edge.

# Details

TODO

# Example

```julia

using TestImages, FileIO, ImageView

img =  testimage("mandril_gray")
img_edges = detect_edges(img, Canny())
img_thin_edges = thin_edges(img_edges, NonmaximaSuppression())

imshow(img)
imshow(img_edges)
imshow(img_thin_edges)
```

# References
TODO.
"""
@with_kw struct NonmaximaSuppression{T₁ <: Union{Real,AbstractGray},
                                     T₂ <: Union{Real,AbstractGray}} <: AbstractEdgeDetectionAlgorithm
    # TO BE DECIDED
    example_keyword_1::T₁ = 0
    example_keyword_2::T₂ = 0
end

function (f::NonmaximaSuppression)(out::GenericGrayImage, img::GenericGrayImage)
    # TODO
end

function (f::NonmaximaSuppression)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    # TODO
end

(f::NonmaximaSuppression)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
export HomogeneousFock

mutable struct HomogeneousFock{D,C} <: AbstractHilbert
    n_sites::Int
    shape::Vector{Int}

    # Constrain on total number of excitations
    n_exc::Int
end

"""
    HomogeneousFock(N, m; excitations=N*m)

Constructs the Hilbert space of `N` identical modes with `m` levels.
If "n_exc" is set, then the total number of excitations is bounded when
generating a state.

!!!note
    If using a sampler that does not respect the symmetries of the system,
    the n_exc will not be respected during sampling.
"""
function HomogeneousFock(n_sites, hilb_dim; excitations = -1)
    n_sites = n_sites isa AbstractGraph ? nv(n_sites) : n_sites

    constrained = excitations >= 0 ? true : false
    if constrained && excitations > n_sites*(hilb_dim-1)
        throw(ErrorException("Constrain is useless: bigger than total number of allowed particles."))
    end
    HomogeneousFock{hilb_dim,constrained}(n_sites, fill(hilb_dim, n_sites),
                                          excitations)
end

Base.similar(hilb::HomogeneousFock, N) =
    HomogeneousFock(N, local_dim(hilb))

@inline nsites(h::HomogeneousFock) = h.n_sites
@inline local_dim(h::HomogeneousFock{D}) where D = D
@inline local_dim(h::HomogeneousFock{D}, i) where D = D
@inline shape(h::HomogeneousFock) = h.shape

@inline spacedimension(h::HomogeneousFock) = local_dim(h)^nsites(h)
@inline indexable(h::HomogeneousFock) = spacedimension(h) != 0
@inline is_homogeneous(h::HomogeneousFock) = true

@inline is_contrained(h::HomogeneousFock{H,C}) where {H,C} = C
@inline constraint_limit(h::HomogeneousFock) = h.n_exc

state(arrT::AbstractArray, T::Type{<:Number}, h::HomogeneousFock, dims::Dims) =
    similar(arrT, T, nsites(h), dims...) .= 0.0

Base.show(io::IO, ::MIME"text/plain", h::HomogeneousFock) =
    print(io, "Hilbert Space with $(nsites(h)) identical sites of dimension $(local_dim(h))")

Base.show(io::IO, h::HomogeneousFock) =
    print(io, "HomogeneousFock($(nsites(h)), $(local_dim(h)))")

## Operations

function flipat!(rng::AbstractRNG, σ::AState, h::HomogeneousFock{N}, i) where N
    @boundscheck checkbounds(σ, i)
    T = eltype(σ)

    @inbounds old_val = σ[i]
    new_val =  floor(rand(rng, T)*(N-1))
    new_val = new_val + (new_val >= old_val)
    @inbounds σ[i]    = new_val

    return old_val, new_val
end

# special case N== 2 to be faster
function flipat!(rng::AbstractRNG, σ::AState, h::HomogeneousFock{2}, i)
    @boundscheck checkbounds(σ, i)
    T = eltype(σ)

    @inbounds old_val = σ[i]
    new_val = old_val == 0.0 ? 1.0 : 0.0
    @inbounds σ[i]    = new_val

    return old_val, new_val
end

function setat!(σ::AState, h::HomogeneousFock, i::Int, val)
    @boundscheck checkbounds(σ, val)

    @inbounds old_val = σ[i]
    @inbounds σ[i] = val

    return old_val
end

set_index!(σ::AState, h::HomogeneousFock, i) = set!(σ, h, i)
function set!(σ::AState, h::HomogeneousFock{N}, val::Integer) where N
    @boundscheck checkbounds_hilbert(h, val)

    val -= 1
    for i=1:nsites(h)
        @inbounds val, σ[i] = divrem(val, N)
    end
    return σ
end

Base.@propagate_inbounds add!(σ::AState, h::HomogeneousFock, val::Integer) =
    set!(σ, h, val+toint(σ, h))

function Random.rand!(rng::AbstractRNG, σ::AbstractArray, h::HomogeneousFock{N, false}) where N
    T = eltype(σ)

    #rand!(rng, σ, 0:(N-1))
    rand!(rng, σ)
    σ .*= N
    σ .= floor.(σ)
    return σ
end

# Specialized for constrained fock spaces
function Random.rand!(rng::AbstractRNG, σ::AState, h::HomogeneousFock{N, true}) where N
    T = eltype(σ)
    n_max = constraint_limit(h)
    n_sites = nsites(h)

    σ .= zero(eltype(σ))

    # add all constraints one by one
    for i=1:constraint_limit(h)
        # select a (non full) site to which it should be added
        site = rand(rng, 1:n_sites)
        i = 0
        @inbounds while σ[site] == N-1
            i += 1
            site = rand(rng, 1:n_sites)
        end
        @inbounds σ[site] += 1
    end

    return σ
end

function toint(σ::AState, h::HomogeneousFock{N}) where N
    tot = 0
    for (i,v)=enumerate(σ)
        tot += Int(v)*N^(i-1)
    end
    return tot + 1
end

function local_index(σ::AState, h::HomogeneousFock, site)
    @boundscheck checkbounds(σ, site)

    @inbounds id = Int(σ[site])+1
    return id
end

function local_index(σ::AState, h::HomogeneousFock{M}, sites::AbstractVector) where M
    @boundscheck checkbounds(σ, sites)

    idx = 1
    for (i,j)=enumerate(sites)
        @inbounds idx += Int(σ[j]) * M^(i-1)
    end
    return idx
end


text:
using Test: @test, @testset
import Infernal
import Rfam

@testset "infernal" begin
    @test isfile(Infernal.path("cmalign"))
    @test success(Infernal.exe("cmfetch", Rfam.cm(), "RF00162"))
end


text:
using bar: bar1
using baz: baz1
function main()
    x = bar1()
    y = baz1()
    @assert(x == 0)
    @assert(y == "foo")
end

main()


text:
struct ConsistencyBars{T<:Real}
    samples::Int
    quantiles::Tuple{T,T}

    @doc """
        ConsistencyBars(; samples::Int=1_000, coverage::Real=0.95)

    Create consistency bars that cover a proportion of `coverage` samples obtained
    by consistency resampling with `samples` resampling steps.

    # References

    Bröcker, J. and Smith, L.A. (2007). *Increasing the reliability of reliability diagrams*.
    Weather and forecasting, 22(3), pp. 651-661.
    """
    function ConsistencyBars(; samples::Int=1_000, coverage::Real=0.95)
        samples > 0 || error("number of samples must be positive")
        zero(coverage) < coverage < one(coverage) ||
            error("coverage must be a number in (0, 1)")
        αo2 = (1 - coverage) / 2
        quantiles = promote(αo2, 1 - αo2)
        return new{eltype(quantiles)}(samples, quantiles)
    end
end


text:
function coordinated_turn_bank(turn_rate, alpha, beta, tas, gamma)

    g0 = GRAVITY_ACCEL
    G = turn_rate * tas / g0

    if abs(gamma) < 1e-8
        phi = G * cos(beta) / (cos(alpha) - G * sin(alpha) * sin(beta))
        phi = atan(phi)
    else
        a = 1 - G * tan(alpha) * sin(beta)
        b = sin(gamma) / cos(beta)
        c = 1 + G^2 * cos(beta)^2

        sq = sqrt(c * (1 - b^2) + G^2 * sin(beta)^2)

        num = (a - b^2) + b * tan(alpha) * sq
        den = a ^ 2 - b^2 * (1 + c * tan(alpha)^2)

        phi = atan(G * cos(beta) / cos(alpha) * num / den)
    end
    return phi
end

function climb_theta(gamma, alpha, beta, phi)
    a = cos(alpha) * cos(beta)
    b = sin(phi) * sin(beta) + cos(phi) * sin(alpha) * cos(beta)
    sq = sqrt(a^2 - sin(gamma)^2 + b^2)
    theta = (a * b + sin(gamma) * sq) / (a^2 - sin(gamma)^2)
    theta = atan(theta)
    return theta
end

function turn_rate_angular_velocity(turn_rate, theta, phi)
    # w = turn_rate * k_h
    # k_h = sin(theta) i_b + sin(phi) * cos(theta) j_b + cos(theta) * sin(phi)
    # w = p * i_b + q * j_b + r * k_b
    p = - turn_rate * sin(theta)
    q = turn_rate * sin(phi) * cos(theta)
    r = turn_rate * cos(theta) * cos(phi)
    return [p, q, r]
end

"""
    body_angular_velocity_to_euler_angles_rates(p, q, r, ψ, θ, ϕ)

Transform body angular velocity (p, q, r) [rad/s] to Euler angles rates
(ψ_dot, θ_dot, ϕ_dot) [rad/s] given the euler angles (θ, ϕ) [rad] using
kinematic angular equations.

# References

- [1] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. Equation (1.4-4) (page 20)
"""
function body_angular_velocity_to_euler_angles_rates(p, q, r, θ, ϕ)

    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    ψ_dot = (q * sϕ + r * cϕ) / cθ
    θ_dot = q * cϕ - r * sϕ
    # ϕ_dot = p + (q * sϕ + r * cϕ) * tan(θ)
    ϕ_dot = p + ψ_dot * sθ

    return [ψ_dot, θ_dot, ϕ_dot]
end

"""
    euler_angles_rates_to_body_angular_velocity(ψ_dot, θ_dot, ϕ_dot, ψ, θ, ϕ)

Transform Euler angles rates (ψ_dot, θ_dot, ϕ_dot) [rad/s] to body angular
velocity (p, q, r) [rad/s] given the euler angles (θ, ϕ) [rad] using
kinematic angular equations.

# References

- [1] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. Equation (1.4-3) (page 20)
"""
function euler_angles_rates_to_body_angular_velocity(ψ_dot, θ_dot, ϕ_dot, θ, ϕ)

    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    p = ϕ_dot - sθ * ψ_dot
    q = cϕ * θ_dot + sϕ*cθ * ψ_dot
    r = -sϕ * θ_dot + cϕ*cθ * ψ_dot

    return [p, q, r]
end

"""
    body_angular_velocity_to_quaternion_rates(p, q, r, q0, q1, q2, q3)

Transform body angular velocity (p, q, r) [rad/s] to quaternion rates [1/s].

- [1] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. Equation (1.8-15) (page 51)

"""
function body_angular_velocity_to_quaternion_rates(p, q, r, q0, q1, q2, q3)

    Ω = [
        0 -p -q -r;
        p  0  r -q;
        q -r  0  p;
        r  q -p  0;
    ]
    q = [q0; q1; q2; q3]
    q_dot = 0.5 * Ω * q

    return q_dot
end

"""
    uvw_to_tasαβ(u, v, w)

Calculate true air speed (TAS), angle of attack (α) and angle of side slip (β)
from velocity expressed in body axis.

# Notes
This function assumes that u, v, w are the body components of the aerodynamic
speed. This is not true in genreal (wind speed different from zero), as u, v, w
represent velocity with respect to an inertial reference frame.

# References

- [1] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. Equation (2.3-6b) (page 78)
"""
function uvw_to_tasαβ(u, v, w)

    tas = sqrt(u*u + v*v + w*w)
    α = atan(w, u)
    β = asin(v / tas)
    return [tas, α, β]
end

"""
    tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

Calculate time derivatives of velcity expressed as TAS, AOA, AOS.

# Notes
Note that tas here is not necessarily true air speed. Could also be inertial speed in the
direction of airspeed. It will concide whith TAS for no wind.

# See also
    uvw_dot_from_tas_α_β(tas, α, β, tas_dot, α_dot, β_dot)

# References
- [1] Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and
 practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Equation 3.33 (page 44).
- [2] Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
 and simulation: dynamics, controls design, and autonomous systems. John Wiley
 & Sons. Equation (2.3-10) (page 81)
"""
function tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

    # [1] 3.31
    tas = sqrt(u*u + v*v + w*w)
    # [1] 3.33a
    tas_dot = (u * u_dot + v * v_dot + w * w_dot) / tas
    # [1] 3.33b
    β_dot = ((u*u + w*w) * v_dot - v * (u * u_dot + w * w_dot)) / (tas^2 * sqrt(u*u + w*w))
    # [2] 2.3.10b
    # β_dot = (v_dot * tas - v * tas_dot) / (tas * sqrt(u*u + w*w))
    # [1] 3.33c
    α_dot = (w_dot * u - w * u_dot) / (u*u + w*w)

    return [tas_dot, α_dot, β_dot]
end

"""
    uvw_dot_from_tas_α_β(tas, α, β, tas_dot, α_dot, β_dot)

Obatain body velocity derivatives given velocity in wind axis and its derivatives.

# See also
    tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

# References
- [1] Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and
 practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Derived from equation 3.32 (page 44).
"""
function uvw_dot_from_tas_α_β(tas, α, β, tas_dot, α_dot, β_dot)
    u_dot = tas_dot * cos(α) * cos(β) - tas * (α_dot * sin(α) * cos(β) + β_dot * cos(α) * sin(β))
    v_dot = tas_dot * sin(β) + tas * β_dot * cos(β)
    w_dot = tas_dot * sin(α) * cos(β) + tas * (α_dot * cos(α) * cos(β) - β_dot * sin(α) * sin(β))
    return [u_dot, v_dot, w_dot]
end


text:
using MJPlayGround
using Base.Test

# write your own tests here
@test 1 == 1

a=2;
a


text:
import Base: show, copy

"""
$(TYPEDEF)
    options for Solver
"""
mutable struct SolverOptions

    ##########################
    ## Satellite Quaternion ##
    ##########################

    sat_att::Bool #if this is true, the 4:7 numbers in the state vector are a quaternion

    ###################
    ## Functionality ##
    ###################
    verbose::Bool # Display solve statistics at each iteration
    benchmark::Bool # Run benchmarks on forward and backward passes
    live_plotting::Bool # Plot state and control trajectories during solve

    ##########################
    ## Convergence Criteria ##
    ##########################
    cost_tolerance::Float64 # dJ < ϵ, cost convergence criteria for unconstrained solve or to enter outerloop for constrained solve
    cost_tolerance_intermediate::Float64 # dJ < ϵ_int, intermediate cost convergence criteria to enter outerloop of constrained solve
    constraint_tolerance::Float64 # max(constraint) < ϵ, constraint convergence criteria
    constraint_tolerance_intermediate::Float64 # max(constraint) < ϵ_int, intermediate constraint convergence criteria
    constraint_tolerance_infeasible::Float64  # infeasible control constraint tolerance
    gradient_norm_tolerance::Float64 # gradient_norm < ϵ, gradient norm convergence criteria
    gradient_norm_tolerance_intermediate::Float64 # gradient_norm_int < ϵ, gradient norm intermediate convergence criteria
    gradient_type::Symbol # type of gradient to evaluate: Augmented Lagrangian ∂L/∂u :AuLa, feedfoward Σ|d|^2 :feedforward, Todorov normalized feedforward :todorov

    ######################
    ## Solve Parameters ##
    ######################
    iterations::Int64 # caps cumulative number of iLQR iterations (see iterations_innerloop, iterations_outerloop)
    dJ_counter_limit::Int64 # restricts the total number of times a forward pass fails, resulting in regularization, before entering an outerloop update

    ###################
    ## Iterative LQR ##
    ###################
    square_root::Bool # use square root method backward pass for numerical conditioning
    iterations_innerloop::Int64 # maximum iterations for inner loop (iLQR)
    line_search_lower_bound::Float64 # forward pass approximate line search lower bound, 0 < line_search_lower_bound < line_search_upper_bound
    line_search_upper_bound::Float64 # forward pass approximate line search upper bound, 0 < line_search_lower_bound < line_search_upper_bound < Inf
    iterations_linesearch::Int64 # maximum number of backtracking steps during forward pass line search

    # Regularization
    bp_reg_initial::Float64 # initial regularization
    bp_reg_increase_factor::Float64 # regularization scaling factor
    bp_reg_max::Float64 # maximum regularization value
    bp_reg_min::Float64 # minimum regularization value
    bp_reg_type::Symbol # type of regularization- control: () + ρI, state: (S + ρI); see "Synthesis and Stabilization of Complex Behaviors through Online Trajectory Optimization"
    bp_reg_fp::Float64 # additive regularization when forward pass reaches max iterations
    eigenvalue_scaling::Float64 # add this multiple of the magnitude of the most negative eigenvalue to Quu decomposition to make positive definite
    eigenvalue_threshold::Float64 # eigenvalues less than this threshold will be increased using the additive eigenvalue scaling
    bp_sqrt_inv_type::Symbol # type of fix for potential inverse failure, :pseudo -pseudo inverse of Wxx, :reg - regularize Wxx to make invertible
    bp_reg_sqrt_initial::Float64 # initial regularization for square root method
    bp_reg_sqrt_increase_factor::Float64 # regularization scaling factor for square root method

    ##########################
    ## Augmented Lagrangian ##
    ##########################
    iterations_outerloop::Int64 # maximum outerloop updates
    dual_min::Float64 # minimum Lagrange multiplier
    dual_max::Float64 # maximum Lagrange multiplier
    penalty_max::Float64 # maximum penalty term
    penalty_initial::Float64 # initial penalty term
    penalty_scaling::Float64 # penalty update multiplier; penalty_scaling > 0
    penalty_scaling_no::Float64 # penalty update multiplier when μ should not be update, typically 1.0 (or 1.0 + ϵ)
    constraint_decrease_ratio::Float64 # update term; 0 < constraint_decrease_ratio < 1
    outer_loop_update_type::Symbol # type of outer loop update (default, momentum, individual, accelerated)
    use_second_order_dual_update::Bool # second order update for Lagrange multipliers once sqrt(cost tolerance | gradient) < desired tolerance
    penalty_update_frequency::Int  # determines how many iterations should pass before the penalty is updated (1 is every iteration)
    constraint_tolerance_second_order_dual_update::Float64 # constraint tolerance for switching to second order dual update
    active_constraint_tolerance::Float64 # numerical tolerance for constraint violation
    use_nesterov::Bool # accelerated gradient descent for dual variables
    use_penalty_burnin::Bool # perform only penalty updates (no dual updates) until constraint_tolerance_intermediate < ϵ_int
    al_type::Symbol  # [:default, :algencan] pick which Augmented Lagrangian formulation to use. ALGENCAN uses a slightly different check for inactive constraints.

    ######################
    ## Infeasible Start ##
    ######################
    R_infeasible::Float64 # regularization term for infeasible controls
    resolve_feasible::Bool # resolve feasible problem after infeasible solve
    feasible_projection::Bool # project infeasible solution into feasible space w/ BP, rollout
    penalty_initial_infeasible::Float64 # initial penalty term for infeasible controls
    penalty_scaling_infeasible::Float64 # penalty update rate for infeasible controls

    ##################
    ## Minimum Time ##
    ##################
    R_minimum_time::Float64 # regularization term for dt
    max_dt::Float64 # maximum allowable dt
    min_dt::Float64 # minimum allowable dt
    minimum_time_tf_estimate::Float64 # initial guess for the length of the minimum time problem (in seconds)
    minimum_time_dt_estimate::Float64 # initial guess for dt of the minimum time problem (in seconds)
    penalty_initial_minimum_time_inequality::Float64 # initial penalty term for minimum time bounds constraints
    penalty_initial_minimum_time_equality::Float64 # initial penalty term for minimum time equality constraints
    penalty_scaling_minimum_time_inequality::Float64 # penalty update rate for minimum time bounds constraints
    penalty_scaling_minimum_time_equality::Float64 # penalty update rate for minimum time equality constraints

    #############################
    ## Solver Numerical Limits ##
    #############################
    max_cost_value::Float64 # maximum cost value, if exceded solve will error
    max_state_value::Float64 # maximum state value, evaluated during rollout, if exceded solve will error
    max_control_value::Float64 # maximum control value, evaluated during rollout, if exceded solve will error

    function SolverOptions(;sat_att=false,verbose=false,benchmark=false,
        live_plotting=false,
        cost_tolerance=1.0e-4,
        cost_tolerance_intermediate=1.0e-3,
        constraint_tolerance=1.0e-3,
        constraint_tolerance_intermediate=sqrt(constraint_tolerance),
        constraint_tolerance_infeasible=1.0e-5,
        gradient_norm_tolerance=1.0e-5,
        gradient_norm_tolerance_intermediate=1.0e-5,
        gradient_type=:todorov,
        iterations=500,
        dJ_counter_limit=10,
        square_root=false,
        iterations_innerloop=250,
        line_search_lower_bound=1.0e-8,
        line_search_upper_bound=10.0,
        iterations_linesearch=20,
        bp_reg_initial=0.0,
        bp_reg_increase_factor=1.6,
        bp_reg_max=1.0e8,
        bp_reg_min=1.0e-8,
        bp_reg_type=:control,
        bp_reg_fp=10.0,
        eigenvalue_scaling=2.0,
        eigenvalue_threshold=1e-8,
        bp_sqrt_inv_type=:pseudo,
        bp_reg_sqrt_initial=1.0e-6,
        bp_reg_sqrt_increase_factor=10.0,
        iterations_outerloop=30,
        dual_min=-1.0e8,
        dual_max=1.0e8,
        penalty_max=1.0e8,
        penalty_initial=1.0,
        penalty_scaling=10.0,
        penalty_scaling_no=1.0,
        constraint_decrease_ratio=0.25,
        outer_loop_update_type=:default,
        use_second_order_dual_update=false,
        penalty_update_frequency=1,
        constraint_tolerance_second_order_dual_update=sqrt(constraint_tolerance),
        active_constraint_tolerance=0.0,
        use_nesterov=false,
        use_penalty_burnin=false,
        al_type=:default,
        R_infeasible=1.0,
        resolve_feasible=true,
        feasible_projection=true,
        penalty_initial_infeasible=1.0,
        penalty_scaling_infeasible=10.0,
        R_minimum_time=1.0,
        max_dt=1.0,
        min_dt=1.0e-3,
        minimum_time_tf_estimate=0.0,
        minimum_time_dt_estimate=0.0,
        penalty_initial_minimum_time_inequality=1.0,
        penalty_initial_minimum_time_equality=1.0,
        penalty_scaling_minimum_time_inequality=1.0,
        penalty_scaling_minimum_time_equality=1.0,
        max_cost_value=1.0e8,
        max_state_value=1.0e8,
        max_control_value=1.0e8)

        new(sat_att,
            verbose,
            benchmark,
            live_plotting,
            cost_tolerance,
            cost_tolerance_intermediate,
            constraint_tolerance,
            constraint_tolerance_intermediate,
            constraint_tolerance_infeasible,
            gradient_norm_tolerance,
            gradient_norm_tolerance_intermediate,
            gradient_type,
            iterations,
            dJ_counter_limit,
            square_root,
            iterations_innerloop,
            line_search_lower_bound,
            line_search_upper_bound,
            iterations_linesearch,
            bp_reg_initial,
            bp_reg_increase_factor,
            bp_reg_max,
            bp_reg_min,
            bp_reg_type,
            bp_reg_fp,
            eigenvalue_scaling,
            eigenvalue_threshold,
            bp_sqrt_inv_type,
            bp_reg_sqrt_initial,
            bp_reg_sqrt_increase_factor,
            iterations_outerloop,
            dual_min,
            dual_max,
            penalty_max,
            penalty_initial,
            penalty_scaling,
            penalty_scaling_no,
            constraint_decrease_ratio,
            outer_loop_update_type,
            use_second_order_dual_update,
            penalty_update_frequency,
            constraint_tolerance_second_order_dual_update,
            active_constraint_tolerance,
            use_nesterov,
            use_penalty_burnin,
            al_type,
            R_infeasible,
            resolve_feasible,
            feasible_projection,
            penalty_initial_infeasible,
            penalty_scaling_infeasible,
            R_minimum_time,
            max_dt,
            min_dt,
            minimum_time_tf_estimate,
            minimum_time_dt_estimate,
            penalty_initial_minimum_time_inequality,
            penalty_initial_minimum_time_equality,
            penalty_scaling_minimum_time_inequality,
            penalty_scaling_minimum_time_equality,
            max_cost_value,
            max_state_value,
            max_control_value)
    end

end

copy(opts::SolverOptions) = SolverOptions(;[name=>getfield(opts,name) for name in fieldnames(typeof(opts))]...)

function Base.:(==)(A::SolverOptions, B::SolverOptions)
    for name in fieldnames(A)
        if getfield(A,name) != getfield(B,name)
            return false
        end
    end
    return true
end


text:
module LawTests

using Test
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData

#e = @macroexpand @law function  l(a)
#  println("a is Any")
#end
#
#display(e)

@law ν(u,x) = (u+1.0)*x[1]
@law dν(du,x) = x[1]*du

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

order = 2
diritag = "boundary"
V = GradConformingFESpace(Float64,model,order,diritag)
U = TrialFESpace(V)

trian = get_triangulation(model)
degree = 1
quad = CellQuadrature(trian,degree)
q = get_coordinates(quad)

x = get_physical_coordinate(trian)
u = zero(U)
du = get_cell_basis(U)

r = ν(u,x)
@test isa(r,CellField)
@test evaluate(r,q) == evaluate(operate(ν,u,x),q)

dr = dν(du,x)
@test isa(dr,CellField)
@test is_basis(dr)
@test is_trial(dr)
@test evaluate(dr,q) == evaluate(operate(dν,du,x),q)

end # module


text:
# Reader
# ======

"""
    BEDReader(input::IO)

Create a data reader of the BED file format.

# Arguments
* `input`: data source
"""
type BEDReader <: Bio.IO.AbstractReader
    state::Ragel.State

    # intermediate values used during parsing
    red::Float32
    green::Float32
    blue::Float32
    block_size_idx::Int
    block_first_idx::Int

    function BEDReader(input::BufferedInputStream)
        return new(Ragel.State(bedparser_start, input), 0, 0, 0, 1, 1)
    end
end

function Bio.IO.stream(reader::BEDReader)
    return reader.state.stream
end

function Intervals.metadatatype(::BEDReader)
    return BEDMetadata
end

function Base.eltype(::Type{BEDReader})
    return BEDInterval
end

function BEDReader(input::IO)
    return BEDReader(BufferedInputStream(input))
end

function IntervalCollection(interval_stream::BEDReader)
    intervals = collect(BEDInterval, interval_stream)
    return IntervalCollection{BEDMetadata}(intervals, true)
end


text:
function gbset(field, value)
    libgb.GxB_Global_Option_set(field, value)
    return nothing
end

function gbget(field)
    return libgb.GxB_Global_Option_get(field)
end

function gbset(A::GBMatrix, field, value)
    libgb.GxB_Matrix_Option_set(A, field, value)
    return nothing
end

function gbget(A::GBMatrix, field)
    return libgb.GxB_Matrix_Option_get(A, field)
end

function gbset(A::GBVector, field, value)
    libgb.GxB_Matrix_Option_set(A, field, value)
    return nothing
end

function gbget(A::GBVector, field)
    return libgb.GxB_Matrix_Option_get(A, field)
end

function format(A::GBVecOrMat)
    t = gbget(A, SPARSITY_STATUS)
    f = gbget(A, FORMAT)
    return (GBSparsity(t), GBFormat(f))
end

const HYPER_SWITCH = libgb.GxB_HYPER_SWITCH
const BITMAP_SWITCH = libgb.GxB_BITMAP_SWITCH
const FORMAT = libgb.GxB_FORMAT
const SPARSITY_STATUS = libgb.GxB_SPARSITY_STATUS
const SPARSITY_CONTROL = libgb.GxB_SPARSITY_CONTROL
const BASE1 = libgb.GxB_PRINT_1BASED
const NTHREADS = libgb.GxB_GLOBAL_NTHREADS
const BURBLE = libgb.GxB_BURBLE

const BYROW = libgb.GxB_BY_ROW
const BYCOL = libgb.GxB_BY_COL

"""
Sparsity options for GraphBLAS. values can be summed to produce additional options.
"""
@cenum GBSparsity::Int32 begin
    GBDENSE = 8 #libgb.GxB_FULL
    GBBITMAP = 4 #libgb.GxB_BITMAP
    GBSPARSE = 2 #libgb.GxB_SPARSE
    GBHYPER = 1 #libgb.GxB_HYPERSPARSE
    GBANYSPARSITY = 15 #libgb.GxB_ANY_SPARSITY
    GBDENSE_OR_BITMAP = 12 #libgb.GxB_FULL + libgb.GXB_BITMAP
    GBSPARSE_OR_HYPER = 3 #libgb.GxB_SPARSE + libgb.GXB_HYPERSPARSE
    #... Probably don't need others
end