text:
# To test a simplified co ax copter.
# Modeling the system as a cube, gimbal and 2 propellers.
# A revolute2 joint between the cube and gimbal
# 1 Rev Joint each between the gimbal and either propeller.

include("../src/plotSol.jl")
include("../src/OrientationConversion.jl")
include("../src/simulate.jl")
#  include("../src/simulateRK.jl")

using JLD
clearconsole()

## Inertial Properties
const mb = 0.755; const mg = 0.215; const mp = 0.02;

const Ib =  [ 0.0009    0.0000   -0.0000;
        0.0000    0.0053    0.0000;
       -0.0000    0.0000    0.0054]#*(10^-3)

const Ip =    [ 0.0029    0.0000    0.0000;
          0.0000    0.5824    0.0000;
          0.0000    0.0000    0.5805]*(10^-4)

const Ig =    [ 0.1435    0.0000    0.0000;
          0.0000    0.0498    0.0000;
          0.0000    0.0000    0.1381]*(10^-3)

## Initialisation
InFrame = InertialFrameAsRB()
Body1 = RigidBody(mb,Ib,2)
Gimbal1 = RigidBody(mg,Ig,3)
Prop11 = RigidBody(mp,Ip,4)
Prop12 = RigidBody(mp,Ip,5)

Box = RigidBody(mb,Ib,6)

Body2 = RigidBody(mb,Ib,7)
Gimbal2 = RigidBody(mg,Ig,8)
Prop21 = RigidBody(mp,Ip,9)
Prop22 = RigidBody(mp,Ip,10)

# CoAxCop1
x0Body1 = [zeros(3);[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Body1,x0Body1)

x0Gimbal1 = [[0.0;0.0;0.143];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Gimbal1,x0Gimbal1)

x0Prop11 = [[0.0;0.0;0.234];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Prop11,x0Prop11)

x0Prop12 = [[0.0;0.0;0.24];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Prop12,x0Prop12)

# Box
x0Box = [[0.08;0.0;0.0];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Box,x0Box)

# CoAxCop2
x0Body2 = [[0.16;0.0;0.0];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Body2,x0Body2)

x0Gimbal2 = [[0.16;0.0;0.143];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Gimbal2,x0Gimbal2)

x0Prop21 = [[0.16;0.0;0.234];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Prop21,x0Prop21)

x0Prop22 = [[0.16;0.0;0.24];[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Prop22,x0Prop22)

## Joint Descriptions
# Joint 1 (Free) between the inertial frame and body.
j1 = Joint(InFrame,Body1,zeros(3),zeros(3))

# CoAxCop1
# Joint 2 (Revolute2) between the body and gimbal
rjBody = [0.0 0.0 0.134][:]
rjGimbal1 = [0.0 0.0 -0.009][:]
j2 = Joint(Body1,Gimbal1,rjBody,rjGimbal1,
     type = "Revolute",jointTorque = [0.0 0.0 0.0][:])

# Joint 3 (Revolute) between the gimbal and prop11
rjGimbal2 = [0.0 0.0 0.091][:]
rjProp1 = [0.0 0.0 0.0][:]
j3 = Joint(Gimbal1, Prop11, rjGimbal2, rjProp1, type = "Revolute",jointTorque = [0.0 0.0 0.001][:])

# Joint 4 (Revolute) between the gimbal and prop12
rjGimbal3 = [0.0 0.0 0.0965][:]
rjProp2 = [0.0 0.0 0.0][:]
j4 = Joint(Gimbal1, Prop12, rjGimbal3, rjProp2, type = "Revolute", jointTorque = [0.0 0.0 -0.001][:])

# Between CACs and Box
# Joint 5 between CoAxCop1 and Box
rjCac1 = [0.04 0.0 0.0][:]
rjBox1 = [-0.04 0.0 0.0][:]
j5 = Joint(Body1, Box, rjCac1, rjBox1, type = "Weld")

# Joint 6 between Box and CoAxCop2
rjBox2 = [0.04 0.0 0.0][:]
rjCac2 = [-0.04 0.0 0.0][:]
j6 = Joint(Box, Body2, rjBox2, rjCac2, type = "Weld")

# CoAxCop2
# Joint 7 (Revolute2) between the body and gimbal
# const rjBody = [0.0 0.0 0.134][:]
# const rjGimbal1 = [0.0 0.0 -0.009][:]
j7 = Joint(Body2,Gimbal2,rjBody,rjGimbal1,
     type = "Revolute")#,jointTorque = [0.0 0.0 0.0][:])

# Joint 8 (Revolute) between the gimbal and prop21
# const rjGimbal2 = [0.0 0.0 0.091][:]
# const rjProp1 = [0.0 0.0 0.0][:]
j8 = Joint(Gimbal2, Prop21, rjGimbal2, rjProp1, type = "Revolute",jointTorque = [0.0 0.0 0.001][:])

# Joint 9 (Revolute) between the gimbal and prop22
# const rjGimbal3 = [0.0 0.0 0.0965][:]
# const rjProp2 = [0.0 0.0 0.0][:]
j9 = Joint(Gimbal2, Prop22, rjGimbal3, rjProp2, type = "Revolute", jointTorque = [0.0 0.0 -0.001][:])

## Simulation
tEnd = 1.0
tSpan = 0.01
const g = MVector{3}([0.0,0.0,-9.806]) # Gravity Vector.
j = (j1,j2,j3,j4,j5,j6,j7,j8,j9) # Tuple of joints
@time tSim, solFinal = simulate(tEnd,tSpan,j...,g=g);
# @time tSim, solFinal = simulateRK45(tEnd, tSpan, j..., g=g);

# save("ModRob1dSecSimRK45NoTol.jld","solFinal", solFinal,"tSim", tSim)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text: tojadencase(sentence::String) = titlecase(sentence, wordsep=isspace)

ext: jl

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Test
using Qlab: savitzky_golay_filter

# test the innput validity checks
# polynomial order greater than window size
@test_throws AssertionError savitzky_golay_filter(1:10, 5, 5)
# window size must be odd
@test_throws AssertionError savitzky_golay_filter(1:10, 10, 2)
# boundary_mode must be one :nearest or :interpolation
@test_throws AssertionError savitzky_golay_filter(1:10, 11, 2; boundary_mode=:silly)

# a window size of 1 with a zeroth order polynomial should return the data
@test savitzky_golay_filter([1.2], 1, 0) == [1.2]
@test savitzky_golay_filter([1.2, 3.4, 5.6], 1, 0) == [1.2, 3.4, 5.6]

# a line smoothed with a polynomial order of 1 should return itself or `boundary_mode==:interpolation`
@test savitzky_golay_filter(1:10, 3, 1; boundary_mode=:interpolation) ≈ 1:10

# take a low order polynomial and smooth with a higher order polynomial to return the exact value
xpts = range(1,10; length=81)
data = 1.2*xpts.^3 + 3.4*xpts.^2 + 5.6*xpts .+ 7.8

@test savitzky_golay_filter(data, 5, 4) ≈ data
@test savitzky_golay_filter(data, 11, 4) ≈ data

# derivatives are easy to analytically calculate
# we need to scale savitzky_golay_filter by the spacing
deriv1 = 3*1.2*xpts.^2 + 2*3.4*xpts .+ 5.6
@test savitzky_golay_filter(data, 5, 4; deriv_order=1) / (xpts[2] - xpts[1]) ≈ deriv1
deriv2 = 2*3*1.2*xpts .+ 2*3.4
@test savitzky_golay_filter(data, 5, 4; deriv_order=2) / (xpts[2] - xpts[1])^2 ≈ deriv2

# test the :nearest boundary condition
# take 2 points and use a 3 point constant smoothing fit
# first point will be average of [1, 1, 2] = 4/3 and second point will be of [1, 2, 2] = 5/3
@test savitzky_golay_filter([1, 2], 3, 0; boundary_mode=:nearest) ≈ [4/3, 5/3]


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiStorageV1beta1CSINodeSpec <: SwaggerModel
    drivers::Any # spec type: Union{ Nothing, Vector{IoK8sApiStorageV1beta1CSINodeDriver} } # spec name: drivers

    function IoK8sApiStorageV1beta1CSINodeSpec(;drivers=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1CSINodeSpec, Symbol("drivers"), drivers)
        setfield!(o, Symbol("drivers"), drivers)
        o
    end
end # type IoK8sApiStorageV1beta1CSINodeSpec

const _property_map_IoK8sApiStorageV1beta1CSINodeSpec = Dict{Symbol,Symbol}(Symbol("drivers")=>Symbol("drivers"))
const _property_types_IoK8sApiStorageV1beta1CSINodeSpec = Dict{Symbol,String}(Symbol("drivers")=>"Vector{IoK8sApiStorageV1beta1CSINodeDriver}")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1CSINodeSpec }) = collect(keys(_property_map_IoK8sApiStorageV1beta1CSINodeSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1CSINodeSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1beta1CSINodeSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1CSINodeSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1CSINodeSpec[property_name]

function check_required(o::IoK8sApiStorageV1beta1CSINodeSpec)
    (getproperty(o, Symbol("drivers")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1beta1CSINodeSpec }, name::Symbol, val)
end


text:
# Simulation example for variable selection using power likelihood
include("varsel.jl")

module Run
using Distributions, HDF5, JLD
using VarSel
using PyPlot
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())

# Settings
ns = [100,1000,5000,10000,50000]  # sample sizes to use
alphas = [Inf,50]  # values of alpha to use
nreps = 10  # number of times to repeat each simulation
n_total = 5*10^4  # total number of MCMC samples
n_keep = 10^4  # number of MCMC samples to record
nburn = int(n_keep/10)  # burn-in
from_file = false  # load previous results from file

# Data distribution
using Skew
a = [0.6,2.7,-3.3,-4.9,-2.5]
Q = [1.0 -0.89 0.93 -0.91 0.98
     -0.89 1.0 -0.94 0.97 -0.91
     0.93 -0.94 1.0 -0.96 0.97
     -0.91 0.97 -0.96 1.0 -0.93
     0.98 -0.91 0.97 -0.93 1.0]
sigma = 1
p = length(a)+1

g(x0) = -1 + 4*(x0 + (1/16)*x0.^2)

function generate_sample(n)
    x = [ones(n) Skew.skewrndNormalized(n,Q,a)']
    x0 = vec(x[:,2])
    y = g(x0) + randn(n)*sigma
    return x,y
end

# Plot data
srand(1)
figure(5,figsize=(8,3.2)); clf(); hold(true)
subplots_adjust(bottom=0.2)
n = 200
xmin,xmax = -4,4
ymin,ymax = -15,15
xs,ys = generate_sample(n)
xs0 = vec(xs[:,2])
plot(xs0,ys,"b.",markersize=3)
xt = linspace(xmin,xmax,5000)
plot(xt,g(xt),"k-")
xlim(xmin,xmax); ylim(ymin,ymax)
title("\$\\mathrm{Data\\,distribution}\$",fontsize=17)
xlabel("\$x_{i 2}\$",fontsize=18)
ylabel("\$y_i\$",fontsize=18)
draw()
savefig("varsel-data-n=$n.png",dpi=150)

# Run simulations
nns = length(ns)
for (i_a,alpha) in enumerate(alphas)
    if from_file
        k_posteriors = load("k_posteriors-alpha=$alpha.jld","k_posteriors")
    else
        k_posteriors = zeros(p+1,nreps,nns)
        for (i_n,n) in enumerate(ns)
            for rep in 1:nreps
                srand(n+rep) # Reset RNG

                # Sample data
                x,y = generate_sample(n)
                x0 = vec(x[:,2])

                # Run sampler
                zeta = (1/n)/(1/n + 1/alpha)
                beta_r,lambda_r,keepers = VarSel.run(x,y,n_total,n_keep,zeta)
                k_r = vec(sum(beta_r.!=0, 1))

                # Compute posterior on k
                edges,counts = hist(k_r[nburn+1:end],-0.5:1:p+0.5)
                k_posteriors[:,rep,i_n] = counts/(n_keep-nburn)

                # traceplot of number of nonzero coefficients
                figure(10); clf(); hold(true)
                plot(k_r + 0.5*rand(length(k_r))-0.25,"k.",markersize=2)
                draw()
                #savefig("traceplot-alpha=$alpha-n=$n-rep=$rep.png",dpi=150)
            end
        end
        save("k_posteriors-alpha=$alpha.jld","k_posteriors",k_posteriors)
    end

    # Plot results
    figure(i_a,figsize=(8,3.2)); clf(); hold(true)
    subplots_adjust(bottom=0.2,right=0.75)
    colors = "bgyrm"
    shapes = "ds^v*"
    plot(0:p,VarSel.prior_k(0:p,p),"ko-",label="prior",markersize=8,linewidth=2)
    for (i_n,n) in enumerate(ns)
        plot(0:p,vec(mean(k_posteriors[:,:,i_n],2)),"$(colors[i_n])$(shapes[i_n])-",label="n = $n",markersize=8,linewidth=2)
    end
    xlim(0,p)
    xticks(0:p)
    if alpha<Inf; title("\$\\mathrm{Coarsened}\$",fontsize=17)
    else; title("\$\\mathrm{Standard}\$",fontsize=17)
    end
    xlabel("\$k \\,\\,\\mathrm{(\\#\\,of\\,nonzero\\,coefficients)}\$",fontsize=17)
    ylabel("\$\\Pi(k|\\mathrm{data})\$",fontsize=15)
    legend(numpoints=1,bbox_to_anchor=(1.02, 0.9), loc=2, borderaxespad=0.)
    draw()
    savefig("varsel-alpha=$alpha.png",dpi=150)
end

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
cd(dirname(@__FILE__))
include("MolDyn.jl")
using Statistics, DataFrames, CSV
datapath = "data/"

num = 100
sidelen = 30.0
datacount = 41

inittemps = collect(range(1.0, 3.0, length=datacount))
meanT, meanP, σT, σP = [Vector{Float64}(undef, datacount) for _ in 1:4]

for i in 1:datacount
    sim = MolDyn.MolSim(num, inittemps[i], sidelen, sidelen)
    MolDyn.alignleft!(sim)
    MolDyn.simulate!(sim, 4000, 0.01)
    # evolve the system to equilibrium
    # collect data
    _, temperature, pressure = MolDyn.simulate!(sim, 40000, 0.001)
    meanT[i] = mean(temperature)
    meanP[i] = mean(pressure)
    σT[i] = std(temperature)
    σP[i] = std(pressure)
end

# save data to file
df = DataFrame(
    "temperature mean" => meanT, "temperature std" => σT,
    "pressure mean" => meanP, "pressure std" => σP)
CSV.write(datapath * "vanderwaals-N$num-L$sidelen.csv", df)

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
using LinearAlgebra: dot
using OnlineStatsBase
using Statistics

function rca(img::AbstractArray, start::CartesianIndex, thresh::Function, nchords = 16)
    function measurechord(ci, slope)
        dx, sx, c = abs(slope[2]), slope[2] > 0 ? 1 : -1, ci.I[2]
        dy, sy, r = -abs(slope[1]), slope[1] > 0 ? 1 : -1, ci.I[1]
        err = dx + dy
        last = (r, c)
        while c >= 1 && c <= size(img, 2) && r >= 1 && r <= size(img, 1) && thresh(img[r, c])
            last = (r, c)
            e2 = 2err
            if e2 >= dy
                err += dy
                c += sx
            end
            if e2 <= dx
                err += dx
                r += sy
            end
        end
        return CartesianIndex(last...)
    end
    @assert thresh(img[start]) "The initial point does not meet the threshold in drawchord(...)"
    bisectors = (((1, 0), (-1, 0)), ((0, 1), (0, -1)), ((1, 1), (-1, -1)), ((-1, 1), (1, -1)))
    rat = map(i -> rationalize(Int32, tan(2π * i / (2 * nchords)), tol = 1.0e-8), 1:nchords÷2) # 9.3 μs
    chords = (
        ((-numerator(r), -denominator(r)) for r in rat)..., #
        ((-denominator(r), numerator(r)) for r in rat)..., #
        ((numerator(r), denominator(r)) for r in rat)..., #
        ((denominator(r), -numerator(r)) for r in rat)..., #
    )
    ci = start # Use successive approximation to find the center
    for (pd, nd) in bisectors
        ci = CartesianIndex((measurechord(ci, pd).I .+ measurechord(ci, nd).I) .÷ 2)
    end
    return CartesianIndex[measurechord(ci, s) for s in chords]
end

"""
    rca(img::AbstractArray, search::Function, measure::Function, nchords = 16)
    rca(blob::Vector{Blob})

Search an image for pixels that meet the `search` threshold.  Then using the the rotation chord algorithm
find the approximate center of the particle and draw 2 `nchords` at equally spaced angles from the center
until the `measure` threshold is no longer met.  The resulting set of perimeter points is stored
as a `Vector{CartesianIndex}`.  It is possible that many regios will meet the `search` threshold so the
actual result is a `Vector{Vector{CartesianIndex}}`.

Note: The measure `threshold` should be more permissive than the `search` threshold.

Example

    rca(img, p->p>0.9, p->p>0.7) #
    bs = blob(img, p->p>0.7)
    rca(bs[1])
"""
function rca(img::AbstractArray, search::Function, measure::Function, nchords = 16)
    function bnds(vci::Vector{CartesianIndex})
        l, r, t, b = size(img, 2), 1, size(img, 1), 1
        for ci in vci
            l = ci.I[2] < l ? ci.I[2] : l
            r = ci.I[2] > r ? ci.I[2] : r
            t = ci.I[1] < t ? ci.I[1] : t
            b = ci.I[1] > b ? ci.I[1] : b
        end
        return CartesianIndices((t:b, l:r))
    end
    res, bounds = Vector{Vector{CartesianIndex}}(), Vector{CartesianIndices}()
    skip = false  # Fast skip over columns that have already been checked...
    for ci in CartesianIndices(img)
        if search(img[ci])
            if !(skip || any(map(cis -> ci in cis, bounds)))
                @assert measure(img[ci]) "The 'measure' threshold must be designed such that a pixel that meets the 'search' threshold must also meet the 'measure' threshold."
                r = rca(img, ci, measure, nchords)
                push!(res, r)
                push!(bounds, bnds(r))
                skip = true
            end
        else
            skip = false
        end
    end
    return res
end

function rca(bs::Vector{Blob})
    function offrca(b)  # RCA offset to the original image coordinates
        offci(bbi, r) =
            CartesianIndex[CartesianIndex(map(i -> ci.I[i] + bbi[i].start - 1, eachindex(ci.I))...) for ci in r]
        return [offci(b.bounds.indices, r) for r in rca(b.mask, p -> p, p -> p)]
    end
    return mapreduce(b -> offrca(b), append!, bs)
end

function colorizedimage(chords::Vector{Vector{CartesianIndex}}, img::AbstractArray)
    ends = RGB(1.0, 0.0, 0.0)
    function drawchords(res::Array, chords::Vector{CartesianIndex}, cl::Color)
        lo2 = length(chords) ÷ 2
        for i = 1:lo2
            drawline(pt->res[pt...]=cl, chords[i], chords[i+lo2])
            res[chords[i]]=ends
            res[chords[i+lo2]]=ends
        end
        return res
    end
    colors =
        convert.(
            RGB,
            distinguishable_colors(
                length(chords) + 2,
                Color[RGB(253 / 255, 255 / 255, 255 / 255), RGB(0, 0, 0), RGB(0 / 255, 168 / 255, 45 / 255)],
                transform = deuteranopic,
            )[3:end],
        )
    res = RGB.(img)
    foreach(z -> drawchords(res, z[1], z[2]), zip(chords, colors))
    return res
end

"""
    metrics(chords::Vector{CartesianIndex})::Dict{Symbol,Float64}

Take a `Vector{CartesianIndex}` from `rca(...)` and compute morphology metrics.
"""
function metrics(chords::Vector{CartesianIndex}, scale = 1.0)::Dict{Symbol,Float64}
    lo2, lo4 = length(chords) ÷ 2, length(chords) ÷ 4
    len(a) = sqrt(dot(a, a))
    leni(i) = len(chords[i].I .- chords[i+lo2].I)
    triarea(a, b) = 0.5 * abs(a[1] * b[2] - a[2] * b[1]) # area of triangle with edges a and b
    s, extrema = Moments(), Extrema()
    dmax, imax = 0.0, -1
    for i = 1:lo2
        d = leni(i)
        if d > dmax
            dmax = d
            imax = i
        end
        fit!(extrema, d)
        fit!(s, d)
    end
    perim = len(chords[1].I .- chords[end].I)
    center = (chords[1].I[1], chords[lo4+1].I[2])
    area = triarea(chords[1].I .- center, chords[end].I .- center)
    yext, xext = Extrema(), Extrema()
    fit!(yext, chords[1].I[1])
    fit!(xext, chords[1].I[2])
    for i = 2:length(chords)
        perim += len(chords[i].I .- chords[i-1].I)
        area += triarea(chords[i].I .- center, chords[i-1].I .- center)
        fit!(yext, chords[i].I[1])
        fit!(xext, chords[i].I[2])
    end
    dperp = leni((imax - 1 + lo4) % lo2 + 1)
    return Dict(
        :DMIN => scale * minimum(extrema),
        :DMAX => scale * maximum(extrema),
        :DAVG => scale * Statistics.mean(s),
        :DSIG => scale * Statistics.var(s),
        :DPERP => dperp,
        :ASPECT => dmax / dperp,
        :PERIM => scale * perim,
        :XCENT => center[2],
        :YCENT => center[1],
        :ORIENT2 => 2π * imax / length(chords) - π,
        :XFERET => scale * (maximum(xext) - minimum(xext)),
        :YFERET => scale * (maximum(yext) - minimum(yext)),
        :AREA => scale^2 * area,
    )
end

function area(chords::Vector{CartesianIndex}, scale = 1.0)
    triarea(a, b) = 0.5 * abs(a[1] * b[2] - a[2] * b[1]) # area of triangle with edges a and b
    center = (chords[1].I[1], chords[length(chords)÷4+1].I[2])
    area = triarea(chords[1].I .- center, chords[end].I .- center)
    for i = 2:length(chords)
        area += triarea(chords[i].I .- center, chords[i-1].I .- center)
    end
    return area
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
m = BabyPOMDP()

for step in stepthrough(m, Starve(), "s,a,o", max_steps=1)
    gfx = render(m, step)
    show(stdout, MIME("text/plain"), gfx)
end


text:
# Copyright 2021 Gravifer
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#=
  Author: Gravifer, X. H. Zhu, H. J. Yan, H. Y. Sun
  Date: 2021-04-10
  Description: TBD
=#

module elliptic_grid_gen

"""
To be done.
"""
function inverted_poisson_jacobi_step(arr) end
function inverted_poisson_jacobi_step!(arr) end

"""
To be done.
"""
function inverted_poisson_jacobi_iterate(arr, ε=1e-10) end
function inverted_poisson_jacobi_iterate!(arr, ε=1e-10) end

"""
To be done.
"""
function dirichlet_boundary_conform end

"""
To be done.
"""
function nuemann_boundary_conform end

end # module elliptic_grid_gen


text:
@testset "DirectPowers" begin

   ×(a,b) = Groups.DirectPower(a,b)

   @testset "Constructors" begin
      G = PermutationGroup(3)

      @test Groups.DirectPowerGroup(G,2) isa AbstractAlgebra.Group
      @test G×G isa AbstractAlgebra.Group
      @test Groups.DirectPowerGroup(G,2) isa Groups.DirectPowerGroup{2, Generic.PermGroup{Int64}}

      @test (G×G)×G == DirectPowerGroup(G, 3)
      @test (G×G)×G == (G×G)×G

      GG = DirectPowerGroup(G,2)
      @test (G×G)() isa GroupElem
      @test (G×G)((G(), G())) isa GroupElem
      @test (G×G)([G(), G()]) isa GroupElem

      @test Groups.DirectPowerGroupElem((G(), G())) == (G×G)()
      @test GG(G(), G()) == (G×G)()

      g = perm"(1,2,3)"

      @test GG(g, g^2) isa GroupElem
      @test GG(g, g^2) isa Groups.DirectPowerGroupElem{2, Generic.perm{Int64}}

      h = GG(g,g^2)

      @test h == GG(h)

      @test GG(g, g^2) isa GroupElem
      @test GG(g, g^2) isa Groups.DirectPowerGroupElem

      @test_throws MethodError GG(g,g,g)
      @test GG(g,g^2) == h

      @test h[1] == g
      @test h[2] == g^2
      h = GG(g, G())
      @test h == GG(g, G())
   end

   @testset "Basic arithmetic" begin
      G = PermutationGroup(3)
      GG = G×G
      i = perm"(1,3)"
      g = perm"(1,2,3)"

      h = GG(g,g^2)
      k = GG(g^3, g^2)

      @test h^2 == GG(g^2,g)
      @test h^6 == GG()

      @test h*h == h^2
      @test h*k == GG(g,g)

      @test h*inv(h) == (G×G)()

      w = GG(g,i)*GG(i,g)
      @test w == GG(perm"(1,2)(3)", perm"(2,3)")
      @test w == inv(w)
      @test w^2 == w*w == GG()
   end

   @testset "elem/parent_types" begin
      G = PermutationGroup(3)
      g = perm"(1,2,3)"

      @test elem_type(G×G) == DirectPowerGroupElem{2, elem_type(G)}
      @test elem_type(G×G×G) == DirectPowerGroupElem{3, elem_type(G)}
      @test parent_type(typeof((G×G)(g,g^2))) == Groups.DirectPowerGroup{2, typeof(G)}
      @test parent(DirectPowerGroupElem((g,g^2,g^3))) == DirectPowerGroup(G,3)

      F = AdditiveGroup(GF(13))

      @test elem_type(F×F) ==
         DirectPowerGroupElem{2, Groups.AddGrpElem{AbstractAlgebra.gfelem{Int}}}
      @test parent_type(typeof((F×F)(1,5))) ==
         Groups.DirectPowerGroup{2, Groups.AddGrp{AbstractAlgebra.GFField{Int}}}
      parent((F×F)(1,5)) == DirectPowerGroup(F,2)

      F = MultiplicativeGroup(GF(13))

      @test elem_type(F×F) ==
         DirectPowerGroupElem{2, Groups.MltGrpElem{AbstractAlgebra.gfelem{Int}}}
      @test parent_type(typeof((F×F)(1,5))) ==
         Groups.DirectPowerGroup{2, Groups.MltGrp{AbstractAlgebra.GFField{Int}}}
      parent((F×F)(1,5)) == DirectPowerGroup(F,2)
   end

   @testset "Additive/Multiplicative groups" begin

      R, x = PolynomialRing(QQ, "x")
      F, a = NumberField(x^3 + x + 1, "a")
      G = PermutationGroup(3)

      @testset "MltGrp basic functionality" begin
         Gr = MltGrp(F)
         @test Gr(a) isa MltGrpElem
         g = Gr(a)
         @test deepcopy(g) isa MltGrpElem
         @test inv(g) == Gr(a^-1)
         @test Gr() == Gr(1)
         @test inv(g)*g == Gr()
      end

      @testset "AddGrp basic functionality" begin
         Gr = AddGrp(F)
         @test Gr(a) isa AddGrpElem
         g = Gr(a)
         @test deepcopy(g) isa AddGrpElem
         @test inv(g) == Gr(-a)
         @test Gr() == Gr(0)
         @test inv(g)*g == Gr()
      end
   end

   @testset "Direct Product of Multiplicative Groups" begin

      R, x = PolynomialRing(QQ, "x")
      F, a = NumberField(x^3 + x + 1, "a")
      FF = Groups.DirectPowerGroup(MltGrp(F),2)

      @test FF([a,1]) isa GroupElem
      @test FF([a,1]) isa DirectPowerGroupElem
      @test FF([a,1]) isa DirectPowerGroupElem{2, MltGrpElem{elem_type(F)}}
      @test_throws DomainError FF(1,0)
      @test_throws DomainError FF([0,1])
      @test_throws DomainError FF([1,0])

      @test MltGrp(F) isa AbstractAlgebra.Group
      @test MltGrp(F) isa MultiplicativeGroup
      @test DirectPowerGroup(MltGrp(F), 2) isa AbstractAlgebra.Group
      @test DirectPowerGroup(MltGrp(F), 2) isa DirectPowerGroup{2, MltGrp{typeof(F)}}

      F, a = NumberField(x^3 + x + 1, "a")
      FF = DirectPowerGroup(MltGrp(F), 2)

      @test FF(a,a+1) == FF([a,a+1])
      @test FF([1,a+1])*FF([a,a]) == FF(a,a^2+a)
      x, y = FF([1,a]), FF([a^2,1])
      @test x*y == FF([a^2, a])
      @test inv(x) == FF([1,-a^2-1])

      @test parent(x) == FF
   end

   @testset "Direct Product of Additive Groups" begin

      R, x = PolynomialRing(QQ, "x")
      F, a = NumberField(x^3 + x + 1, "a")

      # Additive Group
      @test AddGrp(F) isa AbstractAlgebra.Group
      @test AddGrp(F) isa AdditiveGroup
      @test DirectPowerGroup(AddGrp(F), 2) isa AbstractAlgebra.Group
      @test DirectPowerGroup(AddGrp(F), 2) isa DirectPowerGroup{2, AddGrp{typeof(F)}}

      FF = DirectPowerGroup(AdditiveGroup(F), 2)

      @test FF([0,a]) isa AbstractAlgebra.GroupElem
      @test FF(F(0),a) isa DirectPowerGroupElem
      @test FF(0,0) isa DirectPowerGroupElem{2, AddGrpElem{elem_type(F)}}

      @test FF(F(1),a+1) == FF([1,a+1])

      @test FF([F(1),a+1])*FF([a,a]) == FF(1+a,2a+1)

      x, y = FF([1,a]), FF([a^2,1])
      @test x*y == FF(a^2+1, a+1)
      @test inv(x) == FF([F(-1),-a])

      @test parent(x) == FF
   end

   @testset "Misc" begin
      F = GF(5)

      FF = DirectPowerGroup(AdditiveGroup(F),2)
      @test order(FF) == 25

      elts = vec(collect(FF))
      @test length(elts) == 25
      @test all([g*inv(g) == FF() for g in elts])
      @test all(inv(g*h) == inv(h)*inv(g) for g in elts for h in elts)

      FF = DirectPowerGroup(MultiplicativeGroup(F), 3)
      @test order(FF) == 64

      elts = vec(collect(FF))
      @test length(elts) == 64
      @test all([g*inv(g) == FF() for g in elts])
      @test all(inv(g*h) == inv(h)*inv(g) for g in elts for h in elts)

      G = PermutationGroup(3)
      GG = Groups.DirectPowerGroup(G,3)
      @test order(GG) == 216

      @test isa(collect(GG), Vector{Groups.DirectPowerGroupElem{3, elem_type(G)}})
      elts = vec(collect(GG))

      @test length(elts) == 216
      @test all([g*inv(g) == GG() for g in elts])
      @test all(inv(g*h) == inv(h)*inv(g) for g in elts for h in elts)
   end
end


text:
module TestTools

using Printf
using Base: rtoldefault
using Test: @test, record, get_testset, Pass, Fail, Error, Broken
using DiffEqBase: DEProblem, set_u!, step!
using LyapunovExponents: LEProblem, dimension, get_integrator, current_state,
    PhaseTangentDynamics

macro test_nothrow(ex)
    quote
        @test begin
            $(esc(ex))
            true
        end
    end
end

function macro_kwargs(kwargs)
    params = map(
        kw -> begin
            @assert kw isa Expr
            @assert kw.head == :(=)
            Expr(:kw, kw.args...)
        end,
        kwargs)
    return Expr(:parameters, params...)
end

macro test_isapprox_elemwise(args...)
    quote
        record(get_testset(),
               @_test_isapprox_elemwise_result($(esc.(args)...)))
    end
end

macro _test_isapprox_elemwise_result(x, y, kwargs...)
    args = [
        esc(x), esc(y),
        QuoteNode(x), QuoteNode(y), QuoteNode(kwargs),
    ]
    Expr(:call,
         test_isapprox_elemwise,
         esc(macro_kwargs(kwargs)),
         args...)
end

function test_isapprox_elemwise(x, y, xexpr, yexpr, orig_kwargs;
                                skip::Bool = false,
                                broken::Bool = false,
                                io = stdout,
                                rtol::Real = rtoldefault(eltype(x), eltype(y)),
                                atol::Real = 0)
    appx = isapprox.(x, y; rtol=rtol, atol=atol)
    ok = all(appx)

    if !ok
        abserr = abs.(x .- y)
        amp = max.(abs.(x), abs.(y))
        th = atol .+ rtol .* amp
        relerr = abserr ./ amp

        printstyled(io, "Not pairwise isapprox", color=:red)
        println(io, "  ", "rtol=", rtol, "  ", "atol=", atol)
        println(io, "x = ", xexpr)
        println(io, "y = ", yexpr)

        columns = [
            :i => 1:length(x),
            :x => x,
            :y => y,
            :relerr => relerr,
            :abserr => abserr,
            :th => th,
        ]
        table = []
        for (name, value) in columns
            push!(table, (name, sprint.(showcompact, value)))
        end
        mark = Dict(true => "✔", false => "!!")
        push!(table, ("ok?", [mark[x] for x in appx]))
        widths = [max(length(string(name)), maximum(length.(c)))
                  for (name, c) in table]

        for ((name, _), width) in zip(table, widths)
            print(io, "  ")
            print(io, lpad(name, width))
        end
        println(io)
        for i in 1:length(x)
            for ((_, c), width) in zip(table, widths)
                print(io, "  ")
                print(io, lpad(c[i], width))
            end
            println(io)
        end
    end

    orig_expr = Expr(:macrocall,
                     (:@test_isapprox_elemwise).args[1],
                     xexpr, yexpr, orig_kwargs...)

    return if skip
        Broken(:skipped, orig_expr)
    elseif broken
        if ok
            Error(:test_unbroken, orig_expr, ok, nothing)
        else
            Broken(:test, orig_expr)
        end
    else
        if ok
            Pass(:test, orig_expr, nothing, ok)
        else
            Fail(:test, orig_expr, nothing, ok)
        end
    end
end

macro display(ex)
    quote
        println($("$ex ="))
        display($(esc(ex)))
        println()
    end
end

function gen_u_list(rng, u_gen, num_u::Int, u::T) where {T}
    u_list = T[]
    for _ in 1:num_u
        push!(u_list, u_gen(rng, eltype(u), size(u)))
    end
    return u_list
end

function rtol_exponents_default(x, y, n=50)
    return range(log10(rtoldefault(eltype(x), eltype(y))), stop=0.5, length=n)
end

function analyze_diff(x, y, ratio=0.1;
                      rtol_exponents = rtol_exponents_default(x, y))
    for n in rtol_exponents
        rtol = 10^n
        bad = mean(@. abs(x - y) > rtol * max(abs(x), abs(y))) / length(x)
        if bad < ratio
            return rtol, bad
        end
    end
    0, 0
end

macro show_diff_info(x, y)
    quote
        rtol, bad = analyze_diff($(esc(x)), $(esc(y)))
        if bad > 0
            @printf("%.1f%% elements are %s ≉ %s with relative tolerance %g\n",
                    bad * 100, $(string(x)), $(string(y)), rtol)
        end
    end
end

function compare_states(u1, u2; verbose=false, kwargs...)
    local ok
    @test ok = isapprox(u1, u2; kwargs...)
    if !ok
        if verbose
            @display u1
            @display u2
            @display u1 - u2
        end
        @show_diff_info(u1, u2)
    end
    ok
end

function test_same_dynamics(f1!, f2!, u_list, p_list, t_list; kwargs...)
    for u in u_list
        for p in p_list
            for t in t_list
                u1 = similar(u)
                u2 = similar(u)
                f1!(u1, copy(u), p, t)
                f2!(u2, copy(u), p, t)
                ok = compare_states(u1, u2; kwargs...)
                if !ok
                    return
                end
            end
        end
    end
end

function test_same_dynamics(p1::P1, p2::P2, u_list::AbstractArray;
                            t_evolve::Real = one(p1.tspan[1]),
                            kwargs...) where {P1 <: DEProblem,
                                              P2 <: DEProblem}

    for u in u_list
        i1 = get_integrator(p1)
        i2 = get_integrator(p2)
        set_u!(i1, copy(u))
        set_u!(i2, copy(u))
        step!(i1, t_evolve, true)
        step!(i2, t_evolve, true)
        u1 = current_state(i1)
        u2 = current_state(i2)
        ok = compare_states(u1, u2; kwargs...)
        if !ok
            return
        end
    end
end

function test_same_dynamics(p1::P1, p2::P2, num_u::Int;
                            evolve = false,
                            rng = Base.GLOBAL_RNG,
                            u_gen = randn,
                            p_list = [p1.p],
                            t_list = p1.tspan[1]:p1.tspan[1] + 1,
                            t_evolve::Real = one(p1.tspan[1]),
                            kwargs...) where {P1 <: DEProblem,
                                              P2 <: DEProblem}
    u_list = gen_u_list(rng, u_gen, num_u, p1.u0)
    if evolve
        test_same_dynamics(p1, p2, u_list;
                           t_evolve = t_evolve,
                           kwargs...)
    else
        test_same_dynamics(p1.f, p2.f, u_list, p_list, t_list; kwargs...)
    end
end

function test_tangent_dynamics_against_autodiff(
        prob::LEProblem, args...;
        kwargs...)
    @assert ! (prob.tangent_prob.f isa PhaseTangentDynamics)
    prob_ad = LEProblem(prob.phase_prob;
                        # t_attr is unused in the following, but it's
                        # a required argument:
                        t_attr = prob.t_attr,
                        # Q0 is required to match the size of
                        # tangent_prob.u0:
                        Q0 = prob.tangent_prob.u0[:, 2:end])
    prob_ad.tangent_prob.u0 .= prob.tangent_prob.u0
    prob_ad.tangent_prob.f :: PhaseTangentDynamics
    test_same_dynamics(prob.tangent_prob,
                       prob_ad.tangent_prob,
                       args...;
                       t_evolve = prob.t_renorm,
                       t_list = 0:prob.t_renorm,
                       kwargs...)
end

end  # module TestTools

using .TestTools


text:
using Random

correlation = [1 0.8; 0.8 1]

@testset "GaussianCopula" begin
    copula = GaussianCopula(correlation)
    @test copula.correlation == correlation
    @test dimensions(copula) == 2

    @testset "sample" begin
        Random.seed!(8128)

        u = sample(copula, 10^5)

        @test round(cor(u[:, 1], u[:, 2]); digits=2) == 0.79

        Random.seed!()
    end

    @testset "to_standard_normal_space" begin
        Random.seed!(8128)

        u = sample(copula, 10^5)
        s = to_standard_normal_space(copula, u)

        @test isapprox(abs(mean(s[:, 1])), 0.0; atol=0.01)
        @test isapprox(abs(mean(s[:, 2])), 0.0; atol=0.01)

        @test isapprox(std(s[:, 1]), 1.0; atol=0.01)
        @test isapprox(std(s[:, 2]), 1.0; atol=0.01)

        @test round(abs(cor(s[:, 1], s[:, 2])); digits=2) == 0.0

        Random.seed!()
    end

    @testset "to_copula_space" begin
        Random.seed!(8128)

        s = randn(10^5, 2)
        u = to_copula_space(copula, s)

        @test round(cor(u[:, 1], u[:, 2]); digits=2) == 0.79

        Random.seed!()
    end
end


text:
using QuantumLab
using Base.Test

if (!isdefined(:indent))
  indent = ""
end

info("$(now())  $(indent)READING:   h2o.xyz -> h2o::Geometry ")
h2o = readGeometryXYZ("h2o.xyz")

# test BasisSetExchange
# as this takes quite some time we primarily only want to do this during continuous integration (travis-ci) and when we haven-t checked this before
if (!isfile("STO-3G.tx93"))
  info("$(now())  $(indent)OBTAINING:   STO-3G.tx93")
  bseEntries = obtainBasisSetExchangeEntries()
  display(bseEntries)
  stoEntry   = computeBasisSetExchangeEntry("sto-3g",bseEntries)[3]
  downloadBasisSetBasisSetExchange(stoEntry,"STO-3G.tx93")
  @test_throws ErrorException computeBasisSetExchangeEntry("NotDefined",bseEntries)
end

info("$(now())  $(indent)READING:   STO-3G.tx93 -> sto3g::BasisSet")
sto3g = readBasisSetTX93("STO-3G.tx93")

info("$(now())  $(indent)COMPUTING:   sto3g,h2o -> bas")
bas = computeBasis(sto3g,h2o)

info("$(now())  $(indent)COMPUTING:   bas,h2o -> matrixOverlap, matrixKinetic, matrixNuclearAttraction, ERIs")
matrixOverlap = computeMatrixOverlap(bas)
matrixKinetic = computeMatrixKinetic(bas)
matrixNuclearAttraction = computeMatrixNuclearAttraction(bas,h2o)
ERIs = computeTensorElectronRepulsionIntegrals(bas)

info("$(now())  $(indent)COMPUTING:   h2o -> matrixSADguess")
matrixSADguess = computeDensityGuessSAD("HF","STO-3G",h2o)

info("$(now())  $(indent)COMPUTING:   sto3g, h2o -> shells, shells_native")
shells = computeBasisShellsLibInt2(sto3g,h2o)
shells_native = computeBasisShells(sto3g,h2o)

info("$(now())  $(indent)COMPUTING:   (HartreeFock)  shells, h2o, matrixSADguess -> density")
density = evaluateSCF(shells,h2o,mean(matrixSADguess),info=false,detailedinfo=false)[3]

info("$(now())  $(indent)COMPUTING:   density, matrixKinetic, matrixNuclearAttraction, ERIs -> matrixFock")
matrixFock = computeMatrixFock(density,matrixKinetic,matrixNuclearAttraction,ERIs)

densityvirt = inv(matrixOverlap) - density

info("$(now())  $(indent)COMPUTING:   shells, density -> S, T, J, K, Fock")
S = computeMatrixOverlap(shells)
T = computeMatrixKinetic(shells)
J = computeMatrixCoulomb(shells,density)
K = computeMatrixExchange(shells,density)
V = computeMatrixNuclearAttraction(shells,h2o)
fock = computeMatrixFock(T,V,J,K)
moenergies = eigvals(Symmetric(fock),Symmetric(S))


text:
module Datasets

using DataDeps
using HTTP
using InternedStrings
using WordTokenizers

using ..Transformers: Container

export Dataset, Train, Dev, Test
export dataset, datafile, get_batch, get_vocab, batched

include("download_utils.jl")
include("./dataset.jl")

include("translate/wmt.jl")
using .WMT

include("translate/iwslt.jl")
using .IWSLT

include("qa/clozetest.jl")
using .ClozeTest

end


text:
module RocketRaceObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "RaceObservable" begin

    println("Testing: race")

    @testset begin
        @test_throws ErrorException race()
    end

    run_testset([
        (
            source = race(of(1), from(1:5)),
            values = @ts([ 1, c ]),
            source_type = Int
        ),
        (
            source = race(of(1.0) |> async(0), from(1:5)),
            values = @ts([ 1, 2, 3, 4, 5, c ]),
            source_type = Union{Float64, Int}
        ),
        (
            source = race(completed(Int), of(1)),
            values = @ts(c),
            source_type = Int
        ),
        (
            source = race(completed(Int), completed(Int)),
            values = @ts(c),
            source_type = Int
        ),
        (
            source = race(faulted(Float64, "err"), completed(Int)),
            values = @ts(e("err")),
            source_type = Union{Float64, Int}
        ),
        (
            source = race(completed(Int), faulted(Float64, "err")),
            values = @ts(c),
            source_type = Union{Int, Float64}
        ),
        (
            source = race(faulted(Float64, "err1"), faulted(Float64, "err2")),
            values = @ts(e("err1")),
            source_type = Float64
        )
    ])

    @testset begin
        s1 = Subject(Int)
        s2 = Subject(Float64)

        source = race(s1, s2)
        values = Vector{Any}()
        subscription = subscribe!(source, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        ))

        @test eltype(source) === Union{Int, Float64}

        @test values == [ ]

        next!(s1, 1)

        @test values == [ 1 ]

        next!(s2, 2.0)

        @test values == [ 1 ]

        complete!(s2)

        @test values == [ 1 ]

        complete!(s1)

        @test values == [ 1, "completed" ]

        unsubscribe!(subscription)
    end

    @testset begin
        s1 = Subject(Int)
        s2 = Subject(Float64)

        source = race(s1, s2)
        values = Vector{Any}()
        subscription = subscribe!(source, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        ))

        @test eltype(source) === Union{Int, Float64}

        @test values == [ ]

        next!(s2, 2.0)

        @test values == [ 2.0 ]

        next!(s1, 1)

        @test values == [ 2.0 ]

        complete!(s1)

        @test values == [ 2.0 ]

        complete!(s2)

        @test values == [ 2.0, "completed" ]

        unsubscribe!(subscription)
    end

    @testset begin
        s1 = Subject(Int)
        s2 = Subject(Float64)

        source = race(s1, s2)
        values = Vector{Any}()
        subscription = subscribe!(source, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        ))

        @test eltype(source) === Union{Int, Float64}

        @test values == [ ]

        unsubscribe!(subscription)

        next!(s1, 1)

        @test values == [ ]

        next!(s2, 2.0)

        @test values == [ ]

        complete!(s1)

        @test values == [ ]

        complete!(s2)

        @test values == [ ]
    end

end

end


text:
module Tutorial8
using FEniCS

T = 0.05            # final time (increase to 5 for full problem)
num_steps = 50   # number of time steps (#increase to 5000 for full stability)
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create mesh
channel = Rectangle(Point([0.0, 0.0]), Point([2.2, 0.41]))
cylinder = Circle(Point([0.2, 0.2]), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)

# Define boundaries
inflow   = "near(x[0], 0)"
outflow  = "near(x[0], 2.2)"
walls    = "near(x[1], 0) || near(x[1], 0.41)"
cylinder = "on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3"

# Define inflow profile
inflow_profile = ("4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)", "0")

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = FeFunction(V)
u_  = FeFunction(V)
p_n = FeFunction(Q)
p_ = FeFunction(Q)

U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
function epsilon(u)
    return sym(nabla_grad(u))
end
# Define stress tensor
function sigma(u, p)
    return 2*mu*epsilon(u) - p*Identity(len(u))
end

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx  + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx + inner(sigma(U, p_n), epsilon(v))*dx  + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[apply(bc,A1) for bc in bcu]
[apply(bc,A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile("navier_stokes_cylinder/velocity.xdmf")
xdmffile_p = XDMFFile("navier_stokes_cylinder/pressure.xdmf")

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries("navier_stokes_cylinder/velocity_series")
timeseries_p = TimeSeries("navier_stokes_cylinder/pressure_series")

# Save mesh to file (for use in reaction_system.py)
File("navier_stokes_cylinder/cylinder.xml.gz",mesh)

global t = 0
#@showprogress 1 "Solving..." for n=0:num_steps #for native julia progress
for n =0:(num_steps-1)

    # Update current time
    global t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [apply(bc,b1) for bc in bcu]
    solve(A1, vector(u_), b1, "bicgstab", "hypre_amg")

    b2 = assemble(L2)
    [apply(bc,b2) for bc in bcp]
    solve(A2, vector(p_), b2, "bicgstab", "hypre_amg")

    b3 = assemble(L3)
    solve(A3, vector(u_), b3, "cg", "sor")

    write(xdmffile_u,u_, t)
    write(xdmffile_p,p_, t)

    # Save nodal values to file
    store(timeseries_u,vector(u_), t)
    store(timeseries_p,vector(p_), t)
    #update values
    assign(u_n,u_)
    assign(p_n,p_)
end
end#module


text:
@doc raw"""
    LinearDiffusionStateModel(A, B; init)

Returns a linear diffusion process hidden state model ``dX_t = A X_t dt + B dW_t`` with appropriately sized matrices ``A`` and ``B``.

Optional argument `init` stands for the initial condition of the process, which is either
* A vector of length `n` for a fixed (deterministic) initial condition
* A `Distributions.Sampleable` type for a random initial condition

If argument `init` is left out, it is set to either
* a multivariate normal distribution with covariance matrix set to the stationary variance, if it exists
* the zero vector
"""
struct LinearDiffusionStateModel{S, F1<:AbstractMatrix{S}, F2<:AbstractMatrix{S}, TI} <: HiddenStateModel{Vector{S}, ContinuousTime}
    drift_matrix::F1
    diffusion_matrix::F2
    init::TI
    function LinearDiffusionStateModel(A::AbstractMatrix{T}, B::AbstractMatrix{T}, init::TI) where TI<:Union{Distributions.Sampleable, AbstractVector{T}} where T
        n       = size(A, 1)
        m       = size(B, 2)
        size(B, 1) == n || throw(DimensionMismatch("Diffusion matrix has incorrect number of rows."))
        if init isa Distributions.Sampleable && eltype(rand(init)) == eltype(A)
            return new{eltype(A), typeof(A), typeof(B), typeof(init)}(A, B, init)
        elseif eltype(init) == eltype(A)
            return new{eltype(A), typeof(A), typeof(B), typeof(init)}(A, B, init)
        else
            error("Error: initial condition and drift and diffusion matrices have incompatible element types.")
        end
    end
end

function LinearDiffusionStateModel(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
    n       = size(A, 1)
    m       = size(B, 2)
    size(B, 1) == n || throw(DimensionMismatch("Diffusion matrix has incorrect number of rows."))

    mu      = zeros(T, n)
    BB      = B*B'
    Sigma   = find_stationary_variance(A, BB)
    if Sigma == zeros(T, n, n)
        return LinearDiffusionStateModel(A, B, mu)
    else
        return LinearDiffusionStateModel(A, B, Distributions.MvNormal(mu, Sigma))
    end
end

function find_stationary_variance(A::AbstractMatrix{T}, BB::AbstractMatrix{T}) where T
    n       = size(A, 1)
    Sigma   = zeros(T, n, n)
    epsilon = 1.
    temp    = zeros(T, n, n)
    count   = 0
    while epsilon > 1e-4 && count < 10000
        temp   .= Sigma
        Sigma  += (BB + A * Sigma + Sigma * A')/100
        epsilon = LinearAlgebra.norm(Sigma - temp)
        count += 1
    end
    if count == 10000 || Sigma != Sigma
        @warn "Iterative procedure to find stationary variance did not converge within 10000 iterations. Most likely the process does not have a stationary variance. Proceeding with initial variance set to zero."
        return zeros(T, n, n)
    else
        return Sigma/2 + Sigma'/2
    end
end

#####################
### BASIC METHODS ###
#####################

drift(model::LinearDiffusionStateModel)             = model.drift_matrix
diffusion(model::LinearDiffusionStateModel)         = model.diffusion_matrix
initial_condition(model::LinearDiffusionStateModel) = model.init
initialize(model::LinearDiffusionStateModel{T, F1, F2, TI}) where {T, F1, F2, TI<:Distributions.Sampleable}      = rand(model.init)
initialize(model::LinearDiffusionStateModel{T, F1, F2, TI}) where TI<:AbstractVector{T} where {T, F1, F2}        = model.init
state_dim(model::LinearDiffusionStateModel)         = size(drift(model), 1)
noise_dim(model::LinearDiffusionStateModel)         = size(diffusion(model), 2)

function drift_function(state_model::LinearDiffusionStateModel)
    A = drift(state_model)
    f(x) = A*x
    return f
end

function diffusion_function(state_model::LinearDiffusionStateModel)
    B = diffusion(state_model)
    g(x) = B
    return g
end

function Base.show(io::IO, ::MIME"text/plain", model::LinearDiffusionStateModel)
    print(io, "Linear diffusion process model for the hidden state
    type of hidden state:                   ", state_dim(model),"-dimensional vector
    number of independent Brownian motions: ", noise_dim(model),"
    initial condition:                      ", model.init isa Distributions.Sampleable ? "random" : "fixed")
end

function Base.show(io::IO, model::LinearDiffusionStateModel)
    print(io, state_dim(model),"-dimensional linear diffusion with ", noise_dim(model), "-dimensional Brownian motion")
end

function Base.show(io::IO, ::MIME"text/plain", model::Type{LinearDiffusionStateModel{S1, S2, S3, S4}}) where {S1, S2, S3, S4}
    print(io, "LinearDiffusionStateModel{", S1, ",", S2, ",", S3, ",", S4, "}")
end

function Base.show(io::IO, model::Type{LinearDiffusionStateModel{S1, S2, S3, S4}}) where {S1, S2, S3, S4}
    print(io, "LinearDiffusionStateModel")
end

######################
### TIME EVOLUTION ###
######################

(model::LinearDiffusionStateModel)(x::AbstractVector, dt) = x + dt*drift(model)*x + sqrt(dt)*diffusion(model)*randn(eltype(x), noise_dim(model))
(model::LinearDiffusionStateModel)(x::AbstractMatrix, dt) = x + dt*drift(model)*x + sqrt(dt)*diffusion(model)*randn(eltype(x), noise_dim(model), size(x, 2))

function propagate!(x::AbstractVector, model::LinearDiffusionStateModel, dt)
    x .+= dt .* drift(model) * x + sqrt(dt) .* diffusion(model) * randn(eltype(x), noise_dim(model))
    return x
end

function propagate!(x::AbstractMatrix, model::LinearDiffusionStateModel, dt)
    x .+= dt .* drift(model) * x + sqrt(dt) .* diffusion(model) * randn(eltype(x), noise_dim(model), size(x, 2))
    return x
end

# convert linear model into general model
DiffusionStateModel(model::LinearDiffusionStateModel) = DiffusionStateModel(drift_function(model), diffusion_function(model), initial_condition(model))

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
mutable struct Schema
    names::Vector{Symbol}       # column names
    types::Vector{Type}      # Julia types of columns

    function Schema(header::Vector{Symbol}, types::Vector{Type})
        @assert length(header) == length(types) "Sizes mismatch: header = $(length(header)), types = $(length(types))"
        @assert length(header) == length(unique(header)) "Column names must be unique"
        new(header, types)
    end
end

mutable struct Table #<: AbstractDataFrame
    schema::Schema
    data::Vector{Any}

    function Table(s::Schema, d::Vector{Any})
        @assert length(s.names) == length(d)
        for i = 1:length(s.types)
            @assert eltype(d[i]) == s.types[i] "Vector types don't match table's schema."
        end
        new(s,d)
    end
end

struct TableRow
    table::Table
    row::Int
end

struct TableRowIterator
    table::Table
end

check_ds_ts_clash(decimal_separator::Char, thousands_separator::Nothing) = nothing
function check_ds_ts_clash(decimal_separator::Char, thousands_separator::Char)
    @assert thousands_separator != decimal_separator "decimal_separator ($decimal_separator) conflicts with thousands_separator ($thousands_separator)."
end

mutable struct CSVFormat
    dlm::Char
    decimal_separator::Char
    thousands_separator::Union{Nothing, Char}
    missing_str::String
    date_format::Dates.DateFormat

    function CSVFormat(dlm::Char, decimal_separator::Char, thousands_separator::Union{Nothing, Char}, missing_str::String, date_format::Dates.DateFormat)
        check_ds_ts_clash(decimal_separator, thousands_separator)
        new(dlm, decimal_separator, thousands_separator, missing_str, date_format)
    end
end

# Create with default values
CSVFormat(; dlm::Char=';', decimal_separator::Char=',', thousands_separator::Union{Nothing, Char}=nothing, missing_str::String="", date_format::Dates.DateFormat=Dates.ISODateFormat) = CSVFormat(dlm, decimal_separator, thousands_separator, missing_str, date_format)
ismissingstr(str, fm::CSVFormat) = str == fm.missing_str


text:
#####################################################
#forwardprop for different nonlinearities

function forwardprop!(net; nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		BLAS.gemv!('N', 1., net.w[i], net.x[i], 0., net.ax[i])
		BLAS.axpy!(1., net.b[i], net.ax[i])
		nonlinearity[i](net.ax[i], net.x[i+1])
	end
end

function forwardprop_no_rec_linearpart!(net::AE_sparse, i::Int64)
	BLAS.gemv!('N', 1., net.w[i], net.x[i], 0., net.ax[i])
	BLAS.axpy!(1., net.b[i], net.ax[i])
end
function forwardprop_no_rec!(net::AE_sparse, lambdas::Array{Float64, 1}; nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		forwardprop_no_rec_linearpart!(net, i)
		nonlinearity[i](net.ax[i],net.x[i+1],lambdas[i])
		#nonlinearity[i](net.ax[i],net.x[i+1],0)
	end
end
function forwardprop_no_rec!(net::AE_sparse; nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		forwardprop_no_rec_linearpart!(net, i)
		nonlinearity[i](net.ax[i],net.x[i+1],net.b[i])
		#nonlinearity[i](net.ax[i],net.x[i+1],zeros(length(net.b[i]))))
	end
end

#see "g"-function in PhD-Thesis of C. Brito or Rozell2008
function activation_Brito!(input,output,lambda)
	for i in 1:length(input)
		output[i] = clamp(input[i]-lambda,0.,Inf64) #thresholded, linear rectifier
	end
end

#Implementation of 1st part of sparse coding algorithm by C. Brito
#Uses non-local inhibitory connections w^T*w
#epsilon: stop condition for sparsity iteration:
#ratio of update magnitude and quantity magnitude
#r: regression, i.e. implementation of smooth integrator with same fixed point with time constant 1/r
#-> iteration should go on at least 1/r for convergence
function forwardprop_Brito!(net::AE_sparse, lambdas::Array{Float64, 1}; epsilon = 1., r = .1, nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		if nonlinearity[i] != lin!
			error("nonlinearity != linear. Must be linear for sparse coding!")
		else
			#while eps > epsilon
			net.x[i+1][:] = 0.
			net.ax[i][:] = 0.
			n,m = size(net.w[i])
			for k in 1:20
				BLAS.gemv!('N',r,net.w[i],net.x[i]-BLAS.gemv('T',net.w[i],net.x[i+1]),1.0-r,net.ax[i])
				activation_Brito!(net.ax[i],net.x[i+1],lambdas[i])
			end
		end
	end
end

#full implementation with plastic lateral inhibitory connections learned by anti-Hebbian rule
#see Brito's Thesis Chapter 2
#abusing a bit the type AE_sparse by
# - saving/integrating past activations of hidden layer in B-array
# - using A-array for lateral inhibitory connections
function _forwardprop_Brito_full!(net::AE_sparse, lambdas::Array{Float64, 1}, iter::Int64;
										nonlinearity = Array{Function, 1}, iterations = 50, r = 1e-1, memory_decay = 1e-2,
										learningrate_inh = 1e-2, no_recurrence = false)
										#This had worked for sure: iterations = 50, r = 1e-1, memory_decay = 1e-2,
										# #with learningrate 1e-2 and learningrate_inh = 1e-2
	if no_recurrence
		forwardprop_no_rec!(net, lambdas; nonlinearity = [activation_Brito!])
	else
		if iter == -1
			one_over_iter = 0
		else
			one_over_iter = 1. / convert(Float64, iter)
		end
		for i in 1:net.nl
			if nonlinearity[i] != lin!
				error("nonlinearity != linear. Must be linear for sparse coding!")
			else
				net.x[i+1][:] .= 0.
				net.ax[i][:] .= 0.
				n,m = size(net.w[i])
				for j in 1:n
					net.A[i][j,j] = 0. #no self-inhibition
				end
				input = BLAS.gemv('N',net.w[i],net.x[i])
				for k in 1:iterations
					net.ax[i] = r*(input-BLAS.gemv('N',net.A[i],net.x[i+1]))+(1.0-r)*net.ax[i]
					activation_Brito!(net.ax[i],net.x[i+1],lambdas[i])
				end
				if iter != -1 #only learn lateral inhibition weights during overall learning procedure, not during generatehiddenreps!
					if one_over_iter > memory_decay #avoids zero-filling at beginning
						net.B[i][1,:] = (1. - one_over_iter)*net.B[i][1,:] + one_over_iter*net.x[i+1]
					else
						net.B[i][1,:] = (1. - memory_decay)*net.B[i][1,:] + memory_decay*net.x[i+1]
					end
					BLAS.ger!(learningrate_inh,net.x[i+1]-net.B[i][1,:],net.x[i+1],net.A[i])
					clamp!(net.A[i],0.,Inf64) #Dale's law
				end
			end
		end
	end
end
function _forwardprop_Földiak!(net::AE_sparse; no_recurrence = false,
								nonlinearity = relu!, dt = 1e-1, epsilon = 1e-4)
	if no_recurrence
		forwardprop_no_rec!(net; nonlinearity = [nonlinearity])
	else
		net.ax[1] .= 0. # u
		net.x[2] .= 0. # a
		scaling_factor = epsilon/dt
		voltage_incr = scaling_factor*norm(net.ax[1])+1 #+1 to make sure loop is entered
		input_without_recurrence = BLAS.gemv('N',net.w[1],net.x[1])
		while norm(voltage_incr) > scaling_factor*norm(net.ax[1])
			voltage_incr = input_without_recurrence - BLAS.gemv('N',net.A[1],net.x[2]) - net.ax[1]
			BLAS.axpy!(dt, voltage_incr, net.ax[1]) # update membrane potential
			nonlinearity(net.ax[1],net.x[2],net.b[1]) # apply activation function
		end
	end
end
function forwardprop_Brito_full!(net::AE_sparse, lambdas::Array{Float64, 1}, iter::Int64;
									nonlinearity = Array{Function, 1}, no_recurrence = false, Földiak_model = false)
	if Földiak_model
		_forwardprop_Földiak!(net; no_recurrence = no_recurrence)
	else
		_forwardprop_Brito_full!(net, lambdas, iter; nonlinearity = nonlinearity, no_recurrence = no_recurrence)
	end
end

#same for generatehiddenreps (not needed to pay attention on zeros in activities at beginning of training:
#iteration number isn't important)
function forwardprop_Brito_full!(net::AE_sparse, lambdas::Array{Float64, 1};
									nonlinearity = Array{Function, 1}, no_recurrence = false, Földiak_model = false)
	if Földiak_model
		_forwardprop_Földiak!(net; no_recurrence = no_recurrence)
	else
		iter = -1 #codes for iteration number isn't important
		_forwardprop_Brito_full!(net, lambdas, iter; nonlinearity = nonlinearity, no_recurrence = no_recurrence)
	end
end

function getrandompatchpositions(net::AE_sparse_patchy; colorimage = false)
	if colorimage
		im_size = Int(sqrt(length(net.x[1])/3))
	else
		im_size = Int(sqrt(length(net.x[1])))
	end
	nh = length(net.x[2])
	return rand(1:im_size - net.patchsize + 1, nh, 2)
end

function distributeinput!(net::AE_sparse_patchy; rand_patch_pos = false,
	colorimage = false)
	if colorimage
		im_size = Int(sqrt(length(net.x[1])/3))
	else
		im_size = Int(sqrt(length(net.x[1])))
	end
	if rand_patch_pos #patch_positions is globally defined in training.jl
		for i in 1:net.n_AE_sparse_patches
			if colorimage
				net.AE_sparse_patches[i].x[1] = reshape(net.x[1], im_size, im_size, 3)[
					global_patch_positions[i,1]:global_patch_positions[i,1]+net.patchsize-1,
					global_patch_positions[i,2]:global_patch_positions[i,2]+net.patchsize-1,:][:]
			else
				net.AE_sparse_patches[i].x[1] =
					reshape(net.x[1],im_size,im_size)[
					global_patch_positions[i,1]:global_patch_positions[i,1]+net.patchsize-1,
					global_patch_positions[i,2]:global_patch_positions[i,2]+net.patchsize-1][:]
			end
		end
	else
		n_patch = Int(sqrt(net.n_AE_sparse_patches))
		ol = net.patchsize - net.stride
		for i in 1:n_patch
			for j in 1:n_patch
				if colorimage
					net.AE_sparse_patches[(i-1)*n_patch+j].x[1] = reshape(net.x[1],im_size,im_size,3)[(i-1)*net.stride+1 : i*net.patchsize-(i-1)*ol,
					   	(j-1)*net.stride+1 : j*net.patchsize-(j-1)*ol,:][:]
				else
					net.AE_sparse_patches[(i-1)*n_patch+j].x[1] =
						reshape(net.x[1],im_size,im_size)[(i-1)*net.stride+1 : i*net.patchsize-(i-1)*ol,
		              (j-1)*net.stride+1 : j*net.patchsize-(j-1)*ol][:]
				end
			end
		end
	end
end

function forwardprop_Brito_full!(net::AE_sparse_patchy, lambdas, iter;
									nonlinearity = nonlinearity, no_recurrence = false,
									Földiak_model = false, rand_patch_pos = false,
									colorimage = false)
	net.x[end] .= 0
	distributeinput!(net; rand_patch_pos = rand_patch_pos, colorimage = colorimage)
	for i in 1:net.n_AE_sparse_patches
		forwardprop_Brito_full!(net.AE_sparse_patches[i], lambdas, iter; nonlinearity = nonlinearity,
			no_recurrence = no_recurrence, Földiak_model = Földiak_model)
		net.x[2][(i-1)*net.n_out_per_population+1:i*net.n_out_per_population] = deepcopy(net.AE_sparse_patches[i].x[end])
	end
end
#for generatehiddenreps
function forwardprop_Brito_full!(net::AE_sparse_patchy, lambdas;
									nonlinearity = nonlinearity, no_recurrence = false,
									Földiak_model = false, rand_patch_pos = false,
									colorimage = false)
	iter = -1
	forwardprop_Brito_full!(net, lambdas, iter; nonlinearity = nonlinearity,
		no_recurrence = no_recurrence, Földiak_model = Földiak_model, rand_patch_pos = rand_patch_pos,
		colorimage = colorimage)
end


text:
"""
Create a System by parsing power-system parameters from a Matpower, PTI, or JSON file and do
some data checks.

# Examples
```julia
sys = parse_standard_files("case_file.m", configpath = "custom_validation.json",
                    bus_name_formatter = x->string(x["name"]*"-"*string(x["index"])),
                    load_name_formatter = x->strip(join(x["source_id"], "_")))
```

"""
function parse_standard_files(file::String; kwargs...)::System

    # function `parse_file` is in pm_io/common.jl
    data = parse_file(file)

    # Check for at least one bus in input file
    if (length(data["bus"]) < 1)
        @error "There are no buses in this file"
    end

    # in pm2ps_parser.jl
    sys = pm2ps_dict(data; kwargs...)

    return sys

end


text:
using NetCDF
using DataFrames
using DataStructures
using ArgParse
using CSV
using ProfileView

include("netcdf_parser.jl")

function parse_command_line()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Path to folder containing input files"
            arg_type = String
            required = true
        "--slice", "-s"
            help = "Start and stop of a subset of the data"
            arg_type = String
            required = true
        "--output", "-o"
            help = "Path to output folder"
            arg_type = String
            required = true
        "--profile", "-p"
            help = "If profiling the script"
            action = :store_true
    end

    return parse_args(s)
end

print_event(name::String) = println("#" * Dates.format(now(), "HH:MM:SS") * "-" * name)

function main(parsed_args)
    df1 = readnetcdf(parsed_args["input"] * "data1.nc")
    df2 = readnetcdf(parsed_args["input"] * "data2.nc")

    print_event("done_read")

    #= PIPELINE =#
    # 1. join the 2 dataframes
    df = join(df1, df2, on=[:longitude, :latitude, :time], kind=:inner)

    # 2. quick preview on the data
    df_head = DataFrames.head(df, 10)

    print_event("done_head")

    df_head[:time] = Dates.format.(df_head[:time], "yyyy-mm-dd")
    CSV.write(parsed_args["output"] * "head.csv", df_head)

    # 3. subset the data
    slice = split(parsed_args["slice"], ":")
    df = df[parse(Int64, slice[1]) + 1:parse(Int64, slice[2]), :]

    # 4. drop rows with null values
    df = df[(df[:tg] .!= -99.99f0) .& (df[:pp] .!= -999.9f0) .& (df[:rr] .!= -999.9f0), :]

    # 5. drop columns
    delete!(df, [:pp_stderr, :rr_stderr])

    # 6. UDF 1: compute absolute difference between max and min
    function compute_abs_maxmin(column_max, column_min)
        abs.(column_max .- column_min)
    end

    df[:abs_diff] = compute_abs_maxmin(df[:tx], df[:tn])

    # 7. explore the data through aggregations
    # TODO: maybe  could use aggregate(df, <no_column>, [min...])?
    function aggregations(df::DataFrame, aggregations::Array{String})
        dict = OrderedDict()
        exclude = [:longitude, :latitude, :time]
        columns_without_time = [k for k in names(df) if !(k in exclude)]
        # add the columns to the dataframe
        dict["columns"] = columns_without_time

        for function_name in aggregations
            func = getfield(Main, Symbol(function_name))
            dict[function_name] = colwise(func, df[columns_without_time])
        end

        DataFrame(dict)
    end

    df_agg = aggregations(df, ["minimum", "maximum", "mean", "std"])

    print_event("done_agg")

    # need to effectively transpose so it's the same as the other pipelines
    rename!(df_agg, Dict(:columns => :agg, :minimum => :min, :maximum => :max))
    res = OrderedDict()
    aggs = names(df_agg)
    res[aggs[1]] = aggs[2:5]
    for row in eachrow(df_agg)
        res[row[1]] = vec(convert(Array, row[2:5]))
    end
    CSV.write(parsed_args["output"] * "agg.csv", DataFrame(res))

    # 8. compute std per month
    # UDF 2: compute custom year+month format
    df[:year_month] = map(x -> Dates.year(x) * 100 + Dates.month(x), df[:time])

    # groupby cols on on and rename the final columns
    on = [:longitude, :latitude, :year_month]
    cols = [:tg, :tn, :tx, :pp, :rr]
    new_cols = [Symbol(string(k) * "_mean") for k in cols]
    grouped_df = by(df[[on; cols]], on, df -> DataFrame(OrderedDict(zip(new_cols, colwise(mean, df[cols])))))
    delete!(grouped_df, [:longitude, :latitude, :year_month])
    grouped_df = DataFrame(OrderedDict("column" => new_cols, "grouped_sum" => colwise(sum, grouped_df)))

    print_event("done_groupby")

    CSV.write(parsed_args["output"] * "grouped.csv", grouped_df)

end

function run()
    parsed_args = parse_command_line()
    if parsed_args["profile"]
        Profile.clear()
        Profile.init(n = 10^9, delay = 0.01)
        @profile main(parsed_args)
        ProfileView.svgwrite(parsed_args["output"] * "profile.svg")
    else
        main(parsed_args)
    end
end

run()

max_stars_repo_name: radujica/data-analysis-pipelines

max_issues_repo_name: radujica/data-analysis-pipelines

max_forks_repo_name: radujica/data-analysis-pipelines


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

using Markdown
using InteractiveUtils

# ╔═╡ d3b0d688-37c5-11eb-2dba-339e12b29c0c
using DataStructures

# ╔═╡ 1a5162ec-376c-11eb-2960-5d2de63d1513
file = open(f->read(f, String), "day_6.txt")

# ╔═╡ 3d63250e-37bc-11eb-33fc-473c93d2a0b6
begin
	#separing by \n\n
	input = eachmatch(r"((?!\n\n)(.|\s))+",file)

	#getting first match of regex
	input = map(x-> x.match, input)
end

# ╔═╡ 63d8cc20-376c-11eb-0a15-e74649299804
begin
	anwser = 0
	for group in input
		all_c = unique(collect(group))
		filter!(l-> l != '\n' ,all_c)

		anwser+= size(all_c)[1]
	end
end

# ╔═╡ 3e522c16-37bd-11eb-0684-b9a6791b529f
anwser

# ╔═╡ fcd1f052-376f-11eb-2e54-611deedd7c74
md"# Question 2"

# ╔═╡ 29d4e6ea-37c5-11eb-33fa-a3bd8a932ad9
begin
	anwser2 = 0

	for group in input
		yes=[x for x in 'a':'z']

		group = split(group, "\n",keepempty=false)
		for person in group
			yes = yes ∩ collect(person)
		end

		anwser2+= size(yes)[1]
	end
end

# ╔═╡ a34e3080-37c5-11eb-12d2-97fc7249d5b2
anwser2

# ╔═╡ Cell order:
# ╠═d3b0d688-37c5-11eb-2dba-339e12b29c0c
# ╠═1a5162ec-376c-11eb-2960-5d2de63d1513
# ╠═3d63250e-37bc-11eb-33fc-473c93d2a0b6
# ╠═63d8cc20-376c-11eb-0a15-e74649299804
# ╠═3e522c16-37bd-11eb-0684-b9a6791b529f
# ╠═fcd1f052-376f-11eb-2e54-611deedd7c74
# ╠═29d4e6ea-37c5-11eb-33fa-a3bd8a932ad9
# ╠═a34e3080-37c5-11eb-12d2-97fc7249d5b2


text:
"""Loads an edge list file as an adjacency matrix."""
function load_adjacency_matrix(path::String)
    edge_list = load_edge_list(path, UInt16)
    one_indexing!(edge_list)
    adj_matrix = load_adjacency_matrix(edge_list)
    return adj_matrix
end

function load_adjacency_matrix(edge_list::Array)
    n_nodes = size(edge_list, 1)
    adjacency_matrix = zeros(Bool, (n_nodes, n_nodes))
    for (n_s, n_t) in eachrow(edge_list)
        adjacency_matrix[n_s, n_t] = true
    end
    return adjacency_matrix
end

function load_adjacency_matrix(io::IO, n_nodes::Int)
    adjacency_matrix = zeros(Bool, (n_nodes, n_nodes))
    for line in eachline(io)
        n_s, n_t = process_line(line, Int)
        adjacency_matrix[n_s + 1, n_t + 1] = true
    end
    return adjacency_matrix
end

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
function DiscretizedAffineEstimator(mhist::MCEBHistogram, kernel::ContinuousUnivariateDistribution)
    DiscretizedAffineEstimator(mhist, x->pdf(kernel, x))
end

"""
    SincKernel(h)

Implements the `SincKernel` with bandwidth `h` to be used for kernel density estimation
through the `KernelDensity.jl` package. The sinc kernel is defined as follows:
```math
K_{\\text{sinc}}(x) = \\frac{\\sin(x)}{\\pi x}
```
It is not typically used for kernel density estimation, because this kernel is not
a density itself. However, it is particularly well suited to deconvolution problems
and estimation of very smooth densities because its Fourier transform is the following:
```math
K^*_{\\text{sinc}}(t) = \\mathbf 1( t \\in [-1,1])
```
"""
struct SincKernel <: ContinuousUnivariateDistribution
   h::Float64 #Bandwidth
end

_default_bandwidth(::Type{SincKernel}, m::Integer) = 1/sqrt(log(m))
SincKernel(; m = 1_000) = SincKernel(default_bandwidth(SincKernel, m))

cf(a::SincKernel, t) = one(Float64)*(-1/a.h <= t <= 1/a.h)

function pdf(a::SincKernel, t)
   if t==zero(Float64)
       return(one(Float64)/pi/a.h)
   else
       return(sin(t/a.h)/pi/t)
   end
end

"""
    DeLaValleePoussinKernel(h)

Implements the `DeLaValleePoussinKernel` with bandwidth `h` to be used for kernel density estimation
through the `KernelDensity.jl` package. The De La Vallée-Poussin kernel is defined as follows:
```math
K_V(x) = \\frac{\\cos(x)-\\cos(2x)}{\\pi x^2}
```
Its use case is similar to the [`SincKernel`](@ref), however it has the advantage of being integrable
(in the Lebesgue sense). Its Fourier transform is the following:
```math
K^*_V(t) = \\begin{cases}
 1, & \\text{ if } x\\in[-1,1] \\\\
 0, &\\text{ if } |t| \\geq 2 \\\\
 2-|t|,& \\text{ if } |t| \\in [1,2]
 \\end{cases}
```
"""
struct DeLaValleePoussinKernel <: ContinuousUnivariateDistribution
   h::Float64 #Bandwidth
end

_default_bandwidth(a::Type{DeLaValleePoussinKernel}, m::Integer) = 1/sqrt(log(m))

function DeLaValleePoussinKernel(; m = 1_000)
     DeLaValleePoussinKernel(default_bandwidth(DeLaValleePoussinKernel, m))
end

function cf(a::DeLaValleePoussinKernel, t)
   if abs(t * a.h) <= 1
       return(one(Float64))
   elseif abs(t * a.h) <= 2
       return(2*one(Float64) - abs(t * a.h))
   else
       return(zero(Float64))
   end
end

function pdf(a::DeLaValleePoussinKernel, t)
   if t==zero(Float64)
       return(3*one(Float64)/2/pi/a.h)
   else
       return(a.h*(cos(t/a.h)-cos(2*t/a.h))/pi/t^2)
   end
end

"""
    KDEInfinityBandOptions(; kernel=DeLaValleePoussinKernel
                             bandwidth = nothing,
                             a_min,
                             a_max,
                             nboot = 1000)

This struct contains hyperparameters that will be used for constructing a neighborhood
of the marginal density. The steps of the method (and corresponding hyperparameter meanings)
are as follows
* First a kernel density estimate ``\\bar{f}`` of the data is fit with `kernel` as the
kernel and `bandwidth` (the default `bandwidth = nothing` corresponds to automatic
bandwidth selection).
* Second, a Poisson bootstrap with `nboot` replication will be used to estimate a ``L_{\\infty}``
neighborhood ``c_m`` of the true density ``f`` which is such that with probability tending to 1:
```math
\\sup_{x \\in [a_{\\text{min}} , a_{\\text{max}}]} | \\bar{f}(x) - f(x)| \\leq c_m
```
Note that the bound is valid from `a_min` to `a_max`.

## Reference:
  > Paul Deheuvels and Gérard Derzko. Asymptotic certainty bands for kernel density
  > estimators based upon a bootstrap resampling scheme. In Statistical models and methods
  > for biomedical and technical systems, pages 171–186. Springer, 2008
"""
Base.@kwdef struct KDEInfinityBandOptions{T<:Real,
                                         S<:Union{T, Nothing}}
   a_min::T
   a_max::T
   npoints::Integer = 4096
   bandwidth::S = nothing
   kernel = DeLaValleePoussinKernel
   nboot::Integer = 1000
   η_infl::T = 1.1
   n_interp::Integer = 25
end

"""
    KDEInfinityBand

The result of running `StatsBase.fit(opt::KDEInfinityBandOptions, Xs)`. Here `opt` is an instance
of [`KDEInfinityBandOptions`](@ref) and `Xs` is a vector of samples distributed according to
a density ```f``.

## Fields:
* `a_min`,`a_max`, `kernel`: These are the same as the fields in `opt:KDEInfinityBandOptions`.
* `C∞`: This is the Poisson Bootstrap point estimate of ``\\sup_{x \\in [a_{\\text{min}} , a_{\\text{max}}]} | \\bar{f}(x) - f(x)|``
* `fitted_kde`: The fitted `KernelDensity` object.
"""
Base.@kwdef struct KDEInfinityBand{T<:Real}
    C∞::T
    a_min::T
    a_max::T
    fitted_kde
    interp_kde
    kernel
    η_infl::T = 1.1
    n_interp::Integer = 25
end

function _default_bandwidth(kernel, Xs::AbstractVector, ::Nothing)
    _default_bandwidth(kernel, Xs)
end

function _default_bandwidth(kernel, Xs::AbstractVector)
    _default_bandwidth(kernel, length(Xs))
end

function _default_bandwidth(kernel, Xs::AbstractVector, h::Number)
    h
end

function StatsBase.fit(opt::KDEInfinityBandOptions, Xs)
    a_min = opt.a_min
    a_max = opt.a_max
    npoints = opt.npoints
    kernel = opt.kernel
    bandwidth = _default_bandwidth(kernel, Xs, opt.bandwidth)
    nboot = opt.nboot
    tmp_res = certainty_banded_KDE(Xs, a_min, a_max; npoints = npoints,
                         kernel = kernel, bandwidth = bandwidth,
                         nboot = nboot)
    tmp_res = @set tmp_res.η_infl = opt.η_infl
    tmp_res = @set tmp_res.n_interp = opt.n_interp
    tmp_res
end

function certainty_banded_KDE(Xs, a_min, a_max;
                        npoints = 4096,
                        kernel=DeLaValleePoussinKernel,
                        bandwidth = _default_bandwidth(kernel, length(Xs)),
                        nboot = 1_000)

    kernel = kernel(bandwidth)
    h = bandwidth
    m = length(Xs)

    lo, hi = extrema(Xs)
    lo_kde, hi_kde = min(a_min, lo - 6*h), max(a_max, hi + 6*h)
    midpts = range(lo_kde, hi_kde; length = npoints)

    fitted_kde = kde(Xs, KernelDensity.UniformWeights(m), midpts, kernel)
    interp_kde = InterpKDE(fitted_kde)
    midpts_idx = findall( (midpts .>= a_min) .& (midpts .<= a_max) )
    C∞_boot = Vector{Float64}(undef, nboot)

    for k =1:nboot
        # Poisson bootstrap to estimate certainty band
        Z_pois = rand(Poisson(1), m)
        ws =  Weights(Z_pois/sum(Z_pois))
        f_kde_pois =  kde(Xs, ws, midpts, kernel)
        C∞_boot[k] = maximum(abs.(fitted_kde.density[midpts_idx] .- f_kde_pois.density[midpts_idx]))
    end

    C∞ = median(C∞_boot)

    KDEInfinityBand(C∞=C∞,
                    a_min=a_min,
                    a_max=a_max,
                    fitted_kde=fitted_kde,
                    interp_kde=interp_kde,
                    kernel=kernel)
end

function certainty_banded_KDE(Xs, Zs_discr::DiscretizedStandardNormalSamples; kwargs...)
    a_min, a_max = extrema(Zs_discr.mhist.grid)
    certainty_banded_KDE(Xs, a_min, a_max; kwargs...)
end

# plotting of infinite band

@recipe function f(ctband::KDEInfinityBand;
                        show_bands=true)
    y = ctband.fitted_kde.density
    x = ctband.fitted_kde.x
    xlim --> (ctband.a_min, ctband.a_max)
    ylab --> "Density"

    seriestype  -->  :path
    linewidth --> 2
    color --> :purple

    infty_bound = ctband.C∞ * ctband.η_infl

    if show_bands
        _ys_lower, _ys_upper = density_bands_to_ribbons(y, infty_bound)
        ribbons --> (_ys_lower, _ys_lower)
        fillalpha --> 0.2
    end

    xguide --> L"x"
    x,y
end

"""
    set_neighborhood(Zs_discr::DiscretizedStandardNormalSamples,
                          fkde::KDEInfinityBand)

Return a [`DiscretizedStandardNormalSamples`](@ref) instance identical to `Zs_discr` but with
`f_min`,`f_max`,`var_proxy` set based on the fit from `fkde`
(cf.[`KDEInfinityBand`](@ref)).

"""
function set_neighborhood(Zs_discr::DiscretizedStandardNormalSamples,
                          fkde::KDEInfinityBand)

    @unpack n_interp, η_infl = fkde
    grid = Zs_discr.mhist.grid
    fkde_interp = fkde.interp_kde.itp
    n = nobs(Zs_discr)

    f_max = zeros(length(grid) + 1)
    f_min = zeros(length(grid) + 1)
    var_proxy =  zeros(length(grid) + 1)

    C∞ = η_infl * fkde.C∞
    # fill in intermediate points first
    for i in 2:length(grid)
        x_left = grid[i-1]
        x_right = grid[i]
        h = x_right - x_left  # step(grid)...
        interp_res = fkde_interp(range(x_left, x_right, length=n_interp))
        f_max[i] = maximum(interp_res)*h + C∞*h
        f_min[i] = max( minimum(interp_res)*h - C∞*h, 0.0)
        var_proxy[i] = max(maximum(interp_res)*h, C∞*h)
    end

    #fill in first and last element

    # DKW constant with log(n) instead of log(2/alpha)
    C_dkw = sqrt(log(n)/2n)
    f_max[1] = pdf(Zs_discr.mhist)[1] + C_dkw
    f_min[1] = max(pdf(Zs_discr.mhist)[1] - C_dkw, 0.0)
    var_proxy[1] = f_max[1]
    f_max[end] = pdf(Zs_discr.mhist)[end] + C_dkw
    f_min[end] = max(pdf(Zs_discr.mhist)[end] - C_dkw, 0.0)
    var_proxy[end] = f_max[end]

    Zs_discr = @set Zs_discr.var_proxy = var_proxy
    Zs_discr = @set Zs_discr.f_min = f_min
    Zs_discr = @set Zs_discr.f_max = f_max
    Zs_discr
end


text:
using TensorBoardLogger #import the TensorBoardLogger package
using Logging #import Logging package

logger = TBLogger("textlogs", tb_append) #create tensorboard logger

################log scalars example: String################
#using logger interface
with_logger(logger) do
    @info "string/loggerinterface" str = "A computer once beat me at chess, but it was no match for me at kick boxing. *— Emo Philips*"
end
#using explicit function interface
str =
"
All parts should go together without forcing.
You must remember that the parts you are reassembling were disassembled by you.
Therefore, if you can’t get them together again, there must be a reason.
By all means, do not use a hammer.”
*— IBM Manual, 1925*
"
log_text(logger, "string/explicitinterface", str, step = 0)

################log scalar example: y = matrix################
squares =
[["Number" "Square"];
[1 1]; [2 4]; [3 9]; [4 16]; [5 25]; [6 36]; [7 49]; [8 64]; [9 81]; [10 100]]

#using explicit function interface
log_text(logger, "table", squares, step = 0)


text:
using Makie
using MultipleViewGeometry, Test, Random
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleDraw
using MultipleViewGeometry.ModuleMove
using LinearAlgebra
using StaticArrays
using GeometryTypes

# Fix random seed.
Random.seed!(1234)
# Construct two camera matrices and parametrise two planar surfaces.
f = 50
image_width = 640 / 10
image_height = 480 / 10
𝐊₁ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]
𝐑₁ = SMatrix{3,3,Float64,9}(rotxyz(0, 25*(pi/180), 0))
𝐭₁ = [-30.0, 0.0, -5.0]

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

𝐗₁ = [30,15,60]
𝐗₂ = [-30,15,60]
𝐗₃ = [30,-15,60]
𝐗₄ = [-30,-15,60]

scale = 20.0f0
x = Vec3f0(0); baselen = 0.2f0 * scale ; dirlen = 1f0 * scale
# create an array of differently colored boxes in the direction of the 3 axes
rectangles = [
    (HyperRectangle(Vec3f0(x), Vec3f0(dirlen, baselen, baselen)), RGBAf0(1,0,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, dirlen, baselen)), RGBAf0(0,1,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, baselen, dirlen)), RGBAf0(0,0,1,1))
]
meshes = map(GLNormalMesh, rectangles)

scene = mesh(merge(meshes))
scatter!(scene, 𝐗₁[1,:], 𝐗₁[2,:], 𝐗₁[3,:], markersize = 3, color = :red)
scatter!(scene, 𝐗₂[1,:], 𝐗₂[2,:], 𝐗₂[3,:], markersize = 3, color = :green)
scatter!(scene, 𝐗₃[1,:], 𝐗₃[2,:], 𝐗₃[3,:], markersize = 3, color = :black)
scatter!(scene, 𝐗₄[1,:], 𝐗₄[2,:], 𝐗₄[3,:], markersize = 3, color = :blue)

# Convention I: world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
# Convention I: camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
# Convention I: picture_basis = (Point(0.0, 0.0), Vec(1.0, 0.0), Vec(0.0, 1.0))
picture_basis = (Point(0.0, 0.0), Vec(-1.0, 0.0), Vec(0.0, -1.0))

camera₁ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)

relocate!(camera₁, 𝐑₁, 𝐭₁)
draw!(camera₁, scene)

𝐑₁′, 𝐭₁′ = ascertain_pose(camera₁, world_basis... )
𝐊₁′ = obtain_intrinsics(camera₁, CartesianSystem())
#𝐊₁′ = SMatrix{3,3,Float64,9}(f, 0.0, 0.0, 0.0, f, 0.0, 0.0, 0.0 , 1)

𝐏₁ = construct(ProjectionMatrix(), 𝐊₁′, 𝐑₁′, 𝐭₁′)

# Red
ℳ₁ = project(camera₁, 𝐏₁, [SVector(𝐗₁...)])
Base.display(ℳ₁)

# Green
ℳ₂ = project(camera₁, 𝐏₁, [SVector(𝐗₂...)])
Base.display(ℳ₂)

# Black
ℳ₃ = project(camera₁, 𝐏₁, [SVector(𝐗₃...)])
Base.display(ℳ₃)

# Blue
ℳ₄ = project(camera₁, 𝐏₁, [SVector(𝐗₄...)])
Base.display(ℳ₄)

𝐗₁′  = Array([𝐑₁′' -𝐑₁′'*𝐭₁′] * vcat(𝐗₁,1))
𝐗₂′  = Array([𝐑₁′' -𝐑₁′'*𝐭₁′] * vcat(𝐗₂,1))
𝐗₃′  = Array([𝐑₁′' -𝐑₁′'*𝐭₁′] * vcat(𝐗₃,1))
𝐗₄′  = Array([𝐑₁′' -𝐑₁′'*𝐭₁′] * vcat(𝐗₄,1))

scene₂ = mesh(merge(meshes))
scatter!(scene₂, 𝐗₁′[1,:], 𝐗₁′[2,:], 𝐗₁′[3,:], markersize = 3, color = :red)
scatter!(scene₂, 𝐗₂′[1,:], 𝐗₂′[2,:], 𝐗₂′[3,:], markersize = 3, color = :green)
scatter!(scene₂, 𝐗₃′[1,:], 𝐗₃′[2,:], 𝐗₃′[3,:], markersize = 3, color = :black)
scatter!(scene₂, 𝐗₄′[1,:], 𝐗₄′[2,:], 𝐗₄′[3,:], markersize = 3, color = :blue)

relocate!(camera₁, 𝐑₁′, -𝐭₁′)
draw!(camera₁, scene₂)

# Red
hom⁻¹(𝐊₁′*𝐗₁′)
Base.display(hom⁻¹(𝐊₁′*𝐗₁′))

# Green
hom⁻¹(𝐊₁′*𝐗₂′)
Base.display(hom⁻¹(𝐊₁′*𝐗₂′))

# Black
hom⁻¹(𝐊₁′*𝐗₃′)
Base.display(hom⁻¹(𝐊₁′*𝐗₃′))

# Blue
hom⁻¹(𝐊₁′*𝐗₄′)
Base.display(hom⁻¹(𝐊₁′*𝐗₄′))

# Visualize the set of corresponding points
scene = Scene()
M₁ = reshape(reinterpret(Float64,ℳ₁),(2,length(ℳ₁)))
M₂ = reshape(reinterpret(Float64,ℳ₂),(2,length(ℳ₂)))
M₃ = reshape(reinterpret(Float64,ℳ₃),(2,length(ℳ₃)))
M₄ = reshape(reinterpret(Float64,ℳ₄),(2,length(ℳ₄)))

scatter!(scene,M₁[1,:], M₁[2,:], markersize = 1, color = :red, limits = FRect(0, 0, 64, 48.0))
scatter!(scene,M₂[1,:], M₂[2,:], markersize = 1, color = :green, limits = FRect(0, 0, 64, 48.0))
scatter!(scene,M₃[1,:], M₃[2,:], markersize = 1, color = :black, limits = FRect(0, 0, 64, 48.0))
scatter!(scene,M₄[1,:], M₄[2,:], markersize = 1, color = :blue, limits = FRect(0, 0, 64, 48.0))

# scatter!(scene, 𝐗₁′[1,:], 𝐗₁′[2,:], 𝐗₁′[3,:], markersize = 3, color = :red)
# scatter!(scene, 𝐗₂′[1,:], 𝐗₂′[2,:], 𝐗₂′[3,:], markersize = 3, color = :green)
# scatter!(scene, 𝐗₃′[1,:], 𝐗₃′[2,:], 𝐗₃′[3,:], markersize = 3, color = :black)
# scatter!(scene, 𝐗₄′[1,:], 𝐗₄′[2,:], 𝐗₄′[3,:], markersize = 3, color = :blue)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^2*sin(a+b*x+c*x^2), x) == :(-1//2*x*c^-1*cos(a+b*x+c*x^2)+1//4*b*c^-2*cos(a+b*x+c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(x*sin(a+b*x+c*x^2), x) == :(-1//2*c^-1*cos(a+b*x+c*x^2)+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(sin(a+b*x+c*x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*c^-1//2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+(1/2)*2^(1/2)*pi^(1/2)*c^-1//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(x^-2*sin(a+b*x+c*x^2)+-1*b*x^-1*cos(a+b*x+c*x^2), x) == :(-1*x^-1*sin(a+b*x+c*x^2)+2^(1/2)*pi^(1/2)*c^(1/2)*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1*2^(1/2)*pi^(1/2)*c^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(x^2*sin(a+b*x+-1*c*x^2), x) == :((1/2)*x*c^-1*cos(a+b*x+-1*c*x^2)+1//4*b*c^-2*cos(a+b*x+-1*c*x^2)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+1//4*2^(1/2)*pi^(1/2)*c^-3//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1)+-1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1)+1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*cos(a+1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x)))
@test integrate(x*sin(a+b*x+-1*c*x^2), x) == :((1/2)*c^-1*cos(a+b*x+-1*c*x^2)+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1)+1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x)))
@test integrate(sin(a+b*x+-1*c*x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*c^-1//2*cos(a+1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1//2*2^(1/2)*pi^(1/2)*c^-1//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(x^-2*sin(a+b*x+-1*c*x^2)+-1*b*x^-1*cos(a+b*x+-1*c*x^2), x) == :(-1*x^-1*sin(a+b*x+-1*c*x^2)+2^(1/2)*pi^(1/2)*c^(1/2)*cos(a+1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+2^(1/2)*pi^(1/2)*c^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(x^2*sin(1//4+x+x^2), x) == :(1//4*cos(1//4+x+x^2)+-1//2*x*cos(1//4+x+x^2)+1//4*2^(1/2)*pi^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*(1+2x))+1//8*2^(1/2)*pi^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(x*sin(1//4+x+x^2), x) == :(-1//2*cos(1//4+x+x^2)+-1//4*2^(1/2)*pi^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(sin(1//4+x+x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(x^2*sin(a+b*x+c*x^2)^2, x) == :(1//6*x^3+-1//8*x*c^-1*sin(2a+2*b*x+2*c*x^2)+1//16*b*c^-2*sin(2a+2*b*x+2*c*x^2)+1//16*pi^(1/2)*c^-3//2*cos(2a+-1//2*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))+1//16*pi^(1/2)*c^-3//2*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+-1//16*pi^(1/2)*b^2*c^-5//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//16*pi^(1/2)*b^2*c^-5//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))
@test integrate(x*sin(a+b*x+c*x^2)^2, x) == :(1//4*x^2+-1//8*c^-1*sin(2a+2*b*x+2*c*x^2)+-1//8*b*pi^(1/2)*c^-3//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//8*b*pi^(1/2)*c^-3//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate(sin(a+b*x+c*x^2)^2, x) == :((1/2)*x+-1//4*pi^(1/2)*c^-1//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//4*pi^(1/2)*c^-1//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))
@test integrate(x^2*sin(a+b*x+-1*c*x^2)^2, x) == :(1//6*x^3+1//8*x*c^-1*sin(2a+-2*c*x^2+2*b*x)+1//16*b*c^-2*sin(2a+-2*c*x^2+2*b*x)+-1//16*pi^(1/2)*c^-3//2*cos(2a+(1/2)*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))+1//16*pi^(1/2)*c^-3//2*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1)+1//16*pi^(1/2)*b^2*c^-5//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+1//16*pi^(1/2)*b^2*c^-5//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(x*sin(a+b*x+-1*c*x^2)^2, x) == :(1//4*x^2+1//8*c^-1*sin(2a+-2*c*x^2+2*b*x)+1//8*b*pi^(1/2)*c^-3//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+1//8*b*pi^(1/2)*c^-3//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(sin(a+b*x+-1*c*x^2)^2, x) == :((1/2)*x+1//4*pi^(1/2)*c^-1//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+1//4*pi^(1/2)*c^-1//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(x^2*sin(1//4+x+x^2)^2, x) == :(1//6*x^3+1//16*sin(1/2+2x+2*x^2)+-1//8*x*sin(1/2+2x+2*x^2)+-1//16*pi^(1/2)*fresnelc(pi^-1//2*(1+2x))+1//16*pi^(1/2)*fresnels(pi^-1//2*(1+2x)))
@test integrate(x*sin(1//4+x+x^2)^2, x) == :(-1//8*sin(1/2+2x+2*x^2)+1//4*x^2+1//8*pi^(1/2)*fresnelc(pi^-1//2*(1+2x)))
@test integrate(sin(1//4+x+x^2)^2, x) == :((1/2)*x+-1//4*pi^(1/2)*fresnelc(pi^-1//2*(1+2x)))
@test integrate((d+e*x)^2*sin(a+b*x+c*x^2), x) == :(-1//2*e*c^-1*(d+e*x)*cos(a+b*x+c*x^2)+-1//4*e*c^-2*(-1*b*e+2*c*d)*cos(a+b*x+c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*e^2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*e^2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//8*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//8*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate((d+e*x)*sin(a+b*x+c*x^2), x) == :(-1//2*e*c^-1*cos(a+b*x+c*x^2)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//4*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(sin(a+b*x+c*x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*c^-1//2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+(1/2)*2^(1/2)*pi^(1/2)*c^-1//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate((d+e*x)^2*sin(a+b*x+c*x^2)^2, x) == :(1//6*e^-1*(d+e*x)^3+-1//8*e*c^-1*(d+e*x)*sin(2a+2*b*x+2*c*x^2)+-1//16*e*c^-2*(-1*b*e+2*c*d)*sin(2a+2*b*x+2*c*x^2)+-1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//16*pi^(1/2)*c^-3//2*e^2*cos(2a+-1//2*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))+1//16*pi^(1/2)*c^-3//2*e^2*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))
@test integrate(sin(a+b*x+c*x^2)^2*(d+e*x), x) == :(1//4*e^-1*(d+e*x)^2+-1//8*e*c^-1*sin(2a+2*b*x+2*c*x^2)+-1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))
@test integrate(sin(a+b*x+c*x^2)^2, x) == :((1/2)*x+-1//4*pi^(1/2)*c^-1//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//4*pi^(1/2)*c^-1//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.13 (d+e x)^m sin(a+b x+c x^2)^n.jl


text:
# Modelo de separador de corrente
type Splitter
	Splitter()=begin
		new(
			Fracao(Dict{Symbol,Any}(
				:Brief=>"Fracao de massa para a segunda corrente"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>0,
				:PosY=>0.5
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>1,
				:PosY=>0.25
			)),
			Corrente(Dict{Symbol,Any}(
				:Brief=>"Segunda corrente",
				:Symbol=>"_{outx}",
				:PosX=>1,
				:PosY=>0.75
			)),
			[
				:(Fout.P = Fin.P),
				:(Fout.T = Fin.T),
				:(Fout.S = Fin.S),
				:(Fout.H = Fin.H),
				:(Fouts.P = Fin.P),
				:(Fouts.T = Fin.T),
				:(Fouts.S = Fin.S),
				:(Fouts.H = Fin.H),
				:(Fouts.F = Fin.F * y),
				:(Fout.F = Fin.F - Fouts.F),
			],
			[
				"","","","","","","","","","",
			],
			[:y,:Fin,:Fout,:Fouts,]
		)
	end
	y::Fracao
	Fin::Corrente
	Fout::Corrente
	Fouts::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Splitter
function setEquationFlow(in::Splitter)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
end
function atributes(in::Splitter,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/splitter"
	drive!(fields,_)
	return fields
end
Splitter(_::Dict{Symbol,Any})=begin
	newModel=Splitter()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using LinearAlgebra, FastSolversForWeightedTV, CUDA, Flux, Test, BenchmarkTools
CUDA.allowscalar(false)
include("./test_utils.jl")

# flag_gpu = true
flag_gpu = false

# Random data
T = Complex{Float64}
n = (256, 256)
y = randn(T, n..., 2); flag_gpu && (y = y |> gpu)
step = real(T)(1e-5)
rtol = real(T)(1e1*step)

# Norm (2-D)
for g = [mixed_norm(2,2,2; T=T), mixed_norm(2,2,1; T=T), mixed_norm(2,2,Inf; T=T)]

    # Proxy
    λ = real(T)(0.5)*norm(y)^2/g(y)
    x = proxy(y, λ, g)

    # Gradient test
    fun = proxy_objfun(λ, g)
    test_grad(fun, y; step=step, rtol=rtol)

    # Projection test
    ε = real(T)(0.1)*g(y)
    x = project(y, ε, g)
    @test g(x) ≈ ε rtol=rtol

    ## Gradient test
    fun = proj_objfun(ε, g)
    test_grad(fun, y; step=step, rtol=rtol)

end

# Random data
n = (256, 256, 256)
y = randn(T, n..., 3); flag_gpu && (y = y |> gpu)

# Norm (3-D)
for g = [mixed_norm(3,2,2; T=T), mixed_norm(3,2,1; T=T), mixed_norm(3,2,Inf; T=T)]

    # Proxy
    λ = real(T)(0.5)*norm(y)^2/g(y)
    x = proxy(y, λ, g)

    # Gradient test
    fun = proxy_objfun(λ, g)
    test_grad(fun, y; step=step, rtol=rtol)

    # Projection test
    ε = real(T)(0.1)*g(y)
    x = project(y, ε, g)
    @test g(x) ≈ ε rtol=rtol

    ## Gradient test
    fun = proj_objfun(ε, g)
    test_grad(fun, y; step=step, rtol=rtol)

end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
module FactorGraphs
using LightGraphs

import Base:
    eltype, show, ==, Pair,
    Tuple, copy, length, size,
    issubset, zero, getindex
# import Random:
#     randstring, seed!

import LightGraphs:
    AbstractGraph, src, dst, edgetype, nv,
    ne, vertices, edges, is_directed,
    add_vertex!, add_edge!, rem_vertex!, rem_edge!,
    has_vertex, has_edge, inneighbors, outneighbors,
    weights, indegree, outdegree, degree,
    induced_subgraph,
    loadgraph, savegraph, AbstractGraphFormat,
    reverse

import LightGraphs.SimpleGraphs:
    AbstractSimpleGraph, SimpleGraph, SimpleDiGraph,
    SimpleEdge, fadj, badj

export
    FactorGraph
    # addVariable!,
    # addFactor!,
    # LightBayesGraph,
    # filter_edges,
    # filter_vertices,
    # reverse

# import DistributedFactorGraphs: DFGNode
# const AbstractNodeType = DFGNode
import DistributedFactorGraphs: AbstractDFGVariable, AbstractDFGFactor
const AbstractVariableType = AbstractDFGVariable
const AbstractFactorType = AbstractDFGFactor

include("BiMaps.jl")

struct FactorGraph{T <: Integer,V <: AbstractVariableType, F <: AbstractFactorType} <: AbstractGraph{T}
    graph::SimpleGraph{T}
    labels::BiDictMap{T}
    variables::Dict{Symbol,V}
    factors::Dict{Symbol,F}
end

function FactorGraph{T, V, F}(nv::Int=100, nf::Int=100) where {T <: Integer, V <: AbstractVariableType, F <: AbstractFactorType}
    fadjlist = Vector{Vector{T}}()
        sizehint!(fadjlist, nv + nf)
    g = SimpleGraph{T}(0, fadjlist)
    labels = BiDictMap{T}(sizehint=nv+nf)
    variables = Dict{Symbol,V}()
        sizehint!(variables, nv)
    factors = Dict{Symbol,F}()
        sizehint!(factors, nf)
    return FactorGraph{T, V, F}(g, labels, variables, factors)
end

# fg = FactorGraph{Int, AbstractVariableType, AbstractFactorType}()

FactorGraph() = FactorGraph{Int, AbstractVariableType, AbstractFactorType}()
# FactorGraph{V,F}() where {V <: AbstractVariableType, F <: AbstractFactorType} = FactorGraph{Int, V, F}()

function show(io::IO, ::MIME"text/plain", g::FactorGraph)
    dir = is_directed(g) ? "directed" : "undirected"
    print(io, "{$(nv(g)), $(ne(g))} $dir $(eltype(g)) $(typeof(g))")
end

@inline fadj(g::FactorGraph, x...) = fadj(g.graph, x...)
@inline badj(g::FactorGraph, x...) = badj(g.graph, x...)

eltype(g::FactorGraph) = eltype(g.graph)
edgetype(g::FactorGraph) = edgetype(g.graph)
nv(g::FactorGraph) = nv(g.graph)
vertices(g::FactorGraph) = vertices(g.graph)

ne(g::FactorGraph) = ne(g.graph)
edges(g::FactorGraph) = edges(g.graph)

has_vertex(g::FactorGraph, x...) = has_vertex(g.graph, x...)
@inline has_edge(g::FactorGraph, x...) = has_edge(g.graph, x...)

inneighbors(g::FactorGraph, v::Integer) = inneighbors(g.graph, v)
outneighbors(g::FactorGraph, v::Integer) = fadj(g.graph, v)

is_directed(::Type{FactorGraph}) = false
is_directed(::Type{FactorGraph{T,V,F}}) where {T,V,F} = false
is_directed(g::FactorGraph) = false

zero(g::FactorGraph{T,V,F}) where {T,V,F} = FactorGraph{T,V,F}(0,0)

# TODO issubset(g::T, h::T) where T <: FactorGraph = issubset(g.graph, h.graph)

@inline add_edge!(g::FactorGraph, x...) = add_edge!(g.graph, x...)

@inline rem_edge!(g::FactorGraph, x...) = rem_edge!(g.graph, x...)

function addVariable!(g::FactorGraph{T, V, F}, variable::V) where {T, V, F}

    label = variable.label

    haskey(g.labels, label) && (@error "Label already in fg"; return false) #TODO debug error or exception?

    add_vertex!(g.graph) || return false

    g.labels[nv(g.graph)] = label

    push!(g.variables, label=>variable)

    return true
end

function addFactor!(g::FactorGraph{T, V, F}, variableLabels::Vector{Symbol}, factor::F)::Bool where {T, V, F}

    haskey(g.labels, factor.label) && (@error "Label $(factor.label) already in fg"; return false)

    for vlabel in variableLabels
        !haskey(g.labels, vlabel) && (@error "Variable '$(vlabel)' not found in graph when creating Factor '$(factor.label)'"; return false) #TODO debug error or exception?
    end

    add_vertex!(g.graph) || return false

    g.labels[nv(g.graph)] = factor.label

    push!(g.factors, factor.label=>factor)

    # add the edges
    for vlabel in variableLabels
        add_edge!(g, g.labels[vlabel], nv(g.graph)) || return false
    end

    return true
end

function rem_vertex!(g::FactorGraph{T,V,F}, v::Integer) where {T,V,F}
    v in vertices(g) || return false
    lastv = nv(g)

    rem_vertex!(g.graph, v) || return false

    label = g.labels[v]
    delete!(g.variables, label)
    delete!(g.factors, label)

    if v != lastv
        g.labels[v] = g.labels[lastv] #lastSym
    else
        delete!(g.labels, v)
    end

    return true
end

end