1 /***** This code is part of the AFNI software package, which is *****
2 ***** partly in the public domain and partly covered by the GPL. *****
3 ***** See https://afni.nimh.nih.gov/afni for more information. *****/
4
5 #include "mrilib.h"
6
7 static int ARGC ; /* global copies */
8 static char **ARGV ; /* for use in get_ilist_from_args() below */
9
10 static void vstep_print(void) ; /* prototypes for later funcs */
11 static float lhs_legendre( float x, float bot, float top, float n ) ;
12
13 #define IC_POLORT 66666
14
15 #ifdef USE_OMP /* 16 Aug 2021 */
16 #include <omp.h>
17 #include "thd_fitter.c"
18 #include "thd_lasso.c"
19 #endif
20
21 /*------------------------------------------------------------*/
22 /* get an integer list from the args, starting at ARGV[*iarg] */
23
get_ilist_from_args(int * iarg)24 static int * get_ilist_from_args( int *iarg )
25 {
26 char istr[16000] ; int *ilist=NULL , iaa=*iarg ;
27 ENTRY("get_ilist_from_args") ;
28 istr[0] = '\0' ;
29 for( ; iaa < ARGC && isdigit(ARGV[iaa][0]) ; iaa++ ){
30 strcat( istr , ARGV[iaa] ) ;
31 if( istr[ strlen(istr)-1 ] != ',' ) strcat( istr , "," ) ;
32 }
33 ilist = MCW_get_intlist( 999999 , istr ) ;
34 *iarg = iaa ;
35 RETURN(ilist) ;
36 }
37
38 /*------------------------------------------------------------*/
39 /* Get an integer vector from the args */
40
get_intvec_from_args(int * iarg)41 static intvec * get_intvec_from_args( int *iarg )
42 {
43 int *ilist ;
44 intvec *ivec=NULL ;
45 int jj ;
46
47 ilist = get_ilist_from_args( iarg ) ;
48 if( ilist != NULL && ilist[0] > 0 ){
49 MAKE_intvec(ivec,ilist[0]) ;
50 for( jj=1 ; jj <= ilist[0] ; jj++ ){
51 ivec->ar[jj-1] = ilist[jj] ;
52 }
53 }
54 free(ilist) ;
55 return ivec ;
56 }
57
58 /*---------------------------------------------------------------------------*/
59
Tfitter_help(void)60 static void Tfitter_help(void)
61 {
62 printf(
63 "Usage: 3dTfitter [options]\n"
64 "\n"
65 "* At each voxel, assembles and solves a set of linear equations.\n"
66 " ++ The matrix at each voxel may be the same or may be different.\n"
67 " ++ This flexibility (for voxel-wise regressors) is one feature\n"
68 " that makes 3dTfitter different from 3dDeconvolve.\n"
69 " ++ Another distinguishing feature is that 3dTfitter allows for\n"
70 " L2, L1, and L2+L1 (LASSO) regression solvers, and allows you\n"
71 " to impose sign constraints on the solution parameters.\n"
72 "\n"
73 "* Output is a bucket dataset with the beta parameters at each voxel.\n"
74 "\n"
75 "* You can also get output of fitted time series at each voxel, and\n"
76 " the error sum of squares (e.g., for generating statistics).\n"
77 "\n"
78 "* You can also deconvolve with a known kernel function (e.g., an HRF\n"
79 " model in FMRI, or an arterial input function in DSC-MRI, et cetera),\n"
80 " in which case the output dataset is a new time series dataset,\n"
81 " containing the estimate of the source function that, when convolved\n"
82 " with your input kernel function, fits the data (in each voxel).\n"
83 "\n"
84 "* The basic idea is to compute the beta_i so that the following\n"
85 " is approximately true:\n"
86 "\n"
87 " RHS(t) = sum { beta_i * LHS_i(t) }\n"
88 " i>=1\n"
89 "\n"
90 " With the '-FALTUNG' (deconvolution) option, the model expands to be\n"
91 "\n"
92 " RHS(t) = sum { K(j)*S(t-j) } + sum { beta_i * LHS_i(t) }\n"
93 " j>=0 i>=1\n"
94 "\n"
95 " where K() is the user-supplied causal kernel function, and S() is\n"
96 " the source time series to be estimated along with the betas\n"
97 " (which can be thought of as the 'baseline' fit).\n"
98 "\n"
99 "* The model basis functions LHS_i(t) and the kernel function K(t)\n"
100 " can be .1D files (fixed for all voxels) and/or 3D+time datasets\n"
101 " (different for each voxel).\n"
102 "\n"
103 "* The fitting approximation can be done in 4 different ways, minimizing\n"
104 " the errors (differences between RHS(t) and the fitted equation) in\n"
105 " the following ways:\n"
106 " ++ L2 [-l2fit option] = least sum of squares of errors\n"
107 " ++ L1 [-l1fit option] = least sum of absolute values of errors\n"
108 " ++ L2 LASSO = least sum of squares of errors, with an added\n"
109 " [-l2lasso option] L1 penalty on the size of the solution parameters\n"
110 " ++ L2 Square Root LASSO = least square root of the sum of squared errors\n"
111 " [-l2sqrtlasso option] with an added L1 penalty on the solution parameters\n"
112 "\n"
113 "***** Which fitting method is better?\n"
114 " The answer to that question depends strongly on what you are\n"
115 " going to use the results for! And on the quality of the data.\n"
116 "\n"
117 " *************************************************\n"
118 " ***** 3dTfitter is not for the casual user! *****\n"
119 " ***** It has a lot of options which let you *****\n"
120 " ***** control the complex solution process. *****\n"
121 " *************************************************\n"
122 "\n"
123 "----------------------------------\n"
124 "SPECIFYING THE EQUATIONS AND DATA:\n"
125 "----------------------------------\n"
126 " -RHS rset = Specifies the right-hand-side 3D+time dataset.\n"
127 " ('rset' can also be a 1D file with 1 column)\n"
128 " * Exactly one '-RHS' option must be given to 3dTfitter.\n"
129 "\n"
130 " -LHS lset = Specifies a column (or columns) of the left-hand-side matrix.\n"
131 " * More than one 'lset' can follow the '-LHS' option, but each\n"
132 " input filename must NOT start with the '-' character!\n"
133 " * Or you can use multiple '-LHS' options, if you prefer.\n"
134 " * Each 'lset' can be a 3D+time dataset, or a 1D file\n"
135 " with 1 or more columns.\n"
136 " * A 3D+time dataset defines one column in the LHS matrix.\n"
137 " ++ If 'rset' is a 1D file, then you cannot input a 3D+time\n"
138 " dataset with '-LHS'.\n"
139 " ++ If 'rset' is a 3D+time dataset, then the 3D+time dataset(s)\n"
140 " input with '-LHS' must have the same voxel grid as 'rset'.\n"
141 " * A 1D file defines as many columns in the LHS matrix as\n"
142 " are in the file.\n"
143 " ++ For example, you could input the LHS matrix from the\n"
144 " .xmat.1D matrix file output by 3dDeconvolve, if you wanted\n"
145 " to repeat the same linear regression using 3dTfitter,\n"
146 " for some bizarre unfathomable twisted psychotic reason.\n"
147 " (See https://shorturl.at/boxU9 for more details.)\n"
148 " ** If you have a problem where some LHS vectors might be tiny,\n"
149 " causing stability problems, you can choose to omit them\n"
150 " by using the '-vthr' option. By default, only all-zero\n"
151 " vectors will be omitted from the regression.\n"
152 " ** Note that if the scales of the LHS vectors are grossly different\n"
153 " (e.g., 0 < vector#1 < 0.01 and 0 < vector#2 < 1000),\n"
154 " then numerical errors in the calculations might cause the\n"
155 " results to be unreliable. To avoid this problem, you can\n"
156 " scale the vectors (before running 3dTfitter) so that they\n"
157 " have similar magnitudes.\n"
158 " ** Note that if you are fitting a time series dataset that has\n"
159 " nonzero mean, then at least some of your basis vectors\n"
160 " should have nonzero mean, or you won't be able to get a\n"
161 " good fit. If necessary, use '-polort 0' to fit the mean\n"
162 " value of the dataset, so that the zero-mean LHS vectors\n"
163 " can do their work in fitting the fluctuations in the data!\n"
164 " [This means you, HJJ!]\n"
165 " *** Columns are assembled in the order given on the command line,\n"
166 " which means that LHS parameters will be output in that order!\n"
167 " *** If all LHS inputs are 1D vectors AND you are using least\n"
168 " squares fitting without constraints, then 3dDeconvolve would\n"
169 " be more efficient, since each voxel would have the same set\n"
170 " of equations -- a fact that 3dDeconvolve exploits for speed.\n"
171 " ++ But who cares about CPU time? Come on baby, light my fire!\n"
172 "\n"
173 " -polort p = Add 'p+1' Legendre polynomial columns to the LHS matrix.\n"
174 " * These columns are added to the LHS matrix AFTER all other\n"
175 " columns specified by the '-LHS' option, even if the '-polort'\n"
176 " option appears before '-LHS' on the command line.\n"
177 " ** By default, NO polynomial columns will be used.\n"
178 "\n"
179 " -vthr v = The value 'v' (between 0.0 and 0.09, inclusive) defines the\n"
180 " threshold below which LHS vectors will be omitted from\n"
181 " the regression analysis. Each vector's L1 norm (sum of\n"
182 " absolute values) is computed. Any vector whose L1 norm\n"
183 " is less than or equal to 'v' times the largest L1 norm\n"
184 " will not be used in the analysis, and will get 0 weight\n"
185 " in the output. The purpose of this option is to let you\n"
186 " have tiny inputs and have them be ignored.\n"
187 " * By default, 'v' is zero ==> only exactly zero LHS columns\n"
188 " will be ignored in this case.\n"
189 " ** Prior to 18 May 2010, the built-in (and fixed) value of\n"
190 " 'v' was 0.000333. Thus, to get the old results, you should\n"
191 " use option '-vthr 0.000333' -- this means YOU, Rasmus Birn!\n"
192 " * Note that '-vthr' column censoring is done separately for\n"
193 " each voxel's regression problem, so if '-LHS' had any\n"
194 " dataset components (i.e., voxelwise regressors), a different\n"
195 " set of omitted columns could be used betwixt different voxels.\n"
196 "\n"
197 "--------------\n"
198 "DECONVOLUTION:\n"
199 "--------------\n"
200 " -FALTUNG fset fpre pen fac\n"
201 " = Specifies a convolution (German: Faltung) model to be\n"
202 " added to the LHS matrix. Four arguments follow the option:\n"
203 "\n"
204 " -->** 'fset' is a 3D+time dataset or a 1D file that specifies\n"
205 " the known kernel of the convolution.\n"
206 " * fset's time point [0] is the 0-lag point in the kernel,\n"
207 " [1] is the 1-lag into the past point, etc.\n"
208 " ++ Call the data Z(t), the unknown signal S(t), and the\n"
209 " known kernel H(t). The equations being solved for\n"
210 " the set of all S(t) values are of the form\n"
211 " Z(t) = H(0)S(t) + H(1)S(t-1) + ... + H(L)S(t-L) + noise\n"
212 " where L is the last index in the kernel function.\n"
213 " ++++ N.B.: The TR of 'fset' (the source of H) and the TR of the\n"
214 " RHS dataset (the source of Z) MUST be the same, or\n"
215 " the deconvolution results will be revoltingly\n"
216 " meaningless drivel (or worse)!\n"
217 "\n"
218 " -->** 'fpre' is the prefix for the output time series S(t) to\n"
219 " be created -- it will have the same length as the input\n"
220 " 'rset' time series.\n"
221 " ++ If you don't want this time series (why?), set 'fpre'\n"
222 " to be the string 'NULL'.\n"
223 " ++ If you want to see the fit of the model to the data\n"
224 " (a very good idea), use the '-fitts' option, which is\n"
225 " described later.\n"
226 "\n"
227 " -->** 'pen' selects the type of penalty function to be\n"
228 " applied to constrain the deconvolved time series:\n"
229 " ++ The following penalty functions are available:\n"
230 " P0[s] = f^q * sum{ |S(t)|^q }\n"
231 " P1[s] = f^q * sum{ |S(t)-S(t-1)|^q }\n"
232 " P2[s] = f^q * sum{ |2*S(t)-S(t-1)-S(t+1)|^q }\n"
233 " P3[s] = f^q * sum{ |3*S(t)-3*S(t-1)-S(t+1)+S(t-2)|^q }\n"
234 " where S(t) is the deconvolved time series;\n"
235 " where q=1 for L1 fitting, q=2 for L2 fitting;\n"
236 " where f is the value of 'fac' (defined below).\n"
237 " P0 tries to keep S(t) itself small\n"
238 " P1 tries to keep point-to-point fluctuations\n"
239 " in S(t) small (1st derivative)\n"
240 " P2 tries to keep 3 point fluctuations\n"
241 " in S(t) small (2nd derivative)\n"
242 " P3 tries to keep 4 point fluctuations\n"
243 " in S(t) small (3nd derivative)\n"
244 " ++ Higher digits try to make the result function S(t)\n"
245 " smoother. If a smooth result makes sense, then use\n"
246 " the string '012' or '0123' for 'pen'.\n"
247 " ++ In L2 regression, these penalties are analogous to Wiener\n"
248 " (frequency space) deconvolution, with noise spectra\n"
249 " proportional to\n"
250 " P0 ==> fac^2 * 1 (constant in frequency)\n"
251 " P1 ==> fac^2 * freq^2\n"
252 " P2 ==> fac^2 * freq^4\n"
253 " P3 ==> fac^2 * freq^6\n"
254 " However, 3dTfitter does deconvolution in the time\n"
255 " domain, not the frequency domain, and you can choose\n"
256 " to use L2, L1, or LASSO (L2+L1) regression.\n"
257 " ++ The value of 'pen' is a combination of the digits\n"
258 " '0', '1', '2', and/or '3'; for example:\n"
259 " 0 = use P0 only\n"
260 " 1 = use P1 only\n"
261 " 2 = use P2 only\n"
262 " 3 = use P3 only\n"
263 " 01 = use P0+P1 (the sum of these two functions)\n"
264 " 02 = use P0+P2\n"
265 " 12 = use P1+P2\n"
266 " 012 = use P0+P1+P2 (sum of three penalty functions)\n"
267 " 0123 = use P0+P1+P2+P3 (et cetera)\n"
268 " If 'pen' does not contain any of the digits 0..3,\n"
269 " then '01' will be used.\n"
270 "\n"
271 " -->** 'fac' is the positive weight 'f' for the penalty function:\n"
272 " ++ if fac < 0, then the program chooses a penalty factor\n"
273 " for each voxel separately and then scales that by -fac.\n"
274 " ++ use fac = -1 to get this voxel-dependent factor unscaled.\n"
275 " (this is a very reasonable place to start, by the way :-)\n"
276 " ++ fac = 0 is a special case: the program chooses a range\n"
277 " of penalty factors, does the deconvolution regression\n"
278 " for each one, and then chooses the fit it likes best\n"
279 " (as a tradeoff between fit error and solution size).\n"
280 " ++ fac = 0 will be MUCH slower since it solves about 20\n"
281 " problems for each voxel and then chooses what it likes.\n"
282 " setenv AFNI_TFITTER_VERBOSE YES to get some progress\n"
283 " reports, if you want to see what it is doing.\n"
284 " ++ Instead of using fac = 0, a useful alternative is to\n"
285 " do some test runs with several negative values of fac,\n"
286 " [e.g., -1, -2, and -3] and then look at the results to\n"
287 " determine which one is most suitable for your purposes.\n"
288 " ++ It is a good idea to experiment with different fac values,\n"
289 " so you can see how the solution varies, and so you can get\n"
290 " some idea of what penalty level to use for YOUR problems.\n"
291 " ++ SOME penalty has to be applied, since otherwise the\n"
292 " set of linear equations for S(t) is under-determined\n"
293 " and/or ill-conditioned!\n"
294 "\n"
295 " ** If '-LHS' is used with '-FALTUNG', those basis vectors can\n"
296 " be thought of as a baseline to be regressed out at the\n"
297 " same time the convolution model is fitted.\n"
298 " ++ When '-LHS' supplies a baseline, it is important\n"
299 " that penalty type 'pen' include '0', so that the\n"
300 " collinearity between convolution with a constant S(t)\n"
301 " and a constant baseline can be resolved!\n"
302 " ++ Instead of using a baseline here, you could project the\n"
303 " baseline out of a dataset or 1D file using 3dDetrend,\n"
304 " before using 3dTfitter.\n"
305 "\n"
306 " *** At most one '-FALTUNG' option can be used!!!\n"
307 "\n"
308 " *** Consider the time series model\n"
309 " Z(t) = K(t)*S(t) + baseline + noise,\n"
310 " where Z(t) = data time series (in each voxel)\n"
311 " K(t) = kernel (e.g., hemodynamic response function)\n"
312 " S(t) = stimulus time series\n"
313 " baseline = constant, drift, etc.\n"
314 " and * = convolution in time\n"
315 " Then program 3dDeconvolve solves for K(t) given S(t), whereas\n"
316 " 3dTfitter -FALTUNG solves for S(t) given K(t). The difference\n"
317 " between the two cases is that K(t) is presumed to be causal and\n"
318 " have limited support, while S(t) is a full-length time series.\n"
319 "\n"
320 " *** Presumably you know this already, but deconvolution in the\n"
321 " Fourier domain -1\n"
322 " S(t) = F { F[Z] / F[K] }\n"
323 " (where F[] is the Fourier transform) is a bad idea, since\n"
324 " division by small values F[K] will grotesquely amplify the\n"
325 " noise. 3dTfitter does NOT even try to do such a silly thing.\n"
326 "\n"
327 " ****** Deconvolution is a tricky business, so be careful out there!\n"
328 " ++ e.g., Experiment with the different parameters to make\n"
329 " sure the results in your type of problems make sense.\n"
330 " -->>++ Look at the results and the fits with AFNI (or 1dplot)!\n"
331 " Do not blindly assume that the results are accurate.\n"
332 " ++ Also, do not blindly assume that a paper promoting\n"
333 " a new deconvolution method that always works is\n"
334 " actually a good thing!\n"
335 " ++ There is no guarantee that the automatic selection of\n"
336 " of the penalty factor herein will give usable results\n"
337 " for your problem!\n"
338 " ++ You should probably use a mask dataset with -FALTUNG,\n"
339 " since deconvolution can often fail on pure noise\n"
340 " time series.\n"
341 " ++ Unconstrained (no '-cons' options) least squares ('-lsqfit')\n"
342 " is normally the fastest solution method for deconvolution.\n"
343 " This, however, may only matter if you have a very long input\n"
344 " time series dataset (e.g., more than 1000 time points).\n"
345 " ++ For unconstrained least squares deconvolution, a special\n"
346 " sparse matrix algorithm is used for speed. If you wish to\n"
347 " disable this for some reason, set environment variable\n"
348 " AFNI_FITTER_RCMAT to NO before running the program.\n"
349 " ++ Nevertheless, a FALTUNG problem with more than 1000 time\n"
350 " points will probably take a LONG time to run, especially\n"
351 " if 'fac' is chosen to be 0.\n"
352 "\n"
353 "----------------\n"
354 "SOLUTION METHOD:\n"
355 "----------------\n"
356 " -lsqfit = Solve equations via least squares [the default method].\n"
357 " * This is sometimes called L2 regression by mathematicians.\n"
358 " * '-l2fit' and '-L2' are synonyms for this option.\n"
359 "\n"
360 " -l1fit = Solve equations via least sum of absolute residuals.\n"
361 " * This is sometimes called L1 regression by mathematicians.\n"
362 " * '-L1' is a synonym for this option.\n"
363 " * L1 fitting is usually slower than L2 fitting, but\n"
364 " is perhaps less sensitive to outliers in the data.\n"
365 " ++ L1 deconvolution might give nicer looking results\n"
366 " when you expect the deconvolved signal S(t) to\n"
367 " have large-ish sections where S(t) = 0.\n"
368 " [The LASSO solution methods can also have this property.]\n"
369 " * L2 fitting is statistically more efficient when the\n"
370 " noise is KNOWN to be normally (Gaussian) distributed\n"
371 " (and a bunch of other assumptions are also made).\n"
372 " ++ Where such KNOWLEDGE comes from is an interesting question.\n"
373 "\n"
374 " -l2lasso lam [i j k ...]\n"
375 " = Solve equations via least squares with a LASSO (L1) penalty\n"
376 " on the coefficients.\n"
377 " * The positive value 'lam' after the option name is the\n"
378 " weight given to the penalty.\n"
379 " ++ As a rule of thumb, you can try lam = 2 * sigma, where\n"
380 " sigma = standard deviation of noise, but that requires\n"
381 " you to have some idea what the noise level is.\n"
382 " ++ If you enter 'lam' as a negative number, then the code\n"
383 " will CRUDELY estimate sigma and then scale abs(lam) by\n"
384 " that value -- in which case, you can try lam = -2 (or so)\n"
385 " and see if that works well for you.\n"
386 " ++ Or you can use the Square Root LASSO option (next), which\n"
387 " (in theory) does not need to know sigma when setting lam.\n"
388 " ++ If you do not provide lam, or give a value of 0, then a\n"
389 " default value will be used.\n"
390 " * Optionally, you can supply a list of parameter indexes\n"
391 " (after 'lam') that should NOT be penalized in the\n"
392 " the fitting process (e.g., traditionally, the mean value\n"
393 " is not included in the L1 penalty.) Indexes start at 1,\n"
394 " as in 'consign' (below).\n"
395 " ++ If this un-penalized integer list has long stretches of\n"
396 " contiguous entries, you can specify ranges of integers,\n"
397 " as in '1:9' instead of '1 2 3 4 5 6 7 8 9'.\n"
398 " **-->>++ If you want to supply the list of indexes that GET a\n"
399 " L1 penalty, instead of the list that does NOT, you can\n"
400 " put an 'X' character first, as in\n"
401 " -LASSO 0 X 12:41\n"
402 " to indicate that variables 12..41 (inclusive) get the\n"
403 " penalty applied, and the other variables do not. This\n"
404 " inversion might be more useful to you in some cases.\n"
405 " ++ If you also want the indexes to have 1 added to them and\n"
406 " be inverted -- because they came from a 0-based program -- \n"
407 " then use 'X1', as in '-LASSO 0 X1 12:41'.\n"
408 " ++ If you want the indexes to have 1 added to them but NOT\n"
409 " to be inverted, use 'Y1', as in '-LASSO 0 Y1 13:42'.\n"
410 " ++ Note that if you supply an integer list, you MUST supply\n"
411 " a value for lam first, even if that value is 0.\n"
412 " ++ In deconvolution ('-FALTUNG'), all baseline parameters\n"
413 " (from '-LHS' and/or '-polort') are automatically non-penalized,\n"
414 " so there is usually no point to using this un-penalizing feature.\n"
415 " ++ If you are NOT doing deconvolution, then you'll need this\n"
416 " option to un-penalize any '-polort' parameters (if desired).\n"
417 " ** LASSO-ing herein should be considered experimental, and its\n"
418 " implementation is subject to change! You should definitely\n"
419 " play with different 'lam' values to see how well they work\n"
420 " for your particular types of problems. Algorithm is here:\n"
421 " ++ TT Wu and K Lange.\n"
422 " Coordinate descent algorithms for LASSO penalized regression.\n"
423 " Annals of Applied Statistics, 2: 224-244 (2008).\n"
424 " http://arxiv.org/abs/0803.3876\n"
425 " * '-LASSO' is a synonym for this option.\n"
426 "\n"
427 " -lasso_centro_block i j k ...\n"
428 " = Defines a block of coefficients that will be penalized together\n"
429 " with ABS( beta[i] - centromean( beta[i], beta[j] , ... ) )\n"
430 " where the centromean(a,b,...) is computed by sorting the\n"
431 " arguments (a,b,...) and then averaging the central 50%% values.\n"
432 " * The goal is to use LASSO to shrink these coefficients towards\n"
433 " a common value to suppress outliers, rather than the default\n"
434 " LASSO method of shrinking coefficients towards 0, where the\n"
435 " penalty on coefficient beta[i] is just ABS( beta[i] ).\n"
436 " * For example:\n"
437 " -lasso_centro_block 12:26 -lasso_centro_block 27:41\n"
438 " These options define two blocks of coefficients.\n"
439 " -->>*** The intended application of this option is to regularize\n"
440 " (reduce fluctuations) in the 'IM' regression method from\n"
441 " 3dDeconvolve, where each task instance gets a separate\n"
442 " beta fit parameter.\n"
443 " *** That is, the idea is that you run 3dTfitter to get the\n"
444 " 'IM' betas as an alternative to 3dDeconvolve or 3dREMLfit,\n"
445 " since the centromean regularization will damp down wild\n"
446 " fluctuations in the individual task betas.\n"
447 " *** In this example, the two blocks of coefficients correspond\n"
448 " to the beta values for each of two separate tasks.\n"
449 " *** The input '-LHS' matrix is available from 3dDeconvolve's\n"
450 " '-x1D' option.\n"
451 " *** Further details on 'blocks' can be found in this Google Doc\n"
452 " https://shorturl.at/boxU9\n"
453 " including shell commands on how to extract the block indexes\n"
454 " from the header of the matrix file.\n"
455 " *** A 'lam' value for the '-LASSO' option that makes sense is a value\n"
456 " between -1 and -2, but as usual, you'll have to experiment with\n"
457 " your particular data and application.\n"
458 " * If you have more than one block, do NOT let them overlap,\n"
459 " because the program doesn't check for this kind of stoopidity\n"
460 " and then peculiar/bad things will probably happen!\n"
461 " * A block defined here must have at least 5 entries.\n"
462 " In practice, I would recommend at least 12 entries for a\n"
463 " block, or the whole idea of 'shrinking to the centromean'\n"
464 " is silly.\n"
465 " * This option can be abbreviated as '-LCB', since typing\n"
466 " '-lasso_centro_block' correctly is a nontrivial challenge :-)\n"
467 " *** This option is NOT implemented for -l2sqrtlasso :-(\n"
468 " * [New option - 10 Aug 2021 - RWCox]\n"
469 "\n"
470 " -l2sqrtlasso lam [i j k ...]\n"
471 " = Similar to above option, but uses 'Square Root LASSO' instead:\n"
472 " * Approximately speaking, LASSO minimizes E = Q2+lam*L1,\n"
473 " where Q2=sum of squares of residuals and L1=sum of absolute\n"
474 " values of all fit parameters, while Square Root LASSO minimizes\n"
475 " sqrt(Q2)+lam*L1; the method and motivation is described here:\n"
476 " ++ A Belloni, V Chernozhukov, and L Wang.\n"
477 " Square-root LASSO: Pivotal recovery of sparse signals via\n"
478 " conic programming (2010). http://arxiv.org/abs/1009.5689\n"
479 " ++ A coordinate descent algorithm is also used for this optimization\n"
480 " (unlike in the paper above).\n"
481 " ** A reasonable range of 'lam' to use is from 1 to 10 (or so);\n"
482 " I suggest you start with 2 and see how well that works.\n"
483 " ++ Unlike the pure LASSO option above, you do not need to give\n"
484 " give a negative value for lam here -- there is no need for\n"
485 " scaling by sigma -- or so they say.\n"
486 " * The theoretical advantange of Square Root LASSO over\n"
487 " standard LASSO is that a good choice of 'lam' does not\n"
488 " depend on knowing the noise level in the data (that is\n"
489 " what 'Pivotal' means in the paper's title).\n"
490 " * '-SQRTLASSO' is a synonym for this option.\n"
491 "\n"
492 " --------->>**** GENERAL NOTES ABOUT LASSO and SQUARE ROOT LASSO ****<<--------\n"
493 " * LASSO methods are the only way to solve a under-determined\n"
494 " system with 3dTfitter -- one with more vectors on the RHS\n"
495 " than time points. However, a 'solution' to such a problem\n"
496 " doesn't necessarily mean anything -- be careful out there!\n"
497 " * LASSO methods will tend to push small coefficients down\n"
498 " to zero. This feature can be useful when doing deconvolution,\n"
499 " if you expect the result to be zero over large-ish intervals.\n"
500 " ++ L1 regression ('-l1fit') has a similar property, of course.\n"
501 " ++ This difficult-to-estimate bias in the LASSO-computed coefficients\n"
502 " makes it nearly impossible to provide reliable estimates of statistical\n"
503 " significance for the fit (e.g., R^2, F, ...).\n"
504 " * The actual penalty factor lambda used for a given coefficient\n"
505 " is lam scaled by the L2 norm of the corresponding regression\n"
506 " column. The purpose of this is to keep the penalties scale-free:\n"
507 " if a regression column were doubled, then the corresponding fit\n"
508 " coefficient would be cut in half; thus, to keep the same penalty\n"
509 " level, lambda should also be doubled.\n"
510 " * For '-l2lasso', a negative lam additionally means to scale\n"
511 " by the estimate of sigma, as described earlier. This feature\n"
512 " does not apply to Square Root LASSO, however (if you give a\n"
513 " negative lam to '-l2sqrtlasso', its absolute value is used).\n"
514 " -->>** There is no 'best' value of lam; if you are lucky, there is\n"
515 " is a range of lam values that give reasonable results. A good\n"
516 " procedure to follow would be to use several different values of\n"
517 " lam and see how the results vary; for example, the list\n"
518 " lam = -1, -2, -4, -7, -10 might be a good starting point.\n"
519 " * If you don't give ANY numeric value after the LASSO option\n"
520 " (i.e., the next argument on the command line is another option),\n"
521 " then the program will use '-3.1415926536' for the value of lam.\n"
522 " * A tiny value of lam (say 0.01) should give almost the same\n"
523 " results as pure L2 regression.\n"
524 " * Data with a smaller signal-to-noise ratio will probably need\n"
525 " larger values of lam -- you'll have to experiment.\n"
526 " * The number of iterations used for the LASSO solution will be\n"
527 " printed out for the first voxel solved, and for ever 10,000th\n"
528 " one following -- this is mostly for my personal edification.\n"
529 " -->>** Recall: \"3dTfitter is not for the casual user!\"\n"
530 " This statement especially applies when using LASSO, which is a\n"
531 " powerful tool -- and as such, can be dangerous if not used wisely.\n"
532 "\n"
533 "---------------------\n"
534 "SOLUTION CONSTRAINTS:\n"
535 "---------------------\n"
536 " -consign = Follow this option with a list of LHS parameter indexes\n"
537 " to indicate that the sign of some output LHS parameters\n"
538 " should be constrained in the solution; for example:\n"
539 " -consign +1 -3\n"
540 " which indicates that LHS parameter #1 (from the first -LHS)\n"
541 " must be non-negative, and that parameter #3 must be\n"
542 " non-positive. Parameter #2 is unconstrained (e.g., the\n"
543 " output can be positive or negative).\n"
544 " * Parameter counting starts with 1, and corresponds to\n"
545 " the order in which the LHS columns are specified.\n"
546 " * Unlike '-LHS or '-label', only one '-consign' option\n"
547 " can be used.\n"
548 " * Do NOT give the same index more than once after\n"
549 " '-consign' -- you can't specify that an coefficient\n"
550 " is both non-negative and non-positive, for example!\n"
551 " *** Constraints can be used with any of the 4 fitting methods.\n"
552 " *** '-consign' constraints only apply to the '-LHS'\n"
553 " fit parameters. To constrain the '-FALTUNG' output,\n"
554 " use the option below.\n"
555 " * If '-consign' is not used, the signs of the fitted\n"
556 " LHS parameters are not constrained.\n"
557 "\n"
558 " -consFAL c= Constrain the deconvolution time series from '-FALTUNG'\n"
559 " to be positive if 'c' is '+' or to be negative if\n"
560 " 'c' is '-'.\n"
561 " * There is no way at present to constrain the deconvolved\n"
562 " time series S(t) to be positive in some regions and\n"
563 " negative in others.\n"
564 " * If '-consFAL' is not used, the sign of the deconvolved\n"
565 " time series is not constrained.\n"
566 "\n"
567 "---------------\n"
568 "OUTPUT OPTIONS:\n"
569 "---------------\n"
570 " -prefix p = Prefix for the output dataset (LHS parameters) filename.\n"
571 " * Output datasets from 3dTfitter are always in float format.\n"
572 " * If you don't give this option, 'Tfitter' is the prefix.\n"
573 " * If you don't want this dataset, use 'NULL' as the prefix.\n"
574 " * If you are doing deconvolution and do not also give any\n"
575 " '-LHS' options, then this file will not be output, since\n"
576 " it comprises the fit parameters for the '-LHS' vectors.\n"
577 " -->>** If the input '-RHS' file is a .1D file, normally the\n"
578 " output files are written in the AFNI .3D ASCII format,\n"
579 " where each row contains the time series data for one\n"
580 " voxel. If you want to have these files written in the\n"
581 " .1D format, with time represented down the column\n"
582 " direction, be sure to put '.1D' on the end of the prefix,\n"
583 " as in '-prefix Elvis.1D'. If you use '-' or 'stdout' as\n"
584 " the prefix, the resulting 1D file will be written to the\n"
585 " terminal. (See the fun fun fun examples, below.)\n"
586 "\n"
587 " -label lb = Specifies sub-brick labels in the output LHS parameter dataset.\n"
588 " * More than one 'lb' can follow the '-label' option;\n"
589 " however, each label must NOT start with the '-' character!\n"
590 " * Labels are applied in the order given.\n"
591 " * Normally, you would provide exactly as many labels as\n"
592 " LHS columns. If not, the program invents some labels.\n"
593 "\n"
594 " -fitts ff = Prefix filename for the output fitted time series dataset.\n"
595 " * Which is always in float format.\n"
596 " * Which will not be written if this option isn't given!\n"
597 " *** If you want the residuals, subtract this time series\n"
598 " from the '-RHS' input using 3dcalc (or 1deval).\n"
599 "\n"
600 " -errsum e = Prefix filename for the error sums dataset, which\n"
601 " is calculated from the difference between the input\n"
602 " time series and the fitted time series (in each voxel):\n"
603 " * Sub-brick #0 is the sum of squares of differences (L2 sum)\n"
604 " * Sub-brick #1 is the sum of absolute differences (L1 sum)\n"
605 " * The L2 sum value, in particular, can be used to produce\n"
606 " a statistic to measure the significance of a fit model;\n"
607 " cf. the 'Correlation Coefficient Example' far below.\n"
608 "\n"
609 "--------------\n"
610 "OTHER OPTIONS:\n"
611 "--------------\n"
612 " -mask ms = Read in dataset 'ms' as a mask; only voxels with nonzero\n"
613 " values in the mask will be processed. Voxels falling\n"
614 " outside the mask will be set to all zeros in the output.\n"
615 " * Voxels whose time series are all zeros will not be\n"
616 " processed, even if they are inside the mask!\n"
617 "\n"
618 " -quiet = Don't print the fun fun fun progress report messages.\n"
619 " * Why would you want to hide these delightful missives?\n"
620 "\n"
621 "----------------------\n"
622 "ENVIRONMENT VARIABLES:\n"
623 "----------------------\n"
624 " AFNI_TFITTER_VERBOSE = YES means to print out information during\n"
625 " the fitting calculations.\n"
626 " ++ Automatically turned on for 1 voxel -RHS inputs.\n"
627 " AFNI_TFITTER_P1SCALE = number > 0 will scale the P1 penalty by\n"
628 " this value (e.g., to count it more)\n"
629 " AFNI_TFITTER_P2SCALE = number > 0 will scale the P2 penalty by\n"
630 " this value\n"
631 " AFNI_TFITTER_P3SCALE = number > 0 will scale the P3 penalty by\n"
632 " this value\n"
633 " You could set these values on the command line using the AFNI standard\n"
634 " '-Dvariablename=value' command line option.\n"
635 "\n"
636 "------------\n"
637 "NON-Options:\n"
638 "------------\n"
639 "* There is no option to produce statistical estimates of the\n"
640 " significance of the parameter estimates.\n"
641 " ++ 3dTcorrelate might be useful, to compute the correlation\n"
642 " between the '-fitts' time series and the '-RHS' input data.\n"
643 " ++ You can use the '-errsum' option to get around this limitation,\n"
644 " with enough cleverness.\n"
645 "* There are no options for censoring or baseline generation (except '-polort').\n"
646 " ++ You could generate some baseline 1D files using 1deval, perhaps.\n"
647 "* There is no option to constrain the range of the output parameters,\n"
648 " except the semi-infinite ranges provided by '-consign' and/or '-consFAL'.\n"
649 "* This program is NOW parallelized via OpenMP :-) [17 Aug 2021 - RWCox]\n"
650 "\n"
651 "------------------\n"
652 "Contrived Example:\n"
653 "------------------\n"
654 "The datasets 'atm' and 'btm' are assumed to have 99 time points each.\n"
655 "We use 3dcalc to create a synthetic combination of these plus a constant\n"
656 "plus Gaussian noise, then use 3dTfitter to fit the weights of these\n"
657 "3 functions to each voxel, using 4 different methods. Note the use of\n"
658 "the input 1D time series '1D: 99@1' to provide the constant term.\n"
659 "\n"
660 " 3dcalc -a atm+orig -b btm+orig -expr '-2*a+b+gran(100,20)' -prefix 21 -float\n"
661 " 3dTfitter -RHS 21+orig -LHS atm+orig btm+orig '1D: 99@1' -prefix F2u -l2fit\n"
662 " 3dTfitter -RHS 21+orig -LHS atm+orig btm+orig '1D: 99@1' -prefix F1u -l1fit\n"
663 " 3dTfitter -RHS 21+orig -LHS atm+orig btm+orig '1D: 99@1' -prefix F1c -l1fit \\\n"
664 " -consign -1 +3\n"
665 " 3dTfitter -RHS 21+orig -LHS atm+orig btm+orig '1D: 99@1' -prefix F2c -l2fit \\\n"
666 " -consign -1 +3\n"
667 "\n"
668 "In the absence of noise and error, the output datasets should be\n"
669 " #0 sub-brick = -2.0 in all voxels\n"
670 " #1 sub-brick = +1.0 in all voxels\n"
671 " #2 sub-brick = +100.0 in all voxels\n"
672 "\n"
673 "----------------------\n"
674 "Yet More Contrivances:\n"
675 "----------------------\n"
676 "You can input a 1D file for the RHS dataset, as in the example below,\n"
677 "to fit a single time series to a weighted sum of other time series:\n"
678 "\n"
679 " 1deval -num 30 -expr 'cos(t)' > Fcos.1D\n"
680 " 1deval -num 30 -expr 'sin(t)' > Fsin.1D\n"
681 " 1deval -num 30 -expr 'cos(t)*exp(-t/20)' > Fexp.1D\n"
682 " 3dTfitter -quiet -RHS Fexp.1D -LHS Fcos.1D Fsin.1D -prefix -\n"
683 "\n"
684 "* Note the use of the '-' as a prefix to write the results\n"
685 " (just 2 numbers) to stdout, and the use of '-quiet' to hide\n"
686 " the divertingly funny and informative progress messages.\n"
687 "* For the Jedi AFNI Masters out there, the above example can be carried\n"
688 " out on using single complicated command line:\n"
689 "\n"
690 " 3dTfitter -quiet -RHS `1deval -1D: -num 30 -expr 'cos(t)*exp(-t/20)'` \\\n"
691 " -LHS `1deval -1D: -num 30 -expr 'cos(t)'` \\\n"
692 " `1deval -1D: -num 30 -expr 'sin(t)'` \\\n"
693 " -prefix - \n"
694 "\n"
695 " resulting in the single output line below:\n"
696 "\n"
697 " 0.535479 0.000236338\n"
698 "\n"
699 " which are respectively the fit coefficients of 'cos(t)' and 'sin(t)'.\n"
700 "\n"
701 "--------------------------------\n"
702 "Contrived Deconvolution Example:\n"
703 "--------------------------------\n"
704 "(1) Create a 101 point 1D file that is a block of 'activation'\n"
705 " between points 40..50, convolved with a triangle wave kernel\n"
706 " (the '-iresp' input below):\n"
707 " 3dConvolve -input1D -polort -1 -num_stimts 1 \\\n"
708 " -stim_file 1 '1D: 40@0 10@1 950@0' \\\n"
709 " -stim_minlag 1 0 -stim_maxlag 1 5 \\\n"
710 " -iresp 1 '1D: 0 1 2 3 2 1' -nlast 100 \\\n"
711 " | grep -v Result | grep -v '^$' > F101.1D\n"
712 "\n"
713 "(2) Create a 3D+time dataset with this time series in each\n"
714 " voxel, plus noise that increases with voxel 'i' index:\n"
715 " 3dUndump -prefix Fjunk -dimen 100 100 1\n"
716 " 3dcalc -a Fjunk+orig -b F101.1D \\\n"
717 " -expr 'b+gran(0,0.04*(i+1))' \\\n"
718 " -float -prefix F101d\n"
719 " /bin/rm -f Fjunk+orig.*\n"
720 "\n"
721 "(3) Deconvolve, then look what you get by running AFNI:\n"
722 " 3dTfitter -RHS F101d+orig -l1fit \\\n"
723 " -FALTUNG '1D: 0 1 2 3 2 1' F101d_fal1 012 0.0\n"
724 " 3dTfitter -RHS F101d+orig -l2fit \\\n"
725 " -FALTUNG '1D: 0 1 2 3 2 1' F101d_fal2 012 0.0\n"
726 "\n"
727 "(4) View F101d_fal1+orig, F101d_fal2+orig, and F101d+orig in AFNI,\n"
728 " (in Axial image and graph viewers) and see how the fit quality\n"
729 " varies with the noise level and the regression type -- L1 or\n"
730 " L2 regression. Note that the default 'fac' level of 0.0 was\n"
731 " selected in the commands above, which means the program selects\n"
732 " the penalty factor for each voxel, based on the size of the\n"
733 " data time series fluctuations and the quality of the fit.\n"
734 "\n"
735 "(5) Add logistic noise (long tails) to the noise-free 1D time series, then\n"
736 " deconvolve and plot the results directly to the screen, using L1 and L2\n"
737 " and the two LASSO fitting methods:\n"
738 " 1deval -a F101.1D -expr 'a+lran(.5)' > F101n.1D\n"
739 " 3dTfitter -RHS F101n.1D -l1fit \\\n"
740 " -FALTUNG '1D: 0 1 2 3 2 1' stdout 01 -2 | 1dplot -stdin -THICK &\n"
741 " 3dTfitter -RHS F101n.1D -l2fit \\\n"
742 " -FALTUNG '1D: 0 1 2 3 2 1' stdout 01 -2 | 1dplot -stdin -THICK &\n"
743 " 3dTfitter -RHS F101n.1D -l2sqrtlasso 2 \\\n"
744 " -FALTUNG '1D: 0 1 2 3 2 1' stdout 01 -2 | 1dplot -stdin -THICK &\n"
745 " 3dTfitter -RHS F101n.1D -l2lasso -2 \\\n"
746 " -FALTUNG '1D: 0 1 2 3 2 1' stdout 01 -2 | 1dplot -stdin -THICK &\n"
747 " For even more fun, add the '-consfal +' option to the above commands,\n"
748 " to force the deconvolution results to be positive.\n"
749 "\n"
750 " ***N.B.: You can only use 'stdout' as an output filename when\n"
751 " the output will be written as a 1D file (as above)!\n"
752 "\n"
753 "--------------------------------\n"
754 "Correlation Coefficient Example:\n"
755 "--------------------------------\n"
756 "Suppose your initials are HJJ and you want to compute the partial\n"
757 "correlation coefficient of time series Seed.1D with every voxel in\n"
758 "a dataset Rest+orig once a spatially dependent 'artifact' time series\n"
759 "Art+orig has been projected out. You can do this with TWO 3dTfitter\n"
760 "runs, plus 3dcalc:\n"
761 "\n"
762 "(1) Run 3dTfitter with ONLY the artifact time series and get the\n"
763 " error sum dataset\n"
764 " 3dTfitter -RHS Rest+orig -LHS Art+orig -polort 2 -errsum Ebase\n"
765 "\n"
766 "(2) Run 3dTfitter again with the artifact PLUS the seed time series\n"
767 " and get the error sum dataset and also the beta coefficents\n"
768 " 3dTfitter -RHS Rest+orig -LHS Seed.1D Art+orig -polort 2 \\\n"
769 " -errsum Eseed -prefix Bseed\n"
770 "\n"
771 "(3) Compute the correlation coefficient from the amount of variance\n"
772 " reduction between cases 1 and 2, times the sign of the beta\n"
773 " 3dcalc -a Eseed+orig'[0]' -b Ebase+orig'[0]' -c Bseed+orig'[0]' \\\n"
774 " -prefix CorrSeed -expr '(2*step(c)-1)*sqrt(1-a/b)'\n"
775 " 3drefit -fbuc -sublabel 0 'SeedCorrelation' CorrSeed+orig\n"
776 "\n"
777 "More cleverness could be used to compute t- or F-statistics in a\n"
778 "similar fashion, using the error sum of squares between 2 different fits.\n"
779 "(Of course, these are assuming you use the default '-lsqfit' method.)\n"
780 "\n"
781 "--------------------------------\n"
782 "PPI (psycho-physiological interaction) Example:\n"
783 "--------------------------------\n"
784 "Suppose you are running a PPI analysis and want to deconvolve a GAM\n"
785 "signal from the seed time series, hoping (very optimistically) to\n"
786 "convert from the BOLD time series (typical FMRI signal) to a\n"
787 "neurological time series (an impulse signal, say).\n"
788 "\n"
789 "If the BOLD signal at the seed is seed_BOLD.1D and the GAM signal is\n"
790 "GAM.1D, then consider this example for the deconvolution, in order to\n"
791 "create the neuro signal, seed_neuro.1D:\n"
792 "\n"
793 " 3dTfitter -RHS seed_BOLD.1D \\\n"
794 " -FALTUNG GAM.1D seed_neuro.1D 012 -2 \\\n"
795 " -l2lasso -6\n"
796 "\n"
797 "*************************************************************************\n"
798 "** RWCox - Feb 2008, et seq. **\n"
799 "** Created for the glorious purposes of John A Butman, MD, PhD, Poobah **\n"
800 "** But might be useful for some other well-meaning souls out there **\n"
801 "*************************************************************************\n"
802 ) ;
803
804 PRINT_AFNI_OMP_USAGE("3dTfitter",NULL) ;
805 PRINT_COMPILE_DATE ; exit(0) ;
806 }
807
808 /*---------------------------------------------------------------------------*/
809
main(int argc,char * argv[])810 int main( int argc , char *argv[] )
811 {
812 int iarg , ii,jj,kk , nx,ny,nz,nvox , vstep=0 ;
813 THD_3dim_dataset *rhset=NULL ; char *rhsnam="?" ; int rhs1D=0 ;
814 THD_3dim_dataset *lset ; MRI_IMAGE *lim ; int nlset=0 , nlhs=0 ;
815 THD_3dim_dataset *fset=NULL ;
816 RwcPointer_array *dsar ;
817 int ntime , nvar=0 , polort=-1,npol=0 ;
818 char *prefix="Tfitter" ;
819 int meth=2 , nbad=0,ngood=0,nskip=0 ;
820 intvec *convec=NULL , *kvec=NULL ;
821 byte *mask=NULL ; int mnx=0,mny=0,mnz=0 ;
822 float **rvec=NULL , *cvec=NULL ;
823 char **lab=NULL ; int nlab=0 ;
824 int verb=1 ;
825
826 THD_3dim_dataset *fal_set=NULL ; MRI_IMAGE *fal_im=NULL ;
827 char *fal_pre=NULL ; int fal_pencod=3, fal_klen=0 , fal_dcon=0 ;
828 float *fal_kern=NULL , fal_penfac=0.0f ;
829 THD_3dim_dataset *defal_set=NULL ;
830 int nvoff=0 ;
831
832 char *fitts_prefix=NULL; THD_3dim_dataset *fitts_set=NULL;
833 char *ersum_prefix=NULL; THD_3dim_dataset *ersum_set=NULL; /* 23 Jul 2009 */
834 int do_fitts=0 ;
835
836 float vthresh=0.0f ; /* 18 May 2010 */
837
838 intvec *lasso_ivec = NULL ; /* 11 Mar 2011 */
839 int lasso_invert_ivec=0 ; /* 04 Aug 2021 */
840 int lasso_zerobase_ivec=0 ; /* 04 Aug 2021 */
841 float lasso_flam = 8.0f ;
842
843 /*------------------ help the pitifully ignorant user? ------------------*/
844
845 ARGC = argc ; /* 06 Aug 2016 */
846 ARGV = argv ;
847
848 /* help function does not return */
849
850 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) Tfitter_help() ;
851
852 /*------- official startup for the history books -------*/
853
854 mainENTRY("3dTfitter"); machdep();
855 PRINT_VERSION("3dTfitter"); AUTHOR("RWCox") ;
856 AFNI_logger("3dTfitter",argc,argv);
857 AFNI_SETUP_OMP(0) ;
858 (void)COX_clock_time() ;
859
860 /*------------ read command line args ------------*/
861
862 iarg = 1 ; INIT_XTARR(dsar) ;
863 while( iarg < argc ){
864
865 /*-----*/
866
867 if( strcasecmp(argv[iarg],"-vthr" ) == 0 ||
868 strcasecmp(argv[iarg],"-vthresh") == 0 ){ /* 18 May 2010 */
869 if( iarg+1 >= argc )
870 ERROR_exit("Need an argument after '%s'",argv[iarg]);
871 vthresh = (float)strtod(argv[++iarg],NULL) ;
872 if( vthresh < 0.0f ){
873 WARNING_message("-vthr value < 0.00 ==> will use 0.00"); vthresh = 0.00f;
874 } else if( vthresh > 0.09f ){
875 WARNING_message("-vthr value > 0.09 ==> will use 0.09"); vthresh = 0.09f;
876 }
877 iarg++ ; continue ;
878 }
879
880 /*-----*/
881
882 if( strcasecmp(argv[iarg],"-polort") == 0 ){ /* 20 Mar 2008 */
883 char *qpt ;
884 if( polort >= 0 )
885 WARNING_message("you have more than 1 -polort option!") ;
886 if( iarg+1 >= argc )
887 ERROR_exit("Need an argument after '%s'",argv[iarg]);
888 polort = (int)strtod(argv[++iarg],&qpt) ;
889 if( *qpt != '\0' )
890 WARNING_message("Illegal non-numeric value after -polort") ;
891 if( polort < 0 )
892 WARNING_message("-polort value is negative ==> ignoring") ;
893 iarg++ ; continue ;
894 }
895
896 /*-----*/
897
898 if( strcasecmp(argv[iarg],"-faltung") == 0 ){
899 int p0,p1,p2,p3 ;
900 if( fal_set != NULL || fal_im != NULL )
901 ERROR_exit("Can't have two -FALTUNG arguments") ;
902 if( iarg+4 >= argc )
903 ERROR_exit("Need 4 arguments after '%s'",argv[iarg]);
904 iarg++ ;
905 if( strstr(argv[iarg],"1D") != NULL ){
906 fal_im = mri_read_1D(argv[iarg]) ;
907 if( fal_im == NULL )
908 ERROR_exit("Can't read -FALTUNG time series from '%s'",argv[iarg]) ;
909 fal_klen = fal_im->nx ;
910 if( fal_klen < 2 )
911 ERROR_exit("-FALTUNG 1D file '%s' has only 1 time point!",argv[iarg]);
912 if( fal_im->ny > 1 )
913 WARNING_message("Only using first column of -FALTUNG 1D file") ;
914 fal_kern = MRI_FLOAT_PTR(fal_im) ;
915 for( jj=0 ; jj < fal_klen && fal_kern[jj] == 0.0f ; jj++ ) ; /*nada*/
916 for( jj=fal_klen-1 ; jj >= 0 && fal_kern[jj] == 0.0f ; jj-- ) ; /*nada*/
917 if( jj < 0 )
918 ERROR_exit("-FALTUNG 1D file '%s' is all zeros!",argv[iarg]) ;
919 #if 0
920 if( jj < fal_klen-1 ){ /* truncate final zeros */
921 WARNING_message("-FALTUNG 1D file '%s' has %d zeros at its end",
922 argv[iag],fal_klen-1-jj) ;
923 fal_klen = jj+1 ;
924 }
925 #endif
926 } else {
927 fal_set = THD_open_dataset(argv[iarg]) ;
928 CHECK_OPEN_ERROR(fal_set,argv[iarg]) ;
929 fal_klen = DSET_NVALS(fal_set) ;
930 if( fal_klen < 2 )
931 ERROR_exit("-FALTUNG dataset '%s' has only 1 sub-brick!",argv[iarg]);
932 }
933 fal_pre = strdup(argv[++iarg]) ;
934 if( !THD_filename_ok(fal_pre) )
935 ERROR_exit("Illegal filename prefix '%s'",argv[iarg]) ;
936 iarg++ ;
937 p0 = (strchr(argv[iarg],'0') != NULL) ;
938 p1 = (strchr(argv[iarg],'1') != NULL) ;
939 p2 = (strchr(argv[iarg],'2') != NULL) ;
940 p3 = (strchr(argv[iarg],'3') != NULL) ;
941 if( p0==0 && p1==0 && p2==0 && p3==0 ){
942 WARNING_message(
943 "-FALTUNG 'pen' value '%s' illegal: defaulting to '01'",argv[iarg]);
944 p0 = p1 = 1 ;
945 }
946 fal_pencod = p0 + 2*p1 + 4*p2 + 8*p3 ; /* encode in bits */
947 fal_penfac = (float)strtod(argv[++iarg],NULL) ;
948 if( fal_penfac == 0.0f ) fal_penfac = -666.0f ; /* autopen */
949 iarg++ ; continue ;
950 }
951
952 /*-----*/
953
954 if( strncasecmp(argv[iarg],"-consFAL",7) == 0 ){
955 if( argv[iarg][8] == '+' ) fal_dcon = 1 ;
956 else if( argv[iarg][8] == '-' ) fal_dcon = -1 ;
957 else {
958 if( ++iarg >= argc )
959 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
960 fal_dcon = (argv[iarg][0] == '+') ? 1
961 :(argv[iarg][0] == '-') ? -1 : 0 ;
962 if( fal_dcon == 0 )
963 WARNING_message("value after '-consFAL' is not '+' or '-' -- ignoring!");
964 }
965 iarg++ ; continue ;
966 }
967
968 /*-----*/
969
970 if( strcasecmp(argv[iarg],"-mask") == 0 ){
971 THD_3dim_dataset *mset ;
972 if( mask != NULL )
973 ERROR_exit("Can't have two -mask arguments!") ;
974 if( ++iarg >= argc )
975 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
976 mset = THD_open_dataset(argv[iarg]) ;
977 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
978 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
979 mnx = DSET_NX(mset); mny = DSET_NY(mset); mnz = DSET_NZ(mset);
980 mask = THD_makemask( mset, 0, 1.0f,0.0f ); DSET_delete(mset);
981 if( mask == NULL ) ERROR_exit("Can't make mask") ;
982 ii = THD_countmask( mnx*mny*mnz , mask ) ;
983 INFO_message("%d voxels in the mask",ii) ;
984 if( ii < 1 ) ERROR_exit("mask is empty?!") ;
985 iarg++ ; continue ;
986 }
987
988 /*-----*/
989
990 if( strcasecmp(argv[iarg],"-rhs") == 0 ){
991 if( ++iarg >= argc )
992 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
993 if( rhset != NULL )
994 ERROR_exit("Can't have two '%s' options",argv[iarg-1]);
995 rhsnam = malloc(sizeof(char)*(strlen(argv[iarg])+4)) ;
996 strcpy(rhsnam,argv[iarg]) ;
997 if( STRING_HAS_SUFFIX_CASE(rhsnam,"1D") || /* transpose 1D files */
998 strncmp(rhsnam,"1D:",3) == 0 ) strcat(rhsnam,"'");
999 rhset = THD_open_dataset( rhsnam ) ;
1000 if( rhset == NULL )
1001 ERROR_exit("Can't open dataset '%s'",argv[iarg]) ;
1002 rhs1D = DSET_IS_1D(rhset) ; /* 05 Mar 2008 */
1003 iarg++ ; continue ;
1004 }
1005
1006 /*-----*/
1007
1008 if( strcasecmp(argv[iarg],"-lhs") == 0 ){
1009 if( ++iarg >= argc )
1010 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1011 if( argv[iarg][0] == '-' )
1012 ERROR_exit("Illegal argument after '%s'",argv[iarg-1]) ;
1013 for( ; iarg < argc && argv[iarg][0] != '-' ; iarg++ ){
1014 if( strstr(argv[iarg],"1D") != NULL ){
1015 lim = mri_read_1D(argv[iarg]) ;
1016 if( lim == NULL )
1017 ERROR_exit("Can't read 1D file '%s'",argv[iarg]) ;
1018 if( lim->ny > 1 )
1019 INFO_message("1D file '%s' has %d columns",argv[iarg],lim->ny);
1020 ADDTO_XTARR(dsar,lim) ;
1021 ii= XTARR_NUM(dsar)-1 ; XTARR_IC(dsar,ii) = IC_FLIM ;
1022 nvar += lim->ny ;
1023 } else {
1024 lset = THD_open_dataset(argv[iarg]) ;
1025 if( lset == NULL )
1026 ERROR_exit("Can't read dataset '%s'",argv[iarg]) ;
1027 ADDTO_XTARR(dsar,lset) ;
1028 ii= XTARR_NUM(dsar)-1 ; XTARR_IC(dsar,ii) = IC_DSET ;
1029 nvar++ ; nlset++ ;
1030 }
1031 }
1032 continue ;
1033 }
1034
1035 /*-----*/
1036
1037 if( strncasecmp(argv[iarg],"-label",4) == 0 ){
1038 if( ++iarg >= argc )
1039 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1040 if( argv[iarg][0] == '-' )
1041 ERROR_exit("Illegal argument after '%s'",argv[iarg-1]) ;
1042 for( ; iarg < argc && argv[iarg][0] != '-' ; iarg++ ){
1043 lab = (char **)realloc((void *)lab,sizeof(char *)*(nlab+1)) ;
1044 lab[nlab++] = strdup(argv[iarg]) ;
1045 }
1046 continue ;
1047 }
1048
1049 /*-----*/
1050
1051 if( strncasecmp(argv[iarg],"-lsqfit",4) == 0 ||
1052 strncasecmp(argv[iarg],"-l2fit",4) == 0 ||
1053 strcmp (argv[iarg],"-L2") == 0 ){
1054 meth = 2 ; iarg++ ; continue ;
1055 }
1056
1057 /*-----*/
1058
1059 if( strcasecmp(argv[iarg],"-l1fit") == 0 || strcmp(argv[iarg],"-L1") == 0 ){
1060 meth = 1 ; iarg++ ; continue ;
1061 }
1062
1063 /*-----*/
1064
1065 if( strcasecmp(argv[iarg],"-lasso_centro_block") == 0 ||
1066 strcasecmp(argv[iarg],"-LCB") == 0 ){ /* 06 Aug 2021 */
1067 intvec *mblok ;
1068 iarg++ ;
1069 mblok = get_intvec_from_args( &iarg ) ;
1070 if( mblok == NULL ){
1071 WARNING_message("bad -lasso_centro_block (LCB): can't get integer list :(") ;
1072 } else if( mblok->nar < 5 ){
1073 ERROR_message("bad -lasso_centro_block (LCB): not enough entries :(") ;
1074 } else {
1075 THD_lasso_add_centro_block( mblok ) ;
1076 }
1077 continue ; /* iarg was updated in the function above */
1078 }
1079
1080 /*-----*/
1081
1082 #undef ISAL
1083 #define ISAL(s) (isalpha(s[0]) || isalpha(s[1]))
1084 #undef DFAL
1085 #define DFAL 3.1415926536f
1086
1087 if( strcasecmp(argv[iarg],"-l2lasso") == 0 ||
1088 strcasecmp(argv[iarg],"-LASSO" ) == 0 ){
1089 meth = -2 ; ii = 0 ; iarg++ ;
1090 if( iarg >= argc || ISAL(argv[iarg]) ){
1091 lasso_flam = -DFAL ; ii = 1 ;
1092 } else {
1093 lasso_flam = (float)strtod(argv[iarg],NULL) ;
1094 if( lasso_flam == 0.0f ) lasso_flam = -DFAL ;
1095 }
1096 if( lasso_flam < 0.0f ){ THD_lasso_dosigest(1); lasso_flam = -lasso_flam; }
1097 THD_lasso_fixlam(lasso_flam) ; /* cf. thd_lasso.c */
1098
1099 if( iarg < argc-2 && argv[iarg+1][0] == 'X' ){ /* Invert? [04 Aug 2021] */
1100 lasso_invert_ivec = 1 ;
1101 lasso_zerobase_ivec = (argv[iarg+1][1] == '1') ;
1102 iarg++ ;
1103 } else if( iarg < argc-2 && isalpha(argv[iarg+1][0]) ){
1104 lasso_invert_ivec = 0 ;
1105 lasso_zerobase_ivec = (argv[iarg+1][1] == '1') ;
1106 iarg++ ;
1107 }
1108
1109 if( ii == 0 ){ /* flag to scan for index selection */
1110 iarg++ ;
1111 if( iarg < argc && isdigit(argv[iarg][0]) ){ /* get 'free' indexes */
1112 lasso_ivec = get_intvec_from_args( &iarg ) ;
1113 if( lasso_ivec == NULL ){
1114 WARNING_message("illegal LASSO index list :(") ;
1115 }
1116 }
1117 }
1118 continue ; /* iarg was already incremented above */
1119 }
1120
1121 /*-----*/
1122
1123 if( strcasecmp(argv[iarg],"-l2sqrtlasso") == 0 ||
1124 strcasecmp(argv[iarg],"-SQRTLASSO" ) == 0 ||
1125 strcasecmp(argv[iarg],"-LASSOSQRT" ) == 0 ){ /* EXPERIMENTAL */
1126 meth = -1 ; ii = 0 ; iarg++ ;
1127 if( iarg >= argc || ISAL(argv[iarg]) ){
1128 lasso_flam = DFAL ; ii = 1 ;
1129 } else {
1130 lasso_flam = (float)strtod(argv[iarg],NULL) ;
1131 if( lasso_flam == 0.0f ) lasso_flam = DFAL ;
1132 else if( lasso_flam < 0.0f ) lasso_flam = -lasso_flam ;
1133 }
1134 THD_lasso_fixlam(lasso_flam) ; /* cf. thd_lasso.c */
1135 if( ii == 0 ){
1136 iarg++ ;
1137 if( iarg < argc && isdigit(argv[iarg][0]) ){ /* get 'free' indexes */
1138 lasso_ivec = get_intvec_from_args( &iarg ) ;
1139 if( lasso_ivec == NULL ){
1140 WARNING_message("illegal LASSO index list :(") ;
1141 }
1142 }
1143 }
1144 continue ;
1145 }
1146
1147 /*-----*/
1148
1149 if( strncasecmp(argv[iarg],"-prefix",5) == 0 ){
1150 if( ++iarg >= argc )
1151 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1152 prefix = argv[iarg] ;
1153 if( !THD_filename_ok(prefix) )
1154 ERROR_exit("Illegal string after -prefix: '%s'",prefix) ;
1155 iarg++ ; continue ;
1156 }
1157
1158 /*-----*/
1159
1160 if( strncasecmp(argv[iarg],"-fitts",5) == 0 ){
1161 if( ++iarg >= argc )
1162 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1163 fitts_prefix = argv[iarg] ;
1164 if( !THD_filename_ok(fitts_prefix) )
1165 ERROR_exit("Illegal string after -fitts: '%s'",fitts_prefix) ;
1166 do_fitts = 1 ; iarg++ ; continue ;
1167 }
1168
1169 /*-----*/
1170
1171 if( strncasecmp(argv[iarg],"-errsum",5) == 0 ){ /* 23 Jul 2009 */
1172 if( ++iarg >= argc )
1173 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1174 ersum_prefix = argv[iarg] ;
1175 if( !THD_filename_ok(ersum_prefix) )
1176 ERROR_exit("Illegal string after -ersum: '%s'",ersum_prefix) ;
1177 do_fitts = 1 ; iarg++ ; continue ;
1178 }
1179
1180 /*-----*/
1181
1182 if( strncasecmp(argv[iarg],"-consign",7) == 0 ){
1183 char *cpt ; int nvec ;
1184 if( ++iarg >= argc )
1185 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
1186 if( convec != NULL )
1187 ERROR_exit("Can't have 2 -consign options!") ;
1188 MAKE_intvec(convec,1) ;
1189 for( nvec=0 ; iarg < argc ; iarg++ ){
1190 ii = (int)strtod(argv[iarg],&cpt) ;
1191 if( ii == 0 || *cpt != '\0' ) break ; /* bad */
1192 for( jj=0 ; jj < nvec ; jj++ ){
1193 if( abs(convec->ar[jj]) == abs(ii) )
1194 ERROR_exit("Duplicate indexes in -consign!") ;
1195 }
1196 RESIZE_intvec(convec,nvec+1) ;
1197 convec->ar[nvec++] = ii ;
1198 }
1199 if( nvec < 1 )
1200 ERROR_exit("No legal values after -consign?!") ;
1201 continue ;
1202 }
1203
1204 /*-----*/
1205
1206 if( strncasecmp(argv[iarg],"-quiet",2) == 0 ){
1207 verb = 0 ; iarg++ ; continue ;
1208 }
1209
1210 /*-----*/
1211
1212 ERROR_exit("Unknown argument on command line: '%s'",argv[iarg]) ;
1213 }
1214
1215 /*------------ check options for completeness and consistency ----------*/
1216
1217 if( rhset == NULL )
1218 ERROR_exit("No RHS dataset input!?") ;
1219 ntime = DSET_NVALS(rhset) ;
1220 if( ntime < 2 )
1221 ERROR_exit("RHS dataset '%s' has only 1 value per voxel?!",rhsnam) ;
1222
1223 nx = DSET_NX(rhset); ny = DSET_NY(rhset); nz = DSET_NZ(rhset);
1224 nvox = nx*ny*nz;
1225
1226 if( (nvox == 1 || rhs1D) &&
1227 my_getenv("AFNI_TFITTER_VERBOSE") == NULL ) /* 31 Dec 2008 */
1228 AFNI_setenv("AFNI_TFITTER_VERBOSE=YES") ;
1229
1230 if( mask != NULL && (mnx != nx || mny != ny || mnz != nz) )
1231 ERROR_exit("mask and RHS datasets don't match in grid size") ;
1232
1233 if( nvar >= ntime && meth > 0 )
1234 ERROR_exit("too many (%d) LHS time series for %d time points",nvar,ntime) ;
1235
1236 npol = (polort < 0) ? 0 : polort+1 ; /* number of polynomial columns */
1237 if( npol > 0 ){
1238 ADDTO_XTARR(dsar,NULL);
1239 ii = XTARR_NUM(dsar)-1; XTARR_IC(dsar,ii) = IC_POLORT; nvar += npol;
1240 }
1241 nlhs = dsar->num ; /* number of -LHS inputs */
1242
1243 /* deconvolution stuff */
1244
1245 if( fal_klen > 0 ){
1246 if( fal_set != NULL ){
1247 if( DSET_NX(fal_set) != nx ||
1248 DSET_NY(fal_set) != ny || DSET_NZ(fal_set) != nz )
1249 ERROR_exit("-FALTUNG dataset and RHS dataset don't match in grid size");
1250 }
1251 if( fal_klen >= ntime )
1252 ERROR_exit("-FALTUNG kernel size %d must be shorter than dataset %d",
1253 fal_klen , ntime ) ;
1254 if( fal_klen >= ntime/2 )
1255 WARNING_message("-FALTUNG kernel size %d longer than 1/2 dataset %d",
1256 fal_klen , ntime/2 ) ;
1257 if( fal_set != NULL )
1258 fal_kern = (float *)malloc(sizeof(float)*fal_klen) ;
1259 nvoff = ntime ;
1260 } else if( nlhs == 0 ){
1261 ERROR_exit("no -LHS or -polort option given?!") ;
1262 }
1263
1264 if( nvar > 0 )
1265 rvec = (float **)malloc(sizeof(float *)*nvar ) ; /* LHS vectors */
1266
1267 /*--- check LHS inputs and assign ref vectors ---*/
1268
1269 if( nlhs > 0 ){
1270 MAKE_intvec(kvec,nlhs) ;
1271 for( kk=ii=0 ; ii < nlhs ; ii++ ){
1272 if( XTARR_IC(dsar,ii) == IC_FLIM ){ /* 1D image: ref points to data */
1273 float *lar ; int mm ;
1274 lim = (MRI_IMAGE *)XTARR_XT(dsar,ii) ;
1275 jj = lim->nx ; lar = MRI_FLOAT_PTR(lim) ;
1276 if( jj < ntime )
1277 ERROR_exit("LHS 1D file is too short along time axis: %d < %d",jj,ntime) ;
1278 if( jj > ntime )
1279 WARNING_message(
1280 "LHS 1D file too long along time axis: %d > %d: ignoring extra",jj,ntime);
1281 for( mm=0 ; mm < lim->ny ; mm++ ) rvec[kk++] = lar + mm*lim->nx ;
1282 kvec->ar[ii] = -1 ;
1283 } else if( XTARR_IC(dsar,ii) == IC_POLORT ){ /* polort columns */
1284 int mm ;
1285 for( mm=0 ; mm < npol ; mm++,kk++ ){
1286 rvec[kk] = (float *)malloc(sizeof(float)*ntime) ;
1287 for( jj=0 ; jj < ntime ; jj++ )
1288 rvec[kk][jj] = lhs_legendre( (float)jj, 0.0f, ntime-1.0f, mm ) ;
1289 }
1290 kvec->ar[ii] = -1 ;
1291 } else if( XTARR_IC(dsar,ii) == IC_DSET ){ /* dset: create ref vector */
1292 lset = (THD_3dim_dataset *)XTARR_XT(dsar,ii) ;
1293 if( DSET_NX(lset) != nx || DSET_NY(lset) != ny || DSET_NZ(lset) != nz )
1294 ERROR_exit("LHS dataset '%s' doesn't match RHS dataset grid size",
1295 DSET_HEADNAME(lset) ) ;
1296 jj = DSET_NVALS(lset) ;
1297 if( jj < ntime )
1298 ERROR_exit("LHS dataset is too short along time axis: %d < %d",jj,ntime) ;
1299 if( jj > ntime )
1300 WARNING_message(
1301 "LHS dataset too long along time axis: %d > %d: ignoring extra",jj,ntime);
1302 kvec->ar[ii] = kk ; /* index of vector from this dataset */
1303 rvec[kk++] = (float *)malloc(sizeof(float)*(jj+1)) ;
1304 } else {
1305 ERROR_exit("This message should never be seen by mortal eyes!") ;
1306 }
1307 }
1308 }
1309
1310 #if 0
1311 if( verb && nlset == 0 && meth == 2 && convec == NULL )
1312 INFO_message("LHS datasets all 1D files ==> you could use 3dDeconvolve");
1313 #endif
1314
1315 if( nlab < nvar ){
1316 char lll[32] ;
1317 if( verb && strcmp(prefix,"NULL") != 0 )
1318 INFO_message("Making up %d LHS labels (out of %d parameters)",nvar-nlab,nvar);
1319 lab = (char **)realloc((void *)lab,sizeof(char *)*nvar) ;
1320 for( ii=nlab ; ii < nvar ; ii++ ){
1321 sprintf(lll,"Param#%d",ii+1) ; lab[ii] = strdup(lll) ;
1322 }
1323 } else if( nlab > nvar && strcmp(prefix,"NULL") != 0 ){
1324 WARNING_message("Too many (%d) -label strings for %d parameters",nlab,nvar) ;
1325 }
1326
1327 /*--- create constraint vector ---*/
1328
1329 if( convec != NULL ){
1330 if( nvar > 0 ){
1331 cvec = (float *)calloc(sizeof(float),nvar) ;
1332 for( ii=0 ; ii < convec->nar ; ii++ ){
1333 kk = convec->ar[ii] ; jj = abs(kk) ;
1334 if( jj > nvar ) ERROR_exit("Index %d in -consign is too large",jj) ;
1335 cvec[jj-1] = (kk < 0) ? -1.0f : 1.0f ;
1336 }
1337 } else {
1338 WARNING_message("-consign option ignored: no -LHS given!") ;
1339 }
1340 }
1341
1342 /*----- load input datasets -----*/
1343
1344 if( verb && !rhs1D ) INFO_message("loading input datasets into memory") ;
1345
1346 DSET_load(rhset) ; CHECK_LOAD_ERROR(rhset) ;
1347 for( ii=0 ; ii < nlhs ; ii++ ){
1348 if( XTARR_IC(dsar,ii) == IC_DSET ){
1349 lset = (THD_3dim_dataset *)XTARR_XT(dsar,ii) ;
1350 DSET_load(lset) ; CHECK_LOAD_ERROR(lset) ;
1351 }
1352 }
1353
1354 /*------ create output datasets ------*/
1355
1356 if( nvar > 0 && strcmp(prefix,"NULL") != 0 ){ /** coefficients **/
1357 fset = EDIT_empty_copy(rhset) ;
1358 EDIT_dset_items( fset ,
1359 ADN_nvals , nvar ,
1360 ADN_ntt , 0 ,
1361 ADN_func_type , FUNC_BUCK_TYPE ,
1362 ADN_type , HEAD_FUNC_TYPE ,
1363 ADN_datum_all , MRI_float ,
1364 ADN_brick_fac , NULL ,
1365 ADN_prefix , prefix ,
1366 ADN_none ) ;
1367 tross_Copy_History( rhset , fset ) ;
1368 tross_Make_History( "3dTfitter" , argc,argv , fset ) ;
1369
1370 for( jj=0 ; jj < nvar ; jj++ ){ /* create empty bricks to be filled below */
1371 EDIT_substitute_brick( fset , jj , MRI_float , NULL ) ;
1372 EDIT_BRICK_LABEL( fset , jj , lab[jj] ) ;
1373 }
1374 }
1375
1376 if( fitts_prefix != NULL && strcmp(fitts_prefix,"NULL") != 0 ){ /** fitts */
1377 fitts_set = EDIT_empty_copy(rhset) ;
1378 EDIT_dset_items( fitts_set ,
1379 ADN_datum_all , MRI_float ,
1380 ADN_brick_fac , NULL ,
1381 ADN_prefix , fitts_prefix ,
1382 ADN_none ) ;
1383 tross_Copy_History( rhset , fitts_set ) ;
1384 tross_Make_History( "3dTfitter" , argc,argv , fitts_set ) ;
1385
1386 for( jj=0 ; jj < ntime ; jj++ ) /* create empty bricks to be filled below */
1387 EDIT_substitute_brick( fitts_set , jj , MRI_float , NULL ) ;
1388 }
1389
1390 if( ersum_prefix != NULL && strcmp(ersum_prefix,"NULL") != 0 ){ /** ersum */
1391 ersum_set = EDIT_empty_copy(rhset) ; /* 23 Jul 2009 */
1392 EDIT_dset_items( ersum_set ,
1393 ADN_datum_all , MRI_float ,
1394 ADN_brick_fac , NULL ,
1395 ADN_nvals , 2 ,
1396 ADN_ntt , 2 ,
1397 ADN_prefix , ersum_prefix ,
1398 ADN_none ) ;
1399 tross_Copy_History( rhset , ersum_set ) ;
1400 tross_Make_History( "3dTfitter" , argc,argv , ersum_set ) ;
1401
1402 EDIT_substitute_brick( ersum_set , 0 , MRI_float , NULL ) ;
1403 EDIT_substitute_brick( ersum_set , 1 , MRI_float , NULL ) ;
1404 EDIT_BRICK_LABEL ( ersum_set , 0 , "errsum_L2" ) ;
1405 EDIT_BRICK_LABEL ( ersum_set , 1 , "errsum_L1" ) ;
1406 }
1407
1408 if( fal_klen > 0 && strcmp(fal_pre,"NULL") != 0 ){ /** deconvolution **/
1409 defal_set = EDIT_empty_copy(rhset) ;
1410 EDIT_dset_items( defal_set ,
1411 ADN_datum_all , MRI_float ,
1412 ADN_brick_fac , NULL ,
1413 ADN_prefix , fal_pre ,
1414 ADN_none ) ;
1415 tross_Copy_History( rhset , defal_set ) ;
1416 tross_Make_History( "3dTfitter" , argc,argv , defal_set ) ;
1417
1418 for( jj=0 ; jj < ntime ; jj++ ) /* create empty bricks to be filled below */
1419 EDIT_substitute_brick( defal_set , jj , MRI_float , NULL ) ;
1420 }
1421
1422 /*-- Invert LASSO index selection? [04 Aug 2021] --*/
1423
1424 if( lasso_ivec != NULL && lasso_invert_ivec ){ /* 04 Aug 2021 */
1425 intvec *tmp_ivec , *new_ivec=NULL ;
1426
1427 /** DUMP_intvec(lasso_ivec,"LASSO - to be inverted") ; **/
1428
1429 MAKE_intvec( tmp_ivec , nvar ) ;
1430 for( jj=0 ; jj < nvar ; jj++ ) tmp_ivec->ar[jj] = jj+1 ;
1431
1432 for( jj=0 ; jj < lasso_ivec->nar ; jj++ ){
1433 ii = lasso_ivec->ar[jj] ;
1434 if( ii > 0 && ii <= nvar ) tmp_ivec->ar[ii-1] = 0 ; /* un-use it */
1435 }
1436
1437 for( ii=jj=0 ; jj < nvar ; jj++ ){
1438 if( tmp_ivec->ar[jj] > 0 ) ii++ ;
1439 }
1440 if( ii > 0 ){ /* copy survivors (if any) */
1441 MAKE_intvec( new_ivec , ii ) ;
1442 for( ii=jj=0 ; jj < nvar ; jj++ ){
1443 if( tmp_ivec->ar[jj] > 0 ) new_ivec->ar[ii++] = tmp_ivec->ar[jj] ;
1444 }
1445 }
1446 KILL_intvec(tmp_ivec) ; KILL_intvec(lasso_ivec) ; lasso_ivec = new_ivec ;
1447 }
1448
1449 /** if( lasso_ivec != NULL ) DUMP_intvec(lasso_ivec,"final LASSO un-penalized indexes") ; **/
1450
1451 /*-- finalize LASSO setup? --*/
1452
1453 if( meth == -2 || meth == -1 ){
1454 if( nvoff > 0 ){ /* deconvolution */
1455 float *lam = (float *)malloc(sizeof(float)*(nvar+nvoff)) ;
1456 for( ii=0 ; ii < nvar+nvoff ; ii++ )
1457 lam[ii] = (ii < nvoff) ? lasso_flam : 0.0f ;
1458 THD_lasso_setlamvec( nvar+nvoff , lam ) ; free(lam) ;
1459 if( verb && nvar > 0 && lasso_ivec != NULL )
1460 ININFO_message("-FALTUNG ==> All LHS and polort parameters are non-penalized") ;
1461 } else if( lasso_ivec != NULL ){ /* standard regression */
1462 float *lam = (float *)malloc(sizeof(float)*nvar) ;
1463 for( ii=0 ; ii < nvar ; ii++ ) lam[ii] = lasso_flam ;
1464 for( jj=0 ; jj < lasso_ivec->nar ; jj++ ){
1465 ii = lasso_ivec->ar[jj] - 1 ;
1466 if( ii >= 0 && ii < nvar ) lam[ii] = 0.0f ;
1467 }
1468 THD_lasso_setlamvec( nvar , lam ) ; free(lam) ;
1469 }
1470 }
1471
1472 /*------- loop over voxels and process them one at a time ---------*/
1473
1474 if( verb && nvox > 499 ) vstep = nvox / 50 ;
1475 if( vstep > 0 ) fprintf(stderr,"++ voxel loop: ") ;
1476
1477 THD_fitter_do_fitts ( do_fitts ) ; /* 05 Mar 2008 */
1478 THD_fitter_set_vthresh( vthresh ) ; /* 18 May 2010 */
1479
1480
1481 AFNI_OMP_START ;
1482 #pragma omp parallel if( nvox > 1 )
1483 { /* per-thread variables declared here */
1484 int ii , jj , kk ;
1485 float *dvec , *evec , **rrqvec=rvec , *falq_kern=fal_kern ;
1486 THD_3dim_dataset *lset ;
1487 floatvec *bfit ;
1488
1489 #ifdef USE_OMP
1490 #pragma omp master
1491 { int nthr = omp_get_num_threads() ;
1492 if( nthr > 1 ) fprintf(stderr,"[%d OpenMP threads]%s",nthr,(vstep>0)?" ":"\n") ;
1493 }
1494 #endif
1495
1496 /* allocate thread-specific arrays for voxel-dependent data */
1497
1498 #pragma omp critical
1499 { dvec = (float * )malloc(sizeof(float)*ntime) ; /* RHS vector */
1500 evec = (float * )malloc(sizeof(float)*ntime) ; /* RHS vector */
1501
1502 /* This next suff is a little complicated:
1503 rvec = LHS vectors, which are partly pre-loaded earlier
1504 (from .1D files and from -polort basis functions),
1505 and partly need to be loaded for each voxel from
1506 the LHS datasets (if any)
1507 Problem: With OpenMP, the voxelwise loading vectors need
1508 to be distinct to prevent inter-thread clobbers
1509 Solution: Make rrqvec[ii] the same pointer as rvec[ii]
1510 for ii that does NOT come from a LHS dataset
1511 (that is, where rvec[ii] is pre-loaded with data)
1512 and make a new rrqvec[ii] to be loaded from a
1513 LHS dataset when that is necessary
1514 Comment: This is ugly, but the easiest way to retrofit for OpenMP */
1515
1516 if( nlset > 0 && AO_nth > 1 ){ /* have at least one LHS dataset */
1517 /* and have more than 1 thread */
1518 rrqvec = (float **)malloc(sizeof(float *)*nvar ) ;
1519
1520 for( ii=0 ; ii < nvar ; ii++ ){ /* copy all vectors */
1521 rrqvec[ii] = rvec[ii] ;
1522 }
1523 for( ii=0 ; ii < nlhs ; ii++ ){ /* then overwrite a few of them */
1524 kk = kvec->ar[ii] ;
1525 if( kk >= 0 ){
1526 lset = (THD_3dim_dataset *)XTARR_XT(dsar,ii) ;
1527 jj = DSET_NVALS(lset) ;
1528 rrqvec[kk] = (float *)malloc(sizeof(float)*(jj+1)) ;
1529 }
1530 }
1531 }
1532
1533 /* Similarly, if the deconvolution kernel comes
1534 from a dataset, we need a per-voxel kernel vector.
1535 However, since there is only 1 deconvolution kernel
1536 (unlike having multiple LHS vectors), this is simpler. */
1537
1538 if( fal_klen > 0 && fal_set != NULL )
1539 falq_kern = (float *)malloc(sizeof(float)*fal_klen) ;
1540
1541 }
1542
1543 #pragma omp for
1544 for( ii=0 ; ii < nvox ; ii++ ){
1545
1546 if( vstep > 0 && ii%vstep==vstep-1 ) vstep_print() ;
1547
1548 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip */
1549
1550 THD_extract_array( ii , rhset , 0 , dvec ) ; /* get RHS data vector */
1551
1552 for( jj=0 ; jj < ntime && dvec[jj]==0.0f ; jj++ ) ; /*nada*/
1553 if( jj == ntime ){ /*** skip all zero vector ***/
1554 #pragma omp atomic /* ensure nskip ain't updated simultaneously in 2 threads */
1555 nskip++;
1556 continue;
1557 }
1558
1559 for( jj=0 ; jj < ntime ; jj++ ) evec[jj] = dvec[jj] ; /* copy vector */
1560
1561 THD_fitter_voxid(ii) ; /* 10 Sep 2008: for error messages */
1562
1563 for( jj=0 ; jj < nlhs ; jj++ ){ /* get LHS data vectors */
1564 if( XTARR_IC(dsar,jj) == IC_DSET ){ /* out of LHS datasets */
1565 lset = (THD_3dim_dataset *)XTARR_XT(dsar,jj) ;
1566 kk = kvec->ar[jj] ;
1567 THD_extract_array( ii , lset , 0 , rrqvec[kk] ) ;
1568 }
1569 }
1570
1571 /***** do the fitting work *****/
1572
1573 if( fal_klen > 0 ){ /*-- deconvolution --*/
1574
1575 if( fal_set != NULL ) /* get decon kernel if from a 3D+time dataset */
1576 THD_extract_array( ii , fal_set , 0 , falq_kern ) ;
1577
1578 bfit = THD_deconvolve( ntime , dvec ,
1579 0 , fal_klen-1 , falq_kern ,
1580 nvar , rrqvec , meth , cvec , fal_dcon ,
1581 fal_pencod , fal_penfac ) ;
1582
1583 } else { /*-- simple fitting --*/
1584
1585 bfit = THD_fitter( ntime , dvec , nvar , rrqvec , meth , cvec ) ;
1586
1587 }
1588
1589 if( bfit == NULL ){ nbad++; continue; } /*** bad voxel ***/
1590
1591 /** store the results **/
1592
1593 if( nvar > 0 && fset != NULL )
1594 THD_insert_series( ii , fset , nvar , MRI_float , bfit->ar+nvoff , 1 ) ;
1595
1596 if( fal_klen > 0 && defal_set != NULL )
1597 THD_insert_series( ii , defal_set , ntime , MRI_float , bfit->ar , 1 ) ;
1598
1599 if( do_fitts ){
1600 floatvec *fv = THD_retrieve_fitts() ; float *fvar = fv->ar ;
1601 if( fitts_set != NULL )
1602 THD_insert_series( ii , fitts_set , ntime , MRI_float , fvar , 1 ) ;
1603 if( ersum_set != NULL ){ /* 23 Jul 2009 */
1604 float qsum=0.0f , asum=0.0f , val , ee[2] ;
1605 for( jj=0 ; jj < ntime ; jj++ ){
1606 val = fabsf(fvar[jj]-evec[jj]) ; qsum += val*val ; asum += val ;
1607 }
1608 ee[0] = qsum ; ee[1] = asum ;
1609 THD_insert_series( ii , ersum_set , 2 , MRI_float , ee , 1 ) ;
1610 }
1611 }
1612
1613 KILL_floatvec(bfit) ; ngood++ ;
1614
1615 } /* end of loop over voxels */
1616
1617 /* If I was feeling kindly, I'd free the per-thread malloc-ed stuff now. */
1618 /* However, at this moment my kindliness coefficient has drizzled away. */
1619 }
1620 AFNI_OMP_END ;
1621
1622 if( vstep > 0 ) fprintf(stderr," Done!\n") ;
1623 if( nskip > 0 ) WARNING_message("Skipped %d voxels for being all zero",nskip) ;
1624
1625 /*----- clean up and go away -----*/
1626
1627 if( nbad > 0 )
1628 WARNING_message("Fit worked in %d voxels; failed in %d",ngood,nbad) ;
1629 else if( verb )
1630 INFO_message("Fit worked on all %d voxels attempted",ngood) ;
1631
1632 if( fset != NULL ){
1633 if( verb ) ININFO_message("Writing parameter dataset: %s",prefix) ;
1634 DSET_write(fset); DSET_unload(fset);
1635 }
1636
1637 if( defal_set != NULL ){
1638 if( verb ) ININFO_message("Writing FALTUNG dataset: %s",fal_pre) ;
1639 DSET_write(defal_set); DSET_unload(defal_set);
1640 }
1641
1642 if( fitts_set != NULL ){
1643 if( verb ) ININFO_message("Writing fitts dataset: %s",fitts_prefix) ;
1644 DSET_write(fitts_set); DSET_unload(fitts_set);
1645 }
1646
1647 if( ersum_set != NULL ){
1648 if( verb ) ININFO_message("Writing errsum dataset: %s",ersum_prefix) ;
1649 DSET_write(ersum_set); DSET_unload(ersum_set);
1650 }
1651
1652 if( verb )
1653 INFO_message("3dTfitter: CPU time = %.1f s Elapsed = %.1f s",COX_cpu_time(),COX_clock_time()) ;
1654 exit(0);
1655 }
1656
1657 /*---------------------------------------------------------------------------*/
1658
vstep_print(void)1659 static void vstep_print(void)
1660 {
1661 static int nn=0 ;
1662 static char xx[10] = "0123456789" ;
1663 fprintf(stderr , "%c" , xx[nn%10] ) ;
1664 if( nn%10 == 9) fprintf(stderr,".") ;
1665 nn++ ;
1666 }
1667
1668 /*---------------------------------------------------------------------------*/
1669 /* LHS Legendre polynomial */
1670
lhs_legendre(float x,float bot,float top,float n)1671 static float lhs_legendre( float x, float bot, float top, float n )
1672 {
1673 double xx ;
1674 xx = 2.0*(x-bot)/(top-bot) - 1.0 ; /* now in range -1..1 */
1675 return (float)Plegendre(xx,n) ;
1676 }
1677