text:
module ExtendedFESpacesTests

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Algebra
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.CellData

n = 10
mesh = (n,n)
domain = (0,1,0,1) .- 1
order = 1
model = CartesianDiscreteModel(domain, mesh)

Ω = Triangulation(model)

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldcell_to_coods = get_cell_coordinates(Ω)
oldcell_to_is_in = collect1d(lazy_map(is_in,oldcell_to_coods))

Ω_in = Triangulation(model,oldcell_to_is_in)
@test isa(Ω_in,BodyFittedTriangulation)

@test model === get_background_model(Ω_in)

V = TestFESpace(Ω_in,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
test_single_field_fe_space(V)

U = TrialFESpace(V)
test_single_field_fe_space(U)

oldcell_to_is_out = lazy_map(!,oldcell_to_is_in)
Ω_out = Triangulation(model,oldcell_to_is_out)

Γ = InterfaceTriangulation(Ω_in,Ω_out)

degree = 2
dΩ = Measure(Ω,degree)
dΩ_in = Measure(Ω_in,degree)
dΩ_out = Measure(Ω_out,degree)
dΓ = Measure(Γ,degree)

u(x) = x[1]+1

dv = get_fe_basis(V)
du = get_trial_fe_basis(U)
uh = interpolate(u,U)

x = get_cell_points(Ω)
x_in = get_cell_points(Ω_in)
x_out = get_cell_points(Ω_out)
x_Γ = get_cell_points(Γ)

r = dv(x)
test_array(r,collect(r))
r = dv(x_in)
test_array(r,collect(r))
r = dv(x_out)
test_array(r,collect(r))
r = dv.⁺(x_Γ)
test_array(r,collect(r))

r = ∇(dv)(x)
test_array(r,collect(r))
r = ∇(dv)(x_in)
test_array(r,collect(r))
r = ∇(dv)(x_out)
test_array(r,collect(r))
r = ∇(dv).⁺(x_Γ)
test_array(r,collect(r))

r = du(x)
test_array(r,collect(r))
r = du(x_in)
test_array(r,collect(r))
r = du(x_out)
test_array(r,collect(r))
r = du.⁺(x_Γ)
test_array(r,collect(r))

r = ∇(du)(x)
test_array(r,collect(r))
r = ∇(du)(x_in)
test_array(r,collect(r))
r = ∇(du)(x_out)
test_array(r,collect(r))
r = ∇(du).⁺(x_Γ)
test_array(r,collect(r))

r = uh(x)
test_array(r,collect(r))
r = uh(x_in)
test_array(r,collect(r))
r = uh(x_out)
test_array(r,collect(r))
r = uh.⁺(x_Γ)
test_array(r,collect(r))

r = ∇(uh)(x)
test_array(r,collect(r))
r = ∇(uh)(x_in)
test_array(r,collect(r))
r = ∇(uh)(x_out)
test_array(r,collect(r))
r = ∇(uh).⁺(x_Γ)
test_array(r,collect(r))

#using Gridap.Visualization
#writevtk(Ω,"Omega",cellfields=["uh"=>uh])
#writevtk(Ω_in,"Omega_in",cellfields=["uh"=>uh])
#writevtk(Ω_out,"Omega_out",cellfields=["uh"=>uh])
#writevtk(Γ,"gamma",cellfields=["uh"=>uh.plus])

a(u,v) =
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ +
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ_in +
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ_out +
  ∫(jump(v)⊙jump(u) + jump(∇(v))⊙jump(∇(u)))*dΓ

l(v) =
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ +
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ_in +
  ∫(v⊙u + ∇(v)⊙∇(u))*dΩ_out +
  ∫(jump(v))*dΓ

op = AffineFEOperator(a,l,U,V)

V = TestFESpace(Ω_in,ReferenceFE(lagrangian,Float64,2),conformity=:H1)

V = TestFESpace(Ω_in,ReferenceFE(lagrangian,Float64,2),conformity=:H1,constraint=:zeromean)
uh = FEFunction(V,rand(num_free_dofs(V)))
@test sum(∫(uh)*dΩ_in) + 1 ≈ 1

V = TestFESpace(Ω,ReferenceFE(lagrangian,Float64,2),conformity=:H1)

V_in = TestFESpace(Ω_in,ReferenceFE(lagrangian,Float64,2),conformity=:H1)
V = TestFESpace(Ω,ReferenceFE(lagrangian,Float64,2),conformity=:H1)

vh_in = interpolate(V_in) do x
    x[1]
end
vh_in = interpolate(vh_in, V_in)
vh = interpolate(vh_in, V)

#using Gridap.Visualization
#writevtk(Ω,"Ω",cellfields=["vh"=>vh,"vh_in"=>vh_in])

end # module


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: downsample_xform.jl
#
#  Decription:
#       Downsample MFCC frames along the time axis using fixed ratio (from ARGS)
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-24   wm              Initial version
#
################################################################################

using HDF5

const THIS_DIR = dirname(Base.source_path())

function downsample{_Tp<:FloatingPoint}(
    # code_native(downsample, (Int, Vector{Float64}, Vector{Float64}))
    ratio::Int,
    ivec::AbstractVector{_Tp},
    ovec::AbstractVector{_Tp})

    const STRIDE::Int = ratio
    const ILEN::Int = size(ivec, 1)

    opos::Int = 1
    ipos::Int = 1
    @inbounds @simd for ipos = 1:STRIDE:ILEN - STRIDE + 1
        ovec[opos] = 0
        @simd for sub = 0:STRIDE - 1
            ovec[opos] += ivec[ipos + sub]
        end
        ovec[opos] /= STRIDE

        opos += 1
    end
end

function transform(RATIO::Int, MFCC_H5_FNAME::ASCIIString, MFCC_H5_DSET::ASCIIString, OUT_H5_FNAME::ASCIIString, OUT_H5_DSET::ASCIIString)
    fid = h5open(MFCC_H5_FNAME, "r")
    fullX = fid[MFCC_H5_DSET]

    newX = Vector{Float64}[]

    const NCOLS = size(fullX, 2)

    res = Array(eltype(fullX), div(1001, RATIO))

    for coli = 1:NCOLS
        Xi = reshape(fullX[:, coli], (1001, 13))

        v = similar(Xi, eltype(Xi), div(size(Xi)[1], RATIO), size(Xi)[2])
        for i in 1:size(Xi)[2]
            downsample(RATIO, Xi[:, i], res)
            v[:, i] = res
        end

        push!(newX, vec(v))

        if coli % 100 == 0 println(coli) end
    end

    close(fid)

    h5write(OUT_H5_FNAME, OUT_H5_DSET, hcat(newX...))
end

function main()
    if size(ARGS)[1] != 2
        println("Specify integer downsampling RATIO and select input data (1: training, 2: test, 3: both.")
    else
        RATIO = int(ARGS[1])
        SELECT = int(ARGS[2])
        if SELECT == 1 or SELECT == 3
            OUT_H5_FNAME = "trainXscale$(RATIO).h5"
            if isfile(OUT_H5_FNAME)
                rm(OUT_H5_FNAME)
            end
            transform(RATIO, "trainXmfcc.h5", "lid/train/X/mfcc", OUT_H5_FNAME, "lid/train/X/scale$(RATIO)")
        elseif SELECT == 2 or SELECT == 3
            OUT_H5_FNAME = "testXscale$(RATIO).h5"
            if isfile(OUT_H5_FNAME)
                rm(OUT_H5_FNAME)
            end
            transform(RATIO, "testXmfcc.h5", "lid/test/X/mfcc", OUT_H5_FNAME, "lid/test/X/scale$(RATIO)")
        end
    end
end

if ~isinteractive()
    main()
end

ext: jl


text:
using MixedModels
using Test

using MixedModels: dataset, likelihoodratiotest
using MixedModels: pirls!, setβθ!, setθ!, updateL!

include("modelcache.jl")

# explicitly setting theta for these to so that we can do exact textual comparisons
βθ = [0.1955554704948119,  0.05755412761885973, 0.3207843518569843, -1.0582595252774376,
     -2.1047524824609853, -1.0549789653925743,  1.339766125847893,  0.4953047709862237]
gm3 = GeneralizedLinearMixedModel(only(gfms[:verbagg]), dataset(:verbagg), Bernoulli())
pirls!(setβθ!(gm3, βθ))

fm0θ = [1.1656121258575225]
fm0 = updateL!(setθ!(first(models(:sleepstudy)), fm0θ))

fm1θ = [0.9292213288149662, 0.018168393450877257, 0.22264486671069741]
fm1 = updateL!(setθ!(last(models(:sleepstudy)), fm1θ))

fmreθ = [0.32352483854887326, 0.4715395478019364, 0.0,
         0.43705610601403755, 0.016565641868150047, 0.17732248078617097]
# this is a junk model, but it stresses parts of the display code
fmre = LinearMixedModel(@formula(rt_trunc ~ 1+(0+spkr|subj)+(1+load|item)), MixedModels.dataset(:kb07))
updateL!(setθ!(fmre, fmreθ))
fmre.optsum.feval = 1

lrt = likelihoodratiotest(fm0, fm1)

@testset "markdown" begin
    mime = MIME"text/markdown"()
    @test_logs (:warn, "Model has not been fit: results will be nonsense") sprint(show, mime, gm3)
    gm3.optsum.feval = 1
    @testset "lmm" begin
        @test sprint(show, mime, fm0) == """
|             |     Est. |     SE |     z |      p |  σ_subj |
|:----------- | --------:| ------:| -----:| ------:| -------:|
| (Intercept) | 251.4051 | 9.5062 | 26.45 | <1e-99 | 36.0121 |
| days        |  10.4673 | 0.8017 | 13.06 | <1e-38 |         |
| Residual    |  30.8954 |        |       |        |         |
"""
        @test sprint(show, mime, fm1) == """
|             |     Est. |     SE |     z |      p |  σ_subj |
|:----------- | --------:| ------:| -----:| ------:| -------:|
| (Intercept) | 251.4051 | 6.6323 | 37.91 | <1e-99 | 23.7805 |
| days        |  10.4673 | 1.5022 |  6.97 | <1e-11 |  5.7168 |
| Residual    |  25.5918 |        |       |        |         |
"""
    end

    @testset "re without fe" begin
        @test sprint(show, mime, fmre) == """
|             |      Est. |      SE |     z |      p |   σ_subj |   σ_item |
|:----------- | ---------:| -------:| -----:| ------:| --------:| --------:|
| (Intercept) | 2092.3713 | 76.9426 | 27.19 | <1e-99 |          | 349.7858 |
| spkr: old   |           |         |       |        | 377.3837 |          |
| spkr: new   |           |         |       |        | 258.9242 |          |
| load: yes   |           |         |       |        |          | 142.5331 |
| Residual    |  800.3224 |         |       |        |          |          |
"""
    end

    @testset "glmm" begin
        @test sprint(show, mime, gm3) in ("""
|              |    Est. |     SE |     z |      p | σ_subj | σ_item |
|:------------ | -------:| ------:| -----:| ------:| ------:| ------:|
| (Intercept)  |  0.1956 | 0.4052 |  0.48 | 0.6294 | 1.3398 | 0.4953 |
| anger        |  0.0576 | 0.0168 |  3.43 | 0.0006 |        |        |
| gender: M    |  0.3208 | 0.1913 |  1.68 | 0.0935 |        |        |
| btype: scold | -1.0583 | 0.2568 | -4.12 | <1e-04 |        |        |
| btype: shout | -2.1048 | 0.2585 | -8.14 | <1e-15 |        |        |
| situ: self   | -1.0550 | 0.2103 | -5.02 | <1e-06 |        |        |
""","""
|              |    Est. |     SE |     z |      p | σ_subj | σ_item |
|:------------ | -------:| ------:| -----:| ------:| ------:| ------:|
| (Intercept)  |  0.1956 | 0.4052 |  0.48 | 0.6294 | 1.3398 | 0.4953 |
| anger        |  0.0576 | 0.0168 |  3.43 | 0.0006 |        |        |
| gender: M    |  0.3208 | 0.1913 |  1.68 | 0.0935 |        |        |
| btype: scold | -1.0583 | 0.2568 | -4.12 |  <1e-4 |        |        |
| btype: shout | -2.1048 | 0.2585 | -8.14 | <1e-15 |        |        |
| situ: self   | -1.0550 | 0.2103 | -5.02 |  <1e-6 |        |        |
""")
    end

    @testset "lrt" begin

        @test sprint(show, mime, lrt) in ("""
|                                          | model-dof | deviance |  χ² | χ²-dof | P(>χ²) |
|:---------------------------------------- | ---------:| --------:| ---:| ------:|:------ |
| reaction ~ 1 + days + (1 \\| subj)        |         4 |     1794 |     |        |        |
| reaction ~ 1 + days + (1 + days \\| subj) |         6 |     1752 |  42 |      2 | <1e-09 |
""","""
|                                          | model-dof | deviance |  χ² | χ²-dof | P(>χ²) |
|:---------------------------------------- | ---------:| --------:| ---:| ------:|:------ |
| reaction ~ 1 + days + (1 \\| subj)        |         4 |     1794 |     |        |        |
| reaction ~ 1 + days + (1 + days \\| subj) |         6 |     1752 |  42 |      2 | <1e-9  |
""")
    end

    @testset "blockdescription" begin

        @test sprint(show, mime, BlockDescription(gm3)) == """
| rows | subj     | item       | fixed |
|:---- |:-------- |:---------- |:----- |
| 316  | Diagonal |            |       |
| 24   | Dense    | Diag/Dense |       |
| 7    | Dense    | Dense      | Dense |
"""
    end

    @testset "optsum" begin
        fm1.optsum.feval = 1
        fm1.optsum.initial_step = [0.75, 1.0, 0.75]
        fm1.optsum.finitial = 1784.642296192471
        fm1.optsum.final = [0.9292, 0.0182, 0.2226]
        fm1.optsum.fmin =1751.9393444647023
        out =  sprint(show, mime, fm1.optsum)
        @test startswith(out,"""
|                          |                             |
|:------------------------ |:--------------------------- |
| **Initialization**       |                             |
| Initial parameter vector | [1.0, 0.0, 1.0]             |
| Initial objective value  | 1784.642296192471           |
| **Optimizer settings**   |                             |
| Optimizer (from NLopt)   | `LN_BOBYQA`                 |
| Lower bounds             | [0.0, -Inf, 0.0]            |""")
    end

    @testset "varcorr" begin

        @test sprint(show, mime, VarCorr(fm1)) == """
|          | Column      |  Variance |  Std.Dev | Corr. |
|:-------- |:----------- | ---------:| --------:| -----:|
| subj     | (Intercept) | 565.51069 | 23.78047 |       |
|          | days        |  32.68212 |  5.71683 | +0.08 |
| Residual |             | 654.94145 | 25.59182 |       |
"""

        @test sprint(show, mime, VarCorr(gm3)) == """
|      | Column      |  Variance |  Std.Dev |
|:---- |:----------- | ---------:| --------:|
| subj | (Intercept) |  1.794973 | 1.339766 |
| item | (Intercept) |  0.245327 | 0.495305 |
"""
    end
end

@testset "html" begin
    # this is minimal since we're mostly testing that dispatch works
    # the stdlib actually handles most of the conversion
    @test sprint(show, MIME("text/html"), BlockDescription(gm3)) == """

rows	subj	item	fixed

316	Diagonal

24	Dense	Diag/Dense

7	Dense	Dense	Dense

"""
    optsum = sprint(show, MIME("text/html"), fm0.optsum)

    @test occursin("Initialization", optsum)
    @test occursin("LN_BOBYQA", optsum)
end

@testset "latex" begin
    # this is minimal since we're mostly testing that dispatch works
    # the stdlib actually handles most of the conversion

    b = BlockDescription(gm3)

    @test sprint(show, MIME("text/latex"), b) == """
\\begin{tabular}
{l | l | l | l}
rows & subj & item & fixed \\\\
\\hline
316 & Diagonal &  &  \\\\
24 & Dense & Diag/Dense &  \\\\
7 & Dense & Dense & Dense \\\\
\\end{tabular}
"""

    @test sprint(show, MIME("text/xelatex"), b) == sprint(show, MIME("text/latex"), b)

    @test sprint(show, MIME("text/xelatex"), gm3) != sprint(show, MIME("text/latex"), gm3)

    @test startswith(sprint(show, MIME("text/latex"), gm3),"""
\\begin{tabular}
{l | r | r | r | r | r | r}
 & Est. & SE & z & p & \$\\sigma_\\text{subj}\$ & \$\\sigma_\\text{item}\$ \\\\""")

    # not doing the full comparison here because there's a zero-padded exponent
    # that will render differently on different platforms
    @test startswith(sprint(show, MIME("text/latex"), lrt),
                     "\\begin{tabular}\n{l | r | r | r | r | l}\n & model-dof & deviance & \$\\chi^2\$ & \$\\chi^2\$-dof & P(>\$\\chi^2\$) \\\\")

    optsum = sprint(show, MIME("text/latex"), fm0.optsum)

    @test occursin(raw"\textbf{Initialization}", optsum)
    @test occursin(raw"\texttt{LN\_BOBYQA}", optsum)
end

# return these models to their fitted state for the cache
refit!(fm1; progress=false)
refit!(fm0; progress=false)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export BoolWeight

"""
`BoolWeight(x::Bool)`

| Set           |     ``\\oplus``      |  ``\\otimes``  | ``\\bar{0}`` | ``\\bar{1}`` |
|:-------------:|:--------------------:|:--------------:|:------------:|:------------:|
| ``\\{0,1\\}`` |       ``\\lor``      |  ``\\land``    |   ``0``      |    ``1``     |
"""
struct BoolWeight <: Semiring
  x::Bool
end

zero(::Type{BoolWeight}) = BoolWeight(false)
one(::Type{BoolWeight}) = BoolWeight(true)

*(a::BoolWeight, b::BoolWeight) = BoolWeight(a.x && b.x)
+(a::BoolWeight, b::BoolWeight) = BoolWeight(a.x || b.x)

#properties
isidempotent(::Type{W}) where {W <: BoolWeight} = true
iscommulative(::Type{W}) where {W <: BoolWeight} = true
isleft(::Type{W}) where {W <: BoolWeight}= true
isright(::Type{W}) where {W <: BoolWeight}= true
ispath(::Type{W}) where {W <: BoolWeight}= true
iscomplete(::Type{W}) where {W <: BoolWeight}= true


text:
function sumpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CuArrays.zeros(T, dims.us_dims[1], c)
    scatter_add!(Y, X, cluster)
    Y
end

sumpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    sumpool(CuArray{Int64}(cluster), X, c)

function subpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CuArrays.zeros(T, dims.us_dims[1], c)
    scatter_sub!(Y, X, cluster)
    Y
end

subpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    subpool(CuArray{Int64}(cluster), X, c)

function prodpool(cluster::CuArray{Int}, X::CuArray{T},
                  c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CuArrays.ones(T, dims.us_dims[1], c)
    scatter_mul!(Y, X, cluster)
    Y
end

prodpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    prodpool(CuArray{Int64}(cluster), X, c)

function divpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    FT = floattype(T)
    Y = CuArrays.ones(FT, dims.us_dims[1], c)
    scatter_div!(Y, FT.(X), cluster)
    Y
end

divpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    divpool(CuArray{Int64}(cluster), X, c)

function maxpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CuArrays.fill(typemin(T), dims.us_dims[1], c)
    scatter_max!(Y, X, cluster)
    Y
end

maxpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    maxpool(CuArray{Int64}(cluster), X, c)

function minpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CuArrays.fill(typemax(T), dims.us_dims[1], c)
    scatter_min!(Y, X, cluster)
    Y
end

minpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    minpool(CuArray{Int64}(cluster), X, c)

function meanpool(cluster::CuArray{Int}, X::CuArray{T},
                  c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    FT = floattype(T)
    Y = CuArrays.zeros(FT, dims.us_dims[1], c)
    scatter_mean!(Y, FT.(X), cluster)
    Y
end

meanpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    meanpool(CuArray{Int64}(cluster), X, c)

@adjoint function prodpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    prodpool(cluster, X), function (Δ)
        rev_cluster = gather_indices(cluster)
        ∇X = gather(zero(Δ)+Δ, cluster)
        @inbounds for ind = CartesianIndices(cluster)
            ind = Tuple(ind)
            inds = filter(x -> x != ind, rev_cluster[cluster[ind...]])
            for i = 1:size(X, 1)
                ∇X[i, ind...] *= mapreduce(j -> X[i, j...], *, inds; init=one(T))
            end
        end
        (nothing, ∇X)
    end
end

@adjoint function divpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    divpool(cluster, X), function (Δ)
        rev_cluster = gather_indices(cluster)
        ∇X = -gather(zero(Δ)+Δ, cluster)
        ∇X ./= X.^2
        @inbounds for ind = CartesianIndices(cluster)
            ind = Tuple(ind)
            inds = filter(x -> x != ind, rev_cluster[cluster[ind...]])
            for i = 1:size(X, 1)
                ∇X[i, ind...] /= mapreduce(j -> X[i, j...], *, inds; init=one(T))
            end
        end
        (nothing, ∇X)
    end
end

@adjoint function maxpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    max = maxpool(cluster, X)
    max, function (Δ)
       Δu = (X .== gather(max, cluster))
       Δu .*= gather(zero(Δ)+Δ, cluster)
       (nothing, Δu)
    end
end

@adjoint function minpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    min = minpool(cluster, X)
    min, function (Δ)
       Δu = (X .== gather(min, cluster))
       Δu .*= gather(zero(Δ)+Δ, cluster)
       (nothing, Δu)
    end
end


text:
module ArraySlices

    using ..Misc
    using ..Pointers

    export ArraySlice, shift_end

    """
        ArraySlice(array::AbstractArray, interval::UnitRange)

    Creates a new array slice pointing to the hence array shifted by the starting
    value of the interval.
    This is faster than a view, and has better performance when iterating over it.

    Note: It doesn't allow you to use out of bounds indices
    """
    function ArraySlice(arr::AbstractArray, interval::UnitRange)
        assert_interval(interval)
        unsafe_wrap(typeof(arr),
                    Pointer(arr, interval.start - 1),
                    interval.stop - interval.start + 1,
                    own = false)
    end

    """
        shift_start(sl::AbstractArray, i::Integer)

    Shifts the array starting point in a certain number of items.
    Can be called using `<<`
    Usage:
    ```julia-repl
    julia> arr = [1,2,3,4]
    julia> sl = ArraySlice(arr, 2:3)
    2-element Vector{Int64}:
    2
    3
    julia> sl << 1
    1-element Vector{Int64}:
    3
    julia> sl << -1
    3-element Vector{Int64}:
    1
    2
    3
    ```
    """
    @inline Pointers.shift_start(arr::AbstractArray, i::Integer) =
        ArraySlice(arr, (1 + i):length(arr))
    """
        <<(sl::AbstractArray, i::Integer)

    Same as shift_start.
    """
    Base.:<<(arr::AbstractArray, i::Integer) =
        shift_start(arr, i)

    """
        shift_end(sl::AbstractArray, i::Integer)

    Shifts the array ending point in a certain number of items.
    Can be called using `>>`
    Usage:
    ```julia-repl
    julia> arr = [1,2,3,4]
    julia> sl = ArraySlice(arr, 2:3)
    2-element Vector{Int64}:
    2
    3
    julia> sl >> 1
    3-element Vector{Int64}:
    2
    3
    4
    julia> sl >> -1
    1-element Vector{Int64}:
    2
    ```
    """
    @inline shift_end(arr::AbstractArray, i::Integer) =
        ArraySlice(arr, 1:(length(arr) + i))
    """
        >>(sl::AbstractArray, i::Integer)

    Same as array_slice_shift_end.
    """
    Base.:>>(arr::AbstractArray, i::Integer) =
        shift_end(arr, i)

end


text:
"""
Newmark-beta 2nd order ODE solver
"""
struct Newmark <: ODESolver
  nls::NonlinearSolver
  dt::Float64
  γ::Float64
  β::Float64
end

function solve_step!(
  x1::NTuple{3,AbstractVector},
  solver::Newmark,
  op::ODEOperator,
  x0::NTuple{3,AbstractVector},
  t0::Real,
  cache) # -> (uF,tF)

  dt = solver.dt
  γ = solver.γ
  β = solver.β
  t1 = t0+dt
  u0, v0, a0 = x0
  u1, v1, a1 = x1

  if cache === nothing
    newmark_cache = allocate_cache(op,v0,a0)
    nl_cache = nothing
  else
    newmark_cache, nl_cache = cache
  end

  (v,a, ode_cache) = newmark_cache
  ode_cache = update_cache!(ode_cache,op,t1)
  nlop = NewmarkNonlinearOperator(op,t1,dt,γ,β,(u0,v0,a0),newmark_cache)
  nl_cache = solve!(u1,solver.nls,nlop,nl_cache)

  v1 = γ/(β*dt)*(u1-u0) + (1-γ/β)*v0 + dt*(1-γ/(2*β))*a0
  a1 = 1.0/(β*dt^2)*(u1-u0) - 1.0/(β*dt)*v0 - (1-2*β)/(2*β)*a0

  cache = (newmark_cache, nl_cache)
  x1 = (u1,v1,a1)

  return (x1,t1,cache)

end

"""
Nonlinear operator that represents the Newmark nonlinear operator at a
given time step, i.e., A(t,u_n+1,v_n+1,a_n+1)
"""
struct NewmarkNonlinearOperator <: NonlinearOperator
  odeop::ODEOperator
  t1::Float64
  dt::Float64
  γ::Float64
  β::Float64
  x0::NTuple{3,AbstractVector}
  ode_cache
end

function residual!(b::AbstractVector,op::NewmarkNonlinearOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  residual!(b,op.odeop,op.t1,(u1,v1,a1),cache)
end

function jacobian!(A::AbstractMatrix,op::NewmarkNonlinearOperator,x::AbstractVector)
  u1 = x
  u0, v0, a0 = op.x0
  v1, a1, cache = op.ode_cache
  a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
  v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
  z = zero(eltype(A))
  fillstored!(A,z)
  jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
end

function allocate_residual(op::NewmarkNonlinearOperator,x::AbstractVector)
  v1, a1, cache = op.ode_cache
  allocate_residual(op.odeop,x,cache)
end

function allocate_jacobian(op::NewmarkNonlinearOperator,x::AbstractVector)
  v1, a1, cache = op.ode_cache
  allocate_jacobian(op.odeop,x,cache)
end

function zero_initial_guess(op::NewmarkNonlinearOperator)
  x0 = similar(op.u0)
  fill!(x0,zero(eltype(x0)))
  x0
end


text:
for F in (Float32, Float64)
	orig_up = Hermitian(rand(Complex{F}, 25, 25))
	orig_dn = Hermitian(rand(Complex{F}, 25, 25))
	colin = DFW.ColinMatrix(orig_up, orig_dn)
	normal_eig1, normal_eig2 = eigen(orig_up), eigen(orig_dn)
	cache      = DFW.EigCache(colin)
	cached_eig = eigen(colin, cache)
	@test sum(normal_eig1.values) + sum(normal_eig2.values) ≈ sum(cached_eig.values)
	@test Array(normal_eig1) ≈ Array(cached_eig)[1:25, 1:25] ≈ DFW.up(colin)
	@test Array(normal_eig2) ≈ Array(cached_eig)[1:25, 26:50] ≈ DFW.down(colin)
end

for F in (Float32, Float64)
	t = rand(Complex{F}, 50, 50)
	orig       = (t + t')/2
	normal_eig = eigen(orig)
	cache      = DFW.EigCache(orig)
	cached_eig = eigen(orig, cache)
	@test sum(normal_eig.values) ≈ sum(cached_eig.values)
	@test Array(normal_eig) ≈ Array(cached_eig) ≈ orig
end

mat = rand(20, 20)
ats = [DFControl.Atom(:Fe, DFControl.element(:Fe), DFControl.Point3(0.0, 0.0,0.0).*DFControl.angstrom, DFControl.Point3(0.0,0.0,0.0), projections=[DFControl.Projection(DFControl.orbital(:d), 1, 10)]), DFControl.Atom(:Fe, DFControl.element(:Fe), DFControl.Point3(1.0, 0.0,0.0).*DFControl.angstrom, DFControl.Point3(1.0,0.0,0.0), projections=[DFControl.Projection(DFControl.orbital(:d), 11, 20)])]
noncolinmat = convert(DFW.NonColinMatrix, mat)

@test mat[1] == noncolinmat[1]
@test mat[2, 2] == noncolinmat[11, 11]
@test mat[6, 6] == noncolinmat[13, 13]
@test mat[11, 11] == noncolinmat[6, 6]
@test mat[1, 2] == noncolinmat[1, 11]
@test mat[2, 1] == noncolinmat[11, 1]
@test mat[1, 6] == noncolinmat[1, 13]
@test mat[11, 1] == noncolinmat[6, 1]

@test noncolinmat[ats[1]] == [noncolinmat[1:5,1:5] noncolinmat[1:5, 11:15];noncolinmat[11:15, 1:5] noncolinmat[11:15, 11:15]]
@test noncolinmat[ats[1], ats[2]] == [noncolinmat[1:5,6:10] noncolinmat[1:5, 16:20];noncolinmat[11:15, 6:10] noncolinmat[11:15, 16:20]]

@test noncolinmat[ats[1], ats[2], DFW.Up()] == noncolinmat[1:5,6:10]
@test noncolinmat[ats[1], ats[2], DFW.Down()] == noncolinmat[11:15, 16:20]

@test noncolinmat[ats[1], ats[2], DFW.Up(), DFW.Down()] == noncolinmat[1:5,16:20]
@test noncolinmat[ats[1], ats[2], DFW.Down(), DFW.Up()] == noncolinmat[11:15, 6:10]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

cmanpath = joinpath(dirname(@__FILE__), "UnixDomainCM.jl")
include(cmanpath)

npids = addprocs(UnixDomainCM(2))
assert(length(npids) == 2)
test_pids = [remotecall_fetch(myid, x) for x in npids]
assert(npids == test_pids)
rmprocs(npids; waitfor=1.0)

exit(0)

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using FileManager
using Visualization

# potree source CAVA
potree = "C:/Users/marte/Documents/potreeDirectory/pointclouds/CAVA" # replace this path with local potree directory
# FileManager.readhrc(potree)
trie = FileManager.potree2trie(potree)
depth = FileManager.max_depth(trie)
all_files = FileManager.get_all_values(trie)

PC = FileManager.las2pointcloud(all_files...)

# point cloud
Visualization.VIEW(
    [
    Visualization.points(PC.coordinates,PC.rgbs)
    ]
)

# subtrie
sub_trie = FileManager.sub_trie(trie,"r00")
all_files = FileManager.get_all_values(sub_trie)

PC = FileManager.las2pointcloud(all_files...)

# part of point cloud
Visualization.VIEW(
    [
    Visualization.points(PC.coordinates,PC.rgbs)
    ]
)

# cuttrie
cut_trie = deepcopy(trie)
FileManager.cut_trie!(cut_trie,1)
depth = FileManager.max_depth(cut_trie)
all_files = FileManager.get_all_values(cut_trie)

PC = FileManager.las2pointcloud(all_files...)

# first level of detail of point cloud
Visualization.VIEW(
    [
    Visualization.points(PC.coordinates,PC.rgbs)
    ]
)

# potree source CASALETTO
potree = "C:/Users/marte/Documents/potreeDirectory/pointclouds/CASALETTO" # replace this path with local potree directory
trie = FileManager.potree2trie(potree)
depth = FileManager.max_depth(trie)

# cuttrie
cut_trie = deepcopy(trie)
FileManager.cut_trie!(cut_trie,0)
depth = FileManager.max_depth(cut_trie)
all_files = FileManager.get_all_values(cut_trie)

PC = FileManager.las2pointcloud(all_files...)

# root of potree
Visualization.VIEW(
    [
    Visualization.points(PC.coordinates,PC.rgbs)
    ]
)


text:
module str

    export is_unique1, is_unique2

    #=
        Is Unique.

        Implement an algorithm to determine if a string has
        all unique characters. What if you cannot use additional
        data structures?
    =#

    function is_unique1(str::String)::Bool
        return length(Set(str)) == length(str)
    end

    # Instead of a data structure we just use
    # a bit flag to track the characters we've
    # seen thusfar.
    function is_unique2(str::String)::Bool

        seen, offset = Int(0), Int('a')

        for index in firstindex(str):lastindex(str)
            char = str[index]

            isspace(char) && continue

            mask = 1 << (Int(char) - offset)
            seen & mask == mask && return false

            seen |= mask
        end

        return true
    end

    #=
        Check Permutation.

        Given two strings, write a function to decide if one is
        a permutation of the other.
    =#

    export is_permutation1, is_permutation2

    # O(S₁log(S₁) + S₂log(S₂))
    function is_permutation1(s₁::String, s₂::String)

        char_sort(s) = string(sort(collect(s)))

        length(s₁) != length(s₂) && return false

        return char_sort(s₁) == char_sort(s₂)
    end

    # O(S₁ + S₂)
    function is_permutation2(s₁::String, s₂::String)::Bool

        function to_dict(str::String)::Dict{Char, Integer}
            dict = Dict{Char, Integer}()

            for char in str
                if char ∉ keys(dict)
                    dict[char] = 0
                else
                    dict[char] += 1
                end
            end

            return dict
        end

        n₁ = length(s₁); n₂ = length(s₂)

        n₁ != n₂ && return false

        d₁, d₂ = (if (n₁ < n₂)
                    to_dict(s₁), to_dict(s₂)
                 else
                     to_dict(s₂), to_dict(s₁)
                end)

        keys(d₁) != keys(d₂) && return false

        for k in keys(d₁)
            d₁[k] != d₂[k] && return false
        end

        return true
    end

    #=
        URLify.

        Write a method to replacea all spaces in a string with
        '%20'. You may assume that the string has sufficient space
        at the end to hold the additional characters, and that you
        are given the true length of the string.

        EXAMPLE:

        Input: "Mr John Smith    ", 13
        Ouput: "Mr%20%John%20Smith"
    =#

    export urlify1

    # O(N), where N = len(S)
    function urlify1(str::String, len::Integer)::String

        function count_spaces(arr::Array{Char, 1})::Integer
            counter = 0
            for char in arr
                isspace(char) && (counter += 1)
            end
            return counter
        end

        arr = collect(str)
        total_spaces = count_spaces(arr)
        actual, position = length(arr), len

        for index in len:-1:1
            if isspace(arr[position])
                arr[actual] = '0'
                arr[actual-1] = '2'
                arr[actual-2] = '%'

                actual -= 2
            else
                arr[actual] = arr[position]
            end

            actual -= 1; position -= 1
        end

        return String(arr)
    end

    #=
        Palindrome Permutation.

        Given a string, write a function to check if it is
        a permutation of a palindrome. A palindrome is a word
        or phrase that is the same forwards and backwards. A
        permutation is a rearrangement of letters. The palindrome
        does not need to be limited to just dictionary words.

        EXAMPLE:

        Input: Tact Coa
        Output: True (permutations: "taco cat", "atco cta", etc.)
    =#

    export is_pal_perm1

    # O(N), where N = len(S)
    function is_pal_perm1(str::String)::Bool

        function char_counts(str::String)::Array{Integer}
            counts = Dict{Char, Integer}()

            for c in lowercase(str)

                !isletter(c) && continue

                if c ∉ keys(counts)
                    counts[c] = 1
                else
                    counts[c] += 1
                end
            end

            return collect(values(counts))
        end

        isempty(str) && return false

        return count(x -> isodd(x), char_counts(str)) <= 1
    end

    #=
        One Away.

        There are three types of edits that can be performed on
        strings: insert a character, remove a character, or
        replace a character. Given two strings, write a function
        to check if they are one edit (or zero edits) away.

        EXAMPLE:

        pale, ple -> true
        pales, pale -> true
        pale, bale -> true
        pale, bake -> false
    =#

    export one_away1

    function one_away1(s₁::String, s₂)::Bool
        function to_dict(str::String)::Dict{Char, Integer}
            dict = Dict{Char, Integer}()

            for c in lowercase(str)

                !isletter(c) && continue

                if c ∉ keys(dict)
                    dict[c] = 1
                else
                    dict[c] += 1
                end
            end

            return dict
        end

        (isempty(s₁) || isempty(s₂)) && return false

        d₁, d₂ = (if (length(s₁) > length(s₂))
                    to_dict(s₁), to_dict(s₂)
                else
                    to_dict(s₂), to_dict(s₁)
                end)

        δ = 0

        for k in keys(d₁)
            (k ∈ keys(d₂) && d₁[k] == d₂[k]) && continue
            δ += d₁[k]

            δ > 1 && return false
        end

        return true
    end

    #=
        String Compression.

        Implement  a method to perform basic string compression using
        the counts of repeated characters. For example, the string
        "aabcccccaaa" would become a2b1c5a3. If the "compressed" string
        would not become smaller than the original string, your method
        should return the original string. You can assume the string
        has only uppercase and lowercase letters (a-z.)

    =#

    using Printf

    export compress1

    function compress1(str::String)::String
        n = length(str)
        n == 0 && return str

        prev, counter = str[1], 1
        compressed = string()

        for index in 2:n
            current = str[index]

            if prev == current
                counter += 1
            else
                compressed = @sprintf "%s%s%i" compressed prev counter
                counter = 1
            end
            prev = current
        end

        compressed = @sprintf "%s%s%i" compressed prev counter

        l₁, l₂ = length(str), length(compressed)

        return if l₁ <= l₂ str else compressed end
    end

    #=
        String Rotation.

        Assume you have a method isSubstring which checks if one
        word is a substring of another. Given two strings,
        s₁ and s₂, write code to check if s₂ is a rotation
        of s₁ using only one call to isSubstring (e.g,
        "waterbottle" is a rotation of "erbottlewat").
    =#

    # internal package function for is_rotation
    function is_substring(s₁::String, s₂::String)::Bool
        return occursin(s₁, s₂)
    end

    export is_rotation1

    #=
    Brute Force to meet requirments of the "is substring" single
    call constraint.

    O(N) where N = length of s₁
    =#
    function is_rotation1(s₁::String, s₂::String)::Bool

        function get_rotations(str::String, delim::Char='|')::String
            rotations = str * delim

            for index in firstindex(str):lastindex(str)

                prefix = str[index + 1:end]
                suffix = str[1:index]

                rotations *= prefix * suffix * delim
            end

            return rotations
        end

        (length(s₁) == 0 || length(s₂) == 0) && return false

        rotations = get_rotations(s₁)

        return is_substring(s₂, rotations)
    end

end # end module

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
"""
    SOCtoRSOCBridge{T} <: Bridges.Variable.AbstractBridge

Same transformation as [`MOI.Bridges.Constraint.SOCRBridge`](@ref).
"""
struct SOCtoRSOCBridge{T} <: AbstractBridge
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    }
end
function bridge_constrained_variable(
    ::Type{SOCtoRSOCBridge{T}},
    model::MOI.ModelLike,
    set::MOI.SecondOrderCone,
) where {T}
    variables, constraint = MOI.add_constrained_variables(
        model,
        MOI.RotatedSecondOrderCone(MOI.dimension(set)),
    )
    return SOCtoRSOCBridge{T}(variables, constraint)
end

function supports_constrained_variable(
    ::Type{<:SOCtoRSOCBridge},
    ::Type{MOI.SecondOrderCone},
)
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SOCtoRSOCBridge})
    return [(MOI.RotatedSecondOrderCone,)]
end
function MOIB.added_constraint_types(::Type{<:SOCtoRSOCBridge})
    return Tuple{DataType,DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::SOCtoRSOCBridge, ::MOI.NumberOfVariables)
    return length(bridge.variables)
end
function MOI.get(bridge::SOCtoRSOCBridge, ::MOI.ListOfVariableIndices)
    return bridge.variables
end
function MOI.get(
    bridge::SOCtoRSOCBridge,
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.RotatedSecondOrderCone},
)
    return 1
end
function MOI.get(
    bridge::SOCtoRSOCBridge,
    ::MOI.ListOfConstraintIndices{
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    },
)
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::SOCtoRSOCBridge)
    return MOI.delete(model, bridge.variables)
end

# Attributes, Bridge acting as a constraint

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::SOCtoRSOCBridge{T},
) where {T}
    return MOI.SecondOrderCone(length(bridge.variables))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintDual},
    bridge::SOCtoRSOCBridge,
)
    return rotate_result(model, attr, bridge.constraint)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimal,
    bridge::SOCtoRSOCBridge,
    i::IndexInVector,
)
    return rotate_result(model, attr, bridge.variables, i)
end

function MOIB.bridged_function(
    bridge::SOCtoRSOCBridge{T},
    i::IndexInVector,
) where {T}
    return rotate_bridged_function(T, bridge.variables, i)
end
function unbridged_map(
    bridge::SOCtoRSOCBridge{T},
    vis::Vector{MOI.VariableIndex},
) where {T}
    return rotate_unbridged_map(T, bridge.variables, vis)
end


text:
using Pkg
deps = ["Yao", "DelimitedFiles", "FileIO", "Fire", "JLD2", "KrylovKit", "StatsBase"]
extras = ["CUDAnative", "CuArrays"]

USE_CUDA = !("nocuda" in ARGS)

if USE_CUDA
    deps = vcat(deps, extras)
end

for x in deps
    println("Installing $x ...")
    Pkg.add(x)
end

if USE_CUDA
    println("Installing CuYao ...")
    Pkg.clone("https://github.com/QuantumBFS/CuYao.jl")
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text: println("hello World")


text:
export ForwardSensCmeAlgorithm

abstract type ForwardSensCmeAlgorithm end

include("sensoutputs.jl")
include("forwardsensspaceadapters.jl")
include("sparse/forwardsenscmesparse.jl")


text:
using FTCTests
const FTC = FaultTolerantControl
using LinearAlgebra
using Plots
using UnPack
using Debugger
using ReferenceFrameRotations

"""
# Refs
[1] M. Tahavori and A. Hasan, “Fault recoverability for nonlinear systems with application to fault tolerant control of UAVs,” Aerosp. Sci. Technol., vol. 107, p. 106282, 2020, doi: 10.1016/j.ast.2020.106282.

# Notes
- x = state vector, [eta, omega, v, xi]^T
    - eta = Euler angle vector, [phi, theta, psi]^T
    - omega = angular rate vector, [p, q, r]^T
    - v = velocity vector, [u, v, w]^T
    - xi = position vector, [x, y, z]^T
- u = control input vector, thrust vector, [omega_1, omega_2, omega_3, omega_4]^T
- y = output vector, [eta, omega, xi]^T
"""
function compute_minHSV_example(lambda, num::Int; dt=0.01, tf=1.0)
    @assert lambda >= 0.0 && lambda <= 1.0
    Λ = diagm(ones(4))
    Λ[num, num] = lambda
    # @show Λ
    # Quadcopter state space model [1]
    gc = 9.81
    m = 0.65
    l = 0.023
    k = 31.1e-5
    b = 7.5e-7
    J = Diagonal([7.5e-3, 7.4e-3, 1.3e-2])
    Jinv = inv(J)
    A = [zeros(3, 3) Matrix(I, 3, 3) zeros(3, 3) zeros(3, 3);
    zeros(3, 3) zeros(3, 3) zeros(3, 3) zeros(3, 3);
    zeros(3, 3) zeros(3, 3) zeros(3, 3) zeros(3, 3);
    zeros(3, 3) zeros(3, 3) Matrix(I, 3, 3) zeros(3, 3)]
    B = zeros(12, 4)
    for i = 1:12
        for j = 1:4
            if i == 5
                if j == 1
                    B[i, j] = -l*k/J[2, 2]
                elseif j == 3
                    B[i, j] = l*k/J[2, 2]
                end
            elseif i == 4
                if j == 2
                    B[i, j] = -l*k/J[1, 1]
                elseif j == 4
                    B[i, j] = l*k/J[1, 1]
                end
            elseif i == 6
                if j == 1 || j == 3
                    B[i, j] = -b/J[3, 3]
                elseif j == 2 || j == 4
                    B[i, j] = b/J[3, 3]
                end
            elseif i == 9
                if j == 1 || j == 2 || j == 3 || j == 4
                    B[i, j] = k/m
                end
            else
                B[i, j] = 0
            end
        end
    end
    C = [Matrix(I, 3, 3) zeros(3, 3) zeros(3, 3) zeros(3, 3);
    zeros(3, 3) Matrix(I, 3, 3) zeros(3, 3) zeros(3, 3);
    zeros(3, 3) zeros(3, 3) zeros(3, 3) Matrix(I, 3, 3)]

    T(x) = [1 sin(x[1])*tan(x[2]) cos(x[1])*tan(x[2]);
            0 cos(x[1]) -sin(x[1]);
            0 sec(x[2])*sin(x[1]) cos(x[1])*sec(x[2])]
    R(x) = [cos(x[2])*cos(x[3]) sin(x[1])*sin(x[2])*cos(x[3])-cos(x[1])*sin(x[3]) cos(x[1])*sin(x[2])*cos(x[3])+sin(x[1])*sin(x[3]);
            cos(x[2])*sin(x[3]) sin(x[1])*sin(x[2])*sin(x[3])-cos(x[1])*cos(x[3]) cos(x[1])*sin(x[2])*sin(x[3])-sin(x[1])*cos(x[3]);
            -sin(x[2]) sin(x[1])*cos(x[2]) cos(x[1])*cos(x[2])]
    F(x) = vec([
        T(x) * x[4:6]
        Jinv * cross(-x[4:6], J * x[4:6])
        -cross(x[4:6], x[7:9]) + gc * R(x)' * [0; 0; 1]
        R(x) * x[7:9]
    ])

    f(x, u, p, t) = A*x + F(x) + B*(Λ*u).^2
    g(x, u, p, t) = C*x

	# System dimension
	n, m = size(B)
	l = size(C)[1]

    # Initial settings
    x0 = zeros(n)
    u0 = zeros(m)
    pr = zeros(4, 1)

    Wc = FTC.empirical_gramian(f, g, m, n, l; opt=:c, dt=dt, tf=tf, pr=pr, xs=x0, us=u0)
    Wo = FTC.empirical_gramian(f, g, m, n, l; opt=:o, dt=dt, tf=tf, pr=pr, xs=x0, us=u0)
    # minHSV = FTC.min_HSV(Wc, Wo)
    eigvals_Wc = Wc |> LinearAlgebra.eigvals |> minimum |> sqrt
end

# function Dynamics!(multicopter::IslamQuadcopter)
#     @unpack B = multicopter
#     function dynamics!(dx, x, param, t; u, Λ)
#         ν = B * Λ * u
#         f, M = ν[1], ν[2:4]
#         @nested_log FSimZoo.__Dynamics!(multicopter)(dx, x, (), t; f=f, M=M)
#     end
# end

# function test_model()
#     multicopter = IslamQuadcopter()
#     X = State(multicopter)()
#     Λ = diagm(ones(4))
#     u = (multicopter.m * multicopter.g / multicopter.kf) / 4 * ones(4)
#     dX = State(multicopter)()
#     Dynamics!(multicopter)(dX, X, (), 0.0; u=u, Λ=Λ)
#     dX
# end

# function compute_minHSV_Islam(lambda, num::Int; dt=0.01, tf=1.0)
#     @assert lambda >= 0.0 && lambda <= 1.0
#     Λ = diagm(ones(4))
#     Λ[num, num] = lambda
#     @show Λ
#     multicopter = IslamQuadcopter()
#     function f(x, u, param, t)
#         p = x[1:3]
#         v = x[4:6]
#         quat = x[7:10]
#         ω = x[11:13]
#         R = quat_to_dcm(Quaternion(quat...))
#         X = State(multicopter)(p, v, R, ω)
#         dX = State(multicopter)()  # initialisation
#         Dynamics!(multicopter)(dX, X, param, t; u=u, Λ=Λ)
#         dX_quat = dcm_to_quat(DCM(dX.R))
#         dx = [dX.p..., dX.v..., dX_quat.q0, dX_quat.q1, dX_quat.q2, dX_quat.q3, dX.ω...]
#         return dx
#     end
#     g(x, u, p, t) = x[7:10]  # different from [1]
#     n, m, l = 13, 4, 4
#     X0 = State(multicopter)()
#     quat = dcm_to_quat(DCM(X0.R))
#     _quat = [quat.q0, quat.q1, quat.q2, quat.q3]
#     x0 = [X0.p..., X0.v..., _quat..., X0.ω...]
#     # u0 = (multicopter.m * multicopter.g / multicopter.kf) / 4 * ones(4)
#     u0 = zeros(4)
#     pr = zeros(4, 1)

#     Wc = FTC.empirical_gramian(f, g, m, n, l; opt=:c, dt=dt, tf=tf, pr=pr, xs=x0, us=u0, xm=1.0, um=10000.0)
#     Wo = FTC.empirical_gramian(f, g, m, n, l; opt=:o, dt=dt, tf=tf, pr=pr, xs=x0, us=u0, xm=1.0, um=10000.0)
#     _eigvals_Wc = Wc |> LinearAlgebra.eigvals
#     eigvals_Wc = []
#     for eigval in _eigvals_Wc
#         if abs(imag(eigval)) < 1e-6
#             push!(eigvals_Wc, real(eigval))
#         else
#             error("Too large imaginary part")
#         end
#     end
#     min_HSV = eigvals_Wc |> minimum |> sqrt # not min_HSV
#     # minHSV = FTC.min_HSV(Wc, Wo)
# end

function plotting(rotor_idx)
    lambda = 0:0.10:1 |> collect
    HSVs = []
    for i = 1:length(lambda)
        HSVs = push!(HSVs, compute_minHSV_example(lambda[i], rotor_idx))
        # HSVs = push!(HSVs, compute_minHSV_Islam(lambda[i], rotor_idx))
    end
    return plot(lambda,
                HSVs,
                xlabel = "Actuator effectiveness",
                ylabel = "Minimum HSV",
                label = nothing,
               )
end

function main()
    figs = []
    for rotor_idx in 1:4
        fig = plotting(rotor_idx)
        figs = push!(figs, fig)
    end
    plot(figs..., layout=(2, 2))
end


text:
module Multisets

import Base: show, length, getindex, collect, union, intersect, isempty, first
import Base: push!, setindex!, delete!, hash, eltype
import Base: (*), (+), (-)
import Base: issubset, Set, (==)

export Multiset, set_short_show, set_julia_show, set_braces_show, set_key_value_show

"""
A `Multiset` is an unordered collection of things with repetition permitted.
A new `Multiset` container is created with `Multiset{T}()` where `T` is the
type of the objects held in the multiset. If `T` is omitted, it defaults
to `Any`.

A `Multiset` can be created from a collection `list` (such as a `Vector` or
`Set`) with `Multiset(list)`. If an element is repeated in `list` it has
the appropriate multiplicity.
"""
struct Multiset{T} <: AbstractSet{T}
    data::Dict{T,Int}
    function Multiset{T}() where {T}
        d = Dict{T,Int}()
        new(d)
    end
end
Multiset() = Multiset{Any}()
Multiset(x...) = Multiset(collect(x))

function Base.copy(M::Multiset{T}) where {T}
    newMultiset = Multiset{T}()
    for (key, value) in pairs(M)
        newMultiset[key] = value
    end
    return newMultiset
end

function Base.empty!(M::Multiset{T}) where {T}
    for (key, value) in pairs(M)
        M[key] = 0
    end
    clean!(M)
    return M
end

eltype(M::Multiset{T}) where {T} = T

function Multiset(list::AbstractArray{T,d}) where {T,d}
    M = Multiset{T}()
    for x in list
        push!(M, x)
    end
    return M
end

function Multiset(A::Base.AbstractSet{T}) where {T}
    M = Multiset{T}()
    for x in A
        push!(M, x)
    end
    return M
end

"""
`clean!(M)` removes elements of multiplicy 0 from the underlying data
structure supporting `M`.
"""
function clean!(M::Multiset)
    for x in keys(M.data)
        if M[x] == 0
            delete!(M.data, x)
        end
    end
    nothing
end

"""
For a `M[t]` where `M` is a `Multiset` returns the
multiplicity of `t` in `M`. A value of `0` means that
`t` is not a member of `M`.
"""
function getindex(M::Multiset{T}, x)::Int where {T}
    if haskey(M.data, x)
        return M.data[x]
    end
    return 0
end

"""
`push!(M,x,incr)` increases the multiplicity of `x` in `M`
by `incr` (which defaults to 1). `incr` can be negative, but
it is not possible to decrease the multiplicty below 0.
"""
function push!(M::Multiset{T}, x, incr::Int = 1) where {T}
    if haskey(M.data, x)
        M.data[x] += incr
    else
        M.data[x] = incr
    end
    if M.data[x] < 0
        M.data[x] = 0
    end
    return M
end

function setindex!(M::Multiset{T}, m::Int, x) where {T}
    mm = max(m, 0)
    M.data[x] = mm
    return mm
end

function delete!(M::Multiset, x)
    if haskey(M.data, x)
        delete!(M.data, x)
    end
    return M
end

function length(M::Multiset)
    total = 0
    for v in values(M.data)
        total += v
    end
    return total
end

isempty(M::Multiset) = length(M) == 0

function collect(M::Multiset{T}) where {T}
    n = length(M)
    result = Array{T,1}(undef, n)   #  Vector{T}(n)
    i = 0
    for (k, v) in M.data
        for _ = 1:v
            i += 1
            result[i] = k
        end
    end
    try
        sort!(result)
    catch
    end
    return result
end

function braces_string(M::Multiset{T}) where {T}
    elts = collect(M)
    n = length(elts)
    if n == 0
        return "∅"
    end
    str = "{"
    for k = 1:n
        str *= string(elts[k])
        if k < n
            str *= ","
        end
    end
    str *= "}"
    return str
end

short_string(M::Multiset{T}) where {T} = "Multiset{$T} with $(length(M)) elements"

key_value_string(M::Multiset{T}) where {T} =
    "Multiset{$T}($(mapreduce(string, (l,r) -> l == "" ? r : l*", "*r, pairs(M), init="")))"

function julia_string(M::Multiset{T}) where {T}
    elts = collect(M)
    n = length(elts)
    str = "Multiset($T["
    q = ""
    if T <: AbstractString
        q = "\""
    end
    for k = 1:n
        str *= q * string(elts[k]) * q
        if k < n
            str *= ","
        end
    end
    str *= "])"
    return str
end

# This variable controls printing:
# 0 -- {x,y,z}
# 1 -- Multiset{T} with n elements
# 2--  Multiset{T}(x,y,z)

const multi_show_braces = 0
const multi_show_short = 1
const multi_show_julia = 2
const multi_show_key_count = 3
multi_show_flag = multi_show_braces

"""
Set show display mode for multisets, like this:

`Multiset{Int64} with 7 elements`

See also `set_braces_show`,`set_key_value_show` and `set_julia_show`.
"""
set_short_show() = (global multi_show_flag = multi_show_short; nothing)

"""
Set braces display mode for multisets, like this:

`{1,2,2,3,3,3,3}`

See also `set_short_show`, `set_key_value_show` and `set_julia_show`.
"""
set_braces_show() = (global multi_show_flag = multi_show_braces; nothing)

"""
Set Julia style display mode for multisets, like this:

`Multiset(Int64[1,2,2,3,3,3,3])`

See also `set_short_show`, `set_key_value_show` and `set_braces_show`.
"""
set_julia_show() = (global multi_show_flag = multi_show_julia; nothing)

"""
Set key-count style display mode for multisets, like this:

`Multiset{Int64}(1 => 1, 2 => 2, 3 => 4)`

See also `set_short_show`, `set_julia_show` and `set_braces_show`.
"""
set_key_value_show() = (global multi_show_flag = multi_show_key_count; nothing)

function show(io::IO, M::Multiset)
    if multi_show_flag == multi_show_short
        print(io, short_string(M))
    end
    if multi_show_flag == multi_show_braces
        print(io, braces_string(M))
    end
    if multi_show_flag == multi_show_julia
        print(io, julia_string(M))
    end
    if multi_show_flag == multi_show_key_count
        print(io, key_value_string(M))
    end
end

show(M::Multiset) = show(stdout, M)
import Base.Multimedia.display
display(M::Multiset) = print(string(M))

function hash(M::Multiset, h::UInt = UInt(0))
    clean!(M)
    return hash(M.data, h)
end

function Set(M::Multiset{T}) where {T}
    iter = (x for x in keys(M.data) if M.data[x] > 0)
    return Set{T}(iter)
end

function first(M::Multiset{T})::T where {T}
    if length(M) == 0
        error("Multiset must be nonempty")
    end
    clean!(M)
    return first(first(M.data))
end

include("iter.jl")
include("operations.jl")

import Base: keys, values, pairs

function keys(M::Multiset{T}) where {T}
    clean!(M)
    return keys(M.data)
end

function values(M::Multiset{T}) where {T}
    clean!(M)
    return values(M.data)
end

function pairs(M::Multiset{T}) where {T}
    clean!(M)
    return pairs(M.data)
end

end #end of Module


text:
__precompile__(true)
module Hive

using Thrift
using Markdown
using Dates
using Unicode

import Base: close, isready, show, eof, iterate, vcat

export HiveSession, HiveAuth, HiveAuthSASLPlain, Tabular
export eof, close, isready, result, cancel, get_info, show
export catalogs, schemas, tables, tabletypes, functions, columns, execute
export tabular, tabulars, records, columnchunks, columnchunk
export iterate, iteratorsize

export InfoType

# enable logging only during debugging
#using Logging
#const logger = Logging.configure(level=DEBUG)
##const logger = Logging.configure(filename="/tmp/hive$(getpid()).log", level=DEBUG)
#macro logmsg(s)
#    quote
#        debug($(esc(s)))
#    end
#end
macro logmsg(s)
end

# package code goes here
include("HS2/HS2.jl")
using .HS2
include("sess.jl")
include("types.jl")
include("resultset.jl")
include("util.jl")
include("api.jl")

end # module


text:
"""
Calculate the best-efficiency point (BEP) for the pumps in the water network. The BEP flow,
head, and power are added to the pumps' dictionaries.
"""
function calc_pump_bep!(data::Dict{String,<:Any})
    for (pump_id, pump) in data["pump"]
        q_bep = NaN

        # Calculate best efficiency point efficiencies.
        if haskey(pump, "efficiency_curve")
            eff_curve_tuples = pump["efficiency_curve"]
            eff_curve = array_from_tuples(eff_curve_tuples)

            # fit efficiency curve to determine q_bep and eta_bep
            A = hcat(eff_curve[:, 1].^2, eff_curve[:, 1])
            a, b = A \ eff_curve[:, 2]
            q_bep, eta_bep = -0.5*b * inv(a), -0.25*b^2 * inv(a)
        else
            # presume the single-value for efficiency is at the BEP
            eta_bep = pump["efficiency"]
        end

        # Calculate best efficiency point head gains.
        head_curve_tuples = pump["head_curve"]
        head_curve = array_from_tuples(head_curve_tuples)

        if isnan(q_bep)
            # an efficiency curve was not provided, and so we must determine q_bep from the
            # head curve
            if size(head_curve)[1] == 1
                q_bep, g_bep = head_curve[1, 1], head_curve[1, 2]
            else
                A = hcat(head_curve[:, 1].^2, ones(size(head_curve, 1)))
                c, d = A \ head_curve[:, 2]
                q_bep, g_bep = sqrt(-0.25 * d * inv(c)), 0.75 * d
            end
        else
            if size(head_curve)[1] == 1
                # probably rare that there is an efficiency curve and single-point or no head
                # curve -- nonetheless, it is a corner case that may be worth investigating
                # because q_bep would be overdetermined, and the two values may not agree
                #q_bep_alt = head_curve[1,1]
                g_bep = head_curve[1, 2]
            else
                # we could determine q_bep from the head curve and check wether it agrees
                # with the value determined from the efficiency curve -- will presume that
                # the one from the efficiency curve is more accurate
                A = -inv(3.0) * inv(q_bep^2) * head_curve[:, 1].^2 .+ (4.0 * inv(3.0))
                g_bep = A \ head_curve[:, 2]
            end
        end

        # Power at the BEP (1000 and 9.80665 are density and gravity).
        p_bep = _WM._DENSITY * _WM._GRAVITY * inv(eta_bep) * g_bep * q_bep

        # Add BEP values to the pump dictionary.
        pump["q_bep"] = q_bep
        pump["g_bep"] = g_bep
        pump["eta_bep"] = eta_bep
        pump["p_bep"] = p_bep
    end
end


text:
# preparing XRD data for gradient analysis
using Base.Threads

# 1. load data
using HDF5
path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
file = "Bi2O3_19F44_01_all_oned.h5"
datafile = h5open(path * file);
data = datafile["exp"];

stripenames = keys(data);
nstripes = length(stripenames)

temperatures = zeros(nstripes)
dwelltimes = zeros(nstripes)

nspec = length(names(data[stripenames[1]]))
nQ = size(read(data[stripenames[1]]["1"], "integrated_1d"), 1)
X = zeros(nQ, nspec, nstripes)
Y = similar(X)
for j in eachindex(stripenames)
    println("getting data for stripe $j")
    stripe = data[stripenames[j]]
    conditions = split(stripenames[j], "_")
    temperatures[j] = parse(Float64, conditions[4])
    dwelltimes[j] = parse(Float64, conditions[2])
    for i in 1:nspec
        d = stripe["$(i-1)"]
        XY = read(d, "integrated_1d")
        X[:, i, j], Y[:, i, j] = eachcol(XY)
    end
end

# take log10 of dwelltime
@. dwelltimes = log10(dwelltimes)

# 3. get centered position
using JSON
function get_positions(path, nspec, stripenames)
    nstripes = length(stripenames)
    posfile = JSON.parsefile(path * "xrd_centers.json")
    step = 10 # the positions are 10 μm appart
    pos = range(0, length = nspec, step = step) # μm
    pos = collect(pos) ./ 1e3 # mm
    positions = zeros(nspec, nstripes)
    for i in 1:nstripes
        ci = posfile[stripenames[i]] + 1 # add one because python indexing
        ci = convert(Int, ci)
        @. positions[:, i] = pos .- pos[ci]
    end
    return positions
end
positions = get_positions(path, nspec, stripenames)

# 4. interpolate onto same q grid
using Interpolations
function interpolate_xrd(X::AbstractArray, Y::AbstractArray)
    # q = range(extrema(Q)..., length = nQ)
    # q = range(10, 40, length = nQ) # cut off noisy parts
    p = .03 # cut off percentage of Q domain
    minQ = maximum(minimum.(eachcol(X[:, 1, :])))
    maxQ = minimum(maximum.(eachcol(X[:, 1, :])))
    q = range(minQ * (1+p), maxQ * (1-p), length = size(X, 1)) # interpolation grid
    Yq = zeros(length(q), size(Y, 2), size(Y, 3))
    @threads for j in 1:size(Y, 3)
        println("interpolating stripe $j")
        for i in 1:size(Y, 2)
            x, y = @views X[:, i, j], Y[:, i, j]
            dx = diff(x)
            all(x->isapprox(x, dx[1], atol = 1e-5), dx) || throw("q not uniformly spaced")
            x = range(extrema(x)..., length = length(x))
            f = CubicSplineInterpolation(x, y)
            @. Yq[:, i, j] = f.(q)
        end
    end
    return q, Yq
end

q, spectrograms = interpolate_xrd(X, Y)

# 5. exclude bad xrd stripes
using DelimitedFiles
function good_indices(data)
    i = string.(@view readdlm(path * "exclude_xrd.txt", '.')[:, 1])
    findall(!in(i), names(data))
end
ind = good_indices(data)

spectrograms = spectrograms[:, :, ind]
temperatures = temperatures[ind]
dwelltimes = dwelltimes[ind]
positions = positions[:, ind]

path = "SARA/NatCom2020/outer/data/"
f = h5open(path * "Bi2O3_19F44_01_outer_xrd_data.h5", "w")
f["q"] = collect(q)
f["spectrograms"] = spectrograms
f["dimensions"] = ["q", "position", "stripe"]
f["positions"] = positions
f["temperatures"] = temperatures
f["dwelltimes"] = dwelltimes
close(f)
close(datafile)

# path = "SARA/NatCom2020/outer/data/"
# f = h5open(path * "Bi2O3_19F44_01_outer_xrd_data.h5", "r")
# q = read(f, "q")
# spectrograms = read(f, "spectrograms")
# positions = read(f, "positions")
# temperatures = read(f, "temperatures")
# dwelltimes = read(f, "dwelltimes")
# close(f)
#
# using Plots
# plotly()
# i = 2
# heatmap(positions[:, i], q, @views log10.(spectrograms[:, :, i]))
# plot(positions[:, i], coefficients[:, :, i]')
# scatter(positions[:, i], coefficients[2, :, i])
# gui()
# temperatures[i]
# dwelltimes[i]


text:
using Test

include("welcome.jl")

@testset "1. Create the welcome message" begin
    @testset "Welcome message for customer with first letter capitalized" begin
        @test welcome("Formora") == "Welcome to the Tech Palace, FORMORA"
    end

    @testset "Welcome message for customer with only lowercase letters" begin
        @test welcome("oromis") == "Welcome to the Tech Palace, OROMIS"
    end

    @testset "Welcome message for customer with dash in name" begin
        @test welcome("Svit-kona") == "Welcome to the Tech Palace, SVIT-KONA"
    end

    @testset "Welcome message for customer with only uppercase letters" begin
        @test welcome("ARVA") == "Welcome to the Tech Palace, ARVA"
    end

    @testset "Welcome message for customer with non-latin letters" begin
        @test welcome("Andumë") == "Welcome to the Tech Palace, ANDUMË"
    end
end

@testset "2. Add a fancy border" begin
    @testset "Add border with 10 stars per line" begin
        @test add_border("Welcome!", 10) == "**********\nWelcome!\n**********"
    end

    @testset "Add border with 2 stars per line" begin
        @test add_border("Hi", 2) == "**\nHi\n**"
    end
end

@testset "3. Clean up old marketing messages" begin
    @testset "Cleanup message with leading whitespace" begin
        @test clean("                DISCOUNT") == "DISCOUNT"
    end

    @testset "Cleanup message with trailing whitespace" begin
        @test clean("SALE      ") == "SALE"
    end

    @testset "Cleanup message with leading and trailing whitespace" begin
        @test clean("    BUY NOW, SAVE 10%   ") == "BUY NOW, SAVE 10%"
    end
end

max_stars_repo_name: stepinski/julia

max_issues_repo_name: stepinski/julia

max_forks_repo_name: stepinski/julia


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Xpress, Statistics

function read_dat_file(dat_file_path)
    data = readdlm(dat_file_path, comments=true)
    SDP_var_list = []
    Sgen_var_list = []
    Bin_var_list = []
    λ = ""
    i = 1
    while data[i,1] == "VAR_TYPE"
        varname = data[i,3]
        if data[i,2] == "REAL"
            if varname[1:6] == "lambda"
                λ = varname
            elseif varname[1:3] == "Gen"
                push!(Sgen_var_list, varname)
            else
                push!(SDP_var_list, varname)
            end
        elseif data[i,2] == "BOOL"
            push!(Bin_var_list, varname)
        end
        i+=1
    end
    dict_quad_ctr = Dict{String, Dict{Tuple{String,String}, Float64}}()
    dict_linear_ctr = Dict{String, Dict{String,Float64}}()
    dict_bounds_ctr = Dict{String, Dict{String,Float64}}()
    dict_constants_ctr = Dict{String, Float64}()
    dict_Bin_ctr = Dict{String, Dict{String, Float64}}()
    dict_MONO = Dict{String,Array{Any}}()
    for k in i:size(data,1)
        type_elem = data[k,1]
        elem_name = data[k,2]
        var1 = data[k,3]
        var2 = data[k,4]
        value = data[k,5]
        if type_elem == "QUAD"
            if !haskey(dict_quad_ctr, elem_name)
                dict_quad_ctr[elem_name] = Dict((var1,var2)=>value)
            else
                dict_quad_ctr[elem_name][(var1,var2)] = value
            end
        elseif type_elem == "LIN"
            if !haskey(dict_linear_ctr, elem_name)
                dict_linear_ctr[elem_name] = Dict(var2 => value)
            else
                dict_linear_ctr[elem_name][var2] = value
            end
        elseif type_elem == "MONO"
            if !haskey(dict_Bin_ctr, elem_name)
                dict_Bin_ctr[elem_name] = Dict(var1 => value)
            else
                dict_Bin_ctr[elem_name][var1] = value
            end
        elseif type_elem == "LB" || type_elem == "UB"
            if !haskey(dict_bounds_ctr, elem_name)
                dict_bounds_ctr[elem_name] = Dict(type_elem => value)
            else
                dict_bounds_ctr[elem_name][type_elem] = value
            end
        elseif type_elem == "CONST"
            dict_constants_ctr[elem_name] = value
        elseif type_elem == "MONO_DEF"
            if !haskey(dict_MONO, elem_name)
                if value == 1
                    dict_MONO[elem_name] = [var1]
                elseif value == 2
                    dict_MONO[elem_name] = [(var1,var1)]
                end
            else
                if value == 1
                    push!(dict_MONO[elem_name], var1)
                elseif value == 2
                    push!(dict_MONO[elem_name], (var1,var1))
                end
            end
        end
    end
    # println(dict_bounds_ctr)
    return λ, Sgen_var_list, SDP_var_list, Bin_var_list, dict_quad_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO, dict_linear_ctr
end

function read_blocks(DATA_PATH, FORMULATION, INSTANCE_NAME)
    network_name = split(INSTANCE_NAME, '-')[1]
    if FORMULATION == "NO_BLOCKS"
        BLOCKS=readdlm(joinpath(DATA_PATH,"blocks_cholesky", network_name*"_sdp_blocks.txt"))
        VAR = BLOCKS[:,2]
        B1 = ["B1" for i=1:size(BLOCKS,1)]
        BLOCKS = [B1 VAR]
        CLIQUE_TREE = []
    else
        BLOCKS=readdlm(joinpath(DATA_PATH, "blocks_$FORMULATION",network_name*"_sdp_blocks.txt"))
        CLIQUE_TREE = readdlm(joinpath(DATA_PATH, "cliquetree_$FORMULATION",network_name*"_sdp_cliquetree.txt"))
    end
    cliques_dict = Dict{String,Array{String}}()
    for i in 1:size(BLOCKS,1)
        block = BLOCKS[i,1]
        var = BLOCKS[i,2]
            if !haskey(cliques_dict, block)
                cliques_dict[block] = [var]
            else
                push!(cliques_dict[block], var)
            end
    end
    return cliques_dict, CLIQUE_TREE
end

function construct_approximate_solution(mat_var, blocks_dict, SDP_var_list)
  #complex blocks
  # f = open("sol_C.csv", "w")
  eigenvectors_max = Dict{String,Array{Complex{Float64}}}()
  dict_C_blocks = Dict{String, Array{Complex{Float64}}}()
  blocks_dict_C = Dict{String, Array{String}}()
  for (block, list_var) in blocks_dict
    nb_var_R = length(list_var)
    nb_var_C = Int(nb_var_R/2)
    list_var_C = [var[1:end-3] for var in list_var if var[end-1:end]=="Re"]
    # println(list_var_C)
    W = Array{Complex{Float64}}(undef, nb_var_C, nb_var_C)
    for k in 1:nb_var_C
      var1_C =  list_var_C[k]
      var1_Re = var1_C*"_Re"
      var1_Im = var1_C*"_Im"
      for l in k:nb_var_C
        var2_C =  list_var_C[l]
        var2_Re = var2_C*"_Re"
        var2_Im = var2_C*"_Im"
        if k==l
          W[k,l] = JuMP.value(mat_var[(var1_Re,var2_Re)][block]) + JuMP.value(mat_var[(var1_Im,var2_Im)][block])
        else
          W[k,l] = JuMP.value(mat_var[(var1_Re,var2_Re)][block]) + JuMP.value(mat_var[(var1_Im,var2_Im)][block])+ im*(JuMP.value(mat_var[(var1_Im,var2_Re)][block])-JuMP.value(mat_var[(var1_Re,var2_Im)][block]))
          W[l,k] = JuMP.value(mat_var[(var1_Re,var2_Re)][block]) +  JuMP.value(mat_var[(var1_Im,var2_Im)][block])- im*(JuMP.value(mat_var[(var1_Im,var2_Re)][block])-JuMP.value(mat_var[(var1_Re,var2_Im)][block]))
        end
      end
    end
    dict_C_blocks[block] = W
    # println(ishermitian(W))
    # println(dict_C_blocks[block])
    X_block = dict_C_blocks[block]
    K = size(X_block, 1)
    ef = eigen(X_block)
    ef.values[K:K]       #largest k
    eigenvectors_max[block] = ef.vectors[:, K:K]
    # println(eigenvectors_max)
    # norm_eig = norm(ef.vectors[:, K:K])
    # println(norm_eig)
    for i=1:K
      var = list_var_C[i]
      value = eigenvectors_max[block][i]
      # write(f, "$block ; $var ; $value  \n")
    end
    blocks_dict_C[block] = list_var_C
  end
  list_blocks = [block for block in keys(blocks_dict_C)]
  NB_BLOCKS = length(list_blocks)
  sync_pb = Model(with_optimizer(Xpress.Optimizer))
  @variable(sync_pb, θ[1:NB_BLOCKS], lower_bound=0, upper_bound=2*pi)
  obj = @expression(sync_pb, 0.0*θ[1]^2)
  for i in 1:NB_BLOCKS
    block1 = list_blocks[i]
    list_var_1 = blocks_dict_C[block1]
    p1 = eigenvectors_max[block1]
    for j in (i+1):NB_BLOCKS
      block2 = list_blocks[j]
      list_var_2 = blocks_dict_C[block2]
      p2 = eigenvectors_max[block2]
      for k in 1:length(list_var_1)
        var1 = list_var_1[k]
        for l in 1:length(list_var_2)
          var2 = list_var_2[l]
          if var1 == var2
            add_to_expression!(obj, (angle(p1[k])+θ[i]-angle(p2[l])-θ[j])^2)
          end
        end
      end
    end
  end
  @objective(sync_pb, Min, obj)
  optimize!(sync_pb)
  findX = Model(with_optimizer(Xpress.Optimizer))
  X_Re = Dict{String, VariableRef}()
  X_Im = Dict{String, VariableRef}()
  list_var_C = [var[1:end-3] for var in SDP_var_list if var[end-1:end]=="Re"]
  # println(list_var_C)
  for var in list_var_C
    X_Re[var] = @variable(findX, base_name="X_Re_$var")
    X_Im[var] = @variable(findX, base_name="X_Im_$var")
  end

  obj = @expression(findX, 0.0*X_Re["VOLT_1"]^2)
  # f = open("sol_sync_pb.csv", "w")
  for i in 1:NB_BLOCKS
    block = list_blocks[i]
    θ_block = JuMP.value(θ[i])
    p = exp(im*θ_block).*eigenvectors_max[block]
    i = 1
    for var in blocks_dict_C[block]
      value = p[i]
      # write(f, "$block ; $var ; $(abs(eigenvectors_max[block][i])) ; $(angle(eigenvectors_max[block][i])) ; $(abs(value)) ; $(angle(value)) \n")
      add_to_expression!(obj, (X_Re[var]-real(value))^2+(X_Im[var]-imag(value))^2)

      i +=1
    end
  end
  @objective(findX, Min, obj)
  optimize!(findX)
  for var in list_var_C
    value = (JuMP.value(X_Re[var])) + im*(JuMP.value(X_Im[var]))
  end
  return X_Re, X_Im
end

function construct_SDP(blocks_dict, CLIQUE_TREE, Pinput_csv_file, flag, Sgen_var_list, SDP_var_list, Bin_var_list,
   dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO, dict_variables_to_fix, solution_file)

   lines = readlines(Pinput_csv_file)
   NB_BLOCKS = length(blocks_dict)
    #initialize
    m = Model(with_optimizer(Mosek.Optimizer))
    if solution_file == "" #SDP in B&B
        m = Model(with_optimizer(Mosek.Optimizer, LOG=0))
    end
    #variables
    λ_and_Sgen_jumpvar = Dict(varname => @variable(m, base_name="$varname") for varname in Sgen_var_list)
    if flag == "plus"
      λ = "lambda_plus"
      λ_and_Sgen_jumpvar[λ] = @variable(m, base_name="$λ", lower_bound = 0.5, upper_bound=1)
    elseif flag == "minus"
      λ = "lambda_minus"
      λ_and_Sgen_jumpvar[λ] = @variable(m, base_name="$λ", lower_bound = 0, upper_bound=0.5)
   else
     exit("Problem with flag")
   end
   constraints_ref=Dict{String,JuMP.ConstraintRef}()
   for line in lines
       temp = split(line, ';')
       gen_name = split(temp[1])[1]
       Pinput = parse(Float64,temp[2])
       empty = temp[3]
       Pmin = parse(Float64,temp[4])
       Pmax = parse(Float64,temp[5])
       if Pinput > Pmax
           Pinput = Pmax
       elseif Pinput < Pmin
           Pinput = Pmin
       end
       Pgen_var = λ_and_Sgen_jumpvar[gen_name]
       bus = split(gen_name, '_')[end-1]
       ctrname = "_$(bus)_Gen_Activepower_Re"
       if flag == "minus"
         λminus = λ_and_Sgen_jumpvar[λ]
         constraints_ref[ctrname] = @constraint(m, -Pgen_var+ (Pmin+2*(Pinput-Pmin)*λminus) == 0)
      elseif flag == "plus"
        λplus = λ_and_Sgen_jumpvar[λ]
        constraints_ref[ctrname] = @constraint(m, -Pgen_var + (2*Pinput-Pmax+2*(Pmax- Pinput)*λplus) ==0)
      end
   end

   jump_uvar = Dict{String, VariableRef}()
   jumpBinvar = Dict{String, VariableRef}()
    uvar = Dict{String, Tuple{VariableRef, JuMP.GenericAffExpr}}()
    for bin_var in Bin_var_list
      jumpBinvar[bin_var] = @variable(m, base_name = "ξ_$bin_var", lower_bound = 0)
      jump_uvar[bin_var] = @variable(m, base_name = "u_$bin_var", lower_bound = 0, upper_bound=1)
    end

    jumpX = Dict{String,Array{VariableRef,2}}()
    coeff_block = Dict{Tuple{String,String}, Set{String}}()
    mat_var = Dict{Tuple{String,String}, Dict{String,Any}}()

    for (block, var_list) in blocks_dict
      for var1 in var_list
        for var2 in var_list
          mat_var[(var1,var2)] = Dict{String,Any}()
          coeff_block[(var1,var2)] = Set{String}()
        end
      end
    end
    size_block = Dict{String, Int64}()
    i_block = 0
    for (block, var_list) in blocks_dict
      i_block += 1
      i_var1 = 1
      size_block[block] = length(var_list)
      jumpX[block] =  @variable(m,[1:size_block[block],1:size_block[block]], base_name= "X_$block", PSD)
      for var1 in var_list
        i_var2 = 1
        for var2 in var_list
          mat_var[(var1,var2)][block] = jumpX[block][i_var1, i_var2]
          coeff_block[(var1,var2)] = union(coeff_block[(var1,var2)],Set{String}([block]))
          i_var2+=1
        end
        i_var1+=1
      end
    end

    xp = Dict{String,JuMP.GenericAffExpr}()
    objctrnames = union(keys(dict_quad_ctr), keys(dict_Bin_ctr), keys(dict_linear_ctr))
    for objctrname in objctrnames
      xp[objctrname] = 0*jumpX["B1"][1,1]
    end

     for (objctrname, dict_ctr) in dict_quad_ctr
       for (vars, value) in dict_ctr
         var1 = vars[1]
         var2 = vars[2]
         nb_block = length(coeff_block[(var1,var2)])
         min_size = Inf
         min_var = first(mat_var[(var1,var2)])[2]
         for (block, var) in mat_var[(var1,var2)]
           if length(blocks_dict[block]) < min_size
             min_size = length(blocks_dict[block])
             min_var = var
           end
         end
         add_to_expression!(xp[objctrname], value*min_var)
       end
     end

     for (objctrname, dict_ctr) in dict_linear_ctr
       for (var, value) in dict_ctr
         add_to_expression!(xp[objctrname], value*λ_and_Sgen_jumpvar[var])
       end
     end

     for (objctrname, value) in dict_constants_ctr
       add_to_expression!(xp[objctrname], value)
     end

    for (objctrname, dict_ctr) in dict_Bin_ctr
      value_binVar =  Dict{String,Float64}()
      Vkk_binVar = Dict{String, GenericAffExpr}()
       for (mono, value) in dict_ctr
         vars = dict_MONO[mono]
         ξ = ""
         Vkk = ""
         name = ""
         abs2_Vkk = 0*jumpX["B1"][1,1]
           if typeof(vars[1]) == Tuple{String,String} || typeof(vars[1]) == Tuple{SubString{String},SubString{String}}
             var1 = vars[1][1]
             var2 = vars[1][2]
             bin_var = vars[2]
           elseif typeof(vars[1]) == String || typeof(vars[1]) == SubString{String}
              bin_var = vars[1]
              var1 = vars[2][1]
              var2 = vars[2][2]
           end
           value_binVar[bin_var] = value
           nb_block = length(coeff_block[(var1,var2)])
           min_size = Inf
           min_var = first(mat_var[(var1,var2)])[2]
           for (block, variable) in mat_var[(var1,var2)]
             if length(blocks_dict[block]) < min_size
               min_size = length(blocks_dict[block])
               min_var = variable
             end
           end
               if !haskey(Vkk_binVar, bin_var)
                 Vkk_binVar[bin_var] = 1* min_var
               else
               add_to_expression!(Vkk_binVar[bin_var], min_var)
              end
        end
        add_to_expression!(xp[objctrname], sum(value*jumpBinvar[bin_var] for (bin_var,value) in value_binVar))
        for (bin_var, abs2_Vkk) in Vkk_binVar
            ξ = jumpBinvar[bin_var]
            name = "ξ_$bin_var"
            name_VOLTM_ctr = "_$(bin_var[7:end])_Volt_VOLTM_Re"
            ub_abs2_Vkk = dict_bounds_ctr[name_VOLTM_ctr]["UB"]
            lb_abs2_Vkk = dict_bounds_ctr[name_VOLTM_ctr]["LB"]
            if haskey(dict_variables_to_fix, bin_var)
              value = dict_variables_to_fix[bin_var]
              if value == 0
                constraints_ref[name] = @constraint(m, ξ == 0)
                @constraint(m, jump_uvar[bin_var] == 0)
              elseif value ==1
                constraints_ref[name] = @constraint(m, ξ - abs2_Vkk == 0)
                @constraint(m, jump_uvar[bin_var] == 1)
              else
                #McCormick
                @constraint(m, ξ - ub_abs2_Vkk*jump_uvar[bin_var] <= 0)
                @constraint(m, ξ - lb_abs2_Vkk*(jump_uvar[bin_var]-1) - abs2_Vkk <= 0)
                @constraint(m, ξ - lb_abs2_Vkk*jump_uvar[bin_var] >= 0)
                @constraint(m, ξ - ub_abs2_Vkk*(jump_uvar[bin_var]-1) - abs2_Vkk >= 0)
              end
            end
            uvar[bin_var] = (ξ, abs2_Vkk)
          end
        end

    #add objective
    my_timer = @elapsed @objective(m, Min , xp["OBJ"])

    #add constraints
    for (objctrname, exp) in xp
      if objctrname != "OBJ"
        bounds = dict_bounds_ctr[objctrname]
        if haskey(bounds, "LB")
          lb = bounds["LB"]
        else
          lb = "NONE"
        end
        if haskey(bounds, "UB")
          ub = bounds["UB"]
        else
          ub = "NONE"
        end
        if ub == lb
          constraints_ref[objctrname] = @constraint(m, exp == ub)
        elseif lb != "NONE" && ub != "NONE"
          constraints_ref[objctrname] = @constraint(m, lb <= exp <= ub)
        elseif lb == "NONE" && ub != "NONE"
          constraints_ref[objctrname] = @constraint(m,  exp <= ub)
        elseif lb !="NONE" && ub == "NONE"
          constraints_ref[objctrname] = @constraint(m, lb <= exp)
        end
      end
    end
    #println("initial constraints added to model m")
    #constraints linking common terms in blocks
    nb_coupling_constraints=0
    if NB_BLOCKS > 1
      for i in 1:size(CLIQUE_TREE,1)
        B1 = CLIQUE_TREE[i,1]
        B2 = CLIQUE_TREE[i,2]
        vars_B1 = blocks_dict[B1]
        vars_B2 = blocks_dict[B2]
        common_vars = [ var for var in intersect(vars_B1,vars_B2)]
        for i in 1:length(common_vars)
          var1 = common_vars[i]
          for j in i:length(common_vars)
            var2 = common_vars[j]
            JuMPvar1 = mat_var[(var1,var2)][B1]
            JuMPvar2 = mat_var[(var1,var2)][B2]
            constraints_ref[*("cc$(B1)_$(B2)_$var1","_$var2")]=@constraint(m,JuMPvar1-JuMPvar2==0)
            nb_coupling_constraints+=1
          end
        end
      end
    end
    #print(m)
    optimize!(m)

    println("Objective value : ", JuMP.objective_value(m))
    println("Status : ", JuMP.termination_status(m))

    value_bins = Dict{String, Float64}()
    for (varname, tuple) in uvar
        ξ = tuple[1]
        abs2_Vkk = tuple[2]
        u = JuMP.value(ξ)/JuMP.value(abs2_Vkk)
        value_bins[varname] = u
    end
    value_SDP_var = Dict{String, Array{Float64,2}}()
    for (block, X) in jumpX
      value_SDP_var[block] = JuMP.value.(X)
    end

    if solution_file != ""
        X_Re, X_Im = construct_approximate_solution(mat_var, blocks_dict, SDP_var_list)
        isdir("Mosek_solutions") || mkpath("Mosek_solutions")
          f = open(joinpath("Mosek_solutions", solution_file), "w")
          write(f,"#Variable    Value \n")
          for (varname, tuple) in uvar
              ξ = tuple[1]
              abs2_Vkk = tuple[2]
              u = JuMP.value(ξ)/JuMP.value(abs2_Vkk)
              write(f, "$varname    $u \n")
          end
          for (var, val_Re) in X_Re
            value_Re = (JuMP.value(val_Re)) #value_Re = (JuMP.value(X_Re[var]))
            value_Im = (JuMP.value(X_Im[var]))
            write(f, "$(var*"_Re")    $value_Re   \n")
            write(f, "$(var*"_Im")    $value_Im   \n")
          end
          close(f)
      end

  return value_bins, value_SDP_var, JuMP.objective_value(m), primal_status(m)
end

function construct_SDP_refgen(blocks_dict, CLIQUE_TREE, Pinput_csv_file, refbus, Sgen_var_list, SDP_var_list, Bin_var_list,
   dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO, dict_variables_to_fix, solution_file)

   lines = readlines(Pinput_csv_file)
   NB_BLOCKS = length(blocks_dict)
    #initialize
    m = Model(with_optimizer(Mosek.Optimizer))
    if solution_file == "" #SDP in B&B
        m = Model(with_optimizer(Mosek.Optimizer, LOG=0))
    end
    #variables
    Sgen_jumpvar = Dict(varname => @variable(m, base_name="$varname") for varname in Sgen_var_list)

   constraints_ref=Dict{String,JuMP.ConstraintRef}()
   for line in lines
       temp = split(line, ';')
       gen_name = split(temp[1])[1]
       Pinput = parse(Float64,temp[2])
       empty = temp[3]
       Pmin = parse(Float64,temp[4])
       Pmax = parse(Float64,temp[5])
       if Pinput > Pmax
           Pinput = Pmax
       elseif Pinput < Pmin
           Pinput = Pmin
       end
       Pgen_var = Sgen_jumpvar[gen_name]
       bus = split(gen_name, '_')[end-1]
       ctrname = "_$(bus)_Gen_Activepower_Re"
       if gen_name == refbus
         constraints_ref[ctrname] = @constraint(m, Pmin <= Pgen_var <= Pmax)
      else
        constraints_ref[ctrname] = @constraint(m, Pgen_var == Pinput)
      end
   end

   jump_uvar = Dict{String, VariableRef}()
   jumpBinvar = Dict{String, VariableRef}()
    uvar = Dict{String, Tuple{VariableRef, JuMP.GenericAffExpr}}()
    for bin_var in Bin_var_list
      jumpBinvar[bin_var] = @variable(m, base_name = "ξ_$bin_var", lower_bound = 0)
      jump_uvar[bin_var] = @variable(m, base_name = "u_$bin_var", lower_bound = 0, upper_bound=1)
    end

    jumpX = Dict{String,Array{VariableRef,2}}()
    coeff_block = Dict{Tuple{String,String}, Set{String}}()
    mat_var = Dict{Tuple{String,String}, Dict{String,Any}}()

    for (block, var_list) in blocks_dict
      for var1 in var_list
        for var2 in var_list
          mat_var[(var1,var2)] = Dict{String,Any}()
          coeff_block[(var1,var2)] = Set{String}()
        end
      end
    end
    size_block = Dict{String, Int64}()
    i_block = 0
    for (block, var_list) in blocks_dict
      i_block += 1
      i_var1 = 1
      size_block[block] = length(var_list)
      jumpX[block] =  @variable(m,[1:size_block[block],1:size_block[block]], base_name= "X_$block", PSD)
      for var1 in var_list
        i_var2 = 1
        for var2 in var_list
          mat_var[(var1,var2)][block] = jumpX[block][i_var1, i_var2]
          coeff_block[(var1,var2)] = union(coeff_block[(var1,var2)],Set{String}([block]))
          i_var2+=1
        end
        i_var1+=1
      end
    end

    xp = Dict{String,JuMP.GenericAffExpr}()
    objctrnames = union(keys(dict_quad_ctr), keys(dict_Bin_ctr), keys(dict_linear_ctr))
    for objctrname in objctrnames
      xp[objctrname] = 0*jumpX["B1"][1,1]
    end

     for (objctrname, dict_ctr) in dict_quad_ctr
       for (vars, value) in dict_ctr
         var1 = vars[1]
         var2 = vars[2]
         nb_block = length(coeff_block[(var1,var2)])
         min_size = Inf
         min_var = first(mat_var[(var1,var2)])[2]
         for (block, var) in mat_var[(var1,var2)]
           if length(blocks_dict[block]) < min_size
             min_size = length(blocks_dict[block])
             min_var = var
           end
         end
         add_to_expression!(xp[objctrname], value*min_var)
       end
     end

     for (objctrname, dict_ctr) in dict_linear_ctr
       for (var, value) in dict_ctr
         add_to_expression!(xp[objctrname], value*Sgen_jumpvar[var])
       end
     end

     for (objctrname, value) in dict_constants_ctr
       add_to_expression!(xp[objctrname], value)
     end

    for (objctrname, dict_ctr) in dict_Bin_ctr
      value_binVar =  Dict{String,Float64}()
      Vkk_binVar = Dict{String, GenericAffExpr}()
       for (mono, value) in dict_ctr
         vars = dict_MONO[mono]
         ξ = ""
         Vkk = ""
         name = ""
         abs2_Vkk = 0*jumpX["B1"][1,1]
           if typeof(vars[1]) == Tuple{String,String} || typeof(vars[1]) == Tuple{SubString{String},SubString{String}}
             var1 = vars[1][1]
             var2 = vars[1][2]
             bin_var = vars[2]
           elseif typeof(vars[1]) == String || typeof(vars[1]) == SubString{String}
              bin_var = vars[1]
              var1 = vars[2][1]
              var2 = vars[2][2]
           end
           value_binVar[bin_var] = value
           nb_block = length(coeff_block[(var1,var2)])
           min_size = Inf
           min_var = first(mat_var[(var1,var2)])[2]
           for (block, variable) in mat_var[(var1,var2)]
             if length(blocks_dict[block]) < min_size
               min_size = length(blocks_dict[block])
               min_var = variable
             end
           end
               if !haskey(Vkk_binVar, bin_var)
                 Vkk_binVar[bin_var] = 1* min_var
               else
               add_to_expression!(Vkk_binVar[bin_var], min_var)
              end
        end
        add_to_expression!(xp[objctrname], sum(value*jumpBinvar[bin_var] for (bin_var,value) in value_binVar))
        for (bin_var, abs2_Vkk) in Vkk_binVar
            ξ = jumpBinvar[bin_var]
            name = "ξ_$bin_var"
            name_VOLTM_ctr = "_$(bin_var[7:end])_Volt_VOLTM_Re"
            ub_abs2_Vkk = dict_bounds_ctr[name_VOLTM_ctr]["UB"]
            lb_abs2_Vkk = dict_bounds_ctr[name_VOLTM_ctr]["LB"]
            if haskey(dict_variables_to_fix, bin_var)
              value = dict_variables_to_fix[bin_var]
              if value == 0
                constraints_ref[name] = @constraint(m, ξ == 0)
                @constraint(m, jump_uvar[bin_var] == 0)
              elseif value ==1
                constraints_ref[name] = @constraint(m, ξ - abs2_Vkk == 0)
                @constraint(m, jump_uvar[bin_var] == 1)
              else
                #McCormick
                @constraint(m, ξ - ub_abs2_Vkk*jump_uvar[bin_var] <= 0)
                @constraint(m, ξ - lb_abs2_Vkk*(jump_uvar[bin_var]-1) - abs2_Vkk <= 0)
                @constraint(m, ξ - lb_abs2_Vkk*jump_uvar[bin_var] >= 0)
                @constraint(m, ξ - ub_abs2_Vkk*(jump_uvar[bin_var]-1) - abs2_Vkk >= 0)
              end
            end
            uvar[bin_var] = (ξ, abs2_Vkk)
          end
        end

    #add objective
    my_timer = @elapsed @objective(m, Min , xp["OBJ"])

    #add constraints
    for (objctrname, exp) in xp
      if objctrname != "OBJ"
        bounds = dict_bounds_ctr[objctrname]
        if haskey(bounds, "LB")
          lb = bounds["LB"]
        else
          lb = "NONE"
        end
        if haskey(bounds, "UB")
          ub = bounds["UB"]
        else
          ub = "NONE"
        end
        if ub == lb
          constraints_ref[objctrname] = @constraint(m, exp == ub)
        elseif lb != "NONE" && ub != "NONE"
          constraints_ref[objctrname] = @constraint(m, lb <= exp <= ub)
        elseif lb == "NONE" && ub != "NONE"
          constraints_ref[objctrname] = @constraint(m,  exp <= ub)
        elseif lb !="NONE" && ub == "NONE"
          constraints_ref[objctrname] = @constraint(m, lb <= exp)
        end
      end
    end
    #println("initial constraints added to model m")
    #constraints linking common terms in blocks
    nb_coupling_constraints=0
    if NB_BLOCKS > 1
      for i in 1:size(CLIQUE_TREE,1)
        B1 = CLIQUE_TREE[i,1]
        B2 = CLIQUE_TREE[i,2]
        vars_B1 = blocks_dict[B1]
        vars_B2 = blocks_dict[B2]
        common_vars = [ var for var in intersect(vars_B1,vars_B2)]
        for i in 1:length(common_vars)
          var1 = common_vars[i]
          for j in i:length(common_vars)
            var2 = common_vars[j]
            JuMPvar1 = mat_var[(var1,var2)][B1]
            JuMPvar2 = mat_var[(var1,var2)][B2]
            constraints_ref[*("cc$(B1)_$(B2)_$var1","_$var2")]=@constraint(m,JuMPvar1-JuMPvar2==0)
            nb_coupling_constraints+=1
          end
        end
      end
    end
    # print(m)
    optimize!(m)

    println("Objective value : ", JuMP.objective_value(m))
    println("Status : ", JuMP.termination_status(m))

    value_bins = Dict{String, Float64}()
    for (varname, tuple) in uvar
        ξ = tuple[1]
        abs2_Vkk = tuple[2]
        u = JuMP.value(ξ)/JuMP.value(abs2_Vkk)
        value_bins[varname] = u
    end
    value_SDP_var = Dict{String, Array{Float64,2}}()
    for (block, X) in jumpX
      value_SDP_var[block] = JuMP.value.(X)
    end

    if solution_file != ""
        X_Re, X_Im = construct_approximate_solution(mat_var, blocks_dict, SDP_var_list)
        isdir("Mosek_solutions") || mkpath("Mosek_solutions")
          f = open(joinpath("Mosek_solutions", solution_file), "w")
          write(f,"#Variable    Value \n")
          for (varname, tuple) in uvar
              ξ = tuple[1]
              abs2_Vkk = tuple[2]
              u = JuMP.value(ξ)/JuMP.value(abs2_Vkk)
              write(f, "$varname    $u \n")
          end
          for (var, val_Re) in X_Re
            value_Re = (JuMP.value(val_Re)) #value_Re = (JuMP.value(X_Re[var]))
            value_Im = (JuMP.value(X_Im[var]))
            write(f, "$(var*"_Re")    $value_Re   \n")
            write(f, "$(var*"_Im")    $value_Im   \n")
          end
          close(f)
      end

  return value_bins, value_SDP_var, JuMP.objective_value(m), primal_status(m)

end

function solve_SDP(ROPF, flag)
    output_instance_path = ROPF.output_instance_path
    output_decomposition_path = ROPF.output_decomposition_path
    generation_files_path = ROPF.generation_files_path
    FORMULATION = ROPF.decomposition
    INSTANCE_NAME = ROPF.instance_name
    generation = ROPF.generation
    originalSTDOUT = stdout
    outpath = joinpath("Mosek_runs")
    isdir(outpath) || mkpath(outpath)
    outlog = open(joinpath(outpath,"$(INSTANCE_NAME)_$(generation)_$(flag)_$(FORMULATION).log"), "w")
    redirect_stdout(outlog)
    instance_dat_file_path = joinpath(output_instance_path, "$(INSTANCE_NAME).dat")
    Pinput_csv_file = joinpath(generation_files_path,"$(INSTANCE_NAME)_$(generation).csv")
    solution_file = joinpath("$(INSTANCE_NAME)_$(generation)_$(flag).dat")
    λ, Sgen_var_list, SDP_var_list, Bin_var_list, dict_quad_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr,
     dict_MONO, dict_linear_ctr = read_dat_file(instance_dat_file_path)
    cliques_dict, CLIQUE_TREE = read_blocks(output_decomposition_path, FORMULATION, INSTANCE_NAME)
     value_bins, value_SDP_var, obj, statut = construct_SDP(cliques_dict, CLIQUE_TREE, Pinput_csv_file, flag, Sgen_var_list, SDP_var_list, Bin_var_list,
          dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO,  Dict(bin => -1 for bin in Bin_var_list), solution_file)
    close(outlog)
    redirect_stdout(originalSTDOUT)
    return obj, statut
end

function solve_SDP_refgen(ROPF)
    output_instance_path = ROPF.output_instance_path
    output_decomposition_path = ROPF.output_decomposition_path
    generation_files_path = ROPF.generation_files_path
    FORMULATION = ROPF.decomposition
    INSTANCE_NAME = ROPF.instance_name
    generation = ROPF.generation
    originalSTDOUT = stdout
    outpath = joinpath("Mosek_runs")
    isdir(outpath) || mkpath(outpath)
    outlog = open(joinpath(outpath,"$(INSTANCE_NAME)_$(generation)_$(FORMULATION).log"), "w")
    redirect_stdout(outlog)
    instance_dat_file_path = joinpath(output_instance_path, "$(INSTANCE_NAME).dat")
    Pinput_csv_file = joinpath(generation_files_path,"$(INSTANCE_NAME)_$(generation).csv")
    solution_file = joinpath("$(INSTANCE_NAME)_$(generation).dat")
    refbus = readlines(joinpath(generation_files_path, "$(INSTANCE_NAME)_refbus.csv"))[1]
    λ, Sgen_var_list, SDP_var_list, Bin_var_list, dict_quad_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr,
     dict_MONO, dict_linear_ctr = read_dat_file(instance_dat_file_path)
    cliques_dict, CLIQUE_TREE = read_blocks(output_decomposition_path, FORMULATION, INSTANCE_NAME)
     value_bins, value_SDP_var, obj, statut = construct_SDP_refgen(cliques_dict, CLIQUE_TREE, Pinput_csv_file, refbus, Sgen_var_list, SDP_var_list, Bin_var_list,
          dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO,  Dict(bin => -1 for bin in Bin_var_list), solution_file)
    close(outlog)
    redirect_stdout(originalSTDOUT)
    return obj, statut
end


text:
using DataFrames
using CSV
using Plots

function main()

   theme(:ggplot2)
   data = CSV.read("../data/benchmark.csv", DataFrame, header = false)

   plot(log10.(data[!, 1]), color=:blue,
        title="Average run time",
        label="Sequantial",
        xaxis=collect(1:12),
        xlabel="Dataset size (log-10 scale)",
        ylabel="Average Time (log-10 scale)")
   plot!(log10.(data[!, 2]), color=:red, label="Parallel")
   savefig("../images/btime.png")

   plot(data[!, 1], color=:blue,
        title="Average run time",
        label="Sequantial",
        xaxis=collect(1:12),
        xlabel="Dataset size (log-10 scale)",
        ylabel="Average Time (s)")
        plot!(data[!, 2], color=:red, label="Parallel")
   savefig("../images/btime_nolog.png")

   speedup = data[!, 1] ./ data[!, 2]
   plot(speedup, color=:green,
        title="",
        label="",
        xlabel="Dataset size (log-10 scale)",
        ylabel="Speedup")
   savefig("../images/speedup.png")
end

main()

max_stars_repo_name: axmat/Parallel-Kmeans

max_issues_repo_name: axmat/ParallelKMeans

max_forks_repo_name: axmat/ParallelKMeans


text:
using VarStructs

# use syntax sugar :)
using VarStructs: U, UN # union and union with nothing

# declare
@var struct Address
    country::String
    city::UN{String} = nothing
    street::UN{String} = nothing
    unit::UN{ U{String, Int64} } = nothing
    zip::UN{ U{String, Int64} } = nothing
end

@var struct Human
    name::String
    address::UN{Address} = nothing
    id::UN{Int64} = nothing
end

@var struct Animal
    name::String
end

# get instance
h1 = Human(name = "Amin") # don't define all the things
h1.address = Address(country = "Canada") # add things later

h2 =  Human(name = "Someone", address = Address(country = "RandomCountry", street = "Street 1", zip = 1) )
h2.id = 2

a1 = Animal(name = "lion")

# redeclare!
@var struct Animal
    name::String
    id::UN{Int64} = nothing
end

a2 = Animal(name = "lion", id = 1)

# Dispatch!
function introduce_yourself(x::Human)
    println("My name is $(x.name), and I am from $(x.address.country)")
end

function introduce_yourself(x::Animal)
    println("My name is $(x.name), and I live in jungle")
end

introduce_yourself(h1)
introduce_yourself(h2)
introduce_yourself(a1)
introduce_yourself(a2)

# Custom constructor
function Human(n::String)
    @info "Address of human with name $n not found"
    return Human(name = n, address = Address(country = "Earth"))
end

h3 = Human("UnknownPerson")