Merge pull request #185 from bbfrederick/refinedelay

Refinedelay

.gitignore

@@ -19,7 +19,7 @@ v2changes
 
 # any generated test data
 rapidtide/data/examples/dst/*
-rapidtide/tests/tmp/[abcdefghijklmnopqrstuvwxyz]*
+rapidtide/tests/tmp/[0123456789abcdefghijklmnopqrstuvwxyz]*
 rapidtide/tests/*.jpg
 rapidtide/tests/*.json
 rapidtide/tests/*.tsv.gz

CHANGELOG.md

@@ -1,5 +1,22 @@
 # Release history
 
+## Version 3.0alpha1 (12/20/24)
+* (rapidtide) The ``--fixdelay`` option has been split into two options.  ``--initialdelay DELAY`` lets you specify either a float that sets the starting delay for every voxel to that value, or a 3D file specifying the initial delay for each voxel.  ``--nodelayfit`` determines whether the delay can be adjusted from its initial value.  Closes https://github.com/bbfrederick/rapidtide/issues/171.
+* (rapidtide) Reorganized command line options and adjusted the default values.
+* (rapidtide) Help output now shows the filter ranges.
+* (rapidtide, retroglm) Added delay refinement using the ratio of the fit coefficients of the regressor and its time derivative.
+* (rapidtide, retroglm) Fixed a bad shared memory leak.
+* (retroglm) Significantly enhanced logging.
+* (retroglm) Added canary files.
+* (rapidtide) Implemented delay map patching.
+* (rapidtide) Write out individual EV timecourses.
+* (Docker) Cleaned up some internal variables.
+* (Docker) Improved build and testing scripts.
+* (io) Added function to compare nifti files with some tolerance.
+* (docs) Automated more table generation.
+* (package) Merged some dependabot PRs.
+* (package) Fixed a fairly big, but not necessarily impactful bug. mlregress returned R2, not R, so anything referring to the R of a fit was actually squared (R was actually R2, R2 was actually R4).
+
 ## Version 2.9.9.5 (11/15/24)
 * (deployment) New idea - split the README.rst file to remove potentially offending reference stuff.
 

USAGE.rst

@@ -11,7 +11,7 @@ a) Run rapidtide to perform dynamic global signal regression (dGSR) on an fMRI f
 b) Run rapidtide to perform static global signal regression (sGSR) on an fMRI file[1] (this is just global mean regression):
 ::
 
-    rapidtide rapidtide/data/examples/src/sub-RAPIDTIDETEST rapidtide/data/examples/dst/sub-RAPIDTIDETEST_sgsr --fixdelay 0.0 --passes 1
+    rapidtide rapidtide/data/examples/src/sub-RAPIDTIDETEST rapidtide/data/examples/dst/sub-RAPIDTIDETEST_sgsr --nodelayfit --passes 1
 
 
 c) Run tidepool to look at all the interesting maps and timecourses from a):

docs/theoryofoperation.rst

@@ -723,6 +723,45 @@ Finally, if you don't want to do glm filtering at all
 you can shut off the glm filtering with ``--noglm``.
 
 
+Delay refinement
+^^^^^^^^^^^^^^^^
+
+This is new to rapidtide 3.0.  I've added a new method for refining the delay time estimate in
+every voxel based on the filtering step.  To the best of my knowledge, this is something I came up
+with (well not entirely, but this application).
+
+As we remember from freshman physics, you can extrapolate a signal using a Taylor series approximation.
+Which is to say, if you know the value of a function at a time t, and the value of the derivative,
+and the second
+derivative, and so on, you can calculate the signal at another point t + delta t by using a
+weighted sum of those  values.  Neat!  Even more neat is that for sufficiently small values of
+delta t, you can get a pretty good approximation using just the function and it's first derivative.
+
+As always, there are some complications:
+
+    * The mapping between fit coefficient ratio and time delay depends on the function, so it needs
+      to be determined for each regressor.  It's linear for very small delay value, and then the
+      mapping diverges (in a regressor specific way) as the delay increases.
+
+    * As I mentioned, this only works for "small" delay times.  What is small?  For LFO signals
+      in the 0.01 to 0.15 Hz band, this is only really good for about +/-3-5 seconds of offset
+      (and the linear region is only about +/-0.75 seconds (which is why we can't use this method
+      for the initial delay estimation, only for tuning).  The mapping function ends up being
+      sigmoid - you can't really calculate the delay from the ratio when the slope gets close to zero.
+      When that happens depends on the specific regressor, but you can pretty much always do the
+      mapping out to about +/-3.5 seconds.
+
+What is this good for?  Well, one thing I have found is that rapidtide gets much better fits if you
+use a fairly strong spatial smoothing filter (5mm gaussian kernel).  That's great for getting rid
+of a lot of the annoying speckling in the delay maps, but the result is that you lose a lot of fine
+detail in the delay map (which is obvious when you think about it). BUT - we know that delay varies
+relatively smoothly in real brains, so the smoothed delay values, while maybe not exactly right in
+most voxels, aren't far off.  So the delay in any voxel will be within +/-3 seconds of the smoothed
+value in every voxel, so the ratio-of-fit-derivatives method will be able to fit the difference, which
+you can then apply as an offset to find the exact delay in every voxel with much higher spatial resolution.
+Neat, huh?
+
+
 References
 ^^^^^^^^^^
 

docs/usage_rapidtide.rst

@@ -74,8 +74,11 @@ BIDS Outputs:
    :header: "Name", "Extension(s)", "Content", "When present"
    :widths: 30, 10, 30, 20
 
-    "XXX_DONE", "txt", "Run status file", "If rapidtide successfully completed"
-    "XXX_RUNNING", "txt", "Run status file", "If rapidtide was started, but has not yet successfully completed"
+    "XXX_DONE", "txt", "Run status file", "Always if rapidtide has finished without error"
+    "XXX_ISRUNNING", "txt", "Run status file", "Only present if rapidtide is running, or has failed"
+    "XXX_formattedcommandline", "txt", "Command line used to invoke rapidtide, nicely formatted", "Always"
+    "XXX_log", "txt", "Diagnostic log file with a lot of informational output", "Always"
+    "XXX_memusage", "tsv", "Memory usage statistics for performance tuning", "Always"
     "XXX_commandline", "txt", "Raw command line used to invoke rapidtide", "Always"
     "XXX_desc-autocorr_timeseries", "tsv.gz, json", "Autocorrelation of the probe regressor for each pass", "Always"
     "XXX_desc-cleansimdistdata_info", "tsv.gz, json", "Individual sham correlation datapoints after outlier removal", "Present if numnull > 0"
@@ -86,15 +89,21 @@ BIDS Outputs:
     "XXX_desc-corrfitfailreason_info", "nii.gz, json", "Result codes for correlation fit", "Always"
     "XXX_desc-corrfitwindow_info", "nii.gz, json", "The search window for the correlation peak fit", "Present if outputlevel is max"
     "XXX_desc-corrout_info", "nii.gz, json", "Correlation function", "Present if outputlevel is ≥ normal"
-    "XXX_desc-corrtimes_timeseries", "tsv.gz, json", "", "Present if outputlevel is max"
+    "XXX_desc-corrtimes_timeseries", "tsv.gz, json", "Correlation time axis", "Present if outputlevel is max"
+    "XXX_desc-CoV_map", "nii.gz, json", "Voxelwise coefficient of variation of fmri data", "Always"
+    "XXX_desc-delayoffset_hist", "tsv.gz, json", "Histogram of delay offsets calculated from GLM", "Present if refinedelay is enabled"
+    "XXX_desc-delayoffset_map", "nii.gz, json", "Delay offset correction from delay refinement", "Present if refinedelay is enabled"
     "XXX_desc-despeckle_mask", "nii.gz, json", "Voxels that underwent despeckling in the final pass", "Present if despecklepasses > 0 (default) and outputlevel is max"
+    "XXX_desc-EV_timeseries", "tsv.gz, json", "GLM regressor set", "Present if GLM is enabled (default)"
     "XXX_desc-expandedconfounds_timeseries", "tsv.gz, json", "The expanded (via derivatives and powers) set of confound regressors used for prefiltering the data", "Present if doing motion/confound regression"
+    "XXX_desc-filteredglmderivratios_map", "nii.gz, json", "glmderivratios, with outliers patched using median filtered data", "Present if refinedelay is enabled and outputlevel is ≥ normal"
     "XXX_desc-formattedruntimings_info", "tsv", "No description", "Always"
     "XXX_desc-gaussout_info", "nii.gz, json", "Simulated correlation function", "Present if outputlevel is max"
+    "XXX_desc-glmderivratios_map", "nii.gz, json", "Ratio of the first derivative of delayed sLFO to the delayed sLFO", "Present if refinedelay is enabled and outputlevel is ≥ normal"
     "XXX_desc-globallag_hist", "tsv.gz, json", "Histogram of lag times from global lag calculation", "Always"
     "XXX_desc-globalmean_mask", "nii.gz, json", "Voxels used to calculate global mean", "Always"
     "XXX_desc-initialmovingregressor_timeseries", "tsv.gz, json", "The raw and filtered initial probe regressor, at the original sampling resolution", "Always"
-    "XXX_desc-lagtcgenerator_timeseries", "tsv.gz, json", "The lagged timecourse generator", "Present if passes > 1"
+    "XXX_desc-lagtcgenerator_timeseries", "tsv.gz, json", "The lagged timecourse generator", "Always"
     "XXX_desc-lfofilterCleaned_bold", "nii.gz, json", "fMRI data with sLFO signal filtered out", "Present if GLM is enabled (default) and outputlevel is ≥ less"
     "XXX_desc-lfofilterCoeff_map", "nii.gz, json", "Fit coefficient", "Present if GLM is enabled (default) and outputlevel is ≥ normal"
     "XXX_desc-lfofilterCoeffDerivN_map", "nii.gz, json", "Fit coefficient for the Nth temporal derivative", "Present if GLM is enabled (default), glmderivs > 0, and outputlevel is max"
@@ -114,13 +123,16 @@ BIDS Outputs:
     "XXX_desc-lfofilterRemoved_bold", "nii.gz, json", "sLFO signal filtered out of this voxel", "Present if GLM is enabled (default) and outputlevel is ≥ more"
     "XXX_desc-maxcorr_hist", "tsv.gz, json", "Histogram of maximum correlation coefficients", "Always"
     "XXX_desc-maxcorr_map", "nii.gz, json", "Maximum correlation strength", "Always"
+    "XXX_desc-maxcorralt_map", "nii.gz, json", "R value of the GLM fit, with sign", "Present if refinedelay is enabled"
     "XXX_desc-maxcorrsq_map", "nii.gz, json", "Squared maximum correlation strength (proportion of variance explained)", "Always"
     "XXX_desc-maxtime_hist", "tsv.gz, json", "Histogram of maximum correlation times", "Always"
     "XXX_desc-maxtime_map", "nii.gz, json", "Lag time in seconds", "Always"
+    "XXX_desc-maxtimerefined_map", "nii.gz, json", "Lag time in seconds, refined", "Present if refinedelay is enabled"
     "XXX_desc-maxwidth_hist", "tsv.gz, json", "Histogram of correlation peak widths", "Always"
     "XXX_desc-maxwidth_map", "nii.gz, json", "Width of corrrelation peak", "Always"
     "XXX_desc-mean_map", "nii.gz, json", "Voxelwise mean of fmri data", "Always"
-    "XXX_desc-mitimes_timeseries", "tsv.gz, json", "", "Present if outputlevel is max"
+    "XXX_desc-medfiltglmderivratios_map", "nii.gz, json", "Median filtered version of the glmderivratios map", "Present if refinedelay is enabled and outputlevel is ≥ normal"
+    "XXX_desc-mitimes_timeseries", "tsv.gz, json", "Cross mutual information time axis", "Present if outputlevel is max"
     "XXX_desc-movingregressor_timeseries", "tsv.gz, json", "The probe regressor used in each pass, at the time resolution of the data", "Always"
     "XXX_desc-MTT_hist", "tsv.gz, json", "Histogram of correlation peak widths", "Always"
     "XXX_desc-MTT_map", "nii.gz, json", "Mean transit time (estimated)", "Always"
@@ -132,17 +144,17 @@ BIDS Outputs:
     "XXX_desc-plt0p050_mask", "nii.gz, json", "Voxels where the maxcorr value exceeds the p < 0.050 significance level", "Present if numnull > 0"
     "XXX_desc-preprocessedconfounds_timeseries", "tsv.gz, json", "The preprocessed (normalized, filtered, orthogonalized) set of expanded confound regressors used for prefiltering the data", "Present if doing motion/confound regression"
     "XXX_desc-processed_mask", "nii.gz", "No description", "Always"
+    "XXX_desc-ratiotodelayfunc_timeseries", "tsv.gz, json", "The function mapping derivative ratio to delay", "Present if refinedelay is enabled"
     "XXX_desc-refine_mask", "nii.gz, json", "Voxels used for refinement", "Present if passes > 1"
     "XXX_desc-refinedmovingregressor_timeseries", "tsv.gz, json", "The raw and filtered probe regressor produced by the refinement procedure, at the time resolution of the data", "Present if passes > 1"
     "XXX_desc-runoptions_info", "json", "A detailed dump of all internal variables in the program.  Useful for debugging and data provenance.", "Always"
     "XXX_desc-shiftedtcs_bold", "nii.gz, json", "The filtered input fMRI data, in voxels used for refinement, time shifted by the negated delay in every voxel so that the moving blood component is aligned.", "Present if passes > 1 and outputlevel is max"
     "XXX_desc-simdistdata_info", "tsv.gz, json", "Individual sham correlation datapoints", "Present if numnull > 0"
+    "XXX_desc-sLFOamplitude_timeseries", "tsv.gz, json", "Filtered RMS amplitude of the probe regressor, and a linear fit", "Always"
+    "XXX_desc-std_map", "nii.gz, json", "Voxelwise standard deviation of fmri data", "Always"
     "XXX_desc-timepercentile_map", "nii.gz, json", "Percentile ranking of this voxels delay", "Always"
-    "XXX_desc-trimmedcorrtimes_timeseries", "tsv.gz, json", "", "Present if outputlevel is max"
-    "XXX_desc-trimmedmitimes_timeseries", "tsv.gz, json", "", "Present if outputlevel is max"
-    "XXX_formattedcommandline", "txt", "Command line used to invoke rapidtide, nicely formatted", "Always"
-    "XXX_log", "txt", "Diagnostic log file with a lot of informational output", "Always"
-    "XXX_memusage", "tsv", "Memory usage statistics for performance tuning", "Always"
+    "XXX_desc-trimmedcorrtimes_timeseries", "tsv.gz, json", "Trimmed correlation time axis", "Present if outputlevel is max"
+    "XXX_desc-trimmedmitimes_timeseries", "tsv.gz, json", "Trimmed cross mutual information time axis", "Present if outputlevel is max"
 ..
 
 
@@ -174,22 +186,37 @@ CORRFUNCSIZE is the size of the correlation function in TRs at the oversampled T
 The output sizes in TRs (with no motion regression) are as follows:
 
 .. csv-table::  Total image output data size in TRs
-   :header: "Output level", "GLM?", "Number of TRs"
-   :widths: 10, 10, 10
-
-    "min", "No", "16"
-    "less", "No", "16"
-    "normal", "No", "16 + CORRFUNCSIZE"
-    "more", "No", "16 + CORRFUNCSIZE"
-    "max", "No", "17 + CORRFUNCSIZE"
-    "min", "Yes", "24"
-    "less", "Yes", "24 + FMRISIZE"
-    "normal", "Yes", "24 + CORRFUNCSIZE + FMRISIZE"
-    "more", "Yes", "24 + CORRFUNCSIZE + 3*FMRISIZE"
-    "max", "Yes", "25 + 3*CORRFUNCSIZE + 4*FMRISIZE"
+   :header: "Output level", "Passes>1?", "Refine delay?", "GLM?", "Number of TRs"
+   :widths: 10, 10, 10, 10, 20
+
+    "min", "No", "No",  "No", "13"
+    "min", "No", "Yes",  "No", "16"
+    "min", "Yes", "No",  "No", "16"
+    "min", "Yes", "Yes",  "No", "19"
+    "min", "No", "No",  "Yes", "14"
+    "min", "Yes", "No",  "Yes", "17"
+    "less", "No", "No",  "No", "13"
+    "less", "No", "Yes",  "No", "17 + 1*FMRISIZE"
+    "less", "Yes", "No",  "No", "16"
+    "less", "Yes", "Yes",  "No", "20 + 1*FMRISIZE"
+    "normal", "No", "No",  "No", "13 + 1*CORRFUNCSIZE"
+    "normal", "No", "Yes",  "No", "21 + 1*CORRFUNCSIZE + 1*FMRISIZE"
+    "normal", "Yes", "No",  "No", "19 + 1*CORRFUNCSIZE"
+    "normal", "Yes", "Yes",  "No", "27 + 1*CORRFUNCSIZE + 1*FMRISIZE"
+    "more", "No", "No",  "No", "13 + 1*CORRFUNCSIZE"
+    "more", "No", "Yes",  "No", "21 + 1*CORRFUNCSIZE + 2*FMRISIZE"
+    "more", "Yes", "No",  "No", "19 + 1*CORRFUNCSIZE"
+    "more", "Yes", "Yes",  "No", "27 + 1*CORRFUNCSIZE + 2*FMRISIZE"
+    "max", "No", "No",  "No", "13 + 3*CORRFUNCSIZE"
+    "max", "No", "Yes",  "No", "21 + 3*CORRFUNCSIZE + 3*FMRISIZE"
+    "max", "Yes", "No",  "No", "19 + 3*CORRFUNCSIZE"
+    "max", "Yes", "Yes",  "No", "27 + 3*CORRFUNCSIZE + 3*FMRISIZE"
+    "max", "No", "No",  "Yes", "14 + 3*CORRFUNCSIZE + 1*FMRISIZE"
+    "max", "Yes", "No",  "Yes", "20 + 3*CORRFUNCSIZE + 1*FMRISIZE"
 ..
 
-The data size is then this number of TRs times the size of 1 TR worth of data in the input fMRI file.
+The data size is then this number of TRs times the size of 1 TR worth of data in the input fMRI file, (plus the size
+of the various timecourse files and .json sidecars which are much smaller than the image files).
 
 
 As an example, the following table shows the size of the data produced by running a rapidtide analysis on one HCP-YA

pyproject.toml

@@ -43,7 +43,7 @@ dependencies = [
 #dynamic = ['version', 'license', 'keywords']
 dynamic = ['version']
 authors = [
-    {name = "Blaise deB Frederick", email='blaise.frederick@gmail.com' },
+    {name = "Blaise deB. Frederick", email="blaise.frederick@gmail.com"},
     {name = "Taylor Salo"},
     {name = "Daniel M. Drucker, Ph.D."},
     {name = "Jeffrey N Stout"},

rapidtide/RapidtideDataset.py

@@ -1051,6 +1051,8 @@ def setupoverlays(self):
         if self.bidsformat:
             self.funcmaps = [
                 ["lagtimes", "desc-maxtime_map"],
+                ["lagtimesrefined", "desc-maxtimerefined_map"],
+                ["delayoffset", "desc-delayoffset_map"],
                 ["timepercentile", "desc-timepercentile_map"],
                 ["lagstrengths", "desc-maxcorr_map"],
                 ["lagsigma", "desc-maxwidth_map"],

rapidtide/calcnullsimfunc.py

@@ -38,7 +38,7 @@ def _procOneNullCorrelationx(
     thefitter,
     despeckle_thresh=5.0,
     fixdelay=False,
-    fixeddelayvalue=0.0,
+    initialdelayvalue=0.0,
     permutationmethod="shuffle",
     disablethresholds=False,
     rt_floatset=np.float64,
@@ -81,7 +81,7 @@ def getNullDistributionDatax(
     thefitter,
     despeckle_thresh=5.0,
     fixdelay=False,
-    fixeddelayvalue=0.0,
+    initialdelayvalue=0.0,
     numestreps=0,
     nprocs=1,
     alwaysmultiproc=False,
@@ -153,7 +153,7 @@ def nullCorrelation_consumer(inQ, outQ):
                             thefitter,
                             despeckle_thresh=despeckle_thresh,
                             fixdelay=fixdelay,
-                            fixeddelayvalue=fixeddelayvalue,
+                            initialdelayvalue=initialdelayvalue,
                             permutationmethod=permutationmethod,
                             rt_floatset=rt_floatset,
                             rt_floattype=rt_floattype,
@@ -205,7 +205,7 @@ def nullCorrelation_consumer(inQ, outQ):
                 thefitter,
                 despeckle_thresh=despeckle_thresh,
                 fixdelay=fixdelay,
-                fixeddelayvalue=fixeddelayvalue,
+                initialdelayvalue=initialdelayvalue,
                 permutationmethod=permutationmethod,
                 rt_floatset=rt_floatset,
                 rt_floattype=rt_floattype,