text:
using GeoPhyInv
using Statistics
#jl @init_parallel_stencil(2, false, Float32, 4)
#md using Plots; gr();

# ### Medium
medium = Medium(:acou_homo2D, 5); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)

# ### AGeom
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)

# ### SrcWav
tgrid = range(0.0, stop = 2.0, length = 2500); # generate a time grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:p]);
update!(srcwav, [:p], wav);

# ### Plotting the `vp` model
#md p1=heatmap(medium, :vp)
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1)

# ## Acoustic
# Now we have all the required variables to create `SeisForwExpt` object and
# prepare the forward modeling.
# While creating, we switched the `snaps_flag` on, and instruct recording field at
# `tsnaps`.
# Once the `Expt` object is created, do the modeling without *additional* any
# memory allocations using `update!`

pa = SeisForwExpt(
    FdtdAcou(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :p,
    tsnaps = [0.4, 0.5, 0.8],
    tgrid = tgrid,
    rfields = [:p],
    verbose = true,
);

# ### Update `pa` to perform modelling; prints time and allocations
@time update!(pa);

# ### Extracting snaps
snaps1 = pa[:snaps, 1]; # extracting snaps of the first supersource
snaps2 = pa[:snaps, 2]; # second supersource

# ### Plotting pressure snapshots after acoustic simulation

#md function plot_snapshots()
#md     p=[]
#md     for j in 1:2, i in 1:3
#md     	push!(p,heatmap(pa[:snaps, j][i], xlabel="x", ylabel="z", legend=:none))
#md     end
#md     Plots.plot(p..., layout=(2,3), size=(1200,600))
#md end
#md plot_snapshots()

# ## Elastic
# Similarly, we can do the elastic wave-equation modeling
medium = Medium(:elastic_homo2D, 5); # load a simple homogeneous elastic medium from the gallery
update!(medium, [:vp, :vs, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)
# Add source term on `:vz` grid
srcwav = SrcWav(tgrid, ageom, [:vz]);
update!(srcwav, [:vz], wav);
# Perform modelling and extract snaps
pa = SeisForwExpt(
    FdtdElastic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :vz,
    tsnaps = [0.4, 0.5, 0.8],
    rfields = [:vz],
    tgrid = tgrid,
    verbose = true,
)
@time update!(pa);

snaps = pa[:snaps, 1];
snaps = pa[:snaps, 2];

# ### Plotting is vs model
#md p1=heatmap(medium, :vs)
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1)

# ### Plotting snapshots after elastic simulation
#md plot_snapshots() # P and S waves!


text:
@time begin

include("./Initial.jl")
using .Initial
include("./Frame.jl")
using .Frame
include("./DB.jl")
using .DB

time = tick
const frames = [frame(1, deepcopy(nodes))]

function startCalculations()
    insertToDB(deepcopy(nodes), 1)
    for i in 2:totalTicks
        calculateFrame(nodes)
        # push!(frames, frame(i, deepcopy(nodes)))
        insertToDB(deepcopy(nodes), i)

        println("$i frames calculated")
    end
end

connectToDB()
startCalculations()

end

ext: jl

max_stars_repo_name: PseudoDevelopers/Physics-Simulations

max_issues_repo_name: PseudoDevelopers/Physics-Simulations

max_forks_repo_name: PseudoDevelopers/Physics-Simulations


text:
using Documenter, SoftPredicates

makedocs(;
    modules=[SoftPredicates],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zenna/SoftPredicates.jl/blob/{commit}{path}#L{line}",
    sitename="SoftPredicates.jl",
    authors="Zenna Tavares",
    assets=String[],
)

deploydocs(;
    repo="github.com/zenna/SoftPredicates.jl",
)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# NB: If expected_bins not given, assume uniform distribution over bins

# TODO: Add ci()
# TODO: Add one-sided hypothesis tests

immutable ChisqTest <: HypothesisTest
    test_name::ASCIIString
    # Chi-squared
    statistic::Float64
    parameter::Float64
    p_value::Float64
    confidence_interval::Vector{Float64}
    estimate::Float64
    null_value::Float64
    alternative_hypothesis::ASCIIString
    method::ASCIIString
    degrees_of_freedom::Float64
end

function test_statistic{S <: Real,
                        T <: Real}(::Type{ChisqTest},
                                   bins::Vector{S},
                                   expected_bins::Vector{T})
    n = length(bins)
    chisq = 0.0
    for i in 1:n
        if expected_bins[i] < 0.0
            throw(ArgumentError("Negative expected values are not allowed"))
        end
        chisq += (bins[i] - expected_bins[i])^2 / expected_bins[i]
    end
    return chisq
end

function test_statistic{T <: Real}(::Type{ChisqTest},
                                   bins::Vector{T})
    n, s = length(bins), sum(bins)
    val = s / n
    chisq = 0.0
    for i in 1:n
        chisq += (bins[i] - val)^2 / val
    end
    return chisq
end

function p_value{S <: Real,
                 T <: Real}(::Type{ChisqTest},
                            bins::Vector{S},
                            expected_bins::Vector{T})
    chisq = test_statistic(ChisqTest, bins, expected_bins)
    df = length(bins) - 1
    return ccdf(Chisq(df), chisq)
end

function p_value{T <: Real}(::Type{ChisqTest},
                            bins::Vector{T})
    chisq = test_statistic(ChisqTest, bins)
    df = length(bins) - 1
    return ccdf(Chisq(df), chisq)
end

function Base.run{S <: Real,
                  T <: Real}(::Type{ChisqTest},
                             bins::Vector{S},
                             expected_bins::Vector{T})
    n = length(bins)
    df = n - 1.0
    chisq = test_statistic(ChisqTest, bins, expected_bins)
    return ChisqTest("Chi-squared Test",
                     chisq,
                     NaN,
                     ccdf(Chisq(df), chisq),
                     [NaN, NaN],
                     NaN,
                     NaN,
                     "Two-sided",
                     "Chi-squared",
                     df)
end

function Base.run{T <: Real}(::Type{ChisqTest},
                             bins::Vector{T})
    n = length(bins)
    df = n - 1.0
    chisq = test_statistic(ChisqTest, bins)
    return ChisqTest("Chi-squared Test",
                     chisq,
                     NaN,
                     ccdf(Chisq(df), chisq),
                     [NaN, NaN],
                     NaN,
                     NaN,
                     "Two-sided",
                     "Chi-squared",
                     df)
end


text:
str = """
73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450"""

str = join(split(str), "")
l = map(x->parse(Int, x), collect(str))
println(l)

#kernel1 = [1,1,1]
#results1 = conv(l, kernel1)[2:end-2]
#println(results1) # can be used to find max sum

prods = [prod(l[i:i+3]) for i in 1:length(l)-3]
println(maximum(prods))

prods = [prod(l[i:i+12]) for i in 1:length(l)-12]
println(maximum(prods))

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
"""
    Fermi.Integrals

Module to compute and manage molecular integrals.
"""
module Integrals

using Fermi
using Fermi.Libcint
using Fermi.Options
using Fermi.GaussianBasis
using Fermi.Orbitals
using LinearAlgebra
using TensorOperations

import Base: getindex, setindex!, delete!, show

export IntegralHelper, delete!, mo_from_ao!, JKFIT, RIFIT, Chonky, AbstractDFERI, AbstractERI, SparseERI

"""
    Fermi.Integrals.AbstractERI

Abstract type common to all ERI types
"""
abstract type AbstractERI end

"""
    Fermi.Integrals.AbstractDFERI

Abstract type common to all density-fitted ERI types

# Struct tree

**AbstractDFERI** <: AbstractERI
"""
abstract type AbstractDFERI <: AbstractERI end

"""
    Fermi.Integrals.JKFIT

Concrete type representing a density-fitted ERI using a JKFIT auxiliar basis set.

# Fields

| Name       | Type       | Description |
|:-----------|:-----------|:----------------------------------------------------|
|`basisset`  | `BasisSet` | Fermi.GaussianBasis.BasisSet object for the auxiliar JK basis |

# Examples

The JKFIT structure is used to build a IntegralHelper
```julia
julia> jk = Fermi.Integrals.JKFIT()
julia> ints =
```

# Struct tree

JKFIT <: AbstractDFERI <: AbstractERI
"""
struct JKFIT <: AbstractDFERI
    basisset::BasisSet
end

function JKFIT(mol::Molecule = Molecule())

    auxjk = Options.get("jkfit")
    # If aux is auto, determine the aux basis from the basis
    if auxjk == "auto"
        basis = Options.get("basis")
        dunning_name = Regex("cc-pv.z")
        auxjk = occursin(dunning_name, basis) ? basis*"-jkfit" : "cc-pvqz-jkfit"
    end

    return JKFIT(mol, auxjk)
end

function JKFIT(mol::Molecule, basis::String)
    return JKFIT(BasisSet(mol, basis))
end

"""
    Fermi.Integrals.RIFIT

Concrete type representing a density-fitted ERI using a RIFIT auxiliar basis set.

# Fields

| Name       | Type       | Description |
|:-----------|:-----------|:----------------------------------------------------|
|`basisset`  | `BasisSet` | Fermi.GaussianBasis.BasisSet object for the auxiliar RI basis |

# Struct tree

RIFIT <: AbstractDFERI <: AbstractERI
"""
struct RIFIT <: AbstractDFERI
    basisset::BasisSet
end

function RIFIT(mol::Molecule = Molecule())

    auxri = Options.get("rifit")
    # If aux is auto, determine the aux basis from the basis
    if auxri == "auto"
        basis = Options.get("basis")
        dunning_name = Regex("cc-pv.z")
        auxri = occursin(dunning_name, basis) ? basis*"-rifit" : "cc-pvqz-rifit"
    end

    return RIFIT(mol, auxri)
end

function RIFIT(mol::Molecule, basis::String)
    return RIFIT(BasisSet(mol, basis))
end

struct Chonky <:AbstractERI end
struct SparseERI <: AbstractERI end

"""
    IntegralHelper

Manager for integrals computation and storage.
Accesss like a dictionary e.g.,
```
julia> ints = Fermi.Integrals.IntegralHelper()
julia> ints["S"] #Returns overlap matrix
```
A key is associated with each type of integral

    "S"           -> Overlap integral
    "T"           -> Electron kinetic energy integral
    "V"           -> Electron-nuclei attraction integral
    "ERI"         -> Electron repulsion integral

# Fields
    molecule                    Molecule object
    orbitals                    Orbitals used in the integral computation
    basis                       Basis set name
    cache                       Holds computed integrals
    eri_type                    Defines how electron repulsion integrals are handled
"""
struct IntegralHelper{T<:AbstractFloat,E<:AbstractERI,O<:AbstractOrbitals}
    molecule::Molecule
    orbitals::O
    basis::String
    cache::Dict{String,AbstractArray{T}}
    eri_type::E
end

function IntegralHelper(x...;k...)

    precision = Options.get("precision")
    if precision == "single"
        IntegralHelper{Float32}(x...; k...)
    elseif precision == "double"
        IntegralHelper{Float64}(x...; k...)
    else
        throw(FermiException("precision can only be `single` or `double`. Got $precision"))
    end
end

function IntegralHelper{T}(;molecule = Molecule(), orbitals = nothing,
                           basis = Options.get("basis"), eri_type=nothing) where T<:AbstractFloat

    if orbitals === nothing
        bset = BasisSet(molecule, basis)
        orbitals = AtomicOrbitals(bset)
    end

    # If the associated orbitals are AtomicOrbitals and DF is requested, JKFIT is set by default
    if Options.get("df") && orbitals isa AtomicOrbitals && eri_type === nothing
        eri_type = JKFIT(molecule)

    # If df is requested, but the orbitals are not AtomicOrbitals then RIFIT is set by default
    elseif Options.get("df") && eri_type === nothing

        eri_type = RIFIT(molecule)

    # Else, if eri_type is not an AbstractERI object, Chonky is used. This will overwrite invalid entries for eri_type
    elseif !(typeof(eri_type) <: AbstractERI)
        eri_type = orbitals isa AtomicOrbitals ? SparseERI() : Chonky()
    end

    cache = Dict{String, AbstractArray{T}}()
    IntegralHelper(molecule, orbitals, basis, cache, eri_type)
end

function IntegralHelper{T}(bset::BasisSet, eri_type=nothing) where T<:AbstractFloat
    IntegralHelper{T}(
        molecule = bset.molecule,
        orbitals = AtomicOrbitals(bset),
        basis = bset.basis_name,
        eri_type=eri_type
    )
end

"""
    getindex(I::IntegralHelper,entry::String)

Called when `helper["foo"]` syntax is used. If the requested entry already
exists, simply return the entry. If not, compute the requested entry.
"""
function getindex(I::IntegralHelper,entry::String)
    if haskey(I.cache, entry)
        return I.cache[entry]
    else
        compute!(I, entry)
        return I.cache[entry]
    end
end

function setindex!(I::IntegralHelper, A::AbstractArray, key::String)
    I.cache[key] = A
end

function delete!(I::IntegralHelper, keys...)
    for k in keys
        delete!(I.cache, k)
    end
    GC.gc()
end

function delete!(I::IntegralHelper)
    delete!(I, keys(I.cache)...)
end

# Integrals specific to orbital types
include("AtomicIntegrals.jl")
include("ROIntegrals.jl")

# Pretty printing
function string_repr(X::IntegralHelper{T,E,O}) where {T,E,O}
    eri_string = replace("$E", "Fermi.Integrals."=>"")
    orb_string = replace("$O", "Fermi.Orbitals."=>"")
    out = ""
    out = out*" ⇒ Fermi IntegralHelper\n"
    out = out*" ⋅ Data Type:                 $(T)\n"
    out = out*" ⋅ Basis:                     $(X.basis)\n"
    out = out*" ⋅ ERI:                       $(eri_string)\n"
    out = out*" ⋅ Orbitals:                  $(orb_string)\n"
    cache_str = ""
    for k in keys(X.cache)
        cache_str *= k*" "
    end
    out = out*" ⋅ Stored Integrals:          $(cache_str)"
    return out
end

function string_repr(X::AbstractDFERI)
    eri_string = replace("$(typeof(X))", "Fermi.Integrals."=>"")
    return "$(eri_string): $(X.basisset.basis_name)"
end

function show(io::IO, ::MIME"text/plain", X::T) where T<:Union{IntegralHelper,JKFIT,RIFIT}
    print(io, string_repr(X))
end

end # Module


text:
function libsvm2sparse(fname)
    ncols = xrows = yrows = 0
    ydict = Dict{Int32,Int32}()
    info("Sizing $fname...")
    f = open(fname)
    for l in eachline(f)
        ncols += 1
        a = split(l)
        class = parse(Int32,a[1])
        ydict[class] = 1
        for i=2:length(a)
            aa = split(a[i],':')
            idx = parse(Int32,aa[1])
            idx > xrows && (xrows = idx)
        end
    end
    close(f)
    ykeys = sort(collect(keys(ydict)))
    yrows = length(ykeys)
    for i in 1:yrows
        ydict[ykeys[i]] = i
    end
    info("Found $ncols instances with xrows=$xrows yrows=$yrows")

    x = spzeros(Float32,Int32,xrows,ncols)
    y = zeros(Float32,yrows,ncols)

    info("Reading $fname...")
    f = open(fname)
    ncols = 0
    for l in eachline(f)
        ncols += 1
        a = split(l)
        class = parse(Int32,a[1])
        y[ydict[class], ncols] = 1
        for i=2:length(a)
            aa = split(a[i],':')
            idx = parse(Int32,aa[1])
            val = parse(Float32,aa[2])
            x[idx, ncols] = 1
        end
    end
    close(f)
    info("done.")
    return (x,y)
end


text:
using DrWatson
@quickactivate "VMC"
DrWatson.greet()

max_stars_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_issues_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_forks_repo_name: Ved-Mahajan/Variational-Monto-Carlo


text:
@testset "Testing Record API" begin
    # Test some basic `mapleaves` functionality
    mapleaves = CounterTools.mapleaves
    aggregate = CounterTools.aggregate
    Record = CounterTools.Record

    x = 1
    f = x -> 2x
    @test mapleaves(f, 1) == 2
    @test mapleaves(f, (1,2,3)) == (2, 4, 6)
    @test mapleaves(f, ((1, 2), (3, 4))) == ((2, 4), (6, 8))
    @test mapleaves(f, ([1, 2], [3, 4])) == ([2, 4], [6, 8])

    # Test `aggregate` of non-record types
    @test aggregate(1) == 1
    @test aggregate((1,2,3)) == (1,2,3)
    @test aggregate([1,2,3]) == 6
    @test aggregate(((1,2), (3,4), (5,6))) == (9,12)
    @test aggregate([(1,2), (3,4), (5,6)]) == (9,12)

    @test aggregate((x,y) -> max.(x, y), ((10,2), (3,4), (5,6))) == (10,6)

    # Start building some records
    a = Record{:A}((1,2,3))
    @test CounterTools.name(a) == :A
    @test CounterTools.hassubrecord(a) == false
    @test a[1] == 1
    @test a[2] == 2
    @test a[3] == 3

    # mapleaves
    b = mapleaves(f, a)
    @test isa(b, Record)
    @test CounterTools.name(b) == :A
    @test b.data == (2,4,6)

    # Test subtraction
    c = b - a
    @test isa(c, Record)
    @test CounterTools.name(b) == :A
    @test c.data == (1,2,3)

    # Single record aggregation
    @test CounterTools.aggregate(a) == (1,2,3)

    # Dump to devnull to make sure this shows
    show(devnull, a)

    ##### Test over a vector
    data = [
        (1,2,3),
        (4,5,6),
        (7,8,9),
    ]
    a = Record{:A}(data)
    @test CounterTools.aggregate(a) == (12,15,18)

    ##### Work with recursive data structures now
    a = Record{:A}((
        (1,2,3),
        (4,5,6),
    ))

    b = Record{:A}((
        (7,8,9),
        (10,11,12),
    ))

    c = Record{:B}((a,b))

    # Try mapleaves
    d = mapleaves(f, c)
    @test isa(d, Record)
    @test CounterTools.name(d) == :B
    @test CounterTools.hassubrecord(d) == true
    @test isa(d[1], Record)
    @test CounterTools.name(d[1]) == :A

    @test length(d) == 2
    @test length(d[1]) == 2

    @test d[1][1] == (2,4,6)
    @test d[1][2] == (8,10,12)
    @test d[2][1] == (14,16,18)
    @test d[2][2] == (20,22,24)

    e = d - c
    @test e == c
    @test CounterTools.aggregate(e) == (1,2,3) .+ (4,5,6) .+ (7,8,9) .+ (10,11,12)
    @test CounterTools.aggregate(d) == 2 .* CounterTools.aggregate(e)
    show(devnull, e)
end


text:
"""
    getdriver(i::Integer)

Fetch driver by index.
"""
getdriver(i::Integer)::Driver = Driver(GDAL.gdalgetdriver(i))

"""
    getdriver(name::AbstractString)

Fetch a driver based on the short name (such as `GTiff`).
"""
getdriver(name::AbstractString)::Driver = Driver(GDAL.gdalgetdriverbyname(name))

"""
    register(drv::Driver)

Register a driver for use.
"""
function register(drv::Driver)::Nothing
    GDAL.gdalregisterdriver(drv.ptr)
    return nothing
end

"""
    deregister(drv::Driver)

Deregister the passed driver.
"""
function deregister(drv::Driver)::Nothing
    GDAL.gdalderegisterdriver(drv.ptr)
    return nothing
end

"""
    options(drv::Driver)

Return the list of creation options of the driver [an XML string].
"""
options(drv::Driver)::String = GDAL.gdalgetdrivercreationoptionlist(drv.ptr)

driveroptions(name::AbstractString)::String = options(getdriver(name))

"""
    shortname(drv::Driver)

Return the short name of a driver (e.g. `GTiff`).
"""
shortname(drv::Driver)::String = GDAL.gdalgetdrivershortname(drv.ptr)

"""
    longname(drv::Driver)

Return the long name of a driver (e.g. `GeoTIFF`), or empty string.
"""
longname(drv::Driver)::String = GDAL.gdalgetdriverlongname(drv.ptr)

"""
    ndriver()

Fetch the number of registered drivers.
"""
ndriver()::Integer = GDAL.gdalgetdrivercount()

"""
    listdrivers()

Returns a listing of all registered drivers.
"""
listdrivers()::Dict{String,String} = Dict{String,String}([
    shortname(getdriver(i)) => longname(getdriver(i)) for i in 0:(ndriver()-1)
],)

"""
    identifydriver(filename::AbstractString)

Identify the driver that can open a raster file.

This function will try to identify the driver that can open the passed filename
by invoking the Identify method of each registered `Driver` in turn. The
first driver that successful identifies the file name will be returned. If all
drivers fail then `NULL` is returned.
"""
identifydriver(filename::AbstractString)::Driver =
    Driver(GDAL.gdalidentifydriver(filename, C_NULL))

"""
    validate(drv::Driver, options::Vector{<:AbstractString})

Validate the list of creation options that are handled by a drv.

This is a helper method primarily used by `create()` and `copy()` to
validate that the passed in list of creation options is compatible with the
`GDAL_DMD_CREATIONOPTIONLIST` metadata item defined by some drivers.

### Parameters
* `drv`     the handle of the driver with whom the lists of creation option
            must be validated
* `options` the list of creation options. An array of strings, whose last
            element is a `NULL` pointer

### Returns
`true` if the list of creation options is compatible with the `create()` and
`createcopy()` method of the driver, `false` otherwise.

### Additional Remarks
See also: `options(drv::Driver)`

If the `GDAL_DMD_CREATIONOPTIONLIST` metadata item is not defined, this
function will return ``true``. Otherwise it will check that the keys and values
in the list of creation options are compatible with the capabilities declared
by the `GDAL_DMD_CREATIONOPTIONLIST` metadata item. In case of incompatibility
a (non fatal) warning will be emited and ``false`` will be returned.
"""
function validate(
    drv::Driver,
    options::Vector{T},
)::Bool where {T<:AbstractString}
    return Bool(GDAL.gdalvalidatecreationoptions(drv.ptr, options))
end

"""
    copyfiles(drv::Driver, new::AbstractString, old::AbstractString)
    copyfiles(drvname::AbstractString, new::AbstractString, old::AbstractString)

Copy all the files associated with a dataset.
"""
function copyfiles end

function copyfiles(
    drv::Driver,
    new::AbstractString,
    old::AbstractString,
)::Nothing
    result = GDAL.gdalcopydatasetfiles(drv.ptr, new, old)
    @cplerr result "Failed to copy dataset files"
    return nothing
end

function copyfiles(
    drvname::AbstractString,
    new::AbstractString,
    old::AbstractString,
)::Nothing
    copyfiles(getdriver(drvname), new, old)
    return nothing
end

"""
    extensions()

Returns a `Dict{String,String}` of all of the file extensions that can be read
by GDAL,  with their respective drivers' `shortname`s.
"""
function extensions()::Dict{String,String}
    extdict = Dict{String,String}()
    for i in 1:ndriver()
        driver = getdriver(i)
        if !(driver.ptr == C_NULL)
            # exts is a space-delimited list in a String, so split it
            for ext in split(metadataitem(driver, "DMD_EXTENSIONS"))
                extdict[".$ext"] = shortname(driver)
            end
        end
    end
    return extdict
end

"""
    extensiondriver(filename::AbstractString)

Returns a driver shortname that matches the filename extension.

So `extensiondriver("/my/file.tif") == "GTiff"`.
"""
function extensiondriver(filename::AbstractString)::String
    split = splitext(filename)
    extensiondict = extensions()
    ext = split[2] == "" ? split[1] : split[2]
    if !haskey(extensiondict, ext)
        throw(ArgumentError("There are no GDAL drivers for the $ext extension"))
    end
    return extensiondict[ext]
end


text:
const JRLL = [ 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 ]
const JPLL = [ 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 ]

const UTAB = [0, 1, 4, 5,16,17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 85, 256,
              257, 260, 261, 272, 273, 276, 277, 320, 321, 324, 325, 336, 337,
              340, 341, 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045, 1088,
              1089, 1092, 1093, 1104, 1105, 1108, 1109, 1280, 1281, 1284, 1285,
              1296, 1297, 1300, 1301, 1344, 1345, 1348, 1349, 1360, 1361, 1364,
              1365, 4096, 4097, 4100, 4101, 4112, 4113, 4116, 4117, 4160, 4161,
              4164, 4165, 4176, 4177, 4180, 4181, 4352, 4353, 4356, 4357, 4368,
              4369, 4372, 4373, 4416, 4417, 4420, 4421, 4432, 4433, 4436, 4437,
              5120, 5121, 5124, 5125, 5136, 5137, 5140, 5141, 5184, 5185, 5188,
              5189, 5200, 5201, 5204, 5205, 5376, 5377, 5380, 5381, 5392, 5393,
              5396, 5397, 5440, 5441, 5444, 5445, 5456, 5457, 5460, 5461, 16384,
              16385, 16388, 16389, 16400, 16401, 16404, 16405, 16448, 16449,
              16452, 16453, 16464, 16465, 16468, 16469, 16640, 16641, 16644,
              16645, 16656, 16657, 16660, 16661, 16704, 16705, 16708, 16709,
              16720, 16721, 16724, 16725, 17408, 17409, 17412, 17413, 17424,
              17425, 17428, 17429, 17472, 17473, 17476, 17477, 17488, 17489,
              17492, 17493, 17664, 17665, 17668, 17669, 17680, 17681, 17684,
              17685, 17728, 17729, 17732, 17733, 17744, 17745, 17748, 17749,
              20480, 20481, 20484, 20485, 20496, 20497, 20500, 20501, 20544,
              20545, 20548, 20549, 20560, 20561, 20564, 20565, 20736, 20737,
              20740, 20741, 20752, 20753, 20756, 20757, 20800, 20801, 20804,
              20805, 20816, 20817, 20820, 20821, 21504, 21505, 21508, 21509,
              21520, 21521, 21524, 21525, 21568, 21569, 21572, 21573, 21584,
              21585, 21588, 21589, 21760, 21761, 21764, 21765, 21776, 21777,
              21780, 21781, 21824, 21825, 21828, 21829, 21840, 21841, 21844,
              21845]

const CTAB = [0, 1, 256, 257, 2, 3, 258, 259, 512, 513, 768, 769, 514, 515, 770,
              771, 4, 5, 260, 261, 6, 7, 262, 263, 516, 517, 772, 773, 518, 519,
              774, 775, 1024, 1025, 1280, 1281, 1026, 1027, 1282, 1283, 1536,
              1537, 1792, 1793, 1538, 1539, 1794, 1795, 1028, 1029, 1284, 1285,
              1030, 1031, 1286, 1287, 1540, 1541, 1796, 1797, 1542, 1543, 1798,
              1799, 8, 9, 264, 265, 10, 11, 266, 267, 520, 521, 776, 777, 522,
              523, 778, 779, 12, 13, 268, 269, 14, 15, 270, 271, 524, 525, 780,
              781, 526, 527, 782, 783, 1032, 1033, 1288, 1289, 1034, 1035, 1290,
              1291, 1544, 1545, 1800, 1801, 1546, 1547, 1802, 1803, 1036, 1037,
              1292, 1293, 1038, 1039, 1294, 1295, 1548, 1549, 1804, 1805, 1550,
              1551, 1806, 1807, 2048, 2049, 2304, 2305, 2050, 2051, 2306, 2307,
              2560, 2561, 2816, 2817, 2562, 2563, 2818, 2819, 2052, 2053, 2308,
              2309, 2054, 2055, 2310, 2311, 2564, 2565, 2820, 2821, 2566, 2567,
              2822, 2823, 3072, 3073, 3328, 3329, 3074, 3075, 3330, 3331, 3584,
              3585, 3840, 3841, 3586, 3587, 3842, 3843, 3076, 3077, 3332, 3333,
              3078, 3079, 3334, 3335, 3588, 3589, 3844, 3845, 3590, 3591, 3846,
              3847, 2056, 2057, 2312, 2313, 2058, 2059, 2314, 2315, 2568, 2569,
              2824, 2825, 2570, 2571, 2826, 2827, 2060, 2061, 2316, 2317, 2062,
              2063, 2318, 2319, 2572, 2573, 2828, 2829, 2574, 2575, 2830, 2831,
              3080, 3081, 3336, 3337, 3082, 3083, 3338, 3339, 3592, 3593, 3848,
              3849, 3594, 3595, 3850, 3851, 3084, 3085, 3340, 3341, 3086, 3087,
              3342, 3343, 3596, 3597, 3852, 3853, 3598, 3599, 3854, 3855]

@doc raw"""
    pix2xyfRing(resol::Resolution, ipix) :: (Int, Int, Int)

Convert a pixel number into (x, y, face), using RING ordering."""
function pix2xyfRing(resol::Resolution, ipix)
    if ipix < resol.ncap
        # North polar cap
        iring = (1 + isqrt(2 * ipix - 1)) >> 1
        iphi = ipix - 2 * iring * (iring - 1)
        kshift = 0
        nr = iring
        facenum = (iphi - 1) ÷ nr
    elseif ipix < (resol.numOfPixels - resol.ncap)
        # Equatorial region
        ip = ipix - 1 - resol.ncap
        tmp = (resol.order ≥ 0) ? (ip >> (resol.order + 2)) : (ip ÷ resol.nsideTimesFour)
        iring = tmp + resol.nside
        iphi = ip - tmp * resol.nsideTimesFour + 1
        kshift = (iring + resol.nside) & 1
        nr = resol.nside
        ire = iring - resol.nside + 1
        irm = resol.nsideTimesTwo + 2 - ire
        ifm = iphi - ire ÷ 2 + resol.nside - 1
        ifp = iphi - irm ÷ 2 + resol.nside - 1
        if resol.order ≥ 0
            ifm >>= resol.order
            ifp >>= resol.order
        else
            ifm ÷= resol.nside
            ifp ÷= resol.nside
        end

        facenum = if ifp == ifm
            ifp | 4
        else
            if ifp < ifm
                ifp
            else
                ifm + 8
            end
        end
    else
        # South polar cap
        ip = resol.numOfPixels - ipix + 1
        iring = (1 + isqrt(2 * ip - 1)) >> 1
        iphi = 4 * iring + 1 - (ip - 2 * iring * (iring - 1))
        kshift = 0
        nr = iring
        iring = 2 * resol.nsideTimesTwo - iring
        facenum = 8 + (iphi - 1) ÷ nr
    end

    irt = iring - JRLL[facenum + 1] * resol.nside + 1
    ipt = 2 * iphi - JPLL[facenum + 1] * nr - kshift - 1
    ipt ≥ resol.nsideTimesTwo && (ipt -= 8 * resol.nside)

    ((ipt - irt) >> 1, (-ipt - irt) >> 1, facenum)
end

@doc raw"""
    xyf2pixRing(resol::Resolution, ix, iy, facenum) :: Int

Convert (x, y, face) into a pixel number, using RING ordering."""
function xyf2pixRing(resol::Resolution, ix, iy, facenum)
    jr = JRLL[facenum + 1] * resol.nside - ix - iy - 1

    ringinfo = getringinfo(resol, jr, full=false)
    nr = ringinfo.numOfPixels >> 2
    kshift = ringinfo.shifted ? 0 : 1
    jp = (JPLL[facenum + 1] * nr + ix - iy + 1 + kshift) ÷ 2
    jp < 1 && (jp += resol.nsideTimesFour)

    (ringinfo.firstPixIdx - 1) + jp
end

function spreadbits(v::Int)
    (Int64(UTAB[v & 0xff + 1]) |
     Int64(UTAB[(v >> 8) & 0xff + 1] << 16) |
     Int64(UTAB[(v >> 16) & 0xff + 1] << 32) |
     Int64(UTAB[(v >> 24) & 0xff + 1] << 48))
end

function compress_bits(v::Int32)
    raw = (v & 0x5555) | ((v & 0x55550000) >> 15)
    CTAB[raw & 0xff + 1] | (ctab[raw >> 8 + 1] << 4)
end

function compress_bits(v::Int)
    raw = v & 0x5555555555555555
    raw |= raw >> 15
    return (CTAB[raw & 0xff + 1] | (CTAB[(raw >> 8) & 0xff + 1] << 4) |
            (CTAB[(raw >> 32) & 0xff + 1] << 16) |
            (CTAB[(raw >> 40) & 0xff + 1] << 20))
end

@doc raw"""
    pix2xyfNest(resol::Resolution, ipix) :: (Int, Int, Int)

Convert a pixel number into (x, y, face), using NESTED ordering."""
function pix2xyfNest(resol::Resolution, ipix)
    pix = (ipix - 1) & (resol.pixelsPerFace - 1)
    (compress_bits(pix), compress_bits(pix ÷ 2), (ipix - 1) >> (2 * resol.order))
end

@doc raw"""
    xyf2pixNest(resol::Resolution, ix, iy, facenum) :: Int

Convert (x, y, face) into a pixel number, using NESTED ordering."""
function xyf2pixNest(resol::Resolution, ix, iy, facenum)
    facenum << (2 * resol.order) + spreadbits(ix) + 2spreadbits(iy) + 1
end

@doc raw"""
    ring2nest(resol::Resolution, ipix) :: Int

Convert the number of a pixel from RING to NESTED scheme."""
ring2nest(resol::Resolution, ipix) = xyf2pixNest(resol, pix2xyfRing(resol, ipix)...)

@doc raw"""
    nest2ring(resol::Resolution, ipix) :: Int

Convert the number of a pixel from NESTED to RING scheme."""
nest2ring(resol::Resolution, ipix) = xyf2pixRing(resol, pix2xyfNest(resol, ipix)...)


text:
# # BRTs for species distribution forecasting

using SimpleSDMLayers
using EvoTrees
using GBIF
using StatsBase
using StatsPlots

# **Justification for this use case:** Boosted Regression Trees (BRTs) are a
# powerful way to predict the distribution of species. We will see how we can
# get information in and out of layers to use them, and how to use this model to
# predict a new distribution under a climate change scenario. This use-case
# assumes that you have read the manual pages for GBIF integration, future data,
# and pseudo-absences generation.

# We will re-use the data from the pseudo-absences example:

sp = GBIF.taxon("Hypomyces lactifluorum")
observations = occurrences(
    sp, "hasCoordinate" => true, "limit" => 300, "country" => "CA", "country" => "US"
)
while length(observations) < size(observations)
    occurrences!(observations)
end

# We will pick the entire BioClim layers at a 10 minutes resolution, and clip
# them to the observations (this adds a 5 degrees buffer).

layers = [
    clip(layer, observations) for layer in SimpleSDMPredictor(WorldClim, BioClim, 1:19)
];

# To remove the sampling effect, we transform the presences to a grid, and
# generate pseudo-absences using the surface range envelope method.

presences = mask(layers[1], observations, Bool)
absences = rand(SurfaceRangeEnvelope, presences)

# The next step is to extract coordinates at which the species is present or
# pseudo-absent - we can rely on the `replace` method to empty any `false`
# values:

xy_presence = keys(replace(presences, false => nothing));
xy_absence = keys(replace(absences, false => nothing));
xy = vcat(xy_presence, xy_absence);

# With the `xy` list of coordinates, we can get a predictor `X`, and a response
# `y`.

X = hcat([layer[xy] for layer in layers]...);
y = vcat(fill(1.0, length(xy_presence)), fill(0.0, length(xy_absence)));

# To train the model, we will use a random subset representing 70% of the
# dataset:

train_size = floor(Int, 0.7 * length(y));
train_idx = sample(1:length(y), train_size; replace=false);
test_idx = setdiff(1:length(y), train_idx);

# This gives use the training and testing (or evaluation) sets:

Xtrain, Xtest = X[train_idx, :], X[test_idx, :];
Ytrain, Ytest = y[train_idx], y[test_idx];

# In order to fit the tree, we need to define a number of parameters. We will
# use a Gaussian maximum likelihood tree (from `EvoTrees`), which will give us a
# measure of the average prediction, but also the standard deviation. This is
# important in order to communicate uncertainty.

gaussian_tree_parameters = EvoTreeGaussian(;
    loss=:gaussian,
    metric=:gaussian,
    nrounds=100,
    nbins=100,
    λ=0.0,
    γ=0.0,
    η=0.1,
    max_depth=7,
    min_weight=1.0,
    rowsample=0.5,
    colsample=1.0,
)

# We can now fit the BRT. This function has an additional `print_every_n` to
# update on the progress every `n` epochs, but we don't really need it here.

model = fit_evotree(gaussian_tree_parameters, Xtrain, Ytrain; X_eval=Xtest, Y_eval=Ytest)

# The next step is to gather *all* the values of all the layers in a matrix, in
# order to run the full spatial prediction:

all_values = hcat([layer[keys(layer)] for layer in layers]...);

# If the matrix is too big, we could resort to a combination of `clip`, make the
# prediction on each tile, and then use `hcat` and `vcat` to combine them. This
# is not the case here, so we can predict directly:

pred = EvoTrees.predict(model, all_values);

# Once the prediction is done, we can copy the values into a layer.

distribution = similar(layers[1], Float64)
distribution[keys(distribution)] = pred[:, 1]
distribution

# The BRT is able to calculate a measure of relative gain from the different
# variables:

top5_var = importance(model, collect(layernames(WorldClim, BioClim)))[1:5]

# This is an interesting alternative to VIF for variable selection. Let's
# examine how the most important variable relates to the predicted distribution
# score:

most_important_layer = findfirst(isequal(top5_var[1].first), collect(layernames(WorldClim, BioClim)))
histogram(
    layers[most_important_layer][xy_presence]; fill=(0, :teal, 0.2), lc=:teal, frame=:origin, lab="Present"
)
histogram!(
    layers[most_important_layer][xy_absence]; fill=(0, :white, 0.0), frame=:origin, lc=:grey, lab="Absent"
)
xaxis!(layernames(WorldClim, BioClim, most_important_layer))

# It is interesting to notice that despite the importance of this predictor, the
# difference between the presence and absence locations are not as clear as we
# may expect!

# We can similarly extract uncertainty:

uncertainty = similar(layers[1], Float64)
uncertainty[keys(uncertainty)] = pred[:, 2]
uncertainty

# And we can now visualize the prediction, which we force to be in `[0,1]`.

p_dis = plot(rescale(distribution, (0, 1)); c=:bamako, frame=:box)
scatter!(xy_presence; lab="", c=:black, alpha=0.2, msw=0.0, ms=3)

# We can do the same thing for the uncertainty

p_unc = plot(uncertainty; c=:tokyo, frame=:box)

# Of course, this prediction is returing values for the entire range of the
# initial layer, so let's compare the distributions of the prediction score:

histogram(
    distribution[xy_presence]; fill=(0, :teal, 0.2), lc=:teal, frame=:origin, lab="Present"
)
histogram!(
    distribution[xy_absence]; fill=(0, :white, 0.0), frame=:origin, lc=:grey, lab="Absent"
)
xaxis!("Prediction score")

# This looks like a good opportunity to do some thresholding. Note that the
# values are *not* moved back to the unit range, because we'll need the raw
# values for a little surprise later on. We will find the value of the score
# that optimizes Youden's J (Cohen's κ is also a suitable alternative):

cutoff = LinRange(extrema(distribution)..., 500);

obs = y .> 0

tp = zeros(Float64, length(cutoff));
fp = zeros(Float64, length(cutoff));
tn = zeros(Float64, length(cutoff));
fn = zeros(Float64, length(cutoff));

for (i, c) in enumerate(cutoff)
    prd = distribution[xy] .>= c
    tp[i] = sum(prd .& obs)
    tn[i] = sum(.!(prd) .& (.!obs))
    fp[i] = sum(prd .& (.!obs))
    fn[i] = sum(.!(prd) .& obs)
end

# From this, we can calculate a number of validation measures:

tpr = tp ./ (tp .+ fn);
fpr = fp ./ (fp .+ tn);
J = (tp ./ (tp .+ fn)) + (tn ./ (tn .+ fp)) .- 1.0;
ppv = tp ./ (tp .+ fp);

# The ROC-AUC is an overall measure of how good the fit is:

dx = [reverse(fpr)[i] - reverse(fpr)[i - 1] for i in 2:length(fpr)]
dy = [reverse(tpr)[i] + reverse(tpr)[i - 1] for i in 2:length(tpr)]
AUC = sum(dx .* (dy ./ 2.0))

# We can pick the value of the cutoff that maximizes J:

thr_index = last(findmax(J))
τ = cutoff[thr_index]

# Let's have a look at the ROC curve:

plot(fpr, tpr; aspectratio=1, frame=:box, lab="", dpi=600, size=(400, 400))
scatter!([fpr[thr_index]], [tpr[thr_index]]; lab="", c=:black)
plot!([0, 1], [0, 1]; c=:grey, ls=:dash, lab="")
xaxis!("False positive rate", (0, 1))
yaxis!("True positive rate", (0, 1))

# And the precision-recall as well:

plot(tpr, ppv; aspectratio=1, frame=:box, lab="", dpi=600, size=(400, 400))
scatter!([tpr[thr_index]], [ppv[thr_index]]; lab="", c=:black)
plot!([0, 1], [1, 0]; c=:grey, ls=:dash, lab="")
xaxis!("True positive rate", (0, 1))
yaxis!("Positive predictive value", (0, 1))

# We can now map the result using τ as a cutoff for the `distribution` data:

range_mask = broadcast(v -> v >= τ, distribution)

# And finally, plot the whole thing:

plot(distribution; c=:lightgrey, leg=false)
plot!(mask(range_mask, distribution); c=:darkgreen)
scatter!(xy_presence; lab="", c=:orange, alpha=0.5, msw=0.0, ms=2)

# Because our BRT also returns the uncertainty, we can combine both maps into a
# bivariate one, showing both where we expect the species, but also where we are
# uncertain about the prediction:

plot(distribution; leg=false, c=:lightgrey, frame=:grid, xlab="Longitude", ylab="Latitude", grid=false)
bivariate!(mask(range_mask, distribution), mask(range_mask, uncertainty))
p2 = bivariatelegend!(
    mask(range_mask, distribution),
    mask(range_mask, uncertainty);
    inset=(1, bbox(0.04, 0.08, 0.23, 0.23, :center, :left)),
    subplot=2,
    xlab="Prediction",
    ylab="Uncertainty",
    guidefontsize=7,
)

# Now, for the big question - will this range move in the future? To explore
# this, we will get the same variables, but in the future. In order to simplify
# the code, we will limit ourselves to one SSP (585) and one CMIP6 model
# (CanESM5), around 2050:

future_layers = [
    clip(layer, observations) for
    layer in SimpleSDMPredictor(WorldClim, BioClim, CanESM5, SSP585, 1:19; year="2041-2060")
];

# We can get all the future values from this data:

all_future_values = hcat([layer[keys(layer)] for layer in future_layers]...);

# And make a prediction based on our BRT model. This is, of course, assuming
# that BRTs are good at this type of prediction (they're OK).

future_pred = EvoTrees.predict(model, all_future_values);

# As before, we also have a measure of uncertainty. In the interest of keeping
# this vignette small, we will not look at it.

future_distribution = similar(layers[1], Float64)
future_distribution[keys(future_distribution)] = future_pred[:, 1]
future_distribution

# The values in `future_distribution` are in the scale of what the BRT returns,
# so we can compare them with the values of `distribution`:

plot(future_distribution - distribution; clim=(-1, 1), c=:broc, frame=:box)

# This shows the area of predicted gain and loss of presence. Because we have
# thresholded our current distribution, we can look at the predicted ranges of
# suitability:

future_range_mask = broadcast(v -> v >= τ, future_distribution)

# The last step is to get the difference between the future and current masks
# (so +1 is a gain of range, 0 is no change, and -1 is a loss), and to only
# report this for the cells that are both in the current and future data:

range_change = convert(Float32, future_range_mask) - convert(Float32, range_mask)
both_ranges_mask = maximum([future_range_mask, range_mask])

# We can now plot the result, with the brown area being range that becomes
# unfavorable, the green one remaining suitable, and the blue area being newly
# opened range:

plot(distribution; c=:lightgrey, leg=false)
plot!(mask(both_ranges_mask, range_change); c=:roma)


text:
import GR

GR.selntran(0)
GR.settextfontprec(2, 0)
GR.setcharheight(0.024)

GR.text(0.05, 0.85, " !\"#\$%&'()*+,-./")
GR.text(0.05, 0.80, "0123456789:;<=>?")
GR.text(0.05, 0.75, "@ABCDEFGHIJKLMNO")
GR.text(0.05, 0.70, "PQRSTUVWXYZ[\\]^_")
GR.text(0.05, 0.65, "`abcdefghijklmno")
GR.text(0.05, 0.60, "pqrstuvwxyz{|}~ ")

GR.text(0.5, 0.85, " ¡¢£¤¥¦§¨©ª«¬®¯")
GR.text(0.5, 0.80, "°±²³´µ¶·¸¹º»¼½¾¿")
GR.text(0.5, 0.75, "ÀÁÂÃÄÅÆÇÈÉÊËÌÍÎÏ")
GR.text(0.5, 0.70, "ÐÑÒÓÔÕÖ×ØÙÚÛÜÝÞß")
GR.text(0.5, 0.65, "àáâãäåæçèéêëìíîï")
GR.text(0.5, 0.60, "ðñòóôõö÷øùúûüýþÿ")

GR.updatews()


text:
module RandSeis
using SeisIO, SeisIO.Quake
using Random:randexp, randstring, shuffle!
import Dates:now
import SeisIO:InstrumentPosition, code2typ, getbandcode, note!, sμ, μs
export randSeisChannel, randSeisData, randSeisEvent, randSeisHdr, randPhaseCat,
randSeisSrc

include("RandSeis/randseis_utils.jl")
include("RandSeis/randSeisChannel.jl")
include("RandSeis/randSeisData.jl")
include("RandSeis/randSeisEvent.jl")

const fc_vals = Float64[1.0/120.0 1.0/60.0 1.0/30.0 0.2 1.0 1.0 1.0 2.0 4.5 15.0]
const fs_vals = Float64[0.1, 1.0, 2.0, 5.0, 10.0, 20.0, 25.0, 40.0, 50.0, 60.0, 62.5, 80.0, 100.0, 120.0, 125.0, 250.0]
const irregular_units = ["%", "%{cloud_cover}", "{direction_vector}", "Cel", "K", "{none}", "Pa", "T", "V", "W", "m", "m/m", "m/s", "m/s2", "m3/m3", "rad", "rad/s", "rad/s2", "t{SO_2}"]

end


text:
"""
Max-tree morphological representation of an image.

# Details
Let's consider a *thresholding* operation,
```julia
    mask = [val ≥ threshold for val in image]
```
One can identify the connected components (the sets of neighboring true values)
in `mask`. When *image thresholding* is sequentially applied for all possible
thresholds, it generates a collection of connected components that could be
organized into a hierarchical structure called *component tree*. Consider 1D
"image" with values 1, 2 and 3:
```
       2233233312223322
```
The connected components would be
```
    1: AAAAAAAAAAAAAAAA
    2: BBBBBBBB.CCCCCCC
    3: ..DD.EEE....FF..
```
Here, the letters are the labels of the resulting connected components,
and `.` specifies that the pixel value is below the threshold.
In this example, the corresponding *component tree* is:
```
      A
     ⭩ ⭨
    B   C
   ⭩ ⭨   ⭨
  D   E   F
```

A *max-tree* is an efficient representation of the *component tree*.
A connected component ``C`` at threshold level ``t`` is represented by the
single *reference pixel* ``r`` from this level (`image[r] == t`), which is the
parent to all other pixels of ``C`` and also to the *reference pixels* of the
connected components at higher thresholds, which are the children of ``C``.
In our example, the reference pixels (denoted by the letter of the corresponding
component) would be:
```
    1: ........A.......
    2: B........C......
    3: ..D..E......F...
```
I.e.

| Comp | Ref.Pixel |
| ---- | ---------:|
| *A* |  9 |
| *B* |  1 |
| *C* | 10 |
| *D* |  3 |
| *E* |  6 |
| *F* | 13 |

So the whole max-tree could be encoded as a vector of indices of parent pixels:
```
9  1  1  3  1  1  6  6  9  9 10 10 10 13 10 10
```

The *max-tree* is the basis for many morphological operators,
namely connected operators. Unlike morphological openings and closings, these
operators do not require a fixed structuring element, but rather act with a
flexible structuring element that meets a certain criterion.

# See also
[`area_opening`](@ref), [`area_closing`](@ref),
[`diameter_opening`](@ref), [`diameter_closing`](@ref).

# References
1. Salembier, P., Oliveras, A., & Garrido, L. (1998). *Antiextensive Connected Operators for Image and Sequence Processing*. IEEE Transactions on Image Processing, 7(4), 555-570.
   > https://doi.org/10.1109/83.663500
2. Berger, C., Geraud, T., Levillain, R., Widynski, N., Baillard, A., Bertin, E. (2007). *Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging*. In International Conference on Image Processing (ICIP), 41-44.
   > https://doi.org/10.1109/ICIP.2007.4379949
3. Najman, L., & Couprie, M. (2006). *Building the component tree in quasi-linear time*. IEEE Transactions on Image Processing, 15(11), 3531-3539.
   > https://doi.org/10.1109/TIP.2006.877518
4. Carlinet, E., & Geraud, T. (2014). *A Comparative Review of Component Tree Computation Algorithms*. IEEE Transactions on Image Processing, 23(9), 3885-3895.
   > https://doi.org/10.1109/TIP.2014.2336551
"""
struct MaxTree{N}
    """
    If `false`, the tree is constructed from the smallest to the largest values,
    the tree root being the smallest (darkest) pixel.
    Otherwise, if `true`, the direction is from large to small values.
    """
    rev::Bool

    """
    Array of the same shape as the source image.
    Each value is the linear index of the parent pixel in the max-tree.
    """
    parentindices::Array{Int,N}

    """
    The order of the elements in the tree, from root to leaves.
    Each element corresponds to the linear index of the pixel in the image.
    A parent of the pixel always comes before the pixel itself, i.e.
    `parentindices[traverse[1]] == traverse[1]` (the root pixel) and
    `parentindices[traverse[k]] ∈ traverse[1:k-1]` for all `k >= 2`.
    """
    traverse::Vector{Int}

    # create uninitialized MaxTree for image
    MaxTree{N}(image::GenericGrayImage{<:Any, N}, rev::Bool) where N =
        new{N}(rev, Array{Int,N}(undef, size(image)),
               Vector{Int}(undef, length(image)))
end

Base.ndims(::Type{<:MaxTree{N}}) where N = N
Base.ndims(maxtree::MaxTree) = ndims(typeof(maxtree))
Base.length(maxtree::MaxTree) = length(maxtree.parentindices)
Base.size(maxtree::MaxTree) = size(maxtree.parentindices)
Base.:(==)(a::MaxTree, b::MaxTree) =
    (ndims(a) == ndims(b)) && (a.rev == b.rev) &&
    (a.parentindices == b.parentindices) &&
    (a.traverse == b.traverse)

"""
    root_index(maxtree::MaxTree) -> Int

Linear index of the root pixel.
"""
root_index(maxtree::MaxTree) = maxtree.traverse[1]

# Gets the path to the root of the max-tree (defined by `ancestors`) starting from
# `node`. The indices of the nodes along the path are stored into `path`.
# Returns the index of the root
@inline function rootpath!(path::Vector{Int},
                           ancestors::AbstractArray{<:Integer},
                           node::Integer)
    empty!(path)
    p = node
    @label next
    push!(path, p)
    @inbounds q = ancestors[p]
    if p != q
        p = q
        @goto next
    end
    return p
end

# Corrects `maxtree.parentindices`, so that the parent of every pixel
# is a canonical node. I.e. for each connected component at a level l,
# all pixels point to the same representative, which in turn points to the
# representative pixel at the next level.
#
# A canonical node ``p`` is either a root of the max-tree, or
# ``image(parent(p)) != image(p)``
function canonize!(maxtree::MaxTree, image::GenericGrayImage)
    @assert size(maxtree) == size(image)
    parentindices = maxtree.parentindices
    @inbounds for p in maxtree.traverse
        q = parentindices[p]
        if image[q] == image[parentindices[q]]
            parentindices[p] = parentindices[q]
        end
    end
    return maxtree
end

# Checks whether the `pixel` + `offset` is still a valid pixel
# contained inside the image area defined by `axes`.
isvalid_offset(offset::CartesianIndex{N}, pixel::CartesianIndex{N},
               axes::NTuple{N}) where N =
    Base.checkbounds_indices(Bool, axes, (pixel + offset,))

# generates a vector of offsets to the CartesianIndex that defines
# the neighborhood of the pixel in an N-dimensional image.
# `connectivity` is the maximal number of orthogonal steps
# (+1/-1 to a single dimension) that are required to reach a neighbor.
# Note that the corresponding offsets of linear indices would be sorted in ascending order
# ensuring efficient memory access pattern.
function neighbor_cartesian_offsets(::Type{CartesianIndex{N}}, connectivity::Integer) where N
    (connectivity >= 1) || throw(ArgumentError("connectivity should be positive integer"))
    return [ci for ci in CartesianIndices(ntuple(_ -> -1:1, N))
            if 1 <= sum(abs, Tuple(ci)) <= connectivity]
end

# convert offsets of cartesian indices into linear offsets of the given array
linear_offsets(offsets::AbstractVector{<:CartesianIndex},
               arr::AbstractArray) =
    dot.(Tuple.(offsets), Ref(strides(arr)))

"""
    rebuild!(maxtree::MaxTree, image::GenericGrayImage,
             neighbors::AbstractVector{CartesianIndex}) -> maxtree

Rebuilds the `maxtree` for the `image` using `neighbors` as the pixel
connectivity specification.

# Details
The pixels in the connected components generated by the method should be
connected to each other by a path through neighboring pixels.
The pixels ``p_1`` and ``p_2`` are neighbors, if `neighbors` array contains
``d``, such that ``p_2 = p_1 + d``.

# See also
[`MaxTree`](@ref)
"""
function rebuild!(maxtree::MaxTree{N},
                  image::GenericGrayImage{<:Any, N},
                  #=mask::AbstractArray{Bool, N},=#
                  neighbors::AbstractVector{CartesianIndex{N}}) where N
    check_maxtree(maxtree, image, rev=maxtree.rev)
    # pixels need to be sorted according to their gray level.
    sortperm!(maxtree.traverse, vec(image), rev=maxtree.rev)

    neighbor_offsets = linear_offsets(neighbors, maxtree.parentindices)
    p_indices = fill!(maxtree.parentindices, 0) |> vec # 0=uninitialized parent
    p_cis = CartesianIndices(maxtree.parentindices)
    p_axes = axes(maxtree.parentindices)
    r_indices = fill!(similar(p_indices), 0) # temporary array containing current roots
    rootpath = Vector{Int}() # temp array to hold the path to the current root within ancestors

    # traverse the array in reversed order (from highest value to lowest value)
    @inbounds for p in Iterators.Reverse(maxtree.traverse)
        p_ci = p_cis[p]
        p_indices[p] = r_indices[p] = p # it's a new root

        for i in eachindex(neighbors)
            # check that p is in the mask and that it's neighbor is a valid image pixel
            (#=mask[p] && =#isvalid_offset(neighbors[i], p_ci, p_axes)) || continue
            index = p + neighbor_offsets[i] # linear index of the neighbor
            (p_indices[index] == 0) && continue # ignore neighbor without parent (this means image[index] < image[p])
            root = rootpath!(rootpath, r_indices, index) # neighbor's root
            if (root != p) || (length(rootpath) > 1) # attach neighbor's root to the p
                p_indices[root] = p
                r_indices[rootpath] .= p # also attach the ancestral nodes of neighbor to p
            end
        end
    end
    return canonize!(maxtree, image)
end

"""
    MaxTree(image::GenericGrayImage; connectivity=1, rev=false) -> MaxTree

Constructs the *max-tree* of the `image`.

# Arguments
- `connectivity::Integer=1`: defines the pixel neighborhood used to construct
  the connected components. The value is the maximum number of orthogonal steps
  to reach a neighbor. In 2D, it is 1 for a 4-neighborhood and 2 for a
  8-neighborhood. See [`rebuild!](@ref).
- `rev::Bool=false`: if `false`, the max-tree is traversed from the darkest
  (the root node) to the brightest, otherwise it's traversed from the brightest
  (the root) to the darkest.

# Examples
We create a small sample image (Figure 1 from [4]) and build the max-tree.

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Int64,2}|Matrix{Int64}"
julia> image = [15 13 16; 12 12 10; 16 12 14]
3×3 Array{Int64,2}:
 15  13  16
 12  12  10
 16  12  14

julia> mtree = MaxTree(image, connectivity=2)
MaxTree{2}(false, [4 2 4; 8 2 8; 2 2 2], [8, 2, 5, 6, 4, 9, 1, 3, 7])
```
"""
function MaxTree(image::GenericGrayImage;
                 connectivity::Integer=1, rev::Bool=false)
    N = ndims(image)
    # User defined masks are not allowed, as there might be more than one
    # connected component in the mask (and therefore not a single tree that
    # represents the image). Mask here is an image that is 0 on the border
    # and 1 everywhere else
    #=
    mask = np.ones(image.shape)
    for k in range(len(image.shape)):
        np.moveaxis(mask, k, 0)[0] = 0
        np.moveaxis(mask, k, 0)[-1] = 0
    =#
    neighbors = neighbor_cartesian_offsets(CartesianIndex{N}, connectivity)
    return rebuild!(MaxTree{N}(image, rev), #=mask,=# image, neighbors)
end

"""
    areas(maxtree::MaxTree) -> Array{Int}

Computes the areas of all `maxtree` components.

# Returns
The array of the same shape as the original image. The `i`-th element is
the area (in pixels) of the component that is represented by the reference
pixel with index `i`.

# See also
[`diameters`](@ref), [`area_opening`](@ref), [`area_closing`](@ref).
"""
function areas(maxtree::MaxTree)
    areas = fill(1, size(maxtree)) # start with 1-pixel areas
    parentindices = maxtree.parentindices
    @inbounds for p in Iterators.Reverse(maxtree.traverse)
        q = parentindices[p]
        (q != p) && (areas[q] += areas[p])
    end
    return areas
end

"""
    boundingboxes(maxtree::MaxTree) -> Array{NTuple{2, CartesianIndex}}

Computes the minimal bounding boxes of all `maxtree` components.

# Returns
The array of the same shape as the original image. The `i`-th element is
the tuple of the minimal and maximal cartesian indices for the bounding box
of the component that is represented by the reference pixel with index `i`.

# See also
[`diameters`](@ref).
"""
function boundingboxes(maxtree::MaxTree{N}) where N
    # initialize bboxes
    bboxes = [(ci, ci) for ci in CartesianIndices(maxtree.parentindices)]
    @inbounds for p in Iterators.Reverse(maxtree.traverse)
        q = maxtree.parentindices[p]
        bbp, bbq = bboxes[p], bboxes[q]
        bboxes[q] = (min(bbq[1], bbp[1]), max(bbq[2], bbp[2]))
    end
    return bboxes
end

"""
    diameters(maxtree::MaxTree) -> Array{Int}

Computes the "diameters" of all `maxtree` components.

"Diameter" of the *max-tree* connected component is the length of
the widest side of the component's bounding box.

# Returns
The array of the same shape as the original image. The `i`-th element is the
"diameter" of the component that is represented by the reference pixel with
index `i`.

# See also
[`boundingboxes`](@ref), [`areas`](@ref),
[`diameter_opening`](@ref), [`diameter_closing`](@ref).
"""
diameters(maxtree::MaxTree) =
    [maximum(Tuple(bbox[2] - bbox[1])) + 1 for bbox in boundingboxes(maxtree)]

"""
    filter_components!(output::GenericGrayImage, image::GenericGrayImage,
                       maxtree::MaxTree, attrs::AbstractVector,
                       min_attr, all_below_min) -> output

Filters the connected components of the `image` and stores the result in `output`.

The ``output`` is the copy of the ``image`` exluding the connected components,
whose attribute value is below `min_attr`. That is, the pixels of the exluded
component are reset to the value of the reference pixel of its first valid
ancestor (the connected component with the attribute value greater or equal
to `min_attr`).

# Arguments
- `maxtree::MaxTree`: pre-built max-tree representation of the `image`
- `attrs::AbstractVector`: `attrs[i]` is the attribute value for the ``i``-th
   component of the tree (``i`` being the linear index of its *reference pixel*)
- `all_below_min`: the value to fill the `output` if all attributes of all
  components (including the root one) are below `min_attr`

# Details
This function is the basis for [`area_opening`](@ref), [`diameter_opening`](@ref)
and similar transformations.
E.g. for [`area_opening`](@ref) the attribute is the area of the components.
In this case, the max-tree components of the `output` have area no smaller
than `min_attr` pixels.

The method assumes that the attribute values are monotone with respect to the
components hieararchy, i.e. ``attrs[i] <= attrs[maxtree.parentindices[i]]`` for
each `i`.
"""
function filter_components!(output::GenericGrayImage, image::GenericGrayImage,
                            maxtree::MaxTree,
                            attrs::AbstractArray, min_attr,
                            all_below_min = zero(eltype(output)))
    # should have been already checked by higher-level functions
    @assert size(output) == size(image)
    @assert size(image) == size(maxtree)
    @assert length(attrs) == length(maxtree)

    p_root = root_index(maxtree)
    # if p_root is below min_attr, then all others are as well
    (attrs[p_root] < min_attr) && return fill!(output, all_below_min)
    output[p_root] = image[p_root]
    parentindices = maxtree.parentindices

    for p in maxtree.traverse
        q = parentindices[p]
        # this means p is not canonical
        # in other words, it has a parent that has the
        # same image value.
        if image[p] == image[q]
            output[p] = output[q]
            continue
        end
        # if attribute is below the threshold, stop at the lower parent,
        # otherwise continue to reconstruct the image
        output[p] = attrs[p] < min_attr ? output[q] : image[p]
    end

    return output
end

function check_output_image(output::GenericGrayImage, image::GenericGrayImage)
    (size(output) == size(image)) ||
        throw(DimensionMismatch("The sizes of the output and the input image do not match"))
end

# checks if the provided maxtree is compatible with the given options
# or builds the new maxtree if none was given
function check_maxtree(maxtree::MaxTree,
                       image::GenericGrayImage;
                       rev::Bool=false)
    (size(maxtree) == size(image)) ||
        throw(DimensionMismatch("The sizes of the max-tree and the input image do not match"))
    (maxtree.rev == rev) ||
        throw(ArgumentError("The traversal order of the given max-tree is different from the requested one"))
end

"""
    area_opening!(output, image, [maxtree];
                  min_area=64, connectivity=1) -> output

Performs in-place *area opening* of the `image` and stores the result in `output`.
See [`area_opening`](@ref) for the detailed description of the method.
"""
function area_opening!(output::GenericGrayImage, image::GenericGrayImage,
                       maxtree::MaxTree; min_area::Number=64)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=false)
    return filter_components!(output, image, maxtree, areas(maxtree), min_area)
end

"""
    area_opening(image, [maxtree]; min_area=64, connectivity=1) -> Array

Performs an *area opening* of the `image`.

*Area opening* replaces all bright components of an image that
have a surface smaller than `min_area` with the darker value taken from
their first ancestral component (in *max-tree* representation of `image`)
that has the area no smaller than `min_area`.

# Details
Area opening is similar to morphological opening (see [`opening`](@ref)),
but instead of using a fixed structuring element (e.g. disk) it employs
small (less than `min_area`) components of the *max-tree*. Consequently,
the `area_opening` with `min_area = 1` is the identity transformation.

In the binary case, area opening is equivalent to `remove_small_objects`;
this operator is thus extended to gray-level images.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `min_area::Number=64`: the smallest size (in pixels) of the image component
  to keep intact
- `connectivity::Integer=1`: the neighborhood connectivity. The maximum number
  of orthogonal steps to reach a neighbor of the pixel.
  In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An array of the same type and shape as the `image`.

# See also
[`area_opening!`](@ref), [`area_closing`](@ref),
[`diameter_opening`](@ref),
[`MaxTree`](@ref), [`opening`](@ref)

# References
1. Vincent, L. (1993). *Grayscale area openings and closings, their efficient implementation and applications*, Proc. of EURASIP Workshop on Mathematical Morphology and its Applications to Signal Processing, Barcelona, Spain, 22-27
2. Soille, P. (2003). Chapter 6 *Geodesic Metrics* of *Morphological Image Analysis: Principles and Applications*, 2nd edition, Springer.
   > https://doi.org/10.1007/978-3-662-05088-0
3. Salembier, P., Oliveras, A., & Garrido, L. (1998). *Antiextensive Connected Operators for Image and Sequence Processing*. IEEE Transactions on Image Processing, 7(4), 555-570.
   > https://doi.org/10.1109/83.663500
4. Najman, L., & Couprie, M. (2006). *Building the component tree in quasi-linear time*. IEEE Transactions on Image Processing, 15(11), 3531-3539.
   > https://doi.org/10.1109/TIP.2006.877518
5. Carlinet, E., & Geraud, T. (2014). *A Comparative Review of Component Tree Computation Algorithms*. IEEE Transactions on Image Processing, 23(9), 3885-3895.
   > https://doi.org/10.1109/TIP.2014.2336551

# Examples
Creating a test image `f` (quadratic function with a maximum in the center and
4 additional local maxima):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 12;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_aopen = area_opening(f, min_area=8, connectivity=1);
```
The peaks with a surface smaller than 8 are removed.
"""
area_opening

"""
    diameter_opening!(output, image, [maxtree];
                      min_diameter=8, connectivity=1) -> output

Performs in-place *diameter opening* of the `image` and stores the result in `output`.
See [`diameter_opening`](@ref) for the detailed description of the method.
"""
function diameter_opening!(output::GenericGrayImage, image::GenericGrayImage,
                           maxtree::MaxTree; min_diameter=8)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=false)
    return filter_components!(output, image, maxtree, diameters(maxtree), min_diameter)
end

"""
    diameter_opening(image, [maxtree]; min_diameter=8, connectivity=1) -> Array

Performs a *diameter opening* of the `image`.

*Diameter opening* replaces all bright structures of an image that have
the diameter (the widest dimension of their bounding box) smaller than
`min_diameter` with the darker value taken from their first ancestral component
(in *max-tree* representation of `image`) that has the diameter no smaller than
`min_diameter`.

The operator is also called *Bounding Box Opening*. In practice,
the result is similar to a *morphological opening*, but long and thin
structures are not removed.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `min_diameter::Number=8`: the minimal length (in pixels) of the widest
  dimension of the bounding box of the image component to keep intact
- `connectivity::Integer=1`: the neighborhood connectivity. The maximum number
  of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for
  a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An array of the same type and shape as the `image`.

# See also
[`diameter_opening!`](@ref), [`diameter_closing`](@ref),
[`area_opening`](@ref),
[`MaxTree`](@ref), [`opening`](@ref)

# References
1. Walter, T., & Klein, J.-C. (2002). *Automatic Detection of Microaneurysms in Color Fundus Images of the Human Retina by Means of the Bounding Box Closing*. In A. Colosimo, P. Sirabella, A. Giuliani (Eds.), *Medical Data Analysis. Lecture Notes in Computer Science*, vol 2526, 210-220. Springer Berlin Heidelberg.
   > https://doi.org/10.1007/3-540-36104-9_23
2. Carlinet, E., & Geraud, T. (2014). *A Comparative Review of Component Tree Computation Algorithms*. IEEE Transactions on Image Processing, 23(9), 3885-3895.
   > https://doi.org/10.1109/TIP.2014.2336551

# Examples
Creating a test image `f` (quadratic function with a maximum in the center and
4 additional local maxima):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 12;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_dopen = diameter_opening(f, min_diameter=3, connectivity=1);
```
The peaks with a maximal diameter of 2 or less are removed.
For the remaining peaks the widest side of the bounding box is at least 3.
"""
diameter_opening

"""
    area_closing!(output, image, [maxtree];
                  min_area=64, connectivity=1) -> output

Performs in-place *area closing* of the `image` and stores the result in `output`.
See [`area_closing`](@ref) for the detailed description of the method.
"""
function area_closing!(output::GenericGrayImage, image::GenericGrayImage,
                       maxtree::MaxTree; min_area::Number=64)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=true)
    return filter_components!(output, image, maxtree, areas(maxtree), min_area)
end

"""
    area_closing(image, [maxtree]; min_area=64, connectivity=1) -> Array

Performs an *area closing* of the `image`.

*Area closing* replaces all dark components of an image that
have a surface smaller than `min_area` with the brighter value taken from
their first ancestral component (in *max-tree* representation of `image`)
that has the area no smaller than `min_area`.

# Details
*Area closing* is the dual operation to *area opening* (see [`area_opening`](@ref)).
It is similar to morphological closings (see [`closing`](@ref)),
but instead of using a fixed structuring element (e.g. disk) it employs
small (less than `min_area`) components of the *max-tree*. Consequently,
the `area_closing` with `min_area = 1` is the identity transformation.

In the binary case, area closing is equivalent to
`remove_small_holes`; this operator is thus extended to gray-level images.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `min_area::Number=64`: the smallest size (in pixels) of the image component
  to keep intact
- `connectivity::Integer=1`: the neighborhood connectivity. The maximum number
  of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for
  a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An array of the same type and shape as the `image`.

# See also
[`area_closing!`](@ref), [`area_opening`](@ref),
[`diameter_closing`](@ref),
[`MaxTree`](@ref), [`closing`](@ref)

# References
1. Vincent, L. (1993). *Grayscale area openings and closings, their efficient implementation and applications*, Proc. of EURASIP Workshop on Mathematical Morphology and its Applications to Signal Processing, Barcelona, Spain, 22-27
2. Soille, P. (2003). Chapter 6 *Geodesic Metrics* of *Morphological Image Analysis: Principles and Applications*, 2nd edition, Springer.
   > https://doi.org/10.1007/978-3-662-05088-0
3. Salembier, P., Oliveras, A., & Garrido, L. (1998). *Antiextensive Connected Operators for Image and Sequence Processing*. IEEE Transactions on Image Processing, 7(4), 555-570.
   > https://doi.org/10.1109/83.663500
4. Najman, L., & Couprie, M. (2006). *Building the component tree in quasi-linear time*. IEEE Transactions on Image Processing, 15(11), 3531-3539.
   > https://doi.org/10.1109/TIP.2006.877518
5. Carlinet, E., & Geraud, T. (2014). *A Comparative Review of Component Tree Computation Algorithms*. IEEE Transactions on Image Processing, 23(9), 3885-3895.
   > https://doi.org/10.1109/TIP.2014.2336551

# Examples
Creating a test image `f` (quadratic function with a minimum in the center and
4 additional local minima):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 12;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_aclose = area_closing(f, min_area=8, connectivity=1);
```
All small minima are removed, and the remaining minima have at least
a size of 8.
"""
area_closing

"""
    diameter_closing!(output, image, [maxtree];
                      min_diameter=8, connectivity=1) -> output

Performs in-place *diameter closing* of the `image` and stores the result in `output`.
See [`diameter_closing`](@ref) for the detailed description of the method.
"""
function diameter_closing!(output::GenericGrayImage, image::GenericGrayImage,
                           maxtree::MaxTree; min_diameter::Number=8)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=true)
    return filter_components!(output, image, maxtree, diameters(maxtree), min_diameter)
end

"""
    diameter_closing(image, [maxtree]; min_diameter=8, connectivity=1) -> Array

Performs a *diameter closing* of the `image`.

*Diameter closing* replaces all dark structures of an image that have
the diameter (the widest dimension of their bounding box) smaller than
`min_diameter` with the brighter value taken from their first ancestral component
(in *max-tree* representation of `image`) that has the diameter no smaller than
`min_diameter`.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `min_diameter::Number=8`: the minimal length (in pixels) of the widest
  dimension of the bounding box of the image component to keep intact
- `connectivity::Integer=1`: the neighborhood connectivity. The maximum number
  of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for
  a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An array of the same type and shape as the `image`.

# See also
[`diameter_closing!`](@ref), [`diameter_opening`](@ref),
[`area_closing`](@ref),
[`MaxTree`](@ref), [`closing`](@ref)

# References
1. Walter, T., & Klein, J.-C. (2002). *Automatic Detection of Microaneurysms in Color Fundus Images of the Human Retina by Means of the Bounding Box Closing*. In A. Colosimo, P. Sirabella, A. Giuliani (Eds.), *Medical Data Analysis. Lecture Notes in Computer Science*, vol 2526, 210-220. Springer Berlin Heidelberg.
   > https://doi.org/10.1007/3-540-36104-9_23
2. Carlinet, E., & Geraud, T. (2014). *A Comparative Review of Component Tree Computation Algorithms*. IEEE Transactions on Image Processing, 23(9), 3885-3895.
   > https://doi.org/10.1109/TIP.2014.2336551

# Examples
Creating a test image `f` (quadratic function with a minimum in the center and
4 additional local minima):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 12;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_dclose = diameter_closing(f, min_diameter=3, connectivity=1);
```
All small minima with a diameter of 2 or less are removed.
For the remaining minima the widest bounding box side is at least 3.
"""
diameter_closing

# Calculates the local maxima or minima of the image using the max-tree
# representation (maxima if maxtree.rev=false, minima otherwise).
# It's not the most effecient method for extrema calculation, but could be
# useful if the max-tree is already calculated.
# Each minima/maxima region is labeled with the unique id (1 is the global
# maximum/minimum).
function local_extrema!(output::GenericGrayImage,
                        image::GenericGrayImage,
                        maxtree::MaxTree)
    (size(image) == size(maxtree)) || throw(DimensionMismatch())
    (size(output) == size(image)) || throw(DimensionMismatch())

    next_label = oneunit(eltype(output))
    fill!(output, next_label) # initialize output (just has to be non-zero)
    @inbounds for p in Iterators.reverse(maxtree.traverse)
        q = maxtree.parentindices[p]
        # if p is canonical (parent has a different value)
        if image[p] != image[q]
            output[q] = 0
            # if output[p] was the parent of some other canonical
            # pixel, it has been set to zero. Only the leaves
            # (local maxima) are thus > 0.
            if output[p] > 0
                output[p] = next_label
                next_label += 1
            end
        end
    end

    @inbounds for p in Iterators.reverse(maxtree.traverse)
        q = maxtree.parentindices[p]
        # if p is not canonical (parent has the same value)
        if image[p] == image[q]
            # in this case we propagate the value
            output[p] = output[q]
        end
    end

    return output
end

"""
    local_maxima!(output, image, [maxtree]; connectivity=1) -> output

Detects the local maxima of `image` and stores the result in `output`.
See [`local_maxima`](@ref) for the detailed description of the method.
"""
function local_maxima!(output::GenericGrayImage, image::GenericGrayImage,
                       maxtree::MaxTree)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=false)
    return local_extrema!(output, image, maxtree)
end

"""
    local_maxima(image, [maxtree]; connectivity=1) -> Array

Determines and labels all *local maxima* of the `image`.

# Details
The *local maximum* is defined as the connected set of pixels that have the
same value, which is greater than the values of all pixels in direct
neighborhood of the set.

Technically, the implementation is based on the *max-tree* representation
of an image. It's beneficial if the max-tree is already computed,
otherwise [`Images.findlocalmaxima`](@ref) would be more efficient.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `connectivity::Integer=1`: the neighborhood connectivity.
  The maximum number of orthogonal steps to reach a neighbor of the pixel.
  In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An integer array of the same shape as the `image`. Pixels that are not local
maxima have 0 value. Pixels of the same local maximum share the same positive
value (the local maximum id).

# See also
[`MaxTree`](@ref), [`local_maxima!`](@ref), [`local_minima`](@ref),
[`Images.findlocalmaxima`](@ref)

# Examples
Create `f` (quadratic function with a maximum in the center and
4 additional constant maxima):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 10;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:5, 3:5] .= 40; f[3:5, 8:10] .= 60; f[8:10, 3:5] .= 80; f[8:10, 8:10] .= 100;

julia> f_maxima = local_maxima(f); # Get all local maxima of `f`
```
The resulting image contains the 4 labeled local maxima.
"""
local_maxima

"""
    local_minima!(output, image, [maxtree]; connectivity=1) -> output

Detects the local minima of `image` and stores the result in `output`.
See [`local_minima`](@ref) for the detailed description of the method.
"""
function local_minima!(output::GenericGrayImage, image::GenericGrayImage,
                       maxtree::MaxTree; connectivity::Integer=1)
    check_output_image(output, image)
    check_maxtree(maxtree, image, rev=true)
    return local_extrema!(output, image, maxtree)
end

"""
    local_minima(image, [maxtree]; connectivity=1) -> Array

Determines and labels all *local minima* of the `image`.

# Details
The *local minimum* is defined as the connected set of pixels that have the
same value, which is less than the values of all pixels in direct
neighborhood of the set.

Technically, the implementation is based on the *max-tree* representation
of an image. It's beneficial if the max-tree is already computed,
otherwise [`Images.findlocalminima`](@ref) would be more efficient.

# Arguments
- `image::GenericGrayImage`: the ``N``-dimensional input image
- `connectivity::Integer=1`: the neighborhood connectivity. The maximum number
  of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for
  a 4-neighborhood and 2 for a 8-neighborhood.
- `maxtree::MaxTree`: optional pre-built *max-tree*. Note that
  `maxtree` and `connectivity` optional parameters are mutually exclusive.

# Returns
An integer array of the same shape as the `image`. Pixels that are not local
minima have 0 value. Pixels of the same local minimum share the same positive
value (the local minimum id).

# See also
[`MaxTree`](@ref), [`local_minima!`](@ref), [`local_maxima`](@ref),
[`Images.findlocalminima`](@ref)

# Examples
Create `f` (quadratic function with a minimum in the center and
4 additional constant minimum):
```jldoctest; setup = :(using ImageMorphology)
julia> w = 10;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:5, 3:5] .= 40; f[3:5, 8:10] .= 60; f[8:10, 3:5] .= 80; f[8:10, 8:10] .= 100;

julia> f_minima = local_minima(f); # Calculate all local minima of `f`
```
The resulting image contains the labeled local minima.
"""
local_minima

# autogenerate in-place maxtree operations with "connectivity" keyarg
for (fun, rev) in [(:area_opening!, false), (:area_closing!, true),
                   (:diameter_opening!, false), (:diameter_closing!, true),
                   (:local_maxima!, false), (:local_minima!, true)]
    @eval begin
        $fun(output::GenericGrayImage, image::GenericGrayImage;
             connectivity::Integer=1, kwargs...) =
            $fun(output, image,
                 MaxTree(image, connectivity=connectivity, rev=$rev); kwargs...)
    end
end

# autogenerate maxtree operations (simple wrappers of in-place methods)
for fun in [:area_opening, :area_closing,
            :diameter_opening, :diameter_closing,
            :local_maxima, :local_minima]
    fun! = Symbol(fun, "!")
    @eval begin
        $fun(image::GenericGrayImage; kwargs...) =
             $fun!(similar(image), image; kwargs...)
        $fun(image::GenericGrayImage, maxtree::MaxTree; kwargs...) =
             $fun!(similar(image), image, maxtree; kwargs...)
    end
end


text:
using Base: depwarn

Base.@deprecate GLRM(A::AbstractArray, obs::Array{Tuple{Int, Int}, 1}, args...; kwargs...) GLRM(A, args...; obs = obs, kwargs...)

Base.@deprecate ProxGradParams(s::Number,m::Int,c::Float64,ms::Float64) ProxGradParams(s, max_iter=m, abs_tol=c, min_stepsize=ms)

Base.@deprecate expand_categoricals expand_categoricals!

Base.@deprecate errors(g::GLRM) error_metric(g)

Base.@deprecate quadratic QuadLoss

Base.@deprecate logistic LogisticLoss

Base.@deprecate huber HuberLoss

Base.@deprecate LogLoss LogisticLoss

Base.@deprecate l1 L1Loss

Base.@deprecate poisson PoissonLoss

Base.@deprecate ordinal_hinge OrdinalHingeLoss

Base.@deprecate OrdinalHinge OrdinalHingeLoss

Base.@deprecate WeightedHinge WeightedHingeLoss

Base.@deprecate periodic PeriodicLoss

Base.@deprecate quadreg QuadReg

Base.@deprecate constrained_quadreg QuadConstraint

Base.@deprecate onereg OneReg

Base.@deprecate zeroreg ZeroReg

Base.@deprecate nonnegative NonNegConstraint

Base.@deprecate onesparse OneSparseConstraint

Base.@deprecate unitonesparse UnitOneSparseConstraint

Base.@deprecate simplex SimplexConstraint

Base.@deprecate nonneg_onereg NonNegOneReg


text:
@testset "557.reverse-words-in-a-string-iii.jl" begin
    @test reverse_words_557("Let's take LeetCode contest") == "s'teL ekat edoCteeL tsetnoc"
    @test reverse_words_557("God Ding") == "doG gniD"
end


text:
struct MPISum{T}
  data::T
  comm::MPI.Comm
end

MPISum(data::T, comm=MPI.COMM_WORLD) where {T} = MPISum{T}(data, comm)

# Apply the sum
function (A::MPISum)(v)
  return MPI.Allreduce(A.data(v), +, A.comm)
end

max_stars_repo_name: mtfishman/ITensorParallel


text:
# ---
# title: 475. Heaters
# id: problem475
# author: Qling
# date: 2021-04-07
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/heaters/description/>
# hidden: true
# ---
#
# Winter is coming! During the contest, your first job is to design a standard
# heater with a fixed warm radius to warm all the houses.
#
# Every house can be warmed, as long as the house is within the heater's warm
# radius range.
#
# Given the positions of `houses` and `heaters` on a horizontal line, return
# _the minimum radius standard of heaters  so that those heaters could cover all
# houses._
#
# **Notice** that all the `heaters` follow your radius standard, and the warm
# radius will the same.
#
#
#
# **Example 1:**
#
#
#
#     Input: houses = [1,2,3], heaters = [2]
#     Output: 1
#     Explanation: The only heater was placed in the position 2, and if we use the radius 1 standard, then all the houses can be warmed.
#
#
# **Example 2:**
#
#
#
#     Input: houses = [1,2,3,4], heaters = [1,4]
#     Output: 1
#     Explanation: The two heater was placed in the position 1 and 4. We need to use radius 1 standard, then all the houses can be warmed.
#
#
# **Example 3:**
#
#
#
#     Input: houses = [1,5], heaters = [2]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= houses.length, heaters.length <= 3 * 104`
#   * `1 <= houses[i], heaters[i] <= 109`
#
#
## @lc code=start
using LeetCode

function find_radius(houses::Vector{Int}, heaters::Vector{Int})::Int
    sort!(heaters)
    heaters = [-Inf; heaters; Inf]
    res = 0
    for house in houses
        loc = searchsortedfirst(heaters, house)
        res = max(res, min(house - heaters[loc - 1], heaters[loc] - house))
    end

    return convert(Int, res)
end

## @lc code=end


text:
# Data container
immutable BiCGStabData{T}
    v::Vector{T}
    t::Vector{T}
    p::Vector{T}
    phat::Vector{T}
    s::Vector{T}
    shat::Vector{T}
    r::Vector{T}
    rtld::Vector{T}
end

BiCGStabData(n::Int, T::Type) = BiCGStabData(zeros(T, n), zeros(T, n),
                                             zeros(T, n), zeros(T, n),
                                             zeros(T, n), zeros(T, n),
                                             zeros(T, n), zeros(T, n))

# Main solver
function bicgstab!{T<:Real}(A, b::Vector{T}, x::Vector{T}; tol::Float64=1e-6, maxIter::Int64=100,
                                                           tolRho::Float64=1e-40, precon=copy!,
                                                           data=BiCGStabData(length(b), T))
    n = length(b)
    out_int = 0
    n_iter = 0
    converge = false
    bnrm2 = genblas_nrm2(b)
    if bnrm2 == 0.0
        x .= 0.0
        return 1, 0
    end
    A(data.r, x)
    genblas_scal!(-one(T), data.r)
    genblas_axpy!(one(T), b, data.r)
    residual = genblas_nrm2(data.r)/bnrm2
    if residual == 0.0
        return 2, 0
    end
    alpha = 1.0
    omega = 1.0
    data.rtld .= data.r
    rho1 = 0.0
    data.p .= data.r
    for iter = 1:maxIter
        rho = genblas_dot(data.rtld, data.r)
        if abs(rho) <= tolRho
            return -11, iter
        end
        if iter > 1
            beta  = (rho/rho1)*(alpha/omega)
            # p = r + beta*(p - omega*v)
            genblas_scal!(beta, data.p)
            genblas_axpy!(-(beta*omega), data.v, data.p)
            genblas_axpy!(one(T), data.r, data.p)
        end
        precon(data.phat, data.p)
        A(data.v, data.phat)
        alpha = rho/genblas_dot(data.rtld, data.v)
        # s = r - alpha*v
        data.s .= data.r
        genblas_axpy!(-alpha, data.v, data.s)
        # x = x + alpha*phat
        genblas_axpy!(alpha, data.phat, x)
        residual = genblas_nrm2(data.s)/bnrm2
        if residual <= tol
            return 31, iter
        end
        precon(data.shat, data.s)
        A(data.t, data.shat)
        omega = genblas_dot(data.t, data.s)/genblas_dot(data.t, data.t)
        # x = x + omega*shat
        genblas_axpy!(omega, data.shat, x)
        # r = s .- omega.*t
        data.r .= data.s
        genblas_axpy!(-omega, data.t, data.r)
        residual = genblas_nrm2(data.r)/bnrm2
        if residual <= tol
            return 32, iter
        end
        if abs(omega) < 1e-16
            return -12, iter
        end
        rho1 = rho
    end
    return -2, maxIter
end

# API

function bicgstab{T<:Real}(A, b::Vector{T}; tol::Float64=1e-6, maxIter::Int64=100,
                                            tolRho::Float64=1e-40, precon=copy!,
                                            data=BiCGStabData(length(b), T))
    x = zeros(b)
    exit_code, num_iters = bicgstab!(A, b, x, tol=tol, maxIter=maxIter, tolRho=tolRho, precon=precon, data=data)
    return x, exit_code, num_iters
end

export BiCGStabData, bicgstab!, bicgstab


text:
# Please run from the `basic_gppp` directory.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using AbstractGPs, Plots, Random, Stheno

###########################  Define our model  ###########################

# Define a distribution over f₁, f₂, and f₃, where f₃(x) = f₁(x) + f₂(x).
f = @gppp let
    f1 = GP(randn(), SEKernel())
    f2 = GP(SEKernel())
    f3 = f1 + f2
end;

# Randomly sample `N₁` and `N₃` locations at which to observe `f₁` and `f₃` respectively.
rng, N1, N3 = MersenneTwister(123546), 10, 11;
x1 = GPPPInput(:f1, sort(rand(rng, N1) * 10));
x3 = GPPPInput(:f3, sort(rand(rng, N3) * 10));
x = BlockData(x1, x3);

# Generate some toy of `f1` and `f3`, `y1` and `y3` respectively.
fx = f(x);
y = rand(rng, f(x));
y1, y3 = split(x, y);

# Compute the posterior processes.
f_post = posterior(fx, y);

# Define some plotting stuff.
Np, S = 500, 25;
xp_ = range(-2.5, stop=12.5, length=Np);
xp = BlockData(GPPPInput(:f1, xp_), GPPPInput(:f2, xp_), GPPPInput(:f3, xp_));

# Sample jointly from the posterior over each process.
f_samples = rand(rng, f_post(xp, 1e-9), S);
f′1_xp, f′2_xp, f′3_xp = split(xp, f_samples);

# Compute posterior marginals.
ms = marginals(f_post(xp, 1e-9));
f′1_m, f′2_m, f′3_m = split(xp, mean.(ms));
f′1_s, f′2_s, f′3_s = split(xp, std.(ms));

###########################  Plot results  ###########################

gr();
posterior_plot = plot();

# Plot posterior over f1.
plot!(posterior_plot, xp_, f′1_m; ribbon=3f′1_s, color=:red, label="f1", fillalpha=0.3);
plot!(posterior_plot, xp_, f′1_xp; color=:red, label="", alpha=0.2, linewidth=1);
scatter!(posterior_plot, x1.x, y1;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

# Plot posterior over f2.
plot!(posterior_plot, xp_, f′2_m; ribbon=3f′2_s, color=:green, label="f2", fillalpha=0.3);
plot!(posterior_plot, xp_, f′2_xp; color=:green, label="", alpha=0.2, linewidth=1);

# Plot posterior over f3
plot!(posterior_plot, xp_, f′3_m; ribbon=3f′3_s, color=:blue, label="f3", fillalpha=0.3);
plot!(posterior_plot, xp_, f′3_xp; color=:blue, label="", alpha=0.2, linewidth=1);
scatter!(posterior_plot, x3.x, y3;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

display(posterior_plot);


text:
println("输入一个整数")
x=parse(Int,readline())
if x>0
    print(">0")
elseif x==0
    print("=0")
else
    print("<0")
end

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
try
	using RCall, CSV, DataFrames, GZip
catch
	using Pkg
	Pkg.add.(["RCall", "CSV", "DataFrames", "GZip"]);
end

"""
	Estimating and plotting MAP using locfit and ggplot2 in R. It assume your folder contains the posterior estimated through ABCreg
"""
function plotMap(;analysisFolder::String,weak::Bool=true,title::String="Posteriors")

	try
		@eval using RCall
		@eval R"""library(ggplot2);library(locfit);library(data.table)"""

		out = filter(x -> occursin("post",x), readdir(analysisFolder,join=true))
		out          = filter(x -> !occursin(".1.",x),out)

		open(x)      = Array(CSV.read(x,DataFrame,header=false))

		# Control outlier inference. 2Nes non negative values
		flt(x)       = x[x[:,4] .> 0,:]
		posteriors   = flt.(open.(out))

		R"""getmap <- function(df){
				temp = as.data.frame(df)
				d <-locfit(~temp[,1],temp);
				map<-temp[,1][which.max(predict(d,newdata=temp))]
			}"""
		#getmap(x)    = rcopy(R"""matrix(apply($x,2,getmap),nrow=1)""")
		getmap(x)    = rcopy(R"""suppressWarnings(matrix(apply($x,2,getmap),nrow=1))""")

		if !weak
			posteriors = [posteriors[i][:,3:end] for i in eachindex(posteriors)]
			tmp          = getmap.(posteriors)
			maxp         = DataFrame(vcat(tmp...),[:a,:gamNeg,:shape])
			al           = maxp[:,1:1]
			gam          = maxp[:,2:end]
		else
			tmp          = getmap.(posteriors)
			maxp         = DataFrame(vcat(tmp...),[:aw,:as,:a,:gamNeg,:shape])
			al           = maxp[:,1:3]
			gam          = maxp[:,4:end]
		end

		R"""al = as.data.table($al)
			lbls = if(ncol(al) > 1){c(expression(paste('Posterior ',alpha[w])), expression(paste('Posterior ',alpha[s])),expression(paste('Posterior ',alpha)))}else{c(expression(paste('Posterior ',alpha)))}
			clrs = if(ncol(al) > 1){c('#30504f', '#e2bd9a', '#ab2710')}else{c('#ab2710')}

			dal = suppressWarnings(data.table::melt(al))
			pal = suppressWarnings(ggplot(dal) + geom_density(aes(x=value,fill=variable),alpha=0.75) + scale_fill_manual('Posterior distribution',values=clrs ,labels=lbls) + theme_bw() + ggtitle($title) + xlab(expression(alpha)) + ylab(""))
			suppressWarnings(ggsave(pal,filename=paste0($analysisFolder,'/map.png'),dpi=600))
			"""
		CSV.write(analysisFolder * "/map.tsv",maxp,delim='\t',header=true)

		#=RCall.endEmbeddedR();=#

		return(maxp)
	catch
		println("Please install R, ggplot2, data.table and locfit in your system before execute this function")
	end
end


text:
module Excercise_4

# Exercise 4 (JULIA)
# 1. Create two matrices of the same layout and test if addition and subtraction of the matrix works as expected: C = A + B
# 2. Now compare matrix multiplication either this way A * B and this way A .* B. Whats the difference?!
# 3. What about matrix division with "/" or "\"?!
# 4. Create a 3x3 integer matrix A with useful numbers. Now try A+1, A-1, A*2, A/2.
# 5. Now multiply a 3x4 matrix with a suitable (4)vector.

# Part 1
matrix_a = rand(0:10, 4,4)
matrix_b = rand(0:10, 4,4)
println("Matrix A:")
println(matrix_a)
println("Matrix B:")
println(matrix_b)
println("A + B:")
addition = matrix_a + matrix_b
println(addition)
println("A - B:")
subtraction = matrix_a -matrix_b
println(subtraction)

# Part 2

# Matrix multiplication
println("A * B")
multi_1 = matrix_a * matrix_b
println(multi_1)

# Element-wise multiplication
println("A .* B:")
multi_2 = matrix_a .* matrix_b
println(multi_2)
# Part 3

# Part 4

matrix_c = [1 2 3; 3 2 1; 2 3 1]
println("Matrix C:")
println(matrix_c)
println("C+1:")
# println(matrix_c-1)
println("C*2:")
println(matrix_c*2)
println("C/2:")
println(matrix_c/2)

# Part 5

matrix_d = rand(0:10, 3,4)

println("Matrix D:")
println(matrix_d)

end


text:
using Plots
using DifferentialEquations

function generateStencil!(A, N, Δx)
    # Generate the centered difference matrix manually
    for i in 1:N, j in 1:N

    if abs(i-j)<=1
        (A[i,j]+=1)
    end
    if i==j
        (A[i,j]-=3)
    end

    end

    A = A/(Δx^2)
end

# Semi-linear heat equation: $u_t = \delta u + f$
f(t, u)  = zeros(size(u)) #ones(size(x,1)) - .5u
u0_func(x) = if (0.4 < x < 0.6) 1.0 else 0.0 end

# Discretize in time and space
Δx = 0.01
x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!

Δt = 0.1
t = Δt:Δt:1-Δt # Solve only for the interior: the endpoints are known to be zero!

N = length(x)
A = zeros(N,N)
generateStencil!(A, N, Δx)
u0 = u0_func.(x)

function heat_eq(u, p, t)
    du = A * u + f(t, u)
    du
end

function solve_prob(tmax)
    prob = ODEProblem(heat_eq, u0, (Δt, tmax - Δt))
    t = Δt:Δt:tmax-Δt
    sol = solve(prob, tstops=t)

    sol_f(t, x) = sol[t, x]
    plot(1:length(x), 1:length(t), (t, x) -> sol[t, x],
        title="Thermal Evolution of Box centered at 0.5",
        xlabel="Position Bins",
        ylabel="Time Bins",
        st=:heatmap,
        color=:viridis)
    savefig("plots/heat_example.png")
end