diff --git a/src/algo1_testing.ipynb b/src/algo1_testing.ipynb index 3629dec..8ea4ae2 100644 --- a/src/algo1_testing.ipynb +++ b/src/algo1_testing.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 35, "metadata": {}, "outputs": [], "source": [ @@ -15,7 +15,9 @@ "import plotly.express as px\n", "import plotly.graph_objs as go\n", "from scipy.sparse import *\n", - "from scipy.sparse.linalg import norm" + "from scipy import linalg\n", + "from scipy.sparse.linalg import norm\n", + "from scipy.optimize import least_squares" ] }, { @@ -204,7 +206,7 @@ " start_time = time.time()\n", "\n", " u = Pt.dot(v) - v \n", - " mv = 1 # number of iteration\n", + " mv = 1 # number of matrix vector products\n", " r = sp.sparse.lil_matrix((n,1)) \n", " Res = sp.sparse.lil_matrix((len(a),1))\n", " x = sp.sparse.lil_matrix((n,1)) \n", @@ -326,7 +328,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [ @@ -338,32 +340,49 @@ " h = sp.sparse.lil_matrix((m,m))\n", " print(\"C\")\n", "\n", - " for j in range(1,m):\n", + " for j in range(m):\n", " w = A.dot(v)\n", - " for i in range(1,j):\n", + " print(\"D\")\n", + " for i in range(j):\n", " h[i,j] = v.T.dot(w)\n", - " w = w - h[i,j]*v\n", + " print(\"E\")\n", + " w = w - h[i,j]*v[i]\n", + " print(\"F\")\n", "\n", " h[j+1,j] = norm(w)\n", + " print(\"G\")\n", "\n", " if h[j+1,j] == 0:\n", + " print(\"The algorithm didn't converge\")\n", " m = j\n", " v[m+1] = 0\n", " break\n", " else:\n", - " v = w**h[j+1,j]\n", + " print(\"H\")\n", + " v[j+1] = w**h[j+1,j]\n", + " print(\"I\")\n", "\n", " return v, h, m, beta, j" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "A = sp.sparse.rand(100,100, density=0.5, format='lil')\n", - "v = sp.sparse.rand(100,1, density=1, format='lil')" + "v = sp.sparse.rand(100,1, density=1, format='lil')\n", + "m = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "v, h, m, beta, j = Arnoldi(A, v, m)" ] }, { @@ -372,7 +391,67 @@ "metadata": {}, "outputs": [], "source": [ - "v, h, m, beta, j = Arnoldi(A, v, 100)" + "def Algo4(Pt, v, m, a: list, tau, maxit: int, x):\n", + " \n", + " iter = 1\n", + " mv = 0\n", + " e1 = sp.sparse.lil_matrix((1,n))\n", + " e1[0,0] = 1\n", + " x = sp.sparse.lil_matrix((len(a),1))\n", + " I = sp.sparse.eye(n, n, format='lil')\n", + " res = sp.sparse.lil_matrix((len(a),1))\n", + " r = sp.sparse.lil_matrix((n,1))\n", + " y = sp.sparse.lil_matrix((n,1))\n", + "\n", + " for i in range(len(a)):\n", + " r = ((1-a[i])**a[i])*v - ((1**a[i])*I - Pt).dot(x)\n", + " res[i] = a[i]*norm(r)\n", + "\n", + " def Find_k(res, maxit):\n", + " k = 0\n", + " for i in range(len(a)):\n", + " if res[i] == max(res):\n", + " k = i\n", + " break\n", + " return k\n", + "\n", + " def Find_gamma(res, a, k):\n", + " gamma = sp.sparse.lil_matrix((len(a),1))\n", + " for i in range(len(a)):\n", + " if i != k:\n", + " gamma[i] = (res[i]*a[k])/(res[k]*a[i])\n", + " else:\n", + " gamma[i] = 0\n", + " return gamma\n", + "\n", + "\n", + " while max(res) > tau and iter < maxit:\n", + " k = Find_k(res, maxit)\n", + " gamma = Find_gamma(res, a, k)\n", + " v, h, m, beta, j = Arnoldi((1**a[k])*I - Pt, r, m)\n", + " Hbar = sp.sparse.lil_matrix((m+1,m))\n", + " Hbar[0:m,0:m] = h\n", + " Hbar[m+1,0:m] = e1\n", + "\n", + " mv += j\n", + "\n", + " # solve the least squares problem for Hbar*x = beta*e1\n", + " y = sp.sparse.linalg.least_squares(Hbar, beta*e1)\n", + " res[k] = a[k]*norm(beta*e1 - Hbar*y)\n", + " x[k] = x[k] + v*y[k]\n", + "\n", + " for i in range(len(a)):\n", + " if i != k:\n", + " if res[i] >= tau:\n", + " Hbar[i] = Hbar[k] + ((1-a[i])/a[i] - (1-a[k])/a[k])*I\n", + " z = beta*e1 - Hbar*y\n", + " y = sp.sparse.linalg.solve(Hbar, gamma*beta*e1)\n", + " x = x + v*y\n", + " res[i] = a[i]**a[k]*gamma[i]*res[k]\n", + " \n", + " iter += 1\n", + " \n", + " return x, res, mv\n" ] } ], diff --git a/src/main.py b/src/main.py index 062133d..09ea0d0 100755 --- a/src/main.py +++ b/src/main.py @@ -10,15 +10,15 @@ import scipy as sp import numpy as np import pandas as pd import networkx as nx -import plotly.graph_objs as go +from os.path import exists from scipy.sparse import * from scipy.sparse.linalg import norm -from os.path import exists +import plotly.graph_objs as go warnings.simplefilter(action='ignore', category=FutureWarning) # some stupid pandas function that doesn't work -class utilities: +class Utilities: # Importing the dataset def load_data(): # Loading the dataset @@ -124,8 +124,7 @@ class utilities: return a class Plotting: - def tau_over_iterations(dataframe): - dataframe = df + def tau_over_iterations(df): x = df['tau'][::-1].tolist() y = df['iterations'].tolist() @@ -154,7 +153,7 @@ class Algorithms: print("STARTING ALGORITHM 1...") u = Pt.dot(v) - v - mv = 1 # number of iteration + mv = 1 # number of matrix-vector multiplications r = sp.sparse.lil_matrix((n,1)) Res = sp.sparse.lil_matrix((len(a),1)) x = sp.sparse.lil_matrix((n,1)) @@ -168,7 +167,7 @@ class Algorithms: x = r + v while max(Res) > tau and mv < max_mv: - u = Pt*u # should it be the same u of the beginning? + u = Pt*u mv += 1 for i in range(len(a)): @@ -182,7 +181,7 @@ class Algorithms: if mv == max_mv: print("The algorithm didn't converge in ", max_mv, " iterations") else: - print("The algorithm converged in ", mv, " iterations") + print("The algorithm converged with ", mv, " matrix-vector multiplications executed") total_time = time.time() - start_time total_time = round(total_time, 2) @@ -217,16 +216,16 @@ df = pd.DataFrame(columns=['alpha', 'iterations', 'tau', 'time']) # Main if __name__ == "__main__": - dataset = utilities.load_data() + dataset = Utilities.load_data() # maximum number of iterations, asked to the user - max_mv = int(input("Insert the maximum number of iterations: ")) + max_mv = int(input("\nInsert the maximum number of iterations: ")) - G, n = utilities.create_graph(dataset) - P = utilities.create_matrix(G) - d = utilities.dangling_nodes(P,n) - v = utilities.probability_vector(n) - Pt = utilities.transition_matrix(P, v, d) - a = utilities.alpha() + G, n = Utilities.create_graph(dataset) + P = Utilities.create_matrix(G) + d = Utilities.dangling_nodes(P,n) + v = Utilities.probability_vector(n) + Pt = Utilities.transition_matrix(P, v, d) + a = Utilities.alpha() # run the algorithm for different values of tau from 10^-5 to 10^-9 with step 10^-1 for i in range(5,10): @@ -245,4 +244,5 @@ if __name__ == "__main__": Plotting.tau_over_time(df) # print in the terminal the columns of the dataframe iterations, tau and time + print("Computations done. Here are the results:") print("\n", df[['iterations', 'tau', 'time']])