text:
using LinearAlgebra

struct InfoBox
  context::Context
  x::Float64
  y::Float64
  scaling::Float64

  ds_vec::RenderObject
  u1_vec::RenderObject
  u2_vec::RenderObject
  c_vec::RenderObject

  id_text::RenderObject
  s_text::RenderObject
  ds_text::RenderObject
  p_text::RenderObject
  u1_text::RenderObject
  u2_text::RenderObject

  box::RenderObject
end

function InfoBox(context::Context, x::Float64, y::Float64, scaling::Float64)
  ds_vec = vis.make_vector_xy(context, 0.0, 0.0, 0.0, 0.0, [1.0, 0.0, 0.0])
  u1_vec = vis.make_vector_xy(context, 0.0, 0.0, 0.0, 0.0, [1.0, 1.0, 0.0])
  u2_vec = vis.make_vector_xy(context, 0.0, 0.0, 0.0, 0.0, [1.0, 1.0, 0.0])
  c_vec = vis.make_vector_xy(context, 0.0, 0.0, 0.0, 0.0, [1.0, 1.0, 1.0])

  id_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)
  s_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)
  ds_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)
  p_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)
  u1_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)
  u2_text = vis.make_text(context, "test", 0.0, 0.0, 0.5)

  position = RenderData(fill(GLfloat(0.0), 2 * 8), 2, GL_DYNAMIC_DRAW)
  color = RenderData(fill(GLfloat(0.0), 3 * 8), 3, GL_STATIC_DRAW)
  usetext = RenderData(fill(GLfloat(0.0), 1 * 8), 1, GL_STATIC_DRAW)
  texcoord = RenderData(fill(GLfloat(0.0), 2 * 8), 2, GL_STATIC_DRAW)
  elnb = 8

  box = RenderObject(context, [position, color, usetext, texcoord], elnb)

  return InfoBox(context, x, y, scaling, ds_vec, u1_vec, u2_vec, c_vec,
                 id_text, s_text, ds_text, p_text, u1_text, u2_text, box)
end

function render(info::InfoBox)
  xmin = Inf
  xmax = -Inf
  ymin = Inf
  ymax = -Inf
  all_text = [info.id_text, info.s_text, info.ds_text, info.p_text,
              info.u1_text, info.u2_text]
  for text in all_text
    (_xmin, _ymin, _xmax, _ymax) = bounding_box(text)
    xmin = min(xmin, _xmin)
    xmax = max(xmax, _xmax)
    ymin = min(ymin, _ymin)
    ymax = max(ymax, _ymax)

    render(text)
  end
  render(info.ds_vec)
  render(info.u1_vec)
  render(info.u2_vec)
  render(info.c_vec)
  position = GLfloat[xmin, ymin,
                     xmin, ymax,

                     xmin, ymax,
                     xmax, ymax,

                     xmax, ymax,
                     xmax, ymin,

                     xmax, ymin,
                     xmin, ymin]
  update_buffer!(info.box, position, info.context.attributes[1])
  render(info.box)
end

function make_road(context::Context, x::AbstractArray, y::AbstractArray,
                   t::Number)
  @assert length(x) == length(y) > 1
  len = length(x)

  # figure out the gradient
  dpx = similar(x)
  dpy = similar(y)

  for i in 1:len
    upidx = i == len ? len : i + 1
    dwidx = i == 1 ? 1 : i - 1

    dx = x[upidx] - x[dwidx]
    dy = y[upidx] - y[dwidx]

    n = sqrt(dx^2 + dy^2)

    dpx[i] = -dy / n
    dpy[i] = dx / n
  end

  # allocate CPU memory for vertices objects
  points = fill(GLfloat(0), 2 * 2 * len)
  color = fill(GLfloat(0), 2 * 3 * len)
  usetex = fill(GLfloat(0), 2 * len)
  texcoord = fill(GLfloat(0), 2 * 2 * len)
  indices = fill(GLuint(0), 2 * 3 * (len - 1))

  # create points and colors
  for i in 1:len
    points[2 * 2 * (i - 1) + 1] = x[i] + dpx[i] * t / 2
    points[2 * 2 * (i - 1) + 2] = y[i] + dpy[i] * t / 2

    points[2 * 2 * (i - 1) + 3] = x[i] - dpx[i] * t / 2
    points[2 * 2 * (i - 1) + 4] = y[i] - dpy[i] * t / 2

    color[2 * 3 * (i - 1) + 1] = 0.15
    color[2 * 3 * (i - 1) + 2] = 0.15
    color[2 * 3 * (i - 1) + 3] = 0.15

    color[2 * 3 * (i - 1) + 4] = color[2 * 3 * (i - 1) + 1]
    color[2 * 3 * (i - 1) + 5] = color[2 * 3 * (i - 1) + 2]
    color[2 * 3 * (i - 1) + 6] = color[2 * 3 * (i - 1) + 3]
  end

  # select points to form triangles
  for i in 1:(len - 1)
    indices[2 * 3 * (i - 1) + 1] = 2 * (i - 1) + 1 - 1
    indices[2 * 3 * (i - 1) + 2] = 2 * (i - 1) + 2 - 1
    indices[2 * 3 * (i - 1) + 3] = 2 * (i) + 1 - 1

    indices[2 * 3 * (i - 1) + 4] = 2 * (i) + 1 - 1
    indices[2 * 3 * (i - 1) + 5] = 2 * (i) + 2 - 1
    indices[2 * 3 * (i - 1) + 6] = 2 * (i - 1) + 2 - 1
  end

  return RenderObject(context, [RenderData(points, 2, GL_STATIC_DRAW),
                                RenderData(color, 3, GL_STATIC_DRAW),
                                RenderData(usetex, 1, GL_STATIC_DRAW),
                                RenderData(texcoord, 2, GL_STATIC_DRAW)],
                      indices)
end

# dodger blue
function make_car(context::Context,
                  car_color::Array{<: Number, 1}=GLfloat[0.12, 0.56, 1.0])
  width = 1.7
  length = 4.5

  flpos = 0.6 * width / 2
  flw = 0.3 * width
  fll = 0.3 * width

  blpos = 0.6 * width / 2
  blw = 0.3 * width
  bll = 0.3 * width
  position_data = GLfloat[
                          # car
                          -width / 2, -length / 2,
                          -width / 2, length / 2,
                          width / 2, length / 2,
                          width / 2, -length / 2,
                          # left front light
                          -flpos - flw / 2, length / 2 - fll / 2,
                          -flpos - flw / 2, length / 2 + fll / 2,
                          -flpos + flw / 2, length / 2 + fll / 2,
                          -flpos + flw / 2, length / 2 - fll / 2,
                          # right front light
                          flpos - flw / 2, length / 2 - fll / 2,
                          flpos - flw / 2, length / 2 + fll / 2,
                          flpos + flw / 2, length / 2 + fll / 2,
                          flpos + flw / 2, length / 2 - fll / 2,
                          # left front light
                          -blpos - blw / 2, -length / 2 - bll / 2,
                          -blpos - blw / 2, -length / 2 + bll / 2,
                          -blpos + blw / 2, -length / 2 + bll / 2,
                          -blpos + blw / 2, -length / 2 - bll / 2,
                          # right front light
                          blpos - blw / 2, -length / 2 - bll / 2,
                          blpos - blw / 2, -length / 2 + bll / 2,
                          blpos + blw / 2, -length / 2 + bll / 2,
                          blpos + blw / 2, -length / 2 - bll / 2
                         ]
  position = RenderData(position_data, 2, GL_STATIC_DRAW)

  fl_color = GLfloat[1.0, 1.0, 1.0]
  bl_color = GLfloat[1.0, 0.0, 0.0]
  color_data = [repeat(Array{GLfloat}(car_color), 4);
                repeat(fl_color, 8); repeat(bl_color, 8)]
  color = RenderData(color_data, 3, GL_STATIC_DRAW)

  usetex = RenderData(fill(GLfloat(0), 20), 1, GL_STATIC_DRAW)
  texcoord = RenderData(fill(GLfloat(0), 2 * 20), 2, GL_STATIC_DRAW)

  idx1 = GLuint[0, 1, 2,
                0, 2, 3]
  idx = GLuint[]
  offset = GLuint(0)
  for i in 1:5
    append!(idx, idx1 .+ offset)
    offset += 4
  end

  return RenderObject(context, [position, color, usetex, texcoord], idx)
end

function car_lights!(car::RenderObject, on::Union{Bool, Nothing}=nothing)
  @assert car.render_buffers[1].size == 4 * 4 * 2 * 5
  if on == nothing
    car.elnb = (car.elnb == 3 * 2 * 3) ? 3 * 2 * 5 : 3 * 2 * 3 # toggle
  else
    if on == true
      car.elnb = 3 * 2 * 5 # turn on
    else
      car.elnb = 3 * 2 * 3 # turn off
    end
  end
end

# Vectors #####################################################################
function make_vector_thr(context::Context, x::Float64, y::Float64,
                         th::Float64, r::Float64,
                         color::AbstractArray{<: Real}=GLfloat[1.0, 0.0, 0.0])
  th = th - pi / 2

  arrow_r = 0.015
  points = GLfloat[0, 0,
                   0, r,
                   -arrow_r, r - arrow_r * 2,
                   0, r,
                   arrow_r, r - arrow_r * 2,
                   0, r]
  rotM = [cos(th) -sin(th);
          sin(th) cos(th)]
  for i in 1:2:length(points)
    points[i:i+1] = rotM * points[i:i+1]
    points[i] += x
    points[i+1] += y
  end
  color = repeat(Array{GLfloat}(color), 6)

  usetex = fill(GLfloat(0), 6)
  texcoord = fill(GLfloat(0), 2 * 6)

  return RenderObject(context, [RenderData(points, 2, GL_DYNAMIC_DRAW),
                                RenderData(color, 3, GL_STATIC_DRAW),
                                RenderData(usetex, 1, GL_STATIC_DRAW),
                                RenderData(texcoord, 2, GL_STATIC_DRAW)], 6)
end

function make_vector_xy(context::Context, x1::Float64, y1::Float64,
                        x2::Float64, y2::Float64,
                        color::AbstractArray{<: Real}=GLfloat[1.0, 0.0, 0.0])
  dx = x2 - x1
  dy = y2 - y1
  th = atan(dy, dx)
  r = sqrt(dx^2 + dy^2)

  return make_vector_thr(context, x1, y1, th, r, color)
end

function update_vector_thr!(v::RenderObject,
                            x::Float64, y::Float64, th::Float64, r::Float64)
  th = th - pi / 2

  if r < 0.0
    r = -r
    th += pi
  end

  arrow_r = 0.015
  points = GLfloat[0, 0,
                   0, r,
                   -arrow_r, r - arrow_r * 2,
                   0, r,
                   arrow_r, r - arrow_r * 2,
                   0, r]
  rotM = [cos(th) -sin(th);
          sin(th) cos(th)]
  for i in 1:2:length(points)
    points[i:i+1] = rotM * points[i:i+1]
    points[i] += x
    points[i+1] += y
  end

  update_buffer!(v, points, v.context.attributes[1])
end

function update_vector_xy!(v::RenderObject,
                           x1::Float64, y1::Float64, x2::Float64, y2::Float64)
  dx = x2 - x1
  dy = y2 - y1
  r = sqrt(dx^2 + dy^2)
  th = atan(dy, dx)

  update_vector_thr!(v, x1, y1, th, r)
end


text:
include("Graph.jl")

using LinearAlgebra
using SparseArrays

function Degree(G)
    C = zeros(G.n)
    for (ID, u, v, w) in G.E
        C[u] += w
        C[v] += w
    end
    return C
end

function spA(G)
    Is = zeros(Int, G.m*2)
    Js = zeros(Int, G.m*2)
    Vs = zeros(G.m*2)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = w
        Is[G.m + ID] = v
        Js[G.m + ID] = u
        Vs[G.m + ID] = w
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function alpP(G, alp)
    d = Degree(G)
    Is = zeros(Int, G.m*2)
    Js = zeros(Int, G.m*2)
    Vs = zeros(G.m*2)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = alp * w / d[v]
        Is[G.m + ID] = v
        Js[G.m + ID] = u
        Vs[G.m + ID] = alp * w / d[u]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function BetweennessCentrality(G)
    g = getAdjacentList(G)
    C = zeros(G.n)
    p = Array{Array{Int32, 1}, 1}(undef, G.n)
    d = zeros(Int32, G.n)
    S = zeros(Int32, G.n+10)
    sigma = zeros(G.n)
    Q = zeros(Int32, G.n+10)
    delta = zeros(G.n)
    for s = 1 : G.n
        foreach(i -> p[i] = [], 1 : G.n)
        top = 0
        sigma .= 0
        sigma[s] = 1.0
        d .= -1
        d[s] = 0
        front = 1
        rear = 1
        Q[1] = s

        while front <= rear
            v = Q[front]
            front += 1
            top += 1
            S[top] = v
            for w in g[v]
                if d[w] < 0
                    rear += 1
                    Q[rear] = w
                    d[w] = d[v] + 1
                end
                if d[w] == (d[v] + 1)
                    sigma[w] += sigma[v]
                    push!(p[w], v)
                end
            end
        end

        delta .= 0

        while top > 0
            w = S[top]
            top -= 1
            for v in p[w]
                delta[v] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
                if w != s
                    C[w] += delta[w]
                end
            end
        end

    end

    return C
end

function selectPivots(G, eps)
    deg = Degree(G)
    nl = zeros(Int, G.n)
    foreach(i -> nl[i] = i, 1 : G.n)
    sort!(nl, by = x->-deg[x]);
    pivots = []
    kk = round(Int, log2(G.n)/eps^2)
    foreach(i -> push!(pivots, nl[i]), 1 : kk)
    return pivots
end

function ApproxBetweennessCentrality(G; eps = 0.5)
    g = getAdjacentList(G)
    pivots = selectPivots(G, eps)
    C = zeros(G.n)
    p = Array{Array{Int32, 1}, 1}(undef, G.n)
    d = zeros(Int32, G.n)
    S = zeros(Int32, G.n+10)
    sigma = zeros(G.n)
    Q = zeros(Int32, G.n+10)
    delta = zeros(G.n)
    for s in pivots
        foreach(i -> p[i] = [], 1 : G.n)
        top = 0
        sigma .= 0
        sigma[s] = 1.0
        d .= -1
        d[s] = 0
        front = 1
        rear = 1
        Q[1] = s

        while front <= rear
            v = Q[front]
            front += 1
            top += 1
            S[top] = v
            for w in g[v]
                if d[w] < 0
                    rear += 1
                    Q[rear] = w
                    d[w] = d[v] + 1
                end
                if d[w] == (d[v] + 1)
                    sigma[w] += sigma[v]
                    push!(p[w], v)
                end
            end
        end

        delta .= 0

        while top > 0
            w = S[top]
            top -= 1
            for v in p[w]
                delta[v] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
                if w != s
                    C[w] += delta[w]
                end
            end
        end

    end

    return C
end

function ClosenessCentrality(G)
    g = getAdjacentList(G)
    C = zeros(G.n)
    d = zeros(Int32, G.n)
    Q = zeros(Int32, G.n+10)
    for s = 1 : G.n
        d .= -1
        d[s] = 0
        front = 1
        rear = 1
        Q[1] = s

        while front <= rear
            v = Q[front]
            front += 1
            for w in g[v]
                if d[w] < 0
                    rear += 1
                    Q[rear] = w
                    d[w] = d[v] + 1
                end
            end
        end

        C[s] = sum(d)
    end

    foreach(i -> C[i] = 1.0 / C[i], 1 : G.n)

    return C
end

function ForestDistanceClosenessCentrality(G)
    # get Forest Matrix
    IpL = zeros(G.n, G.n)
    for (ID, u, v, w) in G.E
        IpL[u, u] += w
        IpL[v, v] += w
        IpL[u, v] -= w
        IpL[v, u] -= w
    end
    foreach(i -> IpL[i, i] += 1, 1 : G.n)
    W = inv(IpL)

    # calculate FDC
    trace = tr(W)
    C = zeros(G.n)
    foreach(i -> C[i] = G.n / (G.n * W[i, i] + trace - 2), 1 : G.n)
    #foreach(i -> C[i] = W[i, i], 1 : G.n)

    return C
end

function InformationCentrality(G)
    # get pseudo inverse of L
    L = zeros(G.n, G.n)
    for (ID, u, v, w) in G.E
        L[u, u] += w
        L[v, v] += w
        L[u, v] -= w
        L[v, u] -= w
    end
    L .+= (1 / G.n)
    Lp = inv(L)
    Lp .-= (1 / G.n)

    # calculate Information Centrality
    trace = tr(Lp)
    C = zeros(G.n)
    foreach(i -> C[i] = G.n / (G.n * Lp[i, i] + trace), 1 : G.n)
    #foreach(i -> C[i] = Lp[i, i], 1 : G.n)

    return C
end

function PageRank(G; alpha = 0.85)
    aP = alpP(G, alpha)
    C = zeros(G.n)
    C[1] = 1.0
    adC = ((1.0 - alpha) / G.n) * ones(G.n)

    while true
        pC = copy(C)
        C = aP * C + adC
        if norm(C - pC) < 1e-12
            break
        end
    end

    return C
end

function LeaderRank(G)
    n = G.n + 1
    g = Array{Array{Int32, 1}, 1}(undef, n)
    foreach(i -> g[i] = [], 1 : n)
    for (ID, u, v, w) in G.E
        push!(g[u], v)
        push!(g[v], u)
    end
    for i = 1 : n-1
        push!(g[n], i)
        push!(g[i], n)
    end
    d = zeros(Int32, n)
    foreach(i -> d[i] = size(g[i], 1), 1 : n)

    C = zeros(n)
    err = zeros(n)
    pC = zeros(n)
    C[1] = 1.0

    while true
        foreach(i -> pC[i] = C[i], 1 : n)
        for u = 1 : n
            C[u] = 0
            for v in g[u]
                C[u] += (pC[v] / d[v])
            end
        end
        foreach(i -> err[i] = abs(C[i] - pC[i]), 1 : n)
        maxErr = err[argmax(err)]
        if maxErr < 1e-9
            break
        end
    end

    ansC = zeros(G.n)
    foreach(i -> ansC[i] = C[i], 1 : G.n)

    return ansC
end

function LeverageCentrality(G)
    d = Degree(G)
    C = zeros(G.n)
    for (ID, u, v, w) in G.E
        C[u] += ((d[u] - d[v]) / (d[u] + d[v]))
        C[v] += ((d[v] - d[u]) / (d[u] + d[v]))
    end
    foreach(i -> C[i] /= d[i], 1 : G.n)
    return C
end

function EigenvectorCentrality(G)
    C = zeros(G.n)
    C[1] = 1.0
    A = spA(G)

    while true
        pC = copy(C)
        C = A * C
        C ./= C[argmax(C)]
        if norm(C - pC) < 1e-9
            break
        end
    end

    return C
end

function KatzCentrality(G; fac = 0.85)
    A = spA(G)
    getmax(X) = X[argmax(X)]
    x = ones(G.n)
    xc = ones(G.n)
    tt = 0.0
    lmd = 0.0
    while true
        xc = copy(x)
        x = A * x
        tmp = getmax(x)
        x ./= tmp
        if abs(tmp-tt)<1e-6
            lmd = tmp
            break
        end
        tt = tmp
    end
    alp = fac / lmd
    C = zeros(G.n)
    addC = ones(G.n)
    while getmax(addC) > 1e-9
        addC = A * addC
        addC .*= alp
        C = C + addC
    end

    return C
end

function SpanningTreeCentrality(G)
    g = Array{Array{Int32, 1}, 1}(undef, G.n)
    foreach(i -> g[i] = [], 1 : G.n)
    for (ID, u, v, w) in G.E
        push!(g[u], v)
        push!(g[v], u)
    end
    d = zeros(Int32, G.n)
    foreach(i -> d[i] = size(g[i], 1), 1 : G.n)

    L = zeros(G.n-1, G.n-1)
    for (ID, u, v, w) in G.E
        if u < G.n
            L[u, u] += w
        end
        if v < G.n
            L[v, v] += w
        end
        if (u < G.n) && (v < G.n)
            L[u, v] -= w
            L[v, u] -= w
        end
    end
    nL = inv(L)

    I(N) = Diagonal(ones(N))

    C = zeros(G.n)

    for u = 1 : G.n
        if d[u] == 1
            C[u] = 0.0
        else
            U = spzeros(d[u]-1, G.n-1)
            for i = 1 : d[u]-1
                v = g[u][i]
                if u == G.n
                    U[i, v] = 1
                elseif v == G.n
                    U[i, u] = 1
                else
                    U[i, u] = 1
                    U[i, v] = -1
                end
            end
            C[u] = 1.0 - d[u] * det(I(d[u]-1) - U * nL * U')
        end
    end

    return C
end

function STnumber(G)
    L = zeros(G.n-1, G.n-1)
    for (ID, u, v, w) in G.E
        L[u, u] += w
        if v < G.n
            L[v, v] += w
            L[u, v] -= w
            L[v, u] -= w
        end
    end
    return det(L)
end

function reduced(G, sv)
    V = zeros(Int32, G.n-1)
    ori = zeros(Int32, G.n)
    for i = 1 : sv-1
        V[i] = i
        ori[i] = i
    end
    for i = sv+1 : G.n
        V[i-1] = i
        ori[i] = i-1
    end
    edge = Set{Tuple{Int32, Int32}}()
    for (ID, u, v, w) in G.E
        if (u != sv) && (v != sv)
            push!(edge, (ori[u], ori[v]))
        end
    end
    m = length(edge)
    E = Array{Tuple{Int32, Int32, Int32, Float64}, 1}(undef, m)

    ID = 0
    n = G.n-1
    for (u, v) in edge
        ID = ID + 1
        E[ID] = (ID, u, v, 1.0)
    end

    return Graph(n, m, V, E)
end

function STCdp(G)
    C = zeros(G.n)
    all = STnumber(G)
    d = Degree(G)
    for i = 1 : G.n
        C[i] = (all - d[i] * STnumber(reduced(G, i))) / all
    end
    return C
end

ext: jl


text:
using graphtheoryufrj
using Base.Test

include("ratton_first_assignment.jl")
include("ratton_second_assignment.jl")
"""
First assignment: Basics of a graph package (Modeling, BFS, DFS...)
"""
FLAG_RATTON_ASSIGNMENT01 = false
infile_name = "../assets/as_graph.txt"
outfile_name = "parents_q4.txt"
initial_vertex = [1 5 10 50 100 500 1000 5000 10000 25000]

if FLAG_RATTON_ASSIGNMENT01
	# Test the memory representation of the structures of the graph.
	@test test_one(infile_name) == "TEST_ONE_SUCESS"

	# Time to perform 10 bfs on the graph.
	@test test_two(infile_name, initial_vertex) == "TEST_TWO_SUCESS"

	# Time to perform 10 dfs on the graph.
	@test test_three(infile_name, initial_vertex) == "TEST_THREE_SUCESS"

	# Find the parents of given nodes.
	@test test_four(infile_name, outfile_name) == "TEST_FOUR_SUCESS"

	#
	@test test_five(infile_name) == "TEST_FIVE_SUCESS"

	# Empirical distribuiton test.
	@test test_six(infile_name) == "TEST_SIX_SUCESS"

	#
	@test test_seven(infile_name) == "TEST_SEVEN_SUCESS"
end

"""
First assignment: Extending the graph package (Dijkstra, Floyd-Warshall...)
"""
FLAG_RATTON_ASSIGNMENT02 = false
algorithm = prim_mst
initial_vertex_mst = 20
infile_name = ["../assets/grafo_1.txt",
	"../assets/grafo_2.txt",
	"../assets/grafo_3.txt",
	"../assets/grafo_4.txt",
	"../assets/grafo_5.txt"]

if FLAG_RATTON_ASSIGNMENT02
	for i in infile_name
		#
		@test test_one(i) == "TEST_ONE_SUCESS"

		#
		@test test_two(i, algorithm, initial_vertex_mst) == "TEST_TWO_SUCESS"

		#
		@test test_three(i) == "TEST_THREE_SUCESS"

	end

	#
	@test test_colab_one() == "TEST_COLAB_ONE_SUCESS"

	#
	@test test_colab_two() == "TEST_COLAB_TWO_SUCESS"

end


text:
##############
# StateLabel #
##############
struct StateLabel{N}
    label::NTuple{N}
    hash::UInt64
    StateLabel(label::NTuple{N}) where N = new{N}(label, hash(label))
    StateLabel(label::NTuple{N}, h::UInt64) where N = new{N}(label, h)
end

StateLabel(s::StateLabel) = s
# StateLabel{N}(label::NTuple{N}) = StateLabel{N}(label)
StateLabel(i...) = StateLabel(i)

nfactors(::Type{StateLabel{N}}) where N = N
Base.getindex(s::StateLabel, i) = s.label[i]
Base.getindex(s::StateLabel, arr::AbstractArray) = s.label[arr]

Base.copy(s::StateLabel{N}) where N = StateLabel{N}(s.label, s.hash)
Base.hash(s::StateLabel) = s.hash
Base.hash(s::StateLabel, h::UInt64) = hash(hash(s), h)

Base.:(==)(a::StateLabel{N},b::StateLabel{N}) where N = hash(a) == hash(b)

Base.iterate(s::StateLabel) = iterate(s.label)
Base.iterate(s::StateLabel, i) = iterate(s.label, i)
# Base.done(s::StateLabel, i) = done(s.label, i)

Base.first(s::StateLabel) = first(s.label)
Base.last(s::StateLabel) = last(s.label)
# Base.endof(s::StateLabel) = endof(s.label)

Base.length(::StateLabel{N}) where N = N

Base.map(f::Union{Function,DataType}, s::StateLabel) = StateLabel(map(f, s.label))

tensor(a::StateLabel, b::StateLabel) = StateLabel(tuple(a.label..., b.label...))

labelstr(s::StateLabel) = join(map(repr, s.label), ',')
Base.repr(s::StateLabel) = repr(typeof(s)) * "(" * labelstr(s) * ")"
Base.show(io::IO, s::StateLabel) = print(io, repr(s))

##############
# OpLabel #
##############
struct OpLabel{N}
    k::StateLabel{N}
    b::StateLabel{N}
    hash::UInt64
    OpLabel(k::StateLabel{N}, b::StateLabel{N}) where N = new{N}(k, b, hash(k, hash(b)))
    OpLabel(k::StateLabel{N}, b::StateLabel{N}, h::UInt64) where N = new{N}(k, b, h)
end

OpLabel(op::OpLabel) = op
OpLabel(::StateLabel, ::StateLabel) = error("OpLabel can only be constructed if both StateLabels have the same number of factors")
# OpLabel(k::StateLabel{N}, b::StateLabel{N}) where N = OpLabel{N}(k, b)
OpLabel(k, b) = OpLabel(StateLabel(k), StateLabel(b))

klabel(o::OpLabel) = o.k
blabel(o::OpLabel) = o.b

Base.copy(o::OpLabel{N}) where N = OpLabel{N}(o.k, o.b, o.hash)
Base.hash(o::OpLabel) = o.hash
Base.hash(o::OpLabel, h::UInt64) = hash(hash(o), h)

Base.:(==)(a::OpLabel{N}, b::OpLabel{N})  where N = hash(a) == hash(b)

Base.length(::OpLabel{N}) where N = N

Base.adjoint(o::OpLabel) = OpLabel(o.b, o.k)
tensor(o1::OpLabel, o2::OpLabel) = OpLabel(tensor(o1.k, o2.k), tensor(o1.b, o2.b))

ptranspose(k::StateLabel{N}, b::StateLabel{N}, i) where N = OpLabel{N}(setindex(k, b[i], i), setindex(b, k[i], i))
ptranspose(o::OpLabel, i) = ptranspose(o.k, o.b, i)
ptranspose_dual(o::OpLabel, i) = ptranspose(o.b, o.k, i)

traceout(k::StateLabel, b::StateLabel, i) = OpLabel(except(k, i), except(b, i))
traceout(o::OpLabel, i) = traceout(o.k, o.b, i)
traceout_dual(o::OpLabel, i) = traceout(o.b, o.k, i)

Base.repr(o::OpLabel) = repr(typeof(o)) * "(" * ktstr(o.k) * "," * brstr(o.b) * ")"
Base.show(io::IO, o::OpLabel) = print(io, repr(o))

####################
# Helper Functions #
####################
is_sum_x(o::OpLabel, x) = sum(klabel(o))==sum(blabel(o))==x
is_sum_x(s::StateLabel, x) = sum(s) == x

ctpair(k,v) = (k', v')
nzcoeff(k,v) = v!=0
second(t) = t[2]
except(label::StateLabel, i) = StateLabel(label[1:i-1]..., label[i+1:end]...)
setindex(label::StateLabel, x, y) = StateLabel(label[1:y-1]..., x, label[y+1:end]...)

permute(label::StateLabel{N}, perm::Vector) where N = StateLabel{N}(label[perm])
permute(o::OpLabel, perm::Vector) = OpLabel(permute(o.k, perm), permute(o.b, perm))

function switch!(arr, i, j)
    tmp = arr[i]
    arr[i] = arr[j]
    arr[j] = tmp
    return arr
end
switch(label::StateLabel{N}, i, j) where N =  StateLabel{N}(label[switch!([1:N], i, j)])

switch(o::OpLabel, i, j) = OpLabel(switch(o.k, i, j), switch(o.b, i, j))

export StateLabel,
    OpLabel,
    klabel,
    blabel


text:
using ExactWrightFisher

@testset "test the wf likelihood functions with array storage for precomputed coeffs" begin

    function simulate_WF3_data()
        K = 3
        α = ones(K)
        sα = sum(α)
        Pop_size_WF3 = 15
        Ntimes_WF3 = 3
        time_step_WF3 = 0.1
        time_grid_WF3 = [k*time_step_WF3 for k in 0:(Ntimes_WF3-1)]
        Random.seed!(4)
        wfchain_WF3 = Wright_Fisher_K_dim_exact_trajectory(rand(Dirichlet(K,0.3)), time_grid_WF3[1:(end-1)], α)
        wfobs_WF3 = [rand(Multinomial(Pop_size_WF3, wfchain_WF3[:,k])) for k in 1:size(wfchain_WF3,2)] |> l -> hcat(l...)
        time_grid_WF3 = [k*time_step_WF3 for k in 0:(Ntimes_WF3-1)]
        data_WF3 = Dict(zip(time_grid_WF3 , [wfobs_WF3[:,t] for t in 1:size(wfobs_WF3,2)]))
        return data_WF3, wfchain_WF3, α
    end
    data, wfchain_WF3, α = simulate_WF3_data()
    times = data |> keys |> collect |> sort

    current_logw = -0.5*ones(5,5,5)
    current_logw_prime = -0.5*ones(5,5,5)
    current_Λ_max = [2,1,3]
    current_Λ = DualOptimalFiltering.Λ_from_Λ_max(current_Λ_max)
    y = [2,3,1]

    smmax = values(data) |> sum |> sum
    log_ν_ar = Array{Float64}(undef, smmax, smmax)
    log_Cmmi_ar = Array{Float64}(undef, smmax, smmax)
    log_binomial_coeff_ar_offset = Array{Float64}(undef, smmax+1, smmax+1)

    @test_nowarn DualOptimalFiltering.precompute_next_terms_ar!(0, 12, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, 3.1, 0.2)

    @test_nowarn DualOptimalFiltering.WF_loglikelihood_from_adaptive_filtering(α, data, (x, y)-> (x, y); silence = false)

    @test_nowarn DualOptimalFiltering.log_likelihood_WF_fixed_number(α, data, 3; silence = false)

    @test_nowarn DualOptimalFiltering.log_likelihood_WF_fixed_fraction(α, data, 0.99; silence = false)

    @test_nowarn DualOptimalFiltering.log_likelihood_WF_keep_above_threshold(α, data, 0.001; silence = false)

end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Test
using ParseTorrent: pieces2hashes, parsetorrent

@testset "parsetorrent with a file" begin
    @test parsetorrent(read("torrents/bitlove-intro.torrent")) == Dict{Any,Any}(
        "infohash" => "4cb67059ed6bd08362da625b3ae77f6f4a075705",
        "announce" => "http://t.bitlove.org/announce",
        "info" => Dict{Any,Any}(
            "length" => 19211729,
            "name" => "bl001-introduction.webm",
            "pieces" => [
                "90a75dcd4e88d287c7ac5599c108f6036c13c4ce",
                "1ef5468bdff9a4466ad4e446477981cb67d07933",
                "1fa911a663451280953edb723e67611957dc0fe1",
                "2abad6066e29c723f01b0908ec30e0e737514a88",
                "55afda8e14a45e7f797eb47b82b2d0a3b2ca5f36",
                "7e1f49593515ca1b93ad01c3ee050e35f04f5c2e",
                "15b9abb123228002cca6a7d88fc9fc99d24583e1",
                "32704a020d2f121bfc612b7627cd92e2b39ad43c",
                "35bebb2888f7143c2966bb4d5f74e0b875825856",
                "6875f4bb1a9fa631ee35bcd7469b1e8ff37d65a2",
                "cbbeeeadc148ed681b699e88a940f796f51c0915",
                "c69121c81d85055678bf198bb29fc9e504ed8c7f",
                "7e3863c6e1c6a8c824569f1cc0950498dceb03c4",
                "ab4e77dade5f54246559c40915b700a4f734cee0",
                "92c47be2d397afbf06a9e9a573a63a3c683d2aa5",
                "01ad212a1495208b7ffbb173ce5782291695652b",
                "3f6233bf4ea3649c7799a1848f06cade97987525",
                "db37c799e45bd02fc25eacc12e18c6c11b4da3fb",
                "4c73df9307b3939fec3cd5f0df179c50a49c6ca3"
            ],
            "piece length" => 1048576),
            "url-list" => Any["http://spaceboyz.net/~astro/bitlove-show/bl001-introduction.webm"]
        )
end

# @testset "pieces2hashes" begin
#     @test pieces2hashes(Vector{UInt8}(1:100)) == []
# end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6a2c9516-dbda-4005-8d85-c19b5c0f54ee
begin
  using Pkg
  Pkg.activate(Base.Filesystem.homedir() * "/.config/julia/projects/oft")
  using Plots
  using LinearAlgebra
end

# ╔═╡ 4cdf2cb0-568c-47f6-bb13-a4a693e042f6
md"""
## Affine Transformations
An _affine transformation_ is a _linear transformation_ plus _a displacment_.
```math
\begin{array}
  &\mathbb{R}^2 &\to &\mathbb{R}^2& \\
  \begin{bmatrix}
    t_{1}' \\
    t_{2}'
  \end{bmatrix}
  &\mapsto
  &A
  \begin{bmatrix}
    t_{1} \\
    t_{2}
\end{bmatrix}
  + b&
  =
  \begin{bmatrix}
    a_{1,1} & a_{1,2} \\
    a_{2,1} & a_{2,2} \\
  \end{bmatrix}
  \begin{bmatrix}
    t_{1} \\
    t_{2}
  \end{bmatrix}
  +
  \begin{bmatrix}
    b_{1} \\
    b_{2}
  \end{bmatrix}
\end{array}
```
"""

# ╔═╡ f4611f9f-7f42-4164-8135-0fdcbf8f7ca0
begin
  chapeau_xy = [-1 -1 7 4 7 -1
                -5  5 5 3 1  1]
  function draw_flag(A, b, linecolor; new=false)
    if new
      plot()
    end
    xmax = 15
    ymax = xmax
    # chapeau_xy = [-1 -1 7 4 7 -1
    #               -5  5 5 3 1  1]
    transformed = A*chapeau_xy .+ b
    plot!(transformed[1, 1:end],
          transformed[2, 1:end],
          lw=3,
          framestyle=:zerolines,
          aspect_ratio=:equal,
          linecolor=linecolor,
          legend=false,
          xlims=[-xmax, xmax],
          ylims=[-ymax, ymax],
    )
  end
end

# ╔═╡ 5c95082e-75e5-466d-89e5-11b501add71d
md"""
### Translation
"""

# ╔═╡ 6aa15192-3ec7-417f-9847-97248d7b3ca5
begin
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(I(2), [4 3]', :red)
end

# ╔═╡ 434c9a7f-14f1-488e-9580-78b4f250e126
begin
  draw_flag(I(2), [1 2]', :grey; new=true)
  draw_flag(I(2), [-3 4]', :red)
end

# ╔═╡ b7e57e32-d5f1-40ba-b833-01e131ca0e64
md"""
### Scaling
"""

# ╔═╡ cfc871af-0ab3-4f7c-8f3d-a1421a32f3f5
begin
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(2*I(2), [0 0]', :red)
end

# ╔═╡ bd68e9bc-bff8-4820-9e8a-cbcfeee76371
begin
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag((1/3)*I(2), [0 0]', :red)
end

# ╔═╡ 8ca625d1-a595-4b00-8b0c-810c04d5927e
let  # We use let here because we shall repeatedly use the matrix A
  A = [2 0
       0 1/2]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ b34d9ee0-13ea-4b5b-846e-69785abadec2
md"""
### Rotation
"""

# ╔═╡ 6154092f-a161-45a4-8b86-5c3e436c1442
let
  θ = 2π / 3
  A = [cos(θ) -sin(θ)
       sin(θ)  cos(θ)]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ c48341a5-4434-4816-a026-638e59d1344c
let
  θ = -π / 6
  A = [cos(θ) -sin(θ)
       sin(θ)  cos(θ)]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ 2cc2fab1-ca14-4940-9050-e8f7d145a646
md"""
### Flips
"""

# ╔═╡ d9fe42fd-587f-4a1e-bd93-6553208109fe
let
  # horizontal flip
  A = [-1 0
        0 1]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ 99f2817e-f951-4123-9347-79c071329748
let
  # vertical flip
  A = [1  0
       0 -1]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ 516ebe04-e164-424b-832d-c9a91583788c
md"""
### Shearing
"""

# ╔═╡ ae6c8a6b-0ab5-4a0b-b058-2f5e3fecb0af
let
  s = 1.5
  A = [1 s
       0 1]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ fb1068ee-597a-44d5-be66-db42be7d7fb5
let
  s = 1.5
  A = [1 0
       s 1]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ a38ed456-87e3-4572-be89-e4f34d4196f9
let
  s = 1.5
  A = [1 -s
       0 1]
  draw_flag(I(2), [0 0]', :grey; new=true)
  draw_flag(A, [0 0]', :red)
end

# ╔═╡ Cell order:
# ╠═6a2c9516-dbda-4005-8d85-c19b5c0f54ee
# ╟─4cdf2cb0-568c-47f6-bb13-a4a693e042f6
# ╠═f4611f9f-7f42-4164-8135-0fdcbf8f7ca0
# ╟─5c95082e-75e5-466d-89e5-11b501add71d
# ╠═6aa15192-3ec7-417f-9847-97248d7b3ca5
# ╠═434c9a7f-14f1-488e-9580-78b4f250e126
# ╟─b7e57e32-d5f1-40ba-b833-01e131ca0e64
# ╠═cfc871af-0ab3-4f7c-8f3d-a1421a32f3f5
# ╠═bd68e9bc-bff8-4820-9e8a-cbcfeee76371
# ╠═8ca625d1-a595-4b00-8b0c-810c04d5927e
# ╟─b34d9ee0-13ea-4b5b-846e-69785abadec2
# ╠═6154092f-a161-45a4-8b86-5c3e436c1442
# ╠═c48341a5-4434-4816-a026-638e59d1344c
# ╟─2cc2fab1-ca14-4940-9050-e8f7d145a646
# ╠═d9fe42fd-587f-4a1e-bd93-6553208109fe
# ╠═99f2817e-f951-4123-9347-79c071329748
# ╟─516ebe04-e164-424b-832d-c9a91583788c
# ╠═ae6c8a6b-0ab5-4a0b-b058-2f5e3fecb0af
# ╠═fb1068ee-597a-44d5-be66-db42be7d7fb5
# ╠═a38ed456-87e3-4572-be89-e4f34d4196f9


text:
### DCOPF without renewable energy sources
function SCOPF(time,
			   nodes::Nodes,
			   lines::Lines,
			   powerplants::PowerPlants,
			   )

	T = time
	G = powerplants.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	lines_pmax = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		lines_pmax[j, k] = lines.pmax[id]
		lines_pmax[k, j] = lines.pmax[id]
	end

	# Create the B_prime matrix in one step
	B_prime = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		X = lines.reactance[id]				# reactance

		B_prime[j, k] = 1/X					# Fill YBUS with susceptance
		B_prime[k, j] = 1/X					# Symmetric matrix
	end

	for j in 1:size(B_prime)[1]
		B_prime[j, j] = -sum(B_prime[j, :])
	end

	B_prime_s = zeros(length(nodes.id)-1, length(nodes.id)-1, length(lines.id))
	for s in lines.id
		B_temp = B_prime[setdiff(1:end, s), setdiff(1:end, s)]
		B_prime_s[:,:,s] = B_temp
	end

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	K = nodes.id
	S = lines.id

	# Initialise JuMP model: DCOPF
	SCOPF_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(SCOPF_mod, P[T, G] >= 0)
	@variable(SCOPF_mod, Θ[S, T, J])
	@variable(SCOPF_mod, P_inj[T, J])
	@variable(SCOPF_mod, P_flow[S, T, K, J])

	# Fix the voltage angle of the slack bus
	for s in S
		for t in T
			fix(Θ[s, t, slack], 0)
		end
	end

	# Constraints
	# Market clearing/power balance
	@constraint(SCOPF_mod, MarketClearing[t = T],
	    sum(P[t, g] for g in G) == nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(SCOPF_mod, GenerationLimitUp[g = G, t = T],
	    P[t, g] <= powerplants.pmax[g]/Sbase);

	# Lower generation limit
	@constraint(SCOPF_mod, GenerationLimitDown[g = G, t = T],
	    P[t, g] >= powerplants.pmin[g]/Sbase);

	# Ramp-up limit
	@constraint(SCOPF_mod, RampUp[g = G, t = T; t > 1],
		(P[t, g]) - (P[t-1, g]) <= powerplants.rup[g])

	# Ramp-down limit
	@constraint(SCOPF_mod, RampDown[g = G, t = T; t > 1],
		-((P[t, g]) - (P[t-1, g])) <= powerplants.rdn[g])

	### DC power flow constraints
	# Power injection balance
	@constraint(SCOPF_mod, PowerInjectionBal[j = J, t = T],
		sum(P[t, g] for g in getKeyVector(powerplants.node, j)) -	# all generation units connected to a node
			nodes.load[j][t]/Sbase == P_inj[t, j]);

	# Power injection angle
	@constraint(SCOPF_mod, PowerInjectionAng[j = J, t = T, s = S],
		sum(B_prime[j, k] * (Θ[s, t, j] - Θ[s, t, k]) for k in K)
			== P_inj[t, j]);

	@constraint(SCOPF_mod, LinePowerFlow[k = K, j = J, t = T, s = S; j != k &&
		k != lines.from[s] && j != lines.to[s]],
		B_prime[j, k] * (Θ[s, t, j] - Θ[s, t, k]) == P_flow[s, t, j, k]);

	@constraint(SCOPF_mod, LinePowerFlowMax[k = K, j = J, t = T, s = S;
		k != lines.from[s] && j != lines.to[s]],
		P_flow[s, t, j, k] <= lines_pmax[j, k]/Sbase);

	@constraint(SCOPF_mod, LinePowerFlowMin[k = K, j = J, t = T, s = S;
		k != lines.from[s] && j != lines.to[s]],
		P_flow[s, t, j, k] >= -lines_pmax[j, k]/Sbase);

	# Contingency: Line outage s = S
	@constraint(SCOPF_mod, Contingency[k = K, j = J, t = T, s = S;
		k == lines.from[s] && j == lines.to[s]],
		P_flow[s, t, j, k] == 0);

	# Objective function
	@objective(SCOPF_mod, Min,
	    sum(powerplants.mc[g] * P[t, g] * Sbase for g in G, t in T));

	# Initiate optimisation process
	JuMP.optimize!(SCOPF_mod)

	# Export results
	P_opt = JuMP.value.(P)*Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	P_inj_opt = JuMP.value.(P_inj)*Sbase

	return(SCOPF_mod,
	       P_opt,
	       Θ_opt,
	       P_inj_opt,
	       P_flow_opt
	       )
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
const DB_NETCDF = 0
const DB_PDB = 2 # PDB Lite
const DB_TAURUS = 3
const DB_UNKNOWN = 5
const DB_DEBUG = 6
const DB_HDF5X = 7
const DB_PDBP = 1 # PDB Proper
const DB_HDF5_SEC2_OBSOLETE = 0x100
const DB_HDF5_STDIO_OBSOLETE = 0x200
const DB_HDF5_CORE_OBSOLETE = 0x300
const DB_HDF5_MPIO_OBSOLETE = 0x400
const DB_HDF5_MPIOP_OBSOLETE = 0x500
const DB_H5VFD_DEFAULT = 0
const DB_H5VFD_SEC2 = 1
const DB_H5VFD_STDIO = 2
const DB_H5VFD_CORE = 3
const DB_H5VFD_LOG = 4
const DB_H5VFD_SPLIT = 5
const DB_H5VFD_DIRECT = 6
const DB_H5VFD_FAMILY = 7
const DB_H5VFD_MPIO = 8
const DB_H5VFD_MPIP = 9
const DB_H5VFD_SILO = 10
const DB_H5VFD_FIC = 11 # File Image in Core
const DB_NFILES = 256 #Max simultaneously open files
const DB_NFILTERS = 32 #Number of filters defined
const DBAll = 0xffffffffffffffff
const DBNone = 0x0000000000000000
const DBCalc = 0x0000000000000001
const DBMatMatnos = 0x0000000000000002
const DBMatMatlist = 0x0000000000000004
const DBMatMixList = 0x0000000000000008
const DBCurveArrays = 0x0000000000000010
const DBPMCoords = 0x0000000000000020
const DBPVData = 0x0000000000000040
const DBQMCoords = 0x0000000000000080
const DBQVData = 0x0000000000000100
const DBUMCoords = 0x0000000000000200
const DBUMFacelist = 0x0000000000000400
const DBUMZonelist = 0x0000000000000800
const DBUVData = 0x0000000000001000
const DBFacelistInfo = 0x0000000000002000
const DBZonelistInfo = 0x0000000000004000
const DBMatMatnames = 0x0000000000008000
const DBUMGlobNodeNo = 0x0000000000010000
const DBZonelistGlobZoneNo = 0x0000000000020000
const DBMatMatcolors = 0x0000000000040000
const DBCSGMBoundaryInfo = 0x0000000000080000
const DBCSGMZonelist = 0x0000000000100000
const DBCSGMBoundaryNames = 0x0000000000200000
const DBCSGVData = 0x0000000000400000
const DBCSGZonelistZoneNames = 0x0000000000800000
const DBCSGZonelistRegNames = 0x0000000001000000
const DBMMADJNodelists = 0x0000000002000000
const DBMMADJZonelists = 0x0000000004000000
const DBPMGlobNodeNo = 0x0000000008000000
const DBPMGhostNodeLabels = 0x0000000010000000
const DBQMGhostNodeLabels = 0x0000000020000000
const DBQMGhostZoneLabels = 0x0000000040000000
const DBUMGhostNodeLabels = 0x0000000080000000
const DBZonelistGhostZoneLabels = 0x0000000100000000
const DB_COLLINEAR = 130
const DB_NONCOLLINEAR = 131
const DB_CLOBBER = 0
const DB_NOCLOBBER = 1
const DB_READ = 1
const DB_APPEND = 2
const DB_LOCAL = 0
const DB_SUN3 = 10
const DB_SUN4 = 11
const DB_SGI = 12
const DB_RS6000 = 13
const DB_CRAY = 14
const DB_INTEL = 15
const DBOPT_FIRST = 260
const DBOPT_ALIGN = 260
const DBOPT_COORDSYS = 262
const DBOPT_CYCLE = 263
const DBOPT_FACETYPE = 264
const DBOPT_HI_OFFSET = 265
const DBOPT_LO_OFFSET = 266
const DBOPT_LABEL = 267
const DBOPT_XLABEL = 268
const DBOPT_YLABEL = 269
const DBOPT_ZLABEL = 270
const DBOPT_MAJORORDER = 271
const DBOPT_NSPACE = 272
const DBOPT_ORIGIN = 273
const DBOPT_PLANAR = 274
const DBOPT_TIME = 275
const DBOPT_UNITS = 276
const DBOPT_XUNITS = 277
const DBOPT_YUNITS = 278
const DBOPT_ZUNITS = 279
const DBOPT_DTIME = 280
const DBOPT_USESPECMF = 281
const DBOPT_XVARNAME = 282
const DBOPT_YVARNAME = 283
const DBOPT_ZVARNAME = 284
const DBOPT_ASCII_LABEL = 285
const DBOPT_MATNOS = 286
const DBOPT_NMATNOS = 287
const DBOPT_MATNAME = 288
const DBOPT_NMAT = 289
const DBOPT_NMATSPEC = 290
const DBOPT_BASEINDEX = 291 # quad meshes for node and zone
const DBOPT_ZONENUM = 292 # ucd meshes for zone
const DBOPT_NODENUM = 293 # ucd/point meshes for node
const DBOPT_BLOCKORIGIN = 294
const DBOPT_GROUPNUM = 295
const DBOPT_GROUPORIGIN = 296
const DBOPT_NGROUPS = 297
const DBOPT_MATNAMES = 298
const DBOPT_EXTENTS_SIZE = 299
const DBOPT_EXTENTS = 300
const DBOPT_MATCOUNTS = 301
const DBOPT_MATLISTS = 302
const DBOPT_MIXLENS = 303
const DBOPT_ZONECOUNTS = 304
const DBOPT_HAS_EXTERNAL_ZONES = 305
const DBOPT_PHZONELIST = 306
const DBOPT_MATCOLORS = 307
const DBOPT_BNDNAMES = 308
const DBOPT_REGNAMES = 309
const DBOPT_ZONENAMES = 310
const DBOPT_HIDE_FROM_GUI = 311
const DBOPT_TOPO_DIM = 312 # TOPOlogical DIMension
const DBOPT_REFERENCE = 313 # reference object
const DBOPT_GROUPINGS_SIZE = 314 # size of grouping array
const DBOPT_GROUPINGS = 315 # groupings array
const DBOPT_GROUPINGNAMES = 316 # array of size determined by
const DBOPT_ALLOWMAT0 = 317 # Turn off material numer "0" warnings
const DBOPT_MRGTREE_NAME = 318
const DBOPT_REGION_PNAMES = 319
const DBOPT_TENSOR_RANK = 320
const DBOPT_MMESH_NAME = 321
const DBOPT_TV_CONNECTIVITY = 322
const DBOPT_DISJOINT_MODE = 323
const DBOPT_MRGV_ONAMES = 324
const DBOPT_MRGV_RNAMES = 325
const DBOPT_SPECNAMES = 326
const DBOPT_SPECCOLORS = 327
const DBOPT_LLONGNZNUM = 328
const DBOPT_CONSERVED = 329
const DBOPT_EXTENSIVE = 330
const DBOPT_MB_FILE_NS = 331
const DBOPT_MB_BLOCK_NS = 332
const DBOPT_MB_BLOCK_TYPE = 333
const DBOPT_MB_EMPTY_LIST = 334
const DBOPT_MB_EMPTY_COUNT = 335
const DBOPT_MB_REPR_BLOCK_IDX = 336
const DBOPT_MISSING_VALUE = 337
const DBOPT_ALT_ZONENUM_VARS = 338
const DBOPT_ALT_NODENUM_VARS = 339
const DBOPT_GHOST_NODE_LABELS = 340
const DBOPT_GHOST_ZONE_LABELS = 341
const DBOPT_LAST = 499
const DBOPT_H5_FIRST = 500
const DBOPT_H5_VFD = 500
const DBOPT_H5_RAW_FILE_OPTS = 501
const DBOPT_H5_RAW_EXTENSION = 502
const DBOPT_H5_META_FILE_OPTS = 503
const DBOPT_H5_META_EXTENSION = 504
const DBOPT_H5_CORE_ALLOC_INC = 505
const DBOPT_H5_CORE_NO_BACK_STORE = 506
const DBOPT_H5_META_BLOCK_SIZE = 507
const DBOPT_H5_SMALL_RAW_SIZE = 508
const DBOPT_H5_ALIGN_MIN = 509
const DBOPT_H5_ALIGN_VAL = 510
const DBOPT_H5_DIRECT_MEM_ALIGN = 511
const DBOPT_H5_DIRECT_BLOCK_SIZE = 512
const DBOPT_H5_DIRECT_BUF_SIZE = 513
const DBOPT_H5_LOG_NAME = 514
const DBOPT_H5_LOG_BUF_SIZE = 515
const DBOPT_H5_MPIO_COMM = 516
const DBOPT_H5_MPIO_INFO = 517
const DBOPT_H5_MPIP_NO_GPFS_HINTS = 518
const DBOPT_H5_SIEVE_BUF_SIZE = 519
const DBOPT_H5_CACHE_NELMTS = 520
const DBOPT_H5_CACHE_NBYTES = 521
const DBOPT_H5_CACHE_POLICY = 522
const DBOPT_H5_FAM_SIZE = 523
const DBOPT_H5_FAM_FILE_OPTS = 524
const DBOPT_H5_USER_DRIVER_ID = 525
const DBOPT_H5_USER_DRIVER_INFO = 526
const DBOPT_H5_SILO_BLOCK_SIZE = 527
const DBOPT_H5_SILO_BLOCK_COUNT = 528
const DBOPT_H5_SILO_LOG_STATS = 529
const DBOPT_H5_SILO_USE_DIRECT = 530
const DBOPT_H5_FIC_SIZE = 531
const DBOPT_H5_FIC_BUF = 532
const DBOPT_H5_LAST = 599
const DB_TOP = 0 #default--API traps
const DB_NONE = 1 #no errors trapped
const DB_ALL = 2 #all levels trap (traceback)
const DB_ABORT = 3 #abort() is called
const DB_SUSPEND = 4 #suspend error reporting temporarily
const DB_RESUME = 5 #resume normal error reporting
const DB_ALL_AND_DRVR = 6 #DB_ALL + driver error reporting
const E_NOERROR = 0 #No error
const E_BADFTYPE = 1 #Bad file type
const E_NOTIMP = 2 #Callback not implemented
const E_NOFILE = 3 #No data file specified
const E_INTERNAL = 5 #Internal error
const E_NOMEM = 6 #Not enough memory
const E_BADARGS = 7 #Bad argument to function
const E_CALLFAIL = 8 #Low-level function failure
const E_NOTFOUND = 9 #Object not found
const E_TAURSTATE = 10 #Taurus: database state error
const E_MSERVER = 11 #SDX: too many connections
const E_PROTO = 12 #SDX: protocol error
const E_NOTDIR = 13 #Not a directory
const E_MAXOPEN = 14 #Too many open files
const E_NOTFILTER = 15 #Filter(s) not found
const E_MAXFILTERS = 16 #Too many filters
const E_FEXIST = 17 #File already exists
const E_FILEISDIR = 18 #File is actually a directory
const E_FILENOREAD = 19 #File lacks read permission.
const E_SYSTEMERR = 20 #System level error occured.
const E_FILENOWRITE = 21 #File lacks write permission.
const E_INVALIDNAME = 22 # Variable name is invalid
const E_NOOVERWRITE = 23 #Overwrite not permitted
const E_CHECKSUM = 24 #Checksum failed
const E_COMPRESSION = 25 #Compression failed
const E_GRABBED = 26 #Low level driver enabled
const E_NOTREG = 27 #The dbfile pointer is not resitered.
const E_CONCURRENT = 28 #File is opened multiply and not all read-only.
const E_DRVRCANTOPEN = 29 #Driver cannot open the file.
const E_BADOPTCLASS = 30 #Optlist contains options for wrong class
const E_NOTENABLEDINBUILD = 31 #feature not enabled in this build
const E_MAXFILEOPTSETS = 32 #Too many file options sets
const E_NOHDF5 = 33 #HDF5 driver not available
const E_EMPTYOBJECT = 34 #Empty object not currently permitted
const E_OBJBUFF = 35 #No more temp. buffer space for object
const E_NERRORS = 50
const DB_ROWMAJOR = 0
const DB_COLMAJOR = 1
const DB_NOTCENT = 0
const DB_NODECENT = 110
const DB_ZONECENT = 111
const DB_FACECENT = 112
const DB_BNDCENT = 113 # for CSG meshes only
const DB_EDGECENT = 114
const DB_BLOCKCENT = 115 # for 'block-centered' data on multimeshs
const DB_CARTESIAN = 120
const DB_CYLINDRICAL = 121 # x,r; y,theta; z,height; 2D variant is eqiv. to poloar
const DB_SPHERICAL = 122 # x,r; y,theta; z,phi; 2D variant is equiv. to polar
const DB_NUMERICAL = 123
const DB_OTHER = 124
const DB_RECTILINEAR = 100
const DB_CURVILINEAR = 101
const DB_AREA = 140
const DB_VOLUME = 141
const DB_ON = 1000
const DB_ABUTTING = 142
const DB_FLOATING = 143
const DB_VARTYPE_SCALAR = 200
const DB_VARTYPE_VECTOR = 201
const DB_VARTYPE_TENSOR = 202
const DB_VARTYPE_SYMTENSOR = 203
const DB_VARTYPE_ARRAY = 204
const DB_VARTYPE_MATERIAL = 205
const DB_VARTYPE_SPECIES = 206
const DB_VARTYPE_LABEL = 207
const DBCSG_QUADRIC_G = 0x01000000
const DBCSG_SPHERE_PR = 0x02010000
const DBCSG_ELLIPSOID_PRRR = 0x02020000
const DBCSG_PLANE_G = 0x03000000
const DBCSG_PLANE_X = 0x03010000
const DBCSG_PLANE_Y = 0x03020000
const DBCSG_PLANE_Z = 0x03030000
const DBCSG_PLANE_PN = 0x03040000
const DBCSG_PLANE_PPP = 0x03050000
const DBCSG_CYLINDER_PNLR = 0x04000000
const DBCSG_CYLINDER_PPR = 0x04010000
const DBCSG_BOX_XYZXYZ = 0x05000000
const DBCSG_CONE_PNLA = 0x06000000
const DBCSG_CONE_PPA = 0x06010000
const DBCSG_POLYHEDRON_KF = 0x07000000
const DBCSG_HEX_6F = 0x07010000
const DBCSG_TET_4F = 0x07020000
const DBCSG_PYRAMID_5F = 0x07030000
const DBCSG_PRISM_5F = 0x07040000
const DBCSG_QUADRATIC_G = 0x08000000
const DBCSG_CIRCLE_PR = 0x09000000
const DBCSG_ELLIPSE_PRR = 0x09010000
const DBCSG_LINE_G = 0x0A000000
const DBCSG_LINE_X = 0x0A010000
const DBCSG_LINE_Y = 0x0A020000
const DBCSG_LINE_PN = 0x0A030000
const DBCSG_LINE_PP = 0x0A040000
const DBCSG_BOX_XYXY = 0x0B000000
const DBCSG_ANGLE_PNLA = 0x0C000000
const DBCSG_ANGLE_PPA = 0x0C010000
const DBCSG_POLYGON_KP = 0x0D000000
const DBCSG_TRI_3P = 0x0D010000
const DBCSG_QUAD_4P = 0x0D020000
const DBCSG_INNER = 0x7F000000
const DBCSG_OUTER = 0x7F010000
const DBCSG_ON = 0x7F020000
const DBCSG_UNION = 0x7F030000
const DBCSG_INTERSECT = 0x7F040000
const DBCSG_DIFF = 0x7F050000
const DBCSG_COMPLIMENT = 0x7F060000
const DBCSG_XFORM = 0x7F070000
const DBCSG_SWEEP = 0x7F080000
const DB_PREORDER = 0x00000001
const DB_POSTORDER = 0x00000002
const DB_FROMCWR = 0x00000004
const DB_ZONETYPE_BEAM = 10
const DB_ZONETYPE_POLYGON = 20
const DB_ZONETYPE_TRIANGLE = 23
const DB_ZONETYPE_QUAD = 24
const DB_ZONETYPE_POLYHEDRON = 30
const DB_ZONETYPE_TET = 34
const DB_ZONETYPE_PYRAMID = 35
const DB_ZONETYPE_PRISM = 36
const DB_ZONETYPE_HEX = 38
const DB_MAX_H5_OBJ_VALS = 64
const DB_MAX_EXPSTRS = 8 # NO_FORTRAN_DEFINE


text:
module NDS

import ..GreyNumbers: GreyNumber

function dominates(p1::Array{GreyNumber{T}, 1}, p2::Array{GreyNumber{T}, 1})::Bool where {T <: Real}
    n = length(p1)
    notworse = count(i -> p1[i] < p2[i], 1:n)
    better   = count(i -> p1[i] > p2[i], 1:n)
    return (notworse == 0) && (better > 0)
end

function ndsranks(data::Array{GreyNumber{T}, 2})::Array{Int64} where {T <: Real}

    n, p = size(data)

    ranks = zeros(Int64, n)

    mat = convert(Matrix, data)

     for i in 1:n
        for j in 1:n
            if i != j
                if dominates(mat[i,:], mat[j,:])
                    ranks[i] += 1
                end
            end
        end
    end

    return ranks
end

function nds(data::Array{GreyNumber{T}, 2})::Dict where {T <: Real}

    ranks = ndsranks(data)

    bestIndex = sortperm(ranks) |> last

    result = Dict(
        :ranks => ranks,
        :bestindex => bestIndex
    )

    return result
end

end # end of module NDS


text:
module BP

using Polynomials
using Base.Test

type WaveFunction
    N::Int
    int::Float64
    ψ::Vector{Complex{Float64}}
end
WaveFunction(S::SparseMatrixCSC{Complex{Float64},Int64},
             ω::Float64,P::Vector{Complex{Float64}}) =
                WaveFunction(S,ω,P,:landau)
WaveFunction(S::SparseMatrixCSC{Complex{Float64},Int64},
             ω::Float64,P::Vector{Complex{Float64}},gauge::Symbol) =
    WaveFunction(S,ω,P,gauge,1/11,0.001,0.02,0.,0.)
WaveFunction(S::SparseMatrixCSC{Complex{Float64},Int64},
             ω::Float64,P::Vector{Complex{Float64}},gauge::Symbol,
             α::Float64,γ::Float64,κ::Float64) =
                 WaveFunction(S,ω,P,gauge,α,γ,κ, 0., 0.)

function WaveFunction(S::SparseMatrixCSC{Complex{Float64},Int64},
                      ω::Float64,P::Vector{Complex{Float64}},gauge::Symbol,
                      α::Float64,γ::Float64,κ::Float64, m₀::Float64,
                      n₀::Float64)
    N::Int = sqrt(length(P))
    eval(:($(symbol(string("buildham_", gauge, "!")))))(S, N,α,κ,γ,ω,m₀,n₀)
    X = S\P
    return WaveFunction(N,sum(abs2(X)),X)
end

type ExactStates
    N::Int
    gauge::Symbol
    νs::Vector{Float64}
    states::Matrix{Complex{Float64}}
end
ExactStates(nev::Int, gauge::Symbol) = ExactStates(nev, gauge, 45)
ExactStates(nev::Int, gauge::Symbol, N::Int) = ExactStates(nev, gauge,
   N, 1/11, 0.02, 0., 0.)
ExactStates(nev::Int, gauge::Symbol, N::Int, α::Float64, κ::Float64) =
   ExactStates(nev, gauge, N, α, κ, 0., 0.)

function ExactStates(nev::Int, gauge::Symbol, N::Int, α::Float64,
   κ::Float64, m₀::Float64, n₀::Float64)
    M = spzeros(Complex{Float64}, N^2,N^2)
    eval(:($(symbol(string("buildhamexact", gauge, "!")))))(M, N,α,κ,m₀,n₀)
    (d, v, nconv, niter, nmult, resid) = eigs(M, nev=nev, which=:SR,
       ritzvec=true)
    return ExactStates(N, gauge, real(d), v)
end

function getstate(s::ExactStates, ω::Float64)
    i::Int = indmin(abs(s.νs .- ω))
    return reshape(s.states[:,i], (s.N, s.N))
end
function getstate(s::ExactStates, η::Int)
    return reshape(s.states[:,η], (s.N, s.N))
end

type Spectrum
    N::Int
    gauge::Symbol
    pump::Vector{Complex{Float64}}
    νs::Vector{Float64}
    intensity::Vector{Float64}
    states::Vector{WaveFunction}
end
Spectrum(ν::Vector{Float64},P::Vector{Complex{Float64}}) = Spectrum(ν,P,:landau)

function Spectrum(ν::Vector{Float64},P::Vector{Complex{Float64}},gauge::Symbol)
    statevec = Array(WaveFunction, length(ν))
    intvec = Array(Float64, length(ν))
    N::Int = sqrt(length(P))
    A = spzeros(Complex{Float64}, N^2,N^2)
    for (i,ω) in enumerate(ν)
        statevec[i] = WaveFunction(A, ω, P, gauge)
        intvec[i] = statevec[i].int
    end
    return Spectrum(N, gauge, P, ν, intvec, statevec)
end
function Spectrum(ν::Vector{Float64},P::Vector{Complex{Float64}},gauge::Symbol,
                  α::Float64,γ::Float64,κ::Float64)
    statevec = Array(WaveFunction, length(ν))
    intvec = Array(Float64, length(ν))
    N::Int = sqrt(length(P))
    A = spzeros(Complex{Float64}, N^2,N^2)
    for (i,ω) in enumerate(ν)
        statevec[i] = WaveFunction(A, ω, P, gauge, α,γ,κ)
        intvec[i] = statevec[i].int
    end
    return Spectrum(N, gauge, P, ν, intvec, statevec)
end
function Spectrum(ν::Vector{Float64},P::Vector{Complex{Float64}},gauge::Symbol,
                  α::Float64,γ::Float64,κ::Float64, m₀::Float64, n₀::Float64)
    statevec = Array(WaveFunction, length(ν))
    intvec = Array(Float64, length(ν))
    N::Int = sqrt(length(P))
    A = spzeros(Complex{Float64}, N^2,N^2)
    for (i,ω) in enumerate(ν)
        statevec[i] = WaveFunction(A, ω, P, gauge, α,γ,κ, m₀,n₀)
        intvec[i] = statevec[i].int
    end
    return Spectrum(N, gauge, P, ν, intvec, statevec)
end

function getstate(s::Spectrum, ω::Float64)
    i::Int = indmin(abs(s.νs .- ω))
    return reshape(s.states[i].ψ, (s.N, s.N))
end

#Check that matrix is square
function chksquare(A::AbstractMatrix)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square"))
    m
end

#Find radius of ring
function radius(M::Matrix{Float64}, axis::Vector)
    N = chksquare(M)
    isodd(N) || throw(DimensionMismatch("invalid matrix size N=$N. N
       must be odd"))
    length(axis) == N || throw(DimensionMismatch("axis size must match
       matrix dimensions"))
    @test_approx_eq_eps(M, transpose(M), 1e-5)
    half = div(N-1, 2) + 1
    v1 = axis[half:end]
    v2 = M[half:end, half]
    i = indmax(v2)
    r = v1[i]
    return r
end

getm(i::Int64,N::Int64) = div(i-1,N)-div(N-1,2)
getn(i::Int64,N::Int64) = div(N-1,2)-rem(i-1,N)
geti(m::Int,n::Int,N::Int)=(m+div(N-1,2))*N+(div(N-1,2)-n)+1

countentries(N::Int) = N^2 + 8 + 4*(N-2)*3 + (N-2)^2*4

macro hambody(fself, fleft, fright, fup, fdown)
    return quote
        border::Int = div(N-1,2)
        for m in -border:border, n in -border:border
            i  = geti(m,n,N)
            S[i,i] = $fself
        end
        for m in -border+1:border, n in -border:border
            i  = geti(m,n,N)
            S[i,i-N] = $fleft
        end
        for m in -border:border-1, n in -border:border
            i  = geti(m,n,N)
            S[i,i+N] = $fright
        end
        for m in -border:border, n in -border:border-1
            i  = geti(m,n,N)
            S[i,i-1] = $fup
        end
        for m in -border:border, n in -border+1:border
            i  = geti(m,n,N)
            S[i,i+1] = $fdown
        end
    end
end

buildham_landau!(S::SparseMatrixCSC{Complex{Float64},Int},
                 N::Int,α::Float64,κ::Float64,γ::Float64,ω::Float64) =
                     buildham_landau!(S,N,α,κ,γ,ω, 0.,0.)
function buildham_landau!(S::SparseMatrixCSC{Complex{Float64},Int},
   N::Int,α::Float64,κ::Float64,γ::Float64,ω::Float64, m₀::Float64,
   n₀::Float64)
    @hambody(ω + im*γ - 1/2*κ*((n-n₀)^2+(m-m₀)^2), 1, 1,
       exp(-im*2π*α*m), exp(im*2π*α*m))
end
function buildham_symmetric!(S::SparseMatrixCSC{Complex{Float64},Int},
   N::Int,α::Float64,κ::Float64,γ::Float64,ω::Float64, m₀::Float64,
   n₀::Float64)
    @hambody(ω + im*γ - 1/2*κ*((n-n₀)^2+(m-m₀)^2), exp(-im*π*α*n),
       exp(im*π*α*n), exp(-im*π*α*m), exp(im*π*α*m))
end
buildham_exact!(S::SparseMatrixCSC{Complex{Float64},Int},
   N::Int,α::Float64,κ::Float64) = buildham_exact!(S, N,α,κ, 0., 0.)
function buildham_exact!(S::SparseMatrixCSC{Complex{Float64},Int},
   N::Int,α::Float64,κ::Float64, m₀::Float64, n₀::Float64)
    @hambody(1/2*κ*((n-n₀)^2+(m-m₀)^2), -1, -1, -exp(-im*2π*α*m),
       -exp(im*2π*α*m))
end
buildhamexactlandau! = buildham_exact!
function
   buildhamexactsymmetric!(S::SparseMatrixCSC{Complex{Float64},Int},
   N::Int,α::Float64,κ::Float64, m₀::Float64, n₀::Float64)
    @hambody(1/2*κ*((n-n₀)^2+(m-m₀)^2), -exp(-im*π*α*n),
       -exp(im*π*α*n), -exp(-im*π*α*m), -exp(im*π*α*m))
end

function
   genspmat(l::Function,r::Function,u::Function,d::Function,s::Function,
   N::Int,nz::Int,α::Float64)
    iseven(N) && throw(ArgumentError("invalid system size N=$N. N must
       be odd"))
    # Preallocate
    I = Array(Int64,nz)
    J = Array(Int64,nz)
    V = Array(Complex{Float64},nz)
    function setnzelem(i::Int,n::Int,m::Int; pos::ASCIIString = "self")
        if pos=="left"
            k += 1
            J[k] = i-N; I[k] = i; V[k] = l(n,m,α)
        elseif pos=="right"
            k += 1
            J[k] = i+N; I[k] = i; V[k] = r(n,m,α)
        elseif pos=="up"
            k += 1
            J[k] = i-1; I[k] = i; V[k] = u(n,m,α)
        elseif pos=="down"
            k += 1
            J[k] = i+1; I[k] = i; V[k] = d(n,m,α)
        elseif pos=="self"
            k += 1
            J[k] = i; I[k] = i; V[k] = s(n,m,α)
        end
    end
    # maximum value of m or n indices
    maxm = div(N-1,2)
    k = 0
    for i in 1:N^2
        m = getm(i,N)
        n = getn(i,N)
        setnzelem(i,n,m; pos="self")
        #corners
        #top left
        if n==maxm && m==-maxm
            setnzelem(i,n,m; pos="right")
            setnzelem(i,n,m; pos="down")
        #top right
        elseif n==maxm && m==maxm
            setnzelem(i,n,m; pos="left")
            setnzelem(i,n,m; pos="down")
        #bottom right
        elseif n==-maxm && m==maxm
            setnzelem(i,n,m; pos="left")
            setnzelem(i,n,m; pos="up")
        #bottom left
        elseif n==-maxm && m==-maxm
            setnzelem(i,n,m; pos="right")
            setnzelem(i,n,m; pos="up")
        #edges
        #top
        elseif n == maxm
            setnzelem(i,n,m; pos="right")
            setnzelem(i,n,m; pos="left")
            setnzelem(i,n,m; pos="down")
        #right
        elseif m == maxm
            setnzelem(i,n,m; pos="left")
            setnzelem(i,n,m; pos="up")
            setnzelem(i,n,m; pos="down")
        #bottom
        elseif n == -maxm
            setnzelem(i,n,m; pos="left")
            setnzelem(i,n,m; pos="up")
            setnzelem(i,n,m; pos="right")
        #left
        elseif m == -maxm
            setnzelem(i,n,m; pos="down")
            setnzelem(i,n,m; pos="up")
            setnzelem(i,n,m; pos="right")
        else #bulk
            setnzelem(i,n,m; pos="down")
            setnzelem(i,n,m; pos="up")
            setnzelem(i,n,m; pos="right")
            setnzelem(i,n,m; pos="left")
        end
    end
    return sparse(I,J,V)
end

#various pumping schemes
function δpmp(N::Int; A=1., seed=0, σ=0., n0=0, m0=0)
    i = (m0+div(N-1,2)) * N + (div(N-1,2)-n0) + 1
    f = zeros(Complex{Float64}, N^2)
    f[i] = A * one(Complex{Float64})
    f
end
function gausspmp(N::Int; A=1., seed=0, σ=1., n0=0, m0=0)
    x0=m0
    y0=n0
    f = zeros(Complex{Float64}, N,N)
    m = collect(-div(N-1,2):div(N-1,2))
    n = collect(div(N-1,2):-1:-div(N-1,2))
    for c in 1:N, l in 1:N
        x = m[c]
        y = n[l]
        f[l,c] = A*exp(-1/(2σ^2)*((x-x0)^2 + (y-y0)^2))
    end
    reshape(f, N^2)
end
function randpmp(N::Int; A=1., seed=123, σ=0., n0=0, m0=0)
    # seed the RNG #
    srand(seed)
    # generate matrix of random phases in interval [0,2π)
    ϕ = 2π .* rand(N^2)
    A .* exp(im .* ϕ)
end
function homopmp(N::Int; A=1., seed=0, σ=0., n0=0, m0=0)
    A .* ones(Complex{Float64}, N^2)
end

#arbitrary resolution fft
function myfft2(ψr::Matrix{Complex{Float64}}, k1::Float64,
   k2::Float64, xs1::Float64, xs2::Float64, Δx1::Float64, Δx2::Float64)
    (N1,N2) = size(ψr)
    s = zero(Complex{Float64})
    for n2 in 1:N2, n1 in 1:N1
        xn1 = xs1 + (n2-1)*Δx1 #x
        xn2 = xs2 + (n1-1)*Δx2 #y
        cexp = exp(-im*(k1*xn1 + k2*xn2))
        s += ψr[n1,n2]*cexp
    end
    s
end
function myfft2(ψr::Matrix{Complex{Float64}}, k1, k2)
    N1 = length(k2); N2 = length(k1)
    out = Array(Complex{Float64}, N1,N2)
    for j in 1:N2, i in 1:N1
        out[i,j] = myfft2(ψr, k1[j], k2[i], 0., 0., 1., 1.)
    end
    out
end

# maps any momentum to the first Brillouin Zone
function fbz(mom::Float64)
    m = mod(mom, 2π)
    m <= π ? m : m - 2π
end
# Magnetic Brillouin Zone
function mbz(data::Array{Float64,2}, r::Int, q::Int,
   kxmbz::FloatRange{Float64}, k::FloatRange{Float64})
    # data is |ψ(FBZ)|², input
    V = zeros(Float64, (length(k), r)) # |ψ(MBZ)|², output
    for (i,px) in enumerate(kxmbz) # loop over momenta in MBZ
        for j = 0:q-1 # sum over various components
            ptld = px - j * 2π/q
            idx = indmin(abs(k .- fbz(ptld))) # index of \tilde{p_x}
            V[:,i] += data[:, idx]
        end
    end
    V/(4π^2)
end

#comparison to analytics
function compute_hermite_polynomial(n)
    P = Poly([1])
    const x = Poly([0; 1])
    for i = 1:n
        P = 2x*P - polyder(P)
    end
    P
end
function χ(kx0, ky, α, β)
    l = sqrt(2π*α)
    sum = zero(Complex{Float64})
    for j in -20:20 #truncate the sum
        H = polyval(compute_hermite_polynomial(β), kx0/l + j*l)
        sum += exp(-im*ky*j) * exp(-(kx0 + j*l^2)^2/(2l^2)) * H
    end
    sum
end
end #module

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon


text:
#

using Test
using RoME

##

@testset "test Boxes2D generator" begin

##

fg = generateCanonicalFG_Boxes2D!(8, postpose_cb=(g,l)->@show l)

##

@test isapprox( [0;0], getPPE(fg, :x0, :simulated).suggested, atol=1e-8)
@test isapprox( [15;0], getPPE(fg, :x1, :simulated).suggested, atol=1e-8)
@test isapprox( [15;15], getPPE(fg, :x2, :simulated).suggested, atol=1e-8)
@test isapprox( [5;15], getPPE(fg, :x3, :simulated).suggested, atol=1e-8)
@test isapprox( [5;0], getPPE(fg, :x4, :simulated).suggested, atol=1e-8)
@test isapprox( [20;0], getPPE(fg, :x5, :simulated).suggested, atol=1e-8)
@test isapprox( [20;15], getPPE(fg, :x6, :simulated).suggested, atol=1e-8)
@test isapprox( [10;15], getPPE(fg, :x7, :simulated).suggested, atol=1e-8)
@test isapprox( [10;0], getPPE(fg, :x8, :simulated).suggested, atol=1e-8)

##

solveTree!(fg)

##

@test isapprox( [0;0], getPPE(fg, :x0, :simulated).suggested, atol=1e-8)
@test isapprox( [15;0], getPPE(fg, :x1, :simulated).suggested, atol=1e-8)
@test isapprox( [15;15], getPPE(fg, :x2, :simulated).suggested, atol=1e-8)
@test isapprox( [5;15], getPPE(fg, :x3, :simulated).suggested, atol=1e-8)
@test isapprox( [5;0], getPPE(fg, :x4, :simulated).suggested, atol=1e-8)
@test isapprox( [20;0], getPPE(fg, :x5, :simulated).suggested, atol=1e-8)
@test isapprox( [20;15], getPPE(fg, :x6, :simulated).suggested, atol=1e-8)
@test isapprox( [10;15], getPPE(fg, :x7, :simulated).suggested, atol=1e-8)
@test isapprox( [10;0], getPPE(fg, :x8, :simulated).suggested, atol=1e-8)

##

end

@testset "test canonical helix generators" begin

##

# moved to IIF
# tmp = calcHelix_T(0, 3, 25, radius=5, xr_t=t->(1/3)*t)

##

cb(fg_, lp) = @show lp, length(ls(fg_))

fg = generateCanonicalFG_Helix2DSlew!(46, slew_x=2/3, posesperturn=15, radius=10, Qd=diagm( [0.1;0.1;0.05].^2 ), postpose_cb=cb, solverParams=SolverParams(useMsgLikelihoods=false))

## # test slew in x

lastpose = sortDFG(ls(fg))[end]
@test isapprox( getPPE(fg, lastpose, :simulated).suggested , [20,0,1.465088], atol=0.001 )

##

fg = RoME.generateCanonicalFG_Helix2DSpiral!(200, rate_r=0.6, rate_a=6, radius=100, solverParams=SolverParams(graphinit=false))

##
end

@testset "test extending fg with helix generator" begin
##

fg = generateCanonicalFG_Helix2D!(5, posesperturn=15, radius=10, solverParams=SolverParams(graphinit=false))

@test !getSolverParams(fg).graphinit

# check the helix is being constructed in a consistent way, using ppe solveKey :simulated
vars = sortDFG(ls(fg))
ppes = [
  [0.0, 0.0, 1.5707963267948966],
  [0.8645454235739924, 4.067366430758004, 1.151917276019672],
  [3.3086939364114176, 7.431448254773942, 0.7330382545911657],
  [6.909830056250526, 9.510565162951536, 0.31415923447063226],
  [11.045284632676536, 9.945218953682733, -0.10471978645923721],
]

for (i,v) in enumerate(vars)
  @test isapprox(ppes[i], getPPE(fg, v, :simulated).suggested; atol=1e-5 )
end

# check that the graph can be expanded with the same generator function

generateCanonicalFG_Helix2D!(5, dfg=fg, posesperturn=15, radius=10, solverParams=SolverParams(graphinit=false))
# check that no new variables were added
@test length(ls(fg)) == length(vars)
@test !exists(fg, :x5)

##

generateCanonicalFG_Helix2D!(6, dfg=fg, posesperturn=15, radius=10, solverParams=SolverParams(graphinit=false))
# check that only one variable has been added
@test length(ls(fg)) == length(vars) + 1
@test exists(fg, :x5)

# check the simulated location of the new variable is also at the right location
push!(vars, :x5)
push!(ppes, [ 15.0;8.660254037844387;-0.5235988055902416] )

@test isapprox(ppes[end], getPPE(fg, :x5, :simulated).suggested; atol=1e-5 )

##
end

##

# using TensorCast
# using Gadfly
# Gadfly.set_default_plot_size(35cm,25cm)

# ##

# vals = getPPE.(fg, (ls(fg) |> sortDFG), :simulated) .|> x-> x.suggested
# @cast XYT[d,p] := vals[p][d]

# Gadfly.plot(x=XYT[1,:],y=XYT[2,:], Geom.path)

##


text:
using Elemental
using Test
using LinearAlgebra: eigvals

@testset "generel eigenvalues (Schur) with eltype: $elty" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
    n = 10
    A = Elemental.DistMatrix(elty)
    Elemental.gaussian!(A, n, n)
    elvals = Elemental.eigvalsGeneral(A)
    lavals = eigvals(Array(A))
    @test sort(abs.(vec(Array(elvals)))) ≈ sort(abs.(lavals))
end


text:
"""
    AdaptiveParameters <: AbstractAdaptiveParameters

A composite type for the parameters of an adaptive [`AbstractRungeKuttaSolver`](@ref).

# Constructors
```julia
AdaptiveParameters(; δ = 0.0, ϵ = 1e-5, K = 100)
```

# Arguments
- `δ :: Real`    : absolute tolerance.
- `ϵ :: Real`    : relative tolerance.
- `K :: Integer` : maximum number of iterations.

# Methods
- [`show`   ](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct AdaptiveParameters{δ_T, ϵ_T, K_T} <: AbstractAdaptiveParameters
    δ::δ_T
    ϵ::ϵ_T
    K::K_T
end

AdaptiveParameters(; δ = 0.0, ϵ = 1e-5, K = 100) = AdaptiveParameters(δ, ϵ, K)

############################################################################################
#                                         PRINTING                                         #
############################################################################################

"""
    show(io::IO, adaptive::AbstractAdaptiveParameters)

prints a full description of `adaptive` and its contents to a stream `io`.
"""
Base.show(io::IO, adaptive::AbstractAdaptiveParameters) = NSDEBase._show(io, adaptive)

"""
    summary(io::IO, adaptive::AbstractAdaptiveParameters)

prints a brief description of `adaptive` to a stream `io`.
"""
Base.summary(io::IO, adaptive::AbstractAdaptiveParameters) = NSDEBase._summary(io, adaptive)

############################################################################################
#                                         FUNCTIONS                                        #
############################################################################################

function adaptive_step!(
    solution::AbstractRungeKuttaSolution,
    solver::AbstractRungeKuttaSolver,
    adaptive::AbstractAdaptiveParameters,
    cache::AbstractRungeKuttaCache
)
    @↓ u = solution
    @↓ h = solver.stepsize
    @↓ n, m, v = cache
    k = solution.k isa Nothing ? cache.k : solution.k[n]
    @↓ s, b, p, d, q = solver.tableau
    @↓ δ, ϵ, K = solver.adaptive
    zero!(v)
    for i = 1:s
        @. v += (b[i] - d[i]) * k[i]
    end
    @. v /= δ + max(abs(u[n]), abs(u[n+1])) * ϵ
    error = norm(v) / √length(v)
    power = - 1 / (min(p, q) + 1)
    factor = 0.9 * error ^ power
    h *= max(0.5, min((m == 1 ? 2.0 : 1.0), factor))
    if error < 1 || m ≥ K
        n += 1
    else
        m += 1
    end
    @↑ cache = n, m
    @↑ solver.stepsize = h
    return solution
end

function adaptive_step!(
    solution::AbstractRungeKuttaSolution,
    solver::AbstractRungeKuttaSolver,
    adaptive::Nothing,
    cache::AbstractRungeKuttaCache
)
    cache.n += 1
    return solution
end

function adaptive_step!(
    solution::AbstractRungeKuttaSolution,
    solver::AbstractRungeKuttaSolver,
    cache::AbstractRungeKuttaCache
)
    return adaptive_step!(solution, solver, solver.adaptive, cache)
end


text:
#
# Compute Polar Decomposition via SVD
#

struct SVDAlg <: PolarAlg end

function solve!(alg::SVDAlg,
                X::Matrix{T}, U::Matrix{T}, H::Matrix{T}) where {T}

    F = svd(X, full = false)
    PQt = Array{T}(undef, size(U))
    mul!(PQt, F.U, transpose(F.V))
    copyto!(U, PQt)
    copyto!(H, F.V * Diagonal(F.S) * F.Vt)
    H = (1/T(2)) * (H + H')

    return SVDResult(U, H)
end


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

# triple-quote delimited strings
@test """abc""" == "abc"
@test """ab"c""" == "ab\"c"
@test """ab""c""" == "ab\"\"c"
@test """ab"\"c""" == "ab\"\"c"
@test """abc\"""" == "abc\""
n = 3
@test """$n\n""" == "$n\n"
@test """$(n)""" == "3"
@test """$(2n)""" == "6"
@test """$(n+4)""" == "7"
@test """$("string")""" == "string"
a = [3,1,2]
@test """$(a[2])""" == "1"
@test """$(a[3]+7)""" == "9"
@test """$(floor(Int,4.5))""" == "4"
nl = "
"
@test """
     a
     b

     c
     """ == "a$(nl)b$(nl)$(nl)c$(nl)"
@test """
      """ == ""
@test """x
     a
    """ == "x$(nl) a$(nl)"
@test """
     $n
   """ == "  $n$(nl)"
@test """
      a
     b
       c""" == " a$(nl)b$(nl)  c"
# tabs + spaces
@test """
	 a
	 b
	""" == " a$(nl) b$(nl)"
@test """
      a
       """ == "a$(nl) "
s = "   p"
@test """
      $s""" == "$s"
@test """
       $s
      """ == " $s$(nl)"
@test """\t""" == "\t"
@test """
      \t""" == "\t"
@test """
      foo
      \tbar""" == "foo$(nl)\tbar"
@test """
      foo
      \tbar
      """ == "foo$(nl)\tbar$(nl)"
@test """
      foo
      bar\t""" == "foo$(nl)bar\t"
@test """
      $("\n      ")
      """ == "\n      $(nl)"

# issue #34105
let str = "   yoyo"
    interp = """interpolate:\n$str
    next line"""
    @test interp == "interpolate:\n   yoyo\nnext line"
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
"""
    *(a, b)

A more performant version of `a * b` where a, b are matrices
"""
function Base.:(*)(a::Array{DoubleFloat{T},2}, b::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end

    result_rows, result_cols = arows, bcols
    result = reshape(Array{DoubleFloat{T}, 1}(undef, result_rows*result_cols), (result_rows, result_cols))

    for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds result[arow,bcol] = asum
        end
    end

    return result
end

function Base.:(*)(a::Array{DoubleFloat{T},2}, b::Array{T,2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end

    result_rows, result_cols = arows, bcols
    result = reshape(Array{DoubleFloat{T}, 1}(undef, result_rows*result_cols), (result_rows, result_cols))

    for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds result[arow,bcol] = asum
        end
    end

    return result
end

function Base.:(*)(a::Array{T,2}, b::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end

    result_rows, result_cols = arows, bcols
    result = reshape(Array{DoubleFloat{T}, 1}(undef, result_rows*result_cols), (result_rows, result_cols))

    for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds result[arow,bcol] = asum
        end
    end

    return result
end

"""
    mul!(c, a, b)
A more performant version of `c = a * b` where a, b, c are matrices
"""
function LinearAlgebra.mul!(c::Array{DoubleFloat{T},2}, a::Array{DoubleFloat{T},2}, b::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end
    if size(c,1) != arows || size(c,2) != bcols
        throw(DimensionMismatch("result c has dimensions $(size(c)), needs ($arows,$bcols)"))
    end

    @inbounds for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds c[arow,bcol] = asum
        end
    end

    return c
end

function LinearAlgebra.mul!(c::Array{DoubleFloat{T},2}, a::Array{DoubleFloat{T},2}, b::Array{T,2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end
    if size(c,1) != arows || size(c,2) != bcols
        throw(DimensionMismatch("result c has dimensions $(size(c)), needs ($arows,$bcols)"))
    end

    @inbounds for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds c[arow,bcol] = asum
        end
    end

    return c
end

function LinearAlgebra.mul!(c::Array{DoubleFloat{T},2}, a::Array{T,2}, b::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    arows, acols = size(a)
    brows, bcols = size(b)
    if  acols != brows
        throw(DimensionMismatch(string("size(a,1) != size(b,2): ",size(a,1), " != ", size(b,2))))
    end
    if size(c,1) != arows || size(c,2) != bcols
        throw(DimensionMismatch("result c has dimensions $(size(c)), needs ($arows,$bcols)"))
    end

    @inbounds for bcol = 1:bcols
        for arow=1:arows
            asum = zero(DoubleFloat{T})
            for acol=1:acols
                @inbounds asum = fma(a[arow,acol], b[acol,bcol], asum)
            end
            @inbounds c[arow,bcol] = asum
        end
    end

    return c
end


text:
## Test pvalue adjustment methods ##
module Test_pval_pi0_adjustment

using MultipleTesting
using Base.Test

pval1 = [0.0, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.2, 0.4, 0.75, 1.0]
pi0 = 0.4

## compared against:
## qvalue::qvalue(p, pi0, pfdr = TRUE, pi0.method = "bootstrap")

## benjamini_hochberg
m = benjamini_hochberg
t = BenjaminiHochbergAdaptive
ref0 = [0.0, 0.0005, 0.003333333, 0.025, 0.1, 0.166666667, 0.285714286, 0.5, 0.833333333, 1.0]
ref = ref0 .* pi0
println(" ** ", m)
@test_throws MethodError m()
## no integers as input
@test_throws MethodError m([0, 1])
## no valid p-values as input
@test_throws DomainError m([-1.0, 0.5], pi0)
@test_throws DomainError m([0.5, 1.5], pi0)
@test_throws DomainError m([0.5, 0.7], -1.0)
@test_throws DomainError m([0.5, 0.7], 1.5)
## single p-value is returned unchanged
pval = rand(1)
@test m(pval, pi0) == pval .* pi0
@test adjust(pval, t(pi0)) == pval .* pi0
## compare with reference values
@test_approx_eq_eps m(pval1, pi0) ref 1e-6
@test_approx_eq_eps adjust(pval1, t(pi0)) ref 1e-6
## BH Adaptive same as BH for π0 missing or 1
@test_approx_eq_eps m(pval1, 1.0) ref0 1e-6
@test_approx_eq_eps adjust(pval1, t(1.0)) ref0 1e-6
@test_approx_eq_eps adjust(pval1, t()) ref0 1e-6
@test_approx_eq_eps adjust(pval1, t(0.0)) zeros(ref0) 1e-6

# adaptive with pi0 estimator == oracle with estimated pi0
@test adjust(pval1, t(Storey())) == adjust(pval1, t(estimate_pi0(pval1, Storey())))
@test adjust(pval1, t(Storey(0.3))) == adjust(pval1, t(estimate_pi0(pval1, Storey(0.3))))
@test adjust(pval1, t(LeastSlope())) == adjust(pval1, t(estimate_pi0(pval1, LeastSlope())))

## qvalue
m = qValues
ref = [0.0, 0.0002, 0.00133333, 0.01, 0.04, 0.0666667, 0.114286, 0.2, 0.333333, 0.4]
println(" ** ", m)
@test_throws MethodError m()
## no integers as input
@test_throws MethodError m([0, 1])
## no valid p-values as input
@test_throws DomainError m([-1.0, 0.5], pi0)
@test_throws DomainError m([0.5, 1.5], pi0)
@test_throws DomainError m([0.5, 0.7], -1.0)
@test_throws DomainError m([0.5, 0.7], 1.5)
## single p-value is returned unchanged
pval = rand(1)
@test m(pval, pi0) == pval .* pi0
## compare with reference values
@test_approx_eq_eps m(pval1, pi0) ref 1e-6
@test_approx_eq_eps m(pval1, pi0, false) ref 1e-6

## qvalue pfdr
m = qValues
ref = [0.099685, 0.099685, 0.099685, 0.099685, 0.099685, 0.102356, 0.128033, 0.201217, 0.333333, 0.4]
println(" ** ", m)
@test_throws MethodError m()
## no integers as input
@test_throws MethodError m([0, 1])
## no valid p-values as input
@test_throws DomainError m([-1.0, 0.5], pi0)
@test_throws DomainError m([0.5, 1.5], pi0)
@test_throws DomainError m([0.5, 0.7], -1.0)
@test_throws DomainError m([0.5, 0.7], 1.5)
## single p-value is returned unchanged
pval = rand(1)
@test m(pval, pi0) == pval .* pi0
## compare with reference values
@test_approx_eq_eps m(pval1, pi0, true) ref 1e-6

end


text:
module GroundData
"""
Module to obtain data from measurement stations located in various regions of Italy including: Alto Adige, Friuli Venezia Giulia, Lombardia, Trentino and Veneto.
"""
#=
colonna	     |   descrizione
------------------------------------------------------------------------------------
uiid	     |   identificatore univoco delle stazioni di qualunque genere di misura
------------------------------------------------------------------------------------
parameter	 |   tipo di parametro misurato
------------------------------------------------------------------------------------
unit	     |   unita di misura (possibilmente SI)
------------------------------------------------------------------------------------
value	     |   valore misurato
------------------------------------------------------------------------------------
date	     |   anno mese giorno ora (UTM)
------------------------------------------------------------------------------------
longitude	 |   longitudine della stazione
------------------------------------------------------------------------------------
latitude	 |   latitudine della stazione
------------------------------------------------------------------------------------
quote	     |   quota stazione || quota stazione più quota misura per il vento
------------------------------------------------------------------------------------
validation	 |   bool (già segnalato dalla stazione)
------------------------------------------------------------------------------------
note	     |   errori e outlayers? altro?
------------------------------------------------------------------------------------
=#

using CSV
using DataFrames
using Dates
using UUIDs

str = occursin( "GroundData.jl", @__FILE__ ) ? "" : "src\\"
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataAA.jl")
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataER.jl")
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataFVG.jl")
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataL.jl")
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataT.jl")
include("$(@__DIR__)\\$(str)GroundDataRegions\\GroundDataV.jl")

export getGroundData

@enum Region AA=1 FVG=2 L=3 T=4 V=5

const regions_modules = Module[
    GroundDataAA,
    GroundDataFVG,
    GroundDataL,
    GroundDataT,
    GroundDataV
]

Base.convert(::Type{Int64}, s::AbstractString ) = parse( Int64, s )
Base.convert(::Type{Float64}, s::AbstractString ) = isempty(s) || s == "NULL" ? missing : parse( Float64, s )
Base.convert(::Type{DateTime}, s::AbstractString ) = isempty(s) ? missing : DateTime(s)

"""
    pop!( df::DataFrame, columns::Vector{Symbol}, condition::Function  )

Delete rows from `df` where `condition` on the elements of `columns` is true and return the deleted rows as a new `DataFrame`
"""
function pop!( df::DataFrame, columns::Vector{Symbol}, condition::Function  )
    res = filter( x -> any( condition, x[columns] ), df )
    filter!( x -> !any( condition, x[columns] ), df )
    return res
end

"""
    pop!( df::DataFrame, column::Symbol, condition::Function  )

Delete rows from `df` where `condition` on the elements of `column` is true and return the deleted rows as a new `DataFrame`
"""
function pop!( df::DataFrame, column::Symbol, condition::Function  )
    res = filter( x -> condition(x[column]) == true, df )
    filter!( x -> condition(x[column]) == false, df )
    return res
end

"""
    createMap( attributes::AbstractVector, destinations::AbstractVector[; n::Int64=0] )

Create an array of pairs from the tow input vectors, skipping the last `n` pairs
"""
function createMap( attributes::AbstractVector, destinations::AbstractVector )
    min_len = min( length.( [attributes, destinations] )... )
    return [ attributes[i] => destinations[i] for i in 1:min_len ]
end

"""
    standarize(

Generate a dataframe in standard format from `dfSta` and `dfSen` using `bridge` to join the two dataframes and `map` to select the desired columns
and map them to the standard ones.
If `dfSen` is not provided the selection of the column will be done on `dfSta`.
"""
function standardize( map::AbstractVector, dfSta::DataFrame, dfSen::Union{DataFrame, Nothing}=nothing, bridge::Union{ Nothing, Symbol, Pair{Symbol, Symbol} }=nothing )
 # Separate the column that have a mapping from the others
    complete_map = []
    missing_map = []
    for (x,y) in map
        if isnothing(x)
            push!( missing_map, y )
        else
            push!( complete_map, x => y )
        end
    end
 # Check wether there is a second dataframe or there is only one containing all the needed informations
    if isnothing(dfSen) || isnothing(bridge)
        dataframe = select( dfSta, complete_map... )
    else
        dataframe = innerjoin( dfSta, dfSen, on=bridge )
        select!( dataframe, complete_map... )
    end
 # For each column that doesn't have a mapping insert the column as an array of missings
    insertcols!( dataframe, ( missing_map .=> Ref( missings( nrow(dataframe) ) ) )... )
    return dataframe
end

"""
    generateUuidsTable()

Generate a CSV file and return a `DataFrames.DataFrame` containing informations on all the stations of all the regions and their uuids and generate a second file containing
the information on the stations that lack spatial coordinates.
"""
function generateUuidsTable()
    columns = [ :local_id, :name, :longitude, :latitude ]
    df = DataFrame(
        :local_id => Union{Missing, Any}[],
        :name => String[],
        :longitude => Union{Missing, Float64}[],
        :latitude => Union{Missing, Float64}[],
        :type => Symbol[],
        :region => Region[]
    )
    # Obtain data of each "Type" from each "Region", if available
    for r in instances(Region), t in [:METEO, :AIRQUALITY]
        rgn = regions_modules[Integer(r)]
        # Get the data if possible
        station = try
                      rgn.getData( type=t )
                  catch e
                      println( "Skipped $r $t, due to:\n", e )
                      continue
                  end
        stat_info = rgn.getRegionStationsInfo(t)
        map = createMap( stat_info, columns )
        res = standardize( map, station )
        insertcols!( res, :type => t, :region => r )
        append!( df, res )
    end
    # AA e V hanno delle stazioni ( 7 e 2 rispettivamente, che non hanno latitudine e longitudine )
    mdf = pop!( df, [:latitude, :longitude], ismissing )
    CSV.write(split(@__DIR__, "src")[1]*"\\resources\\Ground stations data\\missing_stazioni.csv", mdf)
    disallowmissing!(df)
    unique!( df, [:longitude, :latitude] )
    insertcols!( df, :uuid => [ uuid4() for i in 1:nrow(df) ] )
    CSV.write(split(@__DIR__, "src")[1]*"\\resources\\Ground stations data\\stazioni.csv", df)
    return df
end

#   uuids = generateUuidsTable()

"""
    getGroundData( type::Symbol=METEO, regions::Region... )

Obtain `type` data of `regions` from `source` and save it as `filePath`
"""
function getGroundData( type::Symbol=:METEO, regions::Region... )

 # `s` will be used to regulate the number of columns of the df based on the data type, the conditional assignment
    # checks also wether `type` is valid
    s = type == :METEO ? 0 :
          type == :AIRQUALITY ? 1 : throw( DomainError( type, "`type` must be either `:METEO` OR `:AIRQUALITY`" ) )

    columns = Pair{Symbol, Vector}[
        :uuid                   => String[],
        :parameter              => String[],                    # kind of the measured parameter (es. temperature)
        :unit                   => Union{String, Missing}[],    # unity of measure
        :value                  => Union{Float64, Missing}[],   # measurement value
        :frequency              => Union{String, Missing}[],    # frequency of measurements
        :date                   => Union{DateTime, Missing}[],  # date of the measurement (year, month, day, time) (UTM)
        :longitude              => Float64[],                   # station longitude (longitude of the measurement)
        :latitude               => Float64[],                   # station latitude (latitude of he measurement)
        :height                 => Union{Float64, Missing}[],   # station height (acounting for measurement height for wind measurements)
        :validation             => Any[],                       # wether the measurements have been validated
        :note                   => Any[],                       # errors, outlyers, ecc.
        :rel_measurement_height => String[]                     # relative height of measurement
    ]
 # There are no `:rel_measurement_height` columns for the airquality data so if that is the type of data
    # the created dataframe will lack the specific column
    df = DataFrame( columns[1:end-s] )
    mdf = DataFrame( columns[2:end-s] )
    allowmissing!( mdf, [:longitude, :latitude] )
    for region in regions
     # Use the Region enum number to obtain the corresponding module
        rnum = Integer(region)
        rgn = regions_modules[rnum]
     # Obtain the dataframes that will be used to create the standard format dataframe, while checking for
        # possible problems
        resSta = try
                    rgn.getData( type=type, kind=:STATIONS )
                 catch e
                    println( "Skipped $region, due to:\n", e )
                    continue
                 end
        resSen = try
                    rgn.getData( type=type, kind=:SENSORS )
                 catch e
                    println( "Skipped $region, due to:\n", e )
                    continue
                 end
        # Create the standard dataframe
        attributes = rgn.getRegionAttributes(type)
        map = createMap( attributes, Symbol.(names(df)[2:end]) )
        bridge = rgn.getRegionIds(type)
        res = standardize( map, resSta, resSen, bridge )
        mres = pop!( res, [:latitude, :longitude], ismissing )
        path = split(@__DIR__, "src")[1]*"\\resources\\Ground stations data\\stazioni.csv"
        uuids = isfile(path) ? CSV.read(path, DataFrame) : generateUuidsTable()
        res = innerjoin( uuids[:, [:uuid, :longitude, :latitude]], res, on=[:longitude, :latitude] )
        append!( df, res )
        append!( mdf, mres )
    end
    return df, mdf
end

#   resmt, mresmt = getGroundData( :METEO, AA, FVG, L, T, V )
#   resaq, mresaq = getGroundData( :AIRQUALITY, AA, FVG, L, T, V )

"""
    saveGroundData( path::AbstractString, data::DataFrame; overwrite::Bool=false )

Save `DataFrames.DataFrame` `data` as `path` CSV file.
If `overwrite` is set to `true` adds `data` to the file instead.
"""
function saveGroundData( path::AbstractString, data::DataFrame; overwrite::Bool=false )
    condition = !overwrite && basename(path) in readdir(dirname(path))
    if condition
        # Attributes to confront in order to see if the entry is a duplicate
        attrs = [:uuid, :value, :date, :height, :rel_measurement_height]
        # Preexisting data
        old_data = CSV.read( path, DataFrame )
        filter!( row -> row[attrs] in old_data[:, attrs], data )
    end
    CSV.write( path, data, append=condition )
end

#    path = *( @__DIR__, "\\..\\Dati stazioni")
#    saveGroundData( path*"\\data.csv", resmt )
#    saveGroundData( path*"\\missing_data.csv", mresmt )

end # module

using Plots
using Shapefile
using GeoStats

const shp = Shapefile

stations = shp.Table("C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\Stazioni_Veneto\\Stazioni_Veneto.shp")
comuni = shp.Table("C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis test data\\c0104011_Comuni\\c0104011_Comuni.shp")

resmt, mresmt = getGroundData( :METEO, AA, FVG, L, T, V )
resaq, mresaq = getGroundData( :AIRQUALITY, AA, FVG, L, T, V )

path = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Ground stations data\\25-03-22 - Post incendio"
saveGroundData(path*"\\meteo_data.csv", resmt)
saveGroundData(path*"\\airquality_data.csv", resaq)
saveGroundData(path*"\\airquality_data_missing.csv", mresaq)

unique(resmt.parameter)

# All relative humidity measurements
measurements = resmt[ resmt.parameter .== "Umidità relativa", : ]
# All the measurement's values
values = ( humidity = measurements.value ./ 100, )
# All the coordinates of the measurements
coord = Tuple{Float64, Float64}[ (c[1], c[2]) for c in eachrow(measurements[:, [:longitude, :latitude]]) ]
# Georeference data
D = georef(values, coord)
# Estimation domain
G = CartesianGrid(100, 100)
# Estimation problem
problem = EstimationProblem(D, G, :humidity)
# Solver from the list of solvers
solver = Kriging( :val => ( variogram = GaussianVariogram(range=35.0), ) )
# Solving problem
solution = solve(problem, solver)
# Solution plot
contourf(solution, clabels=true)

=#

max_stars_repo_path: src/Data Gathering/GroundData.jl

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_path: src/Data Gathering/GroundData.jl

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_path: src/Data Gathering/GroundData.jl

max_forks_repo_name: DavideFavaro/Tirocinio