text:
using Distances
using NPZ
using Iterators
using Images

verbose = true

array = convert(Array{Float64,2}, npzread("myo_footprints.npy")')

exp_schedule(t, tf) = (α0 =0.5; αf = 0.01; σ0 = 10.; σf = 1.; (α0 * exp(-t * log(αf/α0)/tf), σ0 * exp(-t * log(σf/σ0)/tf)))

gauss_transform(dmap, d, t, tf, α, σ) = α .* exp(-dmap.^2/σ^2)

tuppush(tup, value) = tuple(push!([i for i=tup], value)...)
tupunshift(tup, value) = tuple(unshift!([i for i=tup], value)...)

@everywhere function indmin_comp(a::Tuple, b::Tuple)
    da, ia = a
    db, ib = a
    mask = da .> db
    da[mask] = db[mask]
    ia[mask] = ib[mask]
    return (da, ia)
end

function findBMU(dist, vect, som, verbose=false)
    typeof(dist)<:PreMetric
    d = size(som)[1]
    mapsize = size(som)[2:end]
    vecsize = size(vect)[2]
    flatmapsize = reduce(*, mapsize)
    flatmap = reshape(som, (d, flatmapsize))
    pd = pairwise(dist, flatmap, vect)
    n = nprocs()
    if n == 1
        dmin, argmin = findmin(pd, 1)
    else
        prev = 0
        chunks = typeof(pd)
        # slower than single proc because of data copy...
        for i=map((x) -> convert(Int, round(x, 0)), linspace(1, flatmapsize, n+1)[2:end])
            #println(prev, i, size(chunks), size(pd))
            push!(chunks, pd[(prev+1):i, :])
            prev = i
        end
        @everywhere f(x) = findmin(x, 1)
        dmin, argmin = reduce(indmin_comp, pmap(f, chunks))
    end
    if vecsize == 1
        return (dmin, ind2sub(mapsize, vec(argmin)))
    else
        indflatmap, indvect = ind2sub((flatmapsize, vecsize), vec(argmin))
        indmap = ind2sub(mapsize, indflatmap)
        return (dmin, indmap)
    end
end

function mgrid(dims...)
    X = Array(Int, tuple(unshift!([i for i=dims], length(dims))...)) # Any[zeros(Int,dims...) for d in 1:length(dims)]
    for d in 1:length(dims)
        base_idx = Any[1:nd for nd in dims]
        for i in 1:dims[d]
            cur_idx = copy(base_idx)
            cur_idx[d] = i
            X[d, cur_idx...] = i
        end
    end
    #@show X
    return X
end

function SOM(input, mapsize, tf, dist, schedule, transform)
    typeof(dist)<:PreMetric
    D, N = size(input)
    true_mapsize = tupunshift(mapsize, D)
    m = minimum(input, 2)
    M = maximum(input, 2)
    som = m .+ M .* rand(true_mapsize)
    steps = [1:N]
    indices = mgrid(mapsize...)
    #get_mapdists(x::Array{Int,length(true_mapsize)}, y::Array{Int,length(true_mapsize)}) = sum(abs(x - y), 1)
    absdiff(x::Int, y::Int) = abs(x - y)
    const broadcast_absdiff = broadcast_function(absdiff)
    const broadcast_mult = broadcast_function(*)
    for t=1:tf
        if t % N == 1
            shuffle!(steps)
        end
        idx = steps[((t-1) % N) + 1]
        inpvec = input[:,[idx]]
        d, bmu = findBMU(dist, inpvec, som)
        bmu = reshape(hcat(bmu...), (bmusize = ones(Int, length(mapsize)+1); bmusize[1] = length(mapsize);tuple(bmusize...)))
        #println(" bmu ", size(bmu), " indices ", size(indices))
        mapdists = float(sum(broadcast_absdiff(indices, bmu), 1))
        α, σ = schedule(t, tf)
        Δmap = transform(mapdists, d, t, tf, α, σ)
        newsize = (shape=ones(Int, length(mapsize)+1); shape[1]=D; tuple(shape...))
        inpvec_bc = reshape(inpvec, newsize)
        som = som + broadcast_mult(inpvec_bc, Δmap)
        if verbose
            println(t,"/",tf, ": input #", idx, ", bmu: ", tuple(bmu...))
        end
    end
    return som
end

function neighbors(point, shape)
    f(x, s) = 1 <= x <= s
    return Iterators.Product([filter!((x)->f(x,s), [p-1, p, p+1]) for (p,s)=zip(point, shape)]...)
end

function umatrix(som, dist; toric=false)
    mapsize = size(som)[2:end]
    umat = zeros(mapsize)
    cartprod = Iterators.Product([[1:i] for i=mapsize]...)
    for point=cartprod
        neighb = neighbors(point, mapsize)
        c=0.
        pvec = som[:,point...]
        for n=neighb
            umat[point...] += evaluate(dist, pvec, som[:,n...])
            c += 1.
        end
        umat[point...] /= c
    end
    #indices = mgrid(mapsize...)
    return umat
end

som = SOM(array, (25,25), 500, CosineDist(), exp_schedule, gauss_transform)
umat = umatrix(som, CosineDist())
NPZ.npzwrite("umatrix.npy", umat)

bmus_distances, bmus = findBMU(CosineDist(), array, som, true)
density = zeros(size(som)[2:end])
#@show bmus
for b=zip(bmus...)
    @show b
    density[b...] += 1
end
NPZ.npzwrite("density.npy", density)

max_stars_repo_name: nathandesdouits/SelfOrganizingMaps

max_issues_repo_name: nathandesdouits/SelfOrganizingMaps

max_forks_repo_name: nathandesdouits/SelfOrganizingMaps


text:
# ---
# title: 1624. Largest Substring Between Two Equal Characters
# id: problem1624
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/largest-substring-between-two-equal-characters/description/>
# hidden: true
# ---
#
# Given a string `s`, return _the length of the longest substring between two
# equal characters, excluding the two characters._ If there is no such substring
# return `-1`.
#
# A **substring** is a contiguous sequence of characters within a string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aa"
#     Output: 0
#     Explanation: The optimal substring here is an empty substring between the two 'a's.
#
# **Example 2:**
#
#
#
#     Input: s = "abca"
#     Output: 2
#     Explanation: The optimal substring here is "bc".
#
#
# **Example 3:**
#
#
#
#     Input: s = "cbzxy"
#     Output: -1
#     Explanation: There are no characters that appear twice in s.
#
#
# **Example 4:**
#
#
#
#     Input: s = "cabbac"
#     Output: 4
#     Explanation: The optimal substring here is "abba". Other non-optimal substrings include "bb" and "".
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 300`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
export PhasePartition
# Thermodynamic states
export ThermodynamicState,
    PhaseDry,
    PhaseDry_ρT,
    PhaseDry_pT,
    PhaseDry_ρθ,
    PhaseDry_pθ,
    PhaseDry_ρp,
    PhaseEquil,
    PhaseEquil_ρTq,
    PhaseEquil_pTq,
    PhaseEquil_ρθq,
    PhaseEquil_pθq,
    PhaseEquil_ρpq,
    PhaseNonEquil,
    PhaseNonEquil_ρTq,
    PhaseNonEquil_ρθq,
    PhaseNonEquil_pθq,
    PhaseNonEquil_ρpq

"""
    ThermodynamicState{FT}

A thermodynamic state, which can be initialized for
various thermodynamic formulations (via its sub-types).
All `ThermodynamicState`'s have access to functions to
compute all other thermodynamic properties.
"""
abstract type ThermodynamicState{FT} end

"""
    PhasePartition

Represents the mass fractions of the moist air mixture.

# Constructors

    PhasePartition(q_tot::Real[, q_liq::Real[, q_ice::Real]])
    PhasePartition(ts::ThermodynamicState)

See also [`PhasePartition_equil`](@ref)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhasePartition{FT <: Real}
    "total specific humidity"
    tot::FT
    "liquid water specific humidity (default: `0`)"
    liq::FT
    "ice specific humidity (default: `0`)"
    ice::FT
    function PhasePartition(tot::FT, liq::FT, ice::FT) where {FT}
        q_tot_safe = max(tot, 0)
        q_liq_safe = max(liq, 0)
        q_ice_safe = max(ice, 0)
        return new{FT}(q_tot_safe, q_liq_safe, q_ice_safe)
    end
end

PhasePartition(q_tot::FT, q_liq::FT) where {FT <: Real} =
    PhasePartition(q_tot, q_liq, zero(FT))
PhasePartition(q_tot::FT) where {FT <: Real} =
    PhasePartition(q_tot, zero(FT), zero(FT))

#####
##### Dry states
#####

"""
    PhaseDry{FT} <: ThermodynamicState

A dry thermodynamic state (`q_tot = 0`).

# Constructors

    PhaseDry(param_set, e_int, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseDry{FT, PS} <: ThermodynamicState{FT}
    "parameter set, used to dispatch planet parameter function calls"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of dry air"
    ρ::FT
end
PhaseDry(param_set::APS, e_int::FT, ρ::FT) where {FT} =
    PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)

"""
    PhaseDry_pT(param_set, p, T)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `p` pressure
 - `T` temperature
"""
function PhaseDry_pT(param_set::APS, p::FT, T::FT) where {FT <: Real}
    e_int = internal_energy(param_set, T)
    ρ = air_density(param_set, T, p)
    return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end

"""
    PhaseDry_ρθ(param_set, ρ, θ_dry)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `θ_dry` dry potential temperature
"""
function PhaseDry_ρθ(param_set::APS, ρ::FT, θ_dry::FT) where {FT <: Real}
    T = air_temperature_given_θρq(param_set, θ_dry, ρ)
    e_int = internal_energy(param_set, T)
    return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end

"""
    PhaseDry_pθ(param_set, p, θ_dry)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `p` pressure
 - `θ_dry` dry potential temperature
"""
function PhaseDry_pθ(param_set::APS, p::FT, θ_dry::FT) where {FT <: Real}
    T = exner_given_pressure(param_set, p) * θ_dry
    e_int = internal_energy(param_set, T)
    ρ = air_density(param_set, T, p)
    return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end

"""
    PhaseDry_ρT(param_set, ρ, T)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `T` temperature
"""
function PhaseDry_ρT(param_set::APS, ρ::FT, T::FT) where {FT <: Real}
    e_int = internal_energy(param_set, T)
    return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end

"""
    PhaseDry_ρp(param_set, ρ, p)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `p` pressure
"""
function PhaseDry_ρp(param_set::APS, ρ::FT, p::FT) where {FT <: Real}
    T = air_temperature_from_ideal_gas_law(param_set, p, ρ)
    e_int = internal_energy(param_set, T)
    return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ)
end

#####
##### Equilibrium states
#####

"""
    PhaseEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic equilibrium (therefore, saturation adjustment
may be needed).

# Constructors

    PhaseEquil(param_set, e_int, ρ, q_tot)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseEquil{FT, PS} <: ThermodynamicState{FT}
    "parameter set, used to dispatch planet parameter function calls"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of air (potentially moist)"
    ρ::FT
    "total specific humidity"
    q_tot::FT
    "temperature: computed via [`saturation_adjustment`](@ref)"
    T::FT
end

"""
    PhaseEquil
Moist thermodynamic phase, given
 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `e_int` internal energy
 - `ρ` density
 - `q_tot` total specific humidity
and, optionally
 - `maxiter` maximum iterations for saturation adjustment
 - `temperature_tol` temperature tolerance for saturation adjustment
 - `sat_adjust` function pointer to particular saturation adjustment method, options include
    - `saturation_adjustment` uses Newtons method with analytic gradients
    - `saturation_adjustment_SecantMethod` uses Secant method
"""
function PhaseEquil(
    param_set::APS,
    e_int::FT,
    ρ::FT,
    q_tot::FT,
    maxiter::Int = 8,
    temperature_tol::FT = FT(1e-1),
    sat_adjust::Function = saturation_adjustment,
) where {FT <: Real}
    phase_type = PhaseEquil
    q_tot_safe = clamp(q_tot, FT(0), FT(1))
    T = sat_adjust(
        param_set,
        e_int,
        ρ,
        q_tot_safe,
        phase_type,
        maxiter,
        temperature_tol,
    )
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot_safe, T)
end

"""
    PhaseEquil_ρθq(param_set, ρ, θ_liq_ice, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` (moist-)air density
 - `θ_liq_ice` liquid-ice potential temperature
 - `q_tot` total specific humidity
 - `temperature_tol` temperature tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function PhaseEquil_ρθq(
    param_set::APS,
    ρ::FT,
    θ_liq_ice::FT,
    q_tot::FT,
    maxiter::Int = 36,
    temperature_tol::FT = FT(1e-1),
) where {FT <: Real}
    phase_type = PhaseEquil
    tol = ResidualTolerance(temperature_tol)
    T = saturation_adjustment_given_ρθq(
        param_set,
        ρ,
        θ_liq_ice,
        q_tot,
        phase_type,
        maxiter,
        tol,
    )
    q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
    e_int = internal_energy(param_set, T, q_pt)
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
end

"""
    PhaseEquil_pθq(param_set, p, θ_liq_ice, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `p` pressure
 - `θ_liq_ice` liquid-ice potential temperature
 - `q_tot` total specific humidity
 - `temperature_tol` temperature tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function PhaseEquil_pθq(
    param_set::APS,
    p::FT,
    θ_liq_ice::FT,
    q_tot::FT,
    maxiter::Int = 30,
    temperature_tol::FT = FT(1e-1),
) where {FT <: Real}
    phase_type = PhaseEquil
    tol = ResidualTolerance(temperature_tol)
    T = saturation_adjustment_given_pθq(
        param_set,
        p,
        θ_liq_ice,
        q_tot,
        phase_type,
        maxiter,
        tol,
    )
    ρ = air_density(param_set, T, p, PhasePartition(q_tot))
    q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
    e_int = internal_energy(param_set, T, q)
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q.tot, T)
end

"""
    PhaseEquil_ρTq(param_set, ρ, T, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `T` temperature
 - `q_tot` total specific humidity
"""
function PhaseEquil_ρTq(
    param_set::APS,
    ρ::FT,
    T::FT,
    q_tot::FT,
) where {FT <: Real}
    phase_type = PhaseEquil
    q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
    e_int = internal_energy(param_set, T, q)
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
end

"""
    PhaseEquil_pTq(param_set, p, T, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `p` pressure
 - `T` temperature
 - `q_tot` total specific humidity
"""
function PhaseEquil_pTq(
    param_set::APS,
    p::FT,
    T::FT,
    q_tot::FT,
) where {FT <: Real}
    phase_type = PhaseEquil
    ρ = air_density(param_set, T, p, PhasePartition(q_tot))
    q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
    e_int = internal_energy(param_set, T, q)
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
end

"""
    PhaseEquil_ρpq(param_set, ρ, p, q_tot, perform_sat_adjust=true)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `p` pressure
 - `q_tot` total specific humidity
 - `perform_sat_adjust` Bool indicating to perform saturation adjustment
"""
function PhaseEquil_ρpq(
    param_set::APS,
    ρ::FT,
    p::FT,
    q_tot::FT,
    perform_sat_adjust = false,
    maxiter::Int = 5,
) where {FT <: Real}

    phase_type = PhaseEquil
    if perform_sat_adjust
        T = saturation_adjustment_ρpq(
            param_set,
            ρ,
            p,
            q_tot,
            phase_type,
            maxiter,
        )
        q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
        e_int = internal_energy(param_set, T, q_pt)
    else
        q_pt = PhasePartition(q_tot)
        T = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt)
        e_int = internal_energy(param_set, T, q_pt)
    end
    return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T)
end

#####
##### Non-equilibrium states
#####

"""
   	 PhaseNonEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic non-equilibrium (therefore, temperature can
be computed directly).

# Constructors

    PhaseNonEquil(param_set, e_int, q::PhasePartition, ρ)

# Fields

$(DocStringExtensions.FIELDS)

"""
struct PhaseNonEquil{FT, PS} <: ThermodynamicState{FT}
    "parameter set, used to dispatch planet parameter function calls"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of air (potentially moist)"
    ρ::FT
    "phase partition"
    q::PhasePartition{FT}
end
function PhaseNonEquil(
    param_set::APS,
    e_int::FT,
    ρ::FT,
    q::PhasePartition{FT} = q_pt_0(FT),
) where {FT}
    return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q)
end

"""
    PhaseNonEquil_ρTq(param_set, ρ, T, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` (moist-)air density
 - `T` temperature
 - `q_pt` phase partition
"""
function PhaseNonEquil_ρTq(
    param_set::APS,
    ρ::FT,
    T::FT,
    q_pt::PhasePartition{FT},
) where {FT <: Real}
    e_int = internal_energy(param_set, T, q_pt)
    return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt)
end

"""
    PhaseNonEquil_ρθq(param_set, ρ, θ_liq_ice, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` (moist-)air density
 - `θ_liq_ice` liquid-ice potential temperature
 - `q_pt` phase partition
and, optionally
 - `potential_temperature_tol` potential temperature for non-linear equation solve
 - `maxiter` maximum iterations for non-linear equation solve
"""
function PhaseNonEquil_ρθq(
    param_set::APS,
    ρ::FT,
    θ_liq_ice::FT,
    q_pt::PhasePartition{FT},
    maxiter::Int = 10,
    potential_temperature_tol::FT = FT(1e-2),
) where {FT <: Real}
    phase_type = PhaseNonEquil
    tol = ResidualTolerance(potential_temperature_tol)
    T = air_temperature_given_θρq_nonlinear(
        param_set,
        θ_liq_ice,
        ρ,
        maxiter,
        tol,
        q_pt,
    )
    e_int = internal_energy(param_set, T, q_pt)
    return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt)
end

"""
    PhaseNonEquil_pθq(param_set, p, θ_liq_ice, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `p` pressure
 - `θ_liq_ice` liquid-ice potential temperature
 - `q_pt` phase partition
"""
function PhaseNonEquil_pθq(
    param_set::APS,
    p::FT,
    θ_liq_ice::FT,
    q_pt::PhasePartition{FT},
) where {FT <: Real}
    T = air_temperature_given_θpq(param_set, θ_liq_ice, p, q_pt)
    ρ = air_density(param_set, T, p, q_pt)
    e_int = internal_energy(param_set, T, q_pt)
    return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt)
end

"""
    PhaseNonEquil_ρpq(param_set, ρ, p, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
 - `ρ` density
 - `p` pressure
 - `q_pt` phase partition
"""
function PhaseNonEquil_ρpq(
    param_set::APS,
    ρ::FT,
    p::FT,
    q_pt::PhasePartition{FT},
) where {FT <: Real}
    T = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt)
    e_int = internal_energy(param_set, T, q_pt)
    return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt)
end


text:
export pdhcat, pdvcat, noconsolidation

Base.setindex!(df::DataFrame, v::AbstractArray, col::Union{Symbol, String}) = setindex!(df, v, [col])

function Base.setindex!(df::DataFrame, v::AbstractArray, cols)
    if isempty(df)
        df_merge = DataFrame(v, copy = false, columns = cols)
        df.pyo = df_merge.pyo
    elseif size(v, 1) == length(df) && !issubset(cols, df.columns)
        ndims(v) == 1 && (v = reshape(v, :, 1))
        df_set = DataFrame(v, copy = false, columns = cols, index = df.index)
        df_merge = df.merge(df_set, left_index = true, right_index = true, copy = false)
        df.pyo = df_merge.pyo
    else
        df[cols] = DataFrame(v)
    end
    return df
end

function pdhcat(dfs...)
    dfs = filter(!isempty, collect(dfs))
    df = dfs[1]
    for dfn in dfs[2:end]
        df = df.merge(dfn, left_index = true, right_index = true, copy = false)
    end
    return df
end

pdvcat(xs...) = pd.concat([xs...], axis = 0, ignore_index = true)

function noconsolidation()
    py"""
def _consolidate_inplace(self):
    pass

def _consolidate(self):
    return self.blocks

from pandas.core.internals import BlockManager
BlockManager._consolidate_inplace = _consolidate_inplace
BlockManager._consolidate = _consolidate
"""
end


text:
function order_updates!(updates)
    function _updt_idx(updt)
        idx = indices(updt.updt_coord)
        @assert length(idx) == 1
        idx[1]
    end
    sort!(updates, by=_updt_idx)
end

function fuse_coord(updates)
    indx = [indices(u.updt_coord)[1] for u in updates]
    @assert length(updates) == length(unique(indx))
    indx
end

function fuse_kernels(updates, cov_mat)
    d = size(cov_mat)[1]
    GaussianRandomWalkMix(2.38^2/d*cov_mat,
                          0.1^2/d*Matrix{Float64}(I, size(cov_mat)...),
                          0.05,
                          [u.t_kernel.pos for u in updates])
end

function fuse_priors(updates)
    prior_list = []
    coord_idx_list = []
    for updt in updates
        push!(prior_list, updt.priors.priors...)
        push!(coord_idx_list, updt.priors.coord_idx...)
    end
    Priors(Tuple(prior_list), Tuple(coord_idx_list))
end

function fuse_aux(updates, aux::Vector{<:UpdtAuxiliary}=[u.aux for u in updates])
    recompute_ODEs = any([a.recompute_ODEs for a in aux])
    UpdtAuxiliary(aux[1].solver, recompute_ODEs)
end


text:
# gas functions
function load_gas_data(set)
    """
    extracts, transforms and loads gas network data
    """
    prod_data = CSV.read("data/$(set[:gas_case])/gas_prod.csv", DataFrame; header=1)
    node_data = CSV.read("data/$(set[:gas_case])/gas_node.csv", DataFrame; header=1)
    pipe_data = CSV.read("data/$(set[:gas_case])/gas_edge.csv", DataFrame; header=1)
    load_fact = CSV.read("data/demand_factor.csv", DataFrame; header=1, limit=set[:T])

    # sets
    E  = collect(1:length(pipe_data[!,:w]))
    E_a = vcat(findall(x->x!=0, pipe_data[!,:b]))
    N  = collect(1:length(node_data[!,:demand]))
    # gas producer data
    c1 = Array(prod_data[!,:c1])
    c2 = Array(prod_data[!,:c2])
    c̀2 = sqrt.(c2)
    ϑ̅ = Array(prod_data[!,:p_max])
    ϑ̲ = Array(prod_data[!,:p_min])
    # node data
    δ = Array(node_data[!,:demand])
    ρ̅ = Array(node_data[!,:presh_max])
    ρ̲ = Array(node_data[!,:presh_min])
    # edge data
    n_s = Array(pipe_data[!,:n_s])
    n_r = Array(pipe_data[!,:n_r])
    w = Array(pipe_data[!,:w])
    s = Array(pipe_data[!,:K_h]) * .5
    b = Array(pipe_data[!,:b])
    κ̅ = Array(pipe_data[!,:kappa_max])
    κ̲ = Array(pipe_data[!,:kappa_min])
    # reference pressure node
    ref = 26
    # node-edge incidence matrix
    A = zeros(size(E,1),size(N,1))
    A_pl = zeros(size(E,1),size(N,1))
    A_mn = zeros(size(E,1),size(N,1))
    for l in E
        A[l,n_s[l]] = 1
        A_pl[l,n_s[l]] = 1
        A[l,n_r[l]] = -1
        A_mn[l,n_r[l]] = -1
    end
    # gas - pressure conversion matrix
    B = zeros(size(N,1),size(E,1))
    for j in E
        B[n_s[j],j] = b[j]
    end
    # initial linepack
    # ψ_o = Array(pipe_data[!,:psi_o])
    ψ_o = zeros(length(E))
    for i in E
        ψ_o[i] = s[i]*(node_data[n_s[i],:presh_init] + node_data[n_r[i],:presh_init]) * 0.05
    end
    # save gas system data
    net_data = Dict(:c1 => c1, :c2 => c2, :c̀2 => c̀2, :ϑ̅ => ϑ̅, :ϑ̲ => ϑ̲,
                    :δ => δ, :k_δ => load_fact[:,1], :ρ̅ => ρ̅, :ρ̲ => ρ̲,
                    :n_s => n_s, :n_r => n_r, :w => w, :ref => ref, :s => s./5,
                    :ψ_o => ψ_o, :κ̅ => κ̅, :κ̲ => κ̲, :b => b, :B => B,
                    :A => A', :A⁺ => A_pl', :A⁻ => A_mn',
                    :E => E, :E_a => E_a, :N => N)
    @info("done loading gas network data")
    return net_data
end

function non_convex_opt(net)
    """
    solves non-convex gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => 15000, "print_level" => 0))
    # variable declaration
    @variable(model, ϱ[1:length(net[:N]),t=1:set[:T]])
    @variable(model, φ[1:length(net[:E]),t=1:set[:T]])
    @variable(model, κ[1:length(net[:E]),t=1:set[:T]])
    @variable(model, ϑ[1:length(net[:N]),t=1:set[:T]])
    @variable(model, φ⁻[1:length(net[:E]),t=1:set[:T]])
    @variable(model, φ⁺[1:length(net[:E]),t=1:set[:T]])
    @variable(model, ψ[1:length(net[:E]),t=1:set[:T]])
    # minimize gas production cost
    @objective(model, Min, sum(sum(net[:c1][i]*ϑ[i,t] + net[:c2][i]*ϑ[i,t]^2 for i in net[:N]) for t in 1:set[:T]))
    # gas variable limits
    @constraint(model, inj_lim_max[i=net[:N],t=1:set[:T]], ϑ[i,t] <= net[:ϑ̅][i])
    @constraint(model, inj_lim_min[i=net[:N],t=1:set[:T]], ϑ[i,t] >= net[:ϑ̲][i])
    @constraint(model, pre_lim_max[i=net[:N],t=1:set[:T]], ϱ[i,t] <= net[:ρ̅][i])
    @constraint(model, pre_lim_min[i=net[:N],t=1:set[:T]], ϱ[i,t] >= net[:ρ̲][i])
    @constraint(model, com_lim_max[i=net[:E],t=1:set[:T]], κ[i,t] <= net[:κ̅][i])
    @constraint(model, com_lim_min[i=net[:E],t=1:set[:T]], κ[i,t] >= net[:κ̲][i])
    @constraint(model, SoC_lim_min[i=net[:E],t=set[:T]  ], ψ[i,t] >= net[:ψ_o][i])
    # gas flow equations
    @NLconstraint(model, w_eq[l=net[:E],t=1:set[:T]],  φ[l,t]*abs(φ[l,t]) - net[:w][l]^2 *((ϱ[ns(l),t] + κ[l,t])^2 - ϱ[nr(l),t]^2) == 0)
    @constraint(model, gas_bal[i=net[:N],t=1:set[:T]],
        sum(net[:A⁺][i,j]*φ⁺[j,t] + net[:A⁻][i,j]*φ⁻[j,t] for j in net[:E])
        == ϑ[i,t] - sum(net[:b][l]*κ[l,t] for l in net[:E] if ns(l) == i) - (Δ(t)*S(t)*ξ[:μ̂])[i])
    @constraint(model, φ_pl[l=net[:E_a],t=1:set[:T]],  φ[l,t] >= 0)
    @constraint(model, φ_avr[l=net[:E],t=1:set[:T]],   φ[l,t] == 1/2*(φ⁺[l,t] + φ⁻[l,t]))
    @constraint(model, ψ_avr[l=net[:E],t=1:set[:T]],   ψ[l,t] == 1/2*net[:s][l]*(ϱ[ns(l),t] + κ[l,t] + ϱ[nr(l),t]))
    @constraint(model, LP_b_t[l=net[:E],t=2:set[:T]],  ψ[l,t] == ψ[l,t-1] +  φ⁺[l,t] - φ⁻[l,t])
    @constraint(model, LP_b_1[l=net[:E],t=1],                   ψ[l,t] == net[:ψ_o][l] +  φ⁺[l,t] - φ⁻[l,t])
    @info("solving the non-convex gas optimization model takes time, wait ...")
    # solve model
    optimize!(model)
    @info("done solving the non-convex gas optimization model in $(solve_time(model)) sec. with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :ϑ => JuMP.value.(ϑ),
    :ϱ => JuMP.value.(ϱ),
    :φ => JuMP.value.(φ),
    :κ => JuMP.value.(κ),
    :ψ => JuMP.value.(ψ),
    :cost => JuMP.objective_value.(model),
    :status => termination_status(model),
    :CPUtime => solve_time(model),
    :model => model)
    return solution
end

function linearization(net,model)
    """
    extracts the senstivity coefficients from the non-convex Weymouth equation
    """
    con_nl = model[:w_eq]
    var_nl = [model[:ϱ];model[:φ];model[:κ]]
    # aux functions
    raw_index(var::MOI.VariableIndex) = var.value
    raw_index(con::NonlinearConstraintIndex) = con.value
    # extract variables
    variables = all_variables(model)
    # extract Jacobian structure
    d = NLPEvaluator(model)
    MOI.initialize(d, [:Jac])
    jac_str = MOI.jacobian_structure(d)
    # extract operating point
    opertaing_point = value.(variables)
    # evaluate Weymouth eq. Jacobian at the operating point
    V = zeros(length(jac_str))
    MOI.eval_constraint_jacobian(d, V, opertaing_point)
    # prepare the result
    I = [js[1] for js in jac_str] # rows, equations
    J = [js[2] for js in jac_str] # cols, variables
    jac = sparse(I, J, V)
    rows = raw_index.(index.(con_nl))
    cols = raw_index.(index.(var_nl))
    # return evaluated Jacobians for each time period (as dictionaries)
    Jacobian = Dict(:ϱ => [], :φ => [], :κ => [])
    function return_jac(net,jac,rows,cols)
        Jac = Matrix(jac[rows, cols])
        Jac_ϱ = Jac[net[:E],net[:N]]
        Jac_φ = Jac[net[:E],length(net[:N]) .+ net[:E]]
        Jac_κ = Jac[net[:E],(length(net[:N]) .+ length(net[:E])) .+ net[:E]]
        return Jac_ϱ, Jac_φ, Jac_κ
    end
    for t in 1:set[:T]
          J_ϱ, J_φ, J_κ  = return_jac(net,jac,rows[:,t],cols[:,t])
          push!(Jacobian[:ϱ],J_ϱ); push!(Jacobian[:φ],J_φ); push!(Jacobian[:κ],J_κ)
    end
    # return operating points for each time period (as dictionaries)
    point = Dict(:ϱ => [], :φ => [], :κ => [])
    for t in 1:set[:T]
          push!(point[:ϱ],opertaing_point[raw_index.(index.(var_nl))][net[:N],t])
          push!(point[:φ],opertaing_point[raw_index.(index.(var_nl))][length(net[:N]) .+ net[:E],t])
          push!(point[:κ],opertaing_point[raw_index.(index.(var_nl))][length(net[:N]) .+ length(net[:E]) .+ net[:E],t])
    end
    # compute linearization coefficients
    lin_res = Dict(:W_0 => [], :W_1 => [], :W_2 => [], :point => point)
    for t in 1:set[:T]
          push!(lin_res[:W_0], inv(Jacobian[:φ][t]) * (Jacobian[:ϱ][t] * point[:ϱ][t] + Jacobian[:κ][t] * point[:κ][t] + Jacobian[:φ][t] * point[:φ][t]) )
          push!(lin_res[:W_1], -inv(Jacobian[:φ][t]) * Jacobian[:ϱ][t])
          push!(lin_res[:W_2], -inv(Jacobian[:φ][t]) * Jacobian[:κ][t])
    end
    @info("done extracting sensetivities and reference pressures")
    return  lin_res
end

function linearized_opt(net,lin_res)
    """
    solves non-convex gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "max_iter" => 15000))
    # variable declaration
    @variable(model, ϱ[1:length(net[:N]),t=1:set[:T]])
    @variable(model, φ[1:length(net[:E]),t=1:set[:T]])
    @variable(model, κ[1:length(net[:E]),t=1:set[:T]])
    @variable(model, ϑ[1:length(net[:N]),t=1:set[:T]])
    @variable(model, φ⁻[1:length(net[:E]),t=1:set[:T]])
    @variable(model, φ⁺[1:length(net[:E]),t=1:set[:T]])
    @variable(model, ψ[1:length(net[:E]),t=1:set[:T]])
    # minimize gas production cost
    @objective(model, Min, sum(sum(net[:c1][i]*ϑ[i,t] + net[:c2][i]*ϑ[i,t]^2 for i in net[:N]) for t in 1:set[:T]))
    # gas variable limits
    @constraint(model, inj_lim_max[i=net[:N],t=1:set[:T]], ϑ[i,t] <= net[:ϑ̅][i])
    @constraint(model, inj_lim_min[i=net[:N],t=1:set[:T]], ϑ[i,t] >= net[:ϑ̲][i])
    @constraint(model, pre_lim_max[i=net[:N],t=1:set[:T]], ϱ[i,t] <= net[:ρ̅][i])
    @constraint(model, pre_lim_min[i=net[:N],t=1:set[:T]], ϱ[i,t] >= net[:ρ̲][i])
    @constraint(model, com_lim_max[i=net[:E],t=1:set[:T]], κ[i,t] <= net[:κ̅][i])
    @constraint(model, com_lim_min[i=net[:E],t=1:set[:T]], κ[i,t] >= net[:κ̲][i])
    @constraint(model, SoC_lim_min[i=net[:E],t=set[:T]  ], ψ[i,t] >= net[:ψ_o][i])
    # gas flow equations
    @constraint(model, w_eq[l=net[:E], t=1:set[:T]],   φ[l,t] .== lin_res[:W_0][t][l] .+ sum(lin_res[:W_1][t][l,n]*ϱ[n,t] for n in net[:N]) .+ sum(lin_res[:W_2][t][l,k]*κ[k,t] for k in net[:E]))
    # @constraint(model, gas_bal[i=net[:N],t=1:set[:T]], sum(net[:A⁺][i,j]*φ⁺[j,t] + net[:A⁻][i,j]*φ⁻[j,t] for j in net[:E]) == ϑ[i,t] - sum(net[:b][l]*κ[l,t] for l in net[:E] if ns(l) == i) - net[:δ][i]*net[:k_δ][t])
    @constraint(model, gas_bal[i=net[:N],t=1:set[:T]],
        sum(net[:A⁺][i,j]*φ⁺[j,t] + net[:A⁻][i,j]*φ⁻[j,t] for j in net[:E])
        == ϑ[i,t] - sum(net[:b][l]*κ[l,t] for l in net[:E] if ns(l) == i) - (Δ(t)*S(t)*ξ[:μ̂])[i])
    @constraint(model, φ_pl[l=net[:E_a],t=1:set[:T]],  φ[l,t] >= 0)
    @constraint(model, φ_avr[l=net[:E],t=1:set[:T]],   φ[l,t] == 1/2*(φ⁺[l,t] + φ⁻[l,t]))
    @constraint(model, ψ_avr[l=net[:E],t=1:set[:T]],   ψ[l,t] == 1/2*net[:s][l]*(ϱ[ns(l),t] + κ[l,t] + ϱ[nr(l),t]))
    @constraint(model, LP_b_t[l=net[:E],t=2:set[:T]],  ψ[l,t] == ψ[l,t-1] +  φ⁺[l,t] - φ⁻[l,t])
    @constraint(model, LP_b_1[l=net[:E],t=1],                   ψ[l,t] == net[:ψ_o][l] +  φ⁺[l,t] - φ⁻[l,t])
    @constraint(model, ref_ϱ[t=1:set[:T]],                  ϱ[net[:ref],t] == lin_res[:point][:ϱ][t][net[:ref]])
    # solve model
    optimize!(model)
    @info("done solving the linearized gas optimization model in $(solve_time(model)) sec. with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :ϑ => JuMP.value.(ϑ),
    :ϱ => JuMP.value.(ϱ),
    :φ => JuMP.value.(φ),
    :κ => JuMP.value.(κ),
    :ψ => JuMP.value.(ψ),
    :φ⁺ => JuMP.value.(φ⁺),
    :φ⁻ => JuMP.value.(φ⁻),
    :cost => JuMP.objective_value.(model),
    :status => termination_status(model),
    :model => model)
    return solution
end

function chance_constrained_gas_opt(net,lin_res,set)
    """
    solves the chance-constrained gas network optimization
    """
    # build model
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    # gas network variable declaration
    P =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    K =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Θ =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁻ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁺ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Ψ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    # auxilary variables for variability control
    s_ϱ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    # auxilary variables for DRO constraints
    y_ϑ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    x_ϑ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    y_ρ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    x_ρ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    y_κ =  [@variable(model, [1:length(net[:E])]) for t = 1:set[:T]]
    x_κ =  [@variable(model, [1:length(net[:E])]) for t = 1:set[:T]]
    # minimize gas injection cost + pressure variability measure
    @objective(model, Min, sum(net[:c1]' * Θ[t] * S(t) * ξ[:μ̂]
                                + tr(Θ[t]' * diagm(net[:c2]) * Θ[t] * S(t) * ξ[:Σ̂] * S(t)')
                           for t in 1:set[:T])
                           # + set[:α_ϱ]*sum(tr((P[t]*S(t)-P[t-1]*S(t-1))*ξ[:Σ̂]*(P[t]*S(t)-P[t-1]*S(t-1))') for t in 2:set[:T])
                           + set[:α_ϱ] * sum(sum(s_ϱ[t]) for t in 2:set[:T])
            )
    # variance constraints
    @constraint(model, pre_var[t=2:set[:T], i=net[:N]],
        [s_ϱ[t][i];(cholesky(ξ[:Σ̂]).U .+ zeros(set[:k_t][set[:T]],set[:k_t][set[:T]]))*(P[t]*S(t) - P[t-1]*S(t-1))[i,:]] in SecondOrderCone())
    # injection deviation constraints
    @constraint(model, inj_lim[t=1:set[:T], i=net[:N]],
        [set[:α_ϑ]*(Θ[t]*S(t)*ξ[:μ̂])[i];
            (cholesky(ξ[:Σ̂]).U .+ zeros(set[:k_t][set[:T]],set[:k_t][set[:T]])) * (Θ[t]*S(t) - Θ[t]*S(t)*ξ[:μ̂]*S(1))[i,:]]
                in SecondOrderCone())
    # linepack flexibility limit constraints
    @constraint(model, lp_lim[t=1:set[:T], i=net[:E]],
        [set[:α_ψ]*(Ψ[t]*S(t)*ξ[:μ̂])[i];
            (cholesky(ξ[:Σ̂]).U .+ zeros(set[:k_t][set[:T]],set[:k_t][set[:T]])) * (Ψ[t]*S(t) - Ψ[t]*S(t)*ξ[:μ̂]*S(1))[i,:]]
                in SecondOrderCone())
    # equivalents of stochastic equations
    @constraint(model, gas_bal[t=1:set[:T]], (net[:A⁺]*Φ⁺[t]*S(t) + net[:A⁻]*Φ⁻[t]*S(t) - Θ[t]*S(t) + net[:B]*K[t]*S(t) + Δ(t)*S(t)) .== 0)
    @constraint(model, lin_wey[t=1:set[:T]], (Φ[t]*S(t) - [lin_res[:W_0][t] zeros(length(net[:E]),set[:k_t][t]-1)]*S(t) - lin_res[:W_1][t]*P[t]*S(t) - lin_res[:W_2][t]*K[t]*S(t)) .== 0)
    @constraint(model, ref_pre[t=1:set[:T]], (P[t][net[:ref],:]' - [lin_res[:point][:ϱ][t][net[:ref]] zeros(1,set[:k_t][t]-1)])*S(t) .== 0)
    @constraint(model, av_flow[t=1:set[:T]], (Φ[t]*S(t) - 1/2*Φ⁺[t]*S(t) - 1/2*Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_def[t=1:set[:T]], (Ψ[t]*S(t) - 1/2*diagm(net[:s])*(K[t] + abs.(net[:A])'*P[t])*S(t)) .== 0)
    @constraint(model, l_p_bal_1[t=1], (Ψ[t]*S(t) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_bal_t[t=2:set[:T]], (Ψ[t]*S(t) - Ψ[t-1]*S(t-1) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    # reformulated chance constraints
    # double-sided gas injection chance constraints
    @constraint(model, ϑ_1[t=1:set[:T], i=1:length(net[:N])],
        [sqrt(set[:ε̅_g]) * (1/2*(net[:ϑ̅][i] - net[:ϑ̲][i]) - x_ϑ[t][i]);vcat(ξ[:Σ̅][2:end,2:end] * (Θ[t]*S(t))[i,2:end],y_ϑ[t][i])] in SecondOrderCone())
    @constraint(model, ϑ_2[t=1:set[:T], i=1:length(net[:N])],
        (Θ[t]*S(t))[i,1] - 1/2*(net[:ϑ̅][i] + net[:ϑ̲][i]) <= y_ϑ[t][i] + x_ϑ[t][i])
    @constraint(model, ϑ_3[t=1:set[:T], i=1:length(net[:N])],
        (Θ[t]*S(t))[i,1] - 1/2*(net[:ϑ̅][i] + net[:ϑ̲][i]) >= - y_ϑ[t][i] - x_ϑ[t][i])
    @constraint(model, ϑ_4[t=1:set[:T], i=1:length(net[:N])],
        1/2*(net[:ϑ̅][i] - net[:ϑ̲][i]) >= x_ϑ[t][i] >= 0)
    @constraint(model, ϑ_5[t=1:set[:T], i=1:length(net[:N])],
        y_ϑ[t][i] >= 0)
    # double-sided pressure chance constraints
    @constraint(model, ρ_1[t=1:set[:T], i=1:length(net[:N])],
        [sqrt(set[:ε̅_g]) * (1/2*(net[:ρ̅][i] - net[:ρ̲][i]) - x_ρ[t][i]);vcat(ξ[:Σ̅][2:end,2:end] * (P[t]*S(t))[i,2:end],y_ρ[t][i])] in SecondOrderCone())
    @constraint(model, ρ_2[t=1:set[:T], i=1:length(net[:N])],
        (P[t]*S(t))[i,1] - 1/2*(net[:ρ̅][i] + net[:ρ̲][i]) <= y_ρ[t][i] + x_ρ[t][i])
    @constraint(model, ρ_3[t=1:set[:T], i=1:length(net[:N])],
        (P[t]*S(t))[i,1] - 1/2*(net[:ρ̅][i] + net[:ρ̲][i]) >= - y_ρ[t][i] - x_ρ[t][i])
    @constraint(model, ρ_4[t=1:set[:T], i=1:length(net[:N])],
        1/2*(net[:ρ̅][i] - net[:ρ̲][i]) >= x_ρ[t][i] >= 0)
    @constraint(model, ρ_5[t=1:set[:T], i=1:length(net[:N])],
        y_ρ[t][i] >= 0)
    # double-sided pressure regulation chance constraints
    @constraint(model, κ_1[t=1:set[:T], i=1:length(net[:E])],
        [sqrt(set[:ε̅_g]) * (1/2*(net[:κ̅][i] - net[:κ̲][i]) - x_κ[t][i]);vcat(ξ[:Σ̅][2:end,2:end] * (K[t]*S(t))[i,2:end],y_κ[t][i])] in SecondOrderCone())
    @constraint(model, κ_2[t=1:set[:T], i=1:length(net[:E])],
        (K[t]*S(t))[i,1] - 1/2*(net[:κ̅][i] + net[:κ̲][i]) <= y_κ[t][i] + x_κ[t][i])
    @constraint(model, κ_3[t=1:set[:T], i=1:length(net[:E])],
        (K[t]*S(t))[i,1] - 1/2*(net[:κ̅][i] + net[:κ̲][i]) >= - y_κ[t][i] - x_κ[t][i])
    @constraint(model, κ_4[t=1:set[:T], i=1:length(net[:E])],
        1/2*(net[:κ̅][i] - net[:κ̲][i]) >= x_κ[t][i] >= 0)
    @constraint(model, κ_5[t=1:set[:T], i=1:length(net[:E])],
        y_κ[t][i] >= 0)
    # chance constraints on linepack and flow variables
    @constraint(model, ψ_T[t=set[:T], i=1:length(net[:E])],
        [(Ψ[t]*S(t)*ξ[:μ̂] - net[:ψ_o])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Ψ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, φ_pl[t=1:set[:T], i=net[:E_a]],
        [(Φ[t]*S(t)*ξ[:μ̂] .- 0)[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Φ[t]*S(t))[i,:]]
            in SecondOrderCone())
    # zero recourse for nodes and edges with no controllable components
    @constraint(model, zero_inj[t=1:set[:T], i=findall(x->x==0, net[:ϑ̅])], Θ[t][i,2:end] .== 0)
    @constraint(model, zero_reg[t=1:set[:T], i=setdiff(findall(x->x==0, net[:κ̅]),findall(x->x<0, net[:κ̲]))],K[t][i,2:end] .== 0)
    # solve model
    optimize!(model)
        @info("done solving chance-constrained gas network optimization in $(solve_time(model)) sec. with status $(termination_status(model))")
    # return solution
    solution = Dict(
    :P => [JuMP.value.(P[t]) for t in 1:set[:T]],
    :Φ => [JuMP.value.(Φ[t]) for t in 1:set[:T]],
    :K => [JuMP.value.(K[t]) for t in 1:set[:T]],
    :Θ => [JuMP.value.(Θ[t]) for t in 1:set[:T]],
    :Φ⁻ => [JuMP.value.(Φ⁻[t]) for t in 1:set[:T]],
    :Φ⁺ => [JuMP.value.(Φ⁺[t]) for t in 1:set[:T]],
    :Ψ => [JuMP.value.(Ψ[t]) for t in 1:set[:T]],
    :cost => sum(net[:c1]' * JuMP.value.(Θ[t]) * S(t) * ξ[:μ̂] + tr(JuMP.value.(Θ[t])' * diagm(net[:c2]) * JuMP.value.(Θ[t]) * S(t) * ξ[:Σ̂] * S(t)') for t in 1:set[:T]),
    :ϱ_std_tot => sum(sum(JuMP.value.(s_ϱ[t])) for t in 2:1:set[:T]) ,
    :status => termination_status(model),
    :obj => JuMP.objective_value(model),
    :model => model)
    return solution
end

function export_results(set,ξ,sol_lin,sol_sto_wo_var,sol_sto_wt_var)
    """
    exports results
    """
    ## wind data
    wind_results = zeros(1+set[:Ns],set[:T])
    wind_results[1,:] =  [sum(W(t)*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    for i in 1:set[:Ns]
        wind_results[i+1,:] = [sum(W(t)*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
    end
    ## gas consumption data
    gas_coms_results = zeros(1+set[:Ns],set[:T])
    gas_coms_results[1,:] =  [sum(Δ(t)*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    for i in 1:set[:Ns]
        gas_coms_results[i+1,:] = [sum(Δ(t)*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
    end
    ## first-stage linepack difference between det and stoch solutions
    δ_Ψ = (sol_sto_wo_var[:Ψ][1] .- sol_lin[:ψ][:,1])./sol_lin[:ψ][:,1] * 100
    δ_Ψ_n = zeros(length(δ_Ψ)); α = 10;
    for i in 1:51 # normalization
        δ_Ψ_n[i] = 1 - 1/2 * (1 - δ_Ψ[i] / sqrt(δ_Ψ[i]^2 + α))
    end
    ## pressure results
    det_pressure = [sum(sol_lin[:ϱ][:,t]) for t in 1:set[:T]]
    sto_pressure_wo_var = zeros(1+set[:Ns],set[:T])
    sto_pressure_wt_var = zeros(1+set[:Ns],set[:T])
    sto_pressure_wo_var[1,:] = [sum(sol_sto_wo_var[:P][t]*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    sto_pressure_wt_var[1,:] = [sum(sol_sto_wt_var[:P][t]*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    for i in 1:set[:Ns]
        sto_pressure_wo_var[i+1,:] = [sum(sol_sto_wo_var[:P][t]*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
        sto_pressure_wt_var[i+1,:] = [sum(sol_sto_wt_var[:P][t]*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
    end
    ## linepack results
    det_linepack = [sum(sol_lin[:ψ][:,t]) for t in 1:set[:T]]
    sto_linepack_wo_var = zeros(1+set[:Ns],set[:T])
    sto_linepack_wt_var = zeros(1+set[:Ns],set[:T])
    sto_linepack_wo_var[1,:] = [sum(sol_sto_wo_var[:Ψ][t]*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    sto_linepack_wt_var[1,:] = [sum(sol_sto_wt_var[:Ψ][t]*S(t)*ξ[:μ̂]) for t in 1:set[:T]]
    for i in 1:set[:Ns]
        sto_linepack_wo_var[i+1,:] = [sum(sol_sto_wo_var[:Ψ][t]*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
        sto_linepack_wt_var[i+1,:] = [sum(sol_sto_wt_var[:Ψ][t]*S(t)*ξ[:ξ̂][:,i]) for t in 1:set[:T]]
    end
    ## presusre regulation
    det_pres_reg = [sum(abs.(sol_lin[:κ][:,t])) for t in 1:set[:T]]
    sto_pres_reg_wo_var = zeros(1+set[:Ns],set[:T])
    sto_pres_reg_wt_var = zeros(1+set[:Ns],set[:T])
    sto_pres_reg_wo_var[1,:] = [sum(abs.(sol_sto_wo_var[:K][t]*S(t)*ξ[:μ̂])) for t in 1:set[:T]]
    sto_pres_reg_wt_var[1,:] = [sum(abs.(sol_sto_wt_var[:K][t]*S(t)*ξ[:μ̂])) for t in 1:set[:T]]
    for i in 1:set[:Ns]
        sto_pres_reg_wo_var[i+1,:] = [sum(abs.(sol_sto_wo_var[:K][t]*S(t)*ξ[:ξ̂][:,i])) for t in 1:set[:T]]
        sto_pres_reg_wt_var[i+1,:] = [sum(abs.(sol_sto_wt_var[:K][t]*S(t)*ξ[:ξ̂][:,i])) for t in 1:set[:T]]
    end

    CSV.write("results/wind_results.csv",  Tables.table([[i for i in 1:set[:T]]' ; wind_results./maximum(wind_results)]'), writeheader=false)
    CSV.write("results/gas_coms_results.csv",  Tables.table([[i for i in 1:set[:T]]' ; gas_coms_results./maximum(gas_coms_results)]'), writeheader=false)
    CSV.write("results/linepack_diff.csv",  Tables.table(δ_Ψ_n'), writeheader=false)

    base_value = max(maximum(det_pressure),maximum(sto_pressure_wo_var),maximum(sto_pressure_wt_var))
    CSV.write("results/det_pressure.csv",  Tables.table([[i for i in 1:set[:T]]' ; det_pressure'./base_value]'), writeheader=false)
    CSV.write("results/sto_pressure_wo_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_pressure_wo_var./base_value]'), writeheader=false)
    CSV.write("results/sto_pressure_wt_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_pressure_wt_var./base_value]'), writeheader=false)

    base_value = max(maximum(det_linepack),maximum(sto_linepack_wo_var),maximum(sto_linepack_wt_var))
    CSV.write("results/det_linepack.csv",  Tables.table([[i for i in 1:set[:T]]' ; det_linepack'./base_value]'), writeheader=false)
    CSV.write("results/sto_linepack_wo_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_linepack_wo_var./base_value]'), writeheader=false)
    CSV.write("results/sto_linepack_wt_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_linepack_wt_var./base_value]'), writeheader=false)

    base_value = max(maximum(det_pres_reg),maximum(sto_pres_reg_wo_var),maximum(sto_pres_reg_wt_var))
    CSV.write("results/det_pres_reg.csv",  Tables.table([[i for i in 1:set[:T]]' ; det_pres_reg'./base_value]'), writeheader=false)
    CSV.write("results/sto_pres_reg_wo_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_pres_reg_wo_var./base_value]'), writeheader=false)
    CSV.write("results/sto_pres_reg_wt_var.csv",  Tables.table([[i for i in 1:set[:T]]' ; sto_pres_reg_wt_var./base_value]'), writeheader=false)
end

function topology_conf(set)
    """
    Identifies sensetivity coefficients and reference pressure
    for different gas network topologies
    """
    # number of valves
    v_num = length(set[:b_valve_edge])
    # number of topology configurations
    C = 2^v_num
    # binary configuration list
    con_list = all_config(v_num)
    # compute senstivities and ref. pressure for each configuration
    output = Dict()
    @info("star enumerating candidate topologies")
    for i in 1:C
        # modify network topology
        net_data = load_gas_data(set)
        for j in 1:v_num
            net_data[:w][set[:b_valve_edge][j]] = net_data[:w][set[:b_valve_edge][j]] * con_list[i,j]
        end
        # run linearization
        sol_non_convex  = non_convex_opt(net_data)
        lin_res         = linearization(net_data,sol_non_convex[:model])
        # save results
        conf_data = Dict(:W_0 => [lin_res[:W_0][t] for t in 1:set[:T]],
                         :W_1 => [lin_res[:W_1][t] for t in 1:set[:T]],
                         :W_2 => [lin_res[:W_2][t] for t in 1:set[:T]],
                         :point => Dict(:ϱ => [lin_res[:point][:ϱ][t] for t in 1:set[:T]])
        )
        push!(output, i => conf_data)
    end
    @info("done enumerating candidate topologies")
    return output
end

function all_config(n)
    """
    Finds the list of possible network topologies
    """
    model = Model(optimizer_with_attributes(Mosek.Optimizer))
    @variable(model, ξ[1:n])
    @constraint(model, con[i=1:n], 0 <= ξ[i] <= 1)
    # create a polyhedron from the sample set
    poly = polyhedron(model, CDDLib.Library(:exact))
    # obtain V-representation of the polyhedron
    vr = vrep(poly)
    vr = MixedMatVRep(vr)
    p = vr.V
    return float(p)
end

function out_of_sample(solution,net,set,ξ)
    """
    Performs out-of-sample analysis
    """
    ϑ = zeros(length(net[:N]),set[:T],set[:Ns])
    ϱ = zeros(length(net[:N]),set[:T],set[:Ns])
    ψ = zeros(length(net[:E]),set[:T],set[:Ns])
    ϕ = zeros(length(net[:E]),set[:T],set[:Ns])
    κ = zeros(length(net[:E]),set[:T],set[:Ns])

    ϑ_mag =  zeros(set[:T])
    ϱ_mag =  zeros(set[:T])
    ψ_mag =  zeros(set[:T])
    ϕ_mag =  zeros(set[:T])
    κ_mag =  zeros(set[:T])

    ϑ_cvar =  zeros(set[:T],set[:Ns])
    ϱ_cvar =  zeros(set[:T],set[:Ns])
    ψ_cvar =  zeros(set[:T],set[:Ns])
    ϕ_cvar =  zeros(set[:T],set[:Ns])
    κ_cvar =  zeros(set[:T],set[:Ns])
    for t in 1:set[:T], s in 1:set[:Ns]
         ϑ[:,t,s] = solution[:Θ][t]*S(t)*ξ[:ξ̂][:,s]
         ϱ[:,t,s] = solution[:P][t]*S(t)*ξ[:ξ̂][:,s]
         ψ[:,t,s] = solution[:Ψ][t]*S(t)*ξ[:ξ̂][:,s]
         ϕ[:,t,s] = solution[:Φ][t]*S(t)*ξ[:ξ̂][:,s]
         κ[:,t,s] = solution[:K][t]*S(t)*ξ[:ξ̂][:,s]
    end
    for t in 1:set[:T], i in 1:48, s in 1:set[:Ns]
        (ϑ[i,t,s] - net[:ϑ̅][i]) >= 0 ? ϑ_mag[t] += (ϑ[i,t,s] - net[:ϑ̅][i]) : NaN
        (net[:ϑ̲][i] - ϑ[i,t,s]) >= 0 ? ϑ_mag[t] += (net[:ϑ̲][i] - ϑ[i,t,s]) : NaN
        (ϑ[i,t,s] - net[:ϑ̅][i]) >= 0 ? ϑ_cvar[t,s] += (ϑ[i,t,s] - net[:ϑ̅][i]) : NaN
        (net[:ϑ̲][i] - ϑ[i,t,s]) >= 0 ? ϑ_cvar[t,s] += (net[:ϑ̲][i] - ϑ[i,t,s]) : NaN

        (ϱ[i,t,s] - net[:ρ̅][i]) >= 0 ? ϱ_mag[t] += (ϱ[i,t,s] - net[:ρ̅][i]) : NaN
        (net[:ρ̲][i] - ϱ[i,t,s]) >= 0 ? ϱ_mag[t] += (net[:ρ̲][i] - ϱ[i,t,s]) : NaN
        (ϱ[i,t,s] - net[:ρ̅][i]) >= 0 ? ϱ_cvar[t,s] += (ϱ[i,t,s] - net[:ρ̅][i]) : NaN
        (net[:ρ̲][i] - ϱ[i,t,s]) >= 0 ? ϱ_cvar[t,s] += (net[:ρ̲][i] - ϱ[i,t,s]) : NaN
    end
    for t in 1:set[:T], i in net[:E_a], s in 1:set[:Ns]
        ϕ[i,t,s] <= 0  ? ϕ_mag[t] += ϕ[i,t,s] : NaN
        ϕ[i,t,s] <= 0  ? ϕ_cvar[t,s] += ϕ[i,t,s] : NaN

        (κ[i,t,s] - net[:κ̅][i]) >= 0  ? κ_mag[t] += (κ[i,t,s] - net[:κ̅][i]) : NaN
        (net[:κ̲][i] - κ[i,t,s]) >= 0  ? κ_mag[t] += (net[:κ̲][i] - κ[i,t,s]) : NaN
        (κ[i,t,s] - net[:κ̅][i]) >= 0  ? κ_cvar[t,s] += (κ[i,t,s] - net[:κ̅][i]) : NaN
        (net[:κ̲][i] - κ[i,t,s]) >= 0  ? κ_cvar[t,s] += (net[:κ̲][i] - κ[i,t,s]) : NaN
    end
    for i in net[:E], s in 1:set[:Ns]
        net[:ψ_o][i] - ψ[i,set[:T],s] >= 0 ? ψ_mag[set[:T]] += (net[:ψ_o][i] - ψ[i,set[:T],s]) : NaN
        net[:ψ_o][i] - ψ[i,set[:T],s] >= 0 ? ψ_cvar[set[:T],s] += (net[:ψ_o][i] - ψ[i,set[:T],s]) : NaN
    end
    ϑ_mag = ϑ_mag./set[:Ns]
    ϱ_mag = ϱ_mag./set[:Ns]/1000
    ϕ_mag = ϱ_mag./set[:Ns]
    κ_mag = ϱ_mag./set[:Ns]/1000
    ψ_mag = ψ_mag./set[:Ns]

    ϑ_cvar = mean(sort(sum(ϑ_cvar,dims=1),dims=2)[950:end])
    ϱ_cvar = mean(sort(sum(ϱ_cvar,dims=1),dims=2)[950:end])/1000
    ψ_cvar = mean(sort(sum(ψ_cvar,dims=1),dims=2)[950:end])
    ϕ_cvar = mean(sort(sum(ϕ_cvar,dims=1),dims=2)[950:end])
    κ_cvar = mean(sort(sum(κ_cvar,dims=1),dims=2)[950:end])/1000

    ### Expected cost
    cost_exp = [net[:c1]'*solution[:Θ][t]*S(t)*ξ[:μ̂] + tr(solution[:Θ][t]'*diagm(net[:c2])*solution[:Θ][t]*S(t)*ξ[:Σ̂]*S(t)') for t in 1:set[:T]]
    mean_lp_deployment = sum(sum(ψ[:,t,s]) for t in 1:set[:T] for s in 1:set[:Ns])*1/set[:Ns]
    mean_compr_deployment = sum(sum(κ[i,t,s]) for i in setdiff(net[:E_a],[51 50]) for t in 1:set[:T] for s in 1:set[:Ns])*1/set[:Ns]
    mean_valve_deployment = sum(sum(abs.(κ[i,t,s])) for i in [51 50] for t in 1:set[:T] for s in 1:set[:Ns])*1/set[:Ns]
    return Dict(:mean_compr_deployment => mean_compr_deployment,
                :mean_lp_deployment => mean_lp_deployment,
                :mean_valve_deployment => mean_valve_deployment,
                :cost_exp => sum(cost_exp),
                :press_violation_exp => sum(ϱ_mag .+ κ_mag),
                :mass_violation_exp => sum(ϑ_mag .+ ϕ_mag .+ ψ_mag),
                :press_violation_wc => sum(ϱ_cvar .+ κ_cvar),
                :mass_violation_wc => sum(ϑ_cvar .+ ϕ_cvar .+ ψ_cvar),
                :first_stage_inj => sum(solution[:Θ][1]*S(1)*ξ[:μ̂]))
end

function cc_opt_MISOCP(net,set,top_data)
    C = Int(2^length(set[:b_valve_edge]))
    # build model
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    # variable declaration
    P =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    F =  [[@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for c = 1:C] for t in 1:set[:T]]
    Z =  [[@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for c = 1:C] for t in 1:set[:T]]
    v =  [@variable(model, [1:1], Bin) for c in 1:C]
    Φ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    K =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Θ =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁻ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁺ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Ψ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    s_ϱ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    # minimize gas production cost
    @objective(model, Min, sum(net[:c1]' * Θ[t] * S(t) * ξ[:μ̂]
                                + tr(Θ[t]' * diagm(net[:c2]) * Θ[t] * S(t) * ξ[:Σ̂] * S(t)')
                           for t in 1:set[:T])
                           + set[:α_ϱ] * sum(sum(s_ϱ[t]) for t in 2:set[:T])
            )
    # variance constraints
    @constraint(model, pre_var[t=2:set[:T], i=net[:N]],
        [s_ϱ[t][i];(cholesky(ξ[:Σ̂]).U .+ zeros(set[:k_t][set[:T]],set[:k_t][set[:T]]))*(P[t]*S(t) - P[t-1]*S(t-1))[i,:]] in SecondOrderCone())
    # equivalents of stochastic equations
    @constraint(model, gas_bal[t=1:set[:T]], (net[:A⁺]*Φ⁺[t]*S(t) + net[:A⁻]*Φ⁻[t]*S(t) - Θ[t]*S(t) + net[:B]*K[t]*S(t) + Δ(t)*S(t)) .== 0)
    @constraint(model, lin_wey_1[t=1:set[:T],c=1:C], (F[t][c]*S(t) - [top_data[c][:W_0][t] zeros(length(net[:E]),set[:k_t][t]-1)]*S(t) - top_data[c][:W_1][t]*P[t]*S(t) - top_data[c][:W_2][t]*K[t]*S(t)) .== 0)
    @constraint(model, lin_wey_2[t=1:set[:T]], (Φ[t] .- sum(Z[t][c] for c in 1:C))*S(t) .== 0)
    @constraint(model, bilin_1[t=1:set[:T],c=1:C,i=net[:E],j=1:t], set[:M̲] <= Z[t][c][i,j] <= set[:M̅])
    @constraint(model, bilin_2[t=1:set[:T],c=1:C,i=1,j=1:t], set[:M̲]*v[c] .<= Z[t][c][i,j])
    @constraint(model, bilin_3[t=1:set[:T],c=1:C,i=1,j=1:t], set[:M̅]*v[c] .>= Z[t][c][i,j])
    @constraint(model, bilin_4[t=1:set[:T],c=1:C,i=1,j=1:t], F[t][c][i,j] .- (1 .- v[c]) * set[:M̅] .<= Z[t][c][i,j])
    @constraint(model, bilin_5[t=1:set[:T],c=1:C,i=1,j=1:t], F[t][c][i,j] .- (1 .- v[c]) * set[:M̲] .>= Z[t][c][i,j])
    @constraint(model, bin_con, sum(v) .== 1)
    @constraint(model, ref_pre[t=1:set[:T]], (P[t][net[:ref],:]' - sum(v[c].*[top_data[c][:point][:ϱ][t][net[:ref]] zeros(1,set[:k_t][t]-1)] for c in 1:C))*S(t) .== 0)
    @constraint(model, av_flow[t=1:set[:T]], (Φ[t]*S(t) - 1/2*Φ⁺[t]*S(t) - 1/2*Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_def[t=1:set[:T]], (Ψ[t]*S(t) - 1/2*diagm(net[:s])*(K[t] + abs.(net[:A])'*P[t])*S(t)) .== 0)
    @constraint(model, l_p_bal_1[t=1], (Ψ[t]*S(t) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_bal_t[t=2:set[:T]], (Ψ[t]*S(t) - Ψ[t-1]*S(t-1) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    # reformulated chance constraints
    @constraint(model, ϑ_max[t=1:set[:T], i=1:length(net[:N])],
        [(net[:ϑ̅] - Θ[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(Θ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ϑ_min[t=1:set[:T], i=1:length(net[:N])],
        [(Θ[t]*S(t)*ξ[:μ̂] .- net[:ϑ̲])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(Θ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ρ_max[t=1:set[:T], i=1:length(net[:N])],
        [(net[:ρ̅] - P[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(P[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ρ_min[t=1:set[:T], i=1:length(net[:N])],
        [(P[t]*S(t)*ξ[:μ̂] - net[:ρ̲])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(P[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, κ_max[t=1:set[:T], i=1:length(net[:E])],
        [(net[:κ̅] - K[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(K[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, κ_min[t=1:set[:T], i=1:length(net[:E])],
        [(K[t]*S(t)*ξ[:μ̂] - net[:κ̲])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(K[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ψ_T[t=set[:T], i=1:length(net[:E])],
        [(Ψ[t]*S(t)*ξ[:μ̂] - net[:ψ_o])[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(Ψ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, φ_pl[t=1:set[:T], i=net[:E_a]],
        [(Φ[t]*S(t)*ξ[:μ̂] .- 0)[i]; z(set[:ε̅_g]*100)*ξ[:Σ̅]*(Φ[t]*S(t))[i,:]]
            in SecondOrderCone())
    # solve model
    optimize!(model)
    @info("done solving topology optimization in $(solve_time(model)) sec. with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :v => [JuMP.value.(v[c]) for c in 1:C],
    :P => [JuMP.value.(P[t]) for t in 1:set[:T]],
    :Φ => [JuMP.value.(Φ[t]) for t in 1:set[:T]],
    :K => [JuMP.value.(K[t]) for t in 1:set[:T]],
    :Θ => [JuMP.value.(Θ[t]) for t in 1:set[:T]],
    :Φ⁻ => [JuMP.value.(Φ⁻[t]) for t in 1:set[:T]],
    :Φ⁺ => [JuMP.value.(Φ⁺[t]) for t in 1:set[:T]],
    :Ψ => [JuMP.value.(Ψ[t]) for t in 1:set[:T]],
    :cost => sum(net[:c1]' * JuMP.value.(Θ[t]) * S(t) * ξ[:μ̂] + tr(JuMP.value.(Θ[t])' * diagm(net[:c2]) * JuMP.value.(Θ[t]) * S(t) * ξ[:Σ̂] * S(t)') for t in 1:set[:T]),
    :ϱ_std_tot => sum(sum(JuMP.value.(α_ϱ[t])) for t in 2:1:set[:T]) ,
    :status => termination_status(model),
    :obj => JuMP.objective_value(model),
    :model => model)
    return solution
end

function cc_opt_MISOCP_2(net,set,top_data)
    C = Int(2^length(set[:b_valve_edge]))
    # build model
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "MSK_DPAR_MIO_TOL_REL_GAP" => 1e-5))
    # model = Model(optimizer_with_attributes(Gurobi.Optimizer))
    # set_optimizer_attribute(model, "MIQCPMethod", 1)
    # variable declaration
    P =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    F =  [[@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for c = 1:C] for t in 1:set[:T]]
    Z =  [[@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for c = 1:C] for t in 1:set[:T]]
    v =  [@variable(model, [1:1], Bin) for c in 1:C]
    Φ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    K =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Θ =  [@variable(model, [1:length(net[:N]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁻ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Φ⁺ = [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    Ψ =  [@variable(model, [1:length(net[:E]), 1:set[:k_t][t]]) for t = 1:set[:T]]
    α_ϱ =  [@variable(model, [1:length(net[:N])]) for t = 1:set[:T]]
    # minimize gas production cost
    @objective(model, Min, sum(net[:c1]' * Θ[t] * S(t) * ξ[:μ̂]
                                + tr(Θ[t]' * diagm(net[:c2]) * Θ[t] * S(t) * ξ[:Σ̂] * S(t)')
                           for t in 1:set[:T])
                           + set[:α_ϱ] * sum(sum(α_ϱ[t]) for t in 2:set[:T])
            )
    # variance constraints
    @constraint(model, pre_var[t=2:set[:T], i=net[:N]],
        [α_ϱ[t][i];(cholesky(ξ[:Σ̂]).U .+ zeros(set[:k_t][set[:T]],set[:k_t][set[:T]]))*(P[t]*S(t) - P[t-1]*S(t-1))[i,:]] in SecondOrderCone())
    # equivalents of stochastic equations
    @constraint(model, gas_bal[t=1:set[:T]],
        (net[:A⁺]*Φ⁺[t]*S(t) + net[:A⁻]*Φ⁻[t]*S(t) - Θ[t]*S(t) + net[:B]*K[t]*S(t) + Δ(t)*S(t)) .== 0)
    @constraint(model, lin_wey[t=1:set[:T],c=1:C],
        (F[t][c]*S(t) - [top_data[c][:W_0][t] zeros(length(net[:E]),set[:k_t][t]-1)]*S(t) - top_data[c][:W_1][t]*P[t]*S(t) - top_data[c][:W_2][t]*K[t]*S(t)) .== 0)
    @constraint(model, lin_wey_2[t=1:set[:T]],
        (Φ[t] .- sum(Z[t][c] for c in 1:C))*S(t) .== 0)
    @constraint(model, lin_wey_3[t=1:set[:T]],
        (P[t][net[:ref],:]' -  sum(v[c].*[top_data[c][:point][:ϱ][t][net[:ref]] zeros(1,set[:k_t][t]-1)] for c in 1:C)) * S(t) .== 0)
    @constraint(model, av_flow[t=1:set[:T]],
        (Φ[t]*S(t) - 1/2*Φ⁺[t]*S(t) - 1/2*Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_def[t=1:set[:T]],
        (Ψ[t]*S(t) - 1/2*diagm(net[:s])*(K[t] + abs.(net[:A])'*P[t])*S(t)) .== 0)
    @constraint(model, l_p_bal_1[t=1],
        (Ψ[t]*S(t) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    @constraint(model, l_p_bal_t[t=2:set[:T]],
        (Ψ[t]*S(t) - Ψ[t-1]*S(t-1) - Φ⁺[t]*S(t) + Φ⁻[t]*S(t)) .== 0)
    # linearization of the bilinear term
    @constraint(model, bilin_1[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] <= set[:M̅])
    @constraint(model, bilin_2[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] >= set[:M̲])
    @constraint(model, bilin_3[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] .<= set[:M̅] * v[c])
    @constraint(model, bilin_4[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] .>= set[:M̲] * v[c])
    @constraint(model, bilin_5[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] .<=  F[t][c][i,j] .- (1 .- v[c]) * set[:M̲])
    @constraint(model, bilin_6[t=1:set[:T],c=1:C,i=net[:E],j=1:t],
        Z[t][c][i,j] .>=  F[t][c][i,j] .- (1 .- v[c]) * set[:M̅])
    # sum up all binaries to one
    @constraint(model, sum(v) .== 1)
    # reformulated chance constraints
    @constraint(model, ϑ_max[t=1:set[:T], i=1:length(net[:N])],
        [(net[:ϑ̅] - Θ[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Θ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ϑ_min[t=1:set[:T], i=1:length(net[:N])],
        [(Θ[t]*S(t)*ξ[:μ̂] .- net[:ϑ̲])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Θ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ρ_max[t=1:set[:T], i=1:length(net[:N])],
        [(net[:ρ̅] - P[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(P[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ρ_min[t=1:set[:T], i=1:length(net[:N])],
        [(P[t]*S(t)*ξ[:μ̂] - net[:ρ̲])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(P[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, κ_max[t=1:set[:T], i=1:length(net[:E])],
        [(net[:κ̅] - K[t]*S(t)*ξ[:μ̂])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(K[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, κ_min[t=1:set[:T], i=1:length(net[:E])],
        [(K[t]*S(t)*ξ[:μ̂] - net[:κ̲])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(K[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, ψ_T[t=set[:T], i=1:length(net[:E])],
        [(Ψ[t]*S(t)*ξ[:μ̂] - net[:ψ_o])[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Ψ[t]*S(t))[i,:]]
            in SecondOrderCone())
    @constraint(model, φ_pl[t=1:set[:T], i=net[:E_a]],
        [(Φ[t]*S(t)*ξ[:μ̂] .- 0)[i]; z(set[:ε̅_g])*ξ[:Σ̅]*(Φ[t]*S(t))[i,:]]
            in SecondOrderCone())
    # solve model
    optimize!(model)
    @info("cc model terminates with status: $(termination_status(model))")
    # return solution
    solution = Dict(
    :v => [JuMP.value.(v[c]) for c in 1:C],
    :P => [JuMP.value.(P[t]) for t in 1:set[:T]],
    :Φ => [JuMP.value.(Φ[t]) for t in 1:set[:T]],
    :K => [JuMP.value.(K[t]) for t in 1:set[:T]],
    :Θ => [JuMP.value.(Θ[t]) for t in 1:set[:T]],
    :Φ⁻ => [JuMP.value.(Φ⁻[t]) for t in 1:set[:T]],
    :Φ⁺ => [JuMP.value.(Φ⁺[t]) for t in 1:set[:T]],
    :Ψ => [JuMP.value.(Ψ[t]) for t in 1:set[:T]],
    :cost => sum(net[:c1]' * JuMP.value.(Θ[t]) * S(t) * ξ[:μ̂] + tr(JuMP.value.(Θ[t])' * diagm(net[:c2]) * JuMP.value.(Θ[t]) * S(t) * ξ[:Σ̂] * S(t)') for t in 1:set[:T]),
    :ϱ_std_tot => sum(sum(JuMP.value.(α_ϱ[t])) for t in 2:1:set[:T]) ,
    :status => termination_status(model),
    :obj => JuMP.objective_value(model),
    :model => model)
    return solution
end


text:
# Julia wrapper for header: /usr/include/libavutil/cast5.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_cast5_alloc,
    av_cast5_init,
    av_cast5_crypt,
    av_cast5_crypt2

function av_cast5_alloc()
    ccall((:av_cast5_alloc, libavutil), Ptr{AVCAST5}, ())
end

function av_cast5_init(ctx, key, key_bits::Integer)
    ccall((:av_cast5_init, libavutil), Cint, (Ptr{AVCAST5}, Ptr{UInt8}, Cint), ctx, key, key_bits)
end

function av_cast5_crypt(ctx, dst, src, count::Integer, decrypt::Integer)
    ccall((:av_cast5_crypt, libavutil), Cvoid, (Ptr{AVCAST5}, Ptr{UInt8}, Ptr{UInt8}, Cint, Cint), ctx, dst, src, count, decrypt)
end

function av_cast5_crypt2(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_cast5_crypt2, libavutil), Cvoid, (Ptr{AVCAST5}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end


text:
#InspectDR: Drawing Smith/polar grids with Cairo layer
#-------------------------------------------------------------------------------

#TODO: Try to collect Smith chart code in this module, if it makes sense.

#==Constants
===============================================================================#
const SMITHLABEL_OFFSET = Float64(2) #Device units from intercept

#==Types
===============================================================================#

#==Base functions
===============================================================================#
function smithlabel(v::Real)
	vint = round(Int, v)
	if vint == v
		return "$vint"
	else
		return "$v"
	end
end

function render_rcircle(ctx::CairoContext, c::DReal)
	pctr = Point2D(c / (1+c), 0)
	radius = 1 / (1+c)
	Cairo.arc(ctx, pctr.x, pctr.y, radius, 0, 2pi)
	Cairo.stroke(ctx)
end
render_rcircles(ctx::CairoContext, clist::Vector{DReal}) =
	(for c in clist; render_rcircle(ctx, c); end)

function render_xcircle(ctx::CairoContext, c::DReal)
	radius = 1/c
	Cairo.arc(ctx, 1, radius, radius, 0, 2pi)
	Cairo.stroke(ctx)
	Cairo.arc(ctx, 1, -radius, radius, 0, 2pi)
	Cairo.stroke(ctx)
end
render_xcircles(ctx::CairoContext, clist::Vector{DReal}) =
	(for c in clist; render_xcircle(ctx, c); end)

#xflip: Mirror x-direction for admittance(Y)-labels
function render_rcirclelabel(ctx::CairoContext, xf::Transform2D, pt::Point2D, xflip::Bool, lbl::DisplayString)
	xscale = xflip? -1: 1
	xalign = xflip? ALIGN_LEFT: ALIGN_RIGHT
	pt = ptmap(xf, Point2D(xscale*pt.x, pt.y))
	render(ctx, lbl, Point2D(pt.x-xscale*SMITHLABEL_OFFSET, pt.y+SMITHLABEL_OFFSET),
		align=ALIGN_TOP|xalign
	)
end
function render_rcirclelabel(ctx::CairoContext, xf::Transform2D, xflip::Bool, refscale::Real, c::DReal)
	xintercept = (c-1)/(c+1)
	pt = Point2D(xintercept, 0)
	cstr = smithlabel(c*refscale)
	tstr = DisplayString("$cstr")
	render_rcirclelabel(ctx, xf, pt, xflip, tstr)
end
render_rcirclelabels(ctx::CairoContext, xf::Transform2D, xflip::Bool, refscale::Real, clist::Vector{DReal}) =
	(for c in clist; render_rcirclelabel(ctx, xf, xflip, refscale::Real, c); end)

#xflip: Mirror x-direction for admittance(Y)-labels
function render_xcirclelabel(ctx::CairoContext, xf::Transform2D, pt::Point2D, xflip::Bool, lbl::DisplayString)
	xscale = xflip? -1: 1
	corr = -xscale*xf.xs/xf.ys #Correct for aspect ratio
	#Compute text angle:
	Θ = atan2(pt.x, pt.y*corr)
	if xflip; Θ+=pi; end

	#Compute text position, including offset:
	pt = Point2D(xscale*pt.x, pt.y)
	vpos = vecmap(xf, pt) #Directional vector of label position
	voffset = (SMITHLABEL_OFFSET/vecnorm(vpos))*vpos
	pt = ptmap(xf, pt)+voffset
	render(ctx, lbl, pt, angle=Θ, align=ALIGN_BOTTOM|ALIGN_HCENTER)
end
function render_xcirclelabel(ctx::CairoContext, xf::Transform2D, xflip::Bool, refscale::Real, c::DReal)
	radius = 1/c
	#Compute intercept with 0dB circle (r=1):
	y = 2*radius/(1+radius*radius)
	x = 1-y*radius

	cstr = smithlabel(c*refscale)
	tstr = DisplayString("$(cstr)j")
	render_xcirclelabel(ctx, xf, Point2D(x, y), xflip, tstr)
	tstr = DisplayString("-$(cstr)j")
	render_xcirclelabel(ctx, xf, Point2D(x, -y), xflip, tstr)
end
render_xcirclelabels(ctx::CairoContext, xf::Transform2D, xflip::Bool, refscale::Real, clist::Vector{DReal}) =
	(for c in clist; render_xcirclelabel(ctx, xf, xflip, refscale, c); end)

#==Extensions to base functionnality
===============================================================================#

#Split complex input data into real/imag components:
function getrealimag{T<:Number}(d::Vector{T})
	x = Array(DReal, length(d))
	y = Array(DReal, length(d))
	for i in 1:length(d)
		x[i] = convert(DReal, real(d[i]))
		y[i] = convert(DReal, imag(d[i]))
	end
	return (x, y)
end

function getrealimag(input::IDataset)
	(x, y) = getrealimag(input.y)
	return IDataset{false}(x, y) #Assume new re/im data is not sorted.
end

function getrealimag(input::IWaveform)
	ds = getrealimag(input.ds)
	return IWaveform(input.id, ds, input.line, input.glyph, getextents(ds))
end

_rescale(inputlist::Vector{IWaveform}, axes::AxesSmith) =
	map((input)->getrealimag(input), inputlist)

_rescale(pt::Point2D, axes::AxesSmith) = pt

#==High-level rendering functions
===============================================================================#
function render_grid(canvas::PCanvas2D, lyt::Layout, grid::GridSmith)
	const ctx = canvas.ctx
	const graphbb = canvas.graphbb
	const MINOR_COLOR = GRID_MINOR_COLOR
	const MAJOR_COLOR = COLOR_BLACK #GRID_MAJOR_COLOR
	const xflip = !grid.zgrid #Flip x-axis for admittance(Y)-grid lines

Cairo.save(ctx) #-----
	setclip(ctx, graphbb)

	#Set transform through Cairo to so it can deal with circles:
	setxfrm(ctx, canvas.xf)
	if xflip
		Cairo.scale(ctx, -1, 1)
	end

	setlinestyle(ctx, :solid, GRID_MINOR_WIDTH)
	Cairo.set_source(ctx, MINOR_COLOR)
		render_rcircles(ctx, grid.minorR)
		render_xcircles(ctx, grid.minorX)
	setlinestyle(ctx, :solid, GRID_MAJOR_WIDTH)
	Cairo.set_source(ctx, MAJOR_COLOR)
		render_rcircles(ctx, grid.majorR)
Cairo.restore(ctx) #-----

Cairo.save(ctx) #-----
	setclip(ctx, graphbb)

	#Draw X=0 line:
	setlinestyle(ctx, :solid, GRID_MINOR_WIDTH)
	Cairo.set_source(ctx, MINOR_COLOR)
	y = ptmap(canvas.xf, Point2D(0, 0)).y
	drawline(ctx::CairoContext, Point2D(graphbb.xmin, y), Point2D(graphbb.xmax, y))

	#Draw labels:
	setfont(ctx, lyt.fntticklabel)
	render_rcirclelabels(ctx, canvas.xf, xflip, grid.ref, grid.labelR)
	render_rcirclelabel(ctx, canvas.xf, Point2D(1,0), xflip, "∞")
	render_xcirclelabels(ctx, canvas.xf, xflip, grid.ref, grid.labelX)

Cairo.restore(ctx) #-----
	nothing
end

function render_axes(canvas::PCanvas2D, lyt::Layout, grid::GridSmith)
	render_graphframe(canvas)

	#Display linear grid:
	#TODO: make it possible to disable?
	grid = gridlines(AxesRect(:lin, :lin), canvas.ext)
	render_ticks(canvas, lyt, grid)
end

#Last line


text:
#=
From http://paulbourke.net/dataformats/ply/

name        type        number of bytes
---------------------------------------
char       character                 1
uchar      unsigned character        1
short      short integer             2
ushort     unsigned short integer    2
int        integer                   4
uint       unsigned integer          4
float      single-precision float    4
double     double-precision float    8
=#
typemap = Dict(
    Char => "char",
    UInt8 => "uchar",
    Int16 => "short",
    UInt16 => "ushort",
    Int32 => "int",
    Int32 => "uint",
    Float32 => "float",
    Float64 => "double",
)

function save(f::Stream{format"PLY_BINARY"}, msh::AbstractMesh)
    io = stream(f)
    points = decompose(Point{3, Float32}, msh)
    faces = decompose(GLTriangleFace, msh)

    n_points = length(points)
    n_faces = length(faces)

    # write the header
    write(io, "ply\n")
    write(io, "format binary_little_endian 1.0\n")
    write(io, "element vertex $n_points\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    write(io, "element face $n_faces\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    write(io, points)

    for f in faces
        write(io, convert(UInt8, 3))
        write(io, raw.(f)...)
    end
    close(io)
end

function save(f::Stream{format"PLY_ASCII"}, msh::AbstractMesh)
    io = stream(f)
    points = coordinates(msh)
    meshfaces = faces(msh)

    n_faces = length(points)
    n_points = length(meshfaces)

    point_type = eltype(points)
    if point_type <: PointMeta
        metafields = point_type.parameters[4]
        metatypes = point_type.parameters[5].parameters
    else
        metafields = []
        metatypes = []
    end

    # write the header
    write(io, "ply\n")
    write(io, "format ascii 1.0\n")
    write(io, "element vertex $n_faces\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    for (metafield, metatype) in zip(metafields, metatypes)
        type = typemap[metatype]
        write(io, "property $type $metafield\n")
    end
    write(io, "element face $n_points\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    for v in points
        metavals = map(f -> getproperty(v, f), metafields)
        print(io, join(Point{3, Float32}(v), " "))
        print(io, " ")
        print(io, join(metavals, " "))
        print(io, "\n")

    end
    for f in meshfaces
        println(io, length(f), " ", join(raw.(ZeroIndex.(f)), " "))
    end
    close(io)
end

function load(fs::Stream{format"PLY_ASCII"}; facetype=GLTriangleFace, pointtype=Point3f0)
    io = stream(fs)
    n_points = 0
    n_faces = 0

    properties = String[]

    # read the header
    line = readline(io)

    while !startswith(line, "end_header")
        if startswith(line, "element vertex")
            n_points = parse(Int, split(line)[3])
        elseif startswith(line, "element face")
            n_faces = parse(Int, split(line)[3])
        elseif startswith(line, "property")
            push!(properties, line)
        end
        line = readline(io)
    end

    faceeltype = eltype(facetype)
    points = Array{pointtype}(undef, n_points)
    #faces = Array{FaceType}(undef, n_faces)
    faces = facetype[]

    # read the data
    for i = 1:n_points
        points[i] = pointtype(parse.(eltype(pointtype), split(readline(io)))) # line looks like: "-0.018 0.038 0.086"
    end

    for i = 1:n_faces
        line = split(readline(io))
        len = parse(Int, popfirst!(line))
        if len == 3
            push!(faces, NgonFace{3, faceeltype}(reinterpret(ZeroIndex{Int}, parse.(Int, line)))) # line looks like: "3 0 1 3"
        elseif len == 4
            push!(faces, convert_simplex(facetype, QuadFace{faceeltype}(reinterpret(ZeroIndex{Int}, parse.(Int, line))))...) # line looks like: "4 0 1 2 3"
        end
    end
    return Mesh(points, faces)
end


text:
# This file was generated, do not modify it. # hide
using PlotlyJS
function random_line()
    a=0.1

    us = 0:0.05π:2π
    ts = -1:0.05:1
    traces=typeof(scatter3d())[]
    marker=attr(
        color="#1f77b4",
        size=12,
        symbol="circle",
        line=attr(color="rgb(0,0,0)", width=0),
    )
    for t in ts
        x=@. (t^2+a)*cos(us)
        y=@. (t^2+a)*sin(us)
        z=t*ones(length(x))
        push!(
            traces,
            scatter3d(
                x=x,y=y,z=z,
                mode="lines",
                marker=marker,
                showlegend=false,
            )
        )
    end

    for u in us
        x = @. (ts^2+a)*cos(u)
        y = @. (ts^2+a)*sin(u)
        z = ts
        push!(
            traces,
            scatter3d(
                x=x,y=y,z=z,
                mode="lines",
                marker=marker,
                showlegend=false
            )
        )
    end

    layout = Layout(
        autosize=true,
        width=500,
        height=500,
        margin=attr(l=0, r=0, b=150, t=0),
    )
    plot(traces,layout)
end

p = random_line()
fdplotly(json(p)) # hide


text:
function ueberpruefeneueFormel(Formel_neu)
A=[1 2; 3 4];
f = [1 2];
m = [-1 -2];

if Formel_neu(f,A,m) == (A[1,1]*f[1]+A[1,2]*f[2]-m[1])^2+(A[2,1]*f[1]+A[2,2]*f[2]-m[2])^2
    println("Eure Umformungen waren korrekt!");
    println("Wir müssen einzeln nach f_1 und f_2 ableiten. Im Ergebnis steht die Ableitung nach f_1 in der ersten und nach f_2 in der zweiten Zeile.");
    println("Das Ziel ist es hierfür wieder eine Kurzschreibweise zu finden. Fasst dafür euer Ergbnis soweit wie möglich zusammen.");
    println("Berechnet dafür transpose(A)*A*f und transpose(A)*m. Schaut, ob euch das bekannt vorkommt.");
    println("transpose(A) nennt man die Transponierte von A.");
    println("Macht all diese Schritte auf einem Papier.");
elseif Formel_neu(f,A,m) == sqrt((A[1,1]*f[1]+A[1,2]*f[2]-m[1])^2+(A[2,1]*f[1]+A[2,2]*f[2]-m[2])^2)
    println("Eure Formel ist noch nicht ganz gerecht. Bedenkt, dass ihr das Quadrat des Abstands berechnet.");
else
    println("Eure Formel stimmt noch nicht. Stellt sicher, dass ihr keine Rechenfehler bei der Matrix-Vektor-Multiplikation gemacht habt.");
end

end


text:
function Gatherer(stYr = 1975,endYr = 1980,windowSize = 5)

#load the data and add to a dictionary the number of countries in the year range
countryYearsNum = Dict{Integer,Integer}()
resultsFile = readdir("./dataTables/")
yrMin = 100000
yrMax = -1
for rf in resultsFile
    fileTmp = open(string("./dataTables/",rf))
    linesTmp = readlines(fileTmp)         #readfile lines
    countryNumTmp = length(split(linesTmp[1],",")) - 1#COLUMN NAME LIST
    yrTmp = parse(Int,((split(rf,"."))[1]))
    countryYearsNum[yrTmp] = countryNumTmp #NUM that receives votes
    close(fileTmp)
    if(yrTmp < yrMin)
        yrMin = yrTmp
    end
    if(yrTmp > yrMax)
        yrMax = yrTmp
    end
end

#sanity check input years
if( endYr < stYr || stYr < 1975 || endYr > 2003 )
    print(string("year range improperly set, for the Gatherer approach end year must be greater than start and the smallest year is ",1975," and largest ",2003," with smallest first"))
    return
end
if( (stYr+windowSize) > endYr)
    print("not enough years between start and end for analysis")
    return
end

#THRESHOLDS
windowConf = Dict() #the threshold for significance in each year of window
yr = stYr
while( (yr+windowSize) <= endYr )
    conf5perc = GathererWindow(yr,yr+windowSize,countryYearsNum)
    yr = yr + windowSize
    windowConf[string(yr-windowSize,"-",yr)] = conf5perc
end

#MATRIX WINDOW AVGS

print(windowConf)
winDicts = windowScores(stYr,endYr,windowSize)
print(winDicts)

collusionDict = countryCollusion(stYr,endYr,windowSize,windowConf,winDicts)

print(collusionDict)
end

function countryCollusion(stYr,endYr,windowSize,windowConf,winDicts)
#FROM THE THRESHOLDS DETERMINE THE COUNTRY PAIRS OF COLLUSION

    collusionDict = Dict()
    yr = stYr
    while( (yr+windowSize) <= endYr )

	threshold = windowConf["$(yr)-$(yr + windowSize)"]

	cntryNames = winDicts["$(yr)-$(yr + windowSize)"]["countries"]
	scoremat = winDicts["$(yr)-$(yr + windowSize)"]["scoremat"]
	scorematAVG = scoremat * (1/(windowSize+1))

	collusionDict["$(yr)-$(yr + windowSize)"] = []
	for ii=1:(size(scoremat)[1])
	    for jj=1:(size(scoremat)[2])
		if(scorematAVG[ii,jj] >= threshold)
		    c1 = cntryNames[ii]
		    c2 = cntryNames[jj]
		    prev = collusionDict["$(yr)-$(yr + windowSize)"]

		    new = append!(prev, ["$(c1)-$(c2)"])
		    collusionDict["$(yr)-$(yr + windowSize)"] = new
		end
	    end
	end
	yr = yr + windowSize
    end

return collusionDict
end

#FOR EACH TIME WINDOW
function windowScores(stYr,endYr,windowSize)

#1-DICTIONARY FOR EACH WINDOW INTERVAL
winDict = Dict()
yr = stYr
while( (yr+windowSize) <= endYr )
    winDict["$(yr)-$(yr+windowSize)"] = Dict()
    winDict["$(yr)-$(yr+windowSize)"]["countries"] = []
    winDict["$(yr)-$(yr+windowSize)"]["scoremat"] = []
    yr = yr + windowSize
end

#2-FILL DICTIONARY WITH ALL POSSIBLE COUNTRY PAIRS IN THOSE YEARS IT COVERS
resultsFile = readdir("./dataTables/")
namesDict = Dict()
scoresDict = Dict()
for rf in resultsFile
    yrTmp = parse(Int,((split(rf,"."))[1]))
    fileTmp = open(string("./dataTables/",rf))                     #each file pipe
    linesTmp = readlines(fileTmp)                            #read each file lines
    tmp = split(linesTmp[1],r",|\n",keep=false) #store the name list
    splice!(tmp,1) #remove first
    namesDict[yrTmp] = tmp
    nameNum = length(namesDict[yrTmp])
    rowNumMat = length(linesTmp)-1
    colNumMat = nameNum
    tmp2 = zeros(rowNumMat,colNumMat)
    for ii=2:length(linesTmp)
	tmp1 = split(linesTmp[ii],r",|\n",keep=false)
	splice!(tmp1,1) #remove first
	tmp2[ii-1,:] = [parse(Int,tmp1[jj]) for jj in 1:length(tmp1)]	#each row of the matrix/data
    end
    scoresDict[string(yrTmp)] = tmp2
    close(fileTmp)
end

yr = stYr
while( (yr+windowSize) <= endYr )
    winInd = 0
    while(winInd <= windowSize)
	new = namesDict[yr + winInd]
        prev = winDict["$(yr)-$(yr+windowSize)"]["countries"]
        winDict["$(yr)-$(yr+windowSize)"]["countries"] = sort(unique(append!(prev,new)))
	winInd = winInd + 1
    end
    yr = yr + windowSize
end

#3-ACCUMULATE POINTS OF EACH PARTICULAR COUNTRY TOWARDS ANOTHER DIFFERENT COUNTRY
yr = stYr
while( (yr+windowSize) <= endYr )
    winInd = 0
    cntryNum = length(winDict["$(yr)-$(yr+windowSize)"]["countries"])#WILL BECOME ARRAAY FOR ROW AND COL SIZE

    winDict["$(yr)-$(yr+windowSize)"]["scoremat"] = zeros(cntryNum,cntryNum)#WILL BECOME ASSYMETRIC!!!
    while(winInd <= windowSize)

        matPrev = winDict["$(yr)-$(yr+windowSize)"]["scoremat"]
	matNew  = zeros(length(matPrev[:,1]),length(matPrev[1,:]))
        namesTotal = winDict["$(yr)-$(yr+windowSize)"]["countries"]
	yrScoresMat = scoresDict[string(yr+winInd)]#the matrix for a PARTICULAR YEAR (might not be nxn)
        yrNames = namesDict[yr+winInd]

	for rowNum=1:length(yrScoresMat[:,1])#change the year nameInds to the total window Inds of names
	    newIndsMat = []
	    for ind=1:length(find(yrScoresMat[rowNum,:]))
		tmpNamesMat = yrNames[find(yrScoresMat[rowNum,:])]
		newIndsMat = sort(append!(find([namesTotal[ii] == tmpNamesMat[ind] for ii in 1:length(namesTotal)]),newIndsMat))
	    end
	    rowNumNew = find([namesTotal[ii] == yrNames[rowNum] for ii in 1:length(namesTotal)])
	    matNew[rowNumNew,newIndsMat] = yrScoresMat[rowNum,find(yrScoresMat[rowNum,:])]
	end

        winDict["$(yr)-$(yr+windowSize)"]["scoremat"] = matPrev + matNew

        winInd = winInd + 1
    end
    yr = yr + windowSize
end

return winDict

end

#THRESHOLD FOR EACH TIME WINDOW
function GathererWindow(stYr,endYr,countryYearsNum)

SCORES = [12,10,8,7,6,5,4,3,2,1]

AVG_SIMULATION = []
iterNum = 250
confInd5perc = max(1,floor(Int,0.05*iterNum))

for ii = 1:iterNum
    ONE_SIMULATION = []
    for yr = stYr:endYr
	NUM = countryYearsNum[yr]#number of countries voting that year
	position = ceil(rand(1,1)*NUM)

	if position[1] <= length(SCORES)
	   score = SCORES[position]
	else
	   score = 0
	end
	append!(ONE_SIMULATION,[score])

    end
    avgSim = mean(ONE_SIMULATION)
    append!(AVG_SIMULATION,[avgSim])

end

sortedAVG_SIMULATION = sort(AVG_SIMULATION,rev=true)
conf5perc = sortedAVG_SIMULATION[confInd5perc]

return conf5perc

end

max_stars_repo_name: mantzaris/eurovision

max_issues_repo_name: mantzaris/eurovision

max_forks_repo_name: mantzaris/eurovision


text:
is_cron(event) = get(ENV, "GITHUB_EVENT_NAME", "") == "schedule"

function handle_cron(event)
    pulls = collect_pulls(ENV["GITHUB_REPOSITORY"])
    repos_versions = map(pull -> repo_and_version_of_pull_request_body(pull.body), pulls)
    filter!(rv -> first(rv) !== nothing, repos_versions)
    unique!(first, repos_versions)  # Send at most one notification per repo.
    for (repo, version) in repos_versions
        maybe_notify(event, repo, version; check_tag=true)
    end
end

function collect_pulls(repo)
    acc = GH.PullRequest[]
    params = (; state="closed", sort="updated", direction="desc", per_page=100)
    kwargs = Dict(:auth => AUTH[], :params => params, :page_limit => 1)
    done = false
    while !done
        get_pulls = retry(; check=(s, e) -> occursin("Server error", e.msg)) do
            GH.pull_requests(repo; kwargs...)
        end
        pulls, pages = get_pulls()
        for pull in pulls
            pull.merged_at === nothing && continue
            if now(UTC) - pull.merged_at < Day(3)
                push!(acc, pull)
            else
                done = true
            end
        end
        if haskey(pages, "next")
            delete!(kwargs, :params)
            kwargs[:start_page] = pages["next"]
        else
            done = true
        end
    end
    return acc
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
using Revise,PlutoUI,ExtendableGrids

# ╔═╡ df8aefc0-5273-11eb-1a67-3552732deae0
using Plots

# ╔═╡ 08fa72d6-5274-11eb-03bd-ef885fdd216d
P=Plots

# ╔═╡ e5cec5c8-5273-11eb-0e0b-03960befcc71
X=collect(0:0.5:10)

# ╔═╡ f088d080-5273-11eb-2085-e5d4734512ff
g1=simplexgrid(X)

# ╔═╡ ff091d22-5273-11eb-2645-89630b592fa0
gridplot(g1,Plotter=P,resolution=(500,100))

# ╔═╡ 2dd685fe-5274-11eb-173c-6f2ccaf41aae
gridplot(g1,map(x->sin(x),g1),Plotter=P)

# ╔═╡ 8d1857ae-5274-11eb-265d-09604d9dba94
g2=simplexgrid(X,X)

# ╔═╡ 9526886c-5274-11eb-33a6-d95a16547472
gridplot(g2,Plotter=P)

# ╔═╡ ae81ae16-5274-11eb-36e6-a1dd77e51abb
gridplot(g2,map((x,y)->(sin(x)*exp(-(y-5)^2)),g2),Plotter=P,colormap=:hot)

# ╔═╡ eb380f5a-5274-11eb-2a7a-87c75706a2ae
p=GridPlotContext(Plotter=P,layout=(1,2),fignum=2,resolution=(800,400))

# ╔═╡ 0114cb24-5275-11eb-3013-0f3b170f25de
begin
	gridplot!(p[1,1],g1,map(x->sin(x),g1),title="1D")
	gridplot!(p[1,2],g2,map((x,y)->(sin(x)*exp(-(y-5)^2)),g2),colormap=:hot,title="2D")
	reveal(p)
end

# ╔═╡ Cell order:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╠═df8aefc0-5273-11eb-1a67-3552732deae0
# ╠═08fa72d6-5274-11eb-03bd-ef885fdd216d
# ╠═e5cec5c8-5273-11eb-0e0b-03960befcc71
# ╠═f088d080-5273-11eb-2085-e5d4734512ff
# ╠═ff091d22-5273-11eb-2645-89630b592fa0
# ╠═2dd685fe-5274-11eb-173c-6f2ccaf41aae
# ╠═8d1857ae-5274-11eb-265d-09604d9dba94
# ╠═9526886c-5274-11eb-33a6-d95a16547472
# ╠═ae81ae16-5274-11eb-36e6-a1dd77e51abb
# ╠═eb380f5a-5274-11eb-2a7a-87c75706a2ae
# ╠═0114cb24-5275-11eb-3013-0f3b170f25de


text:
module Ccolamd

using SparseArrays
using SuiteSparse.CHOLMOD: SuiteSparse_long

const KNOBS = 20
const STATS = 20

function recommended(nnz::SuiteSparse_long, n_row::SuiteSparse_long, n_col::SuiteSparse_long)
    Alen = ccall((:ccolamd_l_recommended, :libccolamd), Csize_t,
        (SuiteSparse_long, SuiteSparse_long, SuiteSparse_long),
        nnz, n_row, n_col)
    if Alen == 0
        error("")
    end
    return Alen
end

function set_defaults(knobs::Vector{Float64})
    if length(knobs) != KNOBS
        error("")
    end
    ccall((:ccolamd_set_defaults, :libccolamd), Nothing, (Ptr{Cdouble},), knobs)
    return knobs
end

function ccolamd!(n_row::SuiteSparse_long, A::Vector{SuiteSparse_long},
    p::Vector{SuiteSparse_long}, knobs::Union{Ptr{Nothing},Vector{Float64}}, stats::Vector{SuiteSparse_long},
    cmember::Union{Ptr{Nothing},Vector{SuiteSparse_long}})

    n_col = length(p) - 1

    if length(stats) != STATS
        error("stats must hcae length $STATS")
    end
    if isa(cmember, Vector) && length(cmember) != n_col
        error("cmember must have length $n_col")
    end

    Alen = recommended(length(A), n_row, n_col)
    resize!(A, Alen)

    for i in eachindex(A)
        A[i] -= 1
    end
    for i in eachindex(p)
        p[i] -= 1
    end
    err = ccall((:ccolamd_l, :libccolamd), SuiteSparse_long,
        (SuiteSparse_long, SuiteSparse_long, SuiteSparse_long, Ptr{SuiteSparse_long},
         Ptr{SuiteSparse_long}, Ptr{Cdouble}, Ptr{SuiteSparse_long}, Ptr{SuiteSparse_long}),
        n_row, n_col, Alen, A, p, knobs, stats, cmember)
    if err == 0
        report(stats)
        error("call to ccolamd return with error code $(stats[4])")
    end

    for i in eachindex(p)
        p[i] += 1
    end

    pop!(p) # remove last zero from pivoting vector
    return p
end

function ccolamd!(n_row, A::Vector{SuiteSparse_long}, p::Vector{SuiteSparse_long}, cmember::Union{Ptr{Nothing},Vector{SuiteSparse_long}})
    n_col = length(p) - 1

    if length(cmember) != n_col
        error("cmember must have length $n_col")
    end

    Alen = recommended(length(A), n_row, n_col)
    resize!(A, Alen)
    stats = zeros(SuiteSparse_long, STATS)
    return ccolamd!(n_row, A, p, C_NULL, stats, cmember)
end

function ccolamd!(n_row, A::Vector{SuiteSparse_long}, p::Vector{SuiteSparse_long})
    n_col = length(p) - 1
    cmember = zeros(SuiteSparse_long, length(p) - 1)
    return ccolamd!(n_row, A, p, cmember)
end

ccolamd(n_row, A::Vector{SuiteSparse_long}, p::Vector{SuiteSparse_long}) =
    ccolamd!(n_row, copy(A), copy(p))

ccolamd(A::SparseMatrixCSC) = ccolamd(size(A, 1), A.rowval, A.colptr)

report(stats::Vector{SuiteSparse_long}) =
    ccall((:ccolamd_l_report, :libccolamd), Nothing, (Ptr{SuiteSparse_long},), stats)

end #module


text:
using JuLIPMaterials: Vec3, Mat3, Ten33, Ten43
using JuLIPMaterials.CLE: onb3D, euclidean, spherical

using Einsum, StaticArrays

using LinearAlgebra

export GreenFunction, IsoGreenFcn3D, grad

# ========== Implementation of the 3D Green's Function =============

abstract type AbstractGreenFunction{T} end

struct GreenFunction3D{T} <: AbstractGreenFunction{T}
   Nquad::Int
   C::Ten43{T}
   remove_singularity::Bool
end

"""
`GreenFunction3D`
construct a green's function type
"""
function GreenFunction(C::Ten43; Nquad = nothing, remove_singularity = true)
   if Nquad == nothing
      error("still need to implement auto-tuning, please provide Nquad")
   end
   return GreenFunction3D(Nquad, C, remove_singularity)
end

GreenFunction(C::Array{T, 4}; kwargs...) where {T} =
      GreenFunction(Ten43{T}(C); kwargs...)

GreenFunction(at::AbstractAtoms; kwargs...) =
      GreenFunction(calculator(at), at; kwargs...)

GreenFunction(calc::AbstractCalculator, at::AbstractAtoms; kwargs...) =
      GreenFunction(elastic_moduli(calc, at); kwargs...)

function (G::GreenFunction3D{T})(x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return eval_green(Vec3(x), G.C, G.Nquad)
end

function grad(G::GreenFunction3D{T}, x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SArray zeros(T, 3, 3, 3)
   end
   return grad_green(Vec3(x), G.C, G.Nquad)
end

"eval_green(x::Vec3, ℂ::Ten43, Nquad::Int)"
function eval_green(x::Vec3{TT}, ℂ::Ten43, Nquad::Int) where TT
   # allocate
   G = @SMatrix zeros(TT, 3, 3)
   zz = @MMatrix zeros(TT, 3, 3)
   # Initialise tensors.
   x̂ = x/norm(x)
   # two vectors orthogonal to x.
   x1, x2 = onb3D(x̂)
   # Integrate
   for ω in range(0.0, step=pi/Nquad, length=Nquad)
      z = cos(ω) * x1 + sin(ω) * x2
      @einsum zz[i,j] = z[α] * ℂ[i,α,j,β] * z[β]
      # Perform integration
      G += inv(zz)
   end
   # Normalise appropriately
   return G / (4*pi*norm(x)*Nquad)
end

function grad_green(x::Vec3{TT}, ℂ::Ten43, Nquad::Int) where TT
   # allocate
   DG = @MArray zeros(TT, 3, 3, 3)
   zz = @MMatrix zeros(TT, 3, 3)
   zT = @MMatrix zeros(TT, 3, 3)
   # Initialise tensors.
   x̂ = x/norm(x)
   # two vectors orthogonal to x.
   x1, x2 = onb3D(x̂)
   # Integrate
   for ω in range(0.0, step=pi/Nquad, length=Nquad)
      z = cos(ω) * x1 + sin(ω) * x2
      @einsum zz[i,j] = z[α] * ℂ[i,α,j,β] * z[β]
      @einsum zT[i,j] = z[α] * ℂ[i,α,j,β] * x̂[β]
      zzinv = inv(zz)
      F = zzinv * (zT + zT') * zzinv
      @einsum DG[i,j,k] = DG[i,j,k] + zzinv[i,j] * x̂[k] - F[i,j] * z[k]
   end
   DG ./= (-4.0 * pi * norm(x)^2 * Nquad)
   return SArray(DG)
end

struct IsoGreenFcn3D{T} <: AbstractGreenFunction{T}
   λ::T
   μ::T
   remove_singularity::Bool
end

IsoGreenFcn3D(λ, μ; remove_singularity = true) =
   IsoGreenFcn3D(λ, μ, remove_singularity)

function (G::IsoGreenFcn3D{T})(x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return eval_greeniso(Vec3(x), G.λ, G.μ)
end

function grad(G::IsoGreenFcn3D{T}, x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SArray zeros(T, 3, 3, 3)
   end
   return grad_greeniso(Vec3{T}(x), G.λ, G.μ)
end

"isotropic CLE Green's function"
function eval_greeniso(x::Vec3{T}, λ::Real, μ::Real) where T
   Id = SMatrix{3,3,T}(I)
   x̂ = x/norm(x)
   return (((λ+3*μ)/(λ+2*μ)/norm(x)) * Id  +
               ((λ+μ)/(λ+2*μ)/norm(x)) * x̂ * x̂') / (8.0*π*μ)
end

"gradient of isotropic CLE Green's function"
function grad_greeniso(x::Vec3{T}, λ::Real, μ::Real) where T
   DG = @MArray zeros(T, 3, 3, 3)
   Id = SMatrix{3,3,T}(I)
   x̂ = x / norm(x)
   for i = 1:3, j = 1:3, k = 1:3
      DG[i,j,k] = (λ+μ) * (Id[i,k] * x̂[j] + Id[j,k] * x̂[i]) - (λ+3*μ) * Id[i,j] * x̂[k] - 3*(λ+μ) * x̂[i] * x̂[j] * x̂[k]
   end
   DG ./= 8 * π * μ*(λ+2*μ) * norm(x)^2
   return DG
end