← Back to branch summary

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

« back to all changes in this revision

Viewing changes to crunch/vbmm2.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:
crunch/vbmm2.cpp
 
1
 
 
2
// vbmm2.cpp
 
3
// VoxBo matrix multiplication
 
4
// Copyright (c) 1998-2002 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
 
 
25
using namespace std;
 
26
 
 
27
 
 
28
 
 
29
// This code provides some simple matrix operations for VoxBo.
 
30
// Some of the functions are optimized to operate in parallel
 
31
// and/or on large matrices that won't fit in memory.  The
 
32
// implementation notes below are strictly for your entertainment.
 
33
 
 
34
#define MIN(a,b) ((a) < (b) ? (a) : (b))
 
35
#define MAX(a,b) ((a) < (b) ? (b) : (a))
 
36
 
 
37
#include <stdio.h>
 
38
#include <math.h>
 
39
#include <fcntl.h>
 
40
#include <unistd.h>
 
41
#include <string.h>
 
42
#include <signal.h>
 
43
#include <ctype.h>
 
44
#include <string>
 
45
#include "vbutil.h"
 
46
#include "vbio.h"
 
47
 
 
48
void do_print(tokenlist &args);
 
49
void do_printsub(tokenlist &args);
 
50
void do_compare(tokenlist &args);
 
51
int do_ident(tokenlist &args);
 
52
int do_zeros(tokenlist &args);
 
53
int do_random(tokenlist &args);
 
54
int do_invert(tokenlist &args);
 
55
int do_add(tokenlist &args);
 
56
int do_subtract(tokenlist &args);
 
57
int do_xyt(tokenlist &args);
 
58
int do_xy(tokenlist &args);
 
59
int do_imxy(tokenlist &args);
 
60
int do_f3(tokenlist &args);
 
61
int do_pinv(tokenlist &args);
 
62
int do_pca(tokenlist &args);
 
63
int do_assemblecols(tokenlist &args);
 
64
int do_assemblerows(tokenlist &args);
 
65
int do_xyz(tokenlist &args);
 
66
void vbmm_help();
 
67
 
 
68
int
 
69
main(int argc,char *argv[])
 
70
{
 
71
  tokenlist args;
 
72
  int err=0;
 
73
  string cmd;
 
74
 
 
75
  args.Transfer(argc-1,argv+1);
 
76
  if (args.size()<1) {
 
77
    vbmm_help();
 
78
    exit(0);
 
79
  }
 
80
  
 
81
  cmd=args[0];
 
82
  args.DeleteFirst();
 
83
 
 
84
  if (cmd=="-xy")
 
85
    err=do_xy(args);
 
86
  else if (cmd=="-imxy")
 
87
    err=do_imxy(args);
 
88
  else if (cmd=="-f3")
 
89
    err=do_f3(args);
 
90
  else if (cmd=="-pinv")
 
91
    err=do_pinv(args);
 
92
  else if (cmd=="-pca")
 
93
    err=do_pca(args);
 
94
  else if (cmd=="-xyz")
 
95
    err=do_xyz(args);
 
96
  else if (cmd=="-xyt")
 
97
    err=do_xyt(args);
 
98
  else if (cmd=="-add")
 
99
    err=do_add(args);
 
100
  else if (cmd=="-subtract")
 
101
    err=do_subtract(args);
 
102
  else if (cmd=="-invert")
 
103
    err=do_invert(args);
 
104
  else if (cmd=="-ident")
 
105
    err=do_ident(args);
 
106
  else if (cmd=="-zeros")
 
107
    err=do_zeros(args);
 
108
  else if (cmd=="-assemblecols")
 
109
    err=do_assemblecols(args);
 
110
  else if (cmd=="-assemblerows")
 
111
    err=do_assemblerows(args);
 
112
  else if (cmd=="-random")
 
113
    err=do_random(args);
 
114
  else if (cmd=="-print")
 
115
    do_print(args);
 
116
  else if (cmd=="-printsub")
 
117
    do_printsub(args);
 
118
  else if (cmd=="-compare")
 
119
    do_compare(args);
 
120
 
 
121
  exit(err);
 
122
}
 
123
 
 
124
// ident args: name size
 
125
 
 
126
int
 
127
do_ident(tokenlist &args)
 
128
{
 
129
  int r;
 
130
 
 
131
  if (args.size() < 2) {
 
132
    printf("[E] vbmm2: need a name and a size to create a matrix.\n");
 
133
    return 5;
 
134
  }
 
135
 
 
136
  r = strtol(args[1]);
 
137
  if (r < 0) {
 
138
    printf("[E] vbmm2: size for identity matrix must be > 0.\n");
 
139
    return 10;
 
140
  }
 
141
 
 
142
  printf("[I] vbmm2: creating identity matrix %s of size %d.\n",args[0].c_str(),r);
 
143
 
 
144
  VBMatrix target(r,r);
 
145
  target.ident();
 
146
  if (target.WriteMAT1(args[0]))
 
147
    printf("[E] vbmm2: identity matrix %s (%dx%d) not created.\n",args[0].c_str(),r,r);
 
148
  else
 
149
    printf("[I] vbmm2: identity matrix %s (%dx%d) created.\n",args[0].c_str(),r,r);
 
150
  return 0;
 
151
}
 
152
 
 
153
// zeros args: name size
 
154
 
 
155
int
 
156
do_zeros(tokenlist &args)
 
157
{
 
158
  if (args.size() < 2) {
 
159
    printf("[E] vbmm2: need a name and a size to create a matrix.\n");
 
160
    return 5;
 
161
  }
 
162
 
 
163
  int r,c;
 
164
  c = r = strtol(args[1]);
 
165
  if (args.size() > 2)         // needed only if it's not square
 
166
    c = strtol(args[2]);
 
167
  if (r < 0 || c < 0) {
 
168
    printf("[E] vbmm2: dimensions for zero matrix must be > 0.\n");
 
169
    return 10;
 
170
  }
 
171
 
 
172
  printf("vbmm2: creating %dx%d zero matrix %s.\n",r,c,args[0].c_str());
 
173
 
 
174
  VBMatrix target(c,r);
 
175
  target.zero();
 
176
  if (target.WriteMAT1(args[0]))
 
177
    printf("vbmm2: zero matrix %s (%dx%d) not created.\n",args[0].c_str(),r,c);
 
178
  else
 
179
    printf("vbmm2: zero matrix %s (%dx%d) created.\n",args[0].c_str(),r,c);
 
180
  return 0;
 
181
}
 
182
 
 
183
int
 
184
do_random(tokenlist &args)
 
185
{
 
186
  if (args.size() < 2) {
 
187
    printf("[E] vbmm2: need a name and a size to create a matrix.\n");
 
188
    return 5;
 
189
  }
 
190
  int r,c;
 
191
  c = r = strtol(args[1]);
 
192
  if (args.size() > 2)         // needed only if it's not square
 
193
    c = strtol(args[2]);
 
194
  if (r < 0 || c < 0) {
 
195
    printf("[E] vbmm2: dimensions for random matrix must be > 0.\n");
 
196
    return 10;
 
197
  }
 
198
 
 
199
  printf("vbmm2: creating random %dx%d matrix %s.\n",r,c,args[0].c_str());
 
200
 
 
201
  VBMatrix target(r,c);
 
202
  target.random();
 
203
  if (target.WriteMAT1(args[0]))
 
204
    printf("[E] vbmm2: error writing random %dx%d matrix %s\n",r,c,args[0].c_str());
 
205
  else
 
206
    printf("[I] vbmm2: random %dx%d matrix %s created\n",r,c,args[0].c_str());
 
207
  return 0;
 
208
}
 
209
 
 
210
// args: in1 in2 out
 
211
 
 
212
int
 
213
do_subtract(tokenlist &args)
 
214
{
 
215
  if (args.size() != 3) {
 
216
    printf("[E] vbmm2: usage: vbmm2 -subtract in1 in2 out\n");
 
217
    return 5;
 
218
  }
 
219
 
 
220
  VBMatrix mat1(args[0]);
 
221
  VBMatrix mat2(args[1]);
 
222
 
 
223
  if (mat1.m == 0 || mat1.n == 0) {
 
224
    printf("[E] vbmm2: first matrix was bad.\n");
 
225
    return 101;
 
226
  }
 
227
  if (mat2.m == 0 || mat2.n == 0) {
 
228
    printf("[E] vbmm2: second matrix was bad.\n");
 
229
    return 102;
 
230
  }
 
231
  if (mat1.m != mat2.m || mat1.n != mat2.n) {
 
232
    fprintf(stderr,"[E] vbmm2: matrix dimensions don't match.\n");
 
233
    return 103;
 
234
  }
 
235
  printf("[I] vbmm2: subtracting matrix %s from matrix %s.\n",args[1].c_str(),args[0].c_str());
 
236
  mat1-=mat2;
 
237
  mat1.WriteMAT1(args[2]);
 
238
  printf("[I] vbmm2: done.\n");
 
239
 
 
240
  return 0;
 
241
}
 
242
 
 
243
int
 
244
do_add(tokenlist &args)
 
245
{
 
246
  if (args.size() != 3) {
 
247
    printf("[E] vbmm2: usage: vbmms -add in1 in2 out\n");
 
248
    return 5;
 
249
  }
 
250
 
 
251
  VBMatrix mat1(args[0]);
 
252
  VBMatrix mat2(args[1]);
 
253
 
 
254
  if (mat1.m == 0 || mat1.n == 0) {
 
255
    printf("[E] vbmm2: first matrix was bad.\n");
 
256
    return 101;
 
257
  }
 
258
  if (mat2.m == 0 || mat2.n == 0) {
 
259
    printf("[E] vbmm2: second matrix was bad.\n");
 
260
    return 102;
 
261
  }
 
262
  if (mat1.m != mat2.m || mat1.n != mat2.n) {
 
263
    fprintf(stderr,"[E] vbmm2: matrix dimensions don't match.\n");
 
264
    return 103;
 
265
  }
 
266
 
 
267
  
 
268
  printf("[I] vbmm2: adding matrix %s and matrix %s.\n",args[0].c_str(),args[1].c_str());
 
269
  mat1+=mat2;
 
270
  mat1.WriteMAT1(args[2]);
 
271
  printf("[I] vbmm2: done.\n");
 
272
 
 
273
  return 0;
 
274
}
 
275
 
 
276
void
 
277
do_compare(tokenlist &args)
 
278
{
 
279
  if (args.size()!=2) {
 
280
    printf("[E] vbmm2: usage: vbmm2 -compare <mat1> <mat2>\n");
 
281
    return;
 
282
  }
 
283
  VBMatrix mat1(args[0]);
 
284
  VBMatrix mat2(args[1]);
 
285
  if (mat1.m!=mat2.m) {
 
286
    printf("[E] vbmm2: matrices have different row count\n");
 
287
    return;
 
288
  }
 
289
  if (mat1.n!=mat2.n) {
 
290
    printf("[E] vbmm2: matrices have different column count\n");
 
291
    return;
 
292
  }
 
293
  int diffs_all=0,diffs_diag=0,diffs_off=0;
 
294
  double totals_all=0.0,totals_diag=0.0,totals_off=0.0;
 
295
  double max_all=0.0,max_diag=0.0,max_off=0.0;
 
296
  double diff;
 
297
  
 
298
  for (int i=0; i<mat1.m; i++) {
 
299
    for (int j=0; j<mat1.n; j++) {
 
300
      diff=fabs(mat1(i,j)-mat2(i,j));
 
301
      if (diff==0.0) continue;
 
302
      diffs_all++;
 
303
      totals_all+=diff;
 
304
      if (diff>max_all) max_all=diff;
 
305
      if (i==j) {
 
306
        diffs_diag++;
 
307
        totals_diag+=diff;
 
308
        if (diff>max_diag) max_diag=diff;
 
309
      }
 
310
      else {
 
311
        diffs_off++;
 
312
        totals_off+=diff;
 
313
        if (diff>max_off) max_off=diff;
 
314
      }
 
315
    }
 
316
  }
 
317
  if (diffs_all)
 
318
    totals_all/=(double)diffs_all;
 
319
  if (diffs_diag)
 
320
    totals_diag/=(double)diffs_diag;
 
321
  if (diffs_off)
 
322
    totals_off/=(double)diffs_off;
 
323
  if (diffs_all==0)
 
324
    printf("[I] vbmm2: matrices are identical\n");
 
325
  else {
 
326
    printf("[I] vbmm2: %d total cells\n",mat1.m*mat1.n);
 
327
    printf("[I] vbmm2:    total: %d different cells, mean abs diff %g, max diff %g\n",diffs_all,totals_all,max_all);
 
328
    printf("[I] vbmm2: diagonal: %d different cells, mean abs diff %g, max diff %g\n",diffs_diag,totals_diag,max_diag);
 
329
    printf("[I] vbmm2: off-diag: %d different cells, mean abs diff %g, max diff %g\n",diffs_off,totals_off,max_off);
 
330
  }
 
331
}
 
332
 
 
333
// xy args: in1 in2 out [col1 col2]
 
334
 
 
335
int
 
336
do_xy(tokenlist &args)
 
337
{
 
338
  int c1,c2;
 
339
  
 
340
  if (args.size() != 3 && args.size() != 5) {
 
341
    printf("[E] vbmm2: usage: vbmm2 -xy in1 in2 out [c1 c2]\n");
 
342
    return 5;
 
343
  }
 
344
  VBMatrix mat1,mat2;
 
345
  mat1.ReadMAT1Header(args[0]);
 
346
  mat2.ReadMAT1Header(args[1]);
 
347
  if (mat1.m==0||mat1.n==0||mat2.m==0||mat2.n==0) {
 
348
    printf("[E] vbmm2: couldn't read matrix headers\n");
 
349
    return 100;
 
350
  }
 
351
  if (mat1.ReadMAT1(args[0])) {
 
352
    printf("[E] vbmm2: first matrix was bad.\n");
 
353
    return 100;
 
354
  }
 
355
 
 
356
  if (args.size()==5) {
 
357
    c1 = strtol(args[3]);
 
358
    c2 = strtol(args[4]);
 
359
  }
 
360
  else {
 
361
    c1=0;
 
362
    c2=mat2.cols-1;
 
363
  }
 
364
 
 
365
  // figure out the outfile name for this part, if we're not doing the whole thing
 
366
  string outname=args[2];
 
367
  char tmps[128];
 
368
  if (c1!=0 || c2!=mat2.cols-1) {
 
369
    sprintf(tmps,"_%08d_%08d",c1,c2);
 
370
    outname += tmps;
 
371
  }
 
372
 
 
373
  if (mat2.ReadMAT1(args[1],-1,-1,c1,c2)) {
 
374
    printf("[E] vbmm2: second matrix was bad.\n");
 
375
    return 101;
 
376
  }
 
377
 
 
378
  printf("vbmm2: multiplying matrix %s by matrix %s (cols %d to %d).\n",
 
379
         args[0].c_str(),args[1].c_str(),c1,c2);
 
380
  mat1*=mat2;
 
381
  if (mat1.WriteMAT1(outname))
 
382
    printf("[E] vbmm2: failed!\n");
 
383
  else
 
384
    printf("[I] vbmm2: done.\n");
 
385
 
 
386
  return 0;
 
387
}
 
388
 
 
389
int
 
390
do_xyt(tokenlist &args)
 
391
{
 
392
  int c1,c2;
 
393
  
 
394
  if (args.size() != 3 && args.size() != 5) {
 
395
    printf("[E] vbmm2: usage: vbmm2 -xy in1 in2 out [c1 c2]\n");
 
396
    return 5;
 
397
  }
 
398
  
 
399
  VBMatrix mat1,mat2;
 
400
  mat1.ReadMAT1Header(args[0]);
 
401
  mat2.ReadMAT1Header(args[1]);
 
402
  if (mat1.m==0||mat1.n==0||mat2.m==0||mat2.n==0) {
 
403
    printf("[E] vbmm2: couldn't read matrix headers\n");
 
404
    return 100;
 
405
  }
 
406
 
 
407
  if (mat1.ReadMAT1(args[0])) {
 
408
    printf("[E] vbmm2: first matrix was bad.\n");
 
409
    return 100;
 
410
  }
 
411
  if (args.size()==5) {
 
412
    c1 = strtol(args[3]);
 
413
    c2 = strtol(args[4]);
 
414
  }
 
415
  else {
 
416
    c1 = 0;
 
417
    c2 = mat2.rows-1;
 
418
  }
 
419
 
 
420
  // figure out the outfile name
 
421
  string outname=args[2];
 
422
  char tmps[128];
 
423
  if (c1!=0 || c2!=mat2.rows-1) {
 
424
    sprintf(tmps,"_%05d_%05d",c1,c2);
 
425
    outname += tmps;
 
426
  }
 
427
 
 
428
  if (mat2.ReadMAT1(args[1],c1,c2,-1,-1)) {
 
429
    printf("[E] vbmm2: second matrix was bad\n");
 
430
    return 101;
 
431
  }
 
432
 
 
433
  printf("[I] vbmm2: multiplying matrix %s by matrix %s (cols %d to %d).\n",
 
434
         args[0].c_str(),args[1].c_str(),c1,c2);
 
435
  mat2.transposed=1;
 
436
  mat1*=mat2;
 
437
  if (mat1.WriteMAT1(outname))
 
438
    printf("[E] vbmm2: failed\n");
 
439
  else
 
440
    printf("[I] vbmm2: done\n");
 
441
 
 
442
  return 0;
 
443
}
 
444
 
 
445
int
 
446
do_xyz(tokenlist &args)
 
447
{
 
448
  if (args.size() != 4) {
 
449
    printf("vbmm2: usage: vbmm2 -xyz in1 in2 in3 out\n");
 
450
    return 5;
 
451
  }
 
452
  
 
453
  VBMatrix mat1(args[0]);
 
454
  VBMatrix mat2(args[1]);
 
455
 
 
456
  if (!(mat1.dataValid() && mat2.dataValid())) {
 
457
    printf("[E] vbmm2: bad input matrix\n");
 
458
    return (100);
 
459
  }
 
460
  if (mat1.n != mat2.m) {
 
461
    printf("[E] vbmm2: bad matrix dimensions for xyz.\n");
 
462
    return (104);
 
463
  }
 
464
 
 
465
  printf("vbmm2: multiplying matrix %s by matrix %s.\n",args(0),args(1));
 
466
  mat1*=mat2;
 
467
  printf("vbmm2: multiplying matrix %s by matrix %s.\n",args(1),args(2));
 
468
  VBMatrix mat3(args[2]);
 
469
  if (mat1.n != mat3.m) {
 
470
    printf("[E] vbmm2: bad matrix dimensions for xyz.\n");
 
471
    return (104);
 
472
  }
 
473
  mat1*=mat3;
 
474
  printf("vbmm2: done.\n");
 
475
 
 
476
  mat1.WriteMAT1(args[3]);
 
477
 
 
478
  return 0;
 
479
}
 
480
 
 
481
int
 
482
do_assemblerows(tokenlist &)
 
483
{
 
484
  return (100);  // error!
 
485
  return (0);  // no error!
 
486
}
 
487
 
 
488
int
 
489
do_assemblecols(tokenlist &args)
 
490
{
 
491
  vglob vg(args[0]+"_*_*");
 
492
  if (vg.size() < 1) {
 
493
    printf("[E] vbmm2: no parts found for %s\n",args(0));
 
494
    return (101);
 
495
  }
 
496
  vector<VBMatrix *> mats;
 
497
  int rows=0,cols=0;
 
498
  // first read all the headers
 
499
  for (size_t i=0; i<vg.size(); i++) {
 
500
    VBMatrix tmp;
 
501
    tmp.ReadMAT1Header(vg[i]);
 
502
    if (!tmp.headerValid()) {
 
503
      printf("[E] vbmm2: invalid matrix in assemble list\n");
 
504
      return (102);
 
505
    }
 
506
    if (rows==0)
 
507
      rows=tmp.m;
 
508
    cols+=tmp.n;
 
509
    if (rows != tmp.m) {
 
510
      printf("[E] vbmm2: wrong-sized matrix %s in assemble list\n",vg[i].c_str());
 
511
      return (103);
 
512
    }
 
513
  }
 
514
  if (rows < 1 || cols < 1) {
 
515
    printf("[E] vbmm2: invalid size for assembled matrix: %d x %d\n",rows,cols);
 
516
    return (103);
 
517
  }
 
518
  VBMatrix newmat(rows,cols);
 
519
  
 
520
  int ind=0;
 
521
  for (size_t i=0; i<vg.size(); i++) {
 
522
    VBMatrix tmp(vg[i]);
 
523
    for (int j=0; j<tmp.n; j++) {
 
524
      VB_Vector vv=tmp.GetColumn(j);
 
525
      newmat.SetColumn(ind,vv);
 
526
      ind++;
 
527
    }
 
528
  }
 
529
  // try to unlink the actual files now that we're merged
 
530
  for (size_t i=0; i<vg.size(); i++)
 
531
    unlink(vg[i].c_str());
 
532
  newmat.WriteMAT1(args[0]);
 
533
  return (0);  // no error!
 
534
}
 
535
 
 
536
int
 
537
do_imxy(tokenlist &args)
 
538
{
 
539
  if (args.size() != 3 && args.size() != 5) {
 
540
    printf("[E] vbmm2: usage: vbmm -imxy in1 in2 out\n");
 
541
    return (100);
 
542
  }
 
543
  
 
544
  VBMatrix mat1(args[0]);
 
545
  VBMatrix mat2(args[1]);
 
546
 
 
547
  if (mat1.m <= 0 || mat1.n <= 0) {
 
548
    printf("[E] vbmm2: first matrix was bad.\n");
 
549
    return (101);
 
550
  }
 
551
  if (mat2.m == 0 || mat2.n == 0) {
 
552
    printf("[E] vbmm2: second matrix was bad.\n");
 
553
    return (102);
 
554
  }
 
555
  if (mat1.n != mat2.m) {
 
556
    printf("[E] vbmm2: incompatible matrix dimensions for I-XY.\n");
 
557
    return (103);
 
558
  }
 
559
 
 
560
  printf("vbmm2: I-XYing matrix %s by matrix %s.\n",
 
561
         args[0].c_str(),args[1].c_str());
 
562
  mat1*=mat2;
 
563
  for (int i=0; i<mat1.m; i++) {
 
564
    VB_Vector tmp=mat1.GetRow(i);
 
565
    for (int j=0; j<mat1.n; j++) {
 
566
      tmp[j]*=(double)-1.0;
 
567
      if (i==j)
 
568
        tmp[j]+=(double)1.0;
 
569
    }
 
570
    mat1.SetRow(i,tmp);
 
571
  }
 
572
  printf("vbmm2: done.\n");
 
573
  mat1.WriteMAT1(args[2]);
 
574
  return 0;
 
575
}
 
576
 
 
577
int
 
578
do_f3(tokenlist &args)
 
579
{
 
580
  if (args.size() != 3) {
 
581
    printf("[E] vbmm2: usage: vbmm -f3 v kg out\n");
 
582
    return (100);
 
583
  }
 
584
  printf("[I] vbmm2: creating F3 matrix (V*KG*invert(KGtKG))\n");
 
585
  printf("[I] vbmm2: V: %s\n",args(0));
 
586
  printf("[I] vbmm2: KG: %s\n",args(1));
 
587
  VBMatrix v(args[0]);
 
588
  VBMatrix kg(args[1]);
 
589
  VBMatrix kgt=kg;
 
590
  VBMatrix tmp;
 
591
 
 
592
  if (v.m==0 || kg.m==0 || kgt.m==0) {
 
593
    printf("[E] vbmm2: couldn't read matrices\n");
 
594
    return 100;
 
595
  }
 
596
  if (v.n != kg.m) {
 
597
    printf("[E] vbmm2: incompatible matrix dimensions\n");
 
598
    return 100;
 
599
  }
 
600
 
 
601
  kgt.transposed=1;
 
602
  kgt*=kg;
 
603
  kgt.transposed=0;
 
604
  invert(kgt,tmp);
 
605
  kgt.clear();  // free mem
 
606
  v*=kg;
 
607
  kg.clear();   // free mem
 
608
  v*=tmp;
 
609
  v.WriteMAT1(args[2]);
 
610
  printf("[I] vbmm2: wrote F3 matrix %s\n",args(2));
 
611
 
 
612
  return 0;
 
613
}
 
614
 
 
615
int
 
616
do_invert(tokenlist &args)
 
617
{
 
618
  if (args.size() != 2) {
 
619
    printf("[E] vbmm2: usage: vbmm -invert in out\n");
 
620
    return (100);
 
621
  }
 
622
  
 
623
  VBMatrix mat(args[0]);
 
624
 
 
625
  if (mat.m <= 0 || mat.n <= 0 || mat.m != mat.n) {
 
626
    printf("[E] vbmm2: input matrix for invert was bad.\n");
 
627
    return (101);
 
628
  }
 
629
  VBMatrix target(mat.m,mat.m);
 
630
  printf("vbmm2: inverting matrix %s.\n",args[0].c_str());
 
631
  invert(mat,target);
 
632
  printf("vbmm2: done.\n");
 
633
  target.WriteMAT1(args[1]);
 
634
  return 0;
 
635
}
 
636
 
 
637
// do_pinv() computes the pseudo-inverse, which is
 
638
// inverse(KGtKG) ## KGt
 
639
 
 
640
int
 
641
do_pinv(tokenlist &args)
 
642
{
 
643
  if (args.size() != 2) {
 
644
    printf("[E] vbmm2: usage: vbmm -pinv in out\n");
 
645
    return (100);
 
646
  }
 
647
  
 
648
  VBMatrix mat(args[0]);
 
649
 
 
650
  if (mat.m <= 0 || mat.n <= 0) {
 
651
    printf("[E] vbmm2: input matrix for pinv was bad.\n");
 
652
    return (101);
 
653
  }
 
654
 
 
655
  VBMatrix target(mat.n,mat.m);
 
656
  printf("vbmm2: pinv'ing matrix %s.\n",args[0].c_str());
 
657
  pinv(mat,target);
 
658
  printf("vbmm2: done.\n");
 
659
  target.WriteMAT1(args[1]);
 
660
  return 0;
 
661
}
 
662
 
 
663
// do_pca() calculates the principle components
 
664
 
 
665
int
 
666
do_pca(tokenlist &args)
 
667
{
 
668
  if (args.size() != 2) {
 
669
    printf("[E] vbmm2: usage: vbmm -pca in out\n");
 
670
    return (100);
 
671
  }
 
672
  
 
673
  VBMatrix mat(args[0]);
 
674
 
 
675
  if (mat.m <= 0 || mat.n <= 0) {
 
676
    printf("[E] vbmm2: input matrix for pinv was bad.\n");
 
677
    return (101);
 
678
  }
 
679
 
 
680
  VB_Vector lambdas;
 
681
  VBMatrix pcs;
 
682
  VBMatrix E;
 
683
  printf("vbmm2: pca'ing matrix %s.\n",args[0].c_str());
 
684
  pca(mat,lambdas,pcs,E);
 
685
  lambdas.print();
 
686
  printf("vbmm2: done.\n");
 
687
  pcs.WriteMAT1(args[1]);
 
688
  return 0;
 
689
}
 
690
 
 
691
void
 
692
do_print(tokenlist &args)
 
693
{
 
694
  if (args.size() < 1)
 
695
    return;
 
696
  for (int i=0; i<args.size(); i++) {
 
697
    VBMatrix tmp(args[i]);
 
698
    if (tmp.dataValid())
 
699
      tmp.print();
 
700
    else {
 
701
      printf("[E] vbmm2: couldn't open matrix %s to print.\n",args(i));
 
702
    }
 
703
  }
 
704
}
 
705
 
 
706
void
 
707
do_printsub(tokenlist &args)
 
708
{
 
709
  if (args.size() != 5)
 
710
    return;
 
711
  int r1=strtol(args[1]);
 
712
  int r2=strtol(args[2]);
 
713
  int c1=strtol(args[3]);
 
714
  int c2=strtol(args[4]);
 
715
  VBMatrix mat(args[0],r1,r2,c1,c2);
 
716
  if (!mat.rowdata) {
 
717
    printf("[E] vbmm2: couldn't read data\n");
 
718
    return;
 
719
  }
 
720
  printf("[I] vbmm: read rows %d-%d and cols %d-%d of %s\n",
 
721
         r1,r2,c1,c2,args(0));
 
722
  mat.print();
 
723
}
 
724
 
 
725
void
 
726
vbmm_help()
 
727
{
 
728
  printf("\nVoxBo vbmm2 (v%s)\n",vbversion.c_str());
 
729
  printf("| vbmm does various bits of matrix arithmetic with MAT1 files.\n");
 
730
  printf("| below, in and out refer to the input and output matrices.  c1 and c2\n");
 
731
  printf("| refer to the start and end columns to be produced (for parallelization).\n");
 
732
  printf("| usage:\n");
 
733
  printf("    vbmm -xyt <in1> <in2> <out> <c1> <c2 >      do part of XYt\n");
 
734
  printf("    vbmm -xy <in1> <in2> <out> <c1> <c2>        do part of XY\n");
 
735
  printf("    vbmm -imxy <in1> <in2> <out>                I-XY in core\n");
 
736
  printf("    vbmm -f3 <v> <kg> <out>                     V*KG*invert(KTtKG)\n");
 
737
  printf("    vbmm -xyz <in1> <in2> <in3> <out>           XYZ in core\n");
 
738
  printf("    vbmm -assemblecols <out>                    assemble out from available parts\n");
 
739
  // printf("    vbmm -assemblerows <out>                    assemble out from available parts\n");
 
740
  printf("    vbmm -add <in1> <in2> <out>                 add two matrices\n");
 
741
  printf("    vbmm -subtract <in1> <in2> <out>            subtract in2 from in1\n");
 
742
  printf("    vbmm -invert <in> <out>                     invert in\n"); 
 
743
  printf("    vbmm -pinv <in> <out>                       pseudo-inverse, in core\n");
 
744
  printf("    vbmm -pca <in> <out>                        calculte pca\n"); 
 
745
  printf("    vbmm -ident <name> <size>                   create an identity matrix\n");
 
746
  printf("    vbmm -zeros <name> <cols> [rows]            create a zero matrix\n");
 
747
  printf("    vbmm -random <name> <cols> [rows]           create a matrix of random numbers\n");
 
748
  printf("    vbmm -print <name>                          display a matrix\n");
 
749
  printf("    vbmm -printsub <name> <r1> <r2> <c1> <c2>   display part of a matrix\n");
 
750
  printf("    vbmm -compare <mat1> <mat2>                 compare two matrices\n");
 
751
}