text:
using DataStructures: OrderedDict, DefaultOrderedDict
const _Config = OrderedDict{Symbol,Any}
"""
    Config

Contains a set of configuration for systems. Configuration for a system contains parameter values.

# Examples
```julia-repl
julia> @config :S => :a => 1
Config for 1 system:
  S
    a = 1
```

See also: [`@config`](@ref)
"""
struct Config
    config::_Config
    Config(c...) = new(_Config(c...))
end

Base.getindex(c::Config, i) = c.config[i]
#HACK: handle :0
Base.getindex(c::Config, i::Int) = c.config[Symbol(i)]
Base.length(c::Config) = length(c.config)
Base.iterate(c::Config) = iterate(c.config)
Base.iterate(c::Config, i) = iterate(c.config, i)
Base.eltype(::Type{Config}) = Pair{Symbol,Any}

Base.get(c::Config, k, d) = get(c.config, k, d)
Base.haskey(c::Config, k) = haskey(c.config, k)

Base.:(==)(c::Config, d::Config) = c.config == d.config

Base.merge(f::Function, c::Config, D...) = merge(f, c.config, [d.config for d in D]...) |> Config

Base.show(io::IO, c::Config) = print(io, "<Config>")
Base.show(io::IO, ::MIME"text/plain", c::Config) = begin
    n = length(c)
    if n == 0
        print(io, "Config empty")
    elseif n == 1
        println(io, "Config for $n system:")
    else
        println(io, "Config for $n systems:")
    end
    f((s, C); color) = begin
        b = IOBuffer()
        x = IOContext(b, :color => color)
        print(x, "  ")
        printstyled(x, s, color=:light_magenta)
        K = keys(C)
        l = isempty(K) ? 0 : maximum(length.(string.(K)))
        for (k, v) in C
            println(x)
            print(x, "    ")
            printstyled(x, rpad(canonicalnamelabel(k), l), color=:light_blue)
            printstyled(x, " = ", color=:light_black)
            print(x, labelstring(v))
        end
        String(take!(b))
    end
    color = get(io, :color, false)
    join(io, f.(c; color), '\n')
end

configure(c::Config) = c
configure(c::AbstractDict) = configure(c...)
configure(c::Pair) = _configure(c.first, c.second)
configure(c::Tuple) = configure(c...)
configure(c::Vector) = configure.(c)
configure(c...) = merge(merge, configure.(c)...)
configure(::Nothing) = configure()
configure(c) = error("unrecognized configuration: $c")
configure() = Config()
_configure(k::Symbol, v) = Config(k => _configure(v))
_configure(k::String, v) = begin
    a = Symbol.(split(k, "."))
    n = length(a)
    if n == 2
        Config(a[1] => _configure(a[2] => v))
    elseif n == 1
        Config(a[1] => _configure(v))
    else
        error("unrecognized configuration key string: $k")
    end
end
_configure(S::Type{<:System}, l) = begin
    P = filter(istag(:parameter), [n.info for n in dependency(S).N])
    K = map(v -> (v.name, v.alias), P) |> Iterators.flatten
    C = _configure(l)
    for k in keys(C)
        (k ∉ K) && error("unrecognized parameter: $S => $k")
    end
    _configure(namefor(S), C)
end
_configure(k, v) = _configure(Symbol(k), v)
_configure(v) = _Config(v)
_configure(v::NamedTuple) = _Config(pairs(v))

parameterflatten(c::Config) = begin
    l = OrderedDict()
    for (s, d) in c
        for (k, v) in d
            l[(s, k)] = v
        end
    end
    l
end
parameterkeys(c::Config) = collect(keys(parameterflatten(c)))
parametervalues(c::Config) = collect(values(parameterflatten(c)))
parameterzip(K, V) = begin
    l = DefaultOrderedDict(OrderedDict)
    for ((s, k), v) in zip(K, V)
        l[s][k] = v
    end
    configure(l)
end

#TODO: wait until TOML 0.5 gets support
# using TOML
# loadconfig(c::AbstractString) = configure(TOML.parse(c))

option(c::Config, s::Symbol, k::Symbol) = begin
    v = get(c, s, missing)
    ismissing(v) ? v : get(v, k, missing)
end
option(c::Config, S::Vector{Symbol}, k::Symbol) = option(c, S, [k])
option(c::Config, s::Symbol, K::Vector{Symbol}) = option(c, [s], K)
option(c::Config, S::Vector{Symbol}, K::Vector{Symbol}) = begin
    v = missing
    for (s, k) in Iterators.product(S, K)
        #HACK: support private parameter (i.e. :S => :_a for S._a)
        k = uncanonicalname(k, s)
        v = option(c, s, k)
        !ismissing(v) && break
    end
    v
end

using DataStructures: OrderedSet
"""
    parameters(S; <keyword arguments>) -> Config

Extract a list of parameters defined for system `S`.

# Arguments
- `S::Type{<:System}`: type of system to be inspected.

# Keyword Arguments
- `alias=false`: show alias instead of parameter name.
- `recursive=false`: extract parameters from other systems declared in `S`.
- `exclude=()`: systems excluded in recurisve search.
- `scope=nothing`: evaluation scope; default is `S.name.module`.

# Examples
```julia-repl
julia> @system S(Controller) begin
           a: aaa => 1 ~ preserve(parameter)
       end;

julia> parameters(S)
Config for 1 system:
  S
    a = 1

julia> parameters(S; alias=true)
Config for 1 system:
  S
    aaa = 1

julia> parameters(S; recursive=true)
Config for 3 systems:
  Clock
    init = 0 hr
    step = 1 hr
  Context
  S
    a = 1

julia> parameters(S; recursive=true, exclude=(Context,))
Config for 1 system:
  S
    a = 1
```
"""
parameters(::Type{S}; alias=false, recursive=false, exclude=(), scope=nothing) where {S<:System} = begin
    #HACK: default evaluation scope is the module where S was originally defined
    isnothing(scope) && (scope = S.name.module)
    V = [n.info for n in dependency(S).N]
    P = filter(istag(:parameter), V)
    key = alias ? (v -> isnothing(v.alias) ? v.name : v.alias) : (v -> v.name)
    K = constsof(S) |> keys
    # evaluate only if parameter has no dependency on other variables
    val(v) = val(v, Val(isempty(v.args)))
    val(v, ::Val{false}) = missing
    val(v, ::Val{true}) = begin
        @gensym CS
        l = (:($k = $CS[$(Meta.quot(k))]) for k in K)
        b = @eval scope let $CS = Cropbox.constsof($S), $(l...); $(v.body) end
        u = @eval scope let $CS = Cropbox.constsof($S), $(l...); $(gettag(v, :unit)) end
        unitfy(b, u)
    end
    C = configure(namefor(S) => ((key(v) => val(v) for v in P)...,))
    if recursive
        T = OrderedSet([@eval scope $(v.type) for v in V])
        T = map(collect(T)) do t
            #HACK: not working for dynamic type (i.e. eltype(Vector{<:System}) = Any)
            et = eltype(t)
            et <: System ? et : t <: System ? t : nothing
        end
        filter!(!isnothing, T)
        filter!(t -> !any(t .<: exclude), T)
        X = (S, T..., exclude...) |> Set
        C = configure(parameters.(T; alias, recursive=true, exclude=X, scope)..., C)
    end
    C
end
#TODO: parameters(::System) to show current configurations (need proper handling of alias)

configmultiply(; base=()) = [configure(base)]
configmultiply(patches::Vector; base=()) = configmultiply(patches...; base)
configmultiply(patches...; base=()) = begin
    C = configexpand(patches[1]; base)
    for p in patches[2:end]
        C = [configexpand(p; base=c) for c in C] |> Iterators.flatten |> collect
    end
    C
end
configexpand(patch; base=()) = begin
    P = configure(patch)
    configs = if isempty(P)
        []
    else
        s, C = only(P)
        k, V = only(C)
        #HACK: allow single patch (i.e. `0 => :a => 1` instead of `1:2`)
        reshape([s => k => v for v in V], :)
    end
    configexpand(configs; base)
end
configexpand(configs::Vector; base=()) = configrebase(configs; base)
configrebase(configs::Vector; base=()) = isempty(configs) ? [configure(base)] : [configure(base, c) for c in configs]
configrebase(config; base=()) = configrebase([config]; base)

configreduce(a::Vector, b) = configure.(a, b)
configreduce(a, b::Vector) = configrebase(b; base=a)
configreduce(a::Vector, b::Vector) = configreduce.(a, b)
configreduce(a, b) = configure(a, b)
configreduce(a) = configure(a)
configreduce(a::Vector) = configure.(a)

using MacroTools: @capture
"""
    @config c.. -> Config | Vector{Config}

Construct a set or multiple sets of configuration.

A basic unit of configuration for a system `S` is represented by a pair in the form of `S => pv`. System name `S` is expressed in a symbol. If actual type of system is used, its name will be automatically converted to a symbol.

A parameter name and corresponding value is then represented by another pair in the form of `p => v`. When specifiying multiple parameters, a tuple of pairs like `(p1 => v1, p2 => v2)` or a named tuple like `(p1 = v1, p2 = v2)` can be used. Parameter name must be a symbol and should indicate a variable declared with `parameter` tag as often used by `preserve` state variable. For example, `:S => (:a => 1, :b => 2)` has the same meaning as `S => (a = 1, b = 2)` in the same scope.

Configurations for multiple systems can be concatenated by a tuple. Multiple elements in `c` separated by commas implicitly forms a tuple. For example, `:S => (:a => 1, :b => 2), :T => :x => 1` represents a set of configuration for two systems `S` and `T` with some parameters. When the same names of system or variable appears again during concatenation, it will be overriden by later ones in an order appeared in a tuple. For example, `:S => :a => 1, :S => :a => 2` results into `:S => :a => 2`. Instead of commas, `+` operator can be used in a similar way as `(:S => :a => 1) + (:S => :a => 2)`. Note parentheses placed due to operator precedence.

When multiple sets of configurations are needed, as in `configs` for [`simulate`](@ref), a vector of `Config` is used. This macro supports some convenient ways to construct a vector by composing simpler configurations. Prefix operator `!` allows *expansion* of any iterable placed in the configuration value. Infix operator `*` allows *multiplication* of a vector of configurations with another vector or a single configuration to construct multiple sets of configurations. For example, `!(:S => :a => 1:2)` is expanded into two sets of separate configurations `[:S => :a => 1, :S => :a => 2]`. `(:S => :a => 1:2) * (:S => :b => 0)` is multiplied into `[:S => (a = 1, b = 0), :S => (a = 2, b = 0)]`.

# Examples
```julia-repl
julia> @config :S => (:a => 1, :b => 2)
Config for 1 system:
  S
    a = 1
    b = 2
```

```julia-repl
julia> @config :S => :a => 1, :S => :a => 2
Config for 1 system:
  S
    a = 2
```

```julia-repl
julia> @config !(:S => :a => 1:2)
2-element Vector{Config}:
 <Config>
 <Config>
```

```julia-repl
julia> @config (:S => :a => 1:2) * (:S => :b => 0)
2-element Vector{Config}:
 <Config>
 <Config>
```
"""
macro config(ex)
    @capture(ex, +(P__) | P__)
    P = map(P) do p
        if @capture(p, !x_)
            :(Cropbox.configexpand($(esc(x))))
        elseif @capture(p, *(x__))
            :(Cropbox.configmultiply($(esc.(x)...)))
        else
            :(Cropbox.configreduce($(esc(p))))
        end
    end
    reduce(P) do a, b
        :(Cropbox.configreduce($a, $b))
    end
end

macro config(ex, exs...)
    :(Cropbox.@config(($(esc(ex)), $(esc.(exs)...))))
end

macro config()
    :(Cropbox.@config(()))
end

export Config, parameters, @config


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots

D = 4
k = 1.
μ = 2.3
c = 1.5

V = Constant(randn(D,D))
S1 = Constant(diagm(rand(D)/D))
S2 = Constant(diagm(rand(D)/D))
KL = Constant(randn(D,D))
KL = 0.5*(KL-KL')
R1 = V*S1*inv(V)
R2 = V*S2*inv(V)
Ss = (S1,S2)
RLs = (R1,R2)
QL = KL
for R in RLs
    mul!(QL, R', R, -1/2, 1)
end
#Put them in cMPS form
Ψ = InfiniteCMPS(QL, (R1,R2); gauge = :left)

h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + c * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2)
H = ∫(h, (-Inf,+Inf))

alg1 = LBFGS(; verbosity = 2, maxiter = 1000000, gradtol = 1e-4);
linalg = GMRES(krylovdim = 80; tol = 1e-5)

#Ψ, ρR, E, e, normgrad, numfg, history = groundstate6(H, Ψ, V, Ss; optalg = alg1, linalg = linalg)
αs,fs, dfs1, dfs2 = groundstate6(H, Ψ, V, Ss; optalg = alg1, linalg = linalg)
#αs = (αs[1:end-1] + αs[2:end])/2
#push!(αs,0.1)
display(plot(αs,[dfs1,dfs2]))
gui()

ext: jl


text:
import NLS_Solver

using ForwardDiff

export NLS_ForwardDiff_From_Model2Fit

# See details in NLS_ForwardDiff_From_Model2Fit doc
#
# In short: by export NLS_Solver, we avoid potential version
# compatibility problems. The only inconvenience is that user have to
# use the "NLS_Solver." prefix.
#
export NLS_Solver

# Important: ################
#
# To reduce number of memory allocations it is important to fully
# parameterize the NLS_ForwardDiff_From_Model2Fit structure
#
# NLS_ForwardDiff_From_Model2Fit{MODEL2FIT_TYPE <: Abstract_Model2Fit,
#                                X_ELEMENT_TYPE,
#                                Y_ELEMENT_TYPE,
#                                X_TYPE <: AbstractVector{X_ELEMENT_TYPE},
#                                Y_TYPE <: AbstractVector{Y_ELEMENT_TYPE}} <: NLS_Solver.AbstractNLS
#
# - Without parametrization ( NLS_ForwardDiff_From_Model2Fit <: NLS_Solver.AbstractNLS ):
#
# julia> @btime eval_r($nls,$θ)
#   2.724 μs (44 allocations: 928 bytes)
#
# - With parametrization:
#
# julia> @btime eval_r($nls,$θ)
#   842.590 ns (1 allocation: 160 bytes)
#
@doc raw"""

A wrapper that allows to use the
[NLS_Solver.jl](https://github.com/vincent-picaud/NLS_Solver.jl)
package to solve the nonlinear least squares problem associated to a
model plus its (X,Y) data.

You must construct an instance of nls problem as follows:
```julia
NLS_ForwardDiff_From_Model2Fit(fit_model::MODEL2FIT_TYPE,
                               X::X_TYPE,
                               Y::Y_TYPE)
```

You can then use the `NLS_Solver` package as usual, but with prefixed
`NLS_Solver.`

```@example
nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y)
conf = NLS_Solver.LevenbergMarquardt_Conf()
result = NLS_Solver.solve(nls,θ_init,conf)
```

**CAVEAT:**

Please note that you should **not** explicitly
```@example
using NLS_Solver
```
as this package is reexported from `NLS_Fit`.
Doing so avoids potential version compatibility problems.

"""
struct NLS_ForwardDiff_From_Model2Fit{MODEL2FIT_TYPE <: Abstract_Model2Fit,
                                      X_ELEMENT_TYPE,
                                      Y_ELEMENT_TYPE,
                                      X_TYPE <: AbstractVector{X_ELEMENT_TYPE},
                                      Y_TYPE <: AbstractVector{Y_ELEMENT_TYPE}} <: NLS_Solver.AbstractNLS
    _fit_model::MODEL2FIT_TYPE
    _X::X_TYPE
    _Y::Y_TYPE

    function NLS_ForwardDiff_From_Model2Fit(fit_model::MODEL2FIT_TYPE,
                                            X::X_TYPE,Y::Y_TYPE) where{MODEL2FIT_TYPE <: Abstract_Model2Fit,
                                                                       X_ELEMENT_TYPE,
                                                                       Y_ELEMENT_TYPE,
                                                                       X_TYPE <: AbstractVector{X_ELEMENT_TYPE},
                                                                       Y_TYPE <: AbstractVector{Y_ELEMENT_TYPE}}
        @assert length(X) == length(Y)

        new{MODEL2FIT_TYPE,X_ELEMENT_TYPE,Y_ELEMENT_TYPE,X_TYPE,Y_TYPE}(fit_model,X,Y)
    end
end

NLS_Solver.parameter_size(nls::NLS_ForwardDiff_From_Model2Fit) = parameter_size(nls._fit_model)
NLS_Solver.residue_size(nls::NLS_ForwardDiff_From_Model2Fit)  = length(nls._Y)

function NLS_Solver.eval_r(nls::NLS_ForwardDiff_From_Model2Fit,θ::AbstractVector)
    Y_model = eval_y(nls._fit_model,nls._X,θ)
    @. Y_model = nls._Y - Y_model
    Y_model
end

function NLS_Solver.eval_r_J(nls::NLS_ForwardDiff_From_Model2Fit, θ::AbstractVector)
    r = NLS_Solver.eval_r(nls,θ)
    J = ForwardDiff.jacobian(θ->NLS_Solver.eval_r(nls,θ), θ)
    r,J
end


text:
import VTKView
import Plots
import PyPlot
import GLMakie

include("examples/plotting.jl")

function makefigs()
    Plotter=PyPlot
    scene=plotting_multiscene(Plotter=Plotter)
    save("multiscene_pyplot.png",scene,Plotter=Plotter)

    Plotter=Plots
    scene=plotting_multiscene(Plotter=Plotter)
    save("multiscene_plots.png",scene,Plotter=Plotter)

    Plotter=GLMakie
    scene=plotting_multiscene(Plotter=Plotter)
    save("multiscene_glmakie.png",scene,Plotter=Plotter)

    Plotter=VTKView
    scene=plotting_multiscene(Plotter=Plotter)
    save("multiscene_vtkview.png",scene,Plotter=Plotter)

end


text:
using Dates

abstract type MyAbstractSimulateur end

export Simulateur
"""
Struct Simulateur

"""
struct Simulateur <: MyAbstractSimulateur

    Indicateur_prop::Indicateur
    name_folder::String

    function Simulateur(Ville::City,Individu::Population,rank::Int)

        A = Indicateur(Individu.nb_pop) ;

        if(ParamSimu.vMPI)
            if(ParamSimu.CALCULATEUR)
                name_folder = split(pwd(),"/")[end] ;
            else
                name_folder = "Res_MPI" ;
            end
        else
            name_folder = Dates.format(now(), "yyyy-mm-ddTHHMMSS") ;
            #mkdir("Res/$(lpad(name_folder,1))") ;
        end

        if(rank==0)
            mkdir("Res/$(lpad(name_folder,1))") ;
            mkdir("Res/$(lpad(name_folder,1))/txt_file/") ;
            mkdir("Res/$(lpad(name_folder,1))/json_file/") ;
            def_parametre_input(name_folder,Ville.name_file) ;
        end

        for i_MC = 1:ParamSimu.NB_MC
            StatutCovid = Statut_covid(Individu) ;
            compute_spreading(A,i_MC,StatutCovid,Ville,Individu) ;
        end

        return new(A,name_folder)

    end

end

function compute_spreading(A::Indicateur, iMC::Int64, statcovid::Statut_covid,Ville::City,Individu::Population)

    Jour = DefJournee(Individu.nb_pop,length(JourneeType.Ecolier)) ;

    for i_day=1:ParamSimu.NB_DAY

        whereis_local = Array{Int32}([]) ;
        if(ParamSimu.PRINT_ETAT)
            println("Calcul en cours pour le jour : ", i_day)
        end
        deplacement_v2!(Jour,Ville,Individu,statcovid) ;
        # On calcule la contamination possible de chaque individu
        contamination_sim = copy(statcovid.contamination) ;
        # Contamination au sein du 5ème arrondissement
        contamination_fct!(statcovid,whereis_local,i_day,Jour,Individu.tissu_work,Individu.statut) ;
        # Contamination en dehors du 5ème arrondissement
        incidence = statslieu_contamination(whereis_local,Ville,Individu) ;
        contamination_externe!(statcovid,incidence,Individu.statut,Individu.indiv_in_out) ;
        # On calcule les indicateurs de suivi de propagation de l'épidémie
        compute_indicateur!(A,iMC,i_day,contamination_sim,statcovid.contamination,Individu.indiv_in_out,statcovid.deb_mesure) ;
        if(ParamSimu.PRINT_ETAT)
            println("Nouvelle contamination : ", A.Nnew_cont_MC[i_day+1,iMC])
            println("Avec biais : ", A.Incidence_biais[i_day+1,iMC])
        end
        mise_a_jour!(statcovid,contamination_sim) ;
        miseajour_tracage!(statcovid,contamination_sim,statcovid.contamination) ;

        if((mesures.INTERVENTION) & (statcovid.mesure_on == false))
            mise_en_place!(statcovid,A.Incidence_biais[i_day+1,iMC]) ;
            #println("statut : ", statcovid.mesure_on)
        elseif(statcovid.mesure_on)
            miseajour_mesure!(statcovid,contamination_sim,statcovid.contamination) ;
            #println("Mise à jour du ptit confinement")
        end

    end

    return() ;

end


text:
module B
b = 2

module SubB
b2 = 2.5
end

end


text:
#=

Calculation of atmoshpere heat transport is not simple. Detail must be careful.
The standard version I found is here:
    http://www.atmos.albany.edu/facstaff/brose/classes/ATM623_Spring2015/Notes/Lectures/Lecture13%20--%20Heat%20transport.html

In particular, snow flux is not part of latent heat flux so must be added separately.

=#

using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("MapTransform.jl")
include("constants.jl")
include("CESMReader.jl")

using .CESMReader
using .MapTransform

function findITCZ(x, lat)

    j_ITCZ_south = -1

    if length(x) < 2
        throw(ErrorException("Need at least 2 points."))
    end

    if x[1] > 0.0
        throw(ErrorException("ITCZ not within this range: " * string(lat)))
    end

    for j = 2:length(x)
        if x[j] > 0.0  # Found the north one
            j_ITCZ_south = j - 1
            break
        end
    end

    if j_ITCZ_south == -1
        throw(ErrorException("Unable to find ITCZ."))
    end

    if any( x[(j_ITCZ_south + 1):end] .<= 0 )
        println(x)
        throw(ErrorException("More than two zero flux nodes found."))
    end

    x_south = x[j_ITCZ_south]
    x_north = x[j_ITCZ_south+1]
    lat_south = lat[j_ITCZ_south]
    lat_north = lat[j_ITCZ_south+1]

    return lat_south + (lat_north - lat_south) / (x_north - x_south) * ( 0.0 - x_south)

end

function mreplace(x)
    return replace(x, missing=>NaN)
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Ocn domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global area = mreplace(ds["area"][:])
    global mask = mreplace(ds["mask"][:])
    global lat  = mreplace(ds["yc"][:])

    mask .= 1
    area .*= 4π * Re^2 / sum(area)
end

lat_bnd = collect(Float64, -90:1:90)
r = MapTransform.Relation(
    lat = lat,
    area = area,
    mask = mask,
    lat_bnd = lat_bnd,
)

ITCZ_lat_idx = (91-10):(91+10) # 10S ~ 10N
ITCZ_lat_bnd = lat_bnd[ITCZ_lat_idx]

_proxy = area * 0 .+ 1.0
sum_valid_area = MapTransform.∫∂a(r, _proxy)[end]

println("Sum of valid area: ", sum_valid_area, "; ratio: ", sum_valid_area / sum(area))

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    FSNT, FLNT, FSNS, FLNS, SHFLX, LHFLX, PRECSC, PRECSL = getData(fh, ["FSNT", "FLNT", "FSNS", "FLNS", "SHFLX", "LHFLX", "PRECSC", "PRECSL"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    # FLX converted to (+) => upward, (-) => downward

    FFLX_TOA = ( - FSNT + FLNT )
    FFLX_SFC = ( - FSNS + FLNS )
    HFLX_SFC = ( SHFLX + LHFLX )
    PRECS    = ( PRECSC + PRECSL )

    _data = - ( FFLX_TOA - FFLX_SFC - HFLX_SFC - PRECS * ρ_fw * L_fusion )  # ~ - ∂T/∂t

    Nt = end_t - beg_t + 1

    global NetRad_TOA = zeros(Float64, length(r.lat_bnd)-1, Nt)
    global TFLX_CONV  = zeros(Float64, length(r.lat_bnd)-1, Nt)
    global AHT        = zeros(Float64, length(r.lat_bnd), Nt)

    for t = 1:Nt
        v = view(_data, :, :, t)
        TFLX_CONV[:, t] = MapTransform.transform(r, v)
        AHT[:, t]       = MapTransform.∫∂a(r, v)

        v2 = view(FFLX_TOA, :, :, t)
        NetRad_TOA[:, t] = MapTransform.transform(r, v2)
    end

    global years = Int(Nt/12)
    global ITCZ_lat = zeros(Float64, years)
    global AHT_AM   = zeros(Float64, length(r.lat_bnd), years)

    for y = 1:years
        AHT_AM[:, y] = mean(AHT[:, ((y-1)*12+1):(y*12)], dims=2)[:, 1]
#        ITCZ_lat[y] = findITCZ(AHT_AM[ITCZ_lat_idx, y], ITCZ_lat_bnd)
    end

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "year", years)
    defDim(ds, "lat_bnd", length(r.lat_bnd))
    defDim(ds, "lat",     length(r.lat_bnd)-1)

    for (varname, vardata, vardim, attrib) in [
        ("NetRad_TOA",  NetRad_TOA,                          ("lat", "time"), Dict()),
        ("TFLX_CONV",   TFLX_CONV,                           ("lat", "time"), Dict()),
        ("AHT",        AHT,                                  ("lat_bnd", "time"), Dict()),
        ("TFLX_CONV_MEAN",   mean(TFLX_CONV, dims=2)[:, 1],  ("lat",), Dict()),
        ("AHT_MEAN",         mean(AHT, dims=2)[:, 1],        ("lat_bnd",), Dict()),

        ("AHT_AM",           AHT_AM,                         ("lat_bnd", "year",), Dict()),
#        ("ITCZ_lat",         ITCZ_lat,                       ("year",), Dict()),

        ("lat_bnd",          r.lat_bnd,                      ("lat_bnd",), Dict()),
        ("area",             r.∂a,                           ("lat",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Code up correct formulation of stochastic forcing][Code up correct formulation of stochastic forcing:1]]
function createNoiseRateForcingHIT(D_a, D_s)
   c = zeros(3,3,3,3);
   delta(i,j) = (i==j) ? 1 : 0;
   for i in 1:3
      for j in 1:3
         for k in 1:3
            for l in 1:3
   	    term1 = -(1/3)*sqrt(D_s/5)*delta(i,j)*delta(k,l);
   	    term2 = (1/2)*(sqrt(D_s/5) + sqrt(D_a/3))*delta(i,k)*delta(j,l);
   	    term3 = (1/2)*(sqrt(D_s/5) - sqrt(D_a/3))*delta(i,l)*delta(j,k);
   	    c[i,j,k,l] = term1 + term2 + term3;
   	    end
         end
      end
   end
   return c;
end
# Code up correct formulation of stochastic forcing:1 ends here

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
using Firefly
using Documenter

setup = quote
    using Firefly
    using Random
    Random.seed!(42)
end

DocMeta.setdocmeta!(Firefly, :DocTestSetup, setup; recursive = true)

makedocs(;
    modules = [Firefly],
    authors = "Miles Lucas <mdlucas@hawaii.edu>",
    repo = "https://github.com/juliahci/Firefly.jl/blob/{commit}{path}#L{line}",
    sitename = "Firefly.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliahci.github.io/Firefly.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
    ],
    strict = true
)

deploydocs(;
    repo = "github.com/JuliaHCI/Firefly.jl",
    push_preview = true
)


text:
# This file was generated by JuDoc, do not modify it. # hide
@show round.(randn(), digits=2)


text:
struct RuntimeDevice{D}
    device::D
end
default_device() = RuntimeDevice(default_device(RUNTIME[]))
default_device(::typeof(HSA_rt)) = get_default_agent()

struct RuntimeQueue{Q}
    queue::Q
end
default_queue() = default_queue(default_device())
default_queue(device::RuntimeDevice) = RuntimeQueue(default_queue(RUNTIME[], device))
default_queue(::typeof(HSA_rt), device::RuntimeDevice) =
    get_default_queue(device.device)
get_device(queue::RuntimeQueue{HSAQueue}) = RuntimeDevice(queue.queue.agent)

default_isa(device::RuntimeDevice{HSAAgent}) =
    get_first_isa(device.device)

struct RuntimeEvent{E}
    event::E
end
create_event(exe) = RuntimeEvent(create_event(RUNTIME[], exe))
Base.wait(event::RuntimeEvent, exe) = wait(event.event, exe)

"Tracks the completion and status of a kernel's execution."
struct HSAStatusSignal
    signal::HSASignal
    exe::HSAExecutable
end
create_event(::typeof(HSA_rt), exe) = HSAStatusSignal(HSASignal(), exe.exe)
function Base.wait(event::RuntimeEvent{HSAStatusSignal}; check_exceptions=true, cleanup=true, kwargs...)
    wait(event.event.signal; kwargs...) # wait for completion signal
    unpreserve!(event) # allow kernel-associated objects to be freed
    exe = event.event.exe
    mod = EXE_TO_MODULE_MAP[exe].value
    agent = exe.agent
    ex = nothing
    signal_handle = event.event.signal.signal[].handle
    if haskey(exe.globals, :__global_exception_flag)
        # Check if any wavefront for this kernel threw an exception
        ex_flag = get_global(exe, :__global_exception_flag)
        ex_flag_ptr = Base.unsafe_convert(Ptr{Int64}, ex_flag)
        ex_flag_value = Base.unsafe_load(ex_flag_ptr)
        if ex_flag_value != 0
            ex_strings = String[]
            if check_exceptions && haskey(exe.globals, :__global_exception_ring)
                # Check for and collect any exceptions, and clear their slots
                ex_ring = get_global(exe, :__global_exception_ring)
                ex_ring_ptr = Base.unsafe_convert(Ptr{Ptr{ExceptionEntry}}, ex_ring)
                ex_ring_ptr = unsafe_load(ex_ring_ptr)
                while (ex_ring_value = unsafe_load(ex_ring_ptr)).kern != 1
                    if ex_ring_value.kern == signal_handle
                        push!(ex_strings, unsafe_string(convert(Ptr{UInt8}, ex_ring_value.ptr)))
                        # FIXME: Write rest of entry first, then CAS 0 to kern field
                        unsafe_store!(ex_ring_ptr, ExceptionEntry())
                    end
                    ex_ring_ptr += sizeof(ExceptionEntry)
                end
            end
            unique!(ex_strings)
            ex = KernelException(RuntimeDevice(agent), isempty(ex_strings) ? nothing : join(ex_strings, '\n'))
        end
    end
    if cleanup
        # Clean-up malloc'd data
        for idx in length(mod.metadata):-1:1
            metadata_value = mod.metadata[idx]
            if metadata_value.kern == signal_handle
                @debug "Cleaning up data: $metadata_value"
                Mem.free(metadata_value.buf)
                deleteat!(mod.metadata, idx)
            end
        end
    end
    ex !== nothing && throw(ex)
end

struct RuntimeExecutable{E}
    exe::E
end
create_executable(device, entry, obj; globals=()) =
    RuntimeExecutable(create_executable(RUNTIME[], device, entry, obj; globals=globals))
function create_executable(::typeof(HSA_rt), device, entry, obj; globals=())
    # link with ld.lld
    @assert ld_lld_path != "" "ld.lld was not found; cannot link kernel"
    path_exe = mktemp() do path_o, io_o
        write(io_o, obj)
        flush(io_o)
        path_exe = path_o*".exe"
        run(`$ld_lld_path -shared -o $path_exe $path_o`)
        path_exe
    end
    data = read(path_exe)
    rm(path_exe)

    return HSAExecutable(device.device, data, entry; globals=globals)
end
get_global(exe::RuntimeExecutable, sym::Symbol) =
    get_global(exe.exe, sym)

struct RuntimeKernel{K}
    kernel::K
end
create_kernel(device, exe, entry, args) =
    RuntimeKernel(create_kernel(RUNTIME[], device, exe, entry, args))
create_kernel(::typeof(HSA_rt), device, exe, entry, args) =
    HSAKernelInstance(device.device, exe.exe, entry, args)
launch_kernel(queue, kern, event; kwargs...) =
    launch_kernel(RUNTIME[], queue, kern, event; kwargs...)
function launch_kernel(::typeof(HSA_rt), queue, kern, event;
                       groupsize=nothing, gridsize=nothing)
    signal = event.event isa HSAStatusSignal ? event.event.signal : event.event
    launch!(queue.queue, kern.kernel, signal;
                       workgroup_size=groupsize, grid_size=gridsize)
end
barrier_and!(queue, events::Vector{<:RuntimeEvent}) = barrier_and!(queue, map(x->x.event,events))
barrier_and!(queue, signals::Vector{HSAStatusSignal}) = barrier_and!(queue, map(x->x.signal,signals))
barrier_or!(queue, events::Vector{<:RuntimeEvent}) = barrier_or!(queue, map(x->x.event,events))
barrier_or!(queue, signals::Vector{HSAStatusSignal}) = barrier_or!(queue, map(x->x.signal,signals))


text:
module Intervals

export Interval

export inf, sup

################################################################################
#
# Types, Constructors, and Accessors
#
################################################################################

abstract AbstractInterval{T}

immutable Interval{T} <: AbstractInterval{T}
    inf::T
    sup::T

    # Diagonal dispatch
    function Interval{T}(i::T, s::T)
        if i > s
            new(s, i)
        else
            new(i, s)
        end
    end
end

# try to promote types
function Interval(i,s)
    t = promote(i,s)
    Interval{typeof(t[1])}(t...)
end

@inline inf(i::Interval) = i.inf
@inline sup(i::Interval) = i.sup

################################################################################
#
# Interval Arithmetic
#
################################################################################

@inline function (+)(i1::Interval, i2::Interval)
    Interval(i1.inf+i2.inf, i1.sup+i2.sup)
end

@inline function (-)(i1::Interval, i2::Interval)
    Interval(i1.inf-i2.inf, i1.sup-i2.sup)
end

@inline function (*)(i1::Interval, i2::Interval)
    a = i1.inf
    b = i1.sup
    c = i2.inf
    d = i2.sup
    ac = a*c
    ad = a*d
    bc = b*c
    bd = b*d
    Interval(min(ac,ad,bc,bd), max(ac,ad,bc,bd))
end

@inline function (/)(i1::Interval, i2::Interval)
    a = i1.inf
    b = i1.sup
    c = i2.inf
    d = i2.sup
    ac = a/c
    ad = a/d
    bc = b/c
    bd = b/d
    Interval(min(ac,ad,bc,bd), max(ac,ad,bc,bd))
end

for f in (:acos, :asin, :atan, :cos, :cosh, :exp, :log, :log10, :sin, :sinh,
          :sqrt, :tan, :tanh, :sind, :cosd, :tand, :sinpi, :cospi, :asind,
          :acosd, :atand, :sec, :csc, :cot, :secd, :cscd, :cotd, :asec, :acsc,
          :acot, :asecd, :acscd, :acotd, :sech, :csch, :coth, :asinh, :acosh,
          :atanh, :asech, :acsch, :acoth, :sinc, :cosc, :rad2deg, :log2, :log1p,
          :exp2, :exp10, :expm1, :abs, :abs2, :sign, :isqrt, :cbrt, :erf, :erfc,
          :erfcx, :erfi, :dawson, :erfinv, :erfcinv, :real, :imag, :conj, :cis,
          :factorial, :gamma, :lgamma, :lfact, :digamma, :trigamma, :airyai,
          :airyprime, :airyaiprime, :airybiprime, :besselj0, :besselj1,
          :bessely0, :bessely1, :eta, :zeta)
    @eval begin
        (Base.$f)(i::Interval) = Interval((Base.$f)(i.inf), (Base.$f)(i.sup))
    end
end

################################################################################
#
# Allen's Interval Algebra
#
################################################################################

@inline function (==)(i1::Interval, i2::Interval)
    i1.inf == i2.inf && i1.sup == i2.sup
end

@doc """
i1: _____

i2:           _____
""" ->
@inline function before(i1::Interval, i2::Interval)
    i1.sup < i2.inf # this is guaranteed because of the inner constructor
end

@doc """
i1: _____

i2:      _____
""" ->
@inline function meets(i1::Interval, i2::Interval)
    i1.sup == i2.inf && i1.inf < i2.inf
end

@doc """
i1:       _____

i2: ___________
""" ->
@inline function finishes(i1::Interval, i2::Interval)
    i1.sup == i2.sup && i1.inf > i2.inf
end

@doc """
i1:      _______

i2: ___________________
""" ->
@inline function during(i1::Interval, i2::Interval)
    i1.inf > i2.inf && i1.sup < i2.sup
end

@doc """
i1: ________

i2: _________________
"""->
@inline function starts(i1::Interval, i2::Interval)
    i1.inf == i2.inf && i1.sup < i2.sup
end

@doc """
i1: ________

i2:    ____________
""" ->
@inline function overlaps(i1::Interval, i2::Interval)
    i1.sup > i2.inf && i1.inf < i2.inf
end

end # module


text:
"""
    function addUpdateLocalInteraction_!(p::Program_,platform::String)

Generate the functions related with Local Interaction Updates.
"""
function addUpdateLocalInteraction_!(p::Program_,platform::String)

    if "UpdateLocalInteraction" in keys(p.agent.declaredUpdates)

        #Construct update computation function
        fcompute = loop_[p.neighbors](p,p.agent.declaredUpdates["UpdateLocalInteraction"],platform)
        fcompute = vectorize_(p.agent,fcompute,p,interaction=true)
        fcompute = wrapInFunction_(:locInterCompute_!,fcompute)
        push!(p.declareF.args,fcompute)

        #Wrap both functions in a clean step function
        if inexpr(p.declareF,:cleanLocalInteraction_!)
            f = push!(p.declareF.args,
                :(
                    function locInterStep_!(ARGS_)
                        @platformAdapt cleanLocalInteraction_!(ARGS_)
                        @platformAdapt locInterCompute_!(ARGS_)
                        return
                    end
                )
                )
            push!(p.execInit.args,
                    :(locInterStep_!(ARGS_))
                )
            push!(p.execInloop.args,
                :(locInterStep_!(ARGS_))
                )
            push!(p.execAfter.args,
                    :(locInterStep_!(ARGS_))
                )
        else
            push!(p.execInit.args,
                    :(@platformAdapt locInterCompute_!(ARGS_))
                )
            push!(p.execInloop.args,
                    :(@platformAdapt locInterCompute_!(ARGS_))
                )
            push!(p.execAfter.args,
                    :(@platformAdapt locInterCompute_!(ARGS_))
                )
        end

    end

    return nothing
end


text:
# a read-only wrapper for an AbstractArray
# because julia's type system is not sufficiently expressive, this will
# only work for arrays (needs to be subtype of DenseArray in order
# to call BLAS functions

import ArrayViews: aview, unsafe_aview, Subs

ContiguousArrays{T, N} =   Union{Array{T, N}, ContiguousView{T, N}, UnsafeContiguousView{T, N}}

const AllArrayViews = Union{ArrayView, UnsafeArrayView}

"""
  A read only wrapper type for AbstractArrays. getindex is defined for
  this function but setindex! is not.
"""
struct ROView{T, N, P <: AbstractArray} <: DenseArray{T, N}
  parent::P

  function ROView{T, N, P}(a::ContiguousArrays) where {T, N, P}
    return new(a)
  end
end

"""
  Outer constructor
"""
function ROView(A::ContiguousArrays{T, N}) where {T, N}
  return ROView{T, N, typeof(A)}(A)
end

# define some functions by delegating to parent array
import Base: size, length, getindex, parent, convert, IndexStyle, convert, unsafe_convert, pointer

size(A::ROView) = size(A.parent)
length(A::ROView) = length(A.parent)

# it appears there is no splatting penalty in this case
@inline getindex(A::ROView, idx...) = getindex(A.parent, idx...)

parent(A::ROView) = A.parent

IndexStyle(A::ROView) = IndexStyle(A.parent)

pointer(A::ROView) = pointer(A.parent)

unsafe_convert(::Type{Ptr{T}}, A::ROView) where {T} = pointer(A)

"""
  Create a sview and then make a read only view of it
"""
@inline ro_sview(A::Array, idx...) = ROView(sview(A, idx...))
@inline ro_sview(A::ROView, idx...) = ROView(sview(parent(A), idx...))
@inline ro_sview(A::AllArrayViews, idx...) = ROView(sview(A, idx...))

# enable taking a sview of a ro_sview
# this must return an ro_sview to preserve the read-only property
# up to 5 argument constructors have to be listed explicitly to avoid ambiguity
# problem with Arrayviews (a varargs is less specific than a positional argument
# and ROView us a subtype of DenseArray

@inline aview(A::ROView, i1::Subs) = ROView(aview(parent(A), i1))
@inline aview(A::ROView, i1::Subs, i2::Subs) = ROView(aview(parent(A), i1, i2))
@inline aview(A::ROView, i1::Subs, i2::Subs, i3::Subs) = ROView(aview(parent(A), i1, i2, i3))
@inline aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs) = ROView(aview(parent(A), i1, i2, i3, i4))
@inline aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs) = ROView(aview(parent(A), i1, i2, i3, i4, i5))
@inline aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) = ROView(aview(parent(A), i1, i2, i3, i4, i5, I...))

@inline unsafe_aview(A::ROView, i1::Subs) = ROView(unsafe_aview(parent(A), i1))
@inline unsafe_aview(A::ROView, i1::Subs, i2::Subs) = ROView(unsafe_aview(parent(A), i1, i2))
@inline unsafe_aview(A::ROView, i1::Subs, i2::Subs, i3::Subs) = ROView(unsafe_aview(parent(A), i1, i2, i3))
@inline unsafe_aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs) = ROView(unsafe_aview(parent(A), i1, i2, i3, i4))
@inline unsafe_aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs) = ROView(unsafe_aview(parent(A), i1, i2, i3, i4, i5))
@inline unsafe_aview(A::ROView, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) = ROView(unsafe_aview(parent(A), i1, i2, i3, i4, i5, I...))

# useful aliases
ROVector{T, P} =  ROView{T, 1, P}
ROMatrix{T, P} =  ROView{T, 2, P}
ROArray{T, N, P} =  ROView{T, N, P}
ROContiguousView{T, N, P} =  ROView{T, N, ContiguousView{T, N, P}}

convert(::Type{Array{T}}, x::ROView{T2}) where {T, T2} = ROView(convert(Array{T}, parent(x)))

copy!(dest::AbstractArray{T, N}, src::ROView{T2, N}) where {T, N, T2} = copy!(dest, src.parent)

import Base.reinterpret
"""
  Reinterpret an ArrayViews.UnsafeContiguousView (needed by BLAS wrappers).

  It is not possible, in general, to reinterpert a regular ContiguousView, so
  the small linear algebra functions will have to avoid calling BLAS routines in
  that case
"""
function reinterpret(::Type{T}, a::UnsafeContiguousView{S},
                     dims::NTuple{N, Int}) where {T, S, N}

# T is the output element type
# S is the input element type
# N is the output dimensionality

  # print the same warnings as the Base code
  # TODO: make sure these don't have a performance impact
  if !isbits(T)
    throw(ArgumentError("cannot reinterpret array to non bitstype $(T)"))
  end

  if !isbits(S)
    throw(ArgumentError("cannot reinterpret array from non bitstype $(S)"))
  end

  # ensure size in memory is same
  size_in = prod(size(a))*sizeof(S)
  len_out = prod(dims)
  size_out = len_out*sizeof(T)

  if size_in != size_out
    throw(DimensionMismatch("input and output sizes must be consistent"))
  end

  pnew = Ptr{T}(a.ptr)
  new_arr = UnsafeContiguousView(pnew, dims)

  return new_arr
end

function reinterpret(::Type{T}, a::UnsafeContiguousView{S}) where {T, S}
# Blas mat-vec wrapper uses this

  size_in = sizeof(S)
  size_out = sizeof(T)

  # if this isn't exact, the length check later will error
  first_dim = size(a, 1)*div(size_in, size_out)
  new_dims = (first_dim, size(a)[2:end]...)

  reinterpret(T, a, new_dims)
end

function reinterpret(::Type{T}, a::UnsafeContiguousView{T}) where T

  return a
end

function reinterpret(::Type{T}, a::ContiguousView{T}) where T

  return a
end

function reinterpret(::Type{T}, a::ROView{S}, dims::NTuple{N, Int}) where {T, S, N}
# reinterpret the parent and wrap it in a ROView

  return ROView(reinterpret(T, parent(a), dims))
end

function reinterpret(::Type{T}, a::ROView{T}) where T

  return a
end
#=
function reinterpret{T}(::Type{T}, a::Union{ContiguousView{T}, UnsafeContiguousView{T}})
# it is always possible to reinterpret to the same type

  return a
end
=#

#TODO: reinterpret a ROView as well


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(softplus(-0.8925450194703983 + 0.6110890427800215*softplus(-0.16626426792027038 + 0.44845618413236066*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.1740395626916258*softplus(0.8063063231762757 + 0.4905142959916642*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8017072131372296*softplus(0.9210510999614945 + -0.4935968211971775*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7966067754659663*softplus(-0.7328670881809516 + 0.6102429954796835*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + softplus(-0.7008872510227278 + -0.258173144464545*softplus(-0.16626426792027038 + 0.44845618413236066*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.49550787111793193*softplus(0.8063063231762757 + 0.4905142959916642*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8462943069272191*softplus(0.9210510999614945 + -0.4935968211971775*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.05869131996649157*softplus(-0.7328670881809516 + 0.6102429954796835*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + softplus(-0.2643911691309153 + -0.31660511055293883*softplus(-0.16626426792027038 + 0.44845618413236066*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.09136517532701127*softplus(0.8063063231762757 + 0.4905142959916642*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.36400896662330684*softplus(0.9210510999614945 + -0.4935968211971775*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.4876638787119827*softplus(-0.7328670881809516 + 0.6102429954796835*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + softplus(0.6679332747384503 + -0.6061375406574494*softplus(-0.16626426792027038 + 0.44845618413236066*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7532640510772808*softplus(0.8063063231762757 + 0.4905142959916642*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.9959397069462161*softplus(0.9210510999614945 + -0.4935968211971775*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8544693505755525*softplus(-0.7328670881809516 + 0.6102429954796835*softplus(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*softplus(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*softplus(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*softplus(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
mutable struct Binding
    name::EXPR
    val::Union{Binding,EXPR,SymbolServer.SymStore,Nothing}
    type::Union{Binding,EXPR,SymbolServer.SymStore,Nothing}
    refs::Vector{Any}
    prev::Union{Binding,SymbolServer.SymStore,Nothing}
    next::Union{Binding,SymbolServer.SymStore,Nothing}
end
Binding(x::EXPR) = Binding(CSTParser.get_name(x), x, nothing, [], nothing, nothing)

function Base.show(io::IO, b::Binding)
    printstyled(io, "Binding(", Expr(b.name),
        b.type === nothing ? "" : ":: ",
        b.refs isa Vector ? "($(length(b.refs)) refs))" : ")", color=:blue)
end

hasbinding(x::EXPR) = hasmeta(x) && hasbinding(x.meta)
bindingof(x) = nothing
bindingof(x::EXPR) = bindingof(x.meta)

hasref(x::EXPR) = hasmeta(x) && hasref(x.meta)
refof(x::EXPR) = hasmeta(x) ? x.meta.ref : nothing

function gotoobjectofref(x::EXPR)
    r = refof(x)
    if r isa SymbolServer.SymStore
        return r
    elseif r isa Binding

    end
end

# Note to self, check consistency of marking self-reference of bindings (i.e.
# for, `function f end` we resolve `f` to itself at this stage.)
function mark_bindings!(x::EXPR, state)
    if hasbinding(x)
        return
    end
    if !hasmeta(x)
        x.meta = Meta()
    end
    if is_binary_call(x)
        if kindof(x[2]) === CSTParser.Tokens.EQ
            if CSTParser.is_func_call(x[1])
                name = CSTParser.get_name(x)
                mark_binding!(x)
                mark_sig_args!(x[1])
                if isidentifier(name)
                    setref!(name, bindingof(x))
                end
            elseif is_curly(x[1])
                mark_typealias_bindings!(x)
            elseif !is_getfield(x[1])
                mark_binding!(x[1], x)
            end
        elseif kindof(x[2]) === CSTParser.Tokens.ANON_FUNC
            mark_binding!(x[1], x)
        end
    elseif is_where(x)
        for i = 3:length(x)
            ispunctuation(x[i]) && continue
            mark_binding!(x[i])
        end
    elseif typof(x) === CSTParser.For
        markiterbinding!(x[2])
    elseif typof(x) === CSTParser.Generator
        for i = 3:length(x)
            ispunctuation(x[i]) && continue
            markiterbinding!(x[i])
        end
    elseif typof(x) === CSTParser.Filter
        for i = 1:length(x) - 2
            ispunctuation(x[i]) && continue
            markiterbinding!(x[i])
        end
    elseif typof(x) === CSTParser.Flatten && length(x) === 1 && length(x[1]) >= 3 && length(x[1][1]) >= 3
        for i = 3:length(x[1][1])
            ispunctuation(x[1][1][i]) && continue
            markiterbinding!(x[1][1][i])
        end
        for i = 3:length(x[1])
            ispunctuation(x[1][i]) && continue
            markiterbinding!(x[1][i])
        end
    elseif typof(x) === CSTParser.Do
        if is_tuple(x[3])
            for i in 1:length(x[3])
                ispunctuation(x[3][i]) && continue
                mark_binding!(x[3][i])
            end
        end
    elseif typof(x) === FunctionDef || CSTParser.defines_macro(x)
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.Function, [], nothing, nothing)
        if isidentifier(name)
            setref!(name, bindingof(x))
        end
        mark_sig_args!(CSTParser.get_sig(x))
    elseif CSTParser.defines_module(x)
        x.meta.binding = Binding(x[2], x, CoreTypes.Module, [], nothing, nothing)
        setref!(x[2], bindingof(x))
    elseif typof(x) === CSTParser.Try && length(x) > 3 && isidentifier(x[4])
        mark_binding!(x[4])
        setref!(x[4], bindingof(x[4]))
    elseif CSTParser.defines_datatype(x)
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.DataType, [], nothing, nothing)
        if isidentifier(name)
            setref!(name, bindingof(x))
        end
        mark_parameters(CSTParser.get_sig(x))
        if CSTParser.defines_struct(x) # mark field block
            blocki = typof(x[3]) === CSTParser.Block ? 3 : 4
            for i in 1:length(x[blocki])
                CSTParser.defines_function(x[blocki][i]) && continue
                mark_binding!(x[blocki][i])
            end
        end
    elseif typof(x) === CSTParser.Local
        if length(x) == 2
            if isidentifier(x[2])
                mark_binding!(x[2])
                setref!(x[2], bindingof(x[2]))
            elseif is_tuple(x[2])
                for i = 1:length(x[2])
                    if isidentifier(x[2][i])
                        mark_binding!(x[2][i])
                        setref!(x[2][i], bindingof(x[2][i]))
                    end
                end
            end
        end

    end
end

function mark_binding!(x::EXPR, val=x)
    if is_kwarg(x) || (is_declaration(x) && is_tuple(x[1]))
        mark_binding!(x[1], x)
    elseif is_tuple(x) || is_parameters(x)
        for arg in x
            ispunctuation(arg) && continue
            mark_binding!(arg, val)
        end
    elseif typof(x) === CSTParser.InvisBrackets
        mark_binding!(CSTParser.rem_invis(x), val)
    elseif !(is_unary_call(x) && kindof(x[1]) === CSTParser.Tokens.DECLARATION)
        if !hasmeta(x)
            x.meta = Meta()
        end
        x.meta.binding = Binding(CSTParser.get_name(x), val, nothing, [], nothing, nothing)
    end
    return x
end

function mark_parameters(sig::EXPR)
    signame = CSTParser.rem_where_subtype(sig)
    if is_curly(signame)
        for i = 3:length(signame) - 1
            ispunctuation(signame[i]) && continue
            mark_binding!(signame[i])
        end
    end
    return sig
end

function markiterbinding!(iter::EXPR)
    if is_binary_call(iter) && kindof(iter[2]) in (CSTParser.Tokens.EQ, CSTParser.Tokens.IN, CSTParser.Tokens.ELEMENT_OF)
        mark_binding!(iter[1], iter)
    elseif typof(iter) === CSTParser.Block
        for i = 1:length(iter)
            ispunctuation(iter[i]) && continue
            markiterbinding!(iter[i])
        end
    end
    return iter
end

function mark_sig_args!(x::EXPR)
    if is_call(x) || is_tuple(x)
        if typof(x[1]) === CSTParser.InvisBrackets && is_declaration(x[1][2])
            mark_binding!(x[1][2])
        end
        for i = 2:length(x) - 1
            a = x[i]
            if is_parameters(a)
                for j = 1:length(a)
                    aa = a[j]
                    if !ispunctuation(aa)
                        mark_binding!(aa)
                    end
                end
            elseif !ispunctuation(a)
                mark_binding!(a)
            end
        end
    elseif is_where(x)
        for i in 3:length(x)
            if !ispunctuation(x[i])
                mark_binding!(x[i])
            end
        end
        mark_sig_args!(x[1])
    elseif is_invis_brackets(x)
        mark_sig_args!(CSTParser.rem_invis(x))
    elseif is_binary_call(x)
        if kindof(x[2]) == CSTParser.Tokens.DECLARATION
            mark_sig_args!(x[1])
        else
            mark_binding!(x[1])
            mark_binding!(x[3])
        end
    elseif is_unary_call(x) && typof(x[2]) == CSTParser.InvisBrackets
        mark_binding!(x[2][2])
    end
end

function mark_typealias_bindings!(x::EXPR)
    mark_binding!(x, x)
    setscope!(x, Scope(x))
    for i = 2:length(x[1])
        if isidentifier(x[1][i])
            mark_binding!(x[1][i])
        elseif is_binary_call(x[1][i]) && kindof(x[1][i][2]) === CSTParser.Tokens.ISSUBTYPE && isidentifier(x[1][i][1])
            mark_binding!(x[1][i][1])
        end
    end
    return x
end

rem_wheres_decls(x) = (is_where(x) || is_declaration(x)) ? x[1] : x

function is_in_funcdef(x)
    if !(parentof(x) isa EXPR)
        return false
    elseif is_where(parentof(x)) || typof(parentof(x)) === CSTParser.InvisBrackets
        return is_in_funcdef(parentof(x))
    elseif typof(parentof(x)) === FunctionDef || CSTParser.is_assignment(parentof(x))
        return true
    else
        return false
    end
end

function _in_func_def(x::EXPR)
    # only called in WhereOpCall
    # check 1st arg contains a call (or op call)
    ex = rem_wheres_decls(x[1])

    !(is_call(ex) || (is_binary_call(ex) && kindof(ex[2]) !== CSTParser.Tokens.DOT) || is_unary_call(ex)) && return false

    # check parent is func def
    return is_in_funcdef(x)
end

function add_binding(x, state, scope=state.scope)
    if bindingof(x) isa Binding
        b = bindingof(x)
        b.prev = nothing
        b.next = nothing
        if isidentifier(b.name)
            name = valofid(b.name)
        elseif is_macroname(b.name)
            name = string(Expr(b.name))
        elseif isoperator(b.name)
            name = string(Expr(b.name))
        else
            return
        end
        # check for global marker
        if isglobal(name, scope)
            scope = _get_global_scope(state.scope)
        end

        if CSTParser.defines_macro(x)
            scope.names[string("@", name)] = b
            mn = CSTParser.get_name(x)
            if isidentifier(mn)
                setref!(mn, b)
            end
        else
            if name_is_getfield(b.name)# && b.type == CoreTypes.Function
                # Question: should `b.name` include the getfield itself?
                # Binding name is part of a getfield : `A.name` so is either
                # 1. Overloading of a function
                # 2. Setting of a field or property
                # We only care about 1.
                resolve_ref(parentof(parentof(b.name))[1], scope, state)
                lhs_ref = refof_maybe_getfield(parentof(parentof(b.name))[1])
                if lhs_ref === nothing
                    # We don't know what we're assigning to, do nothing
                else
                    if lhs_ref isa SymbolServer.ModuleStore && haskey(lhs_ref.vals, Symbol(name))
                        # Overloading
                        tls = retrieve_toplevel_scope(b.val)
                        tls === nothing && return # Shouldn't happen
                        if haskey(tls.names, name) && eventually_overloads(tls.names[name], lhs_ref.vals[Symbol(name)], state.server)
                            # Though we're explicitly naming a function for overloading, it has already been imported to the toplevel scope.
                            overload_method(tls, b, VarRef(lhs_ref.name, Symbol(name)))
                            b.prev = tls.names[name]
                            b.prev.next = b
                            tls.names[name] = b
                        elseif isexportedby(name, lhs_ref)
                            tls.names[name] = b
                            b.prev = maybe_lookup(lhs_ref[Symbol(name)], state.server)
                        else
                            overload_method(tls, b, VarRef(lhs_ref.name, Symbol(name)))
                        end
                    elseif lhs_ref isa Binding && lhs_ref.type == CoreTypes.Module
                        if hasscope(lhs_ref.val) && haskey(scopeof(lhs_ref.val).names, name)
                            b.prev = scopeof(lhs_ref.val).names[name]
                            scopeof(lhs_ref.val).names[name] = b
                        end
                    end
                end

            elseif scopehasbinding(scope, name)
                b.prev = scope.names[name]
                scope.names[name] = b
                b.prev.next = b
            else
                scope.names[name] = b
            end
            # hoist binding for inner constructor to parent scope
            if CSTParser.defines_struct(scope.expr) && CSTParser.defines_function(x) && parentof(scope) isa Scope
                return add_binding(x, state, parentof(scope))
            end
        end
        infer_type(b, scope, state)
    elseif bindingof(x) isa SymbolServer.SymStore
        scope.names[valofid(x)] = bindingof(x)
    end
end

name_is_getfield(x) = parentof(x) isa EXPR && parentof(parentof(x)) isa EXPR && is_getfield_w_quotenode(parentof(parentof(x)))

"""
eventually_overloads(b, x, server)

"""
eventually_overloads(b::Binding, ss::SymbolServer.SymStore, server) = (b.val == ss || b.prev == ss) || (b.prev !== nothing && eventually_overloads(b.prev, ss, server))

eventually_overloads(b::Binding, ss::SymbolServer.VarRef, server) = eventually_overloads(b, maybe_lookup(ss, server), server)

eventually_overloads(b, ss, server) = false

function hoist_prev_binding(b, name, scope, state)
    scope === nothing && return
    if scope.modules !== nothing
        if scopehasmodule(scope, Symbol(valofid(parentof(parentof(b.name))[1]))) # this scope (s1) has a module with matching name
            mod = getscopemodule(scope, Symbol(valofid(parentof(parentof(b.name))[1])))
            if mod isa SymbolServer.ModuleStore && haskey(mod, Symbol(name))
                b.prev = maybe_lookup(mod[Symbol(name)], state.server)
            end
        end
        return # We've reached a scope that loads modules, no need to keep searching upwards
    end
    return hoist_prev_binding(b, name, parentof(scope), state)
end

isglobal(name, scope) = false
isglobal(name::String, scope) = scopehasbinding(scope, "#globals") && name in scope.names["#globals"].refs

function mark_globals(x::EXPR, state)
    if typof(x) === CSTParser.Global
        if !scopehasbinding(state.scope, "#globals")
            state.scope.names["#globals"] = Binding(EXPR(CSTParser.IDENTIFIER, EXPR[], 0, 0, "#globals", CSTParser.NoKind, false, nothing, nothing), nothing, nothing, [], nothing, nothing)
        end
        for i = 2:length(x)
            if isidentifier(x[i]) && !scopehasbinding(state.scope, valofid(x[i]))
                push!(state.scope.names["#globals"].refs, valofid(x[i]))
            end
        end

    end
end

function name_extends_imported_method(b::Binding)
    if b.type == CoreTypes.Function && hasparent(b.name) && is_getfield(parentof(b.name))
        if refof_maybe_getfield(parentof(b.name)[1]) !== nothing

        end
    end
end


text:
# 1st version. No type annotations!
struct Point
    x
    y
end

# this is syntactically same as Point
struct Point2
    x::Any
    y::Any
end

# this uses concrete type
struct Point3
    x::UInt8
    y::UInt8
end

Point3(0x01, 0x01) |> sizeof
Point(0x01, 0x01) |> sizeof
Point(Int128(1), Int128(1)) |> sizeof

#=
julia> Point3(0x01, 0x01) |> sizeof
2

julia> Point(0x01, 0x01) |> sizeof
16

julia> Point(Int128(1), Int128(1)) |> sizeof
16
=#

# this one use concrete type as well
struct Point4
    x::Int128
    y::Int128
end

Point4(Int128(1), Int128(1)) |> sizeof
#=
julia> Point4(Int128(1), Int128(1)) |> sizeof
32
=#

#--- parametric type version

# restart REPL and define this
struct Point{T <: Real}
    x::T
    y::T
end

p = Point(0x01, 0x01)
sizeof(p)
p = Point(Int128(1), Int128(1))
sizeof(p)
#=
julia> p = Point(0x01, 0x01)
Point{UInt8}(0x01, 0x01)

julia> sizeof(p)
2

julia> p = Point(Int128(1), Int128(1))
Point{Int128}(1, 1)

julia> sizeof(p)
32
=#

#--- performance test

struct PointAny
    x::Any
    y::Any
end

using Statistics: mean

function center(points::AbstractVector{T}) where T
    return T(
        mean(p.x for p in points),
        mean(p.y for p in points))
end

make_points(T::Type, n) = [T(rand(), rand()) for _ in 1:n]

using BenchmarkTools

points = make_points(PointAny, 100_000);
@btime center($points)

#=
julia> points = make_points(PointAny, 100_000);

julia> @btime center($points)
  5.384 ms (200007 allocations: 3.05 MiB)
PointAny(0.5008816012742725, 0.4987658522813566)
=#

points = make_points(Point, 100_000);
@btime center($points)

#=
julia> @btime center($points)
  207.027 μs (2 allocations: 32 bytes)
Point4{Float64}(0.5001556466497694, 0.4990206734988199)
=#


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

using JuliaFEM
using JuliaFEM.Testing

@testset "1d strain" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0, 0.0],
        2 => [1.0, 1.0, 1.0])
    u = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0, 0.0],
        2 => [1.0, 1.0, 1.0])
    element = Element(Seg2, [1, 2])
    update!(element, "geometry", X)
    detJ = element([0.0], 0.0, Val{:detJ})
    info("detJ = $detJ")
    @test isapprox(detJ, sqrt(3)/2)
    J = element([0.0], 0.0, Val{:Jacobian})
    info("J = $J")
    @test isapprox(J, [0.5 0.5 0.5])
    update!(element, "displacement", u)
    # FIXME
#   gradu = element("displacement", [0.0], 0.0, Val{:Grad})
#   info("1d bar: ∇u = $gradu")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestBLAS

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasReal, BlasComplex

srand(100)
## BLAS tests - testing the interface code to BLAS routines
@testset for elty in [Float32, Float64, ComplexF32, ComplexF64]
    @testset "syr2k!" begin
        U = randn(5,2)
        V = randn(5,2)
        if elty == ComplexF32 || elty == ComplexF64
            U = complex.(U, U)
            V = complex.(V, V)
        end
        U = convert(Array{elty, 2}, U)
        V = convert(Array{elty, 2}, V)
        @test tril(LinearAlgebra.BLAS.syr2k('L','N',U,V)) ≈ tril(U*transpose(V) + V*transpose(U))
        @test triu(LinearAlgebra.BLAS.syr2k('U','N',U,V)) ≈ triu(U*transpose(V) + V*transpose(U))
        @test tril(LinearAlgebra.BLAS.syr2k('L','T',U,V)) ≈ tril(transpose(U)*V + transpose(V)*U)
        @test triu(LinearAlgebra.BLAS.syr2k('U','T',U,V)) ≈ triu(transpose(U)*V + transpose(V)*U)
    end

    if elty in (ComplexF32, ComplexF64)
        @testset "her2k!" begin
            U = randn(5,2)
            V = randn(5,2)
            U = complex.(U, U)
            V = complex.(V, V)
            U = convert(Array{elty, 2}, U)
            V = convert(Array{elty, 2}, V)
            @test tril(LinearAlgebra.BLAS.her2k('L','N',U,V)) ≈ tril(U*V' + V*U')
            @test triu(LinearAlgebra.BLAS.her2k('U','N',U,V)) ≈ triu(U*V' + V*U')
            @test tril(LinearAlgebra.BLAS.her2k('L','C',U,V)) ≈ tril(U'*V + V'*U)
            @test triu(LinearAlgebra.BLAS.her2k('U','C',U,V)) ≈ triu(U'*V + V'*U)
        end
    end

    o4 = fill(elty(1), 4)
    z4 = zeros(elty, 4)

    I4 = Matrix{elty}(I, 4, 4)
    I43 = Matrix{elty}(I, 4, 3)
    L4 = tril(fill(elty(1), 4,4))
    U4 = triu(fill(elty(1), 4,4))
    Z4 = zeros(elty, (4,4))

    elm1 = convert(elty, -1)
    el2 = convert(elty, 2)
    v14 = convert(Vector{elty}, [1:4;])
    v41 = convert(Vector{elty}, [4:-1:1;])

    let n = 10
        @testset "dot products" begin
            if elty <: Real
                x1 = convert(Vector{elty}, randn(n))
                x2 = convert(Vector{elty}, randn(n))
                @test BLAS.dot(x1,x2) ≈ sum(x1.*x2)
                @test_throws DimensionMismatch BLAS.dot(x1,rand(elty, n + 1))
            else
                z1 = convert(Vector{elty}, complex.(randn(n),randn(n)))
                z2 = convert(Vector{elty}, complex.(randn(n),randn(n)))
                @test BLAS.dotc(z1,z2) ≈ sum(conj(z1).*z2)
                @test BLAS.dotu(z1,z2) ≈ sum(z1.*z2)
                @test_throws DimensionMismatch BLAS.dotc(z1,rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.dotu(z1,rand(elty, n + 1))
            end
        end
        @testset "iamax" begin
            if elty <: Real
                x = convert(Vector{elty}, randn(n))
                @test BLAS.iamax(x) == argmax(abs.(x))
            else
                z = convert(Vector{elty}, complex.(randn(n),randn(n)))
                @test BLAS.iamax(z) == argmax(map(x -> abs(real(x)) + abs(imag(x)), z))
            end
        end
        @testset "axp(b)y" begin
            if elty <: Real
                x1 = convert(Vector{elty}, randn(n))
                x2 = convert(Vector{elty}, randn(n))
                α  = rand(elty)
                β  = rand(elty)
                @test BLAS.axpy!(α,copy(x1),copy(x2)) ≈ α*x1 + x2
                @test BLAS.axpby!(α,copy(x1),β,copy(x2)) ≈ α*x1 + β*x2
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.axpby!(α, copy(x1), β, rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 1:n)
                @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 0:div(n,2), copy(x2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 0:(div(n, 2) - 1))
                @test BLAS.axpy!(α,copy(x1),1:n,copy(x2),1:n) ≈ x2 + α*x1
            else
                z1 = convert(Vector{elty}, complex.(randn(n), randn(n)))
                z2 = convert(Vector{elty}, complex.(randn(n), randn(n)))
                α  = rand(elty)
                @test BLAS.axpy!(α, copy(z1), copy(z2)) ≈ z2 + α * z1
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), 1:div(n, 2), copy(z2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 0:div(n,2), copy(z2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 1:div(n,2), copy(z2), 0:(div(n, 2) - 1))
                @test BLAS.axpy!(α,copy(z1),1:n,copy(z2),1:n) ≈ z2 + α*z1
            end
        end
        @testset "nrm2, iamax, and asum for StridedVectors" begin
            a = rand(elty,n)
            b = view(a,2:2:n,1)
            @test BLAS.nrm2(b) ≈ norm(b)
            if elty <: Real
                @test BLAS.asum(b) ≈ sum(abs.(b))
                @test BLAS.iamax(b) ≈ argmax(abs.(b))
            else
                @test BLAS.asum(b) ≈ sum(abs.(real(b))) + sum(abs.(imag(b)))
                @test BLAS.iamax(b) == argmax(map(x -> abs(real(x)) + abs(imag(x)), b))
            end
        end
        # scal
        α = rand(elty)
        a = rand(elty,n)
        @test BLAS.scal(n,α,a,1) ≈ α * a

        @testset "trsv" begin
            A = triu(rand(elty,n,n))
            @testset "Vector and SubVector" for x in (rand(elty, n), view(rand(elty,2n),1:2:2n))
                @test A\x ≈ BLAS.trsv('U','N','N',A,x)
                @test_throws DimensionMismatch BLAS.trsv('U','N','N',A,Vector{elty}(undef,n+1))
            end
        end
        @testset "ger, her, syr" for x in (rand(elty, n), view(rand(elty,2n), 1:2:2n)),
            y in (rand(elty,n), view(rand(elty,3n), 1:3:3n))

            A = rand(elty,n,n)
            α = rand(elty)

            @test BLAS.ger!(α,x,y,copy(A)) ≈ A + α*x*y'
            @test_throws DimensionMismatch BLAS.ger!(α,Vector{elty}(undef,n+1),y,copy(A))

            A = rand(elty,n,n)
            A = A + transpose(A)
            @test issymmetric(A)
            @test triu(BLAS.syr!('U',α,x,copy(A))) ≈ triu(A + α*x*transpose(x))
            @test_throws DimensionMismatch BLAS.syr!('U',α,Vector{elty}(undef,n+1),copy(A))

            if elty <: Complex
                A = rand(elty,n,n)
                A = A + A'
                α = real(α)
                @test triu(BLAS.her!('U',α,x,copy(A))) ≈ triu(A + α*x*x')
                @test_throws DimensionMismatch BLAS.her!('U',α,Vector{elty}(undef,n+1),copy(A))
            end
        end
        @testset "copy" begin
            x1 = convert(Vector{elty}, randn(n))
            x2 = convert(Vector{elty}, randn(n))
            BLAS.copyto!(x2, 1:n, x1, 1:n)
            @test x2 == x1
            @test_throws DimensionMismatch BLAS.copyto!(x2, 1:n, x1, 1:(n - 1))
            @test_throws ArgumentError BLAS.copyto!(x1, 0:div(n, 2), x2, 1:(div(n, 2) + 1))
            @test_throws ArgumentError BLAS.copyto!(x1, 1:(div(n, 2) + 1), x2, 0:div(n, 2))
        end
        # trmv
        A = triu(rand(elty,n,n))
        x = rand(elty,n)
        @test BLAS.trmv('U','N','N',A,x) ≈ A*x
        @testset "symmetric/Hermitian multiplication" begin
            x = rand(elty,n)
            A = rand(elty,n,n)
            Aherm = A + A'
            Asymm = A + transpose(A)
            @testset "symv and hemv" begin
                @test BLAS.symv('U',Asymm,x) ≈ Asymm*x
                offsizevec, offsizemat = Array{elty}.(undef,(n+1, (n,n+1)))
                @test_throws DimensionMismatch BLAS.symv!('U',one(elty),Asymm,x,one(elty),offsizevec)
                @test_throws DimensionMismatch BLAS.symv('U',offsizemat,x)
                if elty <: BlasComplex
                    @test BLAS.hemv('U',Aherm,x) ≈ Aherm*x
                    @test_throws DimensionMismatch BLAS.hemv('U',offsizemat,x)
                    @test_throws DimensionMismatch BLAS.hemv!('U',one(elty),Aherm,x,one(elty),offsizevec)
                end
            end

            @testset "symm error throwing" begin
                Cnn, Cnm, Cmn = Matrix{elty}.(undef,((n,n), (n,n-1), (n-1,n)))
                @test_throws DimensionMismatch BLAS.symm('L','U',Cnm,Cnn)
                @test_throws DimensionMismatch BLAS.symm('R','U',Cmn,Cnn)
                @test_throws DimensionMismatch BLAS.symm!('L','U',one(elty),Asymm,Cnn,one(elty),Cmn)
                @test_throws DimensionMismatch BLAS.symm!('L','U',one(elty),Asymm,Cnn,one(elty),Cnm)
                if elty <: BlasComplex
                    @test_throws DimensionMismatch BLAS.hemm('L','U',Cnm,Cnn)
                    @test_throws DimensionMismatch BLAS.hemm('R','U',Cmn,Cnn)
                    @test_throws DimensionMismatch BLAS.hemm!('L','U',one(elty),Aherm,Cnn,one(elty),Cmn)
                    @test_throws DimensionMismatch BLAS.hemm!('L','U',one(elty),Aherm,Cnn,one(elty),Cnm)
                end
            end
        end
        @testset "trmm error throwing" begin
            Cnn, Cmn, Cnm = Matrix{elty}.(undef,((n,n), (n+1,n), (n,n+1)))
            @test_throws DimensionMismatch BLAS.trmm('L','U','N','N',one(elty),triu(Cnn),Cmn)
            @test_throws DimensionMismatch BLAS.trmm('R','U','N','N',one(elty),triu(Cnn),Cnm)
        end

        #trsm
        A = triu(rand(elty,n,n))
        B = rand(elty,(n,n))
        @test BLAS.trsm('L','U','N','N',one(elty),A,B) ≈ A\B

        #will work for SymTridiagonal,Tridiagonal,Bidiagonal!
        @testset "banded matrix mv" begin
            @testset "gbmv" begin
                TD  = Tridiagonal(rand(elty,n-1),rand(elty,n),rand(elty,n-1))
                x   = rand(elty,n)
                #put TD into the BLAS format!
                fTD = zeros(elty,3,n)
                fTD[1,2:n] = TD.du
                fTD[2,:] = TD.d
                fTD[3,1:n-1] = TD.dl
                @test BLAS.gbmv('N',n,1,1,fTD,x) ≈ TD*x
            end
            #will work for SymTridiagonal only!
            @testset "sbmv" begin
                if elty <: BlasReal
                    ST  = SymTridiagonal(rand(elty,n),rand(elty,n-1))
                    x   = rand(elty,n)
                    #put TD into the BLAS format!
                    fST = zeros(elty,2,n)
                    fST[1,2:n] = ST.ev
                    fST[2,:] = ST.dv
                    @test BLAS.sbmv('U',1,fST,x) ≈ ST*x
                else
                    dv = real(rand(elty,n))
                    ev = rand(elty,n-1)
                    bH = zeros(elty,2,n)
                    bH[1,2:n] = ev
                    bH[2,:] = dv
                    fullH = diagm(0 => dv, -1 => conj(ev), 1 => ev)
                    @test BLAS.hbmv('U',1,bH,x) ≈ fullH*x
                end
            end
        end
    end

    @testset "gemv" begin
        @test all(BLAS.gemv('N', I4, o4) .== o4)
        @test all(BLAS.gemv('T', I4, o4) .== o4)
        @test all(BLAS.gemv('N', el2, I4, o4) .== el2 * o4)
        @test all(BLAS.gemv('T', el2, I4, o4) .== el2 * o4)
        @test_throws DimensionMismatch BLAS.gemv('N',I43,o4)
        o4cp = copy(o4)
        @test_throws DimensionMismatch BLAS.gemv!('T',one(elty),I43,o4,one(elty),o4cp)
        @test_throws DimensionMismatch BLAS.gemv!('C',one(elty),I43,o4,one(elty),o4cp)
        @test all(BLAS.gemv!('N', one(elty), I4, o4, elm1, o4cp) .== z4)
        @test all(o4cp .== z4)
        o4cp[:] = o4
        @test all(BLAS.gemv!('T', one(elty), I4, o4, elm1, o4cp) .== z4)
        @test all(o4cp .== z4)
        @test all(BLAS.gemv('N', U4, o4) .== v41)
        @test all(BLAS.gemv('N', U4, o4) .== v41)
    end
    @testset "gemm" begin
        @test all(BLAS.gemm('N', 'N', I4, I4) .== I4)
        @test all(BLAS.gemm('N', 'T', I4, I4) .== I4)
        @test all(BLAS.gemm('T', 'N', I4, I4) .== I4)
        @test all(BLAS.gemm('T', 'T', I4, I4) .== I4)
        @test all(BLAS.gemm('N', 'N', el2, I4, I4) .== el2 * I4)
        @test all(BLAS.gemm('N', 'T', el2, I4, I4) .== el2 * I4)
        @test all(BLAS.gemm('T', 'N', el2, I4, I4) .== el2 * I4)
        @test all(LinearAlgebra.BLAS.gemm('T', 'T', el2, I4, I4) .== el2 * I4)
        I4cp = copy(I4)
        @test all(BLAS.gemm!('N', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('N', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('T', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('T', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        @test all(BLAS.gemm('N', 'N', I4, U4) .== U4)
        @test all(BLAS.gemm('N', 'T', I4, U4) .== L4)
        @test_throws DimensionMismatch BLAS.gemm!('N','N', one(elty), I4, I4, elm1, Matrix{elty}(I, 5, 5))
        @test_throws DimensionMismatch BLAS.gemm!('N','N', one(elty), I43, I4, elm1, I4)
        @test_throws DimensionMismatch BLAS.gemm!('T','N', one(elty), I43, I4, elm1, I43)
        @test_throws DimensionMismatch BLAS.gemm!('N','T', one(elty), I43, I43, elm1, I43)
        @test_throws DimensionMismatch BLAS.gemm!('T','T', one(elty), I43, I43, elm1, Matrix{elty}(I, 3, 4))
    end
    @testset "gemm compared to (sy)(he)rk" begin
        if eltype(elm1) <: Complex
            @test all(triu(BLAS.herk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
            @test all(tril(BLAS.herk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
            @test all(triu(BLAS.herk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
            @test all(tril(BLAS.herk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
            @test all(triu(BLAS.herk('U', 'C', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
            @test all(tril(BLAS.herk('L', 'C', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
            @test all(triu(BLAS.herk('U', 'C', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
            @test all(tril(BLAS.herk('L', 'C', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
            ans = similar(L4)
            @test all(tril(BLAS.herk('L','C', L4)) .== tril(BLAS.herk!('L', 'C', real(one(elty)), L4, real(zero(elty)), ans)))
            @test all(LinearAlgebra.copytri!(ans, 'L') .== LinearAlgebra.BLAS.gemm('T', 'N', L4, L4))
            @test_throws DimensionMismatch BLAS.herk!('L','N',real(one(elty)),Matrix{elty}(I, 5, 5),real(one(elty)), Matrix{elty}(I, 6, 6))
        else
            @test all(triu(BLAS.syrk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
            @test all(tril(BLAS.syrk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
            @test all(triu(BLAS.syrk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
            @test all(tril(BLAS.syrk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
            @test all(triu(BLAS.syrk('U', 'T', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
            @test all(tril(BLAS.syrk('L', 'T', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
            @test all(triu(BLAS.syrk('U', 'T', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
            @test all(tril(BLAS.syrk('L', 'T', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
            ans = similar(L4)
            @test all(tril(BLAS.syrk('L','T', L4)) .== tril(BLAS.syrk!('L', 'T', one(elty), L4, zero(elty), ans)))
            @test all(LinearAlgebra.copytri!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
            @test_throws DimensionMismatch BLAS.syrk!('L','N',one(elty), Matrix{elty}(I, 5, 5),one(elty), Matrix{elty}(I, 6, 6))
        end
    end
end

@testset "syr for eltype $elty" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
    A = rand(elty, 5, 5)
    @test triu(A[1,:] * transpose(A[1,:])) ≈ BLAS.syr!('U', one(elty), A[1,:], zeros(elty, 5, 5))
    @test tril(A[1,:] * transpose(A[1,:])) ≈ BLAS.syr!('L', one(elty), A[1,:], zeros(elty, 5, 5))
    @test triu(A[1,:] * transpose(A[1,:])) ≈ BLAS.syr!('U', one(elty), view(A, 1, :), zeros(elty, 5, 5))
    @test tril(A[1,:] * transpose(A[1,:])) ≈ BLAS.syr!('L', one(elty), view(A, 1, :), zeros(elty, 5, 5))
end

@testset "her for eltype $elty" for elty in (Complex{Float32}, Complex{Float64})
    A = rand(elty, 5, 5)
    @test triu(A[1,:] * A[1,:]') ≈ BLAS.her!('U', one(real(elty)), A[1,:], zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:]') ≈ BLAS.her!('L', one(real(elty)), A[1,:], zeros(elty, 5, 5))
    @test triu(A[1,:] * A[1,:]') ≈ BLAS.her!('U', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:]') ≈ BLAS.her!('L', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5))
end

struct WrappedArray{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end

Base.size(A::WrappedArray) = size(A.A)
Base.getindex(A::WrappedArray, i::Int) = A.A[i]
Base.getindex(A::WrappedArray{T, N}, I::Vararg{Int, N}) where {T, N} = A.A[I...]
Base.setindex!(A::WrappedArray, v, i::Int) = setindex!(A.A, v, i)
Base.setindex!(A::WrappedArray{T, N}, v, I::Vararg{Int, N}) where {T, N} = setindex!(A.A, v, I...)
Base.unsafe_convert(::Type{Ptr{T}}, A::WrappedArray{T}) where T = Base.unsafe_convert(Ptr{T}, A.A)

@test_deprecated strides(WrappedArray(rand(5)))
@test_deprecated stride(WrappedArray(rand(5)), 1)

Base.stride(A::WrappedArray, i::Int) = stride(A.A, i)

@testset "strided interface blas" begin
    for elty in (Float32, Float64, ComplexF32, ComplexF64)
    # Level 1
        x = WrappedArray(elty[1, 2, 3, 4])
        y = WrappedArray(elty[5, 6, 7, 8])
        BLAS.blascopy!(2, x, 1, y, 2)
        @test y == WrappedArray(elty[1, 6, 2, 8])
        BLAS.scal!(2, elty(2), x, 1)
        @test x == WrappedArray(elty[2, 4, 3, 4])
        @test BLAS.nrm2(1, x, 2) == elty(2)
        @test BLAS.nrm2(x) == BLAS.nrm2(x.A)
        BLAS.asum(x) == elty(13)
        BLAS.axpy!(4, elty(2), x, 1, y, 1)
        @test y == WrappedArray(elty[5, 14, 8, 16])
        BLAS.axpby!(elty(2), x, elty(3), y)
        @test y == WrappedArray(elty[19, 50, 30, 56])
        @test BLAS.iamax(x) == 2
    # Level 2
        A = WrappedArray(elty[1 2; 3 4])
        x = WrappedArray(elty[1, 2])
        y = WrappedArray(elty[3, 4])
        @test BLAS.gemv!('N', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[13, 26])
        @test BLAS.gbmv!('N', 2, 1, 0, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[15, 40])
        @test BLAS.symv!('U', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[25, 60])
        @test BLAS.trmv!('U', 'N', 'N', A, y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[145, 240])
        @test BLAS.trsv!('U', 'N', 'N', A, y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[25,60])
        @test BLAS.ger!(elty(2), x, y, A) isa WrappedArray{elty,2}
        @test A == WrappedArray(elty[51 122; 103 244])
        @test BLAS.syr!('L', elty(2), x, A) isa WrappedArray{elty,2}
        @test A == WrappedArray(elty[53 122; 107 252])
    # Level 3
        A = WrappedArray(elty[1 2; 3 4])
        B = WrappedArray(elty[5 6; 7 8])
        C = WrappedArray(elty[9 10; 11 12])
        BLAS.gemm!('N', 'N', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([47 54; 97 112])
        BLAS.symm!('L', 'U', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([85 98; 173 200])
        BLAS.syrk!('U', 'N', elty(2), A, elty(1), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([95 120; 173 250])
        BLAS.syr2k!('U', 'N', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([163 244; 173 462])
        BLAS.trmm!('L', 'U', 'N', 'N', elty(2), A, B) isa WrappedArray{elty,2}
        @test B == WrappedArray([38 44; 56 64])
        BLAS.trsm!('L', 'U', 'N', 'N', elty(2), A, B) isa WrappedArray{elty,2}
        @test B == WrappedArray([20 24; 28 32])
    end
    for elty in (Float32, Float64)
    # Level 1
        x = WrappedArray(elty[1, 2, 3, 4])
        y = WrappedArray(elty[5, 6, 7, 8])
        @test BLAS.dot(2, x, 1, y, 2) == elty(19)
    # Level 2
        A = WrappedArray(elty[1 2; 3 4])
        x = WrappedArray(elty[1, 2])
        y = WrappedArray(elty[3, 4])
        BLAS.sbmv!('U', 1, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[17,24])
    end
    for elty in (ComplexF32, ComplexF64)
    # Level 1
        x = WrappedArray(elty[1+im, 2+2im, 3+3im, 4+im])
        y = WrappedArray(elty[5-im, 6-2im, 7-3im, 8-im])
        @test BLAS.dotc(2, x, 1, y, 2) == elty(12-26im)
        @test BLAS.dotu(2, x, 1, y, 2) == elty(26+12im)
    # Level 2
        A = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im])
        x = WrappedArray(elty[1+im, 2+2im])
        y = WrappedArray(elty[5-im, 6-2im])
        @test BLAS.hemv!('U', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[7+17im, 30+14im])
        BLAS.hbmv!('U', 1, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1}
        @test y == WrappedArray(elty[13+39im, 54+30im])
        @test BLAS.her!('L', real(elty(2)), x, A) isa WrappedArray{elty,2}
        @test A == WrappedArray(elty[5 2+2im; 11+3im 20])
    # Level 3
        A = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im])
        B = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im])
        C = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im])
        @test BLAS.hemm!('L', 'U', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([3+27im 6+38im; 35+27im 52+36im])
        @test BLAS.herk!('U', 'N', real(elty(2)), A, real(elty(1)), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([23 50+38im; 35+27im 152])
        @test BLAS.her2k!('U', 'N', elty(2), A, B, real(elty(1)), C) isa WrappedArray{elty,2}
        @test C == WrappedArray([63 138+38im; 35+27im 352])
    end
end

end # module TestBLAS

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using CSV, Gnuplot, DataFrames, ColorSchemes
# https://vega.github.io/vega-lite/examples/circle_bubble_health_income.html
#https://vega.github.io/vega-datasets/
#link = "https://vega.github.io/vega-datasets/data/gapminder-health-income.cvs"
# load file
data = CSV.read("gapminder-health-income.csv")
# selecting and plotting
x, y, popsize = data.income, data.health, data.population
mxs = maximum(popsize)
#
# aesthetics
@gp xlab ="income" ylab = "health" "set auto fix" :-
@gp :- title = "gapminder" :-
@gp :- "set offsets graph .1, graph .1, graph .1, graph .1" :-
@gp :- "set key off"
# the actual plot
@gp :- palette(:viridis)
@gp :- x y 2*(popsize./mxs).^(0.4) popsize./1e6 "w p pt 7 ps var lc pal z notit"
@gp :- "set logscale x" "set cblabel 'population x million'"
@gp :- "set logscale cb"
for i in 1:length(x)
    if data.population[i] > 2.5e8
        @gp :- "set label '$(data.country[i])' at $(x[i]),$(y[i]-1) font ',11' front"
    end
end
save(term ="pdf size 5,3", output ="gapminder2.pdf")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples