HyperMap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Equations\n",
    "\n",
    "equations: $\\alpha-1+(\\frac{\\kappa_0}{\\kappa_C,})^{\\gamma-2}\\\\\\\n",
    "P_0- (\\gamma-1)(\\alpha\\kappa_0)^{(\\gamma-1)}\\Gamma(1-\\gamma,\\alpha\\kappa_0)\\\\\\\n",
    "k_{max_obs}-\\alpha\\kappa_C\\\\\\\n",
    "\\bar{k}\\alpha^2-(1-P_0)k_{bar_obs}\\\\\\\n",
    "\\kappa_0(\\gamma-1)-\\bar{k}(\\gamma-2))\\\\\\\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "import mpmath\n",
    "\n",
    "\n",
    "def solve(gamma,k_bar_obs,k_max_obs) :\n",
    "    def equations(p):\n",
    "        alpha, kappa0,kappaC,P0,k_bar = p\n",
    "        return (alpha-1+(kappa0/kappaC)**(gamma-2),\n",
    "            P0- (gamma-1)*(alpha*kappa0)**(gamma-1)*mpmath.gammainc(1-gamma,alpha*kappa0),\n",
    "            k_max_obs-alpha*kappaC,\n",
    "            k_bar*alpha**2-(1-P0)*k_bar_obs,\n",
    "            kappa0*(gamma-1)-k_bar*(gamma-2))\n",
    "    return  fsolve(equations, (1, 1,k_max_obs,0.5,k_bar_obs))\n",
    "\n",
    "def verify(alpha, kappa0,kappaC,P0,k_bar,gamma,k_bar_obs,k_max_obs) :\n",
    "    def equations(p):\n",
    "        alpha, kappa0,kappaC,P0,k_bar = p\n",
    "        return (alpha-1+(kappa0/kappaC)**(gamma-2),\n",
    "            P0- (gamma-1)*(alpha*kappa0)**(gamma-1)*mpmath.gammainc(1-gamma,alpha*kappa0),\n",
    "            k_max_obs-alpha*kappaC,\n",
    "            k_bar*alpha**2-(1-P0)*k_bar_obs,\n",
    "            kappa0*(gamma-1)-k_bar*(gamma-2))\n",
    "    print(equations((alpha, kappa0,kappaC,P0,k_bar)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.66840165, -2.        ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import fsolve\n",
    "import mpmath\n",
    "def eq(t):\n",
    "    x,y=t\n",
    "    return x**3-19,y+2\n",
    "fsolve(eq,(-1,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- exponent\n",
    "\n",
    "$1+\\frac{\\text{节点数}}{\\sum \\ln(\\text{每个节点的度数})}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import powerlaw\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def get_exponent(graphe):\n",
    "    return 1+(graphe.number_of_nodes()/sum(np.log(graphe.degree().values())))\n",
    "\n",
    "def plot_gamma(graphe) :\n",
    "    #fitness = powerlaw.Fit(graphe.degree().values(),xmin =1)\n",
    "    #print fitness.alpha,fitness.xmin\n",
    "    alpha = 1+(graphe.number_of_nodes()/sum(np.log(graphe.degree().values())))\n",
    "\n",
    "    #powerlaw.plot_pdf(graphe.degree().values())\n",
    "    xdata,ydata = powerlaw.pdf(graphe.degree().values())\n",
    "\n",
    "    vala ={}\n",
    "    xvala=[]\n",
    "    yvala=[]\n",
    "    somme = 0\n",
    "    exceptions =[]\n",
    "    last_exception = -1\n",
    "    for value in graphe.degree().values() :\n",
    "        somme+=1\n",
    "        try :\n",
    "            vala[value]+=1\n",
    "        except :\n",
    "            vala[value] =1\n",
    "\n",
    "    for degree in range(1,max(vala.keys())+1) :\n",
    "        try :\n",
    "            yvala.append(float(vala[degree])/somme)\n",
    "            xvala.append(degree)\n",
    "\n",
    "        except :\n",
    "            if len(exceptions) ==0 :\n",
    "                last_exception =degree\n",
    "                #print degree\n",
    "            exceptions.append(degree)\n",
    "    #slope, intercept, r_value, p_value, std_err = stats.linregress(\n",
    "        #np.log(xvala[:last_exception]),np.log(yvala[:last_exception]))\n",
    "\n",
    "    #print slope\n",
    "\n",
    "    x = xvala\n",
    "\n",
    "    plt.plot(x,np.power(x,-alpha),color = \"orange\")\n",
    "    #plt.plot(x,np.power(x,-2.34296138661),color =\"blue\")\n",
    "    #plfit.plfit(graphe.degree().values(),quiet = False,verbose = True,silent = False)\n",
    "    plt.plot(xvala,yvala,color = \"black\")\n",
    "    plt.yscale('log')\n",
    "    plt.xscale('log')\n",
    "    #plt.ylim([0,0.4])\n",
    "    #plt.xlim([0,20])\n",
    "    plt.show()\n",
    "    #xdata,ydata=powerlaw.pdf(graphe.degree().values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 生成随机图\n",
    "    - $\\mu=\\beta \\times \\frac{\\sin(\\frac{\\pi}{\\beta})}{2\\pi\\bar{k}}$；\n",
    "    - $\\kappa $为$PowerLaw(x_\\text{min}=\\kappa_0,\\gamma)$生成的N个点；\n",
    "    - $\\theta$是在$(0,2\\pi)$均匀分布的$N$个点；\n",
    "    - th: N个$\\theta$的array；\n",
    "    - r 代表半径，i,j\\j,i只计算一次。保存在r2中；\n",
    "    - $khi=N\\times d_{th}/(2\\pi\\times\\mu\\kappa\\kappa_2)$\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import powerlaw\n",
    "import numpy as np\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drange(start, stop, step): # 连续抛出几个数\n",
    "    r = start\n",
    "    while r < stop:\n",
    "        yield r\n",
    "        r += step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clust(N, k_bar, kappa0, gamma, start=1.1, stop=10, step=0.1, nb_networks=100):\n",
    "    def mapping_function(beta):\n",
    "        mean, std = clust_beta_given(beta, N, k_bar, kappa0, gamma, nb_networks)\n",
    "        return [beta, mean, std]\n",
    "\n",
    "    return map(mapping_function, drange(start, stop, step))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clust_beta_given(beta, N, k_bar, kappa0, gamma, nb_networks):\n",
    "    def mapping_function(_):\n",
    "        return get_clustering_from_random_network(beta, N, k_bar, kappa0, gamma)\n",
    "\n",
    "    result = map(mapping_function, range(nb_networks))\n",
    "    return np.mean(result), np.std(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approximated_beta(N, k_bar, kappa0, gamma, c_real, start=1.01, step=0.1, nb_networks=100):\n",
    "    path_followed = []\n",
    "\n",
    "    def rec_approximated_beta(beta_min, beta, beta_max):\n",
    "        if beta_max is not None and (beta_max - beta_min) < step:\n",
    "            return beta\n",
    "        c_mean, c_std = clust_beta_given(beta, N, k_bar, kappa0, gamma, nb_networks)\n",
    "        path_followed.append([beta, c_mean, c_std])\n",
    "        print(beta, c_mean, c_std)\n",
    "        if c_mean > c_real:\n",
    "            return rec_approximated_beta(beta_min, (beta + beta_min) / 2, beta)\n",
    "        if beta_max is not None:\n",
    "            return rec_approximated_beta(beta, (beta + beta_max) / 2, beta_max)\n",
    "        return rec_approximated_beta(beta, 2 * beta, None)\n",
    "\n",
    "    beta_optimal = rec_approximated_beta(start, 1 / (1 - c_real), None)\n",
    "    return beta_optimal, path_followed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_clustering_from_random_network(beta, N, k_bar, kappa0, gamma):\n",
    "    \"\"\"\n",
    "    :rtype: float\n",
    "    \"\"\"\n",
    "    graphe = generate_random_network(beta, N, k_bar, kappa0, gamma)# 下方的随机图\n",
    "    if all(v == 0 for v in nx.clustering(graphe).values()):#聚集系数的值\n",
    "        return 0\n",
    "    return nx.average_clustering(graphe, count_zeros=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_random_network(beta, N, k_bar, kappa0, gamma):\n",
    "    mu = beta * math.sin(math.pi / beta) / (2 * math.pi * k_bar)\n",
    "    graphe_random = nx.Graph()\n",
    "    kappa = powerlaw.Power_Law(xmin=kappa0, parameters=[gamma]).generate_random(N)# 按照参数生成一个大小为N的array\n",
    "    theta = np.random.uniform(0, 2 * math.pi, N)\n",
    "    # print 1+(len(kappa)/sum(np.log(kappa))) #1.6\n",
    "    th = np.outer(np.ones(N), theta)#生成N个theta\n",
    "    # print \"th\"\n",
    "    th2 = np.outer(theta, np.ones(N))\n",
    "    # print \"th2\"\n",
    "    kp = np.outer(np.ones(N), kappa)\n",
    "    # print \"kp\"\n",
    "    kp2 = np.outer(kappa, np.ones(N))\n",
    "    # print \"kp2\"\n",
    "    dth = math.pi - np.abs(math.pi - np.abs(th - th2))#角距离\n",
    "    # print \"dth\"\n",
    "    khi = N * dth / (2 * math.pi * mu * kp * kp2)\n",
    "    # print \"khi\"\n",
    "    r = np.random.uniform(0, 1, [N, N])\n",
    "    # print \"r\"\n",
    "    r2 = np.triu(np.log(1 / r - 1) / beta, k=1)#提取不含对角线的上三角矩阵\n",
    "    p2 = np.triu(np.log(khi), k=1)\n",
    "    # print np.where(r2>p2) ie ln(1/r-1) / beta > ln(khi) <=> r < 1/(1+khi^beta)\n",
    "    graphe_random.add_edges_from(np.transpose(np.where(r2 > p2)))\n",
    "    return graphe_random\n",
    "    # print graphe_random.edges()\n",
    "    # print graphe_random.number_of_nodes()\n",
    "    # print graphe_random.number_of_edges()\n",
    "    # print exponent.get_exponent(graphe_random)\n",
    "    # print graphe_random.degree().values()\n",
    "    # print np.mean(graphe_random.degree().values())\n",
    "    # print max(graphe_random.degree().values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 似然函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# coding=utf-8\n",
    "import networkx as nx\n",
    "import math\n",
    "import numpy as np\n",
    "import operator\n",
    "\n",
    "\n",
    "def optimisation_likelihood(graphe, constantes, theta_init=None, subgroup_of_nodes=None):\n",
    "    \"\"\"On maximise la likelihood, en la maximisant pour chaque i separement\n",
    "    theta_init est la configuration initiale de l'optimisation\n",
    "    subgroup_of_nodes est le sous groupe pour lesquels on maximise la likelihood, les autres restent inchanges\"\"\"\n",
    "    kappa0, gamma, beta, Nobs, k_bar = constantes\n",
    "\n",
    "    kappa = {node: max(kappa0, degree - gamma / beta) for node, degree in graphe.degree().items()}\n",
    "    mu = beta * math.sin(math.pi / beta) / (2 * math.pi * k_bar)\n",
    "    a = (nx.to_numpy_matrix(graphe)).A\n",
    "\n",
    "    if theta_init is None:\n",
    "        theta_current = 2 * math.pi * np.random.randint(0, Nobs, Nobs) / Nobs\n",
    "    else:\n",
    "        theta_current = theta_init\n",
    "\n",
    "    def localloglikelihood(theta_angle, node):\n",
    "        \"\"\"maximisation de la likelihood sur node\"\"\"\n",
    "        th2 = theta_current\n",
    "        dth = math.pi - np.abs(math.pi - np.abs(theta_angle - th2))\n",
    "        khi = Nobs * dth / 2 * math.pi * mu * kappa[node] * kappa.values()\n",
    "        p = 1 / (1 + khi ** beta)\n",
    "\n",
    "        lp = np.log(p)\n",
    "        lp[a[node] == 0] = 0\n",
    "        l_p = np.log(1 - p)\n",
    "        l_p[a[node] == 1] = 0\n",
    "        alp = lp + l_p\n",
    "        return np.sum(alp)\n",
    "\n",
    "    theta_test = np.linspace(0, 2 * math.pi, num=Nobs)\n",
    "\n",
    "    def maximize_locally(node):\n",
    "        \"\"\"\"\"\"\n",
    "        lll_vector = map(lambda theta_angle: localloglikelihood(theta_angle, node), theta_test)\n",
    "        theta_optimal = theta_test[np.argmax(lll_vector)]\n",
    "        return theta_optimal\n",
    "\n",
    "    if subgroup_of_nodes is None: subgroup_of_nodes = graphe.nodes()\n",
    "\n",
    "    best_theta = theta_init\n",
    "    best_likelihood = verifyloglikelihood(theta_init, graphe, constantes, silent=False)\n",
    "    nb_retours = 0\n",
    "    while (nb_retours < 2):\n",
    "        for node, degree in sorted(graphe.degree_iter(subgroup_of_nodes), key=operator.itemgetter(1), reverse=True):\n",
    "            theta_new = maximize_locally(node)\n",
    "            if degree == 2:\n",
    "                neighbor, neighbor_degree = max(graphe.degree_iter(graphe.neighbors(node)), key=operator.itemgetter(1))\n",
    "                theta_new2 = theta_current[neighbor]\n",
    "                print(2, theta_new, theta_new2)\n",
    "            if degree == 1:\n",
    "                theta_new1 = theta_current[graphe.neighbors(node)[0]]\n",
    "                print(1, theta_new, theta_new1)\n",
    "\n",
    "            theta_current[node] = theta_new\n",
    "        current_likelihood = verifyloglikelihood(theta_current, graphe, constantes, silent=False)\n",
    "        if current_likelihood < best_likelihood:\n",
    "            nb_retours += 1\n",
    "        else:\n",
    "            best_likelihood = current_likelihood\n",
    "            best_theta = theta_current\n",
    "            nb_retours = 0\n",
    "    return best_theta\n",
    "\n",
    "'''\n",
    "def optimisation_complete_par_etapes_likelihood(graphe, bornes, theta0, constantes):\n",
    "    \"\"\"On maximise la likelihood en commencant par des sous graphes contenant les noeuds principaux,\n",
    "    on optimise sur le sous graphe complet à chaque étape\"\"\"\n",
    "    kappa0, gamma, beta, Nobs, k_bar = constantes\n",
    "\n",
    "    def optimisation_etape_likelihood(nodes_subgraph, theta_global):\n",
    "        \"\"\"On maximise la likelihood sur le sous graphe\"\"\"\n",
    "        mapping_old_new = {}\n",
    "        mapping_new_old = {}\n",
    "        for i in range(len(nodes_subgraph)):\n",
    "            old = nodes_subgraph[i]\n",
    "            mapping_old_new[old] = i\n",
    "            mapping_new_old[i] = old\n",
    "        graphe_intermediaire = graphe.subgraph(nodes_subgraph)\n",
    "        graphe_intermediaire = nx.relabel_nodes(graphe_intermediaire, mapping_old_new)\n",
    "        Nlocal = len(nodes_subgraph)\n",
    "        local_variables = kappa0, gamma, beta, Nlocal, k_bar\n",
    "\n",
    "        theta_init = np.asarray([theta_global[mapping_new_old[node]] for node in graphe_intermediaire])\n",
    "        theta_local_opt = optimisation_likelihood(graphe_intermediaire, local_variables, theta_init=theta_init)\n",
    "\n",
    "        theta_global_opt = theta_global\n",
    "\n",
    "        for node in graphe_intermediaire:\n",
    "            theta_global_opt[mapping_new_old[node]] = theta_local_opt[node]\n",
    "        return theta_global_opt\n",
    "\n",
    "    theta_current = theta0\n",
    "    for borne in bornes:\n",
    "        nodes_subgraph = [node for node, degree in\n",
    "                          sorted(graphe.degree_iter(), key=operator.itemgetter(1), reverse=True) if\n",
    "                          degree >= borne]\n",
    "\n",
    "        theta_current = optimisation_etape_likelihood(nodes_subgraph, theta_current)\n",
    "\n",
    "    return theta_current\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subgraph_mapping(graphe, nodes_subgraph):\n",
    "    mapping_old_new = {}\n",
    "    mapping_new_old = {}\n",
    "    for i in range(len(nodes_subgraph)):\n",
    "        old = nodes_subgraph[i]\n",
    "        mapping_old_new[old] = i\n",
    "        mapping_new_old[i] = old\n",
    "        # new network\n",
    "    graphe_intermediaire = graphe.subgraph(nodes_subgraph)\n",
    "    graphe_intermediaire = nx.relabel_nodes(graphe_intermediaire, mapping_old_new)\n",
    "    return graphe_intermediaire, mapping_old_new, mapping_new_old"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimisation_par_etapes_likelihood(graphe, bornes, borne_lim, theta0, constantes, method=\"degree\"):\n",
    "    \"\"\"On maximise la likelihood en commencant par des sous graphes contenant les noeuds principaux,\n",
    "    on optimise uniquement sur els nouveaux noeuds à chaque étape\"\"\"\n",
    "    kappa0, gamma, beta, Nobs, k_bar = constantes\n",
    "\n",
    "    def optimisation_etape_likelihood(nodes_subgraph, nodes_evaluated, theta_global):\n",
    "        # rename nodes of new network\n",
    "        graphe_intermediaire, mapping_old_new, mapping_new_old = subgraph_mapping(graphe, nodes_subgraph)\n",
    "        # new \"constants\"\n",
    "        nodes_evaluated_local = map(lambda node: mapping_old_new[node], nodes_evaluated)\n",
    "        Nlocal = len(nodes_subgraph)\n",
    "        local_variables = kappa0, gamma, beta, Nlocal, k_bar\n",
    "        # optimize on subgraph\n",
    "        theta_init = np.asarray([theta_global[mapping_new_old[node]] for node in graphe_intermediaire])\n",
    "        theta_local_opt = optimisation_likelihood(graphe_intermediaire, local_variables,\n",
    "                                                  theta_init=theta_init,\n",
    "                                                  subgroup_of_nodes=nodes_evaluated_local)\n",
    "        # update theta values\n",
    "        theta_global_opt = theta_global\n",
    "\n",
    "        for node in graphe_intermediaire:\n",
    "            theta_global_opt[mapping_new_old[node]] = theta_local_opt[node]\n",
    "        return theta_global_opt\n",
    "    '''\n",
    "    if method == \"degree\":\n",
    "        theta_current = theta0\n",
    "        borne_sup = float(\"inf\")\n",
    "        degree_dict = graphe.degree_iter()\n",
    "        for borne_inf in bornes:\n",
    "            nodes_subgraph = [node for node, degree in degree_dict if degree >= borne_inf]\n",
    "            if borne_inf > borne_lim:\n",
    "                \"\"\"on optimise sur tous les noeuds du sous graphe\"\"\"\n",
    "                theta_current = optimisation_etape_likelihood(nodes_subgraph, nodes_subgraph, theta_current)\n",
    "            else:\n",
    "                \"\"\"on optimise uniquement sur les noeuds qui viennent juste d'etre ajoutés\"\"\"\n",
    "                nodes_evaluated = [node for node, degree in degree_dict if borne_sup > degree >= borne_inf]\n",
    "                theta_current = optimisation_etape_likelihood(nodes_subgraph, nodes_evaluated, theta_current)\n",
    "            borne_sup = borne_inf\n",
    "    '''\n",
    "    if method == \"core_number\":\n",
    "        theta_current = theta0\n",
    "\n",
    "        core_dict = nx.core_number(graphe)\n",
    "        bornes = sorted(set(core_dict.values()), reverse=True)\n",
    "        borne_sup = float(\"inf\")\n",
    "        theta_precedent_centre = [theta_current[node] for node, core_num in core_dict.iteritems() if core_num == 27]\n",
    "        theta_courant_centre = [theta_current[node] for node, core_num in core_dict.iteritems() if core_num == 27]\n",
    "        for borne_inf in bornes:\n",
    "            nodes_subgraph = [node for node, core_number in core_dict.iteritems() if core_number >= borne_inf]\n",
    "            if borne_inf > borne_lim:\n",
    "                \"\"\"on optimise sur tous les noeuds du sous graphe\"\"\"\n",
    "                theta_current = optimisation_etape_likelihood(nodes_subgraph, nodes_subgraph, theta_current)\n",
    "                theta_courant_centre = [theta_current[node] for node, core_num in core_dict.iteritems() if\n",
    "                                        core_num == 27]\n",
    "                print(borne_inf, np.mean(np.abs(\n",
    "                    math.pi - np.absolute(np.asarray(theta_courant_centre) - np.asarray(theta_precedent_centre)))), \\\n",
    "                    np.max(np.abs(\n",
    "                        math.pi - np.absolute(np.asarray(theta_courant_centre) - np.asarray(theta_precedent_centre)))))\n",
    "                theta_precedent_centre = theta_courant_centre\n",
    "            else:\n",
    "                \"\"\"on optimise uniquement sur les noeuds qui viennent juste d'etre ajoutés\"\"\"\n",
    "                nodes_evaluated = [node for node, core_number in core_dict.iteritems() if\n",
    "                                   borne_sup > core_number >= borne_inf]\n",
    "                theta_current = optimisation_etape_likelihood(nodes_subgraph, nodes_evaluated, theta_current)\n",
    "            borne_sup = borne_inf\n",
    "\n",
    "    return theta_current\n",
    "    \"\"\"On maximise la likelihood en commencant par des sous graphes contenant les noeuds principaux,\n",
    "    on optimise sur le sous graphe complet à chaque étape\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verifyloglikelihood(theta_vector, graphe, constantes, silent=False):\n",
    "    \"\"\"On calcule somme (aij * logpij + (1-aij)* log(1-pij))\"\"\"\n",
    "    # aij\n",
    "    a = (nx.to_numpy_matrix(graphe)).A\n",
    "\n",
    "    # pij\n",
    "    kappa0, gamma, beta, Nobs, k_bar = constantes\n",
    "    mu = beta * math.sin(math.pi / beta) / (2 * math.pi * k_bar)\n",
    "    kappa = {node: max(kappa0, degree - gamma / beta) for node, degree in graphe.degree().items()}\n",
    "    kp = np.outer(np.ones(Nobs), kappa.values())\n",
    "    kp2 = np.outer(kappa.values(), np.ones(Nobs))\n",
    "    mukpkp2 = 2 * math.pi * mu * kp * kp2\n",
    "    th = np.outer(np.ones(Nobs), theta_vector)\n",
    "    th2 = np.outer(theta_vector, np.ones(Nobs))\n",
    "    dth = math.pi - np.abs(math.pi - np.abs(th - th2))\n",
    "    khi = Nobs * dth / mukpkp2\n",
    "    p = 1 / (1 + khi ** beta)\n",
    "\n",
    "    # log(pij)\n",
    "    lp = np.triu(np.log(p), k=1)\n",
    "    l_p = np.triu(np.log(1 - p), k=1)\n",
    "\n",
    "    # aij*log(pij)\n",
    "    lp[a == 0] = 0\n",
    "    l_p[a == 1] = 0\n",
    "    if not silent: print(np.sum(lp + l_p), np.sum(lp), np.sum(l_p))\n",
    "    return np.sum(lp + l_p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rayon(graphe, constantes):\n",
    "    kappa0, gamma, beta, Nobs, k_bar = constantes\n",
    "    kappa = np.asarray([max(kappa0, degree - gamma / beta) for node, degree in graphe.degree().items()])\n",
    "    mu = beta * math.sin(math.pi / beta) / (2 * math.pi * k_bar)\n",
    "    rayon = 2 * np.log(Nobs / (math.pi * mu * kappa * kappa0))\n",
    "    return rayon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}