MDN-multiVariate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 55,
   "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",
    "import glob, os\n",
    "\n",
    "#os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"0\"  #uncomment for GPU use\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import IPython\n",
    "from IPython.display import HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "import matplotlib \n",
    "import matplotlib.tri as tri\n",
    "\n",
    "import tensorflow as tf\n",
    "import mdn\n",
    "from tensorflow import keras\n",
    "\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from scipy.stats import multivariate_normal\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import binned_statistic\n",
    "import math as m\n",
    "from joblib import Parallel, delayed, dump, load\n",
    "from keras.callbacks import History "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathSave = \".\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "class getData:\n",
    "    def __init__(self, _paraVecIndices, _xparaVecIndices, _tProfBounds, _trainPercent, _pathSave, _noiseLevel=0.):\n",
    "        self.paraVecIndices = _paraVecIndices\n",
    "        self.xparaVecIndices = _xparaVecIndices\n",
    "        self.tProfBounds = _tProfBounds\n",
    "        self.trainPercent = _trainPercent\n",
    "        self.pathSave = _pathSave\n",
    "        self.noiseLevel = _noiseLevel\n",
    "        \n",
    "    def visualizeData(self, picID, fontSize, Dict_baseline = {}, plotTProf=True, plotBaseline=False): \n",
    "        x_data, x_test, x_cv, y_data, y_test, y_cv, paraVec, xparaVec, rProf, startPoint, endPoint, indexer \\\n",
    "        = self.x_data, self.x_test, self.x_cv, self.y_data, self.y_test, self.y_cv, self.paraVec, self.xparaVec, \\\n",
    "            self.rProf, self.startPoint, self.endPoint, self.indexer    \n",
    "                \n",
    "        print('Added noise ' + str(self.noiseLevel))\n",
    "        print(\"X_train shape: \" + str(x_data.shape))\n",
    "        print(\"Y_train shape: \" + str(y_data.shape))\n",
    "        print(\"X_test shape: \" + str(x_test.shape))\n",
    "        print(\"Y_test shape: \" + str(y_test.shape))\n",
    "        print(\"X_cv shape: \" + str(x_cv.shape))\n",
    "        print(\"Y_cv shape: \" + str(y_cv.shape))\n",
    "        print()\n",
    "        \n",
    "        alphaVal = 0.1\n",
    "        tProfBoolean = False\n",
    "\n",
    "        for ind in range(np.size(paraVec)):\n",
    "            fig, ax = plt.subplots()\n",
    "            n, bins, patches = ax.hist(y_data[:, ind], 100, density=1)\n",
    "            ax.set_xlabel(paraVec[ind])\n",
    "            fig = plt.figure(figsize=(14,4), dpi=200)\n",
    "            nCols = np.size(xparaVec)\n",
    "            nRows = 1\n",
    "            for indX in range(nCols):\n",
    "                ax = fig.add_subplot(nRows,nCols,indX+1)\n",
    "                ax.plot(x_data[:,indX], y_data[:, ind],'ro', alpha=alphaVal)\n",
    "                if plotBaseline:\n",
    "                    ax.plot(x_data[:,indX],Dict_baseline['mu_var' + str(indX)] \\\n",
    "                                 +Dict_baseline['variance_var' + str(indX)],'k.')\n",
    "                    ax.plot(x_data[:,indX],Dict_baseldine['mu_var' + str(indX)] \\\n",
    "                                 -Dict_baseline['variance_var' + str(indX)],'k.')\n",
    "                ax.yaxis.labelpad = 0.8\n",
    "                ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))\n",
    "                plt.title(\"train\")\n",
    "                plt.xlabel(xparaVec[indX])\n",
    "                plt.ylabel(paraVec[ind])\n",
    "            plt.tight_layout()\n",
    "            plt.show()\n",
    "            #fig.savefig(pathSave + \"inverseProblemPics/\" + \"Input0_fig\" + str(picID) + \".pdf\", bbox_inches='tight')\n",
    "            \n",
    "        if  tProfBoolean and plotTProf:   \n",
    "            length = int(np.size(paraVec))\n",
    "            fig = plt.figure(figsize=(length,4), dpi=200)\n",
    "            nCols = np.size(paraVec)\n",
    "            nRows = 2\n",
    "            for indY in range(nCols):\n",
    "                vecInterest = y_data[:, indY]\n",
    "                colorsVec= [plt.cm.jet(i) for i in vecInterest]\n",
    "                ax = fig.add_subplot(nRows,nCols,indY+1)\n",
    "                ax.set_prop_cycle('color',colorsVec)\n",
    "                for i, ii in enumerate(indexer):\n",
    "                    rProf = self.rProf[startPoint:endPoint] # profiles['Dict_Rprof_4p5Gyr' + str(ii)]\n",
    "                    rProf = (rProf-rProf.min())/(rProf.max()-rProf.min())\n",
    "                    tProf = x_data[i,np.size(xparaVec):]\n",
    "                    ax.plot(tProf, rProf, linewidth=0.5)\n",
    "                ax.yaxis.labelpad = 0.8\n",
    "                ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))\n",
    "                if indY > 0:\n",
    "                    ax.axes.get_yaxis().set_ticks([])\n",
    "                ax.set_title(paraVec[indY])\n",
    "                ax.set_xlabel(\"T\") #[K]\")\n",
    "                ax.set_ylabel(\"R\") #[km] \")\n",
    "                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +\n",
    "                     ax.get_xticklabels() + ax.get_yticklabels()):\n",
    "                    item.set_fontsize(12)\n",
    "                #ax.set_xlim([250,2250])\n",
    "                #ax.set_ylim([1700,3400])\n",
    "            plt.tight_layout()\n",
    "            plt.show()\n",
    "           # fig.savefig(pathSave + \"inverseProblemPics/\" + \"Input1_fig\" + str(picID) + \".pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "class plotResults:\n",
    "    def __init__(self, data, ax, ax1, alphaVal, title, _y_actual, Dict_prob_paras ,indPara,\\\n",
    "                 Dict_baseline, kernel, KMIX, plotIndividualPDFs=False):\n",
    "        \n",
    "        numParameters = len(data.paraVec)\n",
    "        \n",
    "        if \"Test\" in title:\n",
    "            self.Dict_Test = Dict_prob_paras\n",
    "        \n",
    "        rho_m  = 3500. \n",
    "        g     = 3.7 \n",
    "        alpha_m = 2.5e-5\n",
    "        T_delta = 2000. \n",
    "        D = 1700e+3         \n",
    "        k_diffusive = 1e-6 \n",
    "        R = 8.314   \n",
    "        \n",
    "        def format_func(_val, tick_number):\n",
    "            f = mticker.ScalarFormatter(useOffset=False, useMathText=True)\n",
    "            _g = lambda x,pos : \"${}$\".format(f._formatSciNotation('%1.2e' % x))\n",
    "            fmt = mticker.FuncFormatter(_g)\n",
    "            return \"{}\".format(fmt(_val))\n",
    "\n",
    "        def dimensionalize(_val,_ind,isVariance=False):\n",
    "            _min = data.pMin[_ind]\n",
    "            _max = data.pMax[_ind]\n",
    "            _val = _val*(_max-_min) \n",
    "            if not isVariance:\n",
    "                _val = _val + _min\n",
    "            \n",
    "            if \"Ra\" in title and not isVariance:\n",
    "                _val = np.log10(rho_m * g * alpha_m * T_delta * np.power(D,3.)/(np.power(10.,_val) * k_diffusive))\n",
    "            \n",
    "            if \"ERef\" in title:\n",
    "                _val = _val*(R * T_delta)\n",
    "            \n",
    "            if \"VRef\" in title:\n",
    "                _val = _val*(R * T_delta) /(rho_m * g * D)\n",
    "            \n",
    "            if \"iniTempTop\" in title:\n",
    "                _val = _val*2000 \n",
    "                if not isVariance:\n",
    "                    _val = _val + 250\n",
    "                    \n",
    "            return _val\n",
    "\n",
    "        xx = dimensionalize(_y_actual,indPara)    \n",
    "        variance = []\n",
    "        combined_mdn_std = []\n",
    "        combined_mdn_mean = []\n",
    "        xxSorted = np.sort(xx)        \n",
    "        \n",
    "        xxSize = np.size(xx)\n",
    "        plotEvery = int(np.floor(xxSize*0.1))\n",
    "        xxSorted.shape = (xxSize,1)\n",
    "        Pr_xx = np.zeros((xxSize,xxSize))\n",
    "        Pr_baseline = np.zeros((xxSize,xxSize))\n",
    "        colors = ['r', 'g', 'm']\n",
    "        \n",
    "        def pdf(x):\n",
    "            return 1./np.sqrt(2.*np.pi) * np.exp(-x**2/2.)\n",
    "        #def cdf(x):\n",
    "        #    return (1. + erf(x/np.sqrt(2))) / 2.\n",
    "        def SSkew(x,a,e,w):\n",
    "            t = (x-e) / w\n",
    "            return 2. / w * pdf(t) * cdf(a*t)\n",
    "        def getGamma(x):\n",
    "            return m.gamma(x)\n",
    "        def tf_beta(y, alpha, beta):\n",
    "            Z = np.divide(np.multiply(getGamma(alpha),getGamma(beta)), getGamma(alpha + beta))\n",
    "            result = np.divide(np.multiply(np.power(y,(alpha - 1.)), np.power((1. - y),(beta - 1.))),Z)\n",
    "            return  result\n",
    "             \n",
    "        if plotIndividualPDFs:\n",
    "            fig3 = plt.figure(figsize=(14,18))\n",
    "           # plt.subplots_adjust(top=0.90, bottom=0.08, left=0.10, right=0.95, hspace=0.3,wspace=0.4)\n",
    "            nRows = np.ceil((xxSize/plotEvery)/3)\n",
    "            nCols = 3\n",
    "            plotCounter = 1                    \n",
    "                \n",
    "        if kernel=='gaussian':\n",
    "            for ind,val in enumerate(xxSorted):\n",
    "                index = np.where(val == xx)[0][0]\n",
    "                muIntermediate = []\n",
    "                piIntermediate = []\n",
    "                sigmaIntermediate = []\n",
    "                for i in range(KMIX):\n",
    "                    out_pi = np.asarray(Dict_prob_paras[\"multi_pi\"+str(i)])[index]\n",
    "                    out_sigma = np.asarray(Dict_prob_paras[\"multi_sigma_\" + str(index) + \"_\" + str(i)])\n",
    "                    out_cov = np.matmul(out_sigma,np.transpose(out_sigma))\n",
    "                    out_mu = np.asarray(Dict_prob_paras[\"multi_mu\"+str(i)])[index,:]\n",
    "                    \n",
    "                    muIntermediate.append(dimensionalize(out_mu[indPara],indPara))\n",
    "                    sigmaIntermediate.append(dimensionalize(np.sqrt(out_cov[indPara,indPara]),indPara,True))\n",
    "                    piIntermediate.append(out_pi)\n",
    "\n",
    "                pr = piIntermediate[0] * norm.pdf(xxSorted, muIntermediate[0], sigmaIntermediate[0])\n",
    "                for i in range(1,KMIX):\n",
    "                    pr += piIntermediate[i] * norm.pdf(xxSorted, muIntermediate[i], sigmaIntermediate[i])\n",
    "                pr.shape = (pr.shape[0],)\n",
    "                \n",
    "                mean_mdn = 0\n",
    "                for i in range(KMIX):\n",
    "                    mean_mdn += piIntermediate[i]*muIntermediate[i]\n",
    "                var_mdn = 0\n",
    "                #http://ecee.colorado.edu/~pao/anonftp/cdc02gaussmix.pdf: Equation 11\n",
    "                for i in range(KMIX):\n",
    "                    var_mdn += piIntermediate[i]*(np.power(sigmaIntermediate[i],2) + np.power(muIntermediate[i],2)) \n",
    "                var_mdn -= np.power(mean_mdn,2)\n",
    "                \n",
    "                combined_mdn_std.append(np.power(var_mdn,0.5))\n",
    "                combined_mdn_mean.append(mean_mdn)\n",
    "                \n",
    "                Pr_xx[ind,:] = pr \n",
    "                if  plotIndividualPDFs and ind%plotEvery == 0:\n",
    "                    ax3 = fig3.add_subplot(nRows,nCols,plotCounter) \n",
    "                    plotCounter += 1\n",
    "                    legendStr = []\n",
    "                    for i in range(np.size(muIntermediate)):\n",
    "                        if piIntermediate[i] >= np.max(piIntermediate)*1e-15:\n",
    "                            yPDF = piIntermediate[i] * norm.pdf(xxSorted, muIntermediate[i], sigmaIntermediate[i])\n",
    "                            ax3.plot(xxSorted,yPDF,colors[i]+'--', linewidth=2.5)\n",
    "                            legendStr.append(\"%.4f\" % piIntermediate[i])\n",
    "                            if \"Test\" in title:\n",
    "                                self.Dict_Test[\"forFig_\" + title + \"_pdf_mixture\" + str(i) + \"_case_\" + str(ind)] = yPDF\n",
    "                    ax3.plot(xxSorted, Pr_xx[ind,:],\"k-\", linewidth=2.5)\n",
    "                    ax3.plot([val, val],[0,max(Pr_xx[ind,:])], \"-\", color='grey', linewidth=2.5)\n",
    "                    \n",
    "                    if \"Test\" in title:\n",
    "                        self.Dict_Test[\"forFig_\" + title + \"_trueVal_\" + str(ind)] = val\n",
    "                    \n",
    "                    for item in ([ax3.title, ax3.xaxis.label, ax3.yaxis.label] +\n",
    "                     ax3.get_xticklabels() + ax3.get_yticklabels()):\n",
    "                        item.set_fontsize(20)\n",
    "                    if \"ERef\" in title:\n",
    "                        ax3.set_xticks([1e+5,3e+5,5e+5])\n",
    "                        ax3.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                        ax3.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                    elif \"VRef\" in title:\n",
    "                        ax3.set_xticks([4e-6, 7e-6, 10e-6])\n",
    "                        ax3.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                        ax3.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                    plt.tight_layout()\n",
    "        \n",
    "        if \"Test\" in title:\n",
    "            self.Dict_Test[\"forFig_\" + title + \"_xxSorted\"] = xxSorted\n",
    "            self.Dict_Test[\"forFig_\" + title + \"_Pr_xx\"] = Pr_xx\n",
    "        \n",
    "        for ind,val in enumerate(xxSorted):\n",
    "            index = np.where(val == xx)[0][0]\n",
    "            muBase = dimensionalize(Dict_baseline[\"mu_var\" + str(0) + \"_\" +  str(indPara)][index],indPara)\n",
    "            sigmaBase = dimensionalize(Dict_baseline[\"variance_var\" + str(0) + \"_\" +  str(indPara)][index],indPara,True)\n",
    "            prb = 0\n",
    "            prb = norm.pdf(xxSorted, muBase, sigmaBase)  \n",
    "            prb.shape = (prb.shape[0],)\n",
    "            Pr_baseline[ind,:] = prb \n",
    "            \n",
    "        if \"Test\" in title:\n",
    "            self.Dict_Test[\"forFig_\" + title + \"_averageSTD\"] = np.mean(np.asarray(combined_mdn_std))            \n",
    "\n",
    "        print(\"Average standard deviation: \" + str(np.mean(np.asarray(combined_mdn_std))))\n",
    "        \n",
    "        if \"Ra\" in title:\n",
    "            titlep = \"$\\log(\\eta_{ref})$ [Pa s]\"\n",
    "        elif \"ERef\" in title:\n",
    "            titlep = r\"$E$ [J mol$^{-1}$]\"\n",
    "        elif \"VRef\" in title:\n",
    "            titlep = r\"$V$ [m$^3$ mol$^{-1}$]\"\n",
    "        elif \"Enrichment_cr\" in title:\n",
    "            titlep = \"$\\Lambda$\"\n",
    "        elif \"iniTempTop\" in title:\n",
    "            titlep = \"$T_{ini}$ [K]\"\n",
    "        \n",
    "        x = np.zeros((np.size(xxSorted),np.size(xxSorted)))\n",
    "        y = np.zeros((np.size(xxSorted),np.size(xxSorted)))\n",
    "        for ind,val in enumerate(xxSorted):\n",
    "            x[ind,:] = val\n",
    "            y[:,ind] = xxSorted[ind]\n",
    "        ax.contourf(x,y,Pr_xx)\n",
    "        \n",
    "        if \"Test\" in title:\n",
    "            self.Dict_Test[\"forFig_\" + title + \"x\"] = x\n",
    "            self.Dict_Test[\"forFig_\" + title + \"y\"] = y\n",
    "        \n",
    "        ax.set_xlabel(\"True\")\n",
    "        ax.set_ylabel(\"Predicted\")\n",
    "        ax.set_title(titlep + \"; MDN\")\n",
    "        \n",
    "        ax1.contourf(x,y,Pr_baseline)  \n",
    "        ax1.set_xlabel(\"True\")\n",
    "        ax1.set_ylabel(\"Predicted\")\n",
    "        ax1.set_title(titlep + \"; MP\")\n",
    "        print()\n",
    "        setTicks = True\n",
    "        if setTicks:\n",
    "            if \"Ra\" in title:\n",
    "                ax.set_xticks([19,20,21,22])\n",
    "                ax.set_yticks([19,20,21,22])\n",
    "                ax1.set_xticks([19,20,21,22])\n",
    "                ax1.set_yticks([19,20,21,22])\n",
    "\n",
    "            if \"ERef\" in title:\n",
    "                ax.set_xticks([1e+5,3e+5,5e+5])\n",
    "                ax.set_yticks([1e+5,2e+5,3e+5,4e+5,5e+5])\n",
    "                ax1.set_xticks([1e+5,3e+5,5e+5])\n",
    "                ax1.set_yticks([1e+5,2e+5,3e+5,4e+5,5e+5])\n",
    "                ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax1.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax1.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "            elif \"VRef\" in title:\n",
    "                ax.set_xticks([4e-6, 7e-6, 10e-6])\n",
    "                ax.set_yticks([4e-6, 5e-6, 6e-6, 7e-6, 8e-6, 9e-6, 10e-6])\n",
    "                ax1.set_xticks([4e-6, 7e-6, 10e-6])\n",
    "                ax1.set_yticks([4e-6, 5e-6, 6e-6, 7e-6, 8e-6, 9e-6, 10e-6])\n",
    "                ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax1.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                ax1.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "\n",
    "            elif \"Enrichment_cr\" in title:\n",
    "                ax.set_xticks([1,10,20,30,40,50])\n",
    "                ax.set_yticks([1,10,20,30,40,50])\n",
    "                ax1.set_xticks([1,10,20,30,40,50])\n",
    "                ax1.set_yticks([1,10,20,30,40,50])\n",
    "            elif \"iniTempTop\" in title:\n",
    "                ax.set_xticks([1600,1650,1700,1750,1800])\n",
    "                ax.set_yticks([1600,1650,1700,1750,1800])\n",
    "                ax1.set_xticks([1600,1650,1700,1750,1800])\n",
    "                ax1.set_yticks([1600,1650,1700,1750,1800])\n",
    "\n",
    "            for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +\n",
    "             ax.get_xticklabels() + ax.get_yticklabels()):\n",
    "                item.set_fontsize(20)\n",
    "\n",
    "            for item in ([ax1.title, ax1.xaxis.label, ax1.yaxis.label] +\n",
    "             ax1.get_xticklabels() + ax1.get_yticklabels()):\n",
    "                item.set_fontsize(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "class plotResults_multi:\n",
    "    def __init__(self, data, Dict_prob_paras, KMIX, id_string):\n",
    "        \n",
    "        rho_m  = 3500. \n",
    "        g     = 3.7 \n",
    "        alpha_m = 2.5e-5\n",
    "        T_delta = 2000. \n",
    "        D = 1700e+3         \n",
    "        k_diffusive = 1e-6 \n",
    "        R = 8.314   \n",
    "        \n",
    "        numParameters = len(data.paraVec)\n",
    "        plotParas = [r\"$\\log(\\eta_{ref})$ [Pa s]\", r\"$E$ [J mol$^{-1}$]\", r\"$V$ [m$^3$ mol$^{-1}$]\", \n",
    "                                     \"$\\Lambda$\", \"$T_{ini}$ [K]\"]\n",
    "        \n",
    "        def format_func(_val, tick_number):\n",
    "            f = mticker.ScalarFormatter(useOffset=False, useMathText=True)\n",
    "            _g = lambda x,pos : \"${}$\".format(f._formatSciNotation('%1.2e' % x))\n",
    "            fmt = mticker.FuncFormatter(_g)\n",
    "            return \"{}\".format(fmt(_val))\n",
    "\n",
    "        def dimensionalize(_val,_ind,isVariance=False):\n",
    "            _min = data.pMin[_ind]\n",
    "            _max = data.pMax[_ind]\n",
    "            _val = _val*(_max-_min) \n",
    "            if not isVariance:\n",
    "                _val = _val + _min\n",
    "            \n",
    "            if _ind==0 and not isVariance:\n",
    "                _val = np.log10(rho_m * g * alpha_m * T_delta * np.power(D,3.)/(np.power(10.,_val) * k_diffusive))\n",
    "            if _ind==1:\n",
    "                _val = _val*(R * T_delta)\n",
    "            if _ind==2:\n",
    "                _val = _val*(R * T_delta) /(rho_m * g * D)\n",
    "            if _ind==4:\n",
    "                _val = _val*2000 \n",
    "                if not isVariance:\n",
    "                    _val = _val + 250    \n",
    "                    \n",
    "            return _val\n",
    "        \n",
    "        def multivariate_gaussian(pos, _mu, Sigma):\n",
    "            #https://scipython.com/blog/visualizing-the-bivariate-gaussian-distribution/\n",
    "            \"\"\"Return the multivariate Gaussian distribution on array pos.\n",
    "            pos is an array constructed by packing the meshed arrays of variables\n",
    "            x_1, x_2, x_3, ..., x_k into its _last_ dimension.\n",
    "\n",
    "            \"\"\"\n",
    "            n = _mu.shape[0]\n",
    "            Sigma_det = np.linalg.det(Sigma)\n",
    "            Sigma_inv = np.linalg.inv(Sigma)\n",
    "            N = np.sqrt((2*np.pi)**n * Sigma_det)\n",
    "            # This einsum call calculates (x-mu)T.Sigma-1.(x-mu) in a vectorized\n",
    "            # way across all the input variables.\n",
    "            fac = np.einsum('...k,kl,...l->...', pos-_mu, Sigma_inv, pos-_mu)\n",
    "            return np.exp(-fac / 2) / N\n",
    "        \n",
    "        figl = plt.figure(figsize=(15,3), dpi=320)\n",
    "        pltcount = 1\n",
    "        for paraInd in range(data.y_test.shape[1]):\n",
    "            axl = figl.add_subplot(1,len(data.paraVec),pltcount)\n",
    "            x = dimensionalize(data.y_test[:,paraInd],paraInd)\n",
    "            rv = np.zeros((x.shape[0],1))\n",
    "            for xind, xval in enumerate(x):\n",
    "                rv[xind,0] = Dict_prob_paras[\"log_likelihood\"][xind]\n",
    "            axl.plot(x, rv, '.')\n",
    "            axl.set_xlabel(plotParas[paraInd])\n",
    "            if paraInd==0:\n",
    "                axl.set_ylabel(\"Log-likelihood\")\n",
    "            axl.set_ylim([-2,16])\n",
    "            setTicks = True\n",
    "            if setTicks:\n",
    "                axl.set_yticks([0,4,8,12,16])\n",
    "                if paraInd==0:\n",
    "                    axl.set_xticks([19,20,21,22])\n",
    "                if paraInd==1:\n",
    "                    axl.set_xticks([1e+5,3e+5,5e+5])\n",
    "                    axl.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                if paraInd==2:\n",
    "                    axl.set_xticks([4e-6, 7e-6, 10e-6])\n",
    "                    axl.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                if paraInd==3:\n",
    "                    axl.set_xticks([1,10,20,30,40,50])\n",
    "                if paraInd==4:\n",
    "                    axl.set_xticks([1600,1700,1800])\n",
    "            \n",
    "            for item in ([axl.title, axl.xaxis.label, axl.yaxis.label] +\n",
    "             axl.get_xticklabels() + axl.get_yticklabels()):\n",
    "                item.set_fontsize(18)\n",
    "            pltcount += 1\n",
    "        \n",
    "        plt.tight_layout(w_pad=0.1)\n",
    "        plt.show()\n",
    "            \n",
    "        for indInterest in [-1,9,14]: \n",
    "            figm = plt.figure(figsize=(14,12), dpi=320)\n",
    "            pltcount = 1\n",
    "\n",
    "            for ind_x, d_x in enumerate(data.paraVec):\n",
    "                for ind_y, d_y in enumerate(data.paraVec):\n",
    "                    if ind_x < ind_y:\n",
    "                        ax2 = figm.add_subplot(len(data.paraVec),len(data.paraVec),pltcount)\n",
    "                        if indInterest==-1:\n",
    "                            x = dimensionalize(data.y_test[:,ind_x],ind_x)\n",
    "                            y = dimensionalize(data.y_test[:,ind_y],ind_y)\n",
    "                            z = Dict_prob_paras[\"log_likelihood\"]\n",
    "                            ax2.scatter(x, y, c=z, s=50, cmap=\"Greys\")\n",
    "                            ax2.scatter(x[9], y[9], s=80, facecolors='none', edgecolors='b', linewidth=2.5)\n",
    "                            ax2.scatter(x[14], y[14], s=80, facecolors='none', edgecolors='r', linewidth=2.5)\n",
    "                            ax2.plot(x, y, color='none')\n",
    "                            ax2.relim()\n",
    "                            ax2.autoscale_view()\n",
    "                        else:\n",
    "                            actVal_x = dimensionalize(data.y_test[indInterest,ind_x],ind_x)\n",
    "                            actVal_y = dimensionalize(data.y_test[indInterest,ind_y],ind_y)\n",
    "                            x = dimensionalize(np.linspace(0,1,1000),ind_x)\n",
    "                            y = dimensionalize(np.linspace(0,1,1000),ind_y)\n",
    "                            x_, y_ = np.meshgrid(x, y)\n",
    "                            pos = np.dstack((x_, y_))\n",
    "                            rv=0\n",
    "                            for i in range(KMIX):\n",
    "                                out_pi = np.asarray(Dict_prob_paras[\"multi_pi\"+str(i)])[indInterest]\n",
    "                                out_sigma = np.asarray(Dict_prob_paras[\"multi_sigma_\" + str(indInterest) + \"_\" + str(i)])\n",
    "                                out_mu = np.asarray(Dict_prob_paras[\"multi_mu\"+str(i)])[indInterest,:]\n",
    "                                out_cov = np.matmul(out_sigma, np.transpose(out_sigma))\n",
    "                                cov_matrix = [[out_cov[ind_x,ind_x],out_cov[ind_x,ind_y]], \n",
    "                                              [out_cov[ind_y,ind_x],out_cov[ind_y,ind_y]]]  # scale matrix from diag\n",
    "\n",
    "                                dim_x = data.pMax[ind_x]-data.pMin[ind_x]\n",
    "                                dim_y = data.pMax[ind_y]-data.pMin[ind_y]                                \n",
    "\n",
    "                                if ind_x==1:\n",
    "                                    dim_x *= (R * T_delta)\n",
    "                                if ind_x==2:\n",
    "                                    dim_x *= (R * T_delta) /(rho_m * g * D)\n",
    "                                if ind_x==4:\n",
    "                                    dim_x *= 2000\n",
    "                                if ind_y==1:\n",
    "                                    dim_y *= (R * T_delta)\n",
    "                                if ind_y==2:\n",
    "                                    dim_y *= (R * T_delta) /(rho_m * g * D)\n",
    "                                if ind_y==4:\n",
    "                                    dim_y *= 2000 \n",
    "\n",
    "                                cov_matrix[0][0] *= dim_x*dim_x \n",
    "                                cov_matrix[0][1] *= dim_x*dim_y\n",
    "                                cov_matrix[1][0] *= dim_y*dim_x \n",
    "                                cov_matrix[1][1] *= dim_y*dim_y\n",
    "\n",
    "                                rv += out_pi*multivariate_gaussian(pos,np.asarray([dimensionalize(out_mu[ind_x],ind_x),\n",
    "                                                                  dimensionalize(out_mu[ind_y],ind_y)]), \n",
    "                                                                 np.asarray(cov_matrix))\n",
    "                            ax2.contour(x_, y_, rv)\n",
    "                            if indInterest==9:\n",
    "                                colorUse = \"b\"\n",
    "                            else:\n",
    "                                colorUse = \"r\"\n",
    "                            ax2.plot(actVal_x,actVal_y,colorUse+'X',markersize=9)\n",
    "                        ax2.set_xlabel(plotParas[ind_x])\n",
    "                        ax2.set_ylabel(plotParas[ind_y])\n",
    "                        setTicks = True\n",
    "                        if setTicks:\n",
    "                            if ind_x==0:\n",
    "                                ax2.set_xticks([19,20,21,22])\n",
    "                            if ind_y==0:\n",
    "                                ax2.set_yticks([19,20,21,22])\n",
    "                            if ind_x==1:\n",
    "                                ax2.set_xticks([1e+5,3e+5,5e+5])\n",
    "                                ax2.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                            if ind_y==1:\n",
    "                                ax2.set_yticks([1e+5,3e+5,5e+5])\n",
    "                                ax2.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                            if ind_x==2:\n",
    "                                ax2.set_xticks([4e-6, 7e-6, 10e-6])\n",
    "                                ax2.xaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                            if ind_y==2:\n",
    "                                ax2.set_yticks([4e-6, 7e-6, 10e-6])\n",
    "                                ax2.yaxis.set_major_formatter(plt.FuncFormatter(format_func))\n",
    "                            if ind_x==3:\n",
    "                                ax2.set_xticks([1,10,20,30,40,50])\n",
    "                            if ind_y==3:\n",
    "                                ax2.set_yticks([1,10,20,30,40,50])\n",
    "                            if ind_x==4:\n",
    "                                ax2.set_xticks([1600,1650,1700,1750,1800])\n",
    "                            if ind_y==4:\n",
    "                                ax2.set_yticks([1600,1650,1700,1750,1800])\n",
    "\n",
    "                            for item in ([ax2.title, ax2.xaxis.label, ax2.yaxis.label] +\n",
    "                             ax2.get_xticklabels() + ax2.get_yticklabels()):\n",
    "                                item.set_fontsize(15)\n",
    "                    pltcount += 1\n",
    "                        \n",
    "            plt.tight_layout(w_pad=0.1)\n",
    "            plt.show()\n",
    "            \n",
    "            if indInterest==-1:\n",
    "                vecInterest = Dict_prob_paras[\"log_likelihood\"]\n",
    "                vecInterestTemp = np.asarray(vecInterest)\n",
    "                vecInterest = np.asarray(vecInterest)\n",
    "                vecInterest = np.interp(vecInterest, (vecInterest.min(), vecInterest.max()), (0, 1))\n",
    "\n",
    "                figc, axc = plt.subplots(figsize=(6,0.2), dpi=320)\n",
    "                gradient= np.sort(vecInterest)\n",
    "                gradient = np.vstack((gradient, gradient))\n",
    "                axc.imshow(gradient, aspect='auto', cmap=plt.get_cmap('Greys'))\n",
    "                axc.set_yticklabels([])\n",
    "                axc.set_yticks([])\n",
    "                #ax.set_xscale('log')\n",
    "\n",
    "                my_xticks = axc.get_xticks()\n",
    "                plt.xticks([my_xticks[1], my_xticks[-2]],['{:.1f}'.format(min(vecInterestTemp)), '{:.1f}'.format(max(vecInterestTemp))]\n",
    "                           ,visible=True, rotation=\"horizontal\")\n",
    "                plt.xlabel('Log-likelihood')\n",
    "\n",
    "            #picDirectory = pathSave + \"Pics/\" + id_string + \"/\"\n",
    "            #if not os.path.exists(picDirectory):\n",
    "            #    os.makedirs(picDirectory)\n",
    "            #figm.savefig(picDirectory + \"multiVar_\" + str(indInterest) + \".pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MDN_K:\n",
    "    def __init__(self, data, x_data,y_data, x_test,y_test, x_cv,y_cv, hSize, KMIX, NEPOCH, learnRate, \\\n",
    "                 paraVec, xparaVec, Dict_baseline_train, Dict_baseline_test, picID, repeats, id_string, \\\n",
    "                 kernel='gaussian',activationFunc='tanh',trainOrload='load'):\n",
    "        \n",
    "        print(\"Hidden Layers: \" + str(hSize))\n",
    "        print(\"Number of Mixtures: \" + str(KMIX))\n",
    "        print('Kernel is ' + str(kernel))\n",
    "        print('Activation is ' + str(activationFunc))\n",
    "        print()\n",
    "            \n",
    "        print(\"Training data: \" + str(x_data.shape) + \",\" + str(y_data.shape))\n",
    "        print(\"Test data: \" + str(x_test.shape) + \",\" + str(y_test.shape))\n",
    "        print(\"CV data: \" + str(x_cv.shape) + \",\" + str(y_cv.shape))\n",
    "        print()\n",
    "        \n",
    "        yNumParameters = y_data.shape[1] \n",
    "        graph = tf.get_default_graph()        \n",
    "        l2reg = keras.regularizers.l2(0.00001/x_data.shape[0])\n",
    "        xSize = x_data.shape[1]\n",
    "        \n",
    "        model = keras.Sequential()\n",
    "        model.add(keras.layers.Dense(hSize[0], input_dim=(xSize), activation=activationFunc, kernel_regularizer=l2reg, dtype=\"float64\"))\n",
    "        for h in hSize[1:]:\n",
    "            model.add(keras.layers.Dense(h, activation=activationFunc, kernel_regularizer=l2reg, dtype=\"float64\"))\n",
    "        model.add(mdn.MDN(yNumParameters, KMIX))\n",
    "\n",
    "        class TimeHistory(keras.callbacks.Callback):\n",
    "            def on_train_begin(self, logs={}):\n",
    "                self.times = []\n",
    "            def on_epoch_begin(self, batch, logs={}):\n",
    "                self.epoch_time_start = time.time()\n",
    "            def on_epoch_end(self, batch, logs={}):\n",
    "                self.times.append(time.time() - self.epoch_time_start)\n",
    "\n",
    "        config = tf.ConfigProto(\n",
    "        intra_op_parallelism_threads=1,\n",
    "        allow_soft_placement=True)\n",
    "\n",
    "        config.gpu_options.allow_growth = True\n",
    "        config.gpu_options.per_process_gpu_memory_fraction = 0.8\n",
    "\n",
    "        session = tf.Session(config=config)\n",
    "        name = id_string + \"_multiVariate\"\n",
    "        with session.as_default():\n",
    "            with session.graph.as_default():\n",
    "                pathSaveK = pathSave + \"/Data_Files/\"\n",
    "                if trainOrload=='load':\n",
    "                    model =  keras.models.load_model(pathSaveK + name + \".hdf5\",\\\n",
    "                            custom_objects={'MDN': mdn.MDN, 'mdn_loss_func': mdn.get_mixture_loss_func(yNumParameters, KMIX)})\n",
    "                    model.summary()\n",
    "                    \n",
    "                else:\n",
    "                    model.compile(loss=mdn.get_mixture_loss_func(yNumParameters,KMIX), optimizer=keras.optimizers.Adam(lr=learnRate))\n",
    "                    model.summary()\n",
    "                    \n",
    "                    cp = keras.callbacks.ModelCheckpoint(pathSaveK + name +'.hdf5', \\\n",
    "                         monitor='val_loss', verbose=1, save_best_only=True, save_weights_only=False, mode='auto', period=1)\n",
    "                    csvlogger = tf.keras.callbacks.CSVLogger(pathSaveK + name + '.txt', separator=\",\", append=False)\n",
    "\n",
    "                    time_callback = TimeHistory()\n",
    "                    history = History()\n",
    "                    history = model.fit(x_data, y_data,\n",
    "                        verbose=0,\n",
    "                        epochs=NEPOCH,\n",
    "                        batch_size=32,\n",
    "                        shuffle=True,\n",
    "                        validation_data=(x_cv, y_cv),\n",
    "                        callbacks=[time_callback, cp, keras.callbacks.TerminateOnNaN(), csvlogger])\n",
    "                        \n",
    "                f = open(pathSaveK + name + '.txt', 'r')\n",
    "                lines = f.readlines()                    \n",
    "                loss = []\n",
    "                loss_cv = []                    \n",
    "                for _l in lines[1:]:\n",
    "                    loss.append(float(_l.split(',')[1]))\n",
    "                    loss_cv.append(float(_l.split(',')[2]))\n",
    "                    \n",
    "                y_trainML = model.predict(x_data)\n",
    "                y_testML = model.predict(x_test)\n",
    "                \n",
    "                Dict_Train = {}\n",
    "                Dict_Test = {}\n",
    "                Dict_CV = {}\n",
    "\n",
    "                def generateDict(x_I, yNumParameters):\n",
    "                    Dict_I = {}\n",
    "                    sigNum = int(yNumParameters*(yNumParameters-1)/2 +yNumParameters)\n",
    "                    mus = np.apply_along_axis((lambda a: a[:KMIX*yNumParameters]), 1, x_I)\n",
    "                    sigs = np.apply_along_axis((lambda a: a[KMIX*yNumParameters:KMIX*yNumParameters+KMIX*sigNum]), 1, x_I)\n",
    "                    pis = np.apply_along_axis((lambda a: a[-KMIX:]), 1, x_I)\n",
    "                    for piInd in range(pis.shape[0]):\n",
    "                        pis[piInd,:] = mdn.softmax(pis[piInd,:])\n",
    "                    \n",
    "                    for i in range(KMIX):\n",
    "                        for ii in range(x_I.shape[0]):\n",
    "                            sig_vals = sigs[ii,i*sigNum:(i+1)*sigNum]\n",
    "                            sig_vals_a = np.concatenate([sig_vals, sig_vals[yNumParameters:][::-1]])\n",
    "                            sig_vals_b = np.reshape(sig_vals_a[::-1], [yNumParameters,yNumParameters])\n",
    "                            sig_vals_c = np.tril(sig_vals_b, k=0)\n",
    "                            for r in range(yNumParameters):\n",
    "                                for c in range(yNumParameters):\n",
    "                                    if r==c:\n",
    "                                        sig_vals_c[r,c] = np.exp(sig_vals_c[r,c])\n",
    "                            Dict_I[\"multi_sigma_\" + str(ii) + \"_\" + str(i)] = sig_vals_c\n",
    "                            \n",
    "                        Dict_I[\"multi_mu\" + str(i)] = mus[:,i*yNumParameters:(i+1)*yNumParameters]\n",
    "                        Dict_I[\"multi_pi\" + str(i)] = pis[:,i]\n",
    "                        \n",
    "                    Dict_I[\"log_likelihood\"] = mdn.get_mixture_log_likelihood(y_test,y_testML,KMIX,yNumParameters).eval()\n",
    "                    return Dict_I\n",
    "                \n",
    "                Dict_Train = generateDict(y_trainML, yNumParameters)\n",
    "                Dict_Test = generateDict(y_testML, yNumParameters)\n",
    "        \n",
    "                alphaVal = 0.5\n",
    "                \n",
    "                plotResults_multi(data, Dict_Test, KMIX, id_string)\n",
    "                \n",
    "                plotIndividual = True\n",
    "                if plotIndividual:\n",
    "                    for ind in range(yNumParameters):\n",
    "                        fig = plt.figure(figsize=(14,9), dpi=1000)\n",
    "                        nCols = 3 \n",
    "                        nRows = 2 \n",
    "                        plotCounter = 1        \n",
    "\n",
    "                        ax = fig.add_subplot(nRows,nCols,plotCounter) \n",
    "                        ax.plot(loss, 'r-', label='Training')\n",
    "                        ax.plot(loss,\"b--\", label='Validation')\n",
    "                        plt.title('Loss')\n",
    "                        ax.set_xlabel('Epochs')\n",
    "                        ax.set_ylabel('Negative log-likelihood')\n",
    "                        ax.legend()\n",
    "                        #ax.locator_params(axis='x', nbins=2)\n",
    "                        plotCounter += 1\n",
    "                        ax = fig.add_subplot(nRows,nCols,plotCounter)\n",
    "                        plotCounter += 1\n",
    "                        ax1 = fig.add_subplot(nRows,nCols,plotCounter)\n",
    "                        plotCounter += 2\n",
    "\n",
    "                        strhere = paraVec[ind]\n",
    "                        pR = plotResults(data,ax,ax1, alphaVal,(paraVec[ind] + \", Train\"), y_data[:,ind], Dict_Train, \\\n",
    "                                     ind, Dict_baseline_train, kernel, KMIX, False)\n",
    "\n",
    "                        ax = fig.add_subplot(nRows,nCols,plotCounter)\n",
    "                        plotCounter += 1\n",
    "                        ax1 = fig.add_subplot(nRows,nCols,plotCounter)\n",
    "                        plotCounter += 1\n",
    "                        pRT = plotResults(data,ax,ax1, alphaVal,(paraVec[ind] + \", Test\"), y_test[:,ind], Dict_Test, \\\n",
    "                                ind, Dict_baseline_test, kernel, KMIX, False)\n",
    "                        plt.tight_layout()\n",
    "                        plt.show()\n",
    "                        \n",
    "                #with open(pathSave + \"/resultsDicts/multiVar_noise_dict\" + str(data.noiseLevel) + \".txt\", \"wb\") as fkl:\n",
    "                #    dump(pRT.Dict_Test, fkl)\n",
    "                    \n",
    "                print(\"-------------------------------------------------------------------------------------------------------\")\n",
    "                print(\"-------------------------------------------------------------------------------------------------------\")\n",
    "                print(\"-------------------------------------------------------------------------------------------------------\")\n",
    "                print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "class meanPredictor:\n",
    "    def __init__(self, x_data,y_data,xparaVec, plotMP=False):\n",
    "        Dict_Mean_Predictor = {}\n",
    "        if 'Tprof_4p5Gyr' in xparaVec and np.size(xparaVec)==1:\n",
    "            tProfOnly = True\n",
    "        else:\n",
    "            tProfOnly = False\n",
    "        for i in range(np.size(xparaVec)):\n",
    "            for j in range(y_data.shape[1]):\n",
    "                if tProfOnly or (xparaVec[i] == 'Tprof_4p5Gyr'):\n",
    "                    observable = np.mean(x_data,axis=1)\n",
    "                else:\n",
    "                    observable = x_data[:,i] \n",
    "                para = y_data[:,j]\n",
    "                steps = 10\n",
    "\n",
    "                bin_means, bin_edges, binnumber = binned_statistic(observable, para, statistic='mean',bins=steps)\n",
    "                bin_variance = [np.std(para[np.where(binnumber == bIndex)[0]]) for bIndex in range(1,steps+1)]\n",
    "                removeVec = []\n",
    "                notRemoveVec = []\n",
    "                for i3 in range(np.size(bin_means)):\n",
    "                    if np.isnan(bin_means[i3]) or bin_variance[i3] == 0.0:\n",
    "                        removeVec.append(i3)\n",
    "                    else:\n",
    "                        notRemoveVec.append(i3)\n",
    "\n",
    "                for r in removeVec:\n",
    "                    bin_variance[r] = np.average([bin_variance[notInd] for notInd in notRemoveVec])\n",
    "                    bin_means[r] = np.average([bin_means[notInd] for notInd in notRemoveVec])\n",
    "\n",
    "                digitizer = np.digitize(observable,bin_edges) - 1\n",
    "                mu = np.zeros(np.size(digitizer))\n",
    "                variance = np.zeros(np.size(digitizer))\n",
    "                for ind in range(np.size(digitizer)):\n",
    "                    mu[ind] = bin_means[digitizer[ind]-1]\n",
    "                    variance[ind] = bin_variance[digitizer[ind]-1]\n",
    "                Dict_Mean_Predictor[\"mu_var\" + str(i) + \"_\" + str(j)] = mu\n",
    "                Dict_Mean_Predictor[\"variance_var\" + str(i) + \"_\" + str(j)] = variance\n",
    "                if plotMP:\n",
    "                    plt.figure(figsize=(4,6))\n",
    "                    plt.plot(observable,para, \"ro\", alpha = 0.4)\n",
    "                    plt.plot((bin_edges[1:]+bin_edges[0:-1])/2.,bin_means, \"k-\")\n",
    "                    plt.plot((bin_edges[1:]+bin_edges[0:-1])/2.,bin_means+bin_variance,\"b--\")\n",
    "                    plt.plot((bin_edges[1:]+bin_edges[0:-1])/2.,bin_means-bin_variance,\"b--\")\n",
    "                    plt.legend([_, \"$\\mu$\", \"$\\mu + \\sigma$\", \"$\\mu - \\sigma$\"], loc=\"lower_right\")\n",
    "        self.Dict_Mean_Predictor = Dict_Mean_Predictor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bigLoop(_zipped, repeats=1, trainOrload='load'):\n",
    "    Dict_kls = {}\n",
    "    hSize, KMIX, x_o, x_b, trainPercent, yIndex, kernel, activation, noiseLevel = _zipped\n",
    "    \n",
    "    fontSize = 14\n",
    "    picID = 0 #yIndices  [yIndices[i5]]\n",
    "    data = getData([yIndex], x_o, x_b, trainPercent, pathSave, noiseLevel)\n",
    "\n",
    "    id_string = str(_zipped)\n",
    "\n",
    "    with open(pathSave + \"/Data_Files/processedData\" + id_string + \".txt\", \"rb\") as fkl:\n",
    "        dataDict = load(fkl)\n",
    "\n",
    "    data.x_data = dataDict['x_data']\n",
    "    data.x_test = dataDict['x_test']\n",
    "    data.x_cv = dataDict['x_cv']     \n",
    "    data.y_data = dataDict['y_data']\n",
    "    data.y_test = dataDict['y_test']\n",
    "    data.y_cv = dataDict['y_cv']     \n",
    "    data.paraVec = dataDict['paraVec']\n",
    "    data.xparaVec = dataDict['xparaVec']\n",
    "    data.rProf = dataDict['rProf']\n",
    "    data.startPoint = dataDict['startPoint']\n",
    "    data.endPoint = dataDict['endPoint']    \n",
    "    data.indexer = dataDict['indexer']    \n",
    "    data.pMax = dataDict['pMax']  \n",
    "    data.pMin = dataDict['pMin']  \n",
    "    data.oMax = dataDict['oMax']\n",
    "    data.oMin = dataDict['oMin']\n",
    "\n",
    "    Dict_MP_Test = meanPredictor(data.x_test, data.y_test, data.xparaVec, False)\n",
    "    Dict_MP_Train = meanPredictor(data.x_data, data.y_data, data.xparaVec, False)\n",
    "\n",
    "    #data.visualizeData(data, picID, fontSize, Dict_MP_Train.Dict_Mean_Predictor)\n",
    "\n",
    "    MDN_K(data,data.x_data,data.y_data, data.x_test,data.y_test, data.x_cv, \\\n",
    "          data.y_cv, hSize, KMIX, 50000, 0.00001, \\\n",
    "          data.paraVec, data.xparaVec, Dict_MP_Train.Dict_Mean_Predictor, \\\n",
    "          Dict_MP_Test.Dict_Mean_Predictor, picID, repeats, id_string, kernel, activation, trainOrload)\n",
    "    picID += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "zipped_new = [[[60, 60, 60], 3, [0,1,2,3,4,5,6], [0, 100], 0.9, [0, 4, 5, 6, 7], 'gaussian', 'tanh', 0.0]]\n",
    "\n",
    "runParallel = False\n",
    "trainOrload = 'load'\n",
    "repeats = 1\n",
    "\n",
    "if runParallel:                   \n",
    "    Parallel(n_jobs=1, verbose=10, backend='loky', prefer='processes')(delayed(bigLoop)(_z, repeats, trainOrload) for _z in zipped_new)\n",
    "else:\n",
    "    [bigLoop(_z, repeats, trainOrload) for _z in zipped_new] "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}