Hi all,

I recently downloaded the latest version of PESTICA (v1.3 for AFNI users) and it has been giving me trouble.  Eventually I tracked the problem down to the fact that line 302 of run_pestica.sh deletes the dataset irfretroicor+orig.*, on which irf_retroicor.m has labored for the past few hours.  I thought others might want to know about this.  Also, the same thing appears at line 265 immediately after a call to matlab_retroicor.m, although in that case it's again irfretroicor+orig.* that's getting deleted and not the dataset that's just been created.

I've also noticed that my output IRFs don't really look like the ones in the tutorial, although they seem reasonably consistent across the few subjects I've gotten through the pipeline.  I've attached an example plot, though I suspect that if there is a problem, this won't be enough information to diagnose it.  Is this plot out of the ordinary?

Thanks,
Isaac Schwabacher

Hi Isaac,

there is actually a second pass run of irf_retroicor() that writes out irfretroicor+orig again.  Below, I've listed the code in run_pestica.sh that follows the second pass of irf_retroicor().  It runs 3dAFNItoNIFTI (to "fix" the written file, otherwise there is a file formatting problem that exists I think due to older versions of the AFNI_matlab library), then 3dcopy to give it an appropriate name based on $fname.  Unfortunately, if this later code fails, you do end up without the results of the several hours run of fit_each_physio_peak() (fit_each_physio_peak is what takes a large amount of time).  A top priority is making this more robust to failure in the next release, and reduce the time taken by fit_each_physio_peak.m.  I think the failure is related to bash, but haven't confirmed it, as I don't get this failure any more.  Can you send me your shell output (note redirecting shell output causes MATLAB to halt, so you might have to copy and paste)?  Sometimes the script just _stops_ after running a stage, and I have not figured out why it doesn't continue, which must be due to bash.  The next version will be structured differently.  Again, I'm really sorry that this is happening to you, I need to make this more transparent.

rm -f irfretroicor.nii irfretroicor+orig.????
3dAFNItoNIFTI -float irfretroicor+orig
fname='echo $correpi | sed "s/\+orig//"'
rm -f $fname.irfretroicor+orig.????
3dcopy irfretroicor.nii $fname.irfretroicor
rm irfretroicor+orig.* irfretroicor.nii

Thanks for trying PESTICA Isaac, and I will do whatever I can to help.  As to the IRFs, as long as they are consistent across subjects and all the heart/breathing rates are physiologically likely (i.e. the breathing rate isn't 2Hz), I would be confident that it is working.  Its been essentially a binary decision in all the data sets I've had parallel monitored pulse/respiration to compare with - either the periodicity matches very closely or it doesn't match at all, and when it fails the coupling maps don't look like they should and the IRFs are arbitrary across subjects.
Erik

Sorry, I don't know whats wrong with me, I did indeed put a huge bug in the script and every time I looked at it, I kept thinking, "no, its supposed to delete the temp file".  I even pasted in the bug when I posted here.  You are correct Isaac, this is a bug, and its one that I introduced.  I'm afraid you have to change that first remove line from:
rm -f irfretroicor.nii irfretroicor+orig.????

to:
rm -f irfretroicor.nii

I don't know how I left this in there, I'm really sorry because I know that first off, its a lot of work to get any new tool started, and secondly this one takes a lot of time to run.
Erik

Hi Erik,

Thanks for your reply; I'm glad I was able to help.  We're using PESTICA on subjects for whom we have measured pulse ox data but no respiratory data, and the PESTICA estimates of the peak locations line up extremely well with our measured cardiac data.  For advertising purposes, I've attached an image illustrating this for the very first subject I ran through PESTICA.  The next thing for me is to hack together a script that gets all subjects through stages 1 and 2, then does everybody's interactive stage 3, and finally crunches stages 4 through 6 noninteractively.  It would be really nice if PESTICA could come up with plausible filters and then ask me afterwards whether they were good enough, but I'm having trouble seeing how one would automate that decision.

Thanks for the cool tool,
Isaac Schwabacher

P.S.: Another thing I noticed while I was picking through the code was that the various MATLAB scripts use the MATLAB call if(exist('some_var')==0) to check whether an optional parameter to a function needs to be initialized with a default value.  But if I'm reading the MATLAB documentation right, that test also checks the filesystem for some_var, some_var.m, and various other files in $MATLABPATH.  I imagine that this could lead to behavior that's hard to reproduce if someone has custom libraries installed.

EDIT: While I was trying to hack run_pestica.sh to read a default value for PESTICA_SLICE_TIMING from the dataset header, I had a problem where MATLAB would ignore the commands I gave it on the command line and instead prompt me for input on stdin.  A little digging revealed, to my utter disbelief, that matlab is merely a shell script that calls the executable MATLAB, and the -r option labors to preserve its argument from shell expansion by a laughable procedure involving awk and multiple invocations of /bin/echo.  I quickly replaced -r with <<< to redirect stdin directly, and the problem was solved.  This has the additional benefit of allowing MATLAB to run in the background, which it otherwise stubbornly refuses to do.

Thanks for the advertisement Isaac!  I've uploaded a temporary bug fix as a replacement master script for those of you who have not already hacked the script with Isaac's fix to workaround my bug.  I will be releasing the next version at the end of the year, and this version will have a automated filtering routine that will be easier to simply check for success/failure.  I'm trying to figure out whats taking so long in the cardiac temporal dither fix script, and if all goes well the next version should be faster.  I'll also try to output text files that AFNI and FSL can use as slicewise regressors - the trick is that they have to be slicewise to have the same accurate timing.  There are a few other things in the todo list, but it sounds like these are at the top of user's requests - if there is anything else, please let me know.

Isaac, thank you so much for figuring out an improved way of passing args to MATLAB.  That is truly awesome and has made my weekend!  Again, thanks for trying out PESTICA and I'm glad it is now working for you!
Erik