Verify Laplacian .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from imports import *\n",
    "\n",
    "from sklearn.manifold import TSNE\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.1/2)**(2*0+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 30/30 [00:14<00:00,  1.96it/s]\n"
     ]
    }
   ],
   "source": [
    "from math import factorial as fact\n",
    "\n",
    "def modified_bessel_term(m,v,u):\n",
    "    return (1/(fact(m)*fact(m+v)))*((u/2)**(2*m+v))\n",
    "\n",
    "def modified_bessel_first_kind_approx(v,u,n_terms=100):\n",
    "    assert isinstance(v, int)\n",
    "    tot=0\n",
    "    for m in range(0,n_terms):\n",
    "        tot += modified_bessel_term(m,v,u)\n",
    "    return tot\n",
    "modified_bessel_first_kind_approx(1,0.01,n_terms=2000)\n",
    "\n",
    "res = []\n",
    "for v in tqdm(range(1,31)):\n",
    "    for u in [0.0000000000000001,0.0000001,0.00001,0.0001,0.001,0.01,0.1]:\n",
    "        for n_terms in [1,2,3,4,5,6,10,100,1000]:\n",
    "            mod = modified_bessel_first_kind_approx(v,u,n_terms)\n",
    "            res += [[v,u,n_terms,mod]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>v</th>\n",
       "      <th>u</th>\n",
       "      <th>n_terms</th>\n",
       "      <th>out</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>1</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>2</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>3</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>4</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>5</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>6</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>10</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>100</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-16</td>\n",
       "      <td>1000</td>\n",
       "      <td>5.000000e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>1</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>2</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>3</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>4</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>5</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>6</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>10</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>100</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-07</td>\n",
       "      <td>1000</td>\n",
       "      <td>5.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>1</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>2</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>3</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>4</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>5</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>6</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>10</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>100</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-05</td>\n",
       "      <td>1000</td>\n",
       "      <td>5.000000e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>1</td>\n",
       "      <td>5.000000e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>2</td>\n",
       "      <td>5.000000e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>3</td>\n",
       "      <td>5.000000e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1860</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>10</td>\n",
       "      <td>3.511075e-162</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1861</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>100</td>\n",
       "      <td>3.511075e-162</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1862</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-04</td>\n",
       "      <td>1000</td>\n",
       "      <td>3.511075e-162</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1863</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>1</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1864</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>2</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1865</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>3</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1866</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>4</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1867</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>5</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1868</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>6</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1869</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>10</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1870</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>100</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1871</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-03</td>\n",
       "      <td>1000</td>\n",
       "      <td>3.511075e-132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1872</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>1</td>\n",
       "      <td>3.511075e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1873</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>2</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1874</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>3</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1875</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>4</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1876</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>5</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1877</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>6</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1878</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>10</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1879</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>100</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1880</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-02</td>\n",
       "      <td>1000</td>\n",
       "      <td>3.511077e-102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1881</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>1</td>\n",
       "      <td>3.511075e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1882</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>2</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1883</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>3</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1884</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>4</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1885</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>5</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1886</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>6</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1887</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>10</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1888</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>100</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1889</th>\n",
       "      <td>30</td>\n",
       "      <td>1.000000e-01</td>\n",
       "      <td>1000</td>\n",
       "      <td>3.511358e-72</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1890 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       v             u  n_terms            out\n",
       "0      1  1.000000e-16        1   5.000000e-17\n",
       "1      1  1.000000e-16        2   5.000000e-17\n",
       "2      1  1.000000e-16        3   5.000000e-17\n",
       "3      1  1.000000e-16        4   5.000000e-17\n",
       "4      1  1.000000e-16        5   5.000000e-17\n",
       "5      1  1.000000e-16        6   5.000000e-17\n",
       "6      1  1.000000e-16       10   5.000000e-17\n",
       "7      1  1.000000e-16      100   5.000000e-17\n",
       "8      1  1.000000e-16     1000   5.000000e-17\n",
       "9      1  1.000000e-07        1   5.000000e-08\n",
       "10     1  1.000000e-07        2   5.000000e-08\n",
       "11     1  1.000000e-07        3   5.000000e-08\n",
       "12     1  1.000000e-07        4   5.000000e-08\n",
       "13     1  1.000000e-07        5   5.000000e-08\n",
       "14     1  1.000000e-07        6   5.000000e-08\n",
       "15     1  1.000000e-07       10   5.000000e-08\n",
       "16     1  1.000000e-07      100   5.000000e-08\n",
       "17     1  1.000000e-07     1000   5.000000e-08\n",
       "18     1  1.000000e-05        1   5.000000e-06\n",
       "19     1  1.000000e-05        2   5.000000e-06\n",
       "20     1  1.000000e-05        3   5.000000e-06\n",
       "21     1  1.000000e-05        4   5.000000e-06\n",
       "22     1  1.000000e-05        5   5.000000e-06\n",
       "23     1  1.000000e-05        6   5.000000e-06\n",
       "24     1  1.000000e-05       10   5.000000e-06\n",
       "25     1  1.000000e-05      100   5.000000e-06\n",
       "26     1  1.000000e-05     1000   5.000000e-06\n",
       "27     1  1.000000e-04        1   5.000000e-05\n",
       "28     1  1.000000e-04        2   5.000000e-05\n",
       "29     1  1.000000e-04        3   5.000000e-05\n",
       "...   ..           ...      ...            ...\n",
       "1860  30  1.000000e-04       10  3.511075e-162\n",
       "1861  30  1.000000e-04      100  3.511075e-162\n",
       "1862  30  1.000000e-04     1000  3.511075e-162\n",
       "1863  30  1.000000e-03        1  3.511075e-132\n",
       "1864  30  1.000000e-03        2  3.511075e-132\n",
       "1865  30  1.000000e-03        3  3.511075e-132\n",
       "1866  30  1.000000e-03        4  3.511075e-132\n",
       "1867  30  1.000000e-03        5  3.511075e-132\n",
       "1868  30  1.000000e-03        6  3.511075e-132\n",
       "1869  30  1.000000e-03       10  3.511075e-132\n",
       "1870  30  1.000000e-03      100  3.511075e-132\n",
       "1871  30  1.000000e-03     1000  3.511075e-132\n",
       "1872  30  1.000000e-02        1  3.511075e-102\n",
       "1873  30  1.000000e-02        2  3.511077e-102\n",
       "1874  30  1.000000e-02        3  3.511077e-102\n",
       "1875  30  1.000000e-02        4  3.511077e-102\n",
       "1876  30  1.000000e-02        5  3.511077e-102\n",
       "1877  30  1.000000e-02        6  3.511077e-102\n",
       "1878  30  1.000000e-02       10  3.511077e-102\n",
       "1879  30  1.000000e-02      100  3.511077e-102\n",
       "1880  30  1.000000e-02     1000  3.511077e-102\n",
       "1881  30  1.000000e-01        1   3.511075e-72\n",
       "1882  30  1.000000e-01        2   3.511358e-72\n",
       "1883  30  1.000000e-01        3   3.511358e-72\n",
       "1884  30  1.000000e-01        4   3.511358e-72\n",
       "1885  30  1.000000e-01        5   3.511358e-72\n",
       "1886  30  1.000000e-01        6   3.511358e-72\n",
       "1887  30  1.000000e-01       10   3.511358e-72\n",
       "1888  30  1.000000e-01      100   3.511358e-72\n",
       "1889  30  1.000000e-01     1000   3.511358e-72\n",
       "\n",
       "[1890 rows x 4 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(res, columns=['v','u','n_terms','out'])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Laplace():\n",
    "    def __init__(self):\n",
    "        self.mus = None\n",
    "        self.bs = None\n",
    "        self.ws = None\n",
    "        self.sigmas = None\n",
    "        self.mus = None\n",
    "        self.cuda = None\n",
    "        self.initialized = False\n",
    "        \n",
    "    @property\n",
    "    def shape(self):\n",
    "        assert self.initialized == True\n",
    "        print('ws shape: ', self.ws.shape)\n",
    "        print('mus shape: ', self.mus.shape)\n",
    "        print('sigmas shape: ', self.sigmas.shape)\n",
    "        \n",
    "    def load(self, ws, sigmas, mus):\n",
    "        self.d = mus.shape[3]\n",
    "        self.n_gaussians = ws.shape[1]\n",
    "        self.ws = ws\n",
    "        self.sigmas = sigmas\n",
    "        self.mus = mus\n",
    "        self.initialized = True\n",
    "        print('GMM Loaded and initialized')\n",
    "\n",
    "    def initialize(self, d, n_gaussians, sigma, mu_min, mu_max, cuda, initialization = 'random'):\n",
    "        assert initialization in ['random','orthogonal']\n",
    "        \n",
    "        self.d = d\n",
    "        self.n_gaussians = n_gaussians\n",
    "        self.initialization = initialization\n",
    "        self.cuda = cuda\n",
    "        if initialization == 'random':\n",
    "            self.ws, self.mus, self.sigmas = make_random_GMM(d, n_gaussians, sigma, mu_min, mu_max, cuda)\n",
    "        elif initialization == 'orthogonal':\n",
    "            self.ws, self.mus, self.sigmas = make_orthogonal_GMM(d, n_gaussians, sigma, cuda)\n",
    "        print('GMM Initialized')\n",
    "        self.initialized = True\n",
    "    \n",
    "    def sample(self,n_points = 1):\n",
    "        assert self.initialized == True\n",
    "        print(self.mus.shape)\n",
    "        mus = self.mus.repeat(1,1,n_points,1)\n",
    "        sigmas = self.sigmas.repeat(1,1,n_points,1)\n",
    "        chosen_gaussians = torch.multinomial(self.ws, n_points, replacement=True)\n",
    "\n",
    "        batch_size = self.ws.shape[0]\n",
    "        chosen_gaussians = chosen_gaussians.view(batch_size, 1, n_points,1).repeat(1,1,1,self.d)\n",
    "\n",
    "        ep = torch.randn(batch_size, 1, n_points,self.d)\n",
    "        if self.cuda:\n",
    "            ep = ep.cuda()\n",
    "        print(mus.shape, chosen_gaussians.shape)\n",
    "        u_mu = torch.gather(mus, 1, chosen_gaussians)\n",
    "        u_sigma = torch.gather(sigmas, 1, chosen_gaussians)\n",
    "        u_i = u_mu + ep*u_sigma\n",
    "        return u_i\n",
    "        \n",
    "    def tsne_plot(self, n_points=200):\n",
    "        assert self.initialized == True\n",
    "        plt.figure(figsize=(10,10))\n",
    "        plt.title('Latent space - TSNE')\n",
    "        plt.xlabel('Dimension 1')\n",
    "        plt.ylabel('Dimension 2')\n",
    "        samples = self.sample(n_points)\n",
    "        samples = sample_from_gmm(n_points, self.ws, self.mus, self.sigmas, self.cuda)\n",
    "        samples = samples.squeeze(0).squeeze(0)\n",
    "        \n",
    "        tsne = TSNE(n_components=2, perplexity=20)\n",
    "        gmm_tsne = tsne.fit_transform(samples.cpu().detach().numpy())\n",
    "        plt.scatter(*gmm_tsne.T, c='blue', marker='o', alpha=0.5)\n",
    "        \n",
    "    def pdf(self, u):\n",
    "        #get the probability density function(s) at point(s) u\n",
    "        #for a mixed gaussian model with weights given by self.ws,\n",
    "        #means given by self.mus and diagonal covariances given by self.sigmas**2\n",
    "\n",
    "        #this is the exponent of the multivariate gaussian\n",
    "        #The equation for a multivariate gaussian is :\n",
    "        #  (2*pi)^(-d/2)*det(cov)^(-1/2)*e^((1/2)*(x-mu).T*inverse_cov*(x-mu))\n",
    "        #in our case the calculation is simplified because\n",
    "        #we assume the covariance matrix is diagonal, so the determinant is just the product\n",
    "        #of the elements, and the inverse is just 1/sigma\n",
    "        ws = self.ws.unsqueeze(2).unsqueeze(3)\n",
    "        covs = self.sigmas**2\n",
    "        \n",
    "        #to avoid infinities - use log covariance\n",
    "        log_covs = torch.log(covs)\n",
    "        \n",
    "        log_det_cov = (torch.sum(log_covs, 2)).unsqueeze(2)\n",
    "        exponent = (-1/2)*((u-self.mus)*(1/covs)*(u-self.mus)).sum(dim=2, keepdim=True)\n",
    "                \n",
    "        if self.cuda:\n",
    "            ws = ws.type(torch.cuda.DoubleTensor)\n",
    "            log_det_cov = log_det_cov.type(torch.cuda.DoubleTensor)\n",
    "            exponent = exponent.type(torch.cuda.DoubleTensor)\n",
    "        else:\n",
    "            ws = ws.type(torch.DoubleTensor)\n",
    "            log_det_cov = log_det_cov.type(torch.DoubleTensor)\n",
    "            exponent = exponent.type(torch.DoubleTensor)\n",
    "        \n",
    "        #multiplier = ((2*np.pi)**(-self.d/2))*np.e**((-1/2)*log_det_cov)\n",
    "        multiplier = np.e**((-1/2)*log_det_cov)#((det_cov)**(-1/2))\n",
    "        \n",
    "        if self.cuda:\n",
    "            multiplier = multiplier.type(torch.cuda.DoubleTensor)\n",
    "        else:\n",
    "            multiplier = multiplier.type(torch.DoubleTensor)\n",
    "        pdf = multiplier*np.e**exponent\n",
    "        #since ours is a mixture of gaussians, we sum the pdf over the number of gaussians we have \n",
    "        #multiplied by the weight each gaussian has\n",
    "        pdf = ((ws*pdf).sum(dim=1,keepdim=True))#.type(torch.cuda.DoubleTensor)   \n",
    "        return pdf\n",
    "\n",
    "\n",
    "    \n",
    "d = 30\n",
    "n_gaussians = 6\n",
    "sigma = 1\n",
    "mu_min = -50\n",
    "mu_max = 50\n",
    "cuda = True\n",
    "n_points=1000\n",
    "        \n",
    "q = GMM()\n",
    "q.initialize(d, n_gaussians, sigma, mu_min, mu_max, cuda, initialization='orthogonal')\n",
    "q.tsne_plot(n_points=500)\n",
    "k = q.sample(n_points=500)\n",
    "q.pdf(k)\n",
    "#p = GMM()\n",
    "#p.initialize(d, n_gaussians, sigma, mu_min, mu_max, cuda, initialization='random')\n",
    "#p.tsne_plot(n_points=500)"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}