You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: parcel_cl_nmds.m
    path of file: ./repos/CanlabCore/CanlabCore/Parcellation_tools
    
    the code of the file until where you have to start completion: function [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)
% :Usage:
% ::
%
%     [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)
%
% :Examples:
% ::
%
%    load Parcellation_info/parcellation.mat
%    [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)

    
function [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)
% :Usage:
% ::
%
%     [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)
%
% :Examples:
% ::
%
%    load Parcellation_info/parcellation.mat
%    [parcel_cl_avgs, NMDS, class_clusters] = parcel_cl_nmds(parcel_cl_avgs)
%
%    parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save')
%    parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save', 'savedir', 'Parcellation_info')
%
% :Complete methods example:

%% 1. Dimension reduction
% 
% Clustering of multivariate data is most stable when the data is not sparse, i.e., the dimensionality is low relative to the number of observations. To limit the dimensionality of the data, a spatio-temporal dimension-reduction step is first performed on the [n x v x N] data matrix of  AUC data for n trials x v voxels x N participants (here, n = 48 trials [usually], v = 17,112 voxels , and N = 26 participants).  A temporal data reduction is first performed to identify components with correlated AUC trial time series within each participant, followed by a spatial reduction to identify components with correlated spatial patterns across subjects.  First, the [n x v] matrix of AUC data for each participant was subjected to PCA, using the [v x v] correlation matrices.  Based on the scree plots across subjects, we saved the first 7 eigenvectors (spatial maps).  These eigenvectors explained 74 +- 2.5% (st. dev. across subjects) of the variance in the full dataset. These eigenvectors were scaled by their variances (eigenvalues) and concatenated across subjects to form an [v x N*7] matrix of eigenvectors.   This matrix was subjected to second (across participant) PCA step to identify components with similar spatial maps across participants.  We retained 12 eigenvectors (maps) based on the scree plot, which explained 65% of the variance across individuals. Component scores in this space were used for clustering.  This is a data reduction step, and the results are not expected to depend