text:
function est_mmvmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    #p = size(Y,2)
    #n = size(Y,1)
    model.w_esteps = zeros(Squantities.n,Squantities.n_esteps,model.g)
    model.Ψ = reshape(repeat(Squantities.psi_initial,model.g),(size(Squantities.psi_initial)[1],model.g))
    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        model.niter = 0
        return
    elseif model.status == "failedE"
        model.status = "Failed: E-step for Ψ failed at initialisation" #probably caused by a bad numerical algorithm ... such as finite differences for GH
        model.niter = 0
        return
    end

    model.niter = Squantities.itmax
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        model.niter = niter
        Squantities.Mstep_calc!(model, Mquantities, Squantities) #Need to update this function
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            model.niter = niter
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            model.niter = niter
            break
        elseif model.status == "failedE"
            model.status = "Failed: E-step for Ψ failed after $niter iterations" #probably caused by a bad numerical algorithm ... such as finite differences for GH
            model.niter = niter
            break
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            break
        end

        Mquantities.logL_old = model.logL
        model.logL_history[niter + 1] = model.logL
    end #for
    #model.niter = Squantities.itmax
end

function est_msmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    model.w_esteps = zeros(Squantities.n,Squantities.n_esteps,model.g)
    model.Ψ = reshape(repeat(Squantities.psi_initial,model.g),(size(Squantities.psi_initial)[1],model.g))
    model.beta = zeros(Squantities.p,model.g)

    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        return
    end
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        Squantities.Mstep_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            break
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            break
        end
        model.logL_history[niter+1] = model.logL
        Mquantities.logL_old = model.logL

    end #for
    #if model == "failedModel"
    #    return "failedModel"
    #else
    #    return model
    #end
end

function est_mfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)

    #p = size(Y,2)
    #n = size(Y,1)
    model.beta =  zeros(Squantities.p,model.g)
    #init_para.beta = zeros(p,init_para.g)

    Squantities.logL_calc!(model, Mquantities, Squantities)
    if model.status == "failed"
        model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after initialisation"
        return
    end
    Mquantities.logL_old = model.logL
    model.logL_history[1] = model.logL
    @inbounds for niter in 1:Squantities.itmax
        Squantities.Mstep_calc!(model, Mquantities, Squantities) #Need to update this function
        if model.status == "failed"
            model.status = "Failed: got a complex eigenvalue in the M-steps after $niter iterations"
            return
        end

        Squantities.logL_calc!(model, Mquantities, Squantities)
        if model.status == "failed"
            model.status = "Failed: got a NaN in a log-likelihod or E-step for Ψ after $niter iterations"
            return
        end
        if ((Squantities.conv_measure == "diff") && (abs(model.logL - Mquantities.logL_old) < Squantities.tol))
            model.niter = niter
            model.logL_history[niter + 1] = model.logL
            break
        end

        if ((Squantities.conv_measure == "ratio") && (abs((model.logL - Mquantities.logL_old)/model.logL) < Squantities.tol))
            model.niter = niter
            model.logL_history[niter + 1] = model.logL
            break
        end

        Mquantities.logL_old = model.logL
        model.logL_history[niter + 1] = model.logL
    end #for
end

max_stars_repo_name: john-c-davey/FactorMixtures

max_issues_repo_name: john-c-davey/FactorMixtures

max_forks_repo_name: john-c-davey/FactorMixtures


text:
# list of point position + velocity.  need to figure out what they
# "spell"...maybe can use min distance from origin as a heuristic here?  it
# seems like all the inputs are large but have opposite signs from their
# velocies (on a scale of like 10_000x or something).  or sum of pairwise
# distances?
#
# or just find the relative position/velocity of each pair of point and use that
# to get the overall (sum) pariwise distance as a function of time then find the
# minimum...
#
# ((a_0 - b_0) + (a_v - b_v)*t)^2 = Δ_0^2 + 2 Δ_0 Δ_v t + Δ_v^2 t^2
# derivative of that is 2Δ_0Δ_v + 2Δ_v^2 t => max/min at t = Δ0 / Δv
#
# that's just for one.  the quadratic coefficients sum but need to do the sum
# first though...

mutable struct Star
    x::Tuple{Int,Int}
    v::Tuple{Int,Int}
end

function parseline(line)
    regex = r"position=< *(-?\d*), *(-?\d*)> velocity=< *(-?\d*), *(-?\d*)>"
    x, y, vx, vy = parse.(Int, match(regex, line).captures)
    Star((x,y), (vx, vy))
end

test_input = parseline.(split(strip(
    """position=< 9,  1> velocity=< 0,  2>
position=< 7,  0> velocity=<-1,  0>
position=< 3, -2> velocity=<-1,  1>
position=< 6, 10> velocity=<-2, -1>
position=< 2, -4> velocity=< 2,  2>
position=<-6, 10> velocity=< 2, -2>
position=< 1,  8> velocity=< 1, -1>
position=< 1,  7> velocity=< 1,  0>
position=<-3, 11> velocity=< 1, -2>
position=< 7,  6> velocity=<-1, -1>
position=<-2,  3> velocity=< 1,  0>
position=<-4,  3> velocity=< 2,  0>
position=<10, -3> velocity=<-1,  1>
position=< 5, 11> velocity=< 1, -2>
position=< 4,  7> velocity=< 0, -1>
position=< 8, -2> velocity=< 0,  1>
position=<15,  0> velocity=<-2,  0>
position=< 1,  6> velocity=< 1,  0>
position=< 8,  9> velocity=< 0, -1>
position=< 3,  3> velocity=<-1,  1>
position=< 0,  5> velocity=< 0, -1>
position=<-2,  2> velocity=< 2,  0>
position=< 5, -2> velocity=< 1,  2>
position=< 1,  4> velocity=< 2,  1>
position=<-2,  7> velocity=< 2, -2>
position=< 3,  6> velocity=<-1, -1>
position=< 5,  0> velocity=< 1,  0>
position=<-6,  0> velocity=< 2,  0>
position=< 5,  9> velocity=< 1, -2>
position=<14,  7> velocity=<-2,  0>
position=<-3,  6> velocity=< 2, -1>
"""), '\n'))

input = parseline.(readlines("10.input"))

move(star::Star, t::Int) = Star(star.x .+ t .* star.v, star.v)

# compute the coefficients of the quadratic expression of squared distance for a pair
coefs(a::Star, b::Star) = (sum((a.x .- b.x).^2),
                           sum(2 .* (a.x .- b.x) .* (a.v .- b.v)),
                           sum((a.v .- b.v).^2)

function tmin(input)
    a,b,c = reduce((x,y)->x.+y, coefs(a,b) for a in input for b in input)
    tmin = round(Int, -b / (2c))
end

function draw(stars::Vector{Star})
    mins, maxs = foldl(((minx,maxx),x) -> (min.(minx,x), max.(maxx,x)),
                       s.x for s in stars,
                       init=(typemax(Int), typemin(Int)))
    img = fill(' ', maxs .- mins .+ 1)
    for s in stars
        img[(s.x .- mins .+ 1)...] = '█'
    end
    cols, rows = axes(img)
    println()
    for r in rows
        for c in cols
            print(img[c,r])
        end
        println()
    end
    img
end

star1(input) = draw(move.(input, tmin(input)))

star1(test_input);

star1(input);

# star 2:
tmin(input)

ext: jl


text:
## labelmap

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem) where {T}
    labelmap!(dict, idx, T(elem))
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem::T) where {T}
    if !haskey(dict, elem)
        push!(dict, elem => [idx])
    else
        push!(dict[elem], idx)
    end
    dict
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::AbstractVector{Int}, itr::S) where {T,S}
    for (i, elem) in zip(idx, itr)
        labelmap!(dict, i, elem)
    end
    dict
end

function labelmap(itr::T) where {T}
    dict = Dict{eltype(T),Vector{Int}}()
    for (idx, elem) in enumerate(itr)
        labelmap!(dict, idx, elem)
    end
    dict
end

## labelfreq

function labelfreq!(dict::Dict{T,Int}, elem) where {T}
    labelfreq!(dict, T(elem))
end

function labelfreq!(dict::Dict{T,Int}, elem::T) where {T}
    cnt = get(dict, elem, 0)
    dict[elem] = cnt + 1
    dict
end

function labelfreq!(dict::Dict{T,Int}, iter::AbstractVector) where {T}
    for elem in iter
        labelfreq!(dict, T(elem))
    end
    dict
end

function rarelabel(labels::Array,threshold::Integer,replace_label::Symbol)
    label_freq = labelfreq(labels)
    for i = 1:length(labels)
        if label_freq[labels[i]]<=threshold
            labels[i]=replace_label
        end
    end
    return labels
end

labelfreq(itr::T) where {T} = countmap(itr, alg = :dict)::Dict{eltype(T),Int}
labelfreq(dict::Dict{T,Vector{Int}})  where {T} = Dict(k => length(v) for (k,v) in dict)::Dict{T,Int}

## Errors

labelmap(A::AbstractMatrix)  = throw(ArgumentError("labelmap not supported for matrices"))
labelfreq(A::AbstractMatrix) = throw(ArgumentError("labelfreq not supported for matrices"))

## General Dict

labelenc(x::Dict) = labelenc(label(x))
nlabel(x::Dict) = length(keys(x))
label(x::Dict) = _arrange_label(collect(keys(x)))


text:
using ..Atmos
using ..ConfigTypes
using ..DiagnosticsMachine

@diagnostics_group(
    "AtmosLESDefault",      # name
    AtmosLESConfigType,     # configuration type
    Nothing,                # params type
    (_...) -> nothing,      # initialization function
    NoInterpolation,        # if/when to interpolate
    # various horizontal averages
    u,
    v,
    w,
    rho,
    temp,
    pres,
    thd,
    et,
    ei,
    ht,
    hi,
    w_ht_sgs,
    # moisture related
    qt,
    ql,
    qi,
    qv,
    thv,
    thl,
    w_qt_sgs,
    # for variances and co-variances
    uu,
    vv,
    ww,
    www,
    eiei,
    wu,
    wv,
    wrho,
    wthd,
    wei,
    qtqt,
    thlthl,
    wqt,
    wql,
    wqi,
    wqv,
    wthv,
    wthl,
    qtthl,
    qtei,
    #cld_top, TODO
    #cld_base, TODO
    #cld_cover, TODO
    #lwp, TODO
    #rwp, TODO
)


text:
function raw_to_counts(raw)
    x_values = convert(Array{Float64}, unique(sort(raw[:, 2])))

    shots = Array(Array{UInt32, 1}, length(x_values))
    for i in 1:length(x_values)
        shots[i] = UInt32[]
    end

    for i in 1:(size(raw)[1])
        x_index = searchsorted(x_values, raw[i, 2])[1]
        for val in raw[i, 4:end]
            if val == ""
                break
            end

            counts = convert(UInt32, val)
            if counts & (1 << 31) != 0
                continue
            end

            push!(shots[x_index], counts)
        end
    end

    x_values, shots
end

function last_complete_row_2d(raw, y_column)
    lines_per_y = 0
    for (i, y) in enumerate(slice(raw, :, y_column))
        if y != raw[1, y_column]
            lines_per_y = i - 1
            break
        end
    end

    row_count = size(raw)[1]
    if lines_per_y == 0
       lines_per_y = row_count
    end

    if row_count % lines_per_y == 0
        return row_count, lines_per_y
    end

    for (i, y) in enumerate(slice(raw, row_count:-1:1, y_column))
        if y != raw[row_count, y_column]
            current = row_count + 1 - i
            return current, lines_per_y
        end
    end

    error("Dataset has only one y value, cannot determine whether it is complete.")
end

function raw_to_counts_2d(raw, lines_per_y::Int64)
    @assert(size(raw)[1] % lines_per_y == 0, "Inconsistent number of measurements per iterator value")

    y_values = convert(Array{Float64}, sort(raw[1:lines_per_y:end, 3]))
    x_values = convert(Array{Float64}, unique(sortrows(slice(raw, 1:lines_per_y, 1:2))[:, 2]))

    shots = Array(Array{UInt32, 1}, length(x_values), length(y_values))
    for j in 1:length(y_values)
        for i in 1:length(x_values)
            shots[i, j] = UInt32[]
        end
    end

    for i in 1:(size(raw)[1])
        x_index = searchsorted(x_values, raw[i, 2])[1]
        y_index = searchsorted(y_values, raw[i, 3])[1]
        for val in raw[i, 5:end]
            if val == ""
                break
            end

            counts = convert(UInt32, val)
            if counts & (1 << 31) != 0
                continue
            end

            push!(shots[x_index, y_index], counts)
        end
    end

    x_values, y_values, shots
end

function walk_match(dir, prefix, suffix)
    if isempty(dir)
        dir = "."
    end
    entries = readdir(dir)
    for e in entries
        if !startswith(e, prefix)
            continue
        end

        path = joinpath(dir, e)
        if endswith(path, suffix)
            return path
        end

        # Recurse into child directories.
        if isdir(path)
            p = walk_match(path, "", suffix)
            if !isempty(p)
                return p
            end
        end
    end
    return ""
end

function find_file(datapath, suffix)
    if endswith(datapath, suffix)
        return datapath
    end

    dirname, start = splitdir(datapath)
    path = walk_match(dirname, start, suffix)
    if isempty(path)
        error("No file found in '$datapath' ending in '$suffix'.")
    end
    return path
end

function read_shots_log(datapath)
    raw_to_counts(readcsv(find_file(datapath, "_shots_log.txt")))
end

function read_shots_log_2d(datapath; y_reversed=false)
    data = readcsv(find_file(datapath, "_2D_shots_log.txt"))
    last_row, lines_per_y = last_complete_row_2d(data, 3)

    if y_reversed
        clipped_data = data[last_row:-1:1, :]
    else
        clipped_data = data[1:last_row, :]
    end

    raw_to_counts_2d(clipped_data, lines_per_y)
end

function read_thresholded_2d(datapath)
    # Read in file, skipping header line
    data = readcsv(find_file(datapath, "_plot_0.csv"))[2:end, :]

    last_row, lines_per_y = last_complete_row_2d(data, 2)

    y_values = convert(Array{Float64}, sort(data[1:lines_per_y:end, 2]))
    x_values = convert(Array{Float64}, unique(sort(data[1:lines_per_y, 1])))

    signal = zeros(Float64, length(x_values), length(y_values))
    for i in 1:(size(data)[1])
        x_index = searchsorted(x_values, data[i, 1])[1]
        y_index = searchsorted(y_values, data[i, 2])[1]
        signal[x_index, y_index] += data[i, 3]
    end

    signal ./= (lines_per_y / length(x_values))

    x_values, y_values, signal
end


text:
include("./imports.jl")

@testset ExtendedTestSet "DataAugmentation.jl" begin
    @testset ExtendedTestSet "base.jl" begin
        include("base.jl")
    end
    @testset ExtendedTestSet "testing.jl" begin
        include("testing.jl")
    end
    @testset ExtendedTestSet "buffered.jl" begin
        include("buffered.jl")
    end
    @testset ExtendedTestSet "projective/base.jl" begin
        include("projective/base.jl")
    end
    @testset ExtendedTestSet "projective/compose.jl" begin
        include("projective/compose.jl")
    end
    @testset ExtendedTestSet "projective/crop.jl" begin
        include("projective/crop.jl")
    end
    @testset ExtendedTestSet "projective/affine.jl" begin
        include("projective/affine.jl")
    end
    @testset ExtendedTestSet "projective/warp.jl" begin
        include("projective/warp.jl")
    end
    @testset ExtendedTestSet "items/image.jl" begin
        include("items/image.jl")
    end
    @testset ExtendedTestSet "items/keypoints.jl" begin
        include("items/keypoints.jl")
    end
    @testset ExtendedTestSet "items/mask.jl" begin
        include("items/mask.jl")
    end
    @testset ExtendedTestSet "preprocessing.jl" begin
        include("preprocessing.jl")
    end
    @testset ExtendedTestSet "oneof.jl" begin
        include("oneof.jl")
    end
    @testset ExtendedTestSet "visualization.jl" begin
        include("visualization.jl")
    end
    @testset ExtendedTestSet "rowtransforms.jl" begin
        include("rowtransforms.jl")
    end
end


text:
function propagation(settings::PyObject, data::PyObject, x_obs::Array{Float64, 1}, msap_source, r, x, z, c_bar, rho_0, I_0, beta)

    # Impedance at the observer
    I_0_obs = 409.74

    # Size of inputs
    n_t = size(r)[1]

    # Get type of input variables
    T = eltype(msap_source)
    msap_prop = zeros(T, (n_t, settings.N_f))

    # Apply spherical spreading and characteristic impedance effects to the MSAP
    # Source: Zorumski report 1982 part 1. Chapter 5.1 Equation 1
    z = max.(z, 0)
    if settings.direct_propagation
        msap_prop = msap_source .* ((settings.r_0^2 ./ r.^2) .* (I_0_obs./I_0))
    else
        msap_prop = msap_source
    end

    # Split mean-square acoustic pressure in frequency sub-bands
    if settings.N_b > 1
        msap_sb = split_subbands(settings, msap_prop)
    end

    # Apply atmospheric absorption on sub-bands
    if settings.absorption
        # Calculate average absorption factor between observer and source
        f_abs = LinearInterpolation((data.abs_alt, data.abs_freq), data.abs)
        alpha_f = f_abs.(z.*ones(T, (1, settings.N_b*settings.N_f)), reshape(data.f_sb, (1,settings.N_f*settings.N_b)).*ones(T, (n_t,1)))

        # Calculate absorption (convert dB to Np: 1dB is 0.115Np)
        # Source: Zorumski report 1982 part 1. Chapter 5.1 Equation 14
        msap_prop = msap_sb .* exp.(-2 * 0.115 * alpha_f .* (r .- settings.r_0))
    end

    # Apply ground effects on sub-bands
    if settings.groundeffects
        # Empirical lateral attenuation for microphone on sideline
        if (settings.lateral_attenuation == true) && (x_obs[2] != 0)
            # Lateral attenuation factor
            Lambda = lateral_attenuation(settings, beta, x_obs)

            # Compute elevation angle from "center-line" observer.
            r_cl = sqrt.((x.-x_obs[1]).^2 .+ 1. ^2 .+ z.^2)
            beta_cl = asin.(z./r_cl) * 180. / pi

            # Ground reflection factor for center-line
            G = ground_reflections(settings, data, x_obs, r_cl, beta_cl, c_bar, rho_0)

            # Apply ground effects
            msap_prop = msap_prop .* (G .* Lambda)
        else
            # Compute ground effects factor
            G = ground_reflections(settings, data, x_obs, r, beta, c_bar, rho_0)

            # Multiply subbands with G
            msap_prop = msap_prop .* G
        end
    end

    # Recombine the mean-square acoustic pressure in the frequency sub-bands
    if settings.N_b > 1
        msap_prop = transpose(reshape(sum(reshape(transpose(msap_prop), (settings.N_b, settings.N_f, n_t)), dims=1), (settings.N_f, n_t)))
    end

    return msap_prop
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
using JGUI
using Base.Test

w = window(title="notebook")
nb = notebook(w); push!(w, nb)
connect(nb, "valueChanged", value -> println("On tab $value"))

b1 = button(nb, "one")
b2 = button(nb, "two")
b3 = button(nb, "three")

push!(nb, b1, "one")
push!(nb, b2, "two")
push!(nb, b3, "three")

setValue(nb, 3)
@assert nb[:value] == 3

pop!(nb, b2)
@assert length(nb) == 2
@assert nb[1][:value] == "one"

raise(w)


text:
# EXCEPTIONS

function squareroot(x)
    try
        return sqrt(x)
    catch
        return sqrt(x + 0im)
    end
end

squareroot(2)
squareroot(-2)

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
#=
district_cooling:
- Julia version:
- Author: qgis1
- Date: 2019-03-28
=#
module district_cooling
    using JuMP
    using Cbc

    using CSV
    using DelimitedFiles
    csv_sep = ','

    function district_cooling_simulation(input_folder, output_folder, tech_infos, network_id, print_function = print)
        m = Model(with_optimizer(Cbc.Optimizer, logLevel=1, seconds=600))

        ####### Load distributions #########

        DEM=readdlm(string(input_folder, "\\DEM_time.csv"), csv_sep)
        #DEM=DEM[4381:6570]

        Electricity_price=readdlm(string(input_folder, "\\Electricity_price_time.csv"), csv_sep)
        #Electricity_price=Electricity_price[4381:6570]

        Heat_exchanger_specific=readdlm(string(input_folder, "\\Heat_exchanger_specific_time.csv"), csv_sep)
        #Heat_exchanger_specific=Heat_exchanger_specific[4381:6570]

        Heat_exchanger_specific_cool=readdlm(string(input_folder, "\\Heat_exchanger_specific_time_cool.csv"), csv_sep)
        #Heat_exchanger_specific_cool=Heat_exchanger_specific_cool[4381:6570]

        eta_HP = readdlm(string(input_folder, "\\eta_HP_1.csv"), csv_sep)
        #eta_HP = eta_HP[4381:6570]
        eta_HP_2 = readdlm(string(input_folder, "\\eta_HP_2.csv"), csv_sep)
        #eta_HP_2 = eta_HP_2[4381:6570]
        eta_HP_3 = readdlm(string(input_folder, "\\eta_HP_3.csv"), csv_sep)
        #eta_HP_3 = eta_HP_3[4381:6570]

        eta_HP_waste_heat_1=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_1.csv"), csv_sep)
        #eta_HP_waste_heat_1=eta_HP_waste_heat_1[4381:6570]
        eta_HP_waste_heat_2=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_2.csv"), csv_sep)
        #eta_HP_waste_heat_2=eta_HP_waste_heat_2[4381:6570]
        eta_HP_waste_heat_3=readdlm(string(input_folder, "\\COP_heat_pump_temperature_group_3.csv"), csv_sep)
        #eta_HP_waste_heat_3=eta_HP_waste_heat_3[4381:6570]

        Available_waste_heat_heat_pump_1=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_1.csv"), csv_sep)
        #Available_waste_heat_heat_pump_1=Available_waste_heat_heat_pump_1[4381:6570]
        Available_waste_heat_heat_pump_2=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_2.csv"), csv_sep)
        #Available_waste_heat_heat_pump_2=Available_waste_heat_heat_pump_2[4381:6570]
        Available_waste_heat_heat_pump_3=readdlm(string(input_folder, "\\Waste_heat_heat_pump_available_temperature_group_3.csv"), csv_sep)
        #Available_waste_heat_heat_pump_3=Available_waste_heat_heat_pump_3[4381:6570]

        ####### Technology input data ########

        ## Heat only boiler ##

        # Installed capacity [MW]
        P_max_HOB = tech_infos["P_max_HOB"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB= tech_infos["Tecnical_minimum_HOB"] #0.2
        P_min_HOB = Tecnical_minimum_HOB*P_max_HOB
        # Efficiency [-]
        eta_HOB = tech_infos["eta_HOB"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB= tech_infos["Percentage_ramp_up_down_HOB"] #0.5
        ramp_up_down_HOB = P_max_HOB*Percentage_ramp_up_down_HOB
        # Variable cost [EUR/MWh]
        cost_var_HOB = tech_infos["cost_var_HOB"] #5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB = tech_infos["cost_fuel_HOB"] #10

        ## Heat only boiler 2 ##

        # Installed capacity [MW]
        P_max_HOB_2 = tech_infos["P_max_HOB_2"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB_2= tech_infos["Tecnical_minimum_HOB_2"] #0.3
        P_min_HOB_2 = Tecnical_minimum_HOB_2*P_max_HOB_2
        # Efficiency [-]
        eta_HOB_2 = tech_infos["eta_HOB_2"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_2= tech_infos["Percentage_ramp_up_down_HOB_2"] #0.5
        ramp_up_down_HOB_2 = P_max_HOB_2*Percentage_ramp_up_down_HOB_2
        # Variable cost [EUR/MWh]
        cost_var_HOB_2 = tech_infos["cost_var_HOB_2"] #5.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_2 = tech_infos["cost_fuel_HOB_2"] #12.5

        ## Heat only boiler 3 ##

        # Installed capacity [MW]
        P_max_HOB_3 = tech_infos["P_max_HOB_3"] #2
        # Technical minimum [MW]
        Tecnical_minimum_HOB_3= tech_infos["Tecnical_minimum_HOB_3"] #0.2
        P_min_HOB_3 = Tecnical_minimum_HOB_3*P_max_HOB_3
        # Efficiency [-]
        eta_HOB_3 = tech_infos["eta_HOB_3"] #0.8
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HOB_3= tech_infos["Percentage_ramp_up_down_HOB_3"] #0.6
        ramp_up_down_HOB_3 = P_max_HOB_3*Percentage_ramp_up_down_HOB_3
        # Variable cost [EUR/MWh]
        cost_var_HOB_3 = tech_infos["cost_var_HOB_3"] #3.4
        # Fuel cost [EUR/MWh]
        cost_fuel_HOB_3 = tech_infos["cost_fuel_HOB_3"] #20

        ## Heat pump ##

        # Installed capacity [MW]
        P_max_HP = tech_infos["P_max_HP"] #2
        # Technical minimum [MW]
        Technical_minimum_HP= tech_infos["Technical_minimum_HP"] #0.2
        P_min_HP = Technical_minimum_HP*P_max_HP
        # Efficiency [-]
        eta_HP = eta_HP
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP= tech_infos["Percentage_ramp_up_down_HP"] #0.9
        ramp_up_down_HP = Percentage_ramp_up_down_HP*P_max_HP
        # Variable cost [EUR/MWh]
        cost_var_HP = tech_infos["cost_var_HP"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP = Electricity_price

        ## Heat pump 2 ##

        # Installed capacity [MW]
        P_max_HP_2 = tech_infos["P_max_HP_2"] #2
        # Technical minimum [MW]
        Technical_minimum_HP_2= tech_infos["Technical_minimum_HP_2"] #0.1
        P_min_HP_2 = Technical_minimum_HP_2*P_max_HP_2
        # Efficiency [-]
        eta_HP_2 = eta_HP_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_2= tech_infos["Percentage_ramp_up_down_HP_2"] #0.95
        ramp_up_down_HP_2 = Percentage_ramp_up_down_HP_2*P_max_HP_2
        # Variable cost [EUR/MWh]
        cost_var_HP_2 = tech_infos["cost_var_HP_2"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_2 = Electricity_price

        ## Heat pump 3 ##

        # Installed capacity [MW]
        P_max_HP_3 = tech_infos["P_max_HP_3"] #2
        # Technical minimum [MW]
        Technical_minimum_HP_3= tech_infos["Technical_minimum_HP_3"] #0.1
        P_min_HP_3 = Technical_minimum_HP_3*P_max_HP_3
        # Efficiency [-]
        eta_HP_3 = eta_HP_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_3= tech_infos["Percentage_ramp_up_down_HP_3"] #0.95
        ramp_up_down_HP_3 = Percentage_ramp_up_down_HP_3*P_max_HP_3
        # Variable cost [EUR/MWh]
        cost_var_HP_3 = tech_infos["cost_var_HP_3"] #0.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_3 = Electricity_price

        ## Cogeneration ##

        # Installed capacity [MW]
        P_max_CHP = tech_infos["P_max_CHP"] #0
        # Technical minimum [MW]
        Technical_minimum_CHP= tech_infos["Technical_minimum_CHP"] #0.0
        P_min_CHP = Technical_minimum_CHP*P_max_CHP
        # Efficiency [-]
        eta_CHP_th = tech_infos["eta_CHP_th"] #0.65
        # Power-to-heat ratio [-]
        cb_CHP = tech_infos["cb_CHP"] #0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP= tech_infos["Percentage_ramp_up_down_CHP"] #1
        ramp_up_down_CHP = Percentage_ramp_up_down_CHP*P_max_CHP
        # Variable cost [EUR/MWh]
        cost_var_CHP = tech_infos["cost_var_CHP"] #3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP = tech_infos["cost_fuel_CHP"] #20
        # Price paid to CHP for electrcity production
        CHP_electricity_price=Electricity_price

        ## Cogeneration 2 ##

        # Installed capacity [MW]
        P_max_CHP_2 = tech_infos["P_max_CHP_2"] #0
        # Technical minimum [MW]
        Technical_minimum_CHP_2= tech_infos["Technical_minimum_CHP_2"] #0.4
        P_min_CHP_2 = Technical_minimum_CHP_2*P_max_CHP_2
        # Efficiency [-]
        eta_CHP_th_2 = tech_infos["eta_CHP_th_2"] #0.65
        # Power-to-heat ratio [-]
        cb_CHP_2 = tech_infos["cb_CHP_2"] #0.3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_CHP_2= tech_infos["Percentage_ramp_up_down_CHP_2"] #0.2
        ramp_up_down_CHP_2 = Percentage_ramp_up_down_CHP_2*P_max_CHP_2
        # Variable cost [EUR/MWh]
        cost_var_CHP_2 = tech_infos["cost_var_CHP_2"] #3.9
        # Fuel cost [EUR/MWh]
        cost_fuel_CHP_2 = tech_infos["cost_fuel_CHP_2"] #15
        # Price paid to CHP for electrcity production
        CHP_electricity_price_2=Electricity_price

        ## Absorption Heat pump ##

        # Installed capacity [MW]
        P_max_HP_abs = tech_infos["P_max_HP_abs"] #9
        # Technical minimum [MW]
        Technical_minimum_HP_abs= tech_infos["Technical_minimum_HP_abs"] #0.1
        P_min_HP_abs = Technical_minimum_HP_abs*P_max_HP_abs
        # Efficiency [-]
        eta_HP_abs = tech_infos["eta_HP_abs"] #0.9
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_abs= tech_infos["Percentage_ramp_up_down_HP_abs"] #1
        ramp_up_down_HP_abs = Percentage_ramp_up_down_HP_abs*P_max_HP_abs
        # Variable cost [EUR/MWh]
        cost_var_HP_abs = tech_infos["cost_var_HP_abs"] #0.1

        ## Waste heat - heat exchangers ##

        # Installed capacity [MW]
        P_HEX_capacity = tech_infos["P_HEX_capacity"] #3
        P_max_HEX=P_HEX_capacity*Heat_exchanger_specific
        # Technical minimum [MW]
        P_min_HEX= tech_infos["P_min_HEX"] #0
        # Variable cost [EUR/MWh]
        cost_var_HEX = tech_infos["cost_var_HEX"] #0.25

        ## Waste cooling - heat exchangers ##

        # Installed capacity [MW]
        P_HEX_capacity_cool = tech_infos["P_HEX_capacity_cool"] #2
        P_max_HEX_cool=P_HEX_capacity_cool*Heat_exchanger_specific_cool
        # Technical minimum [MW]
        P_min_HEX_cool= tech_infos["P_min_HEX_cool"] #0
        # Variable cost [EUR/MWh]
        cost_var_HEX_cool = tech_infos["cost_var_HEX_cool"] #0.25

        ## Waste heat - heat pumps - temperature group 1 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_1 = tech_infos["P_max_HP_waste_heat_1"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_1 = tech_infos["Technical_minimum_HP_waste_heat_1"] #0.0
        P_min_HP_waste_heat_1 = Technical_minimum_HP_waste_heat_1*P_max_HP_waste_heat_1
        # Efficiency [-]
        eta_HP_waste_heat_1 = eta_HP_waste_heat_1
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_1= tech_infos["Percentage_ramp_up_down_HP_waste_heat_1"] #0.95
        ramp_up_down_HP_waste_heat_1 = Percentage_ramp_up_down_HP_waste_heat_1*P_max_HP_waste_heat_1
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_1 = tech_infos["cost_var_HP_waste_heat_1"] #0.7
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_1 = Electricity_price

        ## Waste heat - heat pumps - temperature group 2 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_2 = tech_infos["P_max_HP_waste_heat_2"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_2 = tech_infos["Technical_minimum_HP_waste_heat_2"] #0.1
        P_min_HP_waste_heat_2 = Technical_minimum_HP_waste_heat_2*P_max_HP_waste_heat_2
        # Efficiency [-]
        eta_HP_waste_heat_2 = eta_HP_waste_heat_2
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_2= tech_infos["Percentage_ramp_up_down_HP_waste_heat_2"] #0.95
        ramp_up_down_HP_waste_heat_2 = Percentage_ramp_up_down_HP_waste_heat_2*P_max_HP_waste_heat_2
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_2 = tech_infos["cost_var_HP_waste_heat_2"] #1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_2 = Electricity_price

        ## Waste heat - heat pumps - temperature group 3 ##

        # Installed capacity [MW]
        P_max_HP_waste_heat_3 = tech_infos["P_max_HP_waste_heat_3"] #1
        # Technical minimum [MW]
        Technical_minimum_HP_waste_heat_3 = tech_infos["Technical_minimum_HP_waste_heat_3"] #0.1
        P_min_HP_waste_heat_3 = Technical_minimum_HP_waste_heat_3*P_max_HP_waste_heat_3
        # Efficiency [-]
        eta_HP_waste_heat_3 = eta_HP_waste_heat_3
        # Ramp-up-down speed [MW/h]
        Percentage_ramp_up_down_HP_waste_heat_3= tech_infos["Percentage_ramp_up_down_HP_waste_heat_3"] #0.95
        ramp_up_down_HP_waste_heat_3 = Percentage_ramp_up_down_HP_waste_heat_3*P_max_HP_waste_heat_3
        # Variable cost [EUR/MWh]
        cost_var_HP_waste_heat_3 = tech_infos["cost_var_HP_waste_heat_3"] #1.5
        # Fuel cost [EUR/MWh]
        # Choose between fixed electricity price or electricity market price
        cost_fuel_HP_waste_heat_3 = Electricity_price

        ## Thermal energy storage ##

        # Thermal energy storage size [MWh]
        TES_size = tech_infos["TES_size"] #200
        # Minimum percentage of state of charge of thermal energy storage [-]
        SOC_min = tech_infos["SOC_min"] #0.05
        # Thermal energy storage state of charge at the start/end of the year [-]
        TES_start_end = tech_infos["TES_start_end"] #0.05
        # Thermal storage charge/discharge capacity [MW]
        TES_charge_discharge_time = tech_infos["TES_charge_discharge_time"] #5
        TES_in_out_max = TES_size/TES_charge_discharge_time
        TES_loss = tech_infos["TES_loss"] #0.1/(24*7)

        ######## Constraints and variables ########

        ### Variables ###

        # Technologies' operation
        @variable(m, P_max_HOB >= HOB[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_2 >= HOB_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HOB_3 >= HOB_3[1:length(DEM)] >= 0)

        @variable(m, P_max_HP >= HP[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_2 >= HP_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_3 >= HP_3[1:length(DEM)] >= 0)

        @variable(m, P_max_CHP >= CHP[1:length(DEM)] >= 0)
        @variable(m, P_max_CHP_2 >= CHP_2[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_abs >= HP_abs[1:length(DEM)] >= 0)

        @variable(m, HEX[1:length(DEM)] >= 0)

        @variable(m, HEX_cool[1:length(DEM)] >= 0)

        @variable(m, P_max_HP_waste_heat_1 >= HP_waste_heat_1[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_2 >= HP_waste_heat_2[1:length(DEM)] >= 0)
        @variable(m, P_max_HP_waste_heat_3 >= HP_waste_heat_3[1:length(DEM)] >= 0)

        @variable(m, TES_in_out_max >=TES[1:length(DEM)]>= -TES_in_out_max)
        @variable(m, TES_size >= SOC[1:length(DEM)] >= TES_size*SOC_min)

        #@variable(m, maximum(DEM) >= Error[1:length(DEM)] >= 0)
        @variable(m, 0 >= Error[1:length(DEM)] >= 0)

        # Binary variables, on/off

        @variable(m, OnOffHOB[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_2[1:length(DEM)], Bin)
        @variable(m, OnOffHOB_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP[1:length(DEM)], Bin)
        @variable(m, OnOffHP_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_abs[1:length(DEM)], Bin)

        @variable(m, OnOffCHP[1:length(DEM)], Bin)
        @variable(m, OnOffCHP_2[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_1[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_2[1:length(DEM)], Bin)
        @variable(m, OnOffHP_waste_heat_3[1:length(DEM)], Bin)

        @variable(m, OnOffHP_waste_heat_seasonal[1:length(DEM)], Bin)

        ### Constraints ###

        # Supply + TES = DEM
        @constraint(m, HP+HP_2+HP_3+HEX_cool+HP_abs+HP_waste_heat_1+HP_waste_heat_2+HP_waste_heat_3 - TES + Error .== DEM)

        @constraint(m, HP_abs.==(HOB+HOB_2+HOB_3+CHP+CHP_2+HEX)*eta_HP_abs)

        # Thermal storage constraints
        t=2:(length(DEM)-1)
        @constraint(m, SOC[t]-SOC[t .- 1]-TES[t]+SOC[t]*TES_loss .== 0)
        @constraint(m, SOC[length(DEM)-1]+TES[length(DEM)] == TES_size*TES_start_end)
        @constraint(m, SOC[1]-TES[1] == TES_size*TES_start_end)
        @constraint(m, SOC[length(DEM)] == TES_size*TES_start_end)

        # Waste heat constraints - can't utilize more than available
        @constraint(m, HEX .<=P_max_HEX)
        @constraint(m, HEX_cool .<=P_max_HEX_cool)

        @constraint(m, HP_waste_heat_1 .<= Available_waste_heat_heat_pump_1)
        @constraint(m, HP_waste_heat_2 .<= Available_waste_heat_heat_pump_2)
        @constraint(m, HP_waste_heat_3 .<= Available_waste_heat_heat_pump_3)

        # OnOff constraints - technical minimum
        @constraint(m, OnOffHOB*P_min_HOB-HOB.<=0)
        @constraint(m, -OnOffHOB*P_max_HOB+HOB.<=0)
        @constraint(m, OnOffHOB_2*P_min_HOB_2-HOB_2.<=0)
        @constraint(m, -OnOffHOB_2*P_max_HOB_2+HOB_2.<=0)
        @constraint(m, OnOffHOB_3*P_min_HOB_3-HOB_3.<=0)
        @constraint(m, -OnOffHOB_3*P_max_HOB_3+HOB_3.<=0)

        @constraint(m, OnOffHP*P_min_HP-HP.<=0)
        @constraint(m, -OnOffHP*P_max_HP+HP.<=0)
        @constraint(m, OnOffHP_2*P_min_HP_2-HP_2.<=0)
        @constraint(m, -OnOffHP_2*P_max_HP_2+HP_2.<=0)
        @constraint(m, OnOffHP_3*P_min_HP_3-HP_3.<=0)
        @constraint(m, -OnOffHP_3*P_max_HP_3+HP_3.<=0)

        @constraint(m, OnOffHP_abs*P_min_HP_abs-HP_abs.<=0)
        @constraint(m, -OnOffHP_abs*P_max_HP_abs+HP_abs.<=0)

        @constraint(m, OnOffCHP*P_min_CHP-CHP.<=0)
        @constraint(m, -OnOffCHP*P_max_CHP+CHP.<=0)
        @constraint(m, OnOffCHP_2*P_min_CHP_2-CHP_2.<=0)
        @constraint(m, -OnOffCHP_2*P_max_CHP_2+CHP_2.<=0)

        @constraint(m, OnOffHP_waste_heat_1*P_min_HP_waste_heat_1-HP_waste_heat_1.<=0)
        @constraint(m, -OnOffHP_waste_heat_1*P_max_HP_waste_heat_1+HP_waste_heat_1.<=0)
        @constraint(m, OnOffHP_waste_heat_2*P_min_HP_waste_heat_2-HP_waste_heat_2.<=0)
        @constraint(m, -OnOffHP_waste_heat_2*P_max_HP_waste_heat_2+HP_waste_heat_2.<=0)
        @constraint(m, OnOffHP_waste_heat_3*P_min_HP_waste_heat_3-HP_waste_heat_3.<=0)
        @constraint(m, -OnOffHP_waste_heat_3*P_max_HP_waste_heat_3+HP_waste_heat_3.<=0)

        # Ramp-up-down constraint for each technology

        i=2:(length(DEM)-1)

        @constraint(m, HOB[i .- 1]-HOB[i].<=ramp_up_down_HOB*OnOffHOB[i]+(1 .- OnOffHOB[i])*P_min_HOB)
        @constraint(m, HOB[i]-HOB[i .- 1].<=ramp_up_down_HOB*OnOffHOB[i .- 1]+(1 .- OnOffHOB[i .- 1])*P_min_HOB)
        @constraint(m, HOB_2[i .- 1]-HOB_2[i].<=ramp_up_down_HOB_2*OnOffHOB_2[i]+(1 .- OnOffHOB_2[i])*P_min_HOB_2)
        @constraint(m, HOB_2[i]-HOB_2[i .- 1].<=ramp_up_down_HOB_2*OnOffHOB_2[i .- 1]+(1 .- OnOffHOB_2[i .- 1])*P_min_HOB_2)
        @constraint(m, HOB_3[i .- 1]-HOB_3[i].<=ramp_up_down_HOB_3*OnOffHOB_3[i]+(1 .- OnOffHOB_3[i])*P_min_HOB_3)
        @constraint(m, HOB_3[i]-HOB_3[i .- 1].<=ramp_up_down_HOB_3*OnOffHOB_3[i .- 1]+(1 .- OnOffHOB_3[i .- 1])*P_min_HOB_3)

        @constraint(m, HP[i .- 1]-HP[i].<=ramp_up_down_HP*OnOffHP[i]+(1 .- OnOffHP[i])*P_min_HP)
        @constraint(m, HP[i]-HP[i .- 1].<=ramp_up_down_HP*OnOffHP[i .- 1]+(1 .- OnOffHP[i .- 1])*P_min_HP)
        @constraint(m, HP_2[i .- 1]-HP_2[i].<=ramp_up_down_HP_2*OnOffHP_2[i]+(1 .- OnOffHP_2[i])*P_min_HP_2)
        @constraint(m, HP_2[i]-HP_2[i .- 1].<=ramp_up_down_HP_2*OnOffHP_2[i .- 1]+(1 .- OnOffHP_2[i .- 1])*P_min_HP_2)
        @constraint(m, HP_3[i .- 1]-HP_3[i].<=ramp_up_down_HP_3*OnOffHP_3[i]+(1 .- OnOffHP_3[i])*P_min_HP_3)
        @constraint(m, HP_3[i]-HP_3[i .- 1].<=ramp_up_down_HP_3*OnOffHP_3[i .- 1]+(1 .- OnOffHP_3[i .- 1])*P_min_HP_3)

        @constraint(m, HP_abs[i .- 1]-HP_abs[i].<=ramp_up_down_HP_abs*OnOffHP_abs[i]+(1 .- OnOffHP_abs[i])*P_min_HP_abs)
        @constraint(m, HP_abs[i]-HP_abs[i .- 1].<=ramp_up_down_HP_abs*OnOffHP_abs[i .- 1]+(1 .- OnOffHP_abs[i .- 1])*P_min_HP_abs)

        @constraint(m, CHP[i .- 1]-CHP[i].<=ramp_up_down_CHP*OnOffCHP[i]+(1 .- OnOffCHP[i])*P_min_CHP)
        @constraint(m, CHP[i]-CHP[i .- 1].<=ramp_up_down_CHP*OnOffCHP[i .- 1]+(1 .- OnOffCHP[i .- 1])*P_min_CHP)
        @constraint(m, CHP_2[i .- 1]-CHP_2[i].<=ramp_up_down_CHP_2*OnOffCHP_2[i]+(1 .- OnOffCHP_2[i])*P_min_CHP_2)
        @constraint(m, CHP_2[i]-CHP_2[i .- 1].<=ramp_up_down_CHP_2*OnOffCHP_2[i .- 1]+(1 .- OnOffCHP_2[i .- 1])*P_min_CHP_2)

        @constraint(m, HP_waste_heat_1[i .- 1]-HP_waste_heat_1[i].<=ramp_up_down_HP_waste_heat_1*OnOffHP_waste_heat_1[i]+(1 .- OnOffHP_waste_heat_1[i])*P_min_HP_waste_heat_1)
        @constraint(m, HP_waste_heat_1[i]-HP_waste_heat_1[i .- 1].<=ramp_up_down_HP_waste_heat_1*OnOffHP_waste_heat_1[i .- 1]+(1 .- OnOffHP_waste_heat_1[i .- 1])*P_min_HP_waste_heat_1)
        @constraint(m, HP_waste_heat_2[i .- 1]-HP_waste_heat_1[i].<=ramp_up_down_HP_waste_heat_2*OnOffHP_waste_heat_2[i]+(1 .- OnOffHP_waste_heat_2[i])*P_min_HP_waste_heat_2)
        @constraint(m, HP_waste_heat_2[i]-HP_waste_heat_2[i .- 1].<=ramp_up_down_HP_waste_heat_2*OnOffHP_waste_heat_2[i .- 1]+(1 .- OnOffHP_waste_heat_2[i .- 1])*P_min_HP_waste_heat_2)
        @constraint(m, HP_waste_heat_3[i .- 1]-HP_waste_heat_3[i].<=ramp_up_down_HP_waste_heat_3*OnOffHP_waste_heat_3[i]+(1 .- OnOffHP_waste_heat_3[i])*P_min_HP_waste_heat_3)
        @constraint(m, HP_waste_heat_3[i]-HP_waste_heat_3[i .- 1].<=ramp_up_down_HP_waste_heat_3*OnOffHP_waste_heat_3[i .- 1]+(1 .- OnOffHP_waste_heat_3[i .- 1])*P_min_HP_waste_heat_3)

        ###### Objective function ######

        # Minimize total running costs

        @objective(m,Min, sum(Error)*1e10
                          + sum(HOB)*(cost_var_HOB .+ cost_fuel_HOB/eta_HOB)
                          + sum(HOB_2)*(cost_var_HOB_2 .+ cost_fuel_HOB_2/eta_HOB_2)
                          + sum(HOB_3)*(cost_var_HOB_3 .+ cost_fuel_HOB_3/eta_HOB_3)
                          + sum(HP.*(cost_var_HP .+ cost_fuel_HP./eta_HP))
                          + sum(HP_2.*(cost_var_HP_2 .+ cost_fuel_HP_2./eta_HP_2))
                          + sum(HP_3.*(cost_var_HP_3 .+ cost_fuel_HP_3./eta_HP_3))
                          + sum(HP_abs)*cost_var_HP_abs
                          + sum(CHP)*(cost_var_CHP .+ cost_fuel_CHP/eta_CHP_th)
                          + sum(CHP_2)*(cost_var_CHP_2 .+ cost_fuel_CHP_2/eta_CHP_th_2)
                          + sum(HEX)*cost_var_HEX+sum(HEX_cool)*cost_var_HEX_cool
                          + sum(HP_waste_heat_1.*(cost_var_HP_waste_heat_1 .+ cost_fuel_HP_waste_heat_1./eta_HP_waste_heat_1))
                          + sum(HP_waste_heat_2.*(cost_var_HP_waste_heat_2 .+ cost_fuel_HP_waste_heat_2./eta_HP_waste_heat_2))
                          + sum(HP_waste_heat_3.*(cost_var_HP_waste_heat_3 .+ cost_fuel_HP_waste_heat_3./eta_HP_waste_heat_3))
                          - sum(CHP*cb_CHP.*CHP_electricity_price)
                          - sum(CHP_2*cb_CHP_2.*CHP_electricity_price_2))

        optimize!(m)
        status = termination_status(m)
        if primal_status(m) == MOI.NO_SOLUTION
            print_function("No solution found.")
            return
        end

        ##### Result visualization ####

        Result_HOB=value.(HOB)
        Result_HOB_2=value.(HOB_2)
        Result_HOB_3=value.(HOB_3)

        Result_HP=value.(HP)
        Result_HP_2=value.(HP_2)
        Result_HP_3=value.(HP_3)

        Result_HP_abs=value.(HP_abs)

        Result_CHP=value.(CHP)
        Result_CHP_2=value.(CHP_2)

        Result_HEX=value.(HEX)
        Result_HEX_cool=value.(HEX_cool)

        Result_HP_waste_heat_1=value.(HP_waste_heat_1)
        Result_HP_waste_heat_2=value.(HP_waste_heat_2)
        Result_HP_waste_heat_3=value.(HP_waste_heat_3)

        Result_SOC=value.(SOC)

        Result_Error=value.(Error)

        ## Export results ##

        writedlm(string(output_folder, "\\Result_HOB_", network_id, ".csv"), Result_HOB, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_2_", network_id, ".csv"), Result_HOB_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HOB_3_", network_id, ".csv"), Result_HOB_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_", network_id, ".csv"), Result_HP, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_2_", network_id, ".csv"), Result_HP_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_3_", network_id, ".csv"), Result_HP_3, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_abs_", network_id, ".csv"), Result_HP_abs, csv_sep)

        writedlm(string(output_folder, "\\Result_CHP_", network_id, ".csv"), Result_CHP, csv_sep)
        writedlm(string(output_folder, "\\Result_CHP_2_", network_id, ".csv"), Result_CHP_2, csv_sep)

        writedlm(string(output_folder , "\\Result_HEX_", network_id, ".csv"), Result_HEX, csv_sep)

        writedlm(string(output_folder, "\\Result_HEX_cool_", network_id, ".csv"), Result_HEX_cool, csv_sep)

        writedlm(string(output_folder, "\\Result_HP_waste_heat_1_", network_id, ".csv"), Result_HP_waste_heat_1, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_2_", network_id, ".csv"), Result_HP_waste_heat_2, csv_sep)
        writedlm(string(output_folder, "\\Result_HP_waste_heat_3_", network_id, ".csv"), Result_HP_waste_heat_3, csv_sep)

        writedlm(string(output_folder, "\\Result_DEM_", network_id, ".csv"), DEM, csv_sep)

        writedlm(string(output_folder, "\\Result_SOC_", network_id, ".csv"), Result_SOC, csv_sep)

        writedlm(string(output_folder, "\\Result_Error_", network_id, ".csv"), Result_Error, csv_sep)

    end
    export district_cooling_simulation
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
"""
Model to approximate cross products of node voltages

```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc(pm::AbstractWRModel, n::Int, cnd::Int)
    wdc = PowerModels.var(pm, n, cnd, :wdc)
    wdcr = PowerModels.var(pm, n, cnd, :wdcr)

    for (i,j) in PowerModels.ids(pm, n, :buspairsdc)
        JuMP.@constraint(pm.model, wdcr[(i,j)]^2 <= wdc[i]*wdc[j])
        # InfrastructureModels.relaxation_product(pm.model, wdc[i], wdc[j], wdcr_, wdcr_)
    end
end

function constraint_voltage_dc(pm::AbstractWRConicModel, n::Int, cnd::Int)
    wdc = PowerModels.var(pm, n, cnd, :wdc)
    wdcr = PowerModels.var(pm, n, cnd, :wdcr)

    for (i,j) in PowerModels.ids(pm, n, :buspairsdc)
        relaxation_complex_product_conic(pm.model, wdc[i], wdc[j], wdcr[(i,j)])
    end
end

"""
Limits dc branch current

```
p[f_idx] <= wdc[f_bus] * Imax
```
"""
function constraint_dc_branch_current(pm::AbstractWRModel, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    wdc_fr = PowerModels.var(pm, n, cnd, :wdc, f_bus)

    @constraint(pm.model, p_dc_fr <= wdc_fr * ccm_max * p^2)
end


text:
export build_wire_coordinates, get_minimum_intersection_distance

function build_wire_coordinates(wire_directions)
    path = []
    current_x = 0
    current_y = 0
    dx = Dict('L'=> -1, 'R'=> 1, 'U'=> 0, 'D'=> 0)
    dy = Dict('L'=> 0, 'R'=> 0, 'U'=> 1, 'D'=> -1)
    for direction_instruction in wire_directions
        direction = direction_instruction[1]
        count = parse(Int, direction_instruction[2:end])
        step_x = dx[direction]
        step_y = dy[direction]
        for _ = 1:count
            current_x += step_x
            current_y += step_y
            push!(path, (current_x, current_y))
        end
    end
    return path
end

function get_minimum_intersection_distance(intersections)
    minimum_distance = nothing
    for row in eachrow(intersections)
        distance = sum(abs.(row[1]))
        if isnothing(minimum_distance) || minimum_distance > distance
            minimum_distance = distance
        end
    end
    return minimum_distance
end

d03_inputpath = joinpath(@__DIR__, "input.txt")
d03_inputread(path) = readdlm(path, ',')

function d03_answer1(input)
    coords1 = build_wire_coordinates(input[1, :])
    coords2 = build_wire_coordinates(input[2, :])
    intersections = intersect(coords1, coords2)
    return get_minimum_intersection_distance(intersections)
end

function d03_answer2(input)
    return 0
end


text:
module ConjugateGradientTests

using Test, LinearAlgebra, LazyAlgebra, ConjugateGradient

VERBOSE = true
io = (VERBOSE ? stdout : ConjugateGradient.quiet)
dims = (29,11)
nrows, ncols = dims
@testset "least-square fit ($T)" for T in (Float32, Float64)
    # State the problem (least square fit).
    H = rand(T, dims) # direct model
    x = rand(T, ncols) # ground truth
    y = H*x + T(0.01)*rand(T, nrows) # noisy data

    # LHS matrix and RHS vector of the normal equations.
    A = H'*H
    b = H'*y

    # Least-squares solution using Julia linear algebra.
    x0 = A\b

    # Tolerances.
    rtol = 1e-3
    atol = rtol*vnorm2(x0)

    # Solve with conjugate gradient algorithm.
    n = length(b)
    ws = ConjugateGradient.Context(b;
                                   precond = true,
                                   maxiter = n+1,
                                   ftol = (0.0, 0.0),
                                   gtol = (0.0, 0.0),
                                   xtol = (0.0, 0.0))
    x1 = fill!(similar(b), 0)
    VERBOSE && println("\n# Run the conjugate gradient from 0:")
    status = ConjugateGradient.solve!(x1, A, b, ws, io)
    @test status ≥ 0
    VERBOSE && println("rel. err. for x1: ", vnorm2(x1 - x0)/vnorm2(x0))
    @test vnorm2(x1 - x0) ≤ atol

    # Simple diagonal preconditioner.
    q = convert(Vector{T}, diag(A))
    for i in eachindex(q)
        q[i] = 1/q[i]
    end
    M = Diag(q)

    x2 = fill!(similar(b), 0)
    VERBOSE && println("\n# Run the preconditioned conjugate gradient from 0:")
    status = ConjugateGradient.solve!(x2, A, b, M, ws, io)
    @test status ≥ 0
    VERBOSE && println("rel. err. for x1: ", vnorm2(x1 - x0)/vnorm2(x0))
    @test vnorm2(x1 - x0) ≤ atol
end

end # module


text:
using TemporalGPs: posterior_and_lml, predict, predict_marginals

@testset "linear_gaussian_conditionals" begin
    Dlats = [1, 3]
    Dobss = [1, 2]
    # Dlats = [3]
    # Dobss = [2]
    storages = [
        (name="dense storage Float64", val=ArrayStorage(Float64)),
    ]
    Q_types = [
        Val(:dense),
        Val(:diag),
    ]

    @testset "SmallOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))" for
        Dlat in Dlats,
        Dobs in Dobss,
        Q_type in Q_types,
        storage in storages

        println("SmallOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_small_output_lgc(rng, Dlat, Dobs, Q_type, storage.val)

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )

        Q_type == Val(:diag) && @testset "missing data" begin

            # Generate observation with a missing.
            y = conditional_rand(rng, model, rand(rng, x))
            y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
            y_missing .= y
            y_missing[1] = missing

            # Compute logpdf and posterior under model given the missing.
            x_post, lml = posterior_and_lml(x, model, y_missing)

            # Modify the model and compute the logpdf and posterior.
            new_model = SmallOutputLGC(
                model.A[2:end, :], model.a[2:end], model.Q[2:end, 2:end],
            )
            y_new = y[2:end]
            x_post_new, lml_new = posterior_and_lml(x, new_model, y_new)

            # Verify that both things give the same answer.
            @test x_post ≈ x_post_new
            @test lml ≈ lml_new atol=1e-8 rtol=1e-8

            # Check that everything infers and AD gives the right answer.
            @inferred posterior_and_lml(x, model, y_missing)
            x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
            @test x̄[2].Q isa NamedTuple{(:diag, )}
        end
    end

    @testset "LargeOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))" for
        Dlat in Dlats,
        Dobs in Dobss,
        Q_type in Q_types,
        storage in storages

        println("LargeOutputLGC (Dlat=$Dlat, Dobs=$Dobs, Q=$(Q_type), $(storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_large_output_lgc(rng, Dlat, Dobs, Q_type, storage.val)

        @testset "consistency with SmallOutputLGC" begin
            y = rand(rng, TemporalGPs.predict(x, model))
            vanilla_model = TemporalGPs.SmallOutputLGC(model.A, model.a, model.Q)
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y)
            x_large, lml_large = posterior_and_lml(x, model, y)

            @test x_vanilla.m ≈ x_large.m
            @test x_vanilla.P ≈ x_large.P
            @test lml_vanilla ≈ lml_large

            @test predict(x, vanilla_model) ≈ predict(x, model)

            Q_type == Val(:diag) && @testset "missing data" begin

                # Create missing data.
                y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
                y_missing .= y
                y_missing[1] = missing

                # Compute posterior and lml under both SmallOutputLGC and LargeOutputLGC.
                x_post_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_missing)
                x_post_large, lml_large = posterior_and_lml(x, model, y_missing)

                # Check that they give roughly the same answer.
                @test x_post_vanilla ≈ x_post_large
                @test lml_vanilla ≈ lml_large rtol=1e-8 atol=1e-8

                # Check that everything infers and AD gives the right answer.
                @inferred posterior_and_lml(x, model, y_missing)
                x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
                @test x̄[2].Q isa NamedTuple{(:diag, )}
            end
        end

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )
    end

    @testset "ScalarOutputLGC (Dlat=$Dlat, ($storage.name))" for
        Dlat in Dlats,
        storage in [
            (name="dense storage Float64", val=ArrayStorage(Float64)),
            (name="static storage Float64", val=SArrayStorage(Float64)),
        ]

        println("ScalarOutputLGC (Dlat=$Dlat, ($storage.name))")

        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Dlat, storage.val)
        model = random_scalar_output_lgc(rng, Dlat, storage.val)

        @testset "consistency with LGC" begin
            vanilla_model = lgc_from_scalar_output_lgc(model)
            y_vanilla = rand(rng, predict(x, vanilla_model))
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_vanilla)
            x_scalar, lml_scalar = posterior_and_lml(x, model, only(y_vanilla))

            @test x_vanilla.m ≈ x_scalar.m
            @test x_vanilla.P ≈ x_scalar.P
            @test lml_vanilla ≈ lml_scalar
        end

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=storage.val isa SArrayStorage,
        )
    end

    Dmids = [1, 3]

    @testset "BottleneckLGC (Din=$Din, Dmid=$Dmid, Dout=$Dout, Q=$(Q_type))" for
        Din in Dlats,
        Dout in Dobss,
        Dmid in Dmids,
        Q_type in Q_types

        println("BottleneckLGC (Din=$Din, Dmid=$Dmid, Dout=$Dout, Q=$(Q_type))")

        storage = ArrayStorage(Float64)
        rng = MersenneTwister(123456)
        x = random_gaussian(rng, Din, storage)
        model = random_bottleneck_lgc(rng, Din, Dmid, Dout, Q_type, storage)

        @test TemporalGPs.dim_out(model) == Dout
        @test TemporalGPs.dim_in(model) == Din

        test_interface(
            rng, model, x;
            check_adjoints=true,
            check_infers=true,
            check_allocs=false,
        )

        @testset "consistency with SmallOutputLGC" begin
            vanilla_model = small_output_lgc_from_bottleneck(model)

            @test predict(x, vanilla_model) ≈ predict(x, model)
            @test predict_marginals(x, vanilla_model) ≈ predict_marginals(x, model)

            y = rand(rng, predict(x, model))
            x_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y)
            x_bottle, lml_bottle = posterior_and_lml(x, model, y)
            @test x_vanilla.P ≈ x_bottle.P rtol=1e-6
            @test x_vanilla.m ≈ x_bottle.m rtol=1e-6
            @test lml_vanilla ≈ lml_bottle rtol=1e-6

            Q_type == Val(:diag) && @testset "missing data" begin

                # Create missing data.
                y_missing = Vector{Union{Missing, eltype(y)}}(undef, length(y))
                y_missing .= y
                y_missing[1] = missing

                # Compute posterior and lml under both SmallOutputLGC and LargeOutputLGC.
                x_post_vanilla, lml_vanilla = posterior_and_lml(x, vanilla_model, y_missing)
                x_post_large, lml_large = posterior_and_lml(x, model, y_missing)

                # Check that they give roughly the same answer.
                @test x_post_vanilla ≈ x_post_large rtol=1e-8 atol=1e-8
                @test lml_vanilla ≈ lml_large rtol=1e-8 atol=1e-8

                # Check that everything infers and AD gives the right answer.
                @inferred posterior_and_lml(x, model, y_missing)
                x̄ = adjoint_test(posterior_and_lml, (x, model, y_missing))
                @test x̄[2].fan_out.Q isa NamedTuple{(:diag, )}
            end
        end
    end
end


text:
struct Arc
    u::Any
    v::Any
    key::Any
end

Arc(t::Tuple) = Arc(t[1], t[2], t[3])

source(a::Arc) = a.u
target(a::Arc) = a.v
key(a::Arc) = a.key

add_arc!(g::Graph, a::Arc) = add_arc!(g, a.u, a.v, a.key)
rem_arc!(g::Graph, a::Arc) = rem_arc!(g, a.u, a.v, a.key)

function add_arcs_from!(g::Graph, abunch::Vector{Arc})
    keylist = []
    for a in abunch
        res = add_arc!(g, a)
        push!(keylist, res)
    end
    return keylist
end

function rem_arcs_from!(g::Graph, abunch::Vector{Arc})
    removed = 0
    for a in abunch
        remove = try
            rem_arc!(g, a) # if a doesn't comply,
        catch end # just pass.
        if remove
            removed += 1
        end
    end
    return removed
end

function ==(a1::Arc, a2::Arc)
    return source(a1) == source(a2) && target(a1) == target(a2) && key(a1) == key(a2)
end

getindex(g::Graph, a::Arc) = g.fadj[a.u][a.v][a.key]
getindex(g::Graph, u, v, key) = g[Arc(u, v, key)]
getindex(g::Graph, a::Arc, attrib::Symbol) = g[a][attrib]
getindex(g::Graph, u, v, key, attrib::Symbol) = g[Arc(u, v, key)][attrib]

"""
    setindex!(g, val, a, attrib)
setindex! are exclusive for setting arc attributes only.
To delete an attribute, use delete!(g[a], attrib).
"""
function setindex!(g::Graph, val, a::Arc, attrib::Symbol)
    if has_arc(g, a)
        set!(g.fadj[a.u][a.v][a.key], attrib, val)
        set!(g.badj[a.v][a.u][a.key], attrib, val)
        return g[a,attrib]
    else
        return nothing
    end
end

setindex!(g::Graph, val, u, v, key, attrib::Symbol) = setindex!(
    g, val, Arc(u, v, key), attrib
)

"""
    has_arc(g, a)
has_arc checks wether elementary arc a (u,v,key) exists in graph.
"""
function has_arc(g::Graph, a::Arc)
    try
        g[a]
    catch KeyError
        return false
    end
    return true
end

"""
    ArcIter
An iterator of the arcs in the graph.
"""
struct ArcIter
    len::Integer
    ulist::Vector
    vlist::Vector
    keylist::Vector
end

function iterate(ait::ArcIter, state=1)
    if state <= ait.len
        res = Arc(ait.ulist[state], ait.vlist[state], ait.keylist[state])
        state += 1
        return res, state
    else
        return nothing
    end
end

eltype(::ArcIter) = Arc
length(ait::ArcIter) = ait.len
keys(ait::ArcIter) = 1:ait.len

function rand(ait::ArcIter)
    r = rand(1:ait.len)
    res = (ait.ulist[r], ait.vlist[r], ait.keylist[r])
    return res
end

function ==(ait1::ArcIter, ait2::ArcIter)
    return ait1.len == ait2.len && ait1.ulist == ait2.ulist && ait1.vlist ==
    ait2.vlist && ait1.keylist == ait2.keylist
end

"""
    ArcIter(g)
constructs an iterator of all arc in graph.
"""
function ArcIter(g::Graph, ::Colon, ::Colon)
    ulist = []; sizehint!(ulist, g.na)
    vlist = []; sizehint!(vlist, g.na)
    keylist = []; sizehint!(keylist, g.na)
    @inbounds for n in keys(g.fadj), nbrs in keys(g.fadj[n]), key in keys(g.fadj[n][nbrs])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
    end
    return ArcIter(g.na, ulist, vlist, keylist)
end

"""
    ArcIter(g, u)
constructs an iterator of all arc starting from u (single or multiple node).
ArcIter is precise and will raise an error if an element in u is not in graph.
"""
function ArcIter(g::Graph, u, ::Colon)
    len = 0
    ulist = []
    vlist = []
    keylist = []
    @inbounds for n in u, nbrs in keys(g.fadj[n]), key in keys(g.fadj[n][nbrs])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
        len += 1
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    ArcIter(g, :, v)
constructs an iterator of all arc ending at v.
ArcIter is precise and will raise an error if an element in v is not in graph.
"""
function ArcIter(g::Graph, ::Colon, v)
    len = 0
    ulist = []
    vlist = []
    keylist = []
    @inbounds for nbrs in v, n in keys(g.badj[nbrs]), key in keys(g.badj[nbrs][n])
        push!(ulist, n)
        push!(vlist, nbrs)
        push!(keylist, key)
        len += 1
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    ArcIter(g, u, v)
constructs an iterator of all arc inside u-v.
ArcIter is precise and will raise an error if an element in v is not in graph.
"""
function ArcIter(g::Graph, u, v)
    len = length(g.fadj[u][v])
    ulist = []; sizehint!(ulist, len)
    vlist = []; sizehint!(vlist, len)
    keylist = []; sizehint!(keylist, len)
    @inbounds for key in keys(g.fadj[u][v])
        push!(ulist, u)
        push!(vlist, v)
        push!(keylist, key)
    end
    return ArcIter(len, ulist, vlist, keylist)
end

"""
    arcs(g, u, v)
Returns the arcs of a graph g starting from u, ending at v.
The default value of u and v are Colon() to collect every arc.
If u is an iterable container, v must be Colon, vice versa.
"""
arcs(g::Graph, u=:, v=:) = ArcIter(g, u, v)

function getindex(ait::ArcIter, inds)
    if isa(inds, Integer)
        return Arc(ait.ulist[inds], ait.vlist[inds], ait.keylist[inds])
    else
        res = Arc[]
        for i in inds
            push!(res, getindex(ait, i))
        end
        return res
    end
end

firstindex(ait::ArcIter) = 1
lastindex(ait::ArcIter) = ait.len

function has_attrib(g::Graph, attrib::Symbol)
    for a in arcs(g)
        if !(attrib in keys(g[a]))
            return false
        end
    end
    return true
end


text:
facts("Volume Image") do
    fname = joinpath(dirname(@__FILE__), "../../data", "test-3d.dat")
    t = read_VolumeImage(fname)
    t2 = read_VolumeImage(fname)

    context("Reading") do
	t = read_VolumeImage(fname)
	@fact isa(t, EEG.VolumeImage) --> true
    end

    context("Create") do
        # Array
	n = VolumeImage(t.data, t.units, t.x, t.y, t.z, t.t, t.method, t.info, t.coord_system)
	@fact isa(n, EEG.VolumeImage) --> true
        @fact isequal(n.data, t.data) --> true

        # Vector
        d = zeros(length(t.data))
        x = zeros(length(t.data))
        y = zeros(length(t.data))
        z = zeros(length(t.data))
        s = zeros(length(t.data))  # t is already taken
        i = 1
        for xi in 1:length(t.x)
            for yi in 1:length(t.y)
                for zi in 1:length(t.z)
                    for ti in 1:length(t.t)
                        d[i] = float(t.data[xi, yi, zi, ti])
                        x[i] = float(t.x[xi])
                        y[i] = float(t.y[yi])
                        z[i] = float(t.z[zi])
                        s[i] = float(t.t[ti])
                        i += 1
                    end
                end
            end
        end
        n2 = VolumeImage(d, t.units, x, y, z, s, t.method, t.info, t.coord_system)
        @fact isa(n2, EEG.VolumeImage) --> true
        @fact isequal(n2.data, t.data) --> true
        @fact isequal(n2.x, t.x) --> true
        @fact isequal(n2.y, t.y) --> true
        @fact isequal(n2.z, t.z) --> true
        @fact isequal(n2.t, t.t) --> true

    end

    context("Maths") do

	@fact isequal(t, t2) --> true

        @fact t + t --> t * 2
	t4 = t * 4
	@fact t4 / 2 --> t + t
	@fact t4 / 2 --> t4 - t - t

	t2 = read_VolumeImage(fname)
	t2.units = "A/m^3"
	@fact_throws t + t2

    end

    context("Min/Max") do

	@fact maximum(t) --> 33.2692985535
	@fact minimum(t) --> -7.5189352036

	b = t * 2

	@fact b.data[1, 1, 1] --> 2 * b.data[1, 1, 1]
	@fact b.data[2, 2, 2] --> 2 * b.data[2, 2, 2]

	c = mean([b, t])

	@fact size(c.data) --> size(b.data)
	@fact c.data[1, 1, 1] --> (b.data[1, 1, 1] + t.data[1, 1, 1]) / 2

	@fact maximum([b, t]) --> maximum(t) * 2
	@fact minimum([b, t]) --> minimum(b)

	n = normalise([b, t])

	@fact maximum(n) --> less_than_or_equal(1.0)
	@fact minimum(n) --> greater_than_or_equal(-1.0)

    end

    context("Dimension checks") do
    	t2 = deepcopy(t); t2.x = t2.x[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = deepcopy(t); t2.y = t2.y[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = deepcopy(t); t2.z = t2.z[1:3]
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    	t2 = read_VolumeImage(joinpath(dirname(@__FILE__), "../../data", "test-4d.dat"))
    	@fact_throws KeyError EEG.dimensions_equal(t, t2)
    end

    context("Average") do
    	s = mean(t)
    end

    context("Find dipoles") do

    	t = read_VolumeImage(fname)
    	dips = find_dipoles(mean(t))
    	@fact size(dips) --> (9,)

    	dips = EEG.new_dipole_method(mean(t))
    	@fact size(dips) --> (9,)

    	fname = joinpath(dirname(@__FILE__), "../../data", "test-4d.dat")
    	t2 = read_VolumeImage(fname)
    	dips = find_dipoles(mean(t2))
    	@fact size(dips) --> (3,)

    	# Test dipoles are returned in order of size
    	@fact issorted([dip.size for dip in dips], rev = true) --> true
    end

    context("Best Dipole") do

        dips = find_dipoles(mean(t))
        # Plots.pyplot(size=(1400, 400))
        # p = plot(t, c = :inferno)
        # p = plot(p, Talairach(-0.04, 0.01, 0.02))
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips)
        #p = plot(p, Talairach(-0.05, 0, 0.01), c = :red)

        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        # Take closest
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips, maxdist = 0.0015)
        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        # Take only valid dipole
        dists = [euclidean(Talairach(-0.05, 0, 0.01), dip) for dip=dips]
        bd = best_dipole(Talairach(-0.05, 0, 0.01), dips, maxdist = 0.015)
        @fact float(bd.x) --> roughly(-0.0525, atol =  0.001)
        @fact float(bd.y) --> roughly(-0.00378, atol =  0.001)
        @fact float(bd.z) --> roughly(0.0168099, atol =  0.001)

        bd = best_dipole(Talairach(-0.05, 0, 0.01), Dipole[])
        @fact isnan(bd) --> true

    end

    context("Show") do
    	show(t)
    	show(normalise(t))
    	t.info["Regularisation"] = 1.2
    	show(t)
    end

    context("Plotting") do
    	EEG.plot(mean(t))
    	EEG.plot(mean(t), min_val = 0, max_val = 50)
        EEG.plot(mean(t), elp = joinpath(dirname(@__FILE__), "../../data", "test.elp"))
    	p = EEG.plot(mean(t), threshold = 24 )

        context("Overlay dipole") do

            dips = find_dipoles(mean(t))
            p = plot(p, dips)
            p = plot(p, dips[1], c = :red)

        end
    end
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    abstract type AbstractBridge <: MOI.Bridges.AbstractType

Subtype of [`MOI.Bridges.AbstractBridge`](@ref) for constraint
bridges.
"""
abstract type AbstractBridge <: MOI.Bridges.AbstractBridge end

"""
    MOI.supports_constraint(
        BT::Type{<:AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet},
    )::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging
`F`-in-`S` constraints.

## Implementation notes

 * This method depends only on the type of the inputs, not the runtime values.
 * There is a default fallback, so you need only implement this method for
   constraint types that the bridge implements.
"""
function MOI.supports_constraint(
    ::Type{<:AbstractBridge},
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

"""
    concrete_bridge_type(
        BT::Type{<:AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet}
    )::Type

Return the concrete type of the bridge supporting `F`-in-`S` constraints.

This function can only be called if `MOI.supports_constraint(BT, F, S)` is
`true`.

## Example

The [`SplitIntervalBridge`](@ref) bridges a [`MOI.VariableIndex`](@ref)-in-[`MOI.Interval`](@ref)
constraint into a [`MOI.VariableIndex`](@ref)-in-[`MOI.GreaterThan`](@ref) and a
[`MOI.VariableIndex`](@ref)-in-[`MOI.LessThan`](@ref) constraint.

```jldoctest; setup=:(using MathOptInterface; const MOI = MathOptInterface)
julia> MOI.Bridges.Constraint.concrete_bridge_type(
           MOI.Bridges.Constraint.SplitIntervalBridge{Float64},
           MOI.VariableIndex,
           MOI.Interval{Float64},
       )
MathOptInterface.Bridges.Constraint.SplitIntervalBridge{Float64, MathOptInterface.VariableIndex, MathOptInterface.Interval{Float64}, MathOptInterface.GreaterThan{Float64}, MathOptInterface.LessThan{Float64}}
```
"""
function concrete_bridge_type(
    ::Type{BT},
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
) where {BT}
    return BT
end

function concrete_bridge_type(
    b::MOI.Bridges.AbstractBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return concrete_bridge_type(MOI.Bridges.bridge_type(b, F, S), F, S)
end

"""
    bridge_constraint(
        BT::Type{<:AbstractBridge},
        model::MOI.ModelLike,
        func::AbstractFunction,
        set::MOI.AbstractSet,
    )::BT

Bridge the constraint `func`-in-`set` using bridge `BT` to `model` and returns
a bridge object of type `BT`.

## Implementation notes

 * The bridge type `BT` should be a concrete type, that is, all the type
   parameters of the bridge must be set.
"""
function bridge_constraint end


text:
struct GradConformity <: Conformity end
const H1Conformity = GradConformity

function get_face_own_nodes(reffe::GenericLagrangianRefFE{GradConformity},conf::GradConformity)
  p = get_polytope(reffe)
  orders = get_orders(reffe)
  nodes, face_own_nodes = compute_nodes(p,orders)
  face_own_nodes
end

function get_own_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity},conf::GradConformity)
  p = get_polytope(reffe)
  face_own_nodes = get_face_own_nodes(reffe)
  dofs = get_dof_basis(reffe)
  interior_nodes = dofs.nodes[face_own_nodes[end]]
  compute_own_nodes_permutations(p,interior_nodes)
end

function get_face_own_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity},conf::GradConformity)
  own_nodes_permutations = get_own_nodes_permutations(reffe)
  reffaces = reffe.reffe.metadata
  _reffaces = vcat(reffaces...)
  face_own_nodes_permutations = map(get_own_nodes_permutations,_reffaces)

  _compute_face_own_nodes_permutations(
    num_nodes(reffe),
    get_face_own_nodes(reffe),
    face_own_nodes_permutations,
    own_nodes_permutations)
end

function _compute_face_own_nodes_permutations(
  nnodes,
  face_own_nodes,
  face_own_nodes_permutations,
  own_nodes_permutations)

  if length(face_own_nodes_permutations) == 0
    # Vertex degenerated case
    return [own_nodes_permutations,]
  else
    if nnodes == length(face_own_nodes[end])
      # 0-order degenerated case
      _face_own_nodes_permutations = map( (x) -> fill(Int[],length(x)) , face_own_nodes_permutations )
    else
      # Standard case
      _face_own_nodes_permutations = copy(face_own_nodes_permutations)
    end
    push!(_face_own_nodes_permutations,own_nodes_permutations)
    return _face_own_nodes_permutations
  end
end

# API particular to LagrangianRefFE{GradConformity}

"""
    ReferenceFE{N}(reffe::GenericLagrangianRefFE{GradConformity},iface::Integer) where N
"""
function ReferenceFE{N}(reffe::GenericLagrangianRefFE{GradConformity},iface::Integer) where N
  reffaces = reffe.reffe.metadata
  reffaces[N+1][iface]
end

function ReferenceFE{D}(reffe::GenericLagrangianRefFE{GradConformity,D},iface::Integer) where D
  @assert iface==1 "Only one D-face"
  reffe
end

"""
    get_reffaces(
      ::Type{ReferenceFE{d}},
      reffe::GenericLagrangianRefFE{GradConformity}) where d -> Vector{GenericLagrangianRefFE{GradConformity,M,d}}
"""
function get_reffaces(::Type{ReferenceFE{d}},reffe::GenericLagrangianRefFE{GradConformity}) where d
  ftype_to_reffe, _ = _compute_reffes_and_face_types(reffe,Val{d}())
  ftype_to_reffe
  [reffe for reffe in ftype_to_reffe]
end

"""
    get_face_type(reffe::GenericLagrangianRefFE{GradConformity}, d::Integer) -> Vector{Int}
"""
function get_face_type(reffe::GenericLagrangianRefFE{GradConformity}, d::Integer)
  _, iface_to_ftype = _compute_reffes_and_face_types(reffe,Val{d}())
  iface_to_ftype
end

function _compute_reffes_and_face_types(reffe::GenericLagrangianRefFE{GradConformity},::Val{d}) where d
  p = get_polytope(reffe)
  iface_to_reffe = [ ReferenceFE{d}(reffe,iface) for iface in 1:num_faces(p,d) ]
  _find_unique_with_indices(iface_to_reffe)
end

"""
    is_first_order(reffe::GenericLagrangianRefFE{GradConformity}) -> Bool
"""
function is_first_order(reffe::GenericLagrangianRefFE{GradConformity})
  p = get_polytope(reffe)
  r = true
  r = r && num_vertices(p) == num_nodes(reffe)
  r = r && get_vertex_node(reffe) == collect(1:num_nodes(reffe))
  r
end

"""
    is_affine(reffe::GenericLagrangianRefFE{GradConformity}) -> Bool

Query if the `reffe` leads to an afine map
(true only for first order spaces on top of simplices)
"""
function is_affine(reffe::GenericLagrangianRefFE{GradConformity})
  p = get_polytope(reffe)
  is_first_order(reffe) && is_simplex(p)
end

"""
    is_P(reffe::GenericLagrangianRefFE{GradConformity})
"""
function is_P(reffe::GenericLagrangianRefFE{GradConformity})
  is_simplex(get_polytope(reffe))
end

"""
   is_Q(reffe::GenericLagrangianRefFE{GradConformity})
"""
function is_Q(reffe::GenericLagrangianRefFE{GradConformity})
  monomials = get_prebasis(reffe)
  n = length(get_exponents(monomials))
  is_n_cube(get_polytope(reffe)) && (prod(get_orders(reffe).+1) == n)
end

"""
   is_S(reffe::GenericLagrangianRefFE{GradConformity})
"""
function is_S(reffe::GenericLagrangianRefFE{GradConformity})
  is_n_cube(get_polytope(reffe)) && ! is_Q(reffe)
end

function to_dict(reffe::GenericLagrangianRefFE{GradConformity})
  p = get_polytope(reffe)
  b = get_prebasis(reffe)
  dict = Dict{Symbol,Any}()
  dict[:orders] = collect(get_orders(reffe))
  dict[:extrusion] = Array(TensorValues.get_array(get_extrusion(p)))
  if is_S(reffe)
    dict[:space] = "serendipity"
  else
    dict[:space] = "default"
  end
  dict[:value] = string(get_value_type(b))
  dict
end

function from_dict(::Type{<:LagrangianRefFE},dict::Dict{Symbol,Any})
  orders = Tuple(dict[:orders])
  extrusion = Tuple(dict[:extrusion])
  if dict[:value] == "Float64"
    value = Float64
  else
    @notimplemented
  end
  space = dict[:space]
  p = Polytope(extrusion...)
  if space == "default"
    reffe = LagrangianRefFE(value,p,orders)
  elseif space == "serendipity"
    reffe = SerendipityRefFE(value,p,orders)
  else
    @unreachable "unknown space type"
  end
  reffe
end

# Construction of LagrangianRefFE from Polytopes

"""
    LagrangianRefFE(::Type{T},p::Polytope,orders) where T
    LagrangianRefFE(::Type{T},p::Polytope,order::Int) where T

Builds a `LagrangianRefFE` object on top of the given polytope. `T` is the type of
the value of the approximation space (e.g., `T=Float64` for scalar-valued problems,
`T=VectorValue{N,Float64}` for vector-valued problems with `N` components). The arguments `order` or `orders`
are for the polynomial order of the resulting space, which allows isotropic or anisotropic orders respectively
(provided that the cell topology allows the given anisotropic order). The argument `orders` should be an
indexable collection of `D` integers (e.g., a tuple or a vector), being `D` the number of space dimensions.

In order to be able to use this function, the type of the provided polytope `p` has to implement the
following additional methods. They have been implemented for `ExtrusionPolytope` in the library. They
need to be implemented for new polytope types in order to build Lagangian reference elements on top of them.

- [`compute_monomial_basis(::Type{T},p::Polytope,orders) where T`](@ref)
- [`compute_own_nodes(p::Polytope,orders)`](@ref)
- [`compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)`](@ref)

The following methods are also used in the construction of the `LagrangianRefFE` object. A default implementation
of them is available in terms of the three previous methods. However, the user can also implement them for
new polytope types increasing customization possibilities.

- [`compute_nodes(p::Polytope,orders)`](@ref)
- [`compute_own_nodes_permutations(p::Polytope, interior_nodes)`](@ref)
- [`compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T`](@ref)
"""
function LagrangianRefFE(::Type{T},p::Polytope{D},orders;space::Symbol=_default_space(p)) where {T,D}
  if space == :P && is_n_cube(p)
    return _PDiscRefFE(T,p,orders)
  elseif space == :S && is_n_cube(p)
    SerendipityRefFE(T,p,orders)
  else
    if any(orders.==0) && !all(orders.==0)
      cont = map(i -> i == 0 ? DISC : CONT,orders)
      return _cd_lagrangian_ref_fe(T,p,orders,cont)
    else
      return _lagrangian_ref_fe(T,p,orders)
    end
  end
end

function _default_space(p)
  if is_n_cube(p)
    :Q
  else
    :P
  end
end

function _lagrangian_ref_fe(::Type{T},p::Polytope{D},orders) where {T,D}

  prebasis = compute_monomial_basis(T,p,orders)
  nodes, face_own_nodes = compute_nodes(p,orders)
  dofs = LagrangianDofBasis(T,nodes)
  reffaces = compute_lagrangian_reffaces(T,p,orders)

  nnodes = length(dofs.nodes)
  ndofs = length(dofs.dof_to_node)
  metadata = reffaces
  _reffaces = vcat(reffaces...)
  face_nodes = _generate_face_nodes(nnodes,face_own_nodes,p,_reffaces)
  face_own_dofs = _generate_face_own_dofs(face_own_nodes, dofs.node_and_comp_to_dof)
  face_dofs = _generate_face_dofs(ndofs,face_own_dofs,p,_reffaces)

  if all(orders .== 0 ) && D>0
    conf = L2Conformity()
  else
    conf = GradConformity()
  end

  reffe = GenericRefFE(
    ndofs,
    p,
    prebasis,
    dofs,
    conf,
    metadata,
    face_dofs)

  GenericLagrangianRefFE(reffe,face_nodes)

end

function MonomialBasis(::Type{T},p::Polytope,orders) where T
  compute_monomial_basis(T,p,orders)
end

function LagrangianDofBasis(::Type{T},p::Polytope,orders) where T
  nodes, _ = compute_nodes(p,orders)
  LagrangianDofBasis(T,nodes)
end

# Helpers for LagrangianRefFE constructor

function _generate_face_nodes(nnodes,face_to_own_nodes,polytope,reffaces)

    face_to_num_fnodes = map(num_nodes,reffaces)
    push!(face_to_num_fnodes,nnodes)

    face_to_lface_to_own_fnodes = map(get_face_own_nodes,reffaces)
    push!(face_to_lface_to_own_fnodes,face_to_own_nodes)

    face_to_lface_to_face = get_faces(polytope)

  _generate_face_nodes_aux(
    nnodes,
    face_to_own_nodes,
    face_to_num_fnodes,
    face_to_lface_to_own_fnodes,
    face_to_lface_to_face)
end

function _generate_face_dofs(ndofs,face_to_own_dofs,polytope,reffaces)

    face_to_num_fdofs = map(num_dofs,reffaces)
    push!(face_to_num_fdofs,ndofs)

    face_to_lface_to_own_fdofs = map(get_face_own_dofs,reffaces)
    push!(face_to_lface_to_own_fdofs,face_to_own_dofs)

    face_to_lface_to_face = get_faces(polytope)

  _generate_face_nodes_aux(
    ndofs,
    face_to_own_dofs,
    face_to_num_fdofs,
    face_to_lface_to_own_fdofs,
    face_to_lface_to_face)
end

function _generate_face_nodes_aux(
  nnodes,
  face_to_own_nodes,
  face_to_num_fnodes,
  face_to_lface_to_own_fnodes,
  face_to_lface_to_face)

  if nnodes == length(face_to_own_nodes[end])
    face_fnode_to_node = fill(Int[],length(face_to_own_nodes))
    face_fnode_to_node[end] = collect(1:nnodes)
    return face_fnode_to_node
  end

  face_fnode_to_node = Vector{Int}[]
  for (face, nfnodes) in enumerate(face_to_num_fnodes)
    fnode_to_node = zeros(Int,nfnodes)
    lface_to_face = face_to_lface_to_face[face]
    lface_to_own_fnodes = face_to_lface_to_own_fnodes[face]
    for (lface, faceto) in enumerate(lface_to_face)
      own_nodes = face_to_own_nodes[faceto]
      own_fnodes = lface_to_own_fnodes[lface]
      fnode_to_node[own_fnodes] = own_nodes
    end
    push!(face_fnode_to_node,fnode_to_node)
  end

  face_fnode_to_node
end

# Constructors taking Int

function LagrangianRefFE(::Type{T},p::Polytope{D},order::Int;space::Symbol=_default_space(p)) where {T,D}
  orders = tfill(order,Val{D}())
  LagrangianRefFE(T,p,orders;space=space)
end

function MonomialBasis(::Type{T},p::Polytope{D},order::Int) where {D,T}
  orders = tfill(order,Val{D}())
  MonomialBasis(T,p,orders)
end

function LagrangianDofBasis(::Type{T},p::Polytope{D},order::Int) where {T,D}
  orders = tfill(order,Val{D}())
  LagrangianDofBasis(T,p,orders)
end

# Queries needed to be implemented for polytopes in order to use them
# for building LagrangianRefFEs in a seamless way

"""
    compute_monomial_basis(::Type{T},p::Polytope,orders) where T -> MonomialBasis

Returns the monomial basis of value type `T` and order per direction described by `orders`
on top of the polytope `p`.
"""
function compute_monomial_basis(::Type{T},p::Polytope,orders) where T
  @abstractmethod
end

"""
    compute_own_nodes(p::Polytope{D},orders) where D -> Vector{Point{D,Float64}}

Returns the coordinates of the nodes owned by the interior of the polytope
associated with a Lagrangian space with the order per direction described by `orders`.
"""
function compute_own_nodes(p::Polytope,orders)
  @abstractmethod
end

"""
    compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)

Returns a vector or a tuple with the order per direction at the face `face`
of the polytope `p` when restricting the order per direction `orders` to this face.
`iface` is the face id of `face` in the numeration restricted to the face dimension.
"""
function compute_face_orders(p::Polytope,face::Polytope,iface::Int,orders)
  @abstractmethod
end

"""
    compute_nodes(p::Polytope,orders)

When called

    node_coords, face_own_nodes = compute_nodes(p,orders)

Returns `node_coords`, the nodal coordinates of all the Lagrangian nodes associated with the order per direction
`orders`, and `face_own_nodes`, being a vector of vectors indicating which nodes are owned by each of
the faces of the polytope `p`.
"""
function compute_nodes(p::Polytope,orders)
  _compute_nodes(p,orders)
end

"""
    compute_own_nodes_permutations(
      p::Polytope, own_nodes_coordinates) -> Vector{Vector{Int}}

Returns a vector of vectors with the permutations of the nodes owned by the interior of the
polytope.
"""
function compute_own_nodes_permutations(p::Polytope, interior_nodes)
  perms = _compute_node_permutations(p, interior_nodes)
  perms
end

"""
    compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T

Returns a tuple of length `D` being the number of space dimensions.
The entry `d+1` of this tuple contains a vector of `LagrangianRefFE`
one for each face of dimension `d` on the boundary of the polytope.
"""
function compute_lagrangian_reffaces(::Type{T},p::Polytope,orders) where T
  _compute_lagrangian_reffaces(T,p,orders)
end

# Default implementations

function _compute_nodes(p,orders)
  if any( orders .== 0)
    _compute_constant_nodes(p,orders)
  elseif all(orders .== 1)
    _compute_linear_nodes(p)
  else
    _compute_high_order_nodes(p,orders)
  end
end

function _compute_constant_nodes(p,orders)
  @assert all( orders .== 0) "If an order is 0 in some direction, it should be 0 also in the others"
  x = compute_own_nodes(p,orders)
  facenodes = [Int[] for i in 1:num_faces(p)]
  push!(facenodes[end],1)
  x, facenodes
end

function _compute_linear_nodes(p)
  x = get_vertex_coordinates(p)
  facenodes = [Int[] for i in 1:num_faces(p)]
  for i in 1:num_vertices(p)
    push!(facenodes[i],i)
  end
  x, facenodes
end

function _compute_high_order_nodes(p::Polytope{D},orders) where D
  nodes = Point{D,Float64}[]
  facenodes = [Int[] for i in 1:num_faces(p)]
  _compute_high_order_nodes_dim_0!(nodes,facenodes,p)
  for d in 1:(num_dims(p)-1)
    _compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,Val{d}())
  end
  _compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  (nodes, facenodes)
end

function _compute_high_order_nodes_dim_0!(nodes,facenodes,p)
  x = get_vertex_coordinates(p)
  k = 1
  for vertex in 1:num_vertices(p)
    push!(nodes,x[vertex])
    push!(facenodes[vertex],k)
    k += 1
  end
end

@noinline function _compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,::Val{d}) where d
  x = get_vertex_coordinates(p)
  offset = get_offset(p,d)
  k = length(nodes)+1
  for iface in 1:num_faces(p,d)
    face = Polytope{d}(p,iface)
    face_ref_x = get_vertex_coordinates(face)
    face_prebasis = MonomialBasis(Float64,face,1)
    change = inv(evaluate(face_prebasis,face_ref_x))
    face_shapefuns = change_basis(face_prebasis,change)
    face_vertex_ids = get_faces(p,d,0)[iface]
    face_x = x[face_vertex_ids]
    face_orders = compute_face_orders(p,face,iface,orders)
    face_interior_nodes = compute_own_nodes(face,face_orders)
    face_high_x = evaluate(face_shapefuns,face_interior_nodes)*face_x
    for xi in 1:length(face_high_x)
      push!(nodes,face_high_x[xi])
      push!(facenodes[iface+offset],k)
      k += 1
    end
  end
end

function _compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  k = length(nodes)+1
  p_high_x = compute_own_nodes(p,orders)
  for xi in 1:length(p_high_x)
    push!(nodes,p_high_x[xi])
    push!(facenodes[end],k)
    k += 1
  end
end

_compute_node_permutations(::Polytope{0}, interior_nodes) = [[1]]

function _compute_node_permutations(p, interior_nodes)
  vertex_to_coord = get_vertex_coordinates(p)
  lbasis = MonomialBasis(Float64,p,1)
  change = inv(evaluate(lbasis,vertex_to_coord))
  lshapefuns = change_basis(lbasis,change)
  perms = get_vertex_permutations(p)
  map = evaluate(lshapefuns,interior_nodes)
  pvertex_to_coord = similar(vertex_to_coord)
  node_perms = Vector{Int}[]
  tol = 1.0e-10
  for vertex_to_pvertex in perms
    node_to_pnode = fill(INVALID_PERM,length(interior_nodes))
    pvertex_to_coord[vertex_to_pvertex] = vertex_to_coord
    pinterior_nodes = map*pvertex_to_coord
    for node in 1:length(interior_nodes)
      x = interior_nodes[node]
      pnode = findfirst(i->norm(i-x)<tol,pinterior_nodes)
      if pnode != nothing
         node_to_pnode[node] = pnode
      end
    end
    push!(node_perms,node_to_pnode)
  end
  node_perms
end

_compute_lagrangian_reffaces(::Type{T},p::Polytope{0},orders) where T = ()

function _compute_lagrangian_reffaces(::Type{T},p::Polytope{D},orders) where {T,D}
  reffaces = [ LagrangianRefFE{d}[]  for d in 0:D ]
  p0 = Polytope{0}(p,1)
  reffe0 = LagrangianRefFE(T,p0,())
  for vertex in 1:num_vertices(p)
    push!(reffaces[0+1],reffe0)
  end
  offsets = get_offsets(p)
  for d in 1:(num_dims(p)-1)
    offset = offsets[d+1]
    for iface in 1:num_faces(p,d)
      face = Polytope{d}(p,iface)
      face_orders = compute_face_orders(p,face,iface,orders)
      refface = LagrangianRefFE(T,face,face_orders)
      push!(reffaces[d+1],refface)
    end
  end
  tuple(reffaces...)
end

# Particular implementation for ExtrusionPolytope

function LagrangianRefFE(p::ExtrusionPolytope)
  order = 1
  LagrangianRefFE(Float64,p,order)
end

function compute_monomial_basis(::Type{T},p::ExtrusionPolytope{D},orders) where {D,T}
  extrusion = Tuple(p.extrusion)
  terms = _monomial_terms(extrusion,orders)
  MonomialBasis{D}(T,orders,terms)
end

function compute_own_nodes(p::ExtrusionPolytope{D},orders) where D
  extrusion = Tuple(p.extrusion)
  if all(orders .== 0)
    _interior_nodes_order_0(p)
  else
    _interior_nodes(extrusion,orders)
  end
end

function _interior_nodes_order_0(p)
  x = get_vertex_coordinates(p)
  x0 = sum(x) / length(x)
  [x0,]
end

function compute_face_orders(p::ExtrusionPolytope,face::ExtrusionPolytope{D},iface::Int,orders) where D
  d = num_dims(face)
  offset = get_offset(p,d)
  nface = p.dface.nfaces[iface+offset]
  face_orders = _eliminate_zeros(Val{D}(),nface.extrusion,orders)
  Tuple(face_orders)
end

function _eliminate_zeros(::Val{d},a,o) where d
  b = zero(mutable(Point{d,Int}))
  D = num_components(a)
  k = 1
  for i in 1:D
    m = a[i]
    if (m != 0)
      b[k] = o[i]
      k += 1
    end
  end
  Point(b)
end

function compute_nodes(p::ExtrusionPolytope{D},orders) where D
  _nodes, facenodes = _compute_nodes(p,orders)
  if any( orders .== 0)
    return (_nodes, facenodes)
  end
  terms = _coords_to_terms(_nodes,orders)
  nodes = _terms_to_coords(terms,orders)
  (nodes, facenodes)
end

# Helpers for the ExtrusionPolytope-related implementation

function _monomial_terms(extrusion::NTuple{D,Int},orders) where D
  terms = CartesianIndex{D}[]
  if D == 0
    push!(terms,CartesianIndex(()))
    return terms
  end
  _check_orders(extrusion,orders)
  M = mutable(VectorValue{D,Int})
  term = zero(M)
  _orders = M(orders)
  k = 0
  _add_terms!(terms,term,extrusion,_orders,D,k)
  terms
end

function _interior_nodes(extrusion::NTuple{D,Int},orders) where D
  _check_orders(extrusion,orders)
  terms = CartesianIndex{D}[]
  M = mutable(VectorValue{D,Int})
  term = zero(M)
  _orders = M(orders)
  k = 1
  _add_terms!(terms,term,extrusion,_orders,D,k)
  _terms_to_coords(terms,orders)
end

function _check_orders(extrusion,orders)
  D = length(extrusion)
  @assert length(orders) == D "container of orders not long enough"
  _orders = collect(orders)
  if extrusion[D] == HEX_AXIS
    _orders[D] = 0
  end
  for d in (D-1):-1:1
    if (extrusion[d] == HEX_AXIS || d == 1) && _orders[d+1] == 0
      _orders[d] = 0
    end
  end
  nz = _orders[_orders .!= 0]
  if length(nz) > 1
    @assert all(nz .== nz[1]) "The provided anisotropic order is not compatible with polytope topology"
  end
  nothing
end

function _add_terms!(terms,term,extrusion,orders,dim,k)
  _term = copy(term)
  _orders = copy(orders)
  indexbase = 1
  for i in k:(_orders[dim]-k)
    _term[dim] = i + indexbase
    if dim > 1
      if (extrusion[dim] == TET_AXIS) && i != 0
        _orders .-= 1
      end
      _add_terms!(terms,_term,extrusion,_orders,dim-1,k)
    else
      push!(terms,CartesianIndex(Tuple(_term)))
    end
  end
end

function _coords_to_terms(coords::Vector{<:Point{D}},orders) where D
  indexbase = 1
  terms = CartesianIndex{D}[]
  P = Point{D,Int}
  t = zero(mutable(P))
  for x in coords
    for d in 1:D
      t[d] = round(x[d]*orders[d]) + indexbase
    end
    term = CartesianIndex(Tuple(t))
    push!(terms,term)
  end
  terms
end

function  _terms_to_coords(terms::Vector{CartesianIndex{D}},orders) where D
  P = Point{D,Float64}
  indexbase = 1
  nodes = P[]
  x = zero(mutable(P))
  for t in terms
    for d in 1:D
      x[d] = (t[d] - indexbase) / orders[d]
    end
    node = P(x)
    push!(nodes,node)
  end
  nodes
end

function _extract_nonzeros(mask,values)
  b = Int[]
  for (m,n) in zip(mask,values)
    if (m != 0)
      push!(b, n)
    end
  end
  return Tuple(b)
end

# Precomputed instances

"""
    const VERTEX1 = LagrangianRefFE(Float64,VERTEX,1)
"""
const VERTEX1 = LagrangianRefFE(Float64,VERTEX,1)

"""
    const SEG2 = LagrangianRefFE(Float64,SEGMENT,1)
"""
const SEG2 = LagrangianRefFE(Float64,SEGMENT,1)

"""
    const QUAD4 = LagrangianRefFE(Float64,QUAD,1)
"""
const QUAD4 = LagrangianRefFE(Float64,QUAD,1)

"""
    const TRI3 = LagrangianRefFE(Float64,TRI,1)
"""
const TRI3 = LagrangianRefFE(Float64,TRI,1)

"""
    const TET4 = LagrangianRefFE(Float64,TET,1)
"""
const TET4 = LagrangianRefFE(Float64,TET,1)

"""
    const HEX8 = LagrangianRefFE(Float64,HEX,1)
"""
const HEX8 = LagrangianRefFE(Float64,HEX,1)


text:
using PlotlyJS, LightGraphs
import GraphPlot  # for spring_layout

# Generate a random layout
G =  LightGraphs.euclidean_graph(200, 2, cutoff=0.125)[1]
# Position nodes
pos_x, pos_y = GraphPlot.spring_layout(G)

# Create plot points
edge_x = []
edge_y = []

for edge in edges(G)
    push!(edge_x, pos_x[src(edge)])
    push!(edge_x, pos_x[dst(edge)])
    push!(edge_y, pos_y[src(edge)])
    push!(edge_y, pos_y[dst(edge)])
end

#  Color node points by the number of connections.
color_map = [size(neighbors(G, node))[1] for node in 1:200]

# Create edges
edges_trace = scatter(
    mode="lines",
    x=edge_x,
    y=edge_y,
    line=attr(
        width=0.5,
        color="#888"
    ),
)

# Create nodes
nodes_trace = scatter(
    x=pos_x,
    y=pos_y,
    mode="markers",
    text = [string("# of connections: ", connection) for connection in color_map],
    marker=attr(
        showscale=true,
        colorscale=colors.imola,
        color=color_map,
        size=10,
        colorbar=attr(
            thickness=15,
            title="Node Connections",
            xanchor="left",
            titleside="right"
      )
    )
)

# Create Plot
plot(
    [edges_trace, nodes_trace],
    Layout(
        hovermode="closest",
        title="Network Graph made with Julia",
        titlefont_size=16,
        showlegend=false,
        showarrow=false,
        xaxis=attr(showgrid=false, zeroline=false, showticklabels=false),
        yaxis=attr(showgrid=false, zeroline=false, showticklabels=false)
    )
)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module Gaussian1d
using DataModel
using StatsBase
export Gaussian1dModel
type Gaussian1dModel <: AbstractDataModel
    sdTh::Float64
    sdX::Float64
    x::Array{Float64,1}
    postvar::Float64
    postmean::Float64
    function Gaussian1dModel(  nsdTh::Float64,       nsdX::Float64,       nx::Array{Float64,1})

      n=length(nx)
      postvar=1.0/(1.0/nsdTh^2+n/nsdX^2)
      postmean=sum(nx[1:n]/(n+ nsdX^2/nsdTh^2))
      new(nsdTh,nsdX,nx,postvar,postmean)

    end

end

function posterior(dm::Gaussian1dModel,nobs::Int64)
  x=dm.x
  sdX=dm.sdX
  sdTh=dm.sdTh
  postvar=1.0/(1.0/sdTh^2+nobs/sdX^2)
  postmean=sum(x[1:nobs]/(nobs+ sdX^2/sdTh^2))
  return (postmean,postvar)
end

function DataModel.getllik(dm::Gaussian1dModel;nobs=length(dm.x)) ## add nobs fercitily

    function logdensity(beta)
      return Float64[-0.5*(beta[1]-dm.postmean)^2/dm.postvar ]
    end
    return logdensity
end

function DataModel.getgrad(dm::Gaussian1dModel,batchsize;nobs=length(dm.x))
    x=dm.x
    sdX=dm.sdX
    sdTh=dm.sdTh

    postvar=1.0/(1.0/sdTh^2+nobs/sdX^2)
    postmean=sum(x[1:nobs]/(nobs+ sdX^2/sdTh^2))

    B2=sum(x[1:nobs])/sdX^2
    function grad_log_posterior(beta)
      return Float64[-beta[1]/postvar+B2]
    end
    function grad_log_posterior_sub(beta)

        chosen_indices=sample(1:nobs,batchsize,replace=false)

        return Float64[-beta[1]/postvar+nobs/batchsize*sum(x[chosen_indices])/sdX^2]
    end
    if nobs==batchsize
      return grad_log_posterior
    else
      return grad_log_posterior_sub
    end
end
function DataModel.getfullgrad(dm::Gaussian1dModel;nobs=length(dm.x))
  return DataModel.getgrad(dm,nobs,nobs=nobs)
end

function artificial_Gaussian1d(nObs,seed::Int64 =1)
    srand(seed)
    sdTh=1.0
    sdX=1.0
    th=randn()
    x=th+sdX*randn(nObs)
    srand()
    Gaussian1dModel(sdTh,sdX,x)
end
end