text:
# This file was generated, do not modify it.

using MLJ, RDatasets, ScientificTypes, PrettyPrinting

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

y, X = unpack(hitters, ==(:Salary), col->true);

no_miss = .!ismissing.(y)
y = collect(skipmissing(y))
X = X[no_miss, :]
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=424);

Xc = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))
scitype(Xc)

@pipeline RegPipe(std = Standardizer(),
                  hot = OneHotEncoder(),
                  reg = LinearRegressor())

model = RegPipe()
pipe  = machine(model, Xc, y)
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

pipe.model.reg = RidgeRegressor()
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

r  = range(model, :(reg.lambda), lower=1e-2, upper=100_000, scale=:log10)
tm = TunedModel(model=model, ranges=r, tuning=Grid(resolution=50),
                resampling=CV(nfolds=3, rng=4141), measure=rms)
mtm = machine(tm, Xc, y)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.reg.lambda, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

mtm.model.model.reg = LassoRegressor()
mtm.model.ranges = range(model, :(reg.lambda), lower=500, upper=100_000, scale=:log10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.reg.lambda, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

coefs = mtm.fitresult.fitresult.machine.fitresult
round.(coefs, sigdigits=2)

sum(coefs .≈ 0) / length(coefs)

@load ElasticNetRegressor pkg=MLJLinearModels

mtm.model.model.reg = ElasticNetRegressor()
mtm.model.ranges = [range(model, :(reg.lambda), lower=0.1, upper=100, scale=:log10),
                    range(model, :(reg.gamma),  lower=500, upper=10_000, scale=:log10)]
mtm.model.tuning = Grid(resolution=10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
@show round(best_mdl.reg.lambda, sigdigits=4)
@show round(best_mdl.reg.gamma, sigdigits=4)

ŷ = predict(mtm, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using Pkg

Pkg.develop(path = "..")

using Revise
using Publish
using Artifacts, LazyArtifacts

using Metalhead

# override default theme
cp(artifact"flux-theme", "../_flux-theme"; force = true)

p = Publish.Project(Metalhead)

# serve documentation
serve(Metalhead)

ext: jl


text:
using BinaryBuilder

name = "Sundials32"
version = v"5.3.0"

# Collection of sources required to build Sundials
sources = [
    GitSource("https://github.com/LLNL/sundials.git",
              "887af4374af2271db9310d31eaa9b5aeff49e829"),
    DirectorySource("../bundled@5"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sundials*

# Set up CFLAGS
if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 $WORKSPACE/srcdir/patches/Sundials_windows.patch
    # Work around https://github.com/LLNL/sundials/issues/29
    export CFLAGS="${CFLAGS} -DBUILD_SUNDIALS_LIBRARY"
elif [[ "${target}" == powerpc64le-* ]]; then
    export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
fi

# Set up LAPACK
LAPACK_LIBRARIES="-lgfortran ${libdir}/libopenblas.${dlext}"
if [[ "${target}" == i686-* ]] || [[ "${target}" == x86_64-* ]]; then
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} -lquadmath"
elif [[ "${target}" == powerpc64le-* ]]; then
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} -lgomp -ldl -lm -lpthread"
fi

# Fix the SuperLU_MT library name
atomic_patch -p1 $WORKSPACE/srcdir/patches/Sundials_SuperLU_MT.patch

# Use GCC on Apple/FreeBSD
toolchain="$CMAKE_TARGET_TOOLCHAIN"
if [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-freebsd* ]]; then
    toolchain="${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake"
fi

# Set the mangling scheme manually on Apple
if [[ "${target}" == *-apple-* ]]; then
    mangling="-DSUNDIALS_F77_FUNC_CASE=lower -DSUNDIALS_F77_FUNC_UNDERSCORES=one"
fi

# Build
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_TOOLCHAIN_FILE="$toolchain" \
    -DEXAMPLES_ENABLE_C=OFF \
    -DKLU_ENABLE=ON \
    -DKLU_INCLUDE_DIR="${includedir}" \
    -DKLU_LIBRARY_DIR="$libdir" \
    -DLAPACK_ENABLE=ON \
    -DLAPACK_LIBRARIES:STRING="${LAPACK_LIBRARIES}" \
    -DSUPERLUMT_ENABLE=ON \
    -DSUPERLUMT_INCLUDE_DIR="${includedir}" \
    -DSUPERLUMT_LIBRARY_DIR="$libdir" \
    -DSUPERLUMT_LIBRARIES="${libdir}/libopenblas.${dlext}" \
    -DSUPERLUMT_THREAD_TYPE="OpenMP" \
    -DSUNDIALS_INDEX_SIZE=32 \
    -DBUILD_STATIC_LIBS=OFF \
    ${mangling} \
    ..
make -j${nproc}
make install

# Move libraries to ${libdir} on Windows
if [[ "${target}" == *-mingw* ]]; then
    mv ${prefix}/lib/libsundials_*.${dlext} "${libdir}"
fi
"""

# We attempt to build for all defined platforms
platforms = supported_platforms()
platforms = expand_gfortran_versions(platforms)
platforms = [p for p in platforms if !(arch(p) == :powerpc64le)]

products = [
    LibraryProduct("libsundials_arkode", :libsundials_arkode),
    LibraryProduct("libsundials_cvode", :libsundials_cvode),
    LibraryProduct("libsundials_cvodes", :libsundials_cvodes),
    LibraryProduct("libsundials_ida", :libsundials_ida),
    LibraryProduct("libsundials_idas", :libsundials_idas),
    LibraryProduct("libsundials_kinsol", :libsundials_kinsol),
    LibraryProduct("libsundials_nvecmanyvector", :libsundials_nvecmanyvector),
    LibraryProduct("libsundials_nvecserial", :libsundials_nvecserial),
    LibraryProduct("libsundials_sunlinsolband", :libsundials_sunlinsolband),
    LibraryProduct("libsundials_sunlinsoldense", :libsundials_sunlinsoldense),
    LibraryProduct("libsundials_sunlinsolklu", :libsundials_sunlinsolklu),
    LibraryProduct("libsundials_sunlinsollapackband", :libsundials_sunlinsollapackband),
    LibraryProduct("libsundials_sunlinsollapackdense", :libsundials_sunlinsollapackdense),
    LibraryProduct("libsundials_sunlinsolpcg", :libsundials_sunlinsolpcg),
    LibraryProduct("libsundials_sunlinsolspbcgs", :libsundials_sunlinsolspbcgs),
    LibraryProduct("libsundials_sunlinsolspfgmr", :libsundials_sunlinsolspfgmr),
    LibraryProduct("libsundials_sunlinsolspgmr", :libsundials_sunlinsolspgmr),
    LibraryProduct("libsundials_sunlinsolsptfqmr", :libsundials_sunlinsolsptfqmr),
    LibraryProduct("libsundials_sunmatrixband", :libsundials_sunmatrixband),
    LibraryProduct("libsundials_sunmatrixdense", :libsundials_sunmatrixdense),
    LibraryProduct("libsundials_sunmatrixsparse", :libsundials_sunmatrixsparse),
    LibraryProduct("libsundials_sunnonlinsolfixedpoint", :libsundials_sunnonlinsolfixedpoint),
    LibraryProduct("libsundials_sunnonlinsolnewton", :libsundials_sunnonlinsolnewton),
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
    BuildDependency("SuiteSparse32_jll"),
    Dependency("SuperLU_MT_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# filetypes.jl: contains methods needed to parse different file types

"""
    function readXYZ(
        io::IO;
        basis::AbstractVector{AbstractVector{<:Real}},
        spgrp::Integer
        origin::AbstractVector{<:Real}
        prim = false,
        conv = true,
    )

Reads in an XYZ file supplied by the user.

If `warn` is set to `true`, a warning message will be printed to `err` (defaults to `stderr`)
noting that cell vector data is not present in an XYZ file.
"""
function readXYZ(
    io::IO;
    basis::AbstractVector{AbstractVector{<:Real}},
    spgrp::Integer = 0,
    origin::AbstractVector{<:Real} = [0, 0, 0],
    prim = false,
    conv = true,
)
    # Read the file in character by character
    local sites::Vector{SVector{3,Float64}}
    local atoms::Vector{String}
    for (n, ln) in enumerate(eachline(io))
        if n > 2
            # Data starts after the 2nd line
            data = split(ln)
            atoms[n - 2] = data[1]
            sites[n - 2] = SVector{3, Float64}(parse.(Float64, data[2:4]))
        elseif n == 1
            # This is where you get the number of atoms
            natoms = parse(Int, ln)
            atoms = ["" for m = 1:natoms]
            sites = [SVector{3, Float64}(0, 0, 0) for m = 1:natoms]
        end
    end
    # Convert to fractional coordinates if a basis is supplied
    if basis != _basis{3}()
        sites = [inv(hcat(basis...)) * v for v in sites]
    end
    return CrystalStructure{3}(
        basis = basis,
        spgrp = spgrp,
        origin = origin,
        atoms = atoms,
        sites = sites,
        prim = prim,
        conv = conv
    )
end

"""
    readXYZ(filename::AbstractString; warn=false, err=stderr)

Reads in an XYZ file supplied by the user.

If `warn` is set to `true`, a warning message will be printed to `err` (defaults to `stderr`)
noting that cell vector data is not present in an XYZ file.
"""
function readXYZ(
    filename::AbstractString;
    basis::AbstractVector{AbstractVector{<:Real}},
    spgrp::Integer = 0,
    origin::AbstractVector{<:Real} = [0, 0, 0],
    prim = false,
    conv = true,
)
    open(filename) do io
        readXYZ(io, basis=basis, spgrp=spgrp, origin = origin, prim = prim, conv = conv)
    end
end

"""
    writeXYZ(xtal::AbstractCrystal{3}, io::IO)

Writes 3D crystal data to an XYZ file.
"""
function writeXYZ(xtal::AbstractCrystal{3}, io::IO)
    # Write the number of entries
    sz = length(xtal.atoms)
    println(io, string(sz))
    # Write the comment line
    println(io, "Generated by CrystalStructures.jl")
    # Convert basis vectors to a matrix
    hasbasis = (basis == zeros(SMatrix{3,3,Float64}))
    for n in 1:sz
        # Convert the sites from reduced to Cartesian coordinates
        site =  hasbasis ? hcat(xtal.basis...) * xtal.sites[n] : xtal.sites[n]
        @printf(io, "%s\t%f\t%f\t%f\n", xtal.atoms[n], site...)
    end
    return nothing
end

writeXYZ(xtaldata::AbstractCrystalData{3}, io::IO) = writeXYZ(xtaldata.xtal, io)

function writeXYZ(xtal::AbstractCrystal{3}, filename::AbstractString)
    open(filename) do io
        writeXYZ(xtal, io)
    end
end

writeXYZ(xtaldata::AbstractCrystalData{3}, filename::AbstractString) =
    writeXYZ(xtaldata.xtal, filename)

function writeXYZ(xtaldata::CrystalStructureWithData{3}, filename::AbstractString)
    return writeXYZ(xtaldata.xtal, filename)
end

"""
    readXSF(
        io::IO;
        spgrp::Int = 0,
        origin::Vector{Float64} = [0, 0, 0],
        tryconv=false
    )

Reads in an XCrysDen XSF file from an I/O stream and returns a `CrystalStructure{N}`. Since XSF
data does not contain symmetry information, space group data is not supplied by default, but it may
be included manually by setting the `spgrp` parameter.

This function was written to work around some bugs in an in-house program, `bin2xsf`, that
incorrectly placed the `END_DATAGRID_3D` token on the same line as data. This workaround should not
affect the reading of XSF files with correct syntax.

This function may return an `XSFdata` of arbitrary dimensionality and therefore is not type-stable.
"""
function readXSF(
    io::IO;
    spgrp::Int = 0,
    origin::AbstractVector{<:Number} = [0, 0, 0],
    tryconv::Bool=false
)
    # Iterator that loops through each line
    iter = eachline(io)
    # Variables needed
    D = 3
    primvec = Vector{Vector{Float64}}()
    convvec = Vector{Vector{Float64}}()
    atomnums = Vector{Int}()
    primcoord = Vector{Vector{Float64}}()
    convcoord = Vector{Vector{Float64}}()
    datagrid_2D = Dict{String, DataGrid{2, Float64}}()
    datagrid_3D = Dict{String, DataGrid{3, Float64}}()
    # Using this while loop because we want to advance the iterator in subloops
    while !isnothing(iter)
        # Get the current line
        # TODO: when this becomes `nothing` it should break, otherwise the assignment fails
        ln = let i = iterate(iter); isnothing(i) ? break : i[1]; end
        # Check for empty lines or comment lines; advance if present
        # Comments are only allowed between blocks in a valid XSF file
        # No need to check within the loops below
        (isempty(ln) || startswith(ln, "#") || all(isspace, ln)) && continue
        # Check for a keyword
        # TODO: This could likely be refactored, perhaps with subfunctions?
        if contains(ln, "DIM-GROUP")
            # Move to the next line
            ln = iterate(iter)[1]
            # Get the dimensionality
            D = parse(Int, split(ln)[1])
        elseif contains(ln, "PRIMVEC")
            # Store vectors here
            primvec = Vector{Vector{Float64}}(undef, D)
            for n in 1:D
                # Move to the next line
                ln = iterate(iter)[1]
                # Get the vector and store it
                primvec[n] = parse.(Float64, split(ln))
            end
        elseif contains(ln, "CONVVEC")
            # Store vectors here
            convvec = Vector{Vector{Float64}}(undef, D)
            for n in 1:D
                # Move to the next line
                ln = iterate(iter)[1]
                # Get the vector and store it
                convvec[n] = parse.(Float64, split(ln))
            end
        elseif contains(ln, "PRIMCOORD")
            # Move to the next line
            ln = iterate(iter)[1]
            # Get number of atoms
            natoms = parse(Int, split(ln)[1])
            # Store atom numbers and vectors here
            atomnums = Vector{Int}(undef, natoms)
            primcoord = Vector{Vector{Float64}}(undef, natoms)
            # Determine the number of atoms
            for n in 1:natoms
                ln = iterate(iter)[1]
                # Every line will be of the following format:
                # [atomname] [x] [y] [z]
                data = split(ln)
                atomnums[n] = parse(Int, data[1])
                primcoord[n] = parse.(Float64, data[2:end])
            end
            # TODO: convert the coordinates in the primitive basis
            if !isempty(primvec)
                primcoord = [inv(hcat(primvec...))*v for v in primcoord]
            end
        elseif contains(ln, "CONVCOORD") && isempty(primcoord)
            # Only do this if primcoord was never given
            # Move to the next line
            ln = iterate(iter)[1]
            # Get number of atoms
            natoms = parse(Int, split(ln)[1])
            # Store atom numbers and vectors here
            atomnums = Vector{Int}(undef, natoms)
            convcoord = Vector{Vector{Float64}}(undef, natoms)
            # Determine the number of atoms
            for n in 1:natoms
                ln = iterate(iter)[1]
                # same format as primcoord
                data = split(ln)
                atomnums[n] = parse(Int, data[1])
                convcoord[n] = parse.(Float64, data[2:end])
            end
        # TODO: implement ATOMS and DATAGRID_2D
        elseif contains.(ln, "BEGIN_BLOCK_DATAGRID_3D") || contains(ln, "BEGIN_BLOCK_DATAGRID3D")
            # Need to use both to avoid issues in bin2xsf
            # Skip the comment line that comes right after
            iterate(iter)
            # Loop until end of block
            while !contains(ln, "END_BLOCK_DATAGRID")
                ln = iterate(iter)[1]
                if contains(ln, "DATAGRID_3D_")
                    # Get the name of the datagrid
                    name = rstrip(ln[findfirst("3D_", ln).stop + 1:end])
                    # Get the dimensionality
                    ln = iterate(iter)[1]
                    dims = parse.(Int, split(ln)[1:3])
                    # Form datagrid as a vector and reshape later
                    datavec = Vector{Float64}(undef, prod(dims))
                    # Get the origin, but I'm not sure if I'll ever need this.
                    ln = iterate(iter)[1]
                    origin = parse.(Float64, split(ln)[1:3])
                    # Check the vectors match primitive or conventional cell
                    vec = Vector{Vector{Float64}}(undef, 3)
                    for n = 1:3
                        ln = iterate(iter)[1]
                        vec[n] = parse.(Float64, split(ln))
                    end
                    if vec != primvec && vec != convvec
                        @warn "Datagrid \"" * name *
                              "\" does not use either primitive or conventional cell vectors."
                    end
                    # Get the actual data
                    entry = 0
                    # Loop until we're either out of entries or we hit an end block
                    while entry < prod(dims) || !contains(ln, "END_DATAGRID")
                        ln = iterate(iter)[1]
                        # This code should tolerate weird spacing of numbers
                        # and work around the bin2xsf bug
                        data = split(ln)
                        # Try to parse and remove unparseable data
                        vals = filter(!isnothing, tryparse.(Float64, data))
                        sz = length(vals)
                        # Add the new data to the datagrid
                        datavec[entry .+ (1:sz)] = vals
                        # Shift the entry counter over
                        entry = entry + sz
                    end
                    # Reshape the datagrid
                    # TODO: probably need logic to determine when to restrict this
                    # But for now, we assume that all XSF files contain
                    grid = reshape(datavec, dims[1], dims[2], dims[3])[1:end-1, 1:end-1, 1:end-1]
                    # Add it to the dictionary
                    get!(datagrid_3D, name, DataGrid{3,Float64}(vec, grid))
                end
            end
        end
    end
    # Convert the atom numbers to names
    atoms = ELEMENTS[atomnums]
    # Find the right basis
    if isempty(convvec) || !tryconv
        # No conventional vectors found, or specifically requested not to use them
        if isempty(primvec)
            # No cell vectors found
            basis = _basis{D}()
        else
            # Use primitive vectors
            basis = primvec
        end
    else
        basis = convvec
    end
    # Return the CrystalStructure{D}
    return CrystalStructureWithData{D,Float64}(
        CrystalStructure{D}(
            basis = basis,
            spgrp = spgrp,
            # Ignore supplied origin if not 3D
            origin = (D == 3 ? origin : zeros(Float64, D)),
            atoms = atoms,
            sites = (isempty(primcoord) ? convcoord : primcoord),
            conv = basis == convvec,
            prim = basis == primvec,
        ),
        datagrid_3D
    )
end

"""
    readXSF(
        filename::AbstractString;
        spgrp::Int = 0,
        origin::Vector{Float64} = [0, 0, 0],
        tryconv=true
    )

Reads in an XCrysDen XSF file specified by the user and returns a `CrystalStructure{N}`. Since XSF
data does not contain symmetry information, space group data is not supplied by default, but it may
be included manually by setting the `spgrp` parameter.

This function was written to work around some bugs in an in-house program, `bin2xsf`, that
incorrectly placed the `END_DATAGRID_3D` token on the same line as data. This workaround should not
affect the reading of XSF files with correct syntax.

This function may return an `XSFdata` of arbitrary dimensionality and therefore is not type-stable.
"""
function readXSF(
    filename::AbstractString;
    spgrp::Int = 0,
    origin::AbstractVector{<:Number} = [0, 0, 0],
    tryconv::Bool=false
)
    return open(filename) do io
        readXSF(io, spgrp=spgrp, origin=origin, tryconv=tryconv)
    end
end

"""
    readXSF3D(
        io::IO;
        spgrp::Int = 0,
        origin::Vector{Float64} = [0, 0, 0],
        tryconv=false
    ) -> CrystalStructure{3}

Reads in an XCrysDen XSF file from an I/O stream and returns a `CrystalStructure{N}`. Since XSF
data does not contain symmetry information, space group data is not supplied by default, but it may
be included manually by setting the `spgrp` parameter.

This function was written to work around some bugs in an in-house program, `bin2xsf`, that
incorrectly placed the `END_DATAGRID_3D` token on the same line as data. This workaround should not
affect the reading of XSF files with correct syntax.

This function is type-stable, but will only work on XSF files representing 3D crystal data (data of
other dimensionality and molecular data are not supported).
"""
function readXSF3D(
    io::IO;
    spgrp::Int = 0,
    origin::AbstractVector{<:Number} = [0, 0, 0],
    tryconv=false
)
    # TODO: probably have some better logic for advancing the iterator
    # Used to get basis vectors
    function getbasis!(itr)
        v = Vector{Vector{Float64}}(undef, 3)
        for n in 1:3
            # Advance to the next line
            ln = iterate(itr)[1]
            v[n] = [parse(Float64, s) for s in split(ln)[1:3]]
        end
        return v
    end

    # Used to get atomic positions
    function getsites!(itr, natoms)
        a = Vector{String}(undef, natoms)
        v = Vector{Vector{Float64}}(undef, natoms)
        for n in 1:natoms
            # Advance to the next line
            ln = iterate(itr)[1]
            # Split the data to get the atom and 3 vector components
            data = split(ln)
            # Conversion to a strong occurs in here!
            a[n] = ELEMENTS[parse(Int, data[1])]
            v[n] = [parse(Float64, s) for s in data[2:4]]
        end
        return (a,v)
    end

    # Used to get datagrids
    function getgrid!(itr, sz; trim = true)
        # Create a new vector with length the same size as the datagrid
        v = Vector{Float64}(undef, prod(sz))
        entry = 0
        ln = ""
        while entry < prod(sz) || !contains(ln, "END_DATAGRID")
            ln = iterate(itr)[1]
            dt = filter(!isnothing, [tryparse(Float64, s) for s in split(ln)])
            v[entry .+ (1:length(dt))] = dt
            # Move the entry counter over
            entry = entry + length(dt)
        end
        # Reshape and trim the edges of the grid, if needed
        return reshape(v, sz[1], sz[2], sz[3])[1:end - trim, 1:end - trim, 1:end - trim]
    end

    # Iterator for the file, reads line by line
    iter = eachline(io)
    # Data needed to construct the
    primvec = Vector{Vector{Float64}}(undef, 3)
    convvec = Vector{Vector{Float64}}(undef, 3)
    primcoord = Vector{Vector{Float64}}()
    convcoord = Vector{Vector{Float64}}()
    primatoms = Vector{String}()
    convatoms = Vector{String}()
    datagrids = Dict{String,DataGrid{3,Float64}}()
    # The loop
    while !isnothing(iter)
        # Set ln to the next line - break if there is no other line
        ln = let i = iterate(iter); isnothing(i) ? break : i[1]; end
        # Check for empty lines or comment lines; advance if present
        # Comments are only allowed between blocks in a valid XSF file
        # No need to check within the loops below
        (isempty(ln) || startswith(ln, "#") || all(isspace, ln)) && continue
        # Check for keywords
        if contains(ln, "BEGIN_BLOCK_DATAGRID_3D") || contains(ln, "BEGIN_BLOCK_DATAGRID3D")
            # Skip comment line
            iterate(iter)
            # Loop through the whole datagrid block
            while !contains(ln, "END_BLOCK_DATAGRID")
                ln = iterate(iter)[1]
                # Check for a DATAGRID_3D_ keyword
                if contains(ln, "DATAGRID_3D_")
                    # Get the name of the datagrid - this is complicated
                    name = rstrip(ln[findfirst("3D_", ln).stop + 1:end])
                    # Get the grid size
                    ln = iterate(iter)[1]
                    sz = [parse(Int, s) for s in split(ln)]
                    # Get the origin
                    ln = iterate(iter)[1]
                    origin = [parse(Float64, s) for s in split(ln)]
                    # Get the basis vectors for this dataset
                    basis = getbasis!(iter)
                    # Get the actual data
                    grid = getgrid!(iter, sz, trim=true)
                    # Store in the dictionary
                    get!(datagrids, name, DataGrid{3,Float64}(basis, origin, grid))
                end
            end
        elseif contains(ln, "PRIMVEC")
            primvec = getbasis!(iter)
        elseif contains(ln, "CONVVEC")
            convvec = getbasis!(iter)
        elseif contains(ln, "PRIMCOORD")
            ln = iterate(iter)[1]
            natoms = parse(Int, split(ln)[1])
            (primatoms, primcoord) = getsites!(iter, natoms)
        elseif contains(ln, "CONVCOORD")
            ln = iterate(iter)[1]
            natoms = parse(Int, split(ln)[1])
            (convatoms, convcoord) = getsites!(iter, natoms)
        end
    end
    # Choose which basis vectors to use
    # Only try conventional cell data if there are conventional cell vectors/site data
    tryconv = (tryconv && !isempty(convvec) && !isempty(convcoord))
    # Choose which basis vectors to use
    basis = (tryconv ? convvec : primvec)
    # Construct the CrystalStructureWithData{3,Float64}
    return CrystalStructureWithData{3,Float64}(
        CrystalStructure{3}(
            basis = basis,
            spgrp = spgrp,
            origin = origin,
            atoms = tryconv ? convatoms : primatoms,
            sites = [inv(hcat(basis...)) * v for v in (tryconv ? convcoord : primcoord)],
            conv = tryconv,
            prim = (primcoord == convcoord) || !tryconv
        ),
        datagrids
    )
end

function readXSF3D(
    filename::AbstractString;
    spgrp::Int = 0,
    origin::AbstractVector{<:Number} = [0, 0, 0],
    tryconv=false
)
    return open(filename) do io
        readXSF3D(io, spgrp = spgrp, origin = origin, tryconv = tryconv)
    end
end

"""
    function writeXSF3D(
        xtal::CrystalStructureWithData{3},
        io::IO;
        linesize = 4
    )

Writes an XSF file for a 3D crystal structure with data.
"""
function writeXSF3D(
    xtaldata::CrystalStructureWithData{3, T},
    io::IO;
    linesize = 4
) where T<:Real
    # Needed to format the numbers
    # Matches the format in our in-house bin2xsf
    tostr(x::Number) = rpad(lpad(@sprintf("%.14f", x), 20), 22)
    # Start writing the file
    println(io, "# Written by CrystalStructures.jl\nDIM-GROUP\n 3 1\nPRIMVEC")
    # Write in the primitive cell vectors
    if xtaldata.xtal.prim && !isempty(xtaldata.xtal.basis)
        # Write the primitive cell vectors
        for b in xtaldata.xtal.basis
            println(io, join(tostr.(b)))
        end
    else
        error("Writing crystal structures without primitive vectors is unsupported.")
    end
    # Write PRIMCOORD; skip atoms keyword
    println(io, "PRIMCOORD\n" * lpad(length(xtaldata.xtal.sites), 12) * "  1")
    sites_cart = [hcat(xtaldata.xtal.basis...) * v for v in xtaldata.xtal.sites]
    # TODO: try to figure out what happens when atom names are nonstandard
    atomnums = [ELEMENT_LOOKUP[s] for s in xtaldata.xtal.atoms]
    for (n, v) in zip(atomnums, sites_cart)
        println(io, lpad(string(n), 9) * "  " * join(tostr.(v)))
    end
    # Start writing in the 3D datasets
    println(io, "BEGIN_BLOCK_DATAGRID_3D\nWritten by CrystalStructures.jl")
    for (str,grid) in xtaldata.data
        # Print the datagrid name
        println(io, "    DATAGRID_3D_" * str)
        # Print the datagrid size
        # Note: dimensions of the grid must be increased by 1
        println(io, join([lpad(string(n+1), 6) for n in size(grid.data)]))
        # Print the dataset origin position
        println(io, join([rpad(lpad(@sprintf("%f", x), 20), 22) for x in xtaldata.xtal.origin]))
        # Print the cell vector information for the datagrid
        for b in grid.basis
            println(io, join(tostr.(b)))
        end
        # Print the datagrid
        sz = size(grid.data) .+ 1
        # Counter for inserting newlines
        ctr = 0
        for z = 1:sz[3], y = 1:sz[2], x = 1:sz[1]
            # Generate the new coordinate
            coord = ((x, y, z) .- 1) .% (sz .- 1) .+ 1
            # Write the data at that coordinate
            print(io, tostr(grid.data[coord...]))
            # Insert \n when needed
            ctr = ctr + 1
            ctr % linesize == 0 && print(io, "\n")
        end
        # Print the end of the datagrid
        println(io, "\n    END_DATAGRID_3D")
    end
    println(io, "END_BLOCK_DATAGRID_3D")
    return nothing
end

function writeXSF3D(
    xtaldata::CrystalStructureWithData{3},
    filename::AbstractString;
    linesize = 4
)
    open(filename, write=true) do io
        writeXSF3D(xtaldata, io, linesize=linesize)
    end
    return nothing
end

# TODO: finish methods for dealing with VASP files
"""
    readPOSCAR(io::IO) -> CrystalStructure{3}

Reads a VASP POSCAR input file, which contains information a crystal's unit cell and atoms. Space
group and origin information can be provided as parameters.
"""
function readPOSCAR(
    io::IO;
    spgrp::Integer = 0,
    origin::AbstractVector{<:Number} = [0, 0, 0],
    prim = true,        # To be expected
    conv = false,
)
    iter = eachline(io)
    # Skip the first line (comment)
    # Second line is a scaling factor
    ln = iterate(iterate(iter))[1]
    scalar = parse(Float64, strip(ln))
    # Next 3 lines are the set of basis vectors
    basis = Vector{Vector{Float64}}(undef, 3)
    for n in 1:3
        ln = iterate(iter)[1]
        basis[n] = parse.(Float64, split(ln))
    end
    ln = iterate(iter)[1]
    # TODO: how to handle the missing atomic symbols?
    # TODO: complete the function
end

readPOSCAR(filename::AbstractString) = open(readPOSCAR, filename)

"""
    readPOSCAR() -> CrystalStructure{3}

Reads any file named `POSCAR` in the current directory as a VASP POSCAR file.
"""
readPOSCAR() = readPOSCAR("POSCAR")


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

# Custom scene example code

using Pkg; Pkg.activate(joinpath("..", "..", ".."))

using Raytracer
using Raytracer.Interpreter

"""
    FoggyRenderer <: Renderer

A `Renderer` that returns the color of the `Shape` first hit by a given `Ray` mixed with the fog_color depending on distance travelled.

This renderer returns the color stored in the `material` field of the `Shape` first hit by the given `Ray` at the hit point.
To this renderer there is no difference between radiated light and reflected color. There are no shades, diffusions or reflections.
The color is then mixed with the `fog color` using the formula
``(1-\\exp{-t/λ}) fog_color + \\exp{-t/λ} hit_color``
where `t` is the distance to the hit point and `λ` is the luminosity falloff.
If there are no hits this renderer returns the value of its field `fog_color`.

# Fields

- `world::World`: the `World` to render.
- `fog_color::RGB{Float32}`: color of the fog.
- `falloff::Float32`: color falloff
"""
struct FoggyRenderer <: Renderer
	world::World
	fog_color::RGB{Float32}
	falloff::Float32
end

@doc """
    FoggyRenderer(world::World, fog_color::RGB{Float32}, falloff::Float32)

Constructor for a `FoggyRenderer` instance.
""" FoggyRenderer(::World, ::RGB{Float32}, ::Float32)

"""
    FoggyRenderer(world::World; fog_color::RGB{Float32} = BLACK, falloff::Float32 = 1f0)

Constructor for a `FoggyRenderer` instance.

If no color is specified, it will default on [`BLACK`](@ref).
"""
FoggyRenderer(world::World; fog_color::RGB{Float32} = BLACK, falloff::Float32 = 1f0) = FlatRenderer(world, fog_color, falloff)

"""
    (oor::FoggyRenderer)(ray::Ray)

Render a `Ray` and return a `RBG{Float32}`.
"""
function (fr::FoggyRenderer)(ray::Ray)
    hit = ray_intersection(ray, fr.world)
    isnothing(hit) && return fr.fog_color
	hit_color = hit.material.emitted_radiance(hit.surface_point) + hit.material.brdf.pigment(hit.surface_point)
	t = hit.t/fr.falloff
	fr.fog_color * (1 - exp(-t)) + hit_color * exp(-t)
end

scene = open_stream(parse_scene, "test_scene.sl")

image_tracer = ImageTracer(scene.image, scene.camera; scene.tracer.kwargs...)
renderer = FoggyRenderer(scene.world, WHITE, 10f0)
render(image_tracer, renderer, output_file="foggy_renderer_out.pfm")
tonemapping("foggy_renderer_out.pfm", "foggy_renderer_out.jpg", luminosity = 0.5f0)


text:
using LinearAlgebra: length
using ACE, ACEbase, Test, ACE.Testing
using ACE: evaluate, SymmetricBasis, PIBasis, O3
using StaticArrays

##

@info(" Testset for Multiple Properties in a Linear ACEModel")

# construct the 1p-basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)

# generate a configuration
#TODO check if this tests account/test u, and will it work with more things?
nX = 54
Xs = () -> ACE.State(rr = rand(SVector{3, Float64}), u = rand())
cfg = ACEConfig([Xs() for i in 1:nX])

φ = ACE.Invariant()
basis = SymmetricBasis(φ, B1p, O3(), Bsel)

##

BB = evaluate(basis, cfg)

c_m = rand(SVector{7,Float64}, length(BB))

singlProp = [ACE.LinearACEModel(basis, rand(length(BB)), evaluator = :standard) for i in 1:length(c_m[1])]
multiProp = ACE.LinearACEModel(basis, c_m, evaluator = :standard)

@info("set_params!")
c_s = [[c_m[j][i] for j in 1:length(c_m)] for i in 1:length(c_m[1])]

ACE.set_params!(multiProp,c_m)
for i in 1:length(c_m[1])
    ACE.set_params!(singlProp[i],c_s[i])
    print_tf(@test(c_s[i] ≈ singlProp[i].c))
end
println()

##

# We compare a model with a sinlge property, to the corresponding solution
# of the multiple property model. We do this by giving the same parameters
# to both models and then using isapprox

@info("evaluate")

for i in 1:length(c_m[1])
    print_tf(@test(evaluate(singlProp[i],cfg).val ≈ evaluate(multiProp,cfg)[i].val))
end
println()

@info("grad_params")

multiGradP = ACE.grad_params(multiProp,cfg)

for i in 1:length(c_m[1])
    singl = getproperty.(ACE.grad_params(singlProp[i],cfg),:val)
    multi = [getproperty(multiGradP[j][i], :val) for j in 1:length(c_m)]
    print_tf(@test(singl ≈ multi))
end
println()

##

@info("grad_config")

for i in 1:length(c_m[1])
    singl = ACE.grad_config(singlProp[i],cfg)
    multi = ACE.grad_config(multiProp,cfg)[:,i]

    print_tf(@test(singl ≈ multi))
end
println()

##

@info("adjoint_EVAL_D 1 prop")

#we check by contracting the full matrix that adjoint_eval_config works
#this assumes that grad_params_config works for one porperty. This is
#tested elsewhere.
for i in 1:length(c_m[1])

    #find the jacobian
    Jac = ACE.grad_params_config(singlProp[i],cfg)

    #create a random input emulating the pullback input
    w = rand(SVector{3, Float64}, length(Jac[1,:]))
    u = zeros(1)
    w = [ACE.DState(rr = w[j], u = u) for j in 1:length(w)]

    #calculate the adjoint and make a zeros to fill
    grad = ACE.adjoint_EVAL_D(singlProp[i], cfg, w)
    Jgrad = zeros(length(grad))

    #contract w into the jacobian to get the solution
    for j in 1:length(Jac[:,1])
        Jgrad[j] = sum([ACE.contract(Jac[j,:][k], w[k]) for k in 1:length(w)])
    end

    print_tf(@test(grad ≈ Jgrad))
end
println()

@info("adjoint_EVAL_D >2 prop")

#now we check that multiple properties work. For this, like before, we simply
#compare the single property result to each of the multiple properties.

uo = zeros(1)
function wMaker()
    wtmp = rand(SVector{3, Float64}, length(cfg))
    wtmp = [ACE.DState(rr = wtmp[j], u = uo) for j in 1:length(wtmp)]
    return wtmp
end
TDX1 = ACE.DState{NamedTuple{(:rr, :u), Tuple{SVector{3, Float64}, Vector{Float64}}}}
wo = Matrix{TDX1}(undef, (54,7))
wt  = [wMaker() for i in 1:length(c_m[1])]

for i in 1:length(wt)
    for j in 1:length(wt[i])
        wo[j,i] = wt[i][j]
    end
end

multiEval = ACE.adjoint_EVAL_D(multiProp, cfg, wo)

for i in 1:length(c_m[1])
    singl = ACE.adjoint_EVAL_D(singlProp[i], cfg, wo[:,i])
    multi = [multiEval[j][i] for j in 1:length(c_m)]
    print_tf(@test(singl ≈ multi))
end
println()


text:
@testset "Default post correlation Filter" begin
    post_corr_filter = Tracking.get_default_post_corr_filter(EarlyPromptLateCorrelator())
    @test @inferred(post_corr_filter(1.0 + 0.0im)) == 1.0 + 0.0im

    post_corr_filter = Tracking.get_default_post_corr_filter(
        EarlyPromptLateCorrelator(NumAnts(2))
    )
    @test @inferred(post_corr_filter(SVector(1.0 + 0.0im, 2.0 + 0.0im))) == 1.0 + 0.0im
end

@testset "Integration time" begin
    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(GPSL1, max_integration_time, bit_found)
    @test time == 1ms

    max_integration_time = 2ms
    bit_found = true
    time = @inferred Tracking.get_integration_time(GPSL1, max_integration_time, bit_found)
    @test time == 2ms

    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        GalileoE1B,
        max_integration_time,
        bit_found
    )
    @test time == 2ms

    max_integration_time = 10ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        GalileoE1B,
        max_integration_time,
        bit_found
    )
    @test time == 4ms
end

@testset "Number of chips to integrate" begin
    max_integration_time = 1ms
    code_phase = 200
    bit_found = true
    chips = @inferred Tracking.get_num_chips_to_integrate(
        GPSL1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 1023 - 200

    code_phase = 1200
    chips = @inferred Tracking.get_num_chips_to_integrate(
        GPSL1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 2046 - 1200
end

@testset "Number of samples to integrate" begin
    max_integration_time = 1ms
    sample_frequency = 4e6Hz
    current_code_doppler = 0Hz
    current_code_phase = 0
    bit_found = true
    samples = @inferred Tracking.get_num_samples_left_to_integrate(
        GPSL1,
        max_integration_time,
        sample_frequency,
        current_code_doppler,
        current_code_phase,
        bit_found
    )
    @test samples == 4000
end

@testset "Carrier phase delta" begin
    intermediate_frequency = 100Hz
    carrier_doppler = 50Hz
    frequency = @inferred Tracking.get_current_carrier_frequency(
        intermediate_frequency,
        carrier_doppler
    )
    @test frequency == intermediate_frequency + carrier_doppler
end

@testset "Code phase delta" begin
    code_doppler = 1Hz
    frequency = @inferred Tracking.get_current_code_frequency(GPSL1, code_doppler)
    @test frequency == 1023e3Hz + 1Hz
end

@testset "Update phases" begin
    carrier_phase = π / 2
    carrier_frequency = 10Hz
    sample_frequency = 100Hz
    num_samples = 2000
    phase = @inferred Tracking.update_carrier_phase(
        num_samples,
        carrier_frequency,
        sample_frequency,
        carrier_phase
    )
    @test phase == mod(π / 2 + 0.1 * 2000 + 1, 2) - 1

    code_phase = 10
    code_frequency = 10Hz
    sample_frequency = 100Hz
    num_samples = 2000
    bit_found = true
    phase = @inferred Tracking.update_code_phase(
        GPSL1,
        num_samples,
        code_frequency,
        sample_frequency,
        code_phase,
        bit_found
    )
    @test phase == mod(10 + 0.1 * 2000, 1023)
end

@testset "Correlate" begin
    correlator = EarlyPromptLateCorrelator()
    early_late_sample_shift = 2
    start_sample = 1
    num_samples = 4000
    carrier_frequency = 800e3Hz
    code_frequency = 1023e3Hz
    sample_frequency = 4000e3Hz
    start_carrier_phase = 0.25
    @testset "Correlate for satellite $prn" for prn in 1:32
        range = 0:3999
        signal = cis.(2π .* carrier_frequency .* range ./ sample_frequency .+ π ./ 2) .*
            get_code.(GPSL1, code_frequency .* range ./ sample_frequency .+ 0, prn)
        start_code_phase = 0

        correlator_result = @inferred Tracking.correlate(
            GPSL1,
            correlator,
            signal,
            prn,
            early_late_sample_shift,
            start_sample,
            num_samples,
            carrier_frequency,
            code_frequency,
            sample_frequency,
            start_carrier_phase,
            start_code_phase
        )
        difference_early_late = abs(
            get_early(correlator_result) - get_late(correlator_result)
        )
        # WHY??
        @test get_early(correlator_result) ≈ 1952 || get_early(correlator_result) ≈ 2080.0 ||
            get_early(correlator_result) ≈ 1824.0
        @test get_prompt(correlator_result) ≈ 4000
        @test get_late(correlator_result) ≈ 1952 || get_late(correlator_result) ≈ 2080.0 ||
            get_late(correlator_result) ≈ 1824.0
        @test difference_early_late ≈ 0 atol = 4.5e-11
    end

    range = 0:7999
    num_samples = 8000
    prn = 1
    signal = cis.(2π .* carrier_frequency .* range ./ sample_frequency .+ π ./ 2) .*
        get_code.(GPSL1, code_frequency .* range ./ sample_frequency .+ 0, prn)
    start_code_phase = 0
    correlator_result = @inferred Tracking.correlate(
        GPSL1,
        correlator,
        signal,
        prn,
        early_late_sample_shift,
        start_sample,
        num_samples,
        carrier_frequency,
        code_frequency,
        sample_frequency,
        start_carrier_phase,
        start_code_phase
    )
    @test get_early(correlator_result) ≈ 3904
    @test get_prompt(correlator_result) ≈ 8000
    @test get_late(correlator_result) ≈ 3904

end

@testset "Doppler aiding" begin

    init_carrier_doppler = 10Hz
    init_code_doppler = 1Hz
    carrier_freq_update = 2Hz
    code_freq_update = -0.5Hz
    velocity_aiding = 3Hz

    carrier_freq, code_freq = @inferred Tracking.aid_dopplers(
        GPSL1,
        init_carrier_doppler,
        init_code_doppler,
        carrier_freq_update,
        code_freq_update,
        velocity_aiding
    )

    @test carrier_freq == 10Hz + 2Hz + 3Hz
    @test code_freq == 1Hz + (2Hz + 3Hz) / 1540 - 0.5Hz
end

@testset "Signal" begin
    correlator = EarlyPromptLateCorrelator()
    signal = collect(1:10)
    sample = 5
    signal_sample = @inferred Tracking.get_signal(correlator, signal, sample)
    @test signal_sample == 5
    @test @inferred(Tracking.get_num_samples(signal)) == 10

    correlator = EarlyPromptLateCorrelator(NumAnts(2))
    signal = [1, 2] * collect(1:10)'
    sample = 5
    signal_sample = @inferred Tracking.get_signal(correlator, signal, sample)
    @test signal_sample == [5, 10]
    @test typeof(signal_sample) <: SVector

    @test @inferred(Tracking.get_num_samples(signal)) == 10
end

@testset "Tracking" begin

    carrier_doppler = 200Hz
    start_code_phase = 100
    code_frequency = carrier_doppler / 1540 + 1023kHz
    sample_frequency = 4MHz
    prn = 1
    range = 0:3999
    start_carrier_phase = π / 2
    state = TrackingState(GPSL1, carrier_doppler, start_code_phase)

    signal = cis.(
            2π .* carrier_doppler .* range ./ sample_frequency .+ start_carrier_phase
        ) .*
        get_code.(
            GPSL1,
            code_frequency .* range ./ sample_frequency .+ start_code_phase,
            prn
        )

    track_result = @inferred track(signal, state, prn, sample_frequency)

    iterations = 1000
    code_phases = zeros(iterations)
    carrier_phases = zeros(iterations)
    tracked_code_phases = zeros(iterations)
    tracked_carrier_phases = zeros(iterations)
    tracked_code_dopplers = zeros(iterations)
    tracked_carrier_dopplers = zeros(iterations)
    for i = 1:iterations
        carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sample_frequency +
            start_carrier_phase + π) - π
        code_phase = mod(
            code_frequency * 4000 * i / sample_frequency + start_code_phase,
            1023
        )
        signal = cis.(
                2π .* carrier_doppler .* range ./ sample_frequency .+ carrier_phase
            ) .*
            get_code.(
                GPSL1,
                code_frequency .* range ./ sample_frequency .+ code_phase,
                prn
            )
        track_result = @inferred track(
            signal,
            get_state(track_result),
            prn,
            sample_frequency
        )
        comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) /
            sample_frequency + start_carrier_phase + π) - π
        comp_code_phase = mod(
            code_frequency * 4000 * (i + 1) / sample_frequency + start_code_phase,
            1023
        )
        tracked_code_phases[i] = get_code_phase(track_result)
        tracked_carrier_phases[i] = get_carrier_phase(track_result)
        tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz
        tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz
        code_phases[i] = comp_code_phase
        carrier_phases[i] = comp_carrier_phase
    end
    @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4
    @test tracked_carrier_phases[end] ≈ carrier_phases[end] atol = 5e-2

#    using PyPlot
#    pygui(true)
#    figure("carrier_phases")
#    plot(tracked_carrier_phases - carrier_phases)
#    plot(carrier_phases)
#    figure("Code Phases")
#    plot(300 * (tracked_code_phases .- code_phases))
#    figure("Carrier Doppler")
#    plot(tracked_carrier_dopplers)
#    figure("Code Doppler")
#    plot(tracked_code_dopplers)

end


text:
module test_71_Tables_Rotational_First

using ModiaResult
using ModiaResult.DataFrames
using CSV
ModiaResult.@usingModiaPlot

result1 = CSV.File("$(ModiaResult.path)/test/test_71_Tables_Rotational_First.csv")
result2 = DataFrames.DataFrame(result1)

println("\n... test_71_Tables_Rotational_First.jl:")
println("CSV-Table (result1 = CSV.File(fileName)):\n")
ModiaResult.printResultInfo(result1)

println("\nDataFrame-Table (result2 = DataFrame(result1)):\n")
ModiaResult.printResultInfo(result2)

plot(result1, ["damper.w_rel", "inertia3.w"], prefix="result1: ")
plot(result2, ["damper.w_rel", "inertia3.w"], prefix="result2: ", reuse=true)

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
import Base.getindex
export getindex

getindex{T <: Real}(x::Constant, inds::AbstractArray{T, 1}) = Constant(getindex(x.value, inds))
getindex{T <: Real}(x::Constant, rows::AbstractArray{T, 1}, cols::AbstractArray{T, 1}) = Constant(getindex(x.value, rows, cols))

function getindex{T <: Real}(x::AffineExpr, inds::AbstractArray{T, 1})
  if length(inds) == 1
    # select only one entry of x
    ind = convert(Int64, inds[1])
    vars_to_coeffs_map = Dict{UInt64, Constant}()
    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = c[ind, :]
    end
    constant = x.constant[ind, :]
  else
    # number of rows/cols in the coefficient for x in our canonical form
    num_rows_coeff = length(inds)
    num_cols_coeff = x.size[1] * x.size[2]
    indexer = Constant(sparse(1:length(inds), inds, 1.0, num_rows_coeff, num_cols_coeff))

    vars_to_coeffs_map = Dict{UInt64, Constant}()
    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = indexer * c
    end
    constant = indexer * x.constant
  end
  this = AffineExpr(:index, (x,), vars_to_coeffs_map, constant, (length(inds), 1))
  this.evaluate = ()->x.evaluate()[inds]
  return this
end

function getindex{T <: Real}(x::AffineExpr, rows::AbstractArray{T, 1}, cols::AbstractArray{T, 1})
  if length(rows) == 1 && length(cols) == 1
    # select only one entry of x
    ind = x.size[1] * (convert(Int64, cols[1]) - 1) + convert(Int64, rows[1])
    vars_to_coeffs_map = Dict{UInt64, Constant}()
    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = c[ind, :]
    end
    constant = x.constant[ind, :]
  elseif length(cols) == 1 && typeof(rows) == UnitRange{Int64}
    # build a range object to just select the part of a column in x
    ind = x.size[1] * (convert(Int64, cols[1]) - 1)
    ind_range = (ind + rows[1]) : (ind + rows[end])
    vars_to_coeffs_map = Dict{UInt64, Constant}()

    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = c[ind_range, :]
    end
    constant = x.constant[ind_range, :]
  elseif length(rows) == 1 && typeof(cols) == UnitRange{Int64}
    # build a range object to just select part of a row in x
    start_ind = x.size[1] * (cols[1] - 1)
    end_ind = x.size[1] * (cols[end] - 1)
    ind_range = (start_ind + rows[1]) : x.size[1] : (end_ind + rows[1])

    vars_to_coeffs_map = Dict{UInt64, Constant}()
    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = c[ind_range, :]
    end
    constant = x.constant[ind_range, :]
  else
    # number of rows/cols in the coefficient for x in our canonical form
    num_rows_coeff = length(rows) * length(cols)
    num_cols_coeff = x.size[1] * x.size[2]

    # Create the indexing matrix such that indexer * vec(x) = vec(x[rows, cols])
    J = Array{Int64}(num_rows_coeff)
    k = 1
    num_rows = x.size[1]
    for c in cols
      for r in rows
        J[k] = num_rows * (convert(Int64, c) - 1) + convert(Int64, r)
        k += 1
      end
    end
    indexer = Constant(sparse(1:num_rows_coeff, J, 1.0, num_rows_coeff, num_cols_coeff))

    vars_to_coeffs_map = Dict{UInt64, Constant}()
    for (v, c) in x.vars_to_coeffs_map
      vars_to_coeffs_map[v] = indexer * c
    end
    constant = indexer * x.constant
  end

  this = AffineExpr(:index, (x,), vars_to_coeffs_map, constant, (length(rows), length(cols)))
  this.evaluate = ()->x.evaluate()[rows, cols]
  return this
end

getindex(x::AffineOrConstant, ind::Colon) = getindex(x, 1:size(x,1)*size(x,2))
getindex{T <: Real}(x::AffineOrConstant, rows::Colon, cols::AbstractArray{T, 1}) = getindex(x, 1:size(x, 1), cols)
getindex{T <: Real}(x::AffineOrConstant, rows::AbstractArray{T, 1}, col::Colon) = getindex(x, rows, 1:size(x,2))

getindex(x::AffineOrConstant, ind::Real) = getindex(x, ind:ind)
getindex(x::AffineOrConstant, row::Real, col::Real) = getindex(x, row:row, col:col)
getindex{T <: Real}(x::AffineOrConstant, row::Real, cols::AbstractArray{T, 1}) = getindex(x, row:row, cols)
getindex{T <: Real}(x::AffineOrConstant, rows::AbstractArray{T, 1}, col::Real) = getindex(x, rows, col:col)
getindex(x::AffineOrConstant, rows::Colon, col::Real) = getindex(x, 1:size(x,1), col:col)
getindex(x::AffineOrConstant, row::Real, cols::Colon) = getindex(x, row:row, 1:size(x,2))


text:
using Test
using PolynomialBasis
using IntervalArithmetic
# using Revise
using ImplicitDomainQuadrature
IDQ = ImplicitDomainQuadrature

box = IntervalBox([1.0, 2], [3.0, 5.0])
xL, xR = IDQ.corners(box)
@test all(xL .≈ [1.0, 2.0])
@test all(xR .≈ [3.0, 5.0])

b1, b2, b3, b4 = IDQ.split_box(box,2)
testb1 = IntervalBox([1.0, 2], [2, 3.5])
testb2 = IntervalBox([2.0, 2], [3, 3.5])
testb3 = IntervalBox([1.0, 3.5], [2, 5.0])
testb4 = IntervalBox([2.0, 3.5], [3.0, 5])
@test b1 == testb1
@test b2 == testb2
@test b3 == testb3
@test b4 == testb4

box = IntervalBox([0.,0.],[3.,3.])
b1,b2,b3,b4,b5,b6,b7,b8,b9 = IDQ.split_box(box,3)
testb1 = IntervalBox([0.,0.],[1.,1.])
testb2 = IntervalBox([1.,0.],[2.,1.])
testb3 = IntervalBox([2.,0.],[3.,1.])
testb4 = IntervalBox([0.,1.],[1.,2.])
testb5 = IntervalBox([1.,1.],[2.,2.])
testb6 = IntervalBox([2.,1.],[3.,2.])
testb7 = IntervalBox([0.,2.],[1.,3.])
testb8 = IntervalBox([1.,2.],[2.,3.])
testb9 = IntervalBox([2.,2.],[3.,3.])
@test b1 == testb1
@test b2 == testb2
@test b3 == testb3
@test b4 == testb4
@test b5 == testb5
@test b6 == testb6
@test b7 == testb7
@test b8 == testb8
@test b9 == testb9

box = IntervalBox(0..1, 0..2, 0..5)
@test IDQ.min_diam(box) ≈ 1.0

f(x) = x[1] + 1
box = IntervalBox(-1..1, 1)
s = IDQ.interval_arithmetic_sign_search(f, box, 1e-2,0.0,2)
@test !(s == 0 || s == 1 || s == -1)

f(x) = x[1] + 2
box = IntervalBox(-1..1, 1)
@test IDQ.interval_arithmetic_sign_search(f, box, 1e-2,0.0,2) == 1

f(x) = x[1] - 2
box = IntervalBox(-1..1, 1)
@test IDQ.interval_arithmetic_sign_search(f, box, 1e-2,0.0,2) == -1

f(x) = x[1]
box = IntervalBox(-1..1, 1)
@test IDQ.interval_arithmetic_sign_search(f, box, 1e-2,0.0,2) == 0
@test sign(f, [-1.], [1.], tol = 1e-2) == 0

f(x) = x[1] * (x[1] - 0.1) * (x[1] - 1.0)
@test sign(f, [-1.], [1.], tol = 1e-2) == 0

f(x) = (x[1] - 0.5) * (x[1] + 0.5) + 10
@test sign(f, [-1.], [1.], tol = 1e-2) == 1

f(x) = (x[1] - 0.5) * (x[1] + 0.5) - 3
@test sign(f, [-1.], [1.], tol = 1e-2) == -1

f(x) = x[1] * (x[1] - 1e-4)
s = sign(f, [-1.], [1.], tol = 1e-2)
@test !(s == 0 || s == +1 || s == -1)

box = IntervalBox(2..3, 1)
g(x) = (x[1] - 2.5)*(x[1] - 2.6)*(x[1] - 2.9)
@test sign(g, [2.], [3.], tol = 1e-2) == 0

r = 0.5
f2(x) = x[1]^2 + x[2]^2 - r^2
@test sign(f2, [-1.,-1.], [1.,1.], tol = 1e-2) == 0

g2(x) = sin(2pi * x[1]) * cos(2pi * x[2]) + 2
@test sign(g2, [-1.,-1.], [1.,1.], tol = 1e-2) == 1

P = InterpolatingPolynomial(1, 2, 2)
coeffs = 1:9
PolynomialBasis.update!(P, coeffs)
xL = [-1.,-1.]
xR = [1.,1.]
max_coeff, min_coeff = IDQ.extremal_coeffs_in_box(P, xL,xR)
@test max_coeff ≈ 9.0
@test min_coeff ≈ 1.0

xL = [0.5,0.5]
xR = [1.,1.]
max_coeff, min_coeff = IDQ.extremal_coeffs_in_box(P, xL,xR)
@test max_coeff ≈ 9.0
@test min_coeff ≈ 9.0

xL = [0.6,0.6]
xR = [0.7,0.7]
max_coeff, min_coeff = IDQ.extremal_coeffs_in_box(P, xL,xR)
@test max_coeff == -Inf
@test min_coeff == Inf

P = InterpolatingPolynomial(1, 2, 3)
coeffs = (1:16) .- 5.5
PolynomialBasis.update!(P, coeffs)
xL = [-1.,-1.]
xR = [1.,1.]
s = sign(P, xL, xR)
@test s == 0

P = InterpolatingPolynomial(1, 2, 3)
coeffs = 1.0:16.0
PolynomialBasis.update!(P, coeffs)
s = sign(P, xL,xR)
@test s == 1

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
include("startup.jl")

for col in names(tps)
    tps[!, col] = collect(tps[!, col])
end

tps_specific = stack(tps, ["ageMonths", "cogScore", "ECHOTPCoded", brainmeta...], [:subject, :timepoint])
subj_specific = stack(tps, [:mother_HHS_Education, :simple_race], [:subject])
grp = groupby(subj_specific, [:subject, :variable])
subj_specific = DataFrames.combine(grp, :value=> (v-> coalesce(v...))=> :value)

##

tps_specific = vcat(tps_specific, DataFrame(subject   = get(species, :subject),
                                timepoint = get(species, :timepoint),
                                variable  = fill("sample", nsamples(species)),
                                value     = samplenames(species))
)

for s in samples(species)
    subject   = s.subject
    timepoint = s.timepoint
    for f in features(species)
        variable = name(f)
        value = species[f, s]
        push!(tps_specific, (; subject, timepoint, variable, value))
    end
end

CSV.write("data/wrangled/tidy_timepoints_with_brain.csv", tps_specific)
CSV.write("data/wrangled/tidy_subjects.csv", subj_specific)

##

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
fs1()

@testset "LX input" begin
    set_curpath("index.md")
    mkpath(joinpath(F.PATHS[:assets], "index", "code", "output"))
    write(joinpath(F.PATHS[:assets], "index", "code", "s1.jl"), "println(1+1)")
    write(joinpath(F.PATHS[:assets], "index", "code", "output", "s1a.png"), "blah")
    write(joinpath(F.PATHS[:assets], "index", "code", "output", "s1.out"), "blih")
    st = raw"""
        Some string
        \input{julia}{s1.jl}
        Then maybe
        \output{s1.jl}
        Finally img:
        \input{plot:a}{s1.jl}
        done.
        """;

    F.def_GLOBAL_VARS!()
    F.def_GLOBAL_LXDEFS!()

    m = F.convert_md(st)
    h = F.convert_html(m)

    @test occursin("
Some string
$(read(joinpath(F.PATHS[:assets], "index", "code", "s1.jl"), String))
", h)
    @test occursin("Then maybe
$(read(joinpath(F.PATHS[:assets], "index", "code",  "output", "s1.out"), String))
", h)
    @test occursin("Finally img:  done.", h)
end

@testset  "Input MD" begin
    mkpath(joinpath(F.PATHS[:assets], "ccc"))
    fp = joinpath(F.PATHS[:assets], "ccc", "asset1.md")
    write(fp, "blah **blih**")
    st = raw"""
        Some string
        \textinput{ccc/asset1}
        """
    @test isapproxstr(st |> conv, "
Some string blah blih
")
end

fs2()

@testset "LX input" begin
    set_curpath("index.md")
    mkpath(joinpath(F.PATHS[:site], "assets", "index", "code", "output"))
    write(joinpath(F.PATHS[:site], "assets", "index", "code", "s1.jl"), "println(1+1)")
    write(joinpath(F.PATHS[:site], "assets", "index", "code", "output", "s1a.png"), "blah")
    write(joinpath(F.PATHS[:site], "assets", "index", "code", "output", "s1.out"), "blih")
    st = raw"""
        Some string
        \input{julia}{s1.jl}
        Then maybe
        \output{s1.jl}
        Finally img:
        \input{plot:a}{s1.jl}
        done.
        """;

    F.def_GLOBAL_VARS!()
    F.def_GLOBAL_LXDEFS!()

    m = F.convert_md(st)
    h = F.convert_html(m)

    @test occursin("
Some string
$(read(joinpath(F.PATHS[:site], "assets", "index", "code", "s1.jl"), String))
", h)
    @test occursin("Then maybe
$(read(joinpath(F.PATHS[:site], "assets", "index", "code",  "output", "s1.out"), String))
", h)
    @test occursin("Finally img:  done.", h)
end

@testset "Input MD" begin
    mkpath(joinpath(F.PATHS[:site], "assets", "ccc"))
    fp = joinpath(F.PATHS[:site], "assets", "ccc", "asset1.md")
    write(fp, "blah **blih**")
    st = raw"""
        Some string
        \textinput{ccc/asset1}
        """
    @test isapproxstr(st |> conv, "
Some string blah blih
")
end

@testset "Input err" begin
    gotd()
    s = raw"""
        AA
        \input{julia}{foo/baz}
        \input{plot}{foo/baz}
        \textinput{foo/bar}
        """ |> fd2html_td
    @test isapproxstr(s, """

AA

// Couldn't find a file when trying to resolve an input request with relative path: `foo/baz`. //

// Couldn't find an output directory associated with 'foo/baz' when trying to input a plot. //

// Couldn't find a file when trying to resolve an input request with relative path: `foo/bar`. //

        """)

    fs2()
    gotd()
    # table input
    mkpath(joinpath(td, "_assets", "index", "output"))
    write(joinpath(td, "_assets", "index", "output", "foo.csv"), "bar")
    s = raw"""
        @def fd_rpath = "index.md"
        \tableinput{}{./foo.csv}
        """ |> fd2html_td
    @test isapproxstr(s, """

// Table matching '/assets/index/foo.csv' not found. //

        """)
end


text:
"""
    module NumpyIO

Functions to read and write Numpy .npy and .npz files.
"""
module NumpyIO

include("npy.jl")
include("npz.jl")

end


text:
using MetacommunityDynamics
using Dispersal
using DataFrames
using Plots
using Distributions

classiclevinsmodel = (; C,E, kw...) ->
    LevinsColonization{:O}(probability=C) +
    RandomExtinction{:O}(probability=E);

spatiallyexplicitlevinsmodel = (; C,E, kw...) ->
    SpatiallyExplicitLevinsColonization{:O}(probability=C) +
    AbioticExtinction{:O}(:A, baseprobability=E);

paramset = (
    C = 0:0.01:1,
    E = 0:0.01:1,
)

pops = generate(PoissonProcess, 20)
scatter(pops[:,1], pops[:,2], aspectratio=1, ms=[rand(Exponential(3)) for i in 1:20])

popsgrid = pointstogrid(pops, gridsize=100)
areas = generate(StaticEnvironmentalLayer, zeros(size(popsgrid)), Exponential(3), mask=popsgrid)
heatmap(areas)

initocc = rand(Occupancy, 0.8, 5,5)
arrayout = ArrayOutput((O=initocc, A=areas ), tspan=1:3)
sim!(arrayout, rs)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import ForwardDiff
import Zygote

using FMI.FMIImport: fmi2SampleDirectionalDerivative, fmi2GetJacobian, fmi2SetContinuousStates

FMUPaths = [get_model_filename("SpringFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"]),
            get_model_filename("BouncingBall1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])]

t_start = 0.0
t_step = 0.01
t_stop = 5.0
tData = t_start:t_step:t_stop

for FMUPath in FMUPaths
    myFMU = fmiLoad(FMUPath)
    comp = fmiInstantiate!(myFMU; loggingOn=false)

    fmiSetupExperiment(comp, t_start, t_stop)
    fmiEnterInitializationMode(comp)
    fmiExitInitializationMode(comp)

    x0 = fmiGetContinuousStates(comp)
    numStates = length(x0)

    dx = zeros(numStates)
    t = 0.0
    p = []

    # Jacobians for x0
    FD_jac = ForwardDiff.jacobian(x -> FMI.fx(comp, dx, x, p, t), x0)
    ZG_jac = Zygote.jacobian(FMI.fx, comp, dx, x0, p, t)[3]
    fmiSetContinuousStates(comp, x0)
    samp_jac = fmi2SampleDirectionalDerivative(comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)
    auto_jac = fmi2GetJacobian(comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)

    @test (abs.(auto_jac -   FD_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)
    @test (abs.(auto_jac -   ZG_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)
    @test (abs.(auto_jac - samp_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)

    # Jacobians for random x0 / dx
    x0 = x0 + rand(numStates)
    dx = dx + rand(numStates)
    FD_jac = ForwardDiff.jacobian(x -> FMI.fx(comp, dx, x, p, t), x0)
    ZG_jac = Zygote.jacobian(FMI.fx, comp, dx, x0, p, t)[3]
    fmi2SetContinuousStates(comp, x0)
    samp_jac = fmi2SampleDirectionalDerivative(comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)
    auto_jac = fmi2GetJacobian(comp, comp.fmu.modelDescription.derivativeValueReferences, comp.fmu.modelDescription.stateValueReferences)

    @test (abs.(auto_jac -   FD_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)
    @test (abs.(auto_jac -   ZG_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)
    @test (abs.(auto_jac - samp_jac) .< ones(numStates, numStates).*1e-6) == ones(Bool, numStates, numStates)

    fmiUnload(myFMU)
end


text:
macro import_painter_huge()
    :(
        using ComplexVisual:
                CV_Painter, cv_paint, CV_CanvasPainter, CV_2DCanvasPainter,
                CV_FillPainter, CV_ValueMarkPainter,
                CV_GridPainter, CV_LineSegment, CV_LineSegments,
                CV_LinePainter, CV_DirectionPainter,
                CV_PortraitPainter, cv_clear_cache,
                CV_CombiPainter, →, CV_StyledPainter, ↦,
                CV_Math2DCanvasPainter,
                cv_parallel_lines, cv_arc_lines, cv_star_lines,
                cv_star_arc_lines
    )
end

import Base:show

"""
`CV_Painter`: Instances of this type can paint inside objects with math
coordinate systems (e.g. `CV_Math2DCanvas`).
"""
abstract type CV_Painter end

show(io::IO, p::CV_Painter) = cv_show_impl(io, p)

"""
`CV_CanvasPainter`: Instances of this type can paint inside canvases with
math coordinate systems (e.g. `CV_Math2DCanvas`).
"""
abstract type CV_CanvasPainter <: CV_Painter end

"""
`CV_2DCanvasPainter`: Instances of this type can paint inside 2D canvases
with math coordinate systems (e.g. `CV_Math2DCanvas`).
"""
abstract type CV_2DCanvasPainter <: CV_CanvasPainter end

cv_clear_cache(canvas_painter::CV_2DCanvasPainter) = nothing

"""
```
CV_FillPainter <: CV_2DCanvasPainter
    no fields
```

A painter filling the complete canvas.

A `CV_ContextStyle` is typically used to govern the appearance of
the filling operation.
"""
struct CV_FillPainter <: CV_2DCanvasPainter # {{{
end

"""
```
cv_paint(cc, fill_painter)
    cc            CV_2DCanvasContext
    fill_painter  CV_FillPainter
```

fill the complete (math) canvas.
"""
function cv_paint(cc::CV_2DCanvasContext, fill_painter::CV_FillPainter)
    canvas = cc.canvas
    rectangle(cc.ctx,
        real(canvas.corner_ul), imag(canvas.corner_lr),
        real(canvas.corner_lr) - real(canvas.corner_ul),
        imag(canvas.corner_ul) - imag(canvas.corner_lr))
    fill(cc.ctx)
    return nothing
end
# }}}

"""
```
CV_ValueMarkPainter{N<:Number}
    where     CV_TranslateByOffset{N}
    start     Float64
    len       Float64
    vertical  Bool
```

A painter drawing a vertical line-segment
`(where, start)` to `(where, start + len)` (in the math coordinate system).
Here the `where` is the where-value of a `CV_TranslateByOffset`.
If `vertical` is `true` then the real- and imag-coordinates are
swapped for drawing.

This painter is typically used for Sliders to mark the current
slider position.
"""
struct CV_ValueMarkPainter{N<:Number} <: CV_2DCanvasPainter # {{{
    where     :: CV_TranslateByOffset{N}
    start     :: Float64
    len       :: Float64
    vertical  :: Bool

    function CV_ValueMarkPainter(where::CV_TranslateByOffset{N},
            start::Real, len::Real, vertical::Bool=true) where {N<:Number}
        return new{N}(where, Float64(start), Float64(len), vertical)
    end
end

"""
```
cv_paint(cc, mark_painter)
    cc              CV_2DCanvasContext
    mark_painter    CV_ValueMarkPainter
```

draw a horitzonal or vertical mark to the painter's `where` position.
"""
function cv_paint(cc::CV_2DCanvasContext, mark_painter::CV_ValueMarkPainter)
    canvas, ctx = cc.canvas, cc.ctx
    if mark_painter.vertical
        move_to(ctx, mark_painter.start, mark_painter.where.value)
        rel_line_to(ctx, mark_painter.len, 0)
    else
        move_to(ctx, mark_painter.where.value, mark_painter.start)
        rel_line_to(ctx, 0, mark_painter.len)
    end
    stroke(ctx)
    return nothing
end
# }}}

"""
```
CV_GridPainter <: CV_2DCanvasPainter
    reals   Vector{Float64}
    imags   Vector{Float64}
```

A painter for drawing horizontal and vertical grid lines.
"""
struct CV_GridPainter <: CV_2DCanvasPainter  # {{{
    reals   :: Vector{Float64}
    imags   :: Vector{Float64}
end

"""
`CV_GridPainter(reals, imags)`

```
reals    Union{NTuple{N, Real}, AbstractVector{Float64}}
imags    Union{NTuple{M, Real}, AbstractVector{Float64}}
```

construct grid for given real- and imag-values.
"""
function CV_GridPainter(
        reals::Union{NTuple{N, Real}, AbstractVector{Float64}},
        imags::Union{NTuple{M, Real}, AbstractVector{Float64}}) where {M, N}
    return CV_GridPainter(
        [Float64(x) for x in reals],
        [Float64(y) for y in imags])
end

show(io::IO, m::MIME{Symbol("text/plain")}, gp::CV_GridPainter) =
    cv_show_impl(io, m, gp)

"""
```
cv_paint(cc, grid_painter)
    cc              CV_2DCanvasContext
    grid_painter    CV_GridPainter
```

draw horizontal and vertical grid lines.
"""
function cv_paint(cc::CV_2DCanvasContext, grid_painter::CV_GridPainter) # {{{
    canvas, ctx = cc.canvas, cc.ctx
    left, right = real(canvas.corner_ul), real(canvas.corner_lr)
    top, bottom = imag(canvas.corner_ul), imag(canvas.corner_lr)
    for x in grid_painter.reals
        move_to(ctx, x, bottom)
        line_to(ctx, x, top)
        stroke(ctx)
    end
    for y in grid_painter.imags
        move_to(ctx, left, y)
        line_to(ctx, right, y)
        stroke(ctx)
    end
    return nothing
end # }}}
# }}}

"""
`CV_LineSegment`: list/vector of points in complex plane.
"""
const CV_LineSegment = Vector{Complex{Float64}}

"""
`CV_LineSegments`: list/vector of `CV_LineSegment`
"""
const CV_LineSegments = Vector{CV_LineSegment}

"""
```
CV_LinePainter{dtrafoT, strafoT, scutT} <: CV_2DCanvasPainter
    dst_trafo           dtrafoT
    segments            CV_LineSegments
    auto_close_path     Bool
    src_trafo           strafoT
    src_cut_test        scutT
```

Painting `CV_LineSegments` as curves.

All the points of the line segments are transformed according to `dst_trafo`.

The preimage of the line segments are transformed by `src_trafo` (unless
`src_trafo === nothing`) and afterwards they are checked with `src_cut_test`
(unless `src_cut_test === nothing`) if the line-segment needs to be "cut"
(i.e. the line is interrupted there).
"""
struct CV_LinePainter{dtrafoT, strafoT, scutT} <: CV_2DCanvasPainter # {{{
    dst_trafo       :: dtrafoT
    segments        :: CV_LineSegments
    auto_close_path :: Bool
    src_trafo       :: strafoT
    src_cut_test    :: scutT
end

"""
```
CV_LinePainter(dst_trafo, segments, auto_close_path=false)
    dst_trafo
    segments           CV_LineSegments
    auto_close_path    Bool
```

with `src_trafo` and `src_cut_test` both `nothing`.
"""
function CV_LinePainter(dst_trafo, segments::CV_LineSegments,
        auto_close_path::Bool=false)
    return CV_LinePainter(dst_trafo, segments, auto_close_path,
        nothing, nothing)
end

"""
```
CV_LinePainter(segments, auto_close_path=false)
    segments           CV_LineSegments
    auto_close_path    Bool
```

`dst_trafo = identity` and with `src_trafo` and `src_cut_test` both `nothing`.
"""
CV_LinePainter(segments::CV_LineSegments, auto_close_path::Bool=false) =
    CV_LinePainter(identity, segments, auto_close_path)

"""
```
cv_paint(cc, line_painter)
    cc            CV_2DCanvasContext,
    line_painter  CV_LinePainter{dtrafoT, Nothing, Nothing}
```

paint line segments without cut-test.
"""
function cv_paint(cc::CV_2DCanvasContext,
        line_painter::CV_LinePainter{dtrafoT,
                Nothing, Nothing}) where {dtrafoT}  # {{{
    dst_trafo = line_painter.dst_trafo
    ctx = cc.ctx
    for segment in line_painter.segments
        value = dst_trafo(segment[1])
        move_to(ctx, real(value), imag(value))
        for point in segment[2:end]
            value = dst_trafo(point)
            line_to(ctx, real(value), imag(value))
        end
        if line_painter.auto_close_path
            close_path(ctx)
        end
        stroke(ctx)
    end
    return nothing
end # }}}

"""
```
cv_paint(cc, line_painter)
    cc            CV_2DCanvasContext,
    line_painter  CV_LinePainter{dtrafoT, strafoT, scutT})
```

paint line segments with cut-test.
"""
function cv_paint(cc::CV_2DCanvasContext,
        line_painter::CV_LinePainter{dtrafoT,
                strafoT, scutT}) where {dtrafoT, strafoT, scutT}  # {{{

    ctx, trafo = cc.ctx, line_painter.dst_trafo
    src_trafo, cut_test = line_painter.src_trafo, line_painter.src_cut_test

    for segment in line_painter.segments
        from_point, from_value = src_trafo(segment[1]), trafo(segment[1])
        first, cut = true, false
        for point in segment
            to_point, to_value = src_trafo(point), trafo(point)
            if cut_test(from_point, to_point)
                first, cut = true, true
                from_point, from_value = to_point, to_value
            else
                if first
                    move_to(ctx, real(to_value), imag(to_value))
                    first = false
                else
                    line_to(ctx, real(to_value), imag(to_value))
                    from_point, from_value = to_point, to_value
                end
            end
        end
        line_painter.auto_close_path && !cut && close_path(ctx)
        stroke(ctx)
    end

    return nothing
end # }}}
# }}}

"""
```
CV_DirectionPainter{trafoT} <: CV_2DCanvasPainter
    trafo              trafoT
    segments           CV_LineSegments
    auto_close_path    Bool
    every_len          Float64
    pre_gap            Float64
    arrow              ComplexF64
```

Paint triangles (along a curve) to indicate the direction of a curve.
"""
struct CV_DirectionPainter{trafoT} <: CV_2DCanvasPainter # {{{
    trafo           :: trafoT
    segments        :: CV_LineSegments
    auto_close_path :: Bool
    every_len       :: Float64
    pre_gap         :: Float64
    arrow           :: ComplexF64
end

"""
```
CV_DirectionPainter(trafo, segments, auto_close_path=false;
        every_len=2.0, pre_gap=0.0, arrow=0.3*exp(1im*π*8/9))
    trafo
    segments          CV_LineSegments
    auto_close_path   Bool
    every_len         Real
    pre_gap           Real
    arrow             ComplexF64
```

The `arrow` parameter describes the shape of the triangle:

```
      arrow
        *   ^
        |\\  │
        |░\\ │
        |░░\\│
    ────|░░░┼────>
        |░░/│
        |░/ │
        |/  │
        *   │
     conj(arrow)
```

Moves along the curve. Places after `every_len` a arrow/triangle showing in
the direction of the curve. `pre_gap` can be used to prohibit a
triangle at the beginning (with length `pre_gap`).
"""
CV_DirectionPainter(trafo, segments::CV_LineSegments,
        auto_close_path::Bool=false;
        every_len::Real=2.0, pre_gap::Real=0.0,
        arrow::ComplexF64=0.3*exp(1im*π*8/9)) =
    CV_DirectionPainter(trafo, segments, auto_close_path,
        Float64(every_len), Float64(pre_gap), ComplexF64(arrow))

"""
```
cv_paint(cc, ldirp)
    cc       CV_2DCanvasContext,
    ldirp    CV_DirectionPainter
```

paint traingle/arrows along the curve.
"""
function cv_paint(cc::CV_2DCanvasContext,
        ldirp::CV_DirectionPainter) # {{{
    ctx, trafo, arrow = cc.ctx, ldirp.trafo, ldirp.arrow

    needed_len = ldirp.pre_gap + ldirp.every_len
    for segment in ldirp.segments
        wold = trafo(segment[1])
        for point in segment[2:end]
            wnew = trafo(point)
            wdiff = wnew - wold
            needed_len -= abs(wdiff)
            if needed_len ≤ 0
                needed_len = ldirp.every_len

                zcenter = wold + 0.5*wdiff
                move_to(ctx, real(zcenter), imag(zcenter))

                wnorm = wdiff/abs(wdiff)

                z = zcenter + arrow * wnorm
                line_to(ctx, real(z), imag(z))

                z = zcenter + conj(arrow) * wnorm
                line_to(ctx, real(z), imag(z))

                close_path(ctx)
                fill(ctx)
            end
            wold = wnew
        end
    end
    return nothing
end # }}}
# }}}

"""
```
CV_Math2DCanvasPainterCache
    last_canvas
    last_data
    last_color_matrix    Matrix{UInt32}
```

A cache for caching the color of every pixel in a (math) canvas.
`last_canvas` and `last_data` is used to save the dependencies for
the pixel colors in order to make recomputations possible if something
changes.
"""
mutable struct CV_Math2DCanvasPainterCache # {{{
    last_canvas
    last_data
    last_color_matrix :: Matrix{UInt32}
end

"""
`CV_Math2DCanvasPainterCache()`

init cache.
"""
CV_Math2DCanvasPainterCache() = CV_Math2DCanvasPainterCache(
    nothing, nothing, zeros(UInt32, 1, 1))

"""
```
cv_is_in_cache(cache, canvas, data)
    cache    CV_Math2DCanvasPainterCache
    canvas
    data
```

checks if the cache contains color values "matching" to `canvas` and `data`.
"""
cv_is_in_cache(cache::CV_Math2DCanvasPainterCache, canvas, data) = (
    cache.last_canvas === canvas && cache.last_data === data)
# }}}

"""
```
CV_PortraitPainter{CS}
    trafo           trafoT
    colorscheme     CS
    cache_flag      Bool
    cache           CV_Math2DCanvasPainterCache
```

Fill math coordinate system with phase portrait.

For `colorscheme`, please see the package `ComplexPortraits`.

"""
struct CV_PortraitPainter{trafoT, CS} <: CV_2DCanvasPainter  # {{{
    trafo       :: trafoT
    colorscheme :: CS
    cache_flag  :: Bool
    cache       :: CV_Math2DCanvasPainterCache
end

"""
```
CV_PortraitPainter([trafo=identity, [colorscheme=ComplexPortraits.cs_j(),
                                    [cache_flag=true]]])
    trafo         trafoT
    colorscheme   CS
    cache_flag    Bool
```
"""
function CV_PortraitPainter(trafo::trafoT=identity,
        colorscheme::CS=ComplexPortraits.cs_j(),
        cache_flag::Bool=true) where {trafoT, CS}
    return CV_PortraitPainter{trafoT, CS}(trafo, colorscheme, cache_flag,
        CV_Math2DCanvasPainterCache())
end

"""
```
cv_clear_cache(pp)
    pp    CV_PortraitPainter
```

invalidates the cache (by setting `last_canvas=noting`.
"""
function cv_clear_cache(pp::CV_PortraitPainter)
    if pp.cache_flag
        pp.cache.last_canvas = nothing
    end
    return nothing
end

"""
```
cv_paint(cc, portrait_painter)
    cc                   CV_2DCanvasContext
    portrait_painter     CV_PortraitPainter
```

fill the whole (math) canvas with the phase portrait.
"""
function cv_paint(cc::CV_2DCanvasContext{canvasT},
                  portrait_painter::CV_PortraitPainter{CS}
                  )  where {canvasT <: CV_Math2DCanvas, CS}   # {{{
    canvas = cc.canvas
    cache = portrait_painter.cache
    if portrait_painter.cache_flag && cv_is_in_cache(cache, canvas, nothing)
        color_argb32 = cache.last_color_matrix
    else
        color_matrix = ComplexPortraits.portrait(
            canvas.corner_ul, canvas.corner_lr,
            portrait_painter.trafo,
            no_pixels=(canvas.pixel_height, canvas.pixel_width),
            point_color=portrait_painter.colorscheme)
        color_argb32 = map(x -> convert(Colors.ARGB32, x).color, color_matrix) :: Matrix{UInt32}
        if portrait_painter.cache_flag
            cache.last_canvas, cache.last_data = canvas, nothing
            cache.last_color_matrix = color_argb32
        end
    end
    surface = cc.canvas.surface
    Cairo.flush(surface)
    surface.data .= color_argb32'
    Cairo.mark_dirty(surface)
    return nothing
end # }}}
# }}}

"""
```
CV_CombiPainter{T<:CV_Painter, S<:CV_Painter}
    painter1   T
    painter2   S
```

combines two painters.
"""
struct CV_CombiPainter{T<:CV_Painter, S<:CV_Painter} <: CV_Painter # {{{
    painter1 :: T
    painter2 :: S
end

function cv_clear_cache(cpainter::CV_CombiPainter)
    cv_clear_cache(cpainter.painter1)
    cv_clear_cache(cpainter.painter2)
    return nothing
end

function cv_paint(cc::CV_Context, cpainter::CV_CombiPainter)
    cv_paint(cc, cpainter.painter1)
    cv_paint(cc, cpainter.painter2)
    return nothing
end

"""
`→(painter1, painter2)  = CV_CombiPainter(painter1, painter2)`
"""
function →(painter1::T, painter2::S) where {T<:CV_Painter, S<:CV_Painter}
  return CV_CombiPainter(painter1, painter2)
end

"""
`→(nothing, painter2)  = painter2`
"""
→(painter1::Nothing, painter2::S) where {S<:CV_Painter} = painter2

"""
`→(painter1, nothing)  = painter1`
"""
→(painter1::T, painter2::Nothing) where {T<:CV_Painter} = painter1

show(io::IO, m::MIME{Symbol("text/plain")}, p::CV_CombiPainter) =
    cv_show_impl(io, m, p)

# }}}

"""
```
CV_StyledPainter{styleT<:CV_ContextStyle, painterT<:CV_Painter}
    style       styleT
    painter     painterT
```

Painter with Style :-)

This is *the* way to govern the appearance (color, thickness, etc.)
of the painting operations of a painter.
"""
struct CV_StyledPainter{styleT<:CV_ContextStyle,
                        painterT<:CV_Painter} <: CV_Painter   # {{{
    style    :: styleT
    painter  :: painterT
end

"""
`↦(style, painter) = CV_StyledPainter(style, painter)`
"""
↦(style::CV_ContextStyle, painter::CV_Painter) = CV_StyledPainter(
    style, painter)

show(io::IO, m::MIME{Symbol("text/plain")}, p::CV_StyledPainter) =
    cv_show_impl(io, m, p)

function cv_clear_cache(spainter::CV_StyledPainter)
    cv_clear_cache(spainter.painter)
    return nothing
end

function cv_paint(cc::CV_CanvasContext, styled_painter::CV_StyledPainter)
    cv_prepare(cc, styled_painter.style)
    cv_paint(cc, styled_painter.painter)
    return nothing
end # }}}

"""
```
CV_Math2DCanvasPainter <: CV_2DCanvasPainter
    dst_trafo       dtrafoT
    canvas          canvasT <: CV_Math2DCanvas
    src_trafo       strafoT
    src_cut_test    scutT
```
paint transformed (by `dst_trafo`) canvas in coordinate system.

The preimage of the pixel (of the canvas) are transformed by `src_trafo`
(unless `src_trafo === nothing`) and afterwards they are checked with
`src_cut_test` (unless `src_cut_test === nothing`) if the pixels
needs to be "cut" (i.e. left out).
"""
struct CV_Math2DCanvasPainter{dtrafoT, canvasT<:CV_Math2DCanvas,
                              strafoT, scutT} <: CV_2DCanvasPainter # {{{
    dst_trafo       :: dtrafoT
    canvas          :: canvasT
    src_trafo       :: strafoT
    src_cut_test    :: scutT
end

CV_Math2DCanvasPainter(trafo, canvas::CV_Math2DCanvas) =
    CV_Math2DCanvasPainter(trafo, canvas, nothing, nothing)

CV_Math2DCanvasPainter(canvas::CV_Math2DCanvas) =
    CV_Math2DCanvasPainter(identity, canvas)

show(io::IO, m::MIME{Symbol("text/plain")}, p::CV_Math2DCanvasPainter) =
    cv_show_impl(io, m, p)

"""
```
cv_paint(cc, canvas_painter)
    cc                CV_2DCanvasContext{canvasT <: CV_Math2DCanvas}
    canvas_painter    CV_Math2DCanvasPainter{dtrafoT, canvasT,
                            Nothing, Nothing})
```

Implementation without cut-test
"""
function cv_paint(cc::CV_2DCanvasContext{canvasT},
        canvas_painter::CV_Math2DCanvasPainter{dtrafoT, canvasT,
            Nothing, Nothing}) where {canvasT <: CV_Math2DCanvas, dtrafoT} # {{{
    ctx, trafo = cc.ctx, canvas_painter.dst_trafo
    canvas = canvas_painter.canvas
    width, height = canvas.pixel_width, canvas.pixel_height
    wh, hh = width/2, height/2
    punit = 1/canvas.resolution    # units per pxiel (at canvas)

    # missing optimization for:  trafo isa CV_TranslateByOffset
    data = canvas.surface.data
    for row in range(1, height; step=1)
        yb = (row-1-hh)*punit
        yt = yb + punit
        for col in range(1, width; step=1)
            color = ARGB32(data[col, row], Val{true})
            set_source_rgba(ctx,
                red(color), green(color), blue(color), alpha(color))

            xl = (col-1-wh)*punit
            xr = xl + punit

            z = trafo(xl + yb*1im) :: Complex{Float64}
            move_to(ctx, real(z), imag(z))

            z = trafo(xr + yb*1im) :: Complex{Float64}
            line_to(ctx, real(z), imag(z))

            z = trafo(xr + yt*1im) :: Complex{Float64}
            line_to(ctx, real(z), imag(z))

            z = trafo(xl + yt*1im) :: Complex{Float64}
            line_to(ctx, real(z), imag(z))

            close_path(ctx)
            fill(ctx)
        end
    end
    return nothing
end  # }}}

"""
```
cv_paint(cc, canvas_painter)
    cc                CV_2DCanvasContext{canvasT <: CV_Math2DCanvas}
    canvas_painter    CV_Math2DCanvasPainter{dtrafoT, canvasT,
                            strafoT, scut})
```

Implementation with cut-test
"""
function cv_paint(cc::CV_2DCanvasContext{canvasT},
        canvas_painter::CV_Math2DCanvasPainter{dtrafoT, canvasT,
            strafoT, scutT}) where {canvasT <: CV_Math2DCanvas, dtrafoT,
                strafoT, scutT} # {{{
    ctx, trafo, canvas  = cc.ctx, canvas_painter.dst_trafo, canvas_painter.canvas
    src_trafo, cut_test = canvas_painter.src_trafo, canvas_painter.src_cut_test
    width, height = canvas.pixel_width, canvas.pixel_height
    wh, hh = width/2, height/2
    punit = 1/canvas.resolution    # units per pxiel (at canvas)

    # missing optimization for:  trafo isa CV_TranslateByOffset
    data = canvas.surface.data
    for row in range(1, height; step=1)
        yb = (row-1-hh)*punit
        yt = yb + punit
        for col in range(1, width; step=1)
            xl = (col-1-wh)*punit
            xr = xl + punit

            w1 = trafo(xl + yb*1im) :: Complex{Float64}
            z1 = src_trafo(xl + yb*1im) :: Complex{Float64}

            w2 = trafo(xr + yb*1im) :: Complex{Float64}
            z2 = src_trafo(xr + yb*1im) :: Complex{Float64}

            w3 = trafo(xr + yt*1im) :: Complex{Float64}
            z3 = src_trafo(xr + yt*1im) :: Complex{Float64}

            w4 = trafo(xl + yt*1im) :: Complex{Float64}
            z4 = src_trafo(xl + yt*1im) :: Complex{Float64}

            if !cut_test(z1, z2) && !cut_test(z2, z3) &&
                    !cut_test(z3, z4) && !cut_test(z4, z1)
                color = ARGB32(data[col, row], Val{true})
                set_source_rgba(ctx,
                    red(color), green(color), blue(color), alpha(color))
                move_to(ctx, real(w1), imag(w1))
                line_to(ctx, real(w2), imag(w2))
                line_to(ctx, real(w3), imag(w3))
                line_to(ctx, real(w4), imag(w4))
                close_path(ctx)
                fill(ctx)
            end
        end
    end
    return nothing
end # }}}
# }}}

"""
```
cv_parallel_lines(direction;  width=1.0, lines=5, segments=120)
    direction   ComplexF64    direction for lines
    width       Real          space/width the lines cover
    lines       Integer       number of lines
    segments    Integer       how many segments per line
```

returns the `CV_LineSegments`.
"""
function cv_parallel_lines(direction::ComplexF64;
                           width::Real=1.0, lines::Integer=5,
                           segments::Integer=120)  # {{{
    width = Float64(width)
    if !isfinite(direction)
        cv_error("direction must be finite; direction = ", direction)
    end
    len = abs(direction)
    if len == 0
        cv_error("direction must not be 0")
    end
    if !isfinite(width)
        cv_error("width must be finite; width = ", width)
    end
    if !(width > 0)
        cv_error("width must be positive; width = ", width)
    end
    if !(lines > 0)
        cv_error("lines must be positive; lines = ", lines)
    end
    if !(segments > 0)
        cv_error("segments must be positive; segments = ", segments)
    end

    perp_dir = 1im*direction/len
    perp_parts = (lines == 1) ? (0.0:0.0) : range(-0.5*width, stop=0.5*width, length=lines)

    line_segs = CV_LineSegments()
    seg_parts = collect(0:segments) / segments
    for perp_part in perp_parts
        push!(line_segs, (perp_part*perp_dir - direction/2) .+ (seg_parts .* direction))
    end
    return line_segs
end # }}}

"""
```
cv_arc_lines(ϕ_start, ϕ_end, radii; segments=120)
    ϕ_start    Real                start angle
    ϕ_end      Real                end angle
    radii      NTuple{N, Real}     radii for the arcs
    segments   Integer             how many segments per line
```

create a `CV_LineSegments` with arcs starting at the angle `ϕ_start`
and ending at the angle `ϕ_end`. In `radii` one can give radii at which
the lines are created.
"""
function cv_arc_lines(ϕ_start::Real, ϕ_end::Real, radii::NTuple{N, Real};
                      segments::Integer=120) where {N} # {{{
    ϕ_start, ϕ_end = Float64(ϕ_start), Float64(ϕ_end)

    arc = exp.(1im .* range(ϕ_start, stop=ϕ_end, length=1+segments))
    line_segs = [Float64(radius) * arc for radius in radii]
    return line_segs
end # }}}

"""
```
cv_star_lines(r_start, r_end, angles; segments=120)
    r_start     Real                start radius
    r_end       Real                end radius
    angles      NTuple{N, Real}     lines at this angles
    segments    Integer             how many segments per line
```
create a `CV_LineSegments` with star lines starting the radius `r_start`
and ending at the radius `r_end`. In `angles` one can give angles at
which the lines are created.
"""
function cv_star_lines(r_start::Real, r_end::Real, angles::NTuple{N, Real};
                       segments::Integer=120) where {N}  # {{{
    r_start, r_end = Float64(r_start), Float64(r_end)

    line = collect(range(r_start, stop=r_end, length=1+segments))

    line_segs = [exp(1im*Float64(ϕ)) * line for ϕ in angles]
    return line_segs
end # }}}

"""
```
cv_star_arc_lines(radii, angle; star_segments=120, arc_segments=120)
    radii           NTuple{N, Real}
    angles          NTuple{M, Real};
    star_segments   Integer
    arc_segments    Integer
```

combination of `cv_star_lines` and `cv_arc_lines`.
"""
function cv_star_arc_lines(radii::NTuple{N, Real}, angles::NTuple{M, Real};
                star_segments::Integer=120,
                arc_segments::Integer=120) where {M, N}
    r_start, r_end = minimum(radii), maximum(radii)
    ϕ_start, ϕ_end = minimum(angles), maximum(angles)

    return (cv_star_lines(r_start, r_end, angles; segments=star_segments),
            cv_arc_lines(ϕ_start, ϕ_end, radii; segments=arc_segments))
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
function epistatic_score(arnet::ArNet, arvar::ArVar, seqid::Int; pc::Float64=0.1,min_separation::Int=1)
    @extract arnet : H J p0 idxperm
    @extract arvar : Z M N q

    1 ≤ seqid ≤ M || error("seqid=$seqid should be in the interval [1,...,$M]")

    Da = zeros(q,N)
    Dab = zeros(q,q,N,N)

    xori = Z[:,seqid]
    xmut = [copy(xori) for _ in 1:Threads.nthreads()]
    arlike = [zeros(N) for _ in 1:Threads.nthreads()]
    ppc = (1-pc) * p0 + pc * ones(q)/q
    #_outputarnet!(arlike,xmut, J, H, ppc, N, q)
    # E0 =  sum(log.(arlike))
    # E0=0.0
    @inbounds for i in 1:N
        Threads.@threads for a in 1:q
            xmut[Threads.threadid()][i] = a
            _outputarnet!(arlike[Threads.threadid()],xmut[Threads.threadid()], J, H, ppc, N, q)
            Da[a,i] = -sum(log.(arlike[Threads.threadid()]))
            xmut[Threads.threadid()][i] = xori[i]
        end
    end
    @inbounds for i in 1:N-1
        Threads.@threads for a in 1:q
            xmut[Threads.threadid()][i] = a
            for j in i+1:N
                 for b in 1:q
                    xmut[Threads.threadid()][i] = a
                    xmut[Threads.threadid()][j] = b
                    _outputarnet!(arlike[Threads.threadid()],xmut[Threads.threadid()],J,H,ppc,N,q)
                    Dab[b,a,j,i] = -sum(log.(arlike[Threads.threadid()]))
                    xmut[Threads.threadid()][i] = xori[i]
                    xmut[Threads.threadid()][j] = xori[j]
                end
            end
        end
    end

    Jret = zeros(q,q,N,N)
    @inbounds for i in 1:N-1
        for j in i+1:N
            for a in 1:q
                for b in 1:q
                    Jret[b,a,j,i] = Dab[b,a,j,i] - Da[b,j] - Da[a,i]
                    #Jret[a,b,i,j] = Jret[b,a,j,i]
                end
            end
        end
    end
    Jzsg = zsg(Jret)
    FN = compute_APC(Jzsg, N, q)
    score = compute_ranking(FN, min_separation)

    permtuple=Tuple{Int,Int,Float64}[]
    sizehint!(permtuple,length(permtuple))
    for s in score
        si,sj,val= idxperm[s[1]],idxperm[s[2]],s[3]
        if si > sj
            si,sj = sj,si
        end
        push!(permtuple,(si,sj,val))
    end
    return permtuple
end

function zsg(J::Array{Float64,4})
    q,q,N,N = size(J)
    Jzsg = zeros(q,q,N,N)
    @inbounds for i in 1:N-1
        for j in i+1:N
            Jzsg[:,:,j,i] .= J[:,:,j,i] - repeat(sum(J[:,:,j,i], dims=1)/q, q, 1) - repeat(sum(J[:,:,j,i], dims=2)/q, 1, q) .+ sum(J[:,:,j,i])/q^2
            #Jzsg[:,:,i,j] .= Jzsg[:,:,j,i]'
        end
    end
    Jzsg
end

function compute_APC(J::Array{Float64,4},N,q)
    FN = fill(0.0, N,N)
    @inbounds for i=1:N-1
        for j=i+1:N
            FN[j,i] = norm(J[1:q-1,1:q-1,j,i],2)
            FN[i,j] = FN[j,i]
        end
    end
    FN = correct_APC(FN)
    return FN
end

function correct_APC(S::Matrix)
    N = size(S, 1)
    Si = sum(S, dims=1)
    Sj = sum(S, dims=2)
    Sa = sum(S) * (1 - 1/N)

    S -= (Sj * Si) / Sa
    return S
end

function compute_ranking(S::Matrix{Float64}, min_separation::Int = 5)

    N = size(S, 1)
    R = Array{Tuple{Int,Int,Float64}}(undef, div((N-min_separation)*(N-min_separation+1), 2))
    counter = 0
    for i = 1:N-min_separation, j = i+min_separation:N
        counter += 1
        R[counter] = (i, j, S[j,i])
    end

    sort!(R, by=x->x[3], rev=true)
    return R

end

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
module DistQuads # DistQuads

using Reexport,
      FastGaussQuadrature
@reexport using Distributions
export DistQuad, E, mean, var
import Statistics: mean, var
import LinearAlgebra: dot
import SpecialFunctions: gamma

struct DistQuad{D}
    x
    w
    d::D
end
nodes(dq::DistQuad) = dq.x
weights(dq::DistQuad) = dq.w
distribution(dq::DistQuad) = dq.d

DistQuad(d; N = 32) = DistQuad(d, N)
function DistQuad(d::Distributions.Beta, N)
    gj = FastGaussQuadrature.JacobiRec(N, d.α-1, d.β-1)
    G = gamma(d.α)*gamma(d.β)/gamma(d.α+d.β)
    w = gj[2]/((2.0^(d.α+d.β - 1.0))*G)
    n = (1 .- gj[1])./2
    DistQuad(n, w, d)
end

function DistQuad(d::Distributions.Exponential, N)
    gl = gausslaguerre(N)
    w = gl[2]
    n = gl[1]./d.θ
    DistQuad(n, w, d)
end

function DistQuad(d::Distributions.Gamma, N)
    gl = gausslaguerre(N, d.α-1)
    n = gl[1].*d.θ
    w = gl[2]./gamma(d.α)
    DistQuad(n, w, d)
end

function _DistQuad_Normal_kernel(d, N)
    gh = gausshermite(N)
    w = gh[2]./sqrt(pi)
    n = gh[1].*sqrt(2).*d.σ.+d.μ
    n, w
end
function DistQuad(d::Distributions.Normal, N)
    n, w = _DistQuad_Normal_kernel(d, N)
    DistQuad(n, w, d)
end

function DistQuad(d::Distributions.LogNormal, N)
    n, w = _DistQuad_Normal_kernel(d, N)
    DistQuad(exp.(n), w, d)
end

E(f, dq::DistQuad) = dot(f.(nodes(dq)), weights(dq))
E(f, d) = E(f, DistQuad(d))

mean(dq::DistQuad) = E(identity, dq)
var(dq::DistQuad) = E(x->(x-mean(dq))^2, dq)

end # module


text:
const mac_gdrive = "/Volumes/GoogleDrive/My Drive"
const linux_gdrive = "/home/beatriz/mainen.flipping.5ht@gmail.com"
const files_dir = "Flipping/Datasets/Pharmacology/SwitchingData/Results"
const figs_dir = "Flipping/Datasets/Pharmacology/SwitchingData/Results/figures/"

const columns_types = Dict(
    :Protocol => Union{String,Missing},
    )

const drug_colors = Dict(
    "Altanserin" => RGB(0.6,0.149,0.561),
    "SB242084" => RGB(0.902,0.604,0.737),
    "Citalopram" => RGB(0.1216,0.4667,0.7059),#RGB(0.0902,0.7451,0.8118)
    "Methysergide" => RGB(0.1725,0.6275,0.1725),
    "Way_100135" => RGB(1.0,0.498,0.0549),
    "Optogenetic" => RGB(0.0,0.925,1.0),
    "SB242084_opt" => RGB(0.765,0.01,0.01),
    "training" => RGB(0.498,0.498,0.498),
    "PreVehicle" => RGB(0.498,0.498,0.498),
    "PostVehicle" => RGB(0.498,0.498,0.498),
    "Saline" => RGB(0.498,0.498,0.498),
    "None" => RGB(0.498,0.498,0.498),
    )
const protocol_colors = Dict(
    "High" => RGB(248/255,214/255,21/255),
    "Med." => RGB(242/255,151/255,37/255),
    "Low" => RGB(221/255,43/255,41/255),
    "1.0" => RGB(248/255,214/255,21/255),
    "0.75" => RGB(242/255,151/255,37/255),
    "0.5" => RGB(221/255,43/255,41/255),
    )

const color_series = [
    RGB(0.6,0.149,0.561),
    RGB(0.902,0.604,0.737),
    RGB(0.1216,0.4667,0.7059),
    RGB(0.1725,0.6275,0.1725),
    RGB(1.0,0.498,0.0549),
]

const Treatment_ord = [
    "None",
    "Saline",
    "PostVehicle",
    "PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram"
]

const Treatment_dict = Dict(x=>i for (i,x) in enumerate(Treatment_ord))

const Plotting_position = Dict("Optogenetic" => 1,
                                "Citalopram" => 2,
                                "Methysergide" => 3,
                                "SB242084" => 4,
                                "Altanserin" => 5,
                                "Way_100135" => 6,
                                "SB242084_opt" => 7
                                )


text:
using AudioIO

# Give PortAudio time to load
play([0])
sleep(2)

println("""
                    *
                *   *
            *   *   *
        *   *   *   *
    *   *   *   *   *
*   *   *   *   *   *
""")
wave = SinOsc(440) * LinRamp(0.0, 1.0, 2.0)
play(wave)
sleep(2)
stop(wave)

println("""
                    *
                *   *   *
            *   *   *   *   *
        *   *   *   *   *   *   *
    *   *   *   *   *   *   *   *   *
*   *   *   *   *   *   *   *   *   *   *
""")
wave = SinOsc(440) * LinRamp([0.0, 1.0, 0.0], [2.0, 2.0])
play(wave)
sleep(4)
stop(wave)