tau estim.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "X1 = np.random.multivariate_normal([0,0], 1e-3*np.diag(np.ones(2)),33)\n",
    "X2 = np.random.multivariate_normal([0.2,0.8], 1e-3*np.diag(np.ones(2)),33)\n",
    "X3 = np.random.multivariate_normal([0.8,0.5], 1e-3*np.diag(np.ones(2)),33)\n",
    "X = np.concatenate([X1,np.concatenate([X2,X3])])\n",
    "np.random.shuffle(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x114a6d390>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGGhJREFUeJzt3X1sXFeZx/Hfk+kEJlDq0BhpM0mbrDZNacmmpqZ01T+g\nBTYFltRK0relEkgVlVgVbSmylCyopAUpBQvYP7arpQjEire+RpbZBhmJBiF1SRVHbhslbdiQviQT\npJo27mqbKZ04z/7hGWc8vnfmznhe7DPfj1Rp5s7xzLk4/Ob4Oeeea+4uAEBYlnS6AwCA5iPcASBA\nhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAE6r1MfvGLFCl+zZk2nPh4AFqUDBw782d17\na7XrWLivWbNGY2Njnfp4AFiUzOzlJO0oywBAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDC\nHQAC1LGLmNBZw+M5DY0e0cnJvFb2ZDS4ab0G+rKd7haAJiHcu9DweE47dh9UvjAlScpN5rVj90FJ\nqhrwUV8IkviSABYgwr0LDY0emQn2knxhSkOjR2KDOeoL4a6Hn5nVJumXBIDWo+behU5O5us6LkV/\nIUQpfUkA6CzCvQtdkEnXdVyqHvzzaQugNQj3LmQWffx/3ypoeDwX+drKnkzi96/Wdng8p2vuf1Jr\ntz+ha+5/MvbzAMwP4d5FSsF66nQh8vWzLu3YfTAycAc3rVcmnar5GZl0amaiNerzd+w+qNxkXq5z\nNXoCHmg+JlS7ROWEaJx8YUo7Rw7FroD5l93P6XTh7KyfMUkuKVtjtUwjE7kAGkO4d4mkE6KSNJkv\naDI/PbovrYrZOXJI/7Dxr+SaXdMxSZ+9+iJ9c2BDzfdtZCIXQGMoy3SJ+QboZL6gn+17Zc4XhEva\n+8JEoveIq8XXU88HkAzh3iWaEaAeczzpF0dU3b5ajR5A4wj3LtHKAHUp0cqXgb6sdm3ZoGxPRqbp\nGv2uLRuotwMtYO5x47HW6u/vd26Q3V5X3PvrmVp6o0qTp1Ey6ZS2XpnVfz37p5nPWb4sra9/5nIC\nHGgSMzvg7v212iUauZvZ9WZ2xMyOmtn2iNcvMrO9ZjZuZs+Z2aca6TRaa+fmy+eURdJLTOlUzML3\nCpl0Sp+9+iJlY0o8+cKUfrrvlVlfIKdOFzT42LMsdwTarGa4m1lK0gOSPinpMkm3mtllFc2+JukR\nd++TdIukf292RzE/pU2/8oUppYpXMWV7Mhq6caOGtm2MDeySlJl2bdmgbw5sqLvEU5hytiQA2izJ\nUsirJB1192OSZGYPSbpB0uGyNi7pPcXHF0g62cxOYn6Gx3MafPRZFc5OF1Sm3JVeYnN2dazmOzdt\nnCmtNBLULHcE2itJuGclHS97fkLShyva7JT0azP7kqR3Sfp4U3qHptg5cmgm2EsKZ107dj8nyWqu\nf+/JpGfVzBsJ6iVmWrv9CbYFBtqkWatlbpX0Y3dfJelTkn5iZnPe28zuMLMxMxubmEi2NhrzFzeJ\nmi+crRnsmXRKOzdfPutYI8sqp9zZcgBooyThnpO0uuz5quKxcrdLekSS3P33kt4paUXlG7n7g+7e\n7+79vb29jfUYdWk0RKstVRzctF7JpmClJREN2RYYaL0k4b5f0jozW2tmSzU9YTpS0eYVSR+TJDN7\nv6bDnaF5h5X2k4kTFbzSdKi/eP+n9dT26zTQl52zk6M0veVAtYBfviytl+7/tOJW2lKDB1qrZs3d\n3c+Y2Z2SRiWlJP3I3Q+Z2X2Sxtx9RNJXJP3AzL6s6cnVz3unFtBjRrX9ZNIp080fWq3HD+RmtSld\nMVpaXZObzM9a214qq2y9MqueZenIHSYz6ZS+/pnpUs7KnoxyEUHOlgNAayXaOMzd90jaU3HsnrLH\nhyVd09yuYb6q3llp2/Tql/6L3xt5X9TyHSQrv6XzhSn9bN8rs47H7Qw5uGn9nN0o2XIAaD12hQxY\n3Kg525OZCd+Bvuycmvo19z9Zc6K1MvBLwf7U9utmHR/oy2rs5df1i6ePa8pdKTNtvXLuZwJoLvaW\nCVijG3U1Wg+P+rnh8ZweP5DTlJ9bY//4gRyrZYAWY+QesPKLjqJuvFGuVGM/OZnXErOZMI4St79M\nVB2dG3QAnUG4By6q7FKp8i5NUcFeXlO/9tLe2InYStygA+gMwh2xq2pSZjrrHjnij5qIjfoSYbUM\n0BmEO2JH0Wfd9eL9n458LclfBBKrZYBOYUIVLb39HTfoADqDkTtaPrpOOsoH0DyEO+paVVOufIUN\nuz0CCwvhDkn1j64rV9iUtiUovReAzqLmjoZUW78OoPMIdzSE9evAwka4oyGtXGEDYP4IdzSk0X1r\nALQHE6poSKMrbAC0B+GOhrF+HVi4KMsAQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0A\nAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgBKFu5ldb2ZHzOyomW2P\naXOTmR02s0Nm9vPmdhMAUI+at9kzs5SkByR9QtIJSfvNbMTdD5e1WSdph6Rr3P2Umb2vVR0GANSW\nZOR+laSj7n7M3d+W9JCkGyrafEHSA+5+SpLc/dXmdhMAUI8k4Z6VdLzs+YnisXKXSLrEzJ4ys31m\ndn3UG5nZHWY2ZmZjExMTjfUYAFBTsyZUz5O0TtJHJd0q6Qdm1lPZyN0fdPd+d+/v7e1t0kcDACol\nCfecpNVlz1cVj5U7IWnE3Qvu/qKkP2g67AEAHZAk3PdLWmdma81sqaRbJI1UtBnW9KhdZrZC02Wa\nY03sJwCgDjVXy7j7GTO7U9KopJSkH7n7ITO7T9KYu48UX/t7MzssaUrSoLu/1sqOA0AnDI/nNDR6\nRCcn81rZk9HgpvUa6Kuchuw8c/eOfHB/f7+PjY115LMBdI9mhvHweE47dh9UvjA1cyyTTmnXlg1t\nC3gzO+Du/bXacYUqgGCVwjg3mZdLyk3mtWP3QQ2PV04bJjM0emRWsEtSvjClodEjTehtcxHuAILV\n7DA+OZmv63gn1ay5A8BilSSM6ynbrOzJKBfxnit7Ms3pcBMxcgcQrLjQLR2vt2wzuGm9MunUrGOZ\ndEqDm9Y3td/NwIQqgGBFTYCWLF+WliSdOl2Y81q2J6Ontl8XOaqX1NHVMkknVCnLAAhWKXTv/eWh\nOSEeFeolJyfzc74YSqP6XVs26Knt17Wu001CWQZA0Ab6slq2tL5x7MqezKJaGROFcAcQvKhJ0Dil\nGvpiWhkThbIMgKANj+dkkuJmF3syab3rHefN1NCvvbRXQ6NHYtsvxJUxUQh3AEGrFtTplGnn5stn\navPVJmClc6P6xbAFAeEOIGjVyihD2zZKkq65/0mdnMxriZmmYlYQZstWy0RNtEpaUAFPuAMI1vB4\nLjaws8XySnlQxwW7pJkVMtfc/2TsROtCCncmVAEEqVRiiQrsUnklakVMFCu+n7R4tiBg5A5gQWv0\nQqJqwb31yqwG+rL68sPPJOqDF99voC+7aLYgINwBLFhRFxINPvqsZFJhymeORdW8q42kn3juT9r7\nwkTsRGuU0vsNblofue3vQtuCgLIMgAUravRdOOszwV4SdXFRtZH0qdOFuta+l7/fQF9Wu7ZsULYn\nI9N07b6d+7knxcgdwIJVTx27su3gpvW6K2HZpZbKkflAX3bBhXklwh3AghVX345rW26gL9uUcM9G\n1PTrWefeqTXxhDuABSuqvp1eYrNq7lJrat5xt8+L21BMmrvOvZ62zUbNHcCCFVXfHrpxo4a2baxZ\n8270VnolpRU1lerZUKyTm48xcgewoMXVt2uNfOcboHtfmIg8Xs86906uiWfkDiBI8w3QuFp/rbs7\nNdq22Qh3AEGab4CWX5Varp5b7XXytnyEO4AgRQVrPUpXpVaqZ517J9fEcw9VAMGqXIa4bOkS/c+r\nb9b1Hv968xULak0791AF0PXKJ2NLyxIrrXvfu3T01TdjtyJYiNv5JkFZBkBXiNtI7NjE6ap7zCym\n+6aWI9wBdIW41TPV9nCv9bMLGeEOoCvErZ5JmTX8swsZ4Q6gK8QtS7z1w6urrqpZiNv5JsGEKoCu\nUJoQjdrEq//i984c71mWlrv0Rr6wYG9+nQRLIQFgEUm6FDJRWcbMrjezI2Z21My2V2m31czczGp+\nMACgdWqGu5mlJD0g6ZOSLpN0q5ldFtHufEn/LOnpZncSAFCfJCP3qyQddfdj7v62pIck3RDR7huS\nviXprSb2DwDQgCThnpV0vOz5ieKxGWb2QUmr3f2JJvYNANCgeS+FNLMlkr4r6SsJ2t5hZmNmNjYx\nEb1XMgBg/pKEe07S6rLnq4rHSs6X9AFJvzWzlyRdLWkkalLV3R9093537+/t7W281wCAqpKE+35J\n68xsrZktlXSLpJHSi+7+hruvcPc17r5G0j5Jm92ddY4A0CE1w93dz0i6U9KopOclPeLuh8zsPjPb\n3OoOAgDql+gKVXffI2lPxbF7Ytp+dP7dAgDMB3vLAECACHcACBDhDgABItwBIECEOwAEiHAHgAAR\n7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEO\nAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQ\nIMIdAAJEuANAgBKFu5ldb2ZHzOyomW2PeP1uMztsZs+Z2W/M7OLmdxUAkFTNcDezlKQHJH1S0mWS\nbjWzyyqajUvqd/e/lfSYpG83u6MAgOSSjNyvknTU3Y+5+9uSHpJ0Q3kDd9/r7qeLT/dJWtXcbgIA\n6pEk3LOSjpc9P1E8Fud2Sb+KesHM7jCzMTMbm5iYSN5LAEBdmjqhama3SeqXNBT1urs/6O797t7f\n29vbzI8GAJQ5L0GbnKTVZc9XFY/NYmYfl/RVSR9x9780p3sAgEYkGbnvl7TOzNaa2VJJt0gaKW9g\nZn2Svi9ps7u/2vxuAgDqUTPc3f2MpDsljUp6XtIj7n7IzO4zs83FZkOS3i3pUTN7xsxGYt4OANAG\nScoycvc9kvZUHLun7PHHm9wvAMA8cIUqAASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwB\nIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQAC\nlOgG2d1keDynodEjOjmZ18qejAY3rddAX7bT3QKAujByLzM8ntOO3QeVm8zLJeUm87rr4Wd0xb2/\n1vB4rtPdA4DECPcyQ6NHlC9MzTk+mS9ox+6DBDyARWPRlmVaUT45OZmPfS1fmNLQ6JG6P4MyD4BO\nWJThXiqflEbZucm8duw+KEnzCs6VPRnlqgR8tfCPCnFJLeknANSyKMsyUeWT0sh6PgY3rVcmnYp9\nfWVPJvJ4VK1+x+6D2jlyqCX9BIBaFuXIPW4EXW1knURpNH3vLw/p1OnCrNcy6dTMaLxS3JdNVP2+\nGf0EgFoW5cg9bgQddzypUmll8nRBy5el1ZNJyyRlezLatWVDbCml3rCebz8BoJZFOXIf3LR+Vi1b\nqj6yjlNeJ78gk9abb59RYcolSadOF5RJp/S9m6+oWR+Pq9UvX5bWW4Wz8+4nANTL3L0jH9zf3+9j\nY2MN//x8V6FUTsrWsnxZWl//zOWRnxH1Xpl0Sru2bJAkVssAaBozO+Du/TXbLdZwr1fll8Hpt8/M\nqavXkk6ZhrZtjA14QhxAqxHuZeodpVeT7cnoqe3XNaFX5/DFACCppOGeaELVzK43syNmdtTMtke8\n/g4ze7j4+tNmtqb+LrdO3JWnjWj2Spe4ZZRcDQtgPmpOqJpZStIDkj4h6YSk/WY24u6Hy5rdLumU\nu/+Nmd0i6VuSbm5FhxvRzECOWumSZOQd16bamn1G7wAalWS1zFWSjrr7MUkys4ck3SCpPNxvkLSz\n+PgxSf9mZuadqvlUiFvNYiZF9bAnk9abfzmjwtnZL6ZTpsFN6/W14YP6xdPHNeUuM8kklZpGXYVa\n7YraVq3ZB9DdkpRlspKOlz0/UTwW2cbdz0h6Q9KFzehgM0RdeZpJpyKDXZreKOzd7zxPmfS5/3mW\nL0traNtGjb38un667xVNFX/Y/Vywl1RehVptdN6qNfsAultb17mb2R2S7pCkiy66qCWfEbfHi+lc\nAptJW6/Mau8LE7F7yZw6XZBJuu3qi/TNgQ0z71tt75ly5SPvaqPz7918RVPW7ANAuSThnpO0uuz5\nquKxqDYnzOw8SRdIeq3yjdz9QUkPStOrZRrpcDVR5Y/BR5/VlPus0bW79PD+47r5Q6v1+IFc7GSr\nS/rZvlf04sT/6b//+Lrq6XD5yDuuLLSyJzNTumG1DIBmShLu+yWtM7O1mg7xWyT9Y0WbEUmfk/R7\nSdskPdmJentU+aOybj5zfMq194UJ7dqyQXc/8syc0kqJS3rqj6/X1Y/KkXfUFbUm6dpLeyVN1+YJ\ncwDNVLPmXqyh3ylpVNLzkh5x90Nmdp+ZbS42+6GkC83sqKS7Jc1ZLtkO9U5CnpzMa6AvGxvs9Via\nsth9aAb6stp6ZVZW1t4lPX4gx5JHAC2RqObu7nsk7ak4dk/Z47ck3djcrtWv1n7sUe2bJbVkib69\nLX5zsb0vTMwp67DkEUCrLMpdIeNErYpJLzEtsbltS8sapemlj/NVa5/2uL8qcpN5Ru8Ami6ocB/o\ny2rXlg3K9mRmSiRDN27Ud2+6YlaAl5Y1lkbMOzdfrnTUN0CdqpWFqv2VwBWpAJqtK/aWKRd3pWjl\n8Wsv7Y1cSZNJL1G+cDbyvavtO1Nrf5tW7FkDIDxJ95ZZlPu5N6rWvVejat+lK1FTZrr1w6tj18ab\nVHVteum973r4mcjXuSIVQDMFVZap5d5fJr+n6fB4To8fyM1ciTrlrscP5GInbF21b3o90JdVlitS\nAbRB14T78Hgudv/2qFFz3JYBKYuuzceFdqW4rRC4IhVAM3VNWabaSpaoUXNcmWTKXZl0quHtArgi\nFUA7dE24V6tpRwVz3Jr5bDGM5xPOXJEKoNW6Jtzjwronk44M2mo34SacASx0XVNzj6t179x8eWT7\nqDXzldsKAMBC1TUj90Zq3YzQASxWXRPuEmENoHt0TVkGALoJ4Q4AASLcASBAhDsABIhwB4AAEe4A\nECDCHQACRLgDQIA6dicmM5uQ9HJHPrx+KyT9udOd6JBuPXfOu7sspvO+2N17azXqWLgvJmY2luS2\nViHq1nPnvLtLiOdNWQYAAkS4A0CACPdkHux0BzqoW8+d8+4uwZ03NXcACBAjdwAIEOFexsyuN7Mj\nZnbUzLZHvP4OM3u4+PrTZram/b1svgTnfbeZHTaz58zsN2Z2cSf62Qq1zr2s3VYzczMLYkVFkvM2\ns5uKv/dDZvbzdvexFRL8W7/IzPaa2Xjx3/unOtHPpnB3/psuTaUk/VHSX0taKulZSZdVtPknSf9R\nfHyLpIc73e82nfe1kpYVH38xhPNOeu7FdudL+p2kfZL6O93vNv3O10kal7S8+Px9ne53m877QUlf\nLD6+TNJLne53o/8xcj/nKklH3f2Yu78t6SFJN1S0uUHSfxYfPybpY2ZmbexjK9Q8b3ff6+6ni0/3\nSVrV5j62SpLfuSR9Q9K3JL3Vzs61UJLz/oKkB9z9lCS5+6tt7mMrJDlvl/Se4uMLJJ1sY/+ainA/\nJyvpeNnzE8VjkW3c/YykNyRd2JbetU6S8y53u6RftbRH7VPz3M3sg5JWu/sT7exYiyX5nV8i6RIz\ne8rM9pnZ9W3rXeskOe+dkm4zsxOS9kj6Unu61nxddQ9VzI+Z3SapX9JHOt2XdjCzJZK+K+nzHe5K\nJ5yn6dLMRzX9l9rvzGyDu092tFetd6ukH7v7d8zs7yT9xMw+4O5nO92xejFyPycnaXXZ81XFY5Ft\nzOw8Tf/Z9lpbetc6Sc5bZvZxSV+VtNnd/9KmvrVarXM/X9IHJP3WzF6SdLWkkQAmVZP8zk9IGnH3\ngru/KOkPmg77xSzJed8u6RFJcvffS3qnpvedWXQI93P2S1pnZmvNbKmmJ0xHKtqMSPpc8fE2SU96\nceZlEat53mbWJ+n7mg72EGqvJVXP3d3fcPcV7r7G3ddoer5hs7uPdaa7TZPk3/qwpkftMrMVmi7T\nHGtnJ1sgyXm/IuljkmRm79d0uE+0tZdNQrgXFWvod0oalfS8pEfc/ZCZ3Wdmm4vNfijpQjM7Kulu\nSbFL5xaLhOc9JOndkh41s2fMrPL/EItSwnMPTsLzHpX0mpkdlrRX0qC7L+q/UhOe91ckfcHMnpX0\nC0mfX6wDOK5QBYAAMXIHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABOj/AUTaN8wJ\ng+K5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115092650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.scatter(X[:,0],X[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from cvxpy import *\n",
    "from numba import jit\n",
    "from scipy.stats import multivariate_normal\n",
    "from sklearn.cluster import KMeans\n",
    "class EMtau(object):\n",
    "    # A EM-like algorithm that penalizes the columns of the posterior matrix tau\n",
    "    # resulting to a \"sparse\" estimation of the weights i.e. a estimation of the number of clusters\n",
    "    \n",
    "    def __init__(self, kmax=2, n_iter=10):\n",
    "        self.kmax = kmax\n",
    "        self.n_iter = n_iter\n",
    "        self.fista_iter = 200\n",
    "        self.lambd = 1\n",
    "        self.lambda_pen = 10\n",
    "    \n",
    "    def tau_gen(self, X):\n",
    "        densities = np.array([multivariate_normal(self.centers[k], self.covars[k]).pdf(X) for k in range(self.kmax)]).T*self.pi\n",
    "        s = densities.sum(axis=1)\n",
    "        return (densities.T/(densities.sum(axis=1))).T\n",
    "    \n",
    "    def covar_gen(self, X, i):\n",
    "        a = (X-self.centers[i])*(np.sqrt(self.tau[:,i]).reshape(-1,1))\n",
    "        return a.T.dot(a)/(self.N*self.pi[i])\n",
    "    @jit\n",
    "    def f_gradient(self, xi, X):\n",
    "        #give gradient of f on xi\n",
    "        temp = np.zeros([self.kmax, self.N])\n",
    "        for i in range(self.kmax):\n",
    "            temp[i] = -np.log(multivariate_normal(self.centers[i], self.covars[i]).pdf(X)) - np.log(self.pi[i]/xi[:,i]) +1 + self.lambda_pen*xi[:,i]/(np.linalg.norm(xi, axis=0)[i])\n",
    "        return temp.T\n",
    "\n",
    "    @jit()\n",
    "    def simplex_proj_numba(self, y):\n",
    "        dim = len(y)\n",
    "        u = np.flip(np.sort(y),0)\n",
    "        maxi = 0\n",
    "        lambd = 0\n",
    "        for i in range(dim):\n",
    "            crit = u[i]+1./(i+1)*(1-u[:i+1].sum())\n",
    "            if crit > 0 and i > maxi:\n",
    "                maxi = i\n",
    "        s = u[:maxi+1].sum()\n",
    "        lambd = 1./(maxi+1)*(1.-s)\n",
    "        res = np.zeros(dim)\n",
    "        for j in range(dim):\n",
    "            res[j] = max(y[j]+lambd, 0)\n",
    "        return res\n",
    "    @jit\n",
    "    def proj_simplex_row_matrix(self, Y):\n",
    "        res = np.zeros(Y.shape)\n",
    "        for i in range(Y.shape[0]):\n",
    "            res[i] = self.simplex_proj_numba(Y[i])\n",
    "        return res\n",
    "    #@jit\n",
    "    def expectation(self, X):\n",
    "        # Estimation of tau with penalization on columns\n",
    "        t_current = 1\n",
    "        xi = np.copy(self.tau)\n",
    "        tau_current = np.copy(xi)\n",
    "        tau_next = np.ones([self.N, self.kmax])\n",
    "        i=0\n",
    "        # The thresholds given are purely experimental\n",
    "        while  i < self.fista_iter:\n",
    "            tau_current = np.copy(tau_next)         \n",
    "            f_grad = self.f_gradient(xi, X)\n",
    "            xi_next = xi-self.lambd*f_grad\n",
    "            xi_next = np.nan_to_num(xi_next)\n",
    "            tau_next = self.proj_simplex_row_matrix(xi_next)         \n",
    "            t_next = (1+np.sqrt(1+4*t_current*2))/2\n",
    "            xi = tau_next + (t_current - 1)/t_next*(tau_next - tau_current)\n",
    "            # We set the negative values to 1e-20 to avoid errors with log in the gradient\n",
    "            #xi[xi <= 0] = 1e-20\n",
    "            t_current = t_next\n",
    "            #print np.linalg.norm(tau_next,axis=1)\n",
    "            i+=1\n",
    "        self.tau = tau_next\n",
    "    \n",
    "    def maximization(self, X):\n",
    "        self.pi = self.tau.sum(axis=0)/self.N\n",
    "        self.centers = np.array([(X*(self.tau[:,i].reshape(-1,1))).sum(axis=0)/(self.N*self.pi[i]) for i in range(self.kmax)])\n",
    "        self.covars = np.array([self.covar_gen(X, i) for i in range(self.kmax)])\n",
    "    \n",
    "    def fit(self, X):\n",
    "        #init\n",
    "        self.N = X.shape[0]\n",
    "        kmeans = KMeans(self.kmax)\n",
    "        kmeans.fit(X)\n",
    "        self.centers = kmeans.cluster_centers_\n",
    "        self.pi = np.array([1.0*len(kmeans.labels_[kmeans.labels_==i])/self.N for i in range(self.kmax)])\n",
    "        self.covars = np.array([1e-20*np.diag(np.ones(X.shape[1]))*np.cov((X[kmeans.labels_==i]- kmeans.cluster_centers_[i]).T) for i in range(self.kmax)])\n",
    "        self.tau = self.tau_gen(X)\n",
    "        \n",
    "        #algorithm starts\n",
    "        for i in range(self.n_iter):\n",
    "            print \"iter: \", i\n",
    "            #self.tau[self.tau<=0] = 1e-20\n",
    "            self.expectation(X)\n",
    "            self.maximization(X)\n",
    "        return self.pi, self.centers, self.covars\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "emtau = EMtau(kmax=5, n_iter=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter:  0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:65: RuntimeWarning: divide by zero encountered in log\n",
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:65: RuntimeWarning: invalid value encountered in divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter:  1\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "array must not contain infs or NaNs",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-79-7a5a82c9c7df>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0memtau\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-77-6aad125c4d14>\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m     95\u001b[0m             \u001b[0;32mprint\u001b[0m \u001b[0;34m\"iter: \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     96\u001b[0m             \u001b[0;31m#self.tau[self.tau<=0] = 1e-20\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 97\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpectation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     98\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaximization\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     99\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcenters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcovars\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-77-6aad125c4d14>\u001b[0m in \u001b[0;36mexpectation\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m     63\u001b[0m         \u001b[0;32mwhile\u001b[0m  \u001b[0mi\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfista_iter\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     64\u001b[0m             \u001b[0mtau_current\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau_next\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 65\u001b[0;31m             \u001b[0mf_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_gradient\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     66\u001b[0m             \u001b[0mxi_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlambd\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mf_grad\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m             \u001b[0mxi_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnan_to_num\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxi_next\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, mean, cov, allow_singular, seed)\u001b[0m\n\u001b[1;32m    347\u001b[0m         return multivariate_normal_frozen(mean, cov,\n\u001b[1;32m    348\u001b[0m                                           \u001b[0mallow_singular\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mallow_singular\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 349\u001b[0;31m                                           seed=seed)\n\u001b[0m\u001b[1;32m    350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    351\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_process_parameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdim\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmean\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcov\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, mean, cov, allow_singular, seed)\u001b[0m\n\u001b[1;32m    597\u001b[0m         self.dim, self.mean, self.cov = self._dist._process_parameters(\n\u001b[1;32m    598\u001b[0m                                                             None, mean, cov)\n\u001b[0;32m--> 599\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcov_info\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_PSD\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcov\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallow_singular\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mallow_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    600\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    601\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mlogpdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, M, cond, rcond, lower, check_finite, allow_singular)\u001b[0m\n\u001b[1;32m    146\u001b[0m         \u001b[0;31m# Note that eigh takes care of array conversion, chkfinite,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    147\u001b[0m         \u001b[0;31m# and assertion that the matrix is square.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 148\u001b[0;31m         \u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meigh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlower\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlower\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcheck_finite\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_finite\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    149\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    150\u001b[0m         \u001b[0meps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_eigvalsh_to_eps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrcond\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/linalg/decomp.pyc\u001b[0m in \u001b[0;36meigh\u001b[0;34m(a, b, lower, eigvals_only, overwrite_a, overwrite_b, turbo, eigvals, type, check_finite)\u001b[0m\n\u001b[1;32m    286\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m     \"\"\"\n\u001b[0;32m--> 288\u001b[0;31m     \u001b[0ma1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_asarray_validated\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcheck_finite\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_finite\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    289\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0ma1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0ma1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    290\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'expected square matrix'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/_lib/_util.pyc\u001b[0m in \u001b[0;36m_asarray_validated\u001b[0;34m(a, check_finite, sparse_ok, objects_ok, mask_ok, as_inexact)\u001b[0m\n\u001b[1;32m    226\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'masked arrays are not supported'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    227\u001b[0m     \u001b[0mtoarray\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray_chkfinite\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcheck_finite\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 228\u001b[0;31m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtoarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    229\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mobjects_ok\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    230\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'O'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/numpy/lib/function_base.pyc\u001b[0m in \u001b[0;36masarray_chkfinite\u001b[0;34m(a, dtype, order)\u001b[0m\n\u001b[1;32m   1213\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchar\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtypecodes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'AllFloat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1214\u001b[0m         raise ValueError(\n\u001b[0;32m-> 1215\u001b[0;31m             \"array must not contain infs or NaNs\")\n\u001b[0m\u001b[1;32m   1216\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1217\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: array must not contain infs or NaNs"
     ]
    }
   ],
   "source": [
    "emtau.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "emtau.tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "gmm = GaussianMixture(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussianMixture(covariance_type='full', init_params='kmeans', max_iter=100,\n",
       "        means_init=None, n_components=5, n_init=1, precisions_init=None,\n",
       "        random_state=None, reg_covar=1e-06, tol=0.001, verbose=0,\n",
       "        verbose_interval=10, warm_start=False, weights_init=None)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.09076937,  0.1276379 ,  0.33333333,  0.24256396,  0.20569543])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm.weights_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00871185, -0.00181355],\n",
       "       [ 0.19765643,  0.75945015],\n",
       "       [ 0.84238601,  0.57172857],\n",
       "       [ 0.761032  ,  0.44362419],\n",
       "       [ 0.18991157,  0.85071262]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm.means_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tau = [[Variable() for _ in range(10) ] for _ in range(20)]\n",
    "constraints = [sum_entries(bmat(ligne_n)) == 1 for ligne_n in tau ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Expression(AFFINE, UNKNOWN, (1, 1))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum_entries(bmat(tau[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Expression(AFFINE, UNKNOWN, (20, 10))"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bmat(tau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CVXOPT', 'MOSEK', 'ECOS_BB', 'SCS', 'ECOS', 'LS']\n"
     ]
    }
   ],
   "source": [
    "print installed_solvers()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "    \n",
    "    def expectation(self, X):\n",
    "        # Estimation of tau with penalization on columns\n",
    "        # We use CVXPY, a more efficient implementation would be to solve directly the optimization problem\n",
    "        # with the appropriate procedure\n",
    "        t_current = 1\n",
    "        xi = super(EMtau, self).tau_gen(X)\n",
    "        tau_current = np.copy(xi)\n",
    "        tau_next = np.ones([self.N, self.kmax])\n",
    "        i=0\n",
    "        # The thresholds given are purely experimental\n",
    "        while np.linalg.norm(tau_current-tau_next) > 1e-5 and i < self.fista_iter:\n",
    "            tau_current = np.copy(tau_next)\n",
    "            #tau = [[Variable() for _ in range(self.kmax) ] for _ in range(self.N)]\n",
    "            #constraints = [\n",
    "                #sum_entries(bmat(ligne_n)) == 1 for ligne_n in tau \n",
    "                 #         ]+[\n",
    "             #   item >=0 for sublist in tau for item in sublist \n",
    "                #        ]+[\n",
    "               # item <=1 for sublist in tau for item in sublist \n",
    "             #   ]             \n",
    "            f_grad = self.f_gradient(xi, X)\n",
    "            #import pdb; pdb.set_trace()\n",
    "            xi_next = xi-self.lambd*f_grad\n",
    "            tau_next = self.proj_simplex_row_matrix(xi_next)\n",
    "            #prob = Problem(Minimize(norm(bmat(tau) - xi_next,\"fro\" )**2), constraints)\n",
    "            #prob.solve(solver=\"SCS\")\n",
    "            #tau_next_temp = np.array(bmat(tau).value)\n",
    "            \n",
    "            t_next = (1+np.sqrt(1+4*t_current*2))/2\n",
    "            xi = tau_next + (t_current - 1)/t_next*(tau_next - tau_current)\n",
    "            # We set the negative values to 1e-20 to avoid errors with log in the gradient\n",
    "            xi[xi <= 0] = 1e-20\n",
    "            t_current = t_next\n",
    "            #print np.linalg.norm(tau_next,axis=1)\n",
    "            i+=1\n",
    "        print np.linalg.norm(tau_next, axis=1)\n",
    "        return tau_next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5],\n",
       "       [ 0.5,  0.5]])"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/(2*np.ones([10,2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e-20,   0.00000000e+00],\n",
       "       [  0.00000000e+00,   1.00000000e-20]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1e-20*np.diag(np.ones(X.shape[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}