text:
module scalogram
    using PyCall

"""
Function that uses Pythons' PyWavelets package to compute the scalogram (or spectrogram) of a phase measurement.

Inputs:
        phi::Array{Float64}   -- the vector of measured phases
        sweepRange::Float64   -- the frequency sweep range of the reflectometer (GHz)
        sweepTime::Float64    -- the frequency sweep time of the reflectometer (seconds)
        df::Float64           -- the Δf between each frequency measurement (GHz)

Optional inputs:
        cmor::Int64           -- parameter controlling the width of the wavelet. Default is 5.

Output:
        Wxx   -- the 2-D array of spectrogram values
        freq2 -- the spectrogram frequency values (MHz)
"""
    function WxxCalc(phi::Array{Float64},sweepRange::Number,sweepTime::Float64, df::Float64; cmor::Int64=5)
        fs=1/(1e6/((sweepRange)/sweepTime))
        s0=50
        smax=400
        wave=string("cmor", string(cmor), "-1")
        pywt = pyimport("pywt")
        Wxx,freq2=pywt.cwt(phi,scales,wave,sampling_period=1)
        freq2 = freq2/(df*1e9)
        Wxx = abs.(Wxx)

        Wxx, freq2
    end

"""
Function that calculates the maximum spectrogram frequency from a spectrogram.

Input:
        Wxx::Array{Float64}        -- the 2-D array of spectrogram values from WxxCalc
        freq                       -- vector of reflectometer frequency meausurements (GHz)
        datafreq::Array{Float64}  -- the spectrogram frequency values from WxxCalc

Optional Inputs:
         maxDelay::Number   -- The maximum spectrogram frequency(group delay) to be considered.
         minDelay::Number   -- The minimum spectrogram frequency (group delay) to be considered.
"""
    function max_calc(Wxx::Array{Float64},freq,datafreq::Array{Float64}; maxDelay::Number=1e10, minDelay::Number=0.0)
        local maxfreq = length(datafreq)
        local minfreq = 1
        for j in 1:1:length(datafreq)
            if datafreq[j] < maxDelay
                maxfreq=j
                break
            end
        end
        for j in 1:1:length(datafreq)
            if datafreq[j] < minDelay
                minfreq = j
                break
            end
        end
        maxvec = zeros(length(freq))
        for k in 1:1:length(freq)
            step = findall(x->x==maximum(Wxx[maxfreq:minfreq,k]),Wxx[maxfreq:minfreq,k])[1]
            maxvec[k] = datafreq[maxfreq+step]
        end
        maxvec
    end
end


text:
#Updated implementation inspired by https://groups.google.com/forum/#!topic/julia-users/UPNpP91QGbE

immutable ExampleFlags
    flag::Int64;
    function ExampleFlags(f::Int64)
        @assert 1<=f && f<=2
        return new(f);
    end
end

const flag1 = ExampleFlags(1);
const flag2 = ExampleFlags(2);

function foo(input::FlaggedType)
    if input.flag == flag1
        print("Flag is flag1\n");
    elseif infput.flag == flag2
        print("Flag is flag2\n");
    end
end

ext: jl

max_stars_repo_name: ReidAtcheson/EnumsInJulia

max_issues_repo_name: ReidAtcheson/EnumsInJulia

max_forks_repo_name: ReidAtcheson/EnumsInJulia


text:
# Script simulating the CIR process with different parameters,
# and writing the result to cir_examples.jld

using CIRProc

include("utils.jl")

# general settings
datafile = "cir_examples.jld"

#maxtau_prctile = 0.99

# process base parameters
mu_tau = 0.4
sig_tau = 0.2
theta_tau = 2
dt = 0.005
tmax = 10
n = int(div(tmax, dt))
ts = [(0:(n-1)) * dt]

## different mu_tau
println("-- Simulations for different mu_tau")
mu_taus = [sqrt(0.2^2/2), 0.40, 1.00]
nmu_tau = length(mu_taus)
tau_traj = fill(NaN, nmu_tau, n)
for i in 1:nmu_tau
    p = CIRProcess(mu_taus[i], sig_tau, theta_tau)
    tau_traj[i,:] = CIRProc.simcir(ExactCIRSim(p, dt), n)
end
println("-- Writing data to $datafile/var_mu_tau")
@writedata(datafile, "var_mu_tau",
    mu_taus, sig_tau, theta_tau, ts, tau_traj)
println()

## different sig_tau
println("-- Simulations for different sig_tau")
sig_taus = [0.05, 0.1, 0.2, sqrt(2*mu_tau^2)] # needs to be smaller sqrt(2 mu_tau^2)
nsig_tau = length(sig_taus)
tau_traj = fill(NaN, nsig_tau, n)
for i in 1:nsig_tau
    p = CIRProcess(mu_tau, sig_taus[i], theta_tau)
    tau_traj[i,:] = CIRProc.simcir(ExactCIRSim(p, dt), n)
end
println("-- Writing data to $datafile/var_mu_tau")
@writedata(datafile, "var_sig_tau",
    mu_tau, sig_taus, theta_tau, ts, tau_traj)
println()

## theta_tau
println("-- Simulations for different theta_tau")
theta_taus = [0.25, 1.0, 4.0, 16.0]
ntheta_tau = length(theta_taus)
tau_traj = fill(NaN, ntheta_tau, n)
for i in 1:ntheta_tau
    p = CIRProcess(mu_tau, sig_tau, theta_taus[i])
    tau_traj[i,:] = CIRProc.simcir(ExactCIRSim(p, dt), n)
end
println("-- Writing data to $datafile/var_theta_tau")
@writedata(datafile, "var_theta_tau",
    mu_tau, sig_tau, theta_taus, ts, tau_traj)
println()


text:
using Metatheory

@testset "Type Assertions in Ematcher" begin
    some_theory = @theory begin
        a * b => b * a
        a::Number * b::Number => matched(a,b)
        a::Int64 * b::Int64 => specific(a,b)
        a * (b * c) => (a * b) * c
    end

    G = EGraph(:(2*3))

    @test true == areequal(G, some_theory, :(2 * 3), :(matched(2,3)))
    @test true == areequal(G, some_theory, :(matched(2,3)), :(specific(3,2)))
end


text:
"""
    struct QCurlGradMonomialBasis{...} <: AbstractArray{Monomial}

This type implements a multivariate vector-valued polynomial basis
spanning the space needed for Raviart-Thomas reference elements on n-cubes.
The type parameters and fields of this `struct` are not public.
This type fully implements the [`Field`](@ref) interface, with up to first order
derivatives.
"""
struct QCurlGradMonomialBasis{D,T} <: AbstractVector{Monomial}
  qgrad::QGradMonomialBasis{D,T}
  function QCurlGradMonomialBasis(::Type{T},order::Int,terms::CartesianIndices{D},perms::Matrix{Int}) where {D,T}
    qgrad = QGradMonomialBasis(T,order,terms,perms)
    new{D,T}(qgrad)
  end
end

Base.size(a::QCurlGradMonomialBasis) = (length(a.qgrad),)
# @santiagobadia : Not sure we want to create the monomial machinery
Base.getindex(a::QCurlGradMonomialBasis,i::Integer) = Monomial()
Base.IndexStyle(::QCurlGradMonomialBasis) = IndexLinear()

"""
    QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}

Returns a `QCurlGradMonomialBasis` object. `D` is the dimension
of the coordinate space and `T` is the type of the components in the vector-value.
The `order` argument has the following meaning: the divergence of the  functions in this basis
is in the Q space of degree `order`.
"""
function QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
  @assert T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
  _order = order+1
  _t = tfill(_order,Val{D-1}())
  t = (_order+1,_t...)
  terms = CartesianIndices(t)
  perms = _prepare_perms(D)
  QCurlGradMonomialBasis(T,order,terms,perms)
end

return_type(::QCurlGradMonomialBasis{D,T}) where {D,T} = T

function return_cache(f::QCurlGradMonomialBasis,x::AbstractVector{<:Point})
  return_cache(f.qgrad,x)
end

function evaluate!(cache,f::QCurlGradMonomialBasis,x::AbstractVector{<:Point})
  evaluate!(cache,f.qgrad,x)
end

function return_cache(
  fg::FieldGradientArray{N,<:QCurlGradMonomialBasis},
  x::AbstractVector{<:Point}) where N

  f = fg.fa
  return_cache(FieldGradientArray{N}(f.qgrad),x)
end

function evaluate!(
  cache,
  fg::FieldGradientArray{N,<:QCurlGradMonomialBasis},
  x::AbstractVector{<:Point}) where N

  f = fg.fa
  evaluate!(cache,FieldGradientArray{N}(f.qgrad),x)
end

"""
    num_terms(f::QCurlGradMonomialBasis{D,T}) where {D,T}
"""
num_terms(f::QCurlGradMonomialBasis{D,T}) where {D,T} = length(f.qgrad.terms)*D

get_order(f::QCurlGradMonomialBasis{D,T}) where {D,T} = get_order(f.qgrad)


text:
module test_octo_structured

using Test
using Octo.Adapters.SQL # from to_sql

# IN
@test to_sql([IN ("Germany", "France", "UK")]) == "IN ('Germany', 'France', 'UK')"

struct Customer
end
Schema.model(Customer, table_name="Customers")

customers = from(Customer)

@test to_sql([FROM Customer]) == "FROM Customers"

# IN
@test to_sql([customers.Country IN ("Germany", "France", "UK")]) == "Country IN ('Germany', 'France', 'UK')"
@test to_sql([customers.Country NOT IN ("Germany", "France", "UK")]) == "Country NOT IN ('Germany', 'France', 'UK')"

@test to_sql([WHERE true]) == "WHERE true"
@test to_sql([1,2,3]) == "1 2 3"
@test to_sql([(1,2,3)]) == "1, 2, 3"

struct Department
end
Schema.model(Department, table_name="departments")
struct Employee
end
Schema.model(Employee, table_name="employees")
d = from(Department)
em = from(Employee)
sub = from([SELECT d.department_id FROM d WHERE d.location_id == 1800])
q = [SELECT (em.first_name, em.last_name, em.department_id) FROM em WHERE em.department_id IN sub]
@test to_sql(q) == "SELECT first_name, last_name, department_id FROM employees WHERE department_id IN (SELECT department_id FROM departments WHERE location_id = 1800)"

# https://dba.stackexchange.com/questions/21749/count-with-set-difference-and-union
struct Tbl1
end
Schema.model(Tbl1, table_name="tbl1")

struct Tbl2
end
Schema.model(Tbl2, table_name="tbl2")

struct Tbl3
end
Schema.model(Tbl3, table_name="tbl3")

tbl1 = from(Tbl1, :tbl1)
tbl2 = from(Tbl2, :tbl2)
tbl3 = from(Tbl3, :tbl3)
subexcept = from([SELECT tbl2.id FROM tbl2 UNION SELECT tbl3.id FROM tbl3])
sub = from([SELECT tbl1.id FROM tbl1 EXCEPT subexcept], :t)
@test to_sql([SELECT COUNT(*) FROM sub]) == "SELECT COUNT(*) FROM (SELECT tbl1.id FROM tbl1 EXCEPT (SELECT tbl2.id FROM tbl2 UNION SELECT tbl3.id FROM tbl3)) AS t"

sub = from([SELECT tbl1.id FROM tbl1 EXCEPT SELECT tbl2.id FROM tbl2 EXCEPT SELECT tbl3.id FROM tbl3], :t)
@test to_sql([SELECT COUNT(*) FROM sub]) == "SELECT COUNT(*) FROM (SELECT tbl1.id FROM tbl1 EXCEPT SELECT tbl2.id FROM tbl2 EXCEPT SELECT tbl3.id FROM tbl3) AS t"

# Issue #13
struct User
end
Schema.model(User, table_name="users")
users = from(User)
@test to_sql([INSERT INTO users (users.name, users.email) VALUES ("Jick", "Jick@dd.com")]) == "INSERT INTO users (name, email) VALUES ('Jick', 'Jick@dd.com')"
users = from(User, :users)
@test to_sql([INSERT INTO users (users.name, users.email) VALUES ("Jick", "Jick@dd.com")]) == "INSERT INTO users (name, email) VALUES ('Jick', 'Jick@dd.com')"
@test to_sql([INSERT INTO users (:name, :email) VALUES ("Jick", "Jick@dd.com")]) == "INSERT INTO users (name, email) VALUES ('Jick', 'Jick@dd.com')"

end # module test_octo_structured


text:
module SixDOF

using LinearAlgebra: norm, cross
using Plots

export MassProp, State
export Atmosphere
#export sixdof!
# ------------ GENERIC MODULES -------------------------------------------------
import Dierckx
import CSV
import DataFrames
import JLD
import Dates
using PyPlot
using LinearAlgebra: cross, I

# ------------ FLOW CODES ------------------------------------------------------
# FLOWVLM https://github.com/byuflowlab/FLOWVLM
import FLOWVLM
const vlm = FLOWVLM

# FLOWVPM https://github.com/byuflowlab/FLOWVPM.jl
try                     # Load FLOWVPM if available
    import FLOWVPM
catch e                 # Otherwise load a dummy version of FLOWVPM
    @warn("FLOWVPM module not found. Using dummy module instead.")
    include("FLOWUnsteady_dummy_FLOWVPM.jl")
end
const vpm = FLOWVPM

# GeometricTools https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
const gt = GeometricTools

import FLOWNoise
const noise = FLOWNoise

# BPM https://github.com/byuflowlab/BPM.jl
import BPM

# ------------ GLOBAL VARIABLES ------------------------------------------------
const module_path = splitdir(@__FILE__)[1]                # Path to this module
const def_data_path = joinpath(module_path, "../data/")   # Default path to data folder

# ------------ HEADERS ---------------------------------------------------------
# Load modules
#for module_name in ["vehicle", "vehicle_vlm",
#                    "maneuver", "rotor",
#                    "simulation_types", "simulation", "utils",
#                    "processing", "monitors",
#                    "noise_wopwop", "noise_bpm"]
                    include("FLOWUnsteady_vehicle.jl")
                    include("FLOWUnsteady_vehicle_vlm.jl")
                    include("FLOWUnsteady_maneuver.jl")
                    include("FLOWUnsteady_rotor.jl")
                    include("FLOWUnsteady_simulation_types.jl")
                    include("FLOWUnsteady_simulation.jl")
                    include("FLOWUnsteady_utils.jl")
                    include("FLOWUnsteady_processing.jl")
                    include("FLOWUnsteady_monitors.jl")
                    include("FLOWUnsteady_noise_wopwop.jl")
                    include("FLOWUnsteady_noise_bpm.jl")
#end

#end # END OF MODULE

# ------ General Structs -------

"""
    State(x, y, z, phi, theta, psi, u, v, w, p, q, r)

State of the aircraft: positions in inertial frame, euler angles,
velocities in body frame, angular velocities in body frame.
"""
struct State{TF}
    x::TF  # position (inertial frame)
    y::TF
    z::TF
    phi::TF  # orientation, euler angles
    theta::TF
    psi::TF
    u::TF  # velocity (body frame)
    v::TF
    w::TF
    p::TF  # angular velocity (body frame)
    q::TF
    r::TF
end

"""
    Control(de, dr, da, df, throttle)

Define the control settings: delta elevator, delta rudder, delta aileron,
delta flaps, and throttle.
"""
#struct Control{TF}
#    de::TF  # elevator
#    dr::TF  # rudder
#    da::TF  # aileron
#    df::TF  # rudder
#    throttle::TF
#end

"""
    MassProp(m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz)

Mass and moments of inertia in the body frame.
Ixx = int(y^2 + z^2, dm)
Ixz = int(xz, dm)

We can model our quadcopter as two thin uniform rods crossed at the origin with a point mass
(motor) at the end of each. With this in mind, it’s clear that the symmetries result in a diagonal
inertia matrix
"""
struct MassProp{TF}
    m::TF
    Ixx::TF
    Iyy::TF
    Izz::TF
end

"""
    Reference(S, b, c)

The reference area equal for each rotor, S is the reference area (propeller cross-section, not
area swept out by the propeller)
"""
#struct Reference{TF}
    #S::TF  # area
    #b::TF  # span
    #c::TF  # chord
#end

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

# --------------- Interfaces ---------------

#abstract type AbstractAtmosphereModel end

"""
inizialization of the state position
"""
#function pos(state)
#     x = [state.x; state.y; state.z]
#    Wi = [0.0, 0.0, 0.0]
#    Wb = [0.0, 0.0, 0.0]
#    return x
#end

"""
    inizialization of the state velocities
"""
#function vel(state)
#     xdot = [state.u; state.v; state.w]
#    rho = 1.225  # sea-level properties
#    asound = 340.3
#    return xdot
#end

"""
    inizialization of the state angles
"""
#function angle(state)
#    theta = [state.phi; state.theta; state.psi]
#    g = 9.81
#    return theta
#end

# ----

#abstract type AbstractAeroModel end

"""

"""
#function anglevel(state)
#    omega = [state.p; state.q; state.r]
    # forces and moments in body frame
#    F = [0.0, 0.0, 0.0]
#    M = [0.0, 0.0, 0.0]
#    return omega
#end

# ----

#abstract type AbstractPropulsionModel end

"""
    Every variable used here is a matrix with the three axis for rows and time
    for columns
"""
function inizialization(N)
    row = 3
    matrix = zeros((row,N))
    return matrix
end
#    @warn "propulsionforces function not implemented for AbstractPropulsionModel"
    # forces and moments in body frame
#    F = [0.0, 0.0, 0.0]
#    M = [0.0, 0.0, 0.0]
#    return F, M
#end

# ----

#abstract type AbstractInertialModel end

"""
    gravityforces(model::AbstractInertialModel, atm::AbstractAtmosphereModel, state::State, control::Control, mp::MassProp, ref::Reference)

Compute the gravitational forces and moments in the body frame.
return F, M
"""
function equations(dt, i, mp, atm, k, b, kd, rotors, c, l, xprev, omegaprev, thetaprev, xdotprev)
    #the cycle for is initialized outside and then the function equations is called
        a = Array{Float64}(undef, 3)
        a = acceleration(i, thetaprev, xdotprev, mp, atm, kd, rotors, k)
        omegadot = Array{Float64}(undef, 3)
        #omegadot = omegadot(i, o, mp, l, k, b)
        I = [mp.Ixx   0.0   0.0;
             0.0     mp.Iyy 0.0;
             0.0      0.0  mp.Izz]
        tau = torque(i, l, b, k)
        #Inv = [1/I[1,1];1/I[2,2];1/I[3,3]]
        #prod = [tau[1]*Inv[1];tau[2]*Inv[2];tau[3]*Inv[3]]
        #diff = [((I[2,2]-I[3,3])/I[1,1])*o[2]*o[3];((I[3,3]-I[1,1])/I[2,2])*o[1]*o[3];
        #((I[1,1]-I[2,2])/I[3,3])*o[1]*o[2]]
        #omegadot = prod - diff
        omegadot = inv(I) * (tau - cross(omegaprev, I * omegaprev))
        omega = Array{Float64}(undef, 3)
        omega = omegaprev + dt * omegadot #i am not sure if they are
        #written well for the matrices
        thetadot = Array{Float64}(undef, 3)
        thetadot = omega2thetadot(thetaprev, omega)
        theta = Array{Float64}(undef, 3)
        theta = thetaprev + dt * thetadot
        xdot = Array{Float64}(undef, 3)
        xdot = xdotprev + dt * a
        x = Array{Float64}(undef, 3)
        x = xprev + dt * xdot
        return a ,omegadot, omega, thetadot, theta, xdot, x, tau
end

#    @warn "gravityforces function not implemented for AbstractInertialModel"
    # forces and moments in body frame
#    F = [0.0, 0.0, 0.0]
#    M = [0.0, 0.0, 0.0]
#    return F, M
#end

# ----

#abstract type AbstractController end

"""
    setcontrol(controller::AbstractController, time, atm::AbstractAtmosphereModel, state::State, lastcontrol::Control, mp::MassProp, ref::Reference)

Compute control state for next time step given current state
return control::Control
"""
#function setcontrol(controller::AbstractController, time, atm, state, mp, ref)
#    @warn "setcontrol function not implemented for AbstractController"
#    control = Control(0.0, 0.0, 0.0, 0.0, 0.0)
#    return control
#end

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

# ------------- helper functions (private) --------------

"""
    inertialtobody(state)

Construct a rotation matrix from inertial frame to body frame

The assumed order of rotation is
1) psi radians about the z axis,
2) theta radians about the y axis,
3) phi radians about the x axis.

This is an orthogonal transformation so its inverse is its transpose.
"""
function inertialtobody(theta)

    R = Array{eltype(theta)}(undef,3, 3)

    cphi, ctht, cpsi = cos.([theta[1], theta[2], theta[3]])
    sphi, stht, spsi = sin.([theta[1], theta[2], theta[3]])

    R[1, 1] = ctht*cpsi
    R[1, 2] = ctht*spsi
    R[1, 3] = -stht

    R[2, 1] = sphi*stht*cpsi - cphi*spsi
    R[2, 2] = sphi*stht*spsi + cphi*cpsi
    R[2, 3] = sphi*ctht

    R[3, 1] = cphi*stht*cpsi + sphi*spsi
    R[3, 2] = cphi*stht*spsi - sphi*cpsi
    R[3, 3] = cphi*ctht

    return R

end

function bodytoinertial(theta)

    Ri = Array{eltype(theta)}(undef,3, 3)

    cphi, ctht, cpsi = cos.([theta[1], theta[2], theta[3]])
    sphi, stht, spsi = sin.([theta[1], theta[2], theta[3]])

    Ri[1, 1] = cphi*cpsi - ctht*sphi*spsi
    Ri[1, 2] = -cpsi*sphi - cphi*ctht*spsi
    Ri[1, 3] = stht*spsi

    Ri[2, 1] = ctht*sphi*cpsi + cphi*spsi
    Ri[2, 2] = cphi*ctht*cpsi - sphi*spsi
    Ri[2, 3] = -cpsi*stht

    Ri[3, 1] = sphi*stht
    Ri[3, 2] = cphi*stht
    Ri[3, 3] = ctht

    return Ri

end

"""
time of simulation giving the the initial and final time and the time step
"""
function time(tinit, tfinal, dt);
    n = (tfinal-tinit)/dt
    n = trunc(Int,n)
    t = tinit
    v = Array{Float64}(undef, n)
    v[1] = t
    for i in 2: n
        v[i] = v[i-1] + dt
    end
return v,n

end

"""
 Total thrust, giving the square rotation velocity and the coefficient k
"""
function thrust(rpm, k, rotors)
sum=0
for i in 1:rotors
    sum  = sum + rpm[i]
end
    T = [0.0; 0.0; k*sum]
    return T
end

"""
Torque, giving the square rotation velocity, the distance of rotors from the center
and the coefficients k and b
"""
function torque(rpm, l, b, k)
    tau = [l*k*(rpm[1]-rpm[3]); l*k*(rpm[2]-rpm[4]); b*(rpm[1]-rpm[2]+rpm[3]-rpm[4])]
    return tau
end

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

# ------- Some Default Interface Implementations -----

"""
    StabilityDeriv(CL0, CLalpha, CLq, CLM, CLdf, CLde, alphas,
        CD0, U0, exp_Re, e, Mcc, CDdf, CDde, CDda, CDdr,
        CYbeta, CYp, CYr, CYda, CYdr, Clbeta,
        Clp, Clr, Clda, Cldr,
        Cm0, Cmalpha, Cmq, CmM, Cmdf, Cmde,
        Cnbeta, Cnp, Cnr, Cnda, Cndr)

Stability derivatives of the aircraft.  Most are self explanatory if you are
familiar with stability derivatives (e.g., CLalpha is dCL/dalpha or the lift curve slope).
Some less familiar ones include
- M: Mach number
- alphas: the angle of attack for stall
- U0: the speed for the reference Reynolds number CD0 was computed at
- exp_Re: the coefficient in the denominator of the skin friction coefficient (0.5 laminar, 0.2 turbulent)
- e: Oswald efficiency factor
- Mcc: crest critical Mach number (when compressibility drag rise starts)

"""
#struct StabilityDeriv{TF} <: AbstractAeroModel
#    CL0::TF
#    CLalpha::TF
#    CLq::TF
#    CLM::TF
#    CLdf::TF
#    CLde::TF
#    alphas::TF  # TODO: should probably do in terms of CLmax

#    CD0::TF
#    U0::TF  # velocity corresponding to Reynolds number of CD0  (TODO: rethink this)
#    exp_Re::TF  # exponent for Reynolds number scaling. typical values: exp_Re = 0.5 laminar, 0.2 turbulent
#    e::TF  # Oswald efficiency factor
#    Mcc::TF  # crest-critical Mach number when compressibility drag rise starts (quartic)
#    CDdf::TF
#    CDde::TF
#    CDda::TF
#    CDdr::TF

#    CYbeta::TF
#    CYp::TF
#    CYr::TF
#    CYda::TF
#    CYdr::TF

#    Clbeta::TF
#    Clp::TF
#    Clr::TF
#    Clda::TF
#    Cldr::TF

#    Cm0::TF
#    Cmalpha::TF
#    Cmq::TF
#    CmM::TF
#    Cmdf::TF
#    Cmde::TF

#    Cnbeta::TF
#    Cnp::TF
#    Cnr::TF
#    Cnda::TF
#    Cndr::TF
#end

"""
Computation of the linear acceleration
"""
function acceleration(rpm, theta, xdot, mp, atm, kd, rotors, k)
    gravity= [0.0; 0.0; -atm.g]
    R = bodytoinertial(theta)
    Tb = thrust(rpm, k, rotors)
    Ti = R * Tb
    Fd = -kd * xdot
    a = gravity + 1/mp.m  * Ti + 1/mp.m  * Fd
return a
end
    # Mach number and dynamic pressure
#    rho, asound = properties(atm, state)
#    Mach = Va / asound
#    qdyn = 0.5 * rho * Va^2

    # rename for convenience
#    p = state.p
#    q = state.q
#    r = state.r
#    de = control.de
#    df = control.df
#    dr = control.dr
#    da = control.da

    # lift
#    CL = sd.CL0 + sd.CLalpha*alpha + sd.CLq*q *ref.c/(2*Va) + sd.CLM*Mach
#        + sd.CLdf*df + sd.CLde*de

#    em = exp(-50*(alpha - sd.alphas))
#    ep = exp(50*(alpha + sd.alphas))
#    sigma = (1 + em + ep)/((1 + em)*(1 + ep))
#    CL = (1- sigma)*CL + sigma * 2 * sign(alpha)*sin(alpha)^2*cos(alpha)

    # drag
#    CDp = sd.CD0*(Va/sd.U0)^sd.exp_Re
#    CDi = CL^2/(pi*(ref.b^2/ref.S)*sd.e)
#    CDc = (Mach < sd.Mcc) ? 0.0 : 20*(Mach - sd.Mcc)^4

#    CD = CDp + CDi + CDc + abs(sd.CDdf*df) + abs(sd.CDde*de) + abs(sd.CDda*da) + abs(sd.CDdr*dr)

    # side force
#    CY = sd.CYbeta*beta + (sd.CYp*p + sd.CYr*r)*ref.b/(2*Va) + sd.CYda*da + sd.CYdr*dr

    # rolling moment
#    Cl = sd.Clbeta*beta + (sd.Clp*p + sd.Clr*r)*ref.b/(2*Va) + sd.Clda*da + sd.Cldr*dr

    # pitching moment
#    Cm = sd.Cm0 + sd.Cmalpha*alpha + sd.Cmq*q * ref.c/(2*Va) + sd.CmM*Mach + sd.Cmdf*df + sd.Cmde*de

    # yawing moment
#    Cn = sd.Cnbeta*beta + (sd.Cnp*p + sd.Cnr*r)*ref.b/(2*Va) + sd.Cnda*da + sd.Cndr*dr

    # transfer forces from wind to body axes
#    Rwb = windtobody(alpha, beta)

#    F = Rwb*[-CD, CY, -CL] * qdyn * ref.S

#    M = Rwb*[Cl*ref.b, Cm*ref.c, Cn*ref.b] * qdyn * ref.S

#    return F, M
#end

#@enum PropType CO=1 COUNTER=-1 COCOUNTER=0

"""
    MotorPropBatteryDataFit(CT2, CT1, CT0, CQ2, CQ1, CQ0, D, num, type,
        R, Kv, i0, voltage)

**Inputs**
- CT2, CT1, CT0: quadratic fit to propeller thrust coefficient of form: CT = CT2*J2 + CT1*J + CT0
- CQ2, CQ1, CQ0: quadratic fit to propeller torque coefficient of form: CQ = CQ2*J2 + CQ1*J + CQ0
- D: propeller diameter
- num: number of propellers
- type: CO (torques add), COUNTER (torques add but with minus sign), COCOUNTER (no torque, they cancel out)
- R: motor resistance
- Kv: motor Kv
- i0: motor no-load current
- voltage: battery voltage
"""
#struct MotorPropBatteryDataFit{TF, TI, PropType} <: AbstractPropulsionModel
    # CT = CT2*J2 + CT1*J + CT0
    # CQ = CQ2*J2 + CQ1*J + CQ0
#    CT2::TF  # prop data fit
#    CT1::TF
#    CT0::TF
#    CQ2::TF
#    CQ1::TF
#    CQ0::TF
#    D::TF  # prop diameter
#    num::TI
#    type::PropType

#    R::TF  # motor resistance
#    Kv::TF  # motor Kv
#    i0::TF  # motor no-load current

#    voltage::TF  # battery voltage
#end

function omegadot(rpm, omega, mp, l, k, b)

    I = [mp.Ixx   0.0   0.0;
         0.0     mp.Iyy 0.0;
         0.0      0.0  mp.Izz]
    tau = torque(rpm, l, b, k)
    #Inv = [1/I[1,1];1/I[2,2];1/I[3,3]]
    #prod = [tau[1]*Inv[1];tau[2]*Inv[2];tau[3]*Inv[3]]
    #diff = [((I[2,2]-I[3,3])/I[1,1])*omega[2]*omega[3];((I[3,3]-I[1,1])/I[2,2])*omega[1]*omega[3];
    #((I[1,1]-I[2,2])/I[3,3])*omega[1]*omega[2]]
    #omegadot = prod - diff
    omegadot = inv(I) .* (tau - cross(omega, I .* omega))
    #it seems there is a problem with cross function because they are float
return omegadot
end

    # density
#    rho, _ = properties(atm, state)

#    D = prop.D

    # determine torque for motor/prop match (quadratic equation)
#    a = rho*D^5/(2*pi)^2 * prop.CQ0
#    b = rho*D^4/(2*pi)*Va * prop.CQ1 + 1.0/(prop.R*prop.Kv)
#    c = rho*D^3*Va^2 * prop.CQ2 - control.throttle*prop.voltage/(prop.R*prop.Kv) + prop.i0/prop.Kv
#    Omega = (-b + sqrt(b^2 - 4*a*c))/(2*a)

    # advance ratio
#    n = Omega/(2*pi)
#    J = Va/(n*D)

    # thrust and torque
#    CT = prop.CT0 + prop.CT1*J + prop.CT2*J^2
#    CQ = prop.CQ0 + prop.CQ1*J + prop.CQ2*J^2

#    T = prop.num * CT * rho * n^2 * D^4
#    Q = prop.num * CQ * rho * n^2 * D^5 * Int(prop.type)

#    return [T, 0, 0], [Q, 0, 0]
#end

"""
    UniformGravitationalField()

Assumes center of mass and center of gravity are coincident.
"""
#struct UniformGravitationalField <: AbstractInertialModel end

function omega2thetadot(theta,omega)
    ct, cp = cos.([theta[2], theta[1]])
    st, sp = sin.([theta[2], theta[1]])

    C = [1.0  0.0  -st;
        0.0   cp  ct*sp
        0.0  -sp  ct*cp]

    thetadot = inv(C)*omega
    return thetadot
end

#    W = mp.m * gravity(atm, state)
#    ct, cp = cos.([state.theta, state.phi])
#    st, sp = sin.([state.theta, state.phi])

#    Fg = W*[-st, ct*sp, ct*cp]
#    Mg = [zero(W), zero(W), zero(W)]  # no gravitational moment

#    return Fg, Mg
#end

"""
    ConstantAtmosphere(Wi, Wb, rho, asound, g)

Constant atmospheric properties.
"""
struct Atmosphere{TF,TV}
    Wi::TV
    Wb::TV
    rho::TF
    asound::TF
    g::TF
end

"""
    ConstantController(de, dr, da, df, throttle)

Just a dummy controller that outputs constant control outputs the whole time.
"""
function plotsx(n, x, t)
    inspectdr()
    x1 = zeros(n)
    y1 = zeros(n)
    z1= zeros(n)
    for c in 1:n
    x1[c] = x[1,c]
    y1[c] = x[2,c]
    z1[c] = x[3,c]
    end
    p = Plots.plot(t, x1, label = "x" ,lw=3)
    Plots.plot!(p, t, y1, label = "y" )
    Plots.plot!(p, t, z1, label = "z" )
    Plots.xlabel!("t (s) ")
end

function plotstheta(n, theta, t)
    inspectdr()
    theta1 = zeros(n)
    theta2 = zeros(n)
    theta3= zeros(n)
    for c in 1:n
    theta1[c] = theta[1,c]
    theta2[c] = theta[2,c]
    theta3[c] = theta[3,c]
    end
    p = Plots.plot(t, theta1, label = "roll(phi)" ,lw=3)
    Plots.plot!(p, t, theta2, label = "pitch(theta)" )
    Plots.plot!(p, t, theta3, label = "yaw(qsi)" )
    Plots.xlabel!("t (s) ")
end

function plotsomega(n, omega, t)
    inspectdr()
    omega1 = zeros(n)
    omega2 = zeros(n)
    omega3= zeros(n)
    for c in 1:n
    omega1[c] = omega[1,c]
    omega2[c] = omega[2,c]
    omega3[c] = omega[3,c]
    end
    p = Plots.plot(t, omega1, label = "ω1" ,lw=3)
    Plots.plot!(p, t, omega2, label = "ω2" )
    Plots.plot!(p, t, omega3, label = "ω3" )
    Plots.xlabel!("t (s) ")
end

#    de::TF
#    dr::TF
#    da::TF
#    df::TF
#    throttle::TF
#end

#function setcontrol(controller::ConstantController, time, atm, state, mp, ref)
#    return Control(controller.de, controller.dr, controller.da, controller.df, controller.throttle)
#end
# --------------------------------------------------------

# ------------- main functions (public) --------------

"""
    sixdof!(ds, s, params, time)

dynamic and kinematic ODEs.  Follows format used in DifferentialEquations package.
- s = x, y, z, phi, theta, psi, u, v, w, p, q, r (same order as State)
- params = control, massproperties, reference, aeromodel, propmodel, inertialmodel, atmmodel
"""
#function sixdof!(ds, s, params, time)

#    x, y, z, phi, theta, psi, u, v, w, p, q, r = s
#    mp, ref, aeromodel, propmodel, inertialmodel, atmmodel, controller = params

    # ---- controller -------
#    state = State(s...)
    # -----------------------

    # --------- forces and moments ---------
    # aerodynamics
#    Fa, Ma = aeroforces(aeromodel, atmmodel, state, control, ref, mp)

    # propulsion
#    Fp, Mp = propulsionforces(propmodel, atmmodel, state, control, ref, mp)

    # weight
#    Fg, Mg = gravityforces(inertialmodel, atmmodel, state, control, ref, mp)

    # total forces and moments
#    F = Fa + Fp + Fg
#    M = Ma + Mp + Mg

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

    # ----- derivative of state --------
#    Vb = [u, v, w]
#    omegab = [p, q, r]

    # linear kinematics
#    Rib = inertialtobody(state)
#    rdot = Rib' * Vb

    # angular kinematics
#    phidot = p + (q*sin(phi) + r*cos(phi))*tan(theta)
#    thetadot = q*cos(phi) - r*sin(phi)
#    psidot = (q*sin(phi) + r*cos(phi))/cos(theta)

    # linear dynamics
#    vdot = F/mp.m - cross(omegab, Vb)

    # angular dynamics
#    I = [mp.Ixx -mp.Ixy -mp.Ixz;
#         -mp.Iyz mp.Iyy -mp.Iyz;
#         -mp.Ixz -mp.Iyz mp.Izz]
#    omegadot = I \ (M - cross(omegab, I*omegab))

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

    # TODO: if we need more efficiency we can avoid allocating and then assigning.
#    ds[1:3] = rdot
#    ds[4] = phidot
#    ds[5] = thetadot
#    ds[6] = psidot
#    ds[7:9] = vdot
#    ds[10:12] = omegadot
#end

end # module


text:
function _monkey_delete(dir::AbstractString, frec::Float64; verb = true)
    frec < rand() || return
    for (root, _, files) in walkdir(dir)
        for file in files
            frec > rand() && continue
            rfile = joinpath(root, file)
            _gwrm(rfile)
            verb && println("monkey deleted!! ", rfile)
        end
    end
end


text:
function convert_comment(cursor::LibClang.CXCursor, jlname::String)
	cxcomment = LibClang.clang_Cursor_getParsedComment(cursor)
	kind = LibClang.clang_Comment_getKind(cxcomment)

	comment = kind == LibClang.CXComment_FullComment ? Markdown.MD(cxcomment) : nothing

	return jlname => Comment(comment, [CodeLocation(cursor)])
end

merge_comments!(comments::Dict{String, Comment}, more::Dict{String, Comment}) = foreach(((k, v),) -> merge_comment!(comments, k => v), collect(more))
function merge_comment!(comments::Dict{String, Comment}, comment::Pair{String, Comment})
	(sym, cmt) = comment

	isempty(sym) && return
	if haskey(comments, sym)
		md = if isnothing(comments[sym].md)
			isnothing(cmt.md) ? nothing : cmt.md
		else
			isnothing(cmt.md) ? comments[sym].md : Markdown.MD([comments[sym].md.content..., Markdown.HorizontalRule(), cmt.md.content...,])
		end
		comments[sym] = Comment(md, [comments[sym].locs..., cmt.locs...])
	else
		comments[sym] = cmt
	end
end

function Markdown.MD(cxcomment::LibClang.CXComment)
	hasDetails = false
	hasParams  = false
	hasReturns = false

	contents = []
	for ind in 1:LibClang.clang_Comment_getNumChildren(cxcomment)
		child = LibClang.clang_Comment_getChild(cxcomment, ind-1)
		kind = LibClang.clang_Comment_getKind(child)

		if !hasDetails && !isempty(contents)
			hasDetails = true
			push!(contents, Markdown.Header("Details", 2))
		end

		if kind == LibClang.CXComment_Paragraph
			Bool(LibClang.clang_Comment_isWhitespace(child)) || push!(contents, Markdown.Paragraph(child))
		elseif kind == LibClang.CXComment_BlockCommand
			para = LibClang.clang_BlockCommandComment_getParagraph(child)
			para = Markdown.Paragraph(para)

			num = LibClang.clang_BlockCommandComment_getNumArgs(child)
			cmd = _string(LibClang.clang_BlockCommandComment_getCommandName, child)
			if cmd == "brief"
				push!(contents, para)
			elseif cmd == "note" || cmd == "warning"
				para = Markdown.Paragraph(["$(uppercase(cmd)):", para.content...])
				push!(contents, para)
			elseif cmd == "sa" || cmd == "see"
				para = isempty(para.content) ? para : Markdown.Paragraph("See also: [`$(strip(first(para.content)))`](@ref)")
				push!(contents, para)
			elseif cmd == "returns" || cmd == "return"
				if !hasReturns
					hasReturns = true

					hdr = Markdown.Header("Returns", 2)
					if !isempty(contents) && contents[end] isa Markdown.Header
						contents[end] = hdr
					else
						push!(contents, hdr)
					end
				end

				push!(contents, para)
			elseif cmd == "li"
				if !isempty(contents) && contents[end] isa Markdown.List
					push!(contents[end].items, para)
				else
					push!(contents, Markdown.List(para))
				end
			else
				@warn "Unhandled block-command comment: $(cmd)"
			end
		elseif kind == LibClang.CXComment_ParamCommand
			if !hasParams
				hasParams = true

				hdr = Markdown.Header("Parameters", 2)
				if !isempty(contents) && contents[end] isa Markdown.Header
					contents[end] = hdr
				else
					push!(contents, hdr)
				end
			end

			_addParamCommand(contents, child)
		elseif kind == LibClang.CXComment_VerbatimBlockCommand
			push!(contents, Markdown.Code(child))
		elseif kind == LibClang.CXComment_VerbatimLine
			text = _string(LibClang.clang_VerbatimLineComment_getText, child)
			push!(contents, Markdown.Paragraph(text))
		else
			error("Unhandled document child: $(kind)")
		end
	end

	if length(contents) > 1 && contents[end] isa Markdown.Header
		contents = contents[1:end-1]
	end

	return Markdown.MD(contents)
end

function Markdown.Paragraph(cxcomment::LibClang.CXComment)
	contents = []

	for ind in 1:LibClang.clang_Comment_getNumChildren(cxcomment)
		child = LibClang.clang_Comment_getChild(cxcomment, ind-1)
		kind = LibClang.clang_Comment_getKind(child)

		if kind == LibClang.CXComment_Text
			text = _string(LibClang.clang_TextComment_getText, child)
			push!(contents, text)
		elseif kind == LibClang.CXComment_InlineCommand
			_addInlineCommand(contents, child)
		elseif kind == LibClang.CXComment_HTMLStartTag || kind == LibClang.CXComment_HTMLEndTag
			# TODO: handle HTML stuff...
		else
			error("Unhandled paragraph child: $(kind)")
		end
	end

	return Markdown.Paragraph(contents)
end

function Markdown.Code(cxcomment::LibClang.CXComment)
	lines = []

	for ind in 1:LibClang.clang_Comment_getNumChildren(cxcomment)
		child = LibClang.clang_Comment_getChild(cxcomment, ind-1)
		kind = LibClang.clang_Comment_getKind(child)

		if kind == LibClang.CXComment_VerbatimBlockLine
			line = _string(LibClang.clang_VerbatimBlockLineComment_getText, child)

			# a line with newlines is probably incorrectly parsed, so only keep before a newline
			push!(lines, first(split(line, '\n', limit=2)))
		else
			error("Unhandled code-block child: $(kind)")
		end
	end

	return Markdown.Code(join(lines, '\n'))
end

function _addInlineCommand(contents, cxcomment)
	num = LibClang.clang_InlineCommandComment_getNumArgs(cxcomment)
	cmd = _string(LibClang.clang_InlineCommandComment_getCommandName, cxcomment)

	text = _string(LibClang.clang_InlineCommandComment_getArgText, cxcomment, 0)
	kind = LibClang.clang_InlineCommandComment_getRenderKind(cxcomment)
	if kind == LibClang.CXCommentInlineCommandRenderKind_Bold
		push!(contents, Markdown.Bold(text))
	elseif kind == LibClang.CXCommentInlineCommandRenderKind_Emphasized
		push!(contents, Markdown.Italic(text))
	elseif kind == LibClang.CXCommentInlineCommandRenderKind_Monospaced
		push!(contents, Markdown.Code(text))  # TODO: maybe need Code(language, text)
	else
		push!(contents, text)
	end
end

function _addParamCommand(contents, cxcomment)
	num = LibClang.clang_Comment_getNumChildren(cxcomment)
	num == 1 || error("Incorrect number of parameter-command comment children")

	para = LibClang.clang_Comment_getChild(cxcomment, 0)
	LibClang.clang_Comment_getKind(para) == LibClang.CXComment_Paragraph || error("Expected a parameter-command comment paragraph")

	param = _string(LibClang.clang_ParamCommandComment_getParamName, cxcomment)
	param = Markdown.Paragraph([Markdown.Code(param), ":", Markdown.Paragraph(para).content...,])

	if !isempty(contents) && contents[end] isa Markdown.List
		push!(contents[end].items, param)
	else
		push!(contents, Markdown.List(param))
	end
end


text:
"""
    arbitraryparam(params)

returns the first field in the struct that is a subtype of `OpenSAFTParam`. errors if it finds none.
"""
function arbitraryparam(params)
    paramstype = typeof(params)
    idx = findfirst(z->z <: OpenSAFTParam,fieldtypes(paramstype))
    if isnothing(idx)
        error("The paramater struct ", paramstype, " must contain at least one OpenSAFTParam")
    end
     return fieldnames(paramstype)[idx] |> z->getfield(params,z)
end

"""
    @comps

This macro is an alias to

    model.icomponents

The caveat is that `model` has to exist in the local namespace.
`model` is expected to be an EoSModel type that contains the `icomponents` field.
`icomponents` is an iterator that goes through all component indices.
"""
macro comps()
    return :($(esc(:(model.icomponents))))
end

"""
    @groups

This macro is an alias to

    model.iflattenedgroups

`iflattenedgroups` is an iterator that goes through all groups in flattenedgroups.
"""
macro groups()
    return :($(esc(:(model.iflattenedgroups))))
end

"""
    @groups(component)

This macro is an alias to

    model.igroups[component]

`igroups[component]` is an iterator that goes through all groups in relevent to a given component.
"""
macro groups(component)
    return :($(esc(:(model.igroups[$(component)]))))
end

"""
    @sites(component)

This macro is an alias to

    model.isites[component]

`isites[component]` is an iterator that goes through all sites relevant to
each group in a GC model, and to each main component in a non-GC model.
"""
macro sites(component)
    return :($(esc(:(model.isites[$(component)]))))
end

"""
    @f(func,a,b,c,...)

This macro is an alias to

    func(model, V, T, z, a, b, c, ...)

where `func` is the name of the function, `model` is the model struct,
`V` is the volume, `T` is the absolute temperature, `z` is an array of number
of moles of each component, and `a`, `b`, `c`, ... are arbitrary parameters
that get passed to `func`.

It is very common for functions that are involved in the models to contain the
`model`, `V`, `T` and `z` parameters, so this macro helps reduce code repetition
as long as the first four parameters in the function are written exactly as above.

"""
macro f(func, args...)
    args = [esc(arg) for arg ∈ args]
    return :($(func)($(esc(:model)),$(esc(:V)),$(esc(:T)),$(esc(:z)),$(args...)))
end

"""
    @newmodelgc name parent paramstype

This is a data type that contains all the information needed to use an EoS model.
It also functions as an identifier to ensure that the right functions are called.

The user is expected to create an outter constructor that takes this signature

    function name(components::Array{String,1})

It should then return name(params::paramtype, groups::GroupParam, sites::SiteParam, idealmodel::IdealModel)

= Fields =
The Struct consists of the following fields:

* components: a string lists of components
* longthcomponents: the number of components
* icomponents: an iterator that goes through the indices corresponding to each component

* allcomponentgroups: a list of groups for each component
* lengthallcomponentgroups: a list containing the number of groups for each component
* allcomponentngroups: a list of the group multiplicity of each group corresponding to each group in allcomponentsgroup
* igroups: an iterable that contains a list of group indices corresponding to flattenedgroups for each component

* flattenedgroups: a list of all unique groups--the parameters correspond to this list
* lengthflattenedgroups: the number of unique groups
* allcomponentnflattenedgroups: the group multiplicities corresponding to each group in flattenedgroups
* iflattenedgroups: an iterator that goes through the indices for each flattenedgroup

* allgroupsites: a list containing a list of all sites corresponding to each group in flattenedgroups
* lengthallgroupsites: a list containing the number of unique sites for each group in flattenedgroups
* allgroupnsites: a list of the site multiplicities corresponding to each group in flattenedgroups
* isites: an iterator that goes through the indices corresponding to each group in flattenedgroups

* params: the Struct paramstype that contains all parameters in the model
* idealmodel: the IdealModel struct that determines which ideal model to use
* absolutetolerance: the absolute tolerance for solvers; the default value is 1E-12
* references: reference for this EoS

See the tutorial or browse the implementations to see how this is used.

"""

macro newmodelgc(name, parent, paramstype)
    quote
    struct $name{T <: IdealModel} <: $parent
        components::Array{String,1}
        lengthcomponents::Int
        icomponents::UnitRange{Int}

        allcomponentgroups::Array{Array{String,1},1}
        lengthallcomponentgroups::Array{Int,1}
        allcomponentngroups::Array{Array{Int,1},1}
        igroups::Array{Array{Int,1},1}

        flattenedgroups::Array{String,1}
        lengthflattenedgroups::Int
        allcomponentnflattenedgroups::Array{Array{Int,1},1}
        iflattenedgroups::UnitRange{Int}

        allgroupsites::Array{Array{String,1},1}
        lengthallgroupsites::Array{Int,1}
        allgroupnsites::Array{Array{Int,1},1}
        isites::Array{UnitRange{Int},1}

        params::$paramstype
        idealmodel::T
        absolutetolerance::Float64
        references::Array{String,1}
    end

    function $name(params::$paramstype, groups::GCParam, sites::SiteParam, idealmodel=BasicIdeal;
                    references=String[],
                    absolutetolerance=1E-12,
                    verbose=false)
        return _newmodelgc($name,params,groups,sites,idealmodel;references=references,absolutetolerance=absolutetolerance,verbose=verbose)
    end
    function $name(params::$paramstype, groups::GCParam, idealmodel=BasicIdeal;
                    references=String[],
                    absolutetolerance=1E-12, verbose=false)

        return _newmodelgc($name,params, groups,idealmodel;
                    references=references,
                    absolutetolerance=absolutetolerance,
                    verbose=verbose)
    end

    function Base.show(io::IO, mime::MIME"text/plain", model::$name)
        return eosshow(io, mime, model)
    end

    function Base.show(io::IO, model::$name)
        return eosshow(io, model)
    end
    end |> esc
end

const IDEALTYPE = Type{T} where T<:IdealModel

function _newmodelgc(eostype,params, groups::GCParam, sites::SiteParam, idealmodel::IDEALTYPE=BasicIdeal;
                    references::Vector{String}=String[],
                    absolutetolerance::Float64=1E-12,
                    verbose::Bool=false)

    #!(idealmodel <: IdealModel) && error("idealmodel ", idealmodel, " has to be a concrete subtype of IdealModel.")
    components = groups.components
    lengthcomponents = length(components)
    icomponents = 1:lengthcomponents

    flattenedgroups = groups.flattenedgroups
    lengthflattenedgroups = length(flattenedgroups)
    allcomponentnflattenedgroups = groups.allcomponentnflattenedgroups
    iflattenedgroups = 1:lengthflattenedgroups

    allcomponentgroups = groups.allcomponentgroups
    lengthallcomponentgroups = [length(allcomponentgroups[i]) for i in icomponents]
    allcomponentngroups = groups.allcomponentngroups
    igroups = [[findfirst(isequal(group), flattenedgroups) for group ∈ componentgroups] for componentgroups ∈ allcomponentgroups]

    allgroupsites = sites.allcomponentsites
    lengthallgroupsites = [length(groupsites) for groupsites ∈ allgroupsites]
    allgroupnsites = sites.allcomponentnsites
    isites = [1:lengthallgroupsites[k] for k ∈ iflattenedgroups]

    verbose && idealmodel != BasicIdeal && @info("Now creating ideal model ", idealmodel, ".")
    idealmodelstruct = idealmodel(components; verbose=verbose)

return eostype{idealmodel}(components, lengthcomponents, icomponents,
       allcomponentgroups, lengthallcomponentgroups, allcomponentngroups, igroups,
       flattenedgroups, lengthflattenedgroups, allcomponentnflattenedgroups, iflattenedgroups,
       allgroupsites, lengthallgroupsites, allgroupnsites, isites,
       params, idealmodelstruct, absolutetolerance, references)
end

function _newmodelgc(eostype,params, groups, idealmodel=BasicIdeal;
                    references=String[],
                    absolutetolerance=1E-12,
                    verbose=false)

    components = groups.components
    sites = SiteParam(components)
    return _newmodelgc(eostype,params, groups, sites, idealmodel; references=references, absolutetolerance=absolutetolerance, verbose=verbose)
end

"""

    @newmodel name parent paramstype

This is exactly the same as the above but for non-GC models.
All group parameters are absent in this struct.
The sites are associated to the main component rather than the groups,
and the respective fieldnames are named correspondingly.
"""
macro newmodel(name, parent, paramstype)
    quote
    struct $name{T <: IdealModel} <: $parent
        components::Array{String,1}
        lengthcomponents::Int
        icomponents::UnitRange{Int}

        allcomponentsites::Array{Array{String,1},1}
        lengthallcomponentsites::Array{Int,1}
        allcomponentnsites::Array{Array{Int,1},1}
        isites::Array{UnitRange{Int},1}

        params::$paramstype
        idealmodel::T
        absolutetolerance::Float64
        references::Array{String,1}
    end

    function $name(params::$paramstype, sites::SiteParam, idealmodel=BasicIdeal;
                    references=String[],
                    absolutetolerance=1E-12, verbose=false)

        return _newmodel($name,params,sites,idealmodel;references=references,absolutetolerance=absolutetolerance,verbose=verbose)
    end
    function $name(params::$paramstype, idealmodel=BasicIdeal;
                    references=String[],
                    absolutetolerance=1E-12,
                    verbose=false)

        return _newmodel($name,params, idealmodel; references=references, absolutetolerance=absolutetolerance)
    end

    function Base.show(io::IO, mime::MIME"text/plain", model::$name)
        return eosshow(io, mime, model)
    end

    function Base.show(io::IO, model::$name)
        return eosshow(io, model)
    end
    end |> esc
end

function _newmodel(eostype,params, sites::SiteParam, idealmodel::IDEALTYPE=BasicIdeal;
                references::Array{String,1}=String[],
                absolutetolerance::Float64=1E-12,
                verbose::Bool=false)

    #!(idealmodel <: IdealModel) && error("idealmodel ", idealmodel, " has to be a concrete subtype of IdealModel.")
    arbparam = arbitraryparam(params)
    components = arbparam.components
    lengthcomponents = length(components)
    icomponents = 1:lengthcomponents
    allcomponentsites = sites.allcomponentsites
    lengthallcomponentsites = [length(componentsites) for componentsites ∈ allcomponentsites]
    allcomponentnsites = sites.allcomponentnsites
    isites = [1:lengthallcomponentsites[i] for i ∈ icomponents]

    verbose && idealmodel != BasicIdeal && @info("Now creating ideal model ", idealmodel, ".")
    idealmodelstruct = idealmodel(components; verbose=verbose)

    return eostype{idealmodel}(components, lengthcomponents, icomponents,
        allcomponentsites, lengthallcomponentsites, allcomponentnsites, isites,
        params, idealmodelstruct, absolutetolerance, references)
end

function _newmodel(eostype,params, idealmodel=BasicIdeal;
                references=String[],
                absolutetolerance=1E-12,
                verbose=false)

    arbparam = arbitraryparam(params)
    components = arbparam.components
    sites = SiteParam(components)
    return _newmodel(eostype,params, sites, idealmodel; references=references, absolutetolerance=absolutetolerance)
end

"""
Even simpler model, primarily for the ideal models.
Contains neither sites nor ideal models.
"""
macro newmodelsimple(name, parent, paramstype)
    quote
    struct $name <: $parent
        components::Array{String,1}
        lengthcomponents::Int
        icomponents::UnitRange{Int}

        params::$paramstype
        absolutetolerance::Float64
        references::Array{String,1}
    end

    function $name(params::$paramstype;
                references::Array{String,1}=String[],
                absolutetolerance::Float64=1E-12, verbose::Bool=false)

        return _newmodelsimple($name,params;references=references,absolutetolerance=absolutetolerance,verbose=verbose)
    end

    function Base.show(io::IO, mime::MIME"text/plain", model::$name)
        return eosshow(io, mime, model)
    end

    function Base.show(io::IO, model::$name)
        return eosshow(io, model)
    end
    end |> esc
end

function _newmodelsimple(eostype,params;
                        references::Array{String,1}=String[],
                        absolutetolerance::Float64=1E-12,
                        verbose::Bool=false)

    arbparam = arbitraryparam(params)
    components = arbparam.components
    lengthcomponents = length(components)
    icomponents = 1:lengthcomponents

    return eostype(components, lengthcomponents, icomponents,
        params, absolutetolerance, references)
end


text:
"""
    BatchJob

Stores a batch job id in order to:

- `describe` a job and its parameters
- check on the `status` of a job
- `wait` for a job to complete
- fetch `log_events`

# Fields
- `id::AbstractString`: jobId
"""
@auto_hash_equals struct BatchJob
    id::AbstractString
end

"""
    submit(
        name::AbstractString,
        definition::JobDefinition,
        queue::JobQueue;
        container::AbstractDict=Dict(),
        parameters::Dict{String,String}=Dict{String, String}(),
        num_jobs::Integer=1,
    ) -> BatchJob

Handles submitting the batch job. Returns a `BatchJob` wrapper for the id.
"""
function submit(
    name::AbstractString,
    definition::JobDefinition,
    queue::JobQueue;
    container::AbstractDict=Dict(),
    parameters::Dict{String,String}=Dict{String, String}(),
    num_jobs::Integer=1,
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    debug(logger, "Submitting job \"$name\"")
    input = OrderedDict(
        "parameters" => parameters,
        "containerOverrides" => container,
    )

    if num_jobs > 1
        # https://docs.aws.amazon.com/batch/latest/userguide/array_jobs.html
        @assert 2 <= num_jobs <= 10_000
        push!(input, "arrayProperties" => Dict("size" => num_jobs))
    end

    debug(logger, "Input: $input")

    response = @mock Batch.submit_job(definition.arn, name, queue.arn, input; aws_config=aws_config)
    job = BatchJob(response["jobId"])

    if num_jobs > 1
        info(logger, "Submitted array job \"$(name)\" ($(job.id), n=$(num_jobs))")
    else
        info(logger, "Submitted job \"$(name)\" ($(job.id))")
    end

    return job
end

"""
    describe(job::BatchJob) -> Dict

Provides details about the AWS batch job.
"""
function describe(job::BatchJob; aws_config::AbstractAWSConfig=global_aws_config())
    response = @mock Batch.describe_jobs([job.id]; aws_config=aws_config)
    isempty(response["jobs"]) && error(logger, "Job $(job.id) not found.")
    debug(logger, "Job $(job.id): $response")
    return first(response["jobs"])
end

"""
    JobDefinition

Returns the job definition corresponding to a batch job.
"""
function JobDefinition(job::BatchJob; aws_config::AbstractAWSConfig=global_aws_config())
    JobDefinition(describe(job)["jobDefinition"]; aws_config=aws_config)
end

"""
    status(job::BatchJob) -> JobState

Returns the current status of a job.
"""
function status(job::BatchJob; aws_config::AbstractAWSConfig=global_aws_config())::JobState
    details = describe(job; aws_config=aws_config)
    return parse(JobState, details["status"])
end

"""
    status_reason(job::BatchJob) -> Union{String, Nothing}

A short, human-readable string to provide additional details about the current status of the
job.
"""
function status_reason(job::BatchJob; aws_config::AbstractAWSConfig=global_aws_config())
    details = describe(job; aws_config=aws_config)
    return get(details, "statusReason", nothing)
end

"""
    wait(
        cond::Function,
        job::BatchJob;
        timeout=600,
        delay=5
    )

Polls the batch job state until it hits one of the conditions in `cond`.
The loop will exit if it hits a `failure` condition and will not catch any excpetions.
The polling interval can be controlled with `delay` and `timeout` provides a maximum
polling time.

# Examples
```julia
julia> wait(state -> state < SUCCEEDED, job)
true
```
"""
function Base.wait(
    cond::Function,
    job::BatchJob;
    timeout=600,
    delay=5,
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    completed = false
    last_state = PENDING
    initial = true

    start_time = time()  # System time in seconds since epoch
    while time() - start_time < timeout
        state = status(job; aws_config=aws_config)

        if state != last_state || initial
            info(logger, "$(job.id) status $state")

            if !cond(state)
                completed = true
                break
            end

            last_state = state
        end

        initial && (initial = false)
        sleep(delay)
    end

    if !completed
        message = "Waiting on job $(job.id) timed out"

        if !initial
            message *= " Last known state $last_state"
        end

        throw(BatchJobError(job.id, message))
    end

    return completed
end

"""
    wait(
        job::BatchJob,
        cond::Vector{JobState}=[RUNNING, SUCCEEDED],
        failure::Vector{JobState}=[FAILED];
        kwargs...,
    )

Polls the batch job state until it hits one of the conditions in `cond`.
The loop will exit if it hits a `failure` condition and will not catch any excpetions.
The polling interval can be controlled with `delay` and `timeout` provides a maximum
polling time.
"""
function Base.wait(
    job::BatchJob,
    cond::Vector{JobState}=[RUNNING, SUCCEEDED],
    failure::Vector{JobState}=[FAILED];
    kwargs...,
)
    wait(job; kwargs...) do state
        if state in cond
            false
        elseif state in failure
            throw(BatchJobError(job.id, "Job $(job.id) hit failure condition $state"))
            false
        else
            true
        end
    end
end

"""
    log_events(job::BatchJob) -> Union{Vector{LogEvent}, Nothing}

Fetches the logStreamName, fetches the CloudWatch logs, and returns a vector of log events.
If the log stream does not currently exist then `nothing` is returned.

NOTES:
- The `logStreamName` isn't available until the job is RUNNING, so you may want to use
  `wait(job)` or `wait(job, [AWSBatch.SUCCEEDED])` prior to calling this function.
"""
function log_events(job::BatchJob)
    job_details = describe(job)

    if haskey(job_details["container"], "logStreamName")
        stream = job_details["container"]["logStreamName"]
    else
        return nothing
    end

    info(logger, "Fetching log events from $stream")
    return log_events("/aws/batch/job", stream)
end


text:
# Thrown when the string passed to `parseroman` is not a valid Roman numeral.
type InvalidRomanError <: Exception
    str::ASCIIString
end
Base.showerror(io::IO, err::InvalidRomanError) =
    print(io, err.str, " is not a valid Roman numeral")

immutable RomanNumeral <: Integer
    val::Int
    str::ASCIIString

    RomanNumeral(int::Integer) = new(int, toroman(int))
    RomanNumeral(str::ASCIIString) = begin
        num = parseroman(str)
        new(num, toroman(num))
    end
end
# This macro allows the rn"MMXV" syntactic sugar
macro rn_str(str)
    RomanNumeral(str)
end

typealias RN RomanNumeral

# Standard functions
# Conversion + promotion
Base.convert(::Type{Bool}, num::RN) = true
Base.convert(::Type{BigInt}, num::RN) = BigInt(num.val)
Base.convert{T<:Integer}(::Type{T}, num::RN) = convert(T, num.val)
Base.convert(::Type{RN}, num::Int) = RN(num)
Base.promote_rule{T<:Integer}(::Type{RN}, ::Type{T}) = T
Base.string(num::RN) = num.str

# IO
Base.print(io::IO, num::RN) = print(io, num.str)
Base.show(io::IO, num::RN) = write(io, num.str)

Base.length(num::RN) = length(num.str)
Base.hash(num::RN) = hash(num.str) $ hash(num.val)


text:
function upload_iris_all_formats_to_s3(bucket_name)
    iris_download_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"
    iris_local_path = download(iris_download_path)
    df = CSV.read(iris_local_path, DataFrames.DataFrame)
    # CSV
    verify_file_in_s3(
        bucket_name,
        "iris.csv",
        iris_download_path,
    )
    #df = collect(read_csv("s3://$(bucket_name)/iris.csv"))
    # Parquet
    Parquet.write_parquet("iris.parquet", df)
    verify_file_in_s3(
        bucket_name,
	"iris.parquet",
	p"iris.parquet"
    )
    # Arrow
    Arrow.write("iris.arrow", df)
    verify_file_in_s3(
        bucket_name,
	"iris.arrow",
	p"iris.arrow"
    )
end

@testset "Read/Write to CSV, Parquet, and Arrow" begin
    run_with_job("Reading/writing CSV/Parquet/Arrow") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_iris_all_formats_to_s3(bucket)
        for filetype in ["csv", "parquet", "arrow"]
            read_func = if filetype == "csv" read_csv
                elseif filetype == "parquet" BanyanDataFrames.read_parquet
                elseif filetype == "arrow" read_arrow
                end
            iris = read_func("s3://$(bucket)/iris.$(filetype)")
            @test nrow(iris) == 150
            @test ncol(iris) == 5
            @test collect(names(iris)) ==
                ["sepal_length", "sepal_width", "petal_length", "petal_width", "species"]
            @test collect(propertynames(iris)) ==
                [:sepal_length, :sepal_width, :petal_length, :petal_width, :species]
            iris_filter = iris[map(sl -> sl == 5.1, iris[:, :sepal_length]), :]
            iris_filter_first_row = first(collect(iris_filter))
            iris_filter_last_row = last(collect(iris_filter))
            @test iris_filter_first_row == DataFrames.DataFrameRow(
                DataFrames.DataFrame(
                    sepal_length = 5.1,
                    sepal_width = 3.5,
                    petal_length = 1.4,
                    petal_width = 0.2,
                    species = "setosa",
                ),
                1,
            )
            @test iris_filter_last_row == DataFrames.DataFrameRow(
                DataFrames.DataFrame(
                    sepal_length = 5.1,
                    sepal_width = 2.5,
                    petal_length = 3.0,
                    petal_width = 1.1,
                    species = "versicolor",
                ),
                1,
            )

            # Write filtered dataset to file
            write_func = if filetype == "csv" write_csv
                elseif filetype == "parquet" BanyanDataFrames.write_parquet
                elseif filetype == "arrow" write_arrow
                end
            write_func(iris_filter, "s3://$(bucket)/iris_filter.$(filetype)")
            # Read file and verify
            iris_new = read_func("s3://$(bucket)/iris_filter.$(filetype)")
            @test nrow(iris_new) == 9
            @test ncol(iris_new) == 5
            iris_new_last_row = last(collect(iris_new))
            @test iris_new_last_row == DataFrames.DataFrameRow(
                DataFrames.DataFrame(
                    sepal_length = 5.1,
                    sepal_width = 2.5,
                    petal_length = 3.0,
                    petal_width = 1.1,
                    species = "versicolor",
                ),
                1,
            )
        end
    end
    try
        s3_delete(get_cluster_s3_bucket_name(get_cluster_name()), "iris_new.csv")
        s3_delete(get_cluster_s3_bucket_name(get_cluster_name()), "iris_new.parquet")
        s3_delete(get_cluster_s3_bucket_name(get_cluster_name()), "iris_new.arrow")
    catch
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using Base.depwarn

Base.@deprecate long(data::Array{Float64},temp::Float64,wave::Float64) tlcorrection(data::Array{Float64},temp::Float64,wave::Float64;correction="long",density=2210.0)


text:
function constraint_converter_losses(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, a, b, c, plmax)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    pconv_dc = PowerModels.var(pm, n, cnd, :pconv_dc, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)

    PowerModels.con(pm, n, cnd, :conv_loss)[i] = @constraint(pm.model, pconv_ac + pconv_dc == a + b*iconv)
end

function constraint_conv_transformer(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    qtf_fr = PowerModels.var(pm, n, cnd, :qconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    phi = PowerModels.var(pm, n, cnd, :phi, acbus)
    phi_vmf = PowerModels.var(pm, n, cnd, :phi_vmf, i)
    va = PowerModels.var(pm, n, cnd, :va, acbus)
    vaf = PowerModels.var(pm, n, cnd, :vaf, i)
    cs = PowerModels.var(pm, n, cnd, :cs_vaf,i)

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gtf, btf, phi, phi_vmf, va, vaf, ptf_fr, ptf_to, qtf_fr, qtf_to, tm, cs)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cnd, cs, va, vaf)
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = c1
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = c2
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = c3
        PowerModels.con(pm, n, cnd, :conv_tf_q_to)[i] = c4
    else
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, ptf_fr + ptf_to == 0)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model, qtf_fr + qtf_to == 0)
        @constraint(pm.model, va == vaf)
        @constraint(pm.model, (1+phi) == (1+phi_vmf))
    end
end

"constraints for a voltage magnitude transformer + series impedance"

function lpac_power_flow_constraints(model, g, b, phi_fr, phi_to, va_fr, va_to, p_fr, p_to, q_fr, q_to, tm, cs)

    c1 = @constraint(model, p_fr ==  g/(tm^2)*(1.0 + 2*phi_fr) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*(va_fr-va_to))
    c2 = @constraint(model, q_fr == -b/(tm^2)*(1.0 + 2*phi_fr) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*(va_fr-va_to))
    c3 = @constraint(model, p_to ==  g*(1.0 + 2*phi_to) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*-(va_fr-va_to))
    c4 = @constraint(model, q_to == -b*(1.0 + 2*phi_to) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*-(va_fr-va_to))
    return c1, c2, c3, c4
end

function constraint_conv_reactor(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)
    ppr_to = -PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qpr_to = -PowerModels.var(pm, n, cnd, :qconv_ac, i)
    phi_vmc = PowerModels.var(pm, n, cnd, :phi_vmc, i)
    phi_vmf = PowerModels.var(pm, n, cnd, :phi_vmf, i)
    vac = PowerModels.var(pm, n, cnd, :vac, i)
    vaf = PowerModels.var(pm, n, cnd, :vaf, i)
    cs = PowerModels.var(pm, n, cnd, :cs_vac,i)

    phi_vmc_ub = JuMP.upper_bound(phi_vmc)
    ppr_to_ub = JuMP.upper_bound(PowerModels.var(pm, n, cnd)[:pconv_ac][i])
    qpr_to_ub = JuMP.upper_bound(PowerModels.var(pm, n, cnd)[:qconv_ac][i])
    Smax = sqrt(ppr_to_ub^2 + qpr_to_ub^2)
    zc = rc + im*xc
    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gc, bc, phi_vmf, phi_vmc, vaf, vac, ppr_fr, ppr_to, qpr_fr, qpr_to, 1, cs)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cnd, cs, vaf, vac)
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = c3
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = c4
        c6 = constraint_conv_capacity_PWL(pm, n, cnd, ppr_to, qpr_to, ppr_to_ub, qpr_to_ub, Smax)
   else
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == 0)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == 0)
        @constraint(pm.model, vac == vaf)
        @constraint(pm.model, (1+phi_vmf) == (1+phi_vmc))

    end
end

function constraint_conv_filter(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, bv, filter)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)
    phi_vmf = PowerModels.var(pm, n, cnd, :phi_vmf, i)

    PowerModels.con(pm, n, cnd, :conv_kcl_p)[i] = @constraint(pm.model, ppr_fr + ptf_to == 0 )
    PowerModels.con(pm, n, cnd, :conv_kcl_q)[i] = @constraint(pm.model, qpr_fr + qtf_to + -bv*filter*(1+2*phi_vmf) == 0)
end

function constraint_converter_current(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    phi_vmc = PowerModels.var(pm, n, cnd, :phi_vmc, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)

    # PowerModels.con(pm, n, cnd, :conv_i_sqrt)[i] = @NLconstraint(pm.model, pconv_ac^2 + qconv_ac^2 <= (1.0+phi_vmc)^2 * (iconv)^2)
    @constraint(pm.model, iconv <= Imax)
end

function variable_converter_filter_voltage(pm::AbstractLPACModel; kwargs...)
    variable_converter_filter_voltage_magnitude(pm; kwargs...)
    variable_converter_filter_voltage_angle_cs(pm; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::AbstractLPACModel; kwargs...)
    variable_converter_internal_voltage_magnitude(pm; kwargs...)
    variable_converter_internal_voltage_angle_cs(pm; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end

function variable_converter_filter_voltage_angle_cs(pm::AbstractLPACModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded = true)
    PowerModels.var(pm, nw, cnd)[:cs_vaf] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_cs_vaf",
    lower_bound =  0, # using the maximum bound allowed in cosine linearization. can be changed to bus angle +/- 10 degree?
    upper_bound =  1,
    start = 0
    )
end

function variable_converter_internal_voltage_angle_cs(pm::AbstractLPACModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded = true)
    PowerModels.var(pm, nw, cnd)[:cs_vac] = @variable(pm.model,
    [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_cs_vac",
    lower_bound =  0, # using the maximum bound allowed in cosine linearization. can be changed to bus angle +/- 10 degree?
    upper_bound =  1,
    start = 0
    )
end

function variable_converter_filter_voltage_magnitude(pm::AbstractLPACModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded = true)
    PowerModels.var(pm, nw, cnd)[:phi_vmf] = JuMP.@variable(pm.model,
            [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_phi_vmf",
            lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd) - 1.0,
            upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) - 1.0,
            start = comp_start_value(ref(pm, nw, :convdc, i), "phi_start", cnd)
        )
end

function variable_converter_internal_voltage_magnitude(pm::AbstractLPACModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded = true)
    PowerModels.var(pm, nw, cnd)[:phi_vmc] = JuMP.@variable(pm.model,
        [i in PowerModels.ids(pm, nw, :convdc)], base_name="$(nw)_$(cnd)_phi_vmc",
        lower_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmin", cnd) - 1.0,
        upper_bound = PowerModels.ref(pm, nw, :convdc, i, "Vmmax", cnd) - 1.0,
        start = comp_start_value(ref(pm, nw, :convdc, i), "phi_start", cnd)
        )
end

function constraint_conv_capacity_PWL(pm::AbstractLPACModel, n::Int, cnd::Int, ppr_to, qpr_to, Umax, Imax, Smax)
    np = 20 #no. of segments, can be passed as an argument later
    l = 0
    for i = 1:np
        a= Smax*sin(l)
        b = Smax*cos(l)
        c6 = @constraint(pm.model, a*ppr_to + b*qpr_to <= Smax^2) #current and voltage bounds to be proper to use Umax*Imax because Umax*Imax == Smax
        l = l + 2*pi/np
    end
end

function constraint_cos_angle_diff_PWL(pm::AbstractLPACModel, n::Int, cnd::Int, cs, va_fr, va_to)
    nb = 20 #no. of segments, can be passed as an argument later
    l = -pi/6
    h = pi/6
    inc = (h-l)/(nb+1)
    a = l + inc
    diff = va_fr - va_to
    for i = 1:nb
        c5 = @constraint(pm.model, cs <= -sin(a)*(diff-a) + cos(a))
        a = a + inc
    end
end

function add_dcconverter_voltage_setpoint(sol, pm::AbstractLPACModel)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmconv", :phi_vmc, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmfilt", :phi_vmf, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vaconv", :vac, status_name="islcc", inactive_status_value = 4)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vafilt", :vaf, status_name="islcc", inactive_status_value = 4)
end


text:
sanitize(s, dynvars::Array{Symbol,1}) = s
sanitize(s::Symbol, dynvars::Array{Symbol,1})= s in dynvars ? :($s(0)) : s
function sanitize(eq::Expr, dynvars::Array{Symbol,1})
    if eq.head == :(=)
        return sanitize( :( $(eq.args[1])==$( eq.args[2] ) ), dynvars )
    else
        if eq.head == :call
            if eq.args[1] in dynvars
                return eq
            else
                return Expr(eq.head, [sanitize(e, dynvars) for e in eq.args]...)
            end
        else
            head = eq.head
            args = eq.args
            return Expr(eq.head, [sanitize(e, dynvars) for e in args]...)
        end
    end
end

function sanitize(a,model::Dolo.ASModel)
    dynvars = get_variables(model)
    return sanitize(a, dynvars)
end

function Base.show(io::IO, m::SModel{ID}) where ID
    println(io,
    """NumericModel
      - name: $(m.name)
      - filename: $(m.name)
    """)
end

function Base.show(io::IO, ::MIME"text/html", model::Model)

    source = """

             Model

            name
            $(model.name)

            filename
            $(model.filename)

    """

    header = "
Type	Equation
\n"

    table_lines = []
    for eq_type in keys(model.equations)
        equations = model.equations[eq_type]
        if eq_type in [:controls_ub, :controls_lb]
            continue
        end
        for i=1:length(equations)
            if i==1
                eqt = eq_type
            else
                eqt = ""
            end
            eq = equations[i]
            eq = sanitize(eq, model)
            eqtex = Dolang.latex(eq)
            eqtex = replace(eqtex, "*"=>" ")
            fmt_line = "
$eqt	\\[$eqtex\\]
"
            push!(table_lines, fmt_line)
        end
    end

    source = string(source, "
", header, table_lines..., "
")
    print(io, source)
end


text:
withenv("FORCE_COLOR" => true) do
    if Sys.islinux()
        include("flicker.jl")
    end
    include("test_custompass.jl")
    include("test_bracketmatch.jl")
    include("test_highlighter.jl")
    include("test_rainbowbrackets.jl")
    include("test_logger.jl")
end


text:
#####
##### `sort`
#####

function frule((_, ẋs, _), ::typeof(partialsort), xs::AbstractVector, k; kw...)
    inds = partialsortperm(xs, k; kw...)
    return xs[inds], ẋs[inds]
end

function rrule(::typeof(partialsort), xs::AbstractVector, k::Union{Integer,OrdinalRange}; kwargs...)
    inds = partialsortperm(xs, k; kwargs...)
    ys = xs[inds]

    function partialsort_pullback(Δys)
        function partialsort_add!(Δxs)
            Δxs[inds] += Δys
            return Δxs
        end

        Δxs = InplaceableThunk(partialsort_add!, @thunk(partialsort_add!(zero(xs))))

        return NoTangent(), Δxs, NoTangent()
    end

    return ys, partialsort_pullback
end

function frule((_, ẋs), ::typeof(sort), xs::AbstractVector; kw...)
    inds = sortperm(xs; kw...)
    return xs[inds], ẋs[inds]
end

function rrule(::typeof(sort), xs::AbstractVector; kwargs...)
    inds = sortperm(xs; kwargs...)
    ys = xs[inds]

    function sort_pullback(ȳ)
        Δys = unthunk(ȳ)
        function sort_add!(Δxs)
            Δxs[inds] += Δys
            return Δxs
        end

        Δxs = InplaceableThunk(sort_add!, @thunk(sort_add!(zero(Δys))))

        return NoTangent(), Δxs
    end
    return ys, sort_pullback
end

#####
##### `sortslices`
#####

function frule((_, ẋ), ::typeof(sortslices), x::AbstractArray; dims::Integer, kw...)
    p = sortperm(collect(eachslice(x; dims=dims)); kw...)
    inds = ntuple(d -> d == dims ? p : (:), ndims(x))
    return x[inds...], ẋ[inds...]
end

function rrule(::typeof(sortslices), x::AbstractArray; dims::Integer, kw...)
    p = sortperm(collect(eachslice(x; dims=dims)); kw...)
    inds = ntuple(d -> d == dims ? p : (:), ndims(x))
    function sortslices_pullback(dy)
        # No actual need to zero this, and if you didn't, then you could widen eltype
        # Also, you could use similar(dy) here not x, same size?
        dx = _zerolike_writeat(x, unthunk(dy), (), inds...)
        return (NoTangent(), ProjectTo(x)(dx))
    end
    return x[inds...], sortslices_pullback
end

#####
##### `unique`
#####

function rrule(::typeof(unique), x::AbstractArray{<:Number}; dims=:)
    axes_x = axes(x)
    y = unique(x; dims=dims)  # accepts only dims=: or dims::Integer
    function unique_pullback(dy_raw)
        dy = unthunk(dy_raw)
        if length(x) == length(y)
            # Short-circuit for the case of all unique, since `mask` is fairly expensive:
            dx = reshape(dy, axes_x)
            return (NoTangent(), ProjectTo(x)(dx))
        end

        if dims isa Colon
            xs, ys = vec(x), y
        else
            xs, ys = collect(eachslice(x; dims=dims)), collect(eachslice(y; dims=dims))
        end
        mask = isequal.(permutedims(ys), xs)  # unique([0.0, -0.0, NaN, NaN])
        mask .= (mask .== cumsum(mask, dims=1) .== true)  # this implements  findfirst(mask; dims=1)
        keep = map(I -> I[1], findall(mask))
        if dims isa Colon
            # The function `_zerolike_writeat` allows second derivatives.
            # Should perhaps eventually be shared with `getindex`.
            dx = reshape(_zerolike_writeat(vec(x), vec(dy), (), keep), axes_x)
        else
            inds = ntuple(d -> d==dims ? keep : (:), length(axes_x))
            dx = _zerolike_writeat(x, dy, (), inds...)
        end
        return (NoTangent(), ProjectTo(x)(dx))
    end
    return y, unique_pullback
end


text:
#=
plot package initialization for Pseudospectra.jl

This file is part of Pseudospectra.jl.

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#
export setpsplotter, getpsplotter, setgs, isheadless, defaultgs

const _currentplotter = Ref{Symbol}(:undef)

const _available_plotters = Symbol[:PyPlot, :Plots,
                                   :Makie,
                                   ]
const _enabled_plotters = Dict{Symbol,Bool}(
    :PyPlot => false, :Plots => false,
    :Makie => false,
)

const _currentgs = Ref{Union{Nothing,GUIState}}(nothing)

function link_pyplot()
    include("PseudospectraMPL.jl")
    _enabled_plotters[:PyPlot] = true
end

function link_plots()
    include("PseudospectraPlots.jl")
    _enabled_plotters[:Plots] = true
end

function link_makie()
    include("PseudospectraMakie.jl")
    _enabled_plotters[:Makie] = true
end

getpsplotter() = _currentplotter[]

"""
    setpsplotter(plotter::Symbol=:default)

Select a plotting package for use with Pseudospectra.

Currently `:Plots` and `:PyPlot` are implemented.
Defaults to `:Plots` unless PyPlot is already imported without Plots.
"""
function setpsplotter(plotter::Symbol=:default)
    if plotter == :default
        if _enabled_plotters[:PyPlot] && !_enabled_plotters[:Plots]
            plotter = :PyPlot
        else
            plotter = :Plots
        end
    end
    if plotter ∉ _available_plotters
        throw(ArgumentError("plotter argument must be in $_available_plotters"))
    end
    if !_enabled_plotters[plotter]
        error("selected or default plotter '$plotter' is not enabled")
    end
    _currentplotter[]=plotter
    nothing
end

"""
    setgs(; headless=false, savefigs=true) => gs

Construct a `GUIState` for subsequent use by Pseudospectra functions.

Assumes plotting package has been chosen via `setpsplotter()`.
"""
function setgs(; headless=false, savefigs=true, fig_id=0)
    if _currentplotter[] ∈ [:default, :Plots]
        gs = PseudospectraPlots.PlotsGUIState(nothing,fig_id,nothing;
                                                   headless=headless,
                                                   savefigs=savefigs)
    elseif _currentplotter[] == :PyPlot
        gs = PseudospectraMPL.MPLGUIState(nothing,fig_id,nothing;
                                               headless=headless,
                                               savefigs=savefigs)
    elseif _currentplotter[] == :Makie
        gs = PseudospectraMakie.MakieGUIState(nothing,fig_id,nothing;
                                               headless=headless,
                                               savefigs=savefigs)
    else
        throw(ErrorException("use `setpsplotter` to establish plotting "
                             * "package for Pseudospectra first"))
    end
    _currentgs[] = gs
    return gs
end

"""
provides a default value for the `GUIState` object used by the Pseudospectra
package.  Attempts to initialize such an object if necessary.
"""
function defaultgs()
    if _currentgs[] === nothing
        if _currentplotter[] == :undef
            setpsplotter()
        end
        _currentgs[] = setgs()
    end
    _currentgs[]
end


text:
using AxisSets
using Documenter

makedocs(;
    modules=[AxisSets],
    authors="Invenia Technical Computing Corporation",
    repo="https://github.com/invenia/AxisSets.jl/blob/{commit}{path}#L{line}",
    sitename="AxisSets.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://invenia.github.io/AxisSets.jl",
        assets=String["assets/invenia.css"],
    ),
    pages=[
        "Home" => "index.md",
        "Example" => "example.md",
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/invenia/AxisSets.jl",
    devbranch = "main",
    push_preview = true,
)


text:
function sleepsort(arr::Vector{<:Real})
    out = Vector{eltype(arr)}(0)
    sizehint!(out, length(arr))
    @sync for x in arr
        @async begin
            sleep(x)
            push!(out, x)
        end
    end
    return out
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", sleepsort(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Batches

export batchsplit, sequential_batch_producer, random_batch_producer

using Base.Iterators
using Logging
using MLDataUtils

"""
    batchsplit(n, bsize)

Return number of batches and (possibly fixed) batch size based on a number of samples `n` and target batch size `bsize`.
"""

function batchsplit(n, bsize)
    nbatches = Int(ceil(n / bsize))
    if nbatches == 1
        bsize = n
    end
    nbatches, bsize
end

"""
    sequential_batch_producer(data...; bsize, epochs = 1, func = [identity, ..], log = false)

Create a channel producing batches of maximal size `bsize` based on `data` for the number of `epochs`.

Function `func` functions might be applied to each batch. Batches are created sequentially. Last batch
might be shorter then all the preceding ones.
"""
function sequential_batch_producer(data...; bsize, epochs = 1, func = [identity for _ in 1:length(data)], log = false)
    @assert length(data) > 0
    @assert (func isa Vector{<:Function}) || (func isa Function)
    if func isa Function
        func = [func for _ in 1:length(data)]
    end
    @assert length(data) == length(func)
    @assert length(unique(nobs.(data))) == 1
    n = nobs(data[1])

    nbatches, bsize = batchsplit(n, bsize)

    log && @info "sequential_batch_producer for $n x $epochs = $(epochs*n) samples, $nbatches batches per epoch"
    nbatches = 0
    Channel() do channel
        for _ in 1:epochs
            last = 0
            while last < n
                range = last+1:min(last+bsize, n)
    #             @info range, length(range)
                last += length(range)
                nbatches += 1
                if all(func .== identity)
                    d = tuple((view(X, (Colon() for _ in 1:ndims(X)-1)..., range) for X in data)...)
                else
                    d = tuple((f(X[(Colon() for _ in 1:ndims(X)-1)..., range]) for (f, X) in zip(func, data))...)
                end
                put!(channel, d)
            end
        end
    log && @info "sequential_batch_producer processed $nbatches batches"
    end
end

"""
    random_batch_producer(data...; bsize, epochs = 1, func = identity, log = false)

    Create a channel producing batches of size `bsize` based on `data` for the number of `epochs`.

Function `func` might be applied to each batch. Batches are sampled randomly (with repetition). When the number of instances `n`
is lower than `bsize`, batch size is reduced to `n`.
"""
function random_batch_producer(data...; bsize, epochs = 1, func = identity, log = false)
    @assert length(data) > 0
    @assert length(unique(nobs.(data))) == 1
    n = nobs(data[1])

    nbatches, bsize = batchsplit(n, bsize)

    data_ = flatten(repeated(RandomBatches(data, bsize, nbatches), epochs))
    log && @info "random_batch_producer for $n x $epochs = $(epochs*n) samples, $nbatches batches per epoch"
    Channel() do channel
        nbatches = 0
        for batch in data_
            put!(channel, func.(getobs.(batch)))
            nbatches += 1
        end
    end
end

end

max_stars_repo_name: dhonza/MLGoodies


text:
function camelize(cur_string::AbstractString)
  cur_string = classify(cur_string)

  cur_string = join(map(x -> lcfirst(x), split(cur_string, "/")), "/")

  cur_string
end


text:
# ## Derivatives

# The easiest way to provide derivatives is via Automatic Differentiation.
# We have imported `Zygote` as `Zyg`.
# For automatic differentiation we need custom adjoints for some `StaticArrays`:
Zyg.@adjoint (T::Type{<:SizedMatrix})(x::AbstractMatrix) = T(x), dv -> (nothing, dv)
Zyg.@adjoint (T::Type{<:SizedVector})(x::AbstractVector) = T(x), dv -> (nothing, dv)
Zyg.@adjoint (T::Type{<:SArray})(x::AbstractArray) = T(x), dv -> (nothing, dv)
# This allows us to define the following methods:

"Return the jacobian of `rbf` at `x` (using Zygote)."
function auto_jac( rbf :: RBFModel, x :: AbstractVector{<:Real} )
    Zyg.jacobian( rbf, x )[1]
end

"Evaluate the model and return the jacobian at the same time."
function eval_and_auto_jac( rbf :: RBFModel, x :: AbstractVector{<:Real} )
    y, back = Zyg._pullback( rbf, x )

    T = eltype(y)   # TODO does this make sense?
    n = length(y)
    jac = zeros(T, n, length(x) )
    for i = 1 : length(x)
        e = [ zeros(T, i -1 ); T(1); zeros(T, n - i )  ]
        jac[i, :] .= back(e)[2]
    end

    return y, jac
end

"Return gradient of output `ℓ` of model `rbf` at point `x` (using Zygote)."
function auto_grad( rbf :: RBFModel, x :: AbstractVector{<:Real}, ℓ :: Int = 1)
    Zyg.gradient( χ -> rbf(χ, ℓ), x )[1]
end

"Evaluate output `ℓ` of the model and return the gradient."
function eval_and_auto_grad( rbf :: RBFModel, x :: AbstractVector{<:Real}, ℓ :: Int = 1 )
    y, back = Zyg._pullback( χ -> rbf(χ, ℓ)[end], x)

    grad = back( one(y) )[2]
    return y, grad
end

# !!! note
#     We need at least `ChainRules@v.0.7.64` to have `auto_grad` etc. work for StaticArrays,
#     see [this issue](https://github.com/FluxML/Zygote.jl/issues/860).

# But we don't need `Zygote`, because we can derive the gradients ourselves.
# Assume that ``φ`` is two times continuously differentiable. \
# What is the gradient of a scalar RBF model?
# Using the chain rule and ``ξ = x - x^j`` we get
# ```math
    # \dfrac{∂}{∂ξ_i} \left( φ(\| ξ \|) \right)
    # =
    # φ\prime ( \| ξ \| ) \cdot
    # \dfrac{∂}{∂ξ_i} ( \| ξ \| )
    # =
    # φ\prime ( \| ξ \| ) \cdot
    # \dfrac{ξ_i}{\|ξ\|}.
# ```
# The right term is always bounded, but not well defined for ``ξ = 0``
# (see [^wild_diss] for details). \
# **That is why we require ``φ'(0) \stackrel{!}= 0``.** \
# We have ``\dfrac{∂}{∂x_i} ξ(x) = 1`` and thus
# ```math
    # ∇r(x) = \sum_{i=1}^N \frac{w_i φ\prime( \| x - x^i \| )}{\| x - x^i \|} (x - x^i) + ∇p(x)
# ```

# We can then implement the formula from above.
# For a fixed center ``x^i`` let ``o`` be the distance vector ``x - x^i``
# and let ``ρ`` be the norm ``ρ = \|o\| = \| x- x^i \|``.
# Then, the gradient of a single kernel is:
function grad( k :: ShiftedKernel, o :: AbstractVector{<:Real}, ρ :: Real )
    ρ == 0 ? zero(k.c) : (df( k.φ, ρ )/ρ) .* o
end

# In terms of `x`:
function grad( k :: ShiftedKernel, x :: AbstractVector{<:Real} )
    o = x - k.c     # offset vector
    ρ = norm2( o )  # distance
    return grad( k, o, ρ )
end

# The jacobian of a vector of kernels follows suit:
function jacT( K :: AbstractVector{<:ShiftedKernel}, x :: AbstractVector{<:Real})
    hcat( ( grad(k,x) for k ∈ K )... )
end
## precalculated offsets and distances, 1 per kernel
function jacT( K :: AbstractVector{<:ShiftedKernel}, offsets :: AbstractVector{<:AbstractVector}, dists :: AbstractVector{<:Real} )
    hcat( ( grad(k,o,ρ) for (k,o,ρ) ∈ zip(K,offsets,dists) )... )
end
jac( K :: AbstractVector{<:ShiftedKernel}, args... ) = transpose( jacT(K, args...) )

# Hence, the gradients of an RBFSum are easy:
function grad( rbf :: RBFSum, x :: AbstractVector{<:Real}, ℓ :: Int = 1 )
    #vec( jacT( rbf.kernels, x) * rbf.weights[:,ℓ] )
    vec( rbf.weights[ℓ,:]'jac( rbf.kernels, x ) )
end

function grad( rbf :: RBFSum, offsets :: AbstractVector{<:AbstractVector}, dists :: AbstractVector{<:Real}, ℓ :: Int)
    return vec( rbf.weights[ℓ,:]'jac( rbf.kernels, offsets, dists ) )
end

# The `grad` method looks very similar for the `PolySum`.
# We obtain the jacobian of the polynomial basis system via
# `PolynomialSystem.jacobian`.
function grad( psum :: PolySum, x :: AbstractVector{<:Real} , ℓ :: Int = 1)
    return vec( psum.weights[ℓ,:]'jacobian( psum.polys, x ))
end

# For the `RBFModel` we simply combine both methods:
function _grad( mod :: RBFModel, x :: AbstractVector{<:Real}, ℓ :: Int = 1 )
    return grad(mod.rbf, x, ℓ) + grad( mod.psum, x, ℓ )
end

function grad( mod :: RBFModel, x :: AbstractVector{<:Real}, ℓ :: Int = 1 )
    return _grad(mod, x, ℓ)
end

grad( mod :: RBFModel, x :: Vector{<:Real}, ℓ :: Int = 1 ) = _type_guard( _grad(mod, x, ℓ), Vector, mod.num_vars )
function grad( mod :: RBFModel, x :: T, ℓ :: Int = 1 ) where T <: Union{SVector, MVector, SizedVector}
    return _type_guard( _grad(mod, x, ℓ), T, mod.num_vars )
end

# We can exploit our custom evaluation methods for "distances":
function _offsets_and_dists( rbf :: RBFSum, x :: AbstractVector{<:Real} )
    offsets = [ x - k.c for k ∈ rbf.kernels ]
    dists = norm2.(offsets)
    return offsets, dists
end

function eval_and_grad( rbf :: RBFSum, offsets :: AbstractVector{<:AbstractVector}, dists :: AbstractVector{<:Real}, ℓ :: Int)
    return eval_at_dist( rbf, dists, ℓ ), grad( rbf, offsets, dists, ℓ)
end

function eval_and_grad( rbf :: RBFSum, x :: AbstractVector{<:Real}, ℓ :: Int = 1)
    offsets, dists = _offsets_and_dists(rbf, x)
    return eval_and_grad( rbf, offsets, dists, ℓ)
end

# For the `PolySum` we use `evaluate_and_jacobian`.
function eval_and_grad( psum :: PolySum, x :: AbstractVector{<:Real}, ℓ :: Int = 1)
    res_p, J_p = evaluate_and_jacobian( psum.polys, x )
    return (psum.weights[ℓ,:]'res_p)[1], vec(psum.weights[ℓ,:]'J_p)
end

# Combine for `RBFModel`:
function eval_and_grad( mod :: RBFModel, x :: AbstractVector{<:Real}, ℓ :: Int = 1 )
    res_rbf, g_rbf = eval_and_grad( mod.rbf, x, ℓ )
    res_polys, g_polys = eval_and_grad( mod.psum, x, ℓ )
    return res_rbf .+ res_polys, g_rbf .+ g_polys
end

# For the jacobian, we use the same trick to save evaluations.
function jac( rbf :: RBFSum, x :: AbstractVector{<:Real} )
    offsets, dists = _offsets_and_dists(rbf, x)
    rbf.weights * jac( rbf.kernels, offsets, dists )
end
jacT(rbf :: RBFSum, args... ) = transpose( jac(rbf, args...) )

function jac( psum :: PolySum, x :: AbstractVector{<:Real} )
    psum.weights * jacobian( psum.polys, x )
end

function _jac( mod :: RBFModel, x :: AbstractVector{<:Real} )
    jac( mod.rbf, x ) + jac( mod.psum, x)
end

jac( mod :: RBFModel, x :: AbstractMatrix{<:Real} ) = _jac(mod,x)
jac( mod :: RBFModel, x :: Vector{<:Real}) = convert(Matrix, _jac(mod,x))
jac( mod :: RBFModel, x :: SVector{<:Real}) = convert( SMatrix{mod.num_outputs, mod_nmu_vars}, _jac(mod,x) )
jac( mod :: RBFModel, x :: MVector{<:Real}) = convert( MMatrix{mod.num_outputs, mod_nmu_vars}, _jac(mod,x) )
jac( mod :: RBFModel, x :: SizedVector{<:Real}) = convert( SizedMatrix{mod.num_outputs, mod_nmu_vars}, _jac(mod,x) )

# As before, define an "evaluate-and-jacobian" function that saves evaluations:
function eval_and_jac( rbf :: RBFSum, x :: AbstractVector{<:Real} )
    offsets, dists = _offsets_and_dists(rbf, x)
    res = eval_at_dist( rbf, dists )
    J = rbf.weights * jac( rbf.kernels, offsets, dists )
    return res, J
end

function eval_and_jac( psum :: PolySum, x :: AbstractVector{<:Real} )
    res_p, J_p = evaluate_and_jacobian( psum.polys, x )
    return vec( psum.weights * res_p ), psum.weights * J_p
end

function eval_and_jac( mod :: RBFModel, x :: AbstractVector{<:Real} )
    res_rbf, J_rbf = eval_and_jac( mod.rbf, x )
    res_polys, J_polys = eval_and_jac( mod.psum, x)
    return res_rbf + res_polys, J_rbf + J_polys
end

## TODO type stable eval_and_grad and eval_and_jac ?

# !!! note
#     Hessians are not yet implemented.

# For the Hessian ``Hr \colon ℝ^n \to ℝ^{n\times n}`` we need the gradients of the
# component functions
# ```math
#     ψ_j(ξ) = \frac{ φ'( \left\| ξ \right\| )}{\|ξ\|} ξ_j
# ```
# Suppose ``ξ ≠ 0``.
# First, using the product rule, we have
# ```math
#    \dfrac{∂}{∂ξ_i}
#    \left(
#    \frac{ φ'( \left\| ξ \right\| )}{\|ξ\|} ξ_j
#    \right) =
#    ξ_j
#    \dfrac{∂}{∂ξ_i}
#    \left(
#    \frac{ φ'( \left\| ξ \right\| )}{\|ξ\|}
#    \right)
#    +
#    \frac{ φ'( \left\| ξ \right\| )}{\|ξ\|}
#    \dfrac{∂}{∂ξ_i}
#    ξ_j
# ```
# The last term is easy because of
# ```math
# \frac{∂}{∂ξ_i} ξ_j
# =
# \begin{cases}
#     1 & \text{if }i = j,\\
#     0 & \text{else.}
# \end{cases}
# ```
# For the first term we find
# ```math
#    \dfrac{∂}{∂ξ_i}
#    \left(
#      \frac{ φ'( \left\| ξ \right\| )}
#       {\|ξ\|}
#    \right)
#    =
#    \frac{
#        φ'\left(\left\| ξ \right\|\right) ∂_i \|ξ\|
#        - \|ξ\| ∂_i φ'\left( \left\| ξ \right\|\right)
#     }{
#         \|ξ\|^2
#     }
#     =
#     \frac{
#         \dfrac{φ'(\|ξ\|)}{\|ξ\|} ξ_i - \|ξ\|φ''(\|ξ\|)\dfrac{ξ_i}{\|ξ\|}
#     }{\|ξ\|^2}
# ```
# Hence, the gradient of ``ψ_j`` is
# ```math
#     ∇ψ_j(ξ)
#     =
#     \left( \frac{φ'(\|ξ\|)}{\|ξ\|^3}
#     -
#     \frac{φ''(\|ξ\|)}{\|ξ\|^2} \right) \cdot ξ
#     -\frac{φ'(\|ξ\|)}{\|ξ\|} e^j,
# ```
# where ``e^j ∈ ℝ^n`` is all zeros, except ``e^j_j = 1``.
# For ``ξ = 0`` the first term vanishes due to L'Hôpital's rule:
# ```math
# ∇ψ_j(0) = φ''(0) e^j.
# ```


text:
function lnd_rhs!(du, u, (parameters, F_accumulated), t)
    """
    Slab layer equation
        d(T_lnd)/dt = - (F_accumulated + G) / (h_lnd * ρ_lnd * c_lnd)
        where
            F_accumulated = F_integrated / Δt_coupler
    """
    @unpack lnd_h, lnd_ρ, lnd_c = parameters
    @unpack T_sfc = du

    @. T_sfc = (-F_accumulated) / (lnd_h * lnd_ρ * lnd_c)
end

# set up domain
function hspace_1D(xlim = (-π, π), npoly = 0, helem = 10)
    FT = Float64

    domain = Domains.IntervalDomain(Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]), periodic = true)
    mesh = Meshes.IntervalMesh(domain; nelems = helem)
    topology = Topologies.IntervalTopology(mesh)

    # Finite Volume Approximation: Gauss-Lobatto with 1pt per element
    quad = Spaces.Quadratures.GL{npoly + 1}()
    space = Spaces.SpectralElementSpace1D(topology, quad)

    return space
end

# init simulation
function lnd_init(; xmin = -1000, xmax = 1000, helem = 20, npoly = 0)

    # construct domain spaces - get only surface layer (NB: z should be zero, not z = first central height)
    space = hspace_1D((xmin, xmax), npoly, helem)
    coords = Fields.coordinate_field(space)
    domain = space

    # initial condition
    T_sfc = map(coords) do coord
        T_sfc = 273.0
    end

    # prognostic variable
    Y = Fields.FieldVector(T_sfc = T_sfc)

    return Y, domain
end


text:
module TestRecords

using RecordArrays
using Test

@testset "fill((a = 1, b = 2), 10)" begin
    @testset "$(sprint(show, align))" for align in [nothing, 64]
        A = RecordArrays.fill((a = 1, b = 2), 10; align = align)
        @test A.a == fill(1, 10)
        @test A.b == fill(2, 10)
        A.a .= 1:10
        A.b .= reverse(1:10)
        @test A[1] == (a = 1, b = 10)
        @test A[2] == (a = 2, b = 9)
        @test A.a == 1:10
        @test A.b == 10:-1:1
    end
end

struct OneValue{T}
    value::T
end

@testset "fill(OneValue{Union{Nothing,Int}}, 10)" begin
    Eltype = OneValue{Union{Nothing,Int}}
    @testset "$(sprint(show, align))" for align in [nothing, 64]
        A = RecordArrays.fill(Eltype(0), 10; align = align)
        @test A[1] === Eltype(0)
        @test A.value[1] === 0
        A[2] = Eltype(2)
        A[3] = Eltype(nothing)
        @test A[2] === Eltype(2)
        @test A[3] === Eltype(nothing)
        @test A.value[2] === 2
        @test A.value[3] === nothing
        A.value[4] = 4
        A.value[5] = nothing
        @test A[4] === Eltype(4)
        @test A[5] === Eltype(nothing)
        @test A.value[4] === 4
        @test A.value[5] === nothing
        vals = [[0, 2, nothing, 4, nothing]; fill(0, 5)]
        @test collect(A) == Eltype.(vals)
        @test collect(A.value) == vals
    end
end

# Maybe stop using `Vector{UInt8}` and use padding objects, so that we can put
# GC-managed objects in the array?
#=
struct Node{T}
    head::T
    tail::Union{Nothing,Node{T}}
end

list() = nothing
list(x, xs...) = Node(x, list(xs...))

@testset "fill(Node{Int}(...), 10)" begin
    @testset "$(sprint(show, align))" for align in [nothing, 64]
        A = RecordArrays.fill(list(1, 2), 10; align = align)
        @test A[1] === list(1, 2)
        @test A.head[1] === 1
        @test A.tail[1] === list(2)
        A[1] = list(1, 2, 3)
        @test A[1] === list(1, 2, 3)
        @test A.tail[1] === list(2, 3)
        @test A[2] === list(1, 2)
    end
end
=#

end  # module


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, opesciSLIM.TimeModeling, SeisIO

## Set up model structure
n = (120,100)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) + 0.4f0
v0 = ones(Float32,n) + 0.4f0
v[:,Int(round(end/3)):end] = 2.5f0
v[:,Int(2*round(end/3)):end] = 3.5f0
rho = ones(Float32,n)

# Slowness squared [s^2/km^2]
m = (1./v).^2
m0 = (1./v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2	# number of sources
model = Model(n,d,o,m)	# to include density call Model(n,d,o,m,rho)
model0 = Model(n,d,o,m0)

## Set up receiver geometry
nxrec = 120
xrec = linspace(50.,1150.,nxrec)
yrec = 0.
zrec = linspace(50.,50.,nxrec)

# receiver sampling and recording time
timeR = 1000.	# receiver recording time [ms]
dtR = 2.	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(linspace(100.,1100.,nsrc))
ysrc = convertToCell(linspace(0.,0.,nsrc))
zsrc = convertToCell(linspace(20.,20.,nsrc))

# source sampling and number of time steps
timeS = 1000.
dtS = 2.	# receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01
wavelet = source(timeS,dtS,f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(n),nsrc,ntComp)

######################## WITHOUT DENSITY ############################################

# Options structures
opt = Options(save_data_to_disk=false,
			  file_path=pwd(),	# path to files
			  file_name="shot_record"	# saves files as file_name_xsrc_ysrc.segy
			  )

opt0 = Options(save_data_to_disk=false,
			  file_path=pwd(),	# path to files
			  file_name="linearized_shot_record"	# saves files as file_name_xsrc_ysrc.segy
			  )

# Setup operators
Pr = joProjection(info,recGeometry)
F = joModeling(info,model;options=opt)
F0 = joModeling(info,model0;options=opt0)
Ps = joProjection(info,srcGeometry)
q = joData(srcGeometry,wavelet)

# Combined operator Pr*F*Ps'
Ffull = joModeling(info,model,srcGeometry,recGeometry)

# Nonlinear modeling
d1 = Pr*F*Ps'*q	# equivalent to d = Ffull*q
qad = Ps*F'*Pr'*d1

# Linearized modeling
J = joJacobian(Pr*F0*Ps',q)	# equivalent to J = joSetupJacobian(Ffull,q)
dD = J*dm
rtm = J'*dD

# Plot results
figure(); imshow(d1.data[1], vmin=-1e1, vmax=1e1)
figure(); imshow(reshape(rtm,model.n)', ColorMap("gray"))

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
# 一个码农的成长例程

struct MaNong
    name::String
    age::UInt8
    skills
    others
end

xiucai = MaNong("xiugai",
0,
["加班", "思路广", "欢乐多", "便宜", "好用", "活好", "话少", "英俊潇洒", "风流倜傥", "诗词若滔滔江水", "文赋如黄河泛滥", "……"],
"www.nagexiucai.com")

mutable struct CodePeasant
    name::String
    age::UInt8
    skills
    github::String
end

cp = CodePeasant("JulialangOrgCN", 1, "致力于Julia编程语言在zh-cn区域的推广【http://julialang.org.cn/】", "julialangorgcn")
println(cp.name)
cp.name = "朱华社区"
println(cp.name)

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia