ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/BioCocoa/BioCocoa/trunk/BCFoundation/BCUtils/BCDataMatrix.m
Revision: 350
Committed: Wed Aug 19 16:45:07 2015 UTC (4 years, 9 months ago) by schristley
File size: 37782 byte(s)
Log Message:
MPI timing
Line File contents
1 //
2 // BCDataMatrix.m
3 // BioCocoa
4 //
5 // Created by Scott Christley on 7/25/08.
6 // Copyright (c) 2003-2009 The BioCocoa Project.
7 // All rights reserved.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions
11 // are met:
12 // 1. Redistributions of source code must retain the above copyright
13 // notice, this list of conditions and the following disclaimer.
14 // 2. Redistributions in binary form must reproduce the above copyright
15 // notice, this list of conditions and the following disclaimer in the
16 // documentation and/or other materials provided with the distribution.
17 // 3. The name of the author may not be used to endorse or promote products
18 // derived from this software without specific prior written permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
21 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
22 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
23 // IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
24 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
25 // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
26 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
27 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
29 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 //
31
32 #import "BCDataMatrix.h"
33 #import "BCInternal.h"
34
35 #if 0
36 #ifdef GNUSTEP
37 // GNU Scientific Library
38 #include <gsl/gsl_matrix.h>
39 #include <gsl/gsl_linalg.h>
40 #else
41 // Accelerate framework on Mac
42 #import <Accelerate/Accelerate.h>
43 #endif
44 #endif
45
46 #if 0
47 #pragma mark == ENCODING STRINGS ==
48 #endif
49
50 char * const BCidEncode = @encode(id);
51 char * const BCintEncode = @encode(int);
52 char * const BCdoubleEncode = @encode(double);
53 char * const BCfloatEncode = @encode(float);
54 char * const BClongEncode = @encode(long);
55 char * const BCboolEncode = @encode(BOOL);
56
57 #if 0
58 #pragma mark == FORMAT STRINGS ==
59 #endif
60
61 NSString * const BCParseColumnNames = @"parseColumnNames";
62 NSString * const BCParseRowNames = @"parseRowNames";
63 NSString * const BCColumnNames = @"columnNames";
64 NSString * const BCRowNames = @"rowNames";
65 NSString * const BCSkipHeaderLines = @"skipHeaderLines";
66 NSString * const BCDataLayout = @"dataLayout";
67 NSString * const BCMatrixFormat = @"matrixFormat";
68 NSString * const BCListFormat = @"listFormat";
69 NSString * const BCSeparatorCharacterSet = @"separatorCharacterSet";
70
71 @implementation BCDataMatrix
72
73 + (BCDataMatrix *)emptyDataMatrixWithRows: (unsigned int)rows andColumns: (unsigned int)cols andEncode: (char *)anEncode
74 {
75 return [[[self alloc] initEmptyDataMatrixWithRows: rows andColumns: cols andEncode: anEncode] autorelease];
76 }
77
78 + (BCDataMatrix *)dataMatrixWithContentsOfFile: (NSString *)aFile andEncode: (char *)anEncode
79 {
80 return [[[self alloc] initWithContentsOfFile: aFile andEncode: anEncode] autorelease];
81 }
82
83 + (BCDataMatrix *)dataMatrixWithContentsOfFile: (NSString *)aFile andEncode: (char *)anEncode andFormat: (NSDictionary *)format
84 {
85 return [[[self alloc] initWithContentsOfFile: aFile andEncode: anEncode andFormat: format] autorelease];
86 }
87
88 + (BCDataMatrix *)dataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix
89 {
90 return [[[self alloc] initDataMatrixWithDataMatrix: aMatrix] autorelease];
91 }
92
93 + (BCDataMatrix *)dataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix andEncode: (char *)anEncode
94 {
95 return [[[self alloc] initDataMatrixWithDataMatrix: aMatrix andEncode: anEncode] autorelease];
96 }
97
98 + (BCDataMatrix *)dataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix andEncode: (char *)anEncode isColumnMajor: (BOOL)aFlag
99 {
100 return [[[self alloc] initDataMatrixWithDataMatrix: aMatrix andEncode: anEncode isColumnMajor: aFlag] autorelease];
101 }
102
103 - (BCDataMatrix *)initEmptyDataMatrixWithRows: (unsigned int)rows andColumns: (unsigned int)cols andEncode: (char *)anEncode
104 {
105 return [self initEmptyDataMatrixWithRows: rows andColumns: cols andEncode: anEncode isColumnMajor: NO];
106 }
107
108 - (BCDataMatrix *)initEmptyDataMatrixWithRows: (unsigned int)rows andColumns: (unsigned int)cols andEncode: (char *)anEncode isColumnMajor: (BOOL)aFlag
109 {
110 [super init];
111
112 numOfRows = rows;
113 numOfCols = cols;
114 encode = anEncode;
115 isColumnMajor = aFlag;
116
117 int i, j;
118 if (!strcmp(encode, BCidEncode)) {
119 dataMatrix = malloc(sizeof(id) * numOfRows * numOfCols);
120 id (*grid)[numOfRows][numOfCols] = dataMatrix;
121 for (i = 0;i < numOfRows; ++i)
122 for (j = 0;j < numOfCols; ++j)
123 (*grid)[i][j] = nil;
124 } else if (!strcmp(encode, BCintEncode)) {
125 dataMatrix = malloc(sizeof(int) * numOfRows * numOfCols);
126 int (*grid)[numOfRows][numOfCols] = dataMatrix;
127 for (i = 0;i < numOfRows; ++i)
128 for (j = 0;j < numOfCols; ++j)
129 (*grid)[i][j] = 0;
130 } else if (!strcmp(encode, BClongEncode)) {
131 dataMatrix = malloc(sizeof(long) * numOfRows * numOfCols);
132 long (*grid)[numOfRows][numOfCols] = dataMatrix;
133 for (i = 0;i < numOfRows; ++i)
134 for (j = 0;j < numOfCols; ++j)
135 (*grid)[i][j] = 0;
136 } else if (!strcmp(encode, BCfloatEncode)) {
137 dataMatrix = malloc(sizeof(float) * numOfRows * numOfCols);
138 float (*grid)[numOfRows][numOfCols] = dataMatrix;
139 for (i = 0;i < numOfRows; ++i)
140 for (j = 0;j < numOfCols; ++j)
141 (*grid)[i][j] = 0.0;
142 } else if (!strcmp(encode, BCdoubleEncode)) {
143 dataMatrix = malloc(sizeof(double) * numOfRows * numOfCols);
144 double (*grid)[numOfRows][numOfCols] = dataMatrix;
145 for (i = 0;i < numOfRows; ++i)
146 for (j = 0;j < numOfCols; ++j)
147 (*grid)[i][j] = 0.0;
148 } else if (!strcmp(encode, BCboolEncode)) {
149 dataMatrix = malloc(sizeof(BOOL) * numOfRows * numOfCols);
150 BOOL (*grid)[numOfRows][numOfCols] = dataMatrix;
151 for (i = 0;i < numOfRows; ++i)
152 for (j = 0;j < numOfCols; ++j)
153 (*grid)[i][j] = NO;
154 } else {
155 // throw exception or something
156 NSLog(@"ERROR: BCDataMatrix unknown encoding %s\n", anEncode);
157 return nil;
158 }
159
160 return self;
161 }
162
163 - (BCDataMatrix *)initWithContentsOfFile: (NSString *)aFile andEncode: (char *)anEncode
164 {
165 return [self initWithContentsOfFile: aFile andEncode: anEncode andFormat: nil];
166 }
167
168 - (BCDataMatrix *)initWithContentsOfFile: (NSString *)aFile andEncode: (char *)anEncode andFormat: (NSDictionary *)format
169 {
170 BOOL parseRowNames, parseColNames, matrixLayout;
171 NSCharacterSet *separator;
172
173 // default format is no row/column names and matrix layout
174 parseRowNames = NO;
175 parseColNames = NO;
176 matrixLayout = YES;
177 separator = [NSCharacterSet whitespaceCharacterSet];
178
179 if (format) {
180 NSString *s = [format objectForKey: BCDataLayout];
181 if ([s isEqualToString: BCListFormat]) matrixLayout = NO;
182 if (matrixLayout) {
183 // we can parse the row/columns in matrix layout
184 s = [format objectForKey: BCParseColumnNames];
185 if (s && [s boolValue]) parseColNames = YES;
186 s = [format objectForKey: BCParseRowNames];
187 if (s && [s boolValue]) parseRowNames = YES;
188 } else {
189 // row/columns names must be given to us for list layout
190 colNames = (NSArray *)[format objectForKey: BCColumnNames];
191 if (!colNames) {
192 NSLog(@"ERROR: List of column names required for list layout");
193 return nil;
194 }
195 colNames = [[NSArray alloc] initWithArray: colNames];
196 numOfCols = [colNames count];
197 rowNames = (NSArray *)[format objectForKey: BCRowNames];
198 if (!rowNames) {
199 NSLog(@"ERROR: List of row names required for list layout");
200 return nil;
201 }
202 rowNames = [[NSArray alloc] initWithArray: rowNames];
203 numOfRows = [rowNames count];
204 }
205 s = [format objectForKey: BCSeparatorCharacterSet];
206 if (s) separator = [NSCharacterSet characterSetWithCharactersInString: s];
207 }
208
209 int skipHeaderLines = [[format objectForKey: BCSkipHeaderLines] intValue];
210 int skippingLines = skipHeaderLines;
211
212 NSStringEncoding enc;
213 NSError *error;
214 NSString *contents = [NSString stringWithContentsOfFile: aFile usedEncoding: &enc error: &error];
215 if (!contents) {
216 NSLog(@"ERROR: Unable to read contents of file: %@", aFile);
217 return nil;
218 }
219
220 if (matrixLayout) {
221 NSUInteger start, end, next;
222 NSUInteger i, j;
223 NSAutoreleasePool *pool = [NSAutoreleasePool new];
224 NSRange range;
225 unsigned stringLength = [contents length];
226 BOOL parseFailure = NO;
227
228 // scan through it once to determine number of rows and columns
229 range.location = 0;
230 range.length = 1;
231 i = 0; j = 0;
232 do {
233 // one line at a time
234 [contents getLineStart:&start end:&next contentsEnd:&end forRange:range];
235 range.location = start;
236 range.length = end-start;
237 NSString *s = [contents substringWithRange: range];
238
239 // skip empty lines
240 if ([s length] == 0) {
241 range.location = next;
242 range.length = 1;
243 continue;
244 }
245
246 // skip any header lines
247 if (skippingLines) {
248 --skippingLines;
249 range.location = next;
250 range.length = 1;
251 continue;
252 }
253
254 NSArray *a = [s componentsSeparatedByCharactersInSet: separator];
255
256 if (i == 0) {
257 if (parseColNames) {
258 NSRange r = NSMakeRange(0, [a count]);
259 if (parseRowNames) {
260 ++r.location;
261 --r.length;
262 }
263 colNames = [[a subarrayWithRange: r] retain];
264 numOfCols = [colNames count];
265 } else {
266 if (parseRowNames) {
267 numOfCols = [a count] - 1;
268 } else {
269 numOfCols = [a count];
270 }
271 }
272 }
273 ++i;
274
275 range.location = next;
276 range.length = 1;
277 } while (next < stringLength);
278
279 int expectCols = numOfCols;
280 if (parseRowNames) {
281 rowNames = [NSMutableArray new];
282 expectCols = numOfCols + 1;
283 }
284 numOfRows = i;
285 if (parseColNames) --numOfRows;
286
287 [self initEmptyDataMatrixWithRows: numOfRows andColumns: numOfCols andEncode: anEncode];
288
289 skippingLines = skipHeaderLines;
290
291 // scan through again to parse the data
292 range.location = 0;
293 range.length = 1;
294 i = 0; j = 0;
295 do {
296 // one line at a time
297 [contents getLineStart:&start end:&next contentsEnd:&end forRange:range];
298 range.location = start;
299 range.length = end-start;
300 NSString *s = [contents substringWithRange: range];
301
302 // skip empty lines
303 if ([s length] == 0) {
304 range.location = next;
305 range.length = 1;
306 continue;
307 }
308
309 // skip any header lines
310 if (skippingLines) {
311 --skippingLines;
312 range.location = next;
313 range.length = 1;
314 continue;
315 }
316
317 if ((i == 0) && (parseColNames)) {
318 // skip the column names
319 parseColNames = NO;
320 } else {
321 NSArray *a = [s componentsSeparatedByCharactersInSet: separator];
322 if ([a count] != expectCols) {
323 NSLog(@"Invalid matrix format for data matrix, expected %d items, got %lu", expectCols, [a count]);
324 NSLog(@"Offending line:");
325 NSLog(@"%@", s);
326 parseFailure = YES;
327 break;
328 }
329
330 if (parseRowNames) [(NSMutableArray *)rowNames addObject: [a objectAtIndex: 0]];
331
332 if (!strcmp(encode, BCidEncode)) {
333 id (*grid)[numOfRows][numOfCols] = dataMatrix;
334 for (j = 0; j < numOfCols; ++j)
335 if (parseRowNames) (*grid)[i][j] = [a objectAtIndex: (j + 1)];
336 else (*grid)[i][j] = [a objectAtIndex: j];
337 } else if (!strcmp(encode, BCintEncode)) {
338 int (*grid)[numOfRows][numOfCols] = dataMatrix;
339 for (j = 0; j < numOfCols; ++j)
340 if (parseRowNames) (*grid)[i][j] = [[a objectAtIndex: (j + 1)] intValue];
341 else (*grid)[i][j] = [[a objectAtIndex: j] intValue];
342 } else if (!strcmp(encode, BClongEncode)) {
343 long (*grid)[numOfRows][numOfCols] = dataMatrix;
344 for (j = 0; j < numOfCols; ++j)
345 if (parseRowNames) (*grid)[i][j] = [[a objectAtIndex: (j + 1)] longValue];
346 else (*grid)[i][j] = [[a objectAtIndex: j] longValue];
347 } else if (!strcmp(encode, BCfloatEncode)) {
348 float (*grid)[numOfRows][numOfCols] = dataMatrix;
349 for (j = 0; j < numOfCols; ++j)
350 if (parseRowNames) (*grid)[i][j] = [[a objectAtIndex: (j + 1)] floatValue];
351 else (*grid)[i][j] = [[a objectAtIndex: j] floatValue];
352 } else if (!strcmp(encode, BCdoubleEncode)) {
353 double (*grid)[numOfRows][numOfCols] = dataMatrix;
354 for (j = 0; j < numOfCols; ++j)
355 if (parseRowNames) (*grid)[i][j] = [[a objectAtIndex: (j + 1)] doubleValue];
356 else (*grid)[i][j] = [[a objectAtIndex: j] doubleValue];
357 } else if (!strcmp(encode, BCboolEncode)) {
358 BOOL (*grid)[numOfRows][numOfCols] = dataMatrix;
359 for (j = 0; j < numOfCols; ++j)
360 if (parseRowNames) (*grid)[i][j] = [[a objectAtIndex: (j + 1)] boolValue];
361 else (*grid)[i][j] = [[a objectAtIndex: j] boolValue];
362 }
363
364 ++i;
365 }
366
367 range.location = next;
368 range.length = 1;
369 } while (next < stringLength);
370
371 [pool release];
372 if (parseFailure) {
373 [self dealloc];
374 return nil;
375 }
376
377 } else {
378 // The list layout has one line for each entry in the format
379 // rowName colName value
380
381 NSUInteger start, end, next;
382 NSUInteger i, j;
383 NSAutoreleasePool *pool = [NSAutoreleasePool new];
384 NSRange range;
385 unsigned stringLength = [contents length];
386 BOOL parseFailure = NO;
387
388 [self initEmptyDataMatrixWithRows: numOfRows andColumns: numOfCols andEncode: anEncode];
389
390 range.location = 0;
391 range.length = 1;
392 do {
393 // one line at a time
394 [contents getLineStart:&start end:&next contentsEnd:&end forRange:range];
395 range.location = start;
396 range.length = end-start;
397 NSString *s = [contents substringWithRange: range];
398
399 // skip empty lines
400 if ([s length] == 0) {
401 range.location = next;
402 range.length = 1;
403 continue;
404 }
405
406 // skip any header lines
407 if (skippingLines) {
408 --skippingLines;
409 range.location = next;
410 range.length = 1;
411 continue;
412 }
413
414
415 NSArray *a = [s componentsSeparatedByCharactersInSet: separator];
416 if ([a count] != 3) {
417 NSLog(@"Invalid list format for data matrix, expected 3 items, got %lu", [a count]);
418 NSLog(@"Offending line:");
419 NSLog(@"%@", s);
420 parseFailure = YES;
421 break;
422 }
423
424 i = [rowNames indexOfObject: [a objectAtIndex: 0]];
425 if (i == NSNotFound) {
426 NSLog(@"Found unknown row name in file: %@", [a objectAtIndex: 0]);
427 NSLog(@"Offending line:");
428 NSLog(@"%@", s);
429 parseFailure = YES;
430 break;
431 }
432
433 j = [rowNames indexOfObject: [a objectAtIndex: 1]];
434 if (j == NSNotFound) {
435 NSLog(@"Found unknown row name in file: %@", [a objectAtIndex: 1]);
436 NSLog(@"Offending line:");
437 NSLog(@"%@", s);
438 parseFailure = YES;
439 break;
440 }
441
442 if (!strcmp(encode, BCidEncode)) {
443 id (*grid)[numOfRows][numOfCols] = dataMatrix;
444 (*grid)[i][j] = [a objectAtIndex: 2];
445 } else if (!strcmp(encode, BCintEncode)) {
446 int (*grid)[numOfRows][numOfCols] = dataMatrix;
447 (*grid)[i][j] = [[a objectAtIndex: 2] intValue];
448 } else if (!strcmp(encode, BClongEncode)) {
449 long (*grid)[numOfRows][numOfCols] = dataMatrix;
450 (*grid)[i][j] = [[a objectAtIndex: 2] longValue];
451 } else if (!strcmp(encode, BCfloatEncode)) {
452 float (*grid)[numOfRows][numOfCols] = dataMatrix;
453 (*grid)[i][j] = [[a objectAtIndex: 2] floatValue];
454 } else if (!strcmp(encode, BCdoubleEncode)) {
455 double (*grid)[numOfRows][numOfCols] = dataMatrix;
456 (*grid)[i][j] = [[a objectAtIndex: 2] doubleValue];
457 } else if (!strcmp(encode, BCboolEncode)) {
458 BOOL (*grid)[numOfRows][numOfCols] = dataMatrix;
459 (*grid)[i][j] = [[a objectAtIndex: 2] boolValue];
460 }
461
462 range.location = next;
463 range.length = 1;
464 } while (next < stringLength);
465
466 [pool release];
467 if (parseFailure) {
468 [self dealloc];
469 return nil;
470 }
471 }
472
473 return self;
474 }
475
476 - (BCDataMatrix *)initDataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix
477 {
478 return [self initDataMatrixWithDataMatrix: aMatrix andEncode: [aMatrix matrixEncoding] isColumnMajor: [aMatrix isColumnMajor]];
479 }
480
481 - (BCDataMatrix *)initDataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix andEncode: (char *)anEncode
482 {
483 return [self initDataMatrixWithDataMatrix: aMatrix andEncode: anEncode isColumnMajor: [aMatrix isColumnMajor]];
484 }
485
486 - (BCDataMatrix *)initDataMatrixWithDataMatrix: (BCDataMatrix *)aMatrix andEncode: (char *)anEncode isColumnMajor: (BOOL)aFlag
487 {
488 // initialize with proper settings
489 [self initEmptyDataMatrixWithRows: [aMatrix numberOfRows] andColumns: [aMatrix numberOfColumns] andEncode: anEncode isColumnMajor: [aMatrix isColumnMajor]];
490
491 // copy data from provided data matrix
492 int i;
493 char *otherEncode = [aMatrix matrixEncoding];
494 if (!strcmp(encode, BCidEncode)) {
495 id *grid2 = dataMatrix;
496 if (!strcmp(otherEncode, BCidEncode)) {
497 id *grid1 = [aMatrix dataMatrix];
498 for (i = 0;i < numOfRows*numOfCols; ++i)
499 grid2[i] = grid1[i];
500 } else if (!strcmp(otherEncode, BCintEncode)) {
501 int *grid1 = [aMatrix dataMatrix];
502 for (i = 0;i < numOfRows*numOfCols; ++i)
503 [grid2[i] setIntValue: grid1[i]];
504 } else if (!strcmp(otherEncode, BClongEncode)) {
505 long *grid1 = [aMatrix dataMatrix];
506 for (i = 0;i < numOfRows*numOfCols; ++i)
507 [grid2[i] setLongValue: grid1[i]];
508 } else if (!strcmp(otherEncode, BCfloatEncode)) {
509 float *grid1 = [aMatrix dataMatrix];
510 for (i = 0;i < numOfRows*numOfCols; ++i)
511 [grid2[i] setFloatValue: grid1[i]];
512 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
513 double *grid1 = [aMatrix dataMatrix];
514 for (i = 0;i < numOfRows*numOfCols; ++i)
515 [grid2[i] setDoubleValue: grid1[i]];
516 } else if (!strcmp(otherEncode, BCboolEncode)) {
517 BOOL *grid1 = [aMatrix dataMatrix];
518 for (i = 0;i < numOfRows*numOfCols; ++i)
519 [grid2[i] setBoolValue: grid1[i]];
520 }
521 } else if (!strcmp(encode, BCintEncode)) {
522 int *grid2 = dataMatrix;
523 if (!strcmp(otherEncode, BCidEncode)) {
524 id *grid1 = [aMatrix dataMatrix];
525 for (i = 0;i < numOfRows*numOfCols; ++i)
526 grid2[i] = [grid1[i] intValue];
527 } else if (!strcmp(otherEncode, BCintEncode)) {
528 int *grid1 = [aMatrix dataMatrix];
529 for (i = 0;i < numOfRows*numOfCols; ++i)
530 grid2[i] = grid1[i];
531 } else if (!strcmp(otherEncode, BClongEncode)) {
532 long *grid1 = [aMatrix dataMatrix];
533 for (i = 0;i < numOfRows*numOfCols; ++i)
534 grid2[i] = grid1[i];
535 } else if (!strcmp(otherEncode, BCfloatEncode)) {
536 float *grid1 = [aMatrix dataMatrix];
537 for (i = 0;i < numOfRows*numOfCols; ++i)
538 grid2[i] = grid1[i];
539 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
540 double *grid1 = [aMatrix dataMatrix];
541 for (i = 0;i < numOfRows*numOfCols; ++i)
542 grid2[i] = grid1[i];
543 } else if (!strcmp(otherEncode, BCboolEncode)) {
544 BOOL *grid1 = [aMatrix dataMatrix];
545 for (i = 0;i < numOfRows*numOfCols; ++i)
546 grid2[i] = grid1[i];
547 }
548 } else if (!strcmp(encode, BClongEncode)) {
549 long *grid2 = dataMatrix;
550 if (!strcmp(otherEncode, BCidEncode)) {
551 id *grid1 = [aMatrix dataMatrix];
552 for (i = 0;i < numOfRows*numOfCols; ++i)
553 grid2[i] = [grid1[i] longValue];
554 } else if (!strcmp(otherEncode, BCintEncode)) {
555 int *grid1 = [aMatrix dataMatrix];
556 for (i = 0;i < numOfRows*numOfCols; ++i)
557 grid2[i] = grid1[i];
558 } else if (!strcmp(otherEncode, BClongEncode)) {
559 long *grid1 = [aMatrix dataMatrix];
560 for (i = 0;i < numOfRows*numOfCols; ++i)
561 grid2[i] = grid1[i];
562 } else if (!strcmp(otherEncode, BCfloatEncode)) {
563 float *grid1 = [aMatrix dataMatrix];
564 for (i = 0;i < numOfRows*numOfCols; ++i)
565 grid2[i] = grid1[i];
566 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
567 double *grid1 = [aMatrix dataMatrix];
568 for (i = 0;i < numOfRows*numOfCols; ++i)
569 grid2[i] = grid1[i];
570 } else if (!strcmp(otherEncode, BCboolEncode)) {
571 BOOL *grid1 = [aMatrix dataMatrix];
572 for (i = 0;i < numOfRows*numOfCols; ++i)
573 grid2[i] = grid1[i];
574 }
575 } else if (!strcmp(encode, BCfloatEncode)) {
576 float *grid2 = dataMatrix;
577 if (!strcmp(otherEncode, BCidEncode)) {
578 id *grid1 = [aMatrix dataMatrix];
579 for (i = 0;i < numOfRows*numOfCols; ++i)
580 grid2[i] = [grid1[i] floatValue];
581 } else if (!strcmp(otherEncode, BCintEncode)) {
582 int *grid1 = [aMatrix dataMatrix];
583 for (i = 0;i < numOfRows*numOfCols; ++i)
584 grid2[i] = grid1[i];
585 } else if (!strcmp(otherEncode, BClongEncode)) {
586 long *grid1 = [aMatrix dataMatrix];
587 for (i = 0;i < numOfRows*numOfCols; ++i)
588 grid2[i] = grid1[i];
589 } else if (!strcmp(otherEncode, BCfloatEncode)) {
590 float *grid1 = [aMatrix dataMatrix];
591 for (i = 0;i < numOfRows*numOfCols; ++i)
592 grid2[i] = grid1[i];
593 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
594 double *grid1 = [aMatrix dataMatrix];
595 for (i = 0;i < numOfRows*numOfCols; ++i)
596 grid2[i] = grid1[i];
597 } else if (!strcmp(otherEncode, BCboolEncode)) {
598 BOOL *grid1 = [aMatrix dataMatrix];
599 for (i = 0;i < numOfRows*numOfCols; ++i)
600 grid2[i] = grid1[i];
601 }
602 } else if (!strcmp(encode, BCdoubleEncode)) {
603 double *grid2 = dataMatrix;
604 if (!strcmp(otherEncode, BCidEncode)) {
605 id *grid1 = [aMatrix dataMatrix];
606 for (i = 0;i < numOfRows*numOfCols; ++i)
607 grid2[i] = [grid1[i] doubleValue];
608 } else if (!strcmp(otherEncode, BCintEncode)) {
609 int *grid1 = [aMatrix dataMatrix];
610 for (i = 0;i < numOfRows*numOfCols; ++i)
611 grid2[i] = grid1[i];
612 } else if (!strcmp(otherEncode, BClongEncode)) {
613 long *grid1 = [aMatrix dataMatrix];
614 for (i = 0;i < numOfRows*numOfCols; ++i)
615 grid2[i] = grid1[i];
616 } else if (!strcmp(otherEncode, BCfloatEncode)) {
617 float *grid1 = [aMatrix dataMatrix];
618 for (i = 0;i < numOfRows*numOfCols; ++i)
619 grid2[i] = grid1[i];
620 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
621 double *grid1 = [aMatrix dataMatrix];
622 for (i = 0;i < numOfRows*numOfCols; ++i)
623 grid2[i] = grid1[i];
624 } else if (!strcmp(otherEncode, BCboolEncode)) {
625 BOOL *grid1 = [aMatrix dataMatrix];
626 for (i = 0;i < numOfRows*numOfCols; ++i)
627 grid2[i] = grid1[i];
628 }
629 } else if (!strcmp(encode, BCboolEncode)) {
630 BOOL *grid2 = dataMatrix;
631 if (!strcmp(otherEncode, BCidEncode)) {
632 id *grid1 = [aMatrix dataMatrix];
633 for (i = 0;i < numOfRows*numOfCols; ++i)
634 grid2[i] = [grid1[i] boolValue];
635 } else if (!strcmp(otherEncode, BCintEncode)) {
636 int *grid1 = [aMatrix dataMatrix];
637 for (i = 0;i < numOfRows*numOfCols; ++i)
638 grid2[i] = grid1[i];
639 } else if (!strcmp(otherEncode, BClongEncode)) {
640 long *grid1 = [aMatrix dataMatrix];
641 for (i = 0;i < numOfRows*numOfCols; ++i)
642 grid2[i] = grid1[i];
643 } else if (!strcmp(otherEncode, BCfloatEncode)) {
644 float *grid1 = [aMatrix dataMatrix];
645 for (i = 0;i < numOfRows*numOfCols; ++i)
646 grid2[i] = grid1[i];
647 } else if (!strcmp(otherEncode, BCdoubleEncode)) {
648 double *grid1 = [aMatrix dataMatrix];
649 for (i = 0;i < numOfRows*numOfCols; ++i)
650 grid2[i] = grid1[i];
651 } else if (!strcmp(otherEncode, BCboolEncode)) {
652 BOOL *grid1 = [aMatrix dataMatrix];
653 for (i = 0;i < numOfRows*numOfCols; ++i)
654 grid2[i] = grid1[i];
655 }
656 }
657
658 if ([aMatrix columnNames]) [self setColumnNames: [NSArray arrayWithArray: [aMatrix columnNames]]];
659 if ([aMatrix rowNames]) [self setRowNames: [NSArray arrayWithArray: [aMatrix rowNames]]];
660
661 // set the data format
662 [self setColumnMajor: aFlag];
663
664 return self;
665 }
666
667 - (void)dealloc
668 {
669 if (rowNames) [rowNames release];
670 if (colNames) [colNames release];
671 if (dataMatrix) free(dataMatrix);
672 [super dealloc];
673 }
674
675 - (unsigned int)numberOfRows { return numOfRows; }
676 - (unsigned int)numberOfColumns { return numOfCols; }
677 - (void *)dataMatrix { return dataMatrix; }
678 - (char *)matrixEncoding { return encode; }
679
680 - (BOOL)isColumnMajor { return isColumnMajor; }
681 - (void)setColumnMajor: (BOOL)aFlag
682 {
683 int i, j;
684
685 if (aFlag == isColumnMajor) return;
686
687 if (!strcmp(encode, BCidEncode)) {
688 if (isColumnMajor) {
689 id (*tmp)[numOfRows][numOfCols] = malloc(sizeof(id) * numOfRows * numOfCols);
690 id (*grid)[numOfCols][numOfRows] = dataMatrix;
691 for (i = 0;i < numOfRows; ++i)
692 for (j = 0;j < numOfCols; ++j)
693 (*tmp)[i][j] = (*grid)[j][i];
694 memcpy(grid, tmp, sizeof(id) * numOfRows * numOfCols);
695 free(tmp);
696 } else {
697 id (*tmp)[numOfCols][numOfRows] = malloc(sizeof(id) * numOfRows * numOfCols);
698 id (*grid)[numOfRows][numOfCols] = dataMatrix;
699 for (i = 0;i < numOfRows; ++i)
700 for (j = 0;j < numOfCols; ++j)
701 (*tmp)[j][i] = (*grid)[i][j];
702 memcpy(grid, tmp, sizeof(id) * numOfRows * numOfCols);
703 free(tmp);
704 }
705 } else if (!strcmp(encode, BCintEncode)) {
706 if (isColumnMajor) {
707 int (*tmp)[numOfRows][numOfCols] = malloc(sizeof(int) * numOfRows * numOfCols);
708 int (*grid)[numOfCols][numOfRows] = dataMatrix;
709 for (i = 0;i < numOfRows; ++i)
710 for (j = 0;j < numOfCols; ++j)
711 (*tmp)[i][j] = (*grid)[j][i];
712 memcpy(grid, tmp, sizeof(int) * numOfRows * numOfCols);
713 free(tmp);
714 } else {
715 int (*tmp)[numOfCols][numOfRows] = malloc(sizeof(int) * numOfRows * numOfCols);
716 int (*grid)[numOfRows][numOfCols] = dataMatrix;
717 for (i = 0;i < numOfRows; ++i)
718 for (j = 0;j < numOfCols; ++j)
719 (*tmp)[j][i] = (*grid)[i][j];
720 memcpy(grid, tmp, sizeof(int) * numOfRows * numOfCols);
721 free(tmp);
722 }
723 } else if (!strcmp(encode, BClongEncode)) {
724 if (isColumnMajor) {
725 long (*tmp)[numOfRows][numOfCols] = malloc(sizeof(long) * numOfRows * numOfCols);
726 long (*grid)[numOfCols][numOfRows] = dataMatrix;
727 for (i = 0;i < numOfRows; ++i)
728 for (j = 0;j < numOfCols; ++j)
729 (*tmp)[i][j] = (*grid)[j][i];
730 memcpy(grid, tmp, sizeof(long) * numOfRows * numOfCols);
731 free(tmp);
732 } else {
733 long (*tmp)[numOfCols][numOfRows] = malloc(sizeof(long) * numOfRows * numOfCols);
734 long (*grid)[numOfRows][numOfCols] = dataMatrix;
735 for (i = 0;i < numOfRows; ++i)
736 for (j = 0;j < numOfCols; ++j)
737 (*tmp)[j][i] = (*grid)[i][j];
738 memcpy(grid, tmp, sizeof(long) * numOfRows * numOfCols);
739 free(tmp);
740 }
741 } else if (!strcmp(encode, BCfloatEncode)) {
742 if (isColumnMajor) {
743 float (*tmp)[numOfRows][numOfCols] = malloc(sizeof(float) * numOfRows * numOfCols);
744 float (*grid)[numOfCols][numOfRows] = dataMatrix;
745 for (i = 0;i < numOfRows; ++i)
746 for (j = 0;j < numOfCols; ++j)
747 (*tmp)[i][j] = (*grid)[j][i];
748 memcpy(grid, tmp, sizeof(float) * numOfRows * numOfCols);
749 free(tmp);
750 } else {
751 float (*tmp)[numOfCols][numOfRows] = malloc(sizeof(float) * numOfRows * numOfCols);
752 float (*grid)[numOfRows][numOfCols] = dataMatrix;
753 for (i = 0;i < numOfRows; ++i)
754 for (j = 0;j < numOfCols; ++j)
755 (*tmp)[j][i] = (*grid)[i][j];
756 memcpy(grid, tmp, sizeof(float) * numOfRows * numOfCols);
757 free(tmp);
758 }
759 } else if (!strcmp(encode, BCdoubleEncode)) {
760 if (isColumnMajor) {
761 double (*tmp)[numOfRows][numOfCols] = malloc(sizeof(double) * numOfRows * numOfCols);
762 double (*grid)[numOfCols][numOfRows] = dataMatrix;
763 for (i = 0;i < numOfRows; ++i)
764 for (j = 0;j < numOfCols; ++j)
765 (*tmp)[i][j] = (*grid)[j][i];
766 memcpy(grid, tmp, sizeof(double) * numOfRows * numOfCols);
767 free(tmp);
768 } else {
769 double (*tmp)[numOfCols][numOfRows] = malloc(sizeof(double) * numOfRows * numOfCols);
770 double (*grid)[numOfRows][numOfCols] = dataMatrix;
771 for (i = 0;i < numOfRows; ++i)
772 for (j = 0;j < numOfCols; ++j)
773 (*tmp)[j][i] = (*grid)[i][j];
774 memcpy(grid, tmp, sizeof(double) * numOfRows * numOfCols);
775 free(tmp);
776 }
777 } else if (!strcmp(encode, BCboolEncode)) {
778 if (isColumnMajor) {
779 BOOL (*tmp)[numOfRows][numOfCols] = malloc(sizeof(BOOL) * numOfRows * numOfCols);
780 BOOL (*grid)[numOfCols][numOfRows] = dataMatrix;
781 for (i = 0;i < numOfRows; ++i)
782 for (j = 0;j < numOfCols; ++j)
783 (*tmp)[i][j] = (*grid)[j][i];
784 memcpy(grid, tmp, sizeof(BOOL) * numOfRows * numOfCols);
785 free(tmp);
786 } else {
787 BOOL (*tmp)[numOfCols][numOfRows] = malloc(sizeof(BOOL) * numOfRows * numOfCols);
788 BOOL (*grid)[numOfRows][numOfCols] = dataMatrix;
789 for (i = 0;i < numOfRows; ++i)
790 for (j = 0;j < numOfCols; ++j)
791 (*tmp)[j][i] = (*grid)[i][j];
792 memcpy(grid, tmp, sizeof(BOOL) * numOfRows * numOfCols);
793 free(tmp);
794 }
795 }
796
797 isColumnMajor = aFlag;
798 }
799
800 - (NSArray *)rowNames { return rowNames; }
801 - (void)setRowNames:(NSArray *)anArray { [rowNames release]; rowNames = [anArray retain]; }
802 - (NSArray *)columnNames { return colNames; }
803 - (void)setColumnNames:(NSArray *)anArray { [colNames release]; colNames = [anArray retain]; }
804
805 #if 0
806 #pragma mark == MATRIX OPERATIONS ==
807 #endif
808
809 - (BCDataMatrix *)dataMatrixFromRowRange: (NSRange)rows andColumnRange: (NSRange)cols
810 {
811 // check valid ranges
812 if ((rows.length == 0) || (cols.length == 0)) return nil;
813 if ((rows.location + rows.length) > numOfRows) return nil;
814 if ((cols.location + cols.length) > numOfCols) return nil;
815
816 BCDataMatrix *newMatrix = [BCDataMatrix emptyDataMatrixWithRows: rows.length
817 andColumns: cols.length andEncode: encode];
818
819 int i, j;
820 if (!strcmp(encode, BCidEncode)) {
821 id (*grid1)[numOfRows][numOfCols] = dataMatrix;
822 id (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
823 for (i = 0;i < rows.length; ++i)
824 for (j = 0;j < cols.length; ++j)
825 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
826 } else if (!strcmp(encode, BCintEncode)) {
827 int (*grid1)[numOfRows][numOfCols] = dataMatrix;
828 int (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
829 for (i = 0;i < rows.length; ++i)
830 for (j = 0;j < cols.length; ++j)
831 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
832 } else if (!strcmp(encode, BClongEncode)) {
833 long (*grid1)[numOfRows][numOfCols] = dataMatrix;
834 long (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
835 for (i = 0;i < rows.length; ++i)
836 for (j = 0;j < cols.length; ++j)
837 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
838 } else if (!strcmp(encode, BCfloatEncode)) {
839 float (*grid1)[numOfRows][numOfCols] = dataMatrix;
840 float (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
841 for (i = 0;i < rows.length; ++i)
842 for (j = 0;j < cols.length; ++j)
843 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
844 } else if (!strcmp(encode, BCdoubleEncode)) {
845 double (*grid1)[numOfRows][numOfCols] = dataMatrix;
846 double (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
847 for (i = 0;i < rows.length; ++i)
848 for (j = 0;j < cols.length; ++j)
849 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
850 } else if (!strcmp(encode, BCboolEncode)) {
851 BOOL (*grid1)[numOfRows][numOfCols] = dataMatrix;
852 BOOL (*grid2)[rows.length][cols.length] = [newMatrix dataMatrix];
853 for (i = 0;i < rows.length; ++i)
854 for (j = 0;j < cols.length; ++j)
855 (*grid2)[i][j] = (*grid1)[i + rows.location][j + cols.location];
856 }
857
858 return newMatrix;
859 }
860
861 - (BCDataMatrix *)dataMatrixFromTranspose
862 {
863 BCDataMatrix *newMatrix = [BCDataMatrix emptyDataMatrixWithRows: numOfCols
864 andColumns: numOfRows andEncode: encode];
865
866 int i, j;
867 if (!strcmp(encode, BCidEncode)) {
868 id (*grid1)[numOfRows][numOfCols] = dataMatrix;
869 id (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
870 for (i = 0;i < numOfRows; ++i)
871 for (j = 0;j < numOfCols; ++j)
872 (*grid2)[j][i] = (*grid1)[i][j];
873 } else if (!strcmp(encode, BCintEncode)) {
874 int (*grid1)[numOfRows][numOfCols] = dataMatrix;
875 int (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
876 for (i = 0;i < numOfRows; ++i)
877 for (j = 0;j < numOfCols; ++j)
878 (*grid2)[j][i] = (*grid1)[i][j];
879 } else if (!strcmp(encode, BClongEncode)) {
880 long (*grid1)[numOfRows][numOfCols] = dataMatrix;
881 long (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
882 for (i = 0;i < numOfRows; ++i)
883 for (j = 0;j < numOfCols; ++j)
884 (*grid2)[j][i] = (*grid1)[i][j];
885 } else if (!strcmp(encode, BCfloatEncode)) {
886 float (*grid1)[numOfRows][numOfCols] = dataMatrix;
887 float (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
888 for (i = 0;i < numOfRows; ++i)
889 for (j = 0;j < numOfCols; ++j)
890 (*grid2)[j][i] = (*grid1)[i][j];
891 } else if (!strcmp(encode, BCdoubleEncode)) {
892 double (*grid1)[numOfRows][numOfCols] = dataMatrix;
893 double (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
894 for (i = 0;i < numOfRows; ++i)
895 for (j = 0;j < numOfCols; ++j)
896 (*grid2)[j][i] = (*grid1)[i][j];
897 } else if (!strcmp(encode, BCboolEncode)) {
898 BOOL (*grid1)[numOfRows][numOfCols] = dataMatrix;
899 BOOL (*grid2)[numOfCols][numOfRows] = [newMatrix dataMatrix];
900 for (i = 0;i < numOfRows; ++i)
901 for (j = 0;j < numOfCols; ++j)
902 (*grid2)[j][i] = (*grid1)[i][j];
903 }
904
905 if (rowNames) [newMatrix setColumnNames: [NSArray arrayWithArray: rowNames]];
906 if (colNames) [newMatrix setRowNames: [NSArray arrayWithArray: colNames]];
907
908 return newMatrix;
909
910 }
911
912 - (BOOL)writeToFile: (NSString *)fileName
913 {
914 return [self writeToFile: fileName withFormat: nil];
915 }
916
917 - (BOOL)writeToFile: (NSString *)fileName withFormat: (NSDictionary *)aFormat
918 {
919 BOOL res = NO;
920
921 if (!aFormat) {
922 FILE *aFile = fopen([fileName UTF8String], "w");
923 if (!aFile) return NO;
924
925 double (*A)[numOfRows][numOfCols] = dataMatrix;
926 int i, j;
927
928 if (colNames) {
929 if (rowNames) fprintf(aFile, "row ");
930 for (i = 0; i < [colNames count]; ++i) {
931 if (i != 0) fprintf(aFile, " ");
932 fprintf(aFile, "%s", [[colNames objectAtIndex: i] UTF8String]);
933 }
934 fprintf(aFile, "\n");
935 }
936
937 for (i = 0; i < numOfRows; ++i) {
938 if (rowNames) fprintf(aFile, "%s ", [[rowNames objectAtIndex: i] UTF8String]);
939 fprintf(aFile, "%lf", (*A)[i][0]);
940 for (j = 1; j < numOfCols; ++j) {
941 fprintf(aFile, " %lf", (*A)[i][j]);
942 }
943 fprintf(aFile, "\n");
944 }
945 //fprintf(aFile, "\n");
946
947 fclose(aFile);
948 res = YES;
949 } else {
950 }
951
952 return res;
953 }
954
955 - (void)prettyPrintMatrix
956 {
957 int i, j;
958 if (!strcmp(encode, BCidEncode)) {
959 } else if (!strcmp(encode, BCintEncode)) {
960 int (*grid1)[numOfRows][numOfCols] = dataMatrix;
961 for (i = 0;i < numOfRows; ++i) {
962 for (j = 0;j < numOfCols; ++j) {
963 if (j != 0) printf(" ");
964 printf("%d", (*grid1)[i][j]);
965 }
966 printf("\n");
967 }
968 } else if (!strcmp(encode, BClongEncode)) {
969 long (*grid1)[numOfRows][numOfCols] = dataMatrix;
970 for (i = 0;i < numOfRows; ++i) {
971 for (j = 0;j < numOfCols; ++j) {
972 if (j != 0) printf(" ");
973 printf("%ld", (*grid1)[i][j]);
974 }
975 printf("\n");
976 }
977 } else if (!strcmp(encode, BCfloatEncode)) {
978 float (*grid1)[numOfRows][numOfCols] = dataMatrix;
979 for (i = 0;i < numOfRows; ++i) {
980 for (j = 0;j < numOfCols; ++j) {
981 if (j != 0) printf(" ");
982 printf("%f", (*grid1)[i][j]);
983 }
984 printf("\n");
985 }
986 } else if (!strcmp(encode, BCdoubleEncode)) {
987 double (*grid1)[numOfRows][numOfCols] = dataMatrix;
988 for (i = 0;i < numOfRows; ++i) {
989 for (j = 0;j < numOfCols; ++j) {
990 if (j != 0) printf(" ");
991 printf("%lf", (*grid1)[i][j]);
992 }
993 printf("\n");
994 }
995 } else if (!strcmp(encode, BCboolEncode)) {
996 BOOL (*grid1)[numOfRows][numOfCols] = dataMatrix;
997 for (i = 0;i < numOfRows; ++i) {
998 for (j = 0;j < numOfCols; ++j) {
999 if (j != 0) printf(" ");
1000 if ((*grid1)[i][j]) printf("YES");
1001 else printf(" NO");
1002 }
1003 printf("\n");
1004 }
1005 }
1006
1007 }
1008
1009 #if 0
1010 #pragma mark == LINEAR ALGEBRA OPERATIONS ==
1011 #endif
1012
1013 #if 0
1014 - (BCDataMatrix *)LUfactorize
1015 {
1016 #ifdef GNUSTEP
1017 // create copy of data matrix
1018 BCDataMatrix *m = [BCDataMatrix dataMatrixWithDataMatrix: self andEncode: BCdoubleEncode isColumnMajor: NO];
1019
1020 // create a GSL matrix suitable for LU factorize
1021 gsl_matrix *gm = gsl_matrix_alloc(numOfRows, numOfCols);
1022 if (!gm) {
1023 printf("ERROR: Could not allocate GSL matrix.\n");
1024 return nil;
1025 }
1026
1027 int i, j;
1028 double (*grid)[numOfRows][numOfCols] = [m dataMatrix];
1029 for (i = 0; i < numOfRows; ++i)
1030 for (j = 0; j < numOfCols; ++j)
1031 gsl_matrix_set(gm, i, j, (*grid)[i][j]);
1032
1033 //gsl_matrix_fprintf(stdout, gm, "%lf");
1034
1035 int s;
1036 gsl_permutation *p = gsl_permutation_alloc(numOfRows);
1037 gsl_linalg_LU_decomp(gm, p, &s);
1038
1039 //gsl_matrix_fprintf(stdout, gm, "%lf");
1040
1041 for (i = 0; i < numOfRows; ++i)
1042 for (j = 0; j < numOfCols; ++j)
1043 (*grid)[i][j] = gsl_matrix_get(gm, i, j);
1044
1045 gsl_permutation_free(p);
1046 gsl_matrix_free(gm);
1047
1048 return m;
1049
1050 #else
1051
1052 // create copy of data matrix suitable for LU factorize
1053 BCDataMatrix *m = [BCDataMatrix dataMatrixWithDataMatrix: self andEncode: BCdoubleEncode isColumnMajor: YES];
1054
1055 __CLPK_integer M = numOfRows;
1056 __CLPK_integer N = numOfCols;
1057 __CLPK_integer INFO;
1058 __CLPK_integer ipiv[3];
1059 double *a = [m dataMatrix];
1060
1061 dgetrf_(&M, &N, a, &M, ipiv, &INFO);
1062
1063 // set result matrix to same data format
1064 [m setColumnMajor: [self isColumnMajor]];
1065
1066 return m;
1067 #endif
1068 }
1069 #endif
1070
1071 - (double)determinant
1072 {
1073 if (numOfCols != numOfRows) {
1074 printf("ERROR: Can only calculate determinant for square matrix.\n");
1075 return nan("");
1076 }
1077
1078 if (!strcmp(encode, BCidEncode)) {
1079 printf("ERROR: determinant only valid for numerical matrixes.\n");
1080 return nan("");
1081 } else if (!strcmp(encode, BCboolEncode)) {
1082 printf("ERROR: determinant only valid for numerical matrixes.\n");
1083 return nan("");
1084 }
1085
1086 BCDataMatrix *m = [self LUfactorize];
1087 double (*grid)[numOfCols][numOfRows] = [m dataMatrix];
1088 double det = 1.0;
1089 int i;
1090 for (i = 0; i < numOfCols; ++i) {
1091 det *= (*grid)[i][i];
1092 }
1093
1094 return det;
1095 }
1096
1097 @end