text:
# Function that loads and clean the data into an edgelist
function load_data(df, s)
    # Separate based on WebScale
    scale = [df[df.WebScale .== u,:] for u in unique(df.WebScale)]
    # select for edge interaction columns only
    edgelist = scale[s][:, [:PredTSN, :PreyTSN]]
    # view all entries that contain "san" in PredTSN
    pred_san = edgelist[startswith.(edgelist.PredTSN, "san"),:]
    # remove "san"
    pred_san.PredTSN = replace.(pred_san.PredTSN, "san" .=> "")
    pred_san
    #view unique predID without san
    unique(pred_san,1)
    edgelist.PredTSN == "102" #manually check if in original df
    # remove the "san" from the dataframe itself
    edgelist.PredTSN = replace.(edgelist.PredTSN, "san" .=> "")
    # Do the same for PreyTSN
    backup2 = copy(edgelist)
    old_ids2 = unique(backup2.PreyTSN)
    # keep the indexes before the change
    prey_indexes = startswith.(edgelist.PreyTSN, "san")
    # remove the "san" from the dataframe itself
    edgelist.PreyTSN = replace.(edgelist.PreyTSN, "san" .=> "")
    edgelist
    #check unique preyIDs without san
    check_uniques2 = [in(new_id2, old_ids2) for new_id2 in edgelist.PreyTSN[prey_indexes]]
    filter(isone, check_uniques2)
    # convert "string" to "integer"
    edgelist.PredTSN = parse.(Int64, edgelist.PredTSN)
    edgelist.PreyTSN = parse.(Int64, edgelist.PreyTSN)

    return edgelist
end

# Function that rename the IDs
function renameIDs(edgelist)
    ##Change PredTSN to newID starting from 1
    # Sort by ID
    sort!(edgelist, :PredTSN)
    # Get unique IDs
    unique_ids = unique(edgelist.PredTSN)
    #create new IDs
    newid_pred = collect(1:length(unique_ids))
    # Create Dict where id => newid
    id_dict = Dict(unique_ids[i] => newid_pred[i] for i in 1:length(unique_ids))
    # Add new ids in dataframe
    edgelist.newPredID = [id_dict[id] for id in edgelist.PredTSN]
    ### Do same for prey
    ## Change PreyTSN to new ID starting from end of newPredID
    # Sort by ID
    sort!(edgelist, :PreyTSN)
    # Get unique IDs
    unique_ids_prey = unique(edgelist.PreyTSN)
    # Create new IDs
    newid_prey = collect(length(id_dict)+1:length(id_dict)+length(unique_ids_prey))
    # Create Dict
    id_dict2 = Dict(unique_ids_prey[i] => newid_prey[i] for i in 1:length(unique_ids_prey))
    # Add new ids in dataframe
    edgelist.newPreyID = [id_dict2[id] for id in edgelist.PreyTSN]
    keys = edgelist[:, [:newPredID, :newPreyID]]
    #sort
    sort!(keys, :newPredID)
    #get unique values
    unique!(keys)

    return keys
end

# Function that creates the adjacency matrix
function adjacencyMatrix(keys)
    all_ids = vcat(unique(keys.newPredID), unique(keys.newPreyID))
    ##Create Adjacency matrix
    #create matrix of zeros with n = colPred, m = colPrey
    M = zeros(Int64, length(all_ids)+1, length(all_ids)+1)
    #assign names for cols
    for y in 2:length(all_ids)+1
        M[y,1] = all_ids[y-1]
    end
    #assign names for rows
    for x in 2:length(all_ids)+1
        M[1,x] = all_ids[x-1]
    end
    # Fill adjacency matrix with ones where there are interactions
    for i in 1:length(all_ids)
        M[keys[i, 1]+1, keys[i, 2]+1] = 1
    end
    return M
end

# Function that computes the modularity
function computeModularity(M)
    # Removing firts row and first column of M matrix to keep only values
    M_wrow = M[setdiff(1:end, 1), :]
    M_wcol = M_wrow[:, setdiff(1:end, 1)]
    # Convert ones and zeros in boolean
    M_bool = convert(Array{Bool}, M_wcol .== 1)
    # Convert matrix in bipartite network
    B = UnipartiteNetwork(M_bool)
    # Calculate modularity
    modules = brim(lp(B)...)
    return Q(modules...), Qr(modules...), B
end

# Function that cumputes the nestedness
function computeNestedness(B)
    # Calculate nestedness
    return ρ(B)
end

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ScenarioBased <: Forecast
        label::String
        scenario_count::Int64
        data::TimeSeries.TimeArray
    end

A Discrete Scenario Based forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `scenario_count::Int64`: Number of scenarios
- `data::TimeSeries.TimeArray`: timestamp - scalingfactor
"""
mutable struct ScenarioBased <: Forecast
    "label of component parameter forecasted"
    label::String
    "Number of scenarios"
    scenario_count::Int64
    "timestamp - scalingfactor"
    data::TimeSeries.TimeArray
end

function ScenarioBased(; label, scenario_count, data, )
    ScenarioBased(label, scenario_count, data, )
end

"""Get ScenarioBased label."""
get_label(value::ScenarioBased) = value.label
"""Get ScenarioBased scenario_count."""
get_scenario_count(value::ScenarioBased) = value.scenario_count
"""Get ScenarioBased data."""
get_data(value::ScenarioBased) = value.data

ext: jl


text:
# ### Load packages

using GeoPhyInv
using Statistics
using Plots

model=GIPh.Gallery.Seismic(:acou_homo1); # load a homogeneous model
GIPh.Models.Seismic_addon!(model, randn_perc=0.01); # add some random noise

# a simple acquisition geometry

acqgeom=GeoPhyInv.Gallery.Geom(model.mgrid,:xwell);

# plot the model and source, receivers
p1=JP.seismic(model)
JP.geom!(acqgeom)
plot(p1)

# generate time grid

tgrid = range(0.0,stop=2.0,length=1000)

# Ricker wavelet

wav = GeoPhyInv.Utils.Wavelets.ricker(10.0, tgrid, tpeak=0.25,);

# distribute the same source wavelet to all the supsersources

acqsrc=GeoPhyInv.Acquisition.Src_fixed(acqgeom.nss,1,[:P],wav,tgrid);

# create `Fdtd.Param` object to prepare forward modelling
# * npw corresponds to the number of independently propagating wavefields (1 in most cases)
# Once the `Param` object is created, do the modelling "without any memory allocations" using `mod!`

pa=GeoPhyInv.Fdtd.Param(npw=1,model=model,
	acqgeom=[acqgeom], acqsrc=[acqsrc],
	sflags=[2], rflags=[1],
	tgridmod=tgrid, verbose=true);

@time GeoPhyInv.Fdtd.mod!(pa);

# plot a record after modelling
pdata=plot(pa.c.data[1].d[1,1])
plot(pdata)

# create new seismic model

model_new=J.Gallery.Seismic(:acou_homo1) # prepare another model
J.Models.Seismic_addon!(model_new, randn_perc=0.01)
J.Models.Seismic_addon!(model_new, constant_pert=0.03) # perturb the model
p2=JP.seismic(model_new) # plot new model
JP.geom!(acqgeom)
plot(p2)

# Now, we the change the model in the `Param` object without memory allocation
# This routine can be used during FWI,
# where medium parameters are itertively updated in the same `Fdtd.Param` object

J.Fdtd.update_model!(pa.c, model_new)

# run modelling now and plot data again
@time GeoPhyInv.Fdtd.mod!(pa);

# plot a record after modelling
plot!(pdata, pa.c.data[1].d[1,1])
plot(pdata)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# test postprocessing of strain and stress

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0])

u = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [-1/3, 0.0, 0.0],
    3 => [-1/3, -1/3, 0.0],
    4 => [0.0, -1/3, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [-1/3, 0.0, 1.0],
    7 => [-1/3, -1/3, 1.0],
    8 => [0.0, -1/3, 1.0])

element = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
update!(element, "geometry", 0.0 => X)
update!(element, "displacement", 0.0 => u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)

body = Problem(Elasticity, "[0,1]³ elastic unit block", 3)
add_elements!(body, element)
postprocess!(body, 0.0, Val{:strain})
postprocess!(body, 0.0, Val{:stress})

geom = element("geometry", 0.0)
strain = element("strain", 0.0)
stress = element("stress", 0.0)

for i in 1:length(geom)
    @test isapprox(strain[i], [-1/3, -1/3, 1.0, 0.0, 0.0, 0.0])
    @test isapprox(stress[i], [0.0, 0.0, 288.0, 0.0, 0.0, 0.0])
end


text:
"""
```
    fit_ellipse(observations::AbstractObservations, f::AbstractFittingMethod)
```

Fits an ellipse to a set of scattered `observations` using a specifed ellipse fitting method `f`.

# Details

If each data point is accompanied by a 2x2 covariance matrix then `observations` is an instance of
[UncertainObservations](@ref), otherwise it is an instance of `[Observations](@ref)`

If `f` is `GuaranteedEllipseFit()` and `observations <: UncertainObservations` then
the supplied covariance matrices are incorporated into the estimation procedure.
If `f` is `GuaranteedEllipseFit()` and `observations <: Observations` (data point covariance matrices
are not provided), then the algorithm assumes homogeneous Gaussian noise and estimates the
level of noise from the data.

If `f` is `DirectEllipseFit()` then data point covariance matrices are ignored even if they
are supplied because the algroithm cannot incorporate such information.
"""
function fit_ellipse(observations::AbstractObservations, method::AbstractFittingMethod)
    normalisation_method = IsotropicScalingTranslation()
    normalise = NormaliseDataContext(observations, normalisation_method)
    # Estimate ellipse in a data-driven normalised coordinate system.
    𝛉′ =  method(normalise(observations))
    # Transform parameters back to the original coordinate system.
    𝛉 = normalise(FromNormalisedSpace(), 𝛉′)
    return 𝛉 / norm(𝛉)
end

function fit_ellipse(data::AbstractMatrix, method::AbstractFittingMethod)
    # Convert the N x 2 matrix to a one-dimensional array of 2D points.
    ℳ = svectors(transpose(data), Val{2}())
    fit_ellipse(Observations(tuple(ℳ)), method)
end

function fit_ellipse(observations::AbstractObservations, method::GuaranteedEllipseFit)
    @unpack optimisation_scheme = method
    @unpack seed = optimisation_scheme
    normalisation_method = IsotropicScalingTranslation()
    normalise = NormaliseDataContext(observations, normalisation_method)
    # Convert observations and seed to a data-driven normalised coordinate system.
    method = @set method.optimisation_scheme.seed = ManualEstimation(normalise(ToNormalisedSpace(), seed.𝛉))
    # Apply the guaranteed ellipse fit.
    𝛉′ =  method(normalise(observations))
    # Transform parameters back to the original coordinate system.
    𝛉 = normalise(FromNormalisedSpace(), 𝛉′)
    return 𝛉 / norm(𝛉)
end


text:
import Oceananigans.BoundaryConditions:
    fill_halo_regions!, fill_top_halo!, fill_bottom_halo!, fill_west_halo!, fill_east_halo!, fill_south_halo!, fill_north_halo!

import Oceananigans.Models.HydrostaticFreeSurfaceModels: fill_horizontal_velocity_halos!

# These filling functions won't work so let's not use them.
 fill_west_halo!(c, bc::CubedSphereExchangeBC, args...; kwargs...) = nothing
 fill_east_halo!(c, bc::CubedSphereExchangeBC, args...; kwargs...) = nothing
fill_south_halo!(c, bc::CubedSphereExchangeBC, args...; kwargs...) = nothing
fill_north_halo!(c, bc::CubedSphereExchangeBC, args...; kwargs...) = nothing

function fill_halo_regions!(field::AbstractCubedSphereField, arch, args...; kwargs...)

    for (i, face_field) in enumerate(faces(field))
        # Fill the top and bottom halos the usual way.
        fill_halo_regions!(face_field, arch, get_face(args, i)...; kwargs...)

        # Deal with halo exchanges.
        fill_west_halo!(face_field, field)
        fill_east_halo!(face_field, field)
        fill_south_halo!(face_field, field)
        fill_north_halo!(face_field, field)
    end

    return nothing
end

function fill_west_halo!(field::CubedSphereFaceField{LX, LY, LZ}, cubed_sphere_field) where {LX, LY, LZ}
    location = (LX, LY, LZ)
    dest_halo = underlying_west_halo(field.data, field.grid, LX)

    exchange_info = field.boundary_conditions.west.condition
    src_face_number = exchange_info.to_face
    src_side = exchange_info.to_side
    src_boundary = cubed_sphere_boundary(cubed_sphere_field, location, src_face_number, src_side)

    if sides_in_the_same_dimension(:west, src_side)
        dest_halo .= src_boundary
    else
        dest_halo .= reverse(permutedims(src_boundary, (2, 1, 3)), dims=2)
    end

    return nothing
end

function fill_east_halo!(field::CubedSphereFaceField{LX, LY, LZ}, cubed_sphere_field) where {LX, LY, LZ}
    location = (LX, LY, LZ)
    dest_halo = underlying_east_halo(field.data, field.grid, LX)

    exchange_info = field.boundary_conditions.east.condition
    src_face_number = exchange_info.to_face
    src_side = exchange_info.to_side
    src_boundary = cubed_sphere_boundary(cubed_sphere_field, location, src_face_number, src_side)

    if sides_in_the_same_dimension(:east, src_side)
        dest_halo .= src_boundary
    else
        dest_halo .= reverse(permutedims(src_boundary, (2, 1, 3)), dims=2)
    end

    return nothing
end

function fill_south_halo!(field::CubedSphereFaceField{LX, LY, LZ}, cubed_sphere_field) where {LX, LY, LZ}
    location = (LX, LY, LZ)
    dest_halo = underlying_south_halo(field.data, field.grid, LY)

    exchange_info = field.boundary_conditions.south.condition
    src_face_number = exchange_info.to_face
    src_side = exchange_info.to_side
    src_boundary = cubed_sphere_boundary(cubed_sphere_field, location, src_face_number, src_side)

    if sides_in_the_same_dimension(:south, src_side)
        dest_halo .= src_boundary
    else
        dest_halo .= reverse(permutedims(src_boundary, (2, 1, 3)), dims=1)
    end

    return nothing
end

function fill_north_halo!(field::CubedSphereFaceField{LX, LY, LZ}, cubed_sphere_field) where {LX, LY, LZ}
    location = (LX, LY, LZ)
    dest_halo = underlying_north_halo(field.data, field.grid, LY)

    exchange_info = field.boundary_conditions.north.condition
    src_face_number = exchange_info.to_face
    src_side = exchange_info.to_side
    src_boundary = cubed_sphere_boundary(cubed_sphere_field, location, src_face_number, src_side)

    if sides_in_the_same_dimension(:north, src_side)
        dest_halo .= src_boundary
    else
        dest_halo .= reverse(permutedims(src_boundary, (2, 1, 3)), dims=1)
    end

    return nothing
end

# Don't worry about this when not on a cubed sphere.
fill_horizontal_velocity_halos!(u, v, arch) = nothing

function fill_horizontal_velocity_halos!(u::CubedSphereField, v::CubedSphereField, arch)

    ## Fill them like they're tracers to get the top and bottom filled.
    fill_halo_regions!(u, arch)
    fill_halo_regions!(v, arch)

    ## Now fill the horizontal halos.

    u_loc = (Face, Center, Center)
    v_loc = (Center, Face, Center)

    for face_number in 1:6, side in (:west, :east, :south, :north)
        exchange_info = getproperty(u.boundary_conditions.faces[face_number], side).condition
        src_face_number = exchange_info.to_face
        src_side = exchange_info.to_side

        if sides_in_the_same_dimension(side, src_side)
            cubed_sphere_halo(u, u_loc, face_number, side) .= cubed_sphere_boundary(u, u_loc, src_face_number, src_side)
            cubed_sphere_halo(v, v_loc, face_number, side) .= cubed_sphere_boundary(v, v_loc, src_face_number, src_side)
        else
            u_sign = (isodd(face_number) && side == :west ) || (iseven(face_number) && side == :east ) ? +1 : -1
            v_sign = (isodd(face_number) && side == :north) || (iseven(face_number) && side == :south) ? +1 : -1

            reverse_dim = src_side in (:west, :east) ? 1 : 2

            cubed_sphere_halo(u, u_loc, face_number, side) .= u_sign * reverse(permutedims(cubed_sphere_boundary(v, v_loc, src_face_number, src_side), (2, 1, 3)), dims=reverse_dim)
            cubed_sphere_halo(v, v_loc, face_number, side) .= v_sign * reverse(permutedims(cubed_sphere_boundary(u, u_loc, src_face_number, src_side), (2, 1, 3)), dims=reverse_dim)
        end
    end

    return nothing
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

VOID = [[]], [[1]]

# example 1

V = [[0,0] [1,0] [2,0] [0,1] [1,1] [2,1] [0,2] [1,2] [2,2]];
FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
pattern = repeat([1,2,-3],outer=4);
model = (V,FV);
W,FW = Lar.extrudeSimplicial(model, pattern);
GL.VIEW([ GL.GLGrid(W,FW, GL.COLORS[7], 0.5) ]);

# example 2

model = Lar.extrudeSimplicial( VOID, ones(10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]); #TODO: check

model = Lar.extrudeSimplicial( model, ones(10) )
GL.VIEW([ GL.GLLines(model..., GL.COLORS[7]),
        GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1)) ]);

model = Lar.extrudeSimplicial( model, ones(10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]);

# example 3

model = Lar.extrudeSimplicial( VOID, repeat([1,-1],outer=10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]); #TODO: check

model = Lar.extrudeSimplicial( model, repeat([1,-1],outer=10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7],0.7),
        GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1)) ]);

# example 4

grid_2d = Lar.simplexGrid([3,3])
GL.VIEW([ GL.GLGrid(grid_2d..., GL.COLORS[7],0.7),
        GL.GLLines(grid_2d...),    # TODO:  check lines
        GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1)) ]);

grid_3d = Lar.simplexGrid([2,3,4])
GL.VIEW([ GL.GLGrid(grid_3d..., GL.COLORS[7]),
        GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1)) ]);  # TODO  check


text:
# reimplements methods: hessfact, hessfact!
# and minor functions from LinAlg: full, A_mul_B! ..., reflectorApply!

import LinearAlgebra
import Base: A_mul_B!, Ac_mul_B!, A_mul_Bc!, copymutable, full
import LinearAlgebra: hessfact, hessfact!, Hessenberg, HessenbergQ
import LinearAlgebra: chkstride1, checksquare

import LinearAlgebra: reflector!, reflectorApply!
import LinearAlgebra.LAPACK: gehrd!

function Hessenberg{T<:BigFloatOrComplex}(A::StridedMatrix{T})
  Hessenberg(gehrd!(A)...)
end

hessfact!{T<:BigFloatOrComplex}(A::StridedMatrix{T}) = Hessenberg(A)
hessfact{T<:BigFloatOrComplex}(A::StridedMatrix{T}) = hessfact!(copy(A))

# Perfom Householder/Hessenberg
# The Householder matrices are unitary but not hermitian in the complex case
# to provide real subdiagonal of the Hessenberg matrix.
# mimic the corresponding LAPACK.gehrd! function.
function gehrd!{T<:BigFloatOrComplex}(ilo::Integer, ihi::Integer, A::StridedMatrix{T})
  chkstride1(A)
  n = checksquare(A)
  τ = zeros(T, max(0, n - 1))

  ilo = max(ilo, 1)
  ihi = min(ihi, n)

  for k = ilo:ihi-1
    x = view(A, k+1:ihi, k)
    τk = reflector!(x)
    τ[k] = τk
    reflectorApply!(x, τk, view(A, k+1:n, k+1:n))
    reflectorApply!(view(A, :, k+1:n), x, τk)
  end
  A, τ
end

"""

  `reflectorApply!(x::AbstractVector, τ::Number, A::StridedMatrix)`
  `reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number)`

Apply reflector `x, τ` from left or right to matrix `A`, see `linalg/generic.jl`.
Matrix is modified by this operation.
"""
function reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number)
  n, m = size(A)
  if length(x) != m
    throw(DimensionMismatch("reflector has length $(length(x)), which must match the second d    imension of matrix A, $m"))
  end
  @inbounds begin
  for j = 1:n
    # dot
    vAj = A[j,1]
    for i = 2:m
      vAj += x[i]*A[j,i]
    end

    vAj = τ*vAj

    # ger
    A[j, 1] -= vAj
    for i = 2:m
      A[j, i] -= x[i]'*vAj
    end
   end
  end
  return A
end

# various multiplications with HessenbergQ
# note the adjoint(τ[k]) - supports the case of complex τ!
function A_mul_B!{T<:BigFloatOrComplex}(HQ::HessenbergQ{T}, A::StridedVecOrMat{T})
  n = size(A, 1)
  τ = HQ.τ
  for k = length(τ):-1:1
    τk = τ[k]'
    if τk != 0
      x = view(HQ.factors, k+1:n, k)
      reflectorApply!(x, τk, view(A, k+1:n, :))
    end
  end
  A
end

function A_mul_B!{T<:BigFloatOrComplex}(A::StridedMatrix{T}, HQ::HessenbergQ{T})
  n = size(A, 2)
  τ = HQ.τ
  for k = 1:length(τ)
    τk = τ[k]
    if τk != 0
      x = view(HQ.factors, k+1:n, k)
      reflectorApply!(view(A, :, k+1:n), x, τk)
    end
  end
  A
end

function Ac_mul_B!{T<:BigFloatOrComplex}(HQ::HessenbergQ{T}, A::StridedVecOrMat{T})
  n = size(A, 1)
  τ = HQ.τ
  for k = 1:length(τ)
    τk = τ[k]
    if τk != 0
      x = view(HQ.factors, k+1:n, k)
      reflectorApply!(x, τk, view(A, k+1:n, :))
    end
  end
  A
end

function A_mul_Bc!{T<:BigFloatOrComplex}(A::StridedMatrix{T}, HQ::HessenbergQ{T})
  n = size(A, 2)
  τ = HQ.τ
  for k = length(τ):-1:1
    τk = τ[k]'
    if τk != 0
      x = view(HQ.factors, k+1:n, k)
      reflectorApply!(view(A, :, k+1:n), x, τk)
    end
  end
  A
end

function copymutable(HQ::HessenbergQ)
  Q = eye(HQ.factors)
  A_mul_B!(HQ, Q)
  Q
end

full{T<:BigFloatOrComplex,S<:AbstractMatrix}(HQ::HessenbergQ{T,S}) = copymutable(HQ)

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
using StatisticalRethinking, Optim
#gr(size=(600,600));

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(3, p), 3) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)

N = 10000
samples = sample(p_grid, Weights(posterior), N);

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

MCMCChains.describe(chn)

MCMCChains.hpd(chn, alpha=0.5)

mode(samples)

mean(samples)

median(samples)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module Particles

using OffsetArrays

export Equation, Simulation, run!
export ExplicitEuler

include("equation.jl")
include("integrators.jl")
include("simulation.jl")

end # module

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
@testset "Constructing empty sequences" begin
    @test LongDNASeq() == LongSequence(DNA)
    @test LongRNASeq() == LongSequence(RNA)
    @test LongAminoAcidSeq() == LongSequence(AminoAcid)
    @test LongCharSeq() == LongSequence(Char)
end

@testset "Constructing uninitialized sequences" begin
    @test isa(LongSequence{DNAAlphabet{2}}(0), LongSequence)
    @test isa(LongSequence{DNAAlphabet{4}}(10), LongSequence)
    @test isa(LongSequence{RNAAlphabet{2}}(0), LongSequence)
    @test isa(LongSequence{RNAAlphabet{4}}(10), LongSequence)
    @test isa(LongSequence{AminoAcidAlphabet}(10), LongSequence)

    @test_throws ArgumentError LongSequence{DNAAlphabet{2}}(-1)
    @test_throws ArgumentError LongSequence{DNAAlphabet{4}}(-1)
    @test_throws ArgumentError LongSequence{RNAAlphabet{2}}(-1)
    @test_throws ArgumentError LongSequence{RNAAlphabet{4}}(-1)
    @test_throws ArgumentError LongSequence{AminoAcidAlphabet}(-1)
end

@testset "Conversion from/to strings" begin
    # Check that sequences in strings survive round trip conversion:
    #   String → LongSequence → String
    function test_string_construction(A::Type, seq::AbstractString)
        @test convert(String, LongSequence{A}(seq)) == uppercase(seq)
    end

    function test_string_parse(A::Type, seq::AbstractString)
        @test parse(LongSequence{A}, seq) == LongSequence{A}(seq)
    end

    for len in [0, 1, 2, 3, 10, 32, 1000, 10000]
        test_string_construction(DNAAlphabet{4}, random_dna(len))
        test_string_construction(DNAAlphabet{4}, SubString(random_dna(len), 1:len))
        test_string_construction(DNAAlphabet{4}, lowercase(random_dna(len)))
        test_string_construction(RNAAlphabet{4}, lowercase(random_rna(len)))
        test_string_construction(RNAAlphabet{4}, random_rna(len))
        test_string_construction(AminoAcidAlphabet, random_aa(len))
        test_string_construction(AminoAcidAlphabet, lowercase(random_aa(len)))

        test_string_parse(DNAAlphabet{4}, random_dna(len))
        test_string_parse(DNAAlphabet{4}, SubString(random_dna(len), 1:len))
        test_string_parse(DNAAlphabet{4}, lowercase(random_dna(len)))
        test_string_parse(RNAAlphabet{4}, lowercase(random_rna(len)))
        test_string_parse(RNAAlphabet{4}, random_rna(len))
        test_string_parse(AminoAcidAlphabet, random_aa(len))
        test_string_parse(AminoAcidAlphabet, lowercase(random_aa(len)))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_string_construction(DNAAlphabet{2}, random_dna(len, probs))
        test_string_construction(DNAAlphabet{2}, lowercase(random_dna(len, probs)))
        test_string_construction(RNAAlphabet{2}, random_rna(len, probs))
        test_string_construction(RNAAlphabet{2}, lowercase(random_rna(len, probs)))

        test_string_parse(DNAAlphabet{2}, random_dna(len, probs))
        test_string_parse(DNAAlphabet{2}, lowercase(random_dna(len, probs)))
        test_string_parse(RNAAlphabet{2}, random_rna(len, probs))
        test_string_parse(RNAAlphabet{2}, lowercase(random_rna(len, probs)))
    end

    # non-standard string literal
    @test isa(dna"ACGTMRWSYKVHDBN-", LongDNASeq)
    @test isa(rna"ACGUMRWSYKVHDBN-", LongRNASeq)
    @test isa(aa"ARNDCQEGHILKMFPSTWYVBJZXOU*-", LongAminoAcidSeq)
    @test isa(char"いろは αβγ 甲乙丙", LongCharSeq)

    # Non-nucleotide characters should throw
    @test_throws Exception LongDNASeq("ACCNNCATTTTTTAGATXATAG")
    @test_throws Exception LongRNASeq("ACCNNCATTTTTTAGATXATAG")
    @test_throws Exception LongAminoAcidSeq("ATGHLMY@ZACAGNM")
end

@testset "Construction from vectors" begin
    function test_vector_construction(A, seq::AbstractString)
        T = eltype(A)
        xs = T[convert(T, c) for c in seq]
        @test LongSequence{A}(xs) == LongSequence{A}(seq)
    end

    for len in [0, 1, 10, 32, 1000, 10000]
        test_vector_construction(DNAAlphabet{4}, random_dna(len))
        test_vector_construction(RNAAlphabet{4}, random_rna(len))
        test_vector_construction(AminoAcidAlphabet, random_aa(len))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_vector_construction(DNAAlphabet{2}, random_dna(len, probs))
        test_vector_construction(RNAAlphabet{2}, random_rna(len, probs))
    end
end

@testset "Encode_copy!" begin
    # Note: Other packages use this function, so we need to test it
    # Even though this is NOT exported or part of the API in a normal sense
    function test_copyto!(dst::LongSequence, doff, src, soff, N)
        BioSequences.copyto!(dst, doff, src, soff, N)
        @test String(dst[doff:doff+N-1]) == String(src[soff:soff+N-1])
    end

    probs = [0.25, 0.25, 0.25, 0.25]
    for len in [0, 1, 2, 10, 100]
        for f in [identity, Vector{Char}, Vector{UInt8}]
            test_copyto!(LongSequence{DNAAlphabet{2}}(len), 1, f(random_dna(len, probs)), 1, len)
            test_copyto!(LongSequence{RNAAlphabet{2}}(len), 1, f(random_rna(len, probs)), 1, len)
            test_copyto!(LongSequence{DNAAlphabet{4}}(len), 1, f(random_dna(len)), 1, len)
            test_copyto!(LongSequence{RNAAlphabet{4}}(len), 1, f(random_rna(len)), 1, len)
            test_copyto!(LongSequence{AminoAcidAlphabet}(len), 1, f(random_aa(len)), 1, len)
            test_copyto!(LongSequence{CharAlphabet}(len), 1, f(random_aa(len)), 1, len)
        end
    end

    for len in [10, 32, 100]
        for f in [identity, Vector{Char}, Vector{UInt8}]
            test_copyto!(LongSequence{DNAAlphabet{2}}(len+7), 5, f(random_dna(len+11, probs)), 3, len)
            test_copyto!(LongSequence{RNAAlphabet{2}}(len+7), 5, f(random_rna(len+11, probs)), 3, len)
            test_copyto!(LongSequence{DNAAlphabet{4}}(len+7), 5, f(random_dna(len+11)), 3, len)
            test_copyto!(LongSequence{RNAAlphabet{4}}(len+7), 5, f(random_rna(len+11)), 3, len)
            test_copyto!(LongSequence{AminoAcidAlphabet}(len+7), 5, f(random_aa(len+11)), 3, len)
            test_copyto!(LongSequence{CharAlphabet}(len+7), 5, f(random_aa(len+11)), 3, len)
        end
    end
end

@testset "Conversion between 2-bit and 4-bit encodings" begin
    function test_conversion(A1, A2, seq)
        @test convert(LongSequence{A1}, LongSequence{A2}(seq)) == LongSequence{A1}(seq)
    end

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "")
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "")
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "")
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "")

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "ACGT")
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "ACGT")
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "ACGU")
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "ACGU")

    test_conversion(DNAAlphabet{2}, DNAAlphabet{4}, "ACGT"^100)
    test_conversion(DNAAlphabet{4}, DNAAlphabet{2}, "ACGT"^100)
    test_conversion(RNAAlphabet{4}, RNAAlphabet{2}, "ACGU"^100)
    test_conversion(RNAAlphabet{2}, RNAAlphabet{4}, "ACGU"^100)

    # ambiguous nucleotides cannot be stored in 2-bit encoding
    EncodeError = BioSequences.EncodeError
    @test_throws EncodeError convert(LongSequence{DNAAlphabet{2}}, dna"AN")
    @test_throws EncodeError convert(LongSequence{RNAAlphabet{2}}, rna"AN")

    # test promotion
    a = LongSequence{DNAAlphabet{2}}("ATCG")
    b = LongSequence{DNAAlphabet{4}}("ATCG")
    c = LongSequence{RNAAlphabet{2}}("AUCG")
    d = LongSequence{RNAAlphabet{4}}("AUCG")

    @test typeof(promote(a, b)) == Tuple{LongSequence{DNAAlphabet{4}},LongSequence{DNAAlphabet{4}}}
    @test typeof(promote(c, d)) == Tuple{LongSequence{RNAAlphabet{4}},LongSequence{RNAAlphabet{4}}}
    @test_throws ErrorException typeof(promote(a, d))
    @test_throws ErrorException typeof(promote(a, b, d))
end

@testset "Conversion between RNA and DNA" begin
    @test convert(LongRNASeq, LongDNASeq("ACGTN")) == rna"ACGUN"
    @test convert(LongDNASeq, LongRNASeq("ACGUN")) == dna"ACGTN"
end

@testset "Conversion to Matrices" begin
    dna = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
    dnathrow = [dna"AAA", dna"TTTAAA", dna"CCC", dna"GGG"]

    rna = [rna"AAA", rna"UUU", rna"CCC", rna"GGG"]
    rnathrow = [rna"AAA", rna"UUU", rna"CCCUUU", rna"GGG"]

    prot = [aa"AMG", aa"AMG", aa"AMG", aa"AMG"]

    sitemajdna = [
        DNA_A  DNA_A  DNA_A
        DNA_T  DNA_T  DNA_T
        DNA_C  DNA_C  DNA_C
        DNA_G  DNA_G  DNA_G
    ]
    seqmajdna = [
        DNA_A  DNA_T  DNA_C  DNA_G
        DNA_A  DNA_T  DNA_C  DNA_G
        DNA_A  DNA_T  DNA_C  DNA_G
    ]
    sitemajrna = [
        RNA_A  RNA_A  RNA_A
        RNA_U  RNA_U  RNA_U
        RNA_C  RNA_C  RNA_C
        RNA_G  RNA_G  RNA_G
    ]
    seqmajrna = [
        RNA_A  RNA_U  RNA_C  RNA_G
        RNA_A  RNA_U  RNA_C  RNA_G
        RNA_A  RNA_U  RNA_C  RNA_G
    ]
    sitemajnucint = [
        0x01  0x01  0x01
        0x08  0x08  0x08
        0x02  0x02  0x02
        0x04  0x04  0x04
    ]
    seqmajnucint = [
        0x01  0x08  0x02  0x04
        0x01  0x08  0x02  0x04
        0x01  0x08  0x02  0x04
    ]
    sitemajaa = [
        AA_A AA_M AA_G
        AA_A AA_M AA_G
        AA_A AA_M AA_G
        AA_A AA_M AA_G
    ]
    seqmajaa = [
        AA_A AA_A AA_A AA_A
        AA_M AA_M AA_M AA_M
        AA_G AA_G AA_G AA_G
    ]

    @test seqmatrix(dna, :site) == sitemajdna
    @test seqmatrix(rna, :site) == sitemajrna
    @test seqmatrix(prot, :site) == sitemajaa
    @test seqmatrix(UInt8, dna, :site) == sitemajnucint
    @test seqmatrix(UInt8, rna, :site) == sitemajnucint

    @test seqmatrix(dna, :seq) == seqmajdna
    @test seqmatrix(rna, :seq) == seqmajrna
    @test seqmatrix(prot, :seq) == seqmajaa
    @test seqmatrix(UInt8, dna, :seq) == seqmajnucint
    @test seqmatrix(UInt8, rna, :seq) == seqmajnucint

    @test seqmatrix([dna"", dna"", dna""], :site) == Matrix{DNA}(undef, (3, 0))
    @test seqmatrix([dna"", dna"", dna""], :seq) == Matrix{DNA}(undef, (0, 3))
    @test seqmatrix([rna"", rna"", rna""], :site) == Matrix{RNA}(undef, (3, 0))
    @test seqmatrix([rna"", rna"", rna""], :seq) == Matrix{RNA}(undef, (0, 3))
    @test seqmatrix(UInt8, [dna"", dna"", dna""], :site) == Matrix{UInt8}(undef, (3, 0))
    @test seqmatrix(UInt8, [dna"", dna"", dna""], :seq) == Matrix{UInt8}(undef, (0, 3))
    @test seqmatrix(UInt8, [rna"", rna"", rna""], :site) == Matrix{UInt8}(undef, (3, 0))
    @test seqmatrix(UInt8, [rna"", rna"", rna""], :seq) == Matrix{UInt8}(undef, (0, 3))

    @test_throws ArgumentError seqmatrix(dnathrow, :site)
    @test_throws ArgumentError seqmatrix(rnathrow, :seq)
    @test_throws ArgumentError seqmatrix(dna, :lol)
    @test_throws MethodError seqmatrix(AminoAcid, dna, :site)
    @test_throws ArgumentError seqmatrix(LongDNASeq[], :site)
    @test_throws ArgumentError seqmatrix(LongDNASeq[], :seq)

end


text:
# DynamicalBilliards.jl constants
SV = DynamicalBilliards.SV
Obstacle, Billiard, AbstractParticle = DynamicalBilliards.Obstacle, DynamicalBilliards.Billiard, DynamicalBilliards.AbstractParticle
using AbstractPlotting: RGBf0, RGBAf0

obcolor(::Obstacle) = RGBf0(0,0.6,0)
obcolor(::Union{DynamicalBilliards.RandomWall, DynamicalBilliards.RandomDisk}) = RGBf0(149/255, 88/255, 178/255)
obcolor(::Union{DynamicalBilliards.SplitterWall, DynamicalBilliards.Antidot, DynamicalBilliards.Ellipse}) = RGBf0(0.8,0.0,0)
obcolor(::DynamicalBilliards.PeriodicWall) = RGBf0(0.8,0.8,0)
obfill(o::DynamicalBilliards.Obstacle) = RGBAf0(obcolor(o).r, obcolor(o).g, obcolor(o).b, 0.5)
obfill(o::Union{DynamicalBilliards.Antidot, DynamicalBilliards.Ellipse}) = RGBAf0(obcolor(o), 0.1)
obls(::Obstacle) = nothing
obls(::Union{DynamicalBilliards.SplitterWall, DynamicalBilliards.Antidot, DynamicalBilliards.Ellipse}) = :dash
obls(::DynamicalBilliards.PeriodicWall) = :dotted
oblw(::Obstacle) = 2.0

function bdplot!(ax, o::T; kwargs...) where {T}
    error("Element of type $T does not have a plotting definition yet.")
end

function bdplot!(ax, o::DynamicalBilliards.Semicircle; kwargs...)
    θ1 = atan(o.facedir[2], o.facedir[1]) + π/2 # start of semicircle
    θ2 = θ1 + π
    θ = range(θ1, θ2; length = 200)
    p = [Point2f0(cos(t)*o.r + o.c[1], sin(t)*o.r + o.c[2]) for t in θ]
    lines!(ax, p; color = obcolor(o), linewidth = oblw(o), linestyle = obls(o),
    scale_plot=false, kwargs...)
end

function bdplot!(ax, w::DynamicalBilliards.Wall; kwargs...)
    lines!(ax, Float32[w.sp[1],w.ep[1]], Float32[w.sp[2],w.ep[2]];
    color = obcolor(w), linewidth = oblw(w), scale_plot=false, kwargs...)
end

function bdplot!(ax, o::DynamicalBilliards.Circular; kwargs...)
    θ = 0:0.01:2π
    p = [Point2f0(cos(t)*o.r + o.c[1], sin(t)*o.r + o.c[2]) for t in θ]
    poly!(ax, p; color = obfill(o), scale_plot=false, kwargs...)
    lines!(ax, p; color = obcolor(o), linewidth = oblw(o), linestyle = obls(o),
    scale_plot=false, kwargs...)
end

function bdplot!(ax, bd::Billiard; kwargs...)
    xmin, ymin, xmax, ymax = DynamicalBilliards.cellsize(bd)
    dx = xmax - xmin; dy = ymax - ymin
    for obst in bd; bdplot!(ax, obst; show_axis = false, kwargs...); end
    if !isinf(xmin) && !isinf(xmax)
        AbstractPlotting.xlims!(ax, xmin - 0.1dx, xmax + 0.1dx)
    end
    if !isinf(ymin) && !isinf(ymax)
        AbstractPlotting.ylims!(ax, ymin - 0.1dy, ymax + 0.1dy)
    end
    remove_axis!(ax)
    return ax
end

function remove_axis!(ax)
    ax.bottomspinevisible = false
    ax.leftspinevisible = false
    ax.topspinevisible = false
    ax.rightspinevisible = false
    ax.xgridvisible = false
    ax.ygridvisible = false
    ax.xticksvisible = false
    ax.yticksvisible = false
    ax.xticklabelsvisible = false
    ax.yticklabelsvisible = false
end

bdplot!(ax, p::AbstractParticle; kwargs...) = bdplot!(ax, [p]; kwargs...)

function _estimate_vr(bd)
    xmin, ymin, xmax, ymax = DynamicalBilliards.cellsize(bd)
    f = max((xmax-xmin), (ymax-ymin))
    isinf(f) && error("cellsize of billiard is infinite")
    vr = Float32(f/25)
end

function bdplot!(ax, bd, ps::Vector{<:AbstractParticle};
    use_cell = true, kwargs...)
    c = haskey(kwargs, :color) ? kwargs[:color] : AbstractPlotting.RGBf0(0,0,0)
    N = length(ps)
    xs, ys = Observable(zeros(Float32, N)), Observable(zeros(Float32, N))
    vxs, vys = Observable(zeros(Float32, N)), Observable(zeros(Float32, N))

    # need heuristic for ms and vr
    ms = 6
    vr = _estimate_vr(bd)
    for i in 1:N
        p = ps[i]
        pos = use_cell ? p.pos + p.current_cell : p.pos
        xs[][i] = pos[1]
        ys[][i] = pos[2]
        θ = atan(p.vel[2], p.vel[1])
        vxs[][i] = vr*cos(θ)
        vys[][i] = vr*sin(θ)
    end
    # important: marker hack for zoom-independent size. Will change in the future
    # to allow something like `markerspace = :display`
    marker = Circle(Point2f0(xs[][1], ys[][1]), Float32(1))

    scatter!(ax, xs, ys; color = c, marker = marker, markersize = ms*px)
    arrows!(ax, xs, ys, vxs, vys;
        normalize = false, arrowsize = 0.01px,
        linewidth  = 2,
    )
    return xs, ys
end


text:
using LowRankModels,  Random
import StatsBase: sample, Weights, std
import LinearAlgebra: norm

# test quadratic loss

## generate data
Random.seed!(1);
m,n,k = 300,300,3;
kfit = k+1
# variance of measurement
sigmasq = .1

# coordinates of covariates
X_real = randn(m,k)
# directions of observations
Y_real = randn(k,n)

XY = X_real*Y_real;
A = XY + sqrt(sigmasq)*randn(m,n)

# and the model
losses = QuadLoss()
rx, ry = QuadReg(.1), QuadReg(.1);
glrm = GLRM(A,losses,rx,ry,kfit)
#scale=false, offset=false, X=randn(kfit,m), Y=randn(kfit,n));

# fit w/o initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
println("After fitting, parameters differ from true parameters by $(norm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE\n")

# initialize
init_svd!(glrm)
XYh = glrm.X' * glrm.Y
println("After initialization with the svd, parameters differ from true parameters by $(norm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE\n")

# fit w/ initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
println("After fitting, parameters differ from true parameters by $(norm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")


text:
function beta(obj::GNSSAttitudeContext, arg0::AbsoluteDate)
    return jcall(obj, "beta", jdouble, (AbsoluteDate,), arg0)
end

function beta_ds(obj::GNSSAttitudeContext)
    return jcall(obj, "betaDS", DerivativeStructure, ())
end

function compute_phi(obj::GNSSAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "computePhi", jdouble, (jdouble, jdouble), arg0, arg1)
end

function get_date(obj::GNSSAttitudeContext)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_delta_ds(obj::GNSSAttitudeContext)
    return jcall(obj, "getDeltaDS", DerivativeStructure, ())
end

function get_mu_rate(obj::GNSSAttitudeContext)
    return jcall(obj, "getMuRate", jdouble, ())
end

function get_orbit_angle_since_midnight(obj::GNSSAttitudeContext)
    return jcall(obj, "getOrbitAngleSinceMidnight", jdouble, ())
end

function get_secured_beta(obj::GNSSAttitudeContext)
    return jcall(obj, "getSecuredBeta", jdouble, ())
end

function get_sv_bcos(obj::GNSSAttitudeContext)
    return jcall(obj, "getSVBcos", jdouble, ())
end

function get_turn_span(obj::GNSSAttitudeContext)
    return jcall(obj, "getTurnSpan", TurnSpan, ())
end

function get_yaw_end(obj::GNSSAttitudeContext, arg0::jdouble)
    return jcall(obj, "getYawEnd", jdouble, (jdouble,), arg0)
end

function get_yaw_start(obj::GNSSAttitudeContext, arg0::jdouble)
    return jcall(obj, "getYawStart", jdouble, (jdouble,), arg0)
end

function in_orbit_plane_absolute_angle(obj::GNSSAttitudeContext, arg0::jdouble)
    return jcall(obj, "inOrbitPlaneAbsoluteAngle", jdouble, (jdouble,), arg0)
end

function in_sun_side(obj::GNSSAttitudeContext)
    return jcall(obj, "inSunSide", jboolean, ())
end

function in_turn_time_range(obj::GNSSAttitudeContext)
    return jcall(obj, "inTurnTimeRange", jboolean, ())
end

function linear_model_still_active(obj::GNSSAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "linearModelStillActive", jboolean, (jdouble, jdouble), arg0, arg1)
end

function nominal_yaw(obj::GNSSAttitudeContext, arg0::AbsoluteDate)
    return jcall(obj, "nominalYaw", TimeStampedAngularCoordinates, (AbsoluteDate,), arg0)
end

function orbit_normal_yaw(obj::GNSSAttitudeContext)
    return jcall(obj, "orbitNormalYaw", TimeStampedAngularCoordinates, ())
end

function set_half_span(obj::GNSSAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "setHalfSpan", void, (jdouble, jdouble), arg0, arg1)
end

function set_up_turn_region(obj::GNSSAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "setUpTurnRegion", jboolean, (jdouble, jdouble), arg0, arg1)
end

function time_since_turn_start(obj::GNSSAttitudeContext)
    return jcall(obj, "timeSinceTurnStart", jdouble, ())
end

function turn_corrected_attitude(obj::GNSSAttitudeContext, arg0::DerivativeStructure)
    return jcall(obj, "turnCorrectedAttitude", TimeStampedAngularCoordinates, (DerivativeStructure,), arg0)
end

function turn_corrected_attitude(obj::GNSSAttitudeContext, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "turnCorrectedAttitude", TimeStampedAngularCoordinates, (jdouble, jdouble), arg0, arg1)
end

function yaw_rate(obj::GNSSAttitudeContext, arg0::jdouble)
    return jcall(obj, "yawRate", jdouble, (jdouble,), arg0)
end


text:
# this script basially only handles `ARGS`

Base.push!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
using VSCodeServer
pop!(LOAD_PATH)

let
    args = [popfirst!(Base.ARGS) for _ in 1:5]
    # load Revise ?
    if "USE_REVISE=true" in args
        try
            @static if VERSION ≥ v"1.5"
                using Revise
            else
                @eval using Revise
                Revise.async_steal_repl_backend()
            end
        catch err
        end
    end

    atreplinit() do repl
        "USE_PLOTPANE=true" in args && Base.Multimedia.pushdisplay(VSCodeServer.InlineDisplay())
    end

    conn_pipeline, telemetry_pipeline = args[1:2]
    VSCodeServer.serve(conn_pipeline; is_dev = "DEBUG_MODE=true" in args, crashreporting_pipename = telemetry_pipeline)
end

max_stars_repo_name: fredrikekre/julia-vscode

max_issues_repo_name: fredrikekre/julia-vscode

max_forks_repo_name: fredrikekre/julia-vscode


text:
using DrWatson
using FileIO, BSON, DataFrames
using PrettyTables
using Statistics
using StatsBase
import PGFPlots

using GenerativeAD.Evaluation: MODEL_ALIAS, DATASET_ALIAS, MODEL_TYPE, apply_aliases!
using GenerativeAD.Evaluation: _prefix_symbol, aggregate_stats_mean_max, aggregate_stats_max_mean
using GenerativeAD.Evaluation: PAT_METRICS, PATN_METRICS, PAC_METRICS, BASE_METRICS, TRAIN_EVAL_TIMES
using GenerativeAD.Evaluation: rank_table, print_rank_table, latex_booktabs, convert_anomaly_class

include("./utils/pgf_boxplot.jl")
include("./utils/ranks.jl")
include("./utils/bayes.jl")

const PAT_METRICS_NAMES = ["\$PR@\\%0.01\$","\$PR@\\%0.1\$","\$PR@\\%1\$","\$PR@\\%5\$","\$PR@\\%10\$","\$PR@\\%20\$"]
const PAC_METRICS_NAMES = ["\$AUC@\\#5\$","\$AUC@\\#10\$","\$AUC@\\#50\$","\$AUC@\\#100\$","\$AUC@\\#500\$","\$AUC@\\#1000\$"]
const PATN_METRICS_NAMES = ["\$PR@\\#5\$","\$PR@\\#10\$","\$PR@\\#50\$","\$PR@\\#100\$","\$PR@\\#500\$","\$PR@\\#1000\$"]

DOWNSAMPLE = Dict(
    zip(["ocsvm", "MAF", "RealNVP"], [100, 100, 100]))

function _tabular_filter!(df)
    # filter out ocsvm_nu
    filter!(x -> (x.modelname != "ocsvm_nu"), df)

    # filter autoencoders others than "aae_full", "wae_full" "vae_full"
    filter!(x -> ~(x.modelname in (["aae", "aae_vamp", "vae", "vae_simple", "wae", "wae_vamp"])), df)

    # filter only default parameters of "ocsvm" - rbf + nu = 0.5, tune gamma
    filter!(x -> (x.modelname != "ocsvm_rbf") || occursin("nu=0.5", x.parameters), df)
    df
end

function _filter_autoencoders!(df)
    autoencoders = Set([k for (k,v) in MODEL_TYPE if v == "autoencoders"])
    filter!(x -> x.modelname in autoencoders, df)
end

function _filter_ensembles!(df)
    filter!(x -> occursin("ignore_nan=true_method=mean", x.parameters), df)
end

df_tabular = load(datadir("evaluation/tabular_eval.bson"))[:df];
df_tabular_clean = load(datadir("evaluation/tabular_clean_val_final_eval.bson"))[:df];
df_tabular_autoencoders = _filter_autoencoders!(copy(df_tabular));
# with autoencoders separate, it is now safe to just filter out models that do not get into the big comparison
_tabular_filter!(df_tabular);

df_tabular_bayes = load(datadir("evaluation_bayes/tabular_eval.bson"))[:df];
df_tabular_bayes_outerjoin = combine_bayes(df_tabular, df_tabular_bayes; outer=true);
bayes_models = Set(unique(df_tabular_bayes.modelname))
df_tabular_bayes_innerjoin = filter(x -> x.modelname in bayes_models, df_tabular_bayes_outerjoin);

df_tabular[:anomaly_class] = -1;
df_tabular_clean[:anomaly_class] = 1; # in order to trigger maxmean during autoagg

df_tabular_ens = load(datadir("evaluation_ensembles/tabular_eval.bson"))[:df];
_filter_ensembles!(df_tabular_ens);

@info "Loaded results from tabular evaluation."

# works for both image and tabular data
function basic_summary_table(df; suffix="", prefix="", downsample=Dict{String, Int}())
    agg_names = ["maxmean", "meanmax"]
    agg_funct = [aggregate_stats_max_mean, aggregate_stats_mean_max]
    for (name, agg) in zip(agg_names, agg_funct)
        for metric in [:auc, :tpr_5]
            val_metric = _prefix_symbol("val", metric)
            tst_metric = _prefix_symbol("tst", metric)

            _, rt = sorted_rank(df, agg, val_metric, tst_metric, downsample)

            rt[end-2, 1] = "\$\\sigma_1\$"
            rt[end-1, 1] = "\$\\sigma_{10}\$"
            rt[end, 1] = "rnk"

            file = "$(projectdir())/paper/tables/$(prefix)_$(metric)_$(metric)_$(name)$(suffix).tex"
            open(file, "w") do io
                print_rank_table(io, rt; backend=:tex)
            end
        end
    end
end

## one table with `ocsvm_rbf` and the other with properly tuned `ocsvm`
basic_summary_table(filter(x->x.modelname != "ocsvm_rbf", df_tabular); prefix="tabular", suffix="", downsample=DOWNSAMPLE)
basic_summary_table(filter(x->x.modelname != "ocsvm", df_tabular); prefix="tabular", suffix="_orbf", downsample=DOWNSAMPLE)

## default/clean parameters table
basic_summary_table(filter(x->x.modelname != "ocsvm_rbf", df_tabular_clean); prefix="tabular", suffix="_clean_default")

## bayes results combined with initial random samples and other methods that were not optimized that way
basic_summary_table(filter(x->x.modelname != "ocsvm_rbf", df_tabular_bayes_outerjoin); prefix="tabular", suffix="_bayes_outerjoin")

## bayes results combined with initial random samples
basic_summary_table(df_tabular_bayes_innerjoin; prefix="tabular", suffix="_bayes_innerjoin")

basic_summary_table(df_tabular_ens; prefix="tabular", suffix="_ensembles")
@info "basic_summary_table_tabular"

# how the choice of creation criterion and size of an ensembles affect results
function ensemble_sensitivity(df, df_ensemble; prefix="tabular", suffix="", downsample=Dict{String, Int}())
    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        ranks = []
        # computing baseline
        models, rt = sorted_rank(df, aggregate_stats_max_mean, val_metric, tst_metric, downsample)
        rt["criterion-size"] = "baseline"
        select!(rt, vcat(["criterion-size"], models))
        push!(ranks, rt[end:end, :])

        # computing different ensembles
        for (cn, c) in zip(["AUC", "TPR@5"], [:val_auc, :val_tpr_5])
            for s in [5, 10]
                dff = filter(x -> startswith(x.parameters, "criterion=$(c)") && endswith(x.parameters, "_size=$(s)"), df_ensemble)
                models, rt_ensemble = sorted_rank(dff, aggregate_stats_max_mean, val_metric, tst_metric, downsample; verbose=false)
                rt_ensemble["criterion-size"] = "$(cn)-$(s)"
                select!(rt_ensemble, vcat(["criterion-size"], models))
                mean_dif = DataFrame([Symbol(model) => round.(mean(rt_ensemble[1:end-3, model] .- rt[1:end-3, model]), digits=2) for model in models])
                mean_dif["criterion-size"] = "$(cn)-$(s)"

                push!(ranks, rt_ensemble[end:end, :])
                push!(ranks, mean_dif)
            end
        end

        df_ranks = reduce(vcat, ranks)
        hl_best_rank = LatexHighlighter(
                        (data, i, j) -> ((i == 1) || (i%2 == 0)) && (data[i,j] == minimum(df_ranks[i, 2:end])),
                        ["color{red}","textbf"])

        hl_best_dif = LatexHighlighter(
                        (data, i, j) -> (i != 1) && (i%2 == 1) && (data[i,j] == maximum(df_ranks[i, 2:end])),
                        ["color{blue}","textbf"])

        f_float = (v, i, j) -> (i != 1) && (i%2 == 1) ? ft_printf("%.2f")(v,i,j) : ft_printf("%.1f")(v,i,j)

        file = "$(projectdir())/paper/tables/$(prefix)_ensemble_size_$(metric)$(suffix).tex"
        open(file, "w") do io
            pretty_table(
                io, df_ranks,
                backend=:latex,
                formatters=f_float,
                highlighters=(hl_best_rank, hl_best_dif),
                nosubheader=true,
                tf=latex_booktabs)
        end
    end
end

ensemble_sensitivity(filter(x->x.modelname != "ocsvm_rbf", df_tabular), df_tabular_ens; prefix="tabular", downsample=DOWNSAMPLE)
@info "ensemble_sensitivity_tabular"

# only meanmax as we don't have bayes for images
function bayes_sensitivity(df, df_bayes; prefix="tabular", suffix="", downsample=Dict{String, Int}())
    ranks = []

    # differences in performance in both metrics
    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        # computing baseline
        models, rt = sorted_rank(df, aggregate_stats_mean_max, val_metric, tst_metric, downsample)
        rt["exp"] = "random-$(mn)"
        select!(rt, vcat(["exp"], models))
        push!(ranks, rt[end:end, :])

        models, rt_bayes = sorted_rank(df_bayes, aggregate_stats_mean_max, val_metric, tst_metric)
        rt_bayes["exp"] = "bayes-$(mn)"
        select!(rt_bayes, vcat(["exp"], models))

        metric_dif = Matrix(rt_bayes[1:end-3, 2:end]) - Matrix(rt[1:end-3, 2:end])
        num_inc = sum(metric_dif .> 0, dims=1)
        num_dec = sum(metric_dif .< 0, dims=1)
        mean_dif = round.(mean(metric_dif, dims=1), digits=2)

        rank_dif = Vector(rt_bayes[end, 2:end]) - Vector(rt[end, 2:end])

        push!(rt_bayes, ["rank. change", rank_dif...])
        push!(rt_bayes, ["avg. change", mean_dif...])
        push!(rt_bayes, ["num. inc", num_inc...])
        push!(rt_bayes, ["num. dec", num_dec...])
        push!(ranks, rt_bayes[end-4:end, :])
    end
    df_ranks = reduce(vcat, ranks)
    # first six rows correspond to AUC and the rest to TPR
    # 1. random ranks          :: minimum
    # 2. bayes  ranks          :: minimum
    # 3. bayes - random ranks  :: minimum
    # 4. bayes - random metric :: maximum
    # 5. bayes > random metric :: maximum
    # 6. bayes < random metric :: maximum

    min_rows = [1,2,3,7,8,9]
    max_rows = setdiff(collect(1:12), min_rows)
    hl_best_rank = LatexHighlighter(
                    (data, i, j) -> (i in min_rows) && (data[i,j] == minimum(df_ranks[i, 2:end])),
                    ["color{red}","textbf"])

    hl_best_dif = LatexHighlighter(
                    (data, i, j) -> (i in max_rows) && (data[i,j] == maximum(df_ranks[i, 2:end])),
                    ["color{blue}","textbf"])

    f_float = (v, i, j) -> (i == 4 || i == 10) ? ft_printf("%.2f")(v,i,j) : (i in [5,6,11,12]) ? ft_printf("%.0f")(v,i,j) : ft_printf("%.1f")(v,i,j)

    file = "$(projectdir())/paper/tables/$(prefix)_bayes_comp$(suffix).tex"
    open(file, "w") do io
        pretty_table(
            io, df_ranks,
            backend=:latex,
            formatters=f_float,
            highlighters=(hl_best_rank, hl_best_dif),
            nosubheader=true,
            tf=latex_booktabs)
    end
end

bayes_sensitivity(
    filter(x->x.modelname != "ocsvm_rbf", df_tabular),
    filter(x->x.modelname != "ocsvm_rbf", df_tabular_bayes_outerjoin); prefix="tabular", suffix="", downsample=DOWNSAMPLE)
@info "bayes_sensitivity_tabular"

# training time rank vs avg rank
# should work for all dataset types as it uses the automatic aggregation
function plot_fiteval_time(df; time_col=:fit_t, suffix="", prefix="tabular", format="pdf", downsample=Dict{String,Int}())
    # add first stage time for encoding and training of encoding
    df["total_fit_t"] = df["fit_t"] .+ df["fs_fit_t"]
    # add all eval times together
    df["total_eval_t"] = df["tr_eval_t"] .+ df["tst_eval_t"] .+ df["val_eval_t"] .+ df["fs_eval_t"]

    # compute ranks from averaged values over each dataset
    # should reduce the noise in the results
    agg_cols = vcat(Symbol.(TRAIN_EVAL_TIMES), [:total_eval_t, :total_fit_t])
    df_time_avg = combine(groupby(df, [:dataset, :modelname]), agg_cols .=> mean .=> agg_cols)
    df_time_avg[_prefix_symbol(time_col, "top_10_std")] = 0.0   # add dummy column
    df_time_avg[_prefix_symbol(time_col, "std")] = 0.0          # add dummy column
    df_time_avg[time_col] .= -df_time_avg[time_col]

    # applying this to a copy of a dataframe is prefered
    apply_aliases!(df_time_avg, col="modelname", d=MODEL_RENAME)
	apply_aliases!(df_time_avg, col="modelname", d=MODEL_ALIAS)
	apply_aliases!(df_time_avg, col="dataset", d=DATASET_ALIAS)
    # this is to ensure how the rank_table is sorted
	sort!(df_time_avg, (:dataset, :modelname))

    rtt = rank_table(df_time_avg, time_col)
    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        df_agg = aggregate_stats_auto(df, val_metric; downsample=downsample)
        apply_aliases!(df_agg, col="modelname", d=MODEL_RENAME)
        apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
        apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)

        # this is to ensure how the rank_table is sorted
        sort!(df_agg, (:dataset, :modelname))

        rt = rank_table(df_agg, tst_metric)

        models = names(rt)[2:end]
        x = Vector(rt[end, 2:end])
        y = Vector(rtt[end, 2:end])
        # labels cannot be shifted uniformly
        a = PGFPlots.Axis([
                PGFPlots.Plots.Scatter(x, y),
                [PGFPlots.Plots.Node(m, xx + 0.5, yy + 0.7) for (m, xx, yy) in zip(models, x, y)]...],
                xlabel="avg. rnk",
                ylabel="avg. \$$(String(time_col))\$ rnk")
        file = "$(projectdir())/paper/figures/$(prefix)_$(time_col)_vs_$(metric)$(suffix).$(format)"
        PGFPlots.save(file, a; include_preamble=false)
    end
end

plot_fiteval_time(filter(x->x.modelname != "ocsvm_rbf", df_tabular); time_col=:total_fit_t, format="tex", downsample=DOWNSAMPLE)
plot_fiteval_time(filter(x->x.modelname != "ocsvm_rbf", df_tabular); time_col=:total_eval_t, format="tex", downsample=DOWNSAMPLE)
@info "plot_fiteval_time_tabular"

# show basic table only for autoencoders
function tabular_autoencoders_investigation(df; split_vamp=false, suffix="")
    apply_aliases!(df, col="modelname", d=MODEL_ALIAS)
    apply_aliases!(df, col="dataset", d=DATASET_ALIAS)

    # split aae_full and wae_full based on used prior
    if split_vamp
        aaefv_mask = (df.modelname .== "aaef") .& occursin.("prior=vamp", df.parameters)
        df[aaefv_mask, :modelname] .=  "aaefv"
        waefv_mask = (df.modelname .== "waef") .& occursin.("prior=vamp", df.parameters)
        df[waefv_mask, :modelname] .=  "waefv"
    else
        apply_aliases!(df, col="modelname", d=Dict("aaev" => "aae", "waev" => "wae"))
    end

    # define further splitting of models based on parameters
    jc_mask = occursin.("jacodeco", df.parameters)
    df[jc_mask, :modelname] .=  df[jc_mask, :modelname] .*"-jc"

    lm_mask = occursin.("latent-mean", df.parameters)
    df[lm_mask, :modelname] .=  "disregard"

    l_mask = occursin.("latent_", df.parameters)
    df[l_mask, :modelname] .=  "disregard"

    ls_mask = occursin.("latent-sampled", df.parameters)
    df[ls_mask, :modelname] .=  "disregard"

    rm_mask = occursin.("reconstruction-mean", df.parameters)
    df[rm_mask, :modelname] .=  df[rm_mask, :modelname] .*"-rm"

    r_mask = occursin.("reconstruction_", df.parameters)
    df[r_mask, :modelname] .=  "disregard"

    rs_mask = occursin.("reconstruction-sampled", df.parameters)
    df[rs_mask, :modelname] .=  df[rs_mask, :modelname] .*"-rs"

    d_mask = occursin.("disc_", df.parameters)
    df[d_mask, :modelname] .=  df[d_mask, :modelname] .*"-d"

    el_mask = occursin.("elbo_", df.parameters)
    df[el_mask, :modelname] .=  df[el_mask, :modelname] .*"-el"

    mse_mask = occursin.("score=mse_", df.parameters)
    df[mse_mask, :modelname] .=  df[mse_mask, :modelname] .*"-rm"

    filter!(x -> (x.modelname != "disregard"), df)

    for metric in [:auc, :tpr_5]
        for (name, agg) in zip(
                    ["maxmean", "meanmax"],
                    [aggregate_stats_max_mean, aggregate_stats_mean_max])
            val_metric = _prefix_symbol("val", metric)
            tst_metric = _prefix_symbol("tst", metric)

            df_agg = agg(df, val_metric)
            sort!(df_agg, (:dataset, :modelname))
            rt = rank_table(df_agg, tst_metric)

            rt[end-2, 1] = "\$\\sigma_1\$"
            rt[end-1, 1] = "\$\\sigma_{10}\$"
            rt[end, 1] = "rnk"

            file = "$(projectdir())/paper/tables/tabular_ae_only_$(metric)_$(metric)_$(name)$(suffix).tex"
            open(file, "w") do io
                print_rank_table(io, rt; backend=:tex)
            end

            # due to its width we opted for box plot
            # recompute ranks for each dataset
            ranks = compute_ranks(rt[1:end-3, 2:end])
            models = names(rt)[2:end]

            #### grouping by method
            mgroups = ones(Int, length(models))
            mgroups[startswith.(models, "avae")] .= 1
            mgroups[startswith.(models, "gano")] .= 2
            mgroups[startswith.(models, "aae-")] .= 3
            mgroups[startswith.(models, "aaev")] .= 4
            mgroups[startswith.(models, "aaef-")] .= 5
            mgroups[startswith.(models, "aaefv")] .= 6
            mgroups[startswith.(models, "vae-")] .= 7
            mgroups[startswith.(models, "vaef")] .= 8
            mgroups[startswith.(models, "vaes")] .= 9
            mgroups[startswith.(models, "wae-")] .= 10
            mgroups[startswith.(models, "waev")] .= 11
            mgroups[startswith.(models, "waef-")] .= 12
            mgroups[startswith.(models, "waefv")] .= 13

            sm = sortperm(mgroups)
            mgroups = mgroups[sm]

            models = models[sm]
            ranks = ranks[sm, :]

            #### grouping by score, 5 types + 1 for other scores
            groups = 6*ones(Int, length(models))
            groups[endswith.(models, "-rm")] .= 1
            groups[endswith.(models, "-rs")] .= 2
            groups[endswith.(models, "-jc")] .= 3
            groups[endswith.(models, "-d")] .= 4
            groups[endswith.(models, "-el")] .= 5
            ####

            # compute statistics for boxplot
            rmin, rmax = maximum(ranks, dims=2), maximum(ranks, dims=2)
            rmean = mean(ranks, dims=2)
            rmedian = median(ranks, dims=2)
            rlowq = quantile.(eachrow(ranks), 0.25)
            rhighq = quantile.(eachrow(ranks), 0.75)

            a = pgf_boxplot_grouped_colorpos(rlowq, rhighq, rmedian, rmin, rmax, models, groups, mgroups; h="12cm", w="8cm")
            file = "$(projectdir())/paper/figures/tabular_ae_only_box_$(metric)_$(name)$(suffix).tex"
            open(file, "w") do f
                write(f, a)
            end
        end
    end
end

tabular_autoencoders_investigation(df_tabular_autoencoders; split_vamp=true, suffix="_split_vamp")
@info "tabular_autoencoders_investigation"

# does crossvalidation matters?
function per_seed_ranks_tabular(df; suffix="", downsample=Dict{String, Int}())
    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        ranks = []
        for (name, agg) in zip(
                ["maxmean", "meanmax"],
                [aggregate_stats_max_mean, aggregate_stats_mean_max])
            models, rt = sorted_rank(df, agg, val_metric, tst_metric, downsample)
            rt["exp"] = name
            select!(rt, vcat(["exp"], models))
            push!(ranks, rt[end:end,:])
        end

        for seed in 1:5
            dff = filter(x -> (x.seed == seed), df)
            # silence the warnings for insufficient number of seeds by setting verbose false
            models, rt = sorted_rank(dff, aggregate_stats_max_mean, val_metric, tst_metric, downsample; verbose=false)
            rt["exp"] = "seed=$seed"
            select!(rt, vcat(["exp"], models))
            push!(ranks, rt[end:end,:])
        end

        df_ranks = reduce(vcat, ranks)
        hl_best_rank = LatexHighlighter(
                        (data, i, j) -> (data[i,j] == minimum(df_ranks[i, 2:end])),
                        ["color{red}","textbf"])

        file = "$(projectdir())/paper/tables/tabular_per_seed_ranks_$(metric)$(suffix).tex"
        open(file, "w") do io
            pretty_table(
                io, df_ranks,
                backend=:latex,
                formatters=ft_printf("%.1f"),
                highlighters=(hl_best_rank),
                nosubheader=true,
                tf=latex_booktabs)
        end
    end
end

per_seed_ranks_tabular(filter(x->x.modelname != "ocsvm_rbf", df_tabular), suffix="", downsample=DOWNSAMPLE)
per_seed_ranks_tabular(filter(x->x.modelname != "ocsvm", df_tabular), suffix="_orbf", downsample=DOWNSAMPLE)
@info "per_seed_ranks_tabular"

# metric sensitivity
# should be both for image and tabular data
function metric_sensitivity(df; prefix="tabular", suffix="", downsample=Dict{String, Int}())
    ranks = []

    # compute rank table for each metric
    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        models, rt = sorted_rank(df, aggregate_stats_auto, val_metric, tst_metric, downsample)
        rt["metric"] = "$(mn)"
        select!(rt, vcat(["metric"], models))
        push!(ranks, rt[end:end, :])
    end

    # each row contains average ranks for each metric
    df_ranks = reduce(vcat, ranks)
    rank_dif = Vector(df_ranks[2, 2:end]) .- Vector(df_ranks[1, 2:end])
    push!(df_ranks, ["rank. change", rank_dif...])

    min_rows = [1,2,3]
    max_rows = [3]
    hl_best_rank = LatexHighlighter(
                    (data, i, j) -> (i in min_rows) && (data[i,j] == minimum(df_ranks[i, 2:end])),
                    ["color{red}","textbf"])

    hl_best_dif = LatexHighlighter(
                    (data, i, j) -> (i in max_rows) && (data[i,j] == maximum(df_ranks[i, 2:end])),
                    ["color{blue}","textbf"])

    file = "$(projectdir())/paper/tables/$(prefix)_metric_comp$(suffix).tex"
    open(file, "w") do io
        pretty_table(
            io, df_ranks,
            backend=:latex,
            formatters=ft_printf("%.1f"),
            highlighters=(hl_best_rank, hl_best_dif),
            nosubheader=true,
            tf=latex_booktabs)
    end
end

# needs dummy anomaly_class column
metric_sensitivity(filter(x->x.modelname != "ocsvm_rbf", df_tabular), suffix="", downsample=DOWNSAMPLE)
metric_sensitivity(filter(x->x.modelname != "ocsvm_rbf", df_tabular_clean), suffix="_clean_default")
@info "metric_sensitivity"

######################################################################################
#######################               IMAGES                ##########################
######################################################################################
# df_images_loo = load(datadir("evaluation/images_leave-one-out_eval.bson"))[:df];
# df_images_loo_clean = load(datadir("evaluation/images_leave-one-out_clean_val_final_eval.bson"))[:df];
# filter!(x -> x.seed == 1, df_images_loo)
# filter!(x -> ~(x.modelname in ["aae_ocsvm", "fAnoGAN-GP"]), df_images_loo)
df_images_loi = load(datadir("evaluation/images_leave-one-in_eval.bson"))[:df];
df_images_loi_clean = load(datadir("evaluation/images_leave-one-in_clean_val_final_eval.bson"))[:df];
df_images_loi_clean[df_images_loi_clean.fit_t .=== nothing, :fit_t] .= 1.0; # ConvSkipGANomaly is missing the fit_t

df_images_mnistc = load(datadir("evaluation/images_mnistc_eval.bson"))[:df];
filter!(x -> x.dataset != "MNIST-C_zigzag", df_images_mnistc); # filter out MNISTC-C_zigzag which we forgot and could not compute in time
df_images_mvtec = load(datadir("evaluation/images_mvtec_eval.bson"))[:df];
select!(df_images_mnistc, Not(:anomaly_class)); # mnistc contains anomaly_class but mvtec does not
df_images_single = vcat(df_images_mnistc, df_images_mvtec);
df_images_single[:anomaly_class] = -1; # add dummy anomaly_class to indicate that we have single class anomaly dataset

df_images_mnistc_clean = load(datadir("evaluation/images_mnistc_clean_val_final_eval.bson"))[:df];
filter!(x -> x.dataset != "MNIST-C_zigzag", df_images_mnistc_clean); # filter out MNISTC-C_zigzag which we forgot and could not compute in time
df_images_mvtec_clean = load(datadir("evaluation/images_mvtec_clean_val_final_eval.bson"))[:df];

select!(df_images_mnistc_clean, Not(:anomaly_class));
df_images_single_clean = vcat(df_images_mnistc_clean, df_images_mvtec_clean);
# df_images_single_clean[:anomaly_class] = -1; # this is needed when concatenating with multiclass datasets
df_images_single_clean[:anomaly_class] = 1; # this is needed when concatenating with multiclass datasets
@info "Loaded image results"

df_images_loi_ens = load(datadir("evaluation_ensembles/images_leave-one-in_eval.bson"))[:df];
df_images_mnistc_ens = load(datadir("evaluation_ensembles/images_mnistc_eval.bson"))[:df];
df_images_mvtec_ens = load(datadir("evaluation_ensembles/images_mvtec_eval.bson"))[:df];
df_images_mvtec_ens[:anomaly_class] = -1;
_filter_ensembles!(df_images_loi_ens);
_filter_ensembles!(df_images_mnistc_ens);
_filter_ensembles!(df_images_mvtec_ens);

df_images = vcat(df_images_loi, df_images_single);
df_images_clean = vcat(df_images_loi_clean, df_images_single_clean);
df_images_ens = vcat(df_images_loi_ens, df_images_mnistc_ens, df_images_mvtec_ens);

# renaming dataset to dataset:ac
for df in [df_images, df_images_clean, df_images_ens]
    apply_aliases!(df, col="dataset", d=DATASET_ALIAS)
    for d in Set(["cifar10", "svhn2", "mnist", "fmnist"])
        mask = (df.dataset .== d)
        df[mask, :dataset] .= df[mask, :dataset] .* ":" .* convert_anomaly_class.(df[mask, :anomaly_class], d)
        df[mask, :anomaly_class] .= 1 # it has to be > 0, because otherwise we get too many warnings from the aggregate_stats_max_mean
    end
end

@info "Loaded ensemble image results"
#= do we have all that we need ?
df_images_loi.modelname |> countmap
df_images_mnistc.modelname |> countmap
df_images_mvtec.modelname |> countmap

df_images_loi_clean.modelname |> countmap
df_images_mnistc_clean.modelname |> countmap
df_images_mvtec_clean.modelname |> countmap

df_images_loi.dataset |> countmap
df_images_mnistc.dataset |> countmap
df_images_mvtec.dataset |> countmap

df_images_loi_clean.dataset |> countmap
df_images_mnistc_clean.dataset |> countmap
df_images_mvtec_clean.dataset |> countmap
=#

basic_summary_table(df_images_loi, prefix="images_multi", suffix="")
basic_summary_table(df_images_loi_clean, prefix="images_multi", suffix="_clean_default")
basic_summary_table(df_images_single; prefix="images_single", suffix="")
basic_summary_table(df_images_single_clean; prefix="images_single", suffix="_clean_default")
@info "basic_summary_table_images separated"

SEMANTIC_IMAGE_ANOMALIES = Set(["CIFAR10", "SVHN2"])

# splits single and multi class image datasets into "statistic" and "semantic" anomalies
_split_image_datasets(df) = (
            filter(x -> x.dataset in SEMANTIC_IMAGE_ANOMALIES, df),
            filter(x -> ~(x.dataset in SEMANTIC_IMAGE_ANOMALIES), df)
        )

df_images_semantic, df_images_stat = _split_image_datasets(df_images);
df_images_semantic_clean, df_images_stat_clean =  _split_image_datasets(df_images_clean);
df_images_semantic_ens, df_images_stat_ens = _split_image_datasets(df_images_ens);

basic_summary_table_autoagg(df_images_semantic, prefix="images_semantic", suffix="_per_ac")
basic_summary_table_autoagg(df_images_stat, prefix="images_stat", suffix="_per_ac")
@info "basic_summary_table_images stat/semantic"

basic_summary_table(df_images_semantic_clean, prefix="images_semantic", suffix="_clean_default_per_ac")
basic_summary_table(df_images_stat_clean, prefix="images_stat", suffix="_clean_default_per_ac")
@info "basic_summary_table_images stat/semantic clean"

metric_sensitivity(df_images_semantic, prefix="images_semantic", suffix="_per_ac")
metric_sensitivity(df_images_stat, prefix="images_stat", suffix="_per_ac")

df_images_stat_clean[:anomaly_class] = 1;       # in order to trigger maxmean agg
df_images_semantic_clean[:anomaly_class] = 1;   # in order to trigger maxmean agg
metric_sensitivity(df_images_semantic_clean, prefix="images_semantic", suffix="_clean_default_per_ac")
metric_sensitivity(df_images_stat_clean, prefix="images_stat", suffix="_clean_default_per_ac")
@info "metric_sensitivity"

ensemble_sensitivity(df_images_semantic, df_images_semantic_ens; prefix="images_semantic", suffix="_per_ac")
ensemble_sensitivity(df_images_stat, df_images_stat_ens; prefix="images_stat", suffix="_per_ac")
@info "ensemble_sensitivity"

plot_fiteval_time(df_images; prefix="images", time_col=:total_fit_t, format="tex", suffix="_per_ac")
plot_fiteval_time(df_images; prefix="images", time_col=:total_eval_t, format="tex", suffix="_per_ac")
@info "plot_fiteval_time_images"

# works for both image and tabular data assuming that tabular has anomaly_class column = -1
function basic_summary_table_autoagg(df; suffix="", prefix="", downsample=Dict{String, Int}())
    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        _, rt = sorted_rank(df, aggregate_stats_auto, val_metric, tst_metric, downsample)

        rt[end-2, 1] = "\$\\sigma_1\$"
        rt[end-1, 1] = "\$\\sigma_{10}\$"
        rt[end, 1] = "rnk"

        file = "$(projectdir())/paper/tables/$(prefix)_$(metric)_$(metric)_autoagg$(suffix).tex"
        open(file, "w") do io
            print_rank_table(io, rt; backend=:tex)
        end
    end
end

# this is just specific figure for the main text
# grouptype = true  - renames models according to model type
# grouptype = false - assumes that the results are filtered to a few models/representatives
function plot_knowledge_combined(df_tab, df_img_stat, df_img_sem; grouptype=true, format="pdf", suffix="")
    df_tab, df_tab_clean = df_tab;
    df_img_stat, df_img_stat_clean = df_img_stat;
    df_img_sem, df_img_sem_clean = df_img_sem;

    # tabulars are single class anomaly datasets
    df_tab[:anomaly_class] = -1
    # not include anomaly_class in DataFrame to avoid a wall of warnings about missing anomaly_class from maxmean agg
    # df_img_stat_clean = select(df_img_stat_clean, Not(:anomaly_class))
    # df_img_sem_clean = select(df_img_sem_clean, Not(:anomaly_class))

    # rename models to their corresponding type
    if grouptype
        for df in [df_tab, df_tab_clean, df_img_stat, df_img_stat_clean, df_img_sem, df_img_sem_clean]
            apply_aliases!(df, col="modelname", d=MODEL_TYPE);
        end
    end

    for (mn, metric) in [collect(zip(["AUC", "TPR@5"],[:auc, :tpr_5]))[1]]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        function _incremental_rank(df, criterions, agg)
            ranks, metric_means = [], []
            for criterion in criterions
                df_agg = agg(df, criterion)
                apply_aliases!(df_agg, col="modelname", d=MODEL_RENAME)
                apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
                sort!(df_agg, [:dataset, :modelname])
                rt = rank_table(df_agg, tst_metric)
                # _meanfinite(a) = mean([x for x in a if !isnan(x)])
                # mm = DataFrame([Symbol(model) => _meanfinite(rt[1:end-3, model]) for model in names(rt)[2:end]])
                mm = DataFrame([Symbol(model) => mean(rt[1:end-3, model]) for model in names(rt)[2:end]])
                # @info("", criterion, rt, mm)
                push!(ranks, rt[end:end, 2:end])
                push!(metric_means, mm)
            end
            vcat(ranks...), vcat(metric_means...)
        end

        function _plot(df_ranks, metric_mean, cnames, models, title)
            a = PGFPlots.Axis([PGFPlots.Plots.Linear(
                        1:length(cnames),
                        df_ranks[:, m]) for m in models],
                ylabel="avg. rnk",
                title=title,
                style="xtick=$(_pgf_array(1:length(cnames))),
                    xticklabels=$(_pgf_array(cnames)),
                    width=5cm, height=3cm, scale only axis=true,
                    x tick label style={rotate=50,anchor=east},
                    title style={at={(current bounding box.south)}, anchor=west}")
            b = PGFPlots.Axis([PGFPlots.Plots.Linear(
                        1:length(cnames),
                        metric_mean[:, m],
                            legendentry=m) for m in models],
                ylabel="avg. $mn",
                legendStyle = "at={(0.3,1.30)}, anchor=west",
                style="width=5cm, height=3cm, scale only axis=true,
                xtick=$(_pgf_array(1:length(cnames))),
                xticklabels={}")
            a, b
        end

        for (ctype, cnames, criterions) in collect(zip(
            ["pat", "pac", "patn"],
            [PAT_METRICS_NAMES, PAC_METRICS_NAMES, PATN_METRICS_NAMES],
            [_prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("val", PATN_METRICS)]))[1:1]

            extended_criterions = vcat(criterions, [val_metric])
            extended_cnames = vcat(["clean"], vcat(cnames, ["\$$(mn)_{val}\$"]))
            titles = ["(tab)", "(stat)", "(semantic)"]

            ab_plots = map(enumerate([(df_tab, df_tab_clean), (df_img_stat, df_img_stat_clean), (df_img_sem, df_img_sem_clean)])) do (i, (df, df_clean))
                ranks_clean, metric_means_clean = _incremental_rank(df_clean, [val_metric], aggregate_stats_max_mean)
                ranks_inc, metric_means_inc = _incremental_rank(df, extended_criterions, aggregate_stats_auto)

                ranks_all, metric_means_all = vcat(ranks_clean, ranks_inc; cols=:intersect), vcat(metric_means_clean, metric_means_inc; cols=:intersect)
                # @info("", ranks_clean, ranks_inc)
                # @info("", metric_means_clean, metric_means_inc)

                # reorder table on tabular data as there is additional class of models (flows)
                # one can do this manually at the end
                if i == 1 && grouptype
                    p = [1,2,4,5,3] # make sure that flows are last
                    models = names(ranks_all)[p]
                    select!(ranks_all, p)
                    select!(metric_means_all, p)
                    a, b = _plot(ranks_all, metric_means_all, extended_cnames, models, titles[i])
                else
                    models = names(ranks_all)
                    a, b = _plot(ranks_all, metric_means_all, extended_cnames, models, titles[i])
                end
                a, b
            end

            a_tab, b_tab = ab_plots[1]
            a_img_stat, b_img_stat = ab_plots[2]
            a_img_sem, b_img_sem = ab_plots[3]

            a_img_stat.ylabel = ""
            b_img_stat.ylabel = ""
            a_img_sem.ylabel = ""
            b_img_sem.ylabel = ""
            if grouptype
                for (bst, bsm) in zip(b_img_stat.plots, b_img_sem.plots)
                    bst.legendentry=nothing
                    bsm.legendentry=nothing
                end
            end

            g = PGFPlots.GroupPlot(3, 2, groupStyle = "vertical sep = 0.5cm, horizontal sep = 1.0cm")
            push!(g, b_tab); push!(g, b_img_stat); push!(g, b_img_sem);
            push!(g, a_tab); push!(g, a_img_stat); push!(g, a_img_sem);

            file = "$(projectdir())/paper/figures/combined_knowledge_rank_$(ctype)_$(metric)$(suffix).$(format)"
            PGFPlots.save(file, g; include_preamble=false)
        end
    end
end

plot_knowledge_combined(
    (df_tabular, df_tabular_clean),
    (df_images_stat, df_images_stat_clean),
    (df_images_semantic, df_images_semantic_clean); format="tex", suffix="_grouptype")

ff(df) = filter(x -> x.modelname in representatives, df)

# one high row with all models
representatives=["ocsvm", "vae_full", "vae", "wae_full" , "wae", "fmgan", "vae_ocsvm", "knn", "aae_full", "aae", "RealNVP", "fAnoGAN", "DeepSVDD"]
df_tab = df_tabular |> ff |> copy, df_tabular_clean |> ff |> copy;
df_img_stat = df_images_stat |> ff |> copy, df_images_stat_clean |> ff |> copy;
df_img_sem = df_images_semantic |> ff |> copy, df_images_semantic_clean |> ff |> copy;

plot_knowledge_combined(df_tab, df_img_stat, df_img_sem;
                    grouptype=false, format="tex", suffix="_repre_new_per_ac_one_all")
@info "plot_combine_knowledge_type"

@info "----------------- DONE ---------------------"


text:
function is_convex_mesh_oriented_outward(eM::eMesh{Tri,T2}) where {T2}
	for k = 1:n_tri(eM)
		p3 = vertex_pos_for_tri_ind(eM, k)
		n̂ = triangleNormal(p3)
		c = centroid(p3)
		c = normalize(c)
	end
	return true
end

@testset "swept_mesh" begin
	eM_box = create_swept_mesh(f_swept_triv, (-1.0, 1.0), 1.0, 4, true)
	@test is_convex_mesh_oriented_outward(eM_box)
	@test 8.0 ≈ volume(eM_box)

	eM_box_check = output_eMesh_box()
	@test is_convex_mesh_oriented_outward(eM_box_check)
	@test 8.0 ≈ volume(eM_box_check)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
using LinearAlgebra
# using ForwardDiff
using NLsolve

# http://www.solid.iei.liu.se/Education/TMHL55/TMHL55_lp1_2010/lecture_notes/plasticity_flow_rule_isotropic_hardening.pdf
# http://mms.ensmp.fr/msi_paris/transparents/Georges_Cailletaud/2013-GC-plas3D.pdf

#####################################
# Viscoplastic material definitions #
#####################################
mutable struct ViscoPlastic <: AbstractMaterial
    # Material parameters
    youngs_modulus :: Float64
    poissons_ratio :: Float64
    yield_stress :: Float64
    potential :: Symbol
    params :: Array{Float64, 1}

    # Internal state variables
    plastic_strain :: Vector{Float64}
    dplastic_strain :: Vector{Float64}
    plastic_multiplier :: Float64
    dplastic_multiplier :: Float64

end

POTENTIAL_FUNCTIONS = [:norton]

function ViscoPlastic(potential::Symbol, params :: Array{Float64,1})
    if !any(POTENTIAL_FUNCTIONS .== potential)
        error("Potential $potential not found!")
    end
    youngs_modulus = 0.0
    poissons_ratio = 0.0
    yield_stress = 0.0
    plastic_strain = zeros(6)
    dplastic_strain = zeros(6)
    plastic_multiplier = 0.0
    dplastic_multiplier = 0.0
    dt = 0.0
    return ViscoPlastic(youngs_modulus, poissons_ratio, yield_stress, potential, params,
                        plastic_strain, dplastic_strain, plastic_multiplier,
                        dplastic_multiplier)
end

function ViscoPlastic()
    youngs_modulus = 0.0
    poissons_ratio = 0.0
    yield_stress = 0.0
    plastic_strain = zeros(6)
    dplastic_strain = zeros(6)
    plastic_multiplier = 0.0
    dplastic_multiplier = 0.0
    dt = 0.0
    potential = :norton
    params = [180.0e3, 0.92]
    return ViscoPlastic(youngs_modulus, poissons_ratio, yield_stress, potential, params,
                        plastic_strain, dplastic_strain, plastic_multiplier,
                        dplastic_multiplier)
end

""" Double contraction
"""
function double_contraction(x,y)
    return sum(x.*y)
end

""" Deviatoric stress tensor
"""
function deviatoric_stress(stress)
    stress_tensor = [   stress[1] stress[4] stress[6];
                        stress[4] stress[2] stress[5];
                        stress[6] stress[5] stress[3]]
    stress_dev = stress_tensor - 1/3 * tr(stress_tensor) * Matrix(1.0I, 3, 3)
    return stress_dev
end

function J_2_stress(stress)
    s = deviatoric_stress(stress)
    return 1/2 * double_contraction(s,s)
end

""" Equivalent stress
"""
function equivalent_stress(stress)
    J_2 = J_2_stress(stress)
    return sqrt(3 * J_2)
end

""" Von mises yield function
"""
function von_mises_yield(stress, stress_y)
    return equivalent_stress(stress) - stress_y
end

""" Norton rule
"""
function norton_plastic_potential(stress, params, f)
    K = params[1]
    n = params[2]
    J = f(stress)
    return K/(n+1) * (J / K) ^ (n + 1)
end

function bingham_plastic_potential(stress, params, f)
    eta = params[1]
    return 0.5 * (f(stress) / eta) ^2
end

""" Analytically derivated norton rule
"""
function dNortondStress(stress, params, f)
    K = params[1]
    n = params[2]
    f_ = f(stress)
    stress_v = equivalent_stress(stress)
    stress_dev_vec = deviatoric_stress(stress)
    sol_mat = (f_ / K) ^ n * 3 / 2 * stress_dev_vec / stress_v
    return [sol_mat[1,1], sol_mat[2,2], sol_mat[3,3],
            2*sol_mat[1,2], 2*sol_mat[2,3], 2*sol_mat[1,3]]
end

function find_root(f, df, x; max_iter=50, norm_acc=1e-9)
    converged = false
    for i=1:max_iter
        dx = -df(x) \ f(x)
        x += dx
        norm(dx) < norm_acc && (converged = true; break)
    end
    converged || error("No convergence in radial return!")
    return x
end

function integrate_material!(material::Material{ViscoPlastic})
    mat = material.properties

    E = mat.youngs_modulus
    nu = mat.poissons_ratio
    potential = mat.potential
    params = mat.params
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

    stress = material.stress
    strain = material.strain
    dstress = material.dstress
    dstrain = material.dstrain
    dt = material.dtime
    D = material.jacobian

    dedt = dstrain ./ dt

    fill!(D, 0.0)
    D[1,1] = D[2,2] = D[3,3] = 2.0*mu + lambda
    D[4,4] = D[5,5] = D[6,6] = mu
    D[1,2] = D[2,1] = D[2,3] = D[3,2] = D[1,3] = D[3,1] = lambda

    dstress[:] .= D * dedt .* dt
    stress_tr = stress + dstress
    f = x -> von_mises_yield(x, mat.yield_stress)
    yield = f(stress_tr)

    if yield <= 0
        fill!(mat.dplastic_strain, 0.0)
        return nothing
    else
        if potential == :norton
            x0 = vec([dstress; 0.0])

            function g!(F, x)
                dsigma = x[1:6]
                F[1:6] = dsigma - D * (dedt - dNortondStress(stress + dsigma, params, f)) .* dt
                F[end] = von_mises_yield(stress + dsigma, mat.yield_stress)
            end

            res = nlsolve(g!, x0)
            dsigma = res.zero[1:6]
            dstrain_pl = dNortondStress(stress + dsigma, params, f)

        end
        mat.dplastic_strain = dstrain_pl
        dstress[:] .= D*(dedt - dstrain_pl) .* dt
        return nothing
    end

    return nothing

end


text:
"""
    getclass(repo::Repo{T},id::Int)

Return the class number of element id in  repo.

# Example
```jldoctest
julia> using Repos

julia> r=Repo(["red","green","blue"],2)
Repository of 3 objects in 2 classes

julia> getclass(r,1)
1
```
"""
function getclass(repo::Repo{T},id::Int) where T
    pos=repo.pos[id]
    for c in 1:repo.C
        start=repo.starts[c]
        cend=start+repo.counts[c]
        if pos<cend
            return c
        end
    end
    error("Repo could not determine class of $id")
    return 0
end

"""
    setclass!(repo::Repo{T},id::Int,cls::Int) where T

Change the class of item number id in repo to cls.

# Example
```
julia> using Repos

julia> r=Repo(["red","green","blue"],2)
Repository of 3 objects in 2 classes

julia> setclass!(r,1,2)

julia> print_repo(r)
Repository of 2 classes
  Class 1
    1 - 3: blue
    2 - 2: green
  Class 2
    1 - 1: red
```
"""
function setclass!(repo::Repo{T},id::Int,cls::Int) where T
    curc::Int=getclass(repo,id)
    while curc<cls
        lastposincur=repo.starts[curc]+repo.counts[curc]-1
        curpos=repo.pos[id]
        _swap(repo,lastposincur,curpos)
        repo.counts[curc]-=1
        curc+=1
        repo.starts[curc]-=1
        repo.counts[curc]+=1
    end
    while curc>cls
        firstposincur=repo.starts[curc]
        curpos=repo.pos[id]
        _swap(repo,firstposincur,curpos)
        repo.counts[curc]-=1
        repo.starts[curc]+=1
        curc-=1
        repo.counts[curc]+=1
    end
    return nothing
end

"""
     class(repo::Repo{T},c::Int)

Return an interface to class c in repository repo.

The returned RepoClass structure is an interface to the object of the respective
class in the repository. The RepoClass will stay up to date when the repository
is updated.

# Example
```jldoctest
julia> using Repos

julia> a=["Golden Delicious","Granny Smith","Red Delicious","Valencia"]
4-element Array{String,1}:
 "Golden Delicious"
 "Granny Smith"
 "Red Delicious"
 "Valencia"

julia> r=Repo(a,2)
Repository of 4 objects in 2 classes

julia> apples=class(r,1)
Class of 4 objects

julia> setclass!(r,4,2)

julia> a[3]="Braeburn"
"Braeburn"

julia> print_repo(apples)
Class of 3 objects
    1 - 1: Golden Delicious
    2 - 2: Granny Smith
    3 - 3: Braeburn
```
"""
function class(repo::Repo{T},c::Int) where T
    return RepoClass{T}(repo,c)
end


text:
__precompile__()

module MUMPSjInv

using Distributed
using SparseArrays
using Printf

mutable struct MUMPSfactorization{T}
    ptr::Int64     # pointer to factorization
    worker::Int64  # id of worker that holds factorization
    n::Int64       # matrix size
    a11::T         # first element (HACK for generating the right parametric type)
    time::Float64  # factorization time
end

const MUMPSlibPath  = abspath(joinpath(splitdir(Base.source_path())[1], "..", "lib", "MUMPS"))

include("MUMPSfuncs.jl")

export solveMUMPS, solveMUMPS!, factorMUMPS, applyMUMPS, applyMUMPS!, destroyMUMPS, MUMPSfactorization

end


text:
using ..TopOpt: ElementFEAInfo, ElementMatrix
using ..TopOpt.TopOptProblems: make_cload, convert

"""
    ElementFEAInfo(sp, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {mat_type}

Constructs an instance of `ElementFEAInfo` from a stiffness **truss** problem `sp` using a Gaussian quadrature order of `quad_order`. The element matrix and vector types will be:
1. `SMatrix` and `SVector` if `mat_type` is `:SMatrix` or `:Static`, the default,
2. `MMatrix` and `MVector` if `mat_type` is `:MMatrix`, or
3. `Matrix` and `Vector` otherwise.

The static matrices and vectors are more performant and GPU-compatible therefore they are used by default.
"""
function ElementFEAInfo(
    sp::TrussProblem, quad_order=1, ::Type{Val{mat_type}}=Val{:Static}
) where {mat_type}
    # weights: self-weight element load vectors, all zeros now
    Kes, weights, cellvalues, facevalues = make_Kes_and_fes(sp, quad_order, Val{mat_type})
    element_Kes = convert(
        Vector{<:ElementMatrix},
        Kes;
        bc_dofs=sp.ch.prescribed_dofs,
        dof_cells=sp.metadata.dof_cells,
    )

    # * concentrated load
    # ? why convert a sparse vector back to a Vector?
    fixedload = Vector(make_cload(sp))
    # * distributed load (if any)
    # assemble_f!(fixedload, sp, dloads)

    cellvolumes = get_cell_volumes(sp, cellvalues)
    cells = sp.ch.dh.grid.cells
    return ElementFEAInfo(
        element_Kes,
        weights,
        fixedload,
        cellvolumes,
        cellvalues,
        facevalues,
        sp.metadata,
        sp.black,
        sp.white,
        sp.varind,
        cells,
    )
end

####################################

function get_cell_volumes(sp::TrussProblem{xdim,T}, cellvalues) where {xdim,T}
    dh = sp.ch.dh
    As = getA(sp)
    cellvolumes = zeros(T, getncells(dh.grid))
    for (i, cell) in enumerate(CellIterator(dh))
        truss_reinit!(cellvalues, cell, As[i])
        cellvolumes[i] = sum(
            Ferrite.getdetJdV(cellvalues, q_point) for
            q_point in 1:Ferrite.getnquadpoints(cellvalues)
        )
    end
    return cellvolumes
end

####################################

"""
    compute_local_axes(end_vert_u, end_vert_v)

# Arguments
`end_vert_u, end_vert_v` = vectors for nodal coordinate

# Outputs
`R` = (ndim x ndim) global_from_local transformation matrix
    Each column of the matrix is the local system's axis coordinate in the world system.
    Note that this matrix has its columns as axes
    So should be used as `R*K*R'` instead of `R'*K*R` as indicated in
    https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf
"""
function compute_local_axes(end_vert_u, end_vert_v)
    @assert length(end_vert_u) == length(end_vert_v)
    @assert length(end_vert_u) == 2 || length(end_vert_u) == 3
    xdim = length(end_vert_u)
    L = norm(end_vert_u - end_vert_v)
    @assert L > eps()
    # by convention, the new x axis is along the element's direction
    # directional cosine of the new x axis in the global world frame
    c_x = (end_vert_v[1] - end_vert_u[1]) / L
    c_y = (end_vert_v[2] - end_vert_u[2]) / L
    R = zeros(xdim, xdim)
    if 3 == xdim
        c_z = (end_vert_v[3] - end_vert_u[3]) / L
        if abs(abs(c_z) - 1.0) < eps()
            R[1, 3] = -c_z
            R[2, 2] = 1
            R[3, 1] = c_z
        else
            # local x_axis = element's vector
            new_x = [c_x, c_y, c_z]
            # local y axis = cross product with global z axis
            new_y = -cross(new_x, [0, 0, 1.0])
            new_y /= norm(new_y)
            new_z = cross(new_x, new_y)
            R[:, 1] = new_x
            R[:, 2] = new_y
            R[:, 3] = new_z
        end
    elseif 2 == xdim
        R = [
            c_x -c_y
            c_y c_x
        ]
    end
    return R
end


text:
module MockDofsTests

using Test
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField
using InplaceArrays.ReferenceFEs
using InplaceArrays.ReferenceFEs: MockDofBasis

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
f = MockField{d}(v)
fx = evaluate(f,x)

b = MockDofBasis(x)
test_dof(b,f,fx)

l = 10
ab = fill(b,l)
af = fill(f,l)
ax = fill(x,l)

afx = evaluate(af,ax)
abf = evaluate(ab,af)
@test abf == afx

end # module


text:
# This is testing GDAL/OGR utility executables.
# In `gdal_utils.jl` we test the utility library functions.

@test success(GDAL.gdalinfo_path() do gdinfo
    run(pipeline(`$gdinfo $utmsmall`, devnull))
end)


text:
module TestFunctionLenses
using Test
using Accessors: test_getset_laws
using Accessors

@testset "os" begin
    path = "hello.md"
    path_new = @set splitext(path)[2] = ".jl"
    @test path_new == "hello.jl"

    path = joinpath("root", "somedir", "some.file")
    path_new = @set splitdir(path)[1] = "otherdir"
    @test path_new == joinpath("otherdir", "some.file")

    test_getset_laws(splitext, "hello.world", ("hi", ".jl"), ("ho", ".md"))
    test_getset_laws(splitdir, joinpath("hello", "world"), ("a", "b"), ("A", "B"))
    test_getset_laws(splitpath, joinpath("hello", "world"), ["some"], ["some", "long", "path"])
    test_getset_laws(dirname, joinpath("hello", "world"), "hi", "ho")
    test_getset_laws(basename, joinpath("hello", "world"), "planet", "earth")
end

@testset "first" begin
    obj = (1, 2.0, '3')
    l = @optic first(_)
    @test l === first
    @test l(obj) === 1
    @test set(obj, l, "1") === ("1", 2.0, '3')
    @test (@set first(obj) = "1") === ("1", 2.0, '3')

    obj2 = (a=((b=1,), 2), c=3)
    @test (@set first(obj2.a).b = '1') === (a=((b='1',), 2), c=3)
end

@testset "last" begin
    obj = (1, 2.0, '3')
    l = @optic last(_)
    @test l === last
    @test set(obj, l, '4') === (1, 2.0, '4')
    @test (@set last(obj) = '4') === (1, 2.0, '4')

    obj2 = (a=(1, (b=2,)), c=3)
    @test (@set last(obj2.a).b = '2') === (a=(1, (b='2',)), c=3)
end

@testset "eltype on Number" begin
    @test @set(eltype(Int) = Float32) === Float32
    @test @set(eltype(1.0) = UInt8)   === UInt8(1)

    @inferred set(Int, eltype, Float32)
    @inferred set(1.2, eltype, Float32)

end

@testset "eltype(::Type{<:Array})" begin
    obj = Vector{Int}
    @inferred set(obj, eltype, Float32)
    obj2 = @set eltype(obj) = Float64
    @test obj2 === Vector{Float64}
end

@testset "eltype(::Array)" begin
    obj = [1, 2, 3]
    @inferred set(obj, eltype, Float32)
    obj2 = @set eltype(obj) = Float64
    @test eltype(obj2) == Float64
    @test obj == obj2
end

@testset "(key|val|el)type(::Type{<:Dict})" begin
    obj = Dict{Symbol, Int}
    @test (@set keytype(obj) = String) === Dict{String, Int}
    @test (@set valtype(obj) = String) === Dict{Symbol, String}
    @test (@set eltype(obj) = Pair{String, Any}) === Dict{String, Any}

    obj2 = Dict{Symbol, Dict{Int, Float64}}
    @test (@set keytype(valtype(obj2)) = String) === Dict{Symbol, Dict{String, Float64}}
    @test (@set valtype(valtype(obj2)) = String) === Dict{Symbol, Dict{Int, String}}
end

@testset "(key|val|el)type(::Dict)" begin
    obj = Dict(1 => 2)
    @test typeof(@set keytype(obj) = Float64) === Dict{Float64, Int}
    @test typeof(@set valtype(obj) = Float64) === Dict{Int, Float64}
    @test typeof(@set eltype(obj) = Pair{UInt, Float64}) === Dict{UInt, Float64}
end

@testset "math" begin
    x = 1
    @test 2.0       === @set real(1) = 2.0
    @test 1.0 + 2im === @set imag(1) = 2.0
    @test 1.0 + 2im === @set imag(1+1im) = 2.0

end

end # module


text:
using Oceananigans, Oceananigans.Units, GLMakie

# Specify domain and mesh / grid
grid = RectilinearGrid(size = (60, 60, 1),
                       x = (0, 1200kilometers),
                       y = (0, 1200kilometers),
                       z = (-5000, 0),
                       halo = (3, 3, 3),
                       topology = (Bounded, Bounded, Bounded))

# Specify coriolis, viscosity / turbulence closure, and boundary condition
coriolis = BetaPlane(f₀=1e-4, β=1e-11)
closure = IsotropicDiffusivity(ν=400)

wind_stress(x, y, t) = - 1e-4 * cos(π * x / 600kilometers)
u_bcs = FieldBoundaryConditions(; top=FluxBoundaryCondition(wind_stress))

model = HydrostaticFreeSurfaceModel(; grid, coriolis, closure,
                                      momentum_advection = UpwindBiasedFifthOrder(),
                                      boundary_conditions = (; u=u_bcs))

# Time-stepping
simulation = Simulation(model, Δt=20minutes, stop_time=3years)

progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

run!(simulation)

fig, ax, pl = heatmap(interior(model.velocities.u)[:, :, 1])
display(fig)


text:
module TestGetters
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO.Getters
using Setfield

@testset "test getters" begin
    p = arbitrary(IAEAParticle{0,1})
    p = @set p.x = 100
    @test x(p) ≈ 100
    p = @set p.E = 10
    p = @set p.weight = 1
    @test energy(p) == 10
    p = @set p.weight = 0.1
    @test energy(p) ≈ 1

    p = arbitrary(EGSParticle{Float32})
    p = @set p.x = 1000f0
    @test x(p) ≈ 1000
    p = @set p.zlast = nothing
    @test p.zlast == nothing
    p = @set p.zlast = 1f0
    @test p.zlast == 1f0
end

end


text:
using CompScienceMeshes
using BEAST
using Test

############################################################################################
#
# Test the sign of the double layer potential
#
############################################################################################

γ = 0.0
K = MWDoubleLayer3D(γ)

τ = Mesh(
    [
        point(-4,0,0),
        point(-6,0,0),
        point(-5,-1,0)
    ],
    [
        index(1,2,3)])

σ = Mesh(
    [
        point(4,0,0),
        point(6,0,0),
        point(5,0,-1)
    ],
    [
        index(1,2,3)
    ]
)

@test numcells(τ) == 1
@test numcells(σ) == 1

T = raviartthomas(τ, boundary(τ))
S = raviartthomas(σ, boundary(σ))

@test numfunctions(T) == 3
@test numfunctions(S) == 3

t = chart(τ, first(cells(τ)))
s = chart(σ, first(cells(σ)))

p = neighborhood(t, [1,1]/2)
q = neighborhood(s, [1,1]/2)

R = norm(cartesian(p)-cartesian(q))

rt = refspace(T)
f = (refspace(T))(p)[3][1]
g = (refspace(S))(q)[3][1]

@test dot(f, point(0,1,0)) > 0
@test dot(g, point(0,0,1)) > 0

@show 1/R^2

z = assemble(K, T, S)

qd = quaddata(K, rt, rt, elements(τ), elements(σ))
zlocal = zeros(promote_type(scalartype(K), scalartype(S), scalartype(T)), 3, 3)
strat = BEAST.quadrule(K, rt, rt, 1, t, 1, s, qd)
strat =  BEAST.DoubleQuadStrategy(qd.tpoints[1,1],  qd.bpoints[1,1])
BEAST.momintegrals!(K, rt, rt, t, s, zlocal, strat)
z = zlocal[3,3]