<jupyter_start><jupyter_text>Demo (2 groups)In this demo, it will be demonstrated how to run our GFA implementation using incomplete synthetic data generated (with 2 groups) from the generative model $\mathbf{X}^{(m)} \sim \mathcal{N}(\mathbf{W}^{(m)}\mathbf{Z}, \mathbf{T}^{(m)^{-1}})$. Different experiments can be run by changing the initial arguments.<jupyter_code>import numpy as np
import time
import pickle
import os
import copy
import argparse
import visualization_syntdata
from models import GFA_DiagonalNoiseModel, GFA_OriginalModel
from utils import GFAtools<jupyter_output><empty_output><jupyter_text>Define intial setting Set the initial arguments to generate the data and run the model. The defaults values are set to apply the GFA extended model (``--noise=diagonal``) to incomplete data (``--scenario=incomplete``) with 2 groups (``--num-groups=2``). A second model will be run with complete data after the missing values are imputed using the median of the corresponding features (``--impMedian=True``). Both models are initialised with 10 latent factors (``--K=10``). The original model can be applied to complete data by setting the following values: ``--scenario=complete``, ``--noise=diagonal`` and ``--impMedian=False``.<jupyter_code>def get_args():
    parser = argparse.ArgumentParser(description="GFA with two groups")
    parser.add_argument("--scenario", nargs='?', default='incomplete', type=str,
                        help='Data scenario (complete or incomplete)')
    parser.add_argument("--noise", nargs='?', default='diagonal', type=str,
                        help='Noise assumption for GFA models (diagonal or spherical)')
    parser.add_argument("--num-groups", nargs='?', default=2, type=int,
                        help='Number of groups')
    parser.add_argument("--K", nargs='?', default=10, type=int,
                        help='number of latent factors to initialise the model')
    parser.add_argument("--impMedian", nargs='?', default=True, type=bool,
                        help='(not) impute median')
    
    args = parser.parse_args("")                                                               
    return args<jupyter_output><empty_output><jupyter_text>Define the parameters to generate incomplete data sets. In this example, $10\%$ of the rows are randomly removed from the first group ($\textbf{X}^{(1)}$) and $20\%$ of the elements are randomly removed from the second group ($\textbf{X}^{(2)}$). The results of the experiments will be saved in ``res_dir``.<jupyter_code># Get arguments 
args = get_args() 
# Define parameters to generate incomplete data sets
infoMiss = {'perc': [10,20], #percentage of missing data 
        'type': ['rows','random'], #type of missing data 
        'ds': [1,2]} #groups that will have missing values        
    
# Make directory to save the results of the experiments         
res_dir = f'results/{args.num_groups}groups/GFA_{args.noise}/{args.K}comps/{args.scenario}'
print(res_dir)
if not os.path.exists(res_dir):
        os.makedirs(res_dir)<jupyter_output><empty_output><jupyter_text>Generate dataThe following functions will be used to generate (complete or incomplete) synthetic data with 2 groups. The number of true latent factors is 4: 2 shared factors between groups (first and second factors), 1 specific to the second group (third factor) and 1 specific to the first group (fourth factor). These are manually specified for visualization purposes. The noise precisions and ARD parameters are also manually specified. The loading matrices are generated from the prior: $p(\mathbf{W}|\boldsymbol{\alpha}) = \prod_{m=1}^{M} \prod_{j=1}^{D_{m}} \prod_{k=1}^{K} \mathcal{N}(w^{(m)}_{jk}|0, (\alpha_{k}^{(m)})^{-1})$. Finally, the groups were generated from $\textbf{X}^{(m)} \sim \mathcal{N}(\mathbf{W}^{(m)}\mathbf{Z}, \mathbf{T}^{(m)^{-1}})$.<jupyter_code>def get_data_2g(args, infoMiss=None):

    """ 
    Generate synthetic data with 2 groups.

    Parameters
    ----------
    args : local namespace 
        Arguments selected to run the model.

    infoMiss : dict | None, optional.
        Parameters selected to generate data with missing values.  

    Returns
    -------
    data : dict
        Training and test data as well as model parameters used 
        to generate the data.
    
    """
    Ntrain = 400; Ntest = 100
    N = Ntrain + Ntest #  total number of samples
    M = args.num_groups  #number of groups
    d = np.array([50, 30]) #number of dimensios in each group
    true_K = 4  # true latent factors
    
    # Specify Z manually
    Z = np.zeros((N, true_K))
    for i in range(0, N):
        Z[i,0] = np.sin((i+1)/(N/20))
        Z[i,1] = np.cos((i+1)/(N/20))
        Z[i,2] = 2 * ((i+1)/N-0.5)    
    Z[:,3] = np.random.normal(0, 1, N)          
    
    # Specify noise precisions manually
    tau = [[] for _ in range(d.size)]
    tau[0] = 5 * np.ones((1,d[0]))[0] 
    tau[1] = 10 * np.ones((1,d[1]))[0]
    
    # Specify alphas manually
    alpha = np.zeros((M, true_K))
    alpha[0,:] = np.array([1,1,1e6,1])
    alpha[1,:] = np.array([1,1,1,1e6])     
    
    #W and X
    W = [[] for _ in range(d.size)]
    X_train = [[] for _ in range(d.size)]
    X_test = [[] for _ in range(d.size)]
    for i in range(0, d.size):
        W[i] = np.zeros((d[i], true_K))
        for t in range(0, true_K):
            #generate W from p(W|alpha)
            W[i][:,t] = np.random.normal(0, 1/np.sqrt(alpha[i,t]), d[i])
        X = np.zeros((N, d[i]))
        for j in range(0, d[i]):
            #generate X from the generative model
            X[:,j] = np.dot(Z,W[i][j,:].T) + \
            np.random.normal(0, 1/np.sqrt(tau[i][j]), N*1)    
        # Get training and test data
        X_train[i] = X[0:Ntrain,:] #Training data
        X_test[i] = X[Ntrain:N,:] #Test data
    
    #latent variables for training the model    
    Z = Z[0:Ntrain,:]

    # Generate incomplete training data
    if args.scenario == 'incomplete':
        X_train, missing_Xtrue = generate_missdata(X_train, infoMiss)
    
    # Store data and model parameters            
    data = {'X_tr': X_train, 'X_te': X_test, 'W': W, 'Z': Z, 'tau': tau, 'alpha': alpha, 'true_K': true_K}
    if args.scenario == 'incomplete':
        data.update({'trueX_miss': missing_Xtrue}) 
    return data 

def generate_missdata(X_train, infoMiss):
    """ 
    Generate missing data in the training data.

    Parameters
    ----------
    X_train : list 
        List of arrays containing the data matrix of each group.

    infoMiss : dict 
        Parameters selected to generate data with missing values.  

    Returns
    -------
    X_miss : list 
        List of arrays containing the training data. The groups 
        specified in infoMiss will have missing values.

    missing_Xtrue : list 
        List of arrays containing the true values removed from the
        groups selected in infoMiss.     
    
    """
    missing_Xtrue = [[] for _ in range(len(infoMiss['ds']))]
    for i in range(len(infoMiss['ds'])):
        g_miss = infoMiss['ds'][i]-1  
        if 'random' in infoMiss['type'][i]: 
            #remove entries randomly
            missing_val =  np.random.choice([0, 1], 
                        size=(X_train[g_miss].shape[0],X_train[g_miss].shape[1]), 
                        p=[1-infoMiss['perc'][i-1]/100, infoMiss['perc'][i-1]/100])
            mask_miss =  np.ma.array(X_train[g_miss], mask = missing_val).mask
            missing_Xtrue[i] = np.where(missing_val==1, X_train[g_miss],0)
            X_train[g_miss][mask_miss] = 'NaN'
        elif 'rows' in infoMiss['type'][i]: 
            #remove rows randomly
            Ntrain = X_train[g_miss].shape[0]
            missing_Xtrue[i] = np.zeros((Ntrain, X_train[g_miss].shape[1]))
            n_rows = int(infoMiss['perc'][i-1]/100 * Ntrain)
            shuf_samples = np.arange(Ntrain)
            np.random.shuffle(shuf_samples)
            missing_Xtrue[i][shuf_samples[0:n_rows],:] = X_train[g_miss][shuf_samples[0:n_rows],:]
            X_train[g_miss][shuf_samples[0:n_rows],:] = 'NaN'
        X_miss = X_train
    return X_miss, missing_Xtrue 

# Generate data
if args.scenario == 'complete':
    data = get_data_2g(args)
else:
    data = get_data_2g(args, infoMiss)

#save file with generated data
data_file = f'{res_dir}/[1]Data.dictionary'
with open(data_file, 'wb') as parameters:
        pickle.dump(data, parameters)<jupyter_output><empty_output><jupyter_text>Run model 1Our GFA extension (``GFA_DiagonalNoiseModel``) is applied by selecting ``--noise=diagonal``. The original model (``GFA_OriginalModel``) is run by selecting ``--noise=spherical``. It should take less than 1 minute to run the code below.<jupyter_code># Initialise the model
X_tr = data['X_tr'] #get training data
params = {'num_groups': args.num_groups,
                      'K': args.K, 'scenario': args.scenario}
if 'diagonal' in args.noise:    
    GFAmodel = GFA_DiagonalNoiseModel(X_tr, params)
else:
    GFAmodel = GFA_OriginalModel(X_tr, params) 
    
#Fit the model
time_start = time.process_time()
GFAmodel.fit(X_tr)
GFAmodel.time_elapsed = time.process_time() - time_start
print(f'Computational time: {float("{:.2f}".format(GFAmodel.time_elapsed))}s')<jupyter_output><empty_output><jupyter_text>Predictions Multi-output predictionsThe test data of $\textbf{X}^{(2)}$ is predicted from the test data of $\textbf{X}^{(1)}$ using the ``PredictGroups`` method in [utils.py](utils.py). The MSE calculated between the true and predicted values (on the test set) should be small and better than chance.<jupyter_code># Predictions (Predict group 2 from group 1) 
# Compute mean squared error (MSE)
obs_ds = np.array([1, 0]) #group 1 was observed
gpred = np.where(obs_ds == 0)[0][0] #get the non-observed group
X_test = data['X_te']
X_pred = GFAtools(X_test, GFAmodel).PredictGroups(obs_ds, args.noise)

#Compute MSE 
MSE = np.mean((X_test[gpred] - X_pred[0]) ** 2)
print(f'MSE: {float("{:.2f}".format(MSE))}')

# Compute MSE - chance level (MSE between test values and train means)
Tr_means = np.ones((X_test[gpred].shape[0], X_test[gpred].shape[1])) * \
    np.nanmean(data['X_tr'][gpred], axis=0)           
MSE_chlev = np.mean((X_test[gpred] - Tr_means) ** 2)
print(f'MSE (chance level): {float("{:.2f}".format(MSE_chlev))}')

# Save file containing model outputs and predictions
res_file = f'{res_dir}/[1]ModelOutput.dictionary'
with open(res_file, 'wb') as parameters:
    pickle.dump(GFAmodel, parameters)<jupyter_output><empty_output><jupyter_text>Missing data predictionThe missing values are predicted using the ``PredictMissing`` method in [utils.py](utils.py). Pearson's correlation is computed between the true and missing values for each group.<jupyter_code># Predict missing values
if args.scenario == 'incomplete':
    Corr_miss = np.zeros((1,len(infoMiss['ds'])))
    missing_pred = GFAtools(data['X_tr'], GFAmodel).PredictMissing(infoMiss)
    missing_true = data['trueX_miss']
    for i in range(len(infoMiss['ds'])):
        Corr_miss[0,i] = np.corrcoef(missing_true[i][missing_true[i] != 0], 
                                        missing_pred[i][np.logical_not(np.isnan(missing_pred[i]))])[0,1]                   
        print('Group: ', infoMiss['ds'][i])
        print('Correlation: ', np.around(Corr_miss[0,i],3))<jupyter_output><empty_output><jupyter_text>Run model 2A new model is run with complete data after the missing values are imputed using the median. The test data of $\textbf{X}^{(2)}$ is predicted from the test data of $\textbf{X}^{(1)}$. The MSE between the test and predicted values is calculated. It should take a few seconds to run the code below.<jupyter_code>if args.impMedian: 
    #ensure scenario and noise were correctly selected
    assert args.scenario == 'incomplete' and args.noise == 'diagonal'
    X_impmed = copy.deepcopy(data['X_tr'])
    g_miss = np.array(infoMiss['ds']) - 1 #group with missing values  
    # Impute median before training the model
    for i in range(g_miss.size):
        for j in range(data['X_tr'][g_miss[i]].shape[1]):
            Xtrain_j = data['X_tr'][g_miss[i]][:,j]
            X_impmed[g_miss[i]][np.isnan(X_impmed[g_miss[i]][:,j]),j] = np.nanmedian(Xtrain_j)

    # Initialise the model
    GFAmodel_median = GFA_DiagonalNoiseModel(X_impmed, params, imputation=True)
    # Fit the model
    time_start = time.process_time()
    GFAmodel_median.fit(X_impmed)
    GFAmodel_median.time_elapsed = time.process_time() - time_start
    print(f'Computational time: {float("{:.2f}".format(GFAmodel_median.time_elapsed))}s')

    # Predictions (Predict group 2 from group 1) 
    obs_ds = np.array([1, 0]) #group 1 was observed 
    gpred = np.where(obs_ds == 0)[0][0] #get the non-observed group
    X_test = data['X_te']
    X_pred = GFAtools(X_test, GFAmodel_median).PredictGroups(obs_ds, args.noise)
    MSE_median = np.mean((X_test[gpred] - X_pred[0]) ** 2) 
    print(f'MSE (after imp. median): {float("{:.2f}".format(MSE_median))}')

    # Save file containing model outputs and predictions
    res_med_file = f'{res_dir}/[1]ModelOutput_median.dictionary'
    with open(res_med_file, 'wb') as parameters:
        pickle.dump(GFAmodel_median, parameters)<jupyter_output><empty_output><jupyter_text>VisualizationPlot and save the true and estimated model parameters in ``res_dir``.<jupyter_code># Import auxiliary functions for visualization
from visualization_syntdata import hinton_diag, find_relfactors, match_factors, plot_loadings, plot_Z
import matplotlib.pyplot as plt

model_1 = GFAmodel #model without imputation
model_2 = GFAmodel_median #model with median imputation

#Concatenate (true and estimated) loadings and alphas across groups
#model_1
W_est_1 = np.zeros((np.sum(model_1.d),model_1.k))
alphas_est_1 = np.zeros((model_1.k, args.num_groups))
#model_2
W_est_2 = np.zeros((np.sum(model_2.d),model_2.k))
alphas_est_2 = np.zeros((model_2.k, args.num_groups))
#true parameters
W_true = np.zeros((np.sum(model_1.d),data['true_K']))
alphas_true = np.zeros((data['true_K'], args.num_groups))
d = 0
for m in range(args.num_groups):
    Dm = model_1.d[m]
    alphas_true[:,m] = data['alpha'][m]
    W_true[d:d+Dm,:] = data['W'][m]
    #model_1
    alphas_est_1[:,m] = model_1.E_alpha[m]
    W_est_1[d:d+Dm,:] = model_1.means_w[m]
    #model_2
    alphas_est_2[:,m] = model_2.E_alpha[m]
    W_est_2[d:d+Dm,:] = model_2.means_w[m]
    d += Dm<jupyter_output><empty_output><jupyter_text>Plot and save loading matrices<jupyter_code>#plot true Ws
W_path = f'{res_dir}/[1]W_true.png'
plot_loadings(W_true, model_1.d, W_path) 

#plot estimated Ws - model 1
if model_1.k == data['true_K']:
    #match true and estimated factors
    match_res = match_factors(W_est_1, W_true)
    W_est_1 = match_res[0]
W_path_1 = f'{res_dir}/[1]W_est.png' 
plot_loadings(W_est_1, model_1.d, W_path_1)

#plot estimated Ws - model 2
if model_2.k == data['true_K']:
    #match true and estimated factors
    match_res = match_factors(W_est_2, W_true)
    W_est_2 = match_res[0]
W_path_2 = f'{res_dir}/[1]W_est_median.png' 
plot_loadings(W_est_2, model_2.d, W_path_2)<jupyter_output><empty_output><jupyter_text>Plot and save latent variables<jupyter_code>#plot true latent variables 
Z_path = f'{res_dir}/[1]Z_true.png'    
plot_Z(data['Z'], Z_path)

#plot estimated latent variables - model 1
Z_path_1 = f'{res_dir}/[1]Z_est.png'
if model_1.k == data['true_K']:
    simcomps = match_res[1]
    plot_Z(model_1.means_z[:, simcomps], Z_path_1, match=True, flip=match_res[2])
else:     
    plot_Z(model_1.means_z, Z_path_1)    
    
#plot estimated latent variables - model 2
Z_path_2 = f'{res_dir}/[1]Z_est_median.png'
if model_2.k == data['true_K']:
    simcomps = match_res[1]
    plot_Z(model_2.means_z[:, simcomps], Z_path_2, match=True, flip=match_res[2])
else:     
    plot_Z(model_2.means_z, Z_path_2)<jupyter_output><empty_output><jupyter_text>Plot and save alphas<jupyter_code>#plot true alphas
alphas_path = f'{res_dir}/[1]alphas_true.png'
hinton_diag(np.negative(alphas_true.T), alphas_path)     

#plot estimated alphas - model 1
alphas_path_1 = f'{res_dir}/[1]alphas_est.png'
if model_1.k == data['true_K']:
    hinton_diag(np.negative(alphas_est_1[simcomps,:].T), alphas_path_1)  
else:
    hinton_diag(np.negative(alphas_est_1.T), alphas_path_1)

#plot estimated alphas - model 2
alphas_path_2 = f'{res_dir}/[1]alphas_est_median.png'
if model_2.k == data['true_K']:
    hinton_diag(np.negative(alphas_est_2[simcomps,:].T), alphas_path_2) 
else:
    hinton_diag(np.negative(alphas_est_2.T), alphas_path_2)<jupyter_output><empty_output><jupyter_text>Plot and save Evidence Lower Bound (ELBO)<jupyter_code># Plot ELBO - model 1
L_path_1 = f'{res_dir}/[1]ELBO.png'    
plt.figure(figsize=(5, 4))
plt.plot(model_1.L[1:])
plt.savefig(L_path_1)
plt.close() 

# Plot ELBO - model 2
L_path_2 = f'{res_dir}/[1]ELBO_median.png'
plt.figure(figsize=(5, 4))
plt.plot(model_2.L[1:])
plt.savefig(L_path_2)
plt.close()<jupyter_output><empty_output><jupyter_text>Run analysis_syntdata.py As the variational approximations for GFA are deterministic and the model converges to a local optimum that depends on the initialisation, the algorithm must be initialised multiple times and the best initialisation is chosen based on the solution with largest variational lower bound. We can do it by running the python script [analysis_syntdata.py](analysis_syntdata.py). The data files and outputs of the models are saved in ``res_dir``. The parameters of both models are only plot for the best initialisation. In addition, a .txt file (``results.txt``) is saved with more information about each initialisation and overall results (e.g. multi-output and missing data prediction across the different runs). Here, we run the models 3 times but you can increase the number of initialisations by setting ``num-runs`` to a different number. It should take less than 2 minutes to run the code below.<jupyter_code>import shutil
shutil.rmtree(res_dir, ignore_errors=True)

%run analysis_syntdata.py --scenario=incomplete --noise=diagonal --num-groups=2 --K=10 --num-runs=3 --impMedian=True<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13578477314074647, 'fraction_numerical': 0.015546612271900325, 'mean_word_length': 3.798941798941799, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '370492', 'license': 'MIT', 'n_tokens_mistral': 6427, 'n_tokens_neox': 6179, 'n_words': 1599, 'path': 'demo_2groups.ipynb', 'repo_name': 'ferreirafabio80/gfa', 'size': '18139'}