AAA_CMB_Optmized-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pylab as plt\n",
    "import operator\n",
    "import os\n",
    "import warnings\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "\n",
    "\n",
    "from PIL import Image\n",
    "from pyshtools.legendre import legendre\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "from lib3 import *\n",
    "from parameters import *\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import mpmath as mp\n",
    "from enum import Enum\n",
    "from time import sleep\n",
    "from numba import jit, prange\n",
    "\n",
    "params = {'legend.fontsize': 'x-large',\n",
    "          'figure.figsize': (10, 5),\n",
    "          'axes.labelsize': 'x-large',\n",
    "          'axes.titlesize': 'x-large',\n",
    "          'xtick.labelsize': 'large',\n",
    "          'ytick.labelsize': 'large'}\n",
    "plt.rcParams.update(params)\n",
    "\n",
    "\n",
    "class Color(Enum):\n",
    "    MINIMIZEPOSITION = 1\n",
    "    FINDNEIGHBORHOOD = 2\n",
    "    EVALUATE_DF_AT_POSITION = 3\n",
    "    FINDBESTFORKRANGE = 4\n",
    "    OPTIMIZE_SPECTRUM = 5\n",
    "    CREATE_GAUSSIAN_BACKGROUND = 6\n",
    "\n",
    "\n",
    "def get_chopped_range(lmax, n=20, lmin=1):\n",
    "    llist = sorted(list( { int(np.round(x,0)) for x in np.linspace(lmin,lmax, n) }))\n",
    "    return llist\n",
    "\n",
    "\n",
    "@jit(parallel=True)\n",
    "def eval_my_gegenbauer(aa, bb, cosksi):\n",
    "    A = np.zeros(cosksi.shape)\n",
    "    for i in prange(len(cosksi)):\n",
    "        # A[i] = sp.eval_gegenbauer(aa, bb, cosksi[i])\n",
    "        A[i] = mp.gegenbauer(aa, bb, cosksi[i])\n",
    "    return A.astype(float)\n",
    "\n",
    "def sinxx(xx):\n",
    "    y = np.zeros(xx.shape)\n",
    "    for i, x in enumerate(xx):\n",
    "        if x == 0:\n",
    "            y[i] = 1.0\n",
    "        else:\n",
    "            y[i] = np.sin(x) / x\n",
    "    return y\n",
    "\n",
    "def optimzeSpectraC(df, x0, opt=False, kkk=None):\n",
    "    def erroC(x, df):\n",
    "        k = df.k - 2\n",
    "        dk = x[0] * k + x[1] * k ** 2 + x[2] * k ** 3\n",
    "        dl = x[3] * df.l + x[4] * df.l ** 2\n",
    "        dm = x[5] * np.abs(df.m) + x[6] * df.m ** 2\n",
    "        df[\"population\"] = x[7] + dk + dl + dm\n",
    "        dfG = pd.DataFrame(df.groupby(['k', \"l\"])[\"CL\", \"population\"].apply(lambda x: x.astype(float).mean()),\n",
    "                           columns=[\"CL\", \"population\"])\n",
    "        dfG[\"CL\"] /= dfG[\"CL\"].max()\n",
    "        err = np.sum((dfG.population - dfG.CL) ** 2)\n",
    "        print(err)\n",
    "        return err\n",
    "\n",
    "    if opt:\n",
    "        if kkk:\n",
    "            df = df[(df.k == kkk[0]) * (df.l == kkk[1])]\n",
    "        xout = minimize(erroC, x0, args=(df), method='nelder-mead', options={'xatol': 1e-8, 'disp': True})\n",
    "        print(\"errC\", xout.fun, xout.x)\n",
    "        return xout.x, df\n",
    "    else:\n",
    "        x = x0\n",
    "    ##########################################################\n",
    "    if kkk:\n",
    "        df = df[(df.k == kkk[0]) * (df.l == kkk[1])]\n",
    "        plt.plot(df.m, df.coeff)\n",
    "        plt.show()\n",
    "    k = df.k - 2\n",
    "    dk = x[0] * k + x[1] * k ** 2 + x[2] * k ** 3\n",
    "    dl = x[3] * df.l + x[4] * df.l ** 2\n",
    "    dm = x[5] * np.abs(df.m) + x[6] * df.m ** 2\n",
    "    df[\"population\"] = x[7] + dk + dl + dm\n",
    "    dfG = pd.DataFrame(df.groupby(['k', \"l\"])[\"CL\", \"population\"].apply(lambda x: x.astype(float).mean()),\n",
    "                       columns=[\"CL\", \"population\"])\n",
    "    dfG[\"CL\"] /= dfG[\"CL\"].max()\n",
    "    dfG[\"difference\"] = (dfG.CL - dfG.population) / dfG.population\n",
    "    dfG[\"zero\"] = 0.0\n",
    "    dfG[\"one\"] = 1.0\n",
    "    ##########################################################\n",
    "    ax = dfG.CL.plot(legend=False)\n",
    "    # dfG.population.plot( ax=ax)\n",
    "\n",
    "    plt.xlabel('Hyperspherical Harmonic k')\n",
    "    plt.title('Power Spectrum C')\n",
    "    plt.ylabel(\"Hyperspherical Harmonic Mode Population\")\n",
    "    # plt.ylim(0, 2)\n",
    "    plt.show()\n",
    "\n",
    "def get_legendre(a):\n",
    "    ii = a[0]\n",
    "    lrange = a[1]\n",
    "    chunck = a[2]\n",
    "    lmax = np.max(lrange)\n",
    "    pp = np.zeros([len(chunck), len(lrange), lmax + 1])\n",
    "    for i, x in enumerate(chunck):\n",
    "        pp[i, :, :] = legendre(lmax, x, csphase=-1)[lrange]\n",
    "    return ii, pp\n",
    "\n",
    "def get_P(lrange, costheta):\n",
    "    lmax = np.max(lrange)\n",
    "    pp = np.zeros([0, len(lrange), lmax + 1])\n",
    "    lencos = len(costheta)\n",
    "    number_of_workers = 10\n",
    "    P = {}\n",
    "    mychunks = chunks(costheta, lencos // number_of_workers, lrange)\n",
    "    with ProcessPoolExecutor(max_workers=number_of_workers) as executor:\n",
    "        responses = executor.map(get_legendre, mychunks)\n",
    "    for response in responses:\n",
    "        i = response[0]\n",
    "        P[i] = response[1]\n",
    "    for key in np.sort(list(P.keys())):\n",
    "        pp = np.concatenate([pp, P[key]], axis=0)\n",
    "    return pp\n",
    "\n",
    "def functionG(a):\n",
    "    l = a[0][0] - 1\n",
    "    k = a[0][1] + l\n",
    "    cosksi = a[1]\n",
    "    A = np.array([float(mp.gegenbauer(1 + l, k - l, x)) for x in cosksi])\n",
    "    print(\"ended {}_{}\".format(k, l))\n",
    "    return (a[0], A)\n",
    "\n",
    "def get_gegenbauerAsyn(cosksi,extra_G_keys):\n",
    "    number_of_workers = 10\n",
    "    multiindex = pd.MultiIndex.from_tuples(extra_G_keys, names=('k', 'l'))\n",
    "    columns = list(np.arange(len(cosksi)))\n",
    "    if len(extra_G_keys) == 0:\n",
    "        return\n",
    "    extraG = pd.DataFrame(data=np.zeros([len(multiindex), len(columns)]), index=multiindex, columns=columns)\n",
    "    extra_G_keys = ((x,cosksi) for x in extra_G_keys)\n",
    "    with ProcessPoolExecutor(max_workers=number_of_workers) as executor:\n",
    "        responses = executor.map(functionG, extra_G_keys )\n",
    "    for response in responses:\n",
    "        extraG.loc[response[0], :] = response[1]\n",
    "    return extraG\n",
    "\n",
    "def plot_aitoff(fcolors, cmap=cm.RdBu_r):\n",
    "    hp.mollview(fcolors.squeeze(), cmap=cmap)\n",
    "    hp.graticule()\n",
    "    plt.show()\n",
    "\n",
    "def plot_aitoff_df(l, m, df, cmap=cm.RdBu_r):\n",
    "    def plot_orth(fcolors):\n",
    "        sigma = np.std(fcolors)\n",
    "        hp.orthview(fcolors, min=-2 * sigma, max=2 * sigma, title='Raw WMAP data', unit=r'$\\Delta$T (mK)')\n",
    "        plt.show()\n",
    "\n",
    "    def get_LM(l, m, df):\n",
    "        if m >= 0:\n",
    "            ind = (l + 1) * (l + 2) - l + m\n",
    "            return df[0, ind, :]\n",
    "        if m < 0:\n",
    "            ind = (l + 1) * (l + 2) - l + m\n",
    "            return df[1, ind, :]\n",
    "\n",
    "    fcolors = get_LM(l, m, df)\n",
    "    plot_aitoff(fcolors, cmap=cmap)\n",
    "    plot_orth(fcolors)\n",
    "\n",
    "def B_l(beam_arcmin, ls1):\n",
    "    theta_fwhm = ((beam_arcmin / 60.0) / 180) * math.pi\n",
    "    theta_s = theta_fwhm / (math.sqrt(8 * math.log(2)))\n",
    "    return np.exp(-2 * (ls1 + 0.5) ** 2 * (math.sin(theta_s / 2.0)) ** 2)\n",
    "\n",
    "def is_invertible(a):\n",
    "    return a.shape[0] == a.shape[1] and np.linalg.matrix_rank(a) == a.shape[0]\n",
    "\n",
    "def chunks(lst, n, args=None):\n",
    "    \"\"\"Yield successive n-sized chunks from lst.\"\"\"\n",
    "    maxn = len(lst)\n",
    "    for i in range(0, maxn, n):\n",
    "        nextstep = i + n\n",
    "        if nextstep > maxn:\n",
    "            nextstep = maxn\n",
    "        if args is not None:\n",
    "            yield (i, args, lst[i:nextstep])\n",
    "        else:\n",
    "            yield (i, lst[i:nextstep])\n",
    "\n",
    "def chg2ang(alpha):\n",
    "    residue = alpha % (2 * np.pi)\n",
    "    return int(residue * 180 / np.pi)\n",
    "\n",
    "def newerr(x, SMICA_LR, newmap):\n",
    "    newnewmap = SMICA_LR.squeeze() - x[1] * newmap.squeeze() + x[0]\n",
    "    err = np.std(newnewmap)\n",
    "    return err\n",
    "\n",
    "\n",
    "class HYPER:\n",
    "    def __init__(self, nside3D, sigma_smica, planck_IQU_SMICA, karray, mypath,\n",
    "                 lambda_k=0.0, lambda_l=0.0, lambda_m=0.0, bandwidth=256,\n",
    "                 loadpriorG=False, savePG=False, longG=False):\n",
    "        self.G = pd.DataFrame()\n",
    "        self.loadpriorG=loadpriorG\n",
    "        self.savePG=savePG\n",
    "        self.extra_G_keys={}\n",
    "        self.extra_P_keys = {}\n",
    "        self.longG = longG\n",
    "        self.nside3D = nside3D\n",
    "        self.bandwidth = bandwidth\n",
    "        self.df = np.zeros([1, 1])\n",
    "        self.results = {}\n",
    "        self.sigma_smica = sigma_smica\n",
    "        self.xx, self.yy, self.zz = hp.pix2vec(nside=nside3D, ipix=np.arange(hp.nside2npix(nside3D)))\n",
    "        self.k1 = min(karray)\n",
    "        self.kmax = max(karray)\n",
    "        self.karray = karray\n",
    "        self.mypath = mypath\n",
    "        self.G_filename = os.path.join(self.mypath, \"G_{}_{}_{}_{}.npy\".format(self.nside3D,\n",
    "                                                                                  chg2ang(lambda_k),\n",
    "                                                                                  chg2ang(lambda_l),\n",
    "                                                                                  chg2ang(lambda_m)))\n",
    "        self.P_filename = os.path.join(self.mypath, \"P_{}_{}_{}_{}_\".format(self.nside3D,\n",
    "                                                                               chg2ang(lambda_k),\n",
    "                                                                               chg2ang(lambda_l),\n",
    "                                                                               chg2ang(lambda_m)))\n",
    "        if loadpriorG:\n",
    "            if os.path.exists(self.G_filename):\n",
    "                self.G = pd.read_pickle(self.G_filename)\n",
    "        self.sinksi = np.zeros([1, 1])\n",
    "        self.cosksi = np.zeros([1, 1])\n",
    "        self.phi = np.zeros([1, 1])\n",
    "        self.x = np.zeros([1, 1])\n",
    "        self.y = np.zeros([1, 1])\n",
    "        self.z = np.zeros([1, 1])\n",
    "        self.costheta = np.zeros([1,1])\n",
    "        self.lambda_k, self.lambda_l, self.lambda_m = (lambda_k, lambda_l, lambda_m)\n",
    "        self.SMICA = self.normalizeFColors(planck_IQU_SMICA, sigma_smica)\n",
    "        self.SMICA_LR = self.SMICA\n",
    "        self.newmap = np.zeros(self.SMICA.shape)\n",
    "        self.change_SMICA_resolution(nside3D, doit=True, bandwidth=bandwidth)\n",
    "        self.change_HSH_center(self.lambda_k, self.lambda_l, self.lambda_m, self.karray,\n",
    "                               self.nside3D, loadpriorG=loadpriorG, doit=True, savePG=savePG)\n",
    "\n",
    "    def change_SMICA_resolution(self, nside3D, doit=False, bandwidth=256):\n",
    "        if (self.nside3D == nside3D) and not doit:\n",
    "            return\n",
    "        filename = \"./img1/SMICA_{}_{}.png\".format(nside3D, bandwidth)\n",
    "        title = \"SMICA_nside3D={}_{}\".format(nside3D, bandwidth)\n",
    "        self.SMICA_LR = self.change_resolution(self.SMICA, nside3D=nside3D, bandwidth=bandwidth,\n",
    "                                               title=title, filename=filename,\n",
    "                                               plotme=True, save=True)\n",
    "        mu, sigma = norm.fit(self.SMICA_LR)\n",
    "        self.SMICA_LR = (self.SMICA_LR - mu) / sigma * self.sigma_smica\n",
    "\n",
    "    def factorMP(self, k, l, m):\n",
    "        m = np.abs(m)\n",
    "        a = (-1) ** m * np.sqrt((2 * l + 1) / (2 * np.pi) * mp.factorial(l - m) / mp.factorial(l + m))\n",
    "        b = (-1) ** k * np.sqrt(\n",
    "            2 * (k + 1) / np.pi * mp.factorial(k - l) / mp.factorial(k + l + 1) * 2 ** (2 * l) * mp.factorial(l) ** 2)\n",
    "        c = float64(a * b)\n",
    "        return c\n",
    "\n",
    "    def normalizeFColors(self, fcolors, sigma_smica):\n",
    "        (mu, sigma) = norm.fit(fcolors)\n",
    "        if sigma == 0.0:\n",
    "            return fcolors\n",
    "        return sigma_smica / sigma * (fcolors - mu)\n",
    "\n",
    "    def getSpectralFiles(self, mypath):\n",
    "        return [join(mypath, f) for f in sorted(listdir(mypath)) if isfile(join(mypath, f))]\n",
    "\n",
    "    def creategiff(self, kmax, nside3D, mypath=\"./img1\", prefix=\"aitoff_\", filename=\"CMB\"):\n",
    "        giffiles = self.getSpectralFiles(mypath)\n",
    "        giffiles = [x for x in giffiles if prefix in x]\n",
    "        images = []\n",
    "        for ff in giffiles:\n",
    "            with Image.open(ff).convert(\"RGB\") as img:\n",
    "                images.append(img)\n",
    "        fname = os.path.join(mypath, filename) + '_{}_{}.gif'.format(kmax, nside3D)\n",
    "        images[0].save(fname, save_all=True, append_images=images[1:], optimize=False, duration=1000, loop=0)\n",
    "\n",
    "    def creategiffMem(self, x):\n",
    "        # err, prime, results, fig\n",
    "        images = [y[3] for y in x]\n",
    "        titles = [y[1] for y in x]\n",
    "        errs = [y[0] for y in x]\n",
    "        for i, xx in enumerate(images):\n",
    "            (lk, ll, lm) = [x for x in titles[i]]\n",
    "            err = errs[i]\n",
    "            filename = \"./img1/aitoff_{}_{}_{}_{}.png\".format(kmax, chg2ang(lk), chg2ang(ll), chg2ang(lm))\n",
    "            # draw the canvas, cache the render\n",
    "            xx.seek(0)\n",
    "            im = Image.open(xx)\n",
    "            im.save(filename)\n",
    "\n",
    "    def change_resolution(self, fcolors, nside3D=1024, filename=None,\n",
    "                          title=None, save=False, plotme=None, bandwidth=256):\n",
    "        # noinspection PyUnresolvedReferences\n",
    "        SMICA_alms = hp.map2alm(fcolors, lmax=bandwidth)\n",
    "        fcolors = hp.alm2map(SMICA_alms, nside=nside3D)\n",
    "        (mu, sigma) = norm.fit(fcolors)\n",
    "        if sigma != 0.0:\n",
    "            fcolors = (fcolors - mu) / sigma * self.sigma_smica\n",
    "        fcolors = np.expand_dims(fcolors, axis=1)\n",
    "        self.plot_aitoff(fcolors.squeeze(), kk=0, ll=0, mm=0, err=0, filename=filename, title=title, save=save,\n",
    "                         plotme=plotme)\n",
    "        return fcolors\n",
    "\n",
    "    def plotHistogram(self, fcolors, nside3D, kmax, plotme=False):\n",
    "        (mu, sigma) = norm.fit(fcolors)\n",
    "        n, bins, patch = plt.hist(fcolors, 600, density=1, facecolor=\"r\", alpha=0.25)\n",
    "        y = norm.pdf(bins, mu, sigma)\n",
    "        plt.plot(bins, y)\n",
    "        plt.xlim(mu - 4 * sigma, mu + 4 * sigma)\n",
    "        plt.xlabel(\"Temperature/K\")\n",
    "        plt.ylabel(\"Frequency\")\n",
    "        plt.title(r\"Histogram {}_{} HU Modeling of Planck SMICA Map \".format(kmax, nside3D), y=1.08)\n",
    "        plt.savefig(\"./img1/Histogram_{}_{}.png\".format(kmax, nside3D), dpi=300)\n",
    "        if plotme:\n",
    "            plt.show()\n",
    "\n",
    "    def plot_aitoff(self, fcolors, kk, ll, mm, err, filename=None, title=None, plotme=False, save=False, nosigma=True):\n",
    "        # noinspection PyUnresolvedReferences\n",
    "        plt.clf()\n",
    "        if title is None:\n",
    "            title = \"{}_{}_{}_{}\".format(kmax, chg2ang(kk), chg2ang(ll), chg2ang(mm))\n",
    "\n",
    "        f = plt.figure()\n",
    "        if nosigma:\n",
    "            hp.mollview(fcolors.squeeze(), title=title, min=-2 * self.sigma_smica,\n",
    "                        max=2 * self.sigma_smica, unit=\"K\", cmap=cm.RdBu_r)\n",
    "        else:\n",
    "            mu, sigma = norm.fit(fcolors.squeeze())\n",
    "            hp.mollview(fcolors.squeeze(), title=title, min=mu - 2 * sigma,\n",
    "                        max=mu + 2 * sigma, unit=\"K\", cmap=cm.RdBu_r)\n",
    "        hp.graticule()\n",
    "        if save:\n",
    "            plt.savefig(filename, format='png')\n",
    "        if plotme:\n",
    "            plt.show()\n",
    "        f.clear()\n",
    "        plt.close(f)\n",
    "\n",
    "    def pickledict(self, test_dict, filename='./errarray.pickle'):\n",
    "        import pickle\n",
    "        with open(filename, 'wb') as handle:\n",
    "            pickle.dump(test_dict, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
    "\n",
    "    def unpickledict(self, filename='./errarray.pickle'):\n",
    "        import pickle\n",
    "        with open(filename, 'rb') as handle:\n",
    "            text_dict = pickle.load(handle)\n",
    "        return text_dict\n",
    "\n",
    "    def plot_Single_CL_From_Image(self, fcolors, nside,ymin=1E-6, xmax=300, log=False):\n",
    "        dl_SMICA, ell = self.get_dl(fcolors, nside)\n",
    "        dl_SMICA = dl_SMICA / np.max(dl_SMICA)\n",
    "        fig = plt.figure(figsize=[12, 12])\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.plot(ell, dl_SMICA)\n",
    "        ax.set_xlabel(\"$ell(ell+1)C_ell/2\\pi \\, \\,(\\mu K^2)$\")\n",
    "        ax.set_ylabel('$\\ell$')\n",
    "        ax.set_title(\"Angular Power Spectra\")\n",
    "        ax.legend(loc=\"upper right\")\n",
    "        if log:\n",
    "            ax.set_yscale(\"log\")\n",
    "        ax.set_xlim(1, xmax)\n",
    "        ax.set_ylim(ymin, 1)\n",
    "        ax.grid()\n",
    "        plt.show()\n",
    "        return cl_SMICA, dl_SMICA, ell\n",
    "\n",
    "    def plot_CL_From_Image(self, fcolors, nside, planck_theory_cl, xmax=3000, ymax=1.0):\n",
    "        dl_SMICA, ell = self.get_dl(fcolors, nside)\n",
    "        dl_SMICA = dl_SMICA / np.max(dl_SMICA)\n",
    "        planck_theory_cl[:, 1] = planck_theory_cl[:, 1] / np.max(planck_theory_cl[10:3000, 1])\n",
    "        fig = plt.figure(figsize=[12, 12])\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.plot(planck_theory_cl[:, 0], planck_theory_cl[:, 1])\n",
    "        ax1 = plt.twinx(ax)\n",
    "        ax1.plot(ell, dl_SMICA)\n",
    "        ax.set_xlabel(\"$ell(ell+1)C_ell/2\\pi \\, \\,(\\mu K^2)$\")\n",
    "        ax.set_ylabel('$\\ell$')\n",
    "        ax.set_title(\"Angular Power Spectra\")\n",
    "        ax.legend(loc=\"upper right\")\n",
    "        # ax1.set_yscale(\"log\")\n",
    "        ax.set_xlim(10, xmax)\n",
    "        ax1.set_xlim(10, xmax)\n",
    "        # ax1.set_ylim(1E-4, ymax)\n",
    "        ax1.grid()\n",
    "        plt.show()\n",
    "        return cl_SMICA, dl_SMICA, ell\n",
    "\n",
    "    def get_dl(self, fcolors, nside):\n",
    "        cl_SMICA = hp.anafast(fcolors)\n",
    "        ell = np.arange(len(cl_SMICA))\n",
    "        pl = hp.sphtfunc.pixwin(nside=nside)\n",
    "        # Deconvolve the beam and the pixel window function\n",
    "        dl_SMICA = cl_SMICA / (B_l(10.0, ell) ** 2 * pl ** 2)\n",
    "        dl_SMICA = (ell * (ell + 1) * dl_SMICA / (2 * math.pi))\n",
    "        return dl_SMICA, ell\n",
    "\n",
    "    def plot_ONLY_CL_From_Image(self, fcolors, nside, smica, nsidesmica, xmax=30):\n",
    "        dl_SMICA_HU, ell = self.get_dl(fcolors.squeeze(), nside=nside)\n",
    "        dl_SMICA, ell1 = self.get_dl(smica.squeeze(), nside=nsidesmica)\n",
    "        dl_SMICA_HU /= np.max(dl_SMICA_HU[0:len(ell) // 2])\n",
    "        dl_SMICA /= np.max(dl_SMICA[0:len(ell1) // 2])\n",
    "        fig = plt.figure(figsize=[12, 12])\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.set_xlabel(\"$ell(ell+1)C_ell/2\\pi \\, \\,(\\mu K^2)$\")\n",
    "        line, = ax.plot(ell1, dl_SMICA)\n",
    "        line.set_label('SMICA')\n",
    "        ax1 = plt.twinx(ax)\n",
    "        line1, = ax1.plot(ell, dl_SMICA_HU)\n",
    "        line1.set_label('HU')\n",
    "        ax1.set_ylabel('$\\ell$')\n",
    "        ax1.set_title(\"Angular Power Spectra {}_{}\".format(self.kmax, self.nside3D))\n",
    "        ax1.legend(loc=\"upper right\")\n",
    "        ax1.set_xlim(1, xmax)\n",
    "        plt.ylim(1E-10, 1)\n",
    "        ax1.grid()\n",
    "        plt.savefig(\"./img1/AngularPowerSpectra_{}_{}.png\".format(kmax, nside3D), dpi=300)\n",
    "        plt.show()\n",
    "        return dl_SMICA, dl_SMICA_HU, ell\n",
    "\n",
    "    def plotSH(self, l, m, ll, lm, pp):\n",
    "        fcolors = self.spherharmm(l, m, self.phi, pp)\n",
    "        hp.mollview(fcolors.squeeze(), title=\"{}_{}_{}_{}\".format(l, m, chg2ang(ll), chg2ang(lm)), cmap=cm.RdBu_r)\n",
    "        hp.graticule()\n",
    "        plt.show()\n",
    "        return fcolors\n",
    "\n",
    "    def plotNewMap(self, newmap, err, filename=None, title=None, plotme=False, save=False, nosigma=True):\n",
    "        err0 = np.round(err, 3)\n",
    "        if filename is None:\n",
    "            filename = \"./img1/aitoff_{}_{}_{}__{}_{}_{}.png\".format(self.kmax, self.nside3D, err0,\n",
    "                                                                     chg2ang(self.lambda_k),\n",
    "                                                                     chg2ang(self.lambda_l),\n",
    "                                                                     chg2ang(self.lambda_m))\n",
    "        if title is None:\n",
    "            title = \"{}_{}_{}__{}_{}_{}\".format(self.kmax, self.nside3D, err0, chg2ang(self.lambda_k),\n",
    "                                                chg2ang(self.lambda_l),\n",
    "                                                chg2ang(self.lambda_m))\n",
    "        self.plot_aitoff(newmap, self.lambda_k, self.lambda_l,\n",
    "                         self.lambda_m, err0, title=title, filename=filename, plotme=plotme, save=save, nosigma=nosigma)\n",
    "\n",
    "    def calcStuff(self, alpha):\n",
    "        err = 0.0\n",
    "        results = 0.0\n",
    "        prime = alpha[0]\n",
    "        kmax = alpha[1]\n",
    "        lk = prime[0]\n",
    "        ll = prime[1]\n",
    "        lm = prime[2]\n",
    "        start_time = time()\n",
    "        try:\n",
    "            results, fcolors, ls, err = self.project4D3d_0(kmax)\n",
    "        except Exception as exception_object:\n",
    "            print('Error with getting map for: {0}_{0}_{0}'.format(lk, ll, lm))\n",
    "        stop_time = time()\n",
    "        print(alpha, err)\n",
    "        return [start_time, stop_time, err, prime, results]\n",
    "\n",
    "    def mindict(self, test_dict):\n",
    "        # Using min() + list comprehension + values()\n",
    "        # Finding min value keys in dictionary\n",
    "        temp = min(test_dict.values())\n",
    "        res = [key for key in test_dict if test_dict[key] == temp]\n",
    "        return res[0], test_dict[res[0]]\n",
    "\n",
    "    def calcError(self, x, karray):\n",
    "        t1 = datetime.now()\n",
    "        self.change_HSH_center(x[0], x[1], x[2], karray, nside3D, loadpriorG=False, doit=True, savePG=False)\n",
    "        _, _, err = self.project4D3d_0(karray)\n",
    "        print(x, err, (datetime.now() - t1).microseconds)\n",
    "        return err\n",
    "\n",
    "    def calcErrorDF(self, x):\n",
    "        err = np.sum((np.dot(x.T, self.df[:, 4:]) - self.SMICA_LR.squeeze()) ** 2) * 1E6\n",
    "        return err\n",
    "\n",
    "    def cleanup(self, mypath, prefix):\n",
    "        filelist = [join(mypath, f) for f in sorted(listdir(mypath)) if\n",
    "                    isfile(join(mypath, f)) and f.startswith(prefix)]\n",
    "        for f in filelist:\n",
    "            if os.path.exists(f):\n",
    "                os.remove(f)\n",
    "            else:\n",
    "                print(\"The file does not exist\")\n",
    "\n",
    "    def get_all_gegenbauer(self):\n",
    "        if not self.loadpriorG:\n",
    "            self.G = pd.DataFrame()\n",
    "        G=[]\n",
    "        for k in np.arange(2,np.max(self.karray)+1):\n",
    "                for l in np.arange(1,k):\n",
    "                    if (1+l,k-l) not in self.G.index:\n",
    "                        G.append( (1+l,k-l) )\n",
    "        if len(G)!=0:\n",
    "            self.G = get_gegenbauerAsyn(self.cosksi, G)\n",
    "        if self.savePG:\n",
    "            self.G.to_pickle(self.G_filename)\n",
    "\n",
    "    def change_HSH_center(self, lk, ll, lm, karray, nside, doit=False, loadpriorG=False, savePG=False):\n",
    "        kmax = max(karray)\n",
    "        if nside != self.nside3D:\n",
    "            self.xx, self.yy, self.zz = hp.pix2vec(nside=nside, ipix=np.arange(hp.nside2npix(nside)))\n",
    "            self.nside3D = nside\n",
    "        if not doit:\n",
    "            array_number_match = len(self.karray) == len(karray)\n",
    "            if array_number_match:\n",
    "                array_match = np.sum(self.karray == karray) != 0\n",
    "            else:\n",
    "                array_match = False\n",
    "            if (self.lambda_k, self.lambda_l, self.lambda_m, self.kmax) == (lk, ll, lm, kmax) \\\n",
    "                    and array_number_match and array_match:\n",
    "                return\n",
    "        self.loadpriorG = loadpriorG\n",
    "        if not loadpriorG:\n",
    "            self.G = pd.DataFrame()\n",
    "        self.get_df_size(karray, aggregLevel=\"M\")\n",
    "        self.lambda_k = lk\n",
    "        self.lambda_l = ll\n",
    "        self.lambda_m = lm\n",
    "        self.G_filename = os.path.join(self.mypath, \"G_{}_{}_{}_{}.npy\".format(self.nside3D,\n",
    "                                                                                  chg2ang(lk),\n",
    "                                                                                  chg2ang(ll),\n",
    "                                                                                  chg2ang(lm)))\n",
    "        self.z = self.zz + lk\n",
    "        self.y = self.yy + ll\n",
    "        self.x = self.xx + lm\n",
    "        self.cosksi = np.cos(self.z)\n",
    "        self.sinksi = np.sin(self.z)\n",
    "        self.costheta = np.cos(self.y)\n",
    "        self.phi = self.x\n",
    "        self.get_all_gegenbauer()\n",
    "        kjump = 20\n",
    "        ibeginning=0\n",
    "        for i in np.arange(kjump, self.kmax, kjump):\n",
    "            newPfile = self.P_filename + \"{}.npy\".format(i)\n",
    "            if not os.path.exists(newPfile):\n",
    "                pp = get_legendre((0,np.arange(ibeginning,i), self.costheta))[1].squeeze()\n",
    "                np.save(newPfile,pp)\n",
    "                ibeginning=i\n",
    "\n",
    "\n",
    "    def createGaussianBackground(self, x0, nside, delta, karray):\n",
    "        xx, yy, zz = hp.pix2vec(nside=nside, ipix=np.arange(hp.nside2npix(nside)))\n",
    "        x, y, z = tuple(map(operator.add, (xx, yy, zz), x0))\n",
    "        theta = y\n",
    "        phi = x\n",
    "        ############################################################\n",
    "        ############################################################\n",
    "        # random amplitudes\n",
    "        wavefunction = np.zeros(xx.shape)\n",
    "        costheta = np.cos(theta)\n",
    "        t1 = datetime.now()\n",
    "        pp = self.get_P(karray, costheta)\n",
    "        print((datetime.now() - t1).seconds)\n",
    "        for ii, l in enumerate(karray):\n",
    "            mu = 1 / np.sqrt(2 * l + 1)\n",
    "            sigma = mu / np.sqrt(2 * l + 1)\n",
    "            amp = np.abs(np.random.normal(mu, sigma, [1, 2 * l + 1]))\n",
    "            mlist = list(set([int(np.round(kk)) for kk in np.linspace(-l, l, 20)]))\n",
    "            for m in mlist:\n",
    "                if m > 0:\n",
    "                    wavefunction += amp[0, l + m] * self.spherharmm(ii, l, m, phi, pp) * np.sin(l * delta)\n",
    "                if m == 0:\n",
    "                    wavefunction += amp[0, l + m] * self.spherharmm(ii, l, m, phi, pp)\n",
    "                if m < 0:\n",
    "                    wavefunction += amp[0, l + m] * self.spherharmm(ii, l, m, phi, pp) * np.cos(l * delta)\n",
    "        ############################################################\n",
    "        return wavefunction\n",
    "\n",
    "    def project4D3d(self, karray):\n",
    "        self.kmax = np.max(karray)\n",
    "        self.calc_hyperharmnano2(karray)\n",
    "        SMICA = self.SMICA_LR\n",
    "        C = np.dot(self.df[:, 4:], SMICA)\n",
    "        B = np.dot(self.df[:, 4:], self.df[:, 4:].T)\n",
    "        results = np.linalg.solve(B, C)\n",
    "        self.newmap = np.dot(results.T, self.df[:, 4:])\n",
    "        sigma = np.std(self.newmap)\n",
    "        self.newmap = self.newmap / sigma * self.sigma_smica\n",
    "        err = (1.0 - np.correlate(self.newmap.squeeze(), self.SMICA_LR.squeeze()) * 1e4)[0]\n",
    "        return results, self.newmap, err\n",
    "\n",
    "    def eval_my_gegenbauer(self, k, l):\n",
    "        A = np.zeros(self.cosksi.shape)\n",
    "        for i, x in enumerate(self.cosksi):\n",
    "            A[i] = mp.gegenbauer(1 + l, k - l, x)\n",
    "            # A[i] = sp.eval_gegenbauer(1 + l, k - l, x)\n",
    "        return A\n",
    "\n",
    "    def get_gegenbauer(self):\n",
    "        multiindex = pd.MultiIndex.from_tuples(list(self.extra_G_keys()), names=('k', 'l'))\n",
    "        columns = list(np.arange(len(self.xx)))\n",
    "        extra_G_keys = [x for x in self.extra_G_keys()]\n",
    "        if len(extra_G_keys) == 0:\n",
    "            return\n",
    "        data = np.zeros([len(multiindex), len(columns)])\n",
    "        extraG = pd.DataFrame(data=data, index=multiindex, columns=columns, dtype=float)\n",
    "        for aa, bb in self.extra_G_keys():\n",
    "            l = aa - 1\n",
    "            k = bb + l\n",
    "            a = eval_my_gegenbauer(aa, bb, self.cosksi)\n",
    "            if len(np.argwhere(np.isnan(a))) == 0:\n",
    "                extraG.loc[(aa, bb), :] = a\n",
    "            else:\n",
    "                print(\"error\", aa, bb, k, l)\n",
    "        self.G = pd.concat([self.G, extraG], axis=0)\n",
    "\n",
    "    def get_P(self, l, costheta):\n",
    "        pp = np.zeros([0, kmax + 1, kmax + 1])\n",
    "        lencos = len(costheta)\n",
    "        number_of_workers = 10\n",
    "        mychunks = chunks(costheta, lencos // number_of_workers, l)\n",
    "        P = {}\n",
    "        with ProcessPoolExecutor(max_workers=number_of_workers) as executor:\n",
    "            responses = executor.map(get_legendre, mychunks)\n",
    "        for response in responses:\n",
    "            i = response[0]\n",
    "            P[i] = response[1]\n",
    "        for key in np.sort(list(P.keys())):\n",
    "            pp = np.concatenate([pp, P[key]], axis=0)\n",
    "        return pp\n",
    "\n",
    "    # def calc_hyperharmnano2(self, karray):\n",
    "    #     kmax = np.max(karray)\n",
    "    #     files = sorted( [x for x in self.getSpectralFiles(self.mypath) if \"/PG_data/P_\" in x] )\n",
    "    #     karray = [eval(x.replace(self.mypath,\"\").replace(\".npy\",\"\").split(\"_\")[-1:][0]) for x in files]\n",
    "    #     filesMap = {x:y for (x,y) in zip(karray, files)}\n",
    "    #     jj = 0\n",
    "    #     G = {}\n",
    "    #     listofk = [x for x in sorted(list(filesMap.keys())) if x <=kmax]\n",
    "    #     for k in listofk:\n",
    "    #         pp = np.load(filesMap[k])\n",
    "    #         # if self.longG:\n",
    "    #         #     b = np.zeros(self.xx.shape)\n",
    "    #         llist = [np.count_nonzero(x, axis=0) for x in pp[0,:,:].squeeze()]\n",
    "    #         llist = [x-1 for x in llist if x<k+1]\n",
    "    #         for ii, l in enumerate(llist):\n",
    "    #             if self.longG:\n",
    "    #                 b = np.zeros(self.xx.shape)\n",
    "    #             if (1+l,k-l) not in self.G.index:\n",
    "    #                 G[1+l,k-l]=1.0\n",
    "    #                 continue\n",
    "    #             a = self.sinksi ** l * self.G.loc[(1 + l, k - l), :].values\n",
    "    #             a1 = (-1) ** k * np.sqrt(2 * (k + 1) / np.pi * mp.factorial(k - l)\n",
    "    #                                      * mp.power(2,2 * l) * mp.factorial(l) ** 2 / mp.factorial( k + l + 1))\n",
    "    #             a1 = float(a1)\n",
    "    #             for m in np.arange(-l,l+1):\n",
    "    #                 if not self.longG:\n",
    "    #                     b = a * a1 * self.spherharmm(ii, l, m, self.phi, pp)\n",
    "    #                     c = b.std()\n",
    "    #                     if c != 0:\n",
    "    #                         self.df[jj, 0] = k\n",
    "    #                         self.df[jj, 1] = l\n",
    "    #                         self.df[jj, 2] = m\n",
    "    #                         self.df[jj, 3] = c\n",
    "    #                         self.df[jj, 4::] = b / c\n",
    "    #                         jj += 1\n",
    "    #                     else:\n",
    "    #                         print(k, l, m)\n",
    "    #                 else:\n",
    "    #                     b += a * a1 * self.spherharmm(ii, l, m, self.phi, pp)\n",
    "    #                     print(k,l,m)\n",
    "    #                     if self.longG:\n",
    "    #                         c = b.std()\n",
    "    #                         if c != 0:\n",
    "    #                             self.df[jj, 0] = k\n",
    "    #                             self.df[jj, 1] = 0\n",
    "    #                             self.df[jj, 2] = 0\n",
    "    #                             self.df[jj, 3] = c\n",
    "    #                             self.df[jj, 4::] = b / c\n",
    "    #                             jj += 1\n",
    "    #                         else:\n",
    "    #                             print(k, \"wavefunction failed\")\n",
    "    #     if len(G.keys())!=0:\n",
    "    #         self.extra_G_keys=G.keys\n",
    "    #         self.get_gegenbauer()\n",
    "    #         print(\"DO IT AGAIN\")\n",
    "    #     print(jj)\n",
    "\n",
    "    def calc_hyperharmnano2(self, karray):\n",
    "        kmax = np.max(karray)\n",
    "        npoints = len(self.xx)\n",
    "        files = sorted( [x for x in self.getSpectralFiles(self.mypath) if \"/PG_data/P_\" in x] )\n",
    "        karray0 = [eval(x.replace(self.mypath,\"\").replace(\".npy\",\"\").split(\"_\")[-1:][0]) for x in files]\n",
    "        filesMap = {x:y for (x,y) in zip(karray0, files)}\n",
    "        G = {}\n",
    "        listofk = [x for x in sorted(list(filesMap.keys())) if x <=kmax]\n",
    "        calcpp = len(listofk)==0 or kmax < 50\n",
    "        if calcpp:\n",
    "            listofk= sorted(karray)\n",
    "            lmax= np.max(karray)\n",
    "            pp = np.zeros([len(self.costheta), lmax+1, lmax+1] )\n",
    "            for i, x in enumerate(self.costheta):\n",
    "                pp[i, :, :] = legendre(lmax, x, csphase=-1)\n",
    "            pp=pp[:,1:,:]\n",
    "        ##############################################################\n",
    "        jj = 0\n",
    "        for k in listofk:\n",
    "            llist = np.arange(1,k)\n",
    "            for ii, l in enumerate(llist):\n",
    "                for m in np.arange(-l, l + 1):\n",
    "                    jj+=1\n",
    "        self.df = np.zeros([jj,4+ npoints])\n",
    "        ##############################################################\n",
    "        jj = 0\n",
    "        for k in listofk:\n",
    "            if not calcpp:\n",
    "                pp = np.load(filesMap[k])\n",
    "                llist = [np.count_nonzero(x, axis=0) for x in pp[0,:,:].squeeze()]\n",
    "                llist = [x-1 for x in llist if x<k+1 and x>1]\n",
    "            else:\n",
    "                llist = np.arange(1,k)\n",
    "            for ii, l in enumerate(llist):\n",
    "                if self.longG:\n",
    "                    b = np.zeros(self.xx.shape)\n",
    "                if (1+l,k-l) not in self.G.index:\n",
    "                    G[1+l,k-l]=1.0\n",
    "                    print(k,l)\n",
    "                    continue\n",
    "                a = self.sinksi ** l * self.G.loc[(1 + l, k - l), :].values\n",
    "                a1 = (-1) ** k * np.sqrt(2 * (k + 1) / np.pi * mp.factorial(k - l)\n",
    "                                         * mp.power(2,2 * l) * mp.factorial(l) ** 2 / mp.factorial( k + l + 1))\n",
    "                a1 = float(a1)\n",
    "                for m in np.arange(-l,l+1):\n",
    "                    if self.longG:\n",
    "                        b += a * a1 * self.spherharmm(ii, l, m, self.phi, pp)\n",
    "                    else:\n",
    "                        if calcpp:\n",
    "                            b = a * a1 * self.spherharmm(l-1, l, m, self.phi, pp)\n",
    "                        else:\n",
    "                            b = a * a1 * self.spherharmm(ii, l, m, self.phi, pp)\n",
    "                        c = b.std()\n",
    "                        self.df[jj, 0] = k\n",
    "                        self.df[jj, 1] = l\n",
    "                        self.df[jj, 2] = m\n",
    "                        self.df[jj, 3] = c\n",
    "                        self.df[jj, 4::] = b / c\n",
    "                        jj += 1\n",
    "                if self.longG:\n",
    "                    c = b.std()\n",
    "                    self.df[jj, 0] = k\n",
    "                    self.df[jj, 1] = l\n",
    "                    self.df[jj, 2] = m\n",
    "                    self.df[jj, 3] = c\n",
    "                    self.df[jj, 4::] = b / c\n",
    "                    jj += 1\n",
    "        if len(G.keys())!=0:\n",
    "            self.extra_G_keys=G.keys\n",
    "            print(G.keys())\n",
    "            get_gegenbauerAsyn(self.cosksi, G.keys())\n",
    "            print(\"DO IT AGAIN\")\n",
    "        print(jj)\n",
    "\n",
    "    def spherharmm(self, ii, l, m, phi, pp):\n",
    "        mm = np.abs(m)\n",
    "        # a = float64((-1) ** mm * np.sqrt((2 * l + 1) * mp.fac(l - mm) / mp.fac(l + mm)))\n",
    "        a = float64(np.sqrt((2 * l + 1) * mp.fac(l - mm) / mp.fac(l + mm)))\n",
    "        if m > 0:\n",
    "            return a * pp[:, ii, mm] * np.cos(mm * phi)\n",
    "        if m == 0:\n",
    "            return a * pp[:, ii, mm]\n",
    "        if m < 0:\n",
    "            return a * pp[:, ii, mm] * np.sin(mm * phi)\n",
    "\n",
    "\n",
    "    def get_df_size(self, karray, aggregLevel=\"M\"):\n",
    "        npoints = len(self.xx)\n",
    "        kmax = np.max(karray)\n",
    "        jj = 0\n",
    "        G = {}\n",
    "        files = sorted([x for x in self.getSpectralFiles(self.mypath) if \"/PG_data/P_{}_{}_{}_{}\".format(self.nside3D,\n",
    "                                                                                  chg2ang(self.lambda_k),\n",
    "                                                                                  chg2ang(self.lambda_l),\n",
    "                                                                                  chg2ang(self.lambda_m)) in x])\n",
    "        karray0 = [eval(x.replace(self.mypath, \"\").replace(\".npy\", \"\").split(\"_\")[-1:][0]) for x in files]\n",
    "        filesMap = {x: y for (x, y) in zip(karray0, files)}\n",
    "        listofk = [x for x in sorted(list(filesMap.keys())) if x <= kmax]\n",
    "        for k in listofk:\n",
    "                pp = np.load(filesMap[k])\n",
    "                llist = [np.count_nonzero(x, axis=0) for x in pp[0,:,:].squeeze()]\n",
    "                pp = None\n",
    "                llist = [x-1 for x in llist if x<k+1 and x>1]\n",
    "                for ii, l in enumerate(llist):\n",
    "                    if (1+l,k-l) not in self.G.index:\n",
    "                        G[1+l,k-l]=1.0\n",
    "                    for m in np.arange(-l,l+1):\n",
    "                        if not self.longG:\n",
    "                            jj+=1\n",
    "                    if self.longG:\n",
    "                        jj += 1\n",
    "        self.df = np.zeros([jj,4+npoints])\n",
    "        self.extra_G_keys = G.keys\n",
    "\n",
    "\n",
    "    def plot_orth(self, fcolors):\n",
    "        sigma = np.std(fcolors)\n",
    "        hp.orthview(fcolors, min=-2 * sigma, max=2 * sigma, title='Raw WMAP data', unit=r'$\\Delta$T (mK)')\n",
    "        plt.show()\n",
    "\n",
    "    def plot_aitoff_df(self, l, m, phi, pp=None, cmap=cm.RdBu_r):\n",
    "        if pp is None:\n",
    "            pp = np.array([legendre(l, x) for x in self.costheta])\n",
    "        fcolors = self.spherharmm(l=l, m=m, phi=phi, pp=pp)\n",
    "        sigma = np.std(fcolors)\n",
    "        title = \"{}_{}\".format(l, m)\n",
    "        hp.mollview(fcolors.squeeze(), title=title, min=-2 * sigma,\n",
    "                    max=2 * sigma, unit=\"K\", cmap=cm.RdBu_r)\n",
    "        hp.graticule()\n",
    "        self.plot_orth(fcolors)\n",
    "\n",
    "    def newerr(self, x):\n",
    "        newnewmap = myHyper.SMICA_LR - x[1] * self.newmap + x[0]\n",
    "        err = np.std(newnewmap)\n",
    "        return err\n",
    "\n",
    "    def optimizeNewMap(self, newmap0, SMICA_LR, xx0, nside3D, xmax=300, nosigma=True, mypath=\"./PG_data/\"):\n",
    "        newmap = xx0[1] * newmap0 - xx0[0]\n",
    "        diffmap = SMICA_LR.squeeze() - newmap\n",
    "        self.plotNewMap(newmap, err=0.0, filename=\"newmap_{}_{}.png\".format(kmax, nside3D), title=\"Optimized  CMB Map\", plotme=True, save=False, nosigma=nosigma)\n",
    "        self.plotNewMap(diffmap, err=0.0, filename=\"diffmap_{}_{}.png\".format(kmax, nside3D), title=\"Optimized Diff Map\", plotme=True, save=False, nosigma=nosigma)\n",
    "        return newmap, diffmap\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lib4 import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myHyper = None\n",
    "# History\n",
    "#################################################################\n",
    "#################################################################\n",
    "#################################################################\n",
    "#################################################################\n",
    "# Position of our Universe (error = -1.3924651085864181 for base 2:11)\n",
    "x0 = [4.92463563, 4.95336686, 0.8165214728204138] #[339.46, 341.1, 104.08]\n",
    "y0 = [np.round((xx+1) / np.pi * 180, 2) for xx in x0]\n",
    "print(y0)\n",
    "# In degrees  [238.34, 63.23, 113.78]\n",
    "# In degrees [120.0, 300.0, 240.0]\n",
    "todo = [Color.EVALUATE_DF_AT_POSITION, Color.OPTIMIZE_SPECTRUM ]\n",
    "# Find Neighborhood\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.EVALUATE_DF_AT_POSITION in todo:\n",
    "    # x0 = [2.11999902, 4.84227338, 4.95303838] # -1.4133892288576106 Best_48_48__121_277_283 June02\n",
    "    # x0 = [4.14376593, 1.09484514, 1.89120155]\n",
    "    # x0 = [4.92463563, 4.95336686, 0.8165214728204138]\n",
    "#     x0 = [4.82349567, 4.88348719, 0.84566045]  # -3.442960596169308 9876\n",
    "    # x0 = [2.13302707, 5.08554078, 4.34175881]  # -0.958704 [122.21, 291.38, 248.76]\n",
    "    # 2.13302707, 5.08554078, 4.34175881 -2.037992624771841 Best_29_64__122_291_248\n",
    "    (lk, ll, lm) = x0\n",
    "    print(\"evaluated at {}\".format(x0))\n",
    "    # Create karray\n",
    "    # karray = list(np.arange(2, 49))\n",
    "    # karray = list(np.arange(2, 10))\n",
    "\n",
    "    # karray = list(np.arange(20, 80, 20))\n",
    "    karray = list(np.arange(2,49))\n",
    "    karray = np.array(sorted(list(set(karray)))).astype(int)\n",
    "\n",
    "    nside3D = 48\n",
    "    bandwidth = 48\n",
    "    # kmax=200\n",
    "    # karray += sorted(list(set([int(k) for k in np.geomspace(49,kmax,20)])))\n",
    "    karray = np.array(sorted(list(set(karray)))).astype(int)\n",
    "    print(len(karray))\n",
    "    kmax = max(karray)\n",
    "    #################################\n",
    "    if myHyper is None:\n",
    "        mypath= \"./PG_data\"\n",
    "        myHyper = HYPER(nside3D, sigma_smica, planck_IQU_SMICA, karray, mypath,\n",
    "                        lk, ll, lm, loadpriorG=True, savePG=True, bandwidth=bandwidth, longG=False)\n",
    "    #################################################################\n",
    "    results, newmap0, err = myHyper.project4D3d_0(karray)\n",
    "    filename = \"./img1/SingleBest_{}_{}__{}_{}_{}.png\".format(kmax, nside3D, chg2ang(lk),\n",
    "                                                              chg2ang(ll), chg2ang(lm))\n",
    "    title = \"Best_{}_{}__{}_{}_{}\".format(kmax, nside3D, chg2ang(lk),\n",
    "                                          chg2ang(ll), chg2ang(lm))\n",
    "\n",
    "    myHyper.plotNewMap(newmap0, err, filename=filename, title=title, plotme=True, save=True)\n",
    "    myHyper.plotHistogram(newmap0.squeeze(), nside3D, kmax, plotme=True)\n",
    "    print(lk, ll, lm, err, title)\n",
    "    df = pd.DataFrame(np.concatenate([results, myHyper.df[:, 0:4]], axis=1),\n",
    "                      columns=[\"coeff\", \"k\", \"l\", \"m\", \"std\"])\n",
    "    filename = \"./PG_data/spectrum_{}_{}__{}_{}_{}\".format(kmax, nside3D, chg2ang(lk), chg2ang(ll), chg2ang(lm))\n",
    "    df.to_pickle(filename)\n",
    "    #################################################################\n",
    "    np.save(\"./PG_data/resultsBest_{}_{}.npy\".format(kmax, nside3D), results, allow_pickle=True)\n",
    "    np.save(\"./img1/SMICA_LR.npy\", myHyper.SMICA_LR)\n",
    "    #################################################################\n",
    "    #################################################################\n",
    "    x01 = np.array([0.00140589, 0.43367364])\n",
    "    x00 = minimize(newerr, x01, args=(myHyper.SMICA_LR, newmap0),\n",
    "                   method='nelder-mead', options={'xatol': 1e-8, 'disp': True})\n",
    "    err = x00.fun\n",
    "    xx0 = x00.x\n",
    "    newmap, diffmap = myHyper.optimizeNewMap(newmap0, myHyper.SMICA_LR,xx0,\n",
    "                                             nside3D=nside3D, bandwidth=bandwidth, nosigma=False)\n",
    "    \n",
    "    #################################################################\n",
    "    #################################################################\n",
    "    #################################################################\n",
    "    np.save(\"./img1/newmap.npy\", newmap)\n",
    "    np.save(\"./img1/diffmap.npy\", diffmap)\n",
    "    ###############################################################\n",
    "    ###############################################################\n",
    "#     myHyper.plot_Single_CL_From_Image(diffmap, nside3D, xmax=xmax)\n",
    "#     myHyper.plot_Single_CL_From_Image(newmap, nside3D, xmax=xmax)\n",
    "    ###############################################################\n",
    "    ###############################################################\n",
    "    print(lk, ll, lm, err, title)\n",
    "    df = pd.DataFrame(np.concatenate([results, myHyper.df[:, 0:4]], axis=1),\n",
    "                      columns=[\"coeff\", \"k\", \"l\", \"m\", \"std\"])\n",
    "    filename = \"./PG_data/spectrum_{}_{}__{}_{}_{}\".format(kmax, nside3D, chg2ang(lk), chg2ang(ll), chg2ang(lm))\n",
    "    df.to_pickle(filename)\n",
    "    #################################################################\n",
    "    np.save(\"./PG_data/resultsBest_{}_{}.png\".format(kmax, nside3D), results, allow_pickle=True)\n",
    "    np.save(\"./img1/SMICA_LR.png\", myHyper.SMICA_LR)\n",
    "    #################################################################\n",
    "    #################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myHyper.plotNewMap(diffmap, err, filename=filename, title=title, plotme=True, save=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.FINDNEIGHBORHOOD in todo:\n",
    "    # these three indices are related to position within the hyperspherical hypersurface.  Don't confuse them with the\n",
    "    # quantum numbers k,l,m\n",
    "    x0 = [0.0, 0.0, 0.0]\n",
    "    (lambda_k, lambda_l, lambda_m) = x0\n",
    "    # Create karray\n",
    "    nside3D = 48\n",
    "    bandwidth = 64\n",
    "    karray = list(np.arange(2, 10))\n",
    "    print(len(karray))\n",
    "    kmax = max(karray)\n",
    "    #################################\n",
    "    if myHyper is None:\n",
    "        myHyper = HYPER(nside3D, sigma_smica, planck_IQU_SMICA, karray,\n",
    "                        lambda_k, lambda_l, lambda_m, loadpriorG=False, savePG=False, bandwidth=bandwidth)\n",
    "    else:\n",
    "        myHyper.savePG = False\n",
    "        myHyper.loadpriorG = False\n",
    "    #################################################################\n",
    "    # these three indices are related to position within the hyperspherical hypersurface.  Don't confuse them with the\n",
    "    # quantum numbers k,l,m\n",
    "    errarray = []\n",
    "    olderr = 1110.0\n",
    "    results = 0.0\n",
    "    x00 = []\n",
    "    nside3D = 48\n",
    "    bandwidth = 48\n",
    "    n = 7\n",
    "    myHyper.change_SMICA_resolution(nside3D, doit=False, bandwidth=bandwidth)\n",
    "    for lk in np.linspace(0, 2 * np.pi, n):\n",
    "        for ll in np.linspace(0, 2 * np.pi, n):\n",
    "            for lm in np.linspace(0, 2 * np.pi, n):\n",
    "                start_time = time()\n",
    "                try:\n",
    "                    myHyper.change_HSH_center(lk, ll, lm, karray, nside3D, loadpriorG=False, doit=True,\n",
    "                                              savePG=False)\n",
    "                    results, fcolors, err = myHyper.project4D3d_0(karray)\n",
    "                    myHyper.plotNewMap(fcolors, err, filename=None, title=None, plotme=False, save=True)\n",
    "                    if olderr > err:\n",
    "                        olderr = err\n",
    "                        err0 = np.round(err, 3)\n",
    "                        filename = \"./img1/Bestf_{}_{}_{}__{}_{}_{}.png\".format(myHyper.kmax, myHyper.nside3D,\n",
    "                                                                                err0,\n",
    "                                                                                chg2ang(lk),\n",
    "                                                                                chg2ang(ll),\n",
    "                                                                                chg2ang(lm))\n",
    "                        myHyper.plotNewMap(fcolors, err, filename=filename, title=None, plotme=True, save=True)\n",
    "                        x00.append((lk, ll, lm, err))\n",
    "                        np.save(\"./img1/x0_{}_{}.npy\".format(kmax, nside3D), x00[-1:], allow_pickle=True)\n",
    "                        print(lk, ll, lm, err)\n",
    "                except Exception as aa:\n",
    "                    print('Error with getting map for: {}_{}_{}'.format(lk, ll, lm))\n",
    "                stop_time = time()\n",
    "    myHyper.creategiff(kmax, nside3D, mypath=\"./img1\", prefix=\"aitoff_\", filename=\"CMB\")\n",
    "    sleep(10)\n",
    "    myHyper.creategiff(kmax, nside3D, mypath=\"./img1\", prefix=\"Bestf\", filename=\"BEST\")\n",
    "    sleep(10)\n",
    "    myHyper.cleanup(mypath=\"./img1\", prefix=\"aitoff_\")\n",
    "    # myHyper.cleanup(mypath=\"./img1\", prefix=\"Bestf\")\n",
    "    print(\"Best Position = \", x00[-1:])\n",
    "    bestposition = \"./img1/x0_{}_{}.npy\".format(kmax, nside3D)\n",
    "    np.save(\"./img1/x0_{}_{}.npy\".format(kmax, nside3D), x00[-1:], allow_pickle=True)\n",
    "    np.save(\"./img1/x00_{}_{}.npy\".format(kmax, nside3D), x00, allow_pickle=True)\n",
    "#################################################################\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.MINIMIZEPOSITION in todo:\n",
    "    # Best Position = [(2.0943951023931953, 4.886921905584122, 4.886921905584122, -0.9543073989083581)]\n",
    "    # x0= x00[-1:][0][0:3]\n",
    "    # x0 = (2.0943951023931953, 4.886921905584122, 4.886921905584122)\n",
    "    # x0 = [2.11999902, 4.84227338, 4.95303838]  # ksi, theta, phi = [121.47, 277.44, 283.79]\n",
    "    # x0 = [2.0943951023931953, 5.235987755982988, 4.1887902047863905] # -0.9512365749096219 June 09\n",
    "    # x0 = [2.13302707, 5.08554078, 4.34175881]  # -0.958704 [122.21, 291.38, 241.71] June9 Minimized\n",
    "    # x0 = [4.92463563, 4.95336686, 0.8165214728204138]\n",
    "    # x0 = [4.82349567, 4.88348719, 0.84566045] #- 3.3740886193323982  696721\n",
    "    x0 = [4.82349567, 4.88348719, 0.84566045] #-3.442960596169308 9876\n",
    "    nside3D = 64\n",
    "    bandwidth = 64\n",
    "    (lambda_k, lambda_l, lambda_m) = x0\n",
    "    # Create karray\n",
    "    karray = list(np.arange(2, 30))\n",
    "    print(len(karray))\n",
    "    kmax = max(karray)\n",
    "    #################################\n",
    "    if myHyper is None:\n",
    "        myHyper = HYPER(nside3D, sigma_smica, planck_IQU_SMICA, karray,\n",
    "                        lambda_k, lambda_l, lambda_m, loadpriorG=False, savePG=False, bandwidth=bandwidth)\n",
    "    else:\n",
    "        myHyper.loadpriorG = False\n",
    "        myHyper.savePG = False\n",
    "    #################################\n",
    "    print(x0)\n",
    "    x0 = minimize(myHyper.calcError, x0, args=(karray), method='nelder-mead',\n",
    "                  options={'xatol': 1e-4, 'disp': True})\n",
    "    np.save(\"./img1/x0_{}_{}.npy\".format(kmax, nside3D), x0.x, allow_pickle=True)\n",
    "    x0 = x0.x\n",
    "    (lk, ll, lm) = x0\n",
    "    print(\"minimized at {}\".format(x0))\n",
    "    filename = \"./img1/Best_{}_{}__{}_{}_{}.png\".format(kmax, nside3D, chg2ang(lk), chg2ang(ll), chg2ang(lm))\n",
    "    title = \"Best_{}_{}__{}_{}_{}\".format(kmax, nside3D, chg2ang(lk), chg2ang(ll), chg2ang(lm))\n",
    "    myHyper.change_HSH_center(lk, ll, lm, karray, nside3D, loadpriorG=False, doit=True, savePG=True)\n",
    "    results, newmap, err = myHyper.project4D3d_0(karray)\n",
    "    np.save(\"./img1/results_{}_{}.npy\".format(kmax, nside3D), results, allow_pickle=True)\n",
    "#################################################################\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.OPTIMIZE_SPECTRUM in todo:\n",
    "    (lk, ll, lm) = x0\n",
    "    # Create karray\n",
    "    karray = list(np.arange(20, 222, 20))\n",
    "    karray = np.array(sorted(list(set(karray)))).astype(int)\n",
    "    print(len(karray))\n",
    "    kmax = max(karray)\n",
    "    nside3D = 128\n",
    "    bandwidth = 1024\n",
    "    myHyper = None\n",
    "    myHyper = HYPER(nside3D, sigma_smica, planck_IQU_SMICA, karray,\n",
    "                    lk, ll, lm, loadpriorG=True, savePG=True, bandwidth=bandwidth, longG=True)\n",
    "    myHyper.SMICA_LR = np.load(\"./img1/diffmap.npy\").T\n",
    "    myHyper.plotNewMap(myHyper.SMICA_LR, 0, filename=None, title=\"DiffMap\", plotme=True, save=False, nosigma=False)\n",
    "    results, newmap, err, packedresults = myHyper.project4D3d_0(karray)\n",
    "    # #################################\n",
    "    myHyper.plot_Single_CL_From_Image(newmap, nside3D, xmax=50,log=True)\n",
    "    myHyper.plot_Single_CL_From_Image(diffmap, nside3D, xmax=300)\n",
    "    myHyper.plot_CL_From_Image(newmap, nside3D, planck_theory_cl, xmax=300, ymax=0)\n",
    "\n",
    "    newnewmap = myHyper.SMICA_LR.squeeze() - xx0[1] * newmap.squeeze() + xx0[0]\n",
    "\n",
    "    myHyper.plot_CL_From_Image(diffmap, nside3D, planck_theory_cl, xmax=300, ymax=0.1)\n",
    "\n",
    "    myHyper.plot_CL_From_Image(myHyper.SMICA_LR.squeeze(), nside3D, planck_theory_cl, xmax=3000, ymax=1)\n",
    "    df = pd.DataFrame(np.concatenate([results, myHyper.df[:, 0:4]], axis=1),\n",
    "                      columns=[\"coeff\", \"k\", \"l\", \"m\", \"std\"])\n",
    "    filename = \"./PG_data/spectrumOptimum_{}_{}__{}_{}_{}\".format(kmax, nside3D, chg2ang(lk), chg2ang(ll),\n",
    "                                                                  chg2ang(lm))\n",
    "    df.to_pickle(filename)\n",
    "\n",
    "    # subtract from SMICA the background\n",
    "    SMICA_LR_NO_BG = myHyper.SMICA_LR.squeeze() - diffmap.squeeze()\n",
    "    myHyper.plot_ONLY_CL_From_Image(newmap, nside3D, SMICA_LR_NO_BG, nsidesmica=nside3D, xmax=200)\n",
    "\n",
    "    df[\"CL\"] = df[\"coeff\"] ** 2\n",
    "    df[\"abscoeff\"] = df.coeff.abs()\n",
    "\n",
    "    fcolors = df.coeff.values\n",
    "    (mu, sigma) = norm.fit(fcolors)\n",
    "    n, bins, patch = plt.hist(fcolors, 600, density=1, facecolor=\"r\", alpha=0.25)\n",
    "    y = norm.pdf(bins, mu, sigma)\n",
    "    plt.plot(bins, y)\n",
    "    plt.xlim(mu - 4 * sigma, mu + 4 * sigma)\n",
    "    plt.xlabel(\"Amplitude\")\n",
    "    plt.ylabel(\"Frequency\")\n",
    "    plt.ylim(0, 0.04)\n",
    "    plt.title(\"Amplitude Histogram \\n HU Modeling of Planck SMICA Map\")\n",
    "    plt.savefig(\"./PG_data/AmplitudeHistogram_{}_{}.png\", dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    numberofmodes = df.groupby(['k'])[\"abscoeff\"].count()\n",
    "    meanAmplitude = df.groupby(['k'])[\"abscoeff\"].mean()\n",
    "    meanEnergy = meanAmplitude ** 2 * meanAmplitude.index * numberofmodes\n",
    "    stdd = df.groupby(['k'])[\"abscoeff\"].std()\n",
    "    one = np.ones(stdd.shape)\n",
    "    valueOfK = stdd.index.values\n",
    "\n",
    "    # calculation of Predicted Amplitudes\n",
    "    predictedAmp = 75 / np.sqrt(valueOfK)\n",
    "    # low pass filtering\n",
    "    predictedAmp = np.exp(-0.01 * (valueOfK - 2)) * predictedAmp\n",
    "    # Frequency dependent freezing phase accounting\n",
    "    delta = 0.01\n",
    "    predictedAmp = sinxx(valueOfK * delta) * predictedAmp\n",
    "    plt.plot(numberofmodes, predictedAmp)\n",
    "    plt.scatter(numberofmodes, meanAmplitude)\n",
    "    plt.title(\"Predicted and Observed Mean Amplitudes \\n under Energy Equipartition per mode k\")\n",
    "    plt.xlabel(\"Accessible Modes for given k\")\n",
    "    plt.ylabel(\"Hyperspherical Mode Amplitude per k\")\n",
    "    plt.ylim(0, 25)\n",
    "    plt.xlim(0, 250)\n",
    "    plt.savefig(\"./PG_data/PredicedtAmp.png\", dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    plt.plot(numberofmodes, meanEnergy)\n",
    "    plt.title(\"Mean Energy per k\")\n",
    "    plt.xlabel(\"Accessible Modes for given k\")\n",
    "    plt.ylabel(\"Mean Energy per k \\n arbitrary units\")\n",
    "    plt.xlim(0, 250)\n",
    "    plt.ylim(0, 4.5E5)\n",
    "    plt.savefig(\"./PG_data/MeanEnergy.png\", dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    plt.scatter(numberofmodes, meanAmplitude)\n",
    "    plt.title(\"Mean Amplitude per k\")\n",
    "    plt.xlabel(\"Accessible Modes for given k\")\n",
    "    plt.ylabel(\"Mean Amplitude per k\")\n",
    "    plt.ylim(0, 22.5)\n",
    "    plt.xlim(0, 250)\n",
    "    plt.savefig(\"./PG_data/MeanAmp.png\", dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    plt.scatter(numberofmodes, stdd)\n",
    "    plt.title(\"STD\")\n",
    "    plt.xlabel(\"Accessible Modes for given k\")\n",
    "    plt.ylabel(\"STD\")\n",
    "    plt.xlim(0, 250)\n",
    "    plt.savefig(\"./PG_data/STD.png\", dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "    plt.plot(numberofmodes, meanAmplitude / stdd)\n",
    "    plt.plot(numberofmodes, np.mean((meanAmplitude / stdd)[3:]) * one)\n",
    "    plt.title(\"MeanAmplitude/std\")\n",
    "    plt.xlabel(\"Accessible Modes for given k\")\n",
    "    plt.ylabel(\"MeanAmplitude/std\")\n",
    "    plt.xlim(0, 250)\n",
    "    plt.ylim(0, 2)\n",
    "    plt.savefig(\"./PG_data/PredicedtAmpPerSTD.png\", dpi=300)\n",
    "    plt.show()\n",
    "    ################################################################\n",
    "    ################################################################\n",
    "    ################################################################\n",
    "    ################################################################\n",
    "    ################################################################\n",
    "    ################################################################\n",
    "    #################################################################\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.FINDBESTFORKRANGE in todo:\n",
    "    # myHyper.phi, myHyper.theta, myHyper.ksi = hp.pix2vec(nside=nside3D, ipix=np.arange(hp.nside2npix(nside=nside3D)))\n",
    "    # myHyper.theta, myHyper.phi = hp.pix2ang(nside=nside3D,\n",
    "    #                                                      ipix=np.arange(hp.nside2npix(nside=nside3D)))\n",
    "    # myHyper.costheta = np.cos(myHyper.theta)\n",
    "    # k = 9\n",
    "    # l = 3\n",
    "    # for m in np.arange(-l+1,l):\n",
    "    #     pp = np.array([legendre(l, x) for x in myHyper.costheta])\n",
    "    #     myHyper.plot_aitoff_df(l,m,myHyper.phi, pp=pp)\n",
    "    #\n",
    "    # x0=np.load(\"./img1/x0_{}_{}.npy\".format(kmax, nside3D), allow_pickle=True)\n",
    "    # x00 = np.load(\"./img1/x00_{}_{}.npy\".format(kmax, nside3D), allow_pickle=True)\n",
    "    # x0 = np.load(\"./img1/results_{}_{}.npy\".format(kmax, nside3D), allow_pickle=True)\n",
    "    # x0 = [5.45066325E+00, 1.36266581E+00, -2.65625000E-04]\n",
    "    (lk, ll, lm) = x0\n",
    "    myHyper.change_SMICA_resolution(nside3D, myHyper.sigma_smica, bandwidth=bandwidth)\n",
    "    olderr = 1110.0\n",
    "\n",
    "    for kk in [x for x in karray if x > 10]:\n",
    "        print(kk)\n",
    "        kkarray = [x for x in karray if x <= kk + 1]\n",
    "        kmax = np.max(kkarray)\n",
    "        results, newmap, err, packedresults = myHyper.project4D3d_0(kkarray)\n",
    "        print(lk, ll, lm, err, \"intermediate result\")\n",
    "        if err < olderr:\n",
    "            olderr = err\n",
    "            filename = \"./img1/Best_{}_{}__{}_{}_{}_{}.png\".format(kmax, nside3D, np.round(lk, 1),\n",
    "                                                                   np.round(ll, 1), np.round(lm, 1),\n",
    "                                                                   np.round(err * 1, 4))\n",
    "            title = \"Best_{}_{}__{}_{}_{}_{}\".format(kmax, nside3D, np.round(lk, 1),\n",
    "                                                     np.round(ll, 1), np.round(lm, 1),\n",
    "                                                     np.round(err * 1, 4))\n",
    "            myHyper.plotNewMap(newmap, err, filename=filename, title=title, plotme=True, save=True)\n",
    "            print(lk, ll, lm, err, title)\n",
    "            filenameResults = \"./PG_data/Results_{}_{}__{}_{}_{}_{}.npy\".format(kmax, nside3D, np.round(lk, 1),\n",
    "                                                                                np.round(ll, 1),\n",
    "                                                                                np.round(lm, 1),\n",
    "                                                                                np.round(err * 1, 4))\n",
    "            np.save(filenameResults, packedresults)\n",
    "            np.save(\"./PG_data/resultsBest_{}_{}.npy\".format(kmax, nside3D), results, allow_pickle=True)\n",
    "            myHyper.plotHistogram(newmap.squeeze(), nside3D, kmax, plotme=True)\n",
    "            myHyper.plotNewMap(newmap, err, filename=filename, title=title, plotme=True, save=False)\n",
    "    # dl_SMICA, dl_SMICA_HU, ell = myHyper.plot_ONLY_CL_From_Image(newmap.squeeze(), bandwidth, xmax=10 * kmax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Color.CREATE_GAUSSIAN_BACKGROUND in todo:\n",
    "    nside3D = 128  # 128\n",
    "    bandwidth = 3000\n",
    "    filename = \"./img1/3DGaussianNoise_{}_{}.png\".format(kmax, nside3D)\n",
    "    title = \"3DGaussianNoise_{}_{}\".format(kmax, nside3D)\n",
    "    (lambda_k, lambda_l, lambda_m) = x0\n",
    "    # Create karray\n",
    "    karray = list(np.arange(2, 50))\n",
    "    # kmax=1\n",
    "    # karray += sorted(list(set([int(k) for k in np.geomspace(11,kmax,30)])))\n",
    "    print(len(karray))\n",
    "    kmax = max(karray)\n",
    "    #################################\n",
    "    myHyper = HYPER(nside3D, sigma_smica, planck_IQU_SMICA, karray,\n",
    "                    lambda_k, lambda_l, lambda_m, loadpriorG=True, savePG=True)\n",
    "    #################################################################\n",
    "    delta = 0.001\n",
    "    lmax = 1000\n",
    "    llist = sorted(list(set([int(k) for k in np.geomspace(1, lmax + 1, 100)])))\n",
    "    wavefunction = myHyper.createGaussianBackground(x0, nside3D, delta, karray=llist)\n",
    "    wavefunction = myHyper.normalizeFColors(wavefunction, myHyper.sigma_smica).squeeze()\n",
    "    myHyper.plotNewMap(wavefunction, 0, filename=filename, title=title, plotme=True, save=True)\n",
    "    myHyper.plot_CL_From_Image(wavefunction, nside3D, planck_theory_cl, xmax=2 * lmax, ymax=0.01)\n",
    "    myHyper.plot_ONLY_CL_From_Image(wavefunction, nside3D, myHyper.SMICA, nsidesmica=1024, xmax=2 * lmax)\n",
    "    np.save(\"./PG_data/3DGaussianNoise_{}_{}.npy\".format(kmax, nside3D), wavefunction, allow_pickle=True)\n",
    "###############################################################\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:Cosmos3020]",
   "language": "python",
   "name": "conda-env-Cosmos3020-py"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}