00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 #define PROGRAM_NAME "plug_deconvolve"
00140 #define PROGRAM_AUTHOR "B. Douglas Ward"
00141 #define PROGRAM_INITIAL "09 September 1998"
00142 #define PROGRAM_LATEST "18 March 2003"
00143
00144
00145
00146 #include "afni.h"
00147 #include "matrix.h"
00148
00149 #define MAX_NAME_LENGTH THD_MAX_NAME
00150 #define MAX_XVARS 250
00151 #define MAX_STIMTS 20
00152 #define MAX_GLT 20
00153 #define MAX_CONSTR 20
00154
00155 #define RA_error DC_error
00156
00157
00158
00159
00160
00161
00162 static void DC_error (char * message)
00163 {
00164 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00165
00166 }
00167
00168
00169
00170
00171 #ifndef ALLOW_PLUGINS
00172 # error "Plugins not properly set up -- see machdep.h"
00173 #endif
00174
00175
00176
00177
00178 static char helpstring[] =
00179 " Purpose: Control DC_Fit, DC_Err, and DC_IRF Deconvolution Functions. \n"
00180 " \n"
00181 " Control Base = Baseline polynomial: None, Constant, Linear, \n"
00182 " Quadratic, Cubic, Quartic, or Quintic \n"
00183 " NFirst = Number of first time series point to use \n"
00184 " NLast = Number of last time series point to use \n"
00185 " IRF = Label of stimulus fnc. to use for IRF plot \n"
00186 " \n"
00187 " Concat Label = Name to use as label for concatenation \n"
00188 " File = File containing list of volume indices of \n"
00189 " starting points for individual runs \n"
00190 " Col # = Column of file which contains concat run list \n"
00191 " \n"
00192 " Censor Label = Name to use as label for censor function \n"
00193 " File = Time series file indicating censored points \n"
00194 " Col # = Column of file which contains censor function \n"
00195 " \n"
00196 " StimFnc Label = Name to use as label for this input stimulus \n"
00197 " File = Time series file representing input stimulus \n"
00198 " Col # = Column of file which contains input stimulus \n"
00199 " MinLag = Minimum time delay for impulse response \n"
00200 " MaxLag = Maximum time delay for impulse response \n"
00201 " NPTR = Number of stim fn. time points per TR \n"
00202 " Base = Is this input stimulus part of baseline model?\n"
00203 " \n"
00204 " GLT Mat Label = Name to use as label for this GLT matrix \n"
00205 " File = File containing the GLT matrix \n"
00206 " # Rows = Number of rows (linear constraints) in GLT \n"
00207 " \n"
00208 " \n"
00209 " For more details, see: \n"
00210 " `Deconvolution Analysis of FMRI Time Series Data' by B. Douglas Ward, \n"
00211 " which is contained in file 3dDeconvolve.ps of the AFNI distribution. \n"
00212 ;
00213
00214
00215
00216
00217 #define NBASE 7
00218 static char * baseline_strings[NBASE] = {"None", "Const", "Linear",
00219 "Quadrtc", "Cubic", "Quartic", "Quintic" };
00220
00221 static char * false_or_true[2] = {"False", "True"};
00222
00223
00224
00225 static char * DC_main( PLUGIN_interface * ) ;
00226
00227 static void DC_Fit (int nt, double to, double dt, float * vec, char ** label) ;
00228 static void DC_Err (int nt, double to, double dt, float * vec, char ** label) ;
00229 static void DC_IRF (int nt, double to, double dt, float * vec, char ** label) ;
00230 static int calculate_results();
00231
00232
00233
00234
00235
00236 static PLUGIN_interface * global_plint = NULL;
00237
00238 static int plug_polort;
00239 static int plug_p;
00240 static int plug_q;
00241 static int plug_qp;
00242 static int plug_NFirst;
00243 static int plug_NLast;
00244 static int plug_IRF;
00245 static int initialize;
00246 static int prev_nt;
00247 static char * IRF_label;
00248
00249 static char * concat_label;
00250 static int concat_column;
00251 static int num_blocks;
00252 static int * block_list;
00253
00254 static int num_censor;
00255 static char * censor_label;
00256 static int censor_column;
00257 static float * censor_array;
00258 static int censor_length;
00259 static int * good_list;
00260
00261 static int num_stimts;
00262 static int stim_base[MAX_STIMTS];
00263 static int stim_column[MAX_STIMTS];
00264 static float * stimulus[MAX_STIMTS];
00265 static int stim_length[MAX_STIMTS];
00266 static int min_lag[MAX_STIMTS];
00267 static int max_lag[MAX_STIMTS];
00268 static int nptr[MAX_STIMTS];
00269 static char * stim_label[MAX_STIMTS];
00270
00271 static matrix xdata;
00272 static matrix x_full;
00273 static matrix xtxinv_full;
00274 static matrix xtxinvxt_full;
00275 static matrix x_base;
00276 static matrix xtxinvxt_base;
00277 static matrix x_rdcd[MAX_STIMTS];
00278 static matrix xtxinvxt_rdcd[MAX_STIMTS];
00279
00280
00281 static int glt_num;
00282 static char * glt_label[MAX_GLT];
00283 static int glt_rows[MAX_GLT];
00284 static char * glt_filename[MAX_GLT];
00285 static matrix cxtxinvct[MAX_GLT];
00286 static matrix glt_cmat[MAX_GLT];
00287 static matrix glt_amat[MAX_GLT];
00288 static vector glt_coef[MAX_GLT];
00289 static vector glt_tcoef[MAX_GLT];
00290
00291
00292
00293
00294
00295
00296
00297 #undef USE_BASIS
00298
00299 #include "Deconvolve.c"
00300
00301
00302
00303
00304
00305
00306
00307
00308 static void initialize_options ()
00309 {
00310 int is;
00311 int iglt;
00312
00313
00314
00315 plug_polort = 1;
00316 plug_p = 0;
00317 plug_q = 0;
00318 plug_qp = 0;
00319 plug_NFirst = 0;
00320 plug_NLast = 32767;
00321 plug_IRF = -1;
00322 initialize = 0;
00323 prev_nt = 0;
00324 IRF_label = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (IRF_label);
00325 strcpy (IRF_label, " ");
00326
00327
00328
00329 concat_label = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (concat_label);
00330 strcpy (concat_label, " ");
00331 concat_column = 0;
00332 num_blocks = 1;
00333 block_list = (int *) malloc (sizeof(int) * 1); MTEST (block_list);
00334 block_list[0] = 0;
00335
00336
00337
00338 num_censor = 0;
00339 censor_label = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (censor_label);
00340 strcpy (censor_label, " ");
00341 censor_column = 0;
00342 censor_array = NULL;
00343 censor_length = 0;
00344 good_list = NULL;
00345
00346
00347
00348 num_stimts = 0;
00349 for (is =0; is < MAX_STIMTS; is++)
00350 {
00351 stim_label[is] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00352 MTEST (stim_label[is]);
00353 sprintf (stim_label[is], "Stim #%d ", is+1);
00354
00355 stim_base[is] = 0;
00356 stim_column[is] = 0;
00357 stimulus[is] = NULL;
00358 stim_length[is] = 0;
00359 min_lag[is] = 0;
00360 max_lag[is] = 0;
00361 nptr[is] = 1;
00362 }
00363
00364
00365
00366 matrix_initialize (&xdata);
00367 matrix_initialize (&x_full);
00368 matrix_initialize (&xtxinv_full);
00369 matrix_initialize (&xtxinvxt_full);
00370 matrix_initialize (&x_base);
00371 matrix_initialize (&xtxinvxt_base);
00372
00373 for (is =0; is < MAX_STIMTS; is++)
00374 {
00375 matrix_initialize (&x_rdcd[is]);
00376
00377 matrix_initialize (&xtxinvxt_rdcd[is]);
00378
00379 }
00380
00381
00382
00383 glt_num = 0;
00384 for (iglt =0; iglt < MAX_GLT; iglt++)
00385 {
00386 glt_label[iglt] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00387 MTEST (glt_label[iglt]);
00388 sprintf (glt_label[iglt], "GLT #%d ", iglt+1);
00389
00390 glt_rows[iglt] = 0;
00391 glt_filename[iglt] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00392 MTEST (glt_filename[iglt]);
00393 strcpy (glt_filename[iglt], " ");
00394
00395 matrix_initialize (&cxtxinvct[iglt]);
00396
00397 matrix_initialize (&glt_cmat[iglt]);
00398
00399 matrix_initialize (&glt_amat[iglt]);
00400
00401 vector_initialize (&glt_coef[iglt]);
00402
00403 vector_initialize (&glt_tcoef[iglt]);
00404
00405 }
00406
00407 }
00408
00409
00410
00411
00412
00413
00414
00415 static void reset_options ()
00416 {
00417 int is;
00418 int iglt;
00419
00420
00421
00422 plug_polort = 1;
00423 plug_p = 0;
00424 plug_q = 0;
00425 plug_qp = 0;
00426 plug_NFirst = 0;
00427 plug_NLast = 32767;
00428 plug_IRF = -1;
00429 initialize = 0;
00430 prev_nt = 0;
00431 strcpy (IRF_label, " ");
00432
00433
00434
00435 strcpy (concat_label, " ");
00436 concat_column = 0;
00437 num_blocks = 1;
00438 if (block_list != NULL) free (block_list);
00439 block_list = (int *) malloc (sizeof(int) * 1); MTEST (block_list);
00440 block_list[0] = 0;
00441
00442
00443
00444 num_censor = 0;
00445 strcpy (censor_label, " ");
00446 censor_column = 0;
00447 if (censor_array != NULL)
00448 {
00449 free (censor_array);
00450 censor_array = NULL;
00451 }
00452 censor_length = 0;
00453 if (good_list != NULL)
00454 {
00455 free (good_list);
00456 good_list = NULL;
00457 }
00458
00459
00460
00461 num_stimts = 0;
00462 for (is =0; is < MAX_STIMTS; is++)
00463 {
00464 sprintf (stim_label[is], "Stim #%d ", is+1);
00465
00466 stim_base[is] = 0;
00467 stim_column[is] = 0;
00468 if (stimulus[is] != NULL)
00469 {
00470 free (stimulus[is]);
00471 stimulus[is] = NULL;
00472 }
00473 stim_length[is] = 0;
00474 min_lag[is] = 0;
00475 max_lag[is] = 0;
00476 nptr[is] = 1;
00477 }
00478
00479
00480
00481 matrix_destroy (&xdata);
00482 matrix_destroy (&x_full);
00483 matrix_destroy (&xtxinv_full);
00484 matrix_destroy (&xtxinvxt_full);
00485 matrix_destroy (&x_base);
00486 matrix_destroy (&xtxinvxt_base);
00487
00488 for (is =0; is < MAX_STIMTS; is++)
00489 {
00490 matrix_destroy (&x_rdcd[is]);
00491
00492 matrix_destroy (&xtxinvxt_rdcd[is]);
00493
00494 }
00495
00496
00497
00498 glt_num = 0;
00499 for (iglt =0; iglt < MAX_GLT; iglt++)
00500 {
00501 sprintf (glt_label[iglt], "GLT #%d ", iglt+1);
00502
00503 glt_rows[iglt] = 0;
00504 strcpy (glt_filename[iglt], " ");
00505
00506 matrix_destroy (&cxtxinvct[iglt]);
00507
00508 matrix_destroy (&glt_cmat[iglt]);
00509
00510 matrix_destroy (&glt_amat[iglt]);
00511
00512 vector_destroy (&glt_coef[iglt]);
00513
00514 vector_destroy (&glt_tcoef[iglt]);
00515
00516 }
00517
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527 DEFINE_PLUGIN_PROTOTYPE
00528
00529 PLUGIN_interface * PLUGIN_init( int ncall )
00530 {
00531 PLUGIN_interface * plint ;
00532 int is;
00533 int iglt;
00534
00535
00536 if( ncall > 0 ) return NULL ;
00537
00538
00539
00540
00541
00542
00543 plint = PLUTO_new_interface ("Deconvolution" ,
00544 "Control DC_Fit, DC_Err, and DC_IRF Deconvolution Functions" ,
00545 helpstring, PLUGIN_CALL_VIA_MENU, DC_main);
00546
00547 global_plint = plint ;
00548
00549 PLUTO_short_choose(plint) ;
00550 PLUTO_short_number(plint) ;
00551
00552 PLUTO_add_hint (plint,
00553 "Control DC_Fit, DC_Err, and DC_IRF Deconvolution Functions");
00554
00555 PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00556
00557 PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;
00558
00559
00560 PLUTO_add_option (plint, "Control", "Control", TRUE);
00561 PLUTO_add_string (plint, "Base", NBASE, baseline_strings, 2);
00562 PLUTO_add_number (plint, "NFirst", -1, 32767, 0, -1, TRUE);
00563 PLUTO_add_number (plint, "NLast", 0, 32767, 0, 32767, TRUE);
00564 PLUTO_add_string( plint, "IRF ", 0, NULL, 1);
00565
00566
00567
00568 PLUTO_add_option (plint, "Concat", "Concat", FALSE);
00569 PLUTO_add_string( plint, "Label", 0, NULL, 1);
00570 PLUTO_add_timeseries (plint, "File");
00571 PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00572
00573
00574
00575 PLUTO_add_option (plint, "Censor", "Censor", FALSE);
00576 PLUTO_add_string( plint, "Label", 0, NULL, 1);
00577 PLUTO_add_timeseries (plint, "File");
00578 PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00579
00580
00581
00582 for (is = 0; is < MAX_STIMTS; is++)
00583 {
00584 PLUTO_add_option (plint, "StimFnc", "StimFnc", FALSE);
00585 PLUTO_add_string( plint, "Label", 0, NULL, 1);
00586 PLUTO_add_timeseries (plint, "File");
00587 PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00588 PLUTO_add_number (plint, "MinLag", 0, 100, 0, 0, TRUE);
00589 PLUTO_add_number (plint, "MaxLag", 0, 100, 0, 0, TRUE);
00590 PLUTO_add_number (plint, "NPTR", 1, 100, 0, 0, TRUE);
00591 PLUTO_add_string (plint, "Base", 2, false_or_true, 0);
00592 }
00593
00594
00595 for (is = 0; is < MAX_GLT; is++)
00596 {
00597 PLUTO_add_option (plint, "GLT Mat", "GLT Mat", FALSE);
00598 PLUTO_add_string( plint, "Label", 0, NULL, 1);
00599 PLUTO_add_string( plint, "File", 0, NULL, 1);
00600 PLUTO_add_number (plint, "# Rows", 1, MAX_CONSTR, 0, 0, TRUE);
00601 }
00602
00603
00604 PLUTO_register_1D_funcstr ("DC_Fit" , DC_Fit);
00605 PLUTO_register_1D_funcstr ("DC_Err" , DC_Err);
00606 PLUTO_register_1D_funcstr ("DC_IRF" , DC_IRF);
00607
00608
00609
00610 initialize_options ();
00611
00612
00613 return plint ;
00614 }
00615
00616
00617
00618
00619 static void show_options ()
00620 {
00621 int ib;
00622 int it;
00623 int is;
00624 int iglt;
00625
00626
00627
00628 printf ("\n\n");
00629 printf ("Program: %s \n", PROGRAM_NAME);
00630 printf ("Author: %s \n", PROGRAM_AUTHOR);
00631 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00632 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00633 printf ("\n");
00634
00635
00636
00637 printf ("\nControls: \n");
00638 printf ("Baseline = %10s \n", baseline_strings[plug_polort+1]);
00639 printf ("NFirst = %10d \n", plug_NFirst);
00640 printf ("NLast = %10d \n", plug_NLast);
00641 printf ("IRF label = %10s \n", IRF_label);
00642
00643
00644
00645 if (num_blocks > 0)
00646 {
00647 printf ("\n");
00648 printf ("Concatenation: Label = %8s ", concat_label);
00649 printf ("Column = %3d \n", concat_column);
00650 for (ib = 0; ib < num_blocks; ib++)
00651 printf ("Run #%d Initial Point = %d \n", ib+1, block_list[ib]);
00652 }
00653
00654
00655
00656 if (num_censor > 0)
00657 {
00658 printf ("\n");
00659 printf ("Censor Function: Label = %8s ", censor_label);
00660 printf ("Column = %3d \n", censor_column);
00661 printf ("Censored Points: ");
00662 for (it = 0; it < censor_length; it++)
00663 {
00664 if (censor_array[it] == 0) printf (" %d", it);
00665 }
00666 printf ("\n");
00667 }
00668
00669
00670
00671 if (num_stimts > 0)
00672 {
00673 printf ("\n");
00674 for (is = 0; is < num_stimts; is++)
00675 {
00676 if (stim_base[is])
00677 printf ("Baseline: Label = %8s ", stim_label[is]);
00678 else
00679 printf ("Stimulus: Label = %8s ", stim_label[is]);
00680 printf ("Column = %3d Min. Lag = %3d Max. Lag = %3d ",
00681 stim_column[is], min_lag[is], max_lag[is]);
00682 printf ("NPTR = %d \n", nptr[is]);
00683 }
00684 }
00685
00686
00687
00688 if (glt_num > 0)
00689 {
00690 printf ("\n");
00691 for (iglt = 0; iglt < glt_num; iglt++)
00692 {
00693 printf ("GLT: Label = %8s ", glt_label[iglt]);
00694 printf ("#Rows = %2d Input File: %s \n",
00695 glt_rows[iglt], glt_filename[iglt]);
00696 }
00697 }
00698
00699 }
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709 static char * DC_main( PLUGIN_interface * plint )
00710 {
00711 char * str;
00712 int is;
00713 int iglt;
00714 MRI_IMAGE * stim;
00715
00716 float * far;
00717 int ipt;
00718
00719
00720
00721 reset_options ();
00722
00723
00724
00725 PLUTO_next_option (plint);
00726 str = PLUTO_get_string (plint);
00727 plug_polort = PLUTO_string_index (str, NBASE, baseline_strings) - 1;
00728 plug_NFirst = PLUTO_get_number (plint);
00729 plug_NLast = PLUTO_get_number (plint);
00730 strcpy (IRF_label, PLUTO_get_string (plint));
00731
00732
00733
00734 str = PLUTO_peek_optiontag (plint);
00735 if (str != NULL)
00736
00737
00738 if (strcmp (str, "Concat") == 0)
00739 {
00740 str = PLUTO_get_optiontag (plint);
00741 strcpy (concat_label, PLUTO_get_string(plint));
00742
00743 stim = PLUTO_get_timeseries(plint) ;
00744
00745 if (stim == NULL || stim->nx < 1
00746 || stim->kind != MRI_float)
00747 return "**************************\n"
00748 "Illegal Concat File Input!\n"
00749 "**************************" ;
00750
00751
00752 concat_column = PLUTO_get_number(plint);
00753
00754 if ((concat_column < 0)
00755 ||(concat_column > stim->ny - 1))
00756 return "**********************************\n"
00757 "Illegal Concat File Column Number!\n"
00758 "**********************************" ;
00759
00760
00761 far = MRI_FLOAT_PTR(stim);
00762 num_blocks = stim->nx;
00763 if (block_list != NULL) free (block_list);
00764 block_list = (int *) malloc (sizeof(int) * num_blocks);
00765 MTEST (block_list);
00766
00767 for (ipt = 0; ipt < num_blocks; ipt++)
00768 block_list[ipt] = floor (far[ipt+concat_column*(stim->nx)] + 0.5);
00769 }
00770
00771
00772
00773 str = PLUTO_peek_optiontag (plint);
00774 if (str != NULL)
00775
00776
00777 if (strcmp (str, "Censor") == 0)
00778 {
00779 str = PLUTO_get_optiontag (plint);
00780 strcpy (censor_label, PLUTO_get_string(plint));
00781
00782 stim = PLUTO_get_timeseries(plint) ;
00783
00784 if (stim == NULL || stim->nx < 3
00785 || stim->kind != MRI_float)
00786 return "**************************\n"
00787 "Illegal Censor File Input!\n"
00788 "**************************" ;
00789
00790
00791
00792 censor_column = PLUTO_get_number(plint);
00793
00794 if ((censor_column < 0)
00795 ||(censor_column > stim->ny - 1))
00796 return "**********************************\n"
00797 "Illegal Censor File Column Number!\n"
00798 "**********************************" ;
00799
00800
00801
00802 if (censor_array != NULL)
00803 {
00804 free (censor_array);
00805 censor_array = NULL;
00806 }
00807 far = MRI_FLOAT_PTR(stim);
00808 censor_length = stim->nx;
00809 censor_array = (float *) malloc (sizeof(float) * (stim->nx));
00810 MTEST (censor_array);
00811
00812 num_censor = 1;
00813 for (ipt = 0; ipt < (stim->nx); ipt++)
00814 censor_array[ipt]
00815 = far[ipt + censor_column*(stim->nx)];
00816 }
00817
00818
00819
00820 do
00821 {
00822 str = PLUTO_get_optiontag(plint);
00823 if (str == NULL) break;
00824 if ((strcmp (str, "StimFnc") != 0) && (strcmp (str, "GLT Mat") != 0))
00825 return "************************\n"
00826 "Illegal optiontag found!\n"
00827 "************************";
00828
00829
00830
00831 if (strcmp (str, "StimFnc") == 0)
00832 {
00833 str = PLUTO_get_string(plint);
00834 if (strlen(str) != 0) strcpy (stim_label[num_stimts], str);
00835
00836 if (strcmp(stim_label[num_stimts], IRF_label) == 0)
00837 plug_IRF = num_stimts;
00838
00839 stim = PLUTO_get_timeseries(plint) ;
00840
00841 if (stim == NULL || stim->nx < 3
00842 || stim->kind != MRI_float)
00843 return "*************************\n"
00844 "Illegal Timeseries Input!\n"
00845 "*************************" ;
00846
00847
00848
00849 stim_column[num_stimts] = PLUTO_get_number(plint);
00850
00851 if ((stim_column[num_stimts] < 0)
00852 ||(stim_column[num_stimts] > stim->ny - 1))
00853 return "********************************\n"
00854 "Illegal Stim File Column Number!\n"
00855 "********************************" ;
00856
00857
00858
00859 if (stimulus[num_stimts] != NULL)
00860 {
00861 free (stimulus[num_stimts]);
00862 stimulus[num_stimts] = NULL;
00863 }
00864 far = MRI_FLOAT_PTR(stim);
00865 stim_length[num_stimts] = stim->nx;
00866 stimulus[num_stimts] = (float *) malloc (sizeof(float) * (stim->nx));
00867 MTEST (stimulus[num_stimts]);
00868
00869 for (ipt = 0; ipt < (stim->nx); ipt++)
00870 stimulus[num_stimts][ipt]
00871 = far[ipt + stim_column[num_stimts]*(stim->nx)];
00872
00873
00874
00875 min_lag[num_stimts] = PLUTO_get_number(plint);
00876 max_lag[num_stimts] = PLUTO_get_number(plint);
00877 nptr[num_stimts] = PLUTO_get_number(plint);
00878 str = PLUTO_get_string (plint);
00879 stim_base[num_stimts] = PLUTO_string_index (str, 2, false_or_true);
00880
00881
00882 if (min_lag[num_stimts] > max_lag[num_stimts])
00883 return "**************************\n"
00884 "Require Min Lag <= Max Lag\n"
00885 "**************************" ;
00886
00887 num_stimts++;
00888 }
00889
00890
00891
00892 if (strcmp (str, "GLT Mat") == 0)
00893 {
00894 str = PLUTO_get_string(plint);
00895 if (strlen(str) != 0) strcpy (glt_label[glt_num], str);
00896
00897 strcpy (glt_filename[glt_num], PLUTO_get_string(plint));
00898
00899 glt_rows[glt_num] = PLUTO_get_number(plint);
00900
00901 glt_num++;
00902 }
00903
00904 }
00905
00906 while (1);
00907
00908
00909
00910 plug_qp = (plug_polort + 1) * num_blocks;
00911 plug_q = plug_qp;
00912 plug_p = plug_qp;
00913 for (is = 0; is < num_stimts; is++)
00914 {
00915 if (stim_base[is]) plug_q += max_lag[is] - min_lag[is] + 1;
00916 plug_p += max_lag[is] - min_lag[is] + 1;
00917 if (plug_p > MAX_XVARS)
00918 {
00919 return "****************************\n"
00920 "Too many parameters in model \n"
00921 "****************************" ;
00922 }
00923 }
00924
00925
00926
00927 if (glt_num > 0)
00928 for (iglt = 0; iglt < glt_num; iglt++)
00929 {
00930 matrix_file_read (glt_filename[iglt],
00931 glt_rows[iglt],
00932 plug_p,
00933 &(glt_cmat[iglt]), 0);
00934 if (glt_cmat[iglt].elts == NULL)
00935 {
00936 return "**************************\n"
00937 "Unable to read GLT matrix \n"
00938 "**************************" ;
00939 }
00940 }
00941
00942
00943
00944 show_options ();
00945
00946
00947
00948 initialize = 1 ;
00949 prev_nt = 0;
00950
00951 return NULL ;
00952 }
00953
00954
00955
00956
00957
00958
00959
00960 static int calculate_results
00961 (
00962 int nt,
00963 double dt,
00964 float * vec,
00965 int * NN,
00966 int * nfit,
00967 float * fit,
00968 char ** label,
00969 float ** fitts,
00970 float ** errts
00971
00972 )
00973
00974 {
00975 float * ts_array;
00976
00977 int N;
00978 int p;
00979 int q;
00980 int qp;
00981 int m;
00982 int n;
00983
00984 vector coef;
00985 vector scoef;
00986 vector tcoef;
00987 float fpart[MAX_STIMTS];
00988 float rpart[MAX_STIMTS];
00989 float ffull;
00990 float rfull;
00991 float mse;
00992
00993 vector y;
00994
00995 int NFirst;
00996 int NLast;
00997 int i, it;
00998 int is;
00999
01000 float rms_min = 0.0;
01001 int ok;
01002 int novar;
01003 float fglt[MAX_GLT];
01004 float rglt[MAX_GLT];
01005
01006 int ib;
01007 int irb;
01008
01009
01010
01011 if (! initialize) return (0);
01012
01013
01014
01015 vector_initialize (&coef);
01016 vector_initialize (&scoef);
01017 vector_initialize (&tcoef);
01018 vector_initialize (&y);
01019
01020
01021
01022 qp = plug_qp;
01023 q = plug_q;
01024 p = plug_p;
01025 *nfit = p;
01026
01027
01028
01029 *fitts = (float *) malloc (sizeof(float) * nt); MTEST (*fitts);
01030 *errts = (float *) malloc (sizeof(float) * nt); MTEST (*errts);
01031
01032
01033
01034 if ((num_censor != 0) && (censor_length < nt))
01035 {
01036 DC_error ("Input censor time series file is too short");
01037 return (0);
01038 }
01039
01040
01041
01042 for (ib = 0; ib < num_blocks; ib++)
01043 if ((block_list[ib] < 0) || (block_list[ib] >= nt))
01044 {
01045 DC_error ("Invalid concatenated runs list");
01046 return (0);
01047 }
01048 if (num_blocks > 1)
01049 for (ib = 1; ib < num_blocks; ib++)
01050 if (block_list[ib] <= block_list[ib-1])
01051 {
01052 DC_error ("Invalid concatenated runs list");
01053 return (0);
01054 }
01055
01056
01057
01058 good_list = (int *) malloc (sizeof(int) * nt); MTEST (good_list);
01059 NFirst = plug_NFirst;
01060 if (NFirst < 0)
01061 for (is = 0; is < num_stimts; is++)
01062 if (NFirst < (max_lag[is]+nptr[is]-1)/nptr[is])
01063 NFirst = (max_lag[is]+nptr[is]-1)/nptr[is];
01064 NLast = plug_NLast;
01065
01066 N = 0;
01067 ib = 0;
01068 for (it = block_list[0]; it < nt; it++)
01069 {
01070 if (ib+1 < num_blocks)
01071 if (it >= block_list[ib+1]) ib++;
01072
01073 irb = it - block_list[ib];
01074
01075 if ((irb >= NFirst) && (irb <= NLast))
01076 if ((num_censor == 0) || (censor_array[it]))
01077 {
01078 good_list[N] = it;
01079 N++;
01080 }
01081 }
01082
01083 if (N == 0)
01084 {
01085 DC_error ("No usable time points?");
01086 return (0);
01087 }
01088 if (N <= p)
01089 {
01090 DC_error ("Insufficient time series data for deconvolution fit");
01091 return (0);
01092 }
01093
01094 *NN = N;
01095
01096
01097
01098 if (nt == prev_nt)
01099 {
01100 ok = 1;
01101 }
01102 else
01103 {
01104
01105 demean_base = (plug_polort >= 0) ;
01106 ok = init_indep_var_matrix (p, qp, plug_polort, nt, N, good_list,
01107 block_list, num_blocks, num_stimts, stimulus,
01108 stim_length, min_lag, max_lag, nptr, stim_base, &xdata);
01109
01110
01111
01112 if (ok)
01113 ok = init_regression_analysis (p, qp, num_stimts, stim_base, min_lag,
01114 max_lag, xdata, &x_full, &xtxinv_full, &xtxinvxt_full,
01115 &x_base, &xtxinvxt_base, x_rdcd, xtxinvxt_rdcd);
01116
01117
01118
01119 if (ok)
01120 if (glt_num > 0)
01121 ok = init_glt_analysis (xtxinv_full, glt_num, glt_cmat, glt_amat,
01122 cxtxinvct);
01123 }
01124
01125 if (ok)
01126 {
01127
01128 vector_create (N, &y);
01129 ts_array = vec;
01130 for (i = 0; i < N; i++)
01131 y.elts[i] = ts_array[good_list[i]];
01132
01133
01134
01135 regression_analysis (N, p, q, num_stimts, min_lag, max_lag,
01136 x_full, xtxinv_full, xtxinvxt_full, x_base,
01137 xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
01138 y, rms_min, &mse, &coef, &scoef, &tcoef,
01139 fpart, rpart, &ffull, &rfull, &novar,
01140 *fitts, *errts);
01141
01142
01143
01144 if (glt_num > 0)
01145 glt_analysis (N, p, x_full, y, mse*(N-p), coef, novar, cxtxinvct,
01146 glt_num, glt_rows, glt_cmat, glt_amat,
01147 glt_coef, glt_tcoef, fglt, rglt);
01148
01149
01150
01151 vector_to_array (coef, fit);
01152
01153
01154
01155 printf ("\nResults for Voxel: \n");
01156 report_results (N, qp, q, p, plug_polort, block_list, num_blocks,
01157 num_stimts, stim_label, stim_base, min_lag, max_lag,
01158 coef, tcoef, fpart, rpart, ffull, rfull, mse,
01159 glt_num, glt_label, glt_rows, glt_coef,
01160 glt_tcoef, fglt, rglt, label);
01161 printf ("%s \n", *label);
01162
01163 prev_nt = nt;
01164 }
01165
01166 else
01167 {
01168 vector_create (p, &coef);
01169 vector_to_array (coef, fit);
01170 strcpy (lbuf, "");
01171 *label = lbuf;
01172 prev_nt = 0;
01173 }
01174
01175
01176
01177 vector_destroy (&y);
01178 vector_destroy (&tcoef);
01179 vector_destroy (&scoef);
01180 vector_destroy (&coef);
01181
01182
01183 if (ok) return (1);
01184 else return (0);
01185 }
01186
01187
01188
01189
01190
01191
01192
01193 static void DC_Fit (int nt, double to, double dt, float * vec, char ** label)
01194
01195 {
01196 int N;
01197 int n;
01198 int ifit;
01199 int nfit;
01200 float fit[MAX_XVARS];
01201 int ok;
01202 float * fitts = NULL;
01203 float * errts = NULL;
01204
01205
01206
01207 ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01208 &fitts, &errts);
01209
01210
01211
01212 if (!ok)
01213 for (n = 0; n < nt; n++) vec[n] = 0.0;
01214
01215
01216
01217 else
01218 {
01219 for (n = 0; n < N; n++)
01220 vec[good_list[n]] = fitts[n];
01221 }
01222
01223
01224
01225 free (fitts); fitts = NULL;
01226 free (errts); errts = NULL;
01227
01228 return;
01229 }
01230
01231
01232
01233
01234
01235
01236
01237 static void DC_Err (int nt, double to, double dt, float * vec, char ** label)
01238
01239 {
01240 int N;
01241 float val;
01242 int n;
01243 int ifit;
01244 int nfit;
01245 float fit[MAX_XVARS];
01246 int ok;
01247 float * fitts = NULL;
01248 float * errts = NULL;
01249
01250
01251
01252 ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01253 &fitts, &errts);
01254
01255
01256
01257 for (n = 0; n < nt; n++) vec[n] = 0.0;
01258
01259
01260
01261 if (ok)
01262 {
01263 for (n = 0; n < N; n++)
01264 vec[good_list[n]] = errts[n];
01265 }
01266
01267
01268
01269 free (fitts); fitts = NULL;
01270 free (errts); errts = NULL;
01271
01272 return;
01273 }
01274
01275
01276
01277
01278
01279
01280
01281 static void DC_IRF (int nt, double to, double dt, float * vec, char ** label)
01282
01283 {
01284 int N;
01285 int nfit;
01286 float fit[MAX_XVARS];
01287 int np;
01288 int ip;
01289 int q;
01290 int it;
01291 int ntdnp;
01292 int ok;
01293 float * fitts = NULL;
01294 float * errts = NULL;
01295
01296
01297
01298 ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01299 &fitts, &errts);
01300
01301
01302
01303 if (!ok || (num_stimts < 1))
01304 for (it = 0; it < nt; it++) vec[it] = 0.0;
01305
01306
01307
01308 else
01309 {
01310 if ((num_stimts == 1) || (plug_IRF < 0) || (plug_IRF >= num_stimts))
01311 plug_IRF = 0;
01312
01313 np = max_lag[plug_IRF] - min_lag[plug_IRF] + 1;
01314
01315 q = plug_qp;
01316 for (ip = 0; ip < plug_IRF; ip++)
01317 q += max_lag[ip] - min_lag[ip] + 1;
01318
01319 if (np == 1)
01320 {
01321 for (it = 0; it < nt; it++)
01322 vec[it] = fit[q];
01323 }
01324 else
01325 {
01326 float r;
01327
01328 ntdnp = nt / (np-1);
01329
01330 vec[0] = fit[q];
01331 for (it = 0; it < nt; it++)
01332 {
01333 ip = it / ntdnp + 1;
01334 if (ip < np)
01335 {
01336 r = (float) it / (float) ntdnp - (ip-1);
01337 vec[it] = r * fit[q+ip] + (1.0-r) * fit[q+ip-1];
01338 }
01339 else
01340 vec[it] = fit[q+np-1];
01341 }
01342 }
01343 }
01344
01345
01346
01347 free (fitts); fitts = NULL;
01348 free (errts); errts = NULL;
01349
01350 return;
01351 }
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368