~ubuntu-branches/ubuntu/raring/voxbo/raring

« back to all changes in this revision

Viewing changes to crunch/vbmakeglm.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-06-06 11:33:11 UTC
  • Revision ID: james.westby@ubuntu.com-20100606113311-v3c13imdkkd5n7ae
Tags: upstream-1.8.5~svn1172
ImportĀ upstreamĀ versionĀ 1.8.5~svn1172

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
 
 
2
// vbmakeglm.cpp
 
3
// hack to create voxbo glm's
 
4
// Copyright (c) 1998-2010 by The VoxBo Development Team
 
5
 
 
6
// VoxBo is free software: you can redistribute it and/or modify it
 
7
// under the terms of the GNU General Public License as published by
 
8
// the Free Software Foundation, either version 3 of the License, or
 
9
// (at your option) any later version.
 
10
// 
 
11
// VoxBo is distributed in the hope that it will be useful, but
 
12
// WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
// General Public License for more details.
 
15
// 
 
16
// You should have received a copy of the GNU General Public License
 
17
// along with VoxBo.  If not, see <http://www.gnu.org/licenses/>.
 
18
// 
 
19
// For general information on VoxBo, including the latest complete
 
20
// source code and binary distributions, manual, and associated files,
 
21
// see the VoxBo home page at: http://www.voxbo.org/
 
22
//
 
23
// original version written by Dan Kimberg
 
24
// written to accomplish the same thing as some original IDL code by
 
25
// Geoff Aguirre
 
26
 
 
27
using namespace std;
 
28
 
 
29
#include <stdio.h>
 
30
#include <iostream>
 
31
#include <fstream>
 
32
#include <sstream>
 
33
 
 
34
#include "vbutil.h"
 
35
#include "vbjobspec.h"
 
36
#include "vbprefs.h"
 
37
#include "vbio.h"
 
38
#include "glmutil.h"
 
39
#include "vbx.h"
 
40
#include <sys/wait.h>
 
41
 
 
42
VBPrefs vbp;
 
43
 
 
44
const int BUFLEN=1024;
 
45
bool f_ignorewarnings,f_validateonly,f_setuponly;
 
46
int32 f_run;
 
47
string glmname;
 
48
 
 
49
class GLMParams {
 
50
public:
 
51
  string name;              // name of the sequence
 
52
  string dirname;           // root directory
 
53
  string stem;              // stem name for analysis
 
54
  vector<string>scanlist;   // list of scans (data files)
 
55
  int lows,highs;           // low and high frequencies to remove
 
56
  string middles;           // middle frequencies to remove
 
57
  int pieces;               // n pieces for matrix jobs
 
58
  string kernelname;        // temporal smoothing kernel
 
59
  double kerneltr;          // sampling interval of kernel file
 
60
  string noisemodel;        // noise model
 
61
  string refname;           // file containing reference function
 
62
  string glmfile;           // filename of .glm file
 
63
  string gmatrix;           // g matrix file
 
64
  string email;             // email address
 
65
  vector<string> contrasts; // list of contrasts to add to contrasts.txt
 
66
  int pri;
 
67
  bool auditflag,meannorm,emailflag,meannormset,driftcorrect,driftcorrectset;
 
68
  double TR;
 
69
  int orderg;
 
70
  bool valid;
 
71
  bool rfxgflag;            // should we just dummy up a G matrix with all 1's?
 
72
  VBSequence seq;
 
73
  void init();
 
74
  int RunCheck();
 
75
  int CreateGLMDir();
 
76
  void FixRelativePaths();
 
77
  void CreateGLMJobs();
 
78
  void CreateGLMJobs2();
 
79
  vector<string> CreateGLMScript();
 
80
  int WriteGLMFile(string filename="");
 
81
  int WriteAndSubmitJobs();
 
82
  void Validate();
 
83
  int createsamplefiles();
 
84
};
 
85
 
 
86
class GLMConfig {
 
87
private:
 
88
  GLMParams global,local,*gp;
 
89
public:
 
90
  GLMConfig(string fname,double tr,double kr);
 
91
  int ParseFile(string filename);
 
92
};
 
93
 
 
94
void vbmakeglm_help();
 
95
void vbmakeglm_sample(FILE *);
 
96
void vbmakeglm_make_sample(const string &fname);
 
97
void vbmakeglm_show_sample();
 
98
 
 
99
int
 
100
main(int argc,char **argv)
 
101
{
 
102
  tokenlist args;
 
103
  vbp.init();
 
104
  vbp.read_jobtypes();
 
105
  f_run=vbp.cores;
 
106
 
 
107
  args.Transfer(argc-1,argv+1);
 
108
  if (args.size() == 0) {
 
109
    vbmakeglm_help();
 
110
    exit(0);
 
111
  }
 
112
  if (args[0]=="-h") {
 
113
    vbmakeglm_show_sample();
 
114
    exit(0);
 
115
  }
 
116
  if (args[0]=="-x") {
 
117
    if (args.size() > 1)
 
118
      vbmakeglm_make_sample(args[1]);
 
119
    else
 
120
      printf("vbmakeglm -x needs a filename to store the sample glm file\n");
 
121
    exit(0);
 
122
  }
 
123
  
 
124
  // okay, let's fly
 
125
 
 
126
  f_ignorewarnings=0;
 
127
  f_validateonly=0;
 
128
  f_setuponly=0;
 
129
  vector<string> filelist;
 
130
  int i;
 
131
  double x_tr=0.0,x_kr=0.0;
 
132
 
 
133
  for (i=0; i<args.size(); i++) {
 
134
    if (args[i]=="-i") {
 
135
      f_ignorewarnings=1;
 
136
      continue;
 
137
    }
 
138
    else if (args[i]=="-v") {
 
139
      f_validateonly=1;
 
140
      continue;
 
141
    }
 
142
    else if (args[i]=="-s") {
 
143
      f_setuponly=1;
 
144
      continue;
 
145
    }
 
146
    else if (args[i]=="--run") {
 
147
      f_run=ncores()-1;
 
148
      if (f_run<1) f_run=1;
 
149
    }
 
150
    else if (args[i].compare(0,6,"--run=")==0) {
 
151
      f_run=strtol(args[i].substr(6));
 
152
      if (f_run<1) f_run=1;
 
153
    }
 
154
    else if (args[i]=="-n" && (i+1)<args.size()) {
 
155
      glmname=args[i+1];
 
156
      i++;
 
157
      continue;
 
158
    }
 
159
    else if (args[i]=="-tr" && (i+1)<args.size()) {
 
160
      x_tr=strtod(args[i+1]);
 
161
      i++;
 
162
      continue;
 
163
    }
 
164
    else if (args[i]=="-kr" && (i+1)<args.size()) {
 
165
      x_kr=strtod(args[i+1]);
 
166
      i++;
 
167
      continue;
 
168
    }
 
169
    else {
 
170
      filelist.push_back(args[i]);
 
171
    }
 
172
  }
 
173
  for (i=0; i<(int)filelist.size(); i++) {
 
174
    GLMConfig *glmc=new GLMConfig(filelist[i],x_tr,x_kr);
 
175
    delete glmc;
 
176
  }
 
177
  exit(0);
 
178
}
 
179
 
 
180
// this is the sole constructor, it chews on the passed arguments
 
181
 
 
182
GLMConfig::GLMConfig(string fname,double tr,double kr)
 
183
{
 
184
  local.init();
 
185
  global.init();
 
186
  gp=&global;
 
187
  global.glmfile=fname;
 
188
  global.TR=tr;
 
189
  global.kerneltr=kr;
 
190
  ParseFile(fname);
 
191
}
 
192
 
 
193
void
 
194
GLMParams::init()
 
195
{
 
196
  name="";
 
197
  dirname="";
 
198
  stem="";
 
199
  scanlist.clear();
 
200
  lows=highs=0;
 
201
  pieces=0;
 
202
  kernelname="";
 
203
  kerneltr=0.0;
 
204
  noisemodel="";
 
205
  refname="";
 
206
  gmatrix="";
 
207
  email="";
 
208
  pri=3;
 
209
  auditflag=1;
 
210
  meannorm=0;
 
211
  driftcorrect=0;
 
212
  emailflag=1;
 
213
  TR=0.0;
 
214
  orderg=0;
 
215
  valid=0;
 
216
  seq.init();
 
217
  rfxgflag=0;
 
218
}
 
219
 
 
220
void
 
221
GLMParams::Validate()
 
222
{
 
223
  int i,errs,warns;
 
224
  int timecount;
 
225
  Tes *mytes;
 
226
  stringstream tmps;
 
227
 
 
228
  timecount=0;
 
229
  errs=0;
 
230
  warns=0;
 
231
  valid=0;
 
232
 
 
233
  printf("[I] vbmakeglm: validating GLM %s...\n",name.c_str());
 
234
 
 
235
  if (!(dirname.size())) {
 
236
    printf("[E] vbmakeglm: no dirname specified for glm %s\n",name.c_str());
 
237
    errs++;
 
238
  }
 
239
 
 
240
  for (i=0; i<(int)scanlist.size(); i++) {
 
241
    mytes=new Tes;
 
242
    mytes->ReadHeader(scanlist[i].c_str());
 
243
    if (!(mytes->header_valid)) {
 
244
      printf("[E] vbmakeglm: couldn't get info on data file %s.\n",scanlist[i].c_str());
 
245
      errs++;
 
246
    }
 
247
    else if ((string)"tes1" != mytes->GetFileFormat().getSignature()) {
 
248
      printf("[E] vbmakeglm: couldn't read 4D data from file %s.\n",scanlist[i].c_str());
 
249
      errs++;
 
250
    }
 
251
    else if (mytes->GetFileFormat().getDimensions() != 4) {
 
252
      printf("[E] vbmakeglm: %s isn't a 4D file.\n",scanlist[i].c_str());
 
253
      errs++;
 
254
    }
 
255
    if (mytes->header_valid) {
 
256
      timecount += mytes->dimt;
 
257
      if (i==0 && TR<1.0) {
 
258
        if (mytes->voxsize[3]>0.0)
 
259
          TR=mytes->voxsize[3];
 
260
      }
 
261
    }
 
262
    delete mytes;
 
263
  }
 
264
  // the number of time points ought to be the same as orderg
 
265
  if (orderg == 0)
 
266
    orderg=timecount;
 
267
 
 
268
  if (lows < 0 || lows > (orderg/2)) {
 
269
    printf("[W] vbmakeglm: removing %d low frequencies is a little suspect\n",lows);
 
270
    warns++;
 
271
  }
 
272
 
 
273
  if (highs < 0 || highs > 20) {
 
274
    printf("[W] vbmakeglm: removing %d high frequencies is a little suspect\n",highs);
 
275
    warns++;
 
276
  }
 
277
 
 
278
  if (noisemodel.size()) {
 
279
    Vec p(noisemodel.c_str());
 
280
    if (p.size() != 3) {
 
281
      printf("[E] vbmakeglm: your 1/f parameter file %s has the wrong number of elements\n",xfilename(noisemodel).c_str());
 
282
      errs++;
 
283
    }
 
284
  }
 
285
 
 
286
  if (kernelname.size()) {
 
287
    Vec IRF;
 
288
    if (IRF.ReadFile(kernelname.c_str())) {
 
289
      tmps.str("");
 
290
      tmps << "vbmakeglm: your HRF kernel file " << kernelname << " doesn't exist.";
 
291
      printErrorMsg(VB_ERROR,tmps.str());
 
292
      errs++;
 
293
    }
 
294
    else if (IRF.size() < 6 || IRF.size() > 20) {
 
295
      tmps.str("");
 
296
      tmps << "vbmakeglm: your HRF kernel file " << kernelname << " has "
 
297
           << IRF.size() << " elements, which seems a little suspect.";
 
298
      printErrorMsg(VB_WARNING,tmps.str());
 
299
      warns++;
 
300
    }
 
301
    else if (kerneltr < 1.0) {
 
302
      string trline=GetHeader(IRF.header,"TR(msecs)");
 
303
      if (trline.size()) {
 
304
        tokenlist tmpl;
 
305
        tmpl.ParseLine(trline);
 
306
        if (tmpl.size()>1)
 
307
          kerneltr=strtod(tmpl[1]);
 
308
      }
 
309
    }
 
310
  }
 
311
 
 
312
  // load the G matrix, check its order and see if it has an intercept of interest
 
313
  if (gmatrix.size()) {
 
314
    VBMatrix gmat(gmatrix);
 
315
    // gmat.MakeInCore();
 
316
    if (gmat.m <=0 || gmat.n <= 0) {
 
317
      tmps.str("");
 
318
      tmps << "vbmakeglm: couldn't read G (design) matrix " << gmatrix << ".";
 
319
      printErrorMsg(VB_ERROR,tmps.str());
 
320
      errs++;
 
321
    }
 
322
    orderg=gmat.m;
 
323
    
 
324
    if (gmat.n==1 && meannorm) {
 
325
      tmps.str("");
 
326
      tmps << "vbmakeglm: you have a single covariate and the mean norm flag set - make sure you don't "
 
327
           << "mean normalize data for second tier (group rfx) analyses.";
 
328
      printErrorMsg(VB_WARNING,tmps.str());
 
329
      warns++;
 
330
    }
 
331
  }
 
332
  
 
333
  if (rfxgflag && meannorm) {
 
334
    tmps.str("");
 
335
    tmps << "vbmakeglm: you have the makerrandfxg flag and the mean norm flag set - make sure you don't "
 
336
         << "mean normalize data for second tier (group rfx) analyses.";
 
337
    printErrorMsg(VB_WARNING,tmps.str());
 
338
    warns++;
 
339
  }
 
340
 
 
341
  if (meannorm && noisemodel.size()==0) {
 
342
    tmps.str("");
 
343
    tmps << "vbmakeglm: you have the mean norm flag set, and no noise model - make sure you really want"
 
344
         << "to mean normalize these data.";
 
345
    printErrorMsg(VB_WARNING,tmps.str());
 
346
    warns++;
 
347
  }
 
348
 
 
349
  if (!meannormset) {
 
350
    tmps.str("");
 
351
    tmps << "vbmakeglm: no meannorm flag set -- defaulting to no mean normalization";
 
352
    printErrorMsg(VB_WARNING,tmps.str());
 
353
  }
 
354
 
 
355
  if (orderg != timecount) {
 
356
    printf("[E] vbmakeglm: orderg (%d) doesn't match the number of timepoints (%d)\n",
 
357
           orderg,timecount);
 
358
    errs++;
 
359
  }
 
360
 
 
361
  // FIXME - should actually check order of G matrix
 
362
  if (gmatrix.size()==0 && !rfxgflag) {
 
363
    printf("[E] vbmakeglm: no valid G (design) matrix specified\n");
 
364
    errs++;
 
365
  }
 
366
 
 
367
  if (!errs)
 
368
    valid=1;
 
369
  if (warns && !f_ignorewarnings)
 
370
    valid=0;
 
371
 
 
372
  if (valid) {
 
373
    tmps.str("");
 
374
    tmps << "vbmakeglm: GLM " << name << " is good to go.";
 
375
    printErrorMsg(VB_INFO,tmps.str());
 
376
  }
 
377
  if (TR<1.0) {
 
378
    TR=2000;
 
379
    tmps.str("");
 
380
    tmps << "vbmakeglm: TR not set, using default of 2000ms.";
 
381
    printErrorMsg(VB_WARNING,tmps.str());
 
382
  }
 
383
  if (kerneltr<1.0 && kernelname.size()) {
 
384
    kerneltr=2000;
 
385
    tmps.str("");
 
386
    tmps << "vbmakeglm: HRF TR (sampling rate) not set, using default of 2000ms.";
 
387
    printErrorMsg(VB_WARNING,tmps.str());
 
388
  }
 
389
  // the following added to double-check on the TR bugs
 
390
  printf("[I] vbmakeglm: your data TR is %f\n",TR);
 
391
  printf("[I] vbmakeglm: your HRF kernel TR (sampling rate) is %f\n",kerneltr);
 
392
}
 
393
 
 
394
int
 
395
GLMParams::CreateGLMDir()
 
396
{
 
397
  string fname;   // temporary variable for wherever we need to build a filename
 
398
  int i;
 
399
  struct stat st;
 
400
  stringstream tmps;
 
401
 
 
402
  tokenlist stemstuff;
 
403
  stemstuff.SetSeparator(" \t\n\r/");
 
404
  stemstuff.ParseLine(dirname.c_str());
 
405
  if (!stemstuff.size()) {
 
406
    tmps << "vbmakeglm: couldn't find parent of GLM directory" << dirname << ".";
 
407
    printErrorMsg(VB_ERROR,tmps.str());
 
408
    return 101;
 
409
  }
 
410
  stem=dirname+(string)"/"+stemstuff[stemstuff.size()-1];
 
411
 
 
412
  // create directory
 
413
  mkdir(dirname.c_str(),0777);
 
414
  if (stat(dirname.c_str(),&st)) {
 
415
    tmps << "vbmakeglm: couldn't get to GLM directory " << dirname << ".";
 
416
    printErrorMsg(VB_ERROR,tmps.str());
 
417
    return 102;
 
418
  }
 
419
  // create sub file
 
420
  fname=stem+".sub";
 
421
  ofstream subfile(fname.c_str());
 
422
  if (!subfile) {
 
423
    tmps << "vbmakeglm: couldn't create sub file " << fname << ".";
 
424
    printErrorMsg(VB_ERROR,tmps.str());
 
425
    return 103;
 
426
  }
 
427
  subfile << ";VB98\n;TXT1\n;\n; Subject list generated by vbmakeglm\n;\n\n";
 
428
  for (i=0; i<(int)scanlist.size(); i++)
 
429
    subfile << scanlist[i] << endl;
 
430
  subfile.close();
 
431
  
 
432
  // create ref
 
433
  if (refname.size()) {
 
434
    if (copyfile(refname,stem+".ref")) {
 
435
      tmps << "vbmakeglm: couldn't copy reference function " << refname << ".";
 
436
      printErrorMsg(VB_ERROR,tmps.str());
 
437
      return 104;
 
438
    }
 
439
  }
 
440
 
 
441
  // create glm file for potential future use
 
442
  if (glmfile.size()) {
 
443
    if (WriteGLMFile(stem+".glm")) {
 
444
      tmps << "vbmakeglm: couldn't write glm configuration file " << (stem+".glm") << ".";
 
445
      printErrorMsg(VB_WARNING,tmps.str());
 
446
      tmps.str("");
 
447
    }
 
448
  }
 
449
 
 
450
  // copy G if necessary
 
451
  if (gmatrix.size() > 0) {
 
452
    if (copyfile(gmatrix,stem+".G")) {
 
453
      printf("[E] vbmakeglm: couldn't copy G matrix file %s.\n",gmatrix.c_str());
 
454
      return 105;
 
455
    }
 
456
  }
 
457
  else if (rfxgflag) {
 
458
    gmatrix=stem+".G";
 
459
    ofstream gstr(gmatrix.c_str(),ios::binary);
 
460
    if (gstr) {
 
461
      gstr << "VB98\nMAT1\n";
 
462
      gstr << "DataType:\tFloat\n";
 
463
      gstr << "VoxDims(XY):\t1\t" << orderg << endl << endl;
 
464
      gstr << "# This G matrix generated automatically by vbmakeglm\n\n";
 
465
      gstr << "Parameter:\t0\tInterest\tEffect\n";
 
466
      gstr << "\x0c\n";
 
467
      float pts[orderg];
 
468
      for (i=0; i<orderg; i++)
 
469
        pts[i]=1.0;
 
470
      if (my_endian() != ENDIAN_BIG)
 
471
        swap(pts,orderg);
 
472
      for (i=0; i<(int)(orderg * sizeof(float)); i++) {
 
473
        gstr << *((unsigned char *)pts+i);
 
474
      }
 
475
      gstr.close();
 
476
    }
 
477
  }
 
478
  createsamplefiles();   // after G matrix, so we can count vars of interest
 
479
  return 0;   // no error!
 
480
}
 
481
 
 
482
int
 
483
GLMParams::createsamplefiles()
 
484
{
 
485
  // first find variables of interest
 
486
  GLMInfo glmi;
 
487
  glmi.stemname=stem;
 
488
  glmi.getcovariatenames();
 
489
  string fname=dirname+"/contrasts.txt";
 
490
  vector<string> interestnames;
 
491
  if (access(fname.c_str(),R_OK) || contrasts.size()) {
 
492
    ofstream ofile(fname.c_str());
 
493
    if (ofile) {
 
494
      ofile << "# contrasts.txt\n";
 
495
      ofile << "# in this file you can define contrasts among your covariates of interest\n";
 
496
      if (glmi.cnames.size()) {
 
497
        ofile << "# your covariates of interest are:\n";
 
498
        for (size_t i=0; i<glmi.cnames.size(); i++) {
 
499
          if (glmi.cnames[i][0]=='I') {
 
500
            ofile << "#   "<<strnum(i)<<": "<<glmi.cnames[i].c_str()+1<<endl;
 
501
            interestnames.push_back(glmi.cnames[i].substr(1));
 
502
          }
 
503
        }
 
504
      }
 
505
      ofile << "# you can specify a complete contrast as follows:\n#\n";
 
506
      ofile << "# <name> <scale> vec";  // 0 0 1 0\n";
 
507
      ofile << " 1";
 
508
      for (size_t i=1; i<interestnames.size(); i++)
 
509
        ofile << " 0";
 
510
      ofile << endl << "#\n";
 
511
      ofile << "# (with one value for each covariate of interest)\n";
 
512
      ofile << "#\n";
 
513
      ofile << "# lines beginning with a '#' are comments!\n";
 
514
      ofile << "#\n";
 
515
      ofile << "# the following simple contrasts are provided for your convenience:\n";
 
516
      ofile << endl;
 
517
      for (size_t i=0; i<interestnames.size(); i++) {
 
518
        ofile << interestnames[i] << " t vec";
 
519
        for (size_t j=0; j<interestnames.size(); j++) {
 
520
          if (j==i) ofile << " 1";
 
521
          else ofile << " 0";
 
522
        }
 
523
        ofile << endl;
 
524
      }
 
525
      if (contrasts.size()) {
 
526
        ofile << "\n# the following contrasts were specified:\n";
 
527
        ofile << endl;
 
528
        for (size_t i=0; i<contrasts.size(); i++) {
 
529
          if (glmi.parsecontrast(contrasts[i]))
 
530
            printf("[W] vbgmakeglm: couldn't parse contrast: %s\n",contrasts[i].c_str());
 
531
          else
 
532
            ofile << contrasts[i] << endl;
 
533
        }
 
534
      }
 
535
 
 
536
      ofile.close();
 
537
    }
 
538
  }
 
539
  fname=dirname+"/averages.txt";
 
540
  if (access(fname.c_str(),R_OK)) {
 
541
    ofstream ofile(fname.c_str());
 
542
    if (ofile) {
 
543
      ofile << "# averages.txt\n";
 
544
      ofile << "# \n";
 
545
      ofile << "# In this file you can specify one or more ways to trial-average your data.\n";
 
546
      ofile << "# You can also block-average or whatever else you need, we just call it\n";
 
547
      ofile << "# trial averaging generically.\n";
 
548
      ofile << "# \n";
 
549
      ofile << "# Each trial average needs a separate section that looks like the following:\n";
 
550
      ofile << "# \n";
 
551
      ofile << "# average <name>\n";
 
552
      ofile << "#   units <time/vols>\n";
 
553
      ofile << "#   interval <ms/vols>\n";
 
554
      ofile << "#   nsamples <n>\n";
 
555
      ofile << "#   tr <ms>\n";
 
556
      ofile << "#   trial <n>...\n";
 
557
      ofile << "#   trialset <first> <interval> <count>\n";
 
558
      ofile << "# end\n";
 
559
      ofile << "# \n";
 
560
      ofile << "# Here's what these parameters mean:\n";
 
561
      ofile << "# \n";
 
562
      ofile << "# units: whether the other parameters are in volumes or seconds\n";
 
563
      ofile << "# interval: interval in time or volumes between samples within the trial\n";
 
564
      ofile << "# nsamples: number of time points to include per trial\n";
 
565
      ofile << "# tr: sets the TR if you're using time units\n";
 
566
      ofile << "#\n";
 
567
      ofile << "# The remaining options are two ways to indicate when the trials begin.\n";
 
568
      ofile << "# If your trials are evenly spaced, use 'trialset,' otherwise use 'trial'\n";
 
569
      ofile << "#\n";
 
570
      ofile << "# trialset: specify the start of the first trial, the interval between trial\n";
 
571
      ofile << "#     onsets, and the trial count\n";
 
572
      ofile << "# trial: each trial line lists one or more start times/vols for a trial\n";
 
573
      ofile << "#     (you can include multiple trial lines to help you keep the file neat)\n";
 
574
      ofile << "#\n";
 
575
      ofile << "# --> for trial and trialsets, the first volume is volume 0 (also time 0)\n";
 
576
      ofile << "# --> both time and volume values can be floating point\n";
 
577
      ofile << "#\n";
 
578
      ofile << "# Total data points for this GLM: "<<orderg<<endl;
 
579
      ofile << "# Your TR in ms: "<<TR<<endl;
 
580
      ofile << "# \n";
 
581
 
 
582
      ofile.close();
 
583
    }
 
584
  }
 
585
  return 0;
 
586
}
 
587
 
 
588
int
 
589
GLMConfig::ParseFile(string filename)
 
590
{
 
591
  ifstream infile;
 
592
  char buf[BUFLEN];
 
593
  tokenlist args;
 
594
  stringstream tmps;
 
595
 
 
596
  infile.open(filename.c_str());
 
597
  if (!infile) {
 
598
    tmps << "vbmakeglm: couldn't open " << filename;
 
599
    printErrorMsg(VB_ERROR,tmps.str());
 
600
    return 1;  // error!
 
601
  }
 
602
  while (infile.getline(buf,BUFLEN,'\n')) {
 
603
    args.ParseLine(buf);
 
604
    if (args.size()) {
 
605
      if (args[0][0]=='#' || args[0][0] == ';')
 
606
        continue;
 
607
      if (args[0]=="name" || args[0]=="glm") {
 
608
        local=global;  // start with a copy of the global
 
609
        gp=&local;
 
610
        gp->name=args.Tail();
 
611
      }
 
612
      else if (args[0]=="dirname") {
 
613
        if (args.size() == 2)
 
614
          gp->dirname=args[1];
 
615
        else {
 
616
          tmps.str("");
 
617
          tmps << "vbmakeglm: dirname takes exactly one argument";
 
618
          printErrorMsg(VB_ERROR,tmps.str());
 
619
        }
 
620
      }
 
621
      else if (args[0]=="scan") {
 
622
        if (args.size() != 2) {
 
623
          tmps.str("");
 
624
          tmps << "vbmakeglm: scan takes exactly one argument";
 
625
          printErrorMsg(VB_ERROR,tmps.str());
 
626
          continue;
 
627
        }
 
628
        gp->scanlist.push_back(args[1]);
 
629
      }
 
630
      else if (args[0]=="lows") {
 
631
        if (args.size() == 2)
 
632
          gp->lows=strtol(args[1]);
 
633
        else {
 
634
          tmps.str("");
 
635
          tmps << "vbmakeglm: lows takes exactly one argument";
 
636
          printErrorMsg(VB_ERROR,tmps.str());
 
637
        }
 
638
      }
 
639
      else if (args[0]=="highs") {
 
640
        if (args.size() == 2)
 
641
          gp->highs=strtol(args[1]);
 
642
        else {
 
643
          tmps.str("");
 
644
          tmps << "vbmakeglm: highs takes exactly one argument";
 
645
          printErrorMsg(VB_ERROR,tmps.str());
 
646
        }
 
647
      }
 
648
      else if (args[0]=="middles") {
 
649
        if (args.size() > 0)
 
650
          gp->middles=args[1];
 
651
      }
 
652
      else if (args[0]=="pieces") {
 
653
        if (args.size() == 2)
 
654
          gp->pieces=strtol(args[1]);
 
655
        else {
 
656
          printf("[E] vbmakeglm: pieces takes exactly one argument\n");
 
657
        }
 
658
      }
 
659
      else if (args[0]=="orderg") {
 
660
        if (args.size() == 2)
 
661
          gp->orderg=strtol(args[1]);
 
662
        else {
 
663
          tmps.str("");
 
664
          tmps << "vbmakeglm: orderg takes exactly one argument";
 
665
          printErrorMsg(VB_ERROR,tmps.str());
 
666
        }
 
667
      }
 
668
      else if (args[0]=="refname") {
 
669
        if (args.size() == 2)
 
670
          gp->refname=args[1];
 
671
        else {
 
672
          printf("[E] vbmakeglm: reference takes exactly one argument\n");
 
673
          printErrorMsg(VB_ERROR,tmps.str());
 
674
        }
 
675
      }
 
676
      else if (args[0]=="kernel") {
 
677
        if (args.size() == 2)
 
678
          gp->kernelname=args[1];
 
679
        else if (args.size()==3) {
 
680
          gp->kernelname=args[1];
 
681
          gp->kerneltr=strtod(args[2]);
 
682
        }
 
683
        else {
 
684
          tmps.str("");
 
685
          tmps << "vbmakeglm: kernel takes one or two arguments";
 
686
          printErrorMsg(VB_ERROR,tmps.str());
 
687
        }
 
688
      }
 
689
      else if (args[0]=="noisemodel") {
 
690
        if (args.size() == 2)
 
691
          gp->noisemodel=args[1];
 
692
        else {
 
693
          tmps.str("");
 
694
          tmps << "vbmakeglm: noisemodel takes exactly one argument";
 
695
          printErrorMsg(VB_ERROR,tmps.str());
 
696
        }
 
697
      }
 
698
      else if (args[0]=="gmatrix") {
 
699
        if (args.size() == 2)
 
700
          gp->gmatrix=args[1];
 
701
        else {
 
702
          tmps.str("");
 
703
          tmps << "vbmakeglm: gmatrix takes exactly one argument"; 
 
704
          printErrorMsg(VB_ERROR,tmps.str());
 
705
        }
 
706
      }
 
707
      else if (args[0]=="makerandfxg") {
 
708
        gp->rfxgflag=1;
 
709
        if (args.size() > 1) {
 
710
          tmps.str("");
 
711
          tmps << "vbmakeglm: arguments to makerandfxg ignored";
 
712
          printErrorMsg(VB_WARNING,tmps.str());
 
713
        }
 
714
      }
 
715
      else if (args[0]=="tr") {
 
716
        if (args.size() == 2)
 
717
          gp->TR=strtol(args[1]);
 
718
        else {
 
719
          tmps.str("");
 
720
          tmps << "vbmakeglm: tr takes exactly one argument";
 
721
          printErrorMsg(VB_ERROR,tmps.str());
 
722
        }
 
723
      }
 
724
      else if (args[0]=="pri") {
 
725
        if (args.size() == 2)
 
726
          gp->pri=strtol(args[1]);
 
727
        else {
 
728
          tmps.str("");
 
729
          tmps << "vbmakeglm: pri takes exactly one argument";
 
730
          printErrorMsg(VB_ERROR,tmps.str());
 
731
        }
 
732
      }
 
733
      else if (args[0]=="audit") {
 
734
        if (args.size() != 2) {
 
735
          tmps.str("");
 
736
          tmps << "vbmakeglm: audit takes exactly one argument";
 
737
          printErrorMsg(VB_ERROR,tmps.str());
 
738
        }
 
739
        if (args[1]=="yes" || args[1]=="1")
 
740
          gp->auditflag=1;
 
741
        else if (args[1]=="no" || args[1]=="0")
 
742
          gp->auditflag=0;
 
743
        else {
 
744
          tmps.str("");
 
745
          tmps << "vbmakeglm: unrecognized value for audit flag (should be yes/no)";
 
746
          printErrorMsg(VB_ERROR,tmps.str());
 
747
        }
 
748
      }
 
749
      else if (args[0]=="meannorm") {
 
750
        if (args.size() != 2) {
 
751
          tmps.str("");
 
752
          tmps << "vbmakeglm: meannorm takes exactly one argument";
 
753
          printErrorMsg(VB_ERROR,tmps.str());
 
754
        }
 
755
        if (args[1]=="yes" || args[1]=="1") {
 
756
          gp->meannorm=1;
 
757
          gp->meannormset=1;
 
758
        }
 
759
        else if (args[1]=="no" || args[1]=="0") {
 
760
          gp->meannorm=0;
 
761
          gp->meannormset=1;
 
762
        }
 
763
      }
 
764
      else if (args[0]=="driftcorrect") {
 
765
        if (args.size() != 2) {
 
766
          tmps.str("");
 
767
          tmps << "vbmakeglm: driftcorrect takes exactly one argument";
 
768
          printErrorMsg(VB_ERROR,tmps.str());
 
769
        }
 
770
        if (args[1]=="yes" || args[1]=="1") {
 
771
          gp->driftcorrect=1;
 
772
          gp->driftcorrectset=1;
 
773
        }
 
774
        else if (args[1]=="no" || args[1]=="0") {
 
775
          gp->driftcorrect=0;
 
776
          gp->driftcorrectset=1;
 
777
        }
 
778
        else {
 
779
          tmps.str("");
 
780
          tmps << "vbmakeglm: unrecognized value for driftcorrect flag (should be yes/no)";
 
781
          printErrorMsg(VB_ERROR,tmps.str());
 
782
        }
 
783
      }
 
784
      else if (args[0]=="email") {
 
785
        gp->emailflag=1;
 
786
        if (args.size() == 2)
 
787
          gp->email=args[1];
 
788
      }
 
789
      else if (args[0]=="go" || args[0]=="end") {
 
790
        if (!glmname.size() || glmname==gp->name) {
 
791
          gp->FixRelativePaths();
 
792
          gp->Validate();
 
793
          if (gp->valid && !f_validateonly) {
 
794
            if (!gp->CreateGLMDir()) {
 
795
              if (f_setuponly)
 
796
                ;
 
797
              else if (f_run) {
 
798
                gp->CreateGLMJobs2();
 
799
                runseq(vbp,gp->seq,f_run);
 
800
              }
 
801
              else {
 
802
                gp->CreateGLMJobs();
 
803
                gp->WriteAndSubmitJobs();
 
804
              }
 
805
            }
 
806
          }
 
807
        }
 
808
        gp=&global;
 
809
      }
 
810
      else if (args[0]=="contrast") {
 
811
        gp->contrasts.push_back(args.Tail());
 
812
      }
 
813
      else if (args[0]=="endx" || args[0]=="noend") {
 
814
        gp=&global;
 
815
      }
 
816
      else if (args[0]=="include") {
 
817
        if (args.size() == 2)
 
818
          ParseFile(args[1]);
 
819
      }
 
820
      else {
 
821
        tmps.str("");
 
822
        tmps << "vbmakeglm: unrecognized keyword " << args[0];
 
823
        printErrorMsg(VB_ERROR,tmps.str());
 
824
      }
 
825
    }
 
826
  }
 
827
  infile.close();
 
828
  return 0; // no error!
 
829
}
 
830
 
 
831
void
 
832
GLMParams::CreateGLMJobs()
 
833
{
 
834
  VBJobSpec js;
 
835
  int rowstart,rowfinish,mergeflag;
 
836
  char tmp[STRINGLEN];
 
837
  int i,jobnum=0;
 
838
 
 
839
  seq.name=name;
 
840
 
 
841
  // set pieces heuristically if the user didn't
 
842
  if (pieces==0) {
 
843
    // FIXME is this any good?
 
844
    int cellcount=600000;
 
845
    pieces=(int)ceil((double)orderg*orderg/cellcount);
 
846
    if (pieces > orderg) pieces=orderg;
 
847
    if (pieces<1) pieces=1;
 
848
  }
 
849
  int rowsperjob=orderg/pieces;
 
850
  // are we going to need the merge jobs?
 
851
  if (rowsperjob < orderg)
 
852
    mergeflag=1;
 
853
  else
 
854
    mergeflag=0;
 
855
 
 
856
  // make the exofilt
 
857
  js.init();
 
858
  js.jobtype="vb_makefilter";
 
859
  js.arguments["outfile"]=stem+".ExoFilt";
 
860
  // sprintf(tmp,"lowflag -lf %d",lows);
 
861
  js.arguments["lowflag"]=(string)"-lf "+strnum(lows);
 
862
  if (middles.size())
 
863
    js.arguments["middleflag"]=(string)"-mf "+middles;
 
864
  else
 
865
    js.arguments["middleflag"]="";
 
866
  sprintf(tmp,"highflag -hf %d",highs);
 
867
  js.arguments["highflag"]=(string)"-hf "+strnum(highs);
 
868
  if (kernelname.size())
 
869
    js.arguments["kernelflag"]=(string)"-k "+kernelname;
 
870
  else
 
871
    js.arguments["kernelflag"]="";
 
872
  js.arguments["ndata"]=strnum(orderg);
 
873
  sprintf(tmp,"tr %f",TR);
 
874
  js.arguments["tr"]=strnum(TR);
 
875
  js.magnitude=0;    // set it to something!
 
876
  js.name="make exofilt";
 
877
  js.jnum=jobnum++;
 
878
  int n_exofilt=js.jnum;
 
879
  seq.addJob(js);
 
880
 
 
881
  // make the noisemodel
 
882
  js.init();
 
883
  js.jobtype="vb_makenoisemodel";
 
884
  js.arguments["outfile"]=stem+".IntrinCor";
 
885
  if (noisemodel.size())
 
886
    js.arguments["noisemodelflag"]="-n "+noisemodel;
 
887
  js.arguments["ndata"]=strnum(orderg);
 
888
  js.arguments["tr"]=strnum(TR);
 
889
  js.magnitude=0;    // set it to something!
 
890
  js.name="make noisemodel";
 
891
  js.jnum=jobnum++;
 
892
  int n_noisemodel=js.jnum;
 
893
  seq.addJob(js);
 
894
 
 
895
  // make the KG matrix
 
896
  js.init();
 
897
  js.jobtype="makematkg";
 
898
  js.arguments["stem"]=stem;
 
899
  js.magnitude=0;    // set it to something!
 
900
  js.name="GLM-KG";
 
901
  js.waitfor.insert(n_exofilt);
 
902
  js.jnum=jobnum++;
 
903
  int n_kg=js.jnum;
 
904
  seq.addJob(js);
 
905
 
 
906
  // make the F1 matrix
 
907
  js.init();
 
908
  js.jobtype="matpinv";
 
909
  js.arguments["in"]=stem+".KG";
 
910
  js.arguments["out"]=stem+".F1";
 
911
  js.magnitude=0;    // set it to something!
 
912
  js.name="GLM-F1";
 
913
  js.waitfor.insert(n_kg);
 
914
  js.jnum=jobnum++;
 
915
  int n_f1=js.jnum;
 
916
  seq.addJob(js);
 
917
 
 
918
  // make the R matrix = I-(KG)(F1)
 
919
  js.init();
 
920
  js.jobtype="matimxy";
 
921
  js.arguments["in1"]=stem+".KG";
 
922
  js.arguments["in2"]=stem+".F1";
 
923
  js.arguments["out"]=stem+".R";
 
924
  js.magnitude=0;    // set it to something!
 
925
  js.name="GLM-R";
 
926
  js.waitfor.insert(n_f1);
 
927
  js.jnum=jobnum++;
 
928
  int n_r=js.jnum;
 
929
  seq.addJob(js);
 
930
 
 
931
  // create a K matrix
 
932
  js.init();
 
933
  js.jobtype="makematk";
 
934
  js.arguments["stem"]=stem;
 
935
  js.magnitude=0;    // set it to something!
 
936
  js.name="GLM-K";
 
937
  js.waitfor.insert(n_exofilt);
 
938
  js.waitfor.insert(n_noisemodel);
 
939
  js.jnum=jobnum++;
 
940
  int n_k=js.jnum;
 
941
  seq.addJob(js);
 
942
 
 
943
  // create an empty V matrix
 
944
  js.init();
 
945
  js.jobtype="matzeros";
 
946
  js.arguments["name"]=stem+".V"; 
 
947
  js.arguments["cols"]=strnum(orderg);
 
948
  js.arguments["rows"]=strnum(orderg);
 
949
  js.magnitude=0;    // set it to something!
 
950
  js.name="GLM-Vcreate";
 
951
  js.jnum=jobnum++;
 
952
  int n_vcreate=js.jnum;
 
953
  seq.addJob(js);
 
954
 
 
955
  // build V matrix (V=KKt)
 
956
  set<int32> n_vpieces;
 
957
  rowstart=0;
 
958
  i=0;
 
959
  int n_v;
 
960
  do {
 
961
    rowfinish=rowstart+rowsperjob;
 
962
    if (rowfinish > orderg-1)
 
963
      rowfinish=orderg-1;
 
964
    js.init();
 
965
    js.jobtype="matxxt";
 
966
    js.arguments["in"]=stem+".K";
 
967
    js.arguments["out"]=stem+".V";
 
968
    js.arguments["col1"]=strnum(rowstart);
 
969
    js.arguments["col2"]=strnum(rowfinish);
 
970
    js.magnitude=0;    // set it to something!
 
971
    sprintf(tmp,"GLM-V%d",i++);
 
972
    js.name=tmp;
 
973
    js.waitfor.insert(n_k);
 
974
    js.waitfor.insert(n_vcreate);
 
975
    js.jnum=jobnum;
 
976
    n_vpieces.insert(js.jnum);
 
977
    n_v=js.jnum;  // in case it's just one piece
 
978
    seq.addJob(js);
 
979
    jobnum++;
 
980
    rowstart+=rowsperjob+1;
 
981
  } while (rowstart < orderg);
 
982
 
 
983
  // merge V
 
984
  if (mergeflag) {
 
985
    js.init();
 
986
    js.jobtype="matmerge";
 
987
    js.arguments["out"]=stem+".V";
 
988
    js.name="GLM-MergeV";
 
989
    js.waitfor=n_vpieces;
 
990
    js.jnum=jobnum++;
 
991
    n_v=js.jnum;
 
992
    seq.addJob(js);
 
993
  }
 
994
  
 
995
  // create F3
 
996
  js.init();
 
997
  js.jobtype="matf3";
 
998
  js.arguments["vmatrix"]=stem+".V";
 
999
  js.arguments["kgmatrix"]=stem+".KG";
 
1000
  js.arguments["outfile"]=stem+".F3";
 
1001
  js.magnitude=0;    // set it to something!
 
1002
  js.name="GLM-makeF3";
 
1003
  js.waitfor.insert(n_kg);
 
1004
  js.waitfor.insert(n_v);
 
1005
  js.jnum=jobnum++;
 
1006
  int n_f3=js.jnum;
 
1007
  seq.addJob(js);
 
1008
  
 
1009
  // create an empty RV matrix
 
1010
  js.init();
 
1011
  js.jobtype="matzeros";
 
1012
  js.arguments["name"]=stem+".RV";
 
1013
  js.arguments["cols"]=strnum(orderg);
 
1014
  js.arguments["rows"]=strnum(orderg);
 
1015
  js.magnitude=0;    // set it to something!
 
1016
  js.name="GLM-RVcreate";
 
1017
  js.jnum=jobnum++;
 
1018
  int n_rvcreate=js.jnum;
 
1019
  seq.addJob(js);
 
1020
 
 
1021
  // build RV matrix
 
1022
  set<int32> n_rvpieces;
 
1023
  rowstart=0;
 
1024
  int n_rv;
 
1025
  do {
 
1026
    rowfinish=rowstart+rowsperjob;
 
1027
    if (rowfinish > orderg-1)
 
1028
      rowfinish=orderg-1;
 
1029
    js.init();
 
1030
    js.jobtype="matxy";
 
1031
    js.arguments["in1"]=stem+".R";
 
1032
    js.arguments["in2"]=stem+".V";
 
1033
    js.arguments["out"]=stem+".RV";
 
1034
    js.arguments["col1"]=strnum(rowstart);
 
1035
    js.arguments["col2"]=strnum(rowfinish);
 
1036
    js.magnitude=0;    // set it to something!
 
1037
    sprintf(tmp,"GLM-RV%d",i);
 
1038
    js.name=tmp;
 
1039
    js.waitfor.insert(n_r);
 
1040
    js.waitfor.insert(n_v);
 
1041
    js.waitfor.insert(n_rvcreate);
 
1042
    js.jnum=jobnum;
 
1043
    n_rv=js.jnum;
 
1044
    seq.addJob(js);
 
1045
    n_rvpieces.insert(js.jnum);
 
1046
    n_rv=js.jnum;    // in case it's one piece
 
1047
    jobnum++;
 
1048
    rowstart+=rowsperjob+1;
 
1049
  } while (rowstart < orderg);
 
1050
 
 
1051
  // merge RV
 
1052
  if (mergeflag) {
 
1053
    js.init();
 
1054
    js.jobtype="matmerge";
 
1055
    js.arguments["out"]=stem+".RV";
 
1056
    js.name="GLM-MergeRV";
 
1057
    js.waitfor=n_rvpieces;
 
1058
    js.jnum=jobnum++;
 
1059
    n_rv=js.jnum;
 
1060
    seq.addJob(js);
 
1061
  }
 
1062
  
 
1063
  // create an empty RVRV matrix
 
1064
  js.init();
 
1065
  js.jobtype="matzeros";
 
1066
  js.arguments["name"]=stem+".RVRV";
 
1067
  js.arguments["cols"]=strnum(orderg);
 
1068
  js.arguments["rows"]=strnum(orderg);
 
1069
  js.magnitude=0;    // set it to something!
 
1070
  js.name="GLM-RVRVcreate";
 
1071
  js.jnum=jobnum++;
 
1072
  int n_rvrvcreate=js.jnum;
 
1073
  seq.addJob(js);
 
1074
 
 
1075
  // build RVRV
 
1076
  set<int32> n_rvrvpieces;
 
1077
  int n_rvrv;
 
1078
  rowstart=0;
 
1079
  do {
 
1080
    rowfinish=rowstart+rowsperjob;
 
1081
    if (rowfinish > orderg-1)
 
1082
      rowfinish=orderg-1;
 
1083
    js.init();
 
1084
    js.jobtype="matxx";
 
1085
    js.arguments["in"]=stem+".RV";
 
1086
    js.arguments["out"]=stem+".RVRV";
 
1087
    js.arguments["col1"]=strnum(rowstart);
 
1088
    js.arguments["col2"]=strnum(rowfinish);
 
1089
    js.magnitude=0;    // set it to something!
 
1090
    sprintf(tmp,"GLM-RVRV%d",i);
 
1091
    js.name=tmp;
 
1092
    js.waitfor.insert(n_rvrvcreate);
 
1093
    js.waitfor.insert(n_rv);
 
1094
    js.jnum=jobnum;
 
1095
    n_rvrvpieces.insert(js.jnum);
 
1096
    n_rvrv=js.jnum;
 
1097
    seq.addJob(js);
 
1098
    jobnum++;
 
1099
    rowstart+=rowsperjob+1;
 
1100
  } while (rowstart < orderg);
 
1101
 
 
1102
  // merge RVRV
 
1103
  if (mergeflag) {
 
1104
    js.init();
 
1105
    js.jobtype="matmerge";
 
1106
    js.arguments["out"]=stem+".RVRV";
 
1107
    js.name="GLM-MergeRVRV";
 
1108
    js.waitfor=n_rvrvpieces;
 
1109
    js.jnum=jobnum++;
 
1110
    n_rvrv=js.jnum;
 
1111
    seq.addJob(js);
 
1112
  }
 
1113
 
 
1114
  // traces
 
1115
  js.init();
 
1116
  js.jobtype="comptraces";
 
1117
  js.arguments["stem"]=stem;
 
1118
  js.magnitude=0;    // set it to something!
 
1119
  js.name="GLM-traces";
 
1120
  js.waitfor.insert(n_rvrv);
 
1121
  js.jnum=jobnum++;
 
1122
  int n_traces=js.jnum;
 
1123
  seq.addJob(js);
 
1124
  
 
1125
  // tes regression steps
 
1126
  set<int32> n_regstep;
 
1127
  for (i=0; i<pieces; i++) {
 
1128
    js.init();
 
1129
    js.jobtype="vbregress";
 
1130
    js.arguments["stem"]=stem;
 
1131
    js.arguments["steps"]=strnum(pieces);
 
1132
    js.arguments["index"]=strnum(i+1);
 
1133
    string flags;
 
1134
    if (this->meannorm)
 
1135
      flags+="-m ";
 
1136
    if (this->driftcorrect)
 
1137
      flags+="-d";
 
1138
    js.arguments["flags"]=flags;
 
1139
    js.magnitude=0;    // set it to something!
 
1140
    sprintf(tmp,"regress(%d/%d)",i+1,pieces);
 
1141
    js.name=tmp;
 
1142
    js.waitfor.insert(n_traces);
 
1143
    js.waitfor.insert(n_f1);
 
1144
    js.waitfor.insert(n_r);
 
1145
    js.waitfor.insert(n_exofilt);
 
1146
    js.jnum=jobnum;
 
1147
    n_regstep.insert(js.jnum);
 
1148
    seq.addJob(js);
 
1149
    jobnum++;
 
1150
  } // for i
 
1151
  
 
1152
  if (pieces>1) {
 
1153
    js.init();
 
1154
    js.jobtype="vbmerge4d";
 
1155
    js.arguments["stem"]=stem;
 
1156
    js.arguments["ext"]="prm";
 
1157
    js.magnitude=0;    // set it to something!
 
1158
    js.name="mergeparams";
 
1159
    js.waitfor=n_regstep;
 
1160
    js.jnum=jobnum;
 
1161
    seq.addJob(js);
 
1162
    jobnum++;
 
1163
 
 
1164
    js.init();
 
1165
    js.jobtype="vbmerge4d";
 
1166
    js.arguments["stem"]=stem;
 
1167
    js.arguments["ext"]="res";
 
1168
    js.magnitude=0;    // set it to something!
 
1169
    js.name="mergeparams";
 
1170
    js.waitfor=n_regstep;
 
1171
    js.jnum=jobnum;
 
1172
    seq.addJob(js);
 
1173
    jobnum++;
 
1174
 
 
1175
    // make these two jobs stand in for regstep in the waitfor chain
 
1176
    n_regstep.clear();
 
1177
    n_regstep.insert(jobnum-2);
 
1178
    n_regstep.insert(jobnum-1);
 
1179
  }
 
1180
  
 
1181
  // create smoothness estimate
 
1182
  js.init();
 
1183
  js.jobtype="vbse";
 
1184
  js.arguments["stem"]=stem;
 
1185
  js.magnitude=0;    // set it to something!
 
1186
  js.name="vbse";
 
1187
  js.waitfor=n_regstep;
 
1188
  js.jnum=jobnum;
 
1189
  int n_se=js.jnum;
 
1190
  seq.addJob(js);
 
1191
  jobnum++;
 
1192
 
 
1193
  // audit?
 
1194
  int n_audit;
 
1195
  if (auditflag) {
 
1196
    js.init();
 
1197
    js.jobtype="vb_auditglm";
 
1198
    js.arguments["glmdir"]=xdirname(stem);
 
1199
    js.magnitude=0;    // set it to something!
 
1200
    js.name="GLM-audit";
 
1201
    js.waitfor=n_regstep;
 
1202
    js.waitfor.insert(n_se);
 
1203
    js.jnum=jobnum++;
 
1204
    n_audit=js.jnum;
 
1205
    seq.addJob(js);
 
1206
  }
 
1207
  
 
1208
  // notify?
 
1209
  if (emailflag && email.size()) {
 
1210
    js.init();
 
1211
    js.jobtype="notify";
 
1212
    js.arguments["email"]=email;
 
1213
    js.arguments["msg"]="Your GLM has been solved.";
 
1214
    js.magnitude=0;    // set it to something!
 
1215
    js.name="Notify";
 
1216
    js.waitfor.insert(n_se);
 
1217
    if (auditflag)
 
1218
      js.waitfor.insert(n_audit);
 
1219
    js.waitfor.insert(n_f3);
 
1220
    js.jnum=jobnum++;
 
1221
    seq.addJob(js);
 
1222
  }
 
1223
}
 
1224
 
 
1225
void
 
1226
GLMParams::CreateGLMJobs2()
 
1227
{
 
1228
  VBJobSpec js;
 
1229
  int rowstart,rowfinish,mergeflag;
 
1230
  char tmp[STRINGLEN];
 
1231
  int i,jobnum=0;
 
1232
 
 
1233
  seq.name=name;
 
1234
 
 
1235
  // set pieces heuristically if the user didn't
 
1236
  if (pieces==0) {
 
1237
    // FIXME is this any good?
 
1238
    int cellcount=600000;
 
1239
    pieces=(int)ceil((double)orderg*orderg/cellcount);
 
1240
    if (pieces > orderg) pieces=orderg;
 
1241
    if (pieces<1) pieces=1;
 
1242
  }
 
1243
  int rowsperjob=orderg/pieces;
 
1244
  // are we going to need the merge jobs?
 
1245
  if (rowsperjob < orderg)
 
1246
    mergeflag=1;
 
1247
  else
 
1248
    mergeflag=0;
 
1249
 
 
1250
  // make the exofilt
 
1251
  js.init();
 
1252
  js.jobtype="shellcommand";
 
1253
  js.arguments["command"]=str(format("vbmakefilter -e %s.ExoFilt -lf %d -hf %d %s %s -t %d %f")
 
1254
                              %stem%lows%highs%(middles.size()?"-mf "+middles:"")
 
1255
                              %(kernelname.size()?"-k "+kernelname:"")
 
1256
                              %orderg%TR);
 
1257
  js.name="make exofilt";
 
1258
  js.jnum=jobnum++;
 
1259
  int n_exofilt=js.jnum;
 
1260
  seq.addJob(js);
 
1261
 
 
1262
  // make the noisemodel
 
1263
  js.init();
 
1264
  js.jobtype="shellcommand";
 
1265
  js.arguments["command"]=str(format("vbmakefilter -i %s.IntrinCor %s -t %d %f")
 
1266
                              %stem%(noisemodel.size()?"-n "+noisemodel:"")%orderg%TR);
 
1267
  js.name="make noisemodel";
 
1268
  js.jnum=jobnum++;
 
1269
  int n_noisemodel=js.jnum;
 
1270
  seq.addJob(js);
 
1271
 
 
1272
  // make the KG matrix
 
1273
  js.init();
 
1274
  js.jobtype="shellcommand";
 
1275
  js.arguments["command"]=str(format("makematkg -m %s")%stem);
 
1276
  js.name="GLM-KG";
 
1277
  js.waitfor.insert(n_exofilt);
 
1278
  js.jnum=jobnum++;
 
1279
  int n_kg=js.jnum;
 
1280
  seq.addJob(js);
 
1281
 
 
1282
  // make the F1 matrix
 
1283
  js.init();
 
1284
  js.jobtype="shellcommand";
 
1285
  js.arguments["command"]=str(format("vbmm2 -pinv %1%.KG %1%.F1")%stem);
 
1286
  js.name="GLM-F1";
 
1287
  js.waitfor.insert(n_kg);
 
1288
  js.jnum=jobnum++;
 
1289
  int n_f1=js.jnum;
 
1290
  seq.addJob(js);
 
1291
 
 
1292
  // make the R matrix = I-(KG)(F1)
 
1293
  js.init();
 
1294
  js.jobtype="shellcommand";
 
1295
  js.arguments["command"]=str(format("vbmm2 -imxy %1%.KG %1%.F1 %1%.R")%stem);
 
1296
  js.name="GLM-R";
 
1297
  js.waitfor.insert(n_f1);
 
1298
  js.jnum=jobnum++;
 
1299
  int n_r=js.jnum;
 
1300
  seq.addJob(js);
 
1301
 
 
1302
  // create a K matrix
 
1303
  js.init();
 
1304
  js.jobtype="shellcommand";
 
1305
  js.arguments["command"]=str(format("makematk -m %s")%stem);
 
1306
  js.name="GLM-K";
 
1307
  js.waitfor.insert(n_exofilt);
 
1308
  js.waitfor.insert(n_noisemodel);
 
1309
  js.jnum=jobnum++;
 
1310
  int n_k=js.jnum;
 
1311
  seq.addJob(js);
 
1312
 
 
1313
  // create an empty V matrix
 
1314
  js.init();
 
1315
  js.jobtype="shellcommand";
 
1316
  js.arguments["command"]=str(format("vbmm2 -zeros %s.V %d %d")%stem%orderg%orderg);
 
1317
  js.name="GLM-Vcreate";
 
1318
  js.jnum=jobnum++;
 
1319
  int n_vcreate=js.jnum;
 
1320
  seq.addJob(js);
 
1321
 
 
1322
  // build V matrix (V=KKt)
 
1323
  set<int32> n_vpieces;
 
1324
  rowstart=0;
 
1325
  i=0;
 
1326
  int n_v;
 
1327
  do {
 
1328
    rowfinish=rowstart+rowsperjob;
 
1329
    if (rowfinish > orderg-1)
 
1330
      rowfinish=orderg-1;
 
1331
    js.init();
 
1332
    js.jobtype="shellcommand";
 
1333
    js.arguments["command"]=str(format("vbmm2 -xyt %1%.K %1%.K %1%.V %2% %3%")%stem%rowstart%rowfinish);
 
1334
    sprintf(tmp,"GLM-V%d",i++);
 
1335
    js.name=tmp;
 
1336
    js.waitfor.insert(n_k);
 
1337
    js.waitfor.insert(n_vcreate);
 
1338
    js.jnum=jobnum;
 
1339
    n_vpieces.insert(js.jnum);
 
1340
    n_v=js.jnum;  // in case it's just one piece
 
1341
    seq.addJob(js);
 
1342
    jobnum++;
 
1343
    rowstart+=rowsperjob+1;
 
1344
  } while (rowstart < orderg);
 
1345
 
 
1346
  // merge V
 
1347
  if (mergeflag) {
 
1348
    js.init();
 
1349
    js.jobtype="shellcommand";
 
1350
    js.arguments["command"]=str(format("vbmm2 -assemblecols %s.V")%stem);
 
1351
    js.name="GLM-MergeV";
 
1352
    js.waitfor=n_vpieces;
 
1353
    js.jnum=jobnum++;
 
1354
    n_v=js.jnum;
 
1355
    seq.addJob(js);
 
1356
  }
 
1357
  
 
1358
  // create F3
 
1359
  js.init();
 
1360
  js.jobtype="shellcommand";
 
1361
  js.arguments["command"]=str(format("vbmm2 -f3 %1%.V %1%.KG %1%.F3")%stem);
 
1362
  js.name="GLM-makeF3";
 
1363
  js.waitfor.insert(n_kg);
 
1364
  js.waitfor.insert(n_v);
 
1365
  js.jnum=jobnum++;
 
1366
  int n_f3=js.jnum;
 
1367
  seq.addJob(js);
 
1368
  
 
1369
  // create an empty RV matrix
 
1370
  js.init();
 
1371
  js.jobtype="shellcommand";
 
1372
  js.arguments["command"]=str(format("vbmm2 -zeros %s.RV %d %d")%stem%orderg%orderg);
 
1373
  js.name="GLM-RVcreate";
 
1374
  js.jnum=jobnum++;
 
1375
  int n_rvcreate=js.jnum;
 
1376
  seq.addJob(js);
 
1377
 
 
1378
  // build RV matrix
 
1379
  set<int32> n_rvpieces;
 
1380
  rowstart=0;
 
1381
  int n_rv;
 
1382
  do {
 
1383
    rowfinish=rowstart+rowsperjob;
 
1384
    if (rowfinish > orderg-1)
 
1385
      rowfinish=orderg-1;
 
1386
    js.init();
 
1387
    js.jobtype="shellcommand";
 
1388
    js.arguments["command"]=str(format("vbmm2 -xyt %1%.R %1%.V %1%.RV %2% %3%")%stem%rowstart%rowfinish);
 
1389
    sprintf(tmp,"GLM-RV%d",i);
 
1390
    js.name=tmp;
 
1391
    js.waitfor.insert(n_r);
 
1392
    js.waitfor.insert(n_v);
 
1393
    js.waitfor.insert(n_rvcreate);
 
1394
    js.jnum=jobnum;
 
1395
    n_rv=js.jnum;
 
1396
    seq.addJob(js);
 
1397
    n_rvpieces.insert(js.jnum);
 
1398
    n_rv=js.jnum;    // in case it's one piece
 
1399
    jobnum++;
 
1400
    rowstart+=rowsperjob+1;
 
1401
  } while (rowstart < orderg);
 
1402
 
 
1403
  // merge RV
 
1404
  if (mergeflag) {
 
1405
    js.init();
 
1406
    js.jobtype="shellcommand";
 
1407
    js.arguments["command"]=str(format("vbmm2 -assemblecols %s.RV")%stem);
 
1408
    js.name="GLM-MergeRV";
 
1409
    js.waitfor=n_rvpieces;
 
1410
    js.jnum=jobnum++;
 
1411
    n_rv=js.jnum;
 
1412
    seq.addJob(js);
 
1413
  }
 
1414
  
 
1415
  // create an empty RVRV matrix
 
1416
  js.init();
 
1417
  js.jobtype="shellcommand";
 
1418
  js.arguments["command"]=str(format("vbmm2 -zeros %s.RVRV %d %d")%stem%orderg%orderg);
 
1419
  js.name="GLM-RVRVcreate";
 
1420
  js.jnum=jobnum++;
 
1421
  int n_rvrvcreate=js.jnum;
 
1422
  seq.addJob(js);
 
1423
 
 
1424
  // build RVRV
 
1425
  set<int32> n_rvrvpieces;
 
1426
  int n_rvrv;
 
1427
  rowstart=0;
 
1428
  do {
 
1429
    rowfinish=rowstart+rowsperjob;
 
1430
    if (rowfinish > orderg-1)
 
1431
      rowfinish=orderg-1;
 
1432
    js.init();
 
1433
    js.jobtype="shellcommand";
 
1434
    js.arguments["command"]=str(format("vbmm2 -xy %1%.RV %1%.RV %1%.RVRV %2% %3%")%stem%rowstart%rowfinish);
 
1435
    sprintf(tmp,"GLM-RVRV%d",i);
 
1436
    js.name=tmp;
 
1437
    js.waitfor.insert(n_rvrvcreate);
 
1438
    js.waitfor.insert(n_rv);
 
1439
    js.jnum=jobnum;
 
1440
    n_rvrvpieces.insert(js.jnum);
 
1441
    n_rvrv=js.jnum;
 
1442
    seq.addJob(js);
 
1443
    jobnum++;
 
1444
    rowstart+=rowsperjob+1;
 
1445
  } while (rowstart < orderg);
 
1446
 
 
1447
  // merge RVRV
 
1448
  if (mergeflag) {
 
1449
    js.init();
 
1450
    js.jobtype="shellcommand";
 
1451
    js.arguments["command"]=str(format("vbmm2 -assemblecols %s.RVRV")%stem);
 
1452
    js.name="GLM-MergeRVRV";
 
1453
    js.waitfor=n_rvrvpieces;
 
1454
    js.jnum=jobnum++;
 
1455
    n_rvrv=js.jnum;
 
1456
    seq.addJob(js);
 
1457
  }
 
1458
 
 
1459
  // traces
 
1460
  js.init();
 
1461
  js.jobtype="shellcommand";
 
1462
  js.arguments["command"]=str(format("comptraces -m %s")%stem);
 
1463
  js.name="GLM-traces";
 
1464
  js.waitfor.insert(n_rvrv);
 
1465
  js.jnum=jobnum++;
 
1466
  int n_traces=js.jnum;
 
1467
  seq.addJob(js);
 
1468
  
 
1469
  // tes regression steps
 
1470
  set<int32> n_regstep;
 
1471
  string flags;
 
1472
  if (this->meannorm)
 
1473
    flags+="-m ";
 
1474
  if (this->driftcorrect)
 
1475
    flags+="-d";
 
1476
  for (i=0; i<pieces; i++) {
 
1477
    js.init();
 
1478
    js.jobtype="shellcommand";
 
1479
    js.arguments["command"]=str(format("vbregress %s -p %d %d %s")
 
1480
                                %stem%(i+1)%pieces%flags);
 
1481
    sprintf(tmp,"regress(%d/%d)",i+1,pieces);
 
1482
    js.name=tmp;
 
1483
    js.waitfor.insert(n_traces);
 
1484
    js.waitfor.insert(n_f1);
 
1485
    js.waitfor.insert(n_r);
 
1486
    js.waitfor.insert(n_exofilt);
 
1487
    js.jnum=jobnum;
 
1488
    n_regstep.insert(js.jnum);
 
1489
    seq.addJob(js);
 
1490
    jobnum++;
 
1491
  } // for i
 
1492
  
 
1493
  if (pieces>1) {
 
1494
    js.init();
 
1495
    js.jobtype="shellcommand";
 
1496
    js.arguments["command"]=str(format("vbmerge4d %1%.prm_part_* -o %1%.prm")%stem);
 
1497
    js.name="mergeparams";
 
1498
    js.waitfor=n_regstep;
 
1499
    js.jnum=jobnum;
 
1500
    seq.addJob(js);
 
1501
    jobnum++;
 
1502
 
 
1503
    js.init();
 
1504
    js.jobtype="shellcommand";
 
1505
    js.arguments["command"]=str(format("rm %s.prm_part_*")%stem);
 
1506
    js.name="mergeparams";
 
1507
    js.waitfor.insert(jobnum-1);
 
1508
    js.jnum=jobnum;
 
1509
    seq.addJob(js);
 
1510
    jobnum++;
 
1511
 
 
1512
    js.init();
 
1513
    js.jobtype="shellcommand";
 
1514
    js.arguments["command"]=str(format("vbmerge4d %1%.res_part_* -o %1%.res")%stem);
 
1515
    js.name="mergeparams";
 
1516
    js.waitfor=n_regstep;
 
1517
    js.jnum=jobnum;
 
1518
    seq.addJob(js);
 
1519
    jobnum++;
 
1520
 
 
1521
    js.init();
 
1522
    js.jobtype="shellcommand";
 
1523
    js.arguments["command"]=str(format("rm %s.res_part_*")%stem);
 
1524
    js.name="mergeparams";
 
1525
    js.waitfor.insert(jobnum-1);
 
1526
    js.jnum=jobnum;
 
1527
    seq.addJob(js);
 
1528
    jobnum++;
 
1529
 
 
1530
    // make these four jobs stand in for regstep in the waitfor chain
 
1531
    n_regstep.clear();
 
1532
    n_regstep.insert(jobnum-4);
 
1533
    n_regstep.insert(jobnum-3);
 
1534
    n_regstep.insert(jobnum-2);
 
1535
    n_regstep.insert(jobnum-1);
 
1536
  }
 
1537
  
 
1538
  // create smoothness estimate
 
1539
  js.init();
 
1540
  js.jobtype="shellcommand";
 
1541
  js.arguments["command"]=str(format("vbse %1%.res %1%.se")%stem);
 
1542
  js.name="vbse";
 
1543
  js.waitfor=n_regstep;
 
1544
  js.jnum=jobnum;
 
1545
  int n_se=js.jnum;
 
1546
  seq.addJob(js);
 
1547
  jobnum++;
 
1548
 
 
1549
  // audit?
 
1550
  int n_audit;
 
1551
  if (auditflag) {
 
1552
    js.init();
 
1553
    js.jobtype="shellcommand";
 
1554
    js.arguments["command"]=str(format("glminfo -r %1% > %1%/audit.txt")%xdirname(stem));
 
1555
    js.name="GLM-audit";
 
1556
    js.waitfor=n_regstep;
 
1557
    js.waitfor.insert(n_se);
 
1558
    js.jnum=jobnum++;
 
1559
    n_audit=js.jnum;
 
1560
    seq.addJob(js);
 
1561
  }
 
1562
  
 
1563
  // notify?
 
1564
  if (emailflag && email.size()) {
 
1565
    js.init();
 
1566
    js.jobtype="notify";
 
1567
    js.arguments["email"]=email;
 
1568
    js.arguments["msg"]="Your GLM has been solved.";
 
1569
    js.magnitude=0;    // set it to something!
 
1570
    js.name="Notify";
 
1571
    js.waitfor.insert(n_se);
 
1572
    if (auditflag)
 
1573
      js.waitfor.insert(n_audit);
 
1574
    js.waitfor.insert(n_f3);
 
1575
    js.jnum=jobnum++;
 
1576
    seq.addJob(js);
 
1577
  }
 
1578
  for (SMI j=seq.specmap.begin(); j!=seq.specmap.end(); j++)
 
1579
    j->second.dirname=dirname;
 
1580
}
 
1581
 
 
1582
vector<string>
 
1583
GLMParams::CreateGLMScript()
 
1584
{
 
1585
  VBJobSpec js;
 
1586
  int rowstart,rowfinish,mergeflag;
 
1587
  int i;
 
1588
  string tmps;
 
1589
 
 
1590
  // set pieces heuristically if the user didn't
 
1591
  if (pieces==0) {
 
1592
    // FIXME is this any good?
 
1593
    int cellcount=600000;
 
1594
    pieces=(int)ceil((double)orderg*orderg/cellcount);
 
1595
    if (pieces > orderg) pieces=orderg;
 
1596
    if (pieces<1) pieces=1;
 
1597
  }
 
1598
  int rowsperjob=orderg/pieces;
 
1599
  // are we going to need the merge jobs?
 
1600
  if (rowsperjob < orderg)
 
1601
    mergeflag=1;
 
1602
  else
 
1603
    mergeflag=0;
 
1604
 
 
1605
  vector<string> commandlist;
 
1606
 
 
1607
  // make the exofilt
 
1608
  tmps=(format("vbmakefilter -e %s.ExoFilt -lf %d -hf %d %s %s -t %d %f")
 
1609
        %stem%lows%highs%(middles.size()?"-mf "+middles:"")
 
1610
        %(kernelname.size()?"-k "+kernelname:"")
 
1611
        %orderg%TR).str();
 
1612
  commandlist.push_back(tmps);
 
1613
  // make the noisemodel
 
1614
  tmps=(format("vbmakefilter -i %s.IntrinCor %s -t %d %f")
 
1615
        %stem%(noisemodel.size()?"-n "+noisemodel:"")%orderg%TR).str();
 
1616
  commandlist.push_back(tmps);
 
1617
  // make the KG matrix
 
1618
  tmps=(format("makematkg -m %s")%stem).str();
 
1619
  commandlist.push_back(tmps);
 
1620
  // make the F1 matrix
 
1621
  tmps=(format("vbmm2 -pinv %1%.KG %1%.F1")%stem).str();
 
1622
  commandlist.push_back(tmps);
 
1623
  // make the R matrix = I-(KG)(F1)
 
1624
  tmps=(format("vbmm2 -imxy %1%.KG %1%.F1 %1%.R")%stem).str();
 
1625
  commandlist.push_back(tmps);
 
1626
  // create a K matrix
 
1627
  tmps=(format("makematk -m %s")%stem).str();
 
1628
  commandlist.push_back(tmps);
 
1629
  // create an empty V matrix
 
1630
  tmps=(format("vbmm2 -zeros %s.V %d %d")%stem%orderg%orderg).str();
 
1631
  commandlist.push_back(tmps);
 
1632
 
 
1633
  // build V matrix (V=KKt)
 
1634
  set<int32> n_vpieces;
 
1635
  rowstart=0;
 
1636
  i=0;
 
1637
  do {
 
1638
    rowfinish=rowstart+rowsperjob;
 
1639
    if (rowfinish > orderg-1)
 
1640
      rowfinish=orderg-1;
 
1641
    js.init();
 
1642
    tmps=(format("vbmm2 -xyt %1%.K %1%.K %1%.V %2% %3%")%stem%rowstart%rowfinish).str();
 
1643
    commandlist.push_back(tmps);
 
1644
    rowstart+=rowsperjob+1;
 
1645
  } while (rowstart < orderg);
 
1646
 
 
1647
  // merge V
 
1648
  if (mergeflag) {
 
1649
    tmps=(format("vbmm2 -assemblecols %s.V")%stem).str();
 
1650
    commandlist.push_back(tmps);
 
1651
  }
 
1652
  
 
1653
  // create F3
 
1654
  tmps=(format("vbmm2 -f3 %1%.V %1%.KG %1%.F3")%stem).str();
 
1655
  commandlist.push_back(tmps);
 
1656
  
 
1657
  // create an empty RV matrix
 
1658
  tmps=(format("vbmm2 -zeros %s.RV %d %d")%stem%orderg%orderg).str();
 
1659
  commandlist.push_back(tmps);
 
1660
 
 
1661
  // build RV matrix
 
1662
  rowstart=0;
 
1663
  do {
 
1664
    rowfinish=rowstart+rowsperjob;
 
1665
    if (rowfinish > orderg-1)
 
1666
      rowfinish=orderg-1;
 
1667
    tmps=(format("vbmm2 -xyt %1%.R %1%.V %1%.RV %2% %3%")%stem%rowstart%rowfinish).str();
 
1668
    commandlist.push_back(tmps);
 
1669
    rowstart+=rowsperjob+1;
 
1670
  } while (rowstart < orderg);
 
1671
 
 
1672
  // merge RV
 
1673
  if (mergeflag) {
 
1674
    tmps=(format("vbmm2 -assemblecols %s.RV")%stem).str();
 
1675
    commandlist.push_back(tmps);
 
1676
  }
 
1677
  
 
1678
  // create an empty RVRV matrix
 
1679
  tmps=(format("vbmm2 -zeros %s.RVRV %d %d")%stem%orderg%orderg).str();
 
1680
  commandlist.push_back(tmps);
 
1681
 
 
1682
  // build RVRV
 
1683
  rowstart=0;
 
1684
  do {
 
1685
    rowfinish=rowstart+rowsperjob;
 
1686
    if (rowfinish > orderg-1)
 
1687
      rowfinish=orderg-1;
 
1688
    tmps=(format("vbmm2 -xy %1%.RV %1%.RV %1%.RVRV %2% %3%")%stem%rowstart%rowfinish).str();
 
1689
    commandlist.push_back(tmps);
 
1690
    rowstart+=rowsperjob+1;
 
1691
  } while (rowstart < orderg);
 
1692
 
 
1693
  // merge RVRV
 
1694
  if (mergeflag) {
 
1695
    tmps=(format("vbmm2 -assemblecols %s.RVRV")%stem).str();
 
1696
    commandlist.push_back(tmps);
 
1697
  }
 
1698
 
 
1699
  // traces
 
1700
  tmps=(format("comptraces -m %s")%stem).str();
 
1701
  commandlist.push_back(tmps);
 
1702
  
 
1703
  // tes regression steps
 
1704
  string flags;
 
1705
  if (this->meannorm)
 
1706
    flags+="-m ";
 
1707
  if (this->driftcorrect)
 
1708
    flags+="-d";
 
1709
  for (i=0; i<pieces; i++) {
 
1710
    tmps=(format("vbregress %s -p %d %d %s")
 
1711
          %stem%(i+1)%pieces%flags).str();
 
1712
    commandlist.push_back(tmps);
 
1713
  }
 
1714
  
 
1715
  if (pieces>1) {
 
1716
    tmps=str((format("vbmerge4d %1%.prm_part_* -o %1%.prm")%stem));
 
1717
    commandlist.push_back(tmps);
 
1718
    tmps=str((format("rm %s.prm_part_*")%stem));
 
1719
    commandlist.push_back(tmps);
 
1720
 
 
1721
    tmps=str((format("vbmerge4d %1%.res_part_* -o %1%.res")%stem));
 
1722
    commandlist.push_back(tmps);
 
1723
    tmps=str((format("rm %s.res_part_*")%stem));
 
1724
    commandlist.push_back(tmps);
 
1725
  }
 
1726
  
 
1727
  // create smoothness estimate
 
1728
  tmps=str((format("vbse %1%.res %1%.se")%stem));
 
1729
  commandlist.push_back(tmps);
 
1730
 
 
1731
  // audit?
 
1732
  if (auditflag) {
 
1733
    tmps=str((format("glminfo -r %1% > %1%/audit.txt")%xdirname(stem)));
 
1734
    commandlist.push_back(tmps);
 
1735
  }
 
1736
  return commandlist;
 
1737
}
 
1738
 
 
1739
int
 
1740
GLMParams::WriteAndSubmitJobs()
 
1741
{
 
1742
  // copy a few things to the sequence
 
1743
  seq.priority=pri;
 
1744
  seq.email=email;
 
1745
  for (SMI j=seq.specmap.begin(); j!=seq.specmap.end(); j++)
 
1746
    j->second.dirname=dirname;
 
1747
  return seq.Submit(vbp);
 
1748
}
 
1749
 
 
1750
int
 
1751
GLMParams::WriteGLMFile(string fname)
 
1752
{
 
1753
  FILE *fp;
 
1754
  if (fname.empty())
 
1755
    fname=stem+".glm";
 
1756
  fp=fopen(fname.c_str(),"w");
 
1757
  if (!fp) {
 
1758
    printf("[E] vbmakeglm: couldn't create glm file %s\n",fname.c_str());
 
1759
    return 103;
 
1760
  }
 
1761
  fprintf(fp,"lows %d\n",lows);
 
1762
  fprintf(fp,"highs %d\n",highs);
 
1763
  if (middles.size())
 
1764
    fprintf(fp,"middles %s\n",middles.c_str());
 
1765
  fprintf(fp,"orderg %d\n",orderg);
 
1766
  fprintf(fp,"pieces %d\n",pieces);
 
1767
  fprintf(fp,"kernel %s\n",kernelname.c_str());
 
1768
  fprintf(fp,"noisemodel %s\n",noisemodel.c_str());
 
1769
  if (rfxgflag)
 
1770
    fprintf(fp,"makerandfxg\n");
 
1771
  else
 
1772
    fprintf(fp,"gmatrix %s\n",gmatrix.c_str());
 
1773
  if (refname.size())
 
1774
    fprintf(fp,"refname %s\n",refname.c_str());
 
1775
  fprintf(fp,"pri %d\n",pri);
 
1776
  fprintf(fp,"audit %s\n",(auditflag ? "yes" : "no"));
 
1777
  fprintf(fp,"meannorm %s\n",(meannorm ? "yes" : "no"));
 
1778
  fprintf(fp,"driftcorrect %s\n",(driftcorrect ? "yes" : "no"));
 
1779
  fprintf(fp,"email %s\n",email.c_str());
 
1780
  fprintf(fp,"\n");
 
1781
  fprintf(fp,"glm %s\n",name.c_str());
 
1782
  fprintf(fp,"dirname %s\n",dirname.c_str());
 
1783
  for (int i=0; i<(int)scanlist.size(); i++)
 
1784
    fprintf(fp,"scan %s\n",scanlist[i].c_str());
 
1785
  fprintf(fp,"end\n");
 
1786
  fclose(fp);
 
1787
  return 0;  
 
1788
}
 
1789
 
 
1790
void
 
1791
GLMParams::FixRelativePaths()
 
1792
{
 
1793
  string path=xgetcwd()+"/";
 
1794
  dirname=xabsolutepath(dirname);
 
1795
  kernelname=xabsolutepath(kernelname);
 
1796
  noisemodel=xabsolutepath(noisemodel);
 
1797
  refname=xabsolutepath(refname);
 
1798
  gmatrix=xabsolutepath(gmatrix);
 
1799
  for (size_t i=0; i<scanlist.size(); i++)
 
1800
    scanlist[i]=xabsolutepath(scanlist[i]);
 
1801
}
 
1802
 
 
1803
void
 
1804
vbmakeglm_help()
 
1805
{
 
1806
  printf("\nVoxBo vbmakeglm (v%s)\n",vbversion.c_str());
 
1807
  printf("usage: vbmakeglm [flags] file [file ...]\n");
 
1808
  printf("flags:\n");
 
1809
  printf("    -i            ignore warnings, queue it anyway\n");
 
1810
  printf("    -v            don't queue, just validate\n");
 
1811
  printf("    -s            setup only (create filters)\n");
 
1812
  printf("    -tr <rate>    set data TR (ms)\n");
 
1813
  printf("    -kr <rate>    set HRF kernel sampling rate (ms)\n");
 
1814
  printf("    -h            i'm lazy, show me a sample .glm file\n");
 
1815
  printf("    -n <name>     only actually submit glms matching name\n");
 
1816
  printf("    -x <filename> i'm really lazy, make me a sample .glm file\n");
 
1817
  printf("    --run=<n>     don't queue, run now on n cores\n");
 
1818
  printf("notes:\n");
 
1819
  printf("  If --run is specified without =n, the default is the total number of\n");
 
1820
  printf("  available cores minus 1 (or 1, if that would be 0).\n");
 
1821
  printf("\n");
 
1822
}
 
1823
 
 
1824
void
 
1825
vbmakeglm_make_sample(const string &fname)
 
1826
{
 
1827
  printf("Creating sample glm file in %s...",fname.c_str());  fflush(stdout);
 
1828
  FILE *fp=fopen(fname.c_str(),"a");
 
1829
  if (fp) {
 
1830
    vbmakeglm_sample(fp);
 
1831
    fclose(fp);
 
1832
    printf("done.\n");
 
1833
  }
 
1834
  else {
 
1835
    printf("failed.\n");
 
1836
  }
 
1837
}
 
1838
 
 
1839
void
 
1840
vbmakeglm_show_sample()
 
1841
{
 
1842
  printf("Here's what a valid .glm file looks like, more or less:\n\n");
 
1843
 
 
1844
  vbmakeglm_sample(stdout);
 
1845
}
 
1846
 
 
1847
void
 
1848
vbmakeglm_sample(FILE *fp)
 
1849
{
 
1850
  fprintf(fp,"###############################################################\n");
 
1851
  fprintf(fp,"# Sample .glm file for VoxBo vbmakeglm.\n");
 
1852
  fprintf(fp,"# At the top, put stuff that's common to all GLMS in this file.\n");
 
1853
  fprintf(fp,"###############################################################\n");
 
1854
  fprintf(fp,"\n");
 
1855
 
 
1856
  fprintf(fp,"# how many low frequencies do you want to filter out?\n");
 
1857
  fprintf(fp,"lows 1\n\n");
 
1858
  fprintf(fp,"# how many high frequencies do you want to filter out?\n");
 
1859
  fprintf(fp,"highs 1\n\n");
 
1860
  fprintf(fp,"# which random frequencies fo you want to filter out? (separate by spaces)\n");
 
1861
  fprintf(fp,"# middles 15 16 21\n\n");
 
1862
  fprintf(fp,"# how many total data points per voxel (typically time points)?\n");
 
1863
  fprintf(fp,"orderg 160\n\n");
 
1864
  fprintf(fp,"# how many pieces do you want to break the matrix operations into?\n");
 
1865
  fprintf(fp,"pieces 10\n\n");
 
1866
  fprintf(fp,"# specify a kernel for exogenous smoothing by convolution\n");
 
1867
  fprintf(fp,"kernel /usr/local/VoxBo/elements/filters/Eigen1.ref\n\n");
 
1868
  fprintf(fp,"# specify a model of intrinsic noise\n");
 
1869
  fprintf(fp,"noisemodel /usr/local/VoxBo/elements/noisemodels/smooth_params.ref\n\n");
 
1870
  fprintf(fp,"# specify the location of your G matrix, or just \"makerandfxg\" for makerandfxg\n");
 
1871
  fprintf(fp,"gmatrix /data/Models/mytask.G\n");
 
1872
  fprintf(fp,"# makerandfxg\n\n");
 
1873
  fprintf(fp,"# the location of a reference function to copy into the GLM directory\n");
 
1874
  fprintf(fp,"refname /data/Models/motor.ref\n\n");
 
1875
  fprintf(fp,"# priority (1 for overnight, 2 for low, 3 for normal, 4 and 5 for various emergencies)\n");
 
1876
  fprintf(fp,"pri 3\n\n");
 
1877
  fprintf(fp,"# do some summary statistics on your GLM?\n");
 
1878
  fprintf(fp,"audit yes\n\n");
 
1879
  fprintf(fp,"# mean normalize?  usually say yes only if you have multiple runs of BOLD data\n");
 
1880
  fprintf(fp,"meannorm yes\n\n");
 
1881
  fprintf(fp,"# correct linear drift?  \n");
 
1882
  fprintf(fp,"driftcorrect yes\n\n");
 
1883
  fprintf(fp,"# your email address here\n");
 
1884
  fprintf(fp,"email nobody@nowhere.com\n\n");
 
1885
  fprintf(fp,"\n");
 
1886
 
 
1887
  fprintf(fp,"###############################################################\n");
 
1888
  fprintf(fp,"# Then create paragraphs of specific stuff for each GLM.\n");
 
1889
  fprintf(fp,"# (you can also override globals here)\n");
 
1890
  fprintf(fp,"###############################################################\n");
 
1891
  fprintf(fp,"\n");
 
1892
 
 
1893
  fprintf(fp,"glm larry-glm1\n");
 
1894
  fprintf(fp,"dirname /data/study/larry/glm1\n");
 
1895
  fprintf(fp,"scan /data/study/larry/larry01/larry01.tes\n");
 
1896
  fprintf(fp,"scan /data/study/larry/larry02/larry02.tes\n");
 
1897
  fprintf(fp,"end\n");
 
1898
  fprintf(fp,"\n");
 
1899
 
 
1900
  fprintf(fp,"glm moe-glm1\n");
 
1901
  fprintf(fp,"dirname /data/study/moe/glm1\n");
 
1902
  fprintf(fp,"scan /data/study/moe/moe01/moe01.tes\n");
 
1903
  fprintf(fp,"scan /data/study/moe/moe02/moe02.tes\n");
 
1904
  fprintf(fp,"end\n");
 
1905
  fprintf(fp,"\n");
 
1906
 
 
1907
  fprintf(fp,"glm shemp-glm1\n");
 
1908
  fprintf(fp,"dirname /data/study/shemp/glm1\n");
 
1909
  fprintf(fp,"scan /data/study/shemp/shemp01/shemp01.tes\n");
 
1910
  fprintf(fp,"scan /data/study/shemp/shemp02/shemp02.tes\n");
 
1911
  fprintf(fp,"end\n");
 
1912
  fprintf(fp,"\n");
 
1913
}