{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genre recognition: feature extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The audio genre recognition pipeline:\n",
    "1. GTZAN\n",
    "2. pre-processing\n",
    "3. unsupervised feature extraction\n",
    "4. classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open questions:\n",
    "* Rescale the dataset ? We need to for the algorithm to converge.\n",
    "    * Rescale $n$ features in [0,1] --> converge. But we need to learn the transform.\n",
    "    * Normalize each sample to unit norm --> converge. But higher objective and less sparse Z. We also loose the generative ability of our model.\n",
    "* Is there a way to programmatically assess convergence ? Easy for us to look at the objective function, but for a machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyper-parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* `m`:  number of atoms in the dictionary, sparse code length.\n",
    "* `ls`: weight of the sparse codes l1 penalty (redundant).\n",
    "* `ld`: weigth of the dictionary l2 penalty.\n",
    "* `le`: weight of the encoder l2 penalty.\n",
    "* `lg`: weight of the Dirichlet energy (via the graph Laplacian).\n",
    "* `rtol`: stopping criterion for inner and outer loops.\n",
    "* `N_inner`: hard limit on inner iterations.\n",
    "* `N_outer`: hard limit on outer iterations.\n",
    "* `Ngenres, Nclips, Nframes`: a way to reduce the size of the dataset.\n",
    "* `noise_std`: standard deviation of the Gaussian noise to be added to the data.\n",
    "* `folder`: relative path to HDF5 files.\n",
    "* `filename_*`: name of the HDF5 file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if 'p' in globals().keys():\n",
    "    # Hyper-parameters passed by the experiment runner.\n",
    "    for key, value in p.items():\n",
    "        globals()[key] = value\n",
    "else:\n",
    "    m = 64  # 64, 128, 512\n",
    "    ls = 1\n",
    "    ld = 10\n",
    "    le = None\n",
    "    lg = 1\n",
    "    rtol = 1e-5  # 1e-3, 1e-5, 1e-7\n",
    "    N_inner = 500\n",
    "    N_outer = 50\n",
    "    Ngenres, Nclips, Nframes = 10, 100, 644\n",
    "    noise_std = 0\n",
    "    folder = 'data'\n",
    "    filename_audio = 'audio.hdf5'\n",
    "    filename_graph = 'graph.hdf5'\n",
    "    filename_features = 'features.hdf5'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os, time\n",
    "import numpy as np\n",
    "import scipy.sparse\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# Import auto-encoder definition.\n",
    "%run -n auto_encoder.ipynb\n",
    "#import auto_encoder\n",
    "\n",
    "# Profiling.\n",
    "%reload_ext memory_profiler\n",
    "%reload_ext line_profiler\n",
    "import objgraph\n",
    "\n",
    "#%load_ext autoreload\n",
    "#%autoreload 2\n",
    "\n",
    "toverall = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def datinfo(X, name='Dataset'):\n",
    "    r\"\"\"Print dataset size and dimensionality\"\"\"\n",
    "    print('{}:\\n'\n",
    "          '  size: N={:,} x n={} -> {:,} floats\\n'\n",
    "          '  dim: {:,} features per clip\\n'\n",
    "          '  shape: {}'\n",
    "          .format(name, np.prod(X.shape[:-1]), X.shape[-1],\n",
    "                  np.prod(X.shape), np.prod(X.shape[2:]), X.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = os.path.join(folder, filename_audio)\n",
    "with h5py.File(filename, 'r') as audio:\n",
    "\n",
    "    # Display HDF5 attributes.\n",
    "    print('Attributes:')\n",
    "    for attr in audio.attrs:\n",
    "        print('  {} = {}'.format(attr, audio.attrs[attr]))\n",
    "    sr = audio.attrs['sr']\n",
    "    labels = audio.attrs['labels']\n",
    "\n",
    "    # Show datasets, their dimensionality and data type.\n",
    "    print('Datasets:')\n",
    "    for dname, dset in audio.items():\n",
    "        print('  {:2}: {:24}, {}'.format(dname, dset.shape, dset.dtype))\n",
    "\n",
    "    # Choose dataset: Xa, Xs.\n",
    "    X = audio.get('Xs')\n",
    "\n",
    "    # Full dataset.\n",
    "    n = X.shape[-1]\n",
    "    datinfo(X, 'Full dataset')\n",
    "    print(type(X))\n",
    "\n",
    "    # Load data into memory as a standard NumPy array.\n",
    "    X = X[:Ngenres,:Nclips,:Nframes,...]\n",
    "    datinfo(X, 'Reduced dataset')\n",
    "    print(type(X))\n",
    "\n",
    "    # Resize in place without memory loading via hyperslab.\n",
    "    # Require chunked datasets.\n",
    "    #X.resize((Ngenres, Nclips, Nframes, 2, n))\n",
    "\n",
    "# Squeeze dataset to a 2D array. The auto-encoder does not\n",
    "# care about the underlying structure of the dataset.\n",
    "X.resize(Ngenres * Nclips * Nframes * 2, n)\n",
    "print('Data: {}, {}'.format(X.shape, X.dtype))\n",
    "\n",
    "# Independently rescale each feature.\n",
    "# To be put in an sklearn Pipeline to avoid transductive learning.\n",
    "X -= np.min(X, axis=0)\n",
    "X /= np.max(X, axis=0)\n",
    "\n",
    "# Independently normalize each sample.\n",
    "#X /= np.linalg.norm(X, axis=1)[:,np.newaxis]\n",
    "\n",
    "# Add Gaussian noise.\n",
    "if noise_std is not 0:\n",
    "    X += np.random.normal(scale=noise_std, size=X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = os.path.join(folder, filename_graph)\n",
    "with h5py.File(filename, 'r') as graph:\n",
    "\n",
    "    # Display HDF5 attributes.\n",
    "    print('Attributes:')\n",
    "    for attr in graph.attrs:\n",
    "        print('  {} = {}'.format(attr, graph.attrs[attr]))\n",
    "\n",
    "    # Show datasets, their dimensionality and data type.\n",
    "    print('Datasets:')\n",
    "    for dname, dset in graph.items():\n",
    "        print('  {:10}: {:10}, {}'.format(dname, dset.shape, dset.dtype))\n",
    "\n",
    "    # Data: Laplacian matrix.\n",
    "    pars = []\n",
    "    for par in ('data', 'indices', 'indptr', 'shape'):\n",
    "        pars.append(graph.get('L_'+par))\n",
    "    L = scipy.sparse.csr_matrix(tuple(pars[:3]), shape=pars[3])\n",
    "\n",
    "if L.shape != (X.shape[0], X.shape[0]):\n",
    "    raise ValueError('Graph size does not correspond to data size.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Size of training data and parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = Ngenres * Nclips * Nframes * 2\n",
    "sizeX = N * n / 2.**20\n",
    "sizeZ = N * m / 2.**20\n",
    "sizeD = n * m / 2.**10\n",
    "sizeE = m * n / 2.**10\n",
    "# 32 bits float\n",
    "print('Size X: {:.1f} M --> {:.1f} MiB'.format(sizeX, sizeX*4))\n",
    "print('Size Z: {:.1f} M --> {:.1f} MiB'.format(sizeZ, sizeZ*4))\n",
    "print('Size D: {:.1f} k --> {:.1f} kiB'.format(sizeD, sizeD*4))\n",
    "print('Size E: {:.1f} k --> {:.1f} kiB'.format(sizeE, sizeE*4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ae = auto_encoder(m=m, ls=ls, ld=ld, le=le, lg=lg, rtol=rtol, xtol=None, N_inner=N_inner, N_outer=N_outer)\n",
    "tstart = time.time()\n",
    "Z = ae.fit_transform(X, L)\n",
    "time_features = time.time() - tstart\n",
    "print('Elapsed time: {:.0f} seconds'.format(time_features))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations:\n",
    "* Memory efficiency:\n",
    "    * m=64, 20 songs: 600 MiB --> 170 MiB (pyul mem optimization) --> 120 MiB (float32) --> 150MiB (graph)\n",
    "    * m=64, 40 songs: 900 MiB --> 170 MiB (pyul mem optimization) --> 150 (float32) MiB\n",
    "    * m=128, 200 songs: 800 MiB (pyul mem optimization)\n",
    "    * m=128, 400 songs: 2 GiB (pyul mem optimization) --> 1 GiB (float32)\n",
    "* Time efficiency:\n",
    "    * m=64, 20 songs: 370s --> 230s (float32) --> 515s (graph)\n",
    "    * m=128, 200 songs: 9048s (pyul mem optim) --> 1779s (CDK, ld=10, 10 outer) --> 1992s (CDK, ld=10, 20 outer) --> 3877s (CDK, ld=100, 15 outer)\n",
    "    * m=512, 200 songs: 8814s (CDK, ld=10, 15 outer)\n",
    "    * m=128, 400 songs: 19636s=5h30 (pyul mem optim)\n",
    "    * m=512, 500 songs: 19995s=5h30 (CDK, ld=10, 20 outer)\n",
    "    * m=512, 1000 songs: 35429s=9h50 (CDK, ld=10, 15 outer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time analysis:\n",
    "1. Use ATLAS or OpenBLAS instead of numpy BLAS implementation.\n",
    "1. Multi-threaded ATLAS or OpenBLAS (may not be worth it if we are memory bandwidth limited).\n",
    "1. Compute with float32, it saves memory bandwidth. CPU is then more efficiently used for matrix multiplication.\n",
    "1. Projection in the L2-ball, not on the sphere. It is a convex constraint.\n",
    "1. PyUNLocBoX: do not evaluate the objective at each iteration (configuration). Profile.\n",
    "1. PyUNLocBoX: does capability check use time ? Only once per outer loop iteration. Profile.\n",
    "1. Multiple threads working on independent sub-problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    %prun Z = ae.fit_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Space analysis:\n",
    "1. Avoid copies in PyUNLocBoX.\n",
    "2. Modify data in place and pass references.\n",
    "3. Store data in float64 ? Or compute in float32 ? 32 bits precision should be enough.\n",
    "4. Store Z as scipy.sparse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    import gc\n",
    "    gc.collect()\n",
    "    objgraph.show_most_common_types()\n",
    "    from pyunlocbox import solvers, functions\n",
    "    %mprun -f ae.fit_transform -f ae._minD -f ae._minZ -f solvers.solve -f solvers.forward_backward._pre -f solvers.forward_backward._fista -f functions.norm_l1._prox -T profile.txt ae.fit_transform(X)\n",
    "    #%mprun -f solvers.solve -f solvers.forward_backward._pre -f solvers.forward_backward._fista -f functions.norm_l1._prox -T profile.txt ae.fit_transform(X)\n",
    "    gc.collect()\n",
    "    objgraph.show_most_common_types()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    from pympler import tracker\n",
    "    tr = tracker.SummaryTracker()\n",
    "    Z = ae.fit_transform(X)\n",
    "    tr.print_diff()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective and convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ret = ae.plot_objective()\n",
    "iterations_inner, iterations_outer = ret[:2]\n",
    "objective_g, objective_h, objective_i, objective_j = ret[2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sparsity = sparse_codes(Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations:\n",
    "* The learned atoms seem to represent harmonies and harmonics.\n",
    "* The atoms themselves look sparse. Should we add some prior knowledge on the dictionary ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if ld is not None:\n",
    "    dictenc(ae.D)\n",
    "    atoms_D = atoms(ae.D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if le is not None:\n",
    "    dictenc(ae.E, enc=True)\n",
    "    atoms(ae.E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will store more Z when the various approximations will be implemented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = os.path.join(folder, filename_features)\n",
    "\n",
    "# Remove existing HDF5 file without warning if non-existent.\n",
    "try:\n",
    "    os.remove(filename)\n",
    "except OSError:\n",
    "    pass\n",
    "\n",
    "# Create HDF5 file and datasets.\n",
    "with h5py.File(filename, 'w') as features:\n",
    "\n",
    "    # Metadata.\n",
    "    features.attrs['sr'] = sr\n",
    "    features.attrs['labels'] = labels\n",
    "\n",
    "    # Data.\n",
    "    features.create_dataset('X', data=X.reshape(Ngenres, Nclips, Nframes, 2, n))\n",
    "    features.create_dataset('Z', data=Z.reshape(Ngenres, Nclips, Nframes, 2, Z.shape[-1]))\n",
    "    if ld is not None:\n",
    "        features.create_dataset('D', data=ae.D)\n",
    "    if le is not None:\n",
    "        features.create_dataset('E', data=ae.E)\n",
    "\n",
    "    # Show datasets, their dimensionality and data type.\n",
    "    print('Datasets:')\n",
    "    for dname, dset in features.items():\n",
    "        print('  {:2}: {:22}, {}'.format(dname, dset.shape, dset.dtype))\n",
    "\n",
    "    # Display HDF5 attributes.\n",
    "    print('Attributes:')\n",
    "    for name, value in features.attrs.items():\n",
    "        print('  {} = {}'.format(name, value))\n",
    "\n",
    "print('Overall time: {:.0f} seconds'.format(time.time() - toverall))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}