Repository: EmilianoMarchese/SurpriseMeMore (language: python)

Issue:
Added examples with short guide to the methods.



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 84ad4de6..6c322b59 100644
--- a/README.md
+++ b/README.md
@@ -1,104 +1,61 @@
-Surprisemess
+SurpriseMeMore
 -------------------------------------------------------------------
 
-Surprisemess is a toolbox for detecting mesoscale structure in networks, released as a python3 module. 
+SurpriseMeMore is a toolbox for detecting mesoscale structure in networks, released as a python3 module. 
 
-NEMtropy provides the user with a state of the art solver for a range variety of Maximum Entropy Networks models derived from the ERGM family.
-This module allows you to solve the desired model and generate a number of randomized graphs from the original one: the so-called _graphs ensemble_.
+SurpriseMeMore provides the user with a variety of solvers, based on the _surprise_ framework, for the detection of mesoscale structures ( e.g. communities, core-periphery) in networks.
 
-NEMtropy builds on the current literature on the matter, improving both in speed of convergence and in the scale of the feasible networks.
-To explore Maximum-Entropy modeling on networks, checkout [Maximum Entropy Hub](https://meh.imtlucca.it/).
-
-The models implemented in NEMtropy are presented in a forthcoming [paper](https://arxiv.org/abs/2101.12625) on arXiv.
+The models implemented in SurpriseMeMore are presented in a forthcoming [paper](https://arxiv.org/abs/2106.05055) on arXiv.
 If you use the module for your scientific research, please consider citing us:
 
 ```
-    @misc{vallarano2021fast,
-          title={Fast and scalable likelihood maximization for Exponential Random Graph Models}, 
-          author={Nicolò Vallarano and Matteo Bruno and Emiliano Marchese and Giuseppe Trapani and Fabio Saracco and Tiziano Squartini and Giulio Cimini and Mario Zanon},
-          year={2021},
-          eprint={2101.12625},
-          archivePrefix={arXiv},
-          primaryClass={physics.data-an}
+    @misc{marchese2021detecting,
+      title={Detecting mesoscale structures by surprise}, 
+      author={Emiliano Marchese and Guido Caldarelli and Tiziano Squartini},
+      year={2021},
+      eprint={2106.05055},
+      archivePrefix={arXiv},
+      primaryClass={physics.soc-ph}
     }
 ```
 
 #### Table Of Contents
-- [Currently Implemented Models](#currently-implemented-models)
+- [Currently Implemented Methods](#currently-implemented-methods)
 - [Installation](#installation)
 - [Dependencies](#dependencies)
-- [A Simple Example](#simple-example)
+- [Some Examples](#some-examples)
 - [Development](#development)
-- [Testing](#testing)
 - [Credits](#credits)
 
-## Currently Implemented Models
-The main feature of NEMtropy is (but not limited to) *network randomization*. 
-The specific kind of network to randomize and property to preserve defines the model you need:
-
-* **UBCM** *Undirected Binary Configuration Model* [[1]](#1)
-* **UECM** *Undirected Enhanced Configuration Model* [[1]](#1)
-* **DBCM** *Directed Binary Configuration Model* [[1]](#1)
-* **DECM** *Directed Enhanced Configuration Model* [[1]](#1)
-* **CReMa** [[2]](#2)
-* **BiCM** *Bipartite Configuration Model* [[3]](#3)
-
-The following table may helps you identify the model that fits your needs in function of the type of network you are working with;
-for in-depth discussion please see the references.
-
-[...] | Undirected Graph | Directed Graph | Bipartite Graph
------ | ---------------- | -------------- | --------------
-**Binary Graph** | *UBCM* | *DBCM* | *BiCM*
-**Weighted Graph** | *UECM*, *CReMa*  | *DECM*, *CReMa* | -
-
-The BiCM module is also available as [a standalone package](https://github.com/mat701/BiCM), find its docs [here](https://bipartite-configuration-model.readthedocs.io/en/latest/). 
-
-_References_
-
-* <a id="1">[1]</a>
-    Squartini, Tiziano, Rossana Mastrandrea, and Diego Garlaschelli.
-    "Unbiased sampling of network ensembles."
-    New Journal of Physics 17.2 (2015): 023052.
-    https://arxiv.org/abs/1406.1197
-* <a id="2">[2]</a>
-    Parisi, Federica, Tiziano Squartini, and Diego Garlaschelli.
-    "A faster horse on a safer trail: generalized inference for the efficient reconstruction of weighted networks."
-    New Journal of Physics 22.5 (2020): 053053.
-    https://arxiv.org/abs/1811.09829
-* <a id="3">[3]</a>
-    Saracco, Fabio, Riccardo Di Clemente, Andrea Gabrielli, and Tiziano Squartini.
-	"Randomizing bipartite networks: the case of the World Trade Web." 
-	Scientific reports 5, no. 1 (2015): 1-18.
-    https://doi.org/10.1038/srep10595
+## Currently Implemented Methods
+The available methods, for both directed and undirected networks, are:
 
+* *Community detection on binary networks* 
+* *Community detection on weighted networks with integer weights* 
+* *Community detection on weighted networks with continuous weights* 
+* *Core-Peryphery detection on binary networks* 
+* *Core-Peryphery detection on weighted networks with integer weights*
 
 Installation
 ------------
-
-NEMtropy can be installed via pip. You can get it from your terminal:
+SurpriseMeMore can be installed via pip. You can get it from your terminal:
 
 ```
-    $ pip install NEMtropy
+    $ pip install surprisememore
 ```
 
 If you already install the package and wish to upgrade it,
 you can simply type from your terminal:
 
 ```
-    $ pip install NEMtropy --upgrade
+    $ pip install surprisememore --upgrade
 ```
 
 Dependencies
 ------------
 
-NEMtropy uses <code>numba</code> and <code>powerlaw</code> libraries. They can be installed via pip by running in your terminal the following command:
-
-```
-    $ pip install numba
-    $ pip install powerlaw
-```
-
-For <code>python3.5</code> users the correct command is the following:
+NEMtropy uses <code>numba</code> library. It is installed automatically with surprisememore.
+If you use <code>python3.5</code> you may incur in an error, we suggest installing numba with the following command:
 
 ```
     $ pip install --prefer-binary numba
@@ -107,57 +64,92 @@ For <code>python3.5</code> users the correct command is the following:
 It avoids an error during the installation of <code>llvmlite</code> due to 
 the absence of its wheel in <code>python3.5</code>.
 
-Simple Example
+Some Examples
 --------------
-As an example we solve the UBCM for zachary karate club network.
+As an example, we run community detection on zachary karate club network.
 
 ```
     import numpy as np
     import networkx as nx
-    from NEMtropy import UndirectedGraph
+    from surprisememore import UndirectedGraph
 
+    from surprisememore import UndirectedGraph
+    import networkx as nx
+    
     G = nx.karate_club_graph()
     adj_kar = nx.to_numpy_array(G)
     graph = UndirectedGraph(adj_kar)
+    
+    graph.run_discrete_community_detection(weighted=False,
+                                           num_sim=2)
+```
+The algorithm will find the best partition by optimizing surprise score
+function. At the end of the optimization process, the optimal partition is
+saved as an attribute of the graph class.
 
-    graph.solve_tool(model="cm_exp",
-                 method="newton",
-                 initial_guess="random")
 ```
+    # optimal partition
+    graph.solution
+    
+    # Surprise of the optimal partition
+    graph.surprise
+    
+    # Log surprise
+    graph.log_surprise
+```
+
+Similarly, we can run the algorithm detecting bimodular structure. In the case
+of zachary karate club, the code snippet is the following.
 
-Given the UBCM model, we can generate ten random copies of zachary's karate club.
+#%% md
 
 ```
-    graph.ensemble_sampler(10, cpu_n=2, output_dir="sample/")
+    from surprisememore import UndirectedGraph
+    import networkx as nx
+    
+    G = nx.karate_club_graph()
+    adj_kar = nx.to_numpy_array(G)
+    graph = UndirectedGraph(adjacency=adj_kar)
 ```
 
-These copies are saved as an edgelist, each edgelist can be converted to an
-adjacency matrix by running the NEMtropy build graph function.
+Here we initialized our SupriseMeMore **UndirectedGraph** object with the adjacency
+matrix. The available options are adjacency matrix or edgelist.
 
-```
-    from NEMtropy.network_functions import build_graph_from_edgelist
+* If you use adjacency matrix, then you have to pass the matrix as a **numpy.ndarray**;
+
+* If you use edgelist, then the edgelist has to be passed as a **list of tuple**:
+    * [(u, v), (u, t), ...] for binary networks;
+    * [(u, v, w1), (u, t, w2), ...] for weighted networks;
+
+For more details about edgelist format you can see [link](https://networkx.org/documentation/stable/reference/classes/generated/networkx.DiGraph.add_weighted_edges_from.html?highlight=add_weighted_edges_from#networkx.DiGraph.add_weighted_edges_from).
 
-    edgelist_ens = np.loadtxt("sample/0.txt")
-    ens_adj = build_graph_from_edgelist(edgelist = edgelist_ens,
-                                    is_directed = False,
-                                    is_sparse = False,
-                                    is_weighted = False)
 ```
+    graph.run_discrete_cp_detection(weighted=False, num_sim=2)
+```
+
+In the previous example I passed two optional arguments to the function: *weighted*
+and *num_sim*. The argument *weighted* specify which version of surprise you want 
+to use: binary or weighted. If the network is binary, you don't need to pass 
+the argument "weighted", the class detects by itself that the graph is binary 
+and use the proper method for community/bimodular detection. Instead, if the 
+network has weights, the default method is the weighted one. To run binary 
+community/bimodular detection you must specify "weighted"=False.
 
-These collection of random adjacency matrices can be used as a null model:
-it is enough to compute the expected value of the selected network feature 
-on the ensemble of matrices and to compare it with its original value.
+The arguments *num_sim* specifies the number of time we run over all the edges 
+of the network during the optimization problem. You can find more detail about the
+algorithm in [1](https://arxiv.org/abs/2106.05055), [2](https://www.nature.com/articles/srep19250).
 
-To learn more, please read the two ipython notebooks in the examples directory:
-one is a study case on a [directed graph](https://github.com/nicoloval/NEMtropy/blob/master/examples/Directed%20Graphs.ipynb), 
-while the other is on an [undirected graph](https://github.com/nicoloval/NEMtropy/blob/master/examples/Undirected%20Graphs.ipynb).
+All the implemented algorithms are heuristic, we suggest running them more 
+than once and pick the best solution (the one with higher log_surprise).
 
-You can find complete documentation about NEMtropy library in [docs](https://nemtropy.readthedocs.io/en/master/index.html).
+To learn more, please read the two ipython notebooks in the examples' directory:
+one is a study case on a [community detection](https://github.com/nicoloval/NEMtropy/blob/master/examples/Community%20Detection.ipynb), 
+while the other is on an [bimodular detection](https://github.com/nicoloval/NEMtropy/blob/master/examples/Mesoscale%20Structure%20Detection.ipynb).
 
 Development
 -----------
 Please work on a feature branch and create a pull request to the development 
-branch. If necessary to merge manually do so without fast forward:
+branch. If necessary to merge manually do so without fast-forward:
 
 ```
     $ git merge --no-ff myfeature
@@ -171,29 +163,14 @@ To build a development environment run:
     $ pip install -e '.[dev]'
 ```
 
-Testing
--------
-If you want to test the package integrity, you can run the following 
-bash command from the tests directory:
-
-```
-    $ bash run_all.sh
-```
-
-__P.S.__ _at the moment there may be some problems with the DECM solver functions_
-
 Credits
 -------
 
-_Authors_:
-
-[Nicolò Vallarano](http://www.imtlucca.it/en/nicolo.vallarano/) (a.k.a. [nicoloval](https://github.com/nicoloval))
+_Author_:
 
 [Emiliano Marchese](https://www.imtlucca.it/en/emiliano.marchese/) (a.k.a. [EmilianoMarchese](https://github.com/EmilianoMarchese))
 
-[Matteo Bruno](https://www.imtlucca.it/en/matteo.bruno/) (BiCM) (a.k.a. [mat701](https://github.com/mat701))
 
 _Acknowledgements:_
 
-The module was developed under the supervision of [Tiziano Squartini](http://www.imtlucca.it/en/tiziano.squartini/), [Fabio Saracco](http://www.imtlucca.it/en/fabio.saracco/), [Mario Zanon](http://www.imtlucca.it/it/mario.zanon/), and [Giulio Cimini](https://www.fisica.uniroma2.it/elenco-telefonico/ciminigi/).
-It was developed at [IMT school of advanced studies Lucca](https://www.imtlucca.it/), and financed by the research project Optimized Reconstruction of Complex networkS - ORCS.
+The module was developed under the supervision of [Tiziano Squartini](http://www.imtlucca.it/en/tiziano.squartini/).
diff --git a/examples/Community Detection.ipynb b/examples/Community Detection.ipynb
new file mode 100644
index 00000000..c2cb6eb4
--- /dev/null
+++ b/examples/Community Detection.ipynb
@@ -0,0 +1,562 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Community Detection"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "The scope of this notebook is to provide some examples of how SurpriseMeMore can\n",
+    "be used to detect communities in binary and weighted network. We use the undirected\n",
+    "*les miserables* character network in the next examples, but all the applications\n",
+    " that we are going to see can be generalized to the case of directed networks."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Important\n",
+    "\n",
+    "All the methods we present are heuristic, we suggest then to run it more than once\n",
+    "to increase the chance to find the best partitioning."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "from surprisememore import Undirected_Graph_Class as ug\n",
+    "import numpy as np\n",
+    "import networkx as nx\n",
+    "import pandas as pd\n",
+    "\n",
+    "aux_path = 'out.moreno_lesmis_lesmis'\n",
+    "prova = np.loadtxt(aux_path,comments='%')\n",
+    "edgelist = pd.DataFrame(prova,columns=['source','target','weight'])\n",
+    "network = nx.from_pandas_edgelist(edgelist,source='source',target='target',edge_attr=True,create_using=nx.Graph)\n",
+    "adjacency_matrix = nx.to_numpy_array(network)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   },
+   "execution_count": 8,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Here we initialize our SupriseMeMore **UndirectedGraph** object using the adjacency\n",
+    "matrix. The available options are adjacency matrix or edgelist.\n",
+    "\n",
+    "* If you use adjacency matrix, then you have to pass the matrix as a **numpy.ndarray**;\n",
+    "\n",
+    "* If you use edgelist, then the edgelist has to be passed as a **list of tuple**:\n",
+    "    * [(u, v), (u, t), ...] for binary networks;\n",
+    "    * [(u, v, w1), (u, t, w2), ...] for weighted networks;\n",
+    "\n",
+    "For more details about edgelist format you can see [link](https://networkx.org/documentation/stable/reference/classes/generated/networkx.DiGraph.add_weighted_edges_from.html?highlight=add_weighted_edges_from#networkx.DiGraph.add_weighted_edges_from)."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Basic Usage"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "outputs": [],
+   "source": [
+    "graph = ug.UndirectedGraph(adjacency=adjacency_matrix)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Now that our graph istance is initalized we can run discrete community\n",
+    "detection by tapping the following command."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:00<00:00, 14570.93it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 30605.68it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[17  4  4  4 14 15 16  4  4 11  0 19 20 18 21 22  2  1  1  0  0  5  3  8\n",
+      " 12 10  3  3  3  3  3 13  5  0  0  0  0  0  0  1  1  1  1 28  2  2  2  2\n",
+      "  2  2  2  8  1  1  1  9 27  1 26  5  9  0  0  6  6  0  0  0  0  0  0  7\n",
+      "  7 24 25 23 29]\n",
+      "349.86182922756353\n",
+      "0.0\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "graph.run_discrete_community_detection()\n",
+    "\n",
+    "# The optimal partitioning is given by\n",
+    "print(graph.solution)\n",
+    "# This is a numpy array where the membership of each node is stored. Node with\n",
+    "# the same membership are in the same cluster.\n",
+    "\n",
+    "# The relative logsurprise is\n",
+    "print(graph.log_surprise)\n",
+    "# and the associated p-value (surprise) is\n",
+    "print(graph.surprise)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Important\n",
+    "\n",
+    "Everytime that you run the algorithm the values *solution*, *log_surprise* and\n",
+    "*surprise* are overwritten."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Community Detection Arguments"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "SurpriseMeMore community detection methods allow the user to pass arguments\n",
+    "specifying the used method to use and other useful options. In what follows we\n",
+    "will brifly discuss these options."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### method\n",
+    "\n",
+    "First we can specify the clustering method: **aglomerative** or **fixed-clusters**.\n",
+    "\n",
+    "**Aglomerative** method starts with all the nodes in different partitions and then start\n",
+    "to merge these until no further improvements are possible.\n",
+    "\n",
+    "**Fixed-clusters** solves the partioning problem with a fixed number of clusters\n",
+    "(defined by the user). If *fixed-clusters* is passed, the number of cluster argument\n",
+    " (*num_clusters*) must be specified.\n",
+    "\n",
+    "The default value is *aglomerative*."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### initial_guess\n",
+    "\n",
+    "The user can pass its own initial guess to the algorithm (pay attention that\n",
+    "has to be a proper initial guess) or use one of the implemented one.\n",
+    "\n",
+    "* *random*: membership is assigned to nodes randomly. If the *method* is aglomerative\n",
+    "it doesn't affect the initial guess;\n",
+    "\n",
+    "* *common-neigh-weak*: nodes are in the same initial partition based on the fraction\n",
+    "of common-neighbours. The condtion of aggregation is more relaxed.\n",
+    "\n",
+    "* *common-neigh-strong*: nodes are in the same initial partition based on the fraction\n",
+    "of common-neighbours. The condtion of aggregation is more strictly.\n",
+    "\n",
+    "The default value is *random*."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### weighted\n",
+    "\n",
+    "This argument has to be used when we initialize a UndirectedGraph (DirectedGraph)\n",
+    "instance that is weighted. In that case, if we want to run binary community detection\n",
+    "we must specify *weighted*=False.\n",
+    "\n",
+    "```\n",
+    "    graph.run_discrete_community_detection('weighted'=False)\n",
+    "```\n",
+    "\n",
+    "The above snippet of code run binary community detection on les miserables graph.\n",
+    "\n",
+    "The *weighted* argument is just for discrete community detection methods,\n",
+    "in the case of enhanced or continuous community detection there is no binary\n",
+    "version of the algorithm.\n",
+    "\n",
+    "The default value is *None*, the algorithm will choose the proper method for\n",
+    "the network:\n",
+    "\n",
+    "* weighted network --> weighted surprise;\n",
+    "\n",
+    "* binary netowrk --> binary surprise;"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### num_sim\n",
+    "\n",
+    "Number of times the algorithm will run over all the links trying to improve the\n",
+    "partioning. If no improvements are detected for 10 times in a row then the algorithm\n",
+    "stops.\n",
+    "\n",
+    "The default value is 2."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "##  Discrete Community Detection"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Binary"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:00<00:00, 22463.96it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 44629.60it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 48051.65it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 8280.96it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 7845.83it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 8074.04it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Aglomerative solution [16  5  5  5 13 18 19 20 21  5  3 17 23  8  3  8  2  1  1  4  1  6  3 12\n",
+      " 15 14  3  3  3  3  3  7  6  0  4  4  0  0  0  1  1  1  1  7  2  2  2  2\n",
+      "  2  2  2 12  1  1 22 10 25  1  4  6 10  0  0  9  9  0  0  0  0  0  0 11\n",
+      " 11  0 24  4 26]\n",
+      "Fixed clusters solution [0 1 1 2 2 0 0 0 3 2 0 3 3 2 3 1 1 1 0 0 0 1 0 1 0 1 0 0 1 3 2 2 3 0 0 1 0\n",
+      " 1 3 1 0 3 1 2 1 1 2 0 0 0 1 1 0 0 2 0 2 2 0 3 2 1 2 3 1 1 0 1 1 0 0 3 2 1\n",
+      " 2 0 2]\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_discrete_community_detection(weighted=False,\n",
+    "                                       initial_guess=\"random\",\n",
+    "                                       num_sim=3)\n",
+    "print(\"Aglomerative solution\", graph.solution)\n",
+    "\n",
+    "# for the fixed clusters one\n",
+    "graph.run_discrete_community_detection(method=\"fixed-clusters\",\n",
+    "                                       num_clusters=4,\n",
+    "                                       weighted=False,\n",
+    "                                       initial_guess=\"random\",\n",
+    "                                       num_sim=3)\n",
+    "print(\"Fixed clusters solution\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Weighted"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "outputs": [],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_discrete_community_detection(weighted=True,\n",
+    "                                       initial_guess=\"random\",\n",
+    "                                       num_sim=3)\n",
+    "print(\"Aglomerative solution\", graph.solution)\n",
+    "# for the fixed clusters one\n",
+    "graph.run_discrete_community_detection(method=\"fixed-clusters\",\n",
+    "                                       num_clusters=3,\n",
+    "                                       weighted=True,\n",
+    "                                       initial_guess=\"random\",\n",
+    "                                       num_sim=3)\n",
+    "print(\"Fixed clusters solution\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Enhanced Community Detection"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "For more details about the differences between enhanced methods and discrete one\n",
+    " read the relative paper (you can find the link in the readme)."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "outputs": [],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_enhanced_community_detection(initial_guess=\"common-neigh-weak\",\n",
+    "                                       num_sim=4)\n",
+    "print(\"Aglomerative solution\", graph.solution)\n",
+    "# for the fixed clusters one\n",
+    "graph.run_enhanced_community_detection(method=\"fixed-clusters\",\n",
+    "                                       num_clusters=4,\n",
+    "                                       initial_guess=\"common-neigh-weak\",\n",
+    "                                       num_sim=4)\n",
+    "print(\"Fixed clusters solution\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Continuous Community Detection"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Continuos community detection has to be used when your graph has continuos weights.\n",
+    "It requires integration then it is way slower."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "outputs": [],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_continuous_community_detection(initial_guess=\"common-neigh-weak\",\n",
+    "                                         num_sim=3)\n",
+    "print(\"Aglomerative solution\", graph.solution)\n",
+    "# for the fixed clusters one\n",
+    "graph.run_enhanced_community_detection(method=\"fixed-clusters\",\n",
+    "                                       num_clusters=5,\n",
+    "                                       initial_guess=\"common-neigh-weak\",\n",
+    "                                       num_sim=4)\n",
+    "print(\"Fixed clusters solution\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "outputs": [],
+   "source": [],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "outputs": [],
+   "source": [],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/examples/Mesoscale Structure Detection.ipynb b/examples/Mesoscale Structure Detection.ipynb
new file mode 100644
index 00000000..58c44074
--- /dev/null
+++ b/examples/Mesoscale Structure Detection.ipynb
@@ -0,0 +1,496 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Mesoscale Structures Detection"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "The scope of this notebook is to provide some examples of how SurpriseMeMore can\n",
+    "be used to detect mesoscale structure (e.g. core periphery) in binary\n",
+    " and weighted network. We use the undirected *les miserables* character network\n",
+    " in the next examples, but all the applications\n",
+    " that we are going to see can be generalized to the case of directed networks."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Important\n",
+    "\n",
+    "All the methods we present are heuristic, we suggest then to run it more than once\n",
+    "to increase the chance to find the best partitioning. The methods usually find the\n",
+    "best solution without running it too many times."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "from surprisememore import Undirected_Graph_Class as ug\n",
+    "import numpy as np\n",
+    "import networkx as nx\n",
+    "import pandas as pd\n",
+    "\n",
+    "aux_path = 'out.moreno_lesmis_lesmis'\n",
+    "prova = np.loadtxt(aux_path,comments='%')\n",
+    "edgelist = pd.DataFrame(prova,columns=['source','target','weight'])\n",
+    "network = nx.from_pandas_edgelist(edgelist,source='source',target='target',edge_attr=True,create_using=nx.Graph)\n",
+    "adjacency_matrix = nx.to_numpy_array(network)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   },
+   "execution_count": 4,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Here we initialize our SupriseMeMore **UndirectedGraph** object using the adjacency\n",
+    "matrix. The available options are adjacency matrix or edgelist.\n",
+    "\n",
+    "* If you use adjacency matrix, then you have to pass the matrix as a **numpy.ndarray**;\n",
+    "\n",
+    "* If you use edgelist, then the edgelist has to be passed as a **list of tuple**:\n",
+    "    * [(u, v), (u, t), ...] for binary networks;\n",
+    "    * [(u, v, w1), (u, t, w2), ...] for weighted networks;\n",
+    "\n",
+    "For more details about edgelist format you can see [link](https://networkx.org/documentation/stable/reference/classes/generated/networkx.DiGraph.add_weighted_edges_from.html?highlight=add_weighted_edges_from#networkx.DiGraph.add_weighted_edges_from)."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Basic Usage"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "outputs": [],
+   "source": [
+    "graph = ug.UndirectedGraph(adjacency=adjacency_matrix)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Now that our graph istance is initalized we can run discrete detection of mesoscale\n",
+    "structures by tapping the following command."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:00<00:00, 2400.76it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 2619.18it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0\n",
+      " 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0 1 1 1\n",
+      " 1 1 1]\n",
+      "217.99342901379407\n",
+      "1.0152452966619884e-218\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "graph.run_discrete_cp_detection()\n",
+    "\n",
+    "# The optimal partitioning is given by\n",
+    "print(graph.solution)\n",
+    "# This is a numpy array where the membership of each node is stored. Node with\n",
+    "# the same membership are in the same cluster.\n",
+    "\n",
+    "# The relative logsurprise is\n",
+    "print(graph.log_surprise)\n",
+    "# and the associated p-value (surprise) is\n",
+    "print(graph.surprise)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Important\n",
+    "\n",
+    "Everytime that you run the algorithm the values *solution*, *log_surprise* and\n",
+    "*surprise* are overwritten."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Mesoscale Detection Arguments"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "SurpriseMeMore mesoscale detection methods allow the user to pass arguments\n",
+    "specifying some aspects of the optimization process. In what follows we\n",
+    "will brifly discuss these options."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### initial_guess\n",
+    "\n",
+    "The user can pass its own initial guess to the algorithm (pay attention that\n",
+    "has to be a proper initial guess) or use one of the implemented one.\n",
+    "\n",
+    "* *random*: membership is assigned to nodes randomly. If the *method* is aglomerative\n",
+    "it doesn't affect the initial guess;\n",
+    "\n",
+    "* *ranked*: The 5% nodes with higher degree/strength are in the core. We suggest\n",
+    " using this option if you are looking for *core-periphery* or *bow-tie structures*.\n",
+    "\n",
+    "* *eigenvector*: nodes 5% nodes with higher eigenvector centrality are in the core.\n",
+    " We suggest using this option if you are looking for *core-periphery* or\n",
+    " *bow-tie structures*.\n",
+    "\n",
+    "The default value is *ranked*."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### weighted\n",
+    "\n",
+    "This argument has to be used when we initialize a UndirectedGraph (DirectedGraph)\n",
+    "instance that is weighted. In that case, if we want to run binary mesoscale detection\n",
+    "we must specify *weighted*=False.\n",
+    "\n",
+    "```\n",
+    "    graph.run_discrete_cp_detection('weighted'=False)\n",
+    "```\n",
+    "\n",
+    "The above snippet of code run binary mesoscale detection on les miserables graph.\n",
+    "\n",
+    "The *weighted* argument is just for discrete mesoscale detection methods,\n",
+    "in the case of enhanced mesoscale detection there is no binary version of\n",
+    "the algorithm.\n",
+    "\n",
+    "The default value is *None*, the algorithm will choose the proper method for\n",
+    "the network:\n",
+    "\n",
+    "* weighted network --> weighted surprise;\n",
+    "\n",
+    "* binary netowrk --> binary surprise;"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### num_sim\n",
+    "\n",
+    "Number of times the algorithm will run over all the links trying to improve the\n",
+    "partioning. If no improvements are detected for 10 times in a row then the algorithm\n",
+    "stops.\n",
+    "\n",
+    "The default value is 2."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "##  Discrete Mesoscale Detection\n"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Binary"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:01<00:00, 159.62it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 5079.40it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 5146.19it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The solution is [1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 1 1 0 1 1 0\n",
+      " 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 1 0 0 1 1 0 0 0 0 0 0 1 1 0\n",
+      " 1 1 1]\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_discrete_cp_detection(weighted=False,\n",
+    "                                initial_guess=\"random\",\n",
+    "                                num_sim=3)\n",
+    "print(\"The solution is\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Weighted"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:00<00:00, 3496.26it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 3920.89it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 2977.28it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The solution is [1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0\n",
+      " 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 1 0 0 1 1 0 0 0 0 0 0 1 1 1\n",
+      " 1 1 1]\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_discrete_cp_detection(weighted=True,\n",
+    "                                initial_guess=\"random\",\n",
+    "                                num_sim=3)\n",
+    "print(\"The solution is\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Enhanced Mesoscale Detection"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "For more details about the differences between enhanced methods and discrete ones\n",
+    " read the relative paper (you can find the link in the readme)."
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 254/254 [00:01<00:00, 252.60it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 2662.42it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 2646.33it/s]\n",
+      "100%|██████████| 254/254 [00:00<00:00, 2712.63it/s]\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The solution is [1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0\n",
+      " 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0 1 1 1\n",
+      " 1 1 1]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# An example of how we can run it\n",
+    "graph.run_enhanced_cp_detection(initial_guess=\"eigenvector\",\n",
+    "                                num_sim=4)\n",
+    "print(\"The solution is\", graph.solution)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "outputs": [],
+   "source": [],
+   "metadata": {
+    "collapsed": false,
+    "pycharm": {
+     "name": "#%%\n"
+    }
+   }
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/examples/out.moreno_lesmis_lesmis b/examples/out.moreno_lesmis_lesmis
new file mode 100644
index 00000000..0675592c
--- /dev/null
+++ b/examples/out.moreno_lesmis_lesmis
@@ -0,0 +1,256 @@
+% sym posweighted
+% 254 77 77
+1 2 1
+2 3 8
+2 4 10
+2 5 1
+2 6 1
+2 7 1
+2 8 1
+2 9 2
+2 10 1
+2 11 5
+3 4 6
+3 11 3
+4 11 3
+11 12 1
+11 13 1
+11 14 1
+11 15 1
+11 16 1
+11 24 9
+11 25 7
+11 26 12
+11 27 31
+11 28 17
+11 29 8
+11 30 2
+11 32 3
+11 33 1
+11 34 2
+11 35 3
+11 36 3
+11 37 2
+11 38 2
+11 39 2
+11 44 3
+11 45 1
+11 49 1
+11 50 2
+11 52 2
+11 56 19
+11 59 4
+11 65 1
+11 69 1
+11 70 1
+11 71 1
+11 72 1
+11 73 1
+13 24 2
+17 18 4
+17 19 4
+17 20 4
+17 21 3
+17 22 3
+17 23 3
+17 24 3
+18 19 4
+18 20 4
+18 21 3
+18 22 3
+18 23 3
+18 24 3
+18 27 1
+18 56 1
+19 20 4
+19 21 3
+19 22 3
+19 23 3
+19 24 3
+20 21 4
+20 22 3
+20 23 3
+20 24 3
+21 22 5
+21 23 4
+21 24 4
+22 23 4
+22 24 4
+23 24 4
+24 25 2
+24 26 1
+24 28 5
+24 30 1
+24 31 1
+24 32 2
+25 26 13
+25 27 4
+25 28 1
+25 42 2
+25 43 1
+25 51 1
+25 69 1
+25 70 1
+25 71 1
+26 27 1
+26 28 5
+26 40 1
+26 41 1
+26 42 3
+26 43 2
+26 49 1
+26 56 2
+26 69 5
+26 70 6
+26 71 4
+26 72 1
+26 76 3
+27 28 1
+27 44 1
+27 50 3
+27 52 2
+27 55 1
+27 56 21
+27 73 2
+28 29 1
+28 30 1
+28 32 1
+28 34 1
+28 44 1
+28 49 1
+28 59 6
+28 69 1
+28 70 2
+28 71 1
+28 72 1
+28 73 1
+29 45 3
+29 46 2
+30 35 2
+30 36 2
+30 37 1
+30 38 1
+30 39 1
+31 32 2
+35 36 3
+35 37 2
+35 38 2
+35 39 2
+36 37 2
+36 38 2
+36 39 2
+37 38 2
+37 39 2
+38 39 2
+40 53 1
+40 56 1
+42 43 2
+42 56 5
+42 58 1
+42 63 1
+42 69 1
+42 70 1
+42 71 1
+42 72 1
+42 76 1
+47 48 1
+47 49 2
+49 56 4
+49 58 1
+49 59 7
+49 60 6
+49 61 1
+49 62 2
+49 63 7
+49 64 5
+49 65 5
+49 66 3
+49 67 1
+49 69 1
+49 70 1
+49 72 1
+49 74 2
+49 75 2
+49 76 1
+49 77 1
+50 51 1
+50 52 9
+50 55 1
+50 56 12
+50 57 1
+52 53 1
+52 54 1
+52 55 2
+52 56 6
+55 56 1
+56 57 1
+56 58 1
+56 59 7
+56 60 5
+56 62 1
+56 63 9
+56 64 1
+56 65 5
+56 66 2
+58 59 1
+58 60 2
+58 62 1
+58 63 2
+58 64 2
+58 65 1
+58 66 1
+58 68 3
+59 60 15
+59 61 4
+59 62 6
+59 63 17
+59 64 4
+59 65 10
+59 66 5
+59 67 3
+59 71 1
+59 77 1
+60 61 2
+60 62 5
+60 63 13
+60 64 5
+60 65 9
+60 66 5
+60 67 1
+61 62 2
+61 63 3
+61 64 2
+61 65 2
+61 66 2
+61 67 1
+62 63 6
+62 64 3
+62 65 6
+62 66 5
+62 67 1
+63 64 6
+63 65 12
+63 66 5
+63 67 2
+63 77 1
+64 65 4
+64 66 5
+64 67 1
+64 77 1
+65 66 7
+65 67 3
+65 77 1
+66 67 2
+66 77 1
+67 77 1
+69 70 6
+69 71 4
+69 72 2
+69 76 3
+70 71 4
+70 72 2
+70 76 3
+71 72 2
+71 76 1
+72 76 1
+74 75 3
diff --git a/requirements.txt b/requirements.txt
index d048a2fb..6827f259 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,4 +3,8 @@ numba==0.52.0
 scipy==1.6.0
 networkx==2.5
 surprisemes==0.1.0
-setuptools==53.0.0
\ No newline at end of file
+setuptools==53.0.0
+tqdm==4.60.0
+pytest==6.2.2
+mpmath>=1.2.1
+sympy~=1.8
\ No newline at end of file
diff --git a/setup.py b/setup.py
index fada7a1e..6533d42a 100644
--- a/setup.py
+++ b/setup.py
@@ -4,19 +4,20 @@
     long_description = fh.read()
 
 setup(
-    name="surprisemes",
+    name="surprisememore",
     author="Emiliano Marchese",
     author_email='emilianomarcheserc@gmail.com',
-    packages=["surprisemes"],
+    packages=["surprisememore"],
     package_dir={'': 'src'},
     version="0.1.0",
-    description="Surprisemes is a python module providing methods, based on"
-                " the surprise framework, to detect mesoscale structures in"
-                " graphs and multigraphs.",
+    description="SurpriseMeMore is a python module providing methods, based on"
+                " the surprise framework, to detect mesoscale structures"
+                "(e.g. communities, core-periphery, bow-tie) in graphs and "
+                "multigraphs.",
     long_description=long_description,
     long_description_content_type="text/markdown",
     license="MIT",
-    url="",
+    url="https://github.com/EmilianoMarchese/SurpriseMeMore",
     download_url="",
     keywords=['community detection', 'core-periphery detection',
               'graphs', "multi-graphs", "weighted graphs", "surprise"],
@@ -31,17 +32,17 @@
         "numpy>=1.17",
         "networkx>=2.4",
         "scipy>=1.4",
-        "numba>=0.47"
+        "numba>=0.47",
+        "mpmath>=1.2",
+        "tqdm>=4.5"
                       ],
     extras_require={
-        "dev": ["pytest==6.0.1",
-                "coverage==5.2.1",
-                "pytest-cov==2.10.1",
-                "flake8==3.8.3",
-                "wheel==0.35.1",
-                "check-manifest==0.44",
-                "setuptools==47.1.0",
-                "twine==3.2.0",
-                "tox==3.20.1"],
+        "dev": ["pytest>=6.0.1",
+                "flake8>=3.8.3",
+                "wheel>=0.36.2",
+                "check-manifest>=0.44",
+                "setuptools>=47.1.0",
+                "twine>=3.2.0",
+                "tox>=3.20.1"],
         },
 )
diff --git a/src/surprisemes/Directed_Graph_Class.py b/src/surprisememore/Directed_Graph_Class.py
similarity index 60%
rename from src/surprisemes/Directed_Graph_Class.py
rename to src/surprisememore/Directed_Graph_Class.py
index 8139a5db..07e4498c 100644
--- a/src/surprisemes/Directed_Graph_Class.py
+++ b/src/surprisememore/Directed_Graph_Class.py
@@ -140,15 +140,40 @@ def clean_edges(self):
         self.edgelist = None
         self.is_initialized = False
 
-    def run_cp_detection(self,
-                         initial_guess="random",
-                         weighted=None,
-                         num_sim=2,
-                         sorting_method="default",
-                         print_output=False):
+    def run_enhanced_cp_detection(self,
+                                  initial_guess="random",
+                                  num_sim=2,
+                                  sorting_method="default",
+                                  print_output=False):
 
         self._initialize_problem_cp(
             initial_guess=initial_guess,
+            enhanced=True,
+            weighted=True,
+            sorting_method=sorting_method)
+
+        sol = solver.solver_cp(
+            adjacency_matrix=self.aux_adj,
+            cluster_assignment=self.init_guess,
+            num_sim=num_sim,
+            sort_edges=self.sorting_function,
+            calculate_surprise=self.surprise_function,
+            correct_partition_labeling=self.partition_labeler,
+            flipping_function=self.flipping_function,
+            print_output=print_output)
+
+        self._set_solved_problem(sol)
+
+    def run_discrete_cp_detection(self,
+                                  initial_guess="random",
+                                  weighted=None,
+                                  num_sim=2,
+                                  sorting_method="default",
+                                  print_output=False):
+
+        self._initialize_problem_cp(
+            initial_guess=initial_guess,
+            enhanced=False,
             weighted=weighted,
             sorting_method=sorting_method)
 
@@ -166,6 +191,7 @@ def run_cp_detection(self,
 
     def _initialize_problem_cp(self,
                                initial_guess,
+                               enhanced,
                                weighted,
                                sorting_method):
 
@@ -178,10 +204,20 @@ def _initialize_problem_cp(self,
                 self.aux_adj = self.adjacency
                 self.method = "binary"
         elif weighted:
-            try:
-                self.aux_adj = self.adjacency_weighted
+            if enhanced:
+                self.method = "enhanced"
+            else:
                 self.method = "weighted"
-            except Exception:
+
+            if hasattr(self, "adjacency_weighted"):
+                self.aux_adj = self.adjacency_weighted
+                cond2 = (self.aux_adj.astype(np.int64).sum() !=
+                         self.aux_adj.sum())
+                if cond2:
+                    raise ValueError("The selected method works for discrete "
+                                     "weights, but the initialised graph has "
+                                     "continuous weights.")
+            else:
                 raise TypeError(
                     "You choose weighted core peryphery detection but the"
                     " graph you initialised is binary.")
@@ -215,6 +251,10 @@ def _initialize_problem_cp(self,
                 x,
                 y,
                 True),
+            "enhanced": lambda x, y: cp.calculate_surprise_logsum_cp_enhanced(
+                x,
+                y,
+                True),
         }
 
         try:
@@ -230,6 +270,11 @@ def _set_initial_guess_cp(self, initial_guess):
         # TODO: Sistemare parte pesata
         if isinstance(initial_guess, str):
             if initial_guess == "random":
+                self.init_guess = np.ones(self.n_nodes, dtype=np.int32)
+                aux_n = int(np.ceil((5 * self.n_nodes) / 100))
+                self.init_guess[:aux_n] = 0
+                np.random.shuffle(self.init_guess[:aux_n])
+            elif initial_guess == "ranked":
                 self.init_guess = np.ones(self.n_nodes, dtype=np.int32)
                 aux_n = int(np.ceil((5 * self.n_nodes) / 100))
                 if self.is_weighted:
@@ -239,7 +284,7 @@ def _set_initial_guess_cp(self, initial_guess):
                     self.init_guess[
                         self.degree_sequence_out.argsort()[-aux_n:]] = 0
             elif initial_guess == "eigenvector":
-                self.initial_guess = ax.eigenvector_init_guess(self.adjacency,
+                self.init_guess = ax.eigenvector_init_guess(self.adjacency,
                                                                False)
             else:
                 raise ValueError("Valid values of initial guess are 'random', "
@@ -254,7 +299,7 @@ def _set_initial_guess_cp(self, initial_guess):
         if np.unique(self.init_guess).shape[0] != 2:
             raise ValueError("The custom initial_guess passed is not valid."
                              " The initial guess for core-periphery detection"
-                             " must have nodes' membership that are 0 and 1."
+                             " must have nodes' membership that are 0 or 1."
                              " Pay attention that at least one node has to "
                              "belong to the core (0) or the periphery (1).")
 
@@ -263,21 +308,121 @@ def _set_initial_guess_cp(self, initial_guess):
                 "The length of the initial guess provided is different from"
                 " the network number of nodes.")
 
-    def run_comunity_detection(self,
-                               method="aglomerative",
-                               initial_guess=None,
-                               weighted=None,
-                               num_sim=None,
-                               num_clusters=2,
-                               prob_mix=0.1,
-                               sorting_method="default",
-                               print_output=False):
+    def run_continuous_community_detection(self,
+                                           method="aglomerative",
+                                           initial_guess="random",
+                                           approx=None,
+                                           num_sim=2,
+                                           num_clusters=None,
+                                           prob_mix=0.1,
+                                           sorting_method="default",
+                                           print_output=False
+                                           ):
+        self._initialize_problem_cd(
+            method=method,
+            num_clusters=num_clusters,
+            initial_guess=initial_guess,
+            enhanced=False,
+            weighted=True,
+            continuous=True,
+            sorting_method=sorting_method)
+
+        if method == "aglomerative":
+            sol = solver.solver_com_det_aglom(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                prob_mix=prob_mix,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=approx,
+                is_directed=True,
+                print_output=print_output)
+        elif method == "fixed-clusters":
+            sol = solver.solver_com_det_divis(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=approx,
+                is_directed=True,
+                print_output=print_output)
+        else:
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
+
+        self._set_solved_problem(sol)
+
+    def run_enhanced_community_detection(self,
+                                         method="aglomerative",
+                                         initial_guess="random",
+                                         num_sim=2,
+                                         num_clusters=None,
+                                         prob_mix=0.1,
+                                         sorting_method="default",
+                                         print_output=False
+                                         ):
 
         self._initialize_problem_cd(
             method=method,
             num_clusters=num_clusters,
             initial_guess=initial_guess,
+            enhanced=True,
+            weighted=True,
+            continuous=False,
+            sorting_method=sorting_method)
+
+        if method == "aglomerative":
+            sol = solver.solver_com_det_aglom(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                prob_mix=prob_mix,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=None,
+                is_directed=True,
+                print_output=print_output)
+        elif method == "fixed-clusters":
+            sol = solver.solver_com_det_divis(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=None,
+                is_directed=True,
+                print_output=print_output)
+        else:
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
+
+        self._set_solved_problem(sol)
+
+    def run_discrete_community_detection(self,
+                                         method="aglomerative",
+                                         initial_guess=None,
+                                         weighted=None,
+                                         num_sim=None,
+                                         num_clusters=2,
+                                         prob_mix=0.1,
+                                         sorting_method="default",
+                                         print_output=False):
+
+        self._initialize_problem_cd(
+            method=method,
+            num_clusters=num_clusters,
+            initial_guess=initial_guess,
+            enhanced=False,
             weighted=weighted,
+            continuous=False,
             sorting_method=sorting_method)
 
         if method == "aglomerative":
@@ -289,10 +434,11 @@ def run_comunity_detection(self,
                 calculate_surprise=self.surprise_function,
                 correct_partition_labeling=self.partition_labeler,
                 prob_mix=prob_mix,
-                flipping_function=self.flipping_function,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=None,
                 is_directed=True,
                 print_output=print_output)
-        elif method == "divisive":
+        elif method == "fixed-clusters":
             sol = solver.solver_com_det_divis(
                 adjacency_matrix=self.aux_adj,
                 cluster_assignment=self.init_guess,
@@ -300,11 +446,12 @@ def run_comunity_detection(self,
                 sort_edges=self.sorting_function,
                 calculate_surprise=self.surprise_function,
                 correct_partition_labeling=self.partition_labeler,
-                flipping_function=self.flipping_function,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=None,
                 is_directed=True,
                 print_output=print_output)
         else:
-            raise ValueError("Method can be 'aglomerative' or 'divisive'.")
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
 
         self._set_solved_problem(sol)
 
@@ -312,7 +459,9 @@ def _initialize_problem_cd(self,
                                method,
                                num_clusters,
                                initial_guess,
+                               enhanced,
                                weighted,
+                               continuous,
                                sorting_method):
 
         self._set_initial_guess_cd(method, num_clusters, initial_guess)
@@ -324,13 +473,27 @@ def _initialize_problem_cd(self,
                 self.aux_adj = self.adjacency
                 self.method = "binary"
         elif weighted:
-            try:
-                self.aux_adj = self.adjacency_weighted
+            if enhanced:
+                self.method = "enhanced"
+            elif continuous:
+                self.method = "continuous"
+            else:
                 self.method = "weighted"
-            except Exception:
+
+            if hasattr(self, "adjacency_weighted"):
+                self.aux_adj = self.adjacency_weighted
+                cond1 = (self.method == "enhanced" or
+                         self.method == "weighted")
+                cond2 = (self.aux_adj.astype(np.int64).sum() !=
+                         self.aux_adj.sum())
+                if cond1 and cond2:
+                    raise ValueError("The selected method works for discrete "
+                                     "weights, but the initialised graph has "
+                                     "continuous weights.")
+            else:
                 raise TypeError(
-                    "You choose weighted comunity detection but the graph"
-                    " you initialised is binary.")
+                    "You choose weighted community detection but the"
+                    " graph you initialised is binary.")
         else:
             self.aux_adj = self.adjacency
             self.method = "binary"
@@ -354,13 +517,14 @@ def _initialize_problem_cd(self,
         surp_fun = {
             "binary": cd.calculate_surprise_logsum_clust_bin_new,
             "weighted": cd.calculate_surprise_logsum_clust_weigh_new,
+            "enhanced": cd.calculate_surprise_logsum_clust_enhanced_new,
+            "continuous": cd.calculate_surprise_logsum_clust_weigh_continuos,
         }
 
         self.surprise_function = surp_fun[self.method]
 
         # self.flipping_function = lambda x: CD.flipping_function_comdet(x)
-        self.flipping_function = lambda x, y: cd.flipping_function_comdet_new(
-            x, y, True)
+        # self.flipping_function = cd.flipping_function_comdet_new
 
         self.partition_labeler = lambda x: cd.labeling_communities(x)
 
@@ -369,8 +533,8 @@ def _set_initial_guess_cd(self,
                               num_clusters,
                               initial_guess):
 
-        if num_clusters is None and method == "divisive":
-            raise ValueError("When 'divisive' is passed as clustering 'method'"
+        if num_clusters is None and method == "fixed-clusters":
+            raise ValueError("When 'fixed-clusters' is passed as clustering 'method'"
                              " the 'num_clusters' argument must be specified.")
 
         if isinstance(initial_guess, str):
@@ -378,15 +542,24 @@ def _set_initial_guess_cd(self,
                 if method == "aglomerative":
                     self.init_guess = np.array(
                         [k for k in np.arange(self.n_nodes, dtype=np.int32)])
-                elif method == "divisive":
+                elif method == "fixed-clusters":
                     self.init_guess = np.random.randint(
                         low=num_clusters,
                         size=self.n_nodes)
-            elif initial_guess == "common-neighbours":
+            elif (initial_guess == "common-neigh-weak") or \
+                     (initial_guess == "common-neighbours"):
+                if method == "aglomerative":
+                    self.init_guess = ax.common_neigh_init_guess_weak(
+                        self.adjacency)
+                elif method == "fixed-clusters":
+                    self.init_guess = ax.fixed_clusters_init_guess_cn(
+                        adjacency=self.adjacency,
+                        n_clust=num_clusters)
+            elif initial_guess == "common-neigh-strong":
                 if method == "aglomerative":
-                    self.init_guess = ax.common_neigh_init_guess(
+                    self.init_guess = ax.common_neigh_init_guess_strong(
                         self.adjacency)
-                elif method == "divisive":
+                elif method == "fixed-clusters":
                     self.init_guess = ax.fixed_clusters_init_guess_cn(
                         adjacency=self.adjacency,
                         n_clust=num_clusters)
@@ -395,7 +568,8 @@ def _set_initial_guess_cd(self,
                     "The 'initial_guess' selected is not a valid."
                     "Initial guess can be an array specifying nodes membership"
                     " or an initialisation method ['common-neighbours',"
-                    " random]. For more details see documentation.")
+                    " 'random', 'common-neigh-weak', 'common-neigh-strong']."
+                    " For more details see documentation.")
 
         elif isinstance(initial_guess, np.ndarray):
             self.init_guess = initial_guess
@@ -407,11 +581,11 @@ def _set_initial_guess_cd(self,
                 "The length of the initial guess provided is different"
                 " from the network number of nodes.")
 
-        if (method == "divisive" and
+        if (method == "fixed-clusters" and
                 np.unique(self.init_guess).shape[0] != num_clusters):
             raise ValueError("The number of clusters of a custom initial guess"
                              " must coincide with 'num_clusters' when the "
-                             " divisive method is applied.")
+                             " fixed-clusters method is applied.")
 
     def _set_solved_problem(self, sol):
         self.solution = sol[0]
diff --git a/src/surprisemes/Undirected_Graph_Class.py b/src/surprisememore/Undirected_Graph_Class.py
similarity index 58%
rename from src/surprisemes/Undirected_Graph_Class.py
rename to src/surprisememore/Undirected_Graph_Class.py
index d89f14a2..bf34c60f 100644
--- a/src/surprisemes/Undirected_Graph_Class.py
+++ b/src/surprisememore/Undirected_Graph_Class.py
@@ -122,15 +122,41 @@ def clean_edges(self):
         self.edgelist = None
         self.is_initialized = False
 
-    def run_cp_detection(self,
-                         initial_guess="random",
-                         weighted=None,
-                         num_sim=2,
-                         sorting_method="default",
-                         print_output=False):
+    def run_enhanced_cp_detection(self,
+                                  initial_guess="ranked",
+                                  num_sim=2,
+                                  sorting_method="default",
+                                  print_output=False):
 
         self._initialize_problem_cp(
             initial_guess=initial_guess,
+            enhanced=True,
+            weighted=True,
+            sorting_method=sorting_method)
+
+        sol = solver.solver_cp(
+            adjacency_matrix=self.aux_adj,
+            cluster_assignment=self.init_guess,
+            num_sim=num_sim,
+            sort_edges=self.sorting_function,
+            calculate_surprise=self.surprise_function,
+            correct_partition_labeling=self.partition_labeler,
+            flipping_function=self.flipping_function,
+            print_output=print_output)
+
+        self._set_solved_problem(sol)
+
+    def run_discrete_cp_detection(
+            self,
+            initial_guess="ranked",
+            weighted=None,
+            num_sim=2,
+            sorting_method="default",
+            print_output=False):
+
+        self._initialize_problem_cp(
+            initial_guess=initial_guess,
+            enhanced=False,
             weighted=weighted,
             sorting_method=sorting_method)
 
@@ -148,6 +174,7 @@ def run_cp_detection(self,
 
     def _initialize_problem_cp(self,
                                initial_guess,
+                               enhanced,
                                weighted,
                                sorting_method):
 
@@ -160,10 +187,20 @@ def _initialize_problem_cp(self,
                 self.aux_adj = self.adjacency
                 self.method = "binary"
         elif weighted:
-            try:
-                self.aux_adj = self.adjacency_weighted
+            if enhanced:
+                self.method = "enhanced"
+            else:
                 self.method = "weighted"
-            except Exception:
+
+            if hasattr(self, "adjacency_weighted"):
+                self.aux_adj = self.adjacency_weighted
+                cond2 = (self.aux_adj.astype(np.int64).sum() !=
+                         self.aux_adj.sum())
+                if cond2:
+                    raise ValueError("The selected method works for discrete "
+                                     "weights, but the initialised graph has "
+                                     "continuous weights.")
+            else:
                 raise TypeError(
                     "You choose weighted core peryphery detection but the"
                     " graph you initialised is binary.")
@@ -198,6 +235,10 @@ def _initialize_problem_cp(self,
                 x,
                 y,
                 False),
+            "enhanced": lambda x, y: cp.calculate_surprise_logsum_cp_enhanced(
+                x,
+                y,
+                False),
         }
 
         try:
@@ -214,7 +255,12 @@ def _set_initial_guess_cp(self, initial_guess):
         if isinstance(initial_guess, str):
             if initial_guess == "random":
                 self.init_guess = np.ones(self.n_nodes, dtype=np.int32)
-                aux_n = int(np.ceil((5*self.n_nodes)/100))
+                aux_n = int(np.ceil((5 * self.n_nodes) / 100))
+                self.init_guess[:aux_n] = 0
+                np.random.shuffle(self.init_guess[:aux_n])
+            elif initial_guess == "ranked":
+                self.init_guess = np.ones(self.n_nodes, dtype=np.int32)
+                aux_n = int(np.ceil((5 * self.n_nodes) / 100))
                 if self.is_weighted:
                     self.init_guess[
                         self.strength_sequence.argsort()[-aux_n:]] = 0
@@ -222,12 +268,13 @@ def _set_initial_guess_cp(self, initial_guess):
                     self.init_guess[
                         self.degree_sequence.argsort()[-aux_n:]] = 0
             elif initial_guess == "eigenvector":
-                self.initial_guess = ax.eigenvector_init_guess(self.adjacency,
-                                                               False)
+                self.init_guess = ax.eigenvector_init_guess(
+                    self.adjacency,
+                    False)
             else:
                 raise ValueError("Valid values of initial guess are 'random', "
-                                 "eigenvector or a custom initial guess ("
-                                 "np.ndarray or list).")
+                                 "'eigenvector', 'ranked, or a custom initial"
+                                 " guess (np.ndarray or list).")
 
         elif isinstance(initial_guess, np.ndarray):
             self.init_guess = initial_guess
@@ -237,7 +284,7 @@ def _set_initial_guess_cp(self, initial_guess):
         if np.unique(self.init_guess).shape[0] != 2:
             raise ValueError("The custom initial_guess passed is not valid."
                              " The initial guess for core-periphery detection"
-                             " must have nodes' membership that are 0 and 1."
+                             " must have nodes' membership that are 0 or 1."
                              " Pay attention that at least one node has to "
                              "belong to the core (0) or the periphery (1).")
 
@@ -246,21 +293,121 @@ def _set_initial_guess_cp(self, initial_guess):
                 "The length of the initial guess provided is different from"
                 " the network number of nodes.")
 
-    def run_comunity_detection(self,
-                               method="aglomerative",
-                               initial_guess="random",
-                               weighted=None,
-                               num_sim=2,
-                               num_clusters=None,
-                               prob_mix=0.1,
-                               sorting_method="default",
-                               print_output=False):
+    def run_continuous_community_detection(self,
+                                           method="aglomerative",
+                                           initial_guess="random",
+                                           approx=None,
+                                           num_sim=2,
+                                           num_clusters=None,
+                                           prob_mix=0.1,
+                                           sorting_method="default",
+                                           print_output=False
+                                           ):
+        self._initialize_problem_cd(
+            method=method,
+            num_clusters=num_clusters,
+            initial_guess=initial_guess,
+            enhanced=False,
+            weighted=True,
+            continuous=True,
+            sorting_method=sorting_method)
+
+        if method == "aglomerative":
+            sol = solver.solver_com_det_aglom(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                prob_mix=prob_mix,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=approx,
+                is_directed=False,
+                print_output=print_output)
+        elif method == "fixed-clusters":
+            sol = solver.solver_com_det_divis(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=approx,
+                is_directed=False,
+                print_output=print_output)
+        else:
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
+
+        self._set_solved_problem(sol)
+
+    def run_enhanced_community_detection(self,
+                                         method="aglomerative",
+                                         initial_guess="random",
+                                         num_sim=2,
+                                         num_clusters=None,
+                                         prob_mix=0.1,
+                                         sorting_method="default",
+                                         print_output=False
+                                         ):
 
         self._initialize_problem_cd(
             method=method,
             num_clusters=num_clusters,
             initial_guess=initial_guess,
+            enhanced=True,
+            weighted=True,
+            continuous=False,
+            sorting_method=sorting_method)
+
+        if method == "aglomerative":
+            sol = solver.solver_com_det_aglom(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                prob_mix=prob_mix,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=None,
+                is_directed=False,
+                print_output=print_output)
+        elif method == "fixed-clusters":
+            sol = solver.solver_com_det_divis(
+                adjacency_matrix=self.aux_adj,
+                cluster_assignment=self.init_guess,
+                num_sim=num_sim,
+                sort_edges=self.sorting_function,
+                calculate_surprise=self.surprise_function,
+                correct_partition_labeling=self.partition_labeler,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=None,
+                is_directed=False,
+                print_output=print_output)
+        else:
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
+
+        self._set_solved_problem(sol)
+
+    def run_discrete_community_detection(self,
+                                         method="aglomerative",
+                                         initial_guess="random",
+                                         weighted=None,
+                                         num_sim=2,
+                                         num_clusters=None,
+                                         prob_mix=0.1,
+                                         sorting_method="default",
+                                         print_output=False):
+
+        self._initialize_problem_cd(
+            method=method,
+            num_clusters=num_clusters,
+            initial_guess=initial_guess,
+            enhanced=False,
             weighted=weighted,
+            continuous=False,
             sorting_method=sorting_method)
 
         if method == "aglomerative":
@@ -272,10 +419,11 @@ def run_comunity_detection(self,
                 calculate_surprise=self.surprise_function,
                 correct_partition_labeling=self.partition_labeler,
                 prob_mix=prob_mix,
-                flipping_function=self.flipping_function,
+                flipping_function=cd.flipping_function_comdet_agl_new,
+                approx=None,
                 is_directed=False,
                 print_output=print_output)
-        elif method == "divisive":
+        elif method == "fixed-clusters":
             sol = solver.solver_com_det_divis(
                 adjacency_matrix=self.aux_adj,
                 cluster_assignment=self.init_guess,
@@ -283,11 +431,12 @@ def run_comunity_detection(self,
                 sort_edges=self.sorting_function,
                 calculate_surprise=self.surprise_function,
                 correct_partition_labeling=self.partition_labeler,
-                flipping_function=self.flipping_function,
+                flipping_function=cd.flipping_function_comdet_div_new,
+                approx=None,
                 is_directed=False,
                 print_output=print_output)
         else:
-            raise ValueError("Method can be 'aglomerative' or 'divisive'.")
+            raise ValueError("Method can be 'aglomerative' or 'fixed-clusters'.")
 
         self._set_solved_problem(sol)
 
@@ -295,7 +444,9 @@ def _initialize_problem_cd(self,
                                method,
                                num_clusters,
                                initial_guess,
+                               enhanced,
                                weighted,
+                               continuous,
                                sorting_method):
 
         self._set_initial_guess_cd(method, num_clusters, initial_guess)
@@ -307,12 +458,27 @@ def _initialize_problem_cd(self,
                 self.aux_adj = self.adjacency
                 self.method = "binary"
         elif weighted:
-            try:
-                self.aux_adj = self.adjacency_weighted
+            if enhanced:
+                self.method = "enhanced"
+            elif continuous:
+                self.method = "continuous"
+            else:
                 self.method = "weighted"
-            except Exception:
-                raise TypeError("You choose weighted comunity detection but "
-                                "the graph you initialised is binary.")
+
+            if hasattr(self, "adjacency_weighted"):
+                self.aux_adj = self.adjacency_weighted
+                cond1 = (self.method == "enhanced" or
+                         self.method == "weighted")
+                cond2 = (self.aux_adj.astype(np.int64).sum() !=
+                         self.aux_adj.sum())
+                if cond1 and cond2:
+                    raise ValueError("The selected method works for discrete "
+                                     "weights, but the initialised graph has "
+                                     "continuous weights.")
+            else:
+                raise TypeError(
+                    "You choose weighted core peryphery detection but the"
+                    " graph you initialised is binary.")
         else:
             self.aux_adj = self.adjacency
             self.method = "binary"
@@ -337,15 +503,14 @@ def _initialize_problem_cd(self,
         surp_fun = {
             "binary": cd.calculate_surprise_logsum_clust_bin_new,
             "weighted": cd.calculate_surprise_logsum_clust_weigh_new,
+            "enhanced": cd.calculate_surprise_logsum_clust_enhanced_new,
+            "continuous": cd.calculate_surprise_logsum_clust_weigh_continuos,
         }
 
         self.surprise_function = surp_fun[self.method]
 
         # self.flipping_function = lambda x: CD.flipping_function_comdet(x)
-        self.flipping_function = lambda x, y: cd.flipping_function_comdet_new(
-            x,
-            y,
-            False)
+        # self.flipping_function = cd.flipping_function_comdet_new
 
         self.partition_labeler = lambda x: cd.labeling_communities(x)
 
@@ -353,8 +518,9 @@ def _set_initial_guess_cd(self,
                               method,
                               num_clusters,
                               initial_guess):
-        if num_clusters is None and method == "divisive":
-            raise ValueError("When 'divisive' is passed as clustering 'method'"
+        if num_clusters is None and method == "fixed-clusters":
+            raise ValueError("When 'fixed-clusters' is passed as clustering"
+                             " 'method'"
                              " the 'num_clusters' argument must be specified.")
 
         if isinstance(initial_guess, str):
@@ -362,15 +528,24 @@ def _set_initial_guess_cd(self,
                 if method == "aglomerative":
                     self.init_guess = np.array(
                         [k for k in np.arange(self.n_nodes, dtype=np.int32)])
-                elif method == "divisive":
+                elif method == "fixed-clusters":
                     self.init_guess = np.random.randint(
                         low=num_clusters,
                         size=self.n_nodes)
-            elif initial_guess == "common-neighbours":
+            elif (initial_guess == "common-neigh-weak") or \
+                    (initial_guess == "common-neighbours"):
+                if method == "aglomerative":
+                    self.init_guess = ax.common_neigh_init_guess_weak(
+                        self.adjacency)
+                elif method == "fixed-clusters":
+                    self.init_guess = ax.fixed_clusters_init_guess_cn(
+                        adjacency=self.adjacency,
+                        n_clust=num_clusters)
+            elif initial_guess == "common-neigh-strong":
                 if method == "aglomerative":
-                    self.init_guess = ax.common_neigh_init_guess(
+                    self.init_guess = ax.common_neigh_init_guess_strong(
                         self.adjacency)
-                elif method == "divisive":
+                elif method == "fixed-clusters":
                     self.init_guess = ax.fixed_clusters_init_guess_cn(
                         adjacency=self.adjacency,
                         n_clust=num_clusters)
@@ -379,7 +554,8 @@ def _set_initial_guess_cd(self,
                     "The 'initial_guess' selected is not a valid."
                     "Initial guess can be an array specifying nodes membership"
                     " or an initialisation method ['common-neighbours',"
-                    " random]. For more details see documentation.")
+                    " 'random', 'common-neigh-weak', 'common-neigh-strong']."
+                    " For more details see documentation.")
 
         elif isinstance(initial_guess, np.ndarray):
             self.init_guess = initial_guess.astype(np.int32)
@@ -391,11 +567,11 @@ def _set_initial_guess_cd(self,
                 "The length of the initial guess provided is different from"
                 " the network number of nodes.")
 
-        if (method == "divisive" and
+        if (method == "fixed-clusters" and
                 np.unique(self.init_guess).shape[0] != num_clusters):
             raise ValueError("The number of clusters of a custom initial guess"
                              " must coincide with 'num_clusters' when the "
-                             " divisive method is applied.")
+                             " fixed-clusters method is applied.")
 
     def _set_solved_problem(self, sol):
         self.solution = sol[0]
diff --git a/src/surprisemes/__init__.py b/src/surprisememore/__init__.py
similarity index 100%
rename from src/surprisemes/__init__.py
rename to src/surprisememore/__init__.py
diff --git a/src/surprisememore/__pycache__/Directed_Graph_Class.cpython-37.pyc b/src/surprisememore/__pycache__/Directed_Graph_Class.cpython-37.pyc
new file mode 100644
index 00000000..784e5c15
Binary files /dev/null and b/src/surprisememore/__pycache__/Directed_Graph_Class.cpython-37.pyc differ
diff --git a/src/surprisememore/__pycache__/Undirected_Graph_Class.cpython-37.pyc b/src/surprisememore/__pycache__/Undirected_Graph_Class.cpython-37.pyc
new file mode 100644
index 00000000..512fc942
Binary files /dev/null and b/src/surprisememore/__pycache__/Undirected_Graph_Class.cpython-37.pyc differ
diff --git a/src/surprisemes/__pycache__/__init__.cpython-37.pyc b/src/surprisememore/__pycache__/__init__.cpython-37.pyc
similarity index 100%
rename from src/surprisemes/__pycache__/__init__.cpython-37.pyc
rename to src/surprisememore/__pycache__/__init__.cpython-37.pyc
diff --git a/src/surprisememore/__pycache__/auxiliary_function.cpython-37.pyc b/src/surprisememore/__pycache__/auxiliary_function.cpython-37.pyc
new file mode 100644
index 00000000..a1858e2c
Binary files /dev/null and b/src/surprisememore/__pycache__/auxiliary_function.cpython-37.pyc differ
diff --git a/src/surprisememore/__pycache__/comdet_functions.cpython-37.pyc b/src/surprisememore/__pycache__/comdet_functions.cpython-37.pyc
new file mode 100644
index 00000000..f127b0ca
Binary files /dev/null and b/src/surprisememore/__pycache__/comdet_functions.cpython-37.pyc differ
diff --git a/src/surprisememore/__pycache__/cp_functions.cpython-37.pyc b/src/surprisememore/__pycache__/cp_functions.cpython-37.pyc
new file mode 100644
index 00000000..b6e41dbc
Binary files /dev/null and b/src/surprisememore/__pycache__/cp_functions.cpython-37.pyc differ
diff --git a/src/surprisememore/__pycache__/solver.cpython-37.pyc b/src/surprisememore/__pycache__/solver.cpython-37.pyc
new file mode 100644
index 00000000..e8f21b13
Binary files /dev/null and b/src/surprisememore/__pycache__/solver.cpython-37.pyc differ
diff --git a/src/surprisememore/auxiliary_function.py b/src/surprisememore/auxiliary_function.py
new file mode 100644
index 00000000..4e7cfd09
--- /dev/null
+++ b/src/surprisememore/auxiliary_function.py
@@ -0,0 +1,942 @@
+import networkx as nx
+import numpy as np
+import scipy
+from numba import jit
+from scipy.sparse import isspmatrix
+from scipy.special import comb
+
+from . import comdet_functions as cd
+from . import cp_functions as cp
+
+
+def compute_neighbours(adj):
+    lista_neigh = []
+    for ii in np.arange(adj.shape[0]):
+        lista_neigh.append(adj[ii, :].nonzero()[0])
+    return lista_neigh
+
+
+@jit(nopython=True)
+def compute_cn(adjacency):
+    """ Computes common neighbours table, each entry i,j of this table is the
+     number of common neighbours between i and j.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: numpy.ndarray
+    :return: Common neighbours table.
+    :rtype: numpy.ndarray
+    """
+    cn_table = np.zeros_like(adjacency)
+    for i in np.arange(adjacency.shape[0]):
+        neighbour_i = (adjacency[i, :] + adjacency[:, i]).astype(np.bool_)
+        for j in np.arange(i + 1, adjacency.shape[0]):
+            neighbour_j = (adjacency[j, :] + adjacency[:, j]).astype(np.bool_)
+            cn_table[i, j] = cn_table[j, i] = np.multiply(neighbour_i,
+                                                          neighbour_j).sum()
+    return cn_table
+
+
+@jit(nopython=True)
+def common_neigh_init_guess_strong(adjacency):
+    """Generates a preprocessed initial guess based on the common neighbours
+     of nodes. It makes a stronger aggregation of nodes based on
+      the common neighbours similarity.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: numpy.ndarray
+    :return: Initial guess for nodes memberships.
+    :rtype: np.array
+    """
+    cn_table = compute_cn(adjacency)
+    memberships = np.array(
+        [k for k in np.arange(adjacency.shape[0], dtype=np.int32)])
+    argsorted = np.argsort(adjacency.astype(np.bool_).sum(axis=1))[::-1]
+    for aux_node1 in argsorted:
+        aux_tmp = memberships == aux_node1
+        memberships[aux_tmp] = memberships[np.argmax(cn_table[aux_node1])]
+    return memberships
+
+
+@jit(nopython=True)
+def common_neigh_init_guess_weak(adjacency):
+    """Generates a preprocessed initial guess based on the common neighbours
+     of nodes. It makes a weaker aggregation of nodes based on
+      the common neighbours similarity.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: numpy.ndarray
+    :return: Initial guess for nodes memberships.
+    :rtype: np.array
+    """
+    cn_table = compute_cn(adjacency)
+    memberships = np.array(
+        [k for k in np.arange(adjacency.shape[0], dtype=np.int32)])
+    degree = (adjacency.astype(np.bool_).sum(axis=1)
+              + adjacency.astype(np.bool_).sum(axis=0))
+    avg_degree = np.mean(degree)
+    argsorted = np.argsort(degree)[::-1]
+    for aux_node1 in argsorted:
+        if degree[aux_node1] >= avg_degree:
+            aux_tmp = memberships == aux_node1
+            memberships[aux_tmp] = memberships[np.argmax(cn_table[aux_node1])]
+    return memberships
+
+
+def eigenvector_init_guess(adjacency, is_directed):
+    """Generates an initial guess for core periphery detection method: nodes
+    with higher eigenvector centrality are in the core.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: np.ndarray
+    :param is_directed: True if the network is directed.
+    :type is_directed: bool
+    :return: Initial guess.
+    :rtype: np.ndarray
+    """
+    # TODO: Vedere come funziona la parte pesata
+
+    n_nodes = adjacency.shape[0]
+    aux_nodes = int(np.ceil((n_nodes * 5) / 100))
+    if is_directed:
+        graph = nx.from_numpy_array(adjacency, create_using=nx.DiGraph)
+        centra = nx.eigenvector_centrality_numpy(graph)
+        centra1 = np.array([centra[key] for key in centra])
+        membership = np.ones_like(centra1, dtype=np.int32)
+        membership[np.argsort(centra1)[::-1][:aux_nodes]] = 0
+
+    else:
+        graph = nx.from_numpy_array(adjacency, create_using=nx.Graph)
+        centra = nx.eigenvector_centrality_numpy(graph)
+        centra1 = np.array([centra[key] for key in centra])
+        membership = np.ones_like(centra1, dtype=np.int32)
+        membership[np.argsort(centra1)[::-1][:aux_nodes]] = 0
+
+    return membership
+
+
+def fixed_clusters_init_guess_cn(adjacency, n_clust):
+    """ Generates an intial guess with a fixed number 'n' of clusters.
+    Nodes are organised in clusters based on the number of common neighbors.
+    The starting members of clusters are the 'n' nodes with higher
+    degrees/strengths.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: numpy.ndarray
+    :param n_clust: Partitions number.
+    :type n_clust: int
+    :return: Initial guess.
+    :rtype: numpy.ndarray
+    """
+
+    aux_memb = np.ones(adjacency.shape[0], dtype=np.int32) * (n_clust - 1)
+
+    cn = compute_cn(adjacency)
+    degree = adjacency.astype(np.bool_).sum(axis=1) + adjacency.astype(
+        np.bool_).sum(axis=0)
+    avg_degree = np.mean(degree)
+    degree_indices_g = np.nonzero(degree > 2)[0]
+    degree_indices_l = np.nonzero(degree <= 2)[0]
+
+    arg_max = np.argmax(degree[degree_indices_g])
+    clust_element = degree_indices_g[arg_max]
+
+    cluster_count = 0
+    while cluster_count != n_clust - 1:
+        aux_memb[clust_element] = cluster_count
+        degree_indices_g = np.delete(degree_indices_g, arg_max)
+        if len(degree_indices_g) == 0:
+            break
+        arg_max = np.argmin(cn[clust_element][degree_indices_g])
+        clust_element = degree_indices_g[arg_max]
+        cluster_count += 1
+
+    if np.unique(aux_memb).shape[0] < n_clust - 1:
+        cluster_count += 1
+        arg_max = np.argmax(degree[degree_indices_l])
+        clust_element = degree_indices_l[arg_max]
+        while cluster_count != n_clust - 1:
+            aux_memb[clust_element] = cluster_count
+            degree_indices_l = np.delete(degree_indices_l, arg_max)
+            if len(degree_indices_l) == 0:
+                raise ValueError(
+                    "The number of clusters is higher thant the nodes number.")
+            arg_max = np.argmin(cn[clust_element][degree_indices_l])
+            clust_element = np.argmin(cn[clust_element][degree_indices_l])
+            cluster_count += 1
+
+    aux = np.nonzero(aux_memb == n_clust - 1)[0]
+    np.random.shuffle(aux)
+    for node in aux:
+        if degree[node] < avg_degree:
+            continue
+        aux_list = np.nonzero(aux_memb != n_clust - 1)[0]
+        node_index = aux_list[np.argmax(cn[node, aux_list])]
+        if isinstance(node_index, np.ndarray):
+            node_index = np.random.choice(node_index)
+        aux_memb[node] = aux_memb[node_index]
+
+    return aux_memb
+
+
+def compute_degree(a, is_directed):
+    """Returns matrix *a* degree sequence.
+
+    :param a: Matrix.
+    :type a:  numpy.ndarray
+    :param is_directed: True if the matrix is directed.
+    :type is_directed: bool
+    :return: Degree sequence.
+    :rtype: numpy.ndarray.
+    """
+    # if the matrix is a numpy array
+    if is_directed:
+        if type(a) == np.ndarray:
+            return np.sum(a > 0, 0), np.sum(a > 0, 1)
+        # if the matrix is a scipy sparse matrix
+        elif isspmatrix(a):
+            return np.sum(a > 0, 0).A1, np.sum(a > 0, 1).A1
+    else:
+        if type(a) == np.ndarray:
+            return np.sum(a > 0, 1)
+        # if the matrix is a scipy sparse matrix
+        elif isspmatrix(a):
+            return np.sum(a > 0, 1).A1
+
+
+def compute_strength(a, is_directed):
+    """Returns matrix *a* strength sequence.
+
+    :param a: Matrix.
+    :type a: numpy.ndarray
+    :param is_directed: True if the matrix is directed.
+    :type is_directed: bool
+    :return: Strength sequence.
+    :rtype: numpy.ndarray
+    """
+    if is_directed:
+        # if the matrix is a numpy array
+        if type(a) == np.ndarray:
+            return np.sum(a, 0), np.sum(a, 1)
+        # if the matrix is a scipy sparse matrix
+        elif isspmatrix(a):
+            return np.sum(a, 0).A1, np.sum(a, 1).A1
+    else:
+        # if the matrix is a numpy array
+        if type(a) == np.ndarray:
+            return np.sum(a, 1)
+        # if the matrix is a scipy sparse matrix
+        elif isspmatrix(a):
+            return np.sum(a, 1).A1
+
+
+def from_edgelist(edgelist, is_sparse, is_directed):
+    """Returns np.ndarray or scipy.sparse matrix from edgelist.
+
+    :param edgelist: List of edges, eache edge must be given as a 2-tuples
+     (u,v).
+    :type edgelist: list or numpy.ndarray
+    :param is_sparse: If true the returned matrix is sparse.
+    :type is_sparse: bool
+    :param is_directed: If true the graph is directed.
+    :type is_directed: bool
+    :return: Adjacency matrix.
+    :rtype: numpy.ndarray or scipy.sparse
+    """
+    # TODO: vedere che tipo di sparse e'
+    if is_directed:
+        g = nx.DiGraph()
+    else:
+        g = nx.Graph()
+    g.add_edges_from(edgelist)
+    if is_sparse:
+        return nx.to_scipy_sparse_matrix(g)
+    else:
+        return nx.to_numpy_array(g)
+
+
+def from_weighted_edgelist(edgelist, is_sparse, is_directed):
+    """Returns np.ndarray or scipy.sparse matrix from edgelist.
+
+    :param edgelist: List of weighted edges, eache edge must be given as a
+     3-tuples (u,v,w).
+    :type edgelist: [type]
+    :param is_sparse: If true the returned matrix is sparse.
+    :type is_sparse: bool
+    :param is_directed: If true the graph is directed.
+    :type is_directed: bool
+    :return: Weighted adjacency matrix.
+    :rtype: numpy.ndarray or scipy.sparse
+    """
+    if is_directed:
+        g = nx.DiGraph()
+    else:
+        g = nx.Graph()
+    g.add_weighted_edges_from(edgelist)
+    if is_sparse:
+        return nx.to_scipy_sparse_matrix(g)
+    else:
+        return nx.to_numpy_array(g)
+
+
+def check_symmetric(a, is_sparse, rtol=1e-05, atol=1e-08):
+    """Checks if the matrix is symmetric.
+
+    :param a: Matrix.
+    :type a: numpy.ndarray or scipy.sparse
+    :param is_sparse: If true the matrix is sparse.
+    :type is_sparse: bool
+    :param rtol: Tuning parameter, defaults to 1e-05.
+    :type rtol: float, optional
+    :param atol: Tuning parameter, defaults to 1e-08.
+    :type atol: float, optional
+    :return: True if the matrix is symmetric.
+    :rtype: bool
+    """
+    if is_sparse:
+        return np.all(np.abs(a - a.T) < atol)
+    else:
+        return np.allclose(a, a.T, rtol=rtol, atol=atol)
+
+
+def check_adjacency(adjacency, is_sparse, is_directed):
+    """Functions checking the _validty_ of the adjacency matrix.
+
+    :param adjacency: Adjacency matrix.
+    :type adjacency: numpy.ndarray or scipy.sparse
+    :param is_sparse: If true the matrix is sparse.
+    :type is_sparse: bool
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :raises TypeError: Matrix not square.
+    :raises ValueError: Negative entries.
+    :raises TypeError: Matrix not symmetric.
+    """
+    if adjacency.shape[0] != adjacency.shape[1]:
+        raise TypeError(
+            "Adjacency matrix must be square. If you are passing an edgelist"
+            " use the positional argument 'edgelist='.")
+    if np.sum(adjacency < 0):
+        raise ValueError(
+            "The adjacency matrix entries must be positive."
+        )
+    if (not check_symmetric(adjacency, is_sparse)) and (not is_directed):
+        raise TypeError(
+            "The adjacency matrix seems to be not symmetric, we suggest to use"
+            " 'DirectedGraphClass'.")
+
+
+@jit(nopython=True, fastmath=True)
+def sumLogProbabilities(nextlogp, logp):
+    if nextlogp == 0:
+        stop = True
+    else:
+        stop = False
+        if nextlogp > logp:
+            common = nextlogp
+            diffexponent = logp - common
+        else:
+            common = logp
+            diffexponent = nextlogp - common
+
+        logp = common + ((np.log10(1 + 10 ** diffexponent)) / np.log10(10))
+
+        if (nextlogp - logp) > -4:
+            stop = True
+
+    return logp, stop
+
+
+@jit(nopython=True, fastmath=True)
+def logc(n, k):
+    if k == n:
+        return 0
+    elif (n > 1000) & (k > 1000):  # Stirling's binomial coeff approximation
+        return logStirFac(n) - logStirFac(k) - logStirFac(n - k)
+    else:
+        t = n - k
+        if t < k:
+            t = k
+        logC = sumRange(t + 1, n) - sumFactorial(n - t)
+        return logC
+
+
+@jit(nopython=True, fastmath=True)
+def logStirFac(n):
+    if n <= 1:
+        return 1.0
+    else:
+        return -n + n * np.log10(n) + np.log10(
+            n * (1 + 4.0 * n * (1.0 + 2.0 * n))) / 6.0 + np.log10(np.pi) / 2.0
+
+
+@jit(nopython=True, fastmath=True)
+def sumRange(xmin, xmax):
+    """[summary]
+
+    :param xmin: [description]
+    :type xmin: [type]
+    :param xmax: [description]
+    :type xmax: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    csum = 0
+    for i in np.arange(xmin, xmax + 1):
+        csum += np.log10(i)
+    return csum
+
+
+@jit(nopython=True, fastmath=True)
+def sumFactorial(n):
+    csum = 0
+    if n > 1:
+        for i in np.arange(2, n + 1):
+            csum += np.log10(i)
+    return csum
+
+
+def shuffled_edges(adjacency_matrix, is_directed):
+    """Shuffles edges randomly.
+
+    :param adjacency_matrix: Matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param is_directed: True if graph is directed.
+    :type is_directed: bool
+    :return: Shuffled edgelist.
+    :rtype: mumpy.ndarray
+    """
+    adj = adjacency_matrix.astype(bool).astype(np.int16)
+    if not is_directed:
+        adj = np.triu(adj)
+    edges = np.stack(adj.nonzero(), axis=-1)
+    np.random.shuffle(edges)
+    shuff_edges = edges.astype(np.int32)
+    return shuff_edges
+
+
+def jaccard_sorted_edges(adjacency_matrix):
+    """Returns edges ordered based on jaccard index.
+
+    :param adjacency_matrix: Matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :return: Ordered edgelist.
+    :rtype: numpy.ndarray
+    """
+    G = nx.from_numpy_matrix(adjacency_matrix)
+    jacc = nx.jaccard_coefficient(G, ebunch=G.edges())
+    jacc_array = []
+    for u, v, p in jacc:
+        jacc_array += [[u, v, p]]
+    jacc_array = np.array(jacc_array)
+    jacc_array = jacc_array[jacc_array[:, 2].argsort()][::-1]
+    sorted_edges = jacc_array[:, 0:2]
+    sorted_edges = sorted_edges.astype(np.int32)
+    return sorted_edges
+
+
+def evaluate_surprise_clust_bin(adjacency_matrix,
+                                cluster_assignment,
+                                is_directed):
+    """Calculates the logarithm of the surprise given the current partitions
+     for a binary network.
+
+    :param adjacency_matrix: Binary adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster links
+        p = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment)
+        p = int(p)
+        # All the possible intracluster links
+        M = cd.calculate_possible_intracluster_links(cluster_assignment,
+                                                     is_directed)
+        # Observed links
+        m = np.sum(adjacency_matrix.astype(bool))
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        F = n * (n - 1)
+    else:
+        # intracluster links
+        p = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment)
+        p = int(p / 2)
+        # All the possible intracluster links
+        M = int(cd.calculate_possible_intracluster_links(cluster_assignment,
+                                                         is_directed))
+        # Observed links
+        m = np.sum(adjacency_matrix.astype(bool)) / 2
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        F = int((n * (n - 1)) / 2)
+
+    surprise = surprise_hypergeometric(F, p, M, m)
+    return surprise
+
+
+def surprise_hypergeometric(F, p, M, m):
+    surprise = 0
+    min_p = min(M, m)
+    for p_loop in np.arange(p, min_p + 1):
+        surprise += (comb(M, p_loop, exact=True) * comb(
+            F - M, m - p_loop,
+            exact=True)) / comb(F,
+                                m,
+                                exact=True)
+    return surprise
+
+
+def evaluate_surprise_clust_weigh(adjacency_matrix,
+                                  cluster_assignment,
+                                  is_directed):
+    """Calculates the logarithm of the surprise given the current partitions
+     for a weighted network.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        w = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment)
+        # intracluster possible links
+        Vi = cd.calculate_possible_intracluster_links(cluster_assignment,
+                                                      is_directed)
+        # Total Weight
+        W = np.sum(adjacency_matrix)
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        V = n * (n - 1)
+        # extracluster links
+        Ve = V - Vi
+    else:
+        # intracluster weights
+        w = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment) / 2
+        # intracluster possible links
+        Vi = cd.calculate_possible_intracluster_links(cluster_assignment,
+                                                      is_directed)
+        # Total Weight
+        W = np.sum(adjacency_matrix) / 2
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        V = int((n * (n - 1)) / 2)
+        # extracluster links
+        Ve = V - Vi
+
+    surprise = surprise_negative_hypergeometric(Vi, w, Ve, W, V)
+    return surprise
+
+
+def surprise_negative_hypergeometric(Vi, w, Ve, W, V):
+    """Computes the negative hypergeometric distribution.
+    """
+    surprise = 0
+    for w_loop in range(w, W):
+        surprise += ((comb(Vi + w_loop - 1, w_loop, exact=True) * comb(
+            Ve + W - w_loop, W - w_loop, exact=True)) /
+                     comb(V + W, W, exact=True))
+    return surprise
+
+
+def evaluate_surprise_cp_bin(adjacency_matrix,
+                             cluster_assignment,
+                             is_directed):
+    """Computes core-periphery binary log-surprise given a certain nodes'
+     partitioning.
+
+    :param adjacency_matrix: Binary adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Core periphery assigments.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise
+    :rtype: float
+    """
+    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
+    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
+
+    if is_directed:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = n_c * (n_c - 1)
+        p_x = n_c * n_x * 2
+
+        l_c = cp.compute_sum(adjacency_matrix, core_nodes, core_nodes)
+        l_x = cp.compute_sum(adjacency_matrix, core_nodes,
+                             periphery_nodes) + cp.compute_sum(
+                                                         adjacency_matrix,
+                                                         periphery_nodes,
+                                                         core_nodes)
+
+        l_t = np.sum(adjacency_matrix)
+        n = n_c + n_x
+        p = n * (n - 1)
+
+    else:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = (n_c * (n_c - 1)) / 2
+        p_x = n_c * n_x
+
+        l_c = cp.compute_sum(adjacency_matrix, core_nodes, core_nodes) / 2
+        l_x = (cp.compute_sum(adjacency_matrix, core_nodes,
+                              periphery_nodes) + cp.compute_sum(
+                                                          adjacency_matrix,
+                                                          periphery_nodes,
+                                                          core_nodes)) / 2
+
+        l_t = np.sum(adjacency_matrix) / 2
+        n = n_c + n_x
+        p = (n * (n - 1)) / 2
+
+    if (p_c + p_x) < (l_c + l_x):
+        return 0
+
+    surprise = surprise_bipartite_cp_bin(p, p_c, p_x, l_t, l_c, l_x)
+    return surprise
+
+
+@jit(forceobj=True)
+def surprise_bipartite_cp_bin(p, p_c, p_x, l, l_c, l_x):
+    surprise = 0
+    aux_first = 0
+    for l_c_loop in range(l_c, p_c + 1):
+        aux_first_temp = aux_first
+        for l_x_loop in range(l_x, p_x + 1):
+            if l_c_loop + l_x_loop > l:
+                continue
+            aux = multihyperprobability(p, p_c, p_x,
+                                        l, l_c_loop,
+                                        l_x_loop)
+            surprise += aux
+            if surprise == 0:
+                break
+            if aux/surprise < 1e-3:
+                break
+
+        aux_first = surprise
+        if aux_first - aux_first_temp:
+            if ((aux_first - aux_first_temp) / aux_first) < 1e-3:
+                # pass
+                break
+
+    return surprise
+
+
+# @jit(nopython=True)
+def multihyperprobability(p, p_c, p_x, l, l_c, l_x):
+    """Computes the logarithm of the Multinomial Hypergeometric
+     distribution."""
+    logh = comb(p_c, l_c, True) * comb(p_x, l_x) + comb(
+        p - p_c - p_x,
+        l - l_c - l_x) - comb(p, l)
+    return logh
+
+
+def evaluate_surprise_cp_enh(adjacency_matrix,
+                             cluster_assignment,
+                             is_directed):
+    """Computes core-periphery weighted surprise given a
+     certain nodes' partitioning.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Core periphery assigments.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise
+    :rtype: float
+    """
+    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
+    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
+
+    if is_directed:
+        n_o = core_nodes.shape[0]
+        n_p = periphery_nodes.shape[0]
+        V_o = n_o * (n_o - 1)
+        V_c = n_o * n_p * 2
+
+        l_o, w_o = cp.compute_sum_enh(adjacency_matrix, core_nodes, core_nodes)
+        l_c, w_c = cp.compute_sum_enh(adjacency_matrix,
+                                      core_nodes,
+                                      periphery_nodes) + cp.compute_sum_enh(
+            adjacency_matrix,
+            periphery_nodes,
+            core_nodes)
+        L = np.sum(adjacency_matrix.astype(bool))
+        W = np.sum(adjacency_matrix)
+        # w_p = W - w_o - w_c
+        n = n_o + n_p
+        V = n * (n - 1)
+
+    else:
+        n_o = core_nodes.shape[0]
+        n_p = periphery_nodes.shape[0]
+        V_o = n_o * (n_o - 1) / 2
+        V_c = n_o * n_p
+
+        l_o, w_o = (cp.compute_sum_enh(adjacency_matrix, core_nodes, core_nodes))
+        l_c1, w_c1 = cp.compute_sum_enh(adjacency_matrix,
+                                        core_nodes,
+                                        periphery_nodes)
+        l_c2, w_c2 = cp.compute_sum_enh(adjacency_matrix,
+                                        periphery_nodes,
+                                        core_nodes)
+
+        l_o = l_o / 2
+        w_o = w_o / 2
+        l_c = (l_c1 + l_c2) / 2
+        w_c = (w_c1 + w_c2) / 2
+
+        L = np.sum(adjacency_matrix.astype(bool)) / 2
+        W = np.sum(adjacency_matrix) / 2
+        # w_p = (W - w_o - w_c) / 2
+        n = n_o + n_p
+        V = n * (n - 1) / 2
+
+    # print("V_o", V_o, "l_o", l_o, "V_c", V_c, "l_c", l_c,
+    #      "w_o", w_o, "w_c", w_c, "V", V, "L", L, "W", W)
+
+    surprise = surprise_bipartite_cp_enh(V_o, l_o, V_c, l_c,
+                                         w_o, w_c, V, L, W)
+    return surprise
+
+
+@jit(forceobj=True)
+def surprise_bipartite_cp_enh(V_o, l_o, V_c, l_c, w_o, w_c, V, L, W):
+    surprise = 0
+
+    min_l_o = min(L, V_o + V_c)
+    aux_first = 0
+    aux_second = 0
+    aux_third = 0
+    for l_o_loop in np.arange(l_o, min_l_o + 1):
+        aux_first_temp = aux_first
+        for l_c_loop in np.arange(l_c, min_l_o + 1 - l_o_loop):
+            aux_second_temp = aux_second
+            for w_o_loop in np.arange(w_o, W + 1):
+                aux_third_temp = aux_third
+                for w_c_loop in np.arange(w_c, W + 1 - w_o_loop):
+                    aux = logmulti_hyperprobability_weightenh(
+                                                            V_o, l_o_loop,
+                                                            V_c, l_c_loop,
+                                                            w_o_loop, w_c_loop,
+                                                            V, L, W)
+                    surprise += aux
+                    # print(surprise)
+                    if surprise == 0:
+                        break
+                    if aux / surprise < 1e-3:
+                        break
+
+                aux_third = surprise
+                if aux_third - aux_third_temp:
+                    if ((aux_third - aux_third_temp) / aux_third) < 1e-3:
+                        # pass
+                        break
+                else:
+                    break
+
+            aux_second = aux_third
+            if aux_second - aux_second_temp:
+                if ((aux_second - aux_second_temp) / aux_second) < 1e-3:
+                    # pass
+                    break
+            else:
+                break
+
+        aux_first = aux_second
+        if aux_first - aux_first_temp:
+            if ((aux_first - aux_first_temp) / aux_first) < 1e-4:
+                # pass
+                break
+        else:
+            break
+
+    return surprise
+
+
+@jit(forceobj=True)
+def logmulti_hyperprobability_weightenh(V_o, l_o, V_c, l_c, w_o, w_c, V, L, W):
+    """Computes the of the Negative Multinomial Hypergeometric
+     distribution."""
+    aux1 = (comb(V_o, l_o, exact=True) * comb(V_c, l_c, exact=True) * comb(
+        V - (V_o + V_c), L - (l_o + l_c), exact=True)) / comb(V, L, exact=True)
+    aux2 = (comb(w_o - 1, l_o - 1, exact=True) * comb(w_c - 1, l_c - 1,
+                                                      exact=True) * comb(
+        W - (w_o + w_c) - 1, L - (l_o + l_c) - 1, exact=True)) / comb(W - 1,
+                                                                      W - L,
+                                                                      exact=True)
+    return aux1 * aux2
+
+
+def evaluate_surprise_community_enh(
+        adjacency_matrix,
+        cluster_assignment,
+        is_directed):
+    """Calculates surprise given the current partitions for a weighted network.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        l_o, w_o = cd.intracluster_links_enh(adjacency_matrix,
+                                             cluster_assignment)
+        # intracluster possible links
+        V_o = cd.calculate_possible_intracluster_links(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        W = np.sum(adjacency_matrix)
+        L = np.sum(adjacency_matrix.astype(bool))
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        V = n * (n - 1)
+        # extracluster links
+        # inter_links = V - V_o
+    else:
+        # intracluster weights
+        l_o, w_o = cd.intracluster_links_enh(adjacency_matrix,
+                                             cluster_assignment)
+        l_o = l_o / 2
+        w_o = w_o / 2
+        # intracluster possible links
+        V_o = cd.calculate_possible_intracluster_links(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        W = np.sum(adjacency_matrix) / 2
+        L = np.sum(adjacency_matrix.astype(bool)) / 2
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        V = int((n * (n - 1)) / 2)
+        # extracluster links
+        # inter_links = V - V_o
+
+    # print("V_0", V_o, "l_0", l_o, "w_0", w_o, "V", V, "L", L, "W", W)
+
+    surprise = surprise_clust_enh(V_o, l_o, w_o, V, L, W)
+    return surprise
+
+
+@jit(forceobj=True)
+def surprise_clust_enh(V_o, l_o, w_o, V, L, W):
+    min_l_loop = min(L, V_o)
+
+    surprise = 0.0
+    aux_first = 0.0
+    # print("l_o", l_o, "min l", min_l_loop,"w_0", w_o, "W-L", W)
+    for l_loop in range(l_o, min_l_loop + 1):
+        aux_first_temp = aux_first
+        for w_loop in range(w_o - l_loop + l_o, W - L + l_o + 1):
+            if w_loop <= 0:
+                continue
+            # print(l_loop,  w_loop)
+            aux = logenhancedhypergeometric(V_o, l_loop, w_loop, V, L, W)
+
+            if np.isnan(aux):
+                break
+            surprise += aux
+            # print(aux, surprise)
+            if surprise == 0:
+                break
+            if aux / surprise <= 1e-3:
+                break
+        aux_first = surprise
+        if aux_first - aux_first_temp:
+            if ((aux_first - aux_first_temp) / aux_first) < 1e-4:
+                # pass
+                break
+        else:
+            break
+
+    return surprise
+
+
+@jit(forceobj=True)
+def logenhancedhypergeometric(V_o, l_o, w_o, V, L, W):
+    if l_o < L:
+        aux1 = (comb(V_o, l_o, True) * comb(V - V_o, L - l_o, True)) / comb(V, L, True)
+        aux2 = (comb(w_o - 1, w_o - l_o, True) * comb(W - w_o - 1, (W - L) - (w_o - l_o), True)) / comb(W - 1, W - L, True)
+    else:
+        aux1 = (comb(V_o, l_o, True) / comb(V, L, True))
+        aux2 = comb(w_o - 1, w_o - L, True)
+    return aux1 * aux2
+
+
+def evaluate_surprise_com_det_continuous(
+        adjacency_matrix,
+        cluster_assignment,
+        is_directed):
+    """Calculates the logarithm of the continuous surprise given
+     the current partitions for a weighted network.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        w = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment)
+        # intracluster possible links
+        poss_intr_links = cd.calculate_possible_intracluster_links(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = np.sum(adjacency_matrix)
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        poss_links = n * (n - 1)
+        # extracluster links
+        # inter_links = poss_links - poss_intr_links
+    else:
+        # intracluster weights
+        w = cd.intracluster_links(adjacency_matrix,
+                                  cluster_assignment) / 2
+        # intracluster possible links
+        poss_intr_links = cd.calculate_possible_intracluster_links(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = np.sum(adjacency_matrix) / 2
+        # Possible links
+        n = adjacency_matrix.shape[0]
+        poss_links = int((n * (n - 1)) / 2)
+        # extracluster links
+        # inter_links = poss_links - poss_intr_links
+
+    surprise = cd.continuous_surprise_clust(
+        V=poss_links,
+        W=tot_weights,
+        w_o=w,
+        V_o=poss_intr_links)
+
+    return surprise
diff --git a/src/surprisememore/comdet_functions.py b/src/surprisememore/comdet_functions.py
new file mode 100644
index 00000000..45b91670
--- /dev/null
+++ b/src/surprisememore/comdet_functions.py
@@ -0,0 +1,888 @@
+from . import auxiliary_function as ax
+
+import numpy as np
+from mpmath import mp, ncdf, log10, quad, beta, mpf, sqrt, pi
+from numba import jit
+
+mp.dps = 100
+mp.pretty = True
+
+
+def calculate_possible_intracluster_links(partitions, is_directed):
+    """Computes the number of possible links, given nodes memberships.
+
+    :param partitions: Nodes memberships.
+    :type partitions: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Total number of possible edges.
+    :rtype: float
+    """
+    _, counts = np.unique(partitions, return_counts=True)
+    nr_poss_intr_clust_links = np.sum(counts * (counts - 1))
+    if is_directed:
+        return nr_poss_intr_clust_links
+    else:
+        return nr_poss_intr_clust_links / 2
+
+
+@jit(nopython=True)
+def calculate_possible_intracluster_links_new(partitions, is_directed):
+    """Computes the number of possible links, given nodes memberships.
+    Faster implementation compiled in "nopython" mode.
+
+    :param partitions: Nodes memberships.
+    :type partitions: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Total number of possible edges.
+    :rtype: float
+    """
+    counts = np.bincount(partitions)
+    nr_poss_intr_clust_links = np.sum(counts * (counts - 1))
+    if is_directed:
+        return nr_poss_intr_clust_links
+    else:
+        return nr_poss_intr_clust_links / 2
+
+
+@jit(nopython=True)
+def intracluster_links_new(adj, clust_labels, partitions):
+    """Computes intracluster links or weights. New implementation
+
+    :param adj: Adjacency matrix.
+    :type adj: numpy.ndarray
+    :param clust_labels: Labels of changed clusters.
+    :type clust_labels: numpy.ndarray
+    :param partitions: Nodes memberships.
+    :type partitions: numpy.ndarray
+    :return: Number of intra-cluster links/weights.                
+    :rtype: numpy.ndarray
+    """
+    # print(clust_labels, clust_labels.shape, partitions)
+    nr_intr_clust_links = np.zeros(clust_labels.shape[0])
+    for ii, lab in enumerate(clust_labels):
+        indices = np.where(partitions == lab)[0]
+        nr_intr_clust_links[ii] = intracluster_links_aux(adj, indices)
+    return nr_intr_clust_links
+
+
+@jit(nopython=True)
+def intracluster_links_aux(adj, indices):
+    """Computes intra-cluster links or weights given nodes indices.
+
+    :param adj: [description]
+    :type adj: [type]
+    :param indices: [description]
+    :type indices: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    n_links = 0
+    for ii in indices:
+        for jj in indices:
+            n_links += adj[ii, jj]
+    return n_links
+
+
+def calculate_surprise_logsum_clust_bin_new(
+        adjacency_matrix,
+        cluster_assignment,
+        mem_intr_link,
+        clust_labels,
+        args,
+        approx,
+        is_directed):
+    """Calculates the logarithm of the surprise given the current partitions
+     for a binary network. New faster implementation reducing the number of
+      redundant computations.
+
+    :param adjacency_matrix: Binary adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.                                                       
+    :type cluster_assignment: numpy.ndarray
+    :param mem_intr_link: Intracluster links per cluster
+    :type mem_intr_link: np.array
+    :param clust_labels: Labels of changed clusters.
+    :type clust_labels: np.array
+    :param args: (Observed links, nodes number, possible links)
+    :type args: (int, int)
+    :param approx:
+    :type approx:
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.                                                                              
+    :rtype: float                                                                                       
+    """
+    if is_directed:
+        # intracluster links
+        if len(clust_labels):
+            int_links = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, int_links):
+                mem_intr_link[0][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[0])
+        int_links = int(p)
+        # All the possible intracluster links                                                           
+        poss_int_links = int(calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed))
+        # Observed links                                                                                
+        obs_links = int(args[0])
+        # Possible link
+        poss_links = int(args[2])
+    else:
+        # intracluster links
+        if len(clust_labels):
+            int_links = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, int_links):
+                mem_intr_link[0][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[0])
+        int_links = int(p / 2)
+        # All the possible intracluster links                                                           
+        poss_int_links = int(calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed))
+        # Observed links                                                                                
+        obs_links = int(args[0] / 2)
+        # Possible links
+        poss_links = int(args[2] / 2)
+
+    if int_links == 0:
+        return 0, mem_intr_link
+
+    if approx == "gaussian":
+        surprise = binary_surp_gauss_approx(poss_links,
+                                            poss_int_links,
+                                            obs_links,
+                                            int_links)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    elif approx == "asymptotic":
+        surprise = asymptot_surp_cd_bin_sum(l_o=int_links,
+                                            V=poss_links,
+                                            L=obs_links,
+                                            V_o=poss_int_links)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    else:
+        surprise = surprise_logsum_clust_bin(
+            poss_links,
+            int_links,
+            poss_int_links,
+            obs_links)
+
+    return surprise, mem_intr_link
+
+
+def asymptot_surp_cd_bin(l_o, V, L, V_o):
+    f = lambda x: integrand_asympt_cd_b(x, V, L, V_o)
+    aux = quad(f, [l_o, L], method="gauss-legendre")
+    return aux
+
+
+def asymptot_surp_cd_bin_sum(l_o, V, L, V_o):
+    f = lambda x: integrand_asympt_cd_b(x, V, L, V_o)
+    aux = 0
+    for l_o_loop in range(l_o, L):
+        aux += f(l_o_loop)
+    return aux
+
+
+def integrand_asympt_cd_b(l_o, V, L, V_o):
+    p = mpf(L / V)
+    p_d = mpf(l_o / V_o)
+    p_c = mpf((L - l_o) / (V - V_o))
+    a = a_l_o(l_o, V, L, V_o, p, p_d, p_c)
+    aux = (a * bernoulli(V, L, p)) / (
+                bernoulli(V_o, l_o, p_d) * bernoulli(V - V_o, L - l_o, p_c))
+    return aux
+
+
+def bernoulli(x, y, z):
+    aux = (mpf(z) ** y) * (mpf(1 - z) ** (x - y))
+    return aux
+
+
+def a_l_o(l_o, V, L, V_o, p, p_d, p_c):
+    aux1 = mpf(V*p*(1-p))**2
+    aux2 = mpf(2) * pi * (mpf(V_o * p_d * (1-p_d))**2) * (mpf((V-V_o) * p_c * (1-p_c))**2)
+    if aux2:
+        aux = sqrt(aux1/aux2)
+    else:
+        aux = 1
+    return aux
+
+
+def binary_surp_gauss_approx(V, V0, L, l0):
+    pi = V0 / V
+    P = ncdf(-(l0 - 1 - L * pi) / (np.sqrt(L * pi * (1 - pi))))
+    return P
+
+
+@jit(nopython=True)
+def surprise_logsum_clust_bin(F, p, M, m):
+    """[summary]
+
+    :param F: [description]
+    :type F: [type]
+    :param p: [description]
+    :type p: [type]
+    :param M: [description]
+    :type M: [type]
+    :param m: [description]
+    :type m: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    # stop = False
+    min_p = min(M, m)
+
+    logP = loghyperprobability(F, p, M, m)
+    for p_loop in range(p, min_p + 1):
+        if p_loop == p:
+            continue
+        nextLogP = loghyperprobability(F, p_loop, M, m)
+        [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+        if stop:
+            break
+
+    return -logP
+
+
+@jit(nopython=True)
+def loghyperprobability(F, p, M, m):
+    """Evaluates logarithmic hypergeometric distribution
+
+    :param F:
+    :type F:
+    :param p:
+    :type p:
+    :param M:
+    :type M:
+    :param m:
+    :type m:
+    :return:
+    :rtype:
+    """
+    logH = ax.logc(M, p) + ax.logc(F - M, m - p) - ax.logc(F, m)
+    return logH
+
+
+def calculate_surprise_logsum_clust_weigh_new(
+        adjacency_matrix,
+        cluster_assignment,
+        mem_intr_link,
+        clust_labels,
+        args,
+        approx,
+        is_directed):
+    """Calculates the logarithm of the surprise given the current partitions
+     for a weighted network. New faster implementation.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param mem_intr_link:
+    :type mem_intr_link:
+    :param clust_labels:
+    :type clust_labels:
+    :param args:
+    :type args:
+    :param approx:
+    :type approx:
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        if len(clust_labels):
+            w = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, w):
+                mem_intr_link[1][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[1])
+        intr_weights = p
+        # intracluster possible links
+        poss_intr_links = calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = args[1]
+        # Possible links
+        poss_links = args[2]
+        # extracluster links
+        inter_links = poss_links - poss_intr_links
+    else:
+        # intracluster weights
+        if len(clust_labels):
+            w = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, w):
+                mem_intr_link[1][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[1])
+        intr_weights = p / 2
+        # intracluster possible links
+        poss_intr_links = calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = args[1] / 2
+        # Possible links
+        poss_links = args[2] / 2
+        # extracluster links
+        inter_links = poss_links - poss_intr_links
+
+    if intr_weights == 0:
+        return 0, mem_intr_link
+
+    if approx == "gaussian":
+        surprise = weighted_suprise_approx(
+            poss_links,
+            tot_weights,
+            intr_weights,
+            poss_intr_links)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    elif approx == "asymptotic":
+        surprise = asymptot_surp_cd_wei_sum(poss_links,
+                                            tot_weights,
+                                            poss_intr_links,
+                                            intr_weights)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    else:
+        surprise = surprise_logsum_clust_weigh(
+            poss_intr_links,
+            intr_weights,
+            inter_links,
+            tot_weights,
+            poss_links)
+    return surprise, mem_intr_link
+
+
+def asymptot_surp_cd_wei(V, W, V_o, w_o):
+    f = lambda x: integrand_asympt_cd_w(x, V, W, V_o)
+    aux = quad(f, [w_o, W], method="gauss-legendre")
+    return aux
+
+
+@jit(forceobj=True)
+def asymptot_surp_cd_wei_sum(V, W, V_o, w_o):
+    aux_surp = 0
+    for w_o_loop in range(int(w_o), int(W)):
+        aux = integrand_asympt_cd_w(w_o_loop, V, W, V_o)
+        aux_surp += aux
+        if aux_surp == 0:
+            break
+        if aux / aux_surp < 1e-3:
+            break
+    return aux_surp
+
+
+def integrand_asympt_cd_w(w_o, V, W, V_o):
+    q = mpf(W/(V + W -1))
+    q_d = mpf(w_o/(V_o + w_o - 1))
+    q_c = mpf((W - w_o)/(V - V_o + W - w_o - 1))
+    C = C_w_o(V, W, V_o, w_o, q, q_d, q_c)
+    aux = (C * geometric(V, W, q))/(geometric(V_o, w_o, q_d)*geometric(V-V_o, W-w_o, q_c))
+    return aux
+
+
+def geometric(x, y, z):
+    aux = (mpf(z)**y) * (mpf(1 - z)**x)
+    return aux
+
+
+def C_w_o(V, W, V_o, w_o, q, q_d, q_c):
+    aux = sqrt(mpf(V*q)**2 / (mpf(2)*pi*(mpf(V_o*q_d)**2)*(mpf((V-V_o)*q_c)**2)))
+    return aux
+
+
+def weighted_suprise_approx(V, W, w_o, V_o):
+    """Gaussian approximation of surprise.
+    """
+    rho = W / V
+    aux = ncdf(-(w_o - 1 - V_o * rho) / (np.sqrt(V_o * rho * (1 + rho))))
+    return aux
+
+
+@jit(nopython=True)
+def surprise_logsum_clust_weigh(Vi, w, Ve, W, V):
+    """
+
+    :param Vi: [description]
+    :type Vi: [type]
+    :param w: [description]
+    :type w: [type]
+    :param Ve: [description]
+    :type Ve: [type]
+    :param W: [description]
+    :type W: [type]
+    :param V: [description]
+    :type V: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    # stop = False
+
+    logP = lognegativehyperprobability(Vi, w, Ve, W, V)
+    for w_loop in range(w, W):
+        if w_loop == w:
+            continue
+        nextLogP = lognegativehyperprobability(Vi, w_loop, Ve, W, V)
+        [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+        if stop:
+            break
+
+    return -logP
+
+
+@jit(nopython=True)
+def lognegativehyperprobability(Vi, w, Ve, W, V):
+    """Evaluates logarithmic hypergeometric distribution
+
+    :param Vi:
+    :type Vi:
+    :param w:
+    :type w:
+    :param Ve:
+    :type Ve:
+    :param W:
+    :type W:
+    :param V:
+    :type V:
+    :return:
+    :rtype:
+    """
+    logh = ax.logc(
+        Vi + w - 1,
+        w) + ax.logc(
+        Ve + W - w,
+        W - w) - ax.logc(
+        V + W,
+        W)
+    return logh
+
+
+@jit(nopython=True)
+def intracluster_links_enh_new(adj, clust_labels, partitions):
+    """Computes intracluster links and weights for enhanced community
+     detection method.
+
+    :param adj: Adjacency matrix.
+    :type adj: numpy.array
+    :param clust_labels:
+    :type clust_labels:
+    :param partitions: Nodes memberships.
+    :type partitions: numpy.array
+    :return: Number of intra-cluster links/weights.
+    :rtype: float
+    """
+    nr_intr_clust_links = np.zeros(clust_labels.shape[0])
+    intr_weight = np.zeros(clust_labels.shape[0])
+    for ii, lab in enumerate(clust_labels):
+        indices = np.where(partitions == lab)[0]
+        aux_l, aux_w = intracluster_links_aux_enh(adj, indices)
+        nr_intr_clust_links[ii] = aux_l
+        intr_weight[ii] = aux_w
+    return nr_intr_clust_links, intr_weight
+
+
+@jit(nopython=True)
+def intracluster_links_aux_enh(adj, indices):
+    """Computes intra-cluster links or weights given nodes indices.
+
+    :param adj: [description]
+    :type adj: [type]
+    :param indices: [description]
+    :type indices: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    weight = 0.0
+    n_links = 0.0
+    for ii in indices:
+        for jj in indices:
+            if adj[ii, jj]:
+                weight += adj[ii, jj]
+                n_links += 1
+    return n_links, weight
+
+
+def calculate_surprise_logsum_clust_enhanced_new(
+        adjacency_matrix,
+        cluster_assignment,
+        mem_intr_link,
+        clust_labels,
+        args,
+        approx,
+        is_directed):
+    """Calculates, for a weighted network, the logarithm of the enhanced
+     surprise given the current partitioning.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param mem_intr_link:
+    :type mem_intr_link:
+    :param clust_labels:
+    :type clust_labels:
+    :param args:
+    :type args:
+    :param approx:
+    :type approx:
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        if len(clust_labels):
+            l_aux, w_aux = intracluster_links_enh_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links, w_int in zip(clust_labels, l_aux, w_aux):
+                mem_intr_link[0][node_label] = nr_links
+                mem_intr_link[1][node_label] = w_int
+
+        l_o = int(mem_intr_link[0].sum())
+        w_o = int(mem_intr_link[1].sum())
+
+        # intracluster possible links
+        V_o = int(calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed))
+        # Total Weight
+        W = int(args[1])
+        L = int(args[0])
+        # Possible links
+        # n = adjacency_matrix.shape[0]
+        V = int(args[2])
+        # extracluster links
+        # inter_links = V - V_o
+    else:
+        # intracluster weights
+        if len(clust_labels):
+            l_aux, w_aux = intracluster_links_enh_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links, w_int in zip(clust_labels, l_aux, w_aux):
+                mem_intr_link[0][node_label] = nr_links
+                mem_intr_link[1][node_label] = w_int
+
+        l_o = int(mem_intr_link[0].sum() / 2)
+        w_o = int(mem_intr_link[1].sum() / 2)
+
+        # intracluster possible links
+        V_o = int(calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed))
+        # Total Weight
+        W = int(args[1] / 2)
+        L = int(args[0] / 2)
+        # Possible links
+        # n = adjacency_matrix.shape[0]
+        V = int(args[2] / 2)
+        # extracluster links
+        # inter_links = V - V_o
+
+    if l_o == 0:
+        return 0, mem_intr_link
+
+    # print("V_0", V_o, "l_0", l_o, "w_0", w_o, "V", V, "L", L, "W", W)
+
+    surprise = surprise_logsum_clust_enh(V_o, l_o, w_o, V, L, W)
+    return surprise, mem_intr_link
+
+
+@jit(nopython=True)
+def surprise_logsum_clust_enh(V_o, l_o, w_o, V, L, W):
+    # stop = False
+    # stop1 = False
+    min_l_loop = min(L, V_o)
+
+    logP = logenhancedhypergeometric(V_o, l_o, w_o, V, L, W)
+    logP1 = logP
+    w_loop = w_o
+
+    for l_loop in range(l_o, min_l_loop + 1):
+        if l_loop == 0:
+            continue
+        for w_loop in range(w_o - l_loop + l_o, W - L + l_o + 1):
+            if w_loop <= 0:
+                continue
+            if (w_loop == w_o) and (l_loop == l_o):
+                continue
+            nextLogP = logenhancedhypergeometric(V_o, l_loop, w_loop, V, L, W)
+            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+            if stop:
+                break
+        nextLogP1 = logenhancedhypergeometric(V_o, l_loop, w_loop, V, L, W)
+        [logP1, stop1] = ax.sumLogProbabilities(nextLogP1, logP1)
+        if stop1:
+            break
+
+    return -logP1
+
+
+@jit(nopython=True)
+def logenhancedhypergeometric(V_o, l_o, w_o, V, L, W):
+    if l_o < L:
+        aux1 = (ax.logc(V_o, l_o) + ax.logc(V - V_o, L - l_o)) - ax.logc(V, L)
+        aux2 = (ax.logc(w_o - 1, w_o - l_o) + ax.logc(W - w_o - 1, (W - L) - (
+                    w_o - l_o))) - ax.logc(W - 1, W - L)
+    else:
+        aux1 = (ax.logc(V_o, L) - ax.logc(V, L))
+        aux2 = (ax.logc(w_o, w_o - L) - ax.logc(W, L))
+    return aux1 + aux2
+
+
+def calculate_surprise_logsum_clust_weigh_continuos(
+        adjacency_matrix,
+        cluster_assignment,
+        mem_intr_link,
+        clust_labels,
+        args,
+        approx,
+        is_directed):
+    """Calculates the logarithm of the continuos surprise given the current
+     partitions for a weighted network. New faster implementation.
+
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Nodes memberships.
+    :type cluster_assignment: numpy.ndarray
+    :param mem_intr_link:
+    :type mem_intr_link:
+    :param clust_labels:
+    :type clust_labels:
+    :param args:
+    :type args:
+    :param approx:
+    :type approx:
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise.
+    :rtype: float
+    """
+    if is_directed:
+        # intracluster weights
+        if len(clust_labels):
+            w = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, w):
+                mem_intr_link[1][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[1])
+        intr_weights = p
+        # intracluster possible links
+        poss_intr_links = calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = args[1]
+        # Possible links
+        poss_links = args[2]
+        # extracluster links
+        # inter_links = poss_links - poss_intr_links
+    else:
+        # intracluster weights
+        if len(clust_labels):
+            w = intracluster_links_new(
+                adj=adjacency_matrix,
+                clust_labels=clust_labels,
+                partitions=cluster_assignment)
+
+            for node_label, nr_links in zip(clust_labels, w):
+                mem_intr_link[1][node_label] = nr_links
+
+        p = np.sum(mem_intr_link[1])
+        intr_weights = p / 2
+        # intracluster possible links
+        poss_intr_links = calculate_possible_intracluster_links_new(
+            cluster_assignment,
+            is_directed)
+        # Total Weight
+        tot_weights = args[1] / 2
+        # Possible links
+        poss_links = args[2] / 2
+        # extracluster links
+        # inter_links = poss_links - poss_intr_links
+
+    if intr_weights == 0:
+        return 0, mem_intr_link
+
+    if approx == "gaussian":
+        surprise = weighted_suprise_approx(
+            poss_links,
+            tot_weights,
+            intr_weights,
+            poss_intr_links)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    elif approx == "asymptotic":
+        surprise = asymptot_surp_cd_wei(
+            poss_links,
+            tot_weights,
+            poss_intr_links,
+            intr_weights)
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+    else:
+        surprise = continuous_surprise_clust(
+            poss_links,
+            tot_weights,
+            intr_weights,
+            poss_intr_links)
+
+        if surprise > 0:
+            surprise = np.float64(-log10(surprise))
+        else:
+            surprise = 0
+
+    return surprise, mem_intr_link
+
+
+def continuous_surprise_clust(V, W, w_o, V_o):
+    f = lambda x: integrand_clust(x, V, W, V_o)
+    aux_surp = quad(f, [w_o, W], method="gauss-legendre")
+    return aux_surp
+
+
+def integrand_clust(w_o, V, W, V_o):
+    aux = W * beta(V, W) / (w_o*beta(V_o, w_o) * (W-w_o) * beta(V-V_o, W-w_o))
+    return aux
+
+
+def labeling_communities(partitions):
+    """Gives labels to communities from 0 to number of communities minus one.
+
+    :param partitions: Nodes memberships.
+    :type partitions: numpy.ndarray
+    :return: Re-labeled nodes memberships.
+    :rtype: numpy.ndarray
+    """
+    different_partitions = np.unique(partitions, return_counts=True)
+    aux_argsort = np.argsort(different_partitions[1])[::-1]
+    ordered_clusters = different_partitions[0][aux_argsort]
+    new_partitioning = partitions.copy()
+    new_partitions = np.array([k for k in range(len(ordered_clusters))])
+    for old_part, new_part in zip(ordered_clusters, new_partitions):
+        # print(old_part,new_part)
+        indices_old = np.where(partitions == old_part)[0]
+        new_partitioning[indices_old] = new_part
+
+    return new_partitioning
+
+
+def flipping_function_comdet_agl_new(
+        calculate_surprise,
+        adj,
+        membership,
+        mem_intr_link,
+        args,
+        surprise,
+        approx,
+        is_directed):
+    list_neigh = ax.compute_neighbours(adj)
+
+    for node, node_label in zip(np.arange(membership.shape[0]), membership):
+        for node2 in list_neigh[node]:
+            new_clust = membership[node2]
+            if node_label != new_clust:
+                aux_membership = membership.copy()
+                aux_membership[node] = new_clust
+                # print(np.array([node_label, new_clust]))
+                temp_surprise, temp_mem_intr_link = calculate_surprise(
+                    adjacency_matrix=adj,
+                    cluster_assignment=aux_membership,
+                    mem_intr_link=mem_intr_link.copy(),
+                    clust_labels=np.array([node_label, new_clust]),
+                    args=args,
+                    approx=approx,
+                    is_directed=is_directed)
+                if temp_surprise > surprise:
+                    membership = aux_membership.copy()
+                    surprise = temp_surprise
+                    mem_intr_link = temp_mem_intr_link
+
+    return membership
+
+
+def flipping_function_comdet_div_new(
+        calculate_surprise,
+        adj,
+        membership,
+        mem_intr_link,
+        args,
+        surprise,
+        approx,
+        is_directed):
+    list_neigh = ax.compute_neighbours(adj)
+
+    for node, node_label in zip(np.arange(membership.shape[0]), membership):
+        if len(np.where(membership == node_label)[0]) == 1:
+            continue
+        for node2 in list_neigh[node]:
+            new_clust = membership[node2]
+            if node_label != new_clust:
+                aux_membership = membership.copy()
+                aux_membership[node] = new_clust
+                # print(np.array([node_label, new_clust]))
+                temp_surprise, temp_mem_intr_link = calculate_surprise(
+                    adjacency_matrix=adj,
+                    cluster_assignment=aux_membership,
+                    mem_intr_link=mem_intr_link.copy(),
+                    clust_labels=np.array([node_label, new_clust]),
+                    args=args,
+                    approx=approx,
+                    is_directed=is_directed)
+                if temp_surprise > surprise:
+                    membership = aux_membership.copy()
+                    surprise = temp_surprise
+                    mem_intr_link = temp_mem_intr_link
+
+    return membership
diff --git a/src/surprisememore/cp_functions.py b/src/surprisememore/cp_functions.py
new file mode 100644
index 00000000..d10ff046
--- /dev/null
+++ b/src/surprisememore/cp_functions.py
@@ -0,0 +1,521 @@
+import numpy as np
+from numba import jit
+
+from . import auxiliary_function as ax
+from . import comdet_functions as cd
+
+
+@jit(nopython=True)
+def compute_sum(adj, row_indices, column_indices):
+    """Computes number of links/total weight given nodes indices.
+
+    :param adj: Adjacency matrix.
+    :type adj: numpy.ndarray
+    :param row_indices: Row indices.
+    :type row_indices: numpy.ndarray
+    :param column_indices: Columns indices.
+    :type column_indices: numpy.ndarray
+    :return: Total links/weights
+    :rtype: float
+    """
+    Sum = 0.0
+    for ii in row_indices:
+        for jj in column_indices:
+            Sum += adj[ii, jj]
+    return Sum
+
+
+def calculate_surprise_logsum_cp_weigh(adjacency_matrix,
+                                       cluster_assignment,
+                                       is_directed):
+    """Computes core-periphery weighted log-surprise given a certain nodes'
+     partitioning.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Core periphery assigments.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise
+    :rtype: float
+    """
+    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
+    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
+
+    if is_directed:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = n_c * (n_c - 1)
+        p_x = n_c * n_x * 2
+
+        w_c = compute_sum(adjacency_matrix, core_nodes, core_nodes)
+        w_x = compute_sum(adjacency_matrix,
+                          core_nodes,
+                          periphery_nodes) + compute_sum(adjacency_matrix,
+                                                         periphery_nodes,
+                                                         core_nodes)
+
+        w = np.sum(adjacency_matrix)
+        w_p = w - w_c - w_x
+        n = n_c + n_x
+        p = n * (n - 1)
+        p_p = p - p_c - p_x
+
+    else:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = n_c * (n_c - 1) / 2
+        p_x = n_c * n_x
+
+        w_c = (compute_sum(adjacency_matrix, core_nodes, core_nodes)) / 2
+        w_x = (compute_sum(adjacency_matrix,
+                           core_nodes,
+                           periphery_nodes) + compute_sum(adjacency_matrix,
+                                                          periphery_nodes,
+                                                          core_nodes)) / 2
+
+        w = np.sum(adjacency_matrix) / 2
+        w_p = (w - w_c - w_x) / 2
+        n = n_c + n_x
+        p = n * (n - 1) / 2
+        p_p = p - p_c - p_x
+
+    surprise = surprise_bipartite_logsum_CP_Weigh(p, p_c, p_x, p_p,
+                                                  w, w_c, w_x, w_p)
+    return surprise
+
+
+@jit(nopython=True)
+def surprise_bipartite_logsum_CP_Weigh(p, p_c, p_x, p_p, w, w_c, w_x, w_p):
+    stop = False
+    first_loop_break = False
+
+    logP = logMultiHyperProbabilityWeight(p, p_c, p_x, p_p, w, w_c, w_x, w_p)
+    for w_c_loop in range(w_c, w + 1):
+        for w_x_loop in range(w_x, w + 1 - w_c_loop):
+            if (w_c_loop == w_c) & (w_x_loop == w_x):
+                continue
+            w_p_loop = w - w_c_loop - w_x_loop
+            nextLogP = logMultiHyperProbabilityWeight(p, p_c, p_x,
+                                                      p_p, w, w_c_loop,
+                                                      w_x_loop, w_p_loop)
+            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+
+            if stop:
+                first_loop_break = True
+                break
+        if first_loop_break:
+            break
+    return -logP
+
+
+@jit(nopython=True)
+def logMultiHyperProbabilityWeight(p, p_c, p_x, p_p, w, w_c, w_x, w_p):
+    """Computes the logarithm of the Negative Multinomial Hypergeometric
+     distribution."""
+    logH = (ax.logc(p_c + w_c - 1, w_c) + ax.logc(p_x + w_x - 1, w_x)
+            + ax.logc(p_p + w_p - 1, w_p) - ax.logc(p + w, w))
+    return logH
+
+
+def calculate_surprise_logsum_cp_bin(adjacency_matrix,
+                                     cluster_assignment,
+                                     is_directed):
+    """Computes core-periphery binary log-surprise given a certain nodes'
+     partitioning.
+
+    :param adjacency_matrix: Binary adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Core periphery assigments.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise
+    :rtype: float
+    """
+    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
+    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
+
+    if is_directed:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = n_c * (n_c - 1)
+        p_x = n_c * n_x * 2
+
+        l_c = compute_sum(adjacency_matrix, core_nodes, core_nodes)
+        l_x = compute_sum(adjacency_matrix, core_nodes,
+                          periphery_nodes) + compute_sum(adjacency_matrix,
+                                                         periphery_nodes,
+                                                         core_nodes)
+
+        l_t = np.sum(adjacency_matrix)
+        n = n_c + n_x
+        p = n * (n - 1)
+
+    else:
+        n_c = core_nodes.shape[0]
+        n_x = periphery_nodes.shape[0]
+        p_c = (n_c * (n_c - 1)) / (2)
+        p_x = n_c * n_x
+
+        l_c = compute_sum(adjacency_matrix, core_nodes, core_nodes) / 2
+        l_x = (compute_sum(adjacency_matrix, core_nodes,
+                           periphery_nodes) + compute_sum(adjacency_matrix,
+                                                          periphery_nodes,
+                                                          core_nodes)) / 2
+
+        l_t = np.sum(adjacency_matrix) / 2
+        n = n_c + n_x
+        p = (n * (n - 1)) / 2
+
+    if (p_c + p_x) < (l_c + l_x):
+        return 0
+
+    surprise = surprise_bipartite_logsum_CP_Bin(p, p_c, p_x, l_t, l_c, l_x)
+    return surprise
+
+
+@jit(nopython=True)
+def surprise_bipartite_logsum_CP_Bin(p, p_c, p_x, l, l_c, l_x):
+    stop = False
+    first_loop_break = False
+
+    # min_l_p = min(l, p_c + p_x)
+
+    logP = logMultiHyperProbability(p, p_c, p_x, l, l_c, l_x)
+    for l_c_loop in range(l_c, p_c + 1):
+        for l_x_loop in range(l_x, p_x + 1):
+            if ((l_c_loop == l_c) & (l_x_loop == l_x) or
+                    (l_c_loop + l_x_loop > l)):
+                continue
+            nextLogP = logMultiHyperProbability(p, p_c, p_x,
+                                                l, l_c_loop,
+                                                l_x_loop)
+            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+
+            if stop:
+                first_loop_break = True
+                break
+        if first_loop_break:
+            break
+
+    return -logP
+
+
+@jit(nopython=True)
+def logMultiHyperProbability(p, p_c, p_x, l, l_c, l_x):
+    """Computes the logarithm of the Multinomial Hypergeometric
+     distribution."""
+    logH = ax.logc(p_c, l_c) + ax.logc(p_x, l_x) + ax.logc(
+        p - p_c - p_x,
+        l - l_c - l_x) - ax.logc(p, l)
+    return logH
+
+
+@jit(nopython=True)
+def compute_sum_enh(adj, row_indices, column_indices):
+    """Computes number number of links and weights given nodes indices.
+
+    :param adj: Adjacency matrix.
+    :type adj: numpy.ndarray
+    :param row_indices: Row indices.
+    :type row_indices: numpy.ndarray
+    :param column_indices: Columns indices.
+    :type column_indices: numpy.ndarray
+    :return: Total links/weights
+    :rtype: float
+    """
+    n_link = 0.0
+    weigth = 0.0
+    for ii in row_indices:
+        for jj in column_indices:
+            if adj[ii, jj]:
+                weigth += adj[ii, jj]
+                n_link += 1
+    return n_link, weigth
+
+
+def calculate_surprise_logsum_cp_enhanced(adjacency_matrix,
+                                          cluster_assignment,
+                                          is_directed):
+    """Computes core-periphery enhanced log-surprise given
+     a certain nodes partitioning.
+
+    :param adjacency_matrix: Weighted adjacency matrix.
+    :type adjacency_matrix: numpy.ndarray
+    :param cluster_assignment: Core periphery assigments.
+    :type cluster_assignment: numpy.ndarray
+    :param is_directed: True if the graph is directed.
+    :type is_directed: bool
+    :return: Log-surprise
+    :rtype: float
+    """
+    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
+    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
+
+    if is_directed:
+        n_o = core_nodes.shape[0]
+        n_p = periphery_nodes.shape[0]
+        V_o = int(n_o * (n_o - 1))
+        V_c = int(n_o * n_p * 2)
+
+        l_o, w_o = compute_sum_enh(adjacency_matrix, core_nodes, core_nodes)
+        l_c1, w_c1 = compute_sum_enh(adjacency_matrix,
+                                     core_nodes,
+                                     periphery_nodes)
+        l_c2, w_c2 = compute_sum_enh(adjacency_matrix,
+                                     periphery_nodes,
+                                     core_nodes)
+
+        l_o = int(l_o)
+        w_o = int(w_o)
+        l_c = int(l_c1 + l_c2)
+        w_c = int(w_c1 + w_c2)
+        L = int(np.sum(adjacency_matrix.astype(bool)))
+        W = int(np.sum(adjacency_matrix))
+        # w_p = W - w_o - w_c
+        n = n_o + n_p
+        V = int(n * (n - 1))
+
+    else:
+        n_o = core_nodes.shape[0]
+        n_p = periphery_nodes.shape[0]
+        V_o = int(n_o * (n_o - 1) / 2)
+        V_c = int(n_o * n_p)
+
+        l_o, w_o = (compute_sum_enh(adjacency_matrix, core_nodes, core_nodes))
+        l_c1, w_c1 = compute_sum_enh(adjacency_matrix,
+                                     core_nodes,
+                                     periphery_nodes)
+        l_c2, w_c2 = compute_sum_enh(adjacency_matrix,
+                                     periphery_nodes,
+                                     core_nodes)
+
+        l_o = int(l_o / 2)
+        w_o = int(w_o / 2)
+        l_c = int((l_c1 + l_c2) / 2)
+        w_c = int((w_c1 + w_c2) / 2)
+
+        L = int(np.sum(adjacency_matrix.astype(bool)) / 2)
+        W = int(np.sum(adjacency_matrix) / 2)
+        # w_p = (W - w_o - w_c) / 2
+        n = n_o + n_p
+        V = int(n * (n - 1) / 2)
+
+    # print("V_o", V_o, "l_o", l_o, "V_c", V_c, "l_c", l_c,
+    #      "w_o", w_o, "w_c", w_c, "V", V, "L", L, "W", W)
+
+    surprise = surprise_bipartite_logsum_CP_enh(V_o, l_o, V_c, l_c,
+                                                w_o, w_c, V, L, W)
+    return surprise
+
+
+def surprise_bipartite_logsum_CP_enh(V_o, l_o, V_c, l_c, w_o, w_c, V, L, W):
+    skip = False
+
+    if (l_o > 0) and (l_o < L) and (l_c > 0) and (l_c < L):
+        logP = log_enh_multy_hyper(V_o, l_o, V_c, l_c, w_o, w_c, V, L, W)
+    elif (l_o == 0) and (l_c == L):
+        logP = log_enh_red_univ_hyper(V_c, w_c,
+                                      V, L, W)
+        skip = True
+
+    elif (l_o == L) and (l_c == 0):
+        logP = log_enh_red_univ_hyper(V_o, w_o,
+                                      V, L, W)
+        skip = True
+    elif (l_o == 0) and (l_c > 0) and (l_c < L):
+        logP = cd.logenhancedhypergeometric(V_o=V_c, l_o=l_c,
+                                            w_o=w_c, V=V,
+                                            L=L, W=W)
+
+    elif (l_o > 0) and (l_o < L) and (l_c == 0):
+        logP = cd.logenhancedhypergeometric(V_o=V_o, l_o=l_o,
+                                            w_o=w_o, V=V,
+                                            L=L, W=W)
+    else:
+        logP = log_enh_red_univ_hyper(V - V_o - V_c,
+                                      W - w_o - w_c,
+                                      V, L, W)
+
+    logP1 = logP
+    logP2 = logP
+    logP3 = logP
+    l_c_loop = l_c
+    w_o_loop = w_o
+    w_c_loop = w_c
+
+    for l_o_loop in range(l_o, V_o + 1):
+        for l_c_loop in range(l_c, V_c + 1):
+            if (l_o_loop + l_c_loop > L) or (l_o_loop + l_c_loop == 0):
+                continue
+
+            if ((l_o_loop > 0) and (l_o_loop < L) and
+                    (l_c_loop > 0) and (l_c_loop < L)):
+                w_o_loop, w_c_loop, logP, logP1 = case_one(
+                    logP, logP1,
+                    l_o_loop, l_c_loop,
+                    V_o, l_o, V_c,
+                    l_c, w_o, w_c,
+                    V, L, W)
+                nextLogP2 = log_enh_multy_hyper(V_o, l_o_loop,
+                                                  V_c, l_c_loop,
+                                                  w_o_loop, w_c_loop,
+                                                  V, L, W)
+                [logP2, stop2] = ax.sumLogProbabilities(nextLogP2, logP2)
+                if stop2:
+                    break
+            elif ((l_o_loop == 0) and (l_c_loop > 0) and
+                  (l_c_loop < L)):
+                w_loop, logP = case_two(logP, l_c_loop, V_c,
+                                        l_c, w_c, V, L, W)
+                logP1 = logP
+                nextLogP2 = cd.logenhancedhypergeometric(V_c, l_c_loop, w_loop,
+                                                      V, L, W)
+
+                [logP2, stop2] = ax.sumLogProbabilities(nextLogP2, logP2)
+                if stop2:
+                    break
+            elif ((l_o_loop > 0) and (l_o_loop < L) and
+                  (l_c_loop == 0)):
+
+                continue
+
+        if ((l_o_loop > 0) and (l_o_loop < L) and
+                (l_c_loop > 0) and (l_c_loop < L)):
+
+            nextLogP3 = log_enh_multy_hyper(V_o, l_o_loop,
+                                            V_c, l_c_loop,
+                                            w_o_loop, w_c_loop,
+                                            V, L, W)
+            [logP3, stop3] = ax.sumLogProbabilities(nextLogP3, logP3)
+            if stop3:
+                break
+        elif ((l_o_loop == 0) and (l_c_loop > 0) and
+              (l_c_loop < L)):
+            continue
+        elif ((l_o_loop > 0) and (l_o_loop < L) and
+              (l_c_loop == 0)):
+            w_loop, logP = case_two(logP, l_o_loop, V_o,
+                                    l_o, w_o, V, L, W)
+            logP1 = logP
+            nextLogP3 = cd.logenhancedhypergeometric(V_o, l_o_loop,
+                                                     w_loop, V, L, W)
+
+            [logP3, stop3] = ax.sumLogProbabilities(nextLogP3, logP3)
+            if stop3:
+                break
+
+    if skip:
+        logP3 = logP
+
+    return -logP3
+
+
+@jit(nopython=True)
+def case_one(logP, logP1, l_o_loop, l_c_loop, V_o,
+             l_o, V_c, l_c, w_o, w_c, V, L, W):
+    w_c_loop = w_c
+    w_o_loop = w_o
+    for w_o_loop in range(w_o - l_o_loop + l_o, W - L + l_o + 1):
+        for w_c_loop in range(w_c - l_c_loop + l_c,
+                              W + L + l_c + 1 - w_o_loop):
+            if ((w_o_loop == w_o) & (w_c_loop == w_c) &
+                    (l_o_loop == l_o) & (l_c_loop == l_c)):
+                continue
+
+            nextLogP = log_enh_multy_hyper(V_o, l_o_loop,
+                                           V_c, l_c_loop,
+                                           w_o_loop, w_c_loop,
+                                           V, L, W)
+            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+            if stop:
+                break
+        nextLogP1 = log_enh_multy_hyper(V_o, l_o_loop,
+                                        V_c, l_c_loop,
+                                        w_o_loop, w_c_loop,
+                                        V, L, W)
+        [logP1, stop1] = ax.sumLogProbabilities(nextLogP1, logP1)
+        if stop1:
+            break
+    return w_o_loop, w_c_loop, logP, logP1
+
+
+@jit(nopython=True)
+def case_two(logP, l_loop, V_o,
+             l_o, w_o, V, L, W):
+    w_loop = w_o
+    for w_loop in range(w_o - l_loop + l_o, W - L + l_o + 1):
+        if (w_loop <= 0) & ((w_loop == w_o) & (l_loop == l_o)):
+            continue
+        nextLogP = cd.logenhancedhypergeometric(V_o, l_loop, w_loop, V, L, W)
+        [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
+        if stop:
+            break
+    return w_loop, logP
+
+
+@jit(nopython=True)
+def log_enh_multy_hyper(V_o, l_o, V_c, l_c, w_o, w_c, V, L, W):
+    """Computes the logarithm of the Negative Multinomial
+     Hypergeometric distribution."""
+    aux1 = (ax.logc(V_o, l_o) + ax.logc(V_c, l_c) + ax.logc(V - (V_o + V_c),
+                                                            L - (l_o + l_c))) - ax.logc(V, L)
+    aux2 = (ax.logc(w_o - 1, l_o - 1) + ax.logc(w_c - 1, l_c - 1) + ax.logc(
+        W - (w_o + w_c) - 1, L - (l_o + l_c) - 1)) - ax.logc(W - 1, W - L)
+    return aux1 + aux2
+
+
+@jit(nopython=True)
+def log_enh_red_univ_hyper(V_o, w_o, V, L, W):
+    aux1 = ax.logc(V_o, L) - ax.logc(V, L)
+    aux2 = ax.logc(w_o - 1, w_o - L)
+    return aux1 + aux2
+
+
+def labeling_core_periphery(adjacency_matrix, cluster_assignment):
+    """Function assigning the core and periphery labels based on link density
+
+    :param adjacency_matrix: [description]
+    :type adjacency_matrix: [type]
+    :param cluster_assignment: [description]
+    :type cluster_assignment: [type]
+    :return: [description]
+    :rtype: [type]
+    """
+    core_nodes = np.where(cluster_assignment == 0)[0]
+    periphery_nodes = np.where(cluster_assignment == 1)[0]
+    l_core = np.sum(
+        adjacency_matrix[np.ix_(list(core_nodes), list(core_nodes))] > 0)
+    l_periphery = np.sum(adjacency_matrix[np.ix_(list(periphery_nodes),
+                                                 list(periphery_nodes))] > 0)
+    core_density = l_core / (len(core_nodes) * (len(core_nodes) - 1))
+    periphery_density = l_periphery / (
+            len(periphery_nodes) * (len(periphery_nodes) - 1))
+
+    if periphery_density > core_density:
+        cluster_assignment_new = 1 - cluster_assignment
+        return cluster_assignment_new
+    return cluster_assignment
+
+
+def flipping_function_cp(comm, n_flipping):
+    """Moves n nodes, randomly selected, from a partition to the other.
+
+    :param comm: Nodes memberships.
+    :type comm: numpy.ndarray
+    :param n_flipping: Number of nodes to flipped.
+    :type n_flipping: int
+    :return: New nodes memberships.
+    :rtype: numpy.ndarray
+    """
+    if np.random.random() > 0.5:
+        comm[np.random.choice(np.where(comm == 1)[0],
+                              replace=False,
+                              size=n_flipping)] = 0
+    else:
+        comm[np.random.choice(np.where(comm == 0)[0],
+                              replace=False,
+                              size=n_flipping)] = 1
+    return comm
diff --git a/src/surprisemes/solver.py b/src/surprisememore/solver.py
similarity index 74%
rename from src/surprisemes/solver.py
rename to src/surprisememore/solver.py
index e1d82c63..6e10c50f 100644
--- a/src/surprisemes/solver.py
+++ b/src/surprisememore/solver.py
@@ -1,4 +1,5 @@
 import numpy as np
+from tqdm import tqdm
 
 from . import comdet_functions as cd
 
@@ -37,7 +38,7 @@ def solver_cp(adjacency_matrix,
     sim = 0
     while sim < num_sim:
         # edges_counter = 0
-        for [u, v] in edges_sorted:
+        for [u, v] in tqdm(edges_sorted):
             # surprise_old = surprise
             cluster_assignment_temp1 = cluster_assignment.copy()
             cluster_assignment_temp2 = cluster_assignment.copy()
@@ -108,6 +109,7 @@ def solver_com_det_aglom(
         correct_partition_labeling,
         prob_mix,
         flipping_function,
+        approx,
         is_directed,
         print_output=False):
     """Community detection solver. It carries out the research of the optimal
@@ -129,6 +131,8 @@ def solver_com_det_aglom(
     :type prob_mix: [type]
     :param flipping_function:
     :type flipping_function:
+    :param approx:
+    :type approx:
     :param is_directed:
     :type is_directed:
     :param print_output: [description], defaults to False
@@ -138,17 +142,24 @@ def solver_com_det_aglom(
     """
     prob_random = (1 - prob_mix) / 2
 
-    obs_links = int(np.sum(adjacency_matrix))
+    obs_links = np.sum(adjacency_matrix.astype(bool))
+    obs_weights = np.sum(adjacency_matrix)
     n_nodes = int(adjacency_matrix.shape[0])
-    poss_links = int(n_nodes * (n_nodes - 1))
-    args = (obs_links, poss_links)
+    poss_links = n_nodes * (n_nodes - 1)
+    args = (obs_links, obs_weights, poss_links)
 
-    mem_intr_link = np.zeros(cluster_assignment.shape[0], dtype=np.int32)
+    cluster_assignment = correct_partition_labeling(
+        cluster_assignment.copy())
+    n_clusters = np.unique(cluster_assignment).shape[0]
+
+    mem_intr_link = np.zeros((2, n_clusters), dtype=np.float64)
     for ii in np.unique(cluster_assignment):
         indices = np.where(cluster_assignment == ii)[0]
-        mem_intr_link[ii] = cd.intracluster_links_aux(
+        l_aux, w_aux = cd.intracluster_links_aux_enh(
             adjacency_matrix,
             indices)
+        mem_intr_link[0][ii] = l_aux
+        mem_intr_link[1][ii] = w_aux
 
     surprise, _ = calculate_surprise(
         adjacency_matrix,
@@ -156,6 +167,7 @@ def solver_com_det_aglom(
         mem_intr_link,
         np.array([]),
         args,
+        approx,
         is_directed)
 
     contatore_break = 0
@@ -166,7 +178,7 @@ def solver_com_det_aglom(
         e_sorted = sort_edges(adjacency_matrix)
         # print(cluster_assignment)
         # print(e_sorted)
-        for [u, v] in e_sorted:
+        for [u, v] in tqdm(e_sorted):
             if cluster_assignment[u] != cluster_assignment[v]:
                 clus_u = cluster_assignment[u]
                 clus_v = cluster_assignment[v]
@@ -181,6 +193,7 @@ def solver_com_det_aglom(
                         mem_intr_link.copy(),
                         np.array([clus_u, clus_v]),
                         args,
+                        approx,
                         is_directed)
                     # print(cluster_assignement_temp, cluster_assignment)
                     # print("prob_1", u, v, temp_surprise, surprise)
@@ -193,6 +206,7 @@ def solver_com_det_aglom(
                         mem_intr_link.copy(),
                         np.array([clus_u, clus_v]),
                         args,
+                        approx,
                         is_directed)
                     # print(cluster_assignement_temp, cluster_assignment)
                     # print("prob_2", u, v, temp_surprise, surprise)
@@ -207,6 +221,7 @@ def solver_com_det_aglom(
                         mem_intr_link.copy(),
                         np.array([clus_u, clus_v]),
                         args,
+                        approx,
                         is_directed)
                     # print(cluster_assignement_temp, cluster_assignment)
                     # print("mixing", u, v, temp_surprise, surprise)
@@ -214,7 +229,7 @@ def solver_com_det_aglom(
                 if temp_surprise > surprise:
                     cluster_assignment = cluster_assignement_temp.copy()
                     surprise = temp_surprise
-                    mem_intr_link = temp_mem_intr_link
+                    mem_intr_link = temp_mem_intr_link.copy()
 
         if surprise > previous_surprise:
             contatore_break = 0
@@ -228,8 +243,14 @@ def solver_com_det_aglom(
         sim += 1
 
     cluster_assignment = flipping_function(
-        adjacency_matrix,
-        cluster_assignment.copy())
+        calculate_surprise=calculate_surprise,
+        adj=adjacency_matrix,
+        membership=cluster_assignment.copy(),
+        mem_intr_link=mem_intr_link.copy(),
+        args=args,
+        surprise=surprise,
+        approx=approx,
+        is_directed=is_directed)
 
     cluster_assignement_proper = correct_partition_labeling(
         cluster_assignment.copy())
@@ -244,6 +265,7 @@ def solver_com_det_divis(
         calculate_surprise,
         correct_partition_labeling,
         flipping_function,
+        approx,
         is_directed,
         print_output=False):
     """Community detection solver. It carries out the research of the optimal
@@ -263,6 +285,8 @@ def solver_com_det_divis(
     :type correct_partition_labeling: [type]
     :param flipping_function:
     :type flipping_function:
+    :param approx:
+    :type approx:
     :param is_directed:
     :type is_directed:
     :param print_output: [description], defaults to False
@@ -270,18 +294,22 @@ def solver_com_det_divis(
     :return: [description]
     :rtype: [type]
     """
-    obs_links = int(np.sum(adjacency_matrix))
+    obs_links = np.sum(adjacency_matrix.astype(bool))
+    obs_weights = np.sum(adjacency_matrix)
     n_nodes = int(adjacency_matrix.shape[0])
-    poss_links = int(n_nodes * (n_nodes - 1))
-    args = (obs_links, poss_links)
+    poss_links = n_nodes * (n_nodes - 1)
+    args = (obs_links, obs_weights, poss_links)
 
     n_clusters = np.unique(cluster_assignment).shape[0]
-    mem_intr_link = np.zeros(n_clusters, dtype=np.int32)
+
+    mem_intr_link = np.zeros((2, n_clusters), dtype=np.float64)
     for ii in np.unique(cluster_assignment):
         indices = np.where(cluster_assignment == ii)[0]
-        mem_intr_link[ii] = cd.intracluster_links_aux(
+        l_aux, w_aux = cd.intracluster_links_aux_enh(
             adjacency_matrix,
             indices)
+        mem_intr_link[0][ii] = l_aux
+        mem_intr_link[1][ii] = w_aux
 
     surprise, _ = calculate_surprise(
         adjacency_matrix,
@@ -289,6 +317,7 @@ def solver_com_det_divis(
         mem_intr_link,
         np.array([]),
         args,
+        approx,
         is_directed)
 
     contatore_break = 0
@@ -299,7 +328,7 @@ def solver_com_det_divis(
         e_sorted = sort_edges(adjacency_matrix)
         # print(cluster_assignment)
         # print(E_sorted)
-        for [u, v] in e_sorted:
+        for [u, v] in tqdm(e_sorted):
             cluster_assignment_temp1 = cluster_assignment.copy()
             cluster_assignment_temp2 = cluster_assignment.copy()
 
@@ -315,9 +344,11 @@ def solver_com_det_divis(
                     mem_intr_link.copy(),
                     np.array([clus_u, clus_v]),
                     args,
+                    approx,
                     is_directed)
 
-                if surprise_temp1 > surprise:
+                aux_n_clus = np.unique(cluster_assignment_temp1).shape[0]
+                if (surprise_temp1 > surprise) and (n_clusters == aux_n_clus):
                     cluster_assignment = cluster_assignment_temp1.copy()
                     surprise = surprise_temp1
                     mem_intr_link = temp_mem_intr_link1
@@ -328,9 +359,11 @@ def solver_com_det_divis(
                     mem_intr_link.copy(),
                     np.array([clus_u, clus_v]),
                     args,
+                    approx,
                     is_directed)
 
-                if surprise_temp2 > surprise:
+                aux_n_clus = np.unique(cluster_assignment_temp2).shape[0]
+                if (surprise_temp2 > surprise) and (n_clusters == aux_n_clus):
                     cluster_assignment = cluster_assignment_temp2.copy()
                     surprise = surprise_temp2
                     mem_intr_link = temp_mem_intr_link2
@@ -349,9 +382,11 @@ def solver_com_det_divis(
                     mem_intr_link.copy(),
                     np.array([clus_u, clus_v]),
                     args,
+                    approx,
                     is_directed)
 
-                if surprise_temp1 > surprise:
+                aux_n_clus = np.unique(cluster_assignment_temp1).shape[0]
+                if (surprise_temp1 > surprise) and (n_clusters == aux_n_clus):
                     cluster_assignment = cluster_assignment_temp1.copy()
                     surprise = surprise_temp1
                     mem_intr_link = temp_mem_intr_link1
@@ -362,9 +397,11 @@ def solver_com_det_divis(
                     mem_intr_link.copy(),
                     np.array([clus_u, clus_v]),
                     args,
+                    approx,
                     is_directed)
 
-                if surprise_temp2 > surprise:
+                aux_n_clus = np.unique(cluster_assignment_temp2).shape[0]
+                if (surprise_temp2 > surprise) and (n_clusters == aux_n_clus):
                     cluster_assignment = cluster_assignment_temp2.copy()
                     surprise = surprise_temp2
                     mem_intr_link = temp_mem_intr_link2
@@ -381,105 +418,14 @@ def solver_com_det_divis(
         sim += 1
 
     cluster_assignment = flipping_function(
-        adjacency_matrix,
-        cluster_assignment.copy())
-
-    cluster_assignement_proper = correct_partition_labeling(
-        cluster_assignment.copy())
-    return cluster_assignement_proper, surprise
-
-
-def solver_com_det_old(
-        adjacency_matrix,
-        cluster_assignment,
-        num_sim,
-        sort_edges,
-        calculate_surprise,
-        correct_partition_labeling,
-        prob_mix,
-        flipping_function,
-        print_output=False):
-    """[summary]
-
-    :param adjacency_matrix: [description]
-    :type adjacency_matrix: [type]
-    :param cluster_assignment: [description]
-    :type cluster_assignment: [type]
-    :param num_sim: [description]
-    :type num_sim: [type]
-    :param sort_edges: [description]
-    :type sort_edges: [type]
-    :param calculate_surprise: [description]
-    :type calculate_surprise: [type]
-    :param correct_partition_labeling: [description]
-    :type correct_partition_labeling: [type]
-    :param prob_mix: [description]
-    :type prob_mix: [type]
-    :param print_output: [description], defaults to False
-    :param flipping_function:
-    :type flipping_function:
-    :type print_output: bool, optional
-    :return: [description]
-    :rtype: [type]
-    """
-    prob_random = (1 - prob_mix) / 2
-    surprise = 0
-
-    contatore_break = 0
-
-    sim = 0
-    while sim < num_sim:
-        previous_surprise = surprise
-        e_sorted = sort_edges(adjacency_matrix)
-        # print(cluster_assignment)
-        # print(e_sorted)
-        for [u, v] in e_sorted:
-            if cluster_assignment[u] != cluster_assignment[v]:
-                cluster_assignement_temp = cluster_assignment.copy()
-                random_number = np.random.uniform()
-                if (random_number > prob_mix) & (
-                        random_number <= (prob_mix + prob_random)):
-                    cluster_assignement_temp[v] = cluster_assignment[u]
-                    temp_surprise = calculate_surprise(
-                        adjacency_matrix,
-                        cluster_assignement_temp)
-                    # print(cluster_assignement_temp, cluster_assignment)
-                    # print("prob_1", u, v, temp_surprise, surprise)
-                elif random_number > (prob_mix + prob_random):
-
-                    cluster_assignement_temp[u] = cluster_assignment[v]
-                    temp_surprise = calculate_surprise(
-                        adjacency_matrix,
-                        cluster_assignement_temp)
-                    # print(cluster_assignement_temp, cluster_assignment)
-                    # print("prob_2", u, v, temp_surprise, surprise)
-                else:
-                    aux_cluster = cluster_assignement_temp[u]
-                    cluster_assignement_temp[
-                        cluster_assignement_temp == aux_cluster] = \
-                        cluster_assignement_temp[v]
-                    temp_surprise = calculate_surprise(
-                        adjacency_matrix,
-                        cluster_assignement_temp)
-                    # print(cluster_assignement_temp, cluster_assignment)
-                    # print("mixing", u, v, temp_surprise, surprise)
-
-                if temp_surprise > surprise:
-                    cluster_assignment = cluster_assignement_temp.copy()
-                    surprise = temp_surprise
-        if surprise > previous_surprise:
-            contatore_break = 0
-        else:
-            contatore_break += 1
-
-        if contatore_break >= 10:
-            break
-        if print_output:
-            print(surprise)
-        sim += 1
-
-    cluster_assignment = flipping_function(adjacency_matrix,
-                                           cluster_assignment.copy())
+        calculate_surprise=calculate_surprise,
+        adj=adjacency_matrix,
+        membership=cluster_assignment.copy(),
+        mem_intr_link=mem_intr_link.copy(),
+        args=args,
+        surprise=surprise,
+        approx=approx,
+        is_directed=is_directed)
 
     cluster_assignement_proper = correct_partition_labeling(
         cluster_assignment.copy())
diff --git a/src/surprisemes/__pycache__/Directed_Graph_Class.cpython-37.pyc b/src/surprisemes/__pycache__/Directed_Graph_Class.cpython-37.pyc
deleted file mode 100644
index c052ca97..00000000
Binary files a/src/surprisemes/__pycache__/Directed_Graph_Class.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/__pycache__/Undirected_Graph_Class.cpython-37.pyc b/src/surprisemes/__pycache__/Undirected_Graph_Class.cpython-37.pyc
deleted file mode 100644
index 7ffb2695..00000000
Binary files a/src/surprisemes/__pycache__/Undirected_Graph_Class.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/__pycache__/auxiliary_function.cpython-37.pyc b/src/surprisemes/__pycache__/auxiliary_function.cpython-37.pyc
deleted file mode 100644
index 77c208fa..00000000
Binary files a/src/surprisemes/__pycache__/auxiliary_function.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/__pycache__/comdet_functions.cpython-37.pyc b/src/surprisemes/__pycache__/comdet_functions.cpython-37.pyc
deleted file mode 100644
index 5ec84595..00000000
Binary files a/src/surprisemes/__pycache__/comdet_functions.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/__pycache__/cp_functions.cpython-37.pyc b/src/surprisemes/__pycache__/cp_functions.cpython-37.pyc
deleted file mode 100644
index 5de18c39..00000000
Binary files a/src/surprisemes/__pycache__/cp_functions.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/__pycache__/solver.cpython-37.pyc b/src/surprisemes/__pycache__/solver.cpython-37.pyc
deleted file mode 100644
index a4ed3857..00000000
Binary files a/src/surprisemes/__pycache__/solver.cpython-37.pyc and /dev/null differ
diff --git a/src/surprisemes/auxiliary_function.py b/src/surprisemes/auxiliary_function.py
deleted file mode 100644
index 168911b3..00000000
--- a/src/surprisemes/auxiliary_function.py
+++ /dev/null
@@ -1,486 +0,0 @@
-import networkx as nx
-import numpy as np
-import scipy
-from numba import jit
-from scipy.sparse import csr_matrix, coo_matrix
-from scipy.special import comb
-
-from . import comdet_functions as cd
-
-
-@jit(nopython=True)
-def compute_cn(adjacency):
-    """ Computes common neighbours table, each entry i,j of this table is the
-     number of common neighbours between i and j.
-
-    :param adjacency: Adjacency matrix.
-    :type adjacency: numpy.ndarray
-    :return: Common neighbours table.
-    :rtype: numpy.ndarray
-    """
-    cn_table = np.zeros_like(adjacency)
-    for i in np.arange(adjacency.shape[0]):
-        neighbour_i = (adjacency[i, :] + adjacency[:, i]).astype(np.bool_)
-        for j in np.arange(i + 1, adjacency.shape[0]):
-            neighbour_j = (adjacency[j, :] + adjacency[:, j]).astype(np.bool_)
-            cn_table[i, j] = cn_table[j, i] = np.multiply(neighbour_i,
-                                                          neighbour_j).sum()
-    return cn_table
-
-
-@jit(nopython=True)
-def common_neigh_init_guess(adjacency):
-    """Generates a preprocessed initial guess based on the common neighbours
-     of nodes.
-
-    :param adjacency: Adjacency matrix.
-    :type adjacency: numpy.ndarray
-    :return: Initial guess for nodes memberships.
-    :rtype: np.array
-    """
-    cn_table = compute_cn(adjacency)
-    memberships = np.array(
-        [k for k in np.arange(adjacency.shape[0], dtype=np.int32)])
-    for ii in np.arange(adjacency.shape[0]):
-        aux_node1 = np.random.choice(memberships)
-        memberships[aux_node1] = memberships[np.argmax(cn_table[aux_node1])]
-    return memberships
-
-
-def eigenvector_init_guess(adjacency, is_directed):
-    """Generates an initial guess for core periphery detection method: nodes
-    with higher eigenvector centrality are in the core.
-
-    :param adjacency: Adjacency matrix.
-    :type adjacency: np.ndarray
-    :param is_directed: True if the network is directed.
-    :type is_directed: bool
-    :return: Initial guess.
-    :rtype: np.ndarray
-    """
-    # TODO: Vedere come funziona la parte pesata
-
-    n_nodes = adjacency.shape[0]
-    aux_nodes = int(np.ceil((n_nodes * 5) / 100))
-    if is_directed:
-        graph = nx.from_numpy_array(adjacency, create_using=nx.DiGraph)
-        centra = nx.eigenvector_centrality_numpy(graph)
-        centra1 = np.array([centra[key] for key in centra])
-        membership = np.zeros_like(centra1, dtype=np.int32)
-        membership[np.argsort(centra1)[::-1][:aux_nodes]] = 1
-
-    else:
-        graph = nx.from_numpy_array(adjacency, create_using=nx.Graph)
-        centra = nx.eigenvector_centrality_numpy(graph)
-        centra1 = np.array([centra[key] for key in centra])
-        membership = np.zeros_like(centra1, dtype=np.int32)
-        print(aux_nodes)
-        print(membership[np.argsort(centra1)[::-1]][:aux_nodes])
-        membership[np.argsort(centra1)[::-1][:aux_nodes]] = 1
-
-    return membership
-
-
-def fixed_clusters_init_guess_cn(adjacency, n_clust):
-    """ Generates an intial guess with a fixed number 'n' of clusters.
-    Nodes are organised in clusters based on the number of common neighbors.
-    The starting members of clusters are the 'n' nodes with higher
-    degrees/strengths.
-
-    :param adjacency: Adjacency matrix.
-    :type adjacency: numpy.ndarray
-    :param n_clust: Partitions number.
-    :type n_clust: int
-    :return: Initial guess.
-    :rtype: numpy.ndarray
-    """
-    aux_memb = np.ones(adjacency.shape[0], dtype=np.int32) * n_clust
-    deg = adjacency.sum(axis=1)
-    aux_args_sort = np.argsort(deg)[::-1][0:n_clust]
-    for memb, index in enumerate(aux_args_sort):
-        aux_memb[index] = memb
-
-    common_neigh = compute_cn(adjacency)
-    aux = np.nonzero(aux_memb == n_clust)[0]
-    np.random.shuffle(aux)
-    for node in aux:
-        aux_list = np.nonzero(aux_memb != n_clust)[0]
-        node_index = aux_list[np.argmax(common_neigh[node, aux_list])]
-        if isinstance(node_index, np.ndarray):
-            node_index = np.random.choice(node_index)
-        aux_memb[node] = aux_memb[node_index]
-
-    return aux_memb
-
-
-def compute_degree(a, is_directed):
-    """Returns matrix *a* degree sequence.
-
-    :param a: Matrix.
-    :type a:  numpy.ndarray
-    :param is_directed: True if the matrix is directed.
-    :type is_directed: bool
-    :return: Degree sequence.
-    :rtype: numpy.ndarray.
-    """
-    # if the matrix is a numpy array
-    if is_directed:
-        if type(a) == np.ndarray:
-            return np.sum(a > 0, 0), np.sum(a > 0, 1)
-        # if the matrix is a scipy sparse matrix
-        elif type(a) in [csr_matrix, coo_matrix]:
-            return np.sum(a > 0, 0).A1, np.sum(a > 0, 1).A1
-    else:
-        if type(a) == np.ndarray:
-            return np.sum(a > 0, 1)
-        # if the matrix is a scipy sparse matrix
-        elif type(a) in [csr_matrix, coo_matrix]:
-            return np.sum(a > 0, 1).A1
-
-
-def compute_strength(a, is_directed):
-    """Returns matrix *a* strength sequence.
-
-    :param a: Matrix.
-    :type a: numpy.ndarray
-    :param is_directed: True if the matrix is directed.
-    :type is_directed: bool
-    :return: Strength sequence.
-    :rtype: numpy.ndarray
-    """
-    if is_directed:
-        # if the matrix is a numpy array
-        if type(a) == np.ndarray:
-            return np.sum(a, 0), np.sum(a, 1)
-        # if the matrix is a scipy sparse matrix
-        elif type(a) in [csr_matrix, coo_matrix]:
-            return np.sum(a, 0).A1, np.sum(a, 1).A1
-    else:
-        # if the matrix is a numpy array
-        if type(a) == np.ndarray:
-            return np.sum(a, 1)
-        # if the matrix is a scipy sparse matrix
-        elif type(a) in [csr_matrix, coo_matrix]:
-            return np.sum(a, 1).A1
-
-
-def from_edgelist(edgelist, is_sparse, is_directed):
-    """Returns np.ndarray or scipy.sparse matrix from edgelist.
-
-    :param edgelist: List of edges, eache edge must be given as a 2-tuples
-     (u,v).
-    :type edgelist: list or numpy.ndarray
-    :param is_sparse: If true the returned matrix is sparse.
-    :type is_sparse: bool
-    :param is_directed: If true the graph is directed.
-    :type is_directed: bool
-    :return: Adjacency matrix.
-    :rtype: numpy.ndarray or scipy.sparse
-    """
-    # TODO: vedere che tipo di sparse e'
-    if is_directed:
-        g = nx.DiGraph()
-    else:
-        g = nx.Graph()
-    g.add_edges_from(edgelist)
-    if is_sparse:
-        return nx.to_scipy_sparse_matrix(g)
-    else:
-        return nx.to_numpy_array(g)
-
-
-def from_weighted_edgelist(edgelist, is_sparse, is_directed):
-    """Returns np.ndarray or scipy.sparse matrix from edgelist.
-
-    :param edgelist: List of weighted edges, eache edge must be given as a
-     3-tuples (u,v,w).
-    :type edgelist: [type]
-    :param is_sparse: If true the returned matrix is sparse.
-    :type is_sparse: bool
-    :param is_directed: If true the graph is directed.
-    :type is_directed: bool
-    :return: Weighted adjacency matrix.
-    :rtype: numpy.ndarray or scipy.sparse
-    """
-    if is_directed:
-        g = nx.DiGraph()
-    else:
-        g = nx.Graph()
-    g.add_weighted_edges_from(edgelist)
-    if is_sparse:
-        return nx.to_scipy_sparse_matrix(g)
-    else:
-        return nx.to_numpy_array(g)
-
-
-def check_symmetric(a, is_sparse, rtol=1e-05, atol=1e-08):
-    """Checks if the matrix is symmetric.
-
-    :param a: Matrix.
-    :type a: numpy.ndarray or scipy.sparse
-    :param is_sparse: If true the matrix is sparse.
-    :type is_sparse: bool
-    :param rtol: Tuning parameter, defaults to 1e-05.
-    :type rtol: float, optional
-    :param atol: Tuning parameter, defaults to 1e-08.
-    :type atol: float, optional
-    :return: True if the matrix is symmetric.
-    :rtype: bool
-    """
-    if is_sparse:
-        return np.all(np.abs(a - a.T) < atol)
-    else:
-        return np.allclose(a, a.T, rtol=rtol, atol=atol)
-
-
-def check_adjacency(adjacency, is_sparse, is_directed):
-    """Functions checking the _validty_ of the adjacency matrix.
-
-    :param adjacency: Adjacency matrix.
-    :type adjacency: numpy.ndarray or scipy.sparse
-    :param is_sparse: If true the matrix is sparse.
-    :type is_sparse: bool
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :raises TypeError: Matrix not square.
-    :raises ValueError: Negative entries.
-    :raises TypeError: Matrix not symmetric.
-    """
-    if adjacency.shape[0] != adjacency.shape[1]:
-        raise TypeError(
-            "Adjacency matrix must be square. If you are passing an edgelist"
-            " use the positional argument 'edgelist='.")
-    if np.sum(adjacency < 0):
-        raise ValueError(
-            "The adjacency matrix entries must be positive."
-        )
-    if (not check_symmetric(adjacency, is_sparse)) and (not is_directed):
-        raise TypeError(
-            "The adjacency matrix seems to be not symmetric, we suggest to use"
-            " 'DirectedGraphClass'.")
-
-
-@jit(nopython=True, fastmath=True)
-def sumLogProbabilities(nextLogP, logP):
-    if nextLogP == 0:
-        stop = True
-    else:
-        stop = False
-        if nextLogP > logP:
-            common = nextLogP
-            diffExponent = logP - common
-        else:
-            common = logP
-            diffExponent = nextLogP - common
-
-        logP = common + ((np.log10(1 + 10 ** diffExponent)) / np.log10(10))
-
-        if (nextLogP - logP) > -4:
-            stop = True
-
-    return logP, stop
-
-
-@jit(nopython=True, fastmath=True)
-def logc(n, k):
-    if k == n:
-        return 0
-    elif (n > 1000) & (k > 1000):  # Stirling's binomial coeff approximation
-        return logStirFac(n) - logStirFac(k) - logStirFac(n - k)
-    else:
-        t = n - k
-        if t < k:
-            t = k
-        logC = sumRange(t + 1, n) - sumFactorial(n - t)
-        return logC
-
-
-@jit(nopython=True, fastmath=True)
-def logStirFac(n):
-    if n <= 1:
-        return 1.0
-    else:
-        return -n + n * np.log10(n) + np.log10(
-            n * (1 + 4.0 * n * (1.0 + 2.0 * n))) / 6.0 + np.log10(np.pi) / 2.0
-
-
-@jit(nopython=True, fastmath=True)
-def sumRange(xmin, xmax):
-    """[summary]
-
-    :param xmin: [description]
-    :type xmin: [type]
-    :param xmax: [description]
-    :type xmax: [type]
-    :return: [description]
-    :rtype: [type]
-    """
-    csum = 0
-    for i in np.arange(xmin, xmax + 1):
-        csum += np.log10(i)
-    return csum
-
-
-@jit(nopython=True, fastmath=True)
-def sumFactorial(n):
-    csum = 0
-    if n > 1:
-        for i in np.arange(2, n + 1):
-            csum += np.log10(i)
-    return csum
-
-
-def shuffled_edges(adjacency_matrix, is_directed):
-    """Shuffles edges randomly.
-
-    :param adjacency_matrix: Matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param is_directed: True if graph is directed.
-    :type is_directed: bool
-    :return: Shuffled edgelist.
-    :rtype: mumpy.ndarray
-    """
-    adj = adjacency_matrix.astype(bool).astype(np.int16)
-    if not is_directed:
-        adj = np.triu(adj)
-    edges = np.stack(adj.nonzero(), axis=-1)
-    np.random.shuffle(edges)
-    shuff_edges = edges.astype(np.int32)
-    return shuff_edges
-
-
-def jaccard_sorted_edges(adjacency_matrix):
-    """Returns edges ordered based on jaccard index.
-
-    :param adjacency_matrix: Matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :return: Ordered edgelist.
-    :rtype: numpy.ndarray
-    """
-    G = nx.from_numpy_matrix(adjacency_matrix)
-    jacc = nx.jaccard_coefficient(G, ebunch=G.edges())
-    jacc_array = []
-    for u, v, p in jacc:
-        jacc_array += [[u, v, p]]
-    jacc_array = np.array(jacc_array)
-    jacc_array = jacc_array[jacc_array[:, 2].argsort()][::-1]
-    sorted_edges = jacc_array[:, 0:2]
-    sorted_edges = sorted_edges.astype(np.int32)
-    return sorted_edges
-
-
-def evaluate_surprise_clust_bin(adjacency_matrix,
-                                cluster_assignment,
-                                is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-     for a binary network.
-
-    :param adjacency_matrix: Binary adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.
-    :rtype: float
-    """
-    if is_directed:
-        # intracluster links
-        p = cd.intracluster_links(adjacency_matrix,
-                                  cluster_assignment)
-        p = int(p)
-        # All the possible intracluster links
-        M = cd.calculate_possible_intracluster_links(cluster_assignment,
-                                                     is_directed)
-        # Observed links
-        m = np.sum(adjacency_matrix.astype(bool))
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        F = n * (n - 1)
-    else:
-        # intracluster links
-        p = cd.intracluster_links(adjacency_matrix,
-                                  cluster_assignment)
-        p = int(p / 2)
-        # All the possible intracluster links
-        M = int(cd.calculate_possible_intracluster_links(cluster_assignment,
-                                                         is_directed))
-        # Observed links
-        m = np.sum(adjacency_matrix.astype(bool)) / 2
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        F = int((n * (n - 1)) / 2)
-
-    surprise = surprise_hypergeometric(F, p, M, m)
-    return surprise
-
-
-def surprise_hypergeometric(F, p, M, m):
-    surprise = 0
-    min_p = min(M, m)
-    for p_loop in np.arange(p, min_p + 1):
-        surprise += (comb(M, p_loop, exact=True) * comb(
-            F - M, m - p_loop,
-            exact=True)) / comb(F,
-                                m,
-                                exact=True)
-    return surprise
-
-
-def evaluate_surprise_clust_weigh(adjacency_matrix,
-                                  cluster_assignment,
-                                  is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-     for a weighted network.
-
-    :param adjacency_matrix: Weighted adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.
-    :rtype: float
-    """
-    if is_directed:
-        # intracluster weights
-        w = cd.intracluster_links(adjacency_matrix,
-                                  cluster_assignment)
-        # intracluster possible links
-        Vi = cd.calculate_possible_intracluster_links(cluster_assignment,
-                                                      is_directed)
-        # Total Weight
-        W = np.sum(adjacency_matrix)
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        V = n * (n - 1)
-        # extracluster links
-        Ve = V - Vi
-    else:
-        # intracluster weights
-        w = cd.intracluster_links(adjacency_matrix,
-                                  cluster_assignment) / 2
-        # intracluster possible links
-        Vi = cd.calculate_possible_intracluster_links(cluster_assignment,
-                                                      is_directed)
-        # Total Weight
-        W = np.sum(adjacency_matrix) / 2
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        V = int((n * (n - 1)) / 2)
-        # extracluster links
-        Ve = V - Vi
-
-    surprise = surprise_negative_hypergeometric(Vi, w, Ve, W, V)
-    return surprise
-
-
-def surprise_negative_hypergeometric(Vi, w, Ve, W, V):
-    """Computes the negative hypergeometric distribution.
-    """
-    surprise = 0
-    for w_loop in range(w, W):
-        surprise += ((comb(Vi + w_loop - 1, w_loop, exact=True) * comb(
-            Ve + W - w_loop, W - w_loop, exact=True)) /
-                     comb(V + W, W, exact=True))
-    return surprise
diff --git a/src/surprisemes/comdet_functions.py b/src/surprisemes/comdet_functions.py
deleted file mode 100644
index 18d55265..00000000
--- a/src/surprisemes/comdet_functions.py
+++ /dev/null
@@ -1,557 +0,0 @@
-import numpy as np
-from numba import jit
-
-from . import auxiliary_function as ax
-
-
-def calculate_possible_intracluster_links(partitions, is_directed):
-    """Computes the number of possible links, given nodes memberships.
-
-    :param partitions: Nodes memberships.
-    :type partitions: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Total number of possible edges.
-    :rtype: float
-    """
-    _, counts = np.unique(partitions, return_counts=True)
-    nr_poss_intr_clust_links = np.sum(counts * (counts - 1))
-    if is_directed:
-        return nr_poss_intr_clust_links
-    else:
-        return nr_poss_intr_clust_links / 2
-
-
-@jit(nopython=True)
-def calculate_possible_intracluster_links_new(partitions, is_directed):
-    """Computes the number of possible links, given nodes memberships.
-    Faster implementation compiled in "nopython" mode.
-
-    :param partitions: Nodes memberships.
-    :type partitions: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Total number of possible edges.
-    :rtype: float
-    """
-    counts = np.bincount(partitions)
-    nr_poss_intr_clust_links = np.sum(counts * (counts - 1))
-    if is_directed:
-        return nr_poss_intr_clust_links
-    else:
-        return nr_poss_intr_clust_links / 2
-
-
-@jit(nopython=True)
-def intracluster_links(adj, partitions):
-    """Computes intracluster links or weights.
-
-    :param adj: Adjacency matrix.
-    :type adj: numpy.ndarray
-    :param partitions: Nodes memberships.
-    :type partitions: numpy.ndarray
-    :return: Number of intra-cluster links/weights.
-    :rtype: float
-    """
-    nr_intr_clust_links = 0
-    clust_labels = np.unique(partitions)
-    for lab in clust_labels:
-        indices = np.where(partitions == lab)[0]
-        nr_intr_clust_links += intracluster_links_aux(adj, indices)
-    return nr_intr_clust_links
-
-
-@jit(nopython=True)
-def intracluster_links_new(adj, clust_labels, partitions):
-    """Computes intracluster links or weights. New implementation
-
-    :param adj: Adjacency matrix.
-    :type adj: numpy.ndarray
-    :param clust_labels: Labels of changed clusters.
-    :type clust_labels: numpy.ndarray
-    :param partitions: Nodes memberships.
-    :type partitions: numpy.ndarray
-    :return: Number of intra-cluster links/weights.                
-    :rtype: numpy.ndarray
-    """
-    # print(clust_labels, clust_labels.shape, partitions)
-    nr_intr_clust_links = np.zeros(clust_labels.shape[0])
-    for ii, lab in enumerate(clust_labels):
-        indices = np.where(partitions == lab)[0]
-        nr_intr_clust_links[ii] = intracluster_links_aux(adj, indices)
-    return nr_intr_clust_links
-
-
-@jit(nopython=True)
-def intracluster_links_aux(adj, indices):
-    """Computes intra-cluster links or weights given nodes indices.
-
-    :param adj: [description]
-    :type adj: [type]
-    :param indices: [description]
-    :type indices: [type]
-    :return: [description]
-    :rtype: [type]
-    """
-    n_links = 0
-    for ii in indices:
-        for jj in indices:
-            n_links += adj[ii, jj]
-    return n_links
-
-
-def calculate_surprise_logsum_clust_bin(adjacency_matrix,
-                                        cluster_assignment,
-                                        is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-    for a binary network.
-
-    :param adjacency_matrix: Binary adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.
-    :rtype: float
-    """
-    if is_directed:
-        # intracluster links
-        p = intracluster_links(
-            adjacency_matrix,
-            cluster_assignment)
-        int_links = int(p)
-        # All the possible intracluster links
-        poss_int_links = calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed)
-        # Observed links
-        obs_links = np.sum(adjacency_matrix.astype(bool))
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        poss_links = n * (n - 1)
-    else:
-        # intracluster links
-        p = intracluster_links(
-            adjacency_matrix,
-            cluster_assignment)
-        int_links = int(p / 2)
-        # All the possible intracluster links
-        poss_int_links = int(calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed))
-        # Observed links
-        obs_links = np.sum(adjacency_matrix.astype(bool)) / 2
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        poss_links = int((n * (n - 1)) / 2)
-
-    surprise = surprise_logsum_clust_bin(
-        poss_links,
-        int_links,
-        poss_int_links,
-        obs_links)
-
-    return surprise
-
-
-def calculate_surprise_logsum_clust_bin_new(
-        adjacency_matrix,
-        cluster_assignment,
-        mem_intr_link,
-        clust_labels,
-        args,
-        is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-     for a binary network. New faster implementation reducing the number of
-      redundant computations.
-
-    :param adjacency_matrix: Binary adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.                                                       
-    :type cluster_assignment: numpy.ndarray
-    :param mem_intr_link: Intracluster links per cluster
-    :type mem_intr_link: np.array
-    :param clust_labels: Labels of changed clusters.
-    :type clust_labels: np.array
-    :param args: (Observed links, nodes number, possible links)
-    :type args: (int, int)
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.                                                                              
-    :rtype: float                                                                                       
-    """
-    if is_directed:
-        # intracluster links
-        if len(clust_labels):
-            int_links = intracluster_links_new(
-                adj=adjacency_matrix,
-                clust_labels=clust_labels,
-                partitions=cluster_assignment)
-
-            for node_label, nr_links in zip(clust_labels, int_links):
-                mem_intr_link[node_label] = nr_links
-
-        p = np.sum(mem_intr_link)
-        int_links = int(p)
-        # All the possible intracluster links                                                           
-        poss_int_links = calculate_possible_intracluster_links_new(
-            cluster_assignment,
-            is_directed)
-        # Observed links                                                                                
-        obs_links = args[0]
-        # Possible link
-        poss_links = args[1]
-    else:
-        # intracluster links
-        if len(clust_labels):
-            int_links = intracluster_links_new(
-                adj=adjacency_matrix,
-                clust_labels=clust_labels,
-                partitions=cluster_assignment)
-
-            for node_label, nr_links in zip(clust_labels, int_links):
-                mem_intr_link[node_label] = nr_links
-
-        p = np.sum(mem_intr_link)
-        int_links = int(p / 2)
-        # All the possible intracluster links                                                           
-        poss_int_links = int(calculate_possible_intracluster_links_new(
-            cluster_assignment,
-            is_directed))
-        # Observed links                                                                                
-        obs_links = args[0] / 2
-        # Possible links
-        poss_links = int(args[1] / 2)
-
-    surprise = surprise_logsum_clust_bin(
-        poss_links,
-        int_links,
-        poss_int_links,
-        obs_links)
-
-    return surprise, mem_intr_link
-
-
-@jit(nopython=True)
-def surprise_logsum_clust_bin(F, p, M, m):
-    """[summary]
-
-    :param F: [description]
-    :type F: [type]
-    :param p: [description]
-    :type p: [type]
-    :param M: [description]
-    :type M: [type]
-    :param m: [description]
-    :type m: [type]
-    :return: [description]
-    :rtype: [type]
-    """
-    stop = False
-    min_p = min(M, m)
-
-    logP = loghyperprobability(F, p, M, m)
-    for p_loop in np.arange(p, min_p + 1):
-        if p_loop == p:
-            continue
-        nextLogP = loghyperprobability(F, p_loop, M, m)
-        [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
-        if stop:
-            break
-
-    return -logP
-
-
-@jit(nopython=True)
-def loghyperprobability(F, p, M, m):
-    """Evaluates logarithmic hypergeometric distribution
-
-    :param F:
-    :type F:
-    :param p:
-    :type p:
-    :param M:
-    :type M:
-    :param m:
-    :type m:
-    :return:
-    :rtype:
-    """
-    logH = ax.logc(M, p) + ax.logc(F - M, m - p) - ax.logc(F, m)
-    return logH
-
-
-def calculate_surprise_logsum_clust_weigh(
-        adjacency_matrix,
-        cluster_assignment,
-        is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-     for a weighted network.
-
-    :param adjacency_matrix: Weighted adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.
-    :rtype: float
-    """
-    if is_directed:
-        # intracluster weights
-        w = intracluster_links(adjacency_matrix,
-                               cluster_assignment)
-        # intracluster possible links
-        poss_intr_links = calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed)
-        # Total Weight
-        tot_weights = np.sum(adjacency_matrix)
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        poss_links = n * (n - 1)
-        # extracluster links
-        inter_links = poss_links - poss_intr_links
-    else:
-        # intracluster weights
-        w = intracluster_links(adjacency_matrix,
-                               cluster_assignment) / 2
-        # intracluster possible links
-        poss_intr_links = calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed)
-        # Total Weight
-        tot_weights = np.sum(adjacency_matrix) / 2
-        # Possible links
-        n = adjacency_matrix.shape[0]
-        poss_links = int((n * (n - 1)) / 2)
-        # extracluster links
-        inter_links = poss_links - poss_intr_links
-
-    surprise = surprise_logsum_clust_weigh(
-        poss_intr_links,
-        w,
-        inter_links,
-        tot_weights,
-        poss_links)
-    return surprise
-
-
-def calculate_surprise_logsum_clust_weigh_new(
-        adjacency_matrix,
-        cluster_assignment,
-        mem_intr_link,
-        clust_labels,
-        args,
-        is_directed):
-    """Calculates the logarithm of the surprise given the current partitions
-     for a weighted network. New faster implementation.
-
-    :param adjacency_matrix: Weighted adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Nodes memberships.
-    :type cluster_assignment: numpy.ndarray
-    :param mem_intr_link:
-    :type mem_intr_link:
-    :param clust_labels:
-    :type clust_labels:
-    :param args:
-    :type args:
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise.
-    :rtype: float
-    """
-    if is_directed:
-        # intracluster weights
-        if len(clust_labels):
-            w = intracluster_links_new(
-                adj=adjacency_matrix,
-                clust_labels=clust_labels,
-                partitions=cluster_assignment)
-
-            for node_label, nr_links in zip(clust_labels, w):
-                mem_intr_link[node_label] = nr_links
-
-        p = np.sum(mem_intr_link)
-        intr_weights = int(p)
-        # intracluster possible links
-        poss_intr_links = calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed)
-        # Total Weight
-        tot_weights = args[0]
-        # Possible links
-        poss_links = args[1]
-        # extracluster links
-        inter_links = poss_links - poss_intr_links
-    else:
-        # intracluster weights
-        if len(clust_labels):
-            w = intracluster_links_new(
-                adj=adjacency_matrix,
-                clust_labels=clust_labels,
-                partitions=cluster_assignment)
-
-            for node_label, nr_links in zip(clust_labels, w):
-                mem_intr_link[node_label] = nr_links
-
-        p = np.sum(mem_intr_link)
-        intr_weights = int(p / 2)
-        # intracluster possible links
-        poss_intr_links = calculate_possible_intracluster_links(
-            cluster_assignment,
-            is_directed)
-        # Total Weight
-        tot_weights = int(args[0] / 2)
-        # Possible links
-        poss_links = int(args[1] / 2)
-        # extracluster links
-        inter_links = poss_links - poss_intr_links
-
-    surprise = surprise_logsum_clust_weigh(
-        poss_intr_links,
-        intr_weights,
-        inter_links,
-        tot_weights,
-        poss_links)
-    return surprise, mem_intr_link
-
-
-@jit(nopython=True)
-def surprise_logsum_clust_weigh(Vi, w, Ve, W, V):
-    """
-
-    :param Vi: [description]
-    :type Vi: [type]
-    :param w: [description]
-    :type w: [type]
-    :param Ve: [description]
-    :type Ve: [type]
-    :param W: [description]
-    :type W: [type]
-    :param V: [description]
-    :type V: [type]
-    :return: [description]
-    :rtype: [type]
-    """
-    stop = False
-
-    logP = lognegativehyperprobability(Vi, w, Ve, W, V)
-    for w_loop in range(w, W):
-        if w_loop == w:
-            continue
-        nextLogP = lognegativehyperprobability(Vi, w_loop, Ve, W, V)
-        [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
-        if stop:
-            break
-
-    return -logP
-
-
-@jit(nopython=True)
-def lognegativehyperprobability(Vi, w, Ve, W, V):
-    """Evaluates logarithmic hypergeometric distribution
-
-    :param Vi:
-    :type Vi:
-    :param w:
-    :type w:
-    :param Ve:
-    :type Ve:
-    :param W:
-    :type W:
-    :param V:
-    :type V:
-    :return:
-    :rtype:
-    """
-    logh = ax.logc(
-        Vi + w - 1,
-        w) + ax.logc(
-        Ve + W - w,
-        W - w) - ax.logc(
-        V + W,
-        W)
-    return logh
-
-
-def labeling_communities(partitions):
-    """Gives labels to communities from 0 to number of communities minus one.
-
-    :param partitions: Nodes memberships.
-    :type partitions: numpy.ndarray
-    :return: Re-labeled nodes memberships.
-    :rtype: numpy.ndarray
-    """
-    different_partitions = np.unique(partitions, return_counts=True)
-    aux_argsort = np.argsort(different_partitions[1])[::-1]
-    ordered_clusters = different_partitions[0][aux_argsort]
-    new_partitioning = partitions.copy()
-    new_partitions = np.array([k for k in range(len(ordered_clusters))])
-    for old_part, new_part in zip(ordered_clusters, new_partitions):
-        # print(old_part,new_part)
-        indices_old = np.where(partitions == old_part)[0]
-        new_partitioning[indices_old] = new_part
-
-    return new_partitioning
-
-
-def flipping_function_comdet(comm):
-    """Changes the membership of a randomly selected node.
-
-    :param comm: Nodes memberships.
-    :type comm: numpy.ndarray
-    :return: New nodes memberships.
-    :rtype: numpy.ndarray
-    """
-    labels_set = np.unique(comm)
-    node_index = np.random.randint(0, comm.shape[0])
-    remaining_labels = labels_set[labels_set != comm[node_index]]
-    if remaining_labels.size != 0:
-        new_label = np.random.choice(remaining_labels)
-        comm[node_index] = new_label
-    return comm
-
-
-def flipping_function_comdet_new(adj,
-                                 membership,
-                                 is_directed):
-    obs_links = int(np.sum(adj.astype(bool)))
-    n_nodes = int(adj.shape[0])
-    poss_links = int(n_nodes * (n_nodes - 1))
-    args = (obs_links, poss_links)
-
-    surprise = calculate_surprise_logsum_clust_bin(
-        adjacency_matrix=adj,
-        cluster_assignment=membership,
-        is_directed=is_directed)
-
-    mem_intr_link = np.zeros(membership.shape[0], dtype=np.int32)
-    # print(np.unique(membership), mem_intr_link, membership)
-    for ii in np.unique(membership):
-        indices = np.where(membership == ii)[0]
-        mem_intr_link[ii] = intracluster_links_aux(adj, indices)
-
-    # print("siamo qui")
-    for node, node_label in zip(np.arange(membership.shape[0]), membership):
-
-        for new_clust in np.unique(membership):
-            if node_label != new_clust:
-                aux_membership = membership.copy()
-                aux_membership[node] = new_clust
-                # print(np.array([node_label, new_clust]))
-                temp_surprise, temp_mem_intr_link = calculate_surprise_logsum_clust_bin_new(
-                    adjacency_matrix=adj,
-                    cluster_assignment=aux_membership,
-                    mem_intr_link=mem_intr_link.copy(),
-                    clust_labels=np.array([node_label, new_clust]),
-                    args=args,
-                    is_directed=is_directed)
-                if temp_surprise > surprise:
-                    membership = aux_membership.copy()
-                    surprise = temp_surprise
-                    mem_intr_link = temp_mem_intr_link
-
-    return membership
diff --git a/src/surprisemes/cp_functions.py b/src/surprisemes/cp_functions.py
deleted file mode 100644
index 0b84b0f8..00000000
--- a/src/surprisemes/cp_functions.py
+++ /dev/null
@@ -1,257 +0,0 @@
-import numpy as np
-from numba import jit
-
-from . import auxiliary_function as ax
-
-
-@jit(nopython=True)
-def compute_sum(adj, row_indices, column_indices):
-    """Computes number of links/total weight given nodes indices.
-
-    :param adj: Adjacency matrix.
-    :type adj: numpy.ndarray
-    :param row_indices: Row indices.
-    :type row_indices: numpy.ndarray
-    :param column_indices: Columns indices.
-    :type column_indices: numpy.ndarray
-    :return: Total links/weights
-    :rtype: float
-    """
-    Sum = 0.0
-    for ii in row_indices:
-        for jj in column_indices:
-            Sum += adj[ii, jj]
-    return Sum
-
-
-def calculate_surprise_logsum_cp_weigh(adjacency_matrix,
-                                       cluster_assignment,
-                                       is_directed):
-    """Computes core-periphery weighted log-surprise given a certain nodes'
-     partitioning.
-
-    :param adjacency_matrix: Weighted adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Core periphery assigments.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise
-    :rtype: float
-    """
-    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
-    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
-
-    if is_directed:
-        n_c = core_nodes.shape[0]
-        n_x = periphery_nodes.shape[0]
-        p_c = n_c * (n_c - 1)
-        p_x = n_c * n_x * 2
-
-        w_c = compute_sum(adjacency_matrix, core_nodes, core_nodes)
-        w_x = compute_sum(adjacency_matrix,
-                          core_nodes,
-                          periphery_nodes) + compute_sum(adjacency_matrix,
-                                                         periphery_nodes,
-                                                         core_nodes)
-
-        w = np.sum(adjacency_matrix)
-        w_p = w - w_c - w_x
-        n = n_c + n_x
-        p = n * (n - 1)
-        p_p = p - p_c - p_x
-
-    else:
-        n_c = core_nodes.shape[0]
-        n_x = periphery_nodes.shape[0]
-        p_c = n_c * (n_c - 1) / 2
-        p_x = n_c * n_x
-
-        w_c = (compute_sum(adjacency_matrix, core_nodes, core_nodes)) / 2
-        w_x = (compute_sum(adjacency_matrix,
-                           core_nodes,
-                           periphery_nodes) + compute_sum(adjacency_matrix,
-                                                          periphery_nodes,
-                                                          core_nodes)) / 2
-
-        w = np.sum(adjacency_matrix) / 2
-        w_p = (w - w_c - w_x) / 2
-        n = n_c + n_x
-        p = n * (n - 1) / 2
-        p_p = p - p_c - p_x
-
-    surprise = surprise_bipartite_logsum_CP_Weigh(p, p_c, p_x, p_p,
-                                                  w, w_c, w_x, w_p)
-    return surprise
-
-
-@jit(nopython=True)
-def surprise_bipartite_logsum_CP_Weigh(p, p_c, p_x, p_p, w, w_c, w_x, w_p):
-    stop = False
-    first_loop_break = False
-
-    logP = logMultiHyperProbabilityWeight(p, p_c, p_x, p_p, w, w_c, w_x, w_p)
-    for w_c_loop in range(w_c, w + 1):
-        for w_x_loop in range(w_x, w + 1 - w_c_loop):
-            if (w_c_loop == w_c) & (w_x_loop == w_x):
-                continue
-            w_p_loop = w - w_c_loop - w_x_loop
-            nextLogP = logMultiHyperProbabilityWeight(p, p_c, p_x,
-                                                      p_p, w, w_c_loop,
-                                                      w_x_loop, w_p_loop)
-            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
-
-            if stop:
-                first_loop_break = True
-                break
-        if first_loop_break:
-            break
-    return -logP
-
-
-@jit(nopython=True)
-def logMultiHyperProbabilityWeight(p, p_c, p_x, p_p, w, w_c, w_x, w_p):
-    """Computes the logarithm of the Negative Multinomial Hypergeometric
-     distribution."""
-    logH = (ax.logc(p_c + w_c - 1, w_c) + ax.logc(p_x + w_x - 1, w_x)
-            + ax.logc(p_p + w_p - 1, w_p) - ax.logc(p + w, w))
-    return logH
-
-
-def calculate_surprise_logsum_cp_bin(adjacency_matrix,
-                                     cluster_assignment,
-                                     is_directed):
-    """Computes core-periphery binary log-surprise given a certain nodes'
-     partitioning.
-
-    :param adjacency_matrix: Binary adjacency matrix.
-    :type adjacency_matrix: numpy.ndarray
-    :param cluster_assignment: Core periphery assigments.
-    :type cluster_assignment: numpy.ndarray
-    :param is_directed: True if the graph is directed.
-    :type is_directed: bool
-    :return: Log-surprise
-    :rtype: float
-    """
-    core_nodes = np.unique(np.where(cluster_assignment == 0)[0])
-    periphery_nodes = np.unique(np.where(cluster_assignment == 1)[0])
-
-    if is_directed:
-        n_c = (core_nodes).shape[0]
-        n_x = (periphery_nodes).shape[0]
-        p_c = n_c * (n_c - 1)
-        p_x = n_c * n_x * 2
-
-        l_c = compute_sum(adjacency_matrix, core_nodes, core_nodes)
-        l_x = compute_sum(adjacency_matrix, core_nodes,
-                          periphery_nodes) + compute_sum(adjacency_matrix,
-                                                         periphery_nodes,
-                                                         core_nodes)
-
-        l_t = np.sum(adjacency_matrix)
-        n = n_c + n_x
-        p = n * (n - 1)
-
-    else:
-        n_c = (core_nodes).shape[0]
-        n_x = (periphery_nodes).shape[0]
-        p_c = (n_c * (n_c - 1)) / (2)
-        p_x = n_c * n_x
-
-        l_c = compute_sum(adjacency_matrix, core_nodes, core_nodes) / 2
-        l_x = (compute_sum(adjacency_matrix, core_nodes,
-                           periphery_nodes) + compute_sum(adjacency_matrix,
-                                                          periphery_nodes,
-                                                          core_nodes)) / 2
-
-        l_t = np.sum(adjacency_matrix) / 2
-        n = n_c + n_x
-        p = (n * (n - 1)) / 2
-
-    if (p_c + p_x) < (l_c + l_x):
-        return 0
-
-    surprise = surprise_bipartite_logsum_CP_Bin(p, p_c, p_x, l_t, l_c, l_x)
-    return surprise
-
-
-@jit(nopython=True)
-def surprise_bipartite_logsum_CP_Bin(p, p_c, p_x, l, l_c, l_x):
-    stop = False
-    first_loop_break = False
-
-    min_l_p = min(l, p_c + p_x)
-
-    logP = logMultiHyperProbability(p, p_c, p_x, l, l_c, l_x)
-    for l_c_loop in range(l_c, min_l_p + 1):
-        for l_x_loop in range(l_x, min_l_p + 1 - l_c_loop):
-            if (l_c_loop == l_c) & (l_x_loop == l_x):
-                continue
-            nextLogP = logMultiHyperProbability(p, p_c, p_x, l, l_c_loop,
-                                                l_x_loop)
-            [logP, stop] = ax.sumLogProbabilities(nextLogP, logP)
-
-            if stop:
-                first_loop_break = True
-                break
-        if first_loop_break:
-            break
-
-    return -logP
-
-
-@jit(nopython=True)
-def logMultiHyperProbability(p, p_c, p_x, l, l_c, l_x):
-    """Computes the logarithm of the Multinomial Hypergeometric
-     distribution."""
-    logH = ax.logc(p_c, l_c) + ax.logc(p_x, l_x) + ax.logc(
-        p - p_c - p_x,
-        l - l_c - l_x) - ax.logc(p, l)
-    return logH
-
-
-def labeling_core_periphery(adjacency_matrix, cluster_assignment):
-    """Function assigning the core and periphery labels based on link density
-
-    :param adjacency_matrix: [description]
-    :type adjacency_matrix: [type]
-    :param cluster_assignment: [description]
-    :type cluster_assignment: [type]
-    :return: [description]
-    :rtype: [type]
-    """
-    core_nodes = np.where(cluster_assignment == 0)[0]
-    periphery_nodes = np.where(cluster_assignment == 1)[0]
-    l_core = np.sum(
-        adjacency_matrix[np.ix_(list(core_nodes), list(core_nodes))] > 0)
-    l_periphery = np.sum(adjacency_matrix[np.ix_(list(periphery_nodes),
-                                                 list(periphery_nodes))] > 0)
-    core_density = l_core / (len(core_nodes) * (len(core_nodes) - 1))
-    periphery_density = l_periphery / (
-            len(periphery_nodes) * (len(periphery_nodes) - 1))
-
-    if periphery_density > core_density:
-        cluster_assignment_new = 1 - cluster_assignment
-        return cluster_assignment_new
-    return cluster_assignment
-
-
-def flipping_function_cp(comm, n_flipping):
-    """Moves n nodes, randomly selected, from a partition to the other.
-
-    :param comm: Nodes memberships.
-    :type comm: numpy.ndarray
-    :param n_flipping: Number of nodes to flipped.
-    :type n_flipping: int
-    :return: New nodes memberships.
-    :rtype: numpy.ndarray
-    """
-    if np.random.random() > 0.5:
-        comm[np.random.choice(np.where(comm == 1)[0],
-                              replace=False,
-                              size=n_flipping)] = 0
-    else:
-        comm[np.random.choice(np.where(comm == 0)[0],
-                              replace=False,
-                              size=n_flipping)] = 1
-    return comm


Test patch:
diff --git a/tests/test_initialise_problem.py b/tests/test_initialise_problem.py
index b94a5a5d..acd8ddcc 100644
--- a/tests/test_initialise_problem.py
+++ b/tests/test_initialise_problem.py
@@ -132,7 +132,7 @@ def test_initial_guess_cd_divis_undir(self):
         graph = surp.UndirectedGraph(adjacency=adjacency_und)
 
         num_clusters = 3
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=num_clusters,
                                     initial_guess=None)
 
@@ -142,7 +142,7 @@ def test_initial_guess_cd_divis_dir(self):
         graph = surp.DirectedGraph(adjacency=adjacency_dir)
 
         num_clusters = 3
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=num_clusters,
                                     initial_guess=None)
 
@@ -152,7 +152,7 @@ def test_initial_guess_cd_divis_undir_common_neigh(self):
         graph = surp.UndirectedGraph(adjacency=adjacency_und)
 
         num_clusters = 3
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=num_clusters,
                                     initial_guess="common-neighbours")
 
@@ -162,7 +162,7 @@ def test_initial_guess_cd_divis_dir_common_neigh(self):
         graph = surp.DirectedGraph(adjacency=adjacency_dir)
 
         num_clusters = 3
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=num_clusters,
                                     initial_guess="common-neighbours")
 
@@ -174,7 +174,7 @@ def test_initial_guess_cd_divis_undir_custom(self):
         init_guess = np.array([k for k in np.arange(graph.n_nodes - 10)] +
                               [k for k in np.arange(10)])
 
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=2,
                                     initial_guess=init_guess)
 
@@ -186,7 +186,7 @@ def test_initial_guess_cd_divis_dir_custom(self):
         init_guess = np.array([k for k in np.arange(graph.n_nodes - 10)] +
                               [k for k in np.arange(10)])
 
-        graph._set_initial_guess_cd(method="divisive",
+        graph._set_initial_guess_cd(method="fixed-clusters",
                                     num_clusters=2,
                                     initial_guess=init_guess)