text:
# November, 2020
# Kevin Dorma
# module for common PSV calculations
# Dependancies?
# GridInterpolations
# next, HEM

module PSVrelief

using DataFrames
using GridInterpolations

export PSVvaporFlux, PSVvaporRate, PSVvaporSize, PSVliquidSize, PSVliquidRate
export thermExpansionRate, poolFireReliefRate, PSVsteamFlux, PSVsteamRate, PSVsteamSize, liquidVaporizeReliefRate
export waterPsat, waterTsat, getKsh, PSVareaOrifice, PSVfindOrifice

# to do
# provide either PSV descriptor or area for rating.
# create function to return the next larger PSV descriptor
# implement proper steam equations
# implement HEM, Homogeneous Equilibrium Method for steam relief

function waterTsat(P_kPa)
    # equation for saturated steam, 1/T = A + B*ln P + C/ln P
    # I fit this myself, 100 to 20000 kPaa
    # pressure is in kPaa
    # we need to replace this with IAPWS IFC97 formulas

    lnP = log(P_kPa)

    # constants for fit of 1/T = A + B*ln P + C/ln P
    tA = 0.00379302
    tB = -0.000220828
    tC = -0.000425693

    invT = tA + tB*lnP + tC/lnP
    return((1.0/invT) - 273.15)
end

function waterPsat(T_c)
    # equation for saturated steam, ln P_Pa = A + B/Tk + C*Tk + D*ln(Tk)
    # I fit this myself, 100 to 20000 kPaa
    # pressure is in kPaa

    T_k = T_c + 273.15

    pA = 116.6408494
    pB = -8572.035364
    pC = 0.013736471
    pD = -14.73621925
    lnPsat = pA + pB/T_k + pC*T_k + pD*log(T_k)
    return(exp(lnPsat)/1000.0)

end

function getKsh(PkPa, State)
    # State is either a temperature in degC or a string
    # if this is a temperature, then look up the value in the table
    # if not, just return 1.0 because this is saturated
    # throw an error if the tempeature is below saturation
    # tables for superheat factors. I had a larger table but it threw weird errors. Smaller seems to be better.

    Ksh_tC = fill(0.0,(10))
    Ksh_pkPa = fill(0.0,(26))
    Ksh_table = fill(0.0, (26,10))    # pressure in rows, temperature in columns

# julia code for Ksh temperature values in degC
    Ksh_tC = [93.33333333	,
		148.8888889	,
		204.4444444	,
		260.0	,
		315.5555556	,
		371.1111111	,
		426.6666667	,
		482.2222222	,
		537.7777778	,
		565.5555556	];

    # julia code for Ksh pressure values in kPa
    Ksh_pkPa = [137.8951817	,
		344.7379543	,
		689.4759087	,
		1034.213863	,
		1378.951817	,
		1723.689772	,
		2068.427726	,
		2413.16568	,
		2757.903635	,
		3102.641589	,
		3447.379543	,
		3792.117498	,
		4136.855452	,
		4826.331361	,
		5515.807269	,
		6205.283178	,
		6894.759087	,
		7584.234995	,
		8273.710904	,
		8963.186813	,
		9652.662721	,
		10342.13863	,
		12065.8284	,
		13789.51817	,
		17236.89772	,
		20684.27726	];

    # Julia code for Ksh, rows are pressure, columns are temperature
    # subcooled values are denoted with KSH = 1.
    # Interpolation near the saturation temperature could give artificially low value for Ksh
    # a better method might be to replace the value 1 for the nearest subcooled temperature
    # with the value that gives 1 when interpolated to the saturation tempeature.
    # Clumsy, but it should work
    Ksh_table = [1	0.99455814	0.987	0.93	0.882	0.841	0.805	0.774	0.745	0.732	;
		1	0.997925224	0.987	0.93	0.882	0.841	0.805	0.774	0.745	0.732	;
		1	1	0.998	0.935	0.885	0.843	0.807	0.775	0.746	0.733	;
		1	1	0.984	0.94	0.888	0.846	0.808	0.776	0.747	0.733	;
		1	1	0.979	0.945	0.892	0.848	0.81	0.777	0.748	0.734	;
		1	1	1	0.951	0.895	0.85	0.812	0.778	0.749	0.735	;
		1	1	1	0.957	0.898	0.852	0.813	0.78	0.75	0.736	;
		1	1	1	0.963	0.902	0.854	0.815	0.781	0.75	0.736	;
		1	1	1	0.963	0.906	0.857	0.816	0.782	0.751	0.737	;
		1	1	1	0.961	0.909	0.859	0.818	0.783	0.752	0.738	;
		1	1	1	0.961	0.914	0.862	0.82	0.784	0.753	0.739	;
		1	1	1	0.962	0.918	0.864	0.822	0.785	0.754	0.74	;
		1	1	1	0.964	0.922	0.867	0.823	0.787	0.755	0.74	;
		1	1	1	1	0.931	0.872	0.827	0.789	0.757	0.742	;
		1	1	1	1	0.942	0.878	0.83	0.792	0.759	0.744	;
		1	1	1	1	0.953	0.883	0.834	0.794	0.76	0.745	;
		1	1	1	1	0.959	0.89	0.838	0.797	0.762	0.747	;
		1	1	1	1	0.962	0.896	0.842	0.8	0.764	0.749	;
		1	1	1	1	0.966	0.903	0.846	0.802	0.766	0.75	;
		1	1	1	1	0.973	0.91	0.85	0.805	0.768	0.752	;
		1	1	1	1	0.982	0.918	0.854	0.808	0.77	0.754	;
		1	1	1	1	0.993	0.926	0.859	0.811	0.772	0.755	;
		1	1	1	1	1	0.94	0.862	0.81	0.77	0.752	;
		1	1	1	1	1	0.952	0.861	0.805	0.762	0.744	;
		1	1	1	1	1	0.951	0.852	0.787	0.74	0.721	;
		1	1	1	1	1	1	0.831	0.753	0.704	0.684	];

        # we will use linear interpolation with P and T.
        # Using ln P and 1/T might be more robust and needs to be investigated.

    linear_Ksh = reshape(Ksh_table,(10*26,1));
    Ksh_grid = GridInterpolations.RectangleGrid(Ksh_pkPa, Ksh_tC);  	# rectangular grid

    Ksh = 1.0                # default value
    if (typeof(State)<:Number) # if the State is some kind of number, do the following
        Ksh = GridInterpolations.interpolate(Ksh_grid,linear_Ksh,[PkPa,State])
        # check if we are subcooled
        pSat = waterPsat(State)
        if (pSat < PkPa)
            error("Temperature is in subcooled region")
        end
    end
    return (Ksh)
end

function PSVareaOrifice(letter)
# psvTable[1,:areaMM2]
# from the PSV designation letter, output the PSV area in MM2
    psvTable = DataFrame(Designation = String[], Flanges = String[], areaIN2 = Float64[], areaMM2 = Float64[])
    push!(psvTable, ("D","1.5D2",0.11,70.9676))
    push!(psvTable, ("E","1.5E2",0.20,126.45136))
    push!(psvTable, ("F","1.5F3",0.31,198.06412))
    push!(psvTable, ("G","1.5G3",0.50,324.51548))
    push!(psvTable, ("H","2H3",0.79,506.4506))
    push!(psvTable, ("J","3J4",1.29,830.32092))
    push!(psvTable, ("K","3K4",1.84,1185.80408))
    push!(psvTable, ("L","4L6",2.85,1840.64148))
    push!(psvTable, ("M","4M6",3.60,2322.576))
    push!(psvTable, ("N","4N6",4.34,2799.9944))
    push!(psvTable, ("P","4P6",6.38,4116.1208))
    push!(psvTable, ("Q","6Q8",11.05,7129.018))
    push!(psvTable, ("R","6R10",16.00,10322.56))
    push!(psvTable, ("T","8T10",26.00,16774.16))

    return(psvTable[psvTable[:,:Designation] .== letter, :areaMM2])

end

function PSVfindOrifice(area)
    # from the required area (mm2), find the index for the next size up
    psvTable = DataFrame(Designation = String[], Flanges = String[], areaIN2 = Float64[], areaMM2 = Float64[])
    push!(psvTable, ("D","1.5D2",0.11,70.9676))
    push!(psvTable, ("E","1.5E2",0.20,126.45136))
    push!(psvTable, ("F","1.5F3",0.31,198.06412))
    push!(psvTable, ("G","1.5G3",0.50,324.51548))
    push!(psvTable, ("H","2H3",0.79,506.4506))
    push!(psvTable, ("J","3J4",1.29,830.32092))
    push!(psvTable, ("K","3K4",1.84,1185.80408))
    push!(psvTable, ("L","4L6",2.85,1840.64148))
    push!(psvTable, ("M","4M6",3.60,2322.576))
    push!(psvTable, ("N","4N6",4.34,2799.9944))
    push!(psvTable, ("P","4P6",6.38,4116.1208))
    push!(psvTable, ("Q","6Q8",11.05,7129.018))
    push!(psvTable, ("R","6R10",16.00,10322.56))
    push!(psvTable, ("T","8T10",26.00,16774.16))

    maxIndex = size(psvTable.Designation)[1]
    i = 1
    while ((psvTable.areaMM2[i] < area) & ( i < maxIndex))
        i += 1
    end

    return (psvTable[i,:Designation])
end

function thermExpansionRate(heat, alpha, heatCap)
    # liquid thermal expansion sizing
    # heat kJ/s,
    # alpha thermal expansion coefficient, 1/K
    # heatCap kJ/kg.K
    # return mass flow(kg/hr) = 3600 * alpha (1/K) * heat (kJ/s) / (heatCap (kJ/kg.K))
    return (3600.0 * alpha * heat / heatCap)
end

function liquidVaporizeReliefRate(heat,latent)
    # heat kJ/s
    # latent is in kJ/kg
    # flow rate kg/h = (heat/1000) * (1/latent) * (1/3600)

    flowRate = heat * (1/latent)* 3600
    return (flowRate)
end

function poolFireReliefRate(wettedAreaM2,latent,prompt)
    # wettedAreaM2 is in m2
    # latent is in kJ/kg
    # prompt fire response is either "prompt" or something else
    # heat Watts = C1 F A^0.82
    # flow rate kg/h = (heat/1000) * (1/latent) * (1/3600)

    C1 = 70900.0   # prompt fire fighting DOES NOT exist

    if (prompt == "prompt")
        C1 = 43200.0   # prompt fire fighting exists
    end

    F = 1.0 # no credits
    Q = C1*F*wettedAreaM2^0.82
    flowRate = (Q/1000) * (1/latent)* 3600
    return (flowRate)
end

function PSVsteamRate(areaMM2, Pkpa, State)
    # Napier equation for steam flow through PSV
    # Good for choked flow of steam through orifice
    # Pkpa is in kPaa
    # our function will return the flow rate in kg/hr
    # flowRate = Area * Flux
    psvFlux = PSVsteamFlux(Pkpa, State) # get the flux through the PSV, kg/hr.mm2
    Wkg = areaMM2*psvFlux
    return (Wkg)
end

function PSVsteamSize(Wkg, Pkpa, State)
    # Napier equation for steam flow through PSV
    # Good for choked flow of steam through orifice
    # our function will return the area in MM2
    # flowRate = Area * Flux
    psvFlux = PSVsteamFlux(Pkpa, State) # get the flux through the PSV, kg/hr.mm2
    Amm2 = Wkg/psvFlux
    return (Amm2)
end

# steam flux needs a field for State: Saturated or a temperature in C
function PSVsteamFlux(Pkpa, State)
    # Napier equation for steam flow through PSV
    # mass flow = area*flux
    # flux will be in kg/(mm2.hr)
    # Good for choked flow of steam through orifice
    # Wlbhr = Ain2 Cnapier Kd Ppsi Ksh Kb Kn
    # Wkghr = (1/2.205)*(1/25.4^2)*Amm2 Cnapier Kd Ppsi Ksh Kb Kn
    # flux = (1/2.205)*(1/25.4^2) Cnapier Kd Ppsi Ksh Kb Kn
    # our function will return the flow in kg/h
    # is the napier constant 51.45 or something else?
    # must agree on the value for Kd, 0.975 or 0.938
    # discharge coefficient is usually 0.975 for steam and vapour but
    # could be different depending on manufacturer
    # flowRate = Area * Flux

    Cnapier = 51.45
    Kd = 0.975
    Ksh = getKsh(Pkpa, State)
    Kb = 1.0
    Kn = 1.0     # low pressure and high pressure correction
    if (Pkpa > 10300.0)
        Kn = (2.7644*Pkpa/100.0 - 1000.0)/(3.3242*Pkpa/100.0 - 1061.0)
    end

    Ppsi = Pkpa * (14.503773800721813/100)
    return (Cnapier*Kd*Ppsi*Ksh*Kb*Kn/(2.205*25.4^2))

end

function PSVvaporFlux(P, Tcelcius, MW, k, Z)
    # vapour equation for PSV rating
    # P in kPaa
    # Tcelcius in celcius,
    # k is ratio of specific heats
    # Z compressibility factor
    # MW mole weight
    # Patm is atm pressure in kPaa
    # return value is flux in kg/mm2.hr
#    coeffSI = 13160.0; choop this out, not in the new version
    T = Tcelcius + 273.15;
    Kd = 0.975 # this could be reduced
    Kb = 1.0 # do not consider backpressure derating
    Kc = 1.0 # no derating for rupture disc
    C = apiC(k)
    rootTerm = sqrt(T*Z/MW);
    coeffTerm = 1.0 / (C * Kd * P * Kb * Kc)
    #    areaMM2 = W*coeffTerm*rootTerm;
    # W = areaMM2 * flux
    # flux = 1.0/(CoeffTerm * rootTerm)

    return (1.0/(coeffTerm*rootTerm))
end

function apiC(k)
    # SI form
    C = 0.03948*sqrt(k*(2.0/(k+1.0))^((k+1.0)/(k-1.0)))   # C coefficient API 520A fig 32
    return (C)
end

function PSVvaporRate(areaMM2, P, Tcelcius, MW, k, Z)
    # vapour equation for PSV rating
    # areaMM2 orifice area mm2
    # W flow rate in kg/h
    # P in kPaa
    # Tcelcius in celcius,
    # k is ratio of specific heats
    # Z compressibility factor
    # MW mole weight
    # Patm is atm pressure in kPaa
    # return value is rated flow rate kg/h
    # W = areaMM2 * flux
    # flux = 1.0/(CoeffTerm * rootTerm)
    psvFlux = PSVvaporFlux(P, Tcelcius, MW, k, Z) # get the flux through the PSV, kg/hr.mm2

    return (areaMM2 * psvFlux)
end

function PSVvaporSize(W, P, Tcelcius, MW, k, Z)
    # vapour equation for PSV sizing
    # W flow rate in kg/h
    # P in kPaa
    # Tcelcius in celcius,
    # k is ratio of specific heats
    # Z compressibility factor
    # MW mole weight
    # Patm is atm pressure in kPaa
    # return value is required area mm2
    # W = areaMM2 * flux
    # flux = 1.0/(CoeffTerm * rootTerm)
    psvFlux = PSVvaporFlux(P, Tcelcius, MW, k, Z) # get the flux through the PSV, kg/hr.mm2

    return (W / psvFlux)
end

function PSVliquidSize(W, P, Pback, d, mu)
    # liquid equation for PSV sizing, PSV requiring capacity certification
    # W flow rate in kg/h
    # P in kPag
    # Pback in kPag
    # d density kg/m3
    # mu viscosity cP or mPa.s
    # Patm is atm pressure in kPaa
    # return value is required area mm2
    coeffSI = 11.78;
    Q = 1000.0*W/(d*60.0);  # litres/minute
    G = d/1000.0; # specific gravity
    P1 = P;
    P2 = Pback;

    Kd = 0.65 # This is for a PSV. If we look for a rupture disk, Kd = 0.62
    Kw = 1.0 # do not consider backpressure derating
    Kc = 1.0 # no derating for rupture disc

    Kv = 1.0 # first pass at viscous correction, this will be iterated
    i = 1
    areaMM2 = 0.0;
    maxIter = 10;
    while ( i < maxIter)      # iterate 10 times, we are not going to check for convergence
        rootTerm = sqrt(G/(P1-P2));
        coeffTerm = coeffSI * Q / (Kd * Kw * Kc * Kv)
        areaMM2 = coeffTerm*rootTerm;
#        diam = sqrt(4.0*areaMM2/pi)/1000;
        R = Q * 18800 * G / (mu*sqrt(areaMM2)); # Reynolds number
        Kv = 1.0 / (0.9935 + 2.878/sqrt(R) + 342.75/(R^1.5))

        i += 1
    end
    return (areaMM2)
end

function PSVliquidRate(areaMM2, P, Pback, d, mu)
    # liquid equation for PSV sizing
    # W flow rate in kg/h
    # areaMM2 PSV area
    # P in kPag
    # Pback in kPag
    # d density kg/m3
    # mu viscosity cP or mPa.s
    # Patm is atm pressure in kPaa
    # return value is flow rate kg/h
    # I have not edited this function yet
    # from API 520
    # A = (11.78 Q) / (Kd Kw Kc Kv) * sqrt(G/DP)
    # Q = (A Kd Kw Kc Kv/11.78) sqqrt(DP/G)
    coeffSI = 11.78;
#    Q = 1000*W/(d*60.0);  # litres/minute
#    A = areaMM2 / (25.4^2)
    G = d/1000; # specific gravity
    P1 = P;
    P2 = Pback;

    Kd = 0.65 # typical value
    Kw = 1.0 # do not consider backpressure derating
    Kc = 1.0 # no derating for rupture disc

    Kv = 1.0 # first pass at viscous correction, this will be iterated
    i = 1
    Q = 0.0 # initialize the variable
    W = 0.0;  # initialize the variable
    maxIter = 10;
    while ( i < maxIter)                  # iterate 10 times, we are not going to check for convergence
        # areamm2 = coeff * root
        rootTerm = sqrt((P1-P2)/G);
        coeffTerm = (Kd * Kw * Kc * Kv) / coeffSI
        Q = areaMM2*coeffTerm*rootTerm;   # litres per minute
#        diam = sqrt(4.0*areaMM2/pi)/1000.0;
        R = Q * 18800.0 * G / (mu*sqrt(areaMM2));  # Reynolds number
        Kv = 1.0 / (0.9935 + 2.878/sqrt(R) + 342.75/(R^1.5))

        i += 1
    end
    return (Q*60*d/1000.0)  # mass flow rate kg/h
end

end # module

max_stars_repo_name: kevindorma/PSVrelief

max_issues_repo_name: kevindorma/PSVrelief

max_forks_repo_name: kevindorma/PSVrelief


text:
export AffineMap

"""
    AffineMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM},
              VN<:AbstractVector{NM}} <: AbstractAffineMap{N, S}

Type that represents an affine transformation ``M⋅X ⊕ v`` of a convex set ``X``.

### Fields

- `M` -- matrix/linear map
- `X` -- convex set
- `v` -- translation vector

### Notes

An affine map is the composition of a linear map and a translation. This type is
parametric in the coefficients of the linear map, `NM`, which may be different from
the numeric type of the wrapped set (`N`). However, the numeric type of the
translation vector should be `NM`.

### Examples

For the examples we create a ``3×2`` matrix, a two-dimensional unit square, and
a three-dimensional vector.
Then we combine them in an `AffineMap`.

```jldoctest constructors
julia> A = [1 2; 1 3; 1 4]; X = BallInf([0, 0], 1); b2 = [1, 2]; b3 = [1, 2, 3];

julia> AffineMap(A, X, b3)
AffineMap{Int64,BallInf{Int64,Array{Int64,1}},Int64,Array{Int64,2},Array{Int64,1}}([1 2; 1 3; 1 4], BallInf{Int64,Array{Int64,1}}([0, 0], 1), [1, 2, 3])
```

For convenience, `A` does not need to be a matrix but we also allow to use
`UniformScaling`s resp. scalars (interpreted as a scaling, i.e., a scaled
identity matrix).
Scaling by ``1`` is ignored and simplified to a pure `Translation`.

```jldoctest constructors
julia> using LinearAlgebra

julia> am = AffineMap(2I, X, b2)
AffineMap{Int64,BallInf{Int64,Array{Int64,1}},Int64,Diagonal{Int64,Array{Int64,1}},Array{Int64,1}}([2 0; 0 2], BallInf{Int64,Array{Int64,1}}([0, 0], 1), [1, 2])

julia> AffineMap(2, X, b2) == am
true

julia> AffineMap(1, X, b2)
Translation{Int64,Array{Int64,1},BallInf{Int64,Array{Int64,1}}}(BallInf{Int64,Array{Int64,1}}([0, 0], 1), [1, 2])
```

Applying a linear map to an `AffineMap` object combines the two maps into a new
`AffineMap` instance.
Again we can make use of the conversion for convenience.

```jldoctest constructors
julia> B = [2 0; 0 2]; am2 = B * am
AffineMap{Int64,BallInf{Int64,Array{Int64,1}},Int64,Array{Int64,2},Array{Int64,1}}([4 0; 0 4], BallInf{Int64,Array{Int64,1}}([0, 0], 1), [2, 4])

julia> 2 * am == am2
true
```

The application of an `AffineMap` to a `ZeroSet` or an `EmptySet` is simplified
automatically.

```jldoctest constructors
julia> AffineMap(A, ZeroSet{Int}(2), b3)
Singleton{Int64,Array{Int64,1}}([1, 2, 3])

julia> AffineMap(A, EmptySet{Int}(2), b3)
EmptySet{Int64}(2)
```
"""
struct AffineMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM},
                 VN<:AbstractVector{NM}} <: AbstractAffineMap{N, S}
    M::MAT
    X::S
    v::VN

    # default constructor with dimension match check
    function AffineMap(M::MAT, X::S, v::VN) where {N<:Real, S<:LazySet{N}, NM,
                                                   MAT<:AbstractMatrix{NM},
                                                   VN<:AbstractVector{NM}}

        @assert dim(X) == size(M, 2) "a matrix of size $(size(M)) cannot be " *
            "applied to a set of dimension $(dim(X))"

        @assert size(M, 1) == length(v) "a map with output dimension " *
            "$(size(M, 1)) is incompatible with a translation vector of " *
            "dimension $(length(v))"

        return new{N, S, NM, MAT, VN}(M, X, v)
    end
end

isoperationtype(::Type{<:AffineMap}) = true
isconvextype(::Type{<:AffineMap{N, S}}) where {N, S} = isconvextype(S)

# convenience constructor from a UniformScaling
function AffineMap(M::UniformScaling, X::LazySet, v::AbstractVector)
    return AffineMap(M.λ, X, v)
end

# convenience constructor from a scalar
function AffineMap(α::N, X::LazySet, v::AbstractVector) where {N<:Real}
    if α == one(N)
        return Translation(X, v)
    end
    return AffineMap(Diagonal(fill(α, length(v))), X, v)
end

# simplification for a LinearMap
for MAP in (:AbstractMatrix, :Real)
    @eval begin
        function LinearMap(map::$MAP, am::AffineMap)
            return AffineMap(map * am.M, am.X, map * am.v)
        end
    end
end

# ZeroSet is "almost absorbing" for the linear map (only the dimension changes)
# such that only the translation vector remains
function AffineMap(M::AbstractMatrix{N}, Z::ZeroSet{N}, v::AbstractVector{N}
                  ) where {N<:Real}
    @assert dim(Z) == size(M, 2) "a matrix of size $(size(M)) cannot be " *
        "applied to a set of dimension $(dim(Z))"

    @assert size(M, 1) == length(v) "a map with output dimension " *
        "$(size(M, 1)) is incompatible with a translation vector of " *
        "dimension $(length(v))"

    return Singleton(v)
end

# EmptySet is absorbing for AffineMap
function AffineMap(M::AbstractMatrix{N}, ∅::EmptySet{N}, v::AbstractVector{N}
                  ) where {N<:Real}
    return ∅
end

# --- AbstractAffineMap interface functions ---

function matrix(am::AffineMap)
    return am.M
end

function vector(am::AffineMap{N}) where {N<:Real}
    return am.v
end

function set(am::AffineMap)
    return am.X
end

# --- LazySet interface functions ---

function concretize(am::AffineMap)
    return affine_map(am.M, concretize(am.X), am.v)
end

ext: jl


text:
# !! auto-generated code, do not edit !!
module diagnostic_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct KeyValue <: Readable
    key::String
    value::String
end
struct DiagnosticStatus <: Readable
    level::UInt8
    name::String
    message::String
    hardware_id::String
    values::Vector{KeyValue}
end
struct DiagnosticArray <: Readable
    header::Header
    status::Vector{DiagnosticStatus}
end
end


text:
##########################################################################################
################################## Computing trajectories ################################
##########################################################################################

"""
	function get_trajectories(
		X;
		state = 0,
		fitness = missing,
		threshold = 0.05,
		time_scaling = 1,
	)

Assume rows of `X` are frequencies and time goes along columns, *i.e.* `X[:,1]` is the frequency series of the state of position `1`.
Also assume that `X[2*i-1,:]` (resp. `X[2*i,:]`) is the frequency of state `1` (resp. `-1`) at position `i`.
Same for `fitness` if provided.
"""
function get_trajectories(
	X;
	fitness = missing,
	threshold = 0.05,
	time_scaling = 1,
)
	@assert ismissing(fitness) || size(fitness) == size(X) "fitness and frequencies have inconsistent sizes"
	trajectories = Vector{Trajectory}(undef, 0)
	for i in 1:size(X, 2)
		pos = Int(ceil(i/2))
		state = iseven(i) ? -1 : 1
		append!(
			trajectories,
			Trajectory(
				X[:,i];
				pos,
				state,
				fitness = ismissing(fitness) ? missing : fitness[:, i],
				threshold,
				time_scaling,
			)
		)
	end
	return trajectories
end

"""
	function Trajectory(
		X::Vector{Float64};
		pos = 0,
		state = 0,
		fitness = nothing,
		threshold = 0.05,
		time_scaling = 1,
	)

## Arguments
- `X`: frequency vector

"""
function Trajectory(
	X::Vector{Float64};
	pos = 0,
	state = 0,
	fitness = missing,
	threshold = 0.05,
	time_scaling = 1,
)
	indices = get_traj_index(X, threshold)
	# Assert non overlapping trajectories
	for i in 1:length(indices)-1
		@assert isempty(intersect(indices[i], indices[i+1])) "Overlapping trajectories. Frequency series $X"
	end

	trajectories = Trajectory[]
	for t in indices
		final_state = if X[t[end]] < threshold
			:lost
		elseif X[t[end]] > 1-threshold
			:fixed
		else
			:missing
		end

		traj = Trajectory(;
			t = t * time_scaling,
			f = copy(X[t]),
			pos = pos,
			state = state,
			final_state = final_state,
			ϕpos = ismissing(fitness) ? ones(Missing, length(t)) : copy(fitness[t])
		)
		push!(trajectories, traj)
	end

	return trajectories
end

function get_traj_index(X, threshold)
	rx = Union{Missing, Float64}[missing, missing, missing] # last 3 values in X
	indices = []
	idx = Int64[]
	running = false
	for (i,x) in enumerate(X)
		shift!(rx, x)
		if is_start(rx, threshold, i)
			@assert !running
			running = true
		end

		running && push!(idx, i-1)

		if running && is_stop(rx, threshold, i, length(X))
			running = false
			push!(indices, idx)
			idx = Int64[]
		end
	end
	return indices
end

function is_start(rx, threshold, i)
	if i == 2
		# Special case: rx = [missing, val1, val2]
		return rx[2] < threshold && rx[3] >= threshold
	elseif i > 2
		return rx[1] < threshold && rx[2] < threshold && rx[3] >= threshold
	end
	return false
end
function is_stop(rx, threshold, i, L)
	if i == L
		return true
	elseif i > 2
		if rx[2] < threshold && rx[3] < threshold
			return true
		elseif rx[2] > 1-threshold && rx[3] > 1-threshold
			return true
		end
	end
	return false
end

function shift!(rx, x)
	rx[1] = rx[2]
	rx[2] = rx[3]
	rx[3] = x
	return nothing
end

##########################################################################################
################################## Fixation probability ##################################
##########################################################################################

"""
	fixation_probability(f::Function, trajectories::Vector{Trajectory})

Compute fraction of trajectories `T` for which `f(T)` is true and that fix.
"""
function fixation_probability(f::Function, trajectories::Vector{Trajectory})
	nlost = count(T -> f(T) && T.final_state == :lost, trajectories)
	nfix = count(T -> f(T) && T.final_state == :fixed, trajectories)
	return nfix / (nfix + nlost)
end

"""
	fixation_probability(trajectories::Vector{Trajectory}, fb::FrequencyBin)
	fixation_probability(trajectories, f, df)
"""
function fixation_probability(trajectories::Vector{Trajectory}, fb::FrequencyBin)
	return fixation_probability(T -> inbin(T, fb), trajectories)
end
fixation_probability(trajectories, f, df) = fixation_probability(trajectories, FrequencyBin(f, df))

"""
	fixation_probability(trajectories, fb::FrequencyBin, ϕmin)
	fixation_probability(trajectories, f, df, ϕmin)
"""
function fixation_probability(trajectories, fb::FrequencyBin, ϕmin)
	return fixation_probability(T -> has_fitness_above(T, ϕmin, fb), trajectories)
end
function fixation_probability(trajectories, f, df, ϕmin)
	return fixation_probability(trajectories, FrequencyBin(f, df), ϕmin)
end


text:
module OPRBF
#
# Module OPRBF implements functions related to radial basis function fitting
# and interpolation.
#

export RBFmat, RBFfit, RBFkernel, RBFeval, RBFcondition, RBFderiv

using LinearAlgebra

# ===================================================================
function RBFmat(x, s; ker=1)
#
# Compute the matrix used in evaluating the RBF fit.
#
# Version:        Changes:
# --------        -------------
# 13.06.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 13.06.2021
#
# Inputs:
# -------
# x               : Reference points, Nd-by-Np matrix.
# s               : Radial basis function scale parameter.
# ker             : Flag selecting the type of kernel function: see
#                   RBFkernel.
#
# Outputs:
# --------
# M               : The matrix, Np-by-Np.

   Np = size(x,2)

   # Set up the system of equations.
   M = zeros(Np,Np)
   for irow = 1:Np
      M[irow,irow] = RBFkernel(0.0, 1.0, s, ker=ker)
      for icol = irow+1:Np
         r = norm(x[:,irow] - x[:,icol])
         M[irow,icol] = RBFkernel(r, 1.0, s, ker=ker)
         M[icol,irow] = M[irow,icol]
      end
   end

   return M

end  # RBFmat
function RBFmat(x1, x2, s; ker=1)
#
# This version of RBFmat generates the matrix between two sets of
# points that may be different.
#
# Version:        Changes:
# --------        -------------
# 24.06.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 24.06.2021
#
# Inputs:
# -------
# x1,x2           : Reference points, Nd-by-N1 and Nd-by-N2 matrices.
# s               : Radial basis function scale parameter.
# ker             : Flag selecting the type of kernel function: see
#                   RBFkernel.
#
# Outputs:
# --------
# M               : The matrix, N1-by-N2.

   N1 = size(x1,2)
   N2 = size(x2,2)

   # Set up the system of equations.
   M = zeros(N1,N2)
   for irow = 1:N1
      for icol = 1:N2
         r = norm(x1[:,irow] - x2[:,icol])
         M[irow,icol] = RBFkernel(r, 1.0, s, ker=ker)
      end
   end

   return M

end  # RBFmat

# ===================================================================
function RBFfit(x, y, s; ker=1, con=false)
#
# Fit a radial basis function of the specified kernel to the given points.
#
# [Note, could also make a version that takes different x and xp points,
#  where Nx > Nxp, performing a least-squares solution at the x points
#  for RBF coefficients at the xp points.]
#
# Version:        Changes:
# --------        -------------
# 10.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 10.05.2021      Checked an example 2D surface-fitting problem.
#
# Inputs:
# -------
# x               : Reference points, Nd-by-Np matrix.
# y               : Scalar values at the reference points, Vector of
#                   length Np.
# s               : Radial basis function scale parameter.
# ker             : Flag selecting the type of kernel function: see
#                   RBFkernel.
# con             : Set to true to compute and return the matrix
#                   condition number.
#
# Outputs:
# --------
# w               : Weights, vector of length Np.

   Np = size(x,2)

   # Set up the system of equations.
   M = RBFmat(x, s, ker=ker)

   if (con)
      c = cond(M)
   else
      c = 0.0
   end

   w = M \ y

   return w, c

end  # RBFfit

# ===================================================================
function RBFkernel(r, w, s; ker=1)
#
# Compute the kernel at the given points.
#
# Version:        Changes:
# --------        -------------
# 10.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 10.05.2021      Checked an example 2D surface-fitting problem.
#
# Inputs:
# -------
# r               : Distances at which to evaluate the kernel function.
#                   Vector of length Np.
# w               : Kernel function weights.
# s               : Kernel function scale parameter.
# ker             : Flag specifying the type of kernel:
#                   1: Gaussian.
#                   2: Multiquadratic.
#                   3: Inverse multiquadratic.
#                   4: Student T with nu = 1.
#
# Outputs:
# --------
# y               : Kernel function values.  Vector of length Np.

   if (ker == 1)      # Gaussian
      y = w .* exp.(-(s .* r) .^ 2)
   elseif (ker == 2)  # Quadratic
      y = w .* sqrt.(1.0 .+ (s .* r) .^ 2)
   elseif (ker == 3)  # Inverse quadratic
      y = w ./ sqrt.(1.0 .+ (s .* r) .^ 2)
   elseif (ker == 4)  # Student T, nu = 1
      y = s .* w ./ (pi*(1.0 .+ (s .* r) .^ 2))
   end

   return y

end  # RBFkernel

# ===================================================================
function RBFeval(x, xp, w, s; ker=1)
#
# Compute the value of the function at the given points.
#
# Version:        Changes:
# --------        -------------
# 10.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 10.05.2021      Checked an example 2D surface-fitting problem.
#
# Inputs:
# -------
# x               : Points at which the function is desired.  Nd-by-
#                   Nx matrix.
# xp              : Points at which the weights are defined.  Nd-by-
#                   Np matrix.
# w               : Kernel function weights.  Vector of length Np.
# s               : Kernel function scale parameter.
# ker             : Flag specifying the type of kernel: see RBFkernel.
#
# Outputs:
# --------
# f               : Function values.  Vector of length Nx.

   Nx = size(x,2)
   Np = size(xp,2)
   Nd = size(x,1)

   f = zeros(Nx)
   for ix in 1:Nx
      r = zeros(Np)
      for id in 1:Nd
         r = r .+ (x[id,ix] .- xp[id,:]).^2
      end
      r = sqrt.(r)
      y = RBFkernel(r, w, s; ker=ker)
      f[ix] = sum(y)

   end

   return f

end  # RBFeval

# ===================================================================
function RBFcondition(x, y, sg; clow=5.0e4, chigh=2.0e5, ker=1)
#
# Iterate on the scale parameter s until the matrix condition number
# lies within the specified bounds.
#
# Version:        Changes:
# --------        -------------
# 12.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 12.05.2021      Checked an example 2D surface-fitting problem.
#
# Inputs:
# -------
# x               : Referance points, Nd-by-Np matrix.
# y               : Scalar values at the reference points, Vector of
#                   length Np.
# sg              : Guess for the radial basis function scale parameter.
# clow            : Lower bound for the matrix condition number.
# chigh           : Upper bound for the matrix condition number.
# ker             : Flag selecting the type of kernel function: see
#                   RBFkernel.
#
# Outputs:
# --------
# w               : Weights, vector of length Np.
# s               : Value of the scale parameter satisfying the
#                   condition number.
# c               : The final condition number.

   s = deepcopy(sg)
   w,c = RBFfit(x, y, s; ker=ker, con=true)

   # First bracket the desired condition number.
   if (c < clow)

      # The condition number is low.  Make the scale parameter smaller,
      # which makes the basis function broader.  (Sorry for the confusion,
      # but this is the typical definition.)
      sl = deepcopy(s)
      sh = deepcopy(s)
      iter = 0
      while (c < clow) && (iter < 10)
         sl = 0.5*sl
         w,c = RBFfit(x, y, sl; ker=ker, con=true)
         iter = iter + 1'
      end

      if (c > chigh)
         flg = false
      else
         # The upper bound now lies within the acceptable range.  If
         # iter reached its limit, well, this is anyways the best that
         # we are going to get, so no additional logic is needed.
         s = deepcopy(sl)
         flg = true
      end

   elseif (c > chigh)

      # The condition number is high.  Make the scale parameter larger,
      # which makes the basis function narrower.
      sl = deepcopy(s)
      sh = deepcopy(s)
      iter = 0
      while (c > chigh) && (iter < 10)
         sh = 2.0*sh
         w,c = RBFfit(x, y, sh; ker=ker, con=true)
         iter = iter + 1
      end

      if (c < clow)
         flg = false
      else
         # The lower bound now lies within the acceptable range.
         s = deepcopy(sh)
         flg = true
      end

   else
      # The initial s lies within the acceptable range.
      flg = true
   end

   # Now, if flg is false, cl and ch bracket the acceptable range.  Do
   # a bisecting search.
   iter = 0
   while (!flg) && (iter < 10)
      s = 0.5*(sh + sl)
      w,c = RBFfit(x, y, s; ker=ker, con=true)
      if (c < clow)
         # The condition number is low.  Make the scale parameter smaller.
         sh = deepcopy(s)
      elseif (c > chigh)
         # The condition number is high.  Make the scale parameter larger.
         sl = deepcopy(s)
      else
         flg = true
      end
      iter = iter + 1
   end

   return w,s,c

end  # RBFcondition

# ===================================================================
function RBFderiv(x, xp, w, s, k; ker=1)
#
# Compute the derivatives of the radial basis functions with respect
# to each of the coordinates, in a specified sequence.
#
# Note that this function does not compute df/dx_1, df/dx_2, ... but
# rather df/dx_1, d^2f/dx_1dx_2, ...  If the former is desired,
# then call this function multiple times with the appropriate values
# of k.
#
# Version:        Changes:
# --------        -------------
# 16.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 16.05.2021      Checked a 2D example problem, comparing with finite
#                 difference using RBFeval.
#
# Inputs:
# -------
# x               : Points at which the function is desired.  Nd-by-
#                   Nx matrix.
# xp              : Points at which the weights are defined.  Nd-by-
#                   Np matrix.
# w               : Kernel function weights.
# s               : Kernel function scale parameter.
# k               : Sequence in which derivatives are to be taken.
#                   Vector of length Ndim, containing the integers
#                   1:Nd in the desired order of differentiation.
#                   Example: d^3f/dx_2 dx_1 dx_3; k = [2, 1, 3]
#                   (It doesn't HAVE to contain all the integers from
#                   1 to Ndim, but there cannot be any repeated
#                   indices.)
# ker             : Flag specifying the type of kernel:
#                   1: Gaussian.
#                   2: Multiquadratic.
#                   3: Inverse multiquadratic.
#                   4: Student T with nu = 1.
#
# Outputs:
# --------
# f               : Values of the RBFs, same as a call to RBFeval.
#                   These are computed here regardless, so might as
#                   well return them as output.
# df              : Matrix of derivatives, dimensions Nk-by-Nx,
#                   containing df/dx_k[1], d^2f/dx_k[1] dx_k[2], ...
#                   d^Nk f/dx_k[1] dx_k[2] ... dx_k[Nk].

   Nx = size(x,2)
   Np = size(xp,2)
   Nd = size(x,1)
   Nk = size(k,1)

   # Build up the formulas recursively.
   f  = zeros(Nx)
   df = zeros(Nk,Nx)

   for ix in 1:Nx
      r = zeros(Np)
      for id in 1:Nd
         r = r .+ (x[id,ix] .- xp[id,:]).^2
      end
      r = sqrt.(r)
      y = RBFkernel(r, w, s; ker=ker)
      if (ker == 1)      # Gaussian
         v  =  s^2
         c  = -2.0
         dc =  0.0
      elseif (ker == 2)  # Quadratic
         v  =  (s^2) ./ (1.0 .+ (s .* r).^2)
         c  =  1.0
         dc = -2.0
      elseif (ker == 3)  # Inverse quadratic
         v  =  (s^2) ./ (1.0 .+ (s .* r).^2)
         c  = -1.0
         dc = -2.0
      elseif (ker == 4)  # Student T.
         v  =  (s^2) ./ (1.0 .+ (s .* r).^2)
         c  = -2.0
         dc = -2.0
      end
      f[ix] = sum(y)
      for ik in 1:Nk
         y = y .* (c * v .* (x[k[ik],ix] .- xp[k[ik],:]))
         df[ik,ix] = sum(y)
         c = c + dc
      end

   end

   return f, df

end # RBFderiv

end  # module.

max_stars_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_issues_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_forks_repo_name: SINTEF-Energy-Wind/OutliersPrimer


text:
export grid, gridrects

function grid(xiter, yiter)
    xs = (x for x in xiter, _ in yiter)
    ys = (y for _ in xiter, y in yiter)
    xs, ys
end

export LinAxis
struct LinAxis
    # rect::Rect
    xlim::Tuple{Float64, Float64}
    ylim::Tuple{Float64, Float64}
end

function LinAxis(xlim::Tuple{Real, Real}, ylim::Tuple{Real, Real}, xmargin, ymargin)
    xrange = xlim[2] - xlim[1]
    xlim_m = (xlim[1] - xmargin * xrange / 2, xlim[2] + xmargin * xrange / 2)
    yrange = ylim[2] - ylim[1]
    ylim_m = (ylim[1] - ymargin * yrange / 2, ylim[2] + ymargin * yrange / 2)
    LinAxis(xlim_m, ylim_m)
end

function (a::LinAxis)(r::Rect, p::Point)
    xscaled = (p.x - a.xlim[1]) / (a.xlim[2] - a.xlim[1])
    yscaled = (p.y - a.ylim[1]) / (a.ylim[2] - a.ylim[1])
    P(r, xscaled, yscaled)
end

function (a::LinAxis)(r::Rect, x, y)
    xscaled = (x - a.xlim[1]) / (a.xlim[2] - a.xlim[1])
    yscaled = (y - a.ylim[1]) / (a.ylim[2] - a.ylim[1])
    P(r, xscaled, yscaled)
end

Base.Broadcast.broadcastable(a::LinAxis) = Ref(a)

function gridrects(r::Rect, nrows, ncols, margins, rowgaps, colgaps, height_ratios, width_ratios)
    w = r.width
    h = r.height

    colswidth = (w - margins[1] - margins[3] - sum(colgaps))
    rowsheight = (h - margins[2] - margins[4] - sum(rowgaps))

    colwidths = width_ratios .* colswidth / sum(width_ratios)
    rowheights = height_ratios .* rowsheight / sum(height_ratios)

    xcenters = margins[1] .+ cumsum(colwidths) .- (colwidths ./ 2) .+ cumsum([0, colgaps...]) .+ r.center.x .- (r.width / 2)
    ycenters = margins[2] .+ cumsum(rowheights) .- (rowheights ./ 2) .+ cumsum([0, rowgaps...]) .+ r.center.y .- (r.height / 2)

    centers = P.(grid(xcenters, ycenters)...)

    Rect.(centers, grid(colwidths, rowheights)..., deg(0))
end

function gridrects(r::Rect, nrows, ncols, margin, rowgap, colgap)
    gridrects(
        r, nrows, ncols,
        [margin for _ in 1:4],
        [rowgap for _ in 1:nrows-1],
        [colgap for _ in 1:ncols-1],
        [1 for _ in 1:nrows],
        [1 for _ in 1:ncols])
end


text:
module Logger
import HttpServer.mimetypes

function print_rgb(r, g, b, t)
    print("\e[1m\e[38;2;$r;$g;$b;249m",t)
end

random_color(msg) = print_rgb(rand(1:256), rand(1:256), rand(1:256), msg)

"""
# simple
assumes that the proper engine middleware is included
"""
function simple(ap::Function, ctx)
  s = time()
  output =  ap(ctx)
  f = time()
  random_color("[$(ctx.data[:request][:method]) $(join(ctx.data[:request][:path], "/") * "/")] -- $(f - s) (s)\n")
  output
end

end


text:
name = tempname()
xdoty_tex = L"x \cdot y"

#= This test fails after updating dvisvgm, can this functionality be tested in a less dependant way?
if Sys.islinux()
    render(xdoty_tex, MIME("image/svg"); name=name, open=false)
    svg = open("$(name).svg") do f
        read(f, String)
    end
    @test svg == raw"""
    <?xml version='1.0' encoding='UTF-8'?>
    <!-- This file was generated by dvisvgm 2.10 -->
    <svg version='1.1' xmlns='http://www.w3.org/2000/svg' xmlns:xlink='http://www.w3.org/1999/xlink' width='25.707544pt' height='9.165631pt' viewBox='-72.000004 -72.000009 25.707544 9.165631'>
    <defs>
    <path id='g0-1' d='M2.75447-3.586549C2.75447-4.002589 2.410161-4.346897 1.994121-4.346897S1.233773-4.002589 1.233773-3.586549S1.578082-2.826201 1.994121-2.826201S2.75447-3.170509 2.75447-3.586549Z'/>
    <path id='g1-120' d='M6.800097-5.853248C6.341018-5.767171 6.168864-5.422862 6.168864-5.150284C6.168864-4.805976 6.441442-4.691206 6.642289-4.691206C7.072674-4.691206 7.373945-5.064207 7.373945-5.451554C7.373945-6.054095 6.685327-6.326672 6.082787-6.326672C5.207669-6.326672 4.719898-5.465901 4.590783-5.193323C4.26082-6.269288 3.371356-6.326672 3.113124-6.326672C1.649813-6.326672 .875118-4.447321 .875118-4.131704C.875118-4.07432 .932503-4.002589 1.032926-4.002589C1.147696-4.002589 1.176388-4.088666 1.20508-4.146051C1.692851-5.738478 2.654046-6.039748 3.070086-6.039748C3.715665-6.039748 3.84478-5.437208 3.84478-5.092899C3.84478-4.777283 3.758703-4.447321 3.586549-3.758703L3.098778-1.793274C2.883585-.932503 2.467546-.143462 1.707197-.143462C1.635466-.143462 1.276811-.143462 .975541-.329963C1.492004-.430386 1.606774-.860772 1.606774-1.032926C1.606774-1.31985 1.391581-1.492004 1.119003-1.492004C.774695-1.492004 .401693-1.190734 .401693-.731656C.401693-.129116 1.075965 .143462 1.692851 .143462C2.381468 .143462 2.869239-.401693 3.170509-.989888C3.400048-.143462 4.117358 .143462 4.648167 .143462C6.111479 .143462 6.886174-1.73589 6.886174-2.051506C6.886174-2.123237 6.828789-2.180622 6.742712-2.180622C6.613596-2.180622 6.59925-2.108891 6.556211-1.994121C6.168864-.731656 5.336785-.143462 4.691206-.143462C4.189089-.143462 3.916511-.516463 3.916511-1.104657C3.916511-1.420273 3.973896-1.649813 4.203435-2.596661L4.705552-4.547744C4.920745-5.408516 5.408516-6.039748 6.068441-6.039748C6.097133-6.039748 6.498827-6.039748 6.800097-5.853248Z'/>
    <path id='g1-121' d='M3.773049 1.606774C3.385702 2.151929 2.826201 2.6397 2.123237 2.6397C1.951083 2.6397 1.262465 2.611008 1.047272 1.951083C1.090311 1.965429 1.162042 1.965429 1.190734 1.965429C1.62112 1.965429 1.908044 1.592428 1.908044 1.262465S1.635466 .817733 1.420273 .817733C1.190734 .817733 .688617 .989888 .688617 1.692851C.688617 2.424507 1.305504 2.926624 2.123237 2.926624C3.557857 2.926624 5.006822 1.606774 5.408516 .014346L6.814443-5.58067C6.828789-5.652401 6.857482-5.738478 6.857482-5.824555C6.857482-6.039748 6.685327-6.18321 6.470134-6.18321C6.341018-6.18321 6.039748-6.125826 5.924979-5.69544L4.86336-1.477658C4.791629-1.219427 4.791629-1.190734 4.67686-1.032926C4.389936-.631233 3.916511-.143462 3.227894-.143462C2.424507-.143462 2.352776-.932503 2.352776-1.31985C2.352776-2.137583 2.740123-3.24224 3.127471-4.275166C3.285279-4.691206 3.371356-4.892053 3.371356-5.178977C3.371356-5.781517 2.94097-6.326672 2.238007-6.326672C.918157-6.326672 .387347-4.246474 .387347-4.131704C.387347-4.07432 .444732-4.002589 .545155-4.002589C.674271-4.002589 .688617-4.059973 .746002-4.26082C1.090311-5.465901 1.635466-6.039748 2.194968-6.039748C2.324084-6.039748 2.567969-6.039748 2.567969-5.566324C2.567969-5.193323 2.410161-4.777283 2.194968-4.232128C1.492004-2.352776 1.492004-1.879352 1.492004-1.535043C1.492004-.172154 2.467546 .143462 3.184855 .143462C3.600895 .143462 4.117358 .014346 4.619475-.516463L4.633821-.502117C4.418628 .344309 4.275166 .90381 3.773049 1.606774Z'/>
    </defs>
    <g id='page1'>
    <use x='-72.000004' y='-65.623905' xlink:href='#g1-120'/>
    <use x='-60.829507' y='-65.623905' xlink:href='#g0-1'/>
    <use x='-53.656444' y='-65.623905' xlink:href='#g1-121'/>
    </g>
    </svg>"""
end
# =#

Latexify._writetex(xdoty_tex; name=name)
tex = open("$(name).tex") do f
    read(f, String)
end
@test tex == replace(raw"""
\documentclass[varwidth=100cm]{standalone}
\usepackage{amssymb}
\usepackage{amsmath}

\begin{document}
{
    \Large
    $x \cdot y$
}
\end{document}
""", "\r\n"=>"\n")

Latexify._writetex(L"\ce{ 2 P_1 &<=>[k_{+}][k_{-}] D_{1}}"; name=name)
tex = open("$(name).tex") do f
    read(f, String)
end
@test tex == replace(raw"""
\documentclass[varwidth=100cm]{standalone}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{mhchem}
\begin{document}
{
    \Large
    $\ce{ 2 P_1 &<=>[k_{+}][k_{-}] D_{1}}$
}
\end{document}
""", "\r\n"=>"\n")


text:
# Abstract type definition for Phylogenetic trees.
abstract Phylogeny

# The type definition for a Phylogenetic Tree with branch lengths.
immutable Phylo <: Phylogeny
	name::String
	edge::Array{Int,2}
	Nnode::Int
	tipLabel::Array{String}
	edgeLength::Array{Float64}
	nodeLabel::Array{String}
	rootEdge::Float64
	Phylo(name, edge,
              Nnode, tipLabel,
              edgeLength, nodeLabel,
              rootEdge) = new(name, edge,
                              Nnode, tipLabel,
                              edgeLength, nodeLabel, rootEdge)
end

# The type definition for a Phylogenetic Tree without branch lengths.
immutable Clado <: Phylogeny
	name::String
	edge::Array{Int,2}
	Nnode::Int
	tipLabel::Array{String}
	nodeLabel::Array{String}
	Clado(name,
              edge,
              Nnode,
              tipLabel,
              nodeLabel) = new(name, edge, Nnode, tipLabel, nodeLabel)
end

# Type definition for a small simple representation of a tree.
immutable ReducedTopology <: Phylogeny
	name::String
	indiesArray::Array{Int}

	function ReducedTopology(phy::Phylogeny)
		children = phy.edge[1:size(phy.edge,1), 2]
		parents = phy.edge[1:size(phy.edge,1), 1]
		IArray = [length(findin(children, i)) != 0 ? parents[findin(children, i)[1]] : -1 for i in 1:max(phy.edge)]
		new(phy.name, IArray)
	end
end

function ReducedTopology(phy::Array{Phylogeny})
	outarray = Array(ReducedTopology, length(phy))
	for i in 1:length(phy)
		outarray[i] = ReducedTopology(phy[i])
	end
	return outarray
end

# equality, etc.

function =={T<:Phylogeny}(p1::T,p2::T)
    for field in names(T)
        if getfield(p1,field) != getfield(p2,field)
            return false
        end
    end
    return true
end

# Julia 0.2 compatibility
function isequal{T<:Phylogeny}(p1::T,p2::T)
    for field in names(T)
        if !isequal(getfield(p1,field),getfield(p2,field))
            return false
        end
    end
    return true
end

function hash{T<:Phylogeny}(p::T)
    h = 0
    for field in names(T)
        h = bitmix(hash(getfield(p,field)),h)
    end
    return h
end


text:
module LCMImageSequence

using
  LCMCore,
  Images,
  ImageView

export
  testfnc

function testfnc()
  warn("Initial commit")
end

end


text:
module Jenks

export JenksClassification, GVF, JenksResult

using Random, StatsBase, Statistics, Printf

include("initialisation.jl")
include("metrics.jl")
include("JenksStruct.jl")
include("JenksOptimization.jl")

end


text:
# CodeTable16
# ===========

const CodeTable16 = CodeTable{16}

const BASE16_CODEIGN = 0x10  # IGNore
const BASE16_CODEEND = 0x11  # END
const BASE16_CODEERR = 0xff  # ERRor

ignorecode(::Type{CodeTable16}) = BASE16_CODEIGN

"""
    CodeTable16(asciicode:::String; casesensitive::Bool=false)

Create a code table for base16.
"""
function CodeTable16(asciicode::String; casesensitive::Bool=false)
    if !isascii(asciicode)
        throw(ArgumentError("the code table must be ASCII"))
    elseif length(asciicode) != 16
        throw(ArgumentError("the code size must be 16"))
    end
    encodeword = Vector{UInt8}(undef, 16)
    decodeword = Vector{UInt8}(undef, 256)
    fill!(decodeword, BASE16_CODEERR)
    for (i, char) in enumerate(asciicode)
        bits = UInt8(i-1)
        code = UInt8(char)
        encodeword[bits+1] = code
        decodeword[code+1] = bits
        if !casesensitive
            if isuppercase(char)
                code = UInt8(lowercase(char))
                decodeword[code+1] = bits
            end
            if islowercase(char)
                code = UInt8(uppercase(char))
                decodeword[code+1] = bits
            end
        end
    end
    # NOTE: The padcode is not used.
    return CodeTable16(encodeword, decodeword, 0x00)
end

@inline function encode(table::CodeTable16, x::UInt8)
    return table.encodeword[Int(x)+1]
end

@inline function decode(table::CodeTable16, x::UInt8)
    return table.decodeword[Int(x)+1]
end

"""
The hexadecimal base16 code table (encoding: uppercase; decoding: case-insensitive).
"""
const BASE16_UPPER = CodeTable16("0123456789ABCDEF")

"""
The hexadecimal base16 code table (encoding: lowercase; decoding: case-insensitive).
"""
const BASE16_LOWER = CodeTable16("0123456789abcdef")


text:
using BinDeps

@BinDeps.setup

libbson = library_dependency(
    "libbson",
    aliases = ["libbson", "libbson-1.0"]
    )

provides(Sources, Dict(
    URI("http://github.com/mongodb/libbson/releases/download/1.0.0/libbson-1.0.0.tar.gz") => libbson
    ))

provides(
    BuildProcess,
    Autotools(libtarget = "libbson-1.0.la"),
    libbson,
    os = :Unix
    )

@osx_only begin
    using Homebrew
    Homebrew.rm("libbson")
    provides(Homebrew.HB, "mongo-c", libbson, os = :Darwin)
end

@BinDeps.install Dict(:libbson => :libbson)


text:
export @profile

"""
    @profile

`@profile <expression>` runs your expression, while activating the CUDA profiler upon first
kernel launch.
"""
macro profile(ex)
    quote
        Profile.enabled[] = true
        try
            $(esc(ex))
        finally
            Profile.started[] && CUDAdrv.Profile.stop()
            Profile.enabled[] = false
            Profile.started[] = false
        end
    end
end

module Profile

using CUDAdrv

const enabled = Ref(false)
const started = Ref(false)

macro launch(ex)
    quote
        if enabled[] && !started[]
            started[] = true
            CUDAdrv.Profile.start()
        end

        $(esc(ex))
    end
end

end


text:
using Luna

a = 13e-6
gas = :H2
pres = 5
flength = 200e-2

τfwhm = 20e-15
λ0 = 800e-9
energy = 1e-6

grid = Grid.EnvGrid(flength, λ0, (160e-9, 3000e-9), 40e-12)

m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff = let m=m
    z -> Modes.Aeff(m, z=z)
end

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),
             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, gas)))

linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas)
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.time_1D(output, range(0,flength,length=5), trange=(-0.5e-12, 2e-12))
Plotting.spec_1D(output, range(0,flength,length=5), λrange=(500e-9, 1400e-9))


text:
"""

Create a trankplot.

$(SIGNATURES)

## Positional arguments
```julia
* `df::DataFrame` # DataFrame holding the draws
* `param::AbstractString` # Symbol to use
```

## Keyword arguments
```julia
* `bins = 40` # Number of bins
* `n_draws = 1000` # Number of draws in df
* `n_chains = 4` # Number of chains appended in df
* `n_eff = 0` # Effective number of samples
```

"""
function trankplot(df::DataFrame, param::AbstractString;
        bins=40, n_draws=1000, n_chains=4, n_eff=0, kwargs...)

    nt_args = (n_draws=n_draws, n_chains=n_chains)
    ranks = rank_vector(df[:, param], nt_args)

    f = Figure()
    if n_eff > 0
        ax = Axis(f[1, 1];
            title="Trankplot of parameter $(param) (n_eff = $(n_eff))")
    else
        ax = Axis(f[1, 1])
    end

    Makie.xlims!(ax, 1, 41)
    colors = [:black, :green, :blue, :red]

    for j in 1:4
        bv = bin_vector(ranks[:, j], 40)
        s = [Meta.parse(string(split(string(cut(bv, 40)[i])[2:4], ":")[1])) for i in 1:40]
        for i in 1:40
            lines!([Float64(i), Float64(i+1)], [s[i], s[i]]; color=colors[j])
            if i == 1
                lines!([Float64(i), Float64(i)], [20, s[i]]; color=colors[j])
                lines!([Float64(i+1), Float64(i+1)], [s[i], s[i+1]]; color=colors[j])
            elseif i < 40
                lines!([Float64(i+1), Float64(i+1)], [s[i], s[i+1]]; color=colors[j])
            elseif i == 40
                lines!([Float64(i+1), Float64(i+1)], [s[i], 20]; color=colors[j])
            end
        end
    end
    f
end

export
    trankplot


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export call_request, cast_request

abstract Request
abstract Response

#= Requests

Straightforward request types are defined here. The surprising aspect here is
that requests do not contain a module. Instead, modules are resolved as soon
as the request is made given a symbol which represents that module. As an
example, an Erlang node may make a call request (:FsExtern, :test_f, [1024]).
Here, the call_request constructor will be invoked in order to lookup FsExtern
module and use this to find the function test_f which will then become part
of the constructed request type.

=#

type CallRequest <: Request
    fun::Function
    args::Tuple
end

type StreamingCallRequest <: Request
    fun::Function
    args::Tuple

    StreamingCallRequest(r::CallRequest) = begin
        new(r.fun, r.args)
    end
end

type CastRequest <: Request
    fun::Function
    args::Tuple
end

type InfoRequest <: Request
    command::Symbol
    options::Array{Any}
end

type ErrorRequest <: Request
    err::Any
end

type CallbackCastRequest <: Request
    fun::Function
    args::Tuple
    callback_mod::Symbol
    callback_fun::Symbol

    CallbackCastRequest(r1::CastRequest, r2::InfoRequest) = begin
        new(r1.fun, r1.args, r2.options[:mod], r2.options[:fun])
    end
end

#= External Constructors

These were created in order to ensure that the inputs to the later constructed
types are correctly formed. In the case of call / cast requests for example,
symbols are mapped to module / function objects.

=#

function call_request (mod_sym::Symbol, fun_sym::Symbol, args::Tuple)
    mod::Module = eval(mod_sym)
    fun::Function = eval(mod, fun_sym)
    CallRequest(fun, args)
end

function cast_request (mod_sym::Symbol, fun_sym::Symbol, args::Tuple)
    mod::Module = eval(mod_sym)
    fun::Function = eval(mod, fun_sym)
    CastRequest(fun, args)
end

max_issues_repo_name: esl/FsServer

max_forks_repo_name: esl/FsServer


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

name = "Dbus"
version = v"1.12.16"

# Collection of sources required to build Dbus
sources = [
    ArchiveSource("https://dbus.freedesktop.org/releases/dbus/dbus-$(version).tar.gz",
                  "54a22d2fa42f2eb2a871f32811c6005b531b9613b1b93a0d269b05e7549fec80"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/dbus-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-xml=expat \
    --with-dbus-user=messagebus \
    --with-system-pid-file=/var/run/dbus.pid \
    --disable-verbose-mode \
    --disable-static \
    --enable-inotify \
    --disable-dnotify \
    --disable-asserts \
    --enable-user-session \
    --with-session-socket-dir=/tmp \
    --with-x
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

products = [
    LibraryProduct("libdbus-1", :libdbus),
    ExecutableProduct("dbus-daemon", :dbus_daemon)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Expat_jll", v"2.2.7"; compat="~2.2.7"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
const JULIACOMPUTING_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/Deprecations.jl",
        destination_repo_name = "JuliaComputing-Deprecations.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/FemtoCleaner.jl",
        destination_repo_name = "JuliaComputing-FemtoCleaner.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/NewPkgEval.jl",
        destination_repo_name =
            "NewPkgEval.jl-9f2e2246-6dce-11e8-3d98-4b291446da6e",
        ),
    ]


text:
module seasonality_exponential_smoothingTests

using Test

using TSeriesForecast: HW, loss, forecast

ϵ = 0.1

### Data

data = [42.2056638559382, 24.6491713296294, 32.6673351386246, 37.2573540134764, 45.2424602720635, 29.3504812742718, 36.3442072775119, 41.7820813647366, 49.2765984345619, 31.2754013898459, 37.8506254947152, 38.83704413019, 51.236900336137, 31.8385516171558, 41.3234212576972, 42.7990033724349, 55.708358364385, 33.4071449171312, 42.3166379718116, 45.1571225674731, 59.5760799633937, 34.8373301595328, 44.8416807241332, 46.971249602779, 60.0190309438716, 38.3711785132276, 46.9758641283558, 50.7337964560878, 61.6468731860574, 39.2995693656296, 52.6712090814052, 54.3323168946667, 66.834358381339, 40.8711884667851, 51.8285357927739, 57.4919099342262, 65.2514698518726, 43.0612082202828, 54.7607571288007, 59.8344749355003, 73.2570274672009, 47.6966237298, 61.0977680206996, 66.0557612187001 ]
α = 0.306342956680807
β = 0.000100001296624933
γ = 0.426290711530531
l0 = 32.2596735425297
b0 = 0.701381297812788
h = 8
m = 4
s0 = [9.6962, -9.31324086158846, -1.6935401189679, 1.31060178041439]

expected_mse = 3.10924508566665 #mse = loss/length(time_serie)
y_expected_fitted = [ 42.6572340404844, 24.21081477348, 32.6661829360246, 36.3720578962127, 45.5378079931788, 27.5187168280733, 36.2148118182051, 40.3371306494425, 49.171340569086, 32.1820914045507, 39.3136533640888, 43.5088715586704, 49.8974573204696, 32.8543378844383, 39.7147361880909, 43.4827966071507, 53.6557621665274, 35.8272959907684, 43.3763205153963, 45.3494719969994, 56.8393233421518, 37.3136607326191, 45.4253309990936, 47.9139961056344, 60.4224926964744, 37.7121150317104, 47.5909961330668, 49.9168074175707, 63.1943244235032, 40.5863696216504, 49.3257081290391, 53.4754641480289, 65.7573910045542, 44.0649207101445, 54.1949995740771, 55.5336851714369, 68.2465899064367, 43.4854450109548, 54.8162133971451, 58.7062827647439, 69.0531133095763, 47.5937707433477, 59.24375503884, 64.2240841042244 ]
y_expected_forecast = [76.0983732770355, 51.6033257926593, 63.9686736144129, 68.3717004584965, 78.9040413079506, 54.4089938235743, 66.7743416453279, 71.1773684894116 ]

model = HW(m, α, β, γ, l0, b0, s0)

### Tests

@testset "Seasonality Exponential Smoothing" begin

    @testset "HW Seasonal struct" begin

        @testset "HW default constructor" begin
            model = HW(4)
            @test model.α ≈ 0
            @test model.β ≈ 0
            @test model.γ ≈ 0
            @test model.l0 ≈ 0
            @test model.b0 ≈ 0
            @test model.m ≈ 4
            @test model.s0 ≈ [0, 0, 0, 0]
        end

        @testset "HW parametric constructor" begin
            α = 0.8
            β = 0.001
            γ = 0.42
            l0 = 42.0
            b0 = 2.3
            m = 4
            s0 = [-9.31, 9.32, -1.24, 1.26]
            model = HW(m, α, β, γ, l0, b0, s0)
            @test model.α ≈ α
            @test model.β ≈ β
            @test model.γ ≈ γ
            @test model.l0 ≈ l0
            @test model.b0 ≈ b0
            @test model.m ≈ m
            @test model.s0 ≈ s0
        end

        @testset "HW parametric constructor takes any Number type" begin
            α = Float32(0.8)
            β = Float16(0.001)
            γ = 0.42
            l0 = 42
            b0 = 2.3
            m = 4
            s0 = [-9.31, 9.32, -1.24, 1.26]
            model = HW(m, α, β, γ, l0, b0, s0)
            @test model.α ≈ α
            @test model.β ≈ β
            @test model.γ ≈ γ
            @test model.l0 ≈ l0
            @test model.b0 ≈ b0
            @test model.m ≈ m
            @test model.s0 ≈ s0
        end

        @testset "Length(S0) == to m" begin
            α = Float32(0.8)
            β = Float16(0.001)
            γ = 0.42
            l0 = 42
            b0 = 2.3
            m = 3
            s0 = [-9.31, 9.32, -1.24, 1.26]
            @test_throws ErrorException model = HW(m, α, β, γ, l0, b0, s0)
        end

    end

    @testset "loss" begin
        @test isapprox(loss(model, data)/length(data), expected_mse, atol=ϵ)
    end

    @testset "forecast" begin
        @test isapprox(forecast(model, data, h), y_expected_forecast, atol=ϵ)
        @test length(forecast(model, data, h)) == h
    end

end

end


text:
#=##############################################################################################
Filename: lift_distribution_template.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this is a template file. Convenience methods are provided to calculate the lift
        distribution of a single wing over the specified alphas.
=###############################################################################################

# initialize parameters
"""
    LiftDistribution <: Parameters

# Fields

* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step
* `cfs::Vector{Vector{Array{TF,2}}}`
* `cms::Vector{Vector{Array{TF,2}}}`
* `wake_function::V4`: additional velocity function; default is nothing
* `plot_directory::String`
* `plot_base_name::String`
* `plot_extension::String`
* `plotstepi::AbstractRange`
* `surfacenames::Vector{String}`

"""
struct LiftDistribution{TF,TR} <: Parameters
    CFs::Array{TF,2}
    CMs::Array{TF,2}
    cfs::Vector{Vector{Array{TF,2}}}
    cms::Vector{Vector{Array{TF,2}}}
    wake_function::Vector{Function}
    plot_directory::String
    plot_base_name::String
    plot_extension::String
    plotstepi::TR
    surfacenames::Vector{String}
end

LiftDistribution(CFs, CMs, cfs, cms, plot_directory, plot_base_name, plot_extension, ploti, surfacenames; wake_function=[nothing]) =
    LiftDistribution(CFs, CMs, cfs, cms, wake_function, plot_directory, plot_base_name, plot_extension, ploti, surfacenames)

# function (rs::RotorSweepParameters{V1, V2, V3, V4, V5})(alphas, omega...)

# end

# rs  = RotorSweepParameter(args...)

# rs(alphas, omega, nblades...)

"""
    lift_distribution_template(ploti, alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip;
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        surfacenames = ["default wing"],
        kwargs...)

# Arguments:

* `vinfs::Vector{Float64}`: freestream velocity for each time step
* `ploti`
* `alphas::Vector{Float64}`
* `wing_b::Float64`
* `wing_TR::Float64`
* `wing_AR::Float64`
* `wing_θroot::Float64`
* `wing_θtip::Float64`

# Keyword Arguments:

* `plot_directory = joinpath(topdirectory, "data","plots")`: directory to which plots are saved
* `plot_base_name = "default"`: first part of saved plot filenames
* `plot_extension = ".pdf"`: saved plots file extension
* `step_symbol`
* `surfacenames`

"""
function lift_distribution_template(ploti, alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ;
            plot_directory=joinpath(topdirectory,"data","plots",TODAY),
            plot_base_name="default",
            plot_extension=".pdf",
            step_symbol=L"\alpha ",
            surfacenames=["default wing"],
            kwargs...)

    # prepare subsystems
    wings = simplewing_system(wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ; kwargs...)
    rotors = CCBladeSystem(
        Vector{CC.Rotor}(undef,0),
        Vector{Vector{CC.Section{Float64, Float64, Float64, nothing}}}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Int64}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Bool}(undef,0)
    )
    inertia = nothing
    nonliftingbodies = nothing
    structures = nothing
    motors = nothing
    batteries = nothing

    # build system struct
    aircraft = Aircraft(wings, rotors, inertia, nonliftingbodies, structures, motors, batteries)

    # compile actions
    actions = [solve_wing_CF_CM, lift_moment_distribution]

    # initialize parameters
    wake_function, CFs, CMs = solve_wing_CF_CM(aircraft, alphas) # wake_function, CFs, CMs
    cfs, cms = lift_moment_distribution(aircraft, alphas) # let step_range be replaced by alphas
    _, _, _, _, _, _, _ = post_plot_lift_moment_distribution(aircraft, alphas)

    # check sizes and instantiate struct
    @assert size(CFs) == size(CMs) "size of CFs and CMs inconsistent"
    @assert size(CFs)[2] == length(alphas) "length of parameter CFs and alphas inconsistent"
    @assert size(CFs)[1] == 3 "first dimension of CFs should be of length 1"
    nspanwisepanels1 = size(aircraft.wing_system.system.surfaces[1])[2]
    @assert length(cfs[1]) == length(surfacenames) "length of cfs is inconsistent: expected $(length(surfacenames)); got $(length(cfs))"
    @assert length(cms[1]) == length(surfacenames) "length of cms is inconsistent: expected $(length(surfacenames)); got $(length(cms))"
    @assert size(cfs[1][1]) == (3, nspanwisepanels1) "size of cfs[1] is inconsistent: expected $((3, nspanwisepanels1)); got $(size(cfs[1]))"

    # prepare plot directory
    if !isdir(plot_directory); mkpath(plot_directory); end

    # build parameters struct
    parameters = LiftDistribution(CFs, CMs, cfs, cms, wake_function, plot_directory, plot_base_name, plot_extension, ploti, surfacenames)

    # build freestream_function
    function freestream_function(aircraft, parameters, environment, alphas, stepi)

        # calculate freestream
        vinf = 1.0 # + ti # arbitrary for lift distribution?
        alpha = alphas[stepi]
        beta = 0.0
        Omega = StaticArrays.@SVector zeros(3)
        freestream = Freestream(vinf, alpha, beta, Omega)

        return freestream
    end

    # build environment_function
    function environment_function(aircraft, parameters, alphas, stepi)
        Environment() # arbitrary
    end

    # compile postactions
    postactions = [post_plot_cf_cm_alpha_sweep, post_plot_lift_moment_distribution]

    # build objective_function
    objective_function(aircraft, parameters, freestream, environment, alphas) = 0.0

    return aircraft, parameters, actions, freestream_function, environment_function, postactions, objective_function, alphas, step_symbol
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Statistics, StatsBase, Random, LinearAlgebra

function ARMAmodel(θ, rndseed=1234)
    Random.seed!(rndseed)
    N = 300
    burnin = 50
    α = θ[1]
    β = θ[2]
    σ = sqrt(θ[3])
    data = zeros(N)
    xlag = 0.0
    flag = 0.0
    @inbounds for i = 1:(N+burnin)
        f = σ*randn()
        x = α*xlag + f - β*flag
        xlag = x
        flag = f
        if i > burnin
            data[i-burnin] = x
        end
    end
    data
end

function auxstat(θ, reps)
    data = [ARMAmodel(θ, rand(1:Int64(1e12))) for i = 1:reps]  # reps draws of data
    auxstat.(data)
end

@views function auxstat(x)
    s = std(x)
    c = cor(x[2:end,:],x[1:end-1])
    n = size(x,1)
    b, varb, u, junk1, rsq  = ols(x[2:end,:], x[1:end-1,:])
    ϕ, varϕ, u, junk2, rsq2 = ols(u[2:end,:], u[1:end-1,:])
    σ = std(u)
    sqrt(300.0).*vcat(s, c, b, varb, rsq, ϕ, varϕ, rsq2, σ, pacf(x,collect(1:4)))
end

function TrueParameters()
    [0.95, 0.5, 1.0]
end

function PriorSupport()
    lb = [-0.99, -0.99, 0.0001]
    ub = [0.99, 0.99, 4.0]
    lb,ub
end

function InSupport(θ)
    lb,ub = PriorSupport()
    all(θ .>= lb) & all(θ .<= ub)
end

function PriorMean()
    lb, ub = PriorSupport()
    (ub - lb) ./ 2.0
end

function PriorDraw()
    lb, ub = PriorSupport()
    (ub - lb) .* rand(3) + lb
end

# prior just checks that we're in the bounds
function Prior(theta)
    lb, ub = PriorSupport()
    a = 0.0
    if(all((theta .>= lb) .& (theta .<= ub)))
        a = 1.0
    end
    return a
end

# from Econometrics
# compute ols coefficients, fitted values, and errors
function lsfit(Y, X)
    beta = X\Y
    fit = X*beta
    errors = Y - fit
    return beta, fit, errors
end

function ols(y, x)
    n,k = size(x)
    b, fit, e = lsfit(y,x)
    df = n-k
    sigsq = (e'*e/df)[1,1]
    xx_inv = inv(x'*x)
    ess = (e' * e)[1,1]
    xe = x.*e
    varb = xx_inv*xe'xe*xx_inv
    seb = sqrt.(diag(varb))
    seb = seb.*(seb.>1e-16) # round off to zero when there are restrictions
    t = b ./ seb
    tss = y .- mean(y)
    tss = (tss'*tss)[1,1]
    rsq = (1.0 - ess / tss)
    return b, varb, e, ess, rsq
end

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
import DataFrames
import Statistics

"""
```julia
chill(::Val{:MeanTemperature}, weatherData::DataFrame)::DataFrame
```

Implements the [Mean Temperature](chill_MeanTemperature.md) model.

# Arguments

- `weatherData`: A `DataFrames.DataFrame` containing the weather data used for
    estimating the chill units. `weatherData` must be formatted to match the
    output of [`BiTWeather.read`](@ref) and contain at least the `:temperature` field.

# Return

- A `DataFrames.DataFrame` containing the chill units indexed by
    year/month/day/hour. The returned columns are: x and y.

For temperate climates, you want to use data from the coldest month of the year
for the calculation. Whereas, for cold climates, you want to use data from the
coldest two months of the year.

# Example

```julia
import DataFrames
import Dates
import BiTWeather

configuration = BiTWeather.Configuration(
    dsn = "meteobridge",
    table = "backinthirty",
    fieldMappings = Dict(
        :dateTime =>    BiTWeather.FieldMapping(name = "DateTime",   units = nothing),
        :temperature => BiTWeather.FieldMapping(name = "TempOutNow", units = "F")
    )
)

year::Int = 2021
rangeStart::Dates.Date = Dates.Date(year, 1, 1)
rangeEnd::Dates.Date = Dates.Date(year, 1, 31)
range::Tuple{Float64, Float64} = (rangeStart, rangeEnd)
fields::Vector{Symbol} = [:temperature]
weatherData::DataFrames.DataFrame = BiTWeather.read(configuration, range, fields)
chillData::DataFrames.DataFrame = BiTWeather.chill(Val(:MeanTemperature), weatherData)
println(chillData)
```
"""
function chill(::Val{:MeanTemperature}, weatherData::DataFrames.DataFrame)::DataFrames.DataFrame

    # The chill calculation requires the :temperature field.
    if !(:temperature in DataFrames.propertynames(weatherData))
        error("weatherData must contain the :temperature field")
    end

    mt::DataFrames.DataFrame = weatherData

    # Determine daily max and min temperatures.
    gdf::DataFrames.GroupedDataFrame = DataFrames.groupby(mt, [:year, :month, :day])
    mt = DataFrames.combine(gdf,
        [:temperature] .=> [maximum, minimum] .=> [:x_max_d, :x_min_d]
    )

    # Determine mean max and mean min temperatures.
    mt = DataFrames.combine(mt,
        [:x_max_d, :x_min_d] .=> Statistics.mean .=> [:x_max, :x_min]
    )

    # Determine mean temperature.
    mt[!, :x] = (mt.x_max .+ mt.x_min) ./ 2.0

    # Determine chill units.
    function y(x::Float64)::Float64
        return (-100.0 * x) + (5600.0 / 3.0)
    end
    mt[!, :y] = y.(mt.x)

    DataFrames.select!(mt, [:x, :y])

    for name in DataFrames.propertynames(mt)
        if eltype(mt[!, name]) == Float64
            mt[!, name] = round.(mt[!, name], digits = 2)
        end
    end

    return mt
end

"""
```julia
chillPlot(::Val{:MeanTemperature}, chillData::DataFrame)::PlotlyJS.Plot
```

Generates a plot of the [Cumulative Chill Unit](chill_CumulativeChillUnit.md) model
data.

# Arguments

- `chillData::DataFrames.DataFrame`: The chill data used in generating the plot.
    `chillData` must be formated to match the output of
    [`BiTWeather.chill`](@ref) called with `Val(:MeanTemperature)`.

# Return

- A `PlotlyJS.Plot` containing a plot of the temperature and chill portion data.

# Example

```julia
import DataFrames
import Dates
import PlotlyJS
import BiTWeather

configuration = BiTWeather.Configuration(
    dsn = "meteobridge",
    table = "backinthirty",
    fieldMappings = Dict(
        :dateTime =>    BiTWeather.FieldMapping(name = "DateTime",   units = nothing),
        :temperature => BiTWeather.FieldMapping(name = "TempOutNow", units = "F")
    )
)

year::Int = 2021
rangeStart::Dates.Date = Dates.Date(year, 1, 1)
rangeEnd::Dates.Date = Dates.Date(year, 1, 31)
range::Tuple{Float64, Float64} = (rangeStart, rangeEnd)
fields::Vector{Symbol} = [:temperature]
weatherData::DataFrames.DataFrame = BiTWeather.read(configuration, range, fields)
chillData::DataFrames.DataFrame = BiTWeather.chill(Val(:MeanTemperature), weatherData)
plot::PlotlyJS.Plot = BiTWeather.chillPlot(Val(:MeanTemperature), chillData)
PlotlyJS.display(plot)
```
"""
function chillPlot(::Val{:MeanTemperature}, chillData::DataFrames.DataFrame)::PlotlyJS.Plot
    return _chillPlot(chillData, :x, :y, "Mean Temperature", "Units")
end

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
@init PyObject_TryConvert_AddRules("juliaaa.As", [
    (Any, PyAs_ConvertRule_tryconvert),
])

PyAs_ConvertRule_tryconvert(o, ::Type{S}) where {S} = begin
    # get the type
    to = PyObject_GetAttrString(o, "type")
    isnull(to) && return -1
    err = PyObject_Convert(to, Type)
    Py_DecRef(to)
    ism1(err) && return -1
    t = takeresult(Type)
    # get the value
    vo = PyObject_GetAttrString(o, "value")
    isnull(vo) && return -1
    err = PyObject_Convert(vo, t)
    Py_DecRef(vo)
    ism1(err) && return -1
    v = takeresult(t)
    # convert
    putresult(tryconvert(S, v))
end


text:
print("Enter an upper bound: ")
lower = 0
input = readline()
upper = parse(Int, input)

if upper < 1
    throw(DomainError)
end

attempts = 1
print("Think of a number, ", lower, "--", upper, ", then press ENTER.")
readline()
const maxattempts = round(Int, ceil(-log(1 / (upper - lower)) / log(2)))
println("I will need at most ", maxattempts, " attempts ",
    "(⌈-log(1 / (", upper, " - ", lower, ")) / log(2)⌉ = ",
    maxattempts, ").\n")
previous = -1
guess = -1

while true
    previous = guess
    guess = lower + round(Int, (upper - lower) / 2, RoundNearestTiesUp)

    if guess == previous || attempts > maxattempts
        println("\nThis is impossible; did you forget your number?")
        exit()
    end

    print("I guess ", guess, ".\n[l]ower, [h]igher, or [c]orrect? ")
    input = chomp(readline())

    while input ∉ ["c", "l", "h"]
        print("Please enter one of \"c\", \"l\", or \"h\". ")
        input = chomp(readline())
    end

    if input == "l"
        upper = guess
    elseif input == "h"
        lower = guess
    else
        break
    end

    attempts += 1
end

println("\nI win after ", attempts, attempts == 1 ? " attempt." : " attempts.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#======================================
Methods using semidefinite programming
======================================#
using .SumOfSquares

"""
    new_sos(backend; kwargs...)

Return a new (empty) sum-of-squares optimization problem for the given backend.

### Input

- `backend` -- the backend (also called JuMP `solver`)
- `kwargs`  -- additional keyword arguments

### Output

An instance of `SOSModel` for the given backend and options.
"""
@inline function new_sos(backend, kwargs...)
    𝑂 = Dict(kwargs)

    if :QUIET ∈ keys(𝑂)
        # for mosek solver
        SOSModel(backend, QUIET=𝑂[:QUIET])
    else
        SOSModel(backend)
    end
end

function _enclose(sose::SumOfSquaresEnclosure, f::Function, dom::Interval_or_IntervalBox;
                  kwargs...)

    _enclose_SumOfSquares(f, dom, sose.order, sose.backend; kwargs...)
end

function _enclose_SumOfSquares(f::Function, dom::Interval, order::Int,
                               backend; kwargs...)

    # polynomial variables
    @polyvar x
    p = f(x)

    # box constraints
    B = @set inf(dom) <= x && x <= sup(dom)

    # ============
    # Upper bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p <= γ, domain=B, maxdegree=order)
    @objective(model, Min, γ)
    optimize!(model)
    upper_bound = objective_value(model)

    # ============
    # Lower bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p >= γ, domain=B, maxdegree=order)
    @objective(model, Max, γ)
    optimize!(model)
    lower_bound = objective_value(model)

    return Interval(lower_bound, upper_bound)
end

function _enclose_SumOfSquares(f::Function, dom::IntervalBox{N}, order::Int,
                               backend; kwargs...) where {N}

    # polynomial variables
    @polyvar x[1:N]
    p = f(x...)

    # box constraints
    Bi =[@set inf(dom[i]) <= x[i] && x[i] <= sup(dom[i]) for i in 1:N]
    B = reduce(intersect, Bi)

    # ============
    # Upper bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p <= γ, domain=B, maxdegree=order)
    @objective(model, Min, γ)
    optimize!(model)
    upper_bound = objective_value(model)

    # ============
    # Lower bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p >= γ, domain=B, maxdegree=order)
    @objective(model, Max, γ)
    optimize!(model)
    lower_bound = objective_value(model)

    return Interval(lower_bound, upper_bound)
end


text:
include("euler/euler.jl")
using .Numbers: get_digits
using Combinatorics: with_replacement_combinations

function compute(n::Int)::Int
    max_len = 2
    result = 0
    powers = Dict(i => i ^ n for i in 0:9)
    while length(string(max_len * 9 ^ n)) >= max_len
        for digits ∈ with_replacement_combinations(0:9, max_len)
            candidate = sum(map(x -> powers[x], digits))
            if sort(get_digits(candidate)) == sort(digits)
                result += candidate
            end
        end
        max_len += 1
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Test
# using Revise
using PolynomialBasis

PB = PolynomialBasis

function allequal(v1,v2)
    return all(v1 .≈ v2)
end

function allequal(v1,v2,tol)
    np = length(v1)
    f = length(v2) == np
    return f && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:np])
end

p = [-1.0 1.0]
@test_throws AssertionError PB.tensor_product_points(0,p)
@test_throws AssertionError PB.tensor_product_points(4,p)

tp = PB.tensor_product_points(1,p)
@test allequal(p,tp)

tp = PB.tensor_product_points(2,p)
testtp = [-1.0 -1.0  1.0  1.0
          -1.0  1.0 -1.0  1.0]
@test allequal(tp,testtp)

tp = PB.tensor_product_points(3,p)
testtp = [-1.0   -1.0   -1.0   -1.0    1.0    1.0    1.0    1.0
          -1.0   -1.0    1.0    1.0   -1.0   -1.0    1.0    1.0
          -1.0    1.0   -1.0    1.0   -1.0    1.0   -1.0    1.0]
@test allequal(tp,testtp)

p = [1.0]
tp = PB.tensor_product_points(1,p')
@test allequal(p,tp)

tp = PB.tensor_product_points(2,p')
testtp = [1.0,1.0]
@test allequal(tp,testtp)

tp = PB.tensor_product_points(3,p')
testtp = [1.0;1.0;1.0]
@test allequal(tp,testtp)

basis = PB.LagrangeBasis(0)
@test_throws AssertionError PB.LagrangeTensorProductBasis(0,0)
@test_throws AssertionError PB.LagrangeTensorProductBasis(4,0)
tpb = PB.LagrangeTensorProductBasis(1,0)
@test allequal(tpb(0.0),[1.0])
@test allequal(PB.gradient(tpb,0.0),[0.0])
@test allequal(PB.gradient(tpb,[0.0]),[0.0])
@test_throws AssertionError PB.gradient(tpb,[0.0,0.0])

tpb = PB.LagrangeTensorProductBasis(1,0)
@test allequal(tpb(0.0),[1.0])
@test allequal(PB.gradient(tpb,0.0),[0.0])
@test allequal(PB.gradient(tpb,[0.0]),[0.0])
@test_throws AssertionError PB.gradient(tpb,[0.0,0.0])

tpb = PB.LagrangeTensorProductBasis(1,1)
@test allequal(tpb(-1.0),[1.0,0.0])
@test allequal(tpb(1.0),[0.0,1.0])
@test_throws AssertionError tpb([1.0,2.0])
@test allequal(tpb([1.0]),[0.0,1.0])
@test allequal(PB.gradient(tpb,-1.0),[-0.5,0.5])
@test allequal(PB.gradient(tpb,1.0),[-0.5,0.5])
@test_throws AssertionError PB.gradient(tpb,[1.0,2.0])
@test allequal(PB.gradient(tpb,[1.0]),[-0.5,0.5])

tpb = PB.LagrangeTensorProductBasis(1,2)
@test allequal(tpb(-1.0),[1.0,0.0,0.0])
@test allequal(tpb(0.0),[0.0,1.0,0.0])
@test allequal(tpb(1.0),[0.0,0.0,1.0])
@test_throws AssertionError tpb([1.0,2.0])
@test allequal(tpb([1.0]),[0.0,0.0,1.0])
@test allequal(PB.gradient(tpb,-1.0),[-1.5,2,-0.5])
@test allequal(PB.gradient(tpb,0.0),[-0.5,0,0.5])
@test allequal(PB.gradient(tpb,1.0),[0.5,-2,1.5])
@test_throws AssertionError PB.gradient(tpb,[1.0,2.0])
@test allequal(PB.gradient(tpb,[1.0]),[0.5,-2,1.5])

@test_throws AssertionError PB.LagrangeTensorProductBasis(1,-1)
@test_throws AssertionError PB.LagrangeTensorProductBasis(1,2,start=2)
@test_throws AssertionError PB.LagrangeTensorProductBasis(1,2,stop=-2)
tpb = PB.LagrangeTensorProductBasis(1,2,stop=0.0)
@test allequal(tpb.points,[-1.0 -0.5 0.0])
@test allequal(tpb(-1.0),[1.0,0.0,0.0])
@test allequal(tpb(-0.5),[0.0,1.0,0.0])
@test allequal(tpb(0.0),[0.0,0.0,1.0])
@test_throws AssertionError tpb([1.0,2.0])
@test allequal(tpb([0.0]),[0.0,0.0,1.0])
@test allequal(PB.gradient(tpb,-1.0),[-3.,4,-1])
@test allequal(PB.gradient(tpb,-0.5),[-1.0,0,1])
@test allequal(PB.gradient(tpb,0.0),[1.0,-4,3])
@test_throws AssertionError PB.gradient(tpb,[1.0,2.0])
@test allequal(PB.gradient(tpb,[0.0]),[1.0,-4,3])

function test_basis_on_points(basis::PB.LagrangeTensorProductBasis{2,T,N}) where {T,N}
    flag = true
    for i in 1:N
        vals = zeros(N)
        vals[i] = 1.0
        p = basis.points[:,i]
        flag = flag && basis(p[1],p[2]) ≈ vals
        flag = flag && basis(p...) ≈ vals
    end
    return flag
end

tp2 = PB.LagrangeTensorProductBasis(2,0)
@test allequal(tp2.points,[0.0,0.0])
@test test_basis_on_points(tp2)
@test_throws MethodError PB.gradient(tp2,1.5,1.0,2.0)
@test allequal(PB.gradient(tp2,1,0.0,0.0),[0.0])
@test allequal(PB.gradient(tp2,2,0.0,0.0),[0.0])
@test allequal(PB.gradient(tp2,2,[0.0,0.0]),[0.0])
@test allequal(PB.gradient(tp2,0.0,0.0),[0.0 0.0])
@test allequal(PB.gradient(tp2,[0.0,0.0]),[0.0 0.0])
@test_throws AssertionError PB.gradient(tp2,[0.0])
@test_throws AssertionError PB.gradient(tp2,[0.0,0.0,0.0])

v1 = [1.0,0.0,0.0]
v2 = [0.0,1.0,0.0]
v3 = [0.0,0.0,1.0]
d1 = [-1.5,2.0,-0.5]
d2 = [-0.5,0.0,0.5]
d3 = [0.5,-2.0,1.5]
tp2 = PB.LagrangeTensorProductBasis(2,2)
@test_throws AssertionError PB.gradient(tp2,0,-1,1)
@test_throws AssertionError PB.gradient(tp2,3,-1,1)
@test allequal(PB.gradient(tp2, 1, -1.0, -1.0), kron(d1,v1))
@test allequal(PB.gradient(tp2, 1, -1.0, +0.0), kron(d1,v2))
@test allequal(PB.gradient(tp2, 1, -1.0, +1.0), kron(d1,v3))
@test allequal(PB.gradient(tp2, 1, +0.0, -1.0), kron(d2,v1))
@test allequal(PB.gradient(tp2, 1, +0.0, +0.0), kron(d2,v2))
@test allequal(PB.gradient(tp2, 1, +0.0, +1.0), kron(d2,v3))
@test allequal(PB.gradient(tp2, 1, +1.0, -1.0), kron(d3,v1))
@test allequal(PB.gradient(tp2, 1, +1.0, +0.0), kron(d3,v2))
@test allequal(PB.gradient(tp2, 1, +1.0, +1.0), kron(d3,v3))

@test allequal(PB.gradient(tp2, 1, [+1.0, +1.0]), kron(d3,v3))

@test allequal(PB.gradient(tp2, 2, -1.0, -1.0), kron(v1,d1))
@test allequal(PB.gradient(tp2, 2, -1.0, +0.0), kron(v1,d2))
@test allequal(PB.gradient(tp2, 2, -1.0, +1.0), kron(v1,d3))
@test allequal(PB.gradient(tp2, 2, +0.0, -1.0), kron(v2,d1))
@test allequal(PB.gradient(tp2, 2, +0.0, +0.0), kron(v2,d2))
@test allequal(PB.gradient(tp2, 2, +0.0, +1.0), kron(v2,d3))
@test allequal(PB.gradient(tp2, 2, +1.0, -1.0), kron(v3,d1))
@test allequal(PB.gradient(tp2, 2, +1.0, +0.0), kron(v3,d2))
@test allequal(PB.gradient(tp2, 2, +1.0, +1.0), kron(v3,d3))

@test allequal(PB.gradient(tp2, 2, [+1.0, +1.0]), kron(v3,d3))

@test allequal(PB.gradient(tp2,1.0,1.0),hcat(kron(d3,v3),kron(v3,d3)))
@test_throws AssertionError PB.gradient(tp2,[1.0,1.0,1.0])
@test allequal(PB.gradient(tp2,[1.0,1.0]),hcat(kron(d3,v3),kron(v3,d3)))

function test_basis_on_points(basis::PB.LagrangeTensorProductBasis{3,N}) where {T,N}
    flag = true
    for i in 1:N
        vals = zeros(N)
        vals[i] = 1.0
        p = basis.points[:,i]
        flag = flag && basis(p[1],p[2],p[3]) ≈ vals
        flag = flag && basis(p) ≈ vals
    end
    return flag
end

tp3 = PB.LagrangeTensorProductBasis(3,2)
@test test_basis_on_points(tp3)

@test_throws AssertionError PB.gradient(tp3,0,-1.0,-1.0,-1.0)
@test_throws AssertionError PB.gradient(tp3,4,-1.0,-1.0,-1.0)
@test allequal(PB.gradient(tp3,1,-1.0,-1.0,-1.0),kron(d1,v1,v1))
@test allequal(PB.gradient(tp3,2,-1.0,0.0,-1.0),kron(v1,d2,v1))
@test allequal(PB.gradient(tp3,3,-1.0,0.0,+1.0),kron(v1,v2,d3))

@test_throws AssertionError PB.gradient(tp3,1,[1.0])
@test_throws AssertionError PB.gradient(tp3,1,[1.0,1.0])
@test_throws AssertionError PB.gradient(tp3,1,[1.0,1.0,1.0,1.0])

@test allequal(PB.gradient(tp3,3,[-1.0,-1.0,-1.0]),kron(v1,v1,d1))

@test allequal(PB.gradient(tp3,1.0,-1.0,0.0),hcat(kron(d3,v1,v2),kron(v3,d1,v2),kron(v3,v1,d2)))
@test_throws AssertionError PB.gradient(tp3,[1.0])
@test_throws AssertionError PB.gradient(tp3,[1.0,2.0,1.0,3.0])
@test allequal(PB.gradient(tp3,[1.0,-1.0,0.0]),hcat(kron(d3,v1,v2),kron(v3,d1,v2),kron(v3,v1,d2)))

tp1 = PB.LagrangeTensorProductBasis(1,2)
@test allequal(PB.hessian(tp1,-1),[1,-2,1])
@test allequal(PB.hessian(tp1,0),[1,-2,1])
@test allequal(PB.hessian(tp1,1),[1,-2,1])
@test allequal(PB.hessian(tp1,[-1]),[1,-2,1])
@test allequal(PB.hessian(tp1,[0]),[1,-2,1])
@test allequal(PB.hessian(tp1,[1]),[1,-2,1])

tp2 = PB.LagrangeTensorProductBasis(2,2)
f(x) = x[1]^2 + 2x[1]*x[2] + x[2]^2
coeffs = mapslices(f,Array(tp2.points),dims=1)
@test all([allequal(coeffs*PB.hessian(tp2,tp2.points[:,i]),[2,2,2]) for i = 1:9])

tp3 = PB.LagrangeTensorProductBasis(2,3)
f(x) = x[1]^3 + 2x[1]^2*x[2] + 18.0
tp3h(x) = [6x[1]+4x[2],4x[1],0.0]
coeffs = mapslices(f,Array(tp3.points),dims=1)
p = tp3.points
@test all([allequal(coeffs*PB.hessian(tp3,p[:,i]),tp3h(p[:,i]),1e3eps()) for i = 1:16])


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using FEMMaterials, Test, ForwardDiff

@testset "Calculate plastic strain" begin
    dummyel = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
    dummyip = first(get_integration_points(dummyel))
    # time = [...]
    # strain = [...]
    # stress = run_simulation(mat, time, strain)
    # @test isapprox(stress, stress_expected)

    n = 1.0
    K = 100.0

    mat = Material(ViscoPlastic, (:norton, [K, n]))
    mat.properties.youngs_modulus = 200.0e3
    mat.properties.poissons_ratio = 0.3
    mat.properties.yield_stress = 100.0

    mat.dtime = 0.25
    mat.dstrain .= 1.0e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]*mat.dtime
    integrate_material!(mat)
    @test isapprox(mat.stress+mat.dstress, [50.0, 0.0, 00.0, 0.0, 0.0, 0.0])

    material_postprocess_increment!(mat, dummyel, dummyip, mat.time+mat.dtime)

    integrate_material!(mat)
    @test isapprox(mat.stress+mat.dstress, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test isapprox(mat.properties.dplastic_multiplier, 0.0; atol=1.0e-12)

    material_postprocess_increment!(mat, dummyel, dummyip, mat.time+mat.dtime)

    integrate_material!(mat)
    @test mat.properties.dplastic_strain[1] > 0
    @test mat.properties.dplastic_strain[2] < 0
    @test mat.properties.dplastic_strain[3] < 0

    stress = mat.stress
    dstress = mat.dstress

    @test Materials.von_mises_yield(stress + dstress, 200) <= 0
end

@testset "Invalid potential is not defined" begin
    @test_throws ErrorException Material(ViscoPlastic, (:notNorton, [0.0, 0.0]))
end


text:
module TabularRL

using Lazy
import Base.getindex, Base.setindex!

import ..AbstractVFunction
import ..AbstractQFunction
import ..update!

export VFunction, QFunction

"""
    VFunction(dims::Integer...; init)
Tabular V Function
"""
mutable struct VFunction{W} <: AbstractVFunction
    values::W
end
VFunction(dims::Integer...; init=zeros) = VFunction(init(dims...))

weights(v::VFunction) = v.values
update!(value::VFunction, s, δ) = value[s] .+= δ

@forward VFunction.values Base.getindex
@forward VFunction.values Base.setindex!

#---------------------#
#
# State-Action Value Functions
#
#---------------------#

"""
    QFunction
Tabular Q Function
"""
mutable struct QFunction{W} <: AbstractQFunction
    values::W
end

"""
    QFunction(num_actions, state_dims...)
Creates a action-state value function of size num_actions x state_dims...
"""
QFunction(num_actions::Integer, state_dims::Integer...; init=zeros) = QFunction(init(num_actions, state_dims...))

@forward QFunction.values getindex
@forward QFunction.values setindex!

end # module TabularRL

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
using Knet, ParticleSwarmOptimization
include("Args.jl")
include("Functions.jl")

function main()
    gpu(false)
    args = parse_commandline()
    println("Parsed args:")

    for (arg,val) in args
        println("  $arg  =>  $val")
    end

    problem = RegularProblem(eval(parse(args["f"])), args["dim"], args["l"], args["h"])
    gbest = pso(problem, args["num"], args["epoch"], args["phi1"], args["phi2"], args["w"], args["vmin"], args["vmax"])

    println("\nBest: $(gbest.pbest_fitness), $(gbest.x)")
end

main()