Metric Nearness.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots, LightGraphs, SparseArrays, SimpleWeightedGraphs\n",
    "using Statistics, BenchmarkTools, LinearAlgebra, ProgressMeter\n",
    "using Distributions, Base.Threads\n",
    "using Base.GC\n",
    "using JuMP, ECOS, CSDP, Ipopt, SDPA, ProxSDP, OSQP, SCS, COSMO\n",
    "plotly();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helpers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     9,
     26,
     30,
     59,
     68
    ]
   },
   "outputs": [],
   "source": [
    "function flipcoin(p)\n",
    "    x = rand()\n",
    "    if x <= p\n",
    "        return 1\n",
    "    else\n",
    "        return 0\n",
    "    end\n",
    "end\n",
    "\n",
    "# Type 1 graphs\n",
    "function genGraph(n,p)\n",
    "    G = randn(n,n)\n",
    "    for i = 1:n\n",
    "        G[i,i] = 0\n",
    "        for j = 1:i-1\n",
    "            if flipcoin(p) == 0\n",
    "                G[i,j] = 0\n",
    "                G[j,i] = 0\n",
    "            else\n",
    "                G[i,j] = G[j,i] \n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return G\n",
    "end  \n",
    "\n",
    "function Labels(A,k)\n",
    "    return Set(partialsortperm(A, 1:k, rev=true))\n",
    "end\n",
    "\n",
    "# Type 2 graphs\n",
    "function genGraphDir(n,p,alpha,k)\n",
    "    d = Dirichlet(alpha)\n",
    "    L = Array{Set{Int},1}(undef,n)\n",
    "    for i = 1:n\n",
    "        L[i] = Labels(rand(d),k)\n",
    "    end\n",
    "    \n",
    "    G = zeros(n,n)\n",
    "    for i = 1:n\n",
    "        G[i,i] = 0\n",
    "        for j = 1:i-1\n",
    "            if flipcoin(p) == 0\n",
    "                G[i,j] = 0\n",
    "                G[j,i] = 0\n",
    "            else\n",
    "                if isempty(intersect(L[i],L[j]))\n",
    "                    G[i,j] = 1\n",
    "                    G[j,i] = 1\n",
    "                else\n",
    "                    G[i,j] = 0.00000001\n",
    "                    G[j,i] = 0.00000001\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return G,L\n",
    "end  \n",
    "\n",
    "function enumerate_paths2(s)\n",
    "    P = Array{Any,1}(undef,size(s.parents, 1))\n",
    "    for v = 1:size(s.parents, 1)\n",
    "        P[v] = LightGraphs.enumerate_paths(s, v)\n",
    "    end\n",
    "    \n",
    "    return P\n",
    "end\n",
    "\n",
    "function calculateSTDEV(x,y)\n",
    "    error = 0.0\n",
    "    for i=1:length(y)\n",
    "        error = error + (x-y[i])^2\n",
    "    end\n",
    "    \n",
    "    return sqrt(error/length(y))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running experiment for standard solvers. \n",
    "\n",
    "Modify the line set_optimizer(model,*.optimizer) where * is the name of the optimizer. For the commerical solvers you must have a liscence and those solvers independetly installed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     1
    ]
   },
   "outputs": [],
   "source": [
    "T = zeros(50)\n",
    "for s = 1:3\n",
    "    n = 100*s\n",
    "    time = zeros(5)\n",
    "    for r = 1:5 \n",
    "        G = 0\n",
    "        model = 0\n",
    "        D = 0\n",
    "        G = genGraph(n,1.1)\n",
    "        \n",
    "        t = @elapsed model = Model()\n",
    "        t += @elapsed @variable(model, D[1:n,1:n]);\n",
    "        t += @elapsed @objective(model, Min, sum((D-G).^2));\n",
    "        t += @elapsed for i = 1:n\n",
    "            for j = 1:n\n",
    "                @constraint(model, D[i,j] == D[j,i])\n",
    "            end\n",
    "        end\n",
    "                \n",
    "        t += @elapsed for i = 1:n\n",
    "            for j = 1:n\n",
    "                for k = 1:n\n",
    "                    @constraint(model, D[i,j] <= D[i,k]+D[k,j])\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "        \n",
    "        t += @elapsed for i = 1:n\n",
    "            @constraint(model, D[i,i] ==  0)\n",
    "        end\n",
    "        \n",
    "        t += @elapsed set_optimizer(model, CSDP.Optimizer)\n",
    "        t += @elapsed optimize!(model) \n",
    "        print(termination_status(model))\n",
    "        time[r] = solve_time(model)\n",
    "        @show(t)\n",
    "        flush(stdout)\n",
    "    end\n",
    "    T[s] = mean(time)\n",
    "    @show(T[s])\n",
    "    flush(stdout)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Our implementation of the Brickell et al. algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "function Bregman3Cycle(D)\n",
    "    (n,n) = size(D)\n",
    "    Z = zeros(n,n,n)\n",
    "    maxD = 2\n",
    "    count = 0\n",
    "    while(maxD > 1e-10)\n",
    "        maxD = 0\n",
    "        for i = 1:n\n",
    "            for j = 1:i-1\n",
    "                for k = 1:j-1\n",
    "                    d = (-1*D[i,k] + D[i,j] + D[j,k])/3 \n",
    "                    c = min(d,Z[i,j,k])\n",
    "                    Z[i,j,k] -= c\n",
    "                    D[i,k] += c\n",
    "                    D[k,i] += c\n",
    "                    D[i,j] -= c\n",
    "                    D[j,i] -= c\n",
    "                    D[j,k] -= c\n",
    "                    D[k,j] -= c\n",
    "                    d = (-1*D[j,i] + D[j,k] + D[k,i])/3 \n",
    "                    c = min(d,Z[j,k,i])\n",
    "                    Z[j,k,i] -= c\n",
    "                    D[i,k] -= c\n",
    "                    D[k,i] -= c\n",
    "                    D[i,j] += c\n",
    "                    D[j,i] += c\n",
    "                    D[j,k] -= c\n",
    "                    D[k,j] -= c\n",
    "                    d = (-1*D[k,j] + D[k,i] + D[i,j])/3 \n",
    "                    c = min(d,Z[k,i,j])\n",
    "                    Z[k,i,j] -= c\n",
    "                    D[i,k] -= c\n",
    "                    D[k,i] -= c\n",
    "                    D[i,j] -= c\n",
    "                    D[j,i] -= c\n",
    "                    D[j,k] += c\n",
    "                    D[k,j] += c\n",
    "                end\n",
    "            end\n",
    "        end \n",
    "        FS =  LightGraphs.floyd_warshall_shortest_paths(SimpleWeightedGraph(D))\n",
    "        maxD = norm(FS.dists-D)\n",
    "        count += 1 \n",
    "    end\n",
    "    return (D,count)\n",
    "end "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Our algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "function BregmanOrig(D)\n",
    "    (n,n) = size(D)\n",
    "    #G = sparse(copy(D))\n",
    "    g = SimpleWeightedGraph(D)\n",
    "    Z = Dict()\n",
    "    Z′ = spzeros(n,n)\n",
    "    maxD = 2\n",
    "    count = 0\n",
    "    while(maxD > 1e-10)\n",
    "        \n",
    "        for p in keys(Z)\n",
    "            z = Z[p]\n",
    "            l = length(p)\n",
    "            u = p[1]\n",
    "            v = p[l]\n",
    "            d = -1*g.weights[u,v]\n",
    "            for i = 1:l-1\n",
    "                d = d + g.weights[p[i], p[i+1]]\n",
    "            end\n",
    "            c = min(d/l,z)\n",
    "            for i = 1:l-1\n",
    "                g.weights[p[i],p[i+1]] -= c\n",
    "                g.weights[p[i+1],p[i]] -= c\n",
    "            end\n",
    "            g.weights[u,v] += c\n",
    "            g.weights[v,u] += c\n",
    "            if z == c\n",
    "                delete!(Z,p)\n",
    "            else\n",
    "                Z[p] -= c\n",
    "            end\n",
    "        end\n",
    "    \n",
    "        for i = 1:n\n",
    "            for j = 1:i-1\n",
    "                c = min(g.weights[j,i] - 1e-14,Z′[j,i])\n",
    "                g.weights[j,i] -= c\n",
    "                g.weights[i,j] -= c\n",
    "                Z′[j,i] -= c\n",
    "                Z′[i,j] -= c             \n",
    "            end\n",
    "        end\n",
    "        \n",
    "        FS = LightGraphs.floyd_warshall_shortest_paths(g)\n",
    "        U = FS.dists\n",
    "        P = enumerate_paths2(FS)\n",
    "        maxD = 0\n",
    "        \n",
    "        for i = 1:n\n",
    "            for j = 1:i-1\n",
    "                if g.weights[j,i] - U[j,i] > 0  \n",
    "                    p = P[j][i] #enumerate_paths2(FS,i,j)\n",
    "                    l = length(p)\n",
    "                    u = p[1]\n",
    "                    v = p[l]\n",
    "                    d = -1*g.weights[u,v]\n",
    "                    for k = 1:l-1\n",
    "                        d = d + g.weights[p[k], p[k+1]]\n",
    "                    end\n",
    "                    if d < 0\n",
    "                        c=d/l\n",
    "                        for k = 1:l-1\n",
    "                            g.weights[p[k],p[k+1]] -= c\n",
    "                            g.weights[p[k+1],p[k]] -= c\n",
    "                        end\n",
    "                        g.weights[p[1],p[l]] += c\n",
    "                        g.weights[p[l],p[1]] += c\n",
    "                        if haskey(Z,p)\n",
    "                            Z[p] = Z[p] - c\n",
    "                        else\n",
    "                            Z[p] = -1*c\n",
    "                        end\n",
    "                        if abs(d) > maxD\n",
    "                            maxD = abs(d)\n",
    "                        end\n",
    "                    end\n",
    "                end   \n",
    "            end\n",
    "        end\n",
    "        FS =  LightGraphs.floyd_warshall_shortest_paths(g)\n",
    "        maxD = norm(FS.dists-g.weights)\n",
    "        \n",
    "        count+=1\n",
    "    end\n",
    "    return length(Z)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The code to run the metric nearness experiemnt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "T = zeros(50);\n",
    "Z = zeros(50);\n",
    "errorbars = zeros(50,1);\n",
    "errorbarsZ = zeros(50,1);\n",
    "T2 = zeros(50);\n",
    "errorbars2 = zeros(50,1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i = 1:16\n",
    "    time = zeros(5)\n",
    "    time2 = zeros(5)\n",
    "    z = zeros(5)\n",
    "    n = 100*i\n",
    "    for j = 1:5 \n",
    "        G = genGraph(100*i,1.1)\n",
    "        #G,L = genGraphDir(100*i,1.1, ones(100),20)\n",
    "        time[j] = @elapsed z[j] = BregmanOrig(copy(G))\n",
    "        @show(time[j])\n",
    "        time2[j] = @elapsed Bregman3Cycle(copy(G))\n",
    "        @show(time2[j])\n",
    "    end\n",
    "    errorbarsZ[i] = calculateSTDEV(Z[i],z)\n",
    "    T[i] = mean(time)\n",
    "    T2[i] = mean(time2)\n",
    "    errorbars[i] = calculateSTDEV(T[i],time)\n",
    "    errorbars2[i] = calculateSTDEV(T2[i],time2)\n",
    "    @show(T[i],T2[i])\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.4.1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}