text:
# Accept either Symbol or String keys, but allways convert to Symbol
const Key = Union{Symbol,AbstractString}

"""
`AbstractGeoStack` objects hold multiple [`AbstractGeoArray`](@ref)
that share spatial bounds.

They are `NamedTuple`-like structures that may either contain `NamedTuple`
of [`AbstractGeoArray`](@ref), string paths that will load [`AbstractGeoArray`](@ref),
or a single path that points to as a file itself containing multiple layers, like
NetCDF or HDF5. Use and syntax is similar or identical for all cases.

`getindex` on a `AbstractGeoStack` generally returns a memory backed standard
[`GeoArray`](@ref). `geoarray[:somelayer] |> plot` plots the layers array,
while `geoarray[:somelayer, Lon(1:100), Band(2)] |> plot` will plot the
subset without loading the whole array.

`getindex` on a `AbstractGeoStack` with a key returns another stack with
getindex applied to all the arrays in the stack.
"""
abstract type AbstractGeoStack{T} end

(::Type{T})(data, keys; kwargs...) where T<:AbstractGeoStack =
    T(NamedTuple{Tuple(Symbol.(keys))}(Tuple(data)); kwargs...)

# Standard fields

refdims(s::AbstractGeoStack) = s.refdims
metadata(s::AbstractGeoStack) = s.metadata
window(s::AbstractGeoStack) = s.window

"""
    childkwargs(s::AbstractGeoStack)

Returns the keyword arguments that will be passed to the child array constructor.
"""
childkwargs(s::AbstractGeoStack) = s.childkwargs

# Interface methods ############################################################

"""
    getsource(s::AbstractGeoStack, [key])

Get the lower-level child object. This can be an `AbstractGeoArray` or
another object with GeoData interface methods defined.

Working with the low-level object can be better performance as we do not have to
processes everything needed to build a full `AbstractGeoArray`.
"""
function getsource end

"""
    modify(f, series::AbstractGeoStack)

Apply function `f` to the data of the child `AbstractGeoArray`s.

`f` must return an idenically sized array.

This method triggers a complete rebuild of all objects,
and disk based objects will be transferred to memory.

This is useful for swapping out array backend for an
entire stack to `CuArray` from CUDA.jl to copy data to a GPU,
and potentially other types like `DAarray` from Distributed.jl.
"""
modify(f, stack::AbstractGeoStack) =
    GeoStack(stack; data=mapdata(A -> modify(f, A), stack))

# Base methods #################################################################

@inline Base.getindex(s::AbstractGeoStack, I...) =
    rebuild(s; data=NamedTuple{keys(s)}(a[I...] for a in values(s)))
# Dict/Array hybrid with dims
@inline Base.getindex(s::AbstractGeoStack, key::Key, I::Vararg{<:Dimension}) =
    getindex(s, key, dims2indices(dims(s, key), I)...)

Base.values(s::AbstractGeoStack) = (s[key] for key in keys(s))
Base.length(s::AbstractGeoStack) = length(keys(s))
Base.keys(s::AbstractGeoStack{<:AbstractString}) = Symbol.(querychild(keys, getsource(s), s))
Base.keys(s::AbstractGeoStack{<:NamedTuple}) = Symbol.(keys(getsource(s)))
Base.names(s::AbstractGeoStack) = keys(s)
Base.map(f, s::AbstractGeoStack) = rebuild(s; data=mapdata(f, s))

mapdata(f, s::AbstractGeoStack) =
    NamedTuple{cleankeys(keys(s))}(map(f, values(s)))

"""
    Base.write(filename::AbstractString, T::Type{<:AbstractGeoArray}, s::AbstractGeoStack)

Save all layers of an `AbstractGeoStack` to separate files, using the backend determined
by `T`.

## Example

```julia
write(filename, GDALarray, A)
```
"""
Base.write(filename::AbstractString, ::Type{T}, s::AbstractGeoStack) where T <: AbstractGeoArray =
    for key in keys(s)
        base, ext = splitext(filename)
        fn = joinpath(string(base, "_", key, ext))
        write(fn, T, s[key])
    end

# Memory-based stacks ######################################################

"""
An [`AbstractGeoStack`](@ref) stored in memory.
"""
abstract type MemGeoStack{T} <: AbstractGeoStack{T} end

data(s::MemGeoStack) = s.data
data(s::MemGeoStack{<:NamedTuple}, key::Key) = data(s)[Symbol(key)]

rebuild(s::T; data=data(s), refdims=refdims(s), window=window(s),
        metadata=metadata(s), childkwargs=childkwargs(s), kwargs...) where T<:MemGeoStack =
    basetypeof(T)(data, refdims, window, metadata, childkwargs)

getsource(s::MemGeoStack{<:NamedTuple}, args...) = data(s, args...)

childdata(f, childobj, ::AbstractGeoStack) = f(childobj)

@inline Base.view(s::MemGeoStack, I...) =
    rebuild(s; data=NamedTuple{keys(s)}(view(a, I...) for a in values(s)))

@inline Base.getindex(s::MemGeoStack, key::Key) = begin
    A = rebuild(data(s)[key]; childkwargs(s)...)
    window_ = maybewindow2indices(dims(A), window(s))
    readwindowed(A, window_)
end
@inline Base.getindex(s::MemGeoStack, key::Key, i1::StandardIndices, I::StandardIndices...) = begin
    A = rebuild(data(s)[key]; childkwargs(s)...)
    window_ = maybewindow2indices(dims(A), window(s))
    readwindowed(A, window_, i1, I...)
end

# Disk-based stacks ######################################################

"""
[`AbstractGeoStack`](@ref)s stored on disk.
"""
abstract type DiskGeoStack{T} <: AbstractGeoStack{T} end

rebuild(s::T; data=filename(s), refdims=refdims(s), window=window(s),
        metadata=metadata(s), childtype=childtype(s), childkwargs=childkwargs(s), kwargs...) where T<:DiskGeoStack =
    basetypeof(T)(data, refdims, window, metadata, childtype, childkwargs)

getsource(s::DiskGeoStack, args...) = filename(s, args...)

childtype(stack::DiskGeoStack) = stack.childtype

"""
    filename(s::DiskGeoStack)

Return the filename field of a `DiskGeoStack`.
This may be a `Vector` of `String`, or a `String`.
"""
filename(s::DiskGeoStack) = s.filename
"""
    filename(s::DiskGeoStack, key)

Return the filename field of a `DiskGeoStack` for a given key.

This will always be a single string. However, in some cases
all keys may have the same filename.
"""
filename(s::DiskGeoStack{<:NamedTuple}, key::Key) = filename(s)[Symbol(key)]
filename(s::DiskGeoStack, key::Key) = filename(s)

# Default to using the same method as GeoStack.
# But a method for withsourcedata dispatching on the
# array object may save time in some cases, like when
# the array size has to be determined before loading the data
# like with mmaped grd files.
withsourcedata(f, A::DiskGeoArray, key...) =
    withsourcedata(f, typeof(A), filename(A), key...)
# By default assume the source object is a single object with
# all info and data.
withsourcedata(f, childtype::Type, path, key...) =
    withsource(f, childtype, path, key...)

withsource(f, childtype::Type, path, key...) = f(path)

# Base methods

Base.getindex(s::DiskGeoStack, key::Key) = begin
    filename_ = filename(s, key)
    withsource(childtype(s), filename_, key) do dataset
        A = childtype(s)(dataset, filename_, key; name=Symbol(key), childkwargs(s)...)
        window_ = maybewindow2indices(A, window(s))
        readwindowed(A, window_)
    end
end
Base.getindex(s::DiskGeoStack, key::Key, i1::StandardIndices, I::StandardIndices...) = begin
    filename_ = filename(s, key)
    withsource(childtype(s), filename_, key) do dataset
        A = childtype(s)(dataset, filename_, key; name=Symbol(key), childkwargs(s)...)
        window_ = maybewindow2indices(dims(A), window(s))
        readwindowed(A, window_, i1, I...)
    end
end

@inline Base.view(s::DiskGeoStack, I...) = rebuild(s; window=I)

# Default dims, metadata and missingval methods
#
# For DiskStack we query the underlying object - avoiding building
# an AbstractGeoArray unless we have to. Examples may be an NCDatasets,
# ArchGDAL, or HDF5 `Dataset` object. These sources will add a
# For MemGeoStack the childobj is just a GeoArray as it is allready loaded in memory.
for func in (:dims, :metadata, :missingval)
    @eval begin
        $func(s::AbstractGeoStack, key::Key) = $func(s[key])
        $func(s::DiskGeoStack, key::Key) =
            withsource(childtype(s), getsource(s, key), key) do sourceobj
                $func(sourceobj, key)
            end
    end
end

# Concrete MemGeoStack implementation ######################################################

"""
    GeoStack(data...; keys, kwargs...)
    GeoStack(data::Union{Vector,Tuple}; keys, kwargs...)
    GeoStack(data::NamedTuple; window=(), metadata=nothing, refdims=(), childkwargs=()) =
    GeoStack(s::AbstractGeoStack; [keys, data, refdims, window, metadata])

A concrete `MemGeoStack` implementation. Holds layers of [`GeoArray`](@ref).

## Argumenst

- `data`: A `NamedTuple` of [`GeoArray`](@ref), or a `Vector`, `Tuple` or splatted arguments
  of [`GeoArray`](@ref). The latter options must pass a `keys` keyword argument.

## Keyword Argumenst

- `keys`: Used as stack keys when a `Tuple` or `Vector` or splat of geoarrays are passed in.
- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
  contained arrays when they are accessed.
- `refdims`: Reference dimensions from earlier subsetting.
- `metadata`: Metadata as a [`StackMetadata`](@ref) object.
- `childkwargs`: A `NamedTuple` of keyword arguments to pass to the constructor.
- `refdims`: `Tuple` of  position `Dimension` the array was sliced from.
"""
struct GeoStack{T,R,W,M,K} <: MemGeoStack{T}
    data::T
    refdims::R
    window::W
    metadata::M
    childkwargs::K
end
GeoStack(data::AbstractGeoArray...; keys=name.(data), kwargs...) =
    GeoStack(NamedTuple{cleankeys(keys)}(data); kwargs...)
GeoStack(data::NamedTuple;
         refdims=(),
         window=(),
         metadata=nothing,
         childkwargs=()) =
    GeoStack(data, refdims, window, metadata, childkwargs)
GeoStack(s::AbstractGeoStack; keys=cleankeys(Base.keys(s)),
         data=NamedTuple{keys}(s[key] for key in keys),
         refdims=refdims(s),
         window=(),
         metadata=metadata(s),
         childkwargs=()) =
    GeoStack(data, refdims, window, metadata, childkwargs)

Base.convert(::Type{GeoStack}, src::AbstractGeoStack) = GeoStack(src)

# Concrete DiskGeoStack implementation ######################################################

"""
    DiskStack(filenames...; keys, kwargs...)
    DiskStack(filenames; keys, kwargs...)
    DiskStack(filenames::NamedTuple;
              window=(),
              metadata=nothing,
              childtype,
              childkwargs=()
              refdims=())

Concrete [`DiskGeoStack`](@ref) implementation. Loads a stack of files lazily from disk.

## Arguments

- `filename`: a NamedTuple of stack keys and `String` filenames.

## Keyword arguments

- `keys`: Used as stack keys when a `Tuple`, `Vector` or splat of filenames are passed in.
- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
  contained arrays when they are accessed.
- `metadata`: Metadata as a [`StackMetadata`](@ref) object.
- `childtype`: The type of the child data. eg. `GDALarray`. Required.
- `childkwargs`: A `NamedTuple` of keyword arguments to pass to the `childtype` constructor.
- `refdims`: `Tuple` of  position `Dimension` the array was sliced from.
"""
struct DiskStack{T,R,W,M,C,K} <: DiskGeoStack{T}
    filename::T
    refdims::R
    window::W
    metadata::M
    childtype::C
    childkwargs::K
end
DiskStack(filenames::NamedTuple; refdims=(), window=(), metadata=nothing, childtype, childkwargs=()) =
    DiskStack(filenames, refdims, window, metadata, childtype, childkwargs)
DiskStack(filenames...; kwargs...) = DiskStack(filenames; kwargs...)
DiskStack(filenames; keys, kwargs...) =
    DiskStack(NamedTuple{cleankeys(keys)}((filenames...,)); kwargs...)

# Other base methods

Base.copy(stack::AbstractGeoStack) =
    rebuild(stack; data=map(copy, getsource(stack)))

"""
    Base.copy!(dst::AbstractGeoStack, src::AbstractGeoStack, [keys=keys(dst)])

Copy all or a subset of layers from one stack to another.

## Example

Copy just the `:sea_surface_temp` and `:humidity` layers from `src` to `dst`.

```julia
copy!(dst::AbstractGeoStack, src::AbstractGeoStack, keys=(:sea_surface_temp, :humidity))
```
"""
Base.copy!(dst::MemGeoStack, src::AbstractGeoStack, keys=keys(dst)) = begin
    symkeys = Symbol.(keys)
    # Check all keys first so we don't copy anything if there is any error
    for key in symkeys
        key in Symbol.(Base.keys(dst)) || throw(ArgumentError("key $key not found in dest keys"))
        key in Symbol.(Base.keys(src)) || throw(ArgumentError("key $key not found in source keys"))
    end
    for key in symkeys
        copy!(dst[key], src[key])
    end
end
"""
    Base.copy!(dst::AbstractArray, src::DiskGeoStack, key::Key)

Copy the stack layer `key` to `dst`, which can be any `AbstractArray`.

## Example

Copy the `:humidity` layer from `stack` to `array`.

```julia
copy!(array, stack, :humidity)
```
"""
Base.copy!(dst::AbstractArray, src::AbstractGeoStack, key) =
    copy!(dst, src[key])

"""
    Base.cat(stacks::AbstractGeoStack...; [keys=keys(stacks[1])], dims)

Concatenate all or a subset of layers for all passed in stacks.

## Keyword Arguments

- `keys`: `Tuple` of `Symbol` for the stack keys to concatenate.
- `dims`: Dimension of child array to concatenate on.

## Example

Concatenate the :sea_surface_temp and :humidity layers in the time dimension:

```julia
cat(stacks...; keys=(:sea_surface_temp, :humidity), dims=Ti)
```
"""
Base.cat(stacks::AbstractGeoStack...; keys=keys(stacks[1]), dims) = begin
    vals = Tuple(cat((s[key] for s in stacks)...; dims=dims) for key in keys)
    GeoStack(stacks[1], data=NamedTuple{keys}(vals))
end
Base.first(s::AbstractGeoStack) = s[first(keys(s))]
Base.last(s::AbstractGeoStack) = s[last(keys(s))]


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for color.jl

@testset "operations" begin
    c1 = RGB(1f0, 2f0, 3f0)
    c2 = RGB(0.4f0, 0.5f0, 0.6f0)

    @testset "addition" begin
        @test c1 + c2 == RGB(1.4f0, 2.5f0, 3.6f0)
    end

    @testset "subtraction" begin
        @test c1 - c2 == RGB(0.6f0, 1.5f0, 2.4f0)
    end

    @testset "scalar multiplication" begin
        a = 2
        @test a * c1 == RGB(2f0, 4f0, 6f0)
        @test c1 * a == RGB(2f0, 4f0, 6f0)
        @test a * c1 == c1 * a
    end

    @testset "elementwise ≈" begin
        a = 0
        for i in 1:10
            a += 0.1f0
        end
        @test RGB(a, 2f0 * a, 3f0 * a) ≈ c1
    end

    @testset "elementwise multiplication" begin
        @test c1 * c2 == RGB(1f0 * 0.4f0, 2f0 * 0.5f0, 3f0 * 0.6f0)
    end
end

@testset "iterations" begin
    r = 1f0
    g = 2f0
    b = 3f0
    c = RGB(r, g, b)

    @testset "indexing properties" begin
        @test length(c) === 3
        @test firstindex(c) === 1
        @test lastindex(c) === 3
    end

    @testset "get index" begin
        # linear indexing
        @test r == c[begin]
        @test r == c[1]
        @test g == c[2]
        @test b == c[end]
        @test b == c[3]
        @test_throws BoundsError c[4]

        # cartesian indexing
        @test r == c[CartesianIndex(1)]
        @test g == c[CartesianIndex(2)]
        @test b == c[CartesianIndex(3)]
        @test_throws BoundsError c[CartesianIndex(4)]
    end

    @testset "iterability" begin
        @test all(i == j for (i, j) in zip((r, g, b), c))
    end

    @testset "splat operator" begin
        cc = RGB(c...)
        @test cc === c
    end
end

@testset "broadcasting" begin
    c1 = RGB(1f0, 2f0, 3f0)
    c2 = RGB(4f0, 5f0, 6f0)
    a = 2

    # testing equivalence to custom defined methods
    @test c1 .+ c2 == c1 + c2
    @test c1 .- c2 == c1 - c2
    @test c1 .* c2 == c1 * c2
    @test a .* c1 == a * c1
    @test c1 .* a == c1 * a

    # broadcasting operators can be applied between any broadcastable type instances
    @test all((0.1f0, 0.2f0, 0.3f0) .+ c1 ≈ RGB(1.1f0, 2.2f0, 3.3f0))
    @test all((1f0, 2f0, 3f0) .== c1)

    # it works for any operator valid for the types of the elements
    @test all(c2 ./ c1 ≈ RGB(4f0, 5f0/2f0, 2f0))
end

@testset "color manipulation" begin
    c = RGB(1f0, 2f0, 3f0)

    @test luminosity(c) ≈ 2f0
    @test clamp(c) ≈ RGB(1f0/2f0, 2f0/3f0, 3f0/4f0)
    @test γ_correction(c, 1f0) ≈ c
    @test γ_correction(c, 1.2f0) ≈ RGB(1f0^(1/1.2), 2f0^(1/1.2), 3f0^(1/1.2))
end

@testset "IO" begin
    endian_f = ENDIAN_BOM == 0x04030201 ? ltoh : ntoh

    io = IOBuffer()
    c = RGB(1f0, 2f0, 3f0)

    @testset "show" begin
        # compact
        show(io, c)
        @test String(take!(io)) == "(1.0 2.0 3.0)"
        # extended
        show(io, "text/plain", c)
        @test String(take!(io)) == "RGB color with eltype Float32\nR: 1.0, G: 2.0, B: 3.0"
    end
end

@testset "other" begin
    @testset "eltype" begin
        @test eltype(RGB{Float32}) == Float32
    end

    @testset "zero" begin
        # from type
        @test zero(RGB{Float32}) == RGB(0f0, 0f0, 0f0)
        # from variable
        c = RGB(1f0, 2f0, 3f0)
        @test zero(c) == RGB(0f0, 0f0, 0f0)
    end

    @testset "one" begin
        # from type
        @test one(RGB{Float32}) == RGB(1f0, 1f0, 1f0)
        # from variable
        c = RGB(1f0, 2f0, 3f0)
        @test one(c) == RGB(1f0, 1f0, 1f0)
    end
end

ext: jl


text:
module Alpaca

    ##############################################################################
    ##
    ## Dependencies
    ##
    ##############################################################################

    using Base
    using Statistics
    using RCall
    using DataFrames
    using Distributions
    using Printf
    import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var,
        RegressionModel, model_response, stderror, confint, fit, CoefTable,
        dof_residual, r2, adjr2, deviance, mss, rss, islinear, response
    using Reexport
    using CategoricalArrays
    @reexport using StatsBase
    @reexport using StatsModels

    ##############################################################################
    ##
    ## Exported methods and types
    ##
    ##############################################################################

    export feglm

    ##############################################################################
    ##
    ## Load files
    ##
    ##############################################################################

    include("RegressionResult.jl")
    include("feglm.jl")

    # instantiate R's alpaca package
    # R"if(\"alpaca\" %in% rownames(installed.packages()) == FALSE) {install.packages(\"alpaca\")}"
    # R"if(\"formula.tools\" %in% rownames(installed.packages()) == FALSE) {install.packages(\"formula.tools\")}"

end


text:
function champs_4(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr,
                  conjmatr, sigma, llambda, Ktaubis, epsilon, Ntaumm)

    h1u, h1v = ftau(t,
		    fft_ubar, fft_vbar,
		    A1, A2,
		    matr, conjmatr,
		    sigma, llambda)

    h1u_fft = fft(h1u, dims=1)
    h1v_fft = fft(h1v, dims=1)

    dtu0u = h1u_fft[1, :] / Ntaumm
    dtu0v = h1v_fft[1, :] / Ntaumm

    h1u_fft[1, :] .= 0
    h1v_fft[1, :] .= 0

    h1u = epsilon * ifft(h1u_fft ./ (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
                              fft_ubar, fft_vbar,
			      dtu0u, dtu0v,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar, fft_vbar,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2
    champv = champv1 .+ champv2

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    dttu0u = champu_fft[1, :] / Ntaumm
    dttu0v = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth1u = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth1v = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu, champv = ftau(t,
                          fft_ubar .+ h1u, fft_vbar .+ h1v,
			  A1, A2, matr,
                          conjmatr, sigma, llambda)
    h2u = champu .- dth1u
    h2v = champv .- dth1v

    h2u_fft = fft(h2u, dims=1)
    h2v_fft = fft(h2v, dims=1)

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    dtu1u = champu_fft[1, :] / Ntaumm
    dtu1v = champv_fft[1, :] / Ntaumm

    h2u_fft[1, :] .= 0.0
    h2v_fft[1, :] .= 0.0

    h2u = epsilon * ifft(h2u_fft ./ (1im * Ktaubis), dims=1)
    h2v = epsilon * ifft(h2v_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = d2uftau(t, fft_ubar, fft_vbar,
                               dtu0u, dtu0v, dtu0u, dtu0v,
                               A1, A2, matr, conjmatr, sigma,
                               llambda, Ntaumm)

    champu2, champv2 = dtduftau(t,
				fft_ubar, fft_vbar,
				dtu0u, dtu0v,
                                A1, A2,
				matr, conjmatr,
				sigma, llambda, Ntaumm)

    champu3, champv3 = duftau(t,
			      fft_ubar, fft_vbar,
			      dttu0u, dttu0v,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu4, champv4 = d2tftau(t,
			       fft_ubar, fft_vbar,
			       A1, A2,
                               matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu = champu1 .+ 2 * champu2 .+ champu3 .+ champu4
    champv = champv1 .+ 2 * champv2 .+ champv3 .+ champv4

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dtth1u = epsilon * ifft(champu_fft / (1im * Ktaubis), dims=1)
    dtth1v = epsilon * ifft(champv_fft / (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
                              dtu1u .+ dth1u,
			      dtu1v .+ dth1v,
			      A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .- dtth1u
    champv = champv1 .+ champv2 .- dtth1v

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth2u = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth2v = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    h3u, h3v = ftau(t,
		    fft_ubar .+ h2u,
		    fft_vbar .+ h2v,
		    A1, A2,
		    matr, conjmatr,
		    sigma, llambda)

    h3u = h3u .- dth2u
    h3v = h3v .- dth2v

    h3u_fft = fft(h3u, dims=1)
    h3v_fft = fft(h3v, dims=1)

    h3u_fft[1, :] .= 0
    h3v_fft[1, :] .= 0

    dtauh3u = ifft(h3u_fft, dims=1)
    dtauh3v = ifft(h3v_fft, dims=1)

    h3u = epsilon * ifft(h3u_fft ./ (1im * Ktaubis), dims=1)
    h3v = epsilon * ifft(h3v_fft ./ (1im * Ktaubis), dims=1)

    C3u = fft_ubar .+ h3u
    C3v = fft_vbar .+ h3v

    ffu, ffv = ftau(t,
		    C3u .+ fft_ug,
		    C3v .+ fft_vg,
		    A1, A2,
		    matr, conjmatr,
                    sigma, llambda)

    champu, champv = ftau(t, C3u, C3v, A1, A2, matr, conjmatr,
			  sigma, llambda)

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champubaru = champu_fft[1, :] / Ntaumm
    champubarv = champv_fft[1, :] / Ntaumm

    champu1, champv1 = duftau(t, fft_ubar, fft_vbar,
                              champubaru, champubarv, A1, A2,
                              matr, conjmatr, sigma, llambda,
                              Ntaumm)

    champu2, champv2 = dtftau(t, fft_ubar, fft_vbar, A1, A2,
                              matr, conjmatr, sigma, llambda, Ntaumm)

    champu_fft = fft(champu1 .+ champu2, dims=1)
    champv_fft = fft(champv1 .+ champv2, dims=1)

    dttu0baru = champu_fft[1, :] / Ntaumm
    dttu0barv = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth1baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth1barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
                              champubaru .+ dth1baru,
			      champubarv .+ dth1barv,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu_fft = fft(champu1 .+ champu2, dims=1)
    champv_fft = fft(champv1 .+ champv2, dims=1)

    dttu1baru = champu_fft[1, :] / Ntaumm
    dttu1barv = champv_fft[1, :] / Ntaumm

    champu1, champv1 = d2uftau(t,
			       fft_ubar, fft_vbar,
			       dtu0u, dtu0v,
                               champubaru, champubarv,
			       A1, A2, matr, conjmatr,
                               sigma, llambda, Ntaumm)

    champu2, champv2 = duftau(t,
			      fft_ubar, fft_vbar,
			      dttu0baru, dttu0barv,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda,
                              Ntaumm)

    champu3, champv3 = dtduftau(t,
				fft_ubar, fft_vbar,
				dtu0u, dtu0v,
                                A1, A2,
				matr, conjmatr,
				sigma, llambda, Ntaumm)

    champu4, champv4 = dtduftau(t,
				fft_ubar, fft_vbar,
				champubaru, champubarv,
                                A1, A2,
				matr, conjmatr,
				sigma, llambda,
                                Ntaumm)

    champu5, champv5 = d2tftau(t,
			       fft_ubar, fft_vbar,
			       A1, A2,
			       matr, conjmatr,
                               sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .+ champu3 .+ champu4 .+ champu5
    champv = champv1 .+ champv2 .+ champv3 .+ champv4 .+ champv5

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    dtttu0baru = champu_fft[1, :] / Ntaumm
    dtttu0barv = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dtth1baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dtth1barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = d3uftau(t,
			       fft_ubar, fft_vbar,
			       dtu0u, dtu0v,
                               dtu0u, dtu0v,
			       champubaru, champubarv,
			       A1, A2,
			       matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu2, champv2 = dtd2uftau(t,
				 fft_ubar, fft_vbar,
				 dtu0u, dtu0v,
                                 dtu0u, dtu0v,
				 A1, A2,
				 matr, conjmatr,
				 sigma, llambda, Ntaumm)

    champu3, champv3 = d2uftau(t,
			       fft_ubar, fft_vbar,
                               dttu0baru, dttu0barv,
			       dtu0u, dtu0v,
			       A1, A2,
                               matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu4, champv4 = dtd2uftau(t,
				 fft_ubar, fft_vbar,
				 dtu0u, dtu0v,
                                 champubaru, champubarv,
				 A1, A2, matr, conjmatr,
                                 sigma, llambda, Ntaumm)

    champu5, champv5 = d2tduftau(t,
				 fft_ubar, fft_vbar,
				 dtu0u, dtu0v,
                                 A1, A2,
				 matr, conjmatr,
				 sigma, llambda, Ntaumm)

    champu6, champv6 = dtduftau(t, fft_ubar, fft_vbar,
                                dttu0baru, dttu0barv,
                                A1, A2,
				matr, conjmatr,
				sigma, llambda,
                                Ntaumm)

    champu7, champv7 = d2uftau(t,
			       fft_ubar, fft_vbar,
                               dttu0u, dttu0v,
			       champubaru, champubarv,
                               A1, A2,
			       matr, conjmatr,
                               sigma, llambda, Ntaumm)

    champu8, champv8 = dtduftau(t,
				fft_ubar, fft_vbar,
                                dttu0u, dttu0v,
				A1, A2,
				matr, conjmatr,
				sigma, llambda, Ntaumm)

    champu9, champv9 = duftau(t,
			      fft_ubar, fft_vbar,
                              dtttu0baru, dtttu0barv,
			      A1, A2,
			      matr, conjmatr,
			      sigma, llambda,
                              Ntaumm)

    champu10, champv10 = d3tftau(t,
				 fft_ubar, fft_vbar,
				 A1, A2,
                                 matr, conjmatr,
				 sigma, llambda, Ntaumm)

    champu11, champv11 = d2tduftau(t,
				   fft_ubar, fft_vbar,
                                   champubaru, champubarv,
				   A1, A2,
				   matr, conjmatr,
				   sigma, llambda, Ntaumm)

    champu = (champu1 .+ champu2 .+ 2 * champu3 .+ 2 * champu4
                     .+ 2 * champu5 .+ 2 * champu6 .+ champu7
		     .+ champu8 .+ champu9 .+ champu10 .+ champu11)

    champv = (champv1 .+ champv2 .+ 2 * champv3 .+ 2 * champv4
             .+ 2 * champv5 .+ 2 * champv6 .+ champv7 .+ champv8
	     .+ champv9 .+ champv10 .+ champv11)

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dttth1baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dttth1barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = d2uftau(t,
			       fft_ubar .+ h1u, fft_vbar .+ h1v,
                               dtu1u .+ dth1u, dtu1v .+ dth1v,
                               champubaru .+ dth1baru,
			       champubarv .+ dth1barv,
                               A1, A2,
			       matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu2, champv2 = duftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
                              dttu1baru .+ dtth1baru,
			      dttu1barv .+ dtth1barv,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu3, champv3 = dtduftau(t,
				fft_ubar .+ h1u,
				fft_vbar .+ h1v,
                                dtu1u .+ dth1u,
				dtu1v .+ dth1v,
				A1, A2,
				matr, conjmatr,
				sigma, llambda, Ntaumm)

    champu4, champv4 = dtduftau(t,
				fft_ubar .+ h1u, fft_vbar .+ h1v,
                                champubaru .+ dth1baru,
                                champubarv .+ dth1barv,
				A1, A2,
				matr, conjmatr,
				sigma, llambda,
                                Ntaumm)

    champu5, champv5 = d2tftau(t,
			       fft_ubar .+ h1u, fft_vbar .+ h1v,
			       A1, A2,
                               matr, conjmatr, sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .+ champu3 .+ champu4 .+ champu5 .- dttth1baru
    champv = champv1 .+ champv2 .+ champv3 .+ champv4 .+ champv5 .- dttth1barv

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dtth2baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dtth2barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
                              champubaru .+ dth1baru,
			      champubarv .+ dth1barv,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu_fft = fft(champu1 .+ champu2 .- dtth1baru, dims=1)
    champv_fft = fft(champv1 .+ champv2 .- dtth1barv, dims=1)

    dttu0baru = champu_fft[1, :] / Ntaumm
    dttu0barv = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth2baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth2barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar .+ h2u,
			      fft_vbar .+ h2v,
                              champubaru .+ dth2baru,
			      champubarv .+ dth2barv,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar .+ h2u,
			      fft_vbar .+ h2v,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .- dtth2baru
    champv = champv1 .+ champv2 .- dtth2barv

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dtC3u = champubaru .+ epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dtC3v = champubarv .+ epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champgu = ffu .- dtauh3u .- dtC3u
    champgv = ffv .- dtauh3v .- dtC3v

    champgu_fft = fft(champgu, dims=1)
    champgv_fft = fft(champgv, dims=1)

    champmoyu = champgu_fft[1, :] / Ntaumm
    champmoyv = champgv_fft[1, :] / Ntaumm

    champgu_fft[1, :] .= 0
    champgv_fft[1, :] .= 0

    ichampgu = ifft(champgu_fft ./ (1im * Ktaubis), dims=1)
    ichampgv = ifft(champgv_fft ./ (1im * Ktaubis), dims=1)

    champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv

end


text:
# Slugify.jl -- A library that simplifies a text to an ASCII subset
# By: Emmanuel Raviart <emmanuel@raviart.com>
#
# Copyright (C) 2015 Emmanuel Raviart
# https://github.com/eraviart/Slugify.jl
#
# This file is part of Slugify.jl.
#
# The Slugify.jl package is licensed under the MIT "Expat" License.

module Slugify

export slugify

ASCII_TRANSLATIONS = Dict{Char,String}(
  ' ' => " ",  # U+00A0 NO-BREAK SPACE
  'À' => "A",  # U+00C0 LATIN CAPITAL LETTER A WITH GRAVE
  'Á' => "A",  # U+00C1 LATIN CAPITAL LETTER A WITH ACUTE
  'Â' => "A",  # U+00C2 LATIN CAPITAL LETTER A WITH CIRCUMFLEX
  'Ã' => "A",  # U+00C3 LATIN CAPITAL LETTER A WITH TILDE
  'Ä' => "A",  # U+00C4 LATIN CAPITAL LETTER A WITH DIAERESIS
  'Å' => "A",  # U+00C5 LATIN CAPITAL LETTER A WITH RING ABOVE
  'Æ' => "Ae",  # U+00C6 LATIN CAPITAL LETTER AE
  'Ç' => "C",  # U+00C7 LATIN CAPITAL LETTER C WITH CEDILLA
  'È' => "E",  # U+00C8 LATIN CAPITAL LETTER E WITH GRAVE
  'É' => "E",  # U+00C9 LATIN CAPITAL LETTER E WITH ACUTE
  'Ê' => "E",  # U+00CA LATIN CAPITAL LETTER E WITH CIRCUMFLEX
  'Ë' => "E",  # U+00CB LATIN CAPITAL LETTER E WITH DIAERESIS
  'Ì' => "I",  # U+00CC LATIN CAPITAL LETTER I WITH GRAVE
  'Í' => "I",  # U+00CD LATIN CAPITAL LETTER I WITH ACUTE
  'Î' => "I",  # U+00CE LATIN CAPITAL LETTER I WITH CIRCUMFLEX
  'Ï' => "I",  # U+00CF LATIN CAPITAL LETTER I WITH DIAERESIS
  'Ð' => "Th",  # U+00D0 LATIN CAPITAL LETTER ETH
  'Ñ' => "N",  # U+00D1 LATIN CAPITAL LETTER N WITH TILDE
  'Ò' => "O",  # U+00D2 LATIN CAPITAL LETTER O WITH GRAVE
  'Ó' => "O",  # U+00D3 LATIN CAPITAL LETTER O WITH ACUTE
  'Ô' => "O",  # U+00D4 LATIN CAPITAL LETTER O WITH CIRCUMFLEX
  'Õ' => "O",  # U+00D5 LATIN CAPITAL LETTER O WITH TILDE
  'Ö' => "O",  # U+00D6 LATIN CAPITAL LETTER O WITH DIAERESIS
  'Ø' => "O",  # U+00D8 LATIN CAPITAL LETTER O WITH STROKE
  'Ù' => "U",  # U+00D9 LATIN CAPITAL LETTER U WITH GRAVE
  'Ú' => "U",  # U+00DA LATIN CAPITAL LETTER U WITH ACUTE
  'Û' => "U",  # U+00DB LATIN CAPITAL LETTER U WITH CIRCUMFLEX
  'Ü' => "U",  # U+00DC LATIN CAPITAL LETTER U WITH DIAERESIS
  'Ý' => "Y",  # U+00DD LATIN CAPITAL LETTER Y WITH ACUTE
  'Þ' => "th",  # U+00DE LATIN CAPITAL LETTER THORN
  'ß' => "ss",  # U+00DF LATIN SMALL LETTER SHARP S
  'à' => "a",  # U+00E0 LATIN SMALL LETTER A WITH GRAVE
  'á' => "a",  # U+00E1 LATIN SMALL LETTER A WITH ACUTE
  'â' => "a",  # U+00E2 LATIN SMALL LETTER A WITH CIRCUMFLEX
  'ã' => "a",  # U+00E3 LATIN SMALL LETTER A WITH TILDE
  'ä' => "a",  # U+00E4 LATIN SMALL LETTER A WITH DIAERESIS
  'å' => "a",  # U+00E5 LATIN SMALL LETTER A WITH RING ABOVE
  'æ' => "ae",  # U+00E6 LATIN SMALL LETTER AE
  'ç' => "c",  # U+00E7 LATIN SMALL LETTER C WITH CEDILLA
  'è' => "e",  # U+00E8 LATIN SMALL LETTER E WITH GRAVE
  'é' => "e",  # U+00E9 LATIN SMALL LETTER E WITH ACUTE
  'ê' => "e",  # U+00EA LATIN SMALL LETTER E WITH CIRCUMFLEX
  'ë' => "e",  # U+00EB LATIN SMALL LETTER E WITH DIAERESIS
  'ì' => "i",  # U+00EC LATIN SMALL LETTER I WITH GRAVE
  'í' => "i",  # U+00ED LATIN SMALL LETTER I WITH ACUTE
  'î' => "i",  # U+00EE LATIN SMALL LETTER I WITH CIRCUMFLEX
  'ï' => "i",  # U+00EF LATIN SMALL LETTER I WITH DIAERESIS
  'ð' => "th",  # U+00F0 LATIN SMALL LETTER ETH
  'ñ' => "n",  # U+00F1 LATIN SMALL LETTER N WITH TILDE
  'ò' => "o",  # U+00F2 LATIN SMALL LETTER O WITH GRAVE
  'ó' => "o",  # U+00F3 LATIN SMALL LETTER O WITH ACUTE
  'ô' => "o",  # U+00F4 LATIN SMALL LETTER O WITH CIRCUMFLEX
  'õ' => "o",  # U+00F5 LATIN SMALL LETTER O WITH TILDE
  'ö' => "o",  # U+00F6 LATIN SMALL LETTER O WITH DIAERESIS
  'ø' => "o",  # U+00F8 LATIN SMALL LETTER O WITH STROKE
  'ù' => "u",  # U+00F9 LATIN SMALL LETTER U WITH GRAVE
  'ú' => "u",  # U+00FA LATIN SMALL LETTER U WITH ACUTE
  'û' => "u",  # U+00FB LATIN SMALL LETTER U WITH CIRCUMFLEX
  'ü' => "u",  # U+00FC LATIN SMALL LETTER U WITH DIAERESIS
  'ý' => "y",  # U+00FD LATIN SMALL LETTER Y WITH ACUTE
  'þ' => "th",  # U+00FE LATIN SMALL LETTER THORN
  'ÿ' => "y",  # U+00FF LATIN SMALL LETTER Y WITH DIAERESIS
  'Œ' => "Oe",  # U+0152 LATIN CAPITAL LIGATURE OE
  'œ' => "oe",  # U+0153 LATIN SMALL LIGATURE OE
  '‘' => "'",  # U+2018 LEFT SINGLE QUOTATION MARK
  '’' => "'",  # U+2019 RIGHT SINGLE QUOTATION MARK
)

function char_to_ascii(unicode_char::Char)
  """Convert an unicode character to a string of several ASCII characters."""
  return get(ASCII_TRANSLATIONS, unicode_char) do
    return unicode_char < Char(0x80) ? string(unicode_char) : ""
  end
end

function slugify(string::AbstractString; separator::Char = '-', transform::Function = lowercase)
  """Simplify a string, converting it to an ASCII subset of 0-9, A-Z, a-z & separators."""
  simplified = join([
    slugify_char(char, separator = separator, transform = transform)
    for char in string
  ], "")
  return join(split(simplified, separator, keepempty=false), separator)
end

function slugify_char(unicode_char::Char; separator::Char = '-', transform::Function = lowercase)
  """Convert an unicode character to a subset of ASCII characters or an empty string.

  The result can be composed of several characters (for example, 'œ' becomes 'oe').
  """
  return map(transform(char_to_ascii(unicode_char))) do char
    return char in "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz" ? char : separator
  end
end

end # module


text:
### A Pluto.jl notebook ###
# v0.17.4

using Markdown
using InteractiveUtils

# ╔═╡ cdb3b623-4b24-49a9-b799-ee349287715c
using PlutoUI

# ╔═╡ 855f2a3b-b7ce-4b3d-9ebd-359f8e13bfc2
bigbreak = HTML("<br/>"^4);

# ╔═╡ b2e97780-801b-11ec-0433-c7952bc1387d
md"""
# Getting started

Topics:
  1. How to print
  2. How to assign variables
  3. Ho to comment
  4. Syntax for basic math
"""

# ╔═╡ 2f8bea58-1b67-4ff0-80b9-70f17c5eb26b
md"""
## How to print

In Julia we usually use `println()` to print
"""

# ╔═╡ feb26c6b-f169-4543-be6e-fd48c4b37f3b
with_terminal() do
	println("I'm excited to learn Julia!")
end

# ╔═╡ f4fbcb5f-8744-45d3-a5a8-acef9efd3a79
md"""
_Note: In Pluto, the terminal is often hidden. But Pluto will show the last value of each cell anyway._
"""

# ╔═╡ 18c5fa07-8e59-4c25-be5c-90fd7406b20a
"I'm excited to learn Julia!"

# ╔═╡ a9debd68-5e03-4c69-962d-87449024b3e8
md"""
## How to assign variables

All we need is a variable name, value, and an equal's sign!
Julia will figure out types for us.
"""

# ╔═╡ 4df6e513-6b77-4e6b-8870-d9cf84ebdfc7
my_answer = 42

# ╔═╡ c6dc20f1-d2bd-4f2f-b536-8a1560e928ed
typeof(my_answer)

# ╔═╡ 71777029-f294-42c6-963c-8b6bc15610bc
my_pi = 3.14159

# ╔═╡ 593e8d96-7a1c-47cc-9e23-7aba58867d5a
typeof(my_pi)

# ╔═╡ c66ebcb1-d3f4-4208-ad86-ac3745b8a386
😺 = "smiley cat!"

# ╔═╡ 3a9a5ec6-4ab8-4c44-a621-76a3c4e08233
typeof(😺)

# ╔═╡ 0b7b3b48-18da-426b-91f9-54f8bac5fd76
md"_To type a smiley cat, use **'\' + tab completion** to select the emoji name and then tab again_
"

# ╔═╡ 4ec52ccd-6f3a-4027-87bd-54487ff75562
#\:smi + <tab> --> select with down arrow + <enter> ---> <tab> + <enter> to complete

# ╔═╡ cba0e0e6-531c-4b90-a36f-480b6c5b9d9a
md"""After assigning a value to a variable, we can reassign a value of a different type to that variable without any issue.

_Note: This is the second point in which Pluto is a bit different. Since these notebooks are reactive, you cannot reassign. (But this has many advantages later).
However, we can simulate a normal Julia enviroment by adding a `begin` `end` block._
"""

# ╔═╡ f4e86b43-263f-4240-8c10-67987bfd1103
begin
	🐶 = 1.0
	T1 = typeof(🐶)

	🐶 = 1
	T2 = typeof(🐶)

	"First 🐶 was of type $T1, then is was of type $T2"
end

# ╔═╡ ee1e0582-c26e-4550-a3a8-07d4eb1919b8
md"_Note: until now, we never had to define the type. Julia allows us to write **generic** code._"

# ╔═╡ 54b80a4f-d8bc-40ad-8020-943b3df9578f
bigbreak

# ╔═╡ 83fc7fb2-df70-4959-9dd2-038c6e6df3e9
md"## How to comment"

# ╔═╡ 31ee93da-94ee-41f3-a04a-9be1f3f686ca
# either use the pound/has key!

# ╔═╡ 5364a296-23d4-4c18-8b0e-741501bd3316
#=

For multi-line comments,
use the '#= =#' sequence.

=#

# ╔═╡ 8f0decf1-bdef-456b-bb08-6633876bcd75
bigbreak

# ╔═╡ a3642d38-9c35-4ebc-bc6a-f8d189fffabf
md"## Syntax for basic math"

# ╔═╡ 180d1eba-dd3f-472f-a141-478b0b88f763
begin
	sum = 3 + 7
	difference = 10 - 3
	product = 20 * 5
	quotient = 100 / 10
	power = 10 ^ 2
	modulus = 101 % 2

	(;sum, difference, product, quotient, power, modulus)
end

# ╔═╡ 9ad0c3a5-29f3-40ab-b967-3717f1af4942
bigbreak

# ╔═╡ d1119526-13ae-4555-ab59-cafce95e9ee7
md"## Exercises"

# ╔═╡ 4102dafc-6980-4cbc-a454-41be99a4b64b
md"""
1.1

Look up docs for the convert function.

_Hint: Use `?convert` to access the help._"""

# ╔═╡ 9eee1016-2809-4c1c-8076-1ed172bbbe2c

# ╔═╡ d4ec0db3-5c68-462e-b53c-d495259d6324
md"1.2

Assign 365 to a variable named `days`. Convert days to a float and assign it to variable `days_float` "

# ╔═╡ 4d45400d-b904-418b-80b6-9dc59bda41a6
begin
	days = 0
	days_float = 0
end

# ╔═╡ f6c062f6-b8b7-49d4-ba25-19a5ba26b6e1
begin
	@assert days == 365
	@assert days_float == 365.0
end

# ╔═╡ ada179c2-a434-4947-a85b-e35f340e0d7c
md"""1.3

See what happens when you execute
```julia
convert(Int64, "1")
```
and
```julia
parse(Int64, "1")
```
"""

# ╔═╡ 7ace7d27-1397-4e58-ab44-ce9492e4d1be

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.32"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─cdb3b623-4b24-49a9-b799-ee349287715c
# ╟─855f2a3b-b7ce-4b3d-9ebd-359f8e13bfc2
# ╟─b2e97780-801b-11ec-0433-c7952bc1387d
# ╟─2f8bea58-1b67-4ff0-80b9-70f17c5eb26b
# ╠═feb26c6b-f169-4543-be6e-fd48c4b37f3b
# ╟─f4fbcb5f-8744-45d3-a5a8-acef9efd3a79
# ╠═18c5fa07-8e59-4c25-be5c-90fd7406b20a
# ╟─a9debd68-5e03-4c69-962d-87449024b3e8
# ╠═4df6e513-6b77-4e6b-8870-d9cf84ebdfc7
# ╠═c6dc20f1-d2bd-4f2f-b536-8a1560e928ed
# ╠═71777029-f294-42c6-963c-8b6bc15610bc
# ╠═593e8d96-7a1c-47cc-9e23-7aba58867d5a
# ╠═c66ebcb1-d3f4-4208-ad86-ac3745b8a386
# ╠═3a9a5ec6-4ab8-4c44-a621-76a3c4e08233
# ╟─0b7b3b48-18da-426b-91f9-54f8bac5fd76
# ╠═4ec52ccd-6f3a-4027-87bd-54487ff75562
# ╟─cba0e0e6-531c-4b90-a36f-480b6c5b9d9a
# ╠═f4e86b43-263f-4240-8c10-67987bfd1103
# ╟─ee1e0582-c26e-4550-a3a8-07d4eb1919b8
# ╟─54b80a4f-d8bc-40ad-8020-943b3df9578f
# ╟─83fc7fb2-df70-4959-9dd2-038c6e6df3e9
# ╠═31ee93da-94ee-41f3-a04a-9be1f3f686ca
# ╠═5364a296-23d4-4c18-8b0e-741501bd3316
# ╟─8f0decf1-bdef-456b-bb08-6633876bcd75
# ╟─a3642d38-9c35-4ebc-bc6a-f8d189fffabf
# ╠═180d1eba-dd3f-472f-a141-478b0b88f763
# ╟─9ad0c3a5-29f3-40ab-b967-3717f1af4942
# ╟─d1119526-13ae-4555-ab59-cafce95e9ee7
# ╟─4102dafc-6980-4cbc-a454-41be99a4b64b
# ╠═9eee1016-2809-4c1c-8076-1ed172bbbe2c
# ╟─d4ec0db3-5c68-462e-b53c-d495259d6324
# ╠═4d45400d-b904-418b-80b6-9dc59bda41a6
# ╠═f6c062f6-b8b7-49d4-ba25-19a5ba26b6e1
# ╟─ada179c2-a434-4947-a85b-e35f340e0d7c
# ╠═7ace7d27-1397-4e58-ab44-ce9492e4d1be
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/0 - Getting started.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/0 - Getting started.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/0 - Getting started.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
using vOptGeneric, CPLEX

if length(ARGS) == 0
    m = parseMOP("example.MOP", solver = CplexSolver())
else
    m = parseMOP(ARGS[1], solver = CplexSolver())
end

solve(m, method=:epsilon)

printX_E(m)

using PyPlot
Y_N = getY_N(m)
f1, f2 = map(x -> x[1], Y_N), map(x -> x[2], Y_N)
xlabel("z1") ; ylabel("z2")
plot(f1,f2,"bx", markersize = "6")
!isinteractive() && show()


text:
include("test_2d.jl")
include("test_3d.jl")


text:
using ANOVADDPTest
using Documenter

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

makedocs(;
    modules=[ANOVADDPTest],
    authors="Iván Gutiérrez <ivangutierrez1988@gmail.com> and contributors",
    repo="https://github.com/igutierrezm/ANOVADDPTest.jl/blob/{commit}{path}#{line}",
    sitename="ANOVADDPTest.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://igutierrezm.github.io/ANOVADDPTest.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/igutierrezm/ANOVADDPTest.jl",
)


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

# Parameters
ex = [1.;0.;0.]

h = 1.
r = .05
b1 = Cylinder(r, h, h, color = RGBA(1., 0., 0.))

vert11 = [0.;0.;h / 2]
vert12 = -vert11

# Initial orientation
phi = pi / 4
q1 = UnitQuaternion(RotX(phi))

# Links
N = 1

origin = Origin{Float64}()
links = [Body(b1) for i = 1:N]

# Constraints
jointb1 = EqualityConstraint(Revolute(origin, links[1], ex; p2=vert11))
if N>1
    constraints = [jointb1;[EqualityConstraint(Revolute(links[i - 1], links[i], ex; p1=vert12, p2=vert11)) for i = 2:N]]
else
    constraints = [jointb1]
end

shapes = [b1]

mech = Mechanism(origin, links, constraints;Δt = 0.01, shapes = shapes)
setPosition!(origin,links[1],p2 = vert11,Δq = q1)
previd = links[1].id
for body in Iterators.drop(mech.bodies, 1)
    global previd
    setPosition!(ConstrainedDynamics.getbody(mech, previd), body, p1 = vert12, p2 = vert11)
    previd = body.id
end

storage = simulate!(mech, 10., record = true)
visualize(mech, storage, shapes)


text:
#=
Enumerated type describing how a better patch was found.

!!! JustBetter is better or equal.

Relation between new and previous diff:

>  Worse
= Just Better
< Just Better

new == 0  Perfect  (TODO within epsilon of zero.)
=#

@enum Betterment NotBetter JustBetter PerfectMatch

function classifyBetterment(newDiff, previousDiff)
    if newDiff >= previousDiff
        betterment = NotBetter
    elseif newDiff == 0
        betterment = PerfectMatch
    else
        betterment = JustBetter
    end
    return betterment
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
"
`CrossoverOperator` represents the method that crosses the parents selected to
create a new generation of individuals.

# Fields
- `_method::Function`: method used for crossover.
- `_probability::Float32`: the probability for a set of nParents to being cross.
- `_nParents::Int8`: number of parents used per cross.
- `_nChildren::Int8`: number of children generated by cross.
- `_varArgs::Array{Any}`: aditional arguments for method.

See also: [`setCrossoverOperator`](@ref)
"
struct CrossoverOperator
    _method::Function
    _probability::Float32
    _nParents::UInt8
    _nChildren::UInt8
    _varArgs::Array{Any}
end # struct

"""
    getMethod(crossoverOp::CrossoverOperator)::Function

Obtains the method associated to crossover.
"""
getMethod(crossoverOp::CrossoverOperator)::Function = crossoverOp._method
# function

"""
    getCrossProbability(crossoverOp::CrossoverOperator)::Float32

Obtains the probability for a set of individuals to be cross.
"""
getCrossProbability(crossoverOp::CrossoverOperator)::Float32 = crossoverOp._probability
# function

"""
    getNParents(crossoverOp::CrossoverOperator)::Int8

Obtains the number of parents per cross.
"""
getNParents(crossoverOp::CrossoverOperator)::UInt8 = crossoverOp._nParents
# function

"""
    getNChildren(crossoverOp::CrossoverOperator)::Int8

Obtains the number of children generated by the crossover method.
"""
getNChildren(crossoverOp::CrossoverOperator)::UInt8 = crossoverOp._nChildren
# function

"""
    getFunctionArgs(crossoverOp::CrossoverOperator)::Array

Obtains the aditional arguments associated to crossover method.
"""
getFunctionArgs(crossoverOp::CrossoverOperator)::Array = crossoverOp._varArgs
# function

"""
    getDefaultCrossover(individualType::DataType)::Tuple{Function, Array}

Returns a default crossover method for each registered individual type. For
all `GAGenotype`, it is [`uniformCross`](@ref),
for both `CGPGenotype` and `STGPGenotype` it is [`subtreeCross`](@ref) with a
probability of 0.2; for `GEPGenotype` it is [`kPointRecombinationCross`](@ref)
with one crosspoint.
"""
function getDefaultCrossoverOp(individualType::DataType)::Tuple{Function, Array}
    method = nothing
    varArgs = []
    if individualType <: GAGenotype
        method = uniformCross
    elseif individualType <: CGPGenotype || individualType <: STGPGenotype
        method = subtreeCross
    elseif individualType <: GEPGenotype
        method = kPointRecombinationCross
    end

    return method, varArgs
end # function

"""
    cross_(crossoverOp::CrossoverOperator, selected::Array{Individual},
           experimentInfo::ExperimentInfo)::Array{Individual}

Performs a crossover method between a set of selected individuals. Regardless of
the number of children generated per number of parents necessary to perform the
cross operation, the population will maintain its size.

# Arguments
- `crossoverOp::CrossoverOperator`: structure containing the crossover method, its
    arguments, if any, the number of parents per cross and the number of children
    per cross.
- `selected::Array{Individual}`: the set of selected individuals to cross.
- `experimentInfo::ExperimentInfo`: the information about the experiment currently
    running.

See also: [`cross`](@ref)
"""

function cross_(crossoverOp::CrossoverOperator, selected::Array{Individual},
                experimentInfo::ExperimentInfo)::Array{Individual}

    rng = experimentInfo._rng # Obtaining random number generator

    nParents = crossoverOp._nParents # Number of individual used per cross
    nChildren = crossoverOp._nChildren # Number of individual generated by cross
    nSelected = length(selected) # Number of selected individuals that are going to be cross

    # Calculates the total crossover operations.
    crossIterations = div(nSelected, nParents)
    # Offspring size is the number of children per cross multiplied by the number
    # of crossover operations.
    sizeOffspring = crossIterations * nChildren

    # Reserves memory for the offspring
    offspring = Array{Individual}(undef, sizeOffspring)

    # An array of ranges for offspring, accesing at iteration i of cross and i + 1,
    # we can access the array at the correct positions.
    # For example: suppose nChildren equals 3, the associated array is [0, 3, 6, ..., sizeOffspring].
    # In the first cross, offspring[indexOff[1]+1:indexOff[1+1]] -> offspring[1:3]
    indexOff = collect(0:nChildren:sizeOffspring)
    # An array of ranges for selected parents. See explication for indexOff above.
    indexSelected = collect(0:nParents:nSelected)

    # Obtains an array with the genotypes of selected individuals
    selectedRep = getGenotype(selected)

    # A method is created depending on the genotype, but with the same interface.
    if experimentInfo._individualType <: GPGenotype
        method = (selectedParents) -> crossoverOp._method(selectedParents..., experimentInfo._GPExperimentInfo, rng, crossoverOp._varArgs...)
    else
        method = (selectedParents) -> crossoverOp._method(selectedParents..., rng, crossoverOp._varArgs...)
    end

    # If nChildren are equal to nParents, then the accesing indexes works as it is supposed to.
    if nChildren == nParents
        @inbounds for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i+1]] = selected[indexSelected[i]+1:indexSelected[i+1]]
            end
        end

    # If nChildren are less than nParents, for the case in which there is no
    # cross due to probability, the children will be equal to the nChildren
    # first parents of that cross.
    elseif nChildren < nParents
        @inbounds for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i+1]] = selected[indexSelected[i]+1:indexSelected[i]+nChildren]
            end
        end
    # If nChildren is greater than nParents, an array with the crossed Parents
    # will be carried over through the iterations, for the case in which there is no
    # cross due to probability, we take the parents for that cross as children,
    # and for the rest, we do a sampling among the crossed parents.
    else
        rest =  nChildren-nParents
        selectedForCross = zeros(Bool, len)
        for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                selectedForCross[indexSelected[i]+1:indexSelected[i+1]] = ones(Bool,nParents)
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i]+nParents] = selected[indexSelected[i]+1:indexSelected[i+1]]
                crossedParents = view(selected, selectedForCross)
                if isempty(crossedParents)
                    map!(x -> selected[indexSelected[i]+x], view(offspring, indexOff[i]+nParents+1:indexOff[i+1]), rand(rng, 1:nParents, rest))
                else
                    map!(x -> crossedParents[x], view(offspring, indexOff[i]+nParents+1:indexOff[i+1]), rand(rng, 1:length(crossedParents), rest))
                end
            end
        end
    end

    return offspring
end # function

"""

function cross_(crossoverOp::CrossoverOperator, selected::Array{Individual},
                experimentInfo::ExperimentInfo)::Array{Individual}

    rng = getRNG(experimentInfo) # Obtaining random number generator
    nParents = getNParents(crossoverOp) # Number of individual used per cross
    nChildren = getNChildren(crossoverOp) # Number of individual generated by cross
    nSelected = length(selected) # Number of selected individuals that are going to be cross

    # Calculates the total crossover operations.
    crossIterations = div(nSelected, nParents)
    # Offspring size is the number of children per cross multiplied by the number
    # of crossover operations.
    sizeOffspring = crossIterations * nChildren

    # Reserves memory for the offspring
    offspring = Array{Individual}(undef, sizeOffspring)

    # An array of ranges for offspring, accesing at iteration i of cross and i + 1,
    # we can access the array at the correct positions.
    # For example: suppose nChildren equals 3, the associated array is [0, 3, 6, ..., sizeOffspring].
    # In the first cross, offspring[indexOff[1]+1:indexOff[1+1]] -> offspring[1:3]
    indexOff = collect(0:nChildren:sizeOffspring)
    # An array of ranges for selected parents. See explication for indexOff above.
    indexSelected = collect(0:nParents:nSelected)

    # Obtains an array with the genotypes of selected individuals
    selectedRep = getGenotype(selected)
    # Type of genotype which determines the arguments that cross method receives
    genotypeType = getIndividualType(experimentInfo)

    crossMethod = getMethod(crossoverOp)
    functionArgs = getFunctionArgs(crossoverOp)

    # A method is created depending on the genotype, but with the same interface.
    if genotypeType <: GPGenotype
        method = (selectedParents) -> crossMethod(selectedParents..., experimentInfo._GPExperimentInfo, rng, functionArgs...)
    else
        method = (selectedParents) -> crossMethod(selectedParents..., rng, functionArgs...)
    end

    # If nChildren are equal to nParents, then the accesing indexes works as it is supposed to.
    if nChildren == nParents
        @inbounds for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i+1]] = selected[indexSelected[i]+1:indexSelected[i+1]]
            end
        end

    # If nChildren are less than nParents, for the case in which there is no
    # cross due to probability, the children will be equal to the nChildren
    # first parents of that cross.
    elseif nChildren < nParents
        @inbounds for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i+1]] = selected[indexSelected[i]+1:indexSelected[i]+nChildren]
            end
        end
    # If nChildren is greater than nParents, an array with the crossed Parents
    # will be carried over through the iterations, for the case in which there is no
    # cross due to probability, we take the parents for that cross as children,
    # and for the rest, we do a sampling among the crossed parents.
    else
        rest =  nChildren-nParents
        selectedForCross = zeros(Bool, len)
        for i=1:crossIterations
            if rand(rng) < crossoverOp._probability
                selectedForCross[indexSelected[i]+1:indexSelected[i+1]] = ones(Bool,nParents)
                offspringRep = method(view(selectedRep, indexSelected[i]+1:indexSelected[i+1]))
                map!(x->Individual(x), view(offspring, indexOff[i]+1:indexOff[i+1]), offspringRep)
            else
                offspring[indexOff[i]+1:indexOff[i]+nParents] = selected[indexSelected[i]+1:indexSelected[i+1]]
                crossedParents = view(selected, selectedForCross)
                if isempty(crossedParents)
                    map!(x -> selected[indexSelected[i]+x], view(offspring, indexOff[i]+nParents+1:indexOff[i+1]), rand(rng, 1:nParents, rest))
                else
                    map!(x -> crossedParents[x], view(offspring, indexOff[i]+nParents+1:indexOff[i+1]), rand(rng, 1:length(crossedParents), rest))
                end
            end
        end
    end

    return offspring
end # function
"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# write the triangular elements into arrays rather into objects

include("../../__lib__/std/datastructs/src/datastructs_module.jl")
include("../../__lib__/std/io/json/src/json_module.jl")
include("../src/geom_module.jl")

module meshserial
using ..geom
using ..json
using ..datastructs
# using Threads

function read_mesh(url)
    # file io
    # read the file
    # create a dictionary from it
    # deserialize it

    iostr = open(url, "r");
    iobuff = read(iostr, String);
    close(iostr);

    # dict = json.parse_json(iobuff)
    dict = json.parse(iobuff)
    delmesh = deserialize_mesh(dict)

    return delmesh
end

function write_mesh(delmesh, url)
    # serialize julia object
    # write into file

    serialized_content = ""
    serialized_content = serialize_mesh(delmesh)

    iostr = open(url, "w");
    write(iostr, serialized_content);
    close(iostr);

    return nothing
end

function serialize_mesh(delmesh) # ::DelaunayMesh
    # julia object -> json format

    # mutable struct DelaunayMesh
    #     vertex_array::Vector{<:AbstractPoint}
    #     triangle_array::datastructs.Stack{DelaunayMeshTriangle}

    # "vertex_array": {
    #     "type": Float64,
    #     "length": n,
    #     "x": [],
    #     "y": []
    # },

    # "triangle_array": {
    #     "type": DelaunayMeshTriangle,
    #     "ptr": triangle_array.ptr,
    #     "data": [tri1, tri2,...]
    # }

    # triangle:
    # alive::Bool
    # index_1::Int
    # index_2::Int
    # index_3::Int
    # parent::Int
    # child_1::Int
    # child_2::Int
    # child_3::Int
    # neighbor_12::Int
    # neighbor_23::Int
    # neighbor_31::Int

    str_v = serialize_vertex_array(delmesh)
    str_t = serialize_triangle_array(delmesh)
    str_triidx = serialize_triangle_index(delmesh)

    serialized_content = """
    {
        $(str_v),

        $(str_t),

        $(str_triidx)
    }
    """

    return serialized_content
end

function serialize_triangle_index(delmesh)

    L = size(delmesh.triangle_index,1)

    str_vec_1 = Vector{String}(undef, L)
    str_vec_2 = Vector{String}(undef, L)
    str_vec_3 = Vector{String}(undef, L)

    for kk = 1:L-1
        str_vec_1[kk] = string(delmesh.triangle_index[kk,1]) * ","
        str_vec_2[kk] = string(delmesh.triangle_index[kk,2]) * ","
        str_vec_3[kk] = string(delmesh.triangle_index[kk,3]) * ","
    end

    str_vec_1[L] = string(delmesh.triangle_index[L,1])
    str_vec_2[L] = string(delmesh.triangle_index[L,2])
    str_vec_3[L] = string(delmesh.triangle_index[L,3])

    str_idx1 = reduce(*, str_vec_1)
    str_idx2 = reduce(*, str_vec_2)
    str_idx3 = reduce(*, str_vec_3)

    str = """
    "triangle_index": {
        "length": $(L),
        "idx1": [$(str_idx1)],
        "idx2": [$(str_idx2)],
        "idx3": [$(str_idx3)]
    }
    """
    return str
end

function serialize_vertex_array(delmesh)

    Lp = length(delmesh.vertex_array)

    x_str_vec = Vector{String}(undef, Lp)
    y_str_vec = Vector{String}(undef, Lp)

    delim = ","

    for kk = 1:Lp-1
        v = delmesh.vertex_array[kk]
        x_str_vec[kk] = string(v.x) * delim
        y_str_vec[kk] = string(v.y) * delim
    end
    v = delmesh.vertex_array[Lp]
    x_str_vec[Lp] = string(v.x)
    y_str_vec[Lp] = string(v.y)

    x_str = reduce(*, x_str_vec)
    y_str = reduce(*, y_str_vec)

    str = """
        "vertex_array": {
            "type": "Float64",
            "length": $(length(delmesh.vertex_array)),
            "x": [$(x_str)],
            "y": [$(y_str)]
        }
    """
    return str
end

function serialize_triangle_array(delmesh)
    # "triangle_array": {
    #     "type": DelaunayMeshTriangle,
    #     "ptr": triangle_array.ptr,
    #     "data": [tri1, tri2,...]
    # }

    Lt = length(delmesh.triangle_array)

    alive       = Vector{String}(undef, Lt)
    index_1     = Vector{String}(undef, Lt)
    index_2     = Vector{String}(undef, Lt)
    index_3     = Vector{String}(undef, Lt)
    parent      = Vector{String}(undef, Lt)
    child_1     = Vector{String}(undef, Lt)
    child_2     = Vector{String}(undef, Lt)
    child_3     = Vector{String}(undef, Lt)
    neighbor_12 = Vector{String}(undef, Lt)
    neighbor_23 = Vector{String}(undef, Lt)
    neighbor_31 = Vector{String}(undef, Lt)

    delim = ","
    for kk = 1:(Lt-1)
        tri = delmesh.triangle_array.data[kk]

        alive[kk]       = "$(tri.alive)" * delim
        index_1[kk]     = "$(tri.index_1)" * delim
        index_2[kk]     = "$(tri.index_2)" * delim
        index_3[kk]     = "$(tri.index_3)" * delim
        parent[kk]      = "$(tri.parent)" * delim
        child_1[kk]     = "$(tri.child_1)" * delim
        child_2[kk]     = "$(tri.child_2)" * delim
        child_3[kk]     = "$(tri.child_3)" * delim
        neighbor_12[kk] = "$(tri.neighbor_12)" * delim
        neighbor_23[kk] = "$(tri.neighbor_23)" * delim
        neighbor_31[kk] = "$(tri.neighbor_31)" * delim

    end

    tri = delmesh.triangle_array.data[Lt]

    alive[Lt]       = "$(tri.alive)"
    index_1[Lt]     = "$(tri.index_1)"
    index_2[Lt]     = "$(tri.index_2)"
    index_3[Lt]     = "$(tri.index_3)"
    parent[Lt]      = "$(tri.parent)"
    child_1[Lt]     = "$(tri.child_1)"
    child_2[Lt]     = "$(tri.child_2)"
    child_3[Lt]     = "$(tri.child_3)"
    neighbor_12[Lt] = "$(tri.neighbor_12)"
    neighbor_23[Lt] = "$(tri.neighbor_23)"
    neighbor_31[Lt] = "$(tri.neighbor_31)"

    s_alive       = reduce(*, alive)
    s_index_1     = reduce(*, index_1)
    s_index_2     = reduce(*, index_2)
    s_index_3     = reduce(*, index_3)
    s_parent      = reduce(*, parent)
    s_child_1     = reduce(*, child_1)
    s_child_2     = reduce(*, child_2)
    s_child_3     = reduce(*, child_3)
    s_neighbor_12 = reduce(*, neighbor_12)
    s_neighbor_23 = reduce(*, neighbor_23)
    s_neighbor_31 = reduce(*, neighbor_31)

    str = """
    "triangle_array":
    {
        "type": "DelaunayMeshTriangle",
        "ptr": $(delmesh.triangle_array.ptr),
        "serialized_data": {
            "alive": [$(s_alive)],
            "index_1": [$(s_index_1)],
            "index_2": [$(s_index_2)],
            "index_3": [$(s_index_3)],
            "parent": [$(s_parent)],
            "child_1": [$(s_child_1)],
            "child_2": [$(s_child_2)],
            "child_3": [$(s_child_3)],
            "neighbor_12": [$(s_neighbor_12)],
            "neighbor_23": [$(s_neighbor_23)],
            "neighbor_31": [$(s_neighbor_31)]
        }
    }
    """
    return str
end

#=
function serialize_triangle(tri)
    str = """
    {
        "alive": $(tri.alive),
        "index_1": $(tri.index_1),
        "index_2": $(tri.index_2),
        "index_3": $(tri.index_3),
        "parent": $(tri.parent),
        "child_1": $(tri.child_1),
        "child_2": $(tri.child_2),
        "child_3": $(tri.child_3),
        "neighbor_12": $(tri.neighbor_12),
        "neighbor_23": $(tri.neighbor_23),
        "neighbor_31": $(tri.neighbor_31)
    }
    """
    return str
end
=#

function deserialize_mesh(dict)
    # json dictionary -> julia object

    # "vertex_array": {
    #     "type": Float64,
    #     "length": n,
    #     "x": [],
    #     "y": []
    # },

    # "triangle_array": {
    #     "type": DelaunayMeshTriangle,
    #     "ptr": triangle_array.ptr,
    #     "data": [tri1, tri2,...]
    # }

    # triangle:
    # alive::Bool
    # index_1::Int
    # index_2::Int
    # index_3::Int
    # parent::Int
    # child_1::Int
    # child_2::Int
    # child_3::Int
    # neighbor_12::Int
    # neighbor_23::Int
    # neighbor_31::Int

    point_array    = deserialize_vertex_array(dict)
    triangle_array = deserialize_triangle_array(dict)

    triangle_index = deserialize_triangle_index(dict)

    delmesh                = geom.DelaunayMesh()
    delmesh.vertex_array   = point_array
    delmesh.triangle_array = triangle_array
    delmesh.triangle_index = triangle_index

    return delmesh
end

function deserialize_triangle_array(dict)
    subdict = dict["triangle_array"]
    Lt = Int64(subdict["ptr"])

    data = Vector{geom.DelaunayMeshTriangle}(undef, Lt)
    dict_for_tri_array = subdict["serialized_data"]

    alive       = Vector{Bool}(dict_for_tri_array["alive"])
    index_1     = Vector{Int64}(dict_for_tri_array["index_1"])
    index_2     = Vector{Int64}(dict_for_tri_array["index_2"])
    index_3     = Vector{Int64}(dict_for_tri_array["index_3"])
    parent      = Vector{Int64}(dict_for_tri_array["parent"])
    child_1     = Vector{Int64}(dict_for_tri_array["child_1"])
    child_2     = Vector{Int64}(dict_for_tri_array["child_2"])
    child_3     = Vector{Int64}(dict_for_tri_array["child_3"])
    neighbor_12 = Vector{Int64}(dict_for_tri_array["neighbor_12"])
    neighbor_23 = Vector{Int64}(dict_for_tri_array["neighbor_23"])
    neighbor_31 = Vector{Int64}(dict_for_tri_array["neighbor_31"])

    for kk = 1:Lt

        data[kk] = geom.DelaunayMeshTriangle(
            alive[kk],
            index_1[kk],
            index_2[kk],
            index_3[kk],
            parent[kk],
            child_1[kk],
            child_2[kk],
            child_3[kk],
            neighbor_12[kk],
            neighbor_23[kk],
            neighbor_31[kk])

    end

    triangle_array = datastructs.Stack{geom.DelaunayMeshTriangle}(Lt)
    triangle_array.ptr = Lt
    triangle_array.data = data

    return triangle_array
end

 #=
function deserialize_triangle_element(tri)

    elem = geom.DelaunayMeshTriangle()

    elem.alive = tri["alive"]
    elem.index_1 = tri["index_1"]
    elem.index_2 = tri["index_2"]
    elem.index_3 = tri["index_3"]
    elem.parent = tri["parent"]
    elem.child_1 = tri["child_1"]
    elem.child_2 = tri["child_2"]
    elem.child_3 = tri["child_3"]
    elem.neighbor_12 = tri["neighbor_12"]
    elem.neighbor_23 = tri["neighbor_23"]
    elem.neighbor_31 = tri["neighbor_31"]

    return elem

end
=#

function deserialize_vertex_array(dict)

    subdict = dict["vertex_array"]
    Lp = subdict["length"]
    xs = subdict["x"]
    ys = subdict["y"]

    x = Vector{Float64}(xs)
    y = Vector{Float64}(ys)

    point_array = map( (x, y) -> geom.MutablePoint(x, y), x, y )

    return point_array
end

function deserialize_triangle_index(dict)
    subdict = dict["triangle_index"]

    Lp = subdict["length"]

    i1 = subdict["idx1"]
    i2 = subdict["idx2"]
    i3 = subdict["idx3"]

    idx1 = Vector{Int}(i1)
    idx2 = Vector{Int}(i2)
    idx3 = Vector{Int}(i3)

    triangle_index = [idx1 idx2 idx3]

    return triangle_index
end

end

module stest
using PyPlot
PyPlot.pygui(true)
using ..geom
using ..meshserial

# iostr = open("bla.txt", "w")
# write(iostr, "this is the first line\n")
# write(iostr, "this is the second line\n")
# close(iostr)

# iostr = open("bla.txt", "r")
# iobuff = read(iostr, String)
# close(iostr)

# x_ = [0.0, 1, 2, 3, 2, 1, 0]
# y_ = [0.0, 0, 1, 1, 3, 1, 1]

# point_array = map( (x, y) -> geom.MutablePoint(x, y), x_, y_ )
# delmesh = geom.build_delaunay_mesh!(point_array)

function test_serialize()
    rectangle = geom.Rectangle(-1.0, -1.0, 1.0, 1.0)
    circle = geom.Circle(0.0, 0.0, 0.4)

    distance_fcn = p -> geom.difference_signed_distance_function(
        geom.signed_distance_function(rectangle, p),
        geom.signed_distance_function(circle, p) )

    node_spacing_fcn 	= geom.uniform_node_spacing
    initial_spacing 	= 0.01
    bounding_box 		= geom.Rectangle(-1.0, -1.0, 1.0, 1.0)
    fixed_points 		= zeros(0,2)
    max_iteration 		= 20

    options = geom.MeshGenerationOptions(
        distance_fcn,
        node_spacing_fcn,
        initial_spacing,
        bounding_box,
        fixed_points;
        max_iteration = max_iteration,
        retriangulation_threshold = 0.05)

    @time (point_array, triangle_index, delmesh) = geom.create_mesh(options; visualize = false);

    url = "saved_mesh.txt"

    @time meshserial.write_mesh(delmesh, url)
    # des = meshserial.read_mesh(url)

    # PyPlot.figure()
    # PyPlot.grid()
    # geom.plot_triangulation(delmesh)

    # return nothing
end

function test_read()

    url = "saved_mesh.txt"

    # meshserial.write_mesh(delmesh, url)
    @time des = meshserial.read_mesh(url) # @time

    # PyPlot.figure()
    # PyPlot.grid()
    # geom.plot_triangulation(des)

    return des
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Routines and types to support Helmholtz decomposition

struct ScalarPotentialCache{DVT,VNT,FT} <: AbstractExtraILMCache
   dv :: DVT
   dvn :: VNT
   ftemp :: FT
   divv_temp :: FT
end

struct VectorPotentialCache{DVT,VNT,ST} <: AbstractExtraILMCache
   dv :: DVT
   dvn :: VNT
   stemp :: ST
   curlv_temp :: ST
end

struct VectorFieldCache{ST,VST,FT,VFT,VPT,SPT} <: AbstractExtraILMCache
   ψtemp :: ST
   vψ :: VST
   ϕtemp :: FT
   vϕ :: VFT
   wcache :: VPT
   dcache :: SPT
end

"""
    ScalarPotentialCache(base_cache::BasicILMCache)

Create a cache for calculations involving the scalar potential field
that induces a vector field on the grid. The base cache must be
in support of vector field data.
"""
function ScalarPotentialCache(base_cache::AbstractBasicCache{N,GCT}) where {N,GCT<:Edges{Primal}}
    dv = zeros_surface(base_cache)
    ftemp = zeros_griddiv(base_cache)
    divv_temp = zeros_griddiv(base_cache)
    dvn = ScalarData(dv)
    ScalarPotentialCache(dv,dvn,ftemp,divv_temp)
end

"""
    VectorPotentialCache(base_cache::BasicILMCache)

Create a cache for calculations involving the vector potential field (i.e. streamfunction)
that induces a vector field on the grid. The base cache must be
in support of vector field data.
"""
function VectorPotentialCache(base_cache::BasicILMCache)
    dv = zeros_surface(base_cache)
    stemp = zeros_gridcurl(base_cache)
    curlv_temp = zeros_gridcurl(base_cache)
    dvn = ScalarData(dv)
    VectorPotentialCache(dv,dvn,stemp,curlv_temp)
end

"""
    VectorFieldCache(base_cache::BasicILMCache)

Create a cache for calculations involving the vector field
derived from the vector and scalar potentials.
The base cache must be in support of vector field data.
"""
function VectorFieldCache(base_cache::BasicILMCache)
    ψtemp = zeros_gridcurl(base_cache)
    vψ = zeros_grid(base_cache)
    ϕtemp = zeros_griddiv(base_cache)
    vϕ = zeros_grid(base_cache)
    wfield = VectorPotentialCache(base_cache)
    dfield = ScalarPotentialCache(base_cache)
    VectorFieldCache(ψtemp,vψ,ϕtemp,vϕ,wfield,dfield)
end

"""
    vectorpotential_from_masked_curlv!(ψ::Nodes{Dual},masked_curlv::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache)

Return the vector potential field `ψ` from the masked curl of the vector field, `masked_curlv` (``\\overline{\\nabla\\times\\mathbf{v}}``),
the jump in the vector field across immersed surface `dv`. It solves

``L_n\\psi = -\\overline{\\nabla\\times\\mathbf{v}} - R_n\\mathbf{n}\\times[\\mathbf{v}]``

and returns ``\\psi``.
"""
function vectorpotential_from_masked_curlv!(ψ::Nodes{Dual},curlv::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache)
    @unpack nrm = base_cache
    @unpack stemp = wcache

    fill!(stemp,0.0)
    _curlv_from_surface_jump!(stemp,dv,base_cache)
    stemp .+= curlv

    vectorpotential_from_curlv!(ψ,stemp,base_cache)

    return ψ

end

function _curlv_from_surface_jump!(curlv::Nodes{Dual},dv::VectorData{N},base_cache::BasicILMCache) where {N}
    regularize_normal_cross!(curlv,dv,base_cache)
end

function _curlv_from_surface_jump!(curlv::Nodes{Dual},dv::VectorData{0},base_cache::BasicILMCache)
    return curlv
end

"""
    vectorpotential_from_curlv!(ψ::Nodes{Dual},curlv::Nodes{Dual},base_cache::BasicILMCache)

Return the vector potential field `ψ` from the curl of the masked vector field, `curlv` (``\\nabla\\times\\overline{\\mathbf{v}}``),
the jump in the vector field across immersed surface `dv`. It solves

``L_n\\psi = -\\nabla\\times\\overline{\\mathbf{v}}``

and returns ``\\psi``.
"""
function vectorpotential_from_curlv!(ψ::Nodes{Dual},curlv::Nodes{Dual},base_cache::BasicILMCache)
    inverse_laplacian!(ψ,curlv,base_cache)
    ψ .*= -1.0
end

"""
    vecfield_from_vectorpotential!(v::Edges{Primal},ψ::Nodes{Dual},base_cache::BasicILMCache)

Return the vector field `v` associated with vector potential field `ψ` (in 2-d, a scalar
streamfunction), via the curl

``\\mathbf{v} = \\nabla\\times\\psi \\mathbf{e}_z``
"""
function vecfield_from_vectorpotential!(v::Edges{Primal},ψ::Nodes{Dual},base_cache::BasicILMCache)
    fill!(v,0.0)
    curl!(v,ψ,base_cache)
    return v
end

"""
    masked_curlv_from_curlv_masked!(masked_curlv::Nodes{Dual},curlv::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache)

Return the masked curl of the vector field ``\\overline{\\nabla\\times\\mathbf{v}}`` (`masked_curlv`) from
the curl of the masked vector field ``\\nabla\\times\\overline{\\mathbf{v}}`` (`curlv`). It obtains
this from

``\\overline{\\nabla\\times\\mathbf{v}} = \\nabla\\times\\overline{\\mathbf{v}} - R_n\\mathbf{n}\\times[\\mathbf{v}]``
"""
function masked_curlv_from_curlv_masked!(masked_curlv::Nodes{Dual},curlv::Nodes{Dual},dv::VectorData{N},base_cache::BasicILMCache,wcache::VectorPotentialCache) where {N}
    @unpack nrm = base_cache

    fill!(masked_curlv,0.0)
    #pointwise_cross!(dvn,nrm,dv)
    #regularize!(masked_curlv,dvn,Rn)
    regularize_normal_cross!(masked_curlv,dv,base_cache)
    masked_curlv .*= -1.0
    masked_curlv .+= curlv
    return masked_curlv
end

function masked_curlv_from_curlv_masked!(masked_curlv::Nodes{Dual},curlv::Nodes{Dual},dv::VectorData{0},base_cache::BasicILMCache,wcache::VectorPotentialCache)
    masked_curlv .= curlv
end

"""
    curlv_masked_from_masked_curlv!(curlv::Nodes{Dual},masked_curlv::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache)

Return the curl of the masked vector field ``\\nabla\\times\\overline{\\mathbf{v}}`` (`curlv`) from
the masked curl of the vector field ``\\overline{\\nabla\\times\\mathbf{v}}`` (`masked_curlv`). It obtains
this from

``\\nabla\\times\\overline{\\mathbf{v}} = \\overline{\\nabla\\times\\mathbf{v}} + R_n\\mathbf{n}\\times[\\mathbf{v}]``
"""
function curlv_masked_from_masked_curlv!(curlv::Nodes{Dual},masked_curlv::Nodes{Dual},dv::VectorData{N},base_cache::BasicILMCache,wcache::VectorPotentialCache) where {N}
    @unpack nrm = base_cache
    #@unpack dvn, Rn = wcache

    fill!(curlv,0.0)
    #pointwise_cross!(dvn,nrm,dv)
    #regularize!(curlv,dvn,Rn)
    regularize_normal_cross!(curlv,dv,base_cache)
    curlv .+= masked_curlv
    return curlv
end

function curlv_masked_from_masked_curlv!(curlv::Nodes{Dual},masked_curlv::Nodes{Dual},dv::VectorData{0},base_cache::BasicILMCache,wcache::VectorPotentialCache)
    curlv .= masked_curlv
end

"""
    scalarpotential_from_masked_divv!(ϕ::Nodes{Primal},masked_divv::Nodes{Primal},dv::VectorData,base_cache::BasicILMCache,dcache::ScalarPotentialCache)

Return the scalar potential field `ϕ` from the masked divergence of the vector field, `masked_divv` (``\\overline{\\nabla\\cdot\\mathbf{v}}``)
the jump in the vector field across immersed surface `dv`. It solves

``L_c\\phi = \\overline{\\nabla\\cdot\\mathbf{v}} + R_c\\mathbf{n}\\cdot[\\mathbf{v}]``

and returns ``\\phi``.
"""
function scalarpotential_from_masked_divv!(ϕ::Nodes{Primal},divv::Nodes{Primal},dv::VectorData,base_cache::BasicILMCache,dcache::ScalarPotentialCache)
  @unpack nrm = base_cache
  @unpack dvn, ftemp = dcache

    fill!(ftemp,0.0)
    _divv_from_surface_jump!(ftemp,dv,base_cache)
    ftemp .+= divv

    scalarpotential_from_divv!(ϕ,ftemp,base_cache)
    return ϕ

end

function _divv_from_surface_jump!(divv::Nodes{Primal},dv::VectorData{N},base_cache::BasicILMCache) where N
  regularize_normal_dot!(divv,dv,base_cache)
end

function _divv_from_surface_jump!(divv::Nodes{Primal},dv::VectorData{0},base_cache::BasicILMCache)
    divv
end

"""
    scalarpotential_from_divv!(ϕ::Nodes{Primal},divv::Nodes{Primal},base_cache::BasicILMCache,dcache::ScalarPotentialCache)

Return the scalar potential field `ϕ` from the divergence of the masked vector field, `divv` (``\\nabla\\cdot\\overline{\\mathbf{v}}``)
the jump in the vector field across immersed surface `dv`. It solves

``L_c\\phi = \\nabla\\cdot\\overline{\\mathbf{v}}``

and returns ``\\phi``.
"""
function scalarpotential_from_divv!(ϕ::Nodes{Primal},divv::Nodes{Primal},base_cache::BasicILMCache)

    inverse_laplacian!(ϕ,divv,base_cache)
    return ϕ

end

"""
    vecfield_from_scalarpotential!(v::Edges{Primal},ϕ::Nodes{Primal},base_cache::BasicILMCache)

Return the vector field `v` associated with scalar potential field `ϕ`, via the gradient

``\\mathbf{v} = \\nabla\\phi``
"""
@inline function vecfield_from_scalarpotential!(v::Edges{Primal},ϕ::Nodes{Primal},base_cache::BasicILMCache)
    fill!(v,0.0)
    grad!(v,ϕ,base_cache)
    return v
end

"""
    masked_divv_from_divv_masked!(masked_divv::Nodes{Primal},divv::Nodes{Primal},dv::VectorData,base_cache::BasicILMCache,dcache::ScalarPotentialCache)

Return the masked divergence of the vector field ``\\overline{\\nabla\\cdot\\mathbf{v}}`` (`masked_divv`) from
the divergence of the masked vector field ``\\nabla\\cdot\\overline{\\mathbf{v}}`` (`divv`). It obtains
this from

``\\overline{\\nabla\\cdot\\mathbf{v}} = \\nabla\\cdot\\overline{\\mathbf{v}} - R_c\\mathbf{n}\\cdot[\\mathbf{v}]``
"""
function masked_divv_from_divv_masked!(masked_divv::Nodes{Primal},divv::Nodes{Primal},dv::VectorData{N},base_cache::BasicILMCache,dcache::ScalarPotentialCache) where N
    @unpack nrm = base_cache
    #@unpack dvn, Rc = dcache

    fill!(masked_divv,0.0)
    regularize_normal_dot!(masked_divv,dv,base_cache)
    #pointwise_dot!(dvn,nrm,dv)
    #regularize!(masked_divv,dvn,Rc)
    masked_divv .*= -1.0
    masked_divv .+= divv
    return masked_divv
end

function masked_divv_from_divv_masked!(masked_divv::Nodes{Primal},divv::Nodes{Primal},dv::VectorData{0},base_cache::BasicILMCache,dcache::ScalarPotentialCache)
    masked_divv .= divv
end

"""
    divv_masked_from_masked_divv!(divv::Nodes{Primal},masked_divv::Nodes{Primal},dv::VectorData,base_cache::BasicILMCache,dcache::ScalarPotentialCache)

Return the divergence of the masked vector field ``\\nabla\\cdot\\overline{\\mathbf{v}}`` (`divv`) from the masked divergence
of the vector field ``\\overline{\\nabla\\cdot\\mathbf{v}}`` (`masked_divv`). It obtains
this from

``\\nabla\\cdot\\overline{\\mathbf{v}} = \\overline{\\nabla\\cdot\\mathbf{v}} + R_c\\mathbf{n}\\cdot[\\mathbf{v}]``
"""
function divv_masked_from_masked_divv!(divv::Nodes{Primal},masked_divv::Nodes{Primal},dv::VectorData{N},base_cache::BasicILMCache,dcache::ScalarPotentialCache) where N
    @unpack nrm = base_cache
    #@unpack dvn, Rc = dcache

    fill!(divv,0.0)
    #pointwise_dot!(dvn,nrm,dv)
    #regularize!(divv,dvn,Rc)
    regularize_normal_dot!(divv,dv,base_cache)
    divv .+= masked_divv
    return divv
end

function divv_masked_from_masked_divv!(divv::Nodes{Primal},masked_divv::Nodes{Primal},dv::VectorData{0},base_cache::BasicILMCache,dcache::ScalarPotentialCache)

    divv .= masked_divv
end

"""
    vecfield_helmholtz!(v::Edges{Primal},curlv::Nodes{Dual},divv::Nodes{Primal},dv::VectorData,vp::Union{Edges{Primal},Nothing},base_cache::BasicILMCache,veccache::VectorFieldCache)

Recover the vector field `v` from the masked curl field `curlv` (``\\overline{\\nabla\\times\\mathbf{v}}``),
divergence field `divv` (``\\overline{\\nabla\\cdot\\mathbf{v}}``), surface
vector jump `dv` (``[\\mathbf{v}]``), and additional irrotational, divergence-free vector field `vp`. It obtains this
from the Helmholtz decomposition,

``\\mathbf{v} = \\nabla\\times\\psi + \\nabla\\phi + \\mathbf{v}_p``

in which ``\\psi`` is the solution of

``L_n\\overline{\\psi} = -\\overline{\\nabla\\times\\mathbf{v}} - R_n\\mathbf{n}\\times[\\mathbf{v}]``

in which ``\\phi`` is the solution of

``L_c\\overline{\\phi} = \\overline{\\nabla\\cdot\\mathbf{v}} + R_c\\mathbf{n}\\cdot[\\mathbf{v}]``

It is important to stress that `curlv` is the masked curl of the vector field, not the
curl of the masked vector field. To get this the former from the latter, use [`masked_curlv_from_curlv_masked!`](@ref).
Similarly for `divv`, which is the masked divergence of the vector field, not the divergence
of the masked vector field. One can use [`masked_divv_from_divv_masked!`](@ref) for this.

To specify the irrotational, divergence-free vector field `vp`, one can also simply provide
a tuple, e.g., `(1.0,0.0)`, to specify a uniform vector field.
"""
function vecfield_helmholtz!(v::Edges{Primal},curlv::Nodes{Dual},divv::Nodes{Primal},dv::VectorData,vp::Union{Edges{Primal},Nothing},base_cache::BasicILMCache,veccache::VectorFieldCache)
    @unpack wcache, dcache, ψtemp, vψ, ϕtemp, vϕ = veccache

    fill!(ψtemp,0.0)
    vectorpotential_from_masked_curlv!(ψtemp,curlv,dv,base_cache,wcache)
    vecfield_from_vectorpotential!(vψ,ψtemp,base_cache)

    fill!(ϕtemp,0.0)
    scalarpotential_from_masked_divv!(ϕtemp,divv,dv,base_cache,dcache)
    vecfield_from_scalarpotential!(vϕ,ϕtemp,base_cache)

    v .= vψ .+ vϕ
    if !isnothing(vp)
        v .+= vp
    end
    return v
end

vecfield_helmholtz!(v,curlv,divv,dv,V::Tuple,base_cache,veccache) =
            vecfield_helmholtz!(v,curlv,divv,dv,vecfield_uniformvecfield!(base_cache.gdata_cache,V...,base_cache),base_cache,veccache)

"""
    vectorpotential_uniformvecfield!(ψ::Nodes{Dual},Vx::Real,Vy::Real,base_cache::BasicILMCache)

Return a vector potential field `ψ` associated with a uniform vector field with the components `Vx` and `Vy`.
"""
function vectorpotential_uniformvecfield!(ψ::Nodes{Dual},Vx::Real,Vy::Real,base_cache::BasicILMCache)
    @unpack gcurl_cache = base_cache
    gcurl_cache .= y_gridcurl(base_cache)
    ψ .= Vx*gcurl_cache
    gcurl_cache .= x_gridcurl(base_cache)
    ψ .-= Vy*gcurl_cache
    return ψ
end

"""
    scalarpotential_uniformvecfield!(ϕ::Nodes{Primal},Vx::Real,Vy::Real,base_cache::BasicILMCache)

Return a vector potential field `ϕ` associated with a uniform vector field with the components `Vx` and `Vy`.
"""
function scalarpotential_uniformvecfield!(ϕ::Nodes{Primal},Vx::Real,Vy::Real,base_cache::BasicILMCache)
    @unpack gdiv_cache = base_cache
    gdiv_cache .= x_griddiv(base_cache)
    ϕ .= Vx*gdiv_cache
    gdiv_cache .= y_griddiv(base_cache)
    ϕ .+= Vy*gdiv_cache
    return ϕ
end

"""
    vecfield_uniformvecfield!(v::Edges{Primal},Vx::Real,Vy::Real,base_cache::BasicILMCache)

Return a vector field `v` filled with the uniform components `Vx` and `Vy`, respectively.
"""
function vecfield_uniformvecfield!(v::Edges{Primal},Vx::Real,Vy::Real,base_cache::BasicILMCache)
    v.u .= Vx
    v.v .= Vy
    return v
end


text:
#---------------------------------------------------------------------
#*	only liquid phase
#*--------------------------------------------------------------------
type gibbs_liq
	gibbs_liq()=begin
		NElem=outers.NElem
		new(
			tank_liq(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of elements",
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Universal gas constant",
				:Unit=>"J/mol/K",
				:Default=>8.314
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Number of elements per component"
			)),(NElem,NComp)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Pressure of standard state",
				:Default=>1,
				:DisplayUnit=>"atm"
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Reference temperature",
				:Default=>298.15
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			fill(energy_mol(Dict{Symbol,Any}(
				:Brief=>"Gibbs free-energy change of formation"
			)),(NComp)),
			fill(energy_mol(Dict{Symbol,Any}(
				:Brief=>"Lagrangian multiplier",
				:Symbol=>"\\lambda"
			)),(NElem)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Activity",
				:Symbol=>"\\hat{a}",
				:Lower=>0
			)),(NComp)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Overall component rate of reaction"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fractional conversion of component",
				:Symbol=>"X",
				:Default=>0
			)),(NComp)),
			fill(flow_mol(Dict{Symbol,Any}(
				:Brief=>"Component molar flow rate"
			)),(NComp)),
			[
				:(Outlet.F*Outlet.z = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z + rate*_base_1._base_1.Tank.V),
				:(Outlet.P = _base_1._base_1.Outletm.P),
				:(Outlet.F = sum(Fi)),
				:(Fi = Outlet.F*Outlet.z),
				:(Outlet.F*Outlet.h = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.h),
				:(sumt(Fi*na) = sumt(_base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z*na)),
				:(G = _base_1._base_1.PP.IdealGasGibbsOfFormation(Outlet.T)),
				:(G([1:_base_1._base_1.NComp]) + sumt(lambda*na(:,[1:_base_1._base_1.NComp])) = -Rg*Outlet.T*ln(activ([1:_base_1._base_1.NComp]))),
				:(Fi([1:_base_1._base_1.NComp]) = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z([1:_base_1._base_1.NComp])*(1 - conv([1:_base_1._base_1.NComp]))),
				:(conv([1:_base_1._base_1.NComp]) = 1),
				:(conv([1:_base_1._base_1.NComp]) = 0),
				:(activ = _base_1._base_1.PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z *exp(_base_1._base_1.PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T)),
			],
			[
				"Outlet stream","Mechanical equilibrium","Steady-state","Component molar flow rate","Energy balance","Element balance","Gibbs free-energy of formation","Lagrangian multiplier","Molar conversion","Molar conversion","Molar conversion","Activity",
			],
			[:NElem,:Rg,:na,:Ps,:To,],
			[:Outlet,:G,:lambda,:activ,:rate,:conv,:Fi,]
		)
	end
	_base_1::tank_liq
	NElem::DanaInteger
	Rg::DanaReal
	na::Array{DanaReal}
	Ps::pressure
	To::temperature
	Outlet::liquid_stream
	G::Array{energy_mol}
	lambda::Array{energy_mol}
	activ::Array{DanaReal}
	rate::Array{reaction_mol}
	conv::Array{DanaReal}
	Fi::Array{flow_mol}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export gibbs_liq
function setEquationFlow(in::gibbs_liq)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	#	"Gibbs free-energy of formation without Cp correction"
	#	G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
	#		+ PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
	addEquation(8)
	if (Outletm.z(i) > 1e-16)
		addEquation(9)
	else
		if (Outlet.z(i) > 0)
			addEquation(10)
			# ?

		else
			addEquation(11)
			# ?

		end
	end
	addEquation(12)
end
function atributes(in::gibbs_liq,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/cstr"
	fields[:Brief]="Model of a generic liquid-phase Gibbs CSTR"
	fields[:Info]="
== Assumptions ==
* thermodynamic equilibrium
* steady-state

== Specify ==
* inlet stream
* number of elements related to components
* matrix of elements by components
* equilibrium temperature
"
	drive!(fields,_)
	return fields
end
gibbs_liq(_::Dict{Symbol,Any})=begin
	newModel=gibbs_liq()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module AQFED

__precompile__(true)
module Math
include("math/erfc.jl")
include("math/normal.jl")
include("math/as241.jl")
include("math/lambertw.jl")
include("math/delbourgo_gregory.jl")
end

__precompile__(true)
module Black
include("black/black.jl")
include("black/iv_solver_common.jl")
include("black/iv_solver_halley.jl")
include("black/iv_solver_householder.jl")
include("black/iv_solver_lisor.jl")
include("black/iv_solver_jaeckel.jl")

impliedVolatility(isCall::Bool, price::T, f::T, strike::T, tte::T,  df::T)  where {T} =  impliedVolatilitySRHalley(isCall, price,f,strike,tte, df, T(0.0), 64, Householder())
export impliedVolatility
end

module Bachelier
include("bachelier/bachelier.jl")
include("bachelier/iv_solver_lfk.jl")
end

__precompile__(true)
module Random
include("random/mixmax.jl")
include("random/well1024a.jl")
include("random/mt64.jl")
include("random/blabla.jl")
include("random/philox.jl")
include("random/chacha.jl")
include("random/chacha_simd.jl")
include("random/mrg63k3a.jl")
include("random/ssobol.jl")
include("random/dsobol.jl")
include("random/rngseq.jl")
end

module TermStructure
include("ts/heston.jl")
include("ts/lv.jl")
include("ts/surface.jl")
include("ts/dividend.jl")
end

module American
include("american/andersenlake.jl")
include("american/andersenlaken.jl")
end

module Basket
include("basket/deelstra.jl")
include("basket/asian.jl")
include("basket/pln.jl")
end

module PLN
include("pln/etore_gobet.jl")
end

module MonteCarlo
include("mc/bb.jl")
include("mc/bbvector.jl")
include("mc/payoffs.jl")
include("mc/black.jl")
include("mc/heston.jl")
include("mc/lv.jl")
end

module Collocation
include("collocation/convex.jl") #utilities to ensure convexity of prices, and build X for collocation
include("collocation/poly.jl")
include("collocation/logpoly.jl")
include("collocation/quadratic_spline.jl")
include("collocation/expbspline.jl")
end
end