From 64984b452919e6d77f12ca2fa20b72678a77e3b6 Mon Sep 17 00:00:00 2001
From: Douglas Greve <greve@nmr.mgh.harvard.edu>
Date: Mon, 21 Feb 2022 18:11:25 -0500
Subject: [PATCH] mri_glmfit. #NF. Can now handle perm sim with PVRs

---
 mri_glmfit/mri_glmfit.cpp | 48 ++++++++++++++++++++++++++++++++++++---
 scripts/mri_glmfit-sim    |  5 ++++
 2 files changed, 50 insertions(+), 3 deletions(-)

diff --git a/mri_glmfit/mri_glmfit.cpp b/mri_glmfit/mri_glmfit.cpp
index 5c67eecda33..a65583a1f15 100644
--- a/mri_glmfit/mri_glmfit.cpp
+++ b/mri_glmfit/mri_glmfit.cpp
@@ -548,8 +548,10 @@ double round(double x);
 #include "dti.h"
 #include "image.h"
 #include "stats.h"
+#include "evschutils.h"
 
 int MRISmaskByLabel(MRI *y, MRIS *surf, LABEL *lb, int invflag);
+int RandPermMatrixAndPVR(MATRIX *X, MRI **pvrs, int npvrs);
 
 static int  parse_commandline(int argc, char **argv);
 static void check_options(void);
@@ -1738,11 +1740,11 @@ int main(int argc, char **argv) {
           SmoothSurfOrVol(surf, mriglm->y, mriglm->mask, SmoothLevel);
       }
       if (!strcmp(csd->simtype,"perm")) {
-        if (!OneSamplePerm) MatrixRandPermRows(mriglm->Xg);
+        if (!OneSamplePerm) RandPermMatrixAndPVR(mriglm->Xg,mriglm->pvr,mriglm->npvr);
         else {
           for (n=0; n < mriglm->y->nframes; n++) {
             if (drand48() > 0.5) m = +1;
-            else                m = -1;
+            else                 m = -1;
             mriglm->Xg->rptr[n+1][1] = m;
           }
           //MatrixPrint(stdout,mriglm->Xg);
@@ -3645,7 +3647,8 @@ static void check_options(void) {
     exit(1);
   }
 
-  if(DoSim && !strcmp(csd->simtype,"perm") && npvr != 0){
+  if(0 && DoSim && !strcmp(csd->simtype,"perm") && npvr != 0){
+    // Modified to allow for pvrs with sim
     printf("ERROR: PVR is not supported with permutation simulations\n");
     exit(1);
   }
@@ -4140,6 +4143,45 @@ int MRIloganize(MATRIX **X, MRI **Ct, MRI **intCt, const MATRIX *t, const double
   return(0);
 }
 
+/*!
+  \fn int RandPermMatrixAndPVR(MATRIX *X, MRI **pvrs, int npvrs)
+  \brief Permutes both the design matrix and any PVRs
+ */
+int RandPermMatrixAndPVR(MATRIX *X, MRI **pvrs, int npvrs)
+{
+  int *NewRowOrder, r;
+  MATRIX *X0;
+
+  NewRowOrder = RandPerm(X->rows, NULL);
+  for (r = 0; r < X->rows; r++) {
+    NewRowOrder[r]++;  // Make one-based
+    printf("%d\n",NewRowOrder[r]);
+  }
+  X0 = MatrixCopy(X, NULL);
+  MatrixReorderRows(X0, NewRowOrder, X);
+
+  int n;
+  for(n=0; n < npvrs; n++){
+    int c,r,f;
+    MRI *pvr = pvrs[n];
+    for(c=0; c < pvr->width; c++){
+      double *vect = (double *) calloc(sizeof(double),X->rows);
+      for(r=0; r < pvr->height; r++){
+	for(s=0; s < pvr->depth; s++){
+	  for(f=0; f < X->rows; f++) vect[f] = MRIgetVoxVal(pvr,c,r,s,f);
+	  for(f=0; f < X->rows; f++) {
+	    int f2 = NewRowOrder[f] - 1;
+	    MRIsetVoxVal(pvr,c,r,s,f, vect[f2]);
+	  }
+	}
+      }
+      free(vect);
+    }
+  }
+
+  free(NewRowOrder);
+  return (0);
+}
 
 
 
diff --git a/scripts/mri_glmfit-sim b/scripts/mri_glmfit-sim
index 6414915a375..2f45174dd8b 100755
--- a/scripts/mri_glmfit-sim
+++ b/scripts/mri_glmfit-sim
@@ -123,6 +123,7 @@ echo $glmfitcmd
 set gd2mtx = dods
 set UseTable = 0;
 set label = ();
+set pvrlist = ()
 while($#glmfitcmd)
   set flag = $glmfitcmd[1]; shift glmfitcmd;
   #echo $flag
@@ -250,6 +251,7 @@ while($#glmfitcmd)
     breaksw
   case "--pvr"
     set IsPVR = 1;
+    set pvrlist = ($pvrlist $glmfitcmd[1]);
     shift glmfitcmd;
     breaksw
   # Ignore these options that have two args
@@ -457,6 +459,9 @@ if($DoSim) then
     if($#X)    set cmd = ($cmd --X $X)
     if($#wls)  set cmd = ($cmd --wls $wls)
     if($#framemask)  set cmd = ($cmd --frame-mask $framemask)
+    foreach pvr ($pvrlist) 
+      set cmd = ($cmd --pvr $pvr)
+    end
     # --y must go after --fsgd
     if(! $UseTable) then
       set cmd = ($cmd --mask $mask)