text:
using StatsBase, Random
Random.seed!(1)

grade = ["A","B","C","D","E"]
weightVect = Weights([0.1,0.2,0.1,0.2,0.4])

N = 10^6
data = sample(grade,weightVect,N)
[count(i->(i==g),data) for g in grade]/N

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
    GrayAtmosphericState{FT,FTA1D,FTA2D,I} <:
        AbstractAtmosphericState{FT,I,FTA1D}

Atmospheric conditions, used to compute optical properties with the gray atmosphere approximation

# Fields
$(DocStringExtensions.FIELDS)
"""
struct GrayAtmosphericState{
    FT <: AbstractFloat,
    FTA1D <: AbstractArray{FT, 1},
    FTA2D <: AbstractArray{FT, 2},
    I <: Int,
} <: AbstractAtmosphericState{FT, I, FTA1D}
    "Layer pressures `[Pa, mb]`; `(nlay, ncol)`"
    p_lay::FTA2D
    "Level pressures `[Pa, mb]`; `(nlay+1, ncol)`"
    p_lev::FTA2D
    "Layer temperatures `[K]`; `(nlay, ncol)`"
    t_lay::FTA2D
    "Level temperatures `[K]`; `(nlay+1,ncol)`"
    t_lev::FTA2D
    "Level Altitude `[m]`; `(nlay+1,ncol)`"
    z_lev::FTA2D
    "Surface temperatures `[K]`; `(ncol)`"
    t_sfc::FTA1D
    "lapse rate"
    α::FT
    "optical thickness parameter"
    d0::FTA1D
    "Number of layers."
    nlay::I
    "Number of columns."
    ncol::I
end
Adapt.@adapt_structure GrayAtmosphericState
#---------------------------------------------------------------

# This functions sets up a model temperature and pressure
# distributions for a gray atmosphere based on a pressure grid
# see Schneider 2004, J. Atmos. Sci. (2004) 61 (12): 1317–1340.
# https://doi.org/10.1175/1520-0469(2004)061<1317:TTATTS>2.0.CO;2

function setup_gray_as_pr_grid(
    nlay::Int,
    lat::FTA1D,
    p0::FT,
    pe::FT,
    param_set::RP.ARP,
    ::Type{DA},
    step = "linear",
) where {FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, DA}
    nlev = Int(nlay + 1)
    ncol = length(lat)
    p_lay = DA{FT}(undef, nlay, ncol) # layer mean pressure
    p_lev = DA{FT}(undef, nlev, ncol) # level pressure
    t_lay = DA{FT}(undef, nlay, ncol) # layer mean temperature
    t_lev = DA{FT}(undef, nlev, ncol) # level temperature
    z_lev = DA{FT}(undef, nlev, ncol) # level altitude
    t_sfc = DA{FT}(undef, ncol)       # surface temperature
    efac = log(p0 / pe) / nlay       # multiplication factor for each pressure step
    d0 = DA{FT}(undef, ncol)       # optical depth (function of latitude)
    Δp = (p0 - pe) / nlay          # Δp for linear pressure distribution
    te = FT(300)                   # global mean surface temperature (K)
    tt = FT(200)                   # skin temp at top of atmosphere (K)
    Δt = FT(60)
    α = FT(3.5)                   # lapse rate of radiative equillibrium
    r_d = RP.R_d(param_set)
    grav_ = RP.grav(param_set)
    args = (p_lev, p_lay, t_lev, t_lay, z_lev, t_sfc, lat, d0, efac, p0, pe, Δp, te, tt, Δt, α, r_d, grav_, nlay)
    device = array_device(p_lev)
    if device === CUDADevice()
        max_threads = 256
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        @cuda threads = (tx) blocks = (bx) setup_gray_as_pr_grid_CUDA!(ncol, args...)
    else # launcing Julia native multithreading kernel
        Threads.@threads for gcol in 1:ncol
            setup_gray_as_pr_grid_kernel!(args..., gcol)
        end
    end
    #------------------------------------------------
    return GrayAtmosphericState{eltype(t_sfc), typeof(t_sfc), typeof(p_lev), Int}(
        p_lay,
        p_lev,
        t_lay,
        t_lev,
        z_lev,
        t_sfc,
        α,
        d0,
        nlay,
        ncol,
    )
end

function setup_gray_as_pr_grid_CUDA!(ncol, args...)
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    if gcol ≤ ncol
        setup_gray_as_pr_grid_kernel!(args..., gcol)
    end
    return nothing
end
# This functions sets up a model temperature and pressure
# distributions for a gray atmosphere based on an altitude grid
# with internal GLL point distribution within each cell
# see Schneider 2004, J. Atmos. Sci. (2004) 61 (12): 1317–1340.
# https://doi.org/10.1175/1520-0469(2004)061<1317:TTATTS>2.0.CO;2

function setup_gray_as_alt_grid(
    lat::FTA1D,
    p0::FT,
    zend::FT,
    ncls::Int,
    poly_order::Int,
    param_set::RP.ARP,
    ::Type{DA},
) where {FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, DA}
    ncol = length(lat)
    nlev = Int(ncls + 1 + (poly_order - 1) * ncls)
    nlay = Int(nlev - 1)
    qm1 = poly_order + 1
    zst = FT(0)
    te = FT(300)                   # global mean surface temperature (K)
    tt = FT(200)                   # skin temp at top of atmosphere (K)
    Δt = FT(60)
    α = FT(3.5)           # lapse rate of radiative equillibrium

    p_lay = DA{FT}(undef, nlay, ncol)
    p_lev = DA{FT}(undef, nlev, ncol)
    t_lay = DA{FT}(undef, nlay, ncol)
    t_lev = DA{FT}(undef, nlev, ncol)
    z_lev = DA{FT}(undef, nlev, ncol)
    t_sfc = DA{FT}(undef, ncol)       # surface temperature
    d0 = DA{FT}(undef, ncol)       # optical depth (function of latitude)

    pts, _ = GaussQuadrature.legendre(FT, poly_order + 1, GaussQuadrature.both)
    pts = (pts .+ FT(1)) ./ FT(2)
    Δz_cell = zend / ncls

    for icol in 1:ncol
        #---------bot_at_1---------------------------------------
        z_lev[1, icol] = zst
        z_lev[nlev, icol] = zend

        for icls in 1:ncls
            for j in 1:(poly_order - 1)
                z_lev[(icls - 1) * poly_order + 1 + j, icol] =
                    z_lev[(icls - 1) * poly_order + 1, icol] + pts[1 + j] * Δz_cell
            end
            z_lev[icls * poly_order + 1, icol] = zst + icls * Δz_cell
        end

        ts = te + Δt * (FT(1 / 3) - sin(lat[icol])^2) # surface temp at a given latitude (K)
        d0[icol] = FT((ts / tt)^4 - 1) # optical depth

        p_lev[1, icol] = p0
        t_lev[1, icol] = tt * (1 + d0[icol] * (p_lev[icol, 1] / p0)^α)^FT(0.25)

        for ilay in 1:nlay
            H = RP.R_d(param_set) * t_lev[ilay, icol] / RP.grav(param_set)
            Δz_lay = abs(z_lev[ilay + 1, icol] - z_lev[ilay, icol])
            p_lev[ilay + 1, icol] = p_lev[ilay, icol] * exp(-Δz_lay / H)
            t_lev[ilay + 1, icol] = tt * (1 + d0[icol] * (p_lev[ilay + 1, icol] / p0)^α)^FT(0.25)

            t_lay[ilay, icol] = FT(0.5) * (t_lev[ilay, icol] + t_lev[ilay + 1, icol])
            p_lay[ilay, icol] = FT(0.5) * (p_lev[ilay, icol] + p_lev[ilay + 1, icol])
        end
        #--------------------------------------------------------
        t_sfc[icol] = t_lev[1, icol]
    end
    return GrayAtmosphericState{FT, DA{FT, 1}, DA{FT, 2}, Int}(
        p_lay,
        p_lev,
        t_lay,
        t_lev,
        z_lev,
        t_sfc,
        α,
        d0,
        nlay,
        ncol,
    )
end

#-------------------------------------------------------------------------
function setup_gray_as_pr_grid_kernel!(
    p_lev::FTA2D,
    p_lay::FTA2D,
    t_lev::FTA2D,
    t_lay::FTA2D,
    z_lev::FTA2D,
    t_sfc::FTA1D,
    lat::FTA1D,
    d0::FTA1D,
    efac::FT,
    p0::FT,
    pe::FT,
    Δp::FT,
    te::FT,
    tt::FT,
    Δt::FT,
    α::FT,
    r_d::FT,
    grav_::FT,
    nlay::Int,
    gcol::Int,
) where {FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, FTA2D <: AbstractArray{FT, 2}}
    ts = te + Δt * (FT(1) / FT(3) - sin(lat[gcol])^2) # surface temp at a given latitude (K)
    d0[gcol] = FT((ts / tt)^FT(4) - FT(1)) # optical depth
    nlev = nlay + 1

    #---bot_at_1------------------------------
    p_lev[1, gcol] = p0
    t_lev[1, gcol] = tt * (FT(1) + d0[gcol] * (p_lev[1, gcol] / p0)^α)^FT(0.25)
    z_lev[1, gcol] = FT(0)

    for ilay in 1:nlay
        #                if step == "linear"
        p_lev[ilay + 1, gcol] = p_lev[ilay, gcol] - Δp
        #                else
        #                    p_lev[ilay+1, gcol] = p_lev[ilay, gcol] * exp(-efac)
        #                end
        p_lay[ilay, gcol] = (p_lev[ilay, gcol] + p_lev[ilay + 1, gcol]) * FT(0.5)

        t_lev[ilay + 1, gcol] = tt * (FT(1) + d0[gcol] * (p_lev[ilay + 1, gcol] / p0)^α)^FT(0.25)
        t_lay[ilay, gcol] = tt * (FT(1) + d0[gcol] * (p_lay[ilay, gcol] / p0)^α)^FT(0.25)

        H = r_d * t_lay[ilay, gcol] / grav_
        Δz_lay = H * log(p_lev[ilay, gcol] / p_lev[ilay + 1, gcol])
        z_lev[ilay + 1, gcol] = Δz_lay + z_lev[ilay, gcol]
    end
    t_sfc[gcol] = t_lev[1, gcol]

    #---------------------------------
end
#-------------------------------------------------------------------------

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
function updatechildren(node::Node)
    return node
end

updatechildren(root::RootNode, node::Node) = RootNode(node)

### Operations ###

function updatechildren(indexsum::IndexSumOperation, A::AbstractTensor, indices::Indices)
    return contractioncheck(A)
end

function updatechildren(add::AddOperation, children::Vararg{Node})
    return +(children...)
end

function updatechildren(op::OuterProductOperation, A::Node, B::Node)
    return A⊗B
end

function updatechildren(indexing::IndexingOperation, A::AbstractTensor, indices::Indices)
    return A[indices.indices...]
end

function updatechildren(comp::ComponentTensorOperation, A::AbstractTensor, indices::Indices)
    return componenttensor(A, indices.indices...)
end

### Trigonometry ###

function updatechildren(si::SineOperation, A::AbstractTensor)
    return Base.sin(A)
end

function updatechildren(co::CosineOperation, A::AbstractTensor)
    return Base.cos(A)
end

function updatechildren(ta::TangentOperation, A::AbstractTensor)
    return Base.tan(A)
end

## Indices
function updatevectorspace(i::FreeIndex, V::AbstractVectorSpace)
    return FreeIndex(V, i.name, i.id)
end

function updatevectorspace(i::FixedIndex, V::AbstractVectorSpace)
    return FixedIndex(V, i.value)
end

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcol

"""
    'feffery'_antdcol(;kwargs...)
    'feffery'_antdcol(children::Any;kwargs...)
    'feffery'_antdcol(children_maker::Function;kwargs...)

An AntdCol component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `className` (String; optional)
- `flex` (String | Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `offset` (Real; optional)
- `order` (Real; optional)
- `pull` (Real; optional)
- `push` (Real; optional)
- `span` (Real; optional)
- `style` (Dict; optional)
"""
function 'feffery'_antdcol(; kwargs...)
        available_props = Symbol[:children, :id, :className, :flex, :loading_state, :offset, :order, :pull, :push, :span, :style]
        wild_props = Symbol[]
        return Component("'feffery'_antdcol", "AntdCol", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdcol(children::Any; kwargs...) = 'feffery'_antdcol(;kwargs..., children = children)
'feffery'_antdcol(children_maker::Function; kwargs...) = 'feffery'_antdcol(children_maker(); kwargs...)

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function text_wrap_tests(fname)
    Drawing(1400, 1400, fname)
    origin()
    setopacity(0.8)
    sethue("black")
    fontface("Impact")
    tiles = Tiler(1000, 1000, 4, 3, margin=50)
    for (pos, n) in tiles
        fsize = 12
        opacity = 1.0
        fontsize(fsize)
        @layer begin
            translate(pos)
            sethue("antiquewhite")
            box(O, tiles.tilewidth-2, tiles.tileheight-2, :fill)
            sethue("black")
            if iseven(n)
                textwrap("Luxor is a city in Upper (southern) Egypt. " ^ 12,
                tiles.tilewidth,
                O - (tiles.tilewidth/2 - 2, tiles.tileheight/2 - 2),
                rightgutter=20)
            else
                # textwrap(s::AbstractString, width::Real, pos::Point, linefunc::Function; rightgutter=5)
               textwrap("Luxor is a city in Upper (southern) Egypt. " ^ 12,
                    tiles.tilewidth,
                    O - (tiles.tilewidth/2 - 2,  tiles.tileheight/2 - 2),
                    # l is the text of a line, p is the position of the start of the line, h is the interline height
                    (n, l, p, h) -> begin
                            text(string(length(l)), p + (tiles.tilewidth - 20, 0))
                            opacity *= 0.8
                            setopacity(opacity)
                        end,
                    rightgutter=40)
            end
        end
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

text_wrap_tests("text-wrap-tests.svg")


text:
# ==================================
# squash
# ==================================

const CANDIDATE_SQUASH_TYPES = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64]

"""
    squash(g::ValGraph)
    squash(g::ValDiGraph)
    squash(g::ValOutDiGraph)

Return a copy of `g` with an `eltype` as small as possible.

This can help with performance. Only the `eltype` is changed, all other
types stay the same.
"""
function squash(g::Union{ValGraph, ValDiGraph, ValOutDiGraph})

    nvg = nv(g)
    # TODO find a more future proof way to extract the unparametrized type
    G = typeof(g).name.wrapper
    for V ∈ CANDIDATE_SQUASH_TYPES

        nvg < typemax(V) && return G{V}(g)
    end
end

# TODO implement destructive squash!

# ==================================
#    reverse
# ==================================

"""
    Reverse

Lazy wrapper for value graphs that reverses the directions of the graphs edges.
If this object is mutated, the wrapped graph is mutated as well.
"""
struct Reverse{V <: Integer, V_VALS, E_VALS, G_VALS, G<:AbstractValGraph{V, V_VALS, E_VALS, G_VALS}} <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}

    graph::G
end

"""
    reverse(g::AbstractValGraph)

Reverse the direction of the edges of `g`. Can also return `Reverse(g)` - a lazy
wrapper around `g`.
"""
reverse(g::AbstractValGraph) = is_directed(g) ? Reverse(g) : g

reverse(rg::Reverse) = rg.graph

## ---------------------------------
##   Reverse
## ---------------------------------

wrapped_graph(g::Reverse) = g.graph
wrapped_graph_type(::Type{<:Reverse{V, V_VALS, E_VALS, G_VALS, G}}) where {V, V_VALS, E_VALS, G_VALS, G} = G
@wrap_graph! Reverse include=[ne] exclude=[has_edge, add_edge!, rem_edge!, get_edgeval, set_edgeval!]

Base.zero(RG::Type{<:Reverse}) = Reverse(zero(wrapped_graph_type(RG)))

has_edge(rg::Reverse, s::Integer, d::Integer) =  has_edge(rg.graph, d, s)
add_edge!(rg::Reverse, s::Integer, d::Integer, values) = add_edge!(rg.graph, d, s, values)
rem_edge!(rg::Reverse, s::Integer, d::Integer) = rem_edge!(rg.graph, d, s)

get_edgeval(rg::Reverse, s::Integer, d::Integer, key::Integer) = get_edgeval(rg.graph, d, s, key)
get_edgeval(rg::Reverse, s::Integer, d::Integer, key::Symbol) = get_edgeval(rg.graph, d, s, key)
get_edgeval(rg::Reverse, s::Integer, d::Integer, ::Colon) = get_edgeval(rg.graph, d, s, :)

set_edgeval!(rg::Reverse, s::Integer, d::Integer, key::Integer, value) = set_edgeval!(rg.graph, d, s, key, value)
set_edgeval!(rg::Reverse, s::Integer, d::Integer, key::Symbol, value) = set_edgeval!(rg.graph, d, s, key, value)
set_edgeval!(rg::Reverse, s::Integer, d::Integer, ::Colon, values) = set_edgeval!(rg.graph, d, s, :, values)

outneighbors(rg::Reverse, u) = inneighbors(wrapped_graph(rg), u)
inneighbors(rg::Reverse, v) = outneighbors(wrapped_graph(rg), v)

outedgevals(rg::Reverse, u, key::Integer) = inedgevals(wrapped_graph(rg), u, key)
outedgevals(rg::Reverse, u, key::Symbol) = inedgevals(wrapped_graph(rg), u, key)
outedgevals(rg::Reverse, u, ::Colon) = inedgevals(wrapped_graph(rg), u, :)

inedgevals(rg::Reverse, v, key::Integer) = outedgevals(wrapped_graph(rg), v, key)
inedgevals(rg::Reverse, v, key::Symbol) = outedgevals(wrapped_graph(rg), v, key)
inedgevals(rg::Reverse, v, ::Colon) = outedgevals(wrapped_graph(rg), v, :)

# TODO implementing edges efficiently and with correct edge type


text:
include("header.jl")

if gpu() >= 0; @testset "cuarray" begin
    for nd in (1,2,3)
        sz = ntuple(i->8, nd)
        a0,b0 = rand(sz...),rand(sz...)
        a1,b1 = ka(a0),ka(b0)
        a2,b2 = Param(a0),Param(b0)
        a3,b3 = Param(a1),Param(b1)
        idx = ntuple(i->2:4, nd)
        @test getindex(a0,idx...) == getindex(a1,idx...)
        @test @gcheck getindex(a2,idx...)
        @test @gcheck getindex(a3,idx...)
        if nd == 1
            @test permutedims(a0) == permutedims(a1)
            @test @gcheck permutedims(a2)
            @test @gcheck permutedims(a3)
        elseif nd == 2
            @test permutedims(a0) == permutedims(a1)
            @test permutedims(a0,(2,1)) == permutedims(a1,(2,1))
            @test permutedims(a0,(1,2)) == permutedims(a1,(1,2))
            @test @gcheck permutedims(a2)
            @test @gcheck permutedims(a2,(2,1))
            @test @gcheck permutedims(a2,(1,2))
            @test @gcheck permutedims(a3)
            @test @gcheck permutedims(a3,(2,1))
            @test @gcheck permutedims(a3,(1,2))
        else
            @test permutedims(a0,(1,3,2)) == permutedims(a1,(1,3,2))
            @test @gcheck permutedims(a2,(1,3,2))
            @test @gcheck permutedims(a3,(1,3,2))
        end
        @test hcat(a0,b0) == hcat(a1,b1)
        @test vcat(a0,b0) == vcat(a1,b1)
        @test @gcheck hcat(a2,b2)
        @test @gcheck vcat(a2,b2)
        @test @gcheck hcat(a3,b3)
        @test @gcheck vcat(a3,b3)
        for i in 1:nd
            @test cat(a0,b0,dims=i) == cat(a1,b1,dims=i)
            @test @gcheck cat(a2,b2,dims=i)
            @test @gcheck cat(a3,b3,dims=i)
        end
        @test setindex!(a0,b0[idx...],idx...) == setindex!(a1,b1[idx...],idx...)

        # https://github.com/denizyuret/Knet.jl/issues/368
        @test argmax(a0) == argmax(a1)
        @test argmin(a0) == argmin(a1)
        @test findmax(a0) == findmax(a1)
        @test findmin(a0) == findmin(a1)
        for i in 1:nd
            @test argmax(a0,dims=i) == argmax(a1,dims=i)
            @test argmin(a0,dims=i) == argmin(a1,dims=i)
            @test findmax(a0,dims=i) == findmax(a1,dims=i)
            @test findmin(a0,dims=i) == findmin(a1,dims=i)
        end
    end
end; end


text:
using CPUTime

function compare_optimizers(functionOrProblem, parameters::Parameters = EMPTY_PARAMS;
    Methods = BlackBoxOptim.SingleObjectiveMethodNames, kwargs...)

    parameters = chain(convert(ParamsDict, parameters), kwargs2dict(kwargs))

    results = Any[]
    for m in Methods
        CPUtic()
        res = nothing
        try
            res = bboptimize(functionOrProblem, parameters; Method = m)
        catch e
            @warn "$e when running $m"
            push!(results, (m, e, nothing, missing, missing))
        end
        if res !== nothing
            push!(results, (m, "ok", best_candidate(res), best_fitness(res), CPUtoq()))
        end
    end

    sorted = sort( results, by = (t) -> t[4] )

    if get(parameters, :TraceMode, :compact) != :silent
        println("\n********************************************************************************")
        #println(describe(evaluator))
        for i in 1:length(sorted)
            println("$(i). $(sorted[i][1]): $(sorted[i][2]), fitness = $(sorted[i][4]), time = $(sorted[i][5])")
        end
        println("********************************************************************************\n")
    end

    return sorted
end

function compare_optimizers(problems::Dict{Any, OptimizationProblem},
                            parameters::Parameters = EMPTY_PARAMS;
                            Methods = BlackBoxOptim.SingleObjectiveMethodNames, kwargs...)

    parameters = chain(convert(ParamsDict, parameters), kwargs2dict(kwargs))

    # Lets create an array where we will save how the methods ranks per problem.
    ranks = zeros(length(Methods), length(problems))
    fitnesses = zeros(Float64, length(Methods), length(problems))
    times = zeros(Float64, length(Methods), length(problems))

    problems = collect(problems)

    for i in 1:length(problems)
        name, p = problems[i]
        res = compare_optimizers(p, parameters; Methods = Methods)
        for j in 1:length(res)
            method, best, fitness, elapsedtime = res[j]
            index = findfirst(Methods, method)
            ranks[index, i] = j
            fitnesses[index, i] = fitness
            times[index, i] = elapsedtime
        end
    end

    avg_ranks = round(mean(ranks, dims=2), digits=2)
    avg_fitness = round(mean(fitnesses, dims=2), digits=3)
    avg_times = round(mean(times, dims=2), digits=2)

    perm = sortperm(avg_ranks[:])
    println("\nBy avg rank:")
    for i in 1:length(methods)
        j = perm[i]
        print("\n$(i). $(methods[j]), avg rank = $(avg_ranks[j]), avg fitness = $(avg_fitness[j]), avg time = $(avg_times[j]), ranks = ")
        showcompact(ranks[j,:][:])
    end

    perm = sortperm(avg_fitness[:])
    println("\n\nBy avg fitness:")
    for i in 1:length(methods)
        j = perm[i]
        print("\n$(i). $(methods[j]), avg rank = $(avg_ranks[j]), avg fitness = $(avg_fitness[j]), avg time = $(avg_times[j]), ranks = ")
        showcompact(ranks[j,:][:])
    end

    return ranks, fitnesses
end

"""
Summarize a vector of float values by stating its mean, std dev and median.
"""
function report_on_values(desc, v, lpad = "", rpad = "", digits = 3)
    println("$(lpad)$(desc): $(round(mean(v), sigdigits=digits)) (std. dev = $(round(std(v), sigdigits=digits)), median = $(round(median(v), sigdigits=digits)))")
end

"""
Report the number of times each key was encountered in a count `dict`.

Returns a percentage dict calculated while iterating over the counted items.
"""
function count_dict_report(dict, desc, lpad = "", rpad = "")
    println(desc, ":")
    total = sum(collect(values(dict))) # FIXME collect() should not be required
    pdict = Dict()
    for (r, c) in dict
        pdict[r] = round(100.0*c/total, digits=2)
        println(lpad, r, ": ", c, " (", pdict[r], "%)", rpad)
    end
    pdict
end

"""
Print a report based on a result dict from one set of repeated runs of
an optimization method. Returns the success rate, i.e. number of times the
termination reason was "Within fitness tolerance...".
"""
function report_from_result_dict(statsdict)
    println("Method: $(statsdict[:method])")
    pdict = count_dict_report(statsdict[:reasoncounts], "  Termination reasons", "    ")
    report_on_values("Fitness", statsdict[:fitnesses], "  ")
    report_on_values("Time", statsdict[:times], "  ")
    report_on_values("Num function evals", statsdict[:numevals], "  ")

    success_rate = round(get(pdict, "Within fitness tolerance of optimum", 0.0), digits=3)
    println("  Success rate: ", success_rate, "%\n")
    success_rate
end

function rank_result_dicts_by(result_dicts, byfunc, desc; rev = false,
        descsummary = "mean", digits = 3, rpad = "")

    ranked = BlackBoxOptim.Utils.assign_ranks_within_tolerance(result_dicts; by = byfunc, tolerance = 1e-3, rev = rev)
    println("Ranked by $(descsummary) $(desc):")
    for (rank, rd, value) in ranked
        println("  $(rank). $(rd[:method]), $(round(value, sigdigits=digits))$(rpad)")
    end
    println("")
end

function report_on_methods_results_on_one_problem(problem, result_dicts, numrepeats, max_time, ftol)
    println("********************************************************************************\n")

    println("Problem: $(name(problem)), num dims = $(numdims(problem))")
    println("  Num repeats per method = ", numrepeats)
    println("  Fitness tolerance = ", ftol, " (a run is a success if it reaches to within this value of true optimum)")
    println("  Max time budget per run = ", max_time, " secs\n")

    rank_result_dicts_by(result_dicts, (d) -> d[:success_rate], "success rate (to reach within $(ftol) of optimum)";
            descsummary = "median", rev = true, rpad = "%")
    rank_result_dicts_by(result_dicts, (d) -> median(d[:fitnesses]), "fitness"; descsummary = "median")
    rank_result_dicts_by(result_dicts, (d) -> median(d[:times]), "time (in seconds)";
            descsummary = "median", rpad = " secs")
    rank_result_dicts_by(result_dicts, (d) -> round(Int, median(d[:numevals])), "num function evals"; descsummary = "median")

    for rd in result_dicts
        report_from_result_dict(rd)
    end
end

function repeated_bboptimize(numrepeats, problem, dim, methods, max_time, ftol = 1e-5, extraparams::Parameters = EMPTY_PARAMS)
    fp = instantiate(problem, dim)
    result_dicts = ParamsDict[]

    # Just so they are declared
    ps = best_so_far = nothing

    params = chain(extraparams, ParamsDict(:FitnessTolerance => ftol))

    for m in methods
        ts, fs, nes = zeros(numrepeats), zeros(numrepeats), zeros(Int, numrepeats)
        rcounts = Dict{String,Int}()

        results = map(1:numrepeats) do i
            result = bboptimize(fp, params; MaxTime = max_time, Method = m)

            # Save key data about this run:
            ts[i]  = elapsed_time(result)
            fs[i]  = best_fitness(result)
            nes[i] = f_calls(result)

            # And count only the general stop reasons
            reason = general_stop_reason(result)
            rcounts[reason] = 1 + get(rcounts, reason, 0)
        end

        rdict = ParamsDict(:method => m, :fitnesses => fs, :times => ts, :numevals => nes, :reasoncounts => rcounts)
        rdict[:success_rate] = report_from_result_dict(rdict)
        push!(result_dicts, rdict)
    end

    report_on_methods_results_on_one_problem(fp, result_dicts, numrepeats, max_time, ftol)
end


text:
mutable struct MipSolver
    inst::Instance
    loglevel::Int         # niveau de verbosité
    bestsol::Solution     # meilleure Solution rencontrée
    # Les attribut spécifiques au modèle
    mip_model::Model  # Le modèle MIP
    mip_x         # vecteur des variables d'atterrissage
    mip_b         # vecteur des variables de précédence (before)
    mip_cost      # variable du cout de la solution
    mip_costs     # variables du cout de chaque avion

    # Le constructeur
    function MipSolver(inst::Instance)
        ln1("MipSolver : constructeur avec $(Args.get("external_mip_solver"))")
        this=new()
        this.inst = inst
        this.loglevel = Args.get("level")
        this.bestsol = Solution(inst)
        # solver = Args.get("external_mip_solver") # NON UTILISÉ !
        # Création et configuration du modèle selon le solveur interne sélectionné
        this.mip_model = new_model(mode=:mip, log_level=this.loglevel)
        return this
    end
end

function buildSeqataModel(sv::MipSolver)
    # Quelques variables locales, pour abréger l'écriture
    planes = sv.inst.planes
    n = length(planes)
    model = sv.mip_model

    # Précalcul du coût de retard/avance pour chaque avion et chaque date
    # d'atterrissage.
    # La matrice des coûts est initialisée avec des coûts très élevés.
    # La vraie valeur de ces coûts sera calculée dans le domaine de validité de
    # la date d'atterrissage (dans l'intervale plane.lb:plane.ub)

    # Le big M (on a intérêt à l'optimiser !)
    # M = 4000
    #
    # Calcul de M minimum possible pour chaque couple de (i,j)
    # Si i est avant j, on veut (contrainte c2_sep) :
    #    xj >= xi + S_ij - M  =>  M >= xi - xj + S_ij
    # Donc on construit M[i,j] comme suit (+ 1 de sécurité)
    #    M := ub[i] - lb[j] + S_ij + 1
    M = Matrix{Int}(undef, sv.inst.nb_planes, sv.inst.nb_planes)
    for p1 in sv.inst.planes, p2 in sv.inst.planes
        M[p1.id,p2.id] = p1.ub - p2.lb + get_sep(sv.inst, p1, p2)
    end
    lg1() && @show maximum(M)

    # dates d'atterrissage effective
    @variable(model, x[1:n])
    # @variable(model, x[1:n], Int)
    sv.mip_x = x

    # coût retard/avance pour chaque avion
    @variable(model, costs[1:n])
    sv.mip_costs = costs

    # before : b[i,j] = 1 ssi i est avant j
    @variable(model, b[1:n,1:n], Bin)
    sv.mip_b = b

    # avec JuMP (contrairement à AMPL mais comme OPL) on peut définir des
    # expressions intermédiaires qui ne sont pas des variables de décision
    # (mais qui les utilisent).
    # Ces expressions pourront être exploitées dans le fonction objectif ou
    # pour faciliter l'affichage :
    @expression(model, total_cost, sum(costs))
    sv.mip_cost = total_cost

    #
    # Objectif Minimiser le retard total
    #
    @objective(model, Min, sv.mip_cost)

    #
    # Les contraintes
    #

    # Contrainte sur les bornes de la date d'atterrissage de chaque avion
    # c1_bounds = Array{VariableRef}(undef, n)
    # for i in 1:n
    #     p = planes[i]
    #     # # Création d'un contrainte anonyme (car le nom est inutile
    #     # # et création multiple => warning sous julia0.6)
    #     # # @constraint(model, c1_bounds[i], p.lb <= x[i] <= p.ub)
    #     # #
    #     # # Ces doubles bound mal sont encore mal supportées par les solveurs
    #     # # cplex, glpk et clp (même avec julia-0.6.1 et julia-0.7)
    #     # #     @constraint(model, p.lb <= x[i] <= p.ub)
    #     # @constraint(model, p.lb <= x[i]) # name: c1_bounds_lb[i]
    #     # @constraint(model, x[i] <= p.ub) # name: c1_bounds_ub[i]

    #     @constraint(model, c1_bounds[i], p.lb <= x[i] <= p.ub)
    # end
    @constraint(model, c1_bounds[p in planes], p.lb <= x[p.id] <= p.ub)

    for p1 in planes, p2 in planes
        if p1==p2 continue end
        # ATTENTION ,  POUR L'INSTANT ON DOIT FAIRE M*xxx et pas xxx*M
        # Mais ce sera possible à partir de la version de julia-0.4
        # @constraint(model, c2_sep[p1.id,p2.id],
        #     x[p2.id] >= x[p1.id] + get_sep(s.inst, p1, p2) - M*(1-b[p1.id, p2.id]))
        @constraint(model,
                    x[p2.id] >= x[p1.id] + get_sep(sv.inst, p1, p2)
                                - M[p1.id, p2.id]*(1-b[p1.id, p2.id]))

        # @constraint(model, c3_before_is_exclusive[p1.id, p2.id],
        #     b[p1.id, p2.id] + b[p2.id, p1.id] == 1)
        @constraint(model,
                    b[p1.id, p2.id] + b[p2.id, p1.id] == 1)
    end

    for p in planes
        # Principe: pour chaque segment de l'avion (paire successive de
        # timecosts), on ajoute la contrainte suivante :
        #    Le coût de l'avion est >= la droite associée au segment
        # L'équation de la droite du segment est de la forme :
        #  c = slope*(t-t1) + c1
        #  c = (c2-c1)/(t2-t1) * (t-t1) + c1
        if p.lb==p.ub
            # MAIS... certains avions dégénérés n'ont qu'un seul timecost avec
            # avec date d'atterrissage fiées : p.lb==p.target==p.ub
            # On les traite à part en imposant leur variable
            # date x[p.id] et leur costs[p.id]
            @constraint(model, costs[p.id] == 0.0)
            @constraint(model, x[p.id] == p.target)
        else
            for j = 1:(length(p.timecosts)-1)
                t1 = p.timecosts[j][1]
                c1 = p.timecosts[j][2]
                t2 = p.timecosts[j+1][1]
                c2 = p.timecosts[j+1][2]
                slope = (c2-c1)/(t2-t1)
                cst=@constraint(model, costs[p.id] >= c1 + slope*(x[p.id]-t1))
                # @show t1, c1, t2, c2, slope
                # @show cst
            end
        end
    end

end

# Résoud le problème complet : calcul l'ordre, le timing (dates d'atterrissage)
# et le cout total de la solution optimale au problème SEQATA (coûts linéaires)
# Cette fonction crée le modèle PLNE, le résoud puis met l'objet solution à
# jour.
#
function solve(sv::MipSolver)
    ln1("BEGIN solve(MipSolver)")

    ln1( "="^60 )
    ln1("Construction du modèle\n")
    buildSeqataModel(sv)

    ln1( "="^60 )
    ln1("Lancement de la résolution\n")

    # On prefixe la méthode solve par le nom du module pour éviter tout conflit avec
    # mes propres méthodes solve()
    # @show sv.mip_model
    if lg1()
        @ms optimize!(sv.mip_model)
    else
        optimize!(sv.mip_model)
    end

    ln1( "="^60 )
    ln1("Exploitation des résultats\n")
    if  JuMP.termination_status(sv.mip_model) != MOI.OPTIMAL
        print("ERREUR : pas de solution pour :\n    ")
        ln1("JuMP.termination_status(sv.mip_model)=", JuMP.termination_status(sv.mip_model))
        ln1( to_s(sv.bestsol) )
        exit(1)
    end

    # Extraction des valeurs entières des variables d'atterrissage
    xvals = round.(Int, value.(sv.mip_x))

    # Il reste maintenant à mettre à jour notre objet solution à partir du
    # résultat de la résolution MIP
    #
    mip_costs = value.(sv.mip_costs)
    mip_cost = round(value(sv.mip_cost), digits=Args.get(:cost_precision))

    # copy!(sv.bestsol, sv.inst.planes, xvals, mip_costs, mip_cost) # ORI
    update_from!(sv.bestsol, xvals, mip_costs, mip_cost)

    # On trie juste la solution par date d'atterrissage croissante des avions
    # pour améliorer la présentation de la solution
    sort!(sv.bestsol)

    ln1("END solve(MipSolver)")
end
# END TYPE DescentSolver


text:
############################################################################################
# Checks
############################################################################################

@testset "Core - construct_flatten" begin
    df = FlattenDefault()
    x = [1.]
    constraint = Unconstrained()

    construct_flatten([1.])
    construct_flatten(Unconstrained(), [1.])
    construct_flatten(df, x)
    construct_flatten(df, constraint, x)
    construct_flatten(df, UnflattenStrict(), x)
    construct_flatten(df, UnflattenStrict(), constraint, x)
end

@testset "Core - Param" begin
    Param(_RNG, .1, Normal())
    param = Param(0.1, Normal())
end

@testset "Core - Checkfinite" begin
    @test _checkfinite(.1)
    @test !_checkfinite(Inf)
    @test _checkfinite( (.1, .2, [.2, 3.], zeros(2,3)))

    @test !_checkfinite([.1, -Inf, 3.])
    @test !_checkfinite([[2. 3. ; 4. 5], [.1, -Inf, 3.]])

end

@testset "Core - Checkprior" begin
    @test _checkprior(Distributions.Normal())
    @test !_checkprior(Inf)
    @test _checkprior([Distributions.Normal(), Distributions.Normal()])
    @test !_checkprior([Inf, Distributions.Normal(), Distributions.Normal()])
    @test !_checkprior([[Distributions.Normal(), Distributions.Normal()], [Inf, Distributions.Normal(), Distributions.Normal()]])
    @test !_checkprior(
        [[[Distributions.Normal(), Distributions.Normal()], [Distributions.Normal(), Distributions.Normal()]],
        [Inf, Distributions.Normal(), Distributions.Normal()]]
    )
    @test _checkprior([
        [[Distributions.Normal(), Distributions.Normal()], [Distributions.Normal(), Distributions.Normal()]],
        [Distributions.Normal(), Distributions.Normal()]
    ])
    @test _checkprior([Distributions.Normal() Distributions.Normal() ; Distributions.Normal() Distributions.Normal()])
    @test !_checkprior([Distributions.Normal() Distributions.Normal() ; Fixed() Distributions.Normal()])
end

@testset "Core - Checksampleable" begin
    @test _checksampleable(Distributions.Normal())
    @test !_checksampleable(Fixed())
    @test _checksampleable([Distributions.Normal(), Distributions.Normal()])
    @test !_checksampleable([Fixed(), Distributions.Normal(), Distributions.Normal()])
    @test !_checksampleable([[Distributions.Normal(), Distributions.Normal()], [Fixed(), Distributions.Normal(), Distributions.Normal()]])
    @test !_checksampleable([
    [[Distributions.Normal(), Distributions.Normal()], [Distributions.Normal(), Distributions.Normal()]],
    [Fixed(), Distributions.Normal(), Distributions.Normal()]
    ])
    @test _checksampleable([
    [[Distributions.Normal(), Distributions.Normal()], [Distributions.Normal(), Distributions.Normal()]],
    [Distributions.Normal(), Distributions.Normal()]
    ])
    @test _checksampleable([Distributions.Normal() Distributions.Normal() ; Distributions.Normal() Distributions.Normal()])
    @test !_checksampleable([Distributions.Normal() Distributions.Normal() ; Fixed() Distributions.Normal()])
end

@testset "Core - Checkparams" begin
    @test _checkparams(Param(.1, Normal()))
    @test !_checkparams(Fixed())
end

@testset "Core - _check" begin
    @test ModelWrappers._check(_RNG, Normal(), .1)
    @test !ModelWrappers._check(_RNG, MvNormal(LinearAlgebra.Diagonal(map(abs2, [1., 2.]))), .1)
end

############################################################################################
# Random
############################################################################################
@testset "Core - sample_constraint" begin
    @test ModelWrappers.sample_constraint(_RNG, Fixed(), nothing) isa Nothing
    @test ModelWrappers.sample_constraint(_RNG, Unconstrained(), 1.0) == 1.0
    @test ModelWrappers.sample_constraint(_RNG, Normal(), 3.) isa Float64
end

@testset "Core - log_prior" begin
    @test ModelWrappers.log_prior(Fixed(), nothing) == 0.0
    @test ModelWrappers.log_prior([Normal(), Normal()], 0.0) == 0.0
    @test ModelWrappers.log_prior([Normal(), Normal()], [0.,0.]) isa Float64
    @test ModelWrappers.log_prior(Normal(), 2.) isa Float64
end

@testset "Core - log_prior_with_transform" begin
    @test ModelWrappers.log_prior_with_transform(Fixed(), nothing) == 0.0
    @test ModelWrappers.log_prior_with_transform([Normal(), Normal()], 0.0) == 0.0
    @test ModelWrappers.log_prior_with_transform([Normal(), Normal()], [0.,0.]) isa Float64
    @test ModelWrappers.log_prior_with_transform(Normal(), 2.) isa Float64
end

@testset "Core - log_abs_det_jac" begin
end

############################################################################################
# Utility
############################################################################################
@testset "Core - _get_constraint" begin
    @test _get_constraint(Param(1., Normal())) isa Distribution
    @test _get_constraint(Param(1., Fixed())) isa AbstractConstraint
    @test _get_constraint(Param(1., Unconstrained())) isa AbstractConstraint
end

@testset "Core - _get_val" begin
    @test _get_val(Param(1., Normal())) isa Float64
    @test _get_val(Param(1., Fixed())) isa Float64
    @test _get_val(Param(1., Unconstrained())) isa Float64
end

@testset "Core - _allparam" begin
end

@testset "Core - _anyparam" begin
end

@testset "Core - _paramnames" begin
    @test length(ModelWrappers._paramnames(:a, 1)) == 1
    @test length(ModelWrappers._paramnames(:a, 2)) == 2
    _names = ModelWrappers._paramnames((:a, :b), (1,2))
    @test _names == ["a", "b1", "b2"]
end

@testset "Core - paramnames" begin
    _names2 = ModelWrappers.paramnames(
        FlattenDefault(),
        (a = Unconstrained(), c = Fixed(), b = [Normal(), Normal()]),
        (a = 1., c = 2., b = [3., 4.]),
    )
    @test _names2 == ["a", "b1", "b2"]
    _names3 = ModelWrappers.paramnames(
        (:a,:c,:b),
        FlattenDefault(),
        (a = Unconstrained(), c = Fixed(), b = [Normal(), Normal()]),
        (a = 1., c = 2., b = [3., 4.]),
    )
    @test all(_names2 .== _names3)
end

@testset "Core - paramcount" begin
    _vals = ModelWrappers.paramcount(
        FlattenDefault(),
        (a = 1., c = 2., b = [3., 4.])
    )
    @test _vals == (1,1,2)
    _vals2 = ModelWrappers.paramcount(
        (:a,:c,:b),
        FlattenDefault(),
        (a = 1., c = 2., b = [3., 4.])
    )
    @test all(_vals .== _vals2)
end


text:
# Autogenerated wrapper script for MagicVLSI_jll for x86_64-linux-gnu
export tclmagic

using Tcl_jll
using GLU_jll
using Cairo_jll
using Tk_jll
JLLWrappers.@generate_wrapper_header("MagicVLSI")
JLLWrappers.@declare_library_product(tclmagic, "tclmagic.so")
function __init__()
    JLLWrappers.@generate_init_header(Tcl_jll, GLU_jll, Cairo_jll, Tk_jll)
    JLLWrappers.@init_library_product(
        tclmagic,
        "lib/magic/tcl/tclmagic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
points(S::Jacobi, n) = points(Chebyshev(domain(S)), n)

struct JacobiTransformPlan{T,CPLAN,CJT} <: AbstractTransformPlan{T}
    chebplan::CPLAN
    cjtplan::CJT
end

JacobiTransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
    JacobiTransformPlan{eltype(chebplan),CPLAN,CJT}(chebplan, cjtplan)

plan_transform(S::Jacobi, v::AbstractVector) =
    JacobiTransformPlan(plan_transform(Chebyshev(), v), plan_cheb2jac(v, S.a, S.b))
*(P::JacobiTransformPlan, vals::AbstractVector) = P.cjtplan*(P.chebplan*vals)

struct JacobiITransformPlan{T,CPLAN,CJT} <: AbstractTransformPlan{T}
    ichebplan::CPLAN
    icjtplan::CJT
end

JacobiITransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
    JacobiITransformPlan{eltype(chebplan),CPLAN,CJT}(chebplan, cjtplan)

plan_itransform(S::Jacobi, v::AbstractVector) =
    JacobiITransformPlan(plan_itransform(Chebyshev(), v), plan_jac2cheb(v, S.a, S.b))
*(P::JacobiITransformPlan, cfs::AbstractVector) = P.ichebplan*(P.icjtplan*cfs)

function coefficients(f::AbstractVector{T}, a::Jacobi, b::Chebyshev) where T
    if domain(a) == domain(b) && (!isapproxinteger(a.a-0.5) || !isapproxinteger(a.b-0.5))
        jac2cheb(f, convert(T,a.a), convert(T,a.b))
    else
        defaultcoefficients(f,a,b)
    end
end
function coefficients(f::AbstractVector{T}, a::Chebyshev, b::Jacobi) where T
    isempty(f) && return f
    if domain(a) == domain(b) && (!isapproxinteger(b.a-0.5) || !isapproxinteger(b.b-0.5))
        cheb2jac(f, convert(T,b.a), convert(T,b.b))
    else
        defaultcoefficients(f,a,b)
    end
end

function coefficients(f::AbstractVector,a::Jacobi,b::Jacobi)
    if domain(a) == domain(b) && (!isapproxinteger(a.a-b.a) || !isapproxinteger(a.b-b.b))
        jac2jac(f,a.a,a.b,b.a,b.b)
    else
        defaultcoefficients(f,a,b)
    end
end


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

# ONE-SAMPLE FINITE-SAMPLE ADJUSTMENT

function adjfactor!(V::Matrix, obj::ParModel, corr::Homoscedastic)
end

function adjfactor!(V::Matrix, obj::ParModel, corr::Heteroscedastic)
    (corr.adj == true) && (n = nobs(obj) ; lmul!(n / (n - 1), V))
end

function adjfactor!(V::Matrix, obj::ParModel, corr::Clustered)
    (corr.adj == true) && (n = corr.nc ; lmul!(n / (n - 1), V))
end

function adjfactor!(V::Matrix, obj::ParModel, corr::CrossCorrelated)
    (corr.adj == true) && (n = nobs(obj) ; lmul!(n / (n - 1), V))
end

#==========================================================================================#

# TWO-SAMPLE FINITE-SAMPLE ADJUSTMENT

function adjfactor!(V::Matrix, obj₁::ParModel, obj₂::ParModel, corr::Heteroscedastic)
    if corr.adj == true

        nonmissing₁ = getnonmissing(obj₁)
        nonmissing₂ = getnonmissing(obj₂)
        n           = sum(nonmissing₁ .* nonmissing₂)

        lmul!(n / (n - 1), V)
    end
end

function adjfactor!(V::Matrix, obj₁::ParModel, obj₂::ParModel, corr::Clustered)
    if corr.adj == true

        corr₁ = getcorr(obj₁)
        corr₂ = getcorr(obj₂)
        touse = corr.nonmissing .* corr₁.nonmissing .* corr₂.nonmissing
        ic    = corr.ic[touse]
        n     = length(unique(ic))

        lmul!(n / (n - 1), V)
    end
end

function adjfactor!(V::Matrix, obj₁::ParModel, obj₂::ParModel, corr::CrossCorrelated)
    if corr.adj == true

        corr₁ = getcorr(obj₁)
        corr₂ = getcorr(obj₂)
        touse = corr.nonmissing .* corr₁.nonmissing .* corr₂.nonmissing
        n     = sum(touse)

        lmul!(n / (n - 1), V)
    end
end


text:
vfname        = "SampleFiles/Restricted/10081701.WVP"
sac_fstr      = "SampleFiles/Restricted/20081701.sac-0*"
eq_wvm1_pat   = "SampleFiles/SUDS/*1991191072247.wvm1.sac"
eq_wvm1_file  = "SampleFiles/SUDS/eq_wvm1.sud"
eq_wvm2_pat   = "SampleFiles/SUDS/*1991191072247.wvm2.sac"
eq_wvm2_file  = "SampleFiles/SUDS/eq_wvm2.sud"
lsm_pat       = "SampleFiles/SUDS/*1992187065409.sac"
lsm_file      = "SampleFiles/SUDS/lsm.sud"
rot_pat       = "SampleFiles/SUDS/*1993258220247.sac"
rot_file      = "SampleFiles/SUDS/rotate.sud"

printstyled("  SUDS\n", color=:light_green)

printstyled("    readsudsevt\n", color=:light_green)
# read into SeisEvent
W = SUDS.readsudsevt(rot_file)
W = SUDS.readsudsevt(lsm_file)
W = SUDS.readsudsevt(eq_wvm1_file)
W = SUDS.readsudsevt(eq_wvm2_file)

printstyled("    in read_data\n", color=:light_green)
if safe_isfile(vfname)
  redirect_stdout(out) do
    SUDS.suds_support()
    S = verified_read_data("suds", vfname, v=3, full=true)
    S = verified_read_data("suds", eq_wvm1_file, v=3, full=true)
    S = verified_read_data("suds", eq_wvm2_file, v=3, full=true)
    S = verified_read_data("suds", lsm_file, v=3, full=true)
    S = verified_read_data("suds", rot_file, v=3, full=true)
    S = verified_read_data("suds", "SampleFiles/SUDS/eq_wvm*sud", v=3, full=true)
  end
end

printstyled("    equivalence to SAC readsuds\n", color=:light_green)
# Volcano-seismic event supplied by W. McCausland
if has_restricted
  S = SUDS.read_suds(vfname, full=true)
  S2 = verified_read_data("sac", sac_fstr, full=true)
  for n = 1:S2.n
    id = S2.id[n]
    if startswith(id, ".")
      id = "OV"*id
    end
    i = findid(id, S)
    if i > 0
      @test isapprox(S.x[i], S2.x[n])
      @test isapprox(Float32(S2.fs[n]), Float32(S.fs[i]))
      @test abs(S.t[i][1,2] - S2.t[i][1,2]) < 2000 # SAC only has ~1 ms precision
    end
  end
else
  printstyled("      skipped volcano-seismic test file\n", color=:red)
end

# SUDS sample files
# from eq_wvm1.sud
S = SUDS.read_suds(eq_wvm1_file)
S1 = verified_read_data("sac", eq_wvm1_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from eq_wvm2.sud
S = SUDS.read_suds(eq_wvm2_file)
S1 = verified_read_data("sac", eq_wvm2_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from lsm.sud
S = SeisData()
SUDS.read_suds!(S, lsm_file)
S1 = verified_read_data("sac", lsm_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

# from rot.sud
S = SUDS.read_suds(rot_file)
S1 = verified_read_data("sac", rot_pat)
for n = 1:S1.n
  i = findid(S1.id[n], S)
  if i > 0
    @test isapprox(S.x[i], S1.x[n])
  else
    @warn(string(S1.id[n], " not found in S; check id conversion!"))
  end
end

printstyled("    data types not yet seen in real samples\n", color=:light_green)
f_out = string("fake_suds.sud")
id_buf = codeunits("UWSPORTSBALL")
nb = Int32(3200)
n_unreadable = Int32(60)
suds_types = (Int32, Complex{Float32}, Float64, Float32, Complex{Float64})

# skeleton: 0x53, 0x00, struct_id, struct_size, nbytes_following_struct
suds_struct_tag = (0x53, 0x00, zero(Int16), zero(Int32), zero(Int32))
suds_unreadable = (0x53, 0x36, Int16(4), n_unreadable, zero(Int32))
suds_5 = (0x53, 0x36, Int16(5), Int32(76), zero(Int32))
suds_readable = (0x53, 0x36, Int16(7), Int32(62), nb)

for (j, data_code) in enumerate([0x32, 0x63, 0x64, 0x66, 0x74])
  T = suds_types[j]

  x = rand(T, div(nb, sizeof(T)))
  io = open(f_out, "w")

  # Unreadable packet
  [write(io, i) for i in suds_unreadable]
  write(io, rand(UInt8, n_unreadable))

  # Packet 5
  [write(io, i) for i in suds_5]
  p = position(io)
  write(io, id_buf)
  write(io, rand(Int16, 2))
  write(io, 45.55)
  write(io, 122.62)
  write(io, 77.1f0)
  write(io, rand(UInt8, 7))
  write(io, 0x76)
  write(io, data_code)
  write(io, rand(UInt8, 3))
  write(io, 1.152f8, 0.0f0, 0.0f0)
  write(io, Int16[2, 32])
  skip(io, 4)
  write(io, 0.0f0, -32767.0f0)

  # Packet 7
  [write(io, i) for i in suds_readable]
  p = position(io)
  write(io, id_buf)
  write(io, d2u(now()))
  write(io, zero(Int16))
  write(io, data_code)
  write(io, 0x00)
  skip(io, 4)
  write(io, Int32(div(nb, sizeof(T))))
  write(io, 50.0f0)
  skip(io, 16)
  write(io, zero(Float64))
  write(io, 0.0f0)
  write(io, x)
  close(io)

  if j < 5
    S = verified_read_data("suds", f_out)
    if data_code == 0x63
      @test isapprox(S.x[1], real(x))
    else
      @test isapprox(S.x[1], x)
    end
  else
    @test_throws ErrorException read_data("suds", f_out)
  end
end
safe_rm(f_out)


text:
@testset "promotion" begin
    @polyvar x y z

    @test typeof(@inferred promote(1, y)) == NTuple{2, Term{Int, Monomial{(y,), 1}}}
    @test typeof(@inferred promote(1, y^2)) == NTuple{2, Term{Int, Monomial{(y,), 1}}}
    @test typeof(@inferred promote(1, 2y)) == NTuple{2, Term{Int, Monomial{(y,), 1}}}
    @test typeof(@inferred promote(1, Polynomial(2y))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(y,), 1}}, Vector{Term{Int, Monomial{(y,), 1}}}}}
    @test typeof(@inferred promote(1, Polynomial([2y]))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(y,), 1}}, Vector{Term{Int, Monomial{(y,), 1}}}}}

    @test typeof(@inferred promote(x, 1)) == NTuple{2, Term{Int, Monomial{(x,), 1}}}
    @test typeof(@inferred promote(x, y)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(y, x)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(x, y^2)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(x, 2y)) == NTuple{2, Term{Int, Monomial{(x, y), 2}}}
    @test typeof(@inferred promote(x, Polynomial(2y))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(x, Polynomial([2y]))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}

    @test typeof(@inferred promote(x^2, 1)) == NTuple{2, Term{Int, Monomial{(x,), 1}}}
    @test typeof(@inferred promote(x^2, y)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(y, x^2)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(x^2, y^2)) == NTuple{2, Monomial{(x, y), 2}}
    @test typeof(@inferred promote(x^2, 2y)) == NTuple{2, Term{Int, Monomial{(x, y), 2}}}
    @test typeof(@inferred promote(x^2, Polynomial(2y))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(x^2, Polynomial([2y]))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}

    @test typeof(@inferred promote(2x, 1)) == NTuple{2, Term{Int, Monomial{(x,), 1}}}
    @test typeof(@inferred promote(2x, y)) == NTuple{2, Term{Int, Monomial{(x, y), 2}}}
    @test typeof(@inferred promote(2x, y^2)) == NTuple{2, Term{Int, Monomial{(x, y), 2}}}
    @test typeof(@inferred promote(2x, 2y)) == NTuple{2, Term{Int, Monomial{(x, y), 2}}}
    @test typeof(@inferred promote(2x, Polynomial(2y))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(2x, Polynomial([2y]))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}

    @test typeof(@inferred promote(Polynomial(2x), 1)) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x,), 1}}, Vector{Term{Int, Monomial{(x,), 1}}}}}
    @test typeof(@inferred promote(Polynomial(2x), y)) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(Polynomial(2x), y^2)) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(Polynomial(2x), 2y)) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(Polynomial(2x), Polynomial(2y))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}
    @test typeof(@inferred promote(Polynomial(2x), Polynomial([2y]))) == NTuple{2, Polynomial{Int, Term{Int, Monomial{(x, y), 2}}, Vector{Term{Int, Monomial{(x, y), 2}}}}}

    @test typeof(@inferred(promote(Monomial{tuple(), 0}(), x))) == NTuple{2, Monomial{(x,), 1}}
    @test typeof(@inferred(promote(x, Monomial{tuple(), 0}()))) == NTuple{2, Monomial{(x,), 1}}
end

@testset "conversion" begin
    @polyvar x y z

    @test convert(Monomial{(x, y)}, x) == Monomial{(x, y)}((1, 0))
    @test_throws InexactError convert(Monomial{(x,)}, y)
    @test_throws InexactError convert(Monomial{(x,)}, x * y)
end

@testset "variable_union_type" begin
    @polyvar x y z

    function test(p)
        @test MultivariatePolynomials.variable_union_type(p) == TypedPolynomials.Variable
        @test MultivariatePolynomials.variable_union_type(typeof(p)) == TypedPolynomials.Variable
    end
    test(x)
    test(y)
    test(z)
    test(x*y*z)
    test(2x)
    test(2y)
    test(2z)
    test(x + y + z)
end


text:
export gradiente

function gradiente(
  nlp::AbstractNLPModel;
  atol = 1e-6,
  rtol = 1e-6,
  max_iter = 1000,
  max_time = 10.0
)
  x = copy(nlp.meta.x0)

  f(x) = obj(nlp, x)
  ∇f(x) = grad(nlp, x)

  ϵ = atol + rtol * norm(∇f(x))
  t₀ = time()

  iter = 0
  Δt = time() - t₀
  resolvido = norm(∇f(x)) < ϵ
  cansado = iter > max_iter || Δt > max_time

  status = :unknown

  α = 1.0
  while !(resolvido || cansado)

    x -= α * ∇f(x)
    α *= 0.99

    iter += 1
    Δt = time() - t₀
    resolvido = norm(∇f(x)) < ϵ
    cansado = iter > max_iter || Δt > max_time
  end

  if resolvido
    status = :first_order
  elseif cansado
    if iter > max_iter
      status = :max_iter
    elseif Δt > max_time
      status = :max_time
    end
  end

  return GenericExecutionStats(
    status,
    nlp,
    solution=x,
    objective=f(x),
    dual_feas=norm(∇f(x)),
    elapsed_time=Δt,
    iter=iter
  )
end


text:
export armijo
export goldenLS

"""
armijo(f,fk,dfk,xk,pk)

Backtracked Armijo linesearch
"""
function armijo(f::Function,fk,dfk,xk,pk;maxIter=10, c1=1e-4,b=0.5)
LS = 1
t  = 1
while LS<=maxIter
    if f(xk+t*pk) <= fk + t*c1*dot(dfk,pk)
        break
    end
    t *= b
    LS += 1
end
if LS>maxIter
	LS= -1
	t = 0.0
end
return t,LS
end

"""
goldenLS(f,xk,pk,a,b)

Golden section linesearch
"""
function goldenLS(f::Function,xk,pk,a,b;atol=1e-2,maxIter=40)
	ratio = (sqrt(5)-1)/2.0
	# make sure points we have an interval
	a,b   = (a<b)? (a,b) :  (b,a)
	# evaluate f at left point
	al    = b - ratio*(b-a)
	fl    = f(xk + al*pk)
	# evaluate f at right point
	ar    = a + ratio*(b-a)
	fr    = f(xk + ar*pk)

	LS    = 0
	while (LS<maxIter && abs(ar-al)>atol)
		if fl<fr # continue search in [a,ar]
			b  = ar
			ar = al
			fr = fl
			al = b - ratio*(b-a)
			fl = f(xk + al*pk)
		else # continue search in [al,b]
			a  = al
			al = ar
			fl = fr
			ar = a + ratio*(b-a)
			fr = f(xk + ar*pk)
		end
		# @printf "iter=%d, [al,ar]=[%1.2f,%1.2f], [fl,fr]=[%1.2e,%1.2e]\n" LS al ar fl fr
		LS += 1
	end
	return (al+ar)/2,LS
end


text:
@testset "Gauss–Lobatto" begin
    # Check error
    @test_throws DomainError gausslobatto(-1)
    @test_throws DomainError gausslobatto(0)
    @test_throws DomainError gausslobatto(1)

    n = 2
    x,w = gausslobatto(n)
    @test x[1] ≈ -1.
    @test x[2] ≈ 1.
    @test w[1] ≈ 1.
    @test w[2] ≈ 1.

    n = 3
    x,w = gausslobatto(n)
    @test x[1] ≈ -1.
    @test abs(x[2])<1e-15
    @test x[3] ≈ 1.
    @test w[1] ≈ 1/3
    @test w[2] ≈ 4/3
    @test w[3] ≈ 1/3

    tol = 1e-14
    n = 42
    x,w = gausslobatto(n)
    @test ( length(x) == n) && ( length(w) == n )
    @test abs(x[37] - 0.922259214258616) < tol
    @test abs(w[37] - 0.029306411216166) < tol
    @test ( x[1] == -1 && x[n] == 1 )

    @test dot( w,(x.^2)) ≈ 2/3
    @test dot( w,exp.(x)) ≈ exp(1)-exp(-1)
end


text:
# Perform some tests on Step1EndGrid

using Test
include("../src/Step1EndGrid.jl")

@testset "Step1EndGrid.jl" begin

    # Setting
    β = 0.9932
    eH = 1.0
    eL = 0.1
    E = [eL; eH]
    # Transition matrix: Pr(j|i) = Π(i,j) (columns sum to 1)
    πHH = 0.925 # e_H | e_H
    πHL = 0.5   # e_H | e_L
    Π = [1 - πHL  πHL
        1 - πHH πHH]
    q = 1.3
    σ = 1.5
    a_lb = -2.0
    a_max = 5
    a_size = 100
    A = range(a_lb, a_max, length = a_size)
    A = [A; ]
    Aprime = repeat(A, 1, e_size)

    # Test 1: check β makes sense
    β = 1.1
    @test_throws AssertionError EndGridSearch(A, E, Π, β, σ, q)

    # Test 2: check transition matrix
    β = 0.9932
    Π = [0.3 0.5
        0.5 0.5]
    @test_throws AssertionError EndGridSearch(A, E, Π, β, σ, q)

    # Test 3: check that the price makes sense
    q = 0.8
    @test_throws AssertionError EndGridSearch(A, E, Π, β, σ, q)

    # Test 4: check that σ is appropriate
    σ = 0.9
    @test_throws AssertionError EndGridSearch(A, E, Π, β, σ, q)

end


text:
# Instant arithmetic
for op in (:+,:*,:%,:/)
    @eval ($op)(x::Instant,y::Instant) = throw(ArgumentError("Operation not defined for Instants"))
end
(+)(x::Instant) = x
(-){T<:Instant}(x::T,y::T) = x.periods - y.periods

# TimeType arithmetic
for op in (:+,:*,:%,:/)
    @eval ($op)(x::TimeType,y::TimeType) = throw(ArgumentError("Operation not defined for TimeTypes"))
end
(+)(x::TimeType) = x
(-){T<:TimeType}(x::T,y::T) = x.instant - y.instant

# TimeType-Year arithmetic
function (+)(dt::DateTime,y::Year)
    oy,m,d = yearmonthday(dt); ny = oy+value(y); ld = daysinmonth(ny,m)
    return DateTime(ny,m,d <= ld ? d : ld,hour(dt),minute(dt),second(dt),millisecond(dt))
end
function (+)(dt::Date,y::Year)
    oy,m,d = yearmonthday(dt); ny = oy+value(y); ld = daysinmonth(ny,m)
    return Date(ny,m,d <= ld ? d : ld)
end
function (-)(dt::DateTime,y::Year)
    oy,m,d = yearmonthday(dt); ny = oy-value(y); ld = daysinmonth(ny,m)
    return DateTime(ny,m,d <= ld ? d : ld,hour(dt),minute(dt),second(dt),millisecond(dt))
end
function (-)(dt::Date,y::Year)
    oy,m,d = yearmonthday(dt); ny = oy-value(y); ld = daysinmonth(ny,m)
    return Date(ny,m,d <= ld ? d : ld)
end

# TimeType-Month arithmetic
# monthwrap adds two months with wraparound behavior (i.e. 12 + 1 == 1)
monthwrap(m1,m2) = (v = mod1(m1+m2,12); return v < 0 ? 12 + v : v)
# yearwrap takes a starting year/month and a month to add and returns
# the resulting year with wraparound behavior (i.e. 2000-12 + 1 == 2001)
yearwrap(y,m1,m2) = y + fld(m1 + m2 - 1,12)

function (+)(dt::DateTime,z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y,m,value(z))
    mm = monthwrap(m,value(z)); ld = daysinmonth(ny,mm)
    return DateTime(ny,mm,d <= ld ? d : ld,hour(dt),minute(dt),second(dt),millisecond(dt))
end
function (+)(dt::Date,z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y,m,value(z))
    mm = monthwrap(m,value(z)); ld = daysinmonth(ny,mm)
    return Date(ny,mm,d <= ld ? d : ld)
end
function (-)(dt::DateTime,z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y,m,-value(z))
    mm = monthwrap(m,-value(z)); ld = daysinmonth(ny,mm)
    return DateTime(ny,mm,d <= ld ? d : ld,hour(dt),minute(dt),second(dt),millisecond(dt))
end
function (-)(dt::Date,z::Month)
    y,m,d = yearmonthday(dt)
    ny = yearwrap(y,m,-value(z))
    mm = monthwrap(m,-value(z)); ld = daysinmonth(ny,mm)
    return Date(ny,mm,d <= ld ? d : ld)
end
(+)(x::Date,y::Week) = return Date(UTD(value(x) + 7*value(y)))
(-)(x::Date,y::Week) = return Date(UTD(value(x) - 7*value(y)))
(+)(x::Date,y::Day)  = return Date(UTD(value(x) + y))
(-)(x::Date,y::Day)  = return Date(UTD(value(x) - y))
(+)(x::DateTime,y::Period)   = return DateTime(UTM(value(x)+toms(y)))
(-)(x::DateTime,y::Period)   = return DateTime(UTM(value(x)-toms(y)))
(+)(y::Period,x::TimeType) = x + y
(-)(y::Period,x::TimeType) = x - y

(.+){T<:TimeType}(x::AbstractArray{T}, y::Period) = reshape(T[i + y for i in x], size(x))
(.-){T<:TimeType}(x::AbstractArray{T}, y::Period) = reshape(T[i - y for i in x], size(x))
(.+){T<:TimeType}(y::Period, x::AbstractArray{T}) = x .+ y
(.-){T<:TimeType}(y::Period, x::AbstractArray{T}) = x .- y


text:
include( joinpath( repPath, "misc.jl" ) )
include( joinpath( repPath, "nodeflux.jl" ) )
include( joinpath( repPath, "nodepop.jl" ) )
include( joinpath( repPath, "transitionflux.jl" ) )
include( joinpath( repPath, "nodecomposition.jl" ) )
include( joinpath( repPath, "careerprogression.jl" ) )
include( joinpath( repPath, "subpopulationpop.jl" ) )
include( joinpath( repPath, "subpopulationage.jl" ) )
include( joinpath( repPath, "initpop.jl" ) )


text:
using Ipopt
using Test

@testset "C API" begin
    # First of all, test that hs071 example works
    include("hs071_test.jl")

    # Test setting some options
    # String option
    println("\nString option")
    addOption(prob, "hessian_approximation", "exact")
    @test_throws ErrorException addOption(prob, "hessian_approximation", "badoption")
    println("\nInt option")
    # Int option
    addOption(prob, "file_print_level", 3) == nothing
    @test_throws ErrorException addOption(prob, "file_print_level", -1)
    # Double option
    println("\nFloat option")
    addOption(prob, "derivative_test_tol", 0.5)
    @test_throws ErrorException addOption(prob, "derivative_test_tol", -1.0)

    # Test opening an output file
    openOutputFile(prob, "blah.txt", 5)
    Ipopt.freeProblem(prob) # Needed before the `rm` on Windows.
    # unlink the output file
    rm("blah.txt")

    # Test that the ipopt binary works
    @test Ipopt.amplexefun("-v") == 0
end

@testset "MathProgBase" begin
    include("MPB_wrapper.jl")
end

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end


text:
module Benitor

function fibonacciTail(n::Int,a,b)
    if n <= 0
        return b
    end
    return fibonacciTail(n-1 , b, a+b)
end

end # module


text:
using InfiniteWarpedMixtureModels
using Documenter

DocMeta.setdocmeta!(InfiniteWarpedMixtureModels, :DocTestSetup, :(using InfiniteWarpedMixtureModels); recursive=true)

makedocs(;
    modules=[InfiniteWarpedMixtureModels],
    authors="Shuichi Miyazawa",
    repo="https://github.com/shu13830/InfiniteWarpedMixtureModels.jl/blob/{commit}{path}#{line}",
    sitename="InfiniteWarpedMixtureModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://shu13830.github.io/InfiniteWarpedMixtureModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/shu13830/InfiniteWarpedMixtureModels.jl",
    devbranch="main",
)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdbadge

"""
    'feffery'_antdbadge(;kwargs...)
    'feffery'_antdbadge(children::Any;kwargs...)
    'feffery'_antdbadge(children_maker::Function;kwargs...)

An AntdBadge component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `className` (String; optional)
- `color` (String; optional)
- `count` (Real; optional)
- `dot` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `nClicks` (Real; optional)
- `offset` (Array of Reals; optional)
- `overflowCount` (Real; optional)
- `showZero` (Bool; optional)
- `size` (a value equal to: 'default', 'small'; optional)
- `status` (a value equal to: 'success', 'processing', 'default', 'error', 'warning'; optional)
- `style` (Dict; optional)
- `text` (String; optional)
- `title` (String; optional)
"""
function 'feffery'_antdbadge(; kwargs...)
        available_props = Symbol[:children, :id, :className, :color, :count, :dot, :loading_state, :nClicks, :offset, :overflowCount, :showZero, :size, :status, :style, :text, :title]
        wild_props = Symbol[]
        return Component("'feffery'_antdbadge", "AntdBadge", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdbadge(children::Any; kwargs...) = 'feffery'_antdbadge(;kwargs..., children = children)
'feffery'_antdbadge(children_maker::Function; kwargs...) = 'feffery'_antdbadge(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
using ITensors
using ITensors.ITensorNetworkMaps
using KrylovKit
using LinearAlgebra

χ = 3
d = 2
l = Index(χ, "l")
s = Index(d, "s")

l0 = addtags(l, "c=0")
l1 = addtags(l, "c=1")
A = randomITensor(l0, l1, s)
A′ = prime(dag(A); inds=(l0, l1))

T = ITensorNetworkMap([A, A′]; input_inds=(l1, l1'), output_inds=(l0, l0'))
v = randomITensor(input_inds(T))
@show T(v)
@show T * v
@show (I - T)(v)
@show (3I - 2T)(v)

D, V = eigsolve(T, v)
@show (T - D[1] * I)(V[1])


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# write commands to gnuplot's pipe
function gnuplot_send(s::AbstractString)
    gin = gnuplot_state.fid[1] # gnuplot STDIN
    w = write(gin, string(s,"\n"))
    # check that data was accepted by the pipe
    if !(w > 0)
        println("Something went wrong writing to gnuplot STDIN.")
        return
    end
    flush(gin)
end