text:
# types.jl - Type definitions and structure
import Base: show, summary, convert, getindex, eltype

# Tipo abstracto para definir contenedores del IPC
"""
    abstract type AbstractCPIBase{T <: AbstractFloat}

Tipo abstracto para representar conjuntos de colecciones de datos del IPC.

Vea también: [`FullCPIBase`](@ref), [`VarCPIBase`](@ref) e [`IndexCPIBase`](@ref).
"""
abstract type AbstractCPIBase{T <: AbstractFloat} end

# Tipos para los vectores de fechas
"""
    const DATETYPE = Union{StepRange{Date, Month}, Vector{Date}}
Tipos posibles para el campo de fechas `dates` de un [`AbstractCPIBase`](@ref).
"""
const DATETYPE = Union{StepRange{Date, Month}, Vector{Date}}

# Tipos para los vectores de códigos y nombres
"""
    const DESCTYPE = Union{Vector{String}, Nothing}
Tipos posibles para los nombres en el campo `names` de un [`FullCPIBase`](@ref).
"""
const DESCTYPE = Union{Vector{String}, Nothing}

"""
    const CODETYPE = Union{Vector{String}, Nothing}
Tipos posibles para los códigos en el campo `codes` de un [`FullCPIBase`](@ref).
"""
const CODETYPE = Union{Vector{String}, Nothing}

# El tipo B representa el tipo utilizado para almacenar los índices base.
# Puede ser un tipo flotante, por ejemplo, Float64 o bien, si los datos
# disponibles empiezan con índices diferentes a 100, un vector, Vector{Float64},
# por ejemplo

"""
    FullCPIBase{T, B} <: AbstractCPIBase{T}

    FullCPIBase(ipc::Matrix{T}, v::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B, codes::CODETYPE, names::DESCTYPE) where {T, B}
    FullCPIBase(df::DataFrame, gb::DataFrame)

Contenedor completo para datos desagregados del IPC de un país. Se representa
por:
- Matriz de índices de precios `ipc` que incluye la fila con los índices del
  número base.
- Matriz de variaciones intermensuales `v`. En las filas contiene los períodos y
  en las columnas contiene los gastos básicos.
- Vector de ponderaciones `w` de los gastos básicos.
- Fechas correspondientes `dates` (por meses).
- Índices base `baseindex`.
- Códigos y nombres de los gastos básicos en `codes` y `names`.

El tipo `T` representa el tipo de datos para representar los valores de punto
flotante. El tipo `B` representa el tipo del campo `baseindex`; por ejemplo,
`Float32` o `Vector{Float32}`.
"""
Base.@kwdef struct FullCPIBase{T, B} <: AbstractCPIBase{T}
    ipc::Matrix{T}
    v::Matrix{T}
    w::Vector{T}
    dates::DATETYPE
    baseindex::B
    codes::CODETYPE
    names::DESCTYPE

    function FullCPIBase(ipc::Matrix{T}, v::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B, codes::CODETYPE=nothing, names::DESCTYPE=nothing) where {T, B}
        @debug "Sizes:" size(ipc) size(v) length(dates) first(dates) last(dates)
        size(ipc, 2) == size(v, 2) || throw(ArgumentError("número de columnas debe coincidir entre matriz de índices y variaciones"))
        size(ipc, 2) == length(w) || throw(ArgumentError("número de columnas debe coincidir con vector de ponderaciones"))
        size(ipc, 1) == size(v, 1) == length(dates) || throw(ArgumentError("número de filas de `ipc` debe coincidir con vector de fechas"))
        new{T, B}(ipc, v, w, dates, baseindex, codes, names)
    end
end

"""
    IndexCPIBase{T, B} <: AbstractCPIBase{T}

    IndexCPIBase(ipc::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B) where {T, B}

Contenedor genérico de índices de precios del IPC de un país. Se representa por:
- Matriz de índices de precios `ipc` que incluye la fila con los índices del númbero base.
- Vector de ponderaciones `w` de los gastos básicos.
- Fechas correspondientes `dates` (por meses).
- Índices base `baseindex`.

El tipo `T` representa el tipo de datos para representar los valores de punto
flotante. El tipo `B` representa el tipo del campo `baseindex`; por ejemplo,
`Float32` o `Vector{Float32}`.
"""
Base.@kwdef struct IndexCPIBase{T, B} <: AbstractCPIBase{T}
    ipc::Matrix{T}
    w::Vector{T}
    dates::DATETYPE
    baseindex::B

    function IndexCPIBase(ipc::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B) where {T, B}
        size(ipc, 2) == length(w) || throw(ArgumentError("número de columnas debe coincidir con vector de ponderaciones"))
        size(ipc, 1) == length(dates) || throw(ArgumentError("número de filas debe coincidir con vector de fechas"))
        new{T, B}(ipc, w, dates, baseindex)
    end
end

"""
    VarCPIBase{T, B} <: AbstractCPIBase{T}

    VarCPIBase(v::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B) where {T, B}

Contenedor genérico para de variaciones intermensuales de índices de precios del
IPC de un país. Se representa por:
- Matriz de variaciones intermensuales `v`. En las filas contiene los períodos y
  en las columnas contiene los gastos básicos.
- Vector de ponderaciones `w` de los gastos básicos.
- Fechas correspondientes `dates` (por meses).
- Índices base `baseindex`.

Este tipo es el utilizado en el contenedor de bases del IPC de un país,
denominado [`CountryStructure`](@ref), ya que con los datos de un
`VarCPIBase` es suficiente para computar cualquier medida de inflación basada en
índices de precios individuales o en un estadístico de resumen de las
variaciones intermensuales (como un percentil, o una media truncada).

El tipo `T` representa el tipo de datos para representar los valores de punto
flotante. El tipo `B` representa el tipo del campo `baseindex`; por ejemplo,
`Float32` o `Vector{Float32}`.

Ver también: [`CountryStructure`](@ref), [`UniformCountryStructure`](@ref),
[`MixedCountryStructure`](@ref)
"""
Base.@kwdef struct VarCPIBase{T, B} <: AbstractCPIBase{T}
    v::Matrix{T}
    w::Vector{T}
    dates::DATETYPE
    baseindex::B

    function VarCPIBase(v::Matrix{T}, w::Vector{T}, dates::DATETYPE, baseindex::B) where {T, B}
        size(v, 2) == length(w) || throw(ArgumentError("número de columnas debe coincidir con vector de ponderaciones"))
        size(v, 1) == length(dates) || throw(ArgumentError("número de filas debe coincidir con vector de fechas"))
        new{T, B}(v, w, dates, baseindex)
    end
end

## Constructores
# Los constructores entre tipos crean copias y asignan nueva memoria

function _getbaseindex(baseindex)
    if length(unique(baseindex)) == 1
        return baseindex[1]
    end
    baseindex
end

"""
    FullCPIBase(df::DataFrame, gb::DataFrame)

Este método constructor devuelve una estructura `FullCPIBase` a partir de los
DataFrames de índices de precios `df` y de descripción de los gastos básicos
`gb`.
- El DataFrame `df` posee la siguiente estructura:
    - Contiene en la primera columna las fechas o períodos de los datos. En las
      siguientes columnas, debe contener los códigos de cada una de las
      categorías o gastos básicos de la estructura del IPC, junto con la serie
      de tiempo con los índices de precios individuales.
    - En las filas del DataFrame contiene los períodos por meses. La primera
      fila del DataFrame se utiliza para obtener el índice base. Si el valor es
      el mismo para todos los gastos básicos, se tomará únicamente este valor
      escalar (por ejemplo 100.0 como un Float64). En algunos casos, es posible
      que no se disponga de la información completa, por lo que los índices base
      podrían ser diferentes entre sí. En este caso, `baseindex` almacenará el
      vector de índices base originales.
    - Un ejemplo de cómo puede verse este DataFrame es el siguiente:
```
121×219 DataFrame
 Row │ Date         _011111  _011121  _011131  _011141  _011142  _011151  _011152 ⋯
     │ Date         Float64  Float64  Float64  Float64  Float64  Float64  Float64 ⋯
─────┼─────────────────────────────────────────────────────────────────────────────
   1 │ 2000-12-01   100.0    100.0    100.0    100.0    100.0    100.0    100.0   ⋯
   2 │ 2001-01-01   100.55   103.23   101.66   106.47   100.36   100.0    102.57
   3 │ 2001-02-01   101.47   104.82   102.73   108.38   101.37   100.0    103.35
   4 │ 2001-03-01   101.44   107.74   104.9    103.76   101.32   100.0    104.27
   5 │ 2001-04-01   101.91   107.28   106.19   107.83   101.82   100.0    104.73  ⋯
   6 │ 2001-05-01   102.77   106.12   106.9    109.16   101.81   100.0    105.21
   7 │ 2001-06-01   103.23   109.04   107.4    112.13   102.72   100.0    105.47
   8 │ 2001-07-01   104.35   112.72   107.96   117.19   105.09   100.0    105.66
  ⋮  │     ⋮          ⋮        ⋮        ⋮        ⋮        ⋮        ⋮        ⋮     ⋱
 114 │ 2010-05-01   218.45   501.39   200.28   477.5    179.0    215.0    164.16  ⋯
 115 │ 2010-06-01   219.28   503.35   203.88   476.26   180.94   214.02   164.97
 116 │ 2010-07-01   219.1    503.78   205.19   478.34   181.78   217.6    165.9
 117 │ 2010-08-01   218.52   507.45   206.87   486.72   181.51   223.76   166.46
 118 │ 2010-09-01   218.9    505.8    206.45   501.23   182.04   228.34   166.04  ⋯
 119 │ 2010-10-01   219.51   504.41   205.78   504.4    182.35   221.98   166.3
 120 │ 2010-11-01   219.11   509.63   205.41   502.88   182.16   217.01   166.34
 121 │ 2010-12-01   218.79   511.38   205.09   506.04   182.14   218.63   165.99
                                                   211 columns and 105 rows omitted
```

- El DataFrame `gb` posee la siguiente estructura:
    - La primera columna contiene los códigos de las columnas del DataFrame
      `df`.
    - La segunda columna contiene el nombre o la descripción de cada una de las
      categorías en las columnas de `df`.
    - Y finalmente, la tercer columna, debe contener las ponderaciones asociadas
      a cada una de las categorías o gastos básicos de las columnas de `df`.
    - Los nombres de las columnas no son tomados en cuenta, solamente el orden y
      los tipos.
    - Un ejemplo de cómo puede verse este DataFrame es el siguiente:
```
218×3 DataFrame
 Row │ Code     GoodOrService                       Weight
     │ String   String                              Float64
─────┼────────────────────────────────────────────────────────
   1 │ _011111  Arroz                               0.483952
   2 │ _011121  Pan                                 2.82638
   3 │ _011131  Pastas frescas y secas              0.341395
   4 │ _011141  Productos de tortillería            1.69133
  ⋮  │   ⋮                     ⋮                     ⋮
 216 │ _093111  Gastos por seguros                  0.236691
 217 │ _093121  Gastos por servicios funerarios     0.289885
 218 │ _094111  Gastos por servicios diversos pa…   0.151793
                                              211 rows omitted
```
"""
function FullCPIBase(df::DataFrame, gb::DataFrame)
    # Obtener matriz de índices de precios
    ipc_mat = Matrix(df[!, 2:end])
    # Matrices de variaciones intermensuales de índices de precios
    v_mat = 100 .* (ipc_mat[2:end, :] ./ ipc_mat[1:end-1, :] .- 1)
    # Ponderación de gastos básicos o categorías
    w = gb[!, 3]
    # Actualización de fechas
    dates = df[2, 1]:Month(1):df[end, 1]
    # Revisión de códigos
    codes_df = names(df)[2:end]
    codes_gb = convert.(String, gb[:, 1])
    (codes_df == codes_gb) || throw(AssertionError("Códigos en las columnas de `df` deben coincidir con códigos en las filas de `gb`"))
    # Nombres
    names_gb = gb[:, 2]
    # Estructura de variaciones intermensuales de base del IPC
    return FullCPIBase(ipc_mat[2:end, :], v_mat, w, dates, _getbaseindex(ipc_mat[1, :]), codes_gb, names_gb)
end

"""
    VarCPIBase(df::DataFrame, gb::DataFrame)

Este constructor devuelve una estructura `VarCPIBase` a partir del DataFrame
de índices de precios `df`, que contiene en las columnas las categorías o gastos
básicos del IPC y en las filas los períodos por meses. Las ponderaciones se obtienen
de la estructura `gb`, en la tercera columna de ponderaciones.

Para conocer la estructura de los DataFrames necesarios, vea también: [`FullCPIBase`](@ref).
"""
function VarCPIBase(df::DataFrame, gb::DataFrame)
    # Obtener estructura completa
    cpi_base = FullCPIBase(df, gb)
    # Estructura de variaciones intermensuales de base del IPC
    VarCPIBase(cpi_base)
end

function VarCPIBase(base::FullCPIBase)
    nbase = deepcopy(base)
    VarCPIBase(nbase.v, nbase.w, nbase.dates, nbase.baseindex)
end

# Obtener VarCPIBase de IndexCPIBase con variaciones intermensuales
VarCPIBase(base::IndexCPIBase) = convert(VarCPIBase, deepcopy(base))

"""
    IndexCPIBase(df::DataFrame, gb::DataFrame)

Este constructor devuelve una estructura `IndexCPIBase` a partir del DataFrame
de índices de precios `df`, que contiene en las columnas las categorías o gastos
básicos del IPC y en las filas los períodos por meses. Las ponderaciones se obtienen
de la estructura `gb`, en la tercera columna de ponderaciones.

Para conocer la estructura de los DataFrames necesarios, vea también: [`FullCPIBase`](@ref).
"""
function IndexCPIBase(df::DataFrame, gb::DataFrame)
    # Obtener estructura completa
    cpi_base = FullCPIBase(df, gb)
    # Estructura de índices de precios de base del IPC
    return IndexCPIBase(cpi_base)
end

function IndexCPIBase(base::FullCPIBase)
    nbase = deepcopy(base)
    IndexCPIBase(nbase.ipc, nbase.w, nbase.dates, nbase.baseindex)
end

# Obtener IndexCPIBase de VarCPIBase con capitalización intermensual
IndexCPIBase(base::VarCPIBase) = convert(IndexCPIBase, deepcopy(base))

## Conversión

# Estos métodos sí crean copias a través de la función `convert` de los campos
convert(::Type{T}, base::VarCPIBase) where {T <: AbstractFloat} =
    VarCPIBase(convert.(T, base.v), convert.(T, base.w), base.dates, convert.(T, base.baseindex))
convert(::Type{T}, base::IndexCPIBase) where {T <: AbstractFloat} =
    IndexCPIBase(convert.(T, base.ipc), convert.(T, base.w), base.dates, convert.(T, base.baseindex))
convert(::Type{T}, base::FullCPIBase) where {T <: AbstractFloat} =
    FullCPIBase(convert.(T, base.ipc), convert.(T, base.v), convert.(T, base.w), base.dates, convert.(T, base.baseindex), base.codes, base.names)

# OJO:
# Estos métodos no crean copias, como se indica en la documentación:
# > If T is a collection type and x a collection, the result of convert(T, x)
# > may alias all or part of x.
# Al convertir de esta forma se muta la matriz de variaciones intermensuales y se
# devuelve el mismo tipo, pero sin asignar nueva memoria

function convert(::Type{IndexCPIBase}, base::VarCPIBase)
    vmat = base.v
    capitalize!(vmat, base.baseindex)
    IndexCPIBase(vmat, base.w, base.dates, base.baseindex)
end

function convert(::Type{VarCPIBase}, base::IndexCPIBase)
    ipcmat = base.ipc
    varinterm!(ipcmat, base.baseindex)
    VarCPIBase(ipcmat, base.w, base.dates, base.baseindex)
end

# Función de conversión de FullCPIBase -> VarCPIBase, no crea copia de los datos
# subyacentes.
function convert(::Type{VarCPIBase}, base::FullCPIBase)
    VarCPIBase(base.v, base.w, base.dates, base.baseindex)
end

# Tipo de flotante del contenedor
eltype(::AbstractCPIBase{T}) where {T} = T

## Métodos para mostrar los tipos

function _formatdate(fecha)
    Dates.format(fecha, dateformat"u-yy")
end

function summary(io::IO, base::AbstractCPIBase)
    field = hasproperty(base, :v) ? :v : :ipc
    periods, ngoods = size(getproperty(base, field))
    print(io, typeof(base), ": ", periods, " períodos × ", ngoods, " gastos básicos ")
    datestart, dateend = _formatdate.((first(base.dates), last(base.dates)))
    print(io, datestart, "-", dateend)
    # print(io, "└─→ ", propertynames(base))
end

function show(io::IO, base::Union{VarCPIBase, IndexCPIBase})
    summary(io, base)
    println(io)
    field = hasproperty(base, :v) ? :v : :ipc
    pretty_table(io, getproperty(base, field);
        cell_first_line_only = true,
        row_names = base.dates,
        show_row_number = true,
        row_name_column_title = "Dates",
        header = (1:length(base.w), base.w),
        crop = :both,
        vcrop_mode = :middle,
        formatters = ft_printf("%0.4f")
    )
end

function show(io::IO, base::FullCPIBase)
    summary(io, base)
    println(io)
    pretty_table(io, base.ipc;
        row_names = base.dates,
        show_row_number = true,
        row_name_column_title = "Dates",
        crop_subheader = true,
        header = (base.codes, base.names, base.w),
        crop = :both,
        vcrop_mode = :middle,
        formatters = ft_printf("%0.2f"),
    )
end

"""
    periods(base::VarCPIBase)

Computa el número de períodos (meses) en las base de variaciones intermensuales.
"""
periods(base::VarCPIBase) = size(base.v, 1)


text:
using MAT, Plots

@views function generate_grf3D(lx=100.0, ly=100.0, lz=100.0, sf=1.0, cl=(10.0, 8.0, 5.0), k_m=100.0, nh=10000, nx=64, ny=64, nz=64;
                               cov_typ="expon", do_reset=true, do_viz=true, do_save=false)
    # Derived numerics
    dx, dy, dz  = lx/nx, ly/ny, lz/nz  # numerical grid step size
    me = 0
    co1, co2, co3 = 0, 0, 0
    # Array allocation
    Yf      = @zeros(nx, ny, nz)
    # Visu init
    if do_viz
        y_sl    = Int(ceil(ny/2))
        X, Y, Z = -lx/2:dx:lx/2, -ly/2:dy:ly/2, -lz/2:dz:lz/2
    end

    if cov_typ=="expon"
        # Generate the 3D exponential covariance function
        wtime_it = grf3D_expon!(Yf, sf, cl, nh, nx, ny, nz, dx, dy, dz; me, co1, co2, co3, do_reset)
    elseif cov_typ=="gauss"
        # Generate the 3D Gaussian covariance function
        wtime_it = grf3D_gauss!(Yf, sf, cl[1]/2.0, nh, k_m, nx, ny, nz, dx, dy, dz; me, co1, co2, co3, do_reset)
    else
        error("Undefined covariance function.")
    end

    # Performance
    A_eff    = 2/1e9*nx*ny*nz*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    T_eff    = A_eff/wtime_it                      # Effective memory throughput [GB/s]
    if (me==0) @printf("T_eff = %1.2f GB/s \n", round(T_eff, sigdigits=2)) end

    # Visualisation
    if do_viz
        display(heatmap(X, Z, Array(Yf)[:,y_sl,:]', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Z[1],Z[end]), c=:hot, title="3D RandomField (y-slice)"))
    end
    if do_save  file = matopen("grf3D_$(cov_typ).mat", "w"); write(file, "grf3D", Array(Yf)); close(file)  end

    return Yf
end

ext: jl


text:
dv = DataVector[1, 2, NA]
dm = DataArray([1 2; 3 4])
dt = datazeros(2, 2, 2)

df = DataFrame(quote
                 A = 1:2
                 B = 2:3
               end)

for el in dv
    @assert ndims(el) == 0
end

for el in dm
    @assert ndims(el) == 0
end

for el in dt
    @assert ndims(el) == 0
end

for row in EachRow(df)
    @assert isa(row, DataFrame)
end

for col in EachCol(df)
    @assert isa(col, AbstractDataVector)
end

@assert isequal(map(x -> min(matrix(x)), EachRow(df)), {1,2})
@assert isequal(map(min, EachCol(df)), DataFrame(quote A = 1; B = 2 end))

# @test_fail for x in df; end # Raises an error


text:
using Debugger

function main()

    a = 5
    println(fact(a))

end

function rec_sum(x)

    # base case
    if x == 0
        return x
    else
        # recursive call
        return x + rec_sum(x - 1)
    end

end

function fact(x)

    println(x)

    if x <= 1
        return 1
    else
        return x * fact(x - 1)
    end

end

function fib(x)

    if x < 2
        return x
    else
        return fib(x - 1) + fib(x - 2)
    end

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


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.Preprocess
using JuliaFEM.Testing

@testset "read ascii mesh" begin
    mesh = """
    COOR_2D
    N1          0.0 0.0
    N2          1.0 0.0
    N3          1.0 1.0
    N4          0.0 1.0
    FINSF

    QUAD4
    E1         N1 N2 N3 N4
    FINSF

    SEG2
    E2         N3 N4
    FINSF

    GROUP_NO NOM=NALL
    N1 N2
    FINSF

    GROUP_MA NOM=BODY1
    E1 E2
    FINSF

    FIN
    """

#   m = parse(mesh, Val{:CODE_ASTER_MAIL})
#   @test m["nodes"]["N1"] == [0.0, 0.0]
#   @test m["elements"]["E1"] == ["QUAD4", ["N1", "N2", "N3", "N4"]]
#   @test m["elements"]["E2"] == ["SEG2", ["N3", "N4"]]
#   @test m["elsets"]["BODY1"] == ["E1", "E2"]
#   @test m["nsets"]["NALL"] == ["N1", "N2"]
end

@testset "parse nodes" begin
    section = """
    N9  2.0 3.0 4.0
    COOR_3D
    N1          0.0 0.0 0.0
    N2          1.0 0.0 0.0
    N3          1.0 1.0 0.0
    N4          0.0 1.0 0.0
    N5          0.0 0.0 1.0
    N6          1.0 0.0 1.0
    N7          1.0 1.0 1.0
    N8          0.0 1.0 1.0
    FINSF
    absdflasdf
    N12 3.0 4.0 5.0 6.0
    N13 3.0 4.0 5.0
    """
    nodes = aster_parse_nodes(section)
    @test nodes[1] == Float64[0.0, 0.0, 0.0]
    @test nodes[8] == Float64[0.0, 1.0, 1.0]
    @test length(nodes) == 8
end

function JuliaFEM.get_mesh(::Type{Val{Symbol("block_2d_1elem_quad4")}})
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_1elem_quad4.med"
    mesh = aster_read_mesh(fn)
    return mesh
end

@testset "test reading aster .med file" begin
    mesh = get_mesh("block_2d_1elem_quad4")
    #=
    info("nodes")
    for (k, v) in mesh.nodes
        info("$k => $v")
    end
    info("node sets")
    for (k, v) in mesh.node_sets
        info("$k => $v")
    end
    info("elements")
    for (k, v) in mesh.elements
        info("$k => $v, type = $(mesh.element_types[k])")
    end
    info("element sets")
    for (k, v) in mesh.element_sets
        info("$k => $v")
    end
    =#
    @test length(mesh.element_sets) == 5
    @test length(mesh.node_sets) == 4
    @test length(mesh.elements) == 5
    @test length(mesh.nodes) == 4
    for elset in [:BLOCK, :TOP, :BOTTOM, :LEFT, :RIGHT]
        @test haskey(mesh.element_sets, elset)
        @test length(mesh.element_sets[elset]) == 1
    end
    for nset in [:TOP_LEFT, :TOP_RIGHT, :BOTTOM_LEFT, :BOTTOM_RIGHT]
        @test haskey(mesh.node_sets, nset)
        @test length(mesh.node_sets[nset]) == 1
    end
end

@testset "test filter by element set" begin
    mesh = get_mesh("block_2d_1elem_quad4")
    mesh2 = filter_by_element_set(mesh, :BLOCK)
    @test haskey(mesh2.element_sets, :BLOCK)
    @test length(mesh2.elements) == 1
end

function calculate_volume(mesh_name, eltype)
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(fn, mesh_name)
    elements = create_elements(mesh; element_type=eltype)
    V = 0.0
    time = 0.0
    for element in elements
        for ip in get_integration_points(element)
            detJ = element(ip, time, Val{:detJ})
            detJ > 0 || warn("negative determinant for element $eltype !")
            V += ip.weight*detJ
        end
    end
    info("volume of $eltype is $V")
    return V
end

@testset "calculate volume for 1 element models" begin
    @test isapprox(calculate_volume("TRIANGLE_TRI3_1", :Tri3), 1/2)
    @test isapprox(calculate_volume("TRIANGLE_TRI6_1", :Tri6), 1/2)
    @test isapprox(calculate_volume("TRIANGLE_TRI7_1", :Tri7), 1/2)
    @test isapprox(calculate_volume("SQUARE_QUAD4_1", :Quad4), 2^2)
    @test isapprox(calculate_volume("SQUARE_QUAD8_1", :Quad8), 2^2)
    @test isapprox(calculate_volume("SQUARE_QUAD9_1", :Quad9), 2^2)
    @test isapprox(calculate_volume("TETRA_TET4_1", :Tet4), 1/6)
    @test isapprox(calculate_volume("TETRA_TET10_1", :Tet10), 1/6)
#   @test isapprox(calculate_volume("TETRA_TET14_1", :Tet14), 1/6)
    @test isapprox(calculate_volume("CUBE_HEX8_1", :Hex8), 2^3)
    @test isapprox(calculate_volume("CUBE_HEX20_1", :Hex20), 2^3)
    @test isapprox(calculate_volume("CUBE_HEX27_1", :Hex27), 2^3)
    @test isapprox(calculate_volume("WEDGE_WEDGE6_1", :Wedge6), 1)
#   @test isapprox(calculate_volume("WEDGE_WEDGE15_1", :Wedge15, 1/2))
#   @test isapprox(calculate_volume("PYRAMID_PYRAMID5_1", :Pyramid5, ?))
#   @test isapprox(calculate_volume("PYRAMID_PYRAMID13_1", :Pyramid13, ?))
end

@testset "get nodal field from aster file" begin
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
    medfile = JuliaFEM.Preprocess.RMEDFile(fn)
    temp = JuliaFEM.Preprocess.aster_read_data(medfile, "TEMP")
    info("temp = $temp")
    # more like functional testing, results are what they are,
    # we're happy to just have some results
    @test true
end


text:
function __oddify(n_eval::Int)
    if n_eval % 2 == 0
        n_eval += 1
    end
    return n_eval
end

function __calculate_gamma(n_eval::Int)
    gamma_k = zeros(n_eval-1)
    gamma_k[1] = 1/2
    for kk = 2:(n_eval-1)
        val = kk + 1
        for ii = 2:(kk-1)
            nom = ii
            den = ii + kk
            val *= (nom / den)
        end
        val *= 1/2
        gamma_k[kk] = val
    end
    return gamma_k
end

function __select_noise_sigma(sigma_k, first_idx)
    noise_sigma = 0.0
    if first_idx != 0
        noise_sigma = sigma_k[first_idx]
    end
    return noise_sigma
end

function __selection_criterion(T_ij, sigma_k, eta = 4.0, window_size = 3)
    n_eval = length(sigma_k) + 1
    first_idx_sigma = 0
    for ii = 1:(n_eval-1-window_size+1)

        min_s = min( sigma_k[ii:(ii+window_size-1)]... )
        max_s = max( sigma_k[ii:(ii+window_size-1)]... )

        if max_s <= eta * min_s
            first_idx_sigma = ii
            break
        end
    end

    first_idx_sign = 0
    break_flag = false
    for ii = 2:(n_eval-2)

        sign_1 = sign( T_ij[1, ii+1] )
        for jj = 2:(n_eval - ii)
            sign_2 = sign( T_ij[jj, ii+1] )
            if sign_1 != sign_2
                first_idx_sign = ii
                break_flag = true
                break
            end
        end
        if break_flag
            break
        end

    end
    first_idx = min(first_idx_sign, first_idx_sigma)
    return first_idx
end

function __calculate_sigma(T_ij, gamma_k, n_eval)
    sigma_k = zeros(n_eval-1)
    for kk = 1:(n_eval-1)
        vi = view( T_ij, 1:(n_eval-kk), kk+1 )
        sum_T_ij = sum( vi.*vi )
        sigma_k[kk] = sqrt( gamma_k[kk] / (n_eval - kk) * sum_T_ij )
    end
    return sigma_k
end

# -------------------------- #
# dependent on the dimension #
# -------------------------- #

function __create_spacing(n_eval, x0::T, h::T) where {T <: Number}
    x_i = zeros(n_eval)

    x_i[ div(n_eval, 2)+1 ] = x0

    start_idx = div(n_eval, 2)
    for ii = start_idx:-1:1
        x_i[ii] = x0 - (start_idx + 1 - ii) * h
    end
    start_idx = div(n_eval, 2)+2
    for ii = start_idx:1:n_eval
        x_i[ii] = x0 + (ii - start_idx + 1) * h
    end

    return x_i
end

function __create_spacing(n_eval, x0::Vector{T}, h::Vector{T}) where {T <: Number}

    n_dim = length(x0)
    x_i = zeros(n_dim, n_eval)

    x_i[ :, div(n_eval, 2)+1 ] = x0

    start_idx = div(n_eval, 2)
    for ii = start_idx:-1:1
        x_i[:, ii] = x0 .- (start_idx + 1 - ii) .* h
    end
    start_idx = div(n_eval, 2)+2
    for ii = start_idx:1:n_eval
        x_i[:, ii] = x0 .+ (ii - start_idx + 1) .* h
    end

    return x_i
end

function __allocate_table(n_eval)
    T_ij = zeros(n_eval, n_eval)
    return T_ij
end

function __eval_function(fcn, x_i::Vector{T}, T_ij) where {T}
    for ii = 1:size(T_ij, 1)
        T_ij[ii, 1] = fcn( x_i[ii] )
    end
    return T_ij
end

function __eval_function(fcn, x_i::Matrix{T}, T_ij) where {T}
    n_eval = size(T_ij, 2)
    for ii = 1:n_eval
        T_ij[ii, 1] = fcn( x_i[:, ii] )
    end
    return T_ij
end

function __calculate_differences(T_ij, n_eval)
    for kk = 1:(n_eval-1)
        for ii = 1:(n_eval-kk)
            T_ij[ii, kk+1] = T_ij[ii+1, kk] - T_ij[ii, kk]
        end
    end

    return T_ij
end

function __difference_table(fcn, x_i::Vector{T}, n_eval::Int) where {T}
    # fcn: R -> R

    # allocate the table
    T_ij = __allocate_table(n_eval)

    # calculate the divided differences
    T_ij = __eval_function(fcn, x_i, T_ij)

    T_ij = __calculate_differences(T_ij, n_eval)

    return T_ij
end

function __difference_table(fcn, x_i::Matrix{T}, n_eval::Int) where {T}
    # fcn: R^n -> R

    # allocate the table
    T_ij = __allocate_table(n_eval)

    # calculate the divided differences
    T_ij = __eval_function(fcn, x_i, T_ij)

    T_ij = __calculate_differences(T_ij, n_eval)

    return T_ij
end

function estimate_noise(
    fcn,
    x0;
    n_eval = 7,
    h = 1e3 * sqrt.(eps.(x0))
    )
    # fcn : R -> R

    n_eval      = __oddify(n_eval)
    x_i         = __create_spacing(n_eval, x0, h)
    gamma_k     = __calculate_gamma(n_eval)
    T_ij        = __difference_table(fcn, x_i, n_eval)
    sigma_k     = __calculate_sigma(T_ij, gamma_k, n_eval)
    first_idx   = __selection_criterion(T_ij, sigma_k)
    noise_sigma = __select_noise_sigma(sigma_k, first_idx)

    return noise_sigma
end

function estimate_noise(
    fcn,
    x0::Vector{T};
    step_direction = normalize( 0.5.-rand(length(x0)) ),
    n_eval::Int = 7,
    h = 1e3 * sqrt.(eps.(x0))
    ) where {T}
    # fcn : R^n -> R

    n_eval      = __oddify(n_eval)                                      # indepedent of dim
    x_i         = __create_spacing(n_eval, x0, step_direction .* h)     # double implementation
    gamma_k     = __calculate_gamma(n_eval)                             # indepedent of dim
    T_ij        = __difference_table(fcn, x_i, n_eval)                  # double implementation
    sigma_k     = __calculate_sigma(T_ij, gamma_k, n_eval)              # indepedent of dim
    first_idx   = __selection_criterion(T_ij, sigma_k)                  # indepedent of dim
    noise_sigma = __select_noise_sigma(sigma_k, first_idx)              # indepedent of dim

    return noise_sigma
end

function noisy_numdiff_fw(fcn, x0)
    noise_est = estimate_noise(fcn, x0)
    h = max( sqrt(noise_est), sqrt.(eps.(x0)) )
    return numdiff_fw(fcn, x0, h)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# const LAYER = Layer
# const LAYER = NormalizedLayer

type NeuralNet{LAYER} <: NetStat
  layers::Vector{LAYER}  # note: this doesn't include input layer!!
  params::NetParams
  solverParams::SolverParams
  inputTransformer::Transformer
  transformedInput::VecF   # so we can avoid allocations
  # costmult::VecF

  # TODO: inner constructor which performs some sanity checking on activation/cost combinations:
  function NeuralNet{LAYER}(layers::Vector{LAYER},
                      params::NetParams,
                      solverParams::SolverParams,
                      inputTransformer::Transformer = IdentityTransformer())

    # do some sanity checking on activation/costmodel combos
    if isa(params.mloss, CrossentropyLoss)
      @assert isa(layers[end].activation, layers[end].nout > 1 ? SoftmaxMapping : SigmoidMapping)
    end

    if isa(layers[end].activation, SoftmaxMapping)
      @assert isa(params.mloss, CrossentropyLoss)
    end

    new(layers, params, solverParams, inputTransformer, zeros(first(layers).nin)) #, zeros(last(layers).nout))
  end
end

# # simple constructor which creates all layers the same for given list of node counts.
# # structure should include neuron counts for all layers, including input and output
# function NeuralNet(structure::AVec{Int};
#                     params = NetParams(),
#                     solverParams = SolverParams(),
#                     activation::Mapping = TanhMapping(),
#                     inputTransformer::Transformer = IdentityTransformer())
#   @assert length(structure) > 1

#   layers = LAYER[]
#   for i in 1:length(structure)-1
#     nin, nout = structure[i:i+1]
#     pDropout = getDropoutProb(params, i==1)
#     push!(layers, layerType(nin, nout, activation, params.updater, pDropout))
#   end

#   NeuralNet(layers, params, solverParams, inputTransformer)
# end

function Base.show(io::IO, net::NeuralNet)
  println(io, "NeuralNet{")
  println(io, "  params: $(net.params)")
  println(io, "  solverParams: $(net.solverParams)")
  println(io, "  layers:")
  for layer in net.layers
    println(io, "    ", layer)
  end
  println(io, "}")
end
Base.print(io::IO, net::NeuralNet) = show(io, net)

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

# produces a vector of output (estimated outputs) from the network
function forward!(net::NeuralNet, input::AVecF, istraining::Bool = false)

  # first transform the input
  transform!(net.inputTransformer, net.transformedInput, input)

  # now feed it forward
  output = net.transformedInput
  for layer in net.layers
    forward!(layer, output, istraining)
    output = layer.a
  end

  # update nextr
  for i in 1:length(net.layers)-1
    net.layers[i].nextr = net.layers[i+1].r
  end

  output
end

# given a vector of errors (target - output), update network weights
function backward!(net::NeuralNet, output, target)

  # update δᵢ starting from the output layer
  layer = net.layers[end]
  # δ = sensitivity!(layer.δ, layer.activation, net.params.mloss, layer.x, output, target)
  updateSensitivities!(layer, net.params.mloss, output, target)

  # now update the remaining sensitivities using bakckprop
  for i in length(net.layers)-1:-1:1
    updateSensitivities!(net.layers[i:i+1]...)
  end

  # now update the weights
  for layer in net.layers
    updateWeights!(layer, net.params.updater)

    # println()
    # @show layer
  end

  # # update our η, μ, etc
  # fit!(net.params)
end

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

# ϕₒ = sensitivity!(outputnode.state.ϕ,
#                   outputnode.mapping,
#                   net.mloss,
#                   input,
#                   output,
#                   target)

# online version... returns the feedforward estimate before updating
function OnlineStats.fit!(net::NeuralNet, input::AVecF, target::AVecF)
  output = forward!(net, input, true)
  backward!(net, output, target)
  output
end

# batch version
function OnlineStats.fit!(net::NeuralNet, input::MatF, target::MatF)
  @assert ncols(input) == net.nin
  @assert ncols(target) == net.nout
  @assert nrows(input) == nrows(target)

  Float64[fit!(net, row(input,i), row(target,i)) for i in 1:nrows(input)]
end

# note: when yEqualsX is true, we are updating an autoencoder (or similar) and so we can
# use net.transformedInput instead of target
function OnlineStats.fit!(net::NetStat, data::DataPoint, yEqualsX::Bool = false)
  fit!(net, data.x, yEqualsX ? net.transformedInput : data.y)
  # fit!(net, data.x, transformY ? transform(net.inputTransformer, data.y) : data.y)
end

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

function cost(net::NeuralNet, input::AVecF, target::AVecF)
  output = forward!(net, input)
  sum(value(net.params.mloss, target, output))
end
totalCost(net::NetStat, data::DataPoint) = cost(net, data.x, data.y)
totalCost(net::NetStat, dataset::DataPoints) = sum([totalCost(net, data) for data in dataset])
totalCost(net::NetStat, sampler::DataSampler) = totalCost(net, DataPoints(sampler))

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

function StatsBase.predict(net::NeuralNet, input::AVecF)
  forward!(net, input)
end

function StatsBase.predict(net::NeuralNet, input::AMatF)
  output = zeros(nrows(input), net.layers[end].nout)
  for i in 1:nrows(input)
    row!(output, i, predict(net, row(input, i)))
  end
  output
end

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

# note: we scale standard random normals by (1/sqrt(nin)) so that the distribution of initial (Σ = wx + b)
#       is also approximately standard normal
_initialWeights(nin::Int, nout::Int, activation::Mapping = IdentityMapping()) = 0.5randn(nout, nin) / sqrt(nin)

# initialWeights(nin::Int, nout::Int, activation::Mapping) = (rand(nout, nin) - 0.5) * 2.0 * sqrt(6.0 / (nin + nout))

# note: we scale standard random normals by (1/sqrt(nin)) so that the distribution of initial (Σ = wx + b)
#       is also approximately standard normal
# initialWeights(nin::Int, nout::Int, activation::Mapping) = randn(nout, nin) / sqrt(nin)

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


text:
using Random
using MLJModelInterface
using Test

mutable struct Foo <: MLJType
    rng::AbstractRNG
    x::Int
    y::Int
end

mutable struct Bar <: MLJType
    rng::AbstractRNG
    x::Int
    y::Int
end

mutable struct Super <: MLJType
    sub::Foo
    z::Int
end

mutable struct Partial <: MLJType
    x::Int
    y::Vector{Int}
    Partial(x) = new(x)
end

@testset "equality for MLJType" begin
    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(8), 1, 2)

    @test f1.rng != f2.rng
    @test f1 == f2
    f1.x = 10
    @test f1 != f2
    b = Bar(MersenneTwister(7), 1, 2)
    @test f2 != b

    @test is_same_except(f1, f2, :x)
    f1.y = 20
    @test f1 != f2
    @test is_same_except(f1, f2, :x, :y)

    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(8), 1, 2)
    s1 = Super(f1, 20)
    s2 = Super(f2, 20)
    @test s1 == s2
    s2.sub.x = 10
    @test f1 != f2

    @test !(f1 == Super(f1, 4))

    @test !(isequal(Foo(MersenneTwister(1), 1, 2),
                    Foo(MersenneTwister(1), 1, 2)))

    p1 = Partial(1)
    p2 = Partial(1)
    p2.y = [1,2]
    @test !(p1 == p2)

end

@testset "in(x, collection) for MLJType" begin
    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(7), 1, 2)
    f3 = Super(f1, 20)

    tv = (f1, f3)
    tk = (f2, f3)
    tw = (f3, f3)
    v = [tv...]
    k = [tk...]
    w = [tw...]
    @test f1 in tv
    @test !(f1 in tk)
    @test !(f1 in tw)
    @test f1 in v
    @test !(f1 in k)
    @test !(f1 in w)
    @test f1 in Set(v)
    @test !(f1 in Set(k))
    @test !(f1 in Set(w))

end

@testset "isrepresented" begin
    m = Foo(MersenneTwister(7), 1, 2)
    m2 = Foo(MersenneTwister(8), 1, 2)
    n = Foo(MersenneTwister(7), 3, 2)
    p = Foo(MersenneTwister(9), 4, 5)
    models = [m, n]

    @test isrepresented(m, nothing) == false
    @test isrepresented(m, models)
    @test isrepresented(m2, models)
    @test !isrepresented(p, models)
end

true


text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames
using RData
using LinearAlgebra

Random.seed!(1)

data_dir = joinpath(TuringModels.project_root, "data")
kline_path = joinpath(data_dir, "Kline2.csv")
df = CSV.read(kline_path, DataFrame; delim=';')

dmat_path = joinpath(data_dir, "islandsDistMatrix.rda")
dmat = load(dmat_path)["islandsDistMatrix"]

df.society = 1:10;
df

# ## Model

using Turing

@model function m13_7(Dmat, society, logpop, total_tools)
    rhosq ~ truncated(Cauchy(0, 1), 0, Inf)
    etasq ~ truncated(Cauchy(0, 1), 0, Inf)
    bp ~ Normal(0, 1)
    a ~ Normal(0, 10)

    ## GPL2
    SIGMA_Dmat = etasq * exp.(-rhosq * Dmat.^2)
    SIGMA_Dmat = SIGMA_Dmat + 0.01I
    SIGMA_Dmat = (SIGMA_Dmat' + SIGMA_Dmat) / 2
    g ~ MvNormal(zeros(10), SIGMA_Dmat)

    log_lambda = a .+ g[society] .+ bp * logpop

    total_tools .~ Poisson.(exp.(log_lambda))
end

chns = sample(
    m13_7(dmat, df.society, df.logpop, df.total_tools),
    NUTS(),
    5000
)

# ## Output

# \defaultoutput{}

# ## Original output

m_13_7_rethinking = """
Inference for Stan model: 6422d8042e9cdd08dae2420ad26842f1.
    4 chains, each with iter=10000; warmup=2000; thin=1;
    post-warmup draws per chain=8000, total post-warmup draws=32000.

            mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
    g[1]   -0.27    0.01  0.44  -1.26  -0.50  -0.24  -0.01   0.53  3278    1
    g[2]   -0.13    0.01  0.43  -1.07  -0.34  -0.10   0.12   0.67  3144    1
    g[3]   -0.17    0.01  0.42  -1.10  -0.37  -0.14   0.07   0.59  3096    1
    g[4]    0.30    0.01  0.37  -0.51   0.12   0.30   0.50   1.02  3145    1
    g[5]    0.03    0.01  0.37  -0.77  -0.15   0.04   0.22   0.72  3055    1
    g[6]   -0.46    0.01  0.38  -1.31  -0.64  -0.42  -0.23   0.19  3306    1
    g[7]    0.10    0.01  0.36  -0.70  -0.07   0.11   0.29   0.78  3175    1
    g[8]   -0.26    0.01  0.37  -1.07  -0.43  -0.24  -0.06   0.40  3246    1
    g[9]    0.23    0.01  0.35  -0.52   0.07   0.25   0.42   0.89  3302    1
    g[10]  -0.12    0.01  0.45  -1.04  -0.36  -0.11   0.13   0.77  5359    1
    a       1.31    0.02  1.15  -0.98   0.61   1.31   2.00   3.69  4389    1
    bp      0.25    0.00  0.11   0.02   0.18   0.24   0.31   0.47  5634    1
    etasq   0.34    0.01  0.53   0.03   0.10   0.20   0.39   1.52  5643    1
    rhosq   1.52    0.09 11.82   0.03   0.16   0.39   0.96   7.96 15955    1
    lp__  925.98    0.03  2.96 919.16 924.20 926.34 928.14 930.67  7296    1
""";


text:
using CmdStan

cmd1 = `echo`
cmd2 = `2`
sa = ["3", "4"]
s = "5"

#r1 = cmd1 * cmd2; println(r1)
#r2 = cmd1 * sa; println(r2)
#r3 = cmd1 * s; println(r3)

cs = Base.AbstractCmd[]

for i in 1:5
  push!(cs, `echo $i`)
end

println()
r = CmdStan.par(cs)
run(r)

println()
r = CmdStan.par(`echo hello`, 6)
run(r)


text:
#=
Simulate a ClaDS2 tree conditionned on tip number
=#

"""
    sim_ClaDS2_ntips(n::Int64,σ::Float64,α::Float64,ε::Float64,λ0::Float64 ;prune_extinct = true, sed = 0.001,max_time = 5, max_simulation_try = 100)

Simulate a tree from the ClaDS model conditionned on the number of tips at present

# Arguments
- `n::Int64`: the goal number of tips at present.
- `σ::Float64`: the stochasticity parameter.
- `α::Float64`: the trend parameter.
- `ε::Float64`: turnover rate (extinction / speciation).
- `λ0::Float64`: initial speciation rate.

# Keyword arguments
- `prune_extinct::Bool`: Should extinct lineages be removed from the output tree? Default to false.
"""
function sim_ClaDS2_ntips(n::Int64,σ::Float64,α::Float64,ε::Float64,λ0::Float64 ;
    prune_extinct = true, sed = 0.001,
    max_time = 5, max_simulation_try = 100)
    # accesory functions that will be called latter

    bernouilli=Bernoulli(ε/(1+ε))
    function dead()
        return (rand(bernouilli)[1] == 1)
    end

    lambda_law = LogNormal(log(α),σ)
    function new_rates(λ)
        λ * rand(lambda_law, 2)
    end

    N = n*max_time                  # maximal number of species before killing the simulation
    s = λ0*sed
    breaked = 0

    while breaked < max_simulation_try

        # initialise the tree to one alive lineage with rate λ0
        current_node = 3
        alive = [1,2]
        rates = new_rates(λ0)
        branches = [[[0;1]]; [[0;2]]]
        branch_lengths = [0., 0.]
        dead_branches = [false,false]
        node_time = 0
        times = [0.]
        time_diffs = [0.]
        n_lineages = [1]
        scale = 1.

        while 0 < length(alive) < N
            time_int = randexp()/(sum(rates[alive])*(1+ε))          # time of the new event

            if time_int==0                                          # to avoid numerical issues
                breaked += 1
                break
            end

            if length(alive) == n
                samp = s * time_int                     # probability to sample the tree in this interval with n species
                node_time += time_int
                u = rand()

                if u < samp*scale                       # meaning the tree is sampled
                    time_int *= rand()                  # the end time is uniformly sampled along the branch

                    push!(times, node_time)
                    push!(time_diffs, time_int)
                    push!(n_lineages,length(alive))

                    for i in alive                      # set all extant branch end to hte present time
                        branch_lengths[i] += time_int
                    end

                    t = build_tree(branches, branch_lengths, rates, extinct = dead_branches, prune_extinct = prune_extinct, root_attributes=[λ0])
                    return t
                else
                    scale /= (1. - samp)                # so that sampling is uniform on all time points wit n species
                end
            else
                node_time += time_int
            end

            individual= splice!(alive, sample(eachindex(alive),Weights(rates[alive])))

            is_dead=dead()                                          # check if the new event is a death or a speciation
            push!(times, node_time)
            push!(time_diffs, time_int)

            for node in alive
                branch_lengths[node] += time_int
            end

            # update internal variables
            if is_dead
                branch_lengths[individual] += time_int
                dead_branches[individual] = true
                push!(n_lineages,length(alive))
            else
                sample_rates = new_rates(rates[individual])
                push!(dead_branches,false,false)
                push!(rates, sample_rates[1], sample_rates[2])
                push!(branch_lengths, 0., 0.);
                push!(alive,current_node,current_node+1)
                push!(n_lineages,length(alive))
                push!(branches,[individual;current_node], [individual;current_node+1])

                current_node = current_node+2
                branch_lengths[individual] += time_int
            end
        end

        if length(alive) > 0
            s *=2
        end
    end

    # if the function finished before a good tree was found, returns an empty tree
    return Tree()
end

#=
Simulate a ClaDS2 tree conditionned on time
=#

function sim_ClaDS2_time(root_age,σ,α,ε,λ0 ; return_if_extinct = true, max_node_number = Inf, make_tree = true, prune_extinct = false, return_if_max = true)

    # accesory functions that will be called latter
    bernouilli=Bernoulli(ε/(1+ε))
    function dead()
        return (rand(bernouilli)[1] == 1)
    end

    lambda_law = LogNormal(log(α),σ)
    function new_rates(λ::Float64)
        λ * rand(lambda_law, 2)
    end

    function aux_ext(time, rate, n_max)
        #println(n_max)
        if n_max == 0
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
        end
        node_time = randexp()/(rate*(1+ε))
        if node_time > time
                #offsprings, branch_length, attributes, n_nodes
                return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max
        else
                is_dead = dead()
                if is_dead
                    print("dead")
                    return Tree(Array{Tree,1}(undef,0), node_time, [rate], 1, false), n_max + 1
                else
                    offspring_rates = new_rates(rate)
                    left_time = time - node_time
                    #print("$depth ;")
                    left_tree = aux_ext(left_time, offspring_rates[1], n_max - 1)
                    #print(left_tree[2])
                    if left_tree[1].n_nodes == -1
                        return left_tree
                    end
                    right_tree = aux_ext(left_time, offspring_rates[2], left_tree[2])
                    if right_tree[1].n_nodes == -1
                        return right_tree
                    end
                    return Tree([left_tree[1], right_tree[1]], node_time, [rate]), right_tree[2]
                end
        end
    end
    tree = aux_ext(root_age, λ0, max_node_number)[1]

    if prune_extinct
        if (tree.n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false)
        else
            return prune_extinct_lineages(tree)
            println("pruned")
        end
    else
        if (tree.n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false)
        else
            return tree
        end
    end
end

#=
Simulate a ClaDS2 tree
=#

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
export hom, hom_check

"""
`hom_check(G,H,d)` checks if the dictionary `d` is a graph homomorphism
from `G` to `H`.
"""
function hom_check(G::SimpleGraph, H::SimpleGraph, d::Dict)::Bool
    for e in G.E
        v, w = e
        x = d[v]
        y = d[w]
        if !H[x, y]
            return false
        end
    end
    return true
end

"""
`hom(G,H)` finds a graph homomorphism from `G` to `H`.
"""
function hom(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}

    m = Model(get_solver())

    @variable(m, P[G.V, H.V], Bin)

    for v in G.V
        @constraint(m, sum(P[v, x] for x in H.V) == 1)
    end

    for e in G.E
        u, v = e
        for x in H.V
            for y in H.V
                if !H[x, y]
                    @constraint(m, P[u, x] + P[v, y] <= 1)
                end
            end
        end
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error("No homomorphism")
    end

    PP = Int.(value.(P))

    result = Dict{S,T}()

    for v in G.V
        for x in H.V
            if PP[v, x] > 0
                result[v] = x
            end
        end
    end

    return result

end

max_stars_repo_name: scheinerman/SimpleGraphAlgorithms

max_issues_repo_name: scheinerman/SimpleGraphAlgorithms

max_forks_repo_name: scheinerman/SimpleGraphAlgorithms


text:
using Weave, Compat
using Base.Test

cleanup = true

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="tex")
result = @compat readstring(open("documents/pyplot_formats.tex"))
ref = @compat readstring(open("documents/pyplot_formats_ref.tex"))
result = replace(result, r"\s*PyObject.*\n", "\n") #Remove PyObjects, because they change
ref = replace(ref, r"\s*PyObject.*\n", "\n")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="github")
result = @compat readstring(open("documents/pyplot_formats.md"))
ref = @compat readstring(open("documents/pyplot_formats_ref.md"))
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="rst", fig_ext=".svg")
result = @compat readstring(open("documents/pyplot_formats.rst"))
ref = @compat readstring(open("documents/pyplot_formats_ref.rst"))
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

if cleanup
    rm("documents/pyplot_formats.tex")
    rm("documents/pyplot_formats.rst")
    rm("documents/pyplot_formats.md")
    rm("documents/figures", recursive = true)
end


text:
const Solution = Dict{String, Union{Int, Bool, Set{Int}}}

#TODO add documentation for BeforeRestart
mutable struct Statistics
    infeasibleStatusPerVariable             ::Dict{String, Int}
    numberOfNodes                           ::Int
    numberOfSolutions                       ::Int
    numberOfInfeasibleSolutions             ::Int
    numberOfSolutionsBeforeRestart          ::Int
    numberOfInfeasibleSolutionsBeforeRestart::Int
    numberOfNodesBeforeRestart              ::Int
    AccumulatedRewardBeforeReset            ::Float32
    solutions                               ::Vector{Union{Nothing,Solution}}
    nodevisitedpersolution                  ::Vector{Int}
    objectives                              ::Union{Nothing, Vector{Union{Nothing,Int}}}
    lastPruning                             ::Union{Nothing, Int}
    lastVar                                 ::Union{Nothing, AbstractIntVar} #last var on which we branched
end

mutable struct Limit
    numberOfNodes       ::Union{Int, Nothing}
    numberOfSolutions   ::Union{Int, Nothing}  #the limit can be triggered by set of non-unique solutions
    searchingTime       ::Union{Int, Nothing}
end

"""
    CPModel(trailer::Trailer)
    CPModel()

The structure storing all the informations needed to solve a specific problem and it also stores the solutions.
The CPModel is the center of the solver and evolve during the solving.
The `AbstractStateRepresentation` used by the RL Agent is created from the CPModel.

The CPModel is always created empty and is filled eather by hand by the user (or automatically thanks to written files)
or filled by an `AbstractModelGenerator`.

"""
mutable struct CPModel
    variables               ::Dict{String, AbstractVar}
    branchable              ::Dict{String, Bool}
    constraints             ::Array{Constraint}
    trailer                 ::Trailer
    objective               ::Union{Nothing, AbstractIntVar}
    objectiveBound          ::Union{Nothing, Int}
    statistics              ::Statistics
    limit                   ::Limit
    knownObjective          ::Union{Nothing,Int64}
    adhocInfo               ::Any
    CPModel(trailer) = new(Dict{String, AbstractVar}(), Dict{String, Bool}(), Constraint[], trailer, nothing, nothing, Statistics(Dict{String, Int}(), 0,0, 0, 0, 0, 0, 0, Solution[],Int[], nothing, nothing, nothing), Limit(nothing, nothing, nothing), nothing)
end

CPModel() = CPModel(Trailer())

"""
    addVariable!(model::CPModel, x::AbstractVar; branchable=true)

Add a variable to the model, throwing an error if `x`'s id is already in the model.
The `branchable` argument allows you to tell if we will be able to branch on that variable.
"""
function addVariable!(model::CPModel, x::AbstractVar; branchable=true)
    # Ensure the id is unique
    @assert !haskey(model.variables, x.id) "The id of the variable must be unique"

    @assert !branchable || typeof(x) <: Union{AbstractIntVar, AbstractBoolVar} "You can only branch on Boolean and Integer variables"

    model.statistics.infeasibleStatusPerVariable[id(x)]=0
    model.branchable[x.id] = branchable
    model.variables[x.id] = x
end

"""
    addObjective!(model::CPModel, objective::AbstractVar)

Add an Objective variable to the model. This variable is the variable that needs to be minimized suring the solving.
"""
function addObjective!(model::CPModel, objective::AbstractVar)
    model.objective = objective
    model.statistics.objectives = Int[]  #initialisation of the Array that will contain the score of every solution
end

function addKnownObjective!(model::CPModel, knownObective::Int64)
    model.knownObjective = knownObective
end

function addConstraint!(model::CPModel,constraint::Constraint)
    push!(model.constraints,constraint)
    for var in variablesArray(constraint)
        if haskey(model.branchable, id(var))
            @assert haskey(model.statistics.infeasibleStatusPerVariable, id(var)) "You forget to add the variable $(id(var)) to the model"
            model.statistics.infeasibleStatusPerVariable[id(var)]+=1
        end
    end
end

"""
    function is_branchable(model::CPModel, x::AbstractVar)

Tell if the variable was set as branchable or not.
"""
is_branchable(model::CPModel, x::AbstractVar) = haskey(model.branchable, x.id) && model.branchable[x.id]

"""
    function branchable_variables(model::CPModel)

Return a dict of all branchable variables.
"""
function branchable_variables(model::CPModel)
    to_return = Dict{String, AbstractVar}()
    for x in values(model.variables)
        if is_branchable(model, x)
            to_return[x.id] = x
        end
    end
    to_return
end

"""
    solutionFound(model::CPModel)

Return a boolean, checking whether a solution was found, i.e. every variable is bound.
"""
function solutionFound(model::CPModel)
    for (k, x) in model.variables
        if !isbound(x)
            return false
        end
    end
    return true
end

"""
    triggerFoundSolution!(model::CPModel)

Add the current solution to `model`, and set new constraints for the objective if needed.
"""
function triggerFoundSolution!(model::CPModel)
    @assert solutionFound(model)

    model.statistics.numberOfSolutions += 1
    model.statistics.numberOfSolutionsBeforeRestart += 1

    # Adding the solution
    solution = Solution()
    for (k, x) in model.variables
        solution[k] = assignedValue(x)
    end
    if !(solution in model.statistics.solutions)   #probably not efficient but necessary
        push!(model.statistics.solutions, solution)
        push!(model.statistics.nodevisitedpersolution,model.statistics.numberOfNodes)
        if !isnothing(model.objective)
            @assert !isnothing(model.statistics.objectives)   "did you used SeaPearl.addObjective! to declare your objective function ? "
            push!(model.statistics.objectives, assignedValue(model.objective))
            tightenObjective!(model)
        end
    end
end
"""
    triggerInfeasible!(constraint::Constraint, model::CPModel)

this function increments by one the statistic infeasibleStatusPerVariable for each variable involved in the constraint. infeasibleStatusPerVariable
keeps in track for each variable the number of times the variable was involved in a constraint that led to an infeasible state during a fixpoint. This statistic
is used by the failure-based variable selection heuristic.
"""
function triggerInfeasible!(constraint::Constraint, model::CPModel)
    for var in variablesArray(constraint)
        if haskey(model.branchable, id(var))
            model.statistics.infeasibleStatusPerVariable[id(var)]+=1
        end
    end
    push!(model.statistics.solutions, nothing)
    push!(model.statistics.nodevisitedpersolution,model.statistics.numberOfNodes)

    if !isnothing(model.objective)
        @assert !isnothing(model.statistics.objectives)   "did you used SeaPearl.addObjective! to declare your objective function ? "
        push!(model.statistics.objectives, nothing)
    end
end

"""
    tightenObjective!(model::CPModel)

Set a new constraint to minimize the objective variable.
"""
function tightenObjective!(model::CPModel)
    model.objectiveBound = assignedValue(model.objective)-1
end

"""
    belowLimits(model::CPModel)

Check if `model`' statistics are still under the limits.
"""
belowLimits(model::CPModel) = belowNodeLimit(model) && belowSolutionLimit(model) && belowTimeLimit(model)
belowNodeLimit(model::CPModel) = isnothing(model.limit.numberOfNodes) || model.statistics.numberOfNodes < model.limit.numberOfNodes
belowSolutionLimit(model::CPModel) = isnothing(model.limit.numberOfSolutions) || model.statistics.numberOfSolutions < model.limit.numberOfSolutions
belowTimeLimit(model::CPModel) = isnothing(model.limit.searchingTime) || peektimer() < model.limit.searchingTime
"""
    Base.isempty(model::CPModel)::Bool

Return a boolean describing if the model is empty or not.
"""
function Base.isempty(model::CPModel)::Bool
    (
        isempty(model.variables)
        && isempty(model.constraints)
        && isempty(model.trailer.prior)
        && isempty(model.trailer.current)
        && isnothing(model.objective)
        && isnothing(model.objectiveBound)
        && isempty(model.statistics.solutions)
        && isempty(model.statistics.nodevisitedpersolution)
        && isempty(model.statistics.infeasibleStatusPerVariable)
        && isnothing(model.statistics.objectives)
        && isnothing(model.statistics.lastPruning)
        && isnothing(model.statistics.lastVar)
        && model.statistics.numberOfNodes == 0
        && model.statistics.numberOfSolutions == 0
        && model.statistics.numberOfInfeasibleSolutions == 0
        && model.statistics.numberOfInfeasibleSolutionsBeforeRestart == 0
        && model.statistics.numberOfSolutionsBeforeRestart == 0
        && model.statistics.numberOfNodesBeforeRestart == 0
        && model.statistics.AccumulatedRewardBeforeReset == 0
        && isnothing(model.limit.numberOfNodes)
        && isnothing(model.limit.numberOfSolutions)
        && isnothing(model.limit.searchingTime)
        && isnothing(model.knownObjective)
    )
end

"""
    Base.empty!(model::CPModel)

Empty the CPModel.
"""
function Base.empty!(model::CPModel)
    empty!(model.variables)
    empty!(model.constraints)
    empty!(model.trailer.prior)
    empty!(model.trailer.current)
    model.objective = nothing
    model.objectiveBound = nothing
    empty!(model.statistics.solutions)
    empty!(model.statistics.nodevisitedpersolution)
    empty!(model.statistics.infeasibleStatusPerVariable)
    model.statistics.objectives = nothing
    model.statistics.lastPruning = nothing
    model.statistics.lastVar = nothing
    model.statistics.numberOfNodes = 0
    model.statistics.numberOfSolutions = 0
    model.statistics.numberOfInfeasibleSolutions = 0
    model.statistics.numberOfInfeasibleSolutionsBeforeRestart = 0
    model.statistics.numberOfSolutionsBeforeRestart = 0
    model.statistics.numberOfNodesBeforeRestart = 0
    model.statistics.AccumulatedRewardBeforeReset = 0
    model.limit.numberOfNodes = nothing
    model.limit.numberOfSolutions = nothing
    model.limit.searchingTime = nothing
    model.knownObjective = nothing

    model
end

"""
    reset_model!(model::CPModel)

Reset a given CPModel instance. Make it possible to reuse the same instance instead of having to
delete the old one and create another one. This is used in `launch_experiment!` in order to be able
to use the same CPModel instance to compare different given heuristics.
"""
function reset_model!(model::CPModel)
    restoreInitialState!(model.trailer)
    model.objectiveBound = nothing
    empty!(model.statistics.solutions)
    empty!(model.statistics.nodevisitedpersolution)
    for (key, value) in model.statistics.infeasibleStatusPerVariable
        model.statistics.infeasibleStatusPerVariable[key]=length(getOnDomainChange(model.variables[key]))  #the degree is reset to the initial value : the number of constraints the variable is involved in.
    end
    if !isnothing(model.objective)
        @assert !isnothing(model.statistics.objectives)   "did you used SeaPearl.addObjective! to declare your objective function ?"
        empty!(model.statistics.objectives)
    end
    model.statistics.lastPruning = nothing
    model.statistics.lastVar = nothing
    model.statistics.numberOfNodes = 0
    model.statistics.numberOfSolutions = 0
    model.statistics.numberOfInfeasibleSolutions = 0
    model.statistics.numberOfInfeasibleSolutionsBeforeRestart = 0
    model.statistics.numberOfSolutionsBeforeRestart = 0
    model.statistics.numberOfNodesBeforeRestart = 0
    model.statistics.AccumulatedRewardBeforeReset = 0

end

"""
restart_search!(model::CPModel)

Usefull when dealing with restart based search : ILDS or RBS. Reset to zero usefull statistics on the search that can be used to define
the restart criteria.
"""
function restart_search!(model::CPModel)
    restoreInitialState!(model.trailer)
    model.statistics.lastPruning = nothing
    model.statistics.lastVar = nothing
    model.statistics.numberOfInfeasibleSolutionsBeforeRestart = 0
    model.statistics.numberOfSolutionsBeforeRestart = 0
    model.statistics.numberOfNodesBeforeRestart = 0
end

"""
    domains_cartesian_product(model::CPModel)

Return the cartesian product of the model variables: |D1|x|D2|x ... x|Dn|
Helps providing insights about what is happening during a search.
"""
function domains_cartesian_product(model::CPModel)
    cart_pdt = 1
    for (id, x) in model.variables
        if typeof(x) <: Union{AbstractIntVar, AbstractBoolVar}
            cart_pdt *= length(x.domain)
        end
    end
    return cart_pdt
end

"""
    nb_boundvariables(model::CPModel)

Return the number of variables that have already been assigned to a value.
Helps providing insights about what is happening during a search.
"""
function nb_boundvariables(model::CPModel)
    nb = 0
    for (id, x) in model.variables
        nb += isbound(x) * 1
    end
    return nb
end


text:
ENV["GKSwstype"] = "100" #src

# # Data Vizualization
#
# > The greatest value of a picture is when it forces us to notice what we never expected to see.
# >
# > \- John Tukey
#
# Data visualization is used throughout the data science workflow.  It's used for everything from identifying simple relationships between variables to understanding how hyper-parameters affect model performance.  This course module will focus on exploratory data analysis using visualizations.
#
# Note that Julia does not have a built-in plotting library, but there are many to choose from:
#
# - [**`Gadfly`**](https://github.com/GiovineItalia/Gadfly.jl)
# - [**`GR`**](https://github.com/jheinen/GR.jl)
# - [**`PGFPlotsX`**](https://github.com/KristofferC/PGFPlotsX.jl)
# - [**`PlotlyJS`**](https://github.com/sglyon/PlotlyJS.jl)
# - [**`Plots`**](https://github.com/JuliaPlots/Plots.jl)
# - [**`PyPlot`**](https://github.com/JuliaPy/PyPlot.jl)
# - [**`UnicodePlots`**](https://github.com/Evizero/UnicodePlots.jl)
# - [**`Winston`**](https://github.com/JuliaGraphics/Winston.jl)
#
# If you are new to Julia, this may be an overwhelming number of options.  We highly recommend the **`Plots`** package, as it is a great option for both newcomers and power users.  This course module uses **`Plots`** and its related packages.

# # What is [Plots.jl](https://github.com/JuliaPlots/Plots.jl)?
#
# ![](https://raw.githubusercontent.com/JuliaPlots/PlotReferenceImages.jl/master/PlotDocs/index/lorenz_attractor.gif)
#
# **`Plots`** defines an interface for plotting that works with multiple plotting libraries "backends", including **`PyPlot`**, **`GR`**, and **`UnicodePlots`**.  It allows you to switch between backends, making it painless to go from an interactive javascript plot to something print-ready.  Let's get started with some of the core **`Plots`** concepts.

using Plots
gr()  # set the backend as GR

# The `plot` function has many different methods.

p1 = plot(rand(10))             # y variable only
p2 = plot(-5:4, rand(10))       # x and y variables
p3 = plot(rand(10, 3))          # multiple y variables (each column)
p4 = plot(-5:4, rand(10, 3))    # x and multiple y variables

plot(p1, p2, p3, p4)            # join plots into subplots

# There are many attributes available to change via keyword arguments.

plot(rand(10), seriestype = :scatter, markersize = 5, markercolor = :black)

# ## Series Types
#
# One of the most important keyword arguments to `plot` is `seriestype`, as it determines how the input data is turned into a plot.  For convenience, a `seriestype` also has a function of the same name:

y = rand(10)

scatter(y)  # same as plot(y, seriestype = :scatter)

# ## Argument Aliases
#
# **Plots** supports "aliases" for keyword arguments for fast interactive creation of plots.  You typically don't need to keep checking the documentation for the exact keyword.  Aliases are intuitive and you'll often be able to guess the correct syntax.

plot(y, markershape = :circle)

#-

plot(y, m = :circle)

# ## Magic Arguments
#
# *Magic arguments* is one of **`Plots`**' unique features.  Certain keyword arguments act like "magic" where by providing a tuple, Plots intelligently maps the items to the corresponding plot element.  For example, by providing `marker = (5, .3, :auto)`, the marker size is set as 5, the marker alpha as .3, and the marker shape is automatically determined.  Note that we can also use the alias `m` rather than `marker`:

scatter(randn(10, 3), m = (5, .3, :auto))

# ## Plotting from datasets
#
# Let's now use what we've learned about **`Plots`** to explore the `mtcars` dataset, extracted from a 1974 issue of of *Motor Trend* Magazine.  The variables included are:
#
# - `Model`: Make and model of the car
# - `MPG`: Miles per gallon
# - `Cyl`: Number of Cylinders
# - `Disp`: Displacement (cubic inches)
# - `HP`: Gross horsepower
# - `DRat`: Rear axle ratio
# - `WT`: Weight (in 1000s lbs.)
# - `QSec`: Quarter mile time
# - `VS`: Engine type, V-shaped (0) or straight (1)
# - `AM`: Transmission, Automatic (0) or manual (1)
# - `Gear`: Number of forward gears
# - `Carb`: Number of carburetors

using RDatasets

mtcars = dataset("datasets", "mtcars")

# First, let's `describe` the data.

describe(mtcars)

# What's the relationship between miles per gallon and horsepower?

scatter(mtcars.MPG, mtcars.HP, xlab="MPG", ylab = "HP")

# The `dataframe.column` syntax isn't always convenient for making complex visualizations.  The **`StatsPlots`** package adds a number of features to **`Plots`** for creating statistical graphics.  One of its main features is the `@df` macro that allows you to reference columns of a data table by name inside of a call to `plot`

using StatsPlots

@df mtcars corrplot([:MPG :Cyl :HP :AM], bins=25, grid=false)

# Now let's look at the distribution of `MPG` across `Cyl` using a [violin plot](https://en.wikipedia.org/wiki/Violin_plot).

@df mtcars violin(:Cyl, :HP, legend = :topleft, xlab="Cyl", ylab = "HP")

# **`Plots`** objects can also be mutated with additional data series.  Let's now overlap [box plots](https://en.wikipedia.org/wiki/Box_plot) of the same data, using the `alpha` keyword argument to make the box plots transparent.

@df mtcars boxplot!(:Cyl, :HP, alpha = .5)

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
macro randtest(expr)
    :(@test begin
        srand(100)
        $expr
    end)
end
macro randshow(expr)
    :(begin
        srand(100)
        @show $expr
    end)
end

const just0   = pc.just(0)
const just5   = pc.just(5)
const just10  = pc.just(10)
const just127 = pc.just(127)

const justf0  = pc.just(0.0)
const justf10 = pc.just(10.0)

const signeds_0to10 = pc.signeds(0, 10)
const signeds_5to10 = pc.signeds(5, 10)
const floats_0to10 = pc.floats(0.0, 10.0)
const bools = pc.bools()

const utf8strings_e0to127_l5to10 = pc.utf8strings(0, 127, signeds_5to10)

const vec_e0to10_l0     = pc.vectors(floats_0to10, just0)
const vec_e0to10_l5     = pc.vectors(signeds_0to10, just5)
const vec_e0to10_l5to10 = pc.vectors(signeds_0to10, signeds_5to10)

const dict_k0to10_v0to10_l0 = pc.dicts(signeds_0to10, signeds_0to10, just0)
const dict_k0to127_v0to10_l5to10 = pc.dicts(utf8strings_e0to127_l5to10, signeds_0to10, signeds_0to10)

const set_e0to10_l0 = pc.sets(signeds_0to10, just0)
const set_e0to127_l5 = pc.sets(utf8strings_e0to127_l5to10, just5)

const tuple_00 = pc.tuples(just0, just0)

between{T}(x::T, lb::T, ub::T) = (lb <= x) && (x <= ub)
length_between(x, lb, ub) = between(length(x), lb, ub)

# ----- Just -----

@test pc.arbitrary(just0) == 0
@test collect(pc.shrink(0, just0)) == []

@test pc.arbitrary(just5) == 5
@test collect(pc.shrink(0, just5)) == []

# ----- Numbers -----

@test pc.forAll(between, (signeds_0to10, just0, just10))
@test collect(pc.shrink(1, signeds_0to10)) == Any[]

@randtest pc.forAll(between, (floats_0to10, justf0, justf10))
@test collect(pc.shrink(5.0, floats_0to10)) == Any[]

@test collect(pc.shrink(true, bools)) == Bool[]

# ----- UTF8Strings -----

@test pc.forAll(length_between, (utf8strings_e0to127_l5to10, just5, just10))
@test collect(pc.shrink("123456789", utf8strings_e0to127_l5to10)) == ["12345", "56789"]
@test collect(pc.shrink("123456", utf8strings_e0to127_l5to10)) == ["12345", "23456"]
@test collect(pc.shrink("12345", utf8strings_e0to127_l5to10)) == []

# ----- Vectors -----

@test pc.arbitrary(vec_e0to10_l0) == Float64[]
@test collect(pc.shrink(Float64[], vec_e0to10_l0)) == []

@randtest pc.arbitrary(vec_e0to10_l5) == [1,1,9,10,2]
@test collect(pc.shrink([0,6,9,10,2], vec_e0to10_l5)) == []

@test collect(pc.shrink([0,3,6,9,10,2], vec_e0to10_l5to10)) == Vector[[0,3,6,9,10], [3,6,9,10,2]]

# ----- Dicts -----

@test pc.arbitrary(dict_k0to10_v0to10_l0) == Dict{Int64, Int64}()

@randtest pc.arbitrary(dict_k0to127_v0to10_l5to10) == Dict("#A6Ix"=>0)

# ----- Sets -----

@test pc.arbitrary(set_e0to10_l0) == Set{Int64}()

@randtest pc.arbitrary(set_e0to127_l5) == Set(UTF8String["7\x120XfuxC&","#\0%\t_","\0#A6Ix","!}\x06*j\x17\x13&","V(j]\x7f"])

# ----- Tuples -----

@test pc.arbitrary(tuple_00) == (0,0)

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
module XKCD

using HTTP, JSON3

export Comic

#-----------------------------------------------------------------------------# Comic
"""
    Comic()         # Most recent comic
    Comic(i::Int)   # i-th issue

Get the metadata for an XKCD comic via the XKCD API: `https://xkcd.com/json.html`.

# Example
    Comic(552)  # "Correlation"
"""
struct Comic
    json::JSON3.Object
end

Comic() = Comic(JSON3.read(HTTP.get("https://xkcd.com/info.0.json").body))
Comic(i::Int) = Comic(JSON3.read(HTTP.get("https://xkcd.com/$i/info.0.json").body))

#-----------------------------------------------------------------------------# show
function Base.show(io::IO, o::Comic)
    println(io, "Comic")
    show(io, o.json)
end

function Base.show(io::IO, ::MIME"text/html", c::Comic)
    show(io, MIME"text/html"(), HTML("""

XKCD.Comic $(c.json.num): $(c.json.title)

            Link to Original

    """))
end

#-----------------------------------------------------------------------------# rand
Base.rand(::Type{Comic}) = Comic(rand(1:Comic().json.num))

end # module


text:
using Test
using GeoParams
@testset "Density.jl" begin

#Set alias for density function
@use GeoParamsAliases density=ρ

#Make sure that structs are isbits
x = ConstantDensity()
@test isbits(x)

x = PT_Density()
@test isbits(x)

x = Compressible_Density()
@test isbits(x)

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);

# Define a linear viscous creep law
x1      =   ConstantDensity(ρ=2900kg/m^3)
@test x1.ρ.val == 2900

x1 = nondimensionalize(x1,CharUnits_GEO)
@test x1.ρ.val ≈ 2.9e-16

x2      =   PT_Density()
@test x2.α.val==3e-5
@test x2.ρ0.val==2900.0

x2 = nondimensionalize(x2,CharUnits_GEO)
@test x2.T0.val≈0.21454659702313156

# Compute with density while specifying P & T (not used in case of )
@test compute_density(x2,1.0,1.0) ≈ 2.8419999999999996e-16
@test compute_density(x1,1.0,1.0) ≈ 2.9e-16

# test to allocations
rho = [0.0]
P   = 1.0;
T   = 1.0

x   = ConstantDensity()
num_alloc = @allocated compute_density!(rho, x, P, T)
num_alloc = @allocated compute_density!(rho, x, P, T)
@show num_alloc
@test num_alloc == 0

#Test allocations using ρ alias
num_alloc = @allocated ρ!(rho, x, P, T)
num_alloc = @allocated ρ!(rho, x, P, T)
@test num_alloc == 0

# This does NOT allocate if I test this with @btime;
#   yet it does while running the test here
x   = PT_Density()
num_alloc = @allocated compute_density!(rho, x, P, T)
num_alloc = @allocated compute_density!(rho, x, P, T)
@show num_alloc
@test num_alloc <= 32

# This does NOT allocate if I test this with @btime;
#   yet it does while running the test here
x   = Compressible_Density()
num_alloc = @allocated compute_density!(rho, x, P, T)
num_alloc = @allocated compute_density!(rho, x, P, T)
@show num_alloc
@test num_alloc <=32

# Read Phase diagram interpolation object
fname   =   "./test_data/Peridotite.in"
PD_data =   PerpleX_LaMEM_Diagram(fname);
@test PD_data.meltFrac(1500,1e7) ≈ 0.24368492372485706
@test PD_data.Rho(1500,1e7) ≈ 3042.836820256982
@test PD_data.meltRho(1500,1e7) ≈ 2662.227167592414
@test PD_data.rockRho(1500,1e7) ≈ 3165.467673917775

@test compute_density(PD_data, 1e7, 1500.0) ≈ 3042.836820256982

# Do the same but non-dimensionalize the result
CharDim  =  GEO_units();
PD_data1 =  PerpleX_LaMEM_Diagram(fname, CharDim=CharDim);

rho_ND   =  PD_data1.Rho(nondimensionalize(1500.0K,CharDim), nondimensionalize(1e8*Pa,CharDim))
Vp_ND    =  PD_data1.Vp( nondimensionalize(1500.0K,CharDim), nondimensionalize(1e8*Pa,CharDim))
Vs_ND    =  PD_data1.Vs( nondimensionalize(1500.0K,CharDim), nondimensionalize(1e8*Pa,CharDim))

# redimensionalize and check with value from original structure that did not use non-dimensionalization
@test   ustrip(dimensionalize(rho_ND,kg/m^3,CharDim)) ≈ PD_data.Rho(1500.0,1e8)
@test   ustrip(dimensionalize(Vp_ND, km/s,  CharDim)) ≈ PD_data.Vp( 1500.0,1e8)
@test   ustrip(dimensionalize(Vs_ND, km/s,  CharDim)) ≈ PD_data.Vs( 1500.0,1e8)

# Test computation of density for the whole computational domain, using arrays
MatParam    =   Vector{AbstractMaterialParamsStruct}(undef, 4)

MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=0,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"));

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=1,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = ConstantDensity(ρ=2900kg/m^3));

MatParam[3] =   SetMaterialParams(Name="UpperCrust", Phase=2,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = PT_Density());
MatParam[4] =   SetMaterialParams(Name="UpperCrust", Phase=3,
                        CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                        Density   = Compressible_Density());

Mat_tup = Tuple(MatParam);  # create a tuple to avoid allocations

MatParam1 = Vector{AbstractMaterialParamsStruct}(undef, 4)
MatParam1[1] = SetMaterialParams(Name="Crust", Phase=0,
                            CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pas)),
                            Density   = ConstantDensity(ρ=2900kg/m^3))
MatParam1[2] = SetMaterialParams(Name="Lower Crust", Phase=1,
                            CreepLaws= (PowerlawViscous(n=5.), LinearViscous(η=1e21Pas)),
                            Density  = Compressible_Density(ρ0=3000kg/m^3))
MatParam1[3] = SetMaterialParams(Name="Lower Crust", Phase=2,
                            CreepLaws= LinearViscous(η=1e21Pas),
                            Density  = ConstantDensity())
MatParam1[4] = SetMaterialParams(Name="Lower Crust", Phase=3,
                            CreepLaws= LinearViscous(η=1e21Pas),
                            Density  = ConstantDensity())
Mat_tup1 = Tuple(MatParam1)

# test computing material properties
Phases              = ones(Int64,400,400)*0;
Phases[:,20:end]   .= 1
Phases[:,200:end]  .= 2
Phases[:,300:end]  .= 3

#Phases .= 2;
rho     = zeros(size(Phases))
T       =  ones(size(Phases))
P       =  ones(size(Phases))*10

# Test computing density when Mat_tup1 is provided as a tuple
compute_density!(rho, Mat_tup1, Phases, P,T)
num_alloc = @allocated compute_density!(rho, Mat_tup1, Phases, P,T)   #      287.416 μs (0 allocations: 0 bytes)
@test sum(rho)/400^2 ≈ 2945.000013499999
@test num_alloc == 0

#Same test using function alias
rho     = zeros(size(Phases))
ρ!(rho, Mat_tup1, Phases, P,T)
num_alloc = @allocated compute_density!(rho, Mat_tup1, Phases, P,T)
@test sum(rho)/400^2 ≈ 2945.000013499999
@test num_alloc == 0

# If we employ a phase diagram many allocations occur:
compute_density!(rho, Mat_tup, Phases, P,T)   #        37.189 ms (1439489 allocations: 26.85 MiB)     - the allocations are from the phase diagram
@test sum(rho)/400^2 ≈ 2920.6151145725003

# test computing material properties when we have PhaseRatios, instead of Phase numbers
PhaseRatio  = zeros(size(Phases)...,length(Mat_tup1));
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz+1)
    PhaseRatio[I] = 1.0
end

compute_density!(rho, Mat_tup1, PhaseRatio, P, T)
num_alloc = @allocated compute_density!(rho, Mat_tup1, PhaseRatio, P,T) #   136.776 μs (0 allocations: 0 bytes)
# using BenchmarkTools
# @btime compute_density!($rho, $Mat_tup1, $PhaseRatio, P=$P,T=$T)
@test sum(rho)/400^2 ≈ 2945.000013499999
@test num_alloc == 0           # for some reason this does indicate allocations but @btime does not

# Test calling the routine with only pressure as input.
# This is ok for Mat_tup1, as it only has constant & P-dependent densities.
# Note, however, that if you have P & T dependent densities and do this it will use 0 as defualt value for T
compute_density!(rho, Mat_tup1, PhaseRatio, P)
@test sum(rho)/400^2 ≈ 2945.000013499999

# In case we only want to compute with T, do this:
#  NOTE that in this example the results are actually wrong (as some functions require P as well)
compute_density!(rho, Mat_tup, PhaseRatio, nothing, T)
@test sum(rho)/400^2 ≈ 2920.6151

#Test computation of density given a single phase and P,T as scalars
Phase, P, T = 0, 1., 1.
@test compute_density(Mat_tup1, Phase, P, T) == 2900.0

end

#= Testing more than one density per field (change ln:205 in MaterialParameters.jl)
MatParam = Vector{AbstractMaterialParamsStruct}(undef, 4)
MatParam[1] = SetMaterialParams(Name="Crust", Phase=0,
                            CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pas)),
                            Density   = (PT_Density(), ConstantDensity()))
MatParam[2] = SetMaterialParams(Name="Lower Crust", Phase=1,
                            CreepLaws= (PowerlawViscous(n=5.), LinearViscous(η=1e21Pas)),
                            Density  = (Compressible_Density(ρ0=3000kg/m^3), ConstantDensity()))
MatParam[3] = SetMaterialParams(Name="Lower Crust", Phase=2,
                            CreepLaws= LinearViscous(η=1e21Pas),
                            Density  = (ConstantDensity(), PT_Density(), Compressible_Density(ρ0=3000kg/m^3)))
MatParam[4] = SetMaterialParams(Name="Lower Crust", Phase=3,
                            CreepLaws= LinearViscous(η=1e21Pas),
                            Density  = (ConstantDensity(), PT_Density()))

Mat_tup = Tuple(MatParam)

P,T = 1.,1.
@btime compute_density($Mat_tup, $P, $T)
=#