acdc_iterations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import ase.io\n",
    "from time import time\n",
    "from copy import deepcopy\n",
    "\n",
    "from aml_storage import Labels, Block, Descriptor\n",
    "\n",
    "from utils.librascal import RascalSphericalExpansion\n",
    "from utils.clebsh_gordan import ClebschGordanReal\n",
    "\n",
    "from utils.acdc import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrates the iterative construction of atom-centered density correlations.\n",
    "The calculations incorporate ideas from the N-body iterative contraction of equivariants [paper](https://doi.org/10.1063/5.0021116) and the optimal radial basis [paper](https://doi.org/10.1063/5.0057229)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covariance-optimal basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by computing the spherical expansion representation for a set of methane structures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "rascal_hypers = {\n",
    "    \"interaction_cutoff\": 3.5,\n",
    "    \"cutoff_smooth_width\": 0.5,\n",
    "    \"max_radial\": 3,\n",
    "    \"max_angular\": 2,\n",
    "    \"gaussian_sigma_type\": \"Constant\",\n",
    "    \"compute_gradients\":  True,\n",
    "}\n",
    "\n",
    "nmax, lmax = rascal_hypers[\"max_radial\"], rascal_hypers[\"max_angular\"]\n",
    "calculator = RascalSphericalExpansion(rascal_hypers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = ase.io.read(\"./data/random-methane-10k.extxyz\", \"::100\")\n",
    "for f in frames:\n",
    "    f.cell = [100,100,100]\n",
    "    f.positions += 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/michele/lavoro/code/aml-storage-examples/utils/librascal.py:92: UserWarning: numpy/forward gradients are currently broken with librascal,please use rascaline instead\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "density_expansion = calculator.compute(frames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, select the C,H block, which is the most significant for this dataset, and project to the optimal covariance basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for each (C,H) block, construct the optimal covariance basis\n",
    "opt_eve = {}\n",
    "opt_eva = {}\n",
    "full_norms = np.zeros(len(frames))\n",
    "for index, block in density_expansion:\n",
    "    if index[\"center_species\"] == 6 and index[\"neighbor_species\"] == 1:\n",
    "        feats = block.values.reshape(-1, len(block.features))\n",
    "        full_norms += (block.values**2).sum(axis=-1).sum(axis=-1)\n",
    "        cov = feats.T @ feats / feats.shape[0]\n",
    "        eva, eve = sp.linalg.eigh(cov)\n",
    "        eva = eva[::-1]; eve = eve[:,::-1]\n",
    "        opt_eva[index[\"spherical_harmonics_l\"]] = eva\n",
    "        opt_eve[index[\"spherical_harmonics_l\"]] = eve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmOUlEQVR4nO3deWxc22Hf8e/hvm/DTSIlUSQlPe0iObKfpKQxmhp1nThGUjTIhgCJ4RcUTZsWBYoURZEWbWGgKIIsTtO+JmmSJnAW1w3sIKiLNk0dSPJ7Gorad1IbKXEb7vtw5vSPM7yc4dNCkUPOcn8f4AEW75BzZkzyx3vOub9rrLWIiIj/5KV7ACIikh4KABERn1IAiIj4lAJARMSnFAAiIj6lABAR8amCdA/gTerr621bW1u6hyEiklV6e3vHrbUNb3tcRgdAW1sboVAo3cMQEckqxpinm3mcpoBERHxKASAi4lMKABERn8rIADDGfMEY8+H09HS6hyIikrMyMgCstd+y1n5QXV2d7qGIiOSsjAwAERHZeQoAERGfUgCIiPiUAkBExKcUACIiPqUAEBHxKQWAiIhPKQBERHxKASAi4lMKABERn8rIAFAXkIjIzsvIAFAXkIjIzsvIABARkZ2nABAR8SkFgIiITykARER8SgEgIuJTCgAREZ9SAIiI+JQCQETEpxQAIiI+pQAQEfEpBYCIiE8pAEREfEoBICLiUxkZAKqDFhHZeRkZAKqDFhHZeRkZACIisvMUACIiPqUAEBHxKQWAiIhPKQBERHxKASAi4lMKABERn1IAiIj4lAJARMSnFAAiIj6lABAR8SkFgIiITykARER8SgEgIuJTCgAREZ9SAIiI+JQCQETEpxQAIiI+pQAQEfEpBYCIiE9lZAAYY75gjPlweno63UMREclZGRkA1tpvWWs/qK6uTvdQRERyVkYGgIiI7DwFgIiITykARER8SgEgIuJTCgAREZ9SAIiI+JQCQETEpxQAIiI+pQAQEfEpBYCIiE8pAEREfEoBICLiUwoAERGfUgCIiPiUAkBExKcUACIiPqUAEBHxKQWAiIhPKQBERHxKASAi4lMKABERn1IAiIj4lAJARMSnFAAiIj6lABAR8SkFgIiITykARER8SgEgIuJTCgAREZ9SAIiI+JQCQETEpxQAIiI+tWsBYIz5jDHmr40x/8kY85ndel4REXm1TQWAMeZ3jDGjxphbGz7+OWPMfWPMI2PML77ly1hgDigBBrc2XBERSZWCTT7ud4GvAr+/9gFjTD7wG8Bncb/QrxhjvgnkA1/Z8Pk/C/y1tfb/GWOagF8GfnJ7QxcRke3YVABYa79jjGnb8OFPAY+stQMAxpg/Ar5orf0K8INv+HKTQPEWxioiIim02TOAV2kBnif8exD49OsebIz5EeBvAzW4s4nXPe4D4AOA/fv3b2N4IiLyJtsJgHdirf0G8I1NPO5D4EOAYDBod3pcIiJ+tZ1dQEPAvoR/t8Y/JiIiWWA7AXAFOGSMOWiMKQJ+DPhmaoYlIiI7bbPbQL8GXAaOGGMGjTFfstauAj8PfBu4C/yJtfb2zg1VRERSabO7gH78NR//C+AvUjoiwBjzBeALnZ2dqf7SIiISl5FVENbab1lrP6iurk73UEREclZGBoCIiOy8XdsGuhVLkSh3Xsxs6XPryotori5J8YhERHJHRgfAw9E5Pv9rf73lzz9YX865jgDnOwKcaw8QqNAFyCIia4y1mXut1eHjp+2v/fG3t/S5g5MLXO4P89HjCeaWVwF4r7mS8x31nO8I8On2OipLClM5XBGRjGCM6bXWBt/6uEwOgGAwaEOh0La+xmo0xs2haS71h7nUP07oySTLqzHy8wwnW6o53xHgfEc9PQdqKS3KT9HIRUTSJ6sDIGEb6JcfPnyY0q+9FInS92yKS/3jXOoPc/35FKsxS1F+Ht0HarwzhNP7aijM1xq5iGSfrA6ANak4A3ibueVVrjyZ4HJ/mIuPxrnzcgZroawon7NtdVzodGcIR/dUkZ9ndnQsIiKpoADYosn5FT56HI5PGYV5NDoHQHVpIe+313Gh050hdDRUYIwCQUQyz2YDIKN3AaVDbXkRnzuxh8+d2APAyMwSl+PrBxcfhfn27REAGiqL4+sH7gxhX11ZOoctIvLOdAbwDqy1PJ9Y9NYPLvWHGZ9bBmBfXSnn2+s53xngXEeAxkpdgyAi6aEpoF1greXh6ByXHrlA+O5AmJklt+X0UGOFu/6go5732+uoKStK82hFxC8UAGkQjVnuvJjhYvwM4crjCRYjUYyBE3ur44EQ4GxbHeXFmn0TkZ2R1QGwk9tAd9PKaozrg1NcjJ8h9D2bJBK1FOQZzuyr4Xx8Qblrfw3FBboGQURSI6sDYE22nQG8zeJKlNDTCbd+8Gicm0PTxCyUFOYRPFDH+fiW0xN7qyjQNQgiskUKgCwwvRjh48cTblH5UZj7I7MAVBYX8On2Os511HOhM8DhxkrydA2CiGyStoFmgerSQj57rInPHmsCYGx2me8OuN1Fl/vH+d93RwEIlBfxfnzL6YWOeg4EynQNgohsm84AMthaod3l/jAX+8cZmXFbTvdWl3AuXllxvjPAnurSNI9URDKJpoByjLWWx+PzXIyfHVzuDzO5EAGg3au9dltOVXst4m8KgBwXi1nuDc96F6V9NBBmfiUKwNE9Vd5Vyp86qNprEb/JjQA41GxDv/KTW/vkkhpoDUJLECoaUjquTBSJ116vldqFnk6yEq+9PtWaXHtdUqgtpyK5LKsDYO06gO69hV/u/YXWrX2RxSmw7i9iatug9Wz8vyA0nYSC3L4ydykS5eqzyXiPUZhrz6eIJtReX+hwtRWnWlV7LZJrsjoA1mxrCmhlAV5eg8Er7r/nV2Bu2B3LL4a9Z9YDofUsVLVADu+smVte5craltP+cFLt9acO1nlnCMf2VGnLqUiWUwBsZC3MDMUDIeT+e9EHUbezhso962HQEnQBUVSemufOQJPzK96W00v94/SPzQNQU1bI+wcD3kVpHQ3l2nIqkmUUAJuxugIjt+KBED9TmHzsjpl8aDqeMHV0FgIdOXuWMDKz5F2Qdqk/zNDUIgCNXu11Pec6Aqq9FskCCoCtmh9PDoShq7DirtCltNadHaxNHbX0QGnN7o5vF6zVXq+V2l3uH2d8bgWA/XVlXqmdaq9FMpMCIFViURi7vx4IgyEYuwfE37f6I8lrCY1HIS+3dtkk1l5fjNdez8Zrrw83VXhnB+8fDFBdpi2nIummANhJS9PuzGAotH62sBB2xwrLoaU7eddRRWN6x5ti0Zjl9otp76Y4n6i9jq8fnG2rpaxIbSMiu00BsJusdWsHiVNHwzch5v5KpubA+hlC61loPgkFuXO17spqjGvPp7w1hL7nrva6MN/VXp/rqOdCR4Azqr0W2RUKgHSLLMLL68lTRzND7lh+Eew5nTx1VL0vZxaYF1ZWCT2Z9NYPEmuvz7bVcS5eaneipZp8bTkVSTkFQCaaHopPG11Z34a6uuSOVTQlB8LerpzZhjq9GOEjr+U0ofa6pIBPHwx4pXZHmiq15VQkBbI6AHLljmBvFY18chvqxIA7ZvKh6VjyNtS6DsjL/qt2x2aXuTzgzg4u9Yd5Gl4AXO31Wqnd+Y6Aaq9FtiirA2BNzp0BbMZ8GIZ6E7ah9sLyjDuW2G/UehZae9zW1Cw3OLngnR1cSqi9bqkpjQeCC4Xmam05FdmMnAiAw6cO269+86tb+tya4hpONpyktCDLu/JjMRh/kLyWMHoHbxtq4NCGbajHID97d95YaxkYn/dum3l5IMxUQu312g6j99sD1JXndp+TyFblRACUHiy1nf+qc8ufX5BXwInACYLNQXqaeuhq7KK8MAfm1Zdm3PqBV2txBRbG3bHCMtjbnbzrqLIpvePdhljMcnd4xms5/fjxxCdqry90BjjbptprkTU5EQAnu07ab/zlN7b0uS/nXxIaCdE73Mvt8G2iNkq+yedo3VGCzUGCTUG6mrqoKqpK8ajTwFqYfPKKbajuL2eq9ycHwp5TWbsNNRKNcWNw2ls/eFXt9YWOerpVey0+lhMBkKo1gIXIAtfGrhEaDtE70svN8ZtEYhEMhiN1Rwg2uUDoaeqhpqRm+wPPBJElGL6RPHU0/dwdyy+C5lPJU0c1+7NyG+pSJMrVp5Neqd31wWlXe12QR8/+2vgOo3pOtVar9lp8QwHwBkurS9wYu0HvSC+hkRDXx66zHG8F7azpdGHQ3EOwKUh9aX3Knz9tZl4mb0MdugqrrvSN8oYN21C7obgivePdgrnlVT5+HPZK7e68dAvo5V7ttautUO215DIFwDtYia5wa/yWmzIa6aVvtI/F+C/Gtqo2b8oo2BSkqTx759M/IRpxC8qJawnhR+6YyXMLyolTR4FDWbcNdcKrvXZTRgMJtdfn2gPxYjvVXktuUQBsQyQW4W74LqGREKHhEH2jfcxF5gBorWhdD4TmIC0VLbs+vh21MJG8DXWwF5an3bHiarf1dC0QWnqgrC69431Hw9NLXB4Y5+Ijt8voxbS7EK+pqtg7OzjfEaC1VrXXkr0UACkUjUW5P3mf0HDIO0uYWXFTC3vK93jrB8HmIPsr9+fWX5KxmDsrSNqGehtszB0PdG7Yhno8a7ahWmt5NuGuQbj4aJzL/WHC8+u11xc63dnBufYADZXZuWgu/qQA2EExG+Ph5EMvDHpHeplYmgCgobTBOzsINgU5WH0wtwIBYHluwzbUj2F+zB0rKHU1FolTR1V70jveTbLW8mBkzpsuelXt9fmOAJ9uD1Bdqi2nkrkUALvIWsvj6cfelFFoJMTYovuFWFdSR09TjztDaApyqPYQeSa75tHfylqYepa8lvDy+vo21KrWDdtQT0Nh5l/VuxqNcfvFjLfD6MqTCZYiMfIMnGip9krtgqq9lgyT1QGQ7V1A1lqezz5PCoSX8y8BqC6uprux25syeq/2PfJz7AYyQHwb6s0N21CfuWN5ha4SO/GeCbVtGb8NdXk1yrVnU15tRWLtdde+Wm/9oGt/LUUFORbyklWyOgDWZMsZwGYMzQ25bafxQHg+6/bkVxRW0NXY5V2tfCxwjMK8HJ1emB3ZsA21FyKuCI6y+uS1hJZuKK5M73jfYmFllStPJrnU79YPbg5NYxNqr9emjFR7LbtNAZDhRuZH3BlCfB3h8bS7GX1pQSlnGs54awgn6k9QlJ+jnTfRVRi7mzx1NP4gftDEt6Em7DqqP5LR21CnFyJ89/F6qd2DEbdzrLKkgPfb10vtDjdV5N66kGQUBUCWGV8cTzpDeDTl9uMX5xdzquGUdx3CqYZTlBRk/vz5li1OxrehJtRaLK1tQ61Kvt1mSxDKA+kd7xuMzi7x3YEJLj1yi8rPJtzZTn1FEefiZwfnOwLsr1PttaSWAiDLTS5NcnX0qldfcW/iHhZLYV4hJ+tPeovKZxrPUFaYw3vWYzGY6E9YS7gCIwnbUOvak9cSmk5AfmZOoT2fWODyQNgLhNHZ5NrrC50BzrWr9lq2TwGQY2ZWZugb6fOmjO6E7xC1UQpMAccCx7zqiq7GLiqLMnvufNtW5uHFteRQmBtxxwpKXrENdW9ah/sq1lr6x+a9Uruk2uuGcq/U7v32ALWqvZZ3pADIcfORea6NXvN2Gt0K32I1tkqeyeNI7RFvDaGnqYfq4up0D3dnWQvTgxu2oV6DqLuoi6qW5Bvp7D0DhZl1n4hYzHLn5Yy3frBWe20MHG2u8m6b+amDASqKteVU3kwB4DOLq4vcGLvhBcKNsRusxFYwGA7VHvLCoKeph0Bp5s6bp8zqMgzfSj5LmHrqjuUVuKmixKmjuvaM2obqaq+nvFK73mfrtdenW6u9HUaqvZZXUQD43HJ0mZtjN5MaT9cK7tqr25OuVm4oa0jzaHfJ3Gjy4vLQVYi4cjhK65IDoaUbSjLnzGkpEqX36aR3lfKNhNrr4IFar9ROtdcCCgDZIBKNcDt829t6em30GvPxX377K/cnNZ7uqciO6oZti0VhNL4NdSjkwmHsXvyggYb3ktcSGo5Ahly0N7sU4cqTCVdq1x/mbrz2uqK4IF57HeBcR4Cjzaq99iMFgLzRamyV+xP3vSmj3tFeZldmAdhbvjcpEForW/2zTXFxCl5cTT5TWJx0x4oqk7ehtgahPDPuF7FWe71Wajcw7sK9tqyQc/Gzg/MdAdrrVXvtBwoAeSfRWJSHUw+9axF6R3qZXHa/+BrLGpMaTw9W5WDB3etYCxMDyWsJw7fAuvsSU3sw4SwhCE0noSD9u3ZeTi/GF5STa6+bq0q8s4PznfW01GTWYrikhgJAtiVmYwxMDXhTRqHhEOGlMACBkoAXBsGmIB01HblXcPcmKwtul9FaIDy/AnPD7lh+sdtllFhrUdWS1gVmay1PwwteqV1i7fWBQJl3hfK5jgD1Faq9zgUKAEkpay1PZ54mBcLIgtt7X1NcQ3djtxcIh2sP52bB3etYCzNDCdtQQ64uO36bUSr3bGhDPQNF6bt4LxazPBid9XYYfTQQZnbZ1V4faar0Su1Ue529FACyo6y1DM0NJTWeDs0NAVBZWEl3U7d3tfLRwFEK8ny2d311BUZuJa8lTLq+J0w+NCduQz2b1m2oq9EYt17MeGcHibXXJ1uqvfWDs211lBb5KNizmAJAdt3w/PD6ovJIL09mngBQVlDmNZ4Gm4IcDxynMEPrGnbU/Pgnt6HGF94prV2/UK016G63WVqTlmGu1V5f7A9zuX+cvmdTrMbitdf7a70pozP7alR7naGyOgCy/X4A4owtjHnXIfSO9HoFdyX5JZxuOO3VV5xqOEVxvg/nnmNRGLuffM+EsXtA/Gey/siG220eTcs21PnlVUJPJ70Oo1svXO11aWE+Z+NbTs93BDi+V7XXmSKrA2CNzgByy8TSBFdHrnpnCQ8mH2CxFOUVcbJhveDudMPp3C64e5OlaXdmsHamMBSCBbf4TlFFvOcoYRtqReOuDzGx9vrio3Eejrra66rE2uvOeg41qvY6XRQAkvGml6fpG+3z1hDuTtwlZmMUmAKO1x/3rlY+03CGiqKKdA83Pax1aweJU0fDNyHmFm2pOZC8ltC8+9tQR2eXuNwf9radrtdeF8dvm+mmjPbVlSoQdokCQLLO3MqcC4T4lNHt8dusWldwd7TuqBcIXY1duV9w9yaRRXfP5cSpoxm3AE9+sbvncuLUUXXrri4wP59Y8ErtLvaHGUuovT7fEeBCp9ty2lSl2uudogCQrLcQWeD62HVvyujm+E0isQgGw+Haw0mNp7UltekebnpNDyXfbvNFH6y6i7+oaEoOhL1dUFS+K8Nytddz8QvSXO319KKrve5oKPdK7VR7nVoKAMk5S6tL3By/6c4Qhnu5Pnadpaj7JddZ0+mtIQSbg9SXZkZFQ9pEI5/chjox4I6ZfGg6tmEbaseu3G4zGrPcfTnjldp9/HiChXjt9bE9Vd4Oo7MH61R7vQ0KAMl5kWiEW+FbXn1F32gfC6tu/rmtqi3pauXm8uY0jzYDzIfjt9tc24baC8uuRI6SmuRKi5YetzV1h63VXrtSu3GuPp1iJRqjIM9wel+NV1vRvV+11+9CASC+sxpb5W74rne1ct9IH7MRt8++paIlqQK7paJFC5KxGIw/SF5LGL3D+jbUwwmBEITGY5C/s3+Vr9VeX3y0Vns9RcxCcUEewbZar7LiVEs1Baq9fi0FgPheNBblweSDpMbT6WV3g/nm8ub1grumIAeqDigQAJZm3PpB4t3VFsbdscIy2NudXGtR2bSjw5lZinDl8QSX4ltO7w27QE+svT7fUc97zZWqvU6gABDZIGZjPJp65E0ZhUZCTCxNANBQ2pC0htBe3a5AgPg21Cev2IbqFnKp3r+h5+gUFOzcRX3huWW+OzDhrSE83lB7vbaofNDntdcKAJG3sNbyeOaxFwa9w72MLo4CUFtcm7SGcKj2kL8aT98ksgTDN5Knjqafu2P5RdB8KnnXUc3+HduG+mIqofa6f5yXG2qvz3e6QNjrs9prBYDIO7LWMjg7mNR4+mL+BQBVRVV0N3V7N8k5UnfEfwV3bzLzMnkb6tBViN+ClPLGT25DLU79hX3WWp6EF7yzg8v9YSbitddtgTKv1M4PtdcKAJEUeDH3wuszCg2HeDb7DIDywnJXcBdfRzhef5zCPB8W3L1ONOIWlBPXEsKuCwqTB43HobVnfeoocCjl21DXaq8vPnKldh8NTHi11+81V3pTRp9ur6OqJLf+v1MAiOyA0YVRr+00NBJiYNrtrS8tKOV0w2lvDeFk/UmK8nVhU5KFieRtqIO9EF+Up7g6ORBaeqCsLqVPv1Z7vXbbzCtPJlhejddet9Z4pXbBA9lfe60AENkF4cXw+hnCSIiHk669tiiviFMNp7w1hFMNpygt8Nc89FvFYhB+mHwjndHbYGPueKBzQxvq8ZRuQ11ejdL3bMq7bea15672uig/j679NW5BuTPA6dbsq71WAIikwfTydNKU0f3J+67gLq+Ak/UnvSmjrsYu/zaevsny3IZtqB/D/Jg7VlgWb0ONB0JLEKr2pOyp55dXufJkwrWc9o9z+8VMUu31Wqndsb1VGV97rQAQyQCzK7PrBXfDvdwO3yZqo+SbfI4FjiUV3FUWVaZ7uJnHWph6lryW8PJ6wjbUfW66yNuGehoKU1MyN7WwwncHJrgcX1TeWHt9Ib7DqDMDa68VACIZaCGywLXRa96U0c3xm6zGXOPpkdoj3tbTnsYeakpq0j3czBRZctciJG1DdYvz5BW6SuzEeybUtqVkG+rozBKXB1yp3aWBcZ5PuF1O9RXF3vrBhc569tWl/8xOASCSBZZWl7gxdsMLhBtjN1iO30z+UO0hb8qop6lHBXdvMjsc334aX0sY6oWI64WirD55LaGlG4q3f7b1fGJ9y+mlhNrr1tpS7wrl8x0BGtNQe60AEMlCK9EVbo3f8tYQro1dYzG+n/5g9UHvOoRgc5DGst2/G1jWiK7C2N3kqaPxB/GDxvUaJV7BXH94W9tQ12qv10rtLveHmVlyW047Gyu8M4T32wPUlO387jAFgEgOiMQi3AnfSWo8nYu4ueh9lfuSCu72VuxN82gz3OJkfBtqQq3F0to21KrktYTW4La2oSbWXl985LacrtVeH99b5ZXafaqtjvIdqL1WAIjkoGgsyr3Je159xdWRq8ysuErnPeV7kgJhX+W+jFuczCixGEz0J6wlXIGRhG2ode3JgdB0AvK3dsHYyqqrvV4rtet7llx7faEjwLmOerr216Sk9jqrA8AY8wXgC52dnV9++PBhuocjkrFiNsbDyYfebTR7R3q9grvG0kZ6mnu8aaOD1QcVCG+zMg8vriWHwtyIO1ZQkrwNtfUsVG3trGtxxdVer9028+Yraq9/pLuFPdVbu3YkqwNgjc4ARN6NtZaB6YGkxtOxRbePvq6kLqnxtLOmUwV3b2MtTA9u2IZ6DaKuY4iqlg1tqKeh8N1/ac8sRfh4YMIrtbs3PMu3fv57ONm6tXtfKwBEBGstz2afeWEQGgkxPD8MQHVxNd2N3V4gHKk9Qn5edlcg7IrVZRi+lXyWMPXUHcsreMU21IPvvA01PLdMTVnRli84UwCIyCsNzQ2tB8JwiMG5QQAqCyvpauryzhKOBo6q4G6z5kaTF5eHrkLE3auAsoC7ajnxdpslVTs6HAWAiGzK8PxwUn3Fk5kngCu4W2s8DTYHORE4QeEWF0F9JxaF0fg21LVrE8buxQ8aaHgveeqo4Qik8OxLASAiWzK+OO5VV4RGQjyacjXOxfnFn2g8LSnY/YucstbiFLy4mnymsDjpjhVVugvUEqeOyrd+4Z8CQERSYnJpkqsjV701hPsT97FYCvMKXcFds7ta+UzDGRXcvQtrYWIgeS1h+BbYqDv+09+E9u/b0pdWAIjIjphZmaFvpM+bMro7cZeojVJgCjhWf8zbdtrV2EVFUerv/JXTVhbcLqPBK9D901Bau6UvowAQkV0xH5mnb7TP23p6K3zLK7h7r+49LxC6m7qpLt7atkZ5NwoAEUmLxdVFro9d93Ya3Ry7yUpsBYPxCu7Wpo3qSlJ71y9xFAAikhGWo8vcHLvprSFcH73OUnQJgI7qDi8Mgk1BGsoa0jza3KAAEJGMFIlGuB2+7QVC30gfC6uuuvlA1QGvAvts81may5vTPNrspAAQkaywGlvl3oQruOsd6aV3tJfZlVkAWipakuorWita1We0CQoAEclK0ViUh1MPvTWE3pFeppanAGgqa0qaMmqralMgvIICQERyQszGGJga8KaMQsMhwkthAOpL69fPEJqCtNe0q+AOBYCI5ChrLU9mnnhnB6HhECMLrrK5triW7qb1grvDtYd9GQgKABHxBWstg3ODSVNGQ3NDAFQWVdLT6O6pHGwO8l7dexTkpf4OXJlmswGQ+++EiOQ0Ywz7Kvexr3IfP3zohwF4Ofdy/QxhJMRfDf4VAOWF5ZxpPONNGR0PHPd1wZ3OAEQk540ujK73GQ2H6J/uB1zj6amGU14gnGw4SXF+cZpHu32aAhIReY2JpQnvFpqh4RAPJh9gsRTlFXGy4aS3hnC64TSlBVu7LWM6KQBERDZpenk6qfH03sQ9YjZGQV4BJwInvK2nXY1dlBeWp3u4b6UAEBHZormVOfpG+7xAuDN+h1W7Sr7J52jdUYLN8cbTpi6qinb27l5boQAQEUmRhcgC18aueVNGN8dvEolFMBiO1B3x1hB6mnqoKalJ93AVACIiO2VpdYmb4ze9rafXx66zHF0GoLOm04VBs7tArb5063f22ioFgIjILlmJrriCu3gg9I32sbi6CEBbVZs3ZRRsCtJU3rTj41EAiIikSSQW4W74rncdwtWRq8xF5gBorWhdD4TmIC0VLSl/fgWAiEiGiMai3J+8n9R4Or08DcCe8j1Jjaf7K/dvu+BOASAikqFiNsajqUdJ9RUTSxMANJQ2EGwK8nOnf46Omo4tfX1VQYiIZKg8k8fh2sMcrj3MTxz9Cay1PJ55vB4Iw727UmKnABARSTNjDO3V7bRXt/OjR36U3ZqZUQCIiGSY3brJjf+KskVEBFAAiIj4lgJARMSndm0NwBiTB/wboAoIWWt/b7eeW0REPmlTZwDGmN8xxowaY25t+PjnjDH3jTGPjDG/+JYv80WgFYgAg1sbroiIpMpmzwB+F/gq8PtrHzDG5AO/AXwW9wv9ijHmm0A+8JUNn/+zwBHgkrX2Pxtjvg78n+0NXUREtmNTAWCt/Y4xpm3Dhz8FPLLWDgAYY/4I+KK19ivAD278GsaYQWAl/s/o657LGPMB8AHA/v37NzM8ERHZgu2sAbQAzxP+PQh8+g2P/wbw68aY7wW+87oHWWs/BD4EMMaMGWOebnF89cD4Fj9XoBqYTvcgMohf3o9ce53Z+nq2O+4Dm3nQri0CW2sXgC+94+c0bPX5jDGhzXRhyKsZYz601n6Q7nFkCr+8H7n2OrP19ezWuLezDXQI2Jfw79b4xyQ3fCvdA8gwfnk/cu11Zuvr2ZVxb7oNNL4G8OfW2hPxfxcAD4Dvx/3ivwL8hLX29s4M9d3oDEBE5M02uw30a8Bl4IgxZtAY8yVr7Srw88C3gbvAn2TKL/+4D9M9ABGRTJbR9wMQEZGdoyoIERGfUgCIiPiU7gcgO8oYUw78R9xFgH9lrf3DNA8prfz0fuTSa82W1/Ku4/TNGYAxptwY83vGmP9ijPnJdI8nlV7X1ZRwfJ8x5v8aY+4YY24bY35hJ57rNd1QPwJ83Vr7ZeCHtvq8qRjfhsflG2P6jDF/vhPP94aerJS9H5t5rcaYGmPM140x94wxd40x51L9fKl6rZt8Pf8k/j18yxjzNWNMSSpfyxtez6ZfSzb9vGV1ALzjN+Su/iLaZb8LfO4Nx1eBf2qtPQa8D/wDY8yxxAcYYxqNMZUbPta52edK6Ib6O8Ax4Mfjz9HK+hXjr60ASaFXju8VfgG3e+0T3uG9eOXzveG9gNS+H5947lf4VeB/WmvfA06z4TVn2Gv9xNff8FwtwD8CgvHt6PnAj214TCZ8H2fNz1tWBwDv9g2527+Ido219jvAxBuOv7TWXo3/71ncL4GWDQ/7PuDPjDHFAMaYLwO//g7P5XVDWWtXgD/CNcAO4t572IXvt7e9FwDGmFbgB4Dfes1DNvVevOH5XvdeQArfj7e9VmNMNfA3gN+OP37FWju14WEZ81o38/8dbtq61LjrkMqAF1t5PTv5fZxNP29ZHQDv+A25q7+IMpVxF/R1AR8lftxa+6e4azr+OD5F9rPA33uHL/2qbqgWXAfU3zXG/CaZc1XmrwD/DIi96uAOvhewu+/HQWAM+K/x6a7fMm6O2JNNr9VaOwT8B+AZ8BKYttb+rw2Pyajv40z/ecvFReDXldT9GvBVY8wPkDm/iHaVMaYC+O/AP7bWzmw8bq3998a1uv4m0GGtndvuc1pr54Gf2e7XSRVjzA8Co9baXmPMZ173uJ14L+JfdzffjwKgG/iH1tqPjDG/Cvwi8C83jCkrXqsxphb3x9xBYAr4U2PMT1lr/2DD82bE93E2/Lz55i9ha+28tfZnrLV/P1NX8HeSMaYQ9834h9bab7zmMd8LnAD+B/BL7/gU2dINdQH4IWPME9zZ4d80xvzBxgflyHsxCAxaa9f++vw6LhCSZNFr/VvAY2vtmLU2gvtr9/zGB2XC68mWn7dcDIBM+eHLGMYYg5sHvmut/eXXPKYLV5/xRdxfEAFjzL99h6e5Ahwyxhw0xhThFue+ub2Rp5619p9ba1uttW24Mf6ltfanEh+TK++FtXYYeG6MORL/0PcDdxIfk2Wv9RnwvjGmLP49/f18clE77a8nq37erLVZ/R/QBtxK+HcBMIA7TSwCrgPH0z3OHX4PvoabE1273eaX4h//C2Av8D2ABW4A1+L/fX7D17gAnEz4dyHw5c0+V/zY53EFgf3Av8jE92LDYz+DKzjc+DU29V685fl2/L3YzGsFzgCh+P/3fwbUZupr3eTr+dfAPeAW8N+A4kz7Ps6mn7es7gIyrqTuM7ibv4wAv2St/W1jzOdxC335wO9Ya/9d2gYpIpKhsjoARERk63JxDUBERDZBASAi4lMKABERn1IAiIj4lAJARMSnFAAiIj6lABAR8SkFgIiITykARER86v8Dl1XM5elxlhMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for l in range(lmax+1):\n",
    "    plt.loglog(opt_eva[l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create compressed features, using labels naming consistent with Nigam 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel_thresh = 1.5e-6\n",
    "blocks = []\n",
    "for l in range(lmax+1):\n",
    "    n_selected = len(np.where(opt_eva[l] > sel_thresh)[0])    \n",
    "    de_block = density_expansion.block(center_species = 6, neighbor_species=1, spherical_harmonics_l = l)\n",
    "    block = Block(\n",
    "        values = de_block.values@opt_eve[l][:,:n_selected],\n",
    "        samples = de_block.samples,\n",
    "        components = Labels([\"m\"],np.asarray(range(-l,l+1), dtype=np.int32).reshape(-1,1)),\n",
    "        features = Labels([\"n\"], np.asarray([[n] for n in range(n_selected)], dtype=np.int32))\n",
    "    )\n",
    "    if rascal_hypers[\"compute_gradients\"]:\n",
    "        de_grad = de_block.gradient(\"positions\")\n",
    "        block.add_gradient(\n",
    "            parameter=\"positions\",\n",
    "            gradient_samples=de_grad[0],\n",
    "            gradient=de_grad[1]@opt_eve[l][:,:n_selected],\n",
    "        )\n",
    "    blocks.append( block )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_nu1 = Descriptor(sparse = Labels(names=[\"sigma\", \"lam\", \"nu\"], \n",
    "                                      values=np.asarray([[ 1, l, 1] for l in range(rascal_hypers[\"max_angular\"]+1)], dtype=np.int32)), \n",
    "                      blocks = blocks\n",
    "                     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "# multiplicity of the features, stored as a descriptor array\n",
    "acdc_full_nu1 = Descriptor(\n",
    "    sparse = acdc_nu1.sparse,\n",
    "    blocks = [\n",
    "        Block(\n",
    "            values = np.ones(shape=(1,1,len(acdc_nu1.block(s).features))),\n",
    "            samples = Labels([\"dummy\"], np.zeros(shape=(1,1), dtype=np.int32)),\n",
    "            components = Labels([\"dummy\"], np.zeros(shape=(1,1), dtype=np.int32)),\n",
    "            features = acdc_nu1.block(s).features\n",
    "        ) for s in acdc_nu1.sparse\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe43da791c0>]"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcUklEQVR4nO3de7yVc97/8dendioTSW2iKeV2GBmT2EzUTDFjIuMwcqaSTFJMmcYxQjGOiaEQNSNyrG50u7ln0M0jlN+uRzrIoWKQSs45hPj8/vhe+7bs9m6vtfe11rUO7+fjsR+tvdZ3res9V2veXa7D9zJ3R0RECl+jpAOIiEg8VOgiIkVChS4iUiRU6CIiRUKFLiJSJMqSWnCbNm28Y8eOSS1eRKQgzZ8//wN3L6/ptcQKvWPHjlRWVia1eBGRgmRm/67tNe1yEREpEip0EZEioUIXESkSKnQRkSKhQhcRKRJ1FrqZNTOzl8zsZTNbamZX1DCmqZk9aGbLzWyemXXMSlpg9Wro2RPWrMnWEqTU6TsmhSqdLfSvgYPdvQuwN3ComXWrNmYQ8LG77wKMB66NNWWKsWNhzhwYMyZbS5BSp++YFCrLZPpcM9sSmAOc5e7zUp7/H+Byd3/RzMqANUC5b+bDKyoqPJPz0Js3hw0bNn2+rAwuuijtjxGp1dVXw8aNmz7frBl89VXu84jUxMzmu3tFTa+ldWGRmTUG5gO7ABNSyzzSDngHwN03mtmnQGvgg2qfMxgYDNChQ4dM/jewciX85S8wfTp8880Pz2/cCFdemdFHidSo+uZH8+ZwzDFwww3J5BHJVFqF7u7fAXub2TbAf5rZz919SaYLc/dJwCQIW+iZvHeHHWDrrUOBN2sWSv3MM2HixExTiNTurLPgjjtCuW/YEL5zbdsmnUokPRmd5eLunwCzgUOrvbQKaA8Q7XJpCXwYQ74fWbsWhgyBuXPDnzpoJXFbuzaU+lFHhVJfuDDpRCLpq3ML3czKgW/d/RMzaw4cwqYHPR8DBgAvAscCz2xu/3l9zZz5w+MJE+L+dJEfvmNffQUHHgivvgpvvQWaR04KQTpb6DsAs81sEfD/gH+5+3+Z2RgzOzIaMxlobWbLgT8DF2YnrkhuNG8ejtd8/z0ce2zNB+RF8k1GZ7nEKdOzXESS8OijcPTR4XjN7bcnnUZk82e56EpRkc046ii44IJwoPSee5JOI7J5KnSROlx5JfTqFbbSFy9OOo1I7VToInUoK4P774dttoG+feGzz5JOJFIzFbpIGtq2hQcfDBe4DRy46UVIIvlAhS6Spl/9Cq65JpzaOH580mlENqVCF8nAyJHwhz/A+eeHCbxE8okKXSQDZvD3v0OnTnD88eHKUpF8oUIXyVDLluGio48/hpNOqnmGRpEkqNBF6qFLF7jtNpg9G0aPTjqNSKBCF6mn006DM84I86jPmpV0GhEVukiD3HILdO0K/fuHUxpFkqRCF2mAZs3C/nSA447TJF6SLBW6SAPtvDNMnQoLFsCf/pR0GillKnSRGBxxRLi37Z13wt13J51GSpUKXSQmY8bAQQeFu2m9/HLSaaQUqdBFYlI1iVerVuGmGJ9+mnQiKTUqdJEYbb89PPQQvPlmOK1Rk3hJLqnQRWLWowdcdx088giMG5d0GiklKnSRLDj33DB3+oUXwnPPJZ1GSoUKXSQLzGDKlHBK4wknwJo1SSeSUqBCF8mSrbeGGTPCwdETT9QkXpJ9KnSRLNprL7j9dnj2WbjkkqTTSLFToYtkWf/+MHgwXHstPPZY0mmkmKnQRXLg5pthn31Cua9YkXQaKVYqdJEcqJrEq1GjcNHRV18lnUiKkQpdJEc6dYJ77oGFC+Gcc5JOI8VIhS6SQ4cfDqNGweTJ4d6kInFSoYvk2BVXwG9+A0OHhq11kbio0EVyrHFjuO8+aN06XE36ySdJJ5JioUIXScB224VJvN5+W5N4SXxU6CIJOfBAuP56ePTR8KdIQ9VZ6GbW3sxmm9krZrbUzIbXMKalmc0ys5ejMQOzE1ekuAwfHu5FetFF4WpSkYZIZwt9IzDS3TsD3YBhZta52phhwCvu3gXoBYwzsy1iTSpShMzgrrtgl13CJF6rVyedSApZnYXu7qvdfUH0eD2wDGhXfRiwlZkZ0AL4iPAPgYjUoWoSr/XrNYmXNExG+9DNrCPQFZhX7aVbgT2A94DFwHB3/76G9w82s0ozq1y3bl39EosUoZ//HO64I8ydfvHFSaeRQpV2oZtZC2AGMMLdP6v2cm9gIbAjsDdwq5ltXf0z3H2Su1e4e0V5eXm9Q4sUo1NPDTeYvv76cLcjkUylVehm1oRQ5tPcfWYNQwYCMz1YDrwJ/Cy+mCKl4aaboKICBgyA5cuTTiOFJp2zXAyYDCxz9xtrGfY28Jto/PbA7sDKuEKKlIqmTcMkXmVlmsRLMpfOFnp3oB9wsJktjH76mNkQMxsSjRkLHGhmi4GngQvc/YMsZRYpajvtBPfeC4sWwbBhSaeRQlJW1wB3nwNYHWPeA34XVyiRUnfYYeEOR2PHQvfuMGhQ0omkEOhKUZE8ddllcMghYSt9wYKk00ghUKGL5KnGjWHaNCgvD/vTP/446USS71ToInmsvDxM4vXOO+HMl+83ubpD5AcqdJE8d8ABMG4czJoF112XdBrJZyp0kQJwzjlw/PHhbkezZyedRvKVCl2kAFRN4rXbbmG+l/feSzqR5CMVukiB2GqrcNHR55+HmRm//TbpRJJvVOgiBWTPPeHOO2HOnDCHukgqFbpIgTn55HBu+rhxMLOmmZWkZKnQRQrQuHGw//4wcCC88UbSaSRfqNBFClDTpvDww9CkCfTtC19+mXQiyQcqdJEC1aFDuJJ0yRIYOhTck04kSVOhixSw3r1h9Gi4++5wWqOUNhW6SIG79FL43e/g7LNh/vyk00iSVOgiBa5qEq/ttw+TeH30UdKJJCkqdJEi0KZNOEi6ahX0769JvEqVCl2kSPzyl3DjjfD443DNNUmnkSSo0EWKyLBhYa6XSy+Fp59OOo3kmgpdpIiYhakBdt8dTjop7IKR0qFCFykyLVrAjBnhYiNN4lVaVOgiRWiPPcJ56c8/DxdckHQayRUVukiROvHEcGOM8ePDtLtS/FToIkXshhugWzc4/XR4/fWk00i2qdBFitgWW4SbTDdtGibx+uKLpBNJNqnQRYpc+/Zw332wdCkMGaJJvIqZCl2kBBxyCFx+Odx7L9xxR9JpJFtU6CIl4pJL4NBDYfhwqKxMOo1kgwpdpEQ0ahS20Nu2DZN4ffhh0okkbip0kRLSunWYxOu996BfP03iVWxU6CIlZv/94aab4Ikn4K9/TTqNxKnOQjez9mY228xeMbOlZja8lnG9zGxhNObZ+KOKSFzOOgtOPjnc7eipp5JOI3FJZwt9IzDS3TsD3YBhZtY5dYCZbQNMBI509z2B4+IOKiLxMYNJk6Bz5zCJ17vvJp1I4lBnobv7andfED1eDywD2lUbdjIw093fjsa9H3dQEYnXT34SJvHasAGOPx6++SbpRNJQGe1DN7OOQFdgXrWXdgNamdn/mtl8M+tfy/sHm1mlmVWuW7euXoFFJD677w5TpsCLL8L55yedRhoq7UI3sxbADGCEu39W7eUyYF/gcKA3cKmZ7Vb9M9x9krtXuHtFeXl5A2KLSFyOOy6cm37zzWGaAClcaRW6mTUhlPk0d59Zw5B3gf9x9y/c/QPgOaBLfDFFJJuuuw4OOAAGDYJXX006jdRXOme5GDAZWObuN9Yy7FGgh5mVmdmWwC8J+9pFpABUTeLVrFm46EiTeBWmdLbQuwP9gIOj0xIXmlkfMxtiZkMA3H0Z8CSwCHgJuMvdl2QttYjE7qc/hfvvh1degcGDNYlXISqra4C7zwEsjXHXA9fHEUpEkvHb38KYMeEm0927w9ChSSeSTOhKURH5kYsvhj59YMQIeOmlpNNIJlToIvIjjRrBPffAjjuGM2A0iVfhUKGLyCa23Tbch3TNGjj1VE3iVShU6CJSo4qKcG76k0/ClVcmnUbSoUIXkVqdeWaYZvfyy+Gf/0w6jdRFhS4itTKD22+HPfcMszO+807SiWRzVOgisllbbhkm8frmm3CQVJN45S8VuojUabfdwiRe8+bBX/6SdBqpjQpdRNJy7LFw7rlwyy3wwANJp5GaqNBFJG3XXhuuID3jDFim2ZryjgpdRNLWpAk8+GC4OUbfvvD550knklQqdBHJSLt2YRKv116DP/5Rk3jlExW6iGTs4INh7NiwL33ChKTTSBUVuojUy4UXwu9/D3/+M8ydm3QaARW6iNRTo0YwdWrYBXP88fDBB0knEhW6iNRbq1ZhEq+1a+GUU+C775JOVNpU6CLSIPvuC7feGuZ6GTs26TSlTYUuIg12xhkwYEC429GTTyadpnSp0EWkwcxg4kTYa6+w6+Xtt5NOVJpU6CISiy23DPvTN24Mk3h9/XXSiUqPCl1EYrPrrvD3v4d7kY4cmXSa0qNCF5FYHXNMKPMJE+C++5JOU1pU6CISu6uvhh49wtQAS5cmnaZ0qNBFJHZVk3httVWYxGv9+qQTlQYVuohkxY47hrle3ngjnNaoSbyyT4UuIlnTqxdcdRU89FC4MYZklwpdRLLq/PPhiCPCgdIXX0w6TXFToYtIVjVqBHffDe3bh0m81q1LOlHxUqGLSNa1agUzZoQyP/lkTeKVLSp0EcmJrl3DuelPPQVXXJF0muKkQheRnBk0CAYODLMyPvFE0mmKT52FbmbtzWy2mb1iZkvNbPhmxu5nZhvN7Nh4Y4pIsZgwAbp0gVNPhX//O+k0xSWdLfSNwEh37wx0A4aZWefqg8ysMXAt8M94I4pIMWne/IdJvI49VpN4xanOQnf31e6+IHq8HlgGtKth6DnADOD9WBOKSNHZZZdw5ktlJZx7btJpikdG+9DNrCPQFZhX7fl2wB+A2+p4/2AzqzSzynU6d0mkpB19NJx3Htx2G0yblnSa4pB2oZtZC8IW+Ah3/6zayzcBF7j795v7DHef5O4V7l5RXl6ecVgRKS5//Sv8+tcweDAsWZJ0msKXVqGbWRNCmU9z95k1DKkAHjCzt4BjgYlmdnRcIUWkOJWVhfleqibx+qz6pqJkJJ2zXAyYDCxz9xtrGuPundy9o7t3BKYDQ939kTiDikhx2mGHMDPjihXhtEZN4lV/6Wyhdwf6AQeb2cLop4+ZDTGzIVnOJyIloGfPsPtl+nS4+eak0xSusroGuPscwNL9QHc/rSGBRKQ0nXcevPBC+HO//aB796QTFR5dKSoiecEM/vEP2GmnMInX+zoBOmMqdBHJG9tsEybx+ugjTeJVHyp0EckrXbrAxInw9NNw2WVJpyksKnQRyTsDB4YzXq66Ch5/POk0hUOFLiJ56ZZbYO+9oV8/eOutpNMUBhW6iOSl5s3D/vTvvw+TeG3YkHSi/KdCF5G8tfPOMHUqzJ8PI0YknSb/qdBFJK8deSRccAHccUcod6mdCl1E8t6VV0KvXjBkCCxenHSa/KVCF5G8V1YG998fzlPv2xc+/TTpRPlJhS4iBaFt2zCJ18qVcPrpmsSrJip0ESkYv/oVXHMNzJwJ48cnnSb/qNBFpKCMHAnHHAPnnw9z5iSdJr+o0EWkoJjBlCnQqVOYxGvt2qQT5Q8VuogUnJYtw0VHn3wCJ50EGzcmnSg/qNBFpCD94hfhBtOzZ8Po0UmnyQ8qdBEpWAMGwB//CFdfDbNmJZ0meSp0ESlof/sb7LMP9O8fTmksZSp0ESlozZqFe5ECHHdcaU/ipUIXkYLXqRPccw8sWAB/+lPSaZKjQheRovD738NFF8Gdd4Z7k5YiFbqIFI0xY+Cgg+Css+Dll5NOk3sqdBEpGlWTeLVqFW6KUWqTeKnQRaSobL89PPQQvPkmnHZaaU3ipUIXkaLTowdcfz088giMG5d0mtxRoYtIURoxIux2ufBCeO65pNPkhgpdRIqSGUyeDP/xH3DCCbBmTdKJsk+FLiJFa+utw0VHn34KJ55Y/JN4qdBFpKjttVe4wfSzz8IllySdJrtU6CJS9Pr1gzPPhGuvhcceSzpN9tRZ6GbW3sxmm9krZrbUzIbXMOYUM1tkZovN7AUz65KduCIi9XPTTbDvvmESrxUrkk6THelsoW8ERrp7Z6AbMMzMOlcb8ybQ0933AsYCk+KNKSLSMFWTeDVqFM5++eqrpBPFr85Cd/fV7r4gerweWAa0qzbmBXf/OPp1LvDTuIOKiDRUx45hEq+FC+Gcc5JOE7+M9qGbWUegKzBvM8MGAU/U8v7BZlZpZpXr1q3LZNEiIrE4/HAYNSqc0jhlStJp4mWe5nWxZtYCeBa4yt1n1jLmIGAi0MPdP9zc51VUVHhlZWWGcUVEGu6776B3b3j+eXjxRdh776QTpc/M5rt7RU2vpbWFbmZNgBnAtM2U+S+Au4Cj6ipzEZEkNW4M990HrVtD377hZtPFIJ2zXAyYDCxz9xtrGdMBmAn0c/fX440oIhK/7bYLk3i9/XbxTOKVzhZ6d6AfcLCZLYx++pjZEDMbEo0ZDbQGJkava1+KiOS9Aw+EG26ARx8Nk3kVurT3ocdN+9BFJB+4h2kBpk+HZ56Bnj2TTrR5Dd6HLiJSrMzgrrtg113DJF6rVyedqP5U6CJS8rbaCmbMgPXrC3sSLxW6iAiw554waVKYO/3ii5NOUz8qdBGRyCmnhBtMV93tqNCo0EVEUowfD/vtBwMGwPLlSafJjApdRCRF06bw8MNQVlZ4k3ip0EVEqtlpJ7j3Xli0CIYOLZyLjlToIiI1OOywcIejf/wjTORVCFToIiK1uOwyOOQQOPtsWLAg6TR1U6GLiNSicWOYNg3Ky8P+9I8/rvs9SVKhi4hsRnl5OEj67rvhzJfvv086Ue1U6CIidejWDcaNg1mz4Lrrkk5TOxW6iEgazj47TAswahTMnp10mpqp0EVE0mAGd94Ju+0Wiv2995JOtCkVuohImlq0CJN4ffFFmJnx22+TTvRjKnQRkQx07hy21OfMgYsuSjrNj6nQRUQydNJJMGxYOFA6s8a7LCdDhS4iUg/jxsH++8PAgfDGG0mnCVToIiL1UDWJV5Mm0LcvfPll0olU6CIi9dahQ7iSdMmSMI960pN4qdBFRBqgd28YPRqmTg0HS5OkQhcRaaBLLw3Ffs45MH9+cjlU6CIiDdS4cZg/ffvtwyReH32UTA4VuohIDNq0genTYdUq6N8/mUm8VOgiIjHZf/9wT9LHH4drrsn98lXoIiIxGjo0XHh06aXw9NO5XbYKXUQkRmYwaRLsvnso9lWrcrdsFbqISMyqJvH68svcTuKlQhcRyYI99gg3l37+ebjggtwsU4UuIpIlJ5wQzk0fPz6cAZNtKnQRkSy64YZwC7vTT4fXX8/usuosdDNrb2azzewVM1tqZsNrGGNm9jczW25mi8xsn+zEFREpLFtsAQ89FCbz6tsXVqyAnj1hzZr4l5XOFvpGYKS7dwa6AcPMrHO1MYcBu0Y/g4HbYk0pIlLA2reH++6DpUvhsMPCzTHGjIl/OXUWuruvdvcF0eP1wDKgXbVhRwFTPZgLbGNmO8SeVkSkQB15ZJiN8Y03wlWkt90WTnFs3jy+ZWS0D93MOgJdgXnVXmoHvJPy+7tsWvqY2WAzqzSzynXr1mUYVUSkcK1cGc5Lb9w4/L7llnDKKfDmm/EtI+1CN7MWwAxghLt/Vp+Fufskd69w94ry8vL6fISISEHaYQdo2TJspTdrBhs2wNZbQ9u28S0jrUI3syaEMp/m7jXdQW8V0D7l959Gz4mISGTtWhgyBObODX/GfWC0rK4BZmbAZGCZu99Yy7DHgLPN7AHgl8Cn7r46vpgiIoUv9YbSEybE//l1FjrQHegHLDazhdFzFwMdANz9duC/gT7AcuBLYGDsSUVEZLPqLHR3nwNYHWMcGBZXKBERyZyuFBURKRIqdBGRIqFCFxEpEip0EZEiYeF4ZgILNlsH/Lueb28DfBBjnLjkay7I32zKlRnlykwx5trJ3Wu8MjOxQm8IM6t094qkc1SXr7kgf7MpV2aUKzOllku7XEREioQKXUSkSBRqoU9KOkAt8jUX5G825cqMcmWmpHIV5D50ERHZVKFuoYuISDUqdBGRIpF3hW5mh5rZa9ENpy+s4fWmZvZg9Pq86C5KVa9dFD3/mpn1znGuP0c30l5kZk+b2U4pr31nZgujn8dynOs0M1uXsvwzUl4bYGZvRD8DcpxrfEqm183sk5TXsrm+ppjZ+2a2pJbXa73heZbXV125TonyLDazF8ysS8prb0XPLzSzyhzn6mVmn6b8fY1OeW2z34Es5zovJdOS6Du1bfRaVtaXmbU3s9lRDyw1s+E1jMnu98vd8+YHaAysAHYGtgBeBjpXGzMUuD16fCLwYPS4czS+KdAp+pzGOcx1ELBl9PisqlzR758nuL5OA26t4b3bAiujP1tFj1vlKle18ecAU7K9vqLP/jWwD7Ckltf7AE8QZhjtBszL9vpKM9eBVcsj3JR9XsprbwFtElpfvYD/auh3IO5c1cYeATyT7fUF7ADsEz3eCni9hv8/ZvX7lW9b6PsDy919pbt/AzxAuAF1qqOAu6PH04HfmJlFzz/g7l+7+5uEudn3z1Uud5/t7l9Gv84l3LUp29JZX7XpDfzL3T9y94+BfwGHJpTrJOD+mJa9We7+HPDRZobUdsPzbK6vOnO5+wvRciF336901ldtGvLdjDtXTr5f7r7a3RdEj9cDy9j03spZ/X7lW6Gnc7Pp/xvj7huBT4HWab43m7lSDSL8K1ylmYWbY881s6NjypRJrr7Rf95NN7OqWwXmxfqKdk11Ap5JeTpb6ysdtWXP5vrKVPXvlwP/NLP5ZjY4gTwHmNnLZvaEme0ZPZcX68vMtiQU44yUp7O+vizsCu4KzKv2Ula/X+ncsUgyYGanAhVAz5Snd3L3VWa2M/CMmS129xU5ijQLuN/dvzazMwn/dXNwjpadjhOB6e7+XcpzSa6vvGZmBxEKvUfK0z2i9bUd8C8zezXags2FBYS/r8/NrA/wCLBrjpadjiOA5909dWs+q+vLzFoQ/gEZ4e6fxfW56ci3LfR0bjb9f2PMrAxoCXyY5nuzmQsz+y0wCjjS3b+uet7dV0V/rgT+l/Avd05yufuHKVnuAvZN973ZzJXiRKr953AW11c6asue+I3QzewXhL/Do9z9w6rnU9bX+8B/Et+uxjq5+2fu/nn0+L+BJmbWhjxYX5HNfb9iX19m1oRQ5tPcfWYNQ7L7/Yr7wEADDyqUEQ4GdOKHAyl7VhszjB8fFH0oerwnPz4oupL4Doqmk6sr4SDQrtWebwU0jR63Ad4gpoNDaebaIeXxH4C5/sNBmDejfK2ix9vmKlc07meEA1SWi/WVsoyO1H6Q73B+fNDqpWyvrzRzdSAcFzqw2vM/AbZKefwCcGgOc7Wt+vsjFOPb0bpL6zuQrVzR6y0J+9l/kov1Ff3vngrctJkxWf1+xbZyY/xL6kM4OrwCGBU9N4aw1QvQDHg4+nK/BOyc8t5R0fteAw7Lca6ngLXAwujnsej5A4HF0Rd6MTAox7muBpZGy58N/CzlvadH63E5MDCXuaLfLweuqfa+bK+v+4HVwLeE/ZSDgCHAkOh1AyZEuRcDFTlaX3Xlugv4OOX7VRk9v3O0rl6O/p5H5TjX2Snfr7mk/INT03cgV7miMacRTpRIfV/W1hdhN5gDi1L+nvrk8vulS/9FRIpEvu1DFxGRelKhi4gUCRW6iEiRUKGLiBQJFbqISJFQoYuIFAkVuohIkfj/0TebzVaRMdsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ln_nu1 = np.asarray([[s[\"lam\"], b.values.shape[-1]] for s,b in acdc_nu1 ])\n",
    "plt.plot(ln_nu1[:,0], ln_nu1[:,1], 'b*-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ACDC utilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [],
   "source": [
    "cg = ClebschGordanReal(lmax*3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "default (all to all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 89.4 ms, sys: 497 µs, total: 89.9 ms\n",
      "Wall time: 88.6 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "acdc_nu2 = cg_increment(acdc_nu1, acdc_nu1, clebsch_gordan=cg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 60.1 ms, sys: 0 ns, total: 60.1 ms\n",
      "Wall time: 59.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "acdc_nu2_builder = cg_increment_builder(acdc_nu1, acdc_nu1, clebsch_gordan=cg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Labels([( 1, 0, 2), ( 1, 1, 2), ( 1, 2, 2), (-1, 1, 2), (-1, 2, 2),\n",
       "        ( 1, 3, 2), (-1, 3, 2), ( 1, 4, 2)],\n",
       "       dtype=[('sigma', '<i4'), ('lam', '<i4'), ('nu', '<i4')])"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acdc_nu2_builder.sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "for i in range(len(acdc_nu2.sparse)):\n",
    "    print(np.linalg.norm(\n",
    "        acdc_nu2.block(acdc_nu2.sparse[1]).values - acdc_nu2_builder.block(acdc_nu2.sparse[1]).values\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "specifying features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 74.2 ms, sys: 0 ns, total: 74.2 ms\n",
      "Wall time: 73.4 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "acdc_nu2 = cg_increment(acdc_nu1, acdc_nu1, acdc_full_nu1, acdc_full_nu1, clebsch_gordan=cg)\n",
    "acdc_full_nu2 = full_product_indices(acdc_full_nu1, acdc_full_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.02 s, sys: 89 ms, total: 1.11 s\n",
      "Wall time: 1.11 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "acdc_nu3 = cg_increment(acdc_nu2, acdc_nu1, acdc_full_nu2, acdc_full_nu1, clebsch_gordan = cg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.607837605262184e-12"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features_norm(acdc_nu3).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.607837605262195e-12"
      ]
     },
     "execution_count": 176,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(norm_full1**3).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## analyze jacobian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requires `compute_gradients=True`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "manually densify features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [],
   "source": [
    "j_nu1 = []\n",
    "for idx, block in acdc_nu1:\n",
    "    smp, grad = block.gradient(\"positions\")\n",
    "    j_nu1.append(grad.reshape((grad.shape[0],-1)))\n",
    "j_nu1 = np.concatenate(j_nu1, axis=-1).reshape(10,12,-1)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [],
   "source": [
    "j_nu2 = []\n",
    "for idx, block in acdc_nu2:\n",
    "    smp, grad = block.gradient(\"positions\")\n",
    "    j_nu2.append(grad.reshape((grad.shape[0],-1)))\n",
    "j_nu2 = np.concatenate(j_nu2, axis=-1).reshape(10,12,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "j_nu3 = []\n",
    "for idx, block in acdc_nu3:\n",
    "    smp, grad = block.gradient(\"positions\")\n",
    "    j_nu3.append(grad.reshape((grad.shape[0],-1)))\n",
    "j_nu3 = np.concatenate(j_nu3, axis=-1).reshape(10,12,-1)\n",
    "\n",
    "sv_nu3 = np.asarray([ sp.linalg.svdvals(j_nu3[i])  for i in range(10)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv_nu1 = np.asarray([ sp.linalg.svdvals(j_nu1[i])  for i in range(10)])\n",
    "sv_nu2 = np.asarray([ sp.linalg.svdvals(j_nu2[i])  for i in range(10)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'CN ($\\\\nu=2)$')"
      ]
     },
     "execution_count": 199,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEJCAYAAAByupuRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh90lEQVR4nO3deXTU9b3/8eebJCKbIrImQMEtaiEHBZcoaiAotHLhakFFEBeWIL/rFVsPLdLba70qv7Zayq3bJITFqwIKWK1YURAuS6kaBAoYf7aKSMoWBBcUISHv3x8ZKCKBECbzncn39TiHk8x3hu+8jkfmNZ/Pd/mYuyMiIuFUL+gAIiISHJWAiEiIqQREREJMJSAiEmIqARGREEsNOsDxat68uXfo0CHoGCIiSWXlypU73L3F4duTrgQ6dOhAUVFR0DFERJKKmW080nZNB4mIhJhKQEQkxFQCIiIhphIQEQkxlYCISIipBEREQkwlIJIsVqyACRMqf4rESNJdJyASSitWQG4u7NsHJ50ECxdCdnbQqaQO0EhAJBksXlxZAPv3V/5cvDjoRFJHqAREkkFOTuUIICWl8mdOTtCJpI6Iy3SQmU0B+gLb3b1TdNtvgH8B9gEfAre7+2fxyCOSdLKzK6eAFi+uLABNBUmMxGskMA3oc9i2N4BO7p4FfACMi1MWkeSUnQ3jxqkAJKbiUgLuvgTYedi21929PPrwL0DbeGQREZF/SpRjAncAf6rqSTMbaWZFZlZUWloax1giInVb4CVgZuOBcuDZql7j7vnu3s3du7Vo8Z3bYYuISA0Fep2Amd1G5QHjXHf3ILOIiIRRYCVgZn2AscBV7v51UDlERMIsLtNBZjYDWAFkmlmJmQ0DHgOaAG+Y2WozeyoeWURE5J/iMhJw90FH2FwYj/cWEZGqBX5gWEREgqMSEBEJMZWAiEiIqQREREJMJSAiEmIqARGREFMJiIiEmEpARCTEVAIiIiGmEhARCTGVgIhIiKkERERCTCUgIhJiKgERkRBTCYiIhJhKQEQkxFQCIiIhphIQEQkxlYCISIipBEREQkwlICISYioBEZEQUwmIiISYSkBEJMRUAiIiIaYSEBEJsbiUgJlNMbPtZrbukG3NzOwNM/tb9Odp8cgiIiL/FK+RwDSgz2HbfgYsdPezgYXRxyIiEkdxKQF3XwLsPGxzf2B69PfpwL/GI4uIiPxTkMcEWrn7lujvW4FWVb3QzEaaWZGZFZWWlsYnnYhICCTEgWF3d8CP8ny+u3dz924tWrSIYzIRkbotyBLYZmZtAKI/tweYRUQklIIsgZeBW6O/3wq8FGAWEZFQitcpojOAFUCmmZWY2TDg/wJXm9nfgF7RxyIiEkep8XgTdx9UxVO58Xh/ERE5soQ4MCwiIsFQCYiIhJhKQEQkxFQCIiIhphIQEQkxlYCISIipBEREQkwlICISYioBEZEQUwmIiISYSkBEJMRUAiIiIaYSEBEJMZWAiEiIqQREREJMJSAiEmIqARGREFMJiIiEmEpARCTEVAIiIiGmEhARCTGVgIhIiKkERERCTCUgIhJiKgERkRBTCYiIhFjgJWBm95jZejNbZ2YzzOzkoDOJiIRFoCVgZhnAvwPd3L0TkALcFGQmEZEwCXwkAKQCDcwsFWgIbA44j4hIQvnyyy8pLS2tlX0HWgLu/g/gEeATYAvwubu/fvjrzGykmRWZWVFt/YcQEUlEkyZNIj09nQceeKBW9n/cJWBmjcwsJRZvbmanAf2BjkA60MjMhhz+OnfPd/du7t6tRYsWsXhrEZGEtHv3biZPnszWrVsBaNeuHQMGDGDo0KG18n7HLAEzq2dmN5vZPDPbDrwPbDGz98zsN2Z21gm8fy9gg7uXunsZMBe47AT2JyKSlNasWcPo0aNJT09nxIgRzJ07F4Drr7+eqVOnctFFF9XK+6ZW4zWLgAXAOGCdu1cAmFkzoAfwKzN70d2fqcH7fwJcamYNgT1ALlBUg/2IiCSlffv20aNHD/785z9Tv359brjhBkaNGkV2dnZc3r86JdAr+i39W9x9JzAHmGNmaTV5c3d/y8xmA+8C5cAqIL8m+xIRSRbr1q1j6dKl3HnnnZx00kl07dqVgQMHMnToUJo1axbXLObux36R2blABvCWu+8+ZHsfd3+tFvN9R7du3byoSIMFEUkue/bs4YUXXiASifDnP/+ZBg0aUFJSErcPfTNb6e7dDt9enWMC/w68BNwFrDOz/oc8/XDsIoqI1E0LFiwgIyODW2+9lR07dvDII4/wySefxP1b/5FUZzpoBNDV3XebWQdgtpl1cPdJgNVqOhGRJLR3717mzJlDq1atyM3NpXPnzvTu3ZuRI0eSk5ODWeJ8dFanBOodmAJy94/NLIfKIvgeKgERkYM++OAD8vPzmTZtGp9++ik333wzubm5tGrVihkzZgQd74iqc53ANjPrcuBBtBD6As2BzrWUS0QkqYwaNYrMzEwmTZpEjx49eOONN/if//mfoGMdU3VGAkOpPHPnIHcvB4aaWaRWUomIJLgPP/yQwsJC7rvvPho3bswVV1xB+/btueOOO2jdunXQ8artmCXg7iVHeW55bOOIiCSusrIyXnrpJSKRCAsWLCAlJYWcnByuueYaBg8eHHS8GqnOSEBEJPS2b99OVlYW27Zto127djzwwAMMGzaM9PT0oKOdkOMqATPr6e5vHvhZW6FE4mbFCli8GHJyIDv7u48ltMrLy/njH//Ixo0bGTNmDC1btmTw4MH07NmTPn36kJISk1uoBa5aF4sdfLHZu+5+4YGftZirSrpYTGJmxQrIzYV9++Ckk+B3v4MxY/75eOFCFUEIffLJJxQUFFBYWMiWLVs4++yzee+990hNTe6JkxpfLFbV/k4wj0jwFi+u/MDfv7/y55w53368eHHQCSXOCgoK6NixIw899BAXXHABL730Up0ogKNJhEVlRIKRk1P5jT8lpfLnj3707cc5OUEnlFpWUlLCL3/5S5YvrzzH5corr2TcuHF89NFHzJs3j379+tXpAgAdGJYwy86unPI59BhA5846JlDH7d+/n/nz5xOJRHjllVdwd9LS0rj88svJzMzkwQcfDDpiXKkEJNyys7/9YX/4Y6lzunfvzl/+8hdatmzJ2LFjGTFiBGeccUbQsQJzvNNBB+4g+mWsg4iIxFpFRQXz589n+PDhlJdXXvM6atQonn/+eTZt2sSECRNCXQBwnGcHJQKdHSQix7Jt2zamTJlCQUEBGzZsoHnz5ixZsoTzzjsv6GiBqersIE0HiUidsnr1ai666CLKy8vJycnh4Ycf5rrrrqN+/fpBR0tIKgERSWqlpaVMmzaNtLQ0xowZQ1ZWFuPHj2fQoEFkZmYGHS/hVfuYgJk1MbPLzezM2gwkInIs7s7ixYsZNGgQbdu2ZezYsSxduhSAevXqcf/996sAqul4RgKLgDXAuWZWCtwRXWdYRCSu7r33Xn7729/StGlT7rzzTkaOHMn5558fdKykVO0SOPSAgpkNAF41s+vdfXOtJBMRofJb/7Jly4hEIvz0pz+lc+fODB48mKysLG644QYaNGgQdMSkdrw3kKsHtAE2Aq9ROTrQmEtEYm7Xrl08/fTT5Ofn895773HKKafQv39/OnfuzIUXXsiFFwZy+7I6p9olYGZbgQbAVmAzsAX4Yy3lEpEQKysrIzMzk9LSUi6++GIKCwu58cYbadSoUdDR6pzjGQmc4+5f1FoSEQmtzz//nGeeeYZly5bx3HPPkZaWxsSJE/n+979Ply5dgo5Xpx2zBMzMvFKVBXDgNbGNJiJ1mbvzzjvvEIlEmDlzJl9//TVdu3Zl165dNGvWLGlX6ko21TlFdJGZ3WVm7Q/daGYnmVlPM5sO3Fo78USkrnrppZe45JJLmDVrFoMHD6aoqIiioiKaNWsWdLRQqc50UB/gDmCGmXUEPgNOBlKA14HfufuqWksoInXCypUriUQidOnShdGjR9O7d2+eeuopBg0axCmnnBJ0vNCqzkLz3wBPAE+YWRrQHNjj7p/FIoCZNQUmA50Ap/L6gxWx2LeIBGv37t3MmDGDSCTCypUradCgwcE1eRs0aEBeXl7ACeW4ThF19zIqzwqKpUnAa+4+wMxOAhrGeP8Sdlo3ODCDBw/m5ZdfpnPnzjz22GMMGTKEU089NehYcohA7yJqZqcCq4EzqntgWXcRleNy+DrCWje41nz11VfMmjWLyZMn88ILL5CRkcHbb79NeXk52dnZmGlV2iDFeo3hWOkIlAJTzWyVmU02s++cCGxmI82syMyKSktL459Sktfh6whr3eCYW7duHXfddRcZGRkMGzaMzz77jJKSEgAuvvhiLrvsMhVAAgu6BFKBC4En3f0C4CvgZ4e/yN3z3b2bu3dr0aJFvDNKMjt8HWGtGxxT27dvp0uXLhQUFNC3b1+WLFnC+vXrueSSS4KOJtVUnesEfnGUp93d/+sE3r8EKHH3t6KPZ3OEEhCpsSOtIyw1VlxcTCQSYevWrcycOZOWLVvywgsvcOWVV3L66acHHU9qoDoHhr86wraGwHDgdKDGJeDuW81sk5lluvv/A3KB92q6P5Ej0rrBJ+Sbb75hzpw5RCIRli5dSlpaGgMHDqS8vJzU1FSuu+66oCPKCajOKaKPHvjdzJoAd1N53cBM4NGq/t5xuAt4Nnpm0EfA7THYp4icIHfHzHjyySf58Y9/zFlnncWvf/1rbrvtNjQtW3dU6xRRM2sG/BgYDEwHLnT3XbEI4O6rge8csRapNTpltEp79+7lxRdfJBKJMGzYMIYMGcItt9xC586d6dmzJ/XqBX0YUWKtOscEfgNcD+QDnd19d62nEqktOmX0iP7+97+Tn5/P1KlT2bFjBx07djz4gd+8eXN69eoVcEKpLdUZCfwE2Av8HBh/+Kle7q7rvSV5HOmU0ZCWwIHpHoCBAweydu1a+vXrx6hRo+jVq5e+9YdEdUogE2jl7ssO3Whm3alcW0AkeRw4ZfTASCCEp4x+9NFHFBQU8MILL7Bq1SqaNGlCQUEB6enpB2/pIOFRnRKYCIw7wvbPo8/9S0wTidSmkJ4yWlZWxiuvvEIkEuH111/HzOjbty87d+6kSZMmdOumw3JhVZ0SaOXuaw/f6O5rzaxD7COJ1LIQnTJaUVFBvXr1KC4u5vrrr6dt27b853/+J8OGDaNt27ZBx5MEUJ0SaHqU57TCs0iCKS8v59VXXyUSidC8eXOmT59OVlYWixYtonv37qSmHtd9I6WOq86RnyIzG3H4RjMbDqyMfSQRqYmSkhLuv/9+OnToQP/+/Vm1ahXnnHPOwedzcnJUAPId1fk/YgzwopkN5p8f+t2AkwBdKigSoP3792Nm1KtXj0mTJvHoo49yzTXX8Nhjj9G3b1996MsxVftW0mbWg8qFXwDWu/ubtZbqKHQraRHYvHkzhYWFTJ48mfz8fHr37s2WLVv45ptv6NixY9DxJAFVdSvpan9NcPdFwKKYphKRaquoqOCNN94gEonw8ssvs3//fnr16kWTJk0AaNOmTcAJJRlprCjxo9s11MjevXupX78+FRUVDBs2jH379vGTn/yEESNGcNZZZwUdT5KcSkDiQ7drOC4VFRW8+eabRCIR3n77bf7+97+TlpbG/PnzOeuss6hfv37QEaWO0HXhEh9a4ataduzYwa9//WvOOeccrr76ahYtWsTAgQPZs2cPAN///vdVABJTKgGJD63wVSV35+uvvwZg9erV/PSnPyUjI4Nnn32WkpISHnnkEU45Rbfoktqh6SCJj1jerqGOHFvYsWMH06dPJz8/n759+/Loo4/Ss2dPiouLOffcc4OOJyGhEpD4icXtGurAsYXly5fzxBNPMHv2bPbt28fll1/OZZddBkC9evVUABJXmg6S5JKkxxY+//zzg7/n5+czb948Ro4cydq1a1m2bBk/+tGPAkwnYaYSkOSSRMcW3J3ly5czdOhQWrduzapVqwD41a9+xebNm/n9739Pp06djrEXkdql6SBJLklwK+ivv/6aKVOmEIlEWLduHU2aNOH222+nadOmALRu3TrYgCKHqPZtIxKFbhshicjd2blzJ6effjq7d+8mPT2dzMxM8vLyuOmmm2jcuHHQESXkTvi2ESLyXV988QXPPvsskUiE8vJy1q5dS+PGjSkuLiYjIyPoeCLHpGMCIjWwfv16RowYQXp6OqNHj6ZevXrcddddVFRUAKgAJGloJCBSTV9++SUATZo0YfXq1Tz33HMMGjSIvLw8unXrdnDRdpFkopGAyDGsWrWKUaNGkZ6ezhNPPAHAgAED2Lx5M5MnT+aiiy5SAUjS0khAYqOOXMV7qGnTpvHEE0/wzjvv0KBBA2688UZ69eoFQP369XUPH6kTVAJy4urAVbwHbNy4ke9973sAzJw5k6+++opJkyZxyy23cNpppwWcTiT2EqIEzCwFKAL+4e59g84jx+lIV/EmUQl8/fXXPP/88wdv2/zxxx/Trl07Zs6cyamnnqqpHqnTEuWYwN1AcdAhpIaS6CreQ23ZsoW7776bjIwMbr/9dnbt2vWtO3Y2bdpUBSB1XuAjATNrC1wLPAT8OOA4UhNJcBXvAd988w3btm07OOUzefJk+vfvT15eHldeeaU+9CV0Ai8B4HfAWKBJVS8ws5HASID27dvHJ5Ucn1jcIbQWvf/+++Tn5zN9+nS6dOnCwoULadOmDdu2bdPVvBJqgU4HmVlfYLu7rzza69w93927uXu3Fi1axCmd1AWvvfYaOTk5nHfeeTz22GP06tWL8ePHH3xeBSBhF/RI4HKgn5n9EDgZOMXMnnH3IQHnkiT2t7/9jYyMDBo2bEhxcTGbNm1iwoQJ3H777bRq1SroeCIJJWFuIGdmOcC9xzo7SDeQkyPZt28ff/jDH4hEIrz55ptMnTqV2267jb1795KWlka9eolyDoRIMKq6gZz+ZYTdihUwYULlzyRUVlbGuHHjaNeuHTfeeCMffvghDz74IH369AEqL+pSAYhULejpoIPcfTGwOOAY4ZKkF3mVlZXx17/+la5du5KamsrChQvJzs4mLy+Pa665hpSUlKAjiiSNhCkBCUCSXeS1ceNGCgoKKCws5IsvvmDz5s2ceuqpLF++nLS0tKDjiSQljZPDLEku8lqzZg0//OEP6dixIxMmTKBbt27MmjXr4Jk9KgCRmtNIIMwS+CKvTZs2sW/fPs4880xSU1NZs2YN//Ef/8GwYcN0rYhIDCXM2UHVpbOD6q79+/fzpz/9iUgkwquvvsoNN9zAjBkzDj6nuX6RmtPykpLQfv/73/Ob3/yGTZs20bp1a8aNG8fw4cMPPq8CEKkdKgEJxP79+1m4cCG5ubmkpKSwZcsWzj33XCZOnEi/fv00zy8SJyoBiastW7YwZcoUCgoK2LhxI6+++io/+MEPeOihh3TzNpEA6OwgiYtPP/2UAQMG0L59e37+859z5plnMmvWLHJzcwFUACIB0UhAas327dspLi7mqquuomnTpmzYsIExY8YwcuRIzj777KDjiQgqAYkxd2fRokVEIhFefPFFmjVrRklJCampqRQVFekbv0iC0XSQxMy8efPIzMwkNzeXBQsW8G//9m8sXryY1NTK7xoqAJHEo5GA1Ji7s2TJEtq3b0/Hjh1p2LAhLVu25Be/+AUDBgzg5JNPDjqiiByDRgJy3Hbu3MnEiRM5//zzycnJ4fHHHwegR48eLFu2jCFDhqgARJKERgJyXEaNGsW0adPYu3cvl156KdOmTWPgwIFBxxKRGtJIIBnFcQ2AXbt28cwzzxx83KBBA4YNG8aaNWtYsWIFt956Kw0bNqz1HCJSOzQSSDZxWAPA3XnrrbeIRCLMmjWLPXv20KVLFzp16sTEiRNj+l4iEiyNBJLNkdYAiKH333+fLl26kJ2dzezZsxk6dCjvvvsunTp1iun7iEhi0Egg2RxYA+DASOAE1wBwd4qKivjss8+4+uqradeuHaeddhqRSIRBgwbRpEmTmMQWkcSkEkg2MVoD4Msvv+S5554jEomwatUqunbtytVXX02jRo1YHOPRhYgkLpVAMsrOPqHjAP/93//N+PHj2b17N1lZWTz++OMMHjw4hgFFJFmoBEJg9+7dzJw5k2uvvZY2bdrQrl07BgwYQF5eHpdccomu5BUJMR0YrsPWrFnD6NGjSU9PZ8SIEcydOxeA6667jqlTp3LppZeqAERCTiOBOqisrIwePXqwfPly6tevzw033EBeXh6XXXZZ0NFEJMGoBOqIdevWsWTJEkaPHk1aWhpdu3ZlwIABDB06lGbNmgUdT0QSlBaaT2J79uxh9uzZRCIRli9fzsknn0xJSQmnn3560NFEJMFUtdC8jgkkqYULF5KRkcHQoUPZvn37wUXaVQAicjwCnQ4ys3bA00ArwIF8d58UZKZEtXfvXubMmUOrVq3Izc2lU6dO9O7dmxEjRtCjRw8d4BWRGgn6mEA58BN3f9fMmgArzewNd38v4FwJ44MPPiA/P59p06bx6aefMmjQIHJzc2nVqhUzZswIOp6IJLlAS8DdtwBbor9/aWbFQAagEqDyts2RSITU1FT69+9PXl7ewYXZRURiIeiRwEFm1gG4AHjrCM+NBEYCtG/fPr7B4ujDDz+ksLCQ++67j8aNG3PFFVfQvn177rjjDlq3bh10PBGpgxLi7CAzawz8L/CQu8892mvr2tlBZWVlvPzyyzz11FMsWLCAlJQU5s2bR+/evYOOJiJ1SFVnBwU+EjCzNGAO8OyxCqCuKS0tJSsri61bt9KuXTseeOAB7rjjDjIyMoKOJiIhEfTZQQYUAsXu/tsgs8RDeXk5r7zyChs2bOCee+6hRYsW3HzzzfTs2ZM+ffqQkpISdEQRCZlAp4PMrDuwFFgLVEQ33+fur1b1d5JxOmjjxo1MnjyZwsJCtmzZwllnnUVxcTGpqYEPxEQkJBJyOsjdlwF1+gT3goIC8vLyAOjTpw9PPvkk1157rQpARBKCPolirKSkhMLCQnJzc+nevTtXXXUV9913H8OHD6dDhw5BxxMR+RaVQAzs37+f+fPnE4lEeOWVV6ioqCA1NZXu3btzzjnn8OCDDwYdUUTkiFQCMXDFFVewYsUKWrZsydixYxkxYgRnnHFG0LFERI5JN5A7ThUVFcyfP5/hw4dTXl4OVF7Z+/zzz7Np0yYmTJigAhCRpJEQF4sdj6DODtq2bRtTpkyhoKCADRs20Lx5c5YsWcJ5550X9ywiIscrIc8OSharV6/moosuory8nJycHB5++GGuu+466tevH3Q0EZETohI4gtLSUqZNm0ZaWhpjxowhKyuL8ePHc9NNN3HuuecGHU9EJGZ0TCDK3Vm0aBE33XQTGRkZjB07liVLlgBQr1497r//fhWAiNQ5KoGoe++9l549ezJ//nzuvPNO1q9fz9y5obqVkYiEUCing9ydZcuWEYlEGDt2LFlZWQwZMoSsrCwGDhxIw4YNg44oIhIXoSqBXbt28fTTTxOJRCguLuaUU06hX79+ZGVlccEFF3DBBRcEHVFEJK5CUwJlZWVkZmZSWlrKxRdfTGFhITfeeCONGjUKOpqISGBCUwJpaWlMnDiR888/X9/4RUSiQlMCAIMHDw46gohIQtHZQSIiIaYSEBEJMZWAiEiIqQREREJMJSAiEmIqARGREFMJiIiEmEpARCTEkm5lMTMrBTbW8K83B3bEME5tS6a8yZQVkitvMmWF5MqbTFnhxPJ+z91bHL4x6UrgRJhZ0ZGWV0tUyZQ3mbJCcuVNpqyQXHmTKSvUTl5NB4mIhJhKQEQkxMJWAvlBBzhOyZQ3mbJCcuVNpqyQXHmTKSvUQt5QHRMQEZFvC9tIQEREDqESEBEJsVCUgJlNMbPtZrYu6CzHYmbtzGyRmb1nZuvN7O6gMx2NmZ1sZm+b2Zpo3l8GnelYzCzFzFaZ2StBZzkWM/vYzNaa2WozKwo6z9GYWVMzm21m75tZsZllB52pKmaWGf1veuDPF2Y2JuhcVTGze6L/vtaZ2QwzOzlm+w7DMQEzuxLYDTzt7p2CznM0ZtYGaOPu75pZE2Al8K/u/l7A0Y7IzAxo5O67zSwNWAbc7e5/CThalczsx0A34BR37xt0nqMxs4+Bbu6e8Bc0mdl0YKm7Tzazk4CG7v5ZwLGOycxSgH8Al7h7TS9ErTVmlkHlv6vz3X2PmT0PvOru02Kx/1CMBNx9CbAz6BzV4e5b3P3d6O9fAsVARrCpquaVdkcfpkX/JOw3CzNrC1wLTA46S11iZqcCVwKFAO6+LxkKICoX+DARC+AQqUADM0sFGgKbY7XjUJRAsjKzDsAFwFsBRzmq6PTKamA78Ia7J3Le3wFjgYqAc1SXA6+b2UozGxl0mKPoCJQCU6NTbZPNrFHQoarpJmBG0CGq4u7/AB4BPgG2AJ+7++ux2r9KIEGZWWNgDjDG3b8IOs/RuPt+d+8CtAUuNrOEnHIzs77AdndfGXSW49Dd3S8EfgD8n+jUZiJKBS4EnnT3C4CvgJ8FG+nYotNW/YAXgs5SFTM7DehPZdGmA43MbEis9q8SSEDRufU5wLPuPjfoPNUVHf4vAvoEHKUqlwP9ovPsM4GeZvZMsJGOLvotEHffDrwIXBxsoiqVACWHjAJnU1kKie4HwLvuvi3oIEfRC9jg7qXuXgbMBS6L1c5VAgkmeqC1ECh2998GnedYzKyFmTWN/t4AuBp4P9BQVXD3ce7e1t07UDkF8Ka7x+wbVayZWaPoyQFEp1auARLyDDd33wpsMrPM6KZcICFPZjjMIBJ4KijqE+BSM2sY/XzIpfJYYUyEogTMbAawAsg0sxIzGxZ0pqO4HLiFym+pB05f+2HQoY6iDbDIzP4KvEPlMYGEP/UySbQClpnZGuBtYJ67vxZwpqO5C3g2+v9CF+DhYOMcXbRYr6bym3XCio6uZgPvAmup/NyO2e0jQnGKqIiIHFkoRgIiInJkKgERkRBTCYiIhJhKQEQkxFQCIiIhphIQEQkxlYCISIipBCRUzKy1mc00sw+jN2V71czOMTM3s0cPed29ZnZ/FftoYGb/G70FcSyzfWfdCzM7ycyWRO8eKRJzKgEJjegl9y8Ci939THfvCoyj8srcvcD1Zta8Gru6A5jr7vtjHHEah913yd33AQuBG2P8XiKASkDCpQdQ5u5PHdjg7mvcfSlQTuWl+PdUYz+DgZeg8j76Znbw5mPR0cWpNQl3lHUv/hB9T5GYUwlImHSicqW2qjwODD7ah3j01sNnuPvHAO7+OdDwkOmaNUDWYX9n6WFLGR7406uaudcBF1XztSLHRfOMIlHu/oWZPQ38O7Cnipc1Bz47bNtWKm+ktwk4N/r40P1ecYK59pvZPjNrEl1tTiRmVAISJuuBAcd4ze+ovFvj1Cqe3wMcvsj3ZiDdzC4Bdrj73w590syWAk2OsK973X3BsUJH1Qe+qeZrRapNJSBh8ibwsJmNdPd8ADPLAg5O/7j7zuhC3sOAKYfvwN13RZfTPNndD3wobwZ+SOUCJd+57feJjgTM7HQqy6XsRPYjciQ6JiCh4ZX3Tb8O6BU9RXQ9MIHDpm+AR6mc9qnK60D3Qx5vBm4GfuTuO2qa7yjrXvQA5tV0vyJHo/UERI6TmV0I3OPut8Tp/eYCP3P3D+LxfhIuGgmIHCd3f5fK1dRierHYkUTPRvqDCkBqi0YCIiIhppGAiEiIqQREREJMJSAiEmIqARGREFMJiIiEmEpARCTE/j/KcDKWU1TrrQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(\n",
    "    sv_nu1[:,0]/sv_nu1[:,-1],\n",
    "    sv_nu2[:,0]/sv_nu2[:,-1],\n",
    "    'r.'\n",
    ")\n",
    "plt.plot([1,8],[1,8],'k--')\n",
    "plt.xlabel(r\"CN ($\\nu=1)$\")\n",
    "plt.ylabel(r\"CN ($\\nu=2)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'CN ($\\\\nu=3)$')"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeiklEQVR4nO3dfZAU9b3v8feXXZ5R5ElYQViemREUZFGJiihoEWI0xJCE+FgaUcooRk9Zam7d5N4bE08ZRVMeEznq0eQYjDdK1BMT5RKQkNrDAQxG3OFBFAQFdgmHhxUEdvneP6Z33F13l1mY7d7Z/ryqpmanp6f7O5b0Z/rXv/79zN0REZF4ahd1ASIiEh2FgIhIjCkERERiTCEgIhJjCgERkRhTCIiIxFhhGDsxs07AMqBjsM/fufsPzWww8ALQC1gNXOvuh5vaVu/evb24uLiFKxYRaVtWr169y9371F8eSggAh4BL3L3SzNoDy83sj8BdwDx3f8HMfgncBPyiqQ0VFxezatWqlq9YRKQNMbMtDS0PpTnI0yqDl+2DhwOXAL8Llj8HfC2MekREJC20awJmVmBma4ByYBGwCdjj7lXBKtuA/mHVIyIiIYaAu1e7+1hgAHAOMCrbz5rZbDNbZWarKioqWqpEEZHYCb13kLvvAZYAE4FTzKzmusQA4ONGPjPf3UvcvaRPny9c1xARkeMUSgiYWR8zOyX4uzNwKZAiHQbfCFa7HngljHpERCQtrN5BRcBzZlZAOnhedPf/MLMy4AUz+zHwN+DpkOoRERFCCgF3/zswroHlH5C+PiAix1JaCkuXwuTJMHFi1NVIGxHWmYCInIjSUpgyBQ4fhg4dYPFiBYHkhIaNEMkHS5emA6C6Ov28dGnUFUkboRAQyQeTJ6fPAAoK0s+TJ0ddkbQRag4SyQcTJ6abgHRNQHJMISCSLyZO1MFfck7NQSIiMaYQEBGJMYWAiEiMKQRERGJMISAiEmMKARGRGFMIiIjEmEJARCTGFAIiIjGmEBARiTGFgIhIjCkERERiTCEgIhJjCgERkRhTCIiIxJhCQEQkxhQCIiIxphAQEYkxhYCISIwpBEREYkwhICISY6GEgJmdbmZLzKzMzN4zs7nB8h+Z2cdmtiZ4TA+jHhERSSsMaT9VwN3u/raZnQSsNrNFwXvz3P1nIdUhIiK1hBIC7r4d2B78vd/MUkD/MPYtIiKNC/2agJkVA+OAFcGi75nZ383sGTPrEXY9IiJxFmoImFk34CXgTnffB/wCGAqMJX2m8HAjn5ttZqvMbFVFRUVY5YqItHmhhYCZtScdAM+7+8sA7r7T3avd/Sjwr8A5DX3W3ee7e4m7l/Tp0yeskkVE2rywegcZ8DSQcvdHai0vqrXaDGBtGPWIiEhaWL2DzgeuBd41szXBsvuBWWY2FnBgM3BLSPWIiAjh9Q5aDlgDb70exv5FRKRhumNYRCTGFAIiIjGmEBARiTGFgIhIjCkERERiTCEgIhJjCgERkRhTCIiIxJhCQEQkxhQCIiIxphAQEYkxhYCISIwpBEREYkwhICISYwoBEZEYUwiIiMSYQkBEJMYUAiIiMaYQEBGJMYWAiEiMKQRERGJMISAiEmMKARGRGFMIiIjEmEJARCTGFAIiIjGmEBARibFQQsDMTjezJWZWZmbvmdncYHlPM1tkZhuD5x5h1CMiImlhnQlUAXe7exI4D7jNzJLAvcBidx8OLA5ei4hISEIJAXff7u5vB3/vB1JAf+BK4LlgteeAr4VRj4iIpIV+TcDMioFxwAqgr7tvD97aAfRt5DOzzWyVma2qqKgIp1ARkRgINQTMrBvwEnCnu++r/Z67O+ANfc7d57t7ibuX9OnTJ4RKRUTiIbQQMLP2pAPgeXd/OVi808yKgveLgPKw6hERkfB6BxnwNJBy90dqvfUqcH3w9/XAK2HUIyIiaYUh7ed84FrgXTNbEyy7H3gQeNHMbgK2AN8MqR4RESGkEHD35YA18vaUMGoQEZEv0h3DIiIxphAQEYkxhYCISIwpBEREYkwhICISYwoBEZEYa3YImFlXMytoiWJERCRcxwwBM2tnZt8xsz+YWTmwDtgezA3wkJkNa/kyRUSkJWRzJrAEGArcB/Rz99Pd/VTgAuA/gX82s2tasEYREWkh2dwxPNXdj9Rf6O67SQ8I91IwOJyIiOSZY4aAux8xs3PSf/rKYEawacA6d3+9Zp0WrlNERFrAMUPAzH4IfBkoNLNFwLmkm4juNbNx7v5AC9coIiItJJvmoG8AY4GOpGf/GuDu+8zsZ6RnB1MIiIjkqWwuDFe5e7W7HwA21cwI5u4HgaMtWp2IiLSobELgsJl1Cf4eX7PQzLqjEBARyWvZNAdNcvdDAO5e+6Dfns9nBRMRkTyUTe+gQ40s3wXsynlFIiISGo0dJCISY80KATO7pPaziIjkt+aeCfys3rOIiOSx420OamzSeBERySO6JiAiEmMKARGRGFMIiIjEWHNDoDJ43p/rQkREJHzNCgF3n1T7WURE8lsozUFm9oyZlZvZ2lrLfmRmH5vZmuAxPYxaRETkc2FdE3iW9EQ09c1z97HB4/WQahERkUDWIWBmJ5nZ+WY2tLk7cfdlwO7mfk5ERFpWc84ElgA3Ar8ys9+bWc8c7P97Zvb3oLmoR2MrmdlsM1tlZqsqKipysFsREYFmhIC7l7j7Te5+PvDvwOtmdtoJ7PsXwFDSs5ZtBx5uYt/zg/2X9OnT5wR2KSIitTV3ALl2ZtYf2AL8ifTZwXFx953BjGVHgX8FzjnebYmIyPHJZlIZAMxsB9CZ9DzDn5D+9f7a8e7YzIrcfXvwcgawtqn1RUQk97IOAWBEzfzCzWVmC4DJQG8z2wb8EJhsZmMBBzYDtxzPtkVE5PgdMwTMzDyt0QCoWaex9919VgOLn86yRhERaSHZXBNYYma3m9nA2gvNrIOZXWJmz6G5hqWtKC2Fn/40/SwSA9k0B00j3TV0gZkNBvYAnYAC4E3gUXf/W4tVKBKW0lKYMgUOH4YOHWDxYpg4MeqqRFpUNhPNfwY8ATxhZu2B3sBBd9/TwrWJhGvp0nQAVFenn5cuVQhIm9ecC8O4+xHSvYJE2p7Jk9NnADVnApMnR12RSItrVgiItGkTJ6abgJYuTQeAzgIkBhQCIrVNnKiDv8SKZhYTEYmxbO4T+J9NvO3u/n9yWI+IiIQom+agTxtY1gX4LtALUAiIiOSpbLqIZkb3NLOTgLmk7xt4gSZG/hTJC6WluhAssZbVheFg7oC7gKuB54Cz3f2/W7IwkRanm8NEjn1h2MweAlYC+4Ex7v4jBYC0CQ3dHCYSM9n0DrobOA34H8AnZrav9qNlyxNpQTU3hxUU6OYwia1smoNGAn3dfXnthWZ2Aem5BUTyk24OE8kqBOYB9zWwfG/w3ldzWpFImHRzmMRcNs1Bfd393foLg2XFOa9IRERCk00InNLEe51zVIeIiEQgmxBYZWY3119oZt8FVue+JBERCUs21wTuBBaa2dV8ftAvATqQniBeRETyVDZ3DO8EvmRmFwOjg8V/cPc/t2hlIiLS4rIeStrdlwBLWrAWEREJmYaSFhGJMYWAiEiMKQRERGJMISAiEmMKARGRGAslBMzsGTMrN7O1tZb1NLNFZrYxeO4RRi0iIvK5sM4EngWm1Vt2L7DY3YcDi4PXIiISolBCwN2XAbvrLb6S9CxlBM9fC6MWiaHSUvjpT9PPDb0WibGsbxZrAX3dfXvw9w6gb4S1SFtVfwrJRx+FO+/UlJIigVZxYdjdHfDG3jez2Wa2ysxWVVRUhFiZ5L36U0i+9JKmlBSpJcoQ2GlmRQDBc3ljK7r7fHcvcfeSPn36hFagtAH1p5C86ipNKSlSS5TNQa8C1wMPBs+vRFiLtFUNTSE5ZoymlBQJWLolpoV3YrYAmAz0BnYCPwR+D7wIDAS2AN909/oXj7+gpKTEV61a1WK1ioi0RWa22t1L6i8P5UzA3Wc18taUMPYvrURpqX6Bi7QyUTYHSZzU76WjXjkirUKr6B0kMVC/l4565Yi0CgoBCUf9XjrqlSPSKqg5SMLRUC8dEYmcQkDCM3GiDv4irYyag0REYkwhICISYwoBEZEY0zUBEZFWZOfOnaxdu5ZUKkVZWRmpVIp27dqxePHiFtmfQkBEJGRHjx5ly5YtmQP9hx9+yOOPP46Zcffdd/P8888D0L17dxKJBGPHjm2xWkIZOyiXNHaQiOSLI0eOsGnTJsrKypg2bRpdunThscce47777uPgwYOZ9fr27UsqlaJHjx6sWrWKvXv3kkwm6devH2aWk1oiHTtIRKQtO3jwIO3ataNjx46sWLGChx56iFQqxYYNG6iqqgJg5cqVlJSUcMYZZ3DrrbeSSCRIJpMkEgl69uyZ2VZJyReO0y1KISAi0gy7d+/mlVdeqdNm/+GHH7Jw4UKuvPJKDh48yLvvvksikeCKK67IHOhHjx4NwNSpU5k6dWrE3+JzCgERkVrcnYqKClKpVJ0D/XXXXce1115LRUUFN954Ix06dGDkyJFMmDCB6667jhEjRgAwefJk1q9fH/G3yJ5CQERiyd3Ztm1b5iA/aNAgZsyYwaeffkrfvp9Ped6tWzcSiQQ110+HDRvGhg0bGDx4MIWF+X8Izf9vICLShOrqaj744AMqKysZN24cAJdddhmlpaVUVlZm1vvOd77DjBkz6NatG08++SSDBg0imUwyYMCAOhdnCwoKGD58eOjfo6UoBESkTaiurqagoACAJ598kiVLllBWVsaGDRs4dOgQ48ePp6ZnYSKRYNSoUZn2+mQySe35y2fPnh3Jd4iCQkBE8s6GDRtYsWJFnTb7Tz/9lG3btgGwZMkSVq5cSSKRYNq0aSQSCcaMGZP5/GOPPRZV6a2OQkByI8ypIzVNZSz84x//yFycrXksXLiQTp068ctf/pJ58+ZRWFjI8OHDGTNmDIlEInM28Jvf/IZ27TQqTjYUAnLiwpw6UtNUtinuzvbt2zMH+W9+85uceuqpPP7449x+++2Z9Tp37syoUaPYtWsXAwYM4I477mD27NkMHTqU9u3bf2G7CoDsKQTkxDU0dWRLHZjD3JfkTM0wCd27d6dnz56sXLmSO+64g1Qqxd69ezPrDRs2jGnTpjFp0iQeeuihTJv9oEGD6hzYi4uLI/gWbZNCQE5czdSRNb/OW3LqyDD3Jcdt9+7dPPHEE5k2+/Xr13Pw4EHmz5/PzTffTLdu3ejcuTNXX311nTtn+/XrB8CZZ57JmWeeGfG3iAeNHSS5oWsCsXL06FHeeeedOhdmy8rKuOaaa/jBD37A7t276dWrF4MGDSKRSGQel1xyCUOHDo26/FjS2EHSssKcOlLTVIZm7969dQ70AwcOzLTVX3DBBRw4cICCggKGDh1KIpFgyJAhAPTs2ZPKykq6du0aZfmSBYWASMy5O+Xl5aRSKSorK7n88ssBuPDCC1m+fHlmvY4dO/Ktb30LSF94XbhwIaeddhrDhw+nY8eOX9iuAiA/qDlIJCbcnR07dlBUVATAvHnzeOmll0ilUuzevRuAgQMHsmXLFiDdl/6zzz7LtNcPHjw4czOW5J9W2xxkZpuB/UA1UNVQkdKC1L7eZq1evZpFixZlmnJSqRSHDx/mwIEDFBYWsmvXLtq1a8fMmTPrtNvXmDt3boTVS1giD4HAxe6+K+oiYkd97vPaoUOH2LBhQ52DfFlZGYsWLaJfv368+eab3H///fTv359kMslNN91EMpmkqqqKwsJCHnjggai/grQCrSUEJArqc58X9u/fz7p16zIH+RtvvJERI0bw29/+luuvvx4AM2PIkCEkEgk+/fRTAObMmcNtt93GySefHGX50sq1hhBw4E0zc+BJd58fdUGxoT73rUrNMAn9+/dn8ODBrFmzhiuuuIKtW7dm1mnfvj3nnXceI0aM4OKLL2bBggUkEglGjBhB586d62zvlFNOCfkbSD6K/MKwmfV394/N7FRgEXC7uy+rt85sYDbAwIEDx9dcuJIc0DWBULk7hw4dolOnTuzbt497770305xTXl4OwAMPPMD9999PeXk5d911V52RLocMGdLgMAkix9LYheHIQ6A2M/sRUOnuP2tsHfUOknzy+uuv895779Vps7/66qt54oknqKqqYsCAAQwZMqTOgf7ss8+uM6mJSC60yt5BZtYVaOfu+4O/LwP+d5Q1iTTHkSNHeP/99+vcNduvXz8eeeQRIN0u/9FHH9GvXz8SiQTXXHMNl156KQCFhYXs2LEjyvJFIr8m0BdYGMzaUwj8xt3/FG1JeUBNOKE7cOAA69evz/Sp/973vgfAJZdcUueGqkGDBnHZZZdlXv/xj3+kqKiIHj16hF6zSDZaVXNQNmLfHKRunS1q7969rF+/nnPOOQeABx98kPnz57N58+bMHLMnn3wye/bswcx4+eWXOXDgQGamKt0lK61Vq2wOkuOgbp0nrOZgbmb89a9/5YUXXsg05Wzfvh2AXbt20atXL3r06MG5557LDTfckGmzHzZsWGbO2a9//euRfQ+RXFAI5Bt162yWPXv2sGLFijpt9qlUiqVLlzJmzBjKysp49tlnSSaTXHbZZZkDfZcuXQC45ZZbuOWWWyL+FiItRyGQbyZOTDcB6ZpARlVVFR988EGdHji33norX/rSlygtLWX69OkA9OrVi2QyycyZMzMDnt1www1897vfzfyyF4kbhUA+iulQyrWHSRg2bBjjx49n48aNjB49msOHD2fW69+/P1/96lcBmDhxIm+99RaJRII+ffp8YZvqcy9xpxCQVmf//v1UVlZSVFTE4cOHmTlzJqlUik2bNnH06FEA7rzzTsaPH8/AgQOZO3cuyWSSZDLJqFGj6gyTcMoppzBp0qSovopIq6feQRK5X//617z99tuZppytW7cyc+ZMXnzxRQAuuugiTj311DrTEDY0TIKINE69gyQyO3bsYO3atXUuzJ588sm8+uqrQHrc+rKyMhKJBJMmTSKZTHLeeedlPv/WW29FVbpIm6cQkJw4evQomzdvrtPVsuau2VtvvZVXXnkFSDfPJBKJOvPMvvHGG/To0YN27dpFUrtInKk5SJql9jAJV1xxBe3bt+cnP/kJP/7xjzl48GBmvaKiIj744AM6depEaWkpBw8eJJlM0rdvX/XEEYmAmoOkWWpmn+rQoQPLli3LNNm8//77VFVVAbBu3TpGjhzJmWeeyZw5c+q02dceJmFiDHsyieQLhYBQXl7Oa6+9Vqef/ZYtW1i0aBFTpkxh3759vPfeeyQSCWbMmFFnzlmAyy+/PDM5uYjkFzUHxYC7U15eXmcawlQqxZw5c7jqqqtYvXo1JSUldOzYkZEjR2YO8rNmzWL48OFRly8iOaDmoBg4evQoW7duzRzkR44cyfTp09m5cydFRUWZ9U466SQSiQTV1dUAjBkzho0bNzJ48GAKCgqiKl9EIqAQyENVVVVs2rSJzz77jLPOOgt358ILL2TNmjWZ+WUBbr75ZqZPn07fvn15/PHHGTFiBIlEgv79+9e5ONuhQweGDRsWxVcRkYgpBFqx6urqzC/zn//85yxfvpyysjI2btzI4cOHueiii1i6dClmxhlnnEFJSQmJRCLzqBkmwcy47bbbovwqItJKKQRaibKyMlauXFnnhqqCggLWrVsHwJtvvsn69etJJBJ85StfIZFIcNZZZ2U+/+STT0ZVuojkMYVAiHbt2lXnIP/+++/z6quv0q5dO+bNm8dTTz1F+/btGTFiBGPHjmX06NG4O2aWWU9EJJfUOyjH3J1PPvkkc6C//vrr6d69Ow8++CD33XdfZr2uXbsyatQo3njjDXr16sWmTZuoqqpi6NChFBYqm0Ukt9Q7KMeqq6vZvHkzvXv3pnv37ixbtox77rmHVCrFvn37MuudffbZXHDBBUyZMoWHH3440/3y9NNPr/PLvvYwCiIiYVEIZGnnzp3Mnz8/8wt//fr1fPbZZyxYsIBvf/vbdO3ala5du3LttddmDvTJZJJTTz0VgAkTJjBhwoSIv4WISF1qDgpUV1fzzjvv1GmzLysrY86cOcydO5ePPvqIQYMGUVxcnDnIJxIJpk6dyqBBg3Jej4hILqk5KLBnz546B/qRI0dy8803U1VVxYQJEzh69CgFBQUMHz6cZDLJwIEDATj99NOprKyka9euEX8DEZHciVUIjBs3jjVr1mRed+rUiRtuuAGAjh078tprr1FcXMywYcPo0KFDnc+amQJARNqcWIXArFmzmDVrVqY5p7i4uM4wCTUTkouIxEWsQuCee+6JugQRkVZFdx+JiMRY5CFgZtPMbL2ZvW9m90Zdj4hInEQaAmZWAPwL8GUgCcwys2SUNYmIxEnUZwLnAO+7+wfufhh4Abgy4ppERGIj6hDoD2yt9XpbsExEREIQdQhkxcxmm9kqM1tVUVERdTkiIm1G1CHwMXB6rdcDgmV1uPt8dy9x95KaiVJEROTERR0CK4HhZjbYzDoA3wZejbgmEZHYiHwAOTObDjwKFADPuPsDx1i/AthynLvrDew6zs9GIZ/qzadaIb/qzadaIb/qzada4cTqHeTuX2hKiTwEwmRmqxoaRa+1yqd686lWyK9686lWyK9686lWaJl6o24OEhGRCCkERERiLG4hMD/qApopn+rNp1ohv+rNp1ohv+rNp1qhBeqN1TUBERGpK25nAiIiUkssQsDMnjGzcjNbG3Utx2Jmp5vZEjMrM7P3zGxu1DU1xcw6mdl/mdk7Qb3/K+qajsXMCszsb2b2H1HXcixmttnM3jWzNWaW+8m1c8jMTjGz35nZOjNLmdnEqGtqjJmNDP6b1jz2mdmdUdfVGDP7fvDva62ZLTCzTjnbdhyag8xsElAJ/MrdR0ddT1PMrAgocve3zewkYDXwNXcvi7i0BpmZAV3dvdLM2gPLgbnu/p8Rl9YoM7sLKAFOdvfLo66nKWa2GShx91bfl93MngP+4u5PBTd/dnH3PRGXdUzBaMYfA+e6+/Heg9RizKw/6X9XSXc/aGYvAq+7+7O52H4szgTcfRmwO+o6suHu29397eDv/UCKVjyonqdVBi/bB49W+8vCzAYAXwGeirqWtsTMugOTgKcB3P1wPgRAYAqwqTUGQC2FQGczKwS6AJ/kasOxCIF8ZWbFwDhgRcSlNCloXlkDlAOL3L011/socA9wNOI6suXAm2a22sxmR11MEwYDFcC/BU1tT5lZ16iLytK3gQVRF9EYd/8Y+BnwEbAd2Ovub+Zq+wqBVsrMugEvAXe6+76o62mKu1e7+1jSAwCeY2atssnNzC4Hyt19ddS1NMMF7n426YmXbguaNlujQuBs4BfuPg74FGj1MwUGzVZXAP836loaY2Y9SM+zMhg4DehqZtfkavsKgVYoaFt/CXje3V+Oup5sBaf/S4BpEZfSmPOBK4J29heAS8zs36MtqWnBr0DcvRxYSHoiptZoG7Ct1lng70iHQmv3ZeBtd98ZdSFNmAp86O4V7n4EeBn4Uq42rhBoZYILrU8DKXd/JOp6jsXM+pjZKcHfnYFLgXWRFtUId7/P3Qe4ezHpJoA/u3vOflHlmpl1DToHEDStXAa0yh5u7r4D2GpmI4NFU4BW2Zmhnlm04qagwEfAeWbWJTg+TCF9rTAnYhECZrYAKAVGmtk2M7sp6pqacD5wLelfqTXd16ZHXVQTioAlZvZ30kODL3L3Vt/1Mk/0BZab2TvAfwF/cPc/RVxTU24Hng/+XxgL/CTacpoWBOulpH9Zt1rB2dXvgLeBd0kft3N253AsuoiKiEjDYnEmICIiDVMIiIjEmEJARCTGFAIiIjGmEBARiTGFgIhIjCkERERiTCEgsWJm/czsBTPbFAzK9rqZjTAzN7OHa633T2b2o0a20dnM3gqGIM5lbV+Y98LMOpjZsmD0SJGcUwhIbAS33C8Elrr7UHcfD9xH+s7cQ8DXzax3Fpu6EXjZ3atzXOKz1Bt3yd0PA4uBb+V4XyKAQkDi5WLgiLv/smaBu7/j7n8Bqkjfiv/9LLZzNfAKpMfRN7PM4GPB2UX34ymuiXkvfh/sUyTnFAISJ6NJz9TWmH8Brm7qIB4MPTzE3TcDuPteoEut5pp3gDPrfeYv9aYyrHlMzbLutcCELNcVaRa1M4oE3H2fmf0KuAM42MhqvYE99ZbtID2Q3lZgVPC69nYvPMG6qs3ssJmdFMw2J5IzCgGJk/eAbxxjnUdJj9b4b428fxCoP8n3J8BpZnYusMvdN9Z+08z+ApzUwLb+yd3/37GKDnQEPstyXZGsKQQkTv4M/MTMZrv7fAAzOxPINP+4++5gIu+bgGfqb8Dd/zuYTrOTu9cclD8BppOeoOQLw36f6JmAmfUiHS5HTmQ7Ig3RNQGJDU+Pmz4DmBp0EX0P+Cn1mm+Ah0k3+zTmTeCCWq8/Ab4DXOXuu463vibmvbgY+MPxblekKZpPQKSZzOxs4Pvufm1I+3sZuNfdN4SxP4kXnQmINJO7v016NrWc3izWkKA30u8VANJSdCYgIhJjOhMQEYkxhYCISIwpBEREYkwhICISYwoBEZEYUwiIiMSYQkBEJMb+P97Z2GmzW4H7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(\n",
    "    sv_nu1[:,0]/sv_nu1[:,-1],\n",
    "    sv_nu3[:,0]/sv_nu3[:,-1],\n",
    "    'r.'\n",
    ")\n",
    "plt.plot([1,8],[1,8],'k--')\n",
    "plt.xlabel(r\"CN ($\\nu=1)$\")\n",
    "plt.ylabel(r\"CN ($\\nu=3)$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Independent feature tracking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This performs the cg iteration while avoiding to compute duplicate / linearly dependent features, and tracks the multiplicity so that no density (correlation) is lost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 8.14 s, sys: 8.48 s, total: 16.6 s\n",
      "Wall time: 2.67 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "cg = ClebschGordanReal(12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.98 s, sys: 2.22 s, total: 4.2 s\n",
      "Wall time: 562 ms\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "cg = ClebschGordanReal(lmax*5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_full1 = features_norm(acdc_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_sel_nu2 = canonical_indices(acdc_full_nu1, acdc_full_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 68.6 ms, sys: 63 µs, total: 68.7 ms\n",
      "Wall time: 73.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "canonical_acdc_nu2 = cg_increment(acdc_nu1, acdc_nu1, acdc_full_nu1, acdc_full_nu1, clebsch_gordan=cg, M=acdc_sel_nu2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_sel_nu3 = canonical_indices(acdc_sel_nu2, acdc_full_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 448 ms, sys: 19 µs, total: 448 ms\n",
      "Wall time: 459 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "canonical_acdc_nu3 = cg_increment_builder(canonical_acdc_nu2, acdc_nu1, acdc_sel_nu2, acdc_full_nu1, clebsch_gordan=cg, \n",
    "                                    M=acdc_sel_nu3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_sel_nu4 = canonical_indices(acdc_sel_nu3, acdc_full_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Labels([(-1, 0, 4), ( 1, 0, 4), (-1, 1, 4), ( 1, 1, 4), (-1, 2, 4),\n",
       "        ( 1, 2, 4), (-1, 3, 4), ( 1, 3, 4), (-1, 4, 4), ( 1, 4, 4),\n",
       "        (-1, 5, 4), ( 1, 5, 4), (-1, 6, 4), ( 1, 6, 4), (-1, 7, 4),\n",
       "        ( 1, 7, 4), ( 1, 8, 4)],\n",
       "       dtype=[('sigma', '<i4'), ('lam', '<i4'), ('nu', '<i4')])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acdc_sel_nu4.sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1,0,4) in acdc_sel_nu4.sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 674 ms, sys: 6.38 ms, total: 680 ms\n",
      "Wall time: 676 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "canonical_acdc_nu4 = cg_increment(canonical_acdc_nu3, acdc_nu1, acdc_sel_nu3, acdc_full_nu1, clebsch_gordan=cg, \n",
    "                                    M=acdc_sel_nu4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_sel_nu5 = canonical_indices(acdc_sel_nu4, acdc_full_nu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.89 s, sys: 311 ms, total: 5.2 s\n",
      "Wall time: 5.19 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "canonical_acdc_nu5 = cg_increment(canonical_acdc_nu4, acdc_nu1, acdc_sel_nu4, acdc_full_nu1, clebsch_gordan=cg, \n",
    "                                    M=acdc_sel_nu5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "compute norms and show that they match the norm of nu=1 features to the power nu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnu2_full = features_count(acdc_nu2)\n",
    "nnu2_can = features_count(canonical_acdc_nu2)\n",
    "norm_full2 = features_norm(acdc_nu2)\n",
    "norm_can2 = features_norm(canonical_acdc_nu2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(122, 3.873834622434604e-06) (72, 3.873834622434604e-06) 3.873834622434606e-06\n"
     ]
    }
   ],
   "source": [
    "print((nnu2_full, norm_full2.sum()), (nnu2_can, norm_can2.sum()), (norm_full1**2).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnu3_full = features_count(acdc_nu3)\n",
    "nnu3_can = features_count(canonical_acdc_nu3)\n",
    "norm_full3 = features_norm(acdc_nu3)\n",
    "norm_can3 = features_norm(canonical_acdc_nu3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2204, 9.607837605262185e-10) (564, 9.60783760526219e-10) 9.607837605262195e-10\n"
     ]
    }
   ],
   "source": [
    "print((nnu3_full, norm_full3.sum()), (nnu3_can, norm_can3.sum()),(norm_full1**3).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnu4_can = features_count(canonical_acdc_nu4)\n",
    "norm_can4 = features_norm(canonical_acdc_nu4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4020, 2.6749188416148777e-13) 2.6749188416148813e-13\n"
     ]
    }
   ],
   "source": [
    "print((nnu4_can, norm_can4.sum()),(norm_full1**4).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnu5_can = features_count(canonical_acdc_nu5)\n",
    "norm_can5 = features_norm(canonical_acdc_nu5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(26568, 8.126066561228681e-17) 8.126066561228692e-17\n"
     ]
    }
   ],
   "source": [
    "print((nnu5_can, norm_can5.sum()),(norm_full1**5).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterative scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_thresh_nu2, thrnorm_nu2 = threshold_indices(acdc_nu1, acdc_nu1, acdc_full_nu1, acdc_full_nu1, \n",
    "                                                      sel_threshold=1e-9, l_threshold=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_thresh_nu2, thrnorm_nu2 = threshold_indices(acdc_nu1, acdc_nu1, sel_threshold=1e-9, l_threshold=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_thr_nu2 = cg_increment(acdc_nu1, acdc_nu1, acdc_full_nu1, acdc_full_nu1, M=acdc_thresh_nu2, clebsch_gordan=cg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_thr2 = features_norm(acdc_thr_nu2)\n",
    "nnu2_thr = features_count(acdc_thr_nu2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( norm_thr2.mean(), norm_thr2.mean()+thrnorm_nu2, norm_full2.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nnu2_full, nnu2_can, nnu2_thr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_zip_nu2, acdc_w_nu2 = compress_features(acdc_thr_nu2, 1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_zip2 = features_norm(acdc_zip_nu2)\n",
    "nnu2_zip = features_count(acdc_zip_nu2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(norm_zip2.mean(), norm_thr2.mean(), norm_full2.mean())\n",
    "print(nnu2_full, nnu2_can, nnu2_thr, nnu2_zip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## zipo New heading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "recompute the compressed features!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_zip_post_nu2 = cg_increment(acdc_nu1, acdc_nu1,  M=acdc_w_nu2, clebsch_gordan=cg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_norm(acdc_zip_post_nu2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_norm(acdc_zip_nu2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_thresh_nu3, thrnorm_nu3 = threshold_indices(acdc_zip_nu2, acdc_nu1, sel_threshold=1e-15, l_threshold=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_nu1.block(1).values.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_zip_nu2.block(1).values.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_full_nu1.block(0).values.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = acdc_full_nu2.block(0).features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(0,0,0,0) in acdc_full_nu2.block(0).features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_full_nu2.sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.int32((0,0,0,0))in  features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acdc_thr_nu3 = cg_increment(acdc_zip_nu2, acdc_nu1, M=acdc_thresh_nu3)"
   ]
  }
 ],
 "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "397px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}