← Back to branch summary

~ubuntu-branches/ubuntu/trusty/3depict/trusty

« back to all changes in this revision

Viewing changes to src/backend/APT/APTClasses.cpp

  • Committer: Package Import Robot
  • Author(s): D Haley
  • Date: 2013-05-17 00:52:39 UTC
  • mfrom: (3.1.4 experimental)
  • Revision ID: package-import@ubuntu.com-20130517005239-7bl4mnhkvrhc2ba6
Tags: 0.0.13-1
Upload to unstable 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/backend/APT/APTClasses.cpp
 
1
/* 
 
2
 * APTClasses.h - Generic APT components code
 
3
 * Copyright (C) 2013  D Haley
 
4
 * 
 
5
 * This program is free software: you can redistribute it and/or modify
 
6
 * it under the terms of the GNU General Public License as published by
 
7
 * the Free Software Foundation, either version 3 of the License, or
 
8
 * (at your option) any later version.
 
9
 * 
 
10
 * This program is distributed in the hope that it will be useful,
 
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 * GNU General Public License for more details.
 
14
 * 
 
15
 * You should have received a copy of the GNU General Public License
 
16
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
17
 */
 
18
 
 
19
#include "APTClasses.h"
 
20
#include "../../common/stringFuncs.h"
 
21
 
 
22
#include "../../common/translation.h"
 
23
 
 
24
 
 
25
using std::pair;
 
26
using std::string;
 
27
using std::vector;
 
28
using std::ifstream;
 
29
using std::make_pair;
 
30
 
 
31
const char *POS_ERR_STRINGS[] = { "",
 
32
                                NTRANS("Memory allocation failure on POS load"),
 
33
                                NTRANS("Error opening pos file"),
 
34
                                NTRANS("Pos file empty"),
 
35
                                NTRANS("Pos file size appears to have non-integer number of entries"),
 
36
                                NTRANS("Error reading from pos file (after open)"),
 
37
                                NTRANS("Error - Found NaN in pos file"),
 
38
                                NTRANS("Pos load aborted by interrupt.")
 
39
};
 
40
 
 
41
enum
 
42
{
 
43
        TEXT_ERR_OPEN=1,
 
44
        TEXT_ERR_ONLY_HEADER,
 
45
        TEXT_ERR_REOPEN,
 
46
        TEXT_ERR_READ_CONTENTS,
 
47
        TEXT_ERR_FORMAT,
 
48
        TEXT_ERR_NUM_FIELDS,
 
49
        TEXT_ERR_ALLOC_FAIL,
 
50
        TEXT_ERR_ENUM_END //not an error, just end of enum
 
51
};
 
52
 
 
53
const char *ION_TEXT_ERR_STRINGS[] = { "",
 
54
                                        NTRANS("Error opening file"),
 
55
                                        NTRANS("No numerical data found"),
 
56
                                        NTRANS("Error re-opening file, after first scan"),
 
57
                                        NTRANS("Unable to read file contents after open"),
 
58
                                        NTRANS("Error interpreting field in file"),
 
59
                                        NTRANS("Incorrect number of fields in file"),
 
60
                                        NTRANS("Unable to allocate memory to store data"),
 
61
                                        };
 
62
//!Create an pos file from a vector of IonHits
 
63
unsigned int IonVectorToPos(const vector<IonHit> &ionVec, const string &filename)
 
64
{
 
65
        std::ofstream CFile(filename.c_str(),std::ios::binary);
 
66
        float floatBuffer[4];
 
67
 
 
68
        if (!CFile)
 
69
                return 1;
 
70
 
 
71
        for (unsigned int ui=0; ui<ionVec.size(); ui++)
 
72
        {
 
73
                ionVec[ui].makePosData(floatBuffer);
 
74
                CFile.write((char *)floatBuffer,4*sizeof(float));
 
75
        }
 
76
        return 0;
 
77
}
 
78
 
 
79
 
 
80
 
 
81
void appendPos(const vector<IonHit> &points, const char *name)
 
82
{
 
83
        std::ofstream posFile(name,std::ios::binary|std::ios::app);     
 
84
 
 
85
        float data[4];  
 
86
        
 
87
        for(unsigned int ui=0; ui< points.size(); ui++)
 
88
        {
 
89
                points[ui].makePosData(data);
 
90
                posFile.write((char *)data, 4*sizeof(float));
 
91
        }
 
92
}
 
93
 
 
94
void getPointsFromIons(const vector<IonHit> &ions, vector<Point3D> &p)
 
95
{
 
96
        p.resize(ions.size());
 
97
#pragma omp parallel for
 
98
        for(size_t ui=0;ui<ions.size();ui++)
 
99
                p[ui] = ions[ui].getPosRef();
 
100
}
 
101
 
 
102
unsigned int LimitLoadPosFile(unsigned int inputnumcols, unsigned int outputnumcols, unsigned int index[], vector<IonHit> &posIons,const char *posFile, size_t limitCount,
 
103
                unsigned int &progress, bool (*callback)(bool),bool strongSampling)
 
104
{
 
105
 
 
106
        //Function is only defined for 4 columns here.
 
107
        ASSERT(outputnumcols == 4);
 
108
        //buffersize must be a power of two and at least outputnumcols*sizeof(float)
 
109
        const unsigned int NUMROWS=1;
 
110
        const unsigned int BUFFERSIZE=inputnumcols * sizeof(float) * NUMROWS;
 
111
        const unsigned int BUFFERSIZE2=outputnumcols * sizeof(float) * NUMROWS;
 
112
        char *buffer=new char[BUFFERSIZE];
 
113
 
 
114
        
 
115
        if(!buffer)
 
116
                return POS_ALLOC_FAIL;
 
117
 
 
118
        char *buffer2=new char[BUFFERSIZE2];
 
119
        if(!buffer2)
 
120
        {
 
121
                delete[] buffer;
 
122
                return POS_ALLOC_FAIL;
 
123
        }
 
124
 
 
125
        //open pos file
 
126
        std::ifstream CFile(posFile,std::ios::binary);
 
127
 
 
128
        if(!CFile)
 
129
        {
 
130
                delete[] buffer;
 
131
                delete[] buffer2;
 
132
                return POS_OPEN_FAIL;
 
133
        }
 
134
        
 
135
        CFile.seekg(0,std::ios::end);
 
136
        size_t fileSize=CFile.tellg();
 
137
 
 
138
        if(!fileSize)
 
139
        {
 
140
                delete[] buffer;
 
141
                delete[] buffer2;
 
142
                return POS_EMPTY_FAIL;
 
143
        }
 
144
        
 
145
        CFile.seekg(0,std::ios::beg);
 
146
        
 
147
        //calculate the number of points stored in the POS file
 
148
        size_t pointCount=0;
 
149
        size_t maxIons;
 
150
        size_t maxCols = inputnumcols * sizeof(float);
 
151
        //regular case
 
152
        
 
153
        if(fileSize % maxCols)
 
154
        {
 
155
                delete[] buffer;
 
156
                delete[] buffer2;
 
157
                return POS_SIZE_MODULUS_ERR;    
 
158
        }
 
159
 
 
160
        maxIons =fileSize/maxCols;
 
161
        limitCount=std::min(limitCount,maxIons);
 
162
 
 
163
        //If we are going to load the whole file, don't use a sampling method to do it.
 
164
        if(limitCount == maxIons)
 
165
        {
 
166
                //Close the file
 
167
                CFile.close();
 
168
                delete[] buffer;
 
169
                delete[] buffer2;
 
170
                //Try opening it using the normal functions
 
171
                return GenericLoadFloatFile(inputnumcols, outputnumcols, index, posIons,posFile,progress, callback);
 
172
        }
 
173
 
 
174
        //Use a sampling method to load the pos file
 
175
        std::vector<size_t> ionsToLoad;
 
176
        try
 
177
        {
 
178
                posIons.resize(limitCount);
 
179
 
 
180
                RandNumGen rng;
 
181
                rng.initTimer();
 
182
                unsigned int dummy;
 
183
                randomDigitSelection(ionsToLoad,maxIons,rng,
 
184
                                limitCount,dummy,callback,strongSampling);
 
185
        }
 
186
        catch(std::bad_alloc)
 
187
        {
 
188
                delete[] buffer;
 
189
                delete[] buffer2;
 
190
                return POS_ALLOC_FAIL;
 
191
        }
 
192
 
 
193
 
 
194
        //sort again    
 
195
        GreaterWithCallback<size_t> g(callback,PROGRESS_REDUCE);
 
196
        std::sort(ionsToLoad.begin(),ionsToLoad.end(),g);
 
197
 
 
198
        unsigned int curProg = PROGRESS_REDUCE; 
 
199
 
 
200
        //TODO: probably not very nice to the disk drive. would be better to
 
201
        //scan ahead for contiguous data regions, and load that where possible.
 
202
        //Or switch between different algorithms based upon ionsToLoad.size()/  
 
203
        std::ios::pos_type  nextIonPos;
 
204
        for(size_t ui=0;ui<ionsToLoad.size(); ui++)
 
205
        {
 
206
                nextIonPos =  ionsToLoad[ui]*maxCols;
 
207
                
 
208
                if(CFile.tellg() !=nextIonPos )
 
209
                        CFile.seekg(nextIonPos);
 
210
 
 
211
                CFile.read(buffer,BUFFERSIZE);
 
212
 
 
213
                for (size_t i = 0; i < outputnumcols; i++) // iterate through floats
 
214
                        memcpy(&(buffer2[i * sizeof(float)]), &(buffer[index[i] * sizeof(float)]), sizeof(float));
 
215
                
 
216
                if(!CFile.good())
 
217
                        return POS_READ_FAIL;
 
218
                posIons[ui].setHit((float*)buffer2);
 
219
                //Data bytes stored in pos files are big
 
220
                //endian. flip as required
 
221
                #ifdef __LITTLE_ENDIAN__
 
222
                        posIons[ui].switchEndian();     
 
223
                #endif
 
224
                
 
225
                if(posIons[ui].hasNaN())
 
226
                {
 
227
                        delete[] buffer;
 
228
                        delete[] buffer2;
 
229
                        return POS_NAN_LOAD_ERROR;      
 
230
                }
 
231
                        
 
232
                pointCount++;
 
233
                if(!curProg--)
 
234
                {
 
235
 
 
236
                        progress= (unsigned int)((float)(CFile.tellg())/((float)fileSize)*100.0f);
 
237
                        curProg=PROGRESS_REDUCE;
 
238
                        if(!(*callback)(false))
 
239
                        {
 
240
                                delete[] buffer;
 
241
                                posIons.clear();
 
242
                                return POS_ABORT_FAIL;
 
243
                                
 
244
                        }
 
245
                }
 
246
                                
 
247
        }
 
248
 
 
249
        delete[] buffer;
 
250
        delete[] buffer2;
 
251
        return 0;
 
252
}
 
253
 
 
254
unsigned int GenericLoadFloatFile(unsigned int inputnumcols, unsigned int outputnumcols, 
 
255
                unsigned int index[], vector<IonHit> &posIons,const char *posFile, 
 
256
                        unsigned int &progress, bool (*callback)(bool))
 
257
{
 
258
        ASSERT(outputnumcols==4); //Due to ionHit.setHit
 
259
        //buffersize must be a power of two and at least sizeof(float)*outputnumCols
 
260
        const unsigned int NUMROWS=512;
 
261
        const unsigned int BUFFERSIZE=inputnumcols * sizeof(float) * NUMROWS;
 
262
        const unsigned int BUFFERSIZE2=outputnumcols * sizeof(float) * NUMROWS;
 
263
 
 
264
        char *buffer=new char[BUFFERSIZE];
 
265
        
 
266
        if(!buffer)
 
267
                return POS_ALLOC_FAIL;
 
268
        
 
269
        char *buffer2=new char[BUFFERSIZE2];
 
270
        if(!buffer2)
 
271
        {
 
272
                delete[] buffer;
 
273
                return POS_ALLOC_FAIL;
 
274
        }
 
275
        //open pos file
 
276
        std::ifstream CFile(posFile,std::ios::binary);
 
277
        
 
278
        if(!CFile)
 
279
        {
 
280
                delete[] buffer;
 
281
                delete[] buffer2;
 
282
                return POS_OPEN_FAIL;
 
283
        }
 
284
        
 
285
        CFile.seekg(0,std::ios::end);
 
286
        size_t fileSize=CFile.tellg();
 
287
        
 
288
        if(!fileSize)
 
289
        {
 
290
                delete[] buffer;
 
291
                delete[] buffer2;
 
292
                return POS_EMPTY_FAIL;
 
293
        }
 
294
        
 
295
        CFile.seekg(0,std::ios::beg);
 
296
        
 
297
        //calculate the number of points stored in the POS file
 
298
        IonHit hit;
 
299
        typedef struct IONHIT
 
300
        {
 
301
                float pos[3];
 
302
                float massToCharge;
 
303
        } IONHIT;
 
304
        size_t pointCount=0;
 
305
        //regular case
 
306
        size_t curBufferSize=BUFFERSIZE;
 
307
        size_t curBufferSize2=BUFFERSIZE2;
 
308
        
 
309
        if(fileSize % (inputnumcols * sizeof(float)))
 
310
        {
 
311
                delete[] buffer;
 
312
                delete[] buffer2;
 
313
                return POS_SIZE_MODULUS_ERR;    
 
314
        }
 
315
        
 
316
        try
 
317
        {
 
318
                posIons.resize(fileSize/(inputnumcols*sizeof(float)));
 
319
        }
 
320
        catch(std::bad_alloc)
 
321
        {
 
322
                delete[] buffer;
 
323
                delete[] buffer2;
 
324
                return POS_ALLOC_FAIL;
 
325
        }
 
326
        
 
327
        
 
328
        while(fileSize < curBufferSize) {
 
329
                curBufferSize = curBufferSize >> 1;
 
330
                curBufferSize2 = curBufferSize2 >> 1;
 
331
        }               
 
332
        
 
333
        //Technically this is dependent upon the buffer size.
 
334
        unsigned int curProg = 10000;   
 
335
        size_t ionP=0;
 
336
        int maxCols = inputnumcols * sizeof(float);
 
337
        int maxPosCols = outputnumcols * sizeof(float);
 
338
        do
 
339
        {
 
340
                //Taking curBufferSize chunks at a time, read the input file
 
341
                while((size_t)CFile.tellg() <= fileSize-curBufferSize)
 
342
                {
 
343
                        CFile.read(buffer,curBufferSize);
 
344
                        if(!CFile.good())
 
345
                        {
 
346
                                delete[] buffer;
 
347
                                delete[] buffer2;
 
348
                                return POS_READ_FAIL;
 
349
                        }
 
350
                        
 
351
                        for (unsigned int j = 0; j < NUMROWS; j++) // iterate through rows
 
352
                        {
 
353
                                for (unsigned int i = 0; i < outputnumcols; i++) // iterate through floats
 
354
                                {
 
355
                                        memcpy(&(buffer2[j * maxPosCols + i * sizeof(float)]), 
 
356
                                                &(buffer[j * maxCols + index[i] * sizeof(float)]), sizeof(float));
 
357
                                }
 
358
                        }
 
359
                        
 
360
                        unsigned int ui;
 
361
                        for(ui=0; ui<curBufferSize2; ui+=(sizeof(IONHIT)))
 
362
                        {
 
363
                                hit.setHit((float*)(buffer2+ui));
 
364
                                //Data bytes stored in pos files are big
 
365
                                //endian. flip as required
 
366
                                #ifdef __LITTLE_ENDIAN__
 
367
                                        hit.switchEndian();     
 
368
                                #endif
 
369
                                
 
370
                                if(hit.hasNaN())
 
371
                                {
 
372
                                        delete[] buffer;
 
373
                                        delete[] buffer2;
 
374
                                        return POS_NAN_LOAD_ERROR;      
 
375
                                }
 
376
                                posIons[ionP] = hit;
 
377
                                ionP++;
 
378
                                
 
379
                                pointCount++;
 
380
                        }       
 
381
                        
 
382
                        if(!curProg--)
 
383
                        {
 
384
                                progress= (unsigned int)((float)(CFile.tellg())/((float)fileSize)*100.0f);
 
385
                                curProg=PROGRESS_REDUCE;
 
386
                                if(!(*callback)(false))
 
387
                                {
 
388
                                        delete[] buffer;
 
389
                                        delete[] buffer2;
 
390
                                        posIons.clear();
 
391
                                        return POS_ABORT_FAIL;
 
392
                                
 
393
                                }
 
394
                        }
 
395
                                
 
396
                }
 
397
 
 
398
                curBufferSize = curBufferSize >> 1 ;
 
399
                curBufferSize2 = curBufferSize2 >> 1 ;
 
400
        }while(curBufferSize2 >= sizeof(IONHIT));
 
401
        
 
402
        ASSERT((unsigned int)CFile.tellg() == fileSize);
 
403
        delete[] buffer;
 
404
        delete[] buffer2;
 
405
        
 
406
        return 0;
 
407
}
 
408
 
 
409
 
 
410
//TODO: Add progress
 
411
unsigned int limitLoadTextFile(unsigned int numColsTotal, unsigned int selectedCols[], 
 
412
                        vector<IonHit> &posIons,const char *textFile, const char *delim, const size_t limitCount,
 
413
                                unsigned int &progress, bool (*callback)(bool),bool strongRandom)
 
414
{
 
415
 
 
416
        ASSERT(numColsTotal);
 
417
        ASSERT(textFile);
 
418
 
 
419
        vector<size_t> newLinePositions;
 
420
        std::vector<std::string> subStrs;
 
421
 
 
422
        //Do a brute force scan through the dataset
 
423
        //to locate newlines.
 
424
        char *buffer;
 
425
        const int BUFFER_SIZE=16384; //This is totally a guess. I don't know what is best.
 
426
 
 
427
 
 
428
        //sort the selected columns into increasing order
 
429
        vector<int> sortedCols;
 
430
        for(unsigned int ui=0;ui<numColsTotal;ui++)
 
431
                sortedCols.push_back(selectedCols[ui]);
 
432
 
 
433
        std::sort(sortedCols.begin(),sortedCols.end());
 
434
        //the last entry in the sorted list tells us how many
 
435
        // entries we need in each line at a minimum
 
436
        unsigned int maxCol=sortedCols[sortedCols.size()-1];
 
437
 
 
438
        ifstream CFile(textFile,std::ios::binary);
 
439
 
 
440
        if(!CFile)
 
441
                return TEXT_ERR_OPEN;
 
442
 
 
443
 
 
444
        //seek to the end of the file
 
445
        //to get the filesize
 
446
        size_t maxPos,curPos;
 
447
        CFile.seekg(0,std::ios::end);
 
448
        maxPos=CFile.tellg();
 
449
 
 
450
        CFile.close();
 
451
 
 
452
        CFile.open(textFile);
 
453
 
 
454
        curPos=0;
 
455
 
 
456
 
 
457
        //Scan through file for end of header.
 
458
        //we define this as the split value being able to generate 
 
459
        //1) Enough data to make interpretable columns
 
460
        //2) Enough columns that can be interpreted.
 
461
        while(!CFile.eof() && curPos < maxPos)
 
462
        {
 
463
                string s;
 
464
                curPos = CFile.tellg();
 
465
                getline(CFile,s);
 
466
 
 
467
                if(!CFile.good())
 
468
                        return TEXT_ERR_READ_CONTENTS;
 
469
 
 
470
                splitStrsRef(s.c_str(),delim,subStrs);  
 
471
                stripZeroEntries(subStrs);
 
472
 
 
473
                //Not enough entries in this line
 
474
                //to be interpretable
 
475
                if(subStrs.size() <maxCol)
 
476
                        continue;
 
477
 
 
478
                bool enoughSubStrs;
 
479
                enoughSubStrs=true;
 
480
                for(unsigned int ui=0;ui<numColsTotal; ui++)
 
481
                {
 
482
                        if(selectedCols[ui] >=subStrs.size())
 
483
                        {
 
484
                                enoughSubStrs=false;
 
485
                                break;  
 
486
                        }
 
487
                }
 
488
 
 
489
                if(!enoughSubStrs)
 
490
                        continue;
 
491
 
 
492
                //Unable to stream
 
493
                bool unStreamable;
 
494
                unStreamable=false;
 
495
                for(unsigned int ui=0; ui<numColsTotal; ui++)
 
496
                {
 
497
                        float f;
 
498
                        if(stream_cast(f,subStrs[selectedCols[ui]]))
 
499
                        {
 
500
                                unStreamable=true;
 
501
                                break;
 
502
                        }
 
503
 
 
504
                }
 
505
 
 
506
                //well, we can stream this, so it assume it is not the header.
 
507
                if(!unStreamable)
 
508
                        break;  
 
509
        }       
 
510
 
 
511
        //could not find any data.. only header.
 
512
        if(CFile.eof() || curPos >=maxPos)
 
513
                return TEXT_ERR_ONLY_HEADER;
 
514
 
 
515
 
 
516
        CFile.close();
 
517
 
 
518
        //Re-open the file in binary mode to find the newlines
 
519
        CFile.open(textFile,std::ios::binary);
 
520
 
 
521
        if(!CFile)
 
522
                return TEXT_ERR_REOPEN;
 
523
 
 
524
 
 
525
        //Jump to beyond the header
 
526
        CFile.seekg(curPos);
 
527
 
 
528
 
 
529
        //keep a beginning of file marker
 
530
        newLinePositions.push_back(curPos);
 
531
        bool seenNumeric=false;
 
532
        buffer = new char[BUFFER_SIZE];
 
533
        while(!CFile.eof() && curPos < maxPos)
 
534
        {
 
535
                size_t bytesToRead;
 
536
 
 
537
                if(!CFile.good())
 
538
                {
 
539
                        delete[] buffer;
 
540
                        return TEXT_ERR_READ_CONTENTS;
 
541
                }
 
542
                //read up to BUFFER_SIZE bytes from the file
 
543
                //but only if they are available
 
544
                bytesToRead = std::min(maxPos-curPos,(size_t)BUFFER_SIZE);
 
545
 
 
546
                CFile.read(buffer,bytesToRead);
 
547
 
 
548
                //check that this buffer contains numeric info  
 
549
                for(unsigned int ui=0;ui<bytesToRead; ui++)
 
550
                {
 
551
                        //Check for a unix-style endline
 
552
                        //or the latter part of a windows endline
 
553
                        //either or, whatever.
 
554
                        if( buffer[ui] == '\n')
 
555
                        {
 
556
                                //Check that we have not hit a run of non-numeric data
 
557
                                if(seenNumeric)
 
558
                                        newLinePositions.push_back(ui+curPos);
 
559
                        } 
 
560
                        else if(buffer[ui] >= '0' && buffer[ui] <='9')
 
561
                                seenNumeric=true;
 
562
                }
 
563
                
 
564
                curPos+=bytesToRead;    
 
565
        
 
566
        }
 
567
        
 
568
        //Don't keep any newline that hits the end of the file, but do keep a zero position
 
569
        if(newLinePositions.size())
 
570
                newLinePositions.pop_back();
 
571
        CFile.close();
 
572
 
 
573
 
 
574
        //OK, so now we know where those pesky endlines are. This gives us jump positions
 
575
        //to new lines in the file. Each component must have some form of numeric data 
 
576
        //preceding it. That numeric data may not be fully parseable, but we assume it is until we know better.
 
577
        //
 
578
        //If it is *not* parseable, just throw an error when we find that out.
 
579
 
 
580
        //If we are going to load the whole file, don't use a sampling method to do it.
 
581
        if(limitCount >=newLinePositions.size())
 
582
        {
 
583
                delete[] buffer;
 
584
 
 
585
                vector<vector<float> > data;
 
586
                vector<string> header;
 
587
        
 
588
                //Just use the non-sampling method to load.     
 
589
                if(loadTextData(textFile,data,header,delim))
 
590
                        return TEXT_ERR_FORMAT;
 
591
 
 
592
                if(data.size() !=4)
 
593
                        return TEXT_ERR_NUM_FIELDS;
 
594
 
 
595
                posIons.resize(data[0].size());
 
596
                for(size_t ui=0;ui<data[0].size();ui++)
 
597
                {
 
598
 
 
599
                        ASSERT(numColsTotal==4);
 
600
                        //This is output specific
 
601
                        //and assumes that we have exactly 4 input cols
 
602
                        //to match ot our output vector of pos ions
 
603
                        posIons[ui].setPos(Point3D(data[selectedCols[0]][ui],
 
604
                                                data[selectedCols[1]][ui],
 
605
                                                data[selectedCols[2]][ui]));
 
606
                        posIons[ui].setMassToCharge(data[selectedCols[3]][ui]);
 
607
                }
 
608
 
 
609
                return 0;       
 
610
        }
 
611
 
 
612
        //Generate some random positions to load
 
613
        std::vector<size_t> dataToLoad;
 
614
        try
 
615
        {
 
616
                posIons.resize(limitCount);
 
617
 
 
618
                RandNumGen rng;
 
619
                rng.initTimer();
 
620
                unsigned int dummy;
 
621
                randomDigitSelection(dataToLoad,newLinePositions.size(),rng,
 
622
                                limitCount,dummy,callback,strongRandom);
 
623
        }
 
624
        catch(std::bad_alloc)
 
625
        {
 
626
                delete[] buffer;
 
627
                return TEXT_ERR_ALLOC_FAIL;
 
628
        }
 
629
 
 
630
 
 
631
        //Sort the data such that we are going to
 
632
        //always jump forwards in the file; better disk access and whatnot.
 
633
        GreaterWithCallback<size_t> g(callback,PROGRESS_REDUCE);
 
634
        std::sort(dataToLoad.begin(),dataToLoad.end(),g);
 
635
 
 
636
        //OK, so we have  a list of newlines
 
637
        //that we can use as entry points for random seek.
 
638
        //We also have some random entry points (sorted).
 
639
        // Now re-open the file in text mode and try to load the
 
640
        // data specified at the offsets
 
641
 
 
642
 
 
643
        //Open file in text mode
 
644
        CFile.open(textFile);
 
645
 
 
646
        if(!CFile)
 
647
        {
 
648
                delete[] buffer;
 
649
                return TEXT_ERR_REOPEN;
 
650
        }
 
651
 
 
652
 
 
653
        //OK, now jump to each of the random positions,
 
654
        //as dictated by the endline position 
 
655
        //and attempt a parsing there.
 
656
 
 
657
        subStrs.clear();
 
658
        for(size_t ui=0;ui<dataToLoad.size();ui++)      
 
659
        {
 
660
 
 
661
                std::ios::pos_type  nextIonPos;
 
662
                //Jump to position immediately after the newline
 
663
                nextIonPos = (newLinePositions[dataToLoad[ui]]+1);
 
664
                
 
665
                if(CFile.tellg() !=nextIonPos )
 
666
                        CFile.seekg(nextIonPos);
 
667
 
 
668
                string s;
 
669
 
 
670
                getline(CFile,s);
 
671
 
 
672
                //Now attempt to scan the line for the selected columns.
 
673
                //split around whatever delimiter we can find
 
674
                splitStrsRef(s.c_str(),delim,subStrs);  
 
675
 
 
676
                if(subStrs.size() < numColsTotal)
 
677
                {
 
678
                        //FIXME: Allow skipping of bad lines
 
679
                        delete[] buffer;
 
680
                        return TEXT_ERR_NUM_FIELDS;
 
681
                }
 
682
 
 
683
                float f[4];
 
684
 
 
685
                for(size_t uj=0;uj<sortedCols.size();uj++)
 
686
                {
 
687
 
 
688
                        if(stream_cast(f[uj],subStrs[sortedCols[uj]]))
 
689
                        {
 
690
                                //FIXME: Allow skipping bad lines
 
691
                                //Can't parse line.. Abort.
 
692
                                delete[] buffer;
 
693
                                return TEXT_ERR_FORMAT;
 
694
                        }
 
695
                }
 
696
 
 
697
                ASSERT(ui<posIons.size());
 
698
                posIons[ui].setHit(f);
 
699
        }
 
700
 
 
701
        delete[] buffer;
 
702
        return 0;
 
703
 
 
704
}
 
705
 
 
706
IonHit::IonHit() 
 
707
{
 
708
        //At this point I deliberately don't initialise the point class
 
709
        //as in DEBUG mode, the point class will catch failure to init
 
710
}
 
711
 
 
712
IonHit::IonHit(const IonHit &obj2) : massToCharge(obj2.massToCharge), pos(obj2.pos)
 
713
{
 
714
}
 
715
 
 
716
IonHit::IonHit(const Point3D &p, float newMass) : massToCharge(newMass), pos(p)
 
717
{
 
718
}
 
719
 
 
720
void IonHit::setMassToCharge(float newMass)
 
721
{
 
722
        massToCharge=newMass;
 
723
}
 
724
 
 
725
float IonHit::getMassToCharge() const
 
726
{
 
727
        return massToCharge;
 
728
}       
 
729
 
 
730
 
 
731
void IonHit::setPos(const Point3D &p)
 
732
{
 
733
        pos=p;
 
734
}
 
735
 
 
736
#ifdef __LITTLE_ENDIAN__
 
737
void IonHit::switchEndian()
 
738
{
 
739
        
 
740
        pos.switchEndian();
 
741
        floatSwapBytes(&(massToCharge));
 
742
}
 
743
#endif
 
744
 
 
745
const IonHit &IonHit::operator=(const IonHit &obj)
 
746
{
 
747
        massToCharge=obj.massToCharge;
 
748
        pos = obj.pos;
 
749
 
 
750
        return *this;
 
751
}
 
752
 
 
753
IonHit IonHit::operator+(const Point3D &obj)
 
754
{
 
755
        //FIXME: I think this is wrong???
 
756
        ASSERT(false);
 
757
        pos.add(obj);   
 
758
        return *this;
 
759
}
 
760
 
 
761
void IonHit::makePosData(float *floatArr) const
 
762
{
 
763
        ASSERT(floatArr);
 
764
        //copy positional information
 
765
        pos.copyValueArr(floatArr);
 
766
 
 
767
        //copy mass to charge data
 
768
        *(floatArr+3) = massToCharge;
 
769
                
 
770
        #ifdef __LITTLE_ENDIAN__
 
771
                floatSwapBytes(floatArr);
 
772
                floatSwapBytes((floatArr+1));
 
773
                floatSwapBytes((floatArr+2));
 
774
                floatSwapBytes((floatArr+3));
 
775
        #endif
 
776
}
 
777
 
 
778
Point3D IonHit::getPos() const
 
779
{
 
780
        return pos;
 
781
}       
 
782
 
 
783
bool IonHit::hasNaN()
 
784
{
 
785
        return (std::isnan(massToCharge) || std::isnan(pos[0]) || 
 
786
                                std::isnan(pos[1]) || std::isnan(pos[2]));
 
787
}
 
788
 
 
789
 
 
790
 
 
791
 
 
792
 
 
793
void getPointSum(const std::vector<IonHit> &points,Point3D &centroid)
 
794
{
 
795
        //TODO: Paralellise me
 
796
        centroid=Point3D(0,0,0);
 
797
        for(unsigned int ui=0;ui<points.size();ui++)
 
798
                centroid+=points[ui].getPos();
 
799
}
 
800
 
 
801
BoundCube getIonDataLimits(const std::vector<IonHit> &points)
 
802
{
 
803
        ASSERT(points.size());
 
804
 
 
805
        BoundCube b;    
 
806
        b.setInverseLimits();   
 
807
#ifndef OPENMP
 
808
        float bounds[3][2];
 
809
        for(unsigned int ui=0;ui<3;ui++)
 
810
        {
 
811
                bounds[ui][0]=std::numeric_limits<float>::max();
 
812
                bounds[ui][1]=-std::numeric_limits<float>::max();
 
813
        }
 
814
        
 
815
        for(unsigned int ui=0; ui<points.size(); ui++)
 
816
        {
 
817
 
 
818
                for(unsigned int uj=0; uj<3; uj++)
 
819
                {
 
820
                        Point3D p;
 
821
                        p=points[ui].getPos();
 
822
                        if(p.getValue(uj) < bounds[uj][0])
 
823
                                bounds[uj][0] = p.getValue(uj);
 
824
                        
 
825
                        if(p.getValue(uj) > bounds[uj][1])
 
826
                                bounds[uj][1] = p.getValue(uj);
 
827
                }
 
828
        }
 
829
 
 
830
        b.setBounds(bounds[0][0],bounds[1][0],
 
831
                        bounds[2][0],bounds[0][1],
 
832
                        bounds[1][1],bounds[2][1]);
 
833
#else
 
834
        // parallel version
 
835
        vector<BoundCube> cubes;
 
836
 
 
837
        unsigned int nT=omp_get_max_threads();
 
838
        cubes.resize(nT);
 
839
        for(unsigned int ui=0;ui<cubes.size();ui++)
 
840
                cube[ui].setInverseLimits();
 
841
 
 
842
        unsigned int tCount=1;
 
843
        #pragma omp parallel for reduction(tCount|+)
 
844
        for(unsigned int ui=0;ui<points.size();ui++)
 
845
        {
 
846
                Point3D p;
 
847
                p=points[ui].getPos();
 
848
                for(unsigned int uj=0;uj<3;uj++)
 
849
                {
 
850
                        b.setBounds(uj,0,std::min(b.getBound(uj,0),p[uj]));
 
851
                        b.setBounds(uj,1,std::min(b.getBound(uj,0),p[uj]));
 
852
                }
 
853
        }
 
854
 
 
855
        for(unsigned int ui=0;ui<std::min(tCount,nT);ui++)
 
856
                b.expand(cubes[ui]);
 
857
 
 
858
#endif
 
859
 
 
860
        return b;
 
861
}