run.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9f1f5868",
   "metadata": {},
   "source": [
    "Source code for Can machine learning algorithms improve upon classical palaeoenvironmental reconstruction models?\n",
    "S.P. Leiden, Netherlands."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd13ec86",
   "metadata": {},
   "source": [
    "# Main functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "180045d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.autograd import Variable\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "import torch.optim as optim\n",
    "import random\n",
    "from uq360.algorithms.blackbox_metamodel import MetamodelRegression\n",
    "from tqdm import tqdm\n",
    "from IPython.display import clear_output\n",
    "new_path = 'xxx'\n",
    "\n",
    "# Generate the communtity list\n",
    "\n",
    "def make_sents(sp_matrix):\n",
    "    sent = []\n",
    "    unfort_list = []\n",
    "    ii =0\n",
    "    for i in range(sp_matrix.shape[0]):  \n",
    "\n",
    "        pd_table_sub = sp_matrix.iloc[i,:]>0\n",
    "        c_v = sp_matrix.iloc[i,:][sp_matrix.iloc[i,:]>0]\n",
    "        if len(c_v)>1:\n",
    "            sent.append(c_v)\n",
    "        else:\n",
    "            unfort_list.append(ii)\n",
    "        ii += 1\n",
    "        print('\\r  %s...'%ii,end = '')\n",
    "    return(sent)\n",
    "\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "\n",
    "# The main function for ecological baesd glove\n",
    "def glove_dic(sent,embed_size,vocab,wf,num_epochs = 5,glove_abuna = True):\n",
    "    # Set parameters\n",
    "    xmax = 2\n",
    "    alpha = 0.75\n",
    "    batch_size = 10000\n",
    "    l_rate = 0.001\n",
    "    num_epochs = num_epochs\n",
    "\n",
    "    # Create vocabulary and word lists\n",
    "    word_list = sent\n",
    "    w_list_size = len(word_list)\n",
    "    vocab_size = len(vocab)\n",
    "\n",
    "    # Create word to index mapping\n",
    "    w_to_i = {word: ind for ind, word in enumerate(vocab)}\n",
    "    comat = np.zeros((vocab_size, vocab_size))\n",
    "    if glove_abuna == True:\n",
    "        # Construct co-occurence matrix\n",
    "        for i_init in range(w_list_size):\n",
    "            cc = sent[i_init]\n",
    "            for i in range(len(cc)):\n",
    "                for j in range(len(cc)):\n",
    "                    if i != j:\n",
    "                        ind = w_to_i[cc.index[i]]\n",
    "                        lind = w_to_i[cc.index[j]]\n",
    "                        comat[ind, lind] += cc.values[i]\n",
    "                        comat[lind,ind] += cc.values[j]\n",
    "    # Non-zero co-occurrences\n",
    "    else:\n",
    "        for i_init in range(w_list_size):\n",
    "            cc = sent[i_init]\n",
    "            for i in range(len(cc)):\n",
    "                for j in range(len(cc)):\n",
    "                    if i != j:\n",
    "                        ind = w_to_i[cc.index[i]]\n",
    "                        lind = w_to_i[cc.index[j]]\n",
    "                        comat[ind, lind] += 1\n",
    "                        comat[lind,ind] += 1\n",
    "\n",
    "    coocs = np.transpose(np.nonzero(comat))\n",
    "    # Weight function\n",
    "    # Set up word vectors and biases\n",
    "    import torch\n",
    "    import torch.optim as optim\n",
    "\n",
    "    # Your existing code...\n",
    "    matrix_size = comat.shape[0]\n",
    "    matrix_size = (comat.shape[0], comat.shape[0])\n",
    "    matrix_size2 = (comat.shape[0],256)\n",
    "\n",
    "    matrix_with_embedding = torch.randn(matrix_size2) * 0.01\n",
    "    #bias_matrix = torch.randn(matrix_size) * 0.01\n",
    "    matrix_with_embedding2 = matrix_with_embedding.clone()\n",
    "    bias_r = torch.randn(comat.shape[0]) * 0.01\n",
    "    bias_l = torch.randn(comat.shape[0]) * 0.01\n",
    "    bias_r.requires_grad = True\n",
    "    bias_l.requires_grad = True\n",
    "    bias_l.to(device)\n",
    "    bias_r.to(device)\n",
    "\n",
    "    matrix_with_embedding.requires_grad = True\n",
    "    matrix_with_embedding2.requires_grad = True\n",
    "    #bias_matrix.requires_grad = True\n",
    "    matrix_with_embedding.to(device)\n",
    "    matrix_with_embedding2.to(device)\n",
    "\n",
    "    #bias_matrix.to(device)\n",
    "    l_rate = 0.001  # Reduced learning rate\n",
    "    coocs_matrix = torch.tensor(comat.copy()).to(torch.float64)\n",
    "\n",
    "    optimizer = optim.AdamW([matrix_with_embedding,matrix_with_embedding2, bias_l,bias_r], lr=l_rate)  # Specify parameters for optimization\n",
    "    add_log = -torch.log(coocs_matrix.detach()) \n",
    "    add_log[torch.isinf(add_log)] = -1e10\n",
    "    for epoch in range(num_epochs):\n",
    "        optimizer.zero_grad()\n",
    "        loss = torch.nansum(torch.mul((torch.matmul(matrix_with_embedding, matrix_with_embedding2.T) + torch.add(bias_l, bias_r.T) + add_log)**2, wf_m(coocs_matrix)))\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        print(\"\\r Average loss for epoch \" + str(epoch + 1) + \": \", loss.item(), end='')\n",
    "        print(\"                                                \", end='')\n",
    "\n",
    "    print(' ')\n",
    "\n",
    "    print(' ')\n",
    "    lem_list =[];rem_list=[];spname_list = []\n",
    "    l_embed = matrix_with_embedding#[:,:,0].Ｔ\n",
    "\n",
    "    for word_ind in range(len(vocab)):\n",
    "        # Create embedding by summing left and right embeddings\n",
    "        rem_list.append(l_embed[word_ind].data.numpy())  \n",
    "        spname_list.append(vocab[word_ind])\n",
    "    #cc = pd.DataFrame(np.array(lem_list)[:,:,0])\n",
    "    #cc['Species'] = spname_list\n",
    "    emb_dic = pd.DataFrame(np.array(rem_list)) #[:,:,0]\n",
    "    #emb_dic['Species'] = spname_list\n",
    "    emb_dic.index = vocab\n",
    "    \n",
    "    return(emb_dic)\n",
    "\n",
    "\n",
    "\n",
    "# Generate the communtity embedding\n",
    "\n",
    "def make_com_emb(sp_matrix,emb_dic,embed_size):\n",
    "    emb_sent = []\n",
    "    for i in range(sp_matrix.shape[0]):  \n",
    "\n",
    "        pd_table_sub = sp_matrix.iloc[i,:]#>0\n",
    "        c_v = pd_table_sub[pd_table_sub>0].index#np.unique(pd_table_sub.Species.dropna())\n",
    "        init_vec = np.zeros([embed_size])\n",
    "        for ii in c_v:\n",
    "            #i_sum += 1\n",
    "            init_vec += np.float32(emb_dic.loc[ii].values)\n",
    "\n",
    "        emb_sent.append(init_vec/len(c_v)) #/i_sum\n",
    "\n",
    "    #print('\\r  %s.../95104'%i,end = '')\n",
    "    return(np.array(emb_sent))\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from lightgbm import LGBMRegressor\n",
    "import lightgbm as lgb\n",
    "\n",
    "from sklearn.ensemble import RandomForestRegressor,ExtraTreesRegressor,GradientBoostingRegressor\n",
    "from sklearn.linear_model import LinearRegression,ARDRegression\n",
    "\n",
    "import argparse\n",
    "import pprint\n",
    "#import gensim\n",
    "from sklearn.model_selection import cross_validate,LeaveOneOut\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import GridSearchCV,train_test_split\n",
    "import pickle\n",
    "import math\n",
    "#import random\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error,accuracy_score,recall_score\n",
    "import random\n",
    "def to_mae(y_true, y_pred):\n",
    "    return mean_absolute_error(y_true, y_pred)\n",
    "def to_rmse(y_true, y_pred):\n",
    "    return np.sqrt(mean_squared_error(y_true, y_pred))\n",
    "def to_mse(y_true, y_pred):\n",
    "    return (mean_squared_error(y_true, y_pred))  \n",
    "# The embedding name for printing importantce indicators\n",
    "def emb_names(table):\n",
    "    new_p_list = []\n",
    "    for i in range(table.shape[1]):\n",
    "        new_p_list.append('em'+str(i))\n",
    "    return(new_p_list)\n",
    "\n",
    "# Keep the fossile data have the same taxon with normal data set\n",
    "def trans_to_core(spec,core):\n",
    "    new_pd = pd.DataFrame(np.zeros([core.shape[0], spec.shape[1]]))\n",
    "    new_pd.columns = spec.columns\n",
    "    for i_index,i_names in enumerate(new_pd.columns):\n",
    "        try:            \n",
    "            new_pd.iloc[:,i_index] = core[i_names].values\n",
    "        except:\n",
    "            continue\n",
    "    return(new_pd)\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "def to_pers(df):\n",
    "    row_sums = df.sum(axis=1)\n",
    "    # Divide each value in the row by the corresponding row sum and multiply by 100\n",
    "    df_percentage = df.div(row_sums, axis=0) \n",
    "    return(df_percentage)\n",
    "\n",
    "def gcd_many(s):\n",
    "    g = 0\n",
    "    for i in range(len(s)):\n",
    "        if i == 0:\n",
    "            g = s[i]\n",
    "        else:\n",
    "            g=math.gcd(g,s[i])\n",
    "    return g\n",
    "\n",
    "\n",
    "def to_num(matrix_init):\n",
    "    matrix = matrix_init.copy()\n",
    "    for i in range(matrix.shape[0]):\n",
    "        sub_raw = matrix.iloc[i,:].values\n",
    "        sub_raw_2 = sub_raw[sub_raw>0]\n",
    "        min_raw  = min(sub_raw_2)\n",
    "        sub_raw = sub_raw*10000#sub_raw/min_raw  ## 0/min_raw  = 0\n",
    "        sub_raw = sub_raw/gcd_many(sub_raw[sub_raw>0].astype('int'))\n",
    "    matrix.iloc[i,:] = sub_raw\n",
    "    return(matrix)\n",
    "\n",
    "\n",
    "# The main function for MEMLM\n",
    "def pre_glove(Biological_abundance_matrix,environmental_values_matrix,test_method = 'loo',no_components = 256, data_name='Test',\n",
    "              add_tr_env = [], cores = [],core_sp = [],stack = True, glove_abuna = True,random_no = 123,\n",
    "              recon=False,pass_va = False,save_vali = True,cor_emb = False, emb_method = 'plus',test_set = [],num_epochs = 300):\n",
    "    \n",
    "    # set the flags\n",
    "    if len(test_set) >0:\n",
    "        recon= False\n",
    "        pass_va = False\n",
    "        cores = []\n",
    "        \n",
    "    if emb_method == 'only':\n",
    "        only_emb = True\n",
    "        only_abun = False\n",
    "        print('using the embedding matrix')\n",
    "    elif emb_method == 'none':\n",
    "        only_abun = True\n",
    "        only_emb = False\n",
    "        print('using the abundance matrix')\n",
    "    else:\n",
    "        only_abun = False\n",
    "        only_emb = False\n",
    "        print('using the embedding + abundance matrix')\n",
    "    # Set the random seeds\n",
    "    random.seed(random_no)\n",
    "    torch.manual_seed(random_no)\n",
    "    torch.cuda.manual_seed(random_no)  # If using CUDA\n",
    "    torch.manual_seed(random_no)\n",
    "    # Set the random seed for NumPy\n",
    "    np.random.seed(random_no)   \n",
    "        \n",
    "        \n",
    "        \n",
    "    pd_table = Biological_abundance_matrix\n",
    "    envs = environmental_values_matrix\n",
    "    #save_vali = save_vali            \n",
    "    if recon ==False:\n",
    "        cores = []\n",
    "    else:\n",
    "        #print(cores.shape)\n",
    "        cor_sp = (set(cores.columns)&set(pd_table.columns))\n",
    "    try:\n",
    "        (envs.shape[1])\n",
    "    except:\n",
    "        envs = pd.DataFrame(envs.T)\n",
    "    \n",
    "    #Statistics of species contained in each sample\n",
    "    print('scale: %s'%pd_table.iloc[0,:].sum())\n",
    "    xmax = 2 \n",
    "    alpha = 0.75\n",
    "    def wf(x):\n",
    "        return x\n",
    "    if glove_abuna == True:\n",
    "        if (pd_table.iloc[0,:].sum() > 1.3):\n",
    "            #pd_table_2 = pd_table\n",
    "            #if np.min(pd_table_2[pd_table_2>0]).min()>1:\n",
    "            def wf(x):\n",
    "                if x < xmax:\n",
    "                    return (x/xmax)**alpha\n",
    "                return 1\n",
    "        #else:\n",
    "            #pd_table_2 = to_num(pd_table)\n",
    "\n",
    "    #pd_table_sub = pd_table.iloc[i,:]*250\n",
    "    if recon == True:\n",
    "        if cor_emb == True:\n",
    "            pd_table_2 = pd_table_2[cor_sp]\n",
    "    if only_abun == False:\n",
    "        sent = make_sents(pd_table)\n",
    "        print('creating the embedding now...')   ###\n",
    "    \n",
    "        embed_size = no_components\n",
    "        vocab = pd_table.columns\n",
    "        #embed_size = 256\n",
    "        emb_dic = glove_dic(sent,embed_size,vocab,wf,glove_abuna = glove_abuna,num_epochs = num_epochs)\n",
    "        print(pd_table.shape,embed_size)\n",
    "        paragraph_list = make_com_emb(pd_table,emb_dic,embed_size)\n",
    "        emb_name_list = list(emb_names(paragraph_list))\n",
    "\n",
    "        #print('Species num(availble in asscoiation): %s' % len(corpus_model.dictionary))\n",
    "    print('Species num: %s' % pd_table.shape[1])\n",
    "    print('Samples num: %s' % pd_table.shape[0])\n",
    "    print(' ')\n",
    "    print('=================================')\n",
    "    print('training the embedding ......')\n",
    "    \n",
    "\n",
    "\n",
    "    \n",
    "    if only_emb == True:  \n",
    "        para_names =  emb_name_list   \n",
    "\n",
    "    elif only_abun == True:\n",
    "        para_names = pd_table.columns\n",
    "        paragraph_list = np.array(pd_table)\n",
    "    else:\n",
    "        print('adding the speiese abunadance now...')\n",
    "        paragraph_list = np.concatenate((np.array(paragraph_list),np.array(pd_table)), axis = 1)    ###\n",
    "        #data_name = '%s_rich'%data_name\n",
    "        para_names =  emb_name_list + list(pd_table.columns)\n",
    "        #print([len(para_names),len(list(emb_names(np.array(paragraph_list)))),len(list(emb_names(np.array(paragraph_list))))])    \n",
    "        #print(paragraph_list.shape)\n",
    "        \n",
    "    if  len(add_tr_env) !=0:\n",
    "        data_name = '%s_added'%data_name\n",
    "        paragraph_list = np.concatenate((np.array(paragraph_list),np.array(add_tr_env)), axis = 1)\n",
    "        para_names = emb_name_list + list(add_tr_env.columns)\n",
    "        #print(paragraph_list.shape)\n",
    "    #print(len(para_names))\n",
    "\n",
    "        #print(paragraph_list.shape,cores.shape)\n",
    "    if pass_va == False :  \n",
    "        if test_method == 'cv':\n",
    "            cs=train_data_pre(paragraph_list,envs,data_name=data_name,test_method='cv',para_names = para_names, stack = stack, save_vali = save_vali)\n",
    "        elif test_method == 'loo':\n",
    "            cs=train_data_pre(paragraph_list,envs,data_name=data_name,test_method='loo',para_names = para_names, stack = stack, save_vali = save_vali)\n",
    "        elif test_method == 'diy':\n",
    "            cs=train_data_pre(paragraph_list,envs,data_name=data_name,test_method='diy',para_names = para_names, stack = stack,test_set=test_set)\n",
    "\n",
    "        else:\n",
    "            print('you shall choose one of these validate method:  cv, loo,diy') \n",
    "    if len(cores) !=0:\n",
    "        print('reconstructing the %s'%data_name)\n",
    "        cores_index = cores.index\n",
    "        if len(cores.columns) > len(pd_table.columns):      \n",
    "            cores = cores[cor_sp]\n",
    "            \n",
    "        elif len(cores.columns) < len(pd_table.columns):\n",
    "            cores = trans_to_core(pd_table,cores)\n",
    "        paragraph_list2 = []\n",
    "        #Generate embedding for each sample\n",
    "        cores.columns = pd_table.columns\n",
    "        if only_abun == False:\n",
    "            paragraph_list2 = make_com_emb(cores,emb_dic,embed_size)\n",
    "            if only_emb == False:\n",
    "                paragraph_list = np.concatenate((np.array(paragraph_list)[:,:no_components],np.array(pd_table)), axis = 1)\n",
    "                para_names = emb_name_list+list(pd_table.columns)\n",
    "                cores = np.concatenate((np.array(paragraph_list2),np.array(cores)), axis = 1)    #\n",
    "            else:\n",
    "                cores = np.array(paragraph_list2)\n",
    "        cs = recon_data(paragraph_list,envs,cores,cores_index,pass_va,data_name=data_name,trastacks = stack,test_method=test_method,para_names = para_names,random_no=random_no)\n",
    "        \n",
    "    print('The operation is over. Have a nice day!')\n",
    "    print('')\n",
    "    print('')\n",
    "\n",
    "    return(cs)\n",
    "\n",
    "def train_data_pre(X,y,data_name,GBR = False,n_estimators = 1000,test_method = 'loo',stack = True,save_vali = False,para_names=[],test_set = []):\n",
    "    if save_vali == True:\n",
    "        op_pirnt = 'on'\n",
    "    else:\n",
    "        op_pirnt = 'off'\n",
    "    print('save validation results option: %s'%op_pirnt)\n",
    "    import warnings\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    acc_tt_list = []\n",
    "    try:\n",
    "        y.shape[1]\n",
    "        init_name = data_name\n",
    "        print('has multi environmental values')\n",
    "        for y_col in range(y.shape[1]):\n",
    "            pd_acc = pd.DataFrame(np.zeros([y.shape[1]*2,5]))\n",
    "            data_name = '%s_%s'%(init_name,y.columns[y_col])\n",
    "            X =np.array(X)\n",
    "            yy =np.array(y.iloc[:,y_col])\n",
    "            xx= X[~pd.isnull(yy)]\n",
    "            print('sub mission:  training %s now...'%data_name)\n",
    "            if test_method== 'loo':\n",
    "                acc = train_data_loo(X=xx,y = yy,GBR = False,data_name=data_name,n_estimators = 1000,stack = stack, save_vali = save_vali,test_method=test_method,para_names = para_names)\n",
    "            elif test_method== 'cv':\n",
    "                acc = train_data_cv(X=xx,y = yy,GBR = False,data_name=data_name,n_estimators = 1000,stack = stack, save_vali = save_vali,test_method=test_method,para_names = para_names)\n",
    "            elif test_method== 'diy':\n",
    "                return(train_data_diy(X=xx,y = yy,test_set = test_set, GBR = False,data_name=data_name,n_estimators = 1000,stack = stack, test_method=test_method,para_names = para_names))\n",
    "\n",
    "        acc_tt_list.append([data_name,*acc])  \n",
    "        return(acc_tt_list)\n",
    "    except:\n",
    "        if test_method == 'split':\n",
    "            return(train_data_split_across(X,y,GBR = False,data_name=data_name,n_estimators = 1000,stack = stack,test_method=test_method,para_names = para_names))\n",
    "        if test_method == 'loo':\n",
    "            return(train_data_loo(X,y,GBR = False,data_name = data_name,n_estimators = 1000,stack = stack, save_vali = save_vali,test_method=test_method,para_names = para_names,random_no=random_no))\n",
    "        elif test_method == 'cv':\n",
    "            return(train_data_cv(X,y,GBR = False,data_name=data_name,n_estimators = 1000,stack = stack, save_vali = save_vali,test_method=test_method,para_names = para_names,random_no=random_no))     \n",
    "        elif test_method== 'diy':\n",
    "            return(train_data_diy(X=xx,y = yy,test_set = test_set, GBR = False,data_name=data_name,n_estimators = 1000,stack = stack, test_method=test_method,para_names = para_names)\n",
    ")\n",
    "def train_data_cv(X,y,data_name,test_method = [],para_names = [],GBR = False,n_estimators = 1000,stack = True,save_vali = False,random_no=0):       \n",
    "    \n",
    "    base_list = np.array(range(X.shape[0]))\n",
    "    x_train_names = para_names\n",
    "    ran_index = np.random.permutation(base_list)\n",
    "    y= np.array(y)\n",
    "    X = np.array(X)\n",
    "    X = X[ran_index]\n",
    "    y =y[ran_index]\n",
    "    y_token = 'multi'\n",
    "    y_test_l = []\n",
    "    pre_l = []\n",
    "    kf = KFold(n_splits=2)\n",
    "    print(ran_index[:20],'ran base')\n",
    "    for train, test in kf.split(base_list):\n",
    "    \n",
    "        try:\n",
    "            x_train, x_test, y_train, y_test = X[train],X[test],y[train],y[test]\n",
    "            #print(train[:20],'train')\n",
    "\n",
    "        except:\n",
    "            x_train, x_test, y_train, y_test = X.loc(axis=0)[train],X.loc(axis=0)[test],y.loc(axis=0)[train],y.loc(axis=0)[test]        \n",
    "        #print(x_train.shape,y_train.shape)\n",
    "        predict3 = make_predict(x_train, x_test, y_train, y_test,test_method,data_name,x_train_names,random_no=random_no)\n",
    "        if stack == True:\n",
    "\n",
    "            predict3 = to_stack_2(x_train,y_train,predict3,random_no=random_no)\n",
    "\n",
    "        y_test_l.append(y_test.reshape(-1,1))\n",
    "        pre_l.append(predict3)\n",
    "    \n",
    "    pd_re = tail(pre_l,y_test_l,predict3,data_name,stack = stack,test_method ='cv',save_vali = save_vali)\n",
    "    return(pd_re) #[choose_predict(acc)] #,y_test, predict3,rfr\n",
    "\n",
    "\n",
    "\n",
    "def train_data_diy(X,y,test_set,data_name,test_method = [],para_names = [],GBR = False,n_estimators = 1000,stack = True):       \n",
    "\n",
    "    x_train_names = para_names\n",
    "    X =np.array(X)\n",
    "    y = np.array(y)\n",
    "    y_token = 'multi'\n",
    "    print('test',X.shape)\n",
    "\n",
    "    \n",
    "    x_train = X\n",
    "    x_test = y\n",
    "    y_test = []\n",
    "    y_train = test_set\n",
    "    #print('MEASURED     RECON       RICHNESS')\n",
    "\n",
    "    predict3 = make_predict(x_train, x_test, y_train, y_test,test_method,data_name,x_train_names,print_importa = False)\n",
    "    if stack == True:\n",
    "        #print(y_train.shape,x_train.shape,np.array(predict3).shape)\n",
    "\n",
    "        predict3 = to_stack_2(x_train,y_train,predict3)\n",
    "        #print('OK')\n",
    "    return(predict3) #[choose_predict(acc)] #,y_test, predict3,rfr\n",
    "\n",
    "def make_predict(x_train, x_test, y_train, y_test,data_name , test_method,x_train_names,n_estimators = 1000, GBR = False, print_importa = True,random_no = 0):\n",
    "    \n",
    "    if len(x_train.shape) == 1:\n",
    "            x_train = x_train.reshape(-1, 1)\n",
    "            x_test = x_test.reshape(-1, 1)\n",
    "    if len(y_test)>0:\n",
    "        if len(y_test.shape) == 1:\n",
    "            y_train = y_train.reshape(-1, 1)\n",
    "            y_test = y_test.reshape(-1, 1)\n",
    "            y_token = 'single'\n",
    "    #print([x_train.shape,x_test.shape])\n",
    "    predict3 = three_models( x_train,y_train,x_test,test_method,data_name, x_train_names, print_importa = True,print_model = False,random_no = random_no)\n",
    "\n",
    "\n",
    "    return(predict3)\n",
    "        \n",
    "def to_stack(x_train,y_train,predict3):\n",
    "    rfr2,best_no, quit_lightGBM = mulit_model_stack2(x_train,y_train)\n",
    "    #print(rfr2,best_no, quit_lightGBM)\n",
    "    acc3_l = []\n",
    "    acc3_l.append(np.array(predict3)[0,:])\n",
    "    acc3_l.append(np.array(predict3)[1,:])\n",
    "    acc3_l.append(np.array(predict3)[2,:])\n",
    "    #print(acc3_l)\n",
    "    #print(predict3)\n",
    "    if quit_lightGBM == True:\n",
    "        predict3.append(rfr2.predict((np.array(acc3_l).T)[:,:1]).reshape(-1,1)\n",
    "                       )\n",
    "    else:\n",
    "        predict3.append(rfr2.predict(np.array(acc3_l).T).reshape(-1,1))\n",
    "    #print(predict3)\n",
    "    return(predict3)\n",
    "\n",
    "def to_stack_2(x_train,y_train,predict3,random_no=0):\n",
    "    rfr2,best_no, quit_lightGBM = mulit_model_stack2(x_train,y_train,test_method = 'cv',random_no=random_no)\n",
    "    predict3 = np.array(predict3).T   # (3,8)\n",
    "    acc3_l = []\n",
    "    acc3_l.append(predict3[:,0])\n",
    "    acc3_l.append(predict3[:,1])\n",
    "    acc3_l.append(predict3[:,2])\n",
    "\n",
    "    if quit_lightGBM == True:\n",
    "                predict3 = np.concatenate((predict3,rfr2.predict(np.array(acc3_l).T)[:,:1].reshape(-1,1)\n",
    "                                          ), axis = 1) \n",
    "    else:\n",
    "        predict3 = np.concatenate((predict3,rfr2.predict(np.array(acc3_l).T).reshape(-1,1)), axis = 1) \n",
    "    #print(predict3)\n",
    "    return(predict3)\n",
    "    \n",
    "    \n",
    "def train_data_loo(X,y,test_method, data_name,para_names,GBR = False,n_estimators = 1000, stack = False,save_vali = False):\n",
    "    y = np.array(y)\n",
    "    x_train_names = para_names\n",
    "    X =np.array(X)\n",
    "    y_token = 'multi'\n",
    "    X_no =range(X.shape[0])\n",
    "    loo = LeaveOneOut()\n",
    "    y_test_l = []\n",
    "    pre_l = []\n",
    "    print('MEASURED     RECON       RICHNESS')\n",
    "\n",
    "    # Loo cross_validate\n",
    "    for train, test in loo.split(X_no):\n",
    "\n",
    "        x_train, x_test, y_train, y_test = X[train],X[test],y[train],y[test]\n",
    "\n",
    "        #x_train, x_test, y_train, y_test = train_test_split( X,y, test_size=0.50, random_state=33)\n",
    "        \n",
    "        predict3 = make_predict(x_train, x_test, y_train, y_test,test_method,data_name,x_train_names)\n",
    "\n",
    "        if stack == True:\n",
    "            predict3 = to_stack(x_train,y_train,predict3)\n",
    "\n",
    "        y_test_l.append(y_test)\n",
    "        pre_l.append(predict3)\n",
    "        print(\"   %.3f      %.3f         %s\" %(pre_l[-1][0],y_test_l[-1][0],len(x_test[x_test>0])))\n",
    "    #acc = []\n",
    "    pd_re = tail(pre_l,y_test_l,predict3,data_name,stack = stack,test_method = 'loo',save_vali = save_vali)\n",
    "    \n",
    "    return(pd_re) #[choose_predict(acc)] #,y_test, predict3,rfr\n",
    "\n",
    "def tail(pre_l,y_test_l,predict3,data_name,stack = False,test_method = 'loo',save_vali = False):\n",
    "    if test_method == 'cv':\n",
    "        predict3 = np.vstack(pre_l)\n",
    "        #print(predict3.shape)\n",
    "    else:\n",
    "        predict3 = np.array(pre_l)\n",
    "    y_test =  np.vstack(y_test_l).reshape(-1, 1)\n",
    "\n",
    "    acc_r2 = []\n",
    "    acc_rmse = []\n",
    "    #print(y_test.shape,predict3.shape)\n",
    "    #if test_method == 'cv':\n",
    "    #predict3 = predict3.T\n",
    "    for i in range(predict3.shape[1]):\n",
    "\n",
    "        try:\n",
    "            acc_r2.append(r2_score(y_test, predict3[:,i].reshape(-1, 1)))\n",
    "            acc_rmse.append(to_rmse(y_test, predict3[:,i].reshape(-1, 1)))\n",
    "        except:\n",
    "            acc_r2.append(r2_score(y_test, predict3[:,i]))\n",
    "            acc_rmse.append(to_rmse(y_test, predict3[:,i]))\n",
    "    try:\n",
    "        acc_r2.append(r2_score(y_test, predict3.mean(1).reshape(-1, 1)))\n",
    "        acc_rmse.append(to_rmse(y_test, predict3.mean(1).reshape(-1, 1)))\n",
    "    except:\n",
    "        acc_r2.append(r2_score(y_test, predict3[:,i].mean(1)))\n",
    "        acc_rmse.append(to_rmse(y_test, predict3[:,i].mean(1)))\n",
    "    acc = [acc_r2,acc_rmse]\n",
    "      \n",
    "    print(\"___________result___________\")\n",
    "    print(\"r2_score\")\n",
    "    if stack == False:\n",
    "        print(\"   rfr       etr          lightGBM\")\n",
    "        print(\"   %.3f     %.3f        %.3f\" %(acc[0][0],acc[0][1],acc[0][2]))\n",
    "        print(\"RMSEP\")\n",
    "\n",
    "        print(\"   rfr       etr          lightGBM\")\n",
    "        print(\"   %.3f     %.3f        %.3f\" %(acc[1][0],acc[1][1],acc[1][2]))\n",
    "        va_cloumns = [\"rfr\"  ,     \"etr\",          \"lightGBM\",  \"true\"]\n",
    "    else:\n",
    "\n",
    "        print(\"   rfr       etr       lightGBM       stack       mean\")\n",
    "        print(\"   %.3f     %.3f        %.3f       %.3f       %.3f\" %(acc[0][0],acc[0][1],acc[0][2],acc[0][3],acc[0][4]))\n",
    "        print(\"RMSEP\")\n",
    "\n",
    "        print(\"   rfr       etr       lightGBM       stack       mean\")\n",
    "        print(\"   %.3f     %.3f        %.3f       %.3f       %.3f\" %(acc[1][0],acc[1][1],acc[1][2],acc[1][3],acc[1][4]))\n",
    "        va_cloumns = [\"rfr\"  ,     \"etr\",          \"lightGBM\",   \"stack\",'mean',\"true\"]\n",
    "    ##pickle.dump([predict3,y_test],open('./pkl/cross_val_%s.pkl'%(data_name),\"wb\"), protocol=3)\n",
    "    #print(predict3.shape,y_test.shape)\n",
    "\n",
    "\n",
    "    if save_vali == True:\n",
    "        try:\n",
    "            va_results = pd.DataFrame(np.concatenate((predict3[:,:],y_test),axis = 1))\n",
    "        except:\n",
    "            va_results = pd.DataFrame(np.concatenate((predict3[:,:,0],y_test),axis = 1))\n",
    "        va_results.columns = va_cloumns[:-1]\n",
    "        va_results.to_csv('vadata_%s_%s.csv'%(test_method,data_name))\n",
    "    #pickle.dump([predict3,y_test],open('./pkl/vali_%s_%s.pkl'%(test_method,data_name),\"wb\"), protocol=3)\n",
    "    #r2_score(y_test.reshape((-1, 1)), predict3[-1].reshape((-1, 1)))\n",
    "    #mo.score(y_test.reshape((-1, 1)), predict3[0]\n",
    "    model_num = predict3.shape[1]\n",
    "    pd_re = pd.DataFrame(acc)\n",
    "    pd_re.columns =  ['rfr','etr','lightGBM','stack','mean']#[:model_num]\n",
    "    pd_re.index = ['r2_score','rmsep']\n",
    "    pd_re.to_csv('varesult_%s_%s.csv'%(test_method,data_name))   ######\n",
    "    return(pd_re)\n",
    "\n",
    "\n",
    "\n",
    "def recon_data(paragraph_list,envs,cores,cores_index,pass_va,para_names,test_method, data_name,GBR =False,trastacks = True,random_no=0):\n",
    "    ## starting recon data process!\n",
    "    y_train = envs\n",
    "    x_train_names = para_names\n",
    "    x_train = np.array(paragraph_list)\n",
    "    x_test = cores\n",
    "    predict3 = []\n",
    "    \n",
    "    if len(y_train.shape) == 1:\n",
    "        y_train = y_train.reshape(-1, 1)\n",
    "        y_token = 'single'\n",
    "    if pass_va == False:\n",
    "        print_importa = False\n",
    "    else:\n",
    "        print_importa = True\n",
    "    predict3 = three_models( x_train,y_train,x_test,test_method,data_name,x_train_names, print_importa = print_importa,print_model = False,recon =True,random_no=random_no)\n",
    "    model_name_list = []                 \n",
    "    ## starting stack!!\n",
    "    f_list = []\n",
    "    if (data_name != 'test'):\n",
    "        for ii in range(3):\n",
    "            if trastacks == True:\n",
    "                model_num = 4\n",
    "                rfr2,best_no, quit_lightGBM = mulit_model_stack2(x_train,y_train)\n",
    "\n",
    "                acc3_l = []\n",
    "                acc3_l.append(np.array(predict3[0][ii]))\n",
    "                acc3_l.append(np.array(predict3)[1][ii])\n",
    "                acc3_l.append(np.array(predict3)[2][ii])        #print(np.array(acc3_l).shape)\n",
    "                if quit_lightGBM == True:\n",
    "                    acc3_l.append(rfr2.predict((np.array(acc3_l).T)[:,:1]).reshape(-1,1))\n",
    "                else:\n",
    "                    acc3_l.append(rfr2.predict(np.array(acc3_l).T).reshape(-1,1)\n",
    "                                   )\n",
    "\n",
    "            else:\n",
    "                model_num = 3\n",
    "                print('without stack process...')\n",
    "\n",
    "            for i in range(model_num):\n",
    "                if i == 3:\n",
    "                    i = -1\n",
    "                model_name = ['rfr','etr','lightGBM','stack'][i]+['mean','up','down'][ii]\n",
    "                #pd.DataFrame(predict3[i]).to_csv('./recon/%s_%s.csv'%(data_name,model_name))\n",
    "                print('created the %s_%s model output'%(data_name,model_name))\n",
    "                model_name_list.append(model_name)\n",
    "                f_list.append(acc3_l[i])\n",
    "            t_pd = pd.DataFrame(f_list).T\n",
    "        ## reconstructions has been done\n",
    "\n",
    "        linear_w =rfr2.coef_\n",
    "        t_pd = pd.DataFrame(f_list).T\n",
    "        t_pd.columns = model_name_list\n",
    "        if model_name_list[-1] =='stack':\n",
    "            #t_pd['stack'] =  t_pd['stack'].apply(lambda x: x.replace('[','').replace(']','')) \n",
    "            t_pd['stack'] =  t_pd['stack'].apply(lambda x: x[0]) \n",
    "            t_pd['stack'] = t_pd['stack'].astype(float)\n",
    "\n",
    "            #linear_w.index = \n",
    "            if quit_lightGBM:\n",
    "                linear_w.columns =  ['rfr','etr']\n",
    "            else:\n",
    "                linear_w.columns =   ['rfr','etr','lightGBM']\n",
    "            try:\n",
    "                for i_index in ['stackmean','stackup','stackdown']:\n",
    "                    t_pd[i_index] =  t_pd[i_index].apply(lambda x: x[0]) \n",
    "                    t_pd[i_index] = t_pd[i_index].astype(float)\n",
    "            except:\n",
    "                print('without uq')\n",
    "        t_pd.index = cores_index\n",
    "        t_pd.to_csv('%s_core.csv'%(data_name))\n",
    "        linear_w = pd.DataFrame(np.array(linear_w).reshape(1,-1))\n",
    "        print('weights:')\n",
    "        print(linear_w)\n",
    "        linear_w.to_csv('%s_weights.csv'%(data_name))\n",
    "    else:\n",
    "        if trastacks == True:\n",
    "            model_num = 4\n",
    "            rfr2,best_no, quit_lightGBM = mulit_model_stack2(x_train,y_train)\n",
    "\n",
    "            acc3_l = []\n",
    "            acc3_l.append(np.array(predict3[0]))\n",
    "            acc3_l.append(np.array(predict3[1]))\n",
    "            acc3_l.append(np.array(predict3[2]))\n",
    "            #print(np.array(acc3_l).shape)\n",
    "            if quit_lightGBM == True:\n",
    "                predict3.append(rfr2.predict((np.array(acc3_l).T)[:,:1]))\n",
    "            else:\n",
    "                predict3.append(rfr2.predict(np.array(acc3_l).T))\n",
    "            t_pd = pd.DataFrame(predict3).T\n",
    "        else:\n",
    "                model_num = 3\n",
    "                print('without stack process...')\n",
    "    return(t_pd)\n",
    "\n",
    "\n",
    "\n",
    "def mulit_model_stack2(x_train,y_train,GBR = False,test_method = 'loo',random_no=0):\n",
    "    \n",
    "    test4 = []\n",
    "    predict4 = []\n",
    "    predict5_l = []\n",
    "    best_no = []\n",
    "    acc_list = []\n",
    "    acc_list2 = []\n",
    "\n",
    "    x_train=np.array(x_train)  #\n",
    "    y_train = np.array(y_train)\n",
    "    kf = KFold(n_splits=5)\n",
    "\n",
    "    \n",
    "    base_list = np.array(range(y_train.shape[0]))\n",
    "    ran_index = np.random.permutation(base_list)\n",
    "    x_train = x_train[ran_index]\n",
    "    y_train =y_train[ran_index]\n",
    "    pre_l = []\n",
    "    for train, test in kf.split(ran_index):\n",
    "        predict5 =[]\n",
    "        \n",
    "        try:\n",
    "            x_train2, x_test2, y_train2, y_test2 = x_train[train],x_train[test],y_train[train],y_train[test]\n",
    "        except:\n",
    "            x_train2, x_test2, y_train2, y_test2 = x_train.loc(axis=0)[train],x_train.loc(axis=0)[test],y_train.loc(axis=0)[train],y_train.loc(axis=0)[test]\n",
    "            \n",
    "        predict3 = three_models( x_train2,y_train2,x_test2, print_importa = False,print_model = False,recon = False,random_no=random_no)\n",
    "        #print('pred',np.array(predict3).shape)\n",
    "        test4.append(y_test2.reshape(-1,1))\n",
    "        predict4.append(np.array(predict3).T)\n",
    "        # predict3 (3, 4)\n",
    "        use_inf = True\n",
    "        if test_method == 'loo':\n",
    "                predict5.append(np.array(predict3)[0,:])\n",
    "                predict5.append(np.array(predict3)[1,:])\n",
    "                predict5.append(np.array(predict3)[2,:])        \n",
    "        else:\n",
    "            if use_inf == True:\n",
    "                predict5.append(np.array(predict3)[0,:])\n",
    "                predict5.append(np.array(predict3)[1,:])\n",
    "                predict5.append(np.array(predict3)[2,:])\n",
    "            else:\n",
    "                \n",
    "                predict5.append((np.array(predict3)[0,:]+np.array(predict3)[1,:]).T/2)\n",
    "                predict5.append((np.array(predict3)[0,:]+np.array(predict3)[2,:]).T/2)\n",
    "                predict5.append((np.array(predict3)[1,:]+np.array(predict3)[2,:]).T/2)\n",
    "                predict5.append((np.array(predict3)[1,:]+np.array(predict3)[2,:]+np.array(predict3)[0,:]).T/3)\n",
    "       \n",
    "        predict5_l.append(np.array(predict5).T)\n",
    "\n",
    "        acc2 = []\n",
    "        predict5= np.array(predict5).T\n",
    "\n",
    "        \n",
    "    rfr2 = LinearRegression() #ARDRegression() LinearRegression() \n",
    "    #n_estimators = 1000, min_sample_leaf10)\n",
    "    predict5_lt = np.vstack(predict5_l)\n",
    "    #print('pred2',predict5_lt.shape)\n",
    "\n",
    "    test_lt = np.vstack(test4)\n",
    "\n",
    "    quit_lightGBM = False\n",
    "\n",
    "    rfr2.fit(np.vstack(predict5_lt),test_lt) \n",
    "    return(rfr2,best_no,quit_lightGBM)\n",
    "\n",
    "def three_models(x_train,y_train,x_test,test_method =[],data_name=[],x_train_names=[],GBR =False, print_importa = True,print_model = False,n_estimators = 1000,random_no = 123,recon = False):\n",
    "    from numpy import random\n",
    "    #print(x_train.shape,'three_models')\n",
    "    predict3 = []\n",
    "    \n",
    "    y_train = np.array(y_train)\n",
    "    \n",
    "    \n",
    "    rfr = RandomForestRegressor(n_estimators = 1000,n_jobs=-2, random_state = random.seed(random_no))#n_estimators , min_sample_leaf10)\n",
    "    # training\n",
    "    rfr.fit(x_train,y_train)\n",
    "    rfr_feature_importances_ = rfr.feature_importances_\n",
    "    if recon & (data_name != 'test'):\n",
    "        rfr = MetamodelRegression(base_model=rfr,base_config=None)\n",
    "        rfr.fit(X=x_train, y=y_train.reshape(-1), base_is_prefitted=True)\n",
    "        #rfr.fit(x_train,y_train) #.reshape((-1, 1)),\n",
    "        # predict and save forecast result\n",
    "    predict3.append(rfr.predict(x_test)) \n",
    "\n",
    "\n",
    "\n",
    "    etr = ExtraTreesRegressor(n_estimators = n_estimators,n_jobs=-2, random_state = random.seed(random_no))\n",
    "    etr.fit(x_train,y_train) #.reshape((-1, 1)),\n",
    "    etr_feature_importances_ = etr.feature_importances_\n",
    "    if recon & (data_name != 'test'):\n",
    "        etr = MetamodelRegression(base_model=etr,base_config=None)\n",
    "        etr.fit(X=x_train, y=y_train.reshape(-1), base_is_prefitted=True)\n",
    "    predict3.append(etr.predict(x_test)) \n",
    "\n",
    "    if GBR == True:\n",
    "        gbr = GradientBoostingRegressor(n_estimators = n_estimators,n_jobs=-2, random_state = random.seed(random_no))\n",
    "        gbr.fit(x_train,y_train.reshape(-1)) #.reshape((-1, 1)),\n",
    "        \n",
    "        predict3.append(gbr.predict(x_test)[0])\n",
    "        predict3 = np.array(predict3)\n",
    "        print('the 3rd model is GBR')\n",
    "\n",
    "        imp_table.columns = ['rfr','etr']\n",
    "    else:\n",
    "        gbm = LGBMRegressor(objective='regression', num_leaves=31, learning_rate=0.05, n_estimators=1000,n_jobs=-2,random_state = random.seed(random_no),verbose = -1)\n",
    "\n",
    "        gbm.fit(x_train,y_train) #.reshape((-1, 1)),\n",
    "        gbm_feature_importances_ =gbm.feature_importances_\n",
    "    if recon & (data_name != 'test'):\n",
    "        gbm = MetamodelRegression(base_model=gbm,base_config=None)\n",
    "        gbm.fit(X=x_train, y=y_train.reshape(-1), base_is_prefitted=True)\n",
    "\n",
    "    predict3.append(gbm.predict(x_test))\n",
    "    imp_table = pd.DataFrame([rfr_feature_importances_,etr_feature_importances_,gbm_feature_importances_]).T\n",
    "    imp_table.columns = ['rfr','etr','LightGBM']\n",
    "    if print_importa == True:\n",
    "        \n",
    "\n",
    "        if  len(x_train_names) >0:\n",
    "            imp_table.index = x_train_names\n",
    "        imp_table = imp_table.sort_values(['rfr'],ascending = False)\n",
    "        imp_table.to_csv('importances_%s_%s.csv'%(test_method,data_name))\n",
    "    if print_model == False:\n",
    "        return(predict3)\n",
    "    else:\n",
    "        return(predict3,rfr,etr,model3)\n",
    "def wf_m(x,    xmax = 2 ,alpha = 0.75):\n",
    "    x = x.clone().detach()\n",
    "    x[x >= xmax] =  100\n",
    "\n",
    "    x[x < xmax] =  (x[x < xmax]/xmax)**alpha\n",
    "    x[x == 100] =  1\n",
    "    return(x)\n",
    "\n",
    "def make_recons_mee(sp_m,env,fos,plot_name,random_no=314,ran_time =999,n_estimators = 1000,num_epochs=10000,no_components=256,emb_method = 'plus'):\n",
    "    if len(env.shape) < 2:\n",
    "        env = pd.DataFrame(env.values.reshape(-1,1))\n",
    "    ran_time = 999\n",
    "    if emb_method=='plus':\n",
    "        no_components = 1\n",
    "    fake_raw = fos.shape[0]\n",
    "    fake_col = ran_time\n",
    "\n",
    "    envss = np.random.rand(env.shape[0], ran_time)\n",
    "    envss = np.hstack([env,envss])\n",
    "\n",
    "    imp_raw = sp_m.shape[-1]\n",
    "    rfr_imp = np.zeros([imp_raw, fake_col])\n",
    "    etr_imp = np.zeros([imp_raw, fake_col])\n",
    "    gbm_imp = np.zeros([imp_raw, fake_col])\n",
    "    rfr_recon = np.zeros([fake_raw, fake_col])\n",
    "    etr_recon = np.zeros([fake_raw, fake_col])\n",
    "    gbm_recon = np.zeros([fake_raw, fake_col])\n",
    "    stack_recon = np.zeros([fake_raw, fake_col])\n",
    "    \n",
    "    x_train = sp_m\n",
    "    x_test = fos\n",
    "    for i_num  in tqdm(range(fake_col), desc='Progress'):\n",
    "        random_no = i_num+100\n",
    "        y_train = pd.DataFrame(envss[:,i_num].reshape(-1,1))\n",
    "        recons = pre_glove(x_train,y_train.iloc[:],cores=x_test,save_vali=False,glove_abuna=True,emb_method='plus',recon=True, num_epochs=num_epochs,no_components=no_components,pass_va=True,data_name='test',random_no=random_no)    \n",
    "        recons = recons.values\n",
    "        rfr_recon[:,i_num] = recons[:,0]\n",
    "        etr_recon[:,i_num] = recons[:,1]\n",
    "        gbm_recon[:,i_num] = recons[:,2]\n",
    "        stack_recon[:,i_num] = recons[:,3]\n",
    "        if i_num % 10 == 0:\n",
    "            print('%s precentage ...'%(int(i_num/10)))\n",
    "        clear_output(wait=True)\n",
    "    pd.DataFrame(rfr_recon).to_csv('rfr_recon_%s.csv'%plot_name)\n",
    "    pd.DataFrame(etr_recon).to_csv('efr_recon_%s.csv'%plot_name)\n",
    "    pd.DataFrame(gbm_recon).to_csv('bgm_recon_%s.csv'%plot_name)\n",
    "    pd.DataFrame(stack_recon).to_csv('stack_recon_%s.csv'%plot_name)\n",
    "    return([])\n",
    "\n",
    "def transfer_pre_sqrt(input_ma):\n",
    "    #for ii in range(input_ma.shape[0]):\n",
    "    #    input_ma.iloc[ii,:] =     input_ma.iloc[ii,:]/input_ma.iloc[ii,:].sum()\n",
    "    input_ma = input_ma.div(input_ma.sum(axis=1), axis=0)\n",
    "    return(np.sqrt(input_ma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d23525a",
   "metadata": {},
   "source": [
    "# validation and reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06641007",
   "metadata": {},
   "outputs": [],
   "source": [
    "#SWAP\n",
    "\n",
    "no_components = 256\n",
    "num_epochs = 10000\n",
    "date = 'xxx'\n",
    "name = 'rlgh3'\n",
    "random_no = 0\n",
    "swap_sp= pd.read_csv(new_path+'\\\\'+'swap_spec.csv',index_col = 0)  #/100\n",
    "swap_ph= pd.read_csv(new_path+'\\\\'+'swap_ph.csv',index_col = 0)\n",
    "rlgh3_spec= pd.read_csv(new_path+'\\\\'+'rlgh3.csv',index_col = 0)\n",
    "rlgh_spec= pd.read_csv(new_path+'\\\\'+'rlgh_spec.csv',index_col = 0)\n",
    "#pre_glove(mp20_spec,mp20_env,cores = mp20_core,no_components=64,data_name='em_64_con_mat_ea_1207',test_method = 'cv',glove_abuna = True,recon=True,pass_va = True)\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh3_spec/100,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh3_spec/100,emb_method = 'only',stack = True,no_components=no_components,data_name='%s_e%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh3_spec/100,emb_method = 'none',stack = True,no_components=1,data_name='%s_a%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = False,recon=True,pass_va = True,random_no=random_no)\n",
    "name = 'rlgh'\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh_spec/100,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh_spec/100,emb_method = 'only',stack = True,no_components=no_components,data_name='%s_e%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh_spec/100,emb_method = 'none',stack = True,no_components=1,data_name='%s_a%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c9bd413",
   "metadata": {},
   "outputs": [],
   "source": [
    "#SMPDSv1\n",
    "\n",
    "no_components = 256 \n",
    "num_epochs = 10000\n",
    "name = 'po19'\n",
    "date =  'xxx'\n",
    "no_components = 256\n",
    "random_no = 0\n",
    "po19_envs = pd.read_csv(new_path+'\\\\'+'Pollen site Climate.csv', index_col = 0).MTCO\n",
    "po19_sp = pd.read_csv(new_path+'\\\\'+'SMPDS_Feb2019.csv',index_col = 0,encoding = 'ISO-8859-1').iloc[:,4:]#.drop(['Betula'],axis = 1)\n",
    "po19_core = pd.read_csv(new_path+'\\\\'+'po19_core.csv',index_col = 0,encoding = 'ISO-8859-1')#.drop(['Betula'],axis = 1)\n",
    "for i in range(po19_core.shape[0]):\n",
    "    po19_core.iloc[i,:] = po19_core.iloc[i,:]/po19_core.iloc[i,:].sum()\n",
    "pre_glove(po19_sp,po19_envs,cores = po19_core,emb_method = 'plus',stack = True,num_epochs= num_epochs,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(po19_sp,po19_envs,cores = po19_core,emb_method = 'only',stack = True,num_epochs= num_epochs,no_components=no_components,data_name='%s_e%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "pre_glove(po19_sp,po19_envs,cores = po19_core,emb_method = 'none',stack = True,num_epochs= num_epochs,no_components=1,data_name='%s_a%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53e074da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NIMBIOS\n",
    "\n",
    "no_components = 256\n",
    "num_epochs = 10000\n",
    "name = 'mos'\n",
    "date = 'xxx'\n",
    "mp20_spec= pd.read_csv(new_path+'\\\\'+'MP20_abund.csv',index_col = 1, encoding ='gbk').iloc[:,1:]\n",
    "mp20_env= pd.read_csv(new_path+'\\\\'+'MP20_env.csv',index_col = 0).MAT/10\n",
    "mp20_core= pd.read_csv(new_path+'\\\\'+'mos.csv',index_col = 0,sep = ';').iloc[:,2:].loc['CON1']\n",
    "\n",
    "pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= 10000,glove_abuna = True,recon=True,pass_va = False)\n",
    "pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'only',stack = True,no_components=no_components,data_name='%s_e%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= 10000,glove_abuna = True,recon=True,pass_va = False)\n",
    "pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'none',stack = True,no_components=1,data_name='%s_a%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = False)\n",
    "#mp20_core = pd.read_csv(new_path+'\\\\'+'lla_core_r.csv',index_col = 0,sep = ';')\n",
    "#name = 'lla'\n",
    "\n",
    "#pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= 10000,glove_abuna = True,recon=True,pass_va = True)\n",
    "#pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'only',stack = True,no_components=no_components,data_name='%s_e%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= 10000,glove_abuna = True,recon=True,pass_va = True)\n",
    "#pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'none',stack = True,no_components=1,data_name='%s_a%s_%s'%(name,no_components,date),test_method = 'cv',glove_abuna = True,recon=True,pass_va = True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e250bc1a",
   "metadata": {},
   "source": [
    "# testing the robustness of a reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5350bc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# swap\n",
    "\n",
    "swap_sp= pd.read_csv(new_path+'\\\\'+'swap_spec.csv',index_col = 0)  #/100\n",
    "swap_ph= pd.read_csv(new_path+'\\\\'+'swap_ph.csv',index_col = 0)\n",
    "rlgh3_spec= pd.read_csv(new_path+'\\\\'+'rlgh3.csv',index_col = 0)\n",
    "rlgh_spec= pd.read_csv(new_path+'\\\\'+'rlgh_spec.csv',index_col = 0)\n",
    "\n",
    "pre_glove(swap_sp/100,swap_ph,cores = rlgh_spec/100,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "emb_method = 'plus'\n",
    "make_recons_mee(swap_sp,swap_ph,rlgh_spec/100,'rlgh_test_comb',emb_method = emb_method)\n",
    "emb_method = 'none'\n",
    "make_recons_mee(swap_sp,swap_ph,rlgh_spec/100,'rlgh_test_abund',emb_method = emb_method)\n",
    "\n",
    "emb_method = 'plus'\n",
    "make_recons_mee(swap_sp,swap_ph,rlgh3_spec/100,'rlgh3_test_comb',emb_method = emb_method)\n",
    "emb_method = 'none'\n",
    "make_recons_mee(swap_sp,swap_ph,rlgh3_spec/100,'rlgh3_test_abund',emb_method = emb_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47ef1c44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NIMBIOS\n",
    "\n",
    "mp20_spec= pd.read_csv(new_path+'\\\\'+'MP20_abund.csv',index_col = 1, encoding ='gbk').iloc[:,1:]\n",
    "mp20_env= pd.read_csv(new_path+'\\\\'+'MP20_env.csv',index_col = 0).MAT/10\n",
    "mp20_core= pd.read_csv(new_path+'\\\\'+'mos.csv',index_col = 0,sep = ';').iloc[:,2:].loc['CON1']\n",
    "mp20_spec = transfer_pre_sqrt(mp20_spec)\n",
    "\n",
    "pre_glove(mp20_spec,mp20_env,cores = mp20_core,emb_method = 'plus',stack = True,no_components=no_components,data_name='%s_ae%s_%s'%(name,no_components,date),test_method = 'cv',num_epochs= num_epochs,glove_abuna = True,recon=True,pass_va = True,random_no=random_no)\n",
    "emb_method = 'plus'\n",
    "make_recons_mee(mp20_spec,mp20_env,cores = mp20_core,'mos_test_comb',emb_method = emb_method)\n",
    "emb_method = 'none'\n",
    "make_recons_mee(mp20_spec,mp20_env,cores = mp20_core,'mos_test_abund',emb_method = emb_method)\n",
    "\n",
    "\n",
    "\n",
    "mp20_core= pd.read_csv(new_path+'\\\\'+'lla.csv',index_col = 0,sep = ';').iloc[:,2:].loc['CON1']\n",
    "mp20_core = transfer_pre_sqrt(mp20_core)\n",
    "\n",
    "emb_method = 'plus'\n",
    "make_recons_mee(mp20_spec,mp20_env,cores = mp20_core,'lla_test_comb',emb_method = emb_method)\n",
    "emb_method = 'none'\n",
    "make_recons_mee(mp20_spec,mp20_env,cores = mp20_core,'lla_test_abund',emb_method = emb_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d03c06c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#SMPDSv1\n",
    "\n",
    "po19_envs = pd.read_csv(new_path+'\\\\'+'Pollen site Climate.csv', index_col = 0).MTCO\n",
    "po19_sp = pd.read_csv(new_path+'\\\\'+'SMPDS_Feb2019.csv',index_col = 0,encoding = 'ISO-8859-1').iloc[:,4:]\n",
    "po19_core = pd.read_csv(new_path+'\\\\'+'po19_core_r.csv',index_col = 0,encoding = 'ISO-8859-1')\n",
    "for i in range(po19_core.shape[0]):\n",
    "    po19_core.iloc[i,:] = po19_core.iloc[i,:]/po19_core.iloc[i,:].sum()\n",
    "po19_sp = transfer_pre_sqrt(po19_sp)\n",
    "po19_core = transfer_pre_sqrt(po19_core)\n",
    "emb_method = 'plus'\n",
    "make_recons_mee(po19_sp,po19_envs,po19_core,'po19_test_comb',emb_method = emb_method)\n",
    "emb_method = 'none'\n",
    "make_recons_mee(po19_sp,po19_envs,po19_core,'po19_test_abund',emb_method = emb_method)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d362fad",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "uq360",
   "language": "python",
   "name": "uq360"
  },
  "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}