text:
export ArduinoSurveillanceUnitInternalTemp, ArduinoSurveillanceUnitInternalTempParams

Base.@kwdef struct ArduinoSurveillanceUnitInternalTempParams <: DeviceParams
  portAdress::String
  commandStart::String = "!"
  commandEnd::String = "*"

  pause_ms::Int = 30
  timeout_ms::Int = 500
  delim::String = "#"
  baudrate::Integer = 9600
  ndatabits::Integer = 8
  parity::SPParity = SP_PARITY_NONE
  nstopbits::Integer = 1
end

ArduinoSurveillanceUnitInternalTempParams(dict::Dict) = params_from_dict(ArduinoSurveillanceUnitInternalTempParams, dict)
Base.@kwdef mutable struct ArduinoSurveillanceUnitInternalTemp <: ArduinoSurveillanceUnit
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::ArduinoSurveillanceUnitInternalTempParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String,Union{Device,Missing}}

  ard::Union{SimpleArduino, Nothing} = nothing # Use composition as multiple inheritance is not supported
end

neededDependencies(::ArduinoSurveillanceUnitInternalTemp) = []
optionalDependencies(::ArduinoSurveillanceUnitInternalTemp) = []

function init(su::ArduinoSurveillanceUnitInternalTemp)
  @info "Initializing ArduinoSurveillanceUnit with ID $(su.deviceID)"
  sp = SerialPort(su.params.portAdress)
  open(sp)
	set_speed(sp, su.params.baudrate)
	set_frame(sp, ndatabits=su.params.ndatabits, parity=su.params.parity, nstopbits=su.params.nstopbits)
	# set_flow_control(sp,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
  sleep(2)
  flush(sp)
  write(sp, "!ConnectionEstablished*#")
  response = readuntil(sp, Vector{Char}(su.params.delim), su.params.timeout_ms);
  @info response
  if (response == "ArduinoSurveillanceV1" || response == "ArduinoSurveillanceV2"  )
    @info "Connection to ArduinoSurveillanceUnit established"
    su.sd = SerialDevice(sp, su.params.pause_ms, su.params.timeout_ms, su.params.delim, su.params.delim)
    su.ard = SimpleArduino(;commandStart = su.params.commandStart, commandEnd = su.params.commandEnd, delim = su.params.delim, sd = sd)
  else
    throw(ScannerConfigurationError(string("Connected to wrong Device", response)))
  end
  su.present = true
end

sendCommand(su::ArduinoSurveillanceUnitInternalTemp, cmdString::String) = sendCommand(su.ard, cmdString)

function getTemperatures(Arduino::ArduinoSurveillanceUnitInternalTemp) # toDo: deprecated
  Temps = sendCommand(Arduino, "GET:TEMP")
  TempDelim = "T"

  temp =  tryparse.(Float64, split(Temps, TempDelim))

  # We hardcode this to four here
  tempFloat = zeros(4)

  for i = 1:min(length(tempFloat), length(temp))
    if temp[i] !== nothing
        tempFloat[i] = temp[i]
    end
  end

  return tempFloat
end


text:
export ReaderCountRWLock, read_lock, read_unlock, is_read_locked

mutable struct ReaderCountRWLock <: Base.AbstractLock
    m::Threads.ReentrantLock
    reader_count::Int
    ReaderCountRWLock() = new(Threads.ReentrantLock(), 0)
end

function read_lock(l::ReaderCountRWLock)
    lock(l.m) do
        l.reader_count += 1
    end
end

function read_lock(f, l::ReaderCountRWLock)
    read_lock(l)
    try
        return f()
    finally
        read_unlock(l)
    end
end

function read_unlock(l::ReaderCountRWLock)
    lock(l.m) do
        l.reader_count -= 1
        if l.reader_count < 0
            error("reader count negative")
        end
    end
end

function Base.lock(l::ReaderCountRWLock)
    while true
        lock(l.m)
        if l.reader_count > 0
            unlock(l.m)
        else
            break
        end
    end
end

function Base.unlock(l::ReaderCountRWLock)
    unlock(l.m)
end

function Base.islocked(l::ReaderCountRWLock)
    islocked(l.m)
end

function is_read_locked(l::ReaderCountRWLock)
    l.reader_count > 0
end

ext: jl


text:
@defstruct Index2OnehotLayer Layer (
  name :: AbstractString = "index2onehot",
  (dim :: Int = -2, dim != 0),
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops)),
  (n_class :: Int = 0, n_class > 0)
)

type Index2OnehotLayerState <: LayerState
  layer :: Index2OnehotLayer
  blobs :: Vector{Blob}

  dims  :: Vector{Int}
end

function setup(backend::Backend, layer::Index2OnehotLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type  = eltype(inputs[1])
  dims = Array(Int, length(inputs))
  blobs = Array(Blob, length(inputs))
  for i = 1:length(inputs)
    total_dim = ndims(inputs[i])
    dim = layer.dim < 0 ? layer.dim + total_dim + 1 : layer.dim
    @assert 1 <= dim <= total_dim
    @assert dim != total_dim # should not operate on the mini-batch dimension

    dims[i] = dim
    blob_size = collect(size(inputs[i]))
    @assert blob_size[dim] == 1
    blob_size[dim] = layer.n_class
    blobs[i] = make_blob(backend, data_type, tuple(blob_size...))
  end
  return Index2OnehotLayerState(layer, blobs, dims)
end

function shutdown(backend::Backend, state::Index2OnehotLayerState)
  map(destroy, state.blobs)
end

function forward(backend::CPUBackend, state::Index2OnehotLayerState, inputs::Vector{Blob})
  for ii = 1:length(inputs)
    erase!(state.blobs[ii])

    input = inputs[ii].data
    output = state.blobs[ii].data
    op_dim = state.dims[ii]

    dim_pre, dim_prob, dim_post = split_dims(input, op_dim)
    @assert dim_prob == 1
    for i = 0:dim_pre-1
      for j = 0:dim_post-1
        @inbounds label = convert(Int, input[i + dim_pre*j + 1])
        idx = i + dim_pre*(label + state.layer.n_class*j) + 1
        @inbounds output[idx] = 1
      end
    end
  end
end

function backward(backend::Backend, state::Index2OnehotLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end


text:
"""
    visualize_sampling(input::Input, output)

Draw a schematic representation of the setup.
"""
function visualize_sampling(input::Input, output)

    Drawing(1500, 1500, "basic_test.png")
    origin()
    background("white")
    setopacity(1)

    setline(8)
    sethue("black")
    polysmooth(box(O, 1000, 720, vertices=true), 10, :stroke)
    fontsize(120)
    fontface("italic")
    textcentred("U")

    setline(5)
    inter = 700/(input.m-1)
    counter = 1

    for i in 10:inter:710
        line(Point(-550,i-360), Point(-500,i-360), :stroke)
        line(Point(500,i-360), Point(550,i-360), :stroke)

        if input.r.state[counter] == 1
            p = Point(-550,i-360)
            sethue("blue")
            circle(p,15,:fill)
            sethue("black")
        end

        counter += 1
    end

    sethue("red")
    fontsize(25)
    fontface("italic")
    counter = countmap(output)
    for i = 1:length(output)
        p = Point(550, output[i]*inter-350-inter)
        circle(p,15,:fill)
        if counter[output[i]] > 1
            sethue("black")
            Luxor.text(string(counter[output[i]]), p, valign=:center, halign=:center)
            sethue("red")
        end
    end

    clipreset()
    finish()
    preview()
end

"""
    visualize_proba(input::Input, output, data)

Draw a schematic representation associated to the event with probability `data`.
"""
function visualize_proba(input::Input, output, data)

    nlist = output_mode_occupation(input.n, input.m)

    if output in nlist
        idx = findfirst(x -> x==output, nlist)
        print("event probability: ", data[idx])
        visualize_sampling(input, output)
    else
        throw(ArgumentError("invalid argument(s)"))
    end

end


text:
# File: fVAI.jl
# Author:
# Date: 29-June-2020

# VALUE ADDED INDEX (VAI)
function fVAI(x)

# --------------------------------------------------------------------------
# x:       returns vector
# vaidx:   Value Added Index (VAI)
# --------------------------------------------------------------------------

    n = size(x,1); # data points
    fVAI = ones(1,n); # pre-allocate memory()
    fVAI(1) = 100; # initialise to $100
    for i = 1:i:end
        fVAI(i+1) = fVAI(i) * (1+x[i]./100); # calculate VAI
    end
    vaidx = fVAI;
    return vaidx
end


text:
using CovarianceRealism
using Test

@testset "CovarianceRealism.jl" begin
    # Write your own tests here.
end


text:
## Things outside of the module won't get run when using the package.

#You can create vectors of random normals.
x = randn(40);
x'
# Adding the semicolon to the end hides the output.
# Vectors are column vectors by default, which take up a good chunk of screen.
# So often it's nicer to tranpose them.

#You index into arrays with brackets, like in R
x[3]
x[end]
x[2:5]
x[end-4:end-1]
x[36:end-1]

#You also index into strings that way.
greeting = "Hello world!"
greeting[3:4]
greeting[7:end-1]

#Or arrays of random uniforms.
u = rand(3,2,2)
u[:,:,1]

#Standard syntax for randomly sampling is to call rand on the thing you want to sample from.
rand(1:4, 10)

using Distributions
gamma_3_3 = Gamma(3, 3)
rand(gamma_3_3, 2, 4)

#Julia Arrays are like the arrays/vectors/matrices of R. They will also automatically promote
[3f0 4.0 2 π]

#The individual elements were
3f0 |> typeof
4.0 |> typeof
2 |> typeof
π |> typeof
#Note that pi is an irrational number! It automatically converts to whatever you want it to be.
2π
2π |> typeof
big(2)π
big(2)π |> typeof
# |> is the pipe operator in Julia.
# However, it isn't as powerful or fancy as magrittr's "%>%".
# You can look at Lazy.jl if you want more piping options.

#In "R" a scalar is really a length 1 vector.
#In julia, they are clearly distinct.
scalar4double = 4.0
zerodim4double = fill(4.0)
vector4double = fill(4.0, 1)
matrix4double = fill(4.0, 1, 1)
dim3array4double = fill(4.0, 1, 1, 1)

#They all print differently, but there's also a fundamental difference under the hood.
isbits(scalar4double)
isbits(zerodim4double)
isbits(vector4double)

#Q: Are only numbers bits?
isbits(gamma_3_3)

#"isbits" is low level, referring to how the memory is handled under the hood.
#"isbits" types are raw, they're handled directly, and an array of them will be one solid chunk of memory.
# This makes them
#If something is not bits, that normally means it
#Why would we ever want a zero dim array?
function f!(a::Number)
    a = 2.1
end
function f!(a::AbstractArray)
    a[:] = 2.1
end
f!(scalar4double)
scalar4double
f!(zerodim4double)
zerodim4double
f!(dim3array4double)
dim3array4double

# In Julia, variables are not copied.
# The first f! replaced the input number with another number, but the original scalar4double remained bound to 4
# The second f! modified the contents of the arrays
# The arrays stayed the same, but the contents changed.
# Be weary!
zerodim2 = zerodim4double
zerodim2[] = 17.0
zerodim2
zerodim4double
# Julia does not copy unless you ask it to!
# "=" binds a new name to the old variable.
zerodim3 = copy(zerodim2)
zerodim3[] = 9.4
zerodim2
zerodim3

# It is convention in Julia that whenever a function has side effects -- whenever it modifies inputs
# to end that function name with a "bang"!
# The mutated arguments are also placed first.
copy!(zerodim2, zerodim3)
zerodim2
fill!(zerodim2, -1e9)
zerodim4double

# Also, adding commas makes a column vector.
[3f0, 4.0, 2, π]
# For a matrix
[ 1 2 3 ; 4 5 6]
reshape(1:6, 2, 3)
# Note the differences! Order, and
# In Julia, ranges are lazy -- even when reshaped.

# Also, the colon operator also has low priority.
# What is this in R?
# using RCall
3*1:5+7
10:3:22
10:3:22 |> collect

# R is a lazy language, which is great.
# Julia is eager, but many functions and operations are implemented to be lazy.
# Array slices can also be lazy
@views x[4:7]
# On Julia 0.7, transposing matrices is also lazy. Not yet on 0.6.

#If Julia can't promote them all to the same type, it'll give up and create a vector that can hold anything.
[ 42im 9//2 "fish" π ]
#This isn't a good idea. Don't that.

#You can also create undefined arrays. Be cautious, they may contain NaNs, which can infect code if you're not careful!
Array{Float64}(4, 12)

#Or fill
fill(1.2, 4, 7)
fill(gamma_3_3, 2, 4)

#Or if you just want zeros
zeros(4,7)
zeros(Int, 4, 7)

#You can also use comprehensions. Remember x earlier?
x'
#Lets use it to create a design matrix!
X = [xᵢ^j for xᵢ ∈ x, j ∈ 0:4]

#Lets create some fake data! First, lets get a true beta:
β_true = [2, 3, -5, -2, 1.2]

#And then add a little noise
y = X * β_true .+ 0.1randn(40)

β_hat = X \ y

# "\" is the left devision operator.
# It just means use the thing on the left to divide the thing on the right
2 \ 4
# If we tranpose both X and y, we can get dimensions to line up for a right division
y' / X'

# This is of course equivalent to what we're all used to:
inv(X' * X) * X' * y

# Realstically, we'd always prefer to just use LM and immediately get convenient summary statistics.
using DataFrames, GLM
df = DataFrame(x = x, x² = X[:,3], x³ = X[:,4], x⁴ = X[:,5], y = y)
OLS = glm(@formula(y ~ x + x² + x³ + x⁴), df, Normal(), IdentityLink())

# There are three ways to define functions.
# For an anonymous function, use [arguments] -> [function body]
# For example, if we want a random positive definite matrix
Σ = randn(20,5) |> x -> x' * x

# If we have more than one argument, we need to wrap them in a parenthesis.
anon = (a, b) -> sin(a) - exp(b)
anon(2, 3)

# We can overwrite anon with anything.
anon = 4

# If we create a proper named function, we cannot replace it with something other than a function.
f(a) = 2a
f(4)
f = 3
f(a) = 3a
f(4)
f(a::Int) = "Hello world!"
f(a, b) = 3a + 2b
f(4)
f(4.0)
f(4, -5)

# Julia uses multiple dispatch.
# Julia compiles a different specialized version of that function for each combination of types you give it.
# It then calls the type matching that combination.
# This also lets you define a different body for the types if you want to.
# That is, Julia compiles 4 separate functions for "f" when we do this:
f(4, -5)
f(4.0, -5)
f(4, -5.0)
f(4.0, -5.0)
# And chooses the correct version to use.
# If we want the compiled function to act differently, we can then give Julia a different function body to use for that type
# Just like we did for "Hellow world!" above
# For example:
f(X, β_hat) #Error, dimensions don't line up!
f(a::Matrix, b::Vector) = b' * a' * a * b
f(X, β_hat)
f(4, -5) #Still works!
# Similar to match.call() in R:
pair = (4,-5)
f(pair...)

# The long form just makes it easier to write multi-line functions.
#@inline function f(x, y)
function f(x, y)
    out = zero(promote_type(typeof(x), eltype(y)))
    for yᵢ ∈ y
        out += x * yᵢ
    end
    out
end

@code_llvm f(e, x)

#This is a tuple. A tuple's type is partly defined by its length.
@code_llvm f(e, (1.2, -9.5, 3.2, 37.2))
(1.2, -9.5, 3.2, 37.2) |> typeof

@code_llvm f(e, 2.3)
@code_llvm f(2, 3)

function g(f, a, b, c)
    f(a, b) + c
end

# f is automatically inlined at -O3, not at -O2.
g(f, 2, 3, 4)
@code_llvm g(f, 2, 3, 4)

# In Julia, you end a scope or expression via "end".
f(√2) #Does not have the irrational type! It's just taking the square root of 2.

# R has "single dispatch", which is useful for things like plot, print, and summary methods.

# A simple example why this is useful in Julia is so that we can use "*" for matrix multiplication.
# There are a lot of different methods for "*"!
methods(*) |> length
#Run just methods(*) if you want to see the full list!

# This also makes code faster!
# Compare
S_big = randn(10_000, 5_000) |> x -> x' * x;
X_big = randn(5_000, 5_000);
@time S_big * X_big;
using Compat.LinearAlgebra #This line is currently unnecessary
@time eigfact(S_big);
# However, we know S is symmetric! What happens if we tell Julia that?
S_big2 = Symmetric(S);
@time S_big2 * X_big;
@time eigfact(S_big2);

#The using Compat.LinearAlgebra line is unnecassary in Julia 0.6
#However, starting in Julia 0.7 you'll need "using LinearAlgebra"
#Saying "using Compat.LinearAlgebra" works on both versions.

# Give Julia more information, and it will automatically pick the best method!

# Multiple dispatch is magical, because
# a) It makes code easier to read. When you want to multiply, you use "*".
# When you want an eigenfactorization, you use "eigfact".
# No "eigfact" here, "eigfact_symmetric" there.
# No map_int here and map_dbl there.

# However, what if we actually want elementise multiplication from "*"?
# In Julia, adding a "." automatically broadcosts.
# A broadcast is like an apply or map statement, except it will line up dimensions.

# Recall
size(X)
size(y)
y' * X # Regular matrix-vector product.
X .* y # Broadcast
# Both X and y have 40 rows, so it lines these up
# and then replicates y across each of X's 5 columns.
X .* β_hat'
# β_hat' is 1 x 5, so it lines up with X's 5 columns, and is replicated down the 40 rows.
X .+ 2
# Here, 2 gets replcicated along both colums and rows.

X .* y' # Does not line up -> error
X .* β_hat # Does not line up -> error

# This works for any function.
f.(X)
f.(X, y)

# If we use multiple ones, they will all fuse.
f.(X, y) .* β_hat' .+ 2

# Too many dots can be annoying to write. So Julia offers the macro "@. "
@. f(X, y) * β_hat + 2

# Whenever you see something with an "@" in Julia, it is a macro.
# Macros edit your code. Think of them as non-standard evaluation in R.
# Because they can do extremely weird things, they always have an "@" to let you know something funky is going on.
# A basic example of non-standard evaluation in R is how a lot of functions know the name of the variables you give them.
# Consider -- straight out of Advanced R:
"
x <- seq(0, 2 * pi, length = 100)
sinx <- sin(x)
plot(x, sinx, type = "l")
"
# Similarly:
@show anon;
# If you ever want to know what it is a macro does to your code, there's a macro for that!
@macroexpand @show anon
@macroexpand @. f(X, y) * β_hat + 2

# We can see that "@." does us the favor of adding "." to every function call for us.

# We can get creative!
bunch_o_gammas = @. Gamma(exp(f(X)), exp(-2y))

# What is the meaning of this? Good question!
mean.(bunch_o_gammas)
# Maybe we want to do parametric bootstrap?
samples = rand.(bunch_o_gammas, 4000)

# Method of Moment's estimator for Gamma
function gammas_MoM(x)
    x̄ = mean(x)
    s² = var(x)
    x̄^2 / s², s² / x̄
end
mom_estimators = gammas_MoM.(samples)
function relative_error( estimate, truth )
    @. ( estimate - truth ) / truth
end
#Relative errors
@. relative_error(mom_estimators, params(bunch_o_gammas))

#We can easily condense the steps
N = 4000
@. relative_error(gammas_MoM(rand(bunch_o_gammas, N)), params(bunch_o_gammas))

### Some more array manipulations and fun with broadcasting!
# Remember Σ from earlier?
Σ
# Let's take N samples from a multivariate normal with that as the covariance matrix, and mean:
μ = 2.* randn(5)

W = chol(Σ) * randn(5, N) .+ μ

mean(W, 2)
cov(W')
# We can also recall the definition from multivariate
S = (W .- mean(W,2)) |> x -> x * x' ./ (size(x,2)-1)
#and see that it is a pretty good estimator
Σ
# The matrix and broadcast syntax makes it easy to weave together array manipulations
# and apply definitions, like we see in Johnson and Wichern for the sample Covariance matrix.

# Remember this?
[ 42im 9//2 "fish" π ]
# What if we do want to hold a bunch of different types?
# One straightforward approach is to build our own container!

abstract type four_objects{A,B,C,D} end
struct FourThings{A,B,C,D} <: four_objects{A,B,C,D}
    a::A
    b::B
    c::C
    d::D
end
function FourThings(a::A, b::B, c::C, d::D) where {A,B,C,D}
    FourThings{A,B,C,D}(a,b,c,d)
end
ft1 = FourThings(42im, 9//2, "fish", π)
isbits(ft1)
ft2 = FourThings(42, 9//2,  4f0, π)
isbits(ft2)
mutable struct FourMutableThings{A,B,C,D} <: four_objects{A,B,C,D}
    a::A
    b::B
    c::C
    d::D
end
ft3 = FourMutableThings(42, 9//2,  4f0, π)
isbits(ft3)
struct MyMatrix{T}
    a::T
    b::T
    c::T
    d::T
    e::T
    f::T
    g::T
    h::T
    i::T
    j::T
    k::T
    l::T
    m::T
    n::T
    o::T
    p::T
end
# Then, we get define our operations. For example, multiplication
function Base.:*(x::MyMatrix{T}, y::four) where {T, A,B,C,D,four <: four_objects{A,B,C,D}}
    (
        x.a*y.a + x.e * y.b + x.i * y.c + x.m * y.d,
        x.b*y.a + x.f * y.b + x.j * y.c + x.n * y.d,
        x.c*y.a + x.g * y.b + x.k * y.c + x.o * y.d,
        x.d*y.a + x.h * y.b + x.l * y.c + x.p * y.d,
    )
end
v16 = randn(Float32, 16);
m16 = reshape(v16, 4, 4);
v4 = randn(Float32, 4);
mm = MyMatrix(v16...)
mm * ft2
mm * ft3
ft4 = FourThings(v4...)
ft5 = FourMutableThings(v4...)
m16 * v4
mm * ft4

# What do you expect their performance to be?
# Built in multiplication, vs our own custom version?
using BenchmarkTools
@benchmark $m16 * $v4
@benchmark $mm * $ft2
@benchmark $mm * $ft3
@benchmark $mm * $ft4
@benchmark $mm * $ft5
# Now compare with Julia launched with -O3
# This is why I like -O3

# More cool features
using ForwardDiff

ForwardDiff.derivative(f, 2.3)
ForwardDiff.derivative(sin, 2.3)
cos(2.3)


text:
import NonUniformRandomVariateGeneration.sampleCategorical

@inline function sampleOne(xs::Vector{Float64}, sumXs::Float64, rpio::RPIO)
  M::Int64 = length(xs)
  if rpio.threads > 1 && M > PARALLEL_COUNTS_THRESHOLD
    # rpio.scratch contains the summed weights for each block
    idx::Int64 = sampleOne(rpio.scratch, sumXs)
    nt::Int64 = rpio.threads
    b::Int64 = div(M, nt)
    start::Int64 = (idx-1)*b + 1
    @inbounds return sampleOne(xs, rpio.scratch[idx], start)
  else
    return sampleCategorical(xs, sumXs)
  end
end

@inline function rpEstimate(model::Model, rpio::RPIO{Particle}) where Particle
  potentials::Vector{Float64} = rpio.potentials
  ks::Vector{Int64} = rpio.ks
  v::Float64 = 0.0
  _sampleParticles(model, rpio)
  b::Int64 = div(rpio.N, rpio.threads)
  for p in 1:rpio.n
    s::Float64 = _evaluatePotentials(model, rpio, p)
    for j in 1:p-1
      @inbounds block::Int64 = min(ceil(ks[j]/b), rpio.threads)
      @inbounds rpio.scratch[block] -= potentials[ks[j]]
      @inbounds s -= potentials[ks[j]]
      @inbounds potentials[ks[j]] = 0.0
    end
    v += log(s / (rpio.N-p+1))
    if p < model.n
      @inbounds ks[p] = sampleOne(potentials, s, rpio)
    end
  end
  return v
end

@inline function rpEstimateDiscrete(model::Model, rpio::RPIO{Int64Particle})
  _sampleParticles(model, rpio)
  ζcounts::Vector{Int64} = _getCounts(rpio)
  potentials::Vector{Float64} = Vector{Float64}(undef, length(ζcounts))

  v::Float64 = 0.0
  for p in 1:rpio.n
    s::Float64 = _countPotentials(potentials, ζcounts, model, rpio, p)
    v += log(s / (rpio.N-p+1))
    if p < model.n
      k::Int64 = sampleOne(potentials, s, rpio)
      @inbounds ζcounts[k] -= 1
    end
  end
  return v
end


text:
using Test
using FixedEffects
using StatsBase
using PooledArrays, CategoricalArrays
import Base: ==

==(x::FixedEffect{R,I}, y::FixedEffect{R,I}) where {R,I} =
    x.refs == y.refs && x.interaction == y.interaction && x.n == y.n

@testset "FixedEffect" begin
    fe1 = FixedEffect(1:10)
    @test sprint(show, fe1) == "Fixed Effects"
    if VERSION <= v"1.5"
        @test sprint(show, MIME("text/plain"), fe1) == """
            Fixed Effects:
              refs (10-element Array{Int64,1}):
                [1, 2, 3, 4, 5, ... ]
              interaction (UnitWeights):
                none"""
    else
        @test sprint(show, MIME("text/plain"), fe1) == """
            Fixed Effects:
              refs (10-element Vector{Int64}):
                [1, 2, 3, 4, 5, ... ]
              interaction (UnitWeights):
                none"""
    end
    fe2 = FixedEffect(1:10, interaction=fill(1.23456789, 10))
    if VERSION <= v"1.5"
        @test sprint(show, MIME("text/plain"), fe2) == """
            Fixed Effects:
              refs (10-element Array{Int64,1}):
                [1, 2, 3, 4, 5, ... ]
              interaction (10-element Array{Float64,1}):
                [1.23457, 1.23457, 1.23457, 1.23457, 1.23457, ... ]"""
    else
        @test sprint(show, MIME("text/plain"), fe2) == """
        Fixed Effects:
          refs (10-element Vector{Int64}):
            [1, 2, 3, 4, 5, ... ]
          interaction (10-element Vector{Float64}):
            [1.23457, 1.23457, 1.23457, 1.23457, 1.23457, ... ]"""
    end

    @test_throws DimensionMismatch FixedEffect(1:10, interaction=fill(1, 5))

    @test size(fe1) == (10,)
    @test length(fe1) == 10
    @test eltype(fe1) == Int
    @test eltype(fe2) == Float64

    @test fe1[:] === fe1
    subfe1 = fe1[[1,2]]
    @test subfe1.refs == fe1.refs[1:2]
    @test subfe1.interaction == uweights(2)
    @test subfe1 == fe1[1:2]
    @test subfe1 == fe1[fe1.refs.<=2]
    @test_throws MethodError fe1[[1 2]]

    subfe2 = fe2[[1,2]]
    @test subfe2.refs == fe2.refs[1:2]
    @test subfe2.interaction == fe2.interaction[1:2]
    @test subfe2 == fe2[1:2]
    @test subfe2 == fe2[fe2.refs.<=2]
    @test_throws MethodError fe2[[1 2]]
end

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects


text:
include("sand_05.jl")
using .Sand_05
using Plots
fig = Plots.plot()

NN = 160 # no of cells
tmax = 6.0 # max time
dt = 2.0 # how often to plot
tp = 0.0#5000 # when to start plotting

r= Sand_05.ring(NN; U=-2.0)#5.5)#5.5)
for t in 0.0:dt:tmax
#for i in 1:3#10
    t = r.t
    if t ≥ tp
        x0_ = 0.5 * (r.X_[1:end-1] + r.X_[2:end])
        p1 = Plots.plot!(x0_[1:end],r.Σ_[1:end],label="t: $t",legend=false)#:bottomright)#,title="Σ")#;legend=false)
        #Plots.scatter!(x0_[1:end],r.Σ_[1:end],label=false)
        display(fig)
        p2 = Plots.plot!(x0_[1:end],20*r.J_[1:end],label=false,linestyle=:dash)#,label="t: $t",legend=:topleft)#,title="J")#;legend=false)
        #Plots.scatter!(x0_[1:end],10*r.J_[1:end],shape=:+,label=false)
        display(fig)
    end
    Sand_05.tstep!(r,dt)
    #Sand_05.step!(r)
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
#Get number of detections for each site
function ndetects(y::Array)
  function sum2(row::Array)
    sum(skipmissing(row))
  end
  mapslices(sum2, y, dims=2)
end

#Simulate random covariate frame based on provided formula
function gen_covs(f::FormulaTerm, n::Int)
  covs = varnames(f).rhs
  if isnothing(covs) return DataFrame(_dummy=ones(n)) end
  nc = length(covs)
  out = DataFrame(reshape(rand(Normal(0,1), n*nc), n, nc), :auto)
  DataFrames.rename!(out, covs)
  out
end

#Replicate missing values in a new simulated y matrix
function rep_missing!(ynew::Array, y::Array)
  na_idx = findall(ismissing.(y))
  if length(na_idx) == 0 return nothing end
  ynew[na_idx] = fill(missing, length(na_idx))
end

#Get appropriate value of K
function check_K(K::Int, y::Array{Union{Missing,Int},2})
  if(K < (maximum(skipmissing(y))))
    error("K should be larger than the maximum observed abundance")
  end
  K
end

function check_K(K::Nothing, y::Array{Union{Missing,Int},2})
  maximum(skipmissing(y)) + 20
end

#Get minimum bound for integrating over possible abundance values K
#Minimum bound = maximum abundance ever observed at a site
function get_Kmin(y::Array{Union{Missing,Int64},2})
  N = size(y, 1)
  kmin = zeros(Int64, N)
  for n in 1:N
    if(all(ismissing.(y[n, :]))) continue end
    kmin[n] = maximum(skipmissing(y[n, :]))
  end
  kmin
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
mutable struct DynamicGenerator{
    M <: Machine,
    S <: Shaft,
    A <: AVR,
    TG <: TurbineGov,
    P <: PSS,
} <: DynamicInjection
    number::Int64
    name::String
    bus::Bus
    ω_ref::Float64
    V_ref::Float64
    P_ref::Float64
    Q_ref::Float64
    machine::M
    shaft::S
    avr::A
    tg::TG
    pss::P
    n_states::Int64
    states::Vector{Symbol}
    ext::Dict{String, Any}
    function DynamicGenerator(
        number::Int64,
        name::String,
        bus::Bus,
        ω_ref::Float64,
        V_ref::Float64,
        P_ref::Float64,
        Q_ref::Float64,
        machine::M,
        shaft::S,
        avr::A,
        tg::TG,
        pss::P,
    ) where {M <: Machine, S <: Shaft, A <: AVR, TG <: TurbineGov, P <: PSS}

        n_states =
            (machine.n_states + shaft.n_states + avr.n_states + tg.n_states + pss.n_states)

        states = vcat(machine.states, shaft.states, avr.states, tg.states, pss.states)

        new{M, S, A, TG, P}(
            number,
            name,
            bus,
            ω_ref,
            V_ref,
            P_ref,
            Q_ref,
            machine,
            shaft,
            avr,
            tg,
            pss,
            n_states,
            states,
            Dict{String, Any}(),
        )
    end
end

get_Sbase(device::DynamicGenerator) = device.machine.s_rated
get_bus(device::DynamicGenerator) = device.bus
get_states(device::DynamicGenerator) = device.states
get_n_states(device::DynamicGenerator) = device.n_states
get_number(device::DynamicGenerator) = device.number
get_name(device::DynamicGenerator) = device.name
get_ω_ref(device::DynamicGenerator) = device.ω_ref
get_V_ref(device::DynamicGenerator) = device.V_ref
get_P_ref(device::DynamicGenerator) = device.P_ref
get_Q_ref(device::DynamicGenerator) = device.Q_ref
get_machine(device::DynamicGenerator) = device.machine
get_shaft(device::DynamicGenerator) = device.shaft
get_avr(device::DynamicGenerator) = device.avr
get_tg(device::DynamicGenerator) = device.tg
get_pss(device::DynamicGenerator) = device.pss
get_ext(device::DynamicGenerator) = device.ext


text:
# Belief Propagation                #
# --------------------------------- #
@rule MvNormalMeanPrecision(:μ, Marginalisation) (m_out::PointMass, m_Λ::PointMass) = MvNormalMeanPrecision(mean(m_out), mean(m_Λ))

@rule MvNormalMeanPrecision(:μ, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_Λ::PointMass) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    return MvNormalMeanCovariance(m_out_mean, m_out_cov + cholinv(mean(m_Λ)))
end

# Variational                       #
# --------------------------------- #
@rule MvNormalMeanPrecision(:μ, Marginalisation) (q_out::Any, q_Λ::Any) = MvNormalMeanPrecision(mean(q_out), mean(q_Λ))

@rule MvNormalMeanPrecision(:μ, Marginalisation) (m_out::PointMass, q_Λ::Any) = MvNormalMeanPrecision(mean(m_out), mean(q_Λ))

@rule MvNormalMeanPrecision(:μ, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, q_Λ::Any) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    return MvNormalMeanCovariance(m_out_mean, m_out_cov + cholinv(mean(q_Λ)))
end


text:
struct ObjectData
    objects::Vector{AbstractWorldObj}
    collidable_centers::Vector{Vector{Float32}}
    collidable_corners::Vector{Array{Float32, 2}}
    collidable_norms::Vector{Array{Float32, 2}}
    collidable_safety_radii::Vector{Float32}
end

function ObjectData(map_data::Dict, road_tile_size::Float32, domain_rand::Bool, grid::Vector{Union{Missing, Dict{String,Any}}},
                    grid_width::Int, grid_height::Int)
    # Create the objects array
    objects = Vector{AbstractWorldObj}()

    # The corners for every object, regardless if collidable or not
    #object_corners = []

    # Arrays for checking collisions with N static objects
    # (Dynamic objects done separately)
    # (N x 2): Object position used in calculating reward
    collidable_centers = Vector{Vector{Float32}}()

    # (N x 2 x 4): 4 corners - (x, z) - for object's boundbox
    collidable_corners = Vector{Matrix{Float32}}()

    # (N x 2 x 2): two 2D norms for each object (1 per axis of boundbox)
    collidable_norms = Vector{Matrix{Float32}}()

    # (N): Safety radius for object used in calculating reward
    collidable_safety_radii = Vector{Float32}()

    # For each object
    for (obj_idx, desc) in enumerate(get(map_data, "objects", []))
        kind = desc["kind"]

        pos = Float32.(desc["pos"])
        x, z = pos[1:2]
        y = length(pos) == 3 ? pos[3] : 0f0

        rotate = Float32(desc["rotate"])
        optional = get(desc, "optional", false)

        pos = road_tile_size * [x, y, z]

        # Load the mesh
        #TODO
        mesh = ObjMesh.get(kind)

        if "height" in keys(desc)
            scale = Float32(desc["height"]) / mesh.max_coords[2]
        else
            scale = Float32(desc["scale"])
        end

        @assert !("height" ∈ keys(desc) && "scale" ∈ keys(desc)) "cannot specify both height and scale"

        static = get(desc, "static", true)

        obj_desc = Dict([
            "kind"=> kind,
            "mesh"=> mesh,
            "pos"=> pos,
            "scale"=> scale,
            "y_rot"=> rotate,
            "optional"=> optional,
            "static"=> static
        ])

        # obj = nothing

        if static
            if kind == "trafficlight"
                obj = TrafficLightObj(obj_desc, domain_rand, SAFETY_RAD_MULT)
            else
                obj = WorldObj(obj_desc, domain_rand, SAFETY_RAD_MULT)
            end
        else
            if kind == "duckiebot"
                obj = DuckiebotObj(obj_desc, domain_rand, SAFETY_RAD_MULT, WHEEL_DIST,
                                   ROBOT_WIDTH, ROBOT_LENGTH)
            elseif kind == "duckie"
                obj = DuckieObj(obj_desc, domain_rand, SAFETY_RAD_MULT, road_tile_size)
            else
                msg = "I do not know what object this is: $kind"
                error(msg)
            end
        end

        objects = vcat(objects, [obj])

        # Compute collision detection information

        # angle = rotate * (math.pi / 180)

	# Find drivable tiles object could intersect with
        possible_tiles = find_candidate_tiles(_obj_corners(obj), road_tile_size)

        # If the object intersects with a drivable tile
        if static && kind != "trafficlight" && _collidable_object(
                grid, grid_width, grid_height, _obj_corners(obj), _obj_norm(obj), possible_tiles, road_tile_size)
            collidable_centers = vcat(collidable_centers, [pos])
            collidable_corners = vcat(collidable_corners, [permutedims(_obj_corners(obj))])
            collidable_norms = vcat(collidable_norms, [_obj_norm(obj)])
            collidable_safety_radii = vcat(collidable_safety_radii, [_safety_radius(obj)])
        end
    end

    ObjectData(objects, collidable_centers, collidable_corners,
               collidable_norms, collidable_safety_radii)
end

function _get_curve(grid::Vector{Union{Missing,Dict{String,Any}}}, i::Int, j::Int,
                    width::Int, height::Int, road_tile_size::Float32)
    ##
    #    Get the Bezier curve control points for a given tile
    ##
    tile = _get_tile(grid, i, j, width, height)
    @assert !isnothing(tile)

    kind = tile["kind"]
    angle = tile["angle"]

    # Each tile will have a unique set of control points,
    # Corresponding to each of its possible turns

    if startswith(kind, "straight")
        pts = [
                [-0.2f0 0f0 -0.50f0;
                 -0.2f0 0f0 -0.25f0;
                 -0.2f0 0f0  0.25f0;
                 -0.2f0 0f0  0.50f0],

                [0.2f0 0f0  0.50f0;
                 0.2f0 0f0  0.25f0;
                 0.2f0 0f0 -0.25f0;
                 0.2f0 0f0 -0.50f0],
              ] .* road_tile_size

    elseif kind == "curve_left"
        pts = [
                [-0.2f0 0f0 -0.5f0;
                 -0.2f0 0f0  0.0f0;
                  0.0f0 0f0  0.2f0;
                  0.5f0 0f0  0.2f0],

                [0.5f0 0f0 -0.2f0;
                 0.3f0 0f0 -0.2f0;
                 0.2f0 0f0 -0.3f0;
                 0.2f0 0f0 -0.5f0],
              ] .* road_tile_size

    elseif kind == "curve_right"
        pts = [
                [-0.2f0 0f0 -0.5f0;
                 -0.2f0 0f0 -0.2f0;
                 -0.3f0 0f0 -0.2f0;
                 -0.5f0 0f0 -0.2f0],

                [-0.5f0 0f0  0.2f0;
                 -0.3f0 0f0  0.2f0;
                  0.3f0 0f0  0.0f0;
                  0.2f0 0f0 -0.5f0],
              ] .* road_tile_size

    # Hardcoded all curves for 3way intersection
    elseif startswith(kind, "3way")
        pts = [
                [-0.2f0 0f0 -0.50f0;
                 -0.2f0 0f0 -0.25f0;
                 -0.2f0 0f0  0.25f0;
                 -0.2f0 0f0  0.50f0],

                [-0.2f0 0f0 -0.5f0;
                 -0.2f0 0f0  0.0f0;
                  0.0f0 0f0  0.2f0;
                  0.5f0 0f0  0.2f0],

                [0.2f0 0f0  0.50f0;
                 0.2f0 0f0  0.25f0;
                 0.2f0 0f0 -0.25f0;
                 0.2f0 0f0 -0.50f0],

                [0.5f0 0f0 -0.2f0;
                 0.3f0 0f0 -0.2f0;
                 0.2f0 0f0 -0.2f0;
                 0.2f0 0f0 -0.5f0],

                [0.2f0 0f0 0.5f0;
                 0.2f0 0f0 0.2f0;
                 0.3f0 0f0 0.2f0;
                 0.5f0 0f0 0.2f0],

                [ 0.5f0 0f0 -0.2f0;
                  0.3f0 0f0 -0.2f0;
                 -0.2f0 0f0  0.0f0;
                 -0.2f0 0f0  0.5f0],
              ] .* road_tile_size

    # Template for each side of 4way intersection
    elseif startswith(kind, "4way")
        pts = [
                [-0.2f0 0f0 -0.5f0;
                 -0.2f0 0f0  0.0f0;
                  0.0f0 0f0  0.2f0;
                  0.5f0 0f0  0.2f0],

                [-0.2f0 0f0 -0.50f0;
                 -0.2f0 0f0 -0.25f0;
                 -0.2f0 0f0  0.25f0;
                 -0.2f0 0f0  0.50f0],

                [-0.2f0 0f0 -0.5f0;
                 -0.2f0 0f0 -0.2f0;
                 -0.3f0 0f0 -0.2f0;
                 -0.5f0 0f0 -0.2f0],
              ] .* road_tile_size
    else
        @assert false kind
    end

    # Rotate and align each curve with its place in global frame
    if startswith(kind, "4way")

        # Generate all four sides' curves,
        # with 3-points template above
	function size_gen(rot::Int)
	    mat = gen_rot_matrix([0f0, 1f0, 0f0], rot * Float32(π) / 2f0)
      	    pts_new = map(x -> x * mat, pts)
            add_vec = [(i-0.5f0) 0f0 (j-0.5f0);] * road_tile_size
       	    return hcat(map(x -> x .+ add_vec, pts_new)...)
	end

	fourway_pts = hcat(map(rot->size_gen(rot), 0:3)...)
        return reshape(fourway_pts, 4, 3, :)

    # Hardcoded each curve; just rotate and shift
    elseif startswith(kind, "3way")
        mat = gen_rot_matrix([0f0, 1f0, 0f0], angle * Float32(π) / 2f0)
        #NOTE: pts is 3D matrix, find a work around if * does not work
        pts_new = map(x -> x * mat, pts)
        add_vec = [(i-0.5f0) 0f0 (j-0.5f0);] * road_tile_size
        pts_new = map(x -> x .+ add_vec, pts_new)
        threeway_pts = hcat(pts_new...)

        return reshape(threeway_pts, 4, 3, :)
    else
        mat = gen_rot_matrix([0f0, 1f0, 0f0], angle * Float32(π) / 2f0)
        pts = map(x -> x * mat, pts)
        add_vec = [(i-0.5f0) 0f0 (j - 0.5f0);] * road_tile_size
        pts = map(x -> x .+ add_vec, pts)
    end

    pts = hcat(pts...)
    return reshape(pts, 4, 3, :)
end
#=
struct tile
    coords::Vector{Float32}
    texture
    color
    kind::String
    angle::Float32
    drivable::Bool
    curves
end
=#
struct Grid
    road_tile_size::Float32
    grid_width::Int
    grid_height::Int
    _grid::Vector
    obj_data::ObjectData
    road_vlist::Matrix{Float32}
    road_tlist::Matrix{Float32}
    ground_vlist::Matrix{Float32}
    drivable_tiles::Vector{Dict}
    mesh::ObjectMesh
    start_tile::Union{Nothing,Dict}
end

function _set_tile!(grid::Vector{Union{Missing,Dict{String,Any}}}, i::Int, j::Int,
                    width::Int, height::Int, tile::Dict{String,Any})
    @assert 1 ≤ i ≤ width
    @assert 1 ≤ j ≤ height
    grid[(j-1)*width + i] = tile
end

_get_tile(grid::Grid, i, j) = _get_tile(grid._grid, i, j, grid.grid_width, grid.grid_height)

function _get_tile(grid::Vector{Union{Missing,Dict{String,Any}}}, i::Int, j::Int, width::Int, height::Int)
    ##
    #Returns nothing if the duckiebot is not in a tile.
    ##
    if 1 ≤ i ≤ width && 1 ≤ j ≤ height
        return grid[(j-1)*width + i]
    end

    return nothing
end

function _init_vlists(road_tile_size)
    # Create the vertex list for our road quad
    # Note: the vertices are centered around the origin so we can easily
    # rotate the tiles about their center
    half_size = road_tile_size / 2f0
    verts = [
        -half_size 0f0 -half_size;
         half_size 0f0 -half_size;
         half_size 0f0  half_size;
        -half_size 0f0  half_size;
    ]
    texCoords = [
        1f0 0f0;
        0f0 0f0;
        0f0 1f0;
        1f0 1f0
    ]

    road_vlist = verts
    road_tlist = texCoords

    # Create the vertex list for the ground quad
    verts = [
        -1 -0.8f0  1;
        -1 -0.8f0 -1;
         1 -0.8f0 -1;
         1 -0.8f0  1
    ]

    ground_vlist = verts

    return road_vlist, road_tlist, ground_vlist
end

function Grid(map_data::Dict, domain_rand)
    if "tile_size" ∉ keys(map_data)
        msg = "Must now include explicit tile_size in the map data."
        throw(KeyError(msg))
    end

    road_tile_size = Float32(map_data["tile_size"])

    road_vlist, road_tlist, ground_vlist = _init_vlists(road_tile_size)

    tiles = map_data["tiles"]
    @assert length(tiles) > 0
    @assert length(tiles[1]) > 0

    # Create the grid
    grid_height = length(tiles)
    grid_width = length(tiles[1])
    _grid = Vector{Union{Missing, Dict{String,Any}}}(missing, grid_width * grid_height)

    # We keep a separate list of drivable tiles
    drivable_tiles = []

    # For each row in the grid
    for (j, row) ∈ enumerate(tiles)
        msg = "each row of tiles must have the same length"
        if length(row) != grid_width
            error(msg)
        end

        # For each tile in this row
        for (i, tile) in enumerate(row)
            tile = strip(tile)

            tile == "empty" && continue

            if '/' ∈ tile
                kind, orient = split(tile, '/')
                kind = strip(kind, ' ')
                orient = strip(orient, ' ')
                angle = findfirst(isequal(orient), ["S", "E", "N", "W"]) - 1
                drivable = true
            elseif '4' ∈ tile
                kind = "4way"
                angle = 2
                drivable = true
            else
                kind = tile
                angle = 0
                drivable = false
            end

            tile = Dict([
                "coords"=> [i, j],
                "kind"=> kind * "",
                "angle"=> angle,
                "drivable"=> drivable
            ])

            _set_tile!(_grid, i, j, grid_width, grid_height, tile)

            if drivable
                tile["curves"] = _get_curve(_grid, i, j, grid_width,
                                            grid_height, road_tile_size)
                push!(drivable_tiles, tile)
            end
        end
    end
    #TODO
    mesh = ObjMesh.get("duckiebot")
    obj_data = ObjectData(map_data, road_tile_size, domain_rand, _grid, grid_width, grid_height)

    # Get the starting tile from the map, if specified
    start_tile = nothing
    if "start_tile" ∈ keys(map_data)
        coords = map_data["start_tile"]
        coords .+= 1
        start_tile = _get_tile(_grid, coords..., grid_width, grid_height)
    end

    Grid(road_tile_size, grid_width, grid_height, _grid, obj_data, road_vlist,
         road_tlist, ground_vlist, drivable_tiles, mesh, start_tile)
end

_get_curve(grid::Grid, i, j) = _get_curve(grid._grid, i, j, grid.grid_width,
                                          grid.grid_height, grid.road_tile_size)

mutable struct Map
    map_name::String
    map_file_path::String
    map_data::Dict       # The parsed content of the map_file
    _grid::Grid
end

function Map(map_name::String, domain_rand::Bool)
    ##
    #Load the map layout from a YAML file
    ##
    map_file_path = get_file_path("src/maps", map_name, "yaml")
    #logger.debug('loading map file "%s"' % self.map_file_path)

    map_data = nothing
    open(map_file_path, "r") do f
        map_data = YAML.load(f)
    end

    _grid = Grid(map_data, domain_rand)

    Map(map_name, map_file_path, map_data, _grid)
end


text:
# using CPLEX
# using GLPK
using Cbc
using Clp
using JuMP

include("seqata_test_model_util.jl")
include("../src/solvers/mip_discret.jl")

# ===========
# Préparation des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:itermax, 300)
Args.set(:outdir, "/tmp") # equivalent à l'option -d/tmp

# ===========
# Construction de l'instance mini
inst = instance_build_mini10()
@test inst.nb_planes == 10

# ===========
lg1("1. Création d'un MipDiscretSolver pour Cbc... ")
Args.set("external_mip_solver", :cbc)
Args.set(:level, 0)
solver = MipDiscretSolver(inst)
Args.set(:level, TEST_LEVEL)

# Le nom du solver Cbc n'est pas standard :
#   @show solver_name(solver.model)
#   => "COIN Branch-and-Cut (Cbc)"
# Donc ceci ne fonctionne pas :
#   @test Symbol(lowercase(solver_name(solver.model))) == :cbc
# Je teste alors la présence de la sous-chaine cbs
@test occursin("cbc", lowercase(solver_name(solver.model)))

ln1(" fait.")

# ===========
lg1("2. Résolution par MipDiscretSolver pour Cbc (cost=700.0 ?) LONG (9mn+) ... ")
solve(solver)
bestsol = solver.bestsol
# @test bestsol.cost == 700.0
@test isapprox(bestsol.cost, 700.0) # approx car cbc trouve 699.999999999996
ln1(" (cost=$(solver.bestsol.cost)) fait.\n")


text:
## auc.jl (c) Area Under the ROC
## (c) 2016 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

"""`AUC(::roc; pfa=1.0, pmiss=1.0, normalize=true)` computes the Area Under the Curve
with a sense `auc → 0` indicating better performance.

Optional parameters `pfa` or `pmiss` limit integration over only part of the ROC curve.
`normalize` indicates comparing the partial ROC to the trivial ROC.
"""
function auc(r::Roc; pfa=1.0, pmiss=1.0, normalize=true)
    if pfa != 1.0 || pmiss != 1.0
        if pfa != 1.0 && pmiss != 1.0 || ! (0.0 < pfa ≤ 1.0) || !(0.0 < pmiss ≤ 1.0)
            throw(DomainError)
        end
        if pfa < 1.0
            i = binsearch(-pfa, -r.pfa) ## find threshold where FAR >= pfa
            i < length(r) || error("Unexpected ROC data")
            norm = normalize ? 1/pfa/(2-pfa) : 1
            a = -dot(r.pmiss[i:end-1] + r.pmiss[i+1:end], diff(r.pfa[i:end])) / 2
            ## compensate for the the fact that pfa > r.pfa[i]
            a -= (r.pfa[i] - pfa) * (r.pmiss[i] - (r.pfa[i] - pfa) * (r.pmiss[i+1] - r.pmiss[i]) / 2(r.pfa[i+1] - r.pfa[i]))
            return norm * a
        else
            i = binsearch(pmiss, r.pmiss)
            i < length(r) || error("Unexpected ROC data")
            norm = normalize ? 1/pmiss/(2-pmiss) : 1
            a = dot(r.pfa[1:i] + r.pfa[2:i+1], diff(r.pmiss[1:i+1])) / 2
            a -= (r.pmiss[i+1] - pmiss) * (r.pfa[i+1] - (r.pmiss[i+1] - pmiss) * (r.pfa[i+1] - r.pfa[i]) / 2(r.pmiss[i+1] - r.pmiss[i]))
            return norm * a
        end
    else
        ## The probability that a random non-target score is higher than a target score
        return -dot(r.pmiss[1:end-1] + r.pmiss[2:end], diff(r.pfa)) / 2
    end
end

auc(tar::Vector, non::Vector; kwargs...) = auc(roc(tar, non); kwargs...)

"""`AUC(::roc; pfa=1.0, pmiss=1.0, normalize=true)` computes the traditional Area Under the Curve
with a sense `AUC → 1` indicating better performance.

Optional parameters `pfa` or `pmiss` limit integration over only part of the ROC curve.
`normalize` indicates comparing the partial ROC to the trivial ROC.

You can also call this function as `AUC(targets::Vector, nontargets::Vector; kwargs...)`
"""
function AUC(args...; pfa=1.0, pmiss=1.0, normalize=true)
    if normalize
        return 1 - auc(args...; normalize=normalize, pfa=pfa, pmiss=pmiss)
    else
        return min(pfa, pmiss) - auc(args...; normalize=normalize, pfa=pfa, pmiss=pmiss)
    end
end


text:
#RSDeltaSigmaPort: Generate plots of the Signal-to-Noise Ratio
#-------------------------------------------------------------------------------

#==SNR plot
===============================================================================#
#TODO: Rename `plotSQNR`?
function plotSNR()
	plot = cons(:plot, linlin, title="SQNR vs. Input Level", legend=true,
		#xyaxes=set(xmin=-120, xmax=0, ymin=0, ymax=120),
		xyaxes=set(ymax=120), #Clip infinite gains only
		labels=set(xaxis="Input Level (dBFS)", yaxis="SQNR (dB)"),
	)
	return plot
end

function plotSNR!(plot, snrinfo, dsm::RealDSM; id::String="simulation", color=:green)
	simglyph = cons(:a, glyph=set(shape=:o, size=1, color=color, fillcolor=color))

	#Access data:
	(amp, SNR) = snrinfo.vs_amp_sim
	(pk_amp, pk_SNR) = snrinfo.peak

	#Convert to waveforms:
	SNR = waveform(amp, SNR)

	infostr = @sprintf("peak SNR = %4.1fdB\n@ A = %4.1f dBFS (OSR = %d)",
		pk_SNR, pk_amp, dsm.OSR #orig. coords: (-25, 85)
	)
	push!(plot,
		cons(:wfrm, SNR, simglyph, line=set(style=:dashdot, color=color, width=2), label="simulation"),
		cons(:atext, infostr, x=-25, y=85, align=:cr),
	)
	if !isnothing(snrinfo.vs_amp_predicted)
		(amp_pred, SNR_pred) = snrinfo.vs_amp_predicted
		SNR_pred = waveform(amp_pred, SNR_pred) #Convert to waveform
		push!(plot,
			cons(:wfrm, SNR_pred, line=set(style=:solid, color=:red, width=2), label="theory"),
		)
	end
	return plot
end

function plotSNR(snrinfo, dsm::RealDSM; id::String="simulation", color=:green)
	plot = plotSNR()
	plotSNR!(plot, snrinfo, dsm, id=id, color=color)
	return plot
end

#Last line


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_sastre_m6(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m6!(A_copy)
end

@inline function exp_sastre_m6!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 Bb3 B3 Ba4 Bb4 B4 Ba7_5 Bb5_3 B5 Ba6 Bb6 B6 Bb7 B7 T2k9
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing Bb3 = x*A+x*B2
    coeff1=0.001052151783051235
    coeff2=0.0004675683454147702
    memslots3 .= coeff1.*memslots1 .+ coeff2.*memslots2
    # Computing B3 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*A+x*B2+x*B3
    coeff1=0.2868706220817633
    coeff2=-0.03289442879547955
    coeff3=1.0
    memslots3 .= coeff1.*memslots1 .+ coeff2.*memslots2 .+ coeff3.*memslots4
    # Computing Bb4 = x*B2+x*B3
    coeff1=0.05317514832355802
    coeff2=1.0
    memslots5 .= coeff1.*memslots2 .+ coeff2.*memslots4
    # Computing B4 with operation: mult
    mul!(memslots6,memslots3,memslots5)
    # Deallocating Ba4 in slot 3
    # Deallocating Bb4 in slot 5
    # Computing Ba7_5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=3.09646797193604
    coeff2=0.772360321294401
    coeff3=0.1673139636901279
    coeff4=7.922322450524197
    coeff5=1.0
    # Smart lincomb recycle B3
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots4 .+ coeff5.*memslots6
    mul!(memslots4,true,I*coeff1,true,true)
    # Deallocating B4 in slot 6
    # Computing Bb5_3 = x*A+x*B2
    coeff1=0.002688394980266927
    coeff2=0.0004675683454147702
    memslots3 .= coeff1.*memslots1 .+ coeff2.*memslots2
    # Computing B5 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Bb5_3 in slot 3
    # Computing Ba6 = x*A+x*B2+x*B5
    coeff1=0.39689859154115
    coeff2=0.02219811707032801
    coeff3=1.0
    memslots3 .= coeff1.*memslots1 .+ coeff2.*memslots2 .+ coeff3.*memslots5
    # Computing Bb6 = x*B2+x*B5
    coeff1=0.0277140002806296
    coeff2=1.0
    memslots6 .= coeff1.*memslots2 .+ coeff2.*memslots5
    # Computing B6 with operation: mult
    mul!(memslots7,memslots3,memslots6)
    # Deallocating Ba6 in slot 3
    # Deallocating Bb6 in slot 6
    # Computing Bb7 = x*A+x*B2+x*B5+x*B6
    coeff1=0.3229486011362678
    coeff2=0.08092036376147299
    coeff3=1.930814505527068
    coeff4=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff1.*memslots1 .+ coeff2.*memslots2 .+ coeff3.*memslots5 .+ coeff4.*memslots7
    # Deallocating B2 in slot 2
    # Deallocating B5 in slot 5
    # Deallocating B6 in slot 7
    # Computing B7 with operation: mult
    mul!(memslots2,memslots4,memslots1)
    # Deallocating Ba7_5 in slot 4
    # Deallocating Bb7 in slot 1
    # Computing T2k9 = x*I+x*B7
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle B7
    memslots2 .= coeff2.*memslots2
    mul!(memslots2,true,I*coeff1,true,true)
    return memslots2 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using HTTP

"""
A helper function for converting an `input_source` to an IO if not already an IO object.
"""
function read_input_data(input_source::AbstractString; kwargs...)
    #TODO this should be done async
    #TODO(MD): we should drop this as soon as we have a better option for downloading from URL
    if startswith(input_source, "https://") || startswith(input_source, "http://")
        HTTP.get(input_source; kwargs...)
    else
        open(input_source, "r")
    end
end
read_input_data(input::IOStream; kwargs...) = input

"""
A helper function for memory-mapping an input if possible.
Otherwise, the input itself is returned.
"""
apply_mmap(input) = Mmap.mmap(input)
apply_mmap(input::HTTP.Messages.Response) = input.body

"""
A helper function that prepares an input_source by removing the leading BOM if available
and memory-mapping it if possible.
"""
function read_mmap_data(input_source::Union{AbstractString,IOStream}; kwargs...)
    input = read_input_data(input_source; kwargs...)
    input_io = IOBuffer(apply_mmap(input))
    consumeBOM!(input_io)
    input_io
end
function read_mmap_data(input::IOBuffer; kwargs...)
    consumeBOM!(input)
    input
end

function read_mmap_data(input::IO; kwargs...)
    buff_io = Parsers.BufferedIO(input)
    consumeBOM!(buff_io)
    buff_io
end


text:
"""
`spline_reparametrize`: Reparametrize a string of points with a specified
distance function using cubic spline interpolation
"""
function spline_reparametrize!(U, dist::TD, pin::Bool) where {TD}
    n_images = length(U);
    d = length(U[1]); #dimension of the each image
    s = zeros(n_images);

    # reference image values
    s_ref = LinRange(0, 1, n_images);

    # compute integrated arc length
    for i in 2:n_images
        s[i] = s[i-1] + dist(U[i], U[i-1]);
    end
    # normalize
    @. s = s/s[end];

    # interpolate back onto the uniform mesh
    for j in 1:d
        # Construct spline through the j-th image
        spl = Spline1D(s, [U[i][j] for i in 1:n_images], k=3);
        u_spl = spl(s_ref);
        # Interpolate the j-th image
        for i in 2:n_images-1
            U[i][j] = u_spl[i];
        end
        if !pin
            U[1][j] = u_spl[1];
            U[end][j] = u_spl[end];
        end
    end
    U
end


text:
function credal_marginal_upper_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    credal_marginal_upper_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    function credal_marginal_upper_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        #println("LITERAL")
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )

        return nothing
    end

    function credal_marginal_upper_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
       # println("AND NODE")

        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end

        return nothing
    end

    function credal_marginal_upper_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
       # println("OR NODE")

        n.pr .= 1e-300
        #pr(n) .= 1e-300

        l = Array{Float64}(undef, num_examples(data), length(n.children))
        u = Array{Float64}(undef, num_examples(data), length(n.children))
        c = Array{Float64}(undef, num_examples(data), length(n.children))

        for i=1:num_examples(data)

            l[i,:] .=  exp.(prob_origin(n).log_thetas)
            u[i,:] .=  exp.(prob_origin(n).log_thetas_u)
        end

        for i=1:num_examples(data)

            c[i,:] .= exp.(getindex.(pr.(n.children),i)) ## TODO log is monotone + broadcast
        end

        optx = Array{Float64}(undef, num_examples(data), length(n.children))
        optx = copy(u)

        for i=1:num_examples(data)

            for j=1:length(n.children)
                if c[i,j] == sort(c[i,:])[1]
                    optx[i,j] = 0.0
                    somma = sum(optx[i,:])
                    optx[i,j] = 1.0-somma ## TODO compatta questo pezzo con find o findall
                    break
                end
            end
        end

        for i=1:num_examples(data)
            pr(n)[i] = log(optx[i,:]'c[i,:])
        end

        return nothing

    end

    ## Pass Up on every node in order
    for n in circuit
        credal_marginal_upper_pass_up_node(n, cache, data)
    end

    return nothing
end

#########
###marginal lower

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

function credal_marginal_lower_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    credal_marginal_lower_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

   function credal_marginal_lower_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )
        return nothing
   end

    function credal_marginal_lower_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end
        return nothing
    end

    function credal_marginal_lower_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
        #println("OR NODE")

        n.pr .= 1e-300
        #pr(n) .= 1e-300

        l = Array{Float64}(undef, num_examples(data), length(n.children))
        c = Array{Float64}(undef, num_examples(data), length(n.children))

        for i=1:num_examples(data)

            l[i,:] .=  exp.(prob_origin(n).log_thetas)
        end

        for i=1:num_examples(data)

            c[i,:] .= exp.(getindex.(pr.(n.children),i)) ## TODO log is monotone + broadcast
        end

        optx = Array{Float64}(undef, num_examples(data), length(n.children))
        optx = copy(l)

        for i=1:num_examples(data)

            for j=1:length(n.children)
                if c[i,j] == sort(c[i,:])[1]
                    optx[i,j] = 0.0
                    somma = sum(optx[i,:])
                    optx[i,j] = 1.0-somma ## TODO compatta questo pezzo con find o findall
                    break
                end
            end
        end

        for i=1:num_examples(data)
            pr(n)[i] = log(optx[i,:]'c[i,:])
        end

        return nothing

    end

      ## Pass Up on every node in order
      for n in circuit
          credal_marginal_lower_pass_up_node(n, cache, data)
      end
      return nothing
end

##### marginal_pass_down ##### might need this for complete observation

# function marginal_pass_down(circuit::DownFlowΔ{O,F}) where {O,F}
#     resize_flows(circuit, flow_length(origin(circuit)))
#     for n in circuit
#         reset_downflow_in_progress(n)
#     end
#     for downflow in downflow_sinks(circuit[end])
#         # initialize root flows to 1
#         downflow.downflow .= one(eltype(F))
#     end
#     for n in Iterators.reverse(circuit)
#         marginal_pass_down_node(n)
#     end
# end

# marginal_pass_down_node(n::DownFlowΔNode) = () # do nothing
# marginal_pass_down_node(n::DownFlowLeaf) = ()

# function marginal_pass_down_node(n::DownFlow⋀Cached)
#     # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
#      # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
#     for c in n.children
#         for sink in downflow_sinks(c)
#             if !sink.in_progress
#                 sink.downflow .= downflow(n)
#                 sink.in_progress = true
#             else
#                 sink.downflow .+= downflow(n)
#             end
#         end
#     end
# end

# function marginal_pass_down_node(n::DownFlow⋁Cached)
#     # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
#     # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
#     for (ind, c) in enumerate(n.children)
#         for sink in downflow_sinks(c)
#             if !sink.in_progress
#                 sink.downflow .= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)) )
#                 sink.in_progress = true
#             else
#                 sink.downflow .+= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)))
#             end
#         end
#     end
# end

# #### marginal_pass_up_down

# function marginal_pass_up_down(circuit::DownFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
#     @assert !(E isa Bool)
#     # marginal_pass_up(origin(circuit), data)
#     marginal_pass_down(circuit)
# end

#############
###Conditional Upper

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

# function conditional_upper_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
#     resize_flows(circuit, num_examples(data))

#     conditional_upper_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

#     # pr() now returns a vector of triples (a,b,c), where:
#     #  * a:Float64 is the minValue,
#     #  * b:Float64 is the maxValue,
#     #  * c:Int8 is the type of sub-circuit (to check if contains query variables)

#     function conditional_upper_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
#         npr = pr(n)
#         nsize = length(npr)
#         val = zeros(nsize)
#         # Type 0 observed, 1 query, 2 marginalize
#         type = zeros( nsize)

#         #copy values from data input
#         #values form data input
#         #   *1,0 observed variable values;
#         #   *2,3 query variable values, where  2:0 and 3:1
#         #   *-1 unobserved variables (marginalize variables)
#         if positive(n)
#             val = feature_matrix(data)[:,variable(n)]
#         else
#             val =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
#         end

#         # override query variable values by 0 and 1
#         val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0

#         val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1

#         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+2.0)] .= 1

#         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+3.0)] .= 1

#         #  override marginalize variable  by 1
#         val[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
#         type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 2

#         # Construct the tuple to  propagate 4 values : conditional message, marginal upper message, marginal lower message, type
#         pr(n) .= [x for x in zip(val,val,val,type)]
#         return nothing
#     end

#     function conditional_upper_pass_up_node(n::UpFlow⋀Cached, ::PlainXData)
#         npr = pr(n)
#         #Initialize values for tuples
#         val_min = ones( length(npr))
#         val_max = ones( length(npr))
#         type = zeros( length(npr))

#         for i=1:length(n.children)
#             #Obtain values from child i -> Message
#             val_min_i = getindex.(pr(n.children[i]),1) # vector a: minValue len(a)=number of queries
#             val_max_i = getindex.(pr(n.children[i]),2) # vector b: maxValue
#             type_i = getindex.(pr(n.children[i]),3)    # vector c: type of subcircuit

#             #...
#         end

#         pr(n) .= [x for x in zip(val_min,val_max,val_max,type)]
#         return nothing
#     end

#     function conditional_upper_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData)
#         npr = pr(n)
#         val_min = zeros( length(npr))
#         val_max = zeros( length(npr))
#         type = zeros( length(npr))
#         for i=1:length(n.children)
#             val_min_i = getindex.(pr(n.children[i]),1)
#             val_max_i = getindex.(pr(n.children[i]),2)
#             type_i = getindex.(pr(n.children[i]),3)

#             #.....

#         end
#         pr(n) .= [x for x in zip(val_min,val_max,val_max,type)]
#         return nothing
#     end

#     for n in circuit
#         conditional_upper_pass_up_node(n,data)
#     end
#     return nothing
# end

#############
###Conditional Lower

#########

function conditional_lower_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}

    resize_flows(circuit, num_examples(data))

    conditional_lower_pass_up_node(n::UpFlowΔNode, ::PlainXData, mu::Array{Float64,1}) = ()

    # pr() now returns a vector of triples (a,b,c,d), where:
    #  * a: Float64 "val_min" is the conditional lower message
    #  * b: Float64 "marginal_lower" is the marginal lower message ,
    #  * c: Float64 "marginal_upper" is the marginal upper message,
    #  * d: Int8 "type" is the type of sub-circuit (to check if contains query variables)
    # NB: type takes values  0 observed, 1 query, 2 marginalize for marginal inference, but for conditional inference we only need to distinguish query/no query

    function conditional_lower_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}
        npr = pr(n)
        nsize = length(npr)
        val_min = zeros(nsize)     # or ones?      # primo el tuple
        marginal_lower = zeros(nsize) # secondo el tuple !!! e non zeros ! nella prog marginale P(empty obs)=1 !!!
        marginal_upper = zeros(nsize) # terzo el tuple
        type = zeros( nsize)

        # values from data input
        #   *1,0 observed variable values;
        #   *2,3 query variable values, where  2:0 and 3:1
        #   *-1 unobserved variables (marginalize variables)

        # copy values from data input: this takes care of entries concerning observed variables, for other variables see changes below
        # terminal are only literal (top terminal too. what about bot ?)

        if positive(n)

            marginal_lower = zeros(Float64,nsize) + feature_matrix(data)[:,variable(n)]
            marginal_upper = zeros(Float64,nsize) + feature_matrix(data)[:,variable(n)]
        else

            marginal_lower =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
            marginal_upper =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
        end

        # values 2,3 in feature matrix have to be read as -1 for marginal propagation

        marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0
        marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0

        marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
        marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0

        # override query variable entries, depending both on the queried status and the literal status

        if positive(n)

            val_min[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val_min[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
        else

            val_min[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val_min[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0

        end

         # type of queried variables (independent both from the query status and the literal status, correct)

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+2.0)] .= 1

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+3.0)] .= 1

         #  override marginalize variables entries  by 1, independently from queried status and literal status
         val_min[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_lower[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_upper[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1

         # type of variables to marginalize
         type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 2

         # storing messages
         pr(n) .= [x for x in zip(val_min,marginal_lower,marginal_upper,type)]

         return nothing
    end

    function conditional_lower_pass_up_node(n::UpFlow⋀Cached, ::PlainXData, mu::Array{Float64,1})

        npr = pr(n)

        #Initialize values for tuples
        val_min = ones( length(npr))
        marginal_lower = ones( length(npr))
        marginal_upper = ones( length(npr))
        type = zeros( length(npr))
        vect = ones( length(npr))

        for i=1:length(n.children) # for CSDD this is 2

            # marginal messages are the product of children's messages

            marginal_lower .*= getindex.(pr(n.children[i]),2)
            marginal_upper .*= getindex.(pr(n.children[i]),3)

            ##### conditional message (val_min) is the product of pi_lower and marginal_lower/upper of the children, depending of the types of the latter
            # val_min for now has all entries set to 1
            # start by giving pi_lower values when the type of the children i (in which we are in the for cycle) is 1
            val_min[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= getindex.(pr(n.children[i]),1)[getindex.(pr(n.children[i]),4) .== ones( length(npr))]

            ##### type. as soon as I see a children with type 1, set type = 1. otherwise it stays 0. IN AND/OR NODES TYPE IS ONLY 1 OR 0 !
            type[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= 1
        end

        # preparing logical conditions which will tell where to take marginal_lower and where marginal_upper

        cond_query_first_child = getindex.(pr(n.children[1]),4) .== ones( length(npr))

        cond_positive_val_second_child = getindex.(pr(n.children[2]),1) .> zeros( length(npr))
        cond_negative_val_second_child = getindex.(pr(n.children[2]),1) .< zeros( length(npr))

        cond_query_second_child = getindex.(pr(n.children[2]),4) .== ones( length(npr))

        cond_positive_val_first_child = getindex.(pr(n.children[1]),1) .> zeros( length(npr))
        cond_negative_val_first_child = getindex.(pr(n.children[1]),1) .< zeros( length(npr))

        # completing conditional messages

        #when the query is in the sub
        val_min[cond_query_second_child .& cond_positive_val_second_child].*= getindex.(pr(n.children[1]),2)[cond_query_second_child .& cond_positive_val_second_child]
        val_min[cond_query_second_child  .& cond_negative_val_second_child].*= getindex.(pr(n.children[1]),3)[cond_query_second_child .& cond_negative_val_second_child]

        #when the query is in the prime
        val_min[cond_query_first_child .& cond_positive_val_first_child].*= getindex.(pr(n.children[2]),2)[cond_query_first_child .& cond_positive_val_first_child]
        val_min[cond_query_first_child .& cond_negative_val_first_child ].*= getindex.(pr(n.children[2]),3)[cond_query_first_child .& cond_negative_val_first_child]

        # storing messages

        pr(n) .= [x for x in zip(val_min,marginal_lower,marginal_upper,type)]

        return nothing
    end

    function conditional_lower_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData, mu::Array{Float64,1})

        npr = pr(n)
        val_min = zeros( length(npr))
        marginal_lower = zeros( length(npr))
        marginal_upper = zeros( length(npr))
        type = zeros( length(npr))

        # type :  enough to set it equal to any of its children, eg the first one

        type = getindex.(pr(n.children[1]),4)

        # preparing coefficients for the three optimization problems, for each instance

        c_co = Array{Float64}(undef, length(npr), length(n.children))
        c_marg_lo = Array{Float64}(undef, length(npr), length(n.children))
        c_marg_up = Array{Float64}(undef, length(npr), length(n.children))

        for i=1:length(npr)
             c_co[i,:] = getindex.(getindex.(pr.(n.children),i),1) # vettore dei coefficenti del i-esimo problema COND  (i-esima instance)
             c_marg_lo[i,:] = getindex.(getindex.(pr.(n.children),i),2) # vettore dei coefficenti del i-esimo problema MARG_LOWER (i-esima instance)
             c_marg_up[i,:] = getindex.(getindex.(pr.(n.children),i),3) # vettore dei coefficenti del i-esimo problema MARG_UPPER  (i-esima instance)
        end

        # preparing boundaries (they're the same for each problem, for each instance)

        bounds = Array{Float64,2}(undef, 2, length(n.children))
        bounds = hcat(prob_origin(n).log_thetas, prob_origin(n).log_thetas_u)

        # solving the optimization problems

        for i=1:length(npr)
            val_min[i] = lp(c_co[i,:], bounds, false, true )
            marginal_lower[i] = lp(c_marg_lo[i,:], bounds, false, true)
            marginal_upper[i] = lp(c_marg_up[i,:], bounds, true, true )
        end

        # storing messages

        pr(n) .= [x for x in zip(val_min,marginal_lower,marginal_upper,type)]
        return nothing
    end

    for n in circuit
        conditional_lower_pass_up_node(n,data,mu)
    end
    return nothing
end

#############
###Conditional upper

#########

function conditional_upper_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}

    resize_flows(circuit, num_examples(data))

    conditional_upper_pass_up_node(n::UpFlowΔNode, ::PlainXData, mu::Array{Float64,1}) = ()

    # pr() now returns a vector of triples (a,b,c,d), where:
    #  * a: Float64 "val_max" is the conditional upper message
    #  * b: Float64 "marginal_lower" is the marginal lower message ,
    #  * c: Float64 "marginal_upper" is the marginal upper message,
    #  * d: Int8 "type" is the type of sub-circuit (to check if contains query variables)
    # NB: type takes values  0 observed, 1 query, 2 marginalize for marginal inference, but for conditional inference we only need to distinguish query/no query

    function conditional_upper_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}
        npr = pr(n)
        nsize = length(npr)
        val_max = zeros(nsize)     # or ones?      # primo el tuple
        marginal_lower = zeros(nsize) # secondo el tuple !!! e non zeros ! nella prog marginale P(empty obs)=1 !!!
        marginal_upper = zeros(nsize) # terzo el tuple
        type = zeros( nsize)

        # values from data input
        #   *1,0 observed variable values;
        #   *2,3 query variable values, where  2:0 and 3:1
        #   *-1 unobserved variables (marginalize variables)

        # copy values from data input: this takes care of entries concerning observed variables, for other variables see changes below
        # terminal are only literal (top terminal too. what about bot ?)

        if positive(n)

            marginal_lower = zeros(Float64,nsize) + feature_matrix(data)[:,variable(n)]
            marginal_upper = zeros(Float64,nsize) + feature_matrix(data)[:,variable(n)]
        else

            marginal_lower =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
            marginal_upper =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
        end

        # values 2,3 in feature matrix have to be read as -1 for marginal propagation

        marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0
        marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0

        marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
        marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0

        # override query variable entries, depending both on the queried status and the literal status

        if positive(n)

            val_max[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val_max[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
        else

            val_max[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val_max[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0

            # marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0
            # marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0

        end

         # type of queried variables (independent both from the query status and the literal status, correct)

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+2.0)] .= 1

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+3.0)] .= 1

         #  override marginalize variables entries  by 1, independently from queried status and literal status
         val_max[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_lower[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_upper[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1

         # type of variables to marginalize
         type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 2

         # storing messages
         pr(n) .= [x for x in zip(val_max,marginal_lower,marginal_upper,type)]

         return nothing
    end

    function conditional_upper_pass_up_node(n::UpFlow⋀Cached, ::PlainXData, mu::Array{Float64,1})

        npr = pr(n)

        #Initialize values for tuples
        val_max = ones( length(npr))
        marginal_lower = ones( length(npr))
        marginal_upper = ones( length(npr))
        type = zeros( length(npr))
        vect = ones( length(npr))

        for i=1:length(n.children) # for CSDD this is 2

            # marginal messages are the product of children's messages

            marginal_lower .*= getindex.(pr(n.children[i]),2)
            marginal_upper .*= getindex.(pr(n.children[i]),3)

            ##### conditional message (val_max) is the product of pi_upper and marginal_lower/upper of the children, depending on the types of the latter
            # val_max for now has all entries set to 1
            # start by giving pi_upper values when the type of the children i (in which we are in the for cycle) is 1
            val_max[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= getindex.(pr(n.children[i]),1)[getindex.(pr(n.children[i]),4) .== ones( length(npr))]

            ##### type. as soon as I see a children with type 1, set type = 1. otherwise it stays 0. IN AND/OR NODES TYPE IS ONLY 1 OR 0 !
            type[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= 1
        end

        # preparing logical conditions which will tell where to take marginal_lower and where marginal_upper

        cond_query_first_child = getindex.(pr(n.children[1]),4) .== ones( length(npr))

        cond_positive_val_second_child = getindex.(pr(n.children[2]),1) .> zeros( length(npr))
        cond_negative_val_second_child = getindex.(pr(n.children[2]),1) .< zeros( length(npr))

        cond_query_second_child = getindex.(pr(n.children[2]),4) .== ones( length(npr))

        cond_positive_val_first_child = getindex.(pr(n.children[1]),1) .> zeros( length(npr))
        cond_negative_val_first_child = getindex.(pr(n.children[1]),1) .< zeros( length(npr))

        # completing conditional messages

        #when the query is in the sub
        val_max[cond_query_second_child .& cond_positive_val_second_child].*= getindex.(pr(n.children[1]),3)[cond_query_second_child .& cond_positive_val_second_child]
        val_max[cond_query_second_child  .& cond_negative_val_second_child].*= getindex.(pr(n.children[1]),2)[cond_query_second_child .& cond_negative_val_second_child]

        #when the query is in the prime
        val_max[cond_query_first_child .& cond_positive_val_first_child].*= getindex.(pr(n.children[2]),3)[cond_query_first_child .& cond_positive_val_first_child]
        val_max[cond_query_first_child .& cond_negative_val_first_child ].*= getindex.(pr(n.children[2]),2)[cond_query_first_child .& cond_negative_val_first_child]

        # storing messages

        pr(n) .= [x for x in zip(val_max,marginal_lower,marginal_upper,type)]

        return nothing
    end

    function conditional_upper_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData, mu::Array{Float64,1})

        npr = pr(n)
        val_max = zeros( length(npr))
        marginal_lower = zeros( length(npr))
        marginal_upper = zeros( length(npr))
        type = zeros( length(npr))

        # type :  enough to set it equal to any of its children, eg the first one

        type = getindex.(pr(n.children[1]),4)

        # preparing coefficients for the three optimization problems, for each instance

        c_co = Array{Float64}(undef, length(npr), length(n.children))
        c_marg_lo = Array{Float64}(undef, length(npr), length(n.children))
        c_marg_up = Array{Float64}(undef, length(npr), length(n.children))

        for i=1:length(npr)
             c_co[i,:] = getindex.(getindex.(pr.(n.children),i),1) # vettore dei coefficenti del i-esimo problema COND  (i-esima instance)
             c_marg_lo[i,:] = getindex.(getindex.(pr.(n.children),i),2) # vettore dei coefficenti del i-esimo problema MARG_LOWER (i-esima instance)
             c_marg_up[i,:] = getindex.(getindex.(pr.(n.children),i),3) # vettore dei coefficenti del i-esimo problema MARG_UPPER  (i-esima instance)
        end

        # preparing boundaries (they're the same for each problem, for each instance)

        bounds = Array{Float64,2}(undef, 2, length(n.children))
        bounds = hcat(prob_origin(n).log_thetas, prob_origin(n).log_thetas_u)

        # solving the optimization problems

        for i=1:length(npr)
            val_max[i] = lp(c_co[i,:], bounds, true, true )
            marginal_lower[i] = lp(c_marg_lo[i,:], bounds, false, true)
            marginal_upper[i] = lp(c_marg_up[i,:], bounds, true, true )
        end

        # storing messages

        pr(n) .= [x for x in zip(val_max,marginal_lower,marginal_upper,type)]
        return nothing
    end

    for n in circuit
        conditional_upper_pass_up_node(n,data,mu)
    end
    return nothing
end

max_stars_repo_name: IDSIA/cjuice

max_issues_repo_name: IDSIA/cjuice

max_forks_repo_name: IDSIA/cjuice


text:
module LaplaceBIE

using SurfaceTopology
using GeometryTypes
using LinearAlgebra
using FastGaussQuadrature

eye(A::AbstractMatrix{T}) where T = Matrix{eltype(A)}(I,size(A))
eye(m::Integer) = Matrix(1.0I,m,m)
diagm(x) = Matrix(Diagonal(x))

function vertexareas(points,topology)
    vareas = zeros(Float64,length(points))
    for face in Faces(topology)
        v1,v2,v3 = face
        area = norm(cross(points[v2]-points[v1],points[v3]-points[v1])) /2
        vareas[v1] += area/3
        vareas[v2] += area/3
        vareas[v3] += area/3
    end
    return vareas
end

# function surfacefield(points,normals,topology,psi,mup,Htime)
#     H = tangentderivatives(points,normals,topology,psi)
#     Hn = normalderivatives(points,normals,topology,H,mup,Htime)
#     H .= Hn .* normals
#     return H
# end

@generated function strquad(q::Function,x1,x2,x3,NP::Val{T}) where T
    t, w = gausslegendre(T)
    @inbounds quote # Need to test if that improves performance
        B = dot(x3-x1,x2-x1)/norm(x2-x1)^2
        C = norm(x3-x1)^2/norm(x2-x1)^2
        hS = norm(cross(x2-x1,x3-x1))

        s = 0.
        for i in 1:$T
            Chi = pi/4*(1 + $t[i])

            R = 1/(cos(Chi) + sin(Chi))
            si = 0.
            @simd for j in 1:$T
                rho = R/2*(1 + $t[j])
                si += q(rho*cos(Chi),rho*sin(Chi))*$w[j]
            end

            s += si*R/2 / sqrt(cos(Chi)^2 + B*sin(2*Chi) + C*sin(Chi)^2) * $w[i]
        end
        s *= pi/4
        return s*hS/norm(x2 - x1)
    end
end

strquad(q::Function,x1,x2,x3,NP::Int64) = strquad(q,x1,x2,x3,Val(NP))

normalderivatives(points,normals,topology,Ht,hmag,H0; eps=0.0001, NP=10) = normalderivatives(points,normals,topology,Ht,hmag,x->H0; eps=0.0001, NP=10)

"""
    normalderivatives(points,normals,topology,P∇ψ,μ,∇ψ0::Function; eps=0.0001, NP=100)

Calculates ∇ψ.n approached to the object surface from *interior* region. To use the function surface properties `points`, `normals` and `topology` are required. To use the function it is required to have a tangential field components (see `tangentialderivatives` and `surfacepotential`), boundary jump condition μ and a gradient of external free field.
"""
function normalderivatives(points,normals,topology,Ht,hmag,H0::Function; eps=0.0001, NP=100)

    ### Tangential field
    for xkey in 1:length(points)
        nx = normals[xkey]
        Ht[xkey] = (eye(3) - nx*nx')*Ht[xkey]
    end

    vareas = vertexareas(points,topology)

    function qs(xi,eta,v1,v2,v3,x,nx,Htx)
        y = (1 - xi - eta)*points[v1] + xi*points[v2] + eta*points[v3]
        Hty = (1 - xi - eta)*Ht[v1] + xi*Ht[v2] + eta*Ht[v3]
        ny =  (1 - xi - eta)*normals[v1] + xi*normals[v2] + eta*normals[v3]
        s = - dot(nx,cross(Hty - Htx,cross(ny,-(y-x)/norm(y-x)^2)))
    end

    Hn = Array{Float64}(undef,length(points))

    for xkey in 1:length(points)

        nx = normals[xkey]
        x = points[xkey] + eps*nx
        Htx = Ht[xkey]

        s = 0
        for ykey in 1:length(points)
            !(xkey==ykey) || continue
            y = points[ykey]
            ny = normals[ykey]
            Hty = Ht[ykey]

            s += dot(nx,-(Hty-Htx)*dot((y-x)/norm(y-x)^3,ny)) * vareas[ykey]
            s += -dot(nx,cross(Hty-Htx,cross(ny,-(y-x)/norm(y-x)^3))) * vareas[ykey]
        end

        ### Making a proper hole
        for (v2,v3) in EdgeRing(xkey,topology) #DoubleVertexVRing(xkey,faces)
            area = norm(cross(points[v2]-x,points[v3]-x))/2

            ny = normals[v2]
            y = points[v2]
            s -= -dot(nx,cross(Ht[v2]-Htx,cross(ny,-(y-x)/norm(y-x)^3))) * area/3

            ny = normals[v3]
            y = points[v3]
            s -= -dot(nx,cross(Ht[v3]-Htx,cross(ny,-(y-x)/norm(y-x)^3))) * area/3

            ### Singular triangle integration
            s += strquad((xi,eta) -> qs(xi,eta,xkey,v2,v3,x,nx,Htx),x,points[v2],points[v3],NP)
        end

        Hn[xkey] = dot(H0(points[xkey]),nx)/hmag + 1/4/pi * (1-hmag)/hmag * s
    end

    return Hn
end

"""
    tangentderivatives(points,normals,topology,ψ)

Returns tangential derivative P∇ψ for a shape defined by `points`, `normals` and `topology`.
"""
function tangentderivatives(points,normals,topology,psi)

    H = HField(points,topology,psi)
    for xkey in 1:length(points)
        nx = normals[xkey]
        H[xkey] = (eye(3) - nx*nx')*H[xkey]
    end

    return H
end

function HField(points,topology,psi)

    H = Array{Point{3,Float64}}(undef,length(points))

    for xkey in 1:length(points)
        x = points[xkey]
        psix = psi[xkey]

        vvec = Point{3,Float64}[]
        dphi = Float64[]

        for ykey in VertexRing(xkey,topology)
            y = points[ykey]
            push!(vvec,y-x)
            push!(dphi,psi[ykey]-psix)
        end

        A, B = vvec, dphi

        A = Array(hcat(A...))
        A = transpose(reshape(A,3,div(length(A),3))) ### This looks unecesarry
        H[xkey] = inv(transpose(A)*A)*transpose(A)*B
    end

    return H
end

### Potential simple
surfacepotential(points,normals,topology,hmag,H0) = surfacepotential(points,normals,topology,hmag,x->dot(H0,x))

"""
    surfacepotential(points,normals,topology,μ,ψ0::Function)

Returns a surface field for a given shape defined by a triangular mesh with `points`, `normals` and `topology` and for a boundary conditions where ∇ψ has a jump in th e normal direction charectarized with μ and for the external free field given by a potential ψ0.
"""
function surfacepotential(points,normals,topology,hmag,ψ::Function)

    A = zeros(Float64,length(points),length(points))
    vareas = vertexareas(points,topology)

    for xkey in 1:length(points)

        x = points[xkey]
        nx = normals[xkey]

        for ykey in 1:length(points)
            if xkey==ykey
                continue
            end

            y = points[ykey]
            ny = normals[ykey]

            A[ykey,xkey] = dot(y-x,ny)/norm(y-x)^3 * vareas[ykey]
        end
    end

    B = zeros(Float64,length(points))
    for xkey in 1:length(points)
        B[xkey] = 2*ψ(points[xkey])/(hmag+1)
    end

    A = A'
    psi = (eye(A)*(1- (hmag-1)/(hmag+1)) - 1/2/pi * (hmag-1)/(hmag+1) * (A - diagm(Float64[sum(A[i,:]) for i in 1:size(A,2)]))) \ B

    return psi
end

# ### A method for calculating the field energy
# fieldenergy(points,normals,topology,psix,mup,H0) = fieldenergy(points,normals,topology,psix,mup,x->H0)
# function fieldenergy(points,normals,topology,psix,mup,H0::Function)
#     vareas = vertexareas(points,topology)
#     Area = sum(vareas)

#     s = 0
#     for xkey in 1:length(points)
#         s += psix[xkey]*dot(H0(points[xkey]),normals[xkey]) * vareas[xkey]
#     end

#     Em = 1/8/pi * (1 - mup) * s
#     return Em
# end

export surfacepotential, tangentderivatives, normalderivatives #, fieldenergy

end # module


text:
##

using ACEatoms, Test
ES = ACEatoms.Electrostatics

##

@info("Total dipole moment test")
X = [0.0 0.0 0.0; 1.0 0.0 0.0]
Q = [-1.0; 1.0]
MU = [0.0 0.0 0.0; 0.0 1.0 -2.0]
tot_dipole = ES.get_dipole(X, Q, MU)
println(@test(isapprox(tot_dipole, [0.208194 1.0 -2.0], atol=1e-6) ))

##

@info("Soft charge-charge interaction energy test")
X = [0.0 0.0; 1.0 0.0]
q1 = +1.0
q2 = -1.0
E = ES.soft_coulomb(X[1, :], q1, X[2, :], q2, 1.0)
println(@test( isapprox(E, -14.39964547, atol=1e-6) ))

##

@info("Soft charge-dipole interaction energy test")
X = [0.0 0.0; 1.0 0.0]
q1 = 1.0
dip = [1.0 0.0]
E = ES.soft_q_μ(X[1, :], q1, X[2, :], dip, 1.0)
println(@test( isapprox(E, -2.99792458, atol=1e-6) ))

##