public static void schmiedeckEtAlScript(boolean useDetrend, boolean useFirstFactor, boolean useTaskCombination, boolean ignoreBetweenInMDS) {
        
        // Setting constants and reserving memory space
        String betweenDataName = "between";
        String withinDataName = "raw_person";
        int anzPerFull = 101, anzVar = 18, corSize = (useFirstFactor || useTaskCombination?9:18);
        
        double[] avgPersonCorWthn = new double[anzPerFull+1]; Statik.setTo(avgPersonCorWthn, -999.0);           
        double[] avgPersonCorBtw = new double[anzPerFull+1]; Statik.setTo(avgPersonCorBtw, -999.0);   
        double[] highestEValCorWthn = new double[anzPerFull+1]; Statik.setTo(highestEValCorWthn, -999.0);
        double[] highestEValCorBtwn = new double[anzPerFull+1]; Statik.setTo(highestEValCorBtwn, -999.0);
        double[] highestEValCorSmooth = new double[anzPerFull+1]; Statik.setTo(highestEValCorBtwn, -999.0);
        double[][] highestEVecCorWthn = new double[anzPerFull+1][corSize]; Statik.setTo(highestEVecCorWthn, -999.0);
        double[][] highestEVecCorBtwn = new double[anzPerFull+1][corSize]; Statik.setTo(highestEVecCorBtwn, -999.0);
        double[][] highestEVecCorSmooth = new double[anzPerFull+1][corSize]; Statik.setTo(highestEVecCorBtwn, -999.0);
        double[][] flatFactorCoefficient = new double[anzPerFull+1][corSize]; Statik.setTo(flatFactorCoefficient, -999.0);
        double[][] preferences = Statik.loadMatrix("FloData"+File.separator+"cogito_preferences.dat", '\t');
        int[] omit = new int[]{};
        int nextOmit = 0;
        Vector<Integer> id = new Vector<Integer>();
        Vector<double[][]> corv = new Vector<double[][]>();
        double[][] wthnAvgCov = new double[corSize][corSize], btwnAvgCov = new double[corSize][corSize];
        double[][][] wthnCorMem = new double[anzPerFull][corSize][corSize];
        double[][][] smoothCorMem = new double[anzPerFull][corSize][corSize];
        double[][] wthnDataSmoothRaise = new double[anzPerFull][corSize];
        double[][] wthnMean = new double[anzPerFull][corSize];
        double[][] wthnReducedAvgCov = new double[corSize][corSize], btwnReducedAvgCov = new double[corSize][corSize];
        double[][][] pretestData = new double[anzPerFull][][];
        double[][] pretestMean = new double[anzPerFull][];
        double[][][] pretestCov = new double[anzPerFull][][];
        double[][][] combinedData = new double[anzPerFull][][];
        double[][] combinedMean = new double[anzPerFull][corSize];
        double[][][] combinedWthnCov = new double[anzPerFull][corSize][corSize];

        /*
         * In this loop, the covariance and correlation matrices "within" and "between" are computed
         * for each participant. The within covariance matrix is based on the participant's longitudinal
         * data, the between covariance matrix is based on the pretest data for all participants with the 
         * same time settings as the current participant. Depending on the script parameters, the data is either raw or detrended,
         * and the two indicators for each variable are either left separate or collapsed into one. 
         * Participants can be omitted from the analysis using the omit vector. All participants with numerical errors
         * are reported to standard output.
         */ 
        System.out.print("Failed IDs = ");
        for (int partNr=1; partNr<=anzPerFull; partNr++) if (nextOmit >= omit.length || partNr != omit[nextOmit]) {
            double[][] wthnDataWithID = Statik.loadMatrix("FloData"+File.separator+withinDataName+"_"+partNr+".dat", ' ', true, ".");
            int[] allButOne = new int[anzVar]; for (int i=0; i<anzVar; i++) allButOne[i] = i+1;
            double[][] wthnData = Statik.submatrix(wthnDataWithID, null, allButOne);
            double[][] wthnDataSmooth = Statik.gaussSmoothening(wthnData, 3.0, 0.0, Model.MISSING, true);
            double[][] wthnDataDetrend = new double[wthnData.length][wthnData[0].length];
            for (int i=0; i<wthnData.length; i++) for (int j=0; j<wthnData[i].length; j++) 
                wthnDataDetrend[i][j] = (wthnData[i][j]==Model.MISSING?Model.MISSING:wthnData[i][j]-wthnDataSmooth[i][j]);
            
            int k = anzVar / corSize;
            for (int i=0; i<corSize; i++) {
                wthnDataSmoothRaise[partNr-1][i] = 0;
                for (int j=0; j<k; j++) wthnDataSmoothRaise[partNr-1][i] += wthnDataSmooth[wthnDataSmooth.length-1][k*i+j]-wthnDataSmooth[0][k*i+j];
                wthnDataSmoothRaise[partNr-1][i] /= (double)k;
            }
            double[][] wthnCov = new double[anzVar][anzVar]; double[] wthnMeanAll = new double[anzVar];
            Statik.covarianceMatrixAndMeans((useDetrend?wthnDataDetrend:wthnData), wthnMeanAll, wthnCov, Model.MISSING);
            double[][] wthnCor = Statik.correlationFromCovariance(wthnCov);
            globalWarning = false;
            if (useFirstFactor) {
                wthnCor = computePairwiseFactorCorrelationSimple(wthnCor);
                if (globalWarning) System.out.print(partNr+"\t");
            }
            if (useTaskCombination) {
                wthnCor = computePairwiseCombinedCorrelationSimple(wthnCor);
                if (globalWarning) System.out.print(partNr+"\t");
            }
            if (!globalWarning) {
                corv.add(wthnCor);
                id.add(partNr);
                wthnAvgCov = Statik.add(wthnAvgCov, wthnCor);
                wthnCorMem[partNr-1] = Statik.copy(wthnCor);
                if (useFirstFactor || useTaskCombination) wthnMean[partNr-1] = Statik.copy(wthnMeanAll); 
                else for (int i=0; i<corSize; i++) wthnMean[partNr-1][i] = (wthnMeanAll[2*i]+wthnMeanAll[2*i+1])/2.0;
            }
            
            double[][] smoothCov = new double[anzVar][anzVar]; double[] smoothMean = new double[anzVar];
            Statik.covarianceMatrixAndMeans(wthnDataSmooth, smoothMean, smoothCov);
            double[][] smoothCor = Statik.correlationFromCovariance(smoothCov);
            globalWarning = false;
            if (useFirstFactor) {
                smoothCor = computePairwiseFactorCorrelationSimple(smoothCor);
                if (globalWarning) System.out.print(partNr+"\t");
            }
            if (useTaskCombination) {
                smoothCor = computePairwiseCombinedCorrelationSimple(smoothCor);
                if (globalWarning) System.out.print(partNr+"\t");
            }
            if (!globalWarning) {
                smoothCorMem[partNr-1] = Statik.copy(smoothCor);
            }
            
            double[][] btwDataWithID = Statik.loadMatrix("FloData"+File.separator+betweenDataName+"_"+partNr+".dat", ' ', true, ".");
            allButOne = new int[anzVar]; for (int i=0; i<anzVar; i++) allButOne[i] = i+1;
            double[][] btwnData = Statik.submatrix(btwDataWithID, null, allButOne);

            pretestData[partNr-1] = btwnData;
            if (useTaskCombination) {
                pretestData[partNr-1] = new double[btwnData.length][9]; for (int i=0; i<btwnData.length; i++) {
                    double[] row = btwnData[i];
                    for (int j=0; j<9; j++) pretestData[partNr-1][i][j] = (Model.isMissing(row[2*j]) || Model.isMissing(row[2*j+1])?Model.MISSING:(row[2*j]+row[2*j+1])/2.0);
                }
            }
            pretestMean[partNr-1] = new double[corSize]; pretestCov[partNr-1] = new double[corSize][corSize];
            Statik.covarianceMatrixAndMeans(pretestData[partNr-1], pretestMean[partNr-1], pretestCov[partNr-1], Model.MISSING);
            globalWarning = false;
            double[][] btwCor = Statik.correlationFromCovariance(pretestCov[partNr-1]);
            if (!globalWarning) {
                if (!ignoreBetweenInMDS) corv.add(btwCor);
                if (!ignoreBetweenInMDS) id.add(partNr+1000);
                btwnAvgCov = Statik.add(btwnAvgCov, btwCor);
            } else System.out.print("*\t");
            
            avgPersonCorWthn[partNr] = 0; for (int i=0; i<wthnCor.length; i++) for (int j=i+1; j<wthnCor.length; j++) avgPersonCorWthn[partNr] += wthnCor[i][j]; 
            avgPersonCorWthn[partNr] /= (wthnCor.length*(wthnCor.length-1)/2);
            avgPersonCorBtw[partNr] = 0; for (int i=0; i<btwCor.length; i++) for (int j=i+1; j<btwCor.length; j++) avgPersonCorBtw[partNr] += btwCor[i][j]; 
            avgPersonCorBtw[partNr] /= (btwCor.length*(btwCor.length-1)/2);
            double[] eval = new double[corSize];
            double[][] evec = Statik.eigenvectorsOfSymmetrical(wthnCor, eval, true);
            highestEValCorWthn[partNr] = eval[0];
            highestEVecCorWthn[partNr] = Statik.multiply((evec[0][0]<0?-1:1), evec[0]);
            evec = Statik.eigenvectorsOfSymmetrical(btwCor, eval, true);
            highestEValCorBtwn[partNr] = eval[0];
            highestEVecCorBtwn[partNr] = Statik.multiply((evec[0][0]<0?-1:1), evec[0]);
            evec = Statik.eigenvectorsOfSymmetrical(smoothCor, eval, true);
            highestEValCorSmooth[partNr] = eval[0];
            highestEVecCorSmooth[partNr] = Statik.multiply((evec[0][0]<0?-1:1), evec[0]);
            if (useFirstFactor || useTaskCombination) {
                RAMModel flatFactor = floVeryFlatFactorModel();
                flatFactor.setDataDistribution(wthnCor, new double[9], 1);
                double[] starting = new double[18]; for (int i=0; i<9; i++) {
                    starting[i] = highestEVecCorWthn[partNr][i]; 
                    starting[i+9] = 1-Math.pow(highestEVecCorWthn[partNr][i],2);
                }
                double[] est = flatFactor.estimateML(starting); 
                flatFactorCoefficient[partNr] = new double[9]; for (int i=0; i<9; i++) flatFactorCoefficient[partNr][i] = est[i];
            }
            
            combinedData[partNr-1] = (useDetrend?wthnDataDetrend:wthnData);
            if (useTaskCombination) {
                combinedData[partNr-1] = new double[wthnData.length][9]; for (int i=0; i<wthnData.length; i++) {
                    double[] row = (useDetrend?wthnDataDetrend:wthnData)[i];
                    for (int j=0; j<9; j++) combinedData[partNr-1][i][j] = (Model.isMissing(row[2*j]) || Model.isMissing(row[2*j+1])?Model.MISSING:(row[2*j]+row[2*j+1])/2.0);
                }
            }
            Statik.covarianceMatrixAndMeans(combinedData[partNr-1], combinedMean[partNr-1], combinedWthnCov[partNr-1], Model.MISSING);
        } else nextOmit++;
        btwnAvgCov = Statik.multiply(1.0/anzPerFull, btwnAvgCov);
        wthnAvgCov = Statik.multiply(1.0/anzPerFull, wthnAvgCov);

        // Preparing computation of residual error when predicting a single value from covariance matrix. 
        double[] residualWthnBiasOnWthnData = new double[corSize], residualWthnR2OnWthnData = new double[corSize], 
                residualAvgWthnBiasOnWthnData = new double[corSize], residualAvgWthnR2OnWthnData = new double[corSize], 
                residualAvgWthnBiasOnBtwnData = new double[corSize], residualAvgWthnR2OnBtwnData = new double[corSize], 
                residualAvgBtwnBiasOnWthnData = new double[corSize], residualAvgBtwnR2OnWthnData = new double[corSize],
                residualAvgBtwnBiasOnBtwnData = new double[corSize], residualAvgBtwnR2OnBtwnData = new double[corSize], 
                residualZeroBiasOnWthnData = new double[corSize], residualZeroR2OnWthnData = new double[corSize],
                residualZeroBiasOnBtwnData = new double[corSize], residualZeroR2OnBtwnData = new double[corSize];
        
        double[] grandMean = new double[corSize]; double[][] combinedBtwnCov = new double[corSize][corSize];
        Statik.covarianceMatrixAndMeans(combinedMean, grandMean, combinedBtwnCov, Model.MISSING);
        double[][] combinedWthnAvgCov = new double[corSize][corSize];
        for (int partNr=1; partNr <= anzPerFull; partNr++) for (int i=0; i<corSize; i++) for (int j=0; j<corSize; j++) combinedWthnAvgCov[i][j] += combinedWthnCov[partNr-1][i][j];
        combinedWthnAvgCov = Statik.multiply(1.0/(double)anzPerFull, combinedWthnAvgCov);
        
        /*
         * The following loop computes all data prediction precision. There are 2x4 conditions: Prediction is either based on combined data, in which the two indicators
         * for each variable are collapsed, or on the full data. The matrix used is either the pretest covariance matrix (between condition), on the average within covariance matrix,
         * on the specific within covariance matrix of this participant, or on the identity matrix. 
         * 
         */
        int anzCond = 8; int[] totMeasures = new int[anzCond];
        double[][] residualBias = new double[anzCond][corSize], residualR2 = new double[anzCond][corSize];
        double[][] ident = Statik.identityMatrix(corSize);
        for (int cond=0; cond<anzCond; cond++) {
            for (int partNr=1; partNr <= anzPerFull; partNr++) {
                double[][] data = null, matrix = null;
                if (cond==0) {data = pretestData[partNr-1]; matrix = pretestCov[partNr-1];}
                if (cond==1) {data = pretestData[partNr-1]; matrix = combinedWthnAvgCov;}
                if (cond==2) {data = pretestData[partNr-1]; matrix = combinedWthnCov[partNr-1];}
                if (cond==3) {data = pretestData[partNr-1]; matrix = ident;}
                if (cond==4) {data = combinedData[partNr-1]; matrix = pretestCov[partNr-1];}
                if (cond==5) {data = combinedData[partNr-1]; matrix = combinedWthnAvgCov;}
                if (cond==6) {data = combinedData[partNr-1]; matrix = combinedWthnCov[partNr-1];}
                if (cond==7) {data = combinedData[partNr-1]; matrix = ident;}
                double[] mean = Statik.meanVector(data,Model.MISSING);
                
                for (int t=0; t<data.length; t++) {
                    boolean hasMissing = false; for (int i=0; i<corSize; i++) if (Model.isMissing(data[t][i])) hasMissing=true;
                    if (!hasMissing) {
                        double[] datarow = Statik.subtract(data[t], mean); 
                        SaturatedRAMModel sat = new SaturatedRAMModel(corSize, true);
                        sat.filter = new int[corSize-1];
                        sat.setParameterToDistribution(new double[corSize], matrix);
                        for (int i=0; i<corSize; i++) {
                            for (int j=0; j<corSize-1; j++) sat.filter[j] = (j<i?j:j+1);
                            sat.setFilter(sat.filter);
                            double[] scores = sat.getAllScores(null, new double[][]{Statik.subvector(datarow, i)})[0];
                            double residual = scores[i]-datarow[i];
                            residualBias[cond][i] += residual; residualR2[cond][i] += residual*residual;
                        }
                        totMeasures[cond]++;
                    }
                }
            }
        }
        System.out.println("\r\n\r\nPrediction Precision bias and R^2");
        for (int cond=0; cond<anzCond; cond++) {
            residualBias[cond] = Statik.multiply(1.0/(double)(totMeasures[cond]), residualBias[cond]);
            residualR2[cond] = Statik.multiply(1.0/(double)(totMeasures[cond]), residualR2[cond]);
        }
        for (int cond=0; cond<anzCond; cond++) System.out.println("Bias Condition "+cond+" : \t"+Statik.matrixToString(residualBias[cond]));
        for (int cond=0; cond<anzCond; cond++) System.out.println("R^2 Condition "+cond+"  : \t"+Statik.matrixToString(residualR2[cond]));
        // END prediction precision
        
        /*
         * The following block prepares the Kullback Leibler Divergences between all "between" and "within" covariance matrices and establishes by simulation the null distributions. 
         */
        System.out.println();
        int anzPer = corv.size();
        double[][][] cor = new double[anzPer][][];
        corv.toArray(cor);
        Integer[] idI = new Integer[anzPer]; id.toArray(idI);
        int[] ids = new int[anzPer]; for (int i=0; i<anzPer; i++) ids[i] = idI[i].intValue();
        
        int anzBtw = 0; for (int i=0; i<anzPer; i++) if (ids[i] > 1000) anzBtw++;
        double[][][] allBtwCov = new double[anzBtw][][], block1BtwCov = new double[anzBtw][][], 
                block2BtwCov = new double[anzBtw][][], block3BtwCov = new double[anzBtw][][]; 
        int k = 0; for (int i=0; i<anzPer; i++) if (ids[i] > 1000) {allBtwCov[k] = cor[i]; block1BtwCov[k] = Statik.submatrix(cor[i], new int[]{0,1,2}, new int[]{0,1,2});
           block2BtwCov[k] = Statik.submatrix(cor[i], new int[]{3,4,5}, new int[]{3,4,5}); block3BtwCov[k] = Statik.submatrix(cor[i], new int[]{6,7,8}, new int[]{6,7,8});
           k++;
        }
        int simTrials = 100;
        System.out.println("********** Starting simulation of null distribution for average KL (within) ************");
        floSimulateAverageKLUnderNull(wthnAvgCov, null, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (between) ************");
        floSimulateAverageKLUnderNull(btwnAvgCov, allBtwCov, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (PS within) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(wthnAvgCov, new int[]{0,1,2}, new int[]{0,1,2}), null, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (PS between) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(btwnAvgCov, new int[]{0,1,2}, new int[]{0,1,2}), block1BtwCov, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (WM within) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(wthnAvgCov, new int[]{3,4,5}, new int[]{3,4,5}), null, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (WM between) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(btwnAvgCov, new int[]{3,4,5}, new int[]{3,4,5}), block2BtwCov, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (EM within) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(wthnAvgCov, new int[]{6,7,8}, new int[]{6,7,8}), null, simTrials);
        System.out.println("********** Starting simulation of null distribution for average KL (EM between) ************");
        floSimulateAverageKLUnderNull(Statik.submatrix(btwnAvgCov, new int[]{6,7,8}, new int[]{6,7,8}), block3BtwCov, simTrials);

        floErgodicSubspace(btwnAvgCov, wthnAvgCov);
        
        if (!ignoreBetweenInMDS) {
            System.out.println("Distances of within covariance matrices to between centroid (5% significant greater than 62.70), then distance to corresponding between: ");
            for (int i=0; i<anzPer; i++) if (ids[i] < 1000) {
                double distToBtwCentroid = Statik.symmetricalKullbackLeiblerNormal(cor[i], btwnAvgCov);
                double distToCorBtw = Statik.symmetricalKullbackLeiblerNormal(cor[i], cor[i+1]);
                System.out.println(ids[i]+"\t"+distToBtwCentroid+"\t"+distToCorBtw);
            }
        }
        
        double[][] distance = new double[anzPer][anzPer];
        for (int i=0; i<anzPer; i++) for (int j=i; j<anzPer; j++) distance[i][j] = distance[j][i] = Statik.symmetricalKullbackLeiblerNormal(cor[i], cor[j]);
        // following two lines just to iron out numerical inprecisions on the diagonal
        for (int i=0; i<anzPer; i++) for (int j=i+1; j<anzPer; j++) distance[i][j] = distance[j][i] = distance[i][j] - (distance[i][i] + distance[j][j])/2.0;
        for (int i=0; i<anzPer; i++) distance[i][i] = 0.0; 

        double[][] btwToWthnDist = new double[4][anzPer];
        int totWthn = 0, totBtwn = 0;
        // averageKLWithin is the average KL of all within covariance matrix. averageKLBetween is the average KL from with to between for each participant.
        double averageKLWithin = 0, averageKLBetween = 0, stdvKLBetween = 0; 
        for (int i=0; i<anzPer; i++) for (int j=0; j<anzPer; j++) 
            if (ids[i]<1000 && ids[j]<1000)        {totWthn++; averageKLWithin += distance[i][j];}
            else if (ids[i]<1000 && ids[j]==ids[i]+1000) {totBtwn++; averageKLBetween += distance[i][j]; 
                                                          stdvKLBetween += distance[i][j]*distance[i][j]; btwToWthnDist[0][i] = distance[i][j];}
        averageKLWithin /= (double)totWthn; averageKLBetween /= (double)totBtwn;
        stdvKLBetween = Math.sqrt( stdvKLBetween/(double)totBtwn - averageKLBetween*averageKLBetween);
        System.out.println("Average KL within: "+averageKLWithin);
        if (!ignoreBetweenInMDS) System.out.println("Average KL within to between: "+averageKLBetween+" +- "+stdvKLBetween);
        
        for (int block = 0; block < 3; block++) {
            int[] selector = new int[]{3*block, 3*block+1, 3*block+2};
            totWthn = 0; totBtwn = 0;
            averageKLWithin = 0; averageKLBetween = 0; stdvKLBetween = 0; 
            for (int i=0; i<anzPer; i++) for (int j=0; j<anzPer; j++) {
                double value = Statik.symmetricalKullbackLeiblerNormal(Statik.submatrix(cor[i],selector,selector), Statik.submatrix(cor[j],selector,selector));
                if (ids[i]<1000 && ids[j]<1000)        {totWthn++; averageKLWithin += value;}
                else if (ids[i]<1000 && ids[j]==ids[i]+1000) {totBtwn++; averageKLBetween += value; 
                                                              stdvKLBetween += value*value; btwToWthnDist[1+block][i] = value;}
            }
            averageKLWithin /= (double)totWthn; averageKLBetween /= (double)totBtwn; 
            stdvKLBetween = Math.sqrt( stdvKLBetween/(double)totBtwn - averageKLBetween*averageKLBetween );
            System.out.println("Average KL within (Block "+block+"): "+averageKLWithin);
            if (!ignoreBetweenInMDS) System.out.println("Average KL within to between (Block "+block+"): "+averageKLBetween+" +- "+stdvKLBetween);
        }
        double[][] btwToWthnDistPruned = new double[totBtwn][4];
        k = 0; for (int i=0; i<anzPer; i++) if (btwToWthnDist[0][i] != 0.0) {for (int j=0; j<4; j++) btwToWthnDistPruned[k][j] = btwToWthnDist[j][i]; k++;} 
        System.out.println("All between/within distances: \r\n All\tPS\tWM\tEM\r\n"+Statik.matrixToString(btwToWthnDistPruned));

        /*
         * The following block starts the multidimensional scaling, using the standard algorithm implemented in Onyx. 
         */
        int dim = 5;
        double[] scree = new double[dim];
        double[][] coordinates = DataPreprocessing.multiDimensionalScaling(distance, scree, dim, 0.001);
        Statik.writeMatrix(coordinates, "toFlo80DimMDS.csv",'\t');
        
        /*
         *  The following block computes the agglomerative clustering, using the standard algorithm implemented in Onyx. 
         */
        AgglomerativeClustering cluster = new AgglomerativeClustering(distance, true);
        int[] membership = cluster.cluster();

        double[][] coordID = new double[anzPer][dim+2]; 
        for (int i=0; i<anzPer; i++) {
            coordID[i][0] = id.elementAt(i); 
            coordID[i][1] = membership[i];
            for (int j=0; j<dim; j++) coordID[i][j+2] = coordinates[i][j];
        } 
        System.out.println("\t\t"+Statik.matrixToString(scree)+"\r\n");
        System.out.println(Statik.matrixToString(coordID));

        System.out.println("Within Cov = \r\n"+Statik.matrixToString(wthnAvgCov));
        System.out.println("Between Cov = \r\n"+Statik.matrixToString(btwnAvgCov));

    } // END script for Schmiedeck et al. manuscript