- "print(__doc__)\n\nfrom distutils.version import LooseVersion\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nfrom sklearn.decomposition import SparseCoder\n\n\ndef ricker_function(resolution, center, width):\n \"\"\"Discrete sub-sampled Ricker (Mexican hat) wavelet\"\"\"\n x = np.linspace(0, resolution - 1, resolution)\n x = ((2 / ((np.sqrt(3 * width) * np.pi ** 1 / 4)))\n * (1 - ((x - center) ** 2 / width ** 2))\n * np.exp((-(x - center) ** 2) / (2 * width ** 2)))\n return x\n\n\ndef ricker_matrix(width, resolution, n_components):\n \"\"\"Dictionary of Ricker (Mexican hat) wavelets\"\"\"\n centers = np.linspace(0, resolution - 1, n_components)\n D = np.empty((n_components, resolution))\n for i, center in enumerate(centers):\n D[i] = ricker_function(resolution, center, width)\n D /= np.sqrt(np.sum(D ** 2, axis=1))[:, np.newaxis]\n return D\n\n\nresolution = 1024\nsubsampling = 3 # subsampling factor\nwidth = 100\nn_components = resolution // subsampling\n\n# Compute a wavelet dictionary\nD_fixed = ricker_matrix(width=width, resolution=resolution,\n n_components=n_components)\nD_multi = np.r_[tuple(ricker_matrix(width=w, resolution=resolution,\n n_components=n_components // 5)\n for w in (10, 50, 100, 500, 1000))]\n\n# Generate a signal\ny = np.linspace(0, resolution - 1, resolution)\nfirst_quarter = y < resolution / 4\ny[first_quarter] = 3.\ny[np.logical_not(first_quarter)] = -1.\n\n# List the different sparse coding methods in the following format:\n# (title, transform_algorithm, transform_alpha,\n# transform_n_nozero_coefs, color)\nestimators = [('OMP', 'omp', None, 15, 'navy'),\n ('Lasso', 'lasso_lars', 2, None, 'turquoise'), ]\nlw = 2\n# Avoid FutureWarning about default value change when numpy >= 1.14\nlstsq_rcond = None if LooseVersion(np.__version__) >= '1.14' else -1\n\nplt.figure(figsize=(13, 6))\nfor subplot, (D, title) in enumerate(zip((D_fixed, D_multi),\n ('fixed width', 'multiple widths'))):\n plt.subplot(1, 2, subplot + 1)\n plt.title('Sparse coding against %s dictionary' % title)\n plt.plot(y, lw=lw, linestyle='--', label='Original signal')\n # Do a wavelet approximation\n for title, algo, alpha, n_nonzero, color in estimators:\n coder = SparseCoder(dictionary=D, transform_n_nonzero_coefs=n_nonzero,\n transform_alpha=alpha, transform_algorithm=algo)\n x = coder.transform(y.reshape(1, -1))\n density = len(np.flatnonzero(x))\n x = np.ravel(np.dot(x, D))\n squared_error = np.sum((y - x) ** 2)\n plt.plot(x, color=color, lw=lw,\n label='%s: %s nonzero coefs,\\n%.2f error'\n % (title, density, squared_error))\n\n # Soft thresholding debiasing\n coder = SparseCoder(dictionary=D, transform_algorithm='threshold',\n transform_alpha=20)\n x = coder.transform(y.reshape(1, -1))\n _, idx = np.where(x != 0)\n x[0, idx], _, _, _ = np.linalg.lstsq(D[idx, :].T, y, rcond=lstsq_rcond)\n x = np.ravel(np.dot(x, D))\n squared_error = np.sum((y - x) ** 2)\n plt.plot(x, color='darkorange', lw=lw,\n label='Thresholding w/ debiasing:\\n%d nonzero coefs, %.2f error'\n % (len(idx), squared_error))\n plt.axis('tight')\n plt.legend(shadow=False, loc='best')\nplt.subplots_adjust(.04, .07, .97, .90, .09, .2)\nplt.show()"
0 commit comments